source: trunk/examples/extended/electromagnetic/TestEm5/src/RunAction.cc@ 1303

Last change on this file since 1303 was 1230, checked in by garnier, 16 years ago

update to geant4.9.3
File size: 11.6 KB
Line 
1//
2// ********************************************************************
3// * License and Disclaimer *
4// * *
5// * The Geant4 software is copyright of the Copyright Holders of *
6// * the Geant4 Collaboration. It is provided under the terms and *
7// * conditions of the Geant4 Software License, included in the file *
8// * LICENSE and available at http://cern.ch/geant4/license . These *
9// * include a list of copyright holders. *
10// * *
11// * Neither the authors of this software system, nor their employing *
12// * institutes,nor the agencies providing financial support for this *
13// * work make any representation or warranty, express or implied, *
14// * regarding this software system or assume any liability for its *
15// * use. Please see the license in the file LICENSE and URL above *
16// * for the full disclaimer and the limitation of liability. *
17// * *
18// * This code implementation is the result of the scientific and *
19// * technical work of the GEANT4 collaboration. *
20// * By using, copying, modifying or distributing the software (or *
21// * any work based on the software) you agree to acknowledge its *
22// * use in resulting scientific publications, and indicate your *
23// * acceptance of all terms of the Geant4 Software license. *
24// ********************************************************************
25//
26// $Id: RunAction.cc,v 1.29 2009/06/18 19:08:18 vnivanch Exp $
27// GEANT4 tag $Name: geant4-09-03-cand-01 $
28//
29//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
30//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
31
32#include "RunAction.hh"
33#include "DetectorConstruction.hh"
34#include "PrimaryGeneratorAction.hh"
35#include "HistoManager.hh"
36
37#include "G4Run.hh"
38#include "G4RunManager.hh"
39#include "G4UnitsTable.hh"
40#include "G4EmCalculator.hh"
41
42#include "Randomize.hh"
43#include <iomanip>
44
45//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
46
47RunAction::RunAction(DetectorConstruction* det, PrimaryGeneratorAction* kin,
48 HistoManager* histo)
49:detector(det), primary(kin), histoManager(histo)
50{ }
51
52//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
53
54RunAction::~RunAction()
55{ }
56
57//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
58
59void RunAction::BeginOfRunAction(const G4Run* aRun)
60{
61 G4cout << "### Run " << aRun->GetRunID() << " start." << G4endl;
62
63 //initialisation
64 EnergyDeposit = EnergyDeposit2 = 0.;
65 TrakLenCharged = TrakLenCharged2 = 0.;
66 TrakLenNeutral = TrakLenNeutral2 = 0.;
67 nbStepsCharged = nbStepsCharged2 = 0.;
68 nbStepsNeutral = nbStepsNeutral2 = 0.;
69 MscProjecTheta = MscProjecTheta2 = 0.;
70 MscThetaCentral = 3*ComputeMscHighland();
71
72 nbGamma = nbElect = nbPosit = 0;
73
74 Transmit[0] = Transmit[1] = Reflect[0] = Reflect[1] = 0;
75
76 MscEntryCentral = 0;
77
78 EnergyLeak[0] = EnergyLeak[1] = EnergyLeak2[0] = EnergyLeak2[1] = 0.;
79
80 histoManager->book();
81
82 // save Rndm status
83 G4RunManager::GetRunManager()->SetRandomNumberStore(true);
84 CLHEP::HepRandom::showEngineStatus();
85}
86
87//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
88
89void RunAction::EndOfRunAction(const G4Run* aRun)
90{
91 // compute mean and rms
92 //
93 G4int TotNbofEvents = aRun->GetNumberOfEvent();
94 if (TotNbofEvents == 0) return;
95
96 G4double EnergyBalance = EnergyDeposit + EnergyLeak[0] + EnergyLeak[1];
97 EnergyBalance /= TotNbofEvents;
98
99 EnergyDeposit /= TotNbofEvents; EnergyDeposit2 /= TotNbofEvents;
100 G4double rmsEdep = EnergyDeposit2 - EnergyDeposit*EnergyDeposit;
101 if (rmsEdep>0.) rmsEdep = std::sqrt(rmsEdep/TotNbofEvents);
102 else rmsEdep = 0.;
103
104 TrakLenCharged /= TotNbofEvents; TrakLenCharged2 /= TotNbofEvents;
105 G4double rmsTLCh = TrakLenCharged2 - TrakLenCharged*TrakLenCharged;
106 if (rmsTLCh>0.) rmsTLCh = std::sqrt(rmsTLCh/TotNbofEvents);
107 else rmsTLCh = 0.;
108
109 TrakLenNeutral /= TotNbofEvents; TrakLenNeutral2 /= TotNbofEvents;
110 G4double rmsTLNe = TrakLenNeutral2 - TrakLenNeutral*TrakLenNeutral;
111 if (rmsTLNe>0.) rmsTLNe = std::sqrt(rmsTLNe/TotNbofEvents);
112 else rmsTLNe = 0.;
113
114 nbStepsCharged /= TotNbofEvents; nbStepsCharged2 /= TotNbofEvents;
115 G4double rmsStCh = nbStepsCharged2 - nbStepsCharged*nbStepsCharged;
116 if (rmsStCh>0.) rmsStCh = std::sqrt(rmsTLCh/TotNbofEvents);
117 else rmsStCh = 0.;
118
119 nbStepsNeutral /= TotNbofEvents; nbStepsNeutral2 /= TotNbofEvents;
120 G4double rmsStNe = nbStepsNeutral2 - nbStepsNeutral*nbStepsNeutral;
121 if (rmsStNe>0.) rmsStNe = std::sqrt(rmsTLCh/TotNbofEvents);
122 else rmsStNe = 0.;
123
124 G4double Gamma = (G4double)nbGamma/TotNbofEvents;
125 G4double Elect = (G4double)nbElect/TotNbofEvents;
126 G4double Posit = (G4double)nbPosit/TotNbofEvents;
127
128 G4double transmit[2];
129 transmit[0] = 100.*Transmit[0]/TotNbofEvents;
130 transmit[1] = 100.*Transmit[1]/TotNbofEvents;
131
132 G4double reflect[2];
133 reflect[0] = 100.*Reflect[0]/TotNbofEvents;
134 reflect[1] = 100.*Reflect[1]/TotNbofEvents;
135
136 G4double rmsMsc = 0., tailMsc = 0.;
137 if (MscEntryCentral > 0) {
138 MscProjecTheta /= MscEntryCentral; MscProjecTheta2 /= MscEntryCentral;
139 rmsMsc = MscProjecTheta2 - MscProjecTheta*MscProjecTheta;
140 if (rmsMsc > 0.) rmsMsc = std::sqrt(rmsMsc);
141 tailMsc = 100.- (100.*MscEntryCentral)/(2*Transmit[1]);
142 }
143
144 EnergyLeak[0] /= TotNbofEvents; EnergyLeak2[0] /= TotNbofEvents;
145 G4double rmsEl0 = EnergyLeak2[0] - EnergyLeak[0]*EnergyLeak[0];
146 if (rmsEl0>0.) rmsEl0 = std::sqrt(rmsEl0/TotNbofEvents);
147 else rmsEl0 = 0.;
148
149 EnergyLeak[1] /= TotNbofEvents; EnergyLeak2[1] /= TotNbofEvents;
150 G4double rmsEl1 = EnergyLeak2[1] - EnergyLeak[1]*EnergyLeak[1];
151 if (rmsEl1>0.) rmsEl1 = std::sqrt(rmsEl1/TotNbofEvents);
152 else rmsEl1 = 0.;
153
154
155 //Stopping Power from input Table.
156 //
157 G4Material* material = detector->GetAbsorberMaterial();
158 G4double length = detector->GetAbsorberThickness();
159 G4double density = material->GetDensity();
160
161 G4ParticleDefinition* particle = primary->GetParticleGun()
162 ->GetParticleDefinition();
163 G4String partName = particle->GetParticleName();
164 G4double energy = primary->GetParticleGun()->GetParticleEnergy();
165
166 G4EmCalculator emCalculator;
167 G4double dEdxTable = 0., dEdxFull = 0.;
168 if (particle->GetPDGCharge()!= 0.) {
169 dEdxTable = emCalculator.GetDEDX(energy,particle,material);
170 dEdxFull = emCalculator.ComputeTotalDEDX(energy,particle,material);
171 }
172 G4double stopTable = dEdxTable/density;
173 G4double stopFull = dEdxFull /density;
174
175 //Stopping Power from simulation.
176 //
177 G4double meandEdx = EnergyDeposit/length;
178 G4double stopPower = meandEdx/density;
179
180 G4cout << "\n ======================== run summary ======================\n";
181
182 G4int prec = G4cout.precision(3);
183
184 G4cout << "\n The run was " << TotNbofEvents << " " << partName << " of "
185 << G4BestUnit(energy,"Energy") << " through "
186 << G4BestUnit(length,"Length") << " of "
187 << material->GetName() << " (density: "
188 << G4BestUnit(density,"Volumic Mass") << ")" << G4endl;
189
190 G4cout.precision(4);
191
192 G4cout << "\n Total energy deposit in absorber per event = "
193 << G4BestUnit(EnergyDeposit,"Energy") << " +- "
194 << G4BestUnit(rmsEdep, "Energy")
195 << G4endl;
196
197 G4cout << "\n -----> Mean dE/dx = " << meandEdx/(MeV/cm) << " MeV/cm"
198 << "\t(" << stopPower/(MeV*cm2/g) << " MeV*cm2/g)"
199 << G4endl;
200
201 G4cout << "\n From formulas :" << G4endl;
202 G4cout << " restricted dEdx = " << dEdxTable/(MeV/cm) << " MeV/cm"
203 << "\t(" << stopTable/(MeV*cm2/g) << " MeV*cm2/g)"
204 << G4endl;
205
206 G4cout << " full dEdx = " << dEdxFull/(MeV/cm) << " MeV/cm"
207 << "\t(" << stopFull/(MeV*cm2/g) << " MeV*cm2/g)"
208 << G4endl;
209
210 G4cout << "\n Leakage : primary = "
211 << G4BestUnit(EnergyLeak[0],"Energy") << " +- "
212 << G4BestUnit(rmsEl0, "Energy")
213 << " secondaries = "
214 << G4BestUnit(EnergyLeak[1],"Energy") << " +- "
215 << G4BestUnit(rmsEl1, "Energy")
216 << G4endl;
217
218 G4cout << " Energy balance : edep + eleak = "
219 << G4BestUnit(EnergyBalance,"Energy")
220 << G4endl;
221
222 G4cout << "\n Total track length (charged) in absorber per event = "
223 << G4BestUnit(TrakLenCharged,"Length") << " +- "
224 << G4BestUnit(rmsTLCh, "Length") << G4endl;
225
226 G4cout << " Total track length (neutral) in absorber per event = "
227 << G4BestUnit(TrakLenNeutral,"Length") << " +- "
228 << G4BestUnit(rmsTLNe, "Length") << G4endl;
229
230 G4cout << "\n Number of steps (charged) in absorber per event = "
231 << nbStepsCharged << " +- " << rmsStCh << G4endl;
232
233 G4cout << " Number of steps (neutral) in absorber per event = "
234 << nbStepsNeutral << " +- " << rmsStNe << G4endl;
235
236 G4cout << "\n Number of secondaries per event : Gammas = " << Gamma
237 << "; electrons = " << Elect
238 << "; positrons = " << Posit << G4endl;
239
240 G4cout << "\n Number of events with the primary particle transmitted = "
241 << transmit[1] << " %" << G4endl;
242
243 G4cout << " Number of events with at least 1 particle transmitted "
244 << "(same charge as primary) = " << transmit[0] << " %" << G4endl;
245
246 G4cout << "\n Number of events with the primary particle reflected = "
247 << reflect[1] << " %" << G4endl;
248
249 G4cout << " Number of events with at least 1 particle reflected "
250 << "(same charge as primary) = " << reflect[0] << " %" << G4endl;
251
252 // compute width of the Gaussian central part of the MultipleScattering
253 //
254 G4cout << "\n MultipleScattering:"
255 << "\n rms proj angle of transmit primary particle = "
256 << rmsMsc/mrad << " mrad (central part only)" << G4endl;
257
258 G4cout << " computed theta0 (Highland formula) = "
259 << ComputeMscHighland()/mrad << " mrad" << G4endl;
260
261 G4cout << " central part defined as +- "
262 << MscThetaCentral/mrad << " mrad; "
263 << " Tail ratio = " << tailMsc << " %" << G4endl;
264
265 G4cout.precision(prec);
266
267 // normalize histograms
268 //
269 G4int ih = 1;
270 G4double binWidth = histoManager->GetBinWidth(ih);
271 G4double unit = histoManager->GetHistoUnit(ih);
272 G4double fac = unit/(TotNbofEvents*binWidth);
273 histoManager->Scale(ih,fac);
274
275 ih = 10;
276 binWidth = histoManager->GetBinWidth(ih);
277 unit = histoManager->GetHistoUnit(ih);
278 fac = unit/(TotNbofEvents*binWidth);
279 histoManager->Scale(ih,fac);
280
281 // save histograms
282 histoManager->save();
283
284 // show Rndm status
285 CLHEP::HepRandom::showEngineStatus();
286}
287
288//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
289
290G4double RunAction::ComputeMscHighland()
291{
292 //compute the width of the Gaussian central part of the MultipleScattering
293 //projected angular distribution.
294 //Eur. Phys. Jour. C15 (2000) page 166, formule 23.9
295
296 G4double t = (detector->GetAbsorberThickness())
297 /(detector->GetAbsorberMaterial()->GetRadlen());
298 if (t < DBL_MIN) return 0.;
299
300 G4ParticleGun* particle = primary->GetParticleGun();
301 G4double T = particle->GetParticleEnergy();
302 G4double M = particle->GetParticleDefinition()->GetPDGMass();
303 G4double z = std::abs(particle->GetParticleDefinition()->GetPDGCharge()/eplus);
304
305 G4double bpc = T*(T+2*M)/(T+M);
306 G4double teta0 = 13.6*MeV*z*std::sqrt(t)*(1.+0.038*std::log(t))/bpc;
307 return teta0;
308}
309
310//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
Note: See TracBrowser for help on using the repository browser.