| Geant4 Cross Reference |
1 // 1
2 // *******************************************
3 // * License and Disclaimer
4 // *
5 // * The Geant4 software is copyright of th
6 // * the Geant4 Collaboration. It is provided
7 // * conditions of the Geant4 Software License
8 // * LICENSE and available at http://cern.ch/
9 // * include a list of copyright holders.
10 // *
11 // * Neither the authors of this software syst
12 // * institutes,nor the agencies providing fin
13 // * work make any representation or warran
14 // * regarding this software system or assum
15 // * use. Please see the license in the file
16 // * for the full disclaimer and the limitatio
17 // *
18 // * This code implementation is the result
19 // * technical work of the GEANT4 collaboratio
20 // * By using, copying, modifying or distri
21 // * any work based on the software) you ag
22 // * use in resulting scientific publicati
23 // * acceptance of all terms of the Geant4 Sof
24 // *******************************************
25 //
26 /// \file electromagnetic/TestEm5/src/Run.cc
27 /// \brief Implementation of the Run class
28 //
29 //
30 //....oooOO0OOooo........oooOO0OOooo........oo
31 //....oooOO0OOooo........oooOO0OOooo........oo
32
33 #include "Run.hh"
34
35 #include "DetectorConstruction.hh"
36 #include "HistoManager.hh"
37 #include "PrimaryGeneratorAction.hh"
38
39 #include "G4EmCalculator.hh"
40 #include "G4SystemOfUnits.hh"
41 #include "G4UnitsTable.hh"
42
43 #include <iomanip>
44
45 //....oooOO0OOooo........oooOO0OOooo........oo
46
47 Run::Run(DetectorConstruction* det) : fDetecto
48
49 //....oooOO0OOooo........oooOO0OOooo........oo
50
51 void Run::SetPrimary(G4ParticleDefinition* par
52 {
53 fParticle = particle;
54 fEkin = energy;
55
56 // compute theta0
57 fMscThetaCentral = 3 * ComputeMscHighland();
58 }
59
60 //....oooOO0OOooo........oooOO0OOooo........oo
61
62 void Run::Merge(const G4Run* run)
63 {
64 const Run* localRun = static_cast<const Run*
65
66 // pass information about primary particle
67 fParticle = localRun->fParticle;
68 fEkin = localRun->fEkin;
69
70 fMscThetaCentral = localRun->fMscThetaCentra
71
72 // accumulate sums
73 //
74 fEnergyDeposit += localRun->fEnergyDeposit;
75 fEnergyDeposit2 += localRun->fEnergyDeposit2
76 fTrakLenCharged += localRun->fTrakLenCharged
77 fTrakLenCharged2 += localRun->fTrakLenCharge
78 fTrakLenNeutral += localRun->fTrakLenNeutral
79 fTrakLenNeutral2 += localRun->fTrakLenNeutra
80 fNbStepsCharged += localRun->fNbStepsCharged
81 fNbStepsCharged2 += localRun->fNbStepsCharge
82 fNbStepsNeutral += localRun->fNbStepsNeutral
83 fNbStepsNeutral2 += localRun->fNbStepsNeutra
84 fMscProjecTheta += localRun->fMscProjecTheta
85 fMscProjecTheta2 += localRun->fMscProjecThet
86
87 fTypes[0] += localRun->fTypes[0];
88 fTypes[1] += localRun->fTypes[1];
89 fTypes[2] += localRun->fTypes[2];
90 fTypes[3] += localRun->fTypes[3];
91
92 fNbGamma += localRun->fNbGamma;
93 fNbElect += localRun->fNbElect;
94 fNbPosit += localRun->fNbPosit;
95
96 fTransmit[0] += localRun->fTransmit[0];
97 fTransmit[1] += localRun->fTransmit[1];
98 fReflect[0] += localRun->fReflect[0];
99 fReflect[1] += localRun->fReflect[1];
100
101 fMscEntryCentral += localRun->fMscEntryCentr
102
103 fEnergyLeak[0] += localRun->fEnergyLeak[0];
104 fEnergyLeak[1] += localRun->fEnergyLeak[1];
105 fEnergyLeak2[0] += localRun->fEnergyLeak2[0]
106 fEnergyLeak2[1] += localRun->fEnergyLeak2[1]
107
108 G4Run::Merge(run);
109 }
110
111 //....oooOO0OOooo........oooOO0OOooo........oo
112
113 void Run::EndOfRun()
114 {
115 // compute mean and rms
116 //
117 G4int TotNbofEvents = numberOfEvent;
118 if (TotNbofEvents == 0) return;
119
120 G4double EnergyBalance = fEnergyDeposit + fE
121 EnergyBalance /= TotNbofEvents;
122
123 fEnergyDeposit /= TotNbofEvents;
124 fEnergyDeposit2 /= TotNbofEvents;
125 G4double rmsEdep = fEnergyDeposit2 - fEnergy
126 if (rmsEdep > 0.)
127 rmsEdep = std::sqrt(rmsEdep / TotNbofEvent
128 else
129 rmsEdep = 0.;
130
131 fTrakLenCharged /= TotNbofEvents;
132 fTrakLenCharged2 /= TotNbofEvents;
133 G4double rmsTLCh = fTrakLenCharged2 - fTrakL
134 if (rmsTLCh > 0.)
135 rmsTLCh = std::sqrt(rmsTLCh / TotNbofEvent
136 else
137 rmsTLCh = 0.;
138
139 fTrakLenNeutral /= TotNbofEvents;
140 fTrakLenNeutral2 /= TotNbofEvents;
141 G4double rmsTLNe = fTrakLenNeutral2 - fTrakL
142 if (rmsTLNe > 0.)
143 rmsTLNe = std::sqrt(rmsTLNe / TotNbofEvent
144 else
145 rmsTLNe = 0.;
146
147 fNbStepsCharged /= TotNbofEvents;
148 fNbStepsCharged2 /= TotNbofEvents;
149 G4double rmsStCh = fNbStepsCharged2 - fNbSte
150 if (rmsStCh > 0.)
151 rmsStCh = std::sqrt(rmsStCh / TotNbofEvent
152 else
153 rmsStCh = 0.;
154
155 fNbStepsNeutral /= TotNbofEvents;
156 fNbStepsNeutral2 /= TotNbofEvents;
157 G4double rmsStNe = fNbStepsNeutral2 - fNbSte
158 if (rmsStNe > 0.)
159 rmsStNe = std::sqrt(rmsStNe / TotNbofEvent
160 else
161 rmsStNe = 0.;
162
163 G4double Gamma = (G4double)fNbGamma / TotNbo
164 G4double Elect = (G4double)fNbElect / TotNbo
165 G4double Posit = (G4double)fNbPosit / TotNbo
166
167 G4double transmit[2];
168 transmit[0] = 100. * fTransmit[0] / TotNbofE
169 transmit[1] = 100. * fTransmit[1] / TotNbofE
170
171 G4double reflect[2];
172 reflect[0] = 100. * fReflect[0] / TotNbofEve
173 reflect[1] = 100. * fReflect[1] / TotNbofEve
174
175 G4double rmsMsc = 0., tailMsc = 0.;
176 if (fMscEntryCentral > 0) {
177 fMscProjecTheta /= fMscEntryCentral;
178 fMscProjecTheta2 /= fMscEntryCentral;
179 rmsMsc = fMscProjecTheta2 - fMscProjecThet
180 if (rmsMsc > 0.) {
181 rmsMsc = std::sqrt(rmsMsc);
182 }
183 if (fTransmit[1] > 0.0) {
184 tailMsc = 100. - (100. * fMscEntryCentra
185 }
186 }
187
188 fEnergyLeak[0] /= TotNbofEvents;
189 fEnergyLeak2[0] /= TotNbofEvents;
190 G4double rmsEl0 = fEnergyLeak2[0] - fEnergyL
191 if (rmsEl0 > 0.)
192 rmsEl0 = std::sqrt(rmsEl0 / TotNbofEvents)
193 else
194 rmsEl0 = 0.;
195
196 fEnergyLeak[1] /= TotNbofEvents;
197 fEnergyLeak2[1] /= TotNbofEvents;
198 G4double rmsEl1 = fEnergyLeak2[1] - fEnergyL
199 if (rmsEl1 > 0.)
200 rmsEl1 = std::sqrt(rmsEl1 / TotNbofEvents)
201 else
202 rmsEl1 = 0.;
203
204 // Stopping Power from input Table.
205 //
206 const G4Material* material = fDetector->GetA
207 G4double length = fDetector->GetAbsorberThic
208 G4double density = material->GetDensity();
209 G4String partName = fParticle->GetParticleNa
210
211 G4EmCalculator emCalculator;
212 G4double dEdxTable = 0., dEdxFull = 0.;
213 if (fParticle->GetPDGCharge() != 0.) {
214 dEdxTable = emCalculator.GetDEDX(fEkin, fP
215 dEdxFull = emCalculator.ComputeTotalDEDX(f
216 }
217 G4double stopTable = dEdxTable / density;
218 G4double stopFull = dEdxFull / density;
219
220 // Stopping Power from simulation.
221 //
222 G4double meandEdx = fEnergyDeposit / length;
223 G4double stopPower = meandEdx / density;
224
225 G4cout << "\n ======================== run s
226
227 G4int prec = G4cout.precision(3);
228
229 G4cout << "\n The run was " << TotNbofEvents
230 << G4BestUnit(fEkin, "Energy") << " t
231 << material->GetName() << " (density:
232 << G4endl;
233
234 G4cout.precision(4);
235
236 G4cout << "\n Total energy deposit in absorb
237 << G4BestUnit(fEnergyDeposit, "Energy
238 << G4endl;
239
240 G4cout << "\n -----> Mean dE/dx = " << meand
241 << "\t(" << stopPower / (MeV * cm2 /
242
243 G4cout << "\n From formulas :" << G4endl;
244 G4cout << " restricted dEdx = " << dEdxTab
245 << "\t(" << stopTable / (MeV * cm2 /
246
247 G4cout << " full dEdx = " << dEdxFul
248 << "\t(" << stopFull / (MeV * cm2 / g
249
250 G4cout << "\n Leakage : primary = " << G4Be
251 << G4BestUnit(rmsEl0, "Energy")
252 << " secondaries = " << G4BestUnit(
253 << G4BestUnit(rmsEl1, "Energy") << G4
254
255 G4cout << " Energy balance : edep + eleak =
256
257 G4cout << "\n Total track length (charged) i
258 << G4BestUnit(fTrakLenCharged, "Lengt
259 << G4endl;
260
261 G4cout << " Total track length (neutral) in
262 << G4BestUnit(fTrakLenNeutral, "Lengt
263 << G4endl;
264
265 G4cout << "\n Number of steps (charged) in a
266 << rmsStCh << G4endl;
267
268 G4cout << " Number of steps (neutral) in abs
269 << rmsStNe << G4endl;
270
271 G4cout << "\n Number of secondaries per even
272 << "; positrons = " << Posit << G4e
273
274 G4cout << "\n Number of events with the prim
275 << G4endl;
276
277 G4cout << " Number of events with at least
278 << "(same charge as primary) = " << t
279
280 G4cout << "\n Number of events with the prim
281 << G4endl;
282
283 G4cout << " Number of events with at least
284 << "(same charge as primary) = " << r
285
286 // compute width of the Gaussian central par
287 //
288 G4cout << "\n MultipleScattering:"
289 << "\n rms proj angle of transmit pr
290 << " mrad (central part only)" << G4e
291
292 G4cout << " computed theta0 (Highland formu
293 << " mrad" << G4endl;
294
295 G4cout << " central part defined as +- " <<
296 << " Tail ratio = " << tailMsc << "
297
298 // gamma process counts
299 //
300 G4cout << "\n Gamma process counts:" << G4en
301 G4cout << " Photoeffect " << fTypes[0] <<
302 G4cout << " Compton " << fTypes[1] <<
303 G4cout << " Conversion " << fTypes[2] <<
304 G4cout << " Rayleigh " << fTypes[3] <<
305
306 // normalize histograms
307 //
308 G4AnalysisManager* analysisManager = G4Analy
309
310 G4int ih = 1;
311 G4double binWidth = analysisManager->GetH1Wi
312 G4double fac = 1. / (TotNbofEvents * binWidt
313 analysisManager->ScaleH1(ih, fac);
314
315 ih = 10;
316 binWidth = analysisManager->GetH1Width(ih);
317 fac = 1. / (TotNbofEvents * binWidth);
318 analysisManager->ScaleH1(ih, fac);
319
320 ih = 12;
321 analysisManager->ScaleH1(ih, 1. / TotNbofEve
322
323 // reset default precision
324 G4cout.precision(prec);
325 }
326
327 //....oooOO0OOooo........oooOO0OOooo........oo
328
329 G4double Run::ComputeMscHighland()
330 {
331 // compute the width of the Gaussian central
332 // projected angular distribution.
333 // Eur. Phys. Jour. C15 (2000) page 166, for
334
335 G4double t =
336 (fDetector->GetAbsorberThickness()) / (fDe
337 if (t < DBL_MIN) return 0.;
338
339 G4double T = fEkin;
340 G4double M = fParticle->GetPDGMass();
341 G4double z = std::abs(fParticle->GetPDGCharg
342
343 G4double bpc = T * (T + 2 * M) / (T + M);
344 G4double teta0 = 13.6 * MeV * z * std::sqrt(
345 return teta0;
346 }
347
348 //....oooOO0OOooo........oooOO0OOooo........oo
349