| Geant4 Cross Reference |
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 /// \file electromagnetic/TestEm5/src/Run.cc
27 /// \brief Implementation of the Run class
28 //
29 //
30 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
31 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
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........oooOO0OOooo........oooOO0OOooo......
46
47 Run::Run(DetectorConstruction* det) : fDetector(det) {}
48
49 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
50
51 void Run::SetPrimary(G4ParticleDefinition* particle, G4double energy)
52 {
53 fParticle = particle;
54 fEkin = energy;
55
56 // compute theta0
57 fMscThetaCentral = 3 * ComputeMscHighland();
58 }
59
60 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
61
62 void Run::Merge(const G4Run* run)
63 {
64 const Run* localRun = static_cast<const Run*>(run);
65
66 // pass information about primary particle
67 fParticle = localRun->fParticle;
68 fEkin = localRun->fEkin;
69
70 fMscThetaCentral = localRun->fMscThetaCentral;
71
72 // accumulate sums
73 //
74 fEnergyDeposit += localRun->fEnergyDeposit;
75 fEnergyDeposit2 += localRun->fEnergyDeposit2;
76 fTrakLenCharged += localRun->fTrakLenCharged;
77 fTrakLenCharged2 += localRun->fTrakLenCharged2;
78 fTrakLenNeutral += localRun->fTrakLenNeutral;
79 fTrakLenNeutral2 += localRun->fTrakLenNeutral2;
80 fNbStepsCharged += localRun->fNbStepsCharged;
81 fNbStepsCharged2 += localRun->fNbStepsCharged2;
82 fNbStepsNeutral += localRun->fNbStepsNeutral;
83 fNbStepsNeutral2 += localRun->fNbStepsNeutral2;
84 fMscProjecTheta += localRun->fMscProjecTheta;
85 fMscProjecTheta2 += localRun->fMscProjecTheta2;
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->fMscEntryCentral;
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........oooOO0OOooo........oooOO0OOooo......
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 + fEnergyLeak[0] + fEnergyLeak[1];
121 EnergyBalance /= TotNbofEvents;
122
123 fEnergyDeposit /= TotNbofEvents;
124 fEnergyDeposit2 /= TotNbofEvents;
125 G4double rmsEdep = fEnergyDeposit2 - fEnergyDeposit * fEnergyDeposit;
126 if (rmsEdep > 0.)
127 rmsEdep = std::sqrt(rmsEdep / TotNbofEvents);
128 else
129 rmsEdep = 0.;
130
131 fTrakLenCharged /= TotNbofEvents;
132 fTrakLenCharged2 /= TotNbofEvents;
133 G4double rmsTLCh = fTrakLenCharged2 - fTrakLenCharged * fTrakLenCharged;
134 if (rmsTLCh > 0.)
135 rmsTLCh = std::sqrt(rmsTLCh / TotNbofEvents);
136 else
137 rmsTLCh = 0.;
138
139 fTrakLenNeutral /= TotNbofEvents;
140 fTrakLenNeutral2 /= TotNbofEvents;
141 G4double rmsTLNe = fTrakLenNeutral2 - fTrakLenNeutral * fTrakLenNeutral;
142 if (rmsTLNe > 0.)
143 rmsTLNe = std::sqrt(rmsTLNe / TotNbofEvents);
144 else
145 rmsTLNe = 0.;
146
147 fNbStepsCharged /= TotNbofEvents;
148 fNbStepsCharged2 /= TotNbofEvents;
149 G4double rmsStCh = fNbStepsCharged2 - fNbStepsCharged * fNbStepsCharged;
150 if (rmsStCh > 0.)
151 rmsStCh = std::sqrt(rmsStCh / TotNbofEvents);
152 else
153 rmsStCh = 0.;
154
155 fNbStepsNeutral /= TotNbofEvents;
156 fNbStepsNeutral2 /= TotNbofEvents;
157 G4double rmsStNe = fNbStepsNeutral2 - fNbStepsNeutral * fNbStepsNeutral;
158 if (rmsStNe > 0.)
159 rmsStNe = std::sqrt(rmsStNe / TotNbofEvents);
160 else
161 rmsStNe = 0.;
162
163 G4double Gamma = (G4double)fNbGamma / TotNbofEvents;
164 G4double Elect = (G4double)fNbElect / TotNbofEvents;
165 G4double Posit = (G4double)fNbPosit / TotNbofEvents;
166
167 G4double transmit[2];
168 transmit[0] = 100. * fTransmit[0] / TotNbofEvents;
169 transmit[1] = 100. * fTransmit[1] / TotNbofEvents;
170
171 G4double reflect[2];
172 reflect[0] = 100. * fReflect[0] / TotNbofEvents;
173 reflect[1] = 100. * fReflect[1] / TotNbofEvents;
174
175 G4double rmsMsc = 0., tailMsc = 0.;
176 if (fMscEntryCentral > 0) {
177 fMscProjecTheta /= fMscEntryCentral;
178 fMscProjecTheta2 /= fMscEntryCentral;
179 rmsMsc = fMscProjecTheta2 - fMscProjecTheta * fMscProjecTheta;
180 if (rmsMsc > 0.) {
181 rmsMsc = std::sqrt(rmsMsc);
182 }
183 if (fTransmit[1] > 0.0) {
184 tailMsc = 100. - (100. * fMscEntryCentral) / (2 * fTransmit[1]);
185 }
186 }
187
188 fEnergyLeak[0] /= TotNbofEvents;
189 fEnergyLeak2[0] /= TotNbofEvents;
190 G4double rmsEl0 = fEnergyLeak2[0] - fEnergyLeak[0] * fEnergyLeak[0];
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] - fEnergyLeak[1] * fEnergyLeak[1];
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->GetAbsorberMaterial();
207 G4double length = fDetector->GetAbsorberThickness();
208 G4double density = material->GetDensity();
209 G4String partName = fParticle->GetParticleName();
210
211 G4EmCalculator emCalculator;
212 G4double dEdxTable = 0., dEdxFull = 0.;
213 if (fParticle->GetPDGCharge() != 0.) {
214 dEdxTable = emCalculator.GetDEDX(fEkin, fParticle, material);
215 dEdxFull = emCalculator.ComputeTotalDEDX(fEkin, fParticle, material);
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 summary ======================\n";
226
227 G4int prec = G4cout.precision(3);
228
229 G4cout << "\n The run was " << TotNbofEvents << " " << partName << " of "
230 << G4BestUnit(fEkin, "Energy") << " through " << G4BestUnit(length, "Length") << " of "
231 << material->GetName() << " (density: " << G4BestUnit(density, "Volumic Mass") << ")"
232 << G4endl;
233
234 G4cout.precision(4);
235
236 G4cout << "\n Total energy deposit in absorber per event = "
237 << G4BestUnit(fEnergyDeposit, "Energy") << " +- " << G4BestUnit(rmsEdep, "Energy")
238 << G4endl;
239
240 G4cout << "\n -----> Mean dE/dx = " << meandEdx / (MeV / cm) << " MeV/cm"
241 << "\t(" << stopPower / (MeV * cm2 / g) << " MeV*cm2/g)" << G4endl;
242
243 G4cout << "\n From formulas :" << G4endl;
244 G4cout << " restricted dEdx = " << dEdxTable / (MeV / cm) << " MeV/cm"
245 << "\t(" << stopTable / (MeV * cm2 / g) << " MeV*cm2/g)" << G4endl;
246
247 G4cout << " full dEdx = " << dEdxFull / (MeV / cm) << " MeV/cm"
248 << "\t(" << stopFull / (MeV * cm2 / g) << " MeV*cm2/g)" << G4endl;
249
250 G4cout << "\n Leakage : primary = " << G4BestUnit(fEnergyLeak[0], "Energy") << " +- "
251 << G4BestUnit(rmsEl0, "Energy")
252 << " secondaries = " << G4BestUnit(fEnergyLeak[1], "Energy") << " +- "
253 << G4BestUnit(rmsEl1, "Energy") << G4endl;
254
255 G4cout << " Energy balance : edep + eleak = " << G4BestUnit(EnergyBalance, "Energy") << G4endl;
256
257 G4cout << "\n Total track length (charged) in absorber per event = "
258 << G4BestUnit(fTrakLenCharged, "Length") << " +- " << G4BestUnit(rmsTLCh, "Length")
259 << G4endl;
260
261 G4cout << " Total track length (neutral) in absorber per event = "
262 << G4BestUnit(fTrakLenNeutral, "Length") << " +- " << G4BestUnit(rmsTLNe, "Length")
263 << G4endl;
264
265 G4cout << "\n Number of steps (charged) in absorber per event = " << fNbStepsCharged << " +- "
266 << rmsStCh << G4endl;
267
268 G4cout << " Number of steps (neutral) in absorber per event = " << fNbStepsNeutral << " +- "
269 << rmsStNe << G4endl;
270
271 G4cout << "\n Number of secondaries per event : Gammas = " << Gamma << "; electrons = " << Elect
272 << "; positrons = " << Posit << G4endl;
273
274 G4cout << "\n Number of events with the primary particle transmitted = " << transmit[1] << " %"
275 << G4endl;
276
277 G4cout << " Number of events with at least 1 particle transmitted "
278 << "(same charge as primary) = " << transmit[0] << " %" << G4endl;
279
280 G4cout << "\n Number of events with the primary particle reflected = " << reflect[1] << " %"
281 << G4endl;
282
283 G4cout << " Number of events with at least 1 particle reflected "
284 << "(same charge as primary) = " << reflect[0] << " %" << G4endl;
285
286 // compute width of the Gaussian central part of the MultipleScattering
287 //
288 G4cout << "\n MultipleScattering:"
289 << "\n rms proj angle of transmit primary particle = " << rmsMsc / mrad
290 << " mrad (central part only)" << G4endl;
291
292 G4cout << " computed theta0 (Highland formula) = " << ComputeMscHighland() / mrad
293 << " mrad" << G4endl;
294
295 G4cout << " central part defined as +- " << fMscThetaCentral / mrad << " mrad; "
296 << " Tail ratio = " << tailMsc << " %" << G4endl;
297
298 // gamma process counts
299 //
300 G4cout << "\n Gamma process counts:" << G4endl;
301 G4cout << " Photoeffect " << fTypes[0] << G4endl;
302 G4cout << " Compton " << fTypes[1] << G4endl;
303 G4cout << " Conversion " << fTypes[2] << G4endl;
304 G4cout << " Rayleigh " << fTypes[3] << G4endl;
305
306 // normalize histograms
307 //
308 G4AnalysisManager* analysisManager = G4AnalysisManager::Instance();
309
310 G4int ih = 1;
311 G4double binWidth = analysisManager->GetH1Width(ih);
312 G4double fac = 1. / (TotNbofEvents * binWidth);
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. / TotNbofEvents);
322
323 // reset default precision
324 G4cout.precision(prec);
325 }
326
327 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
328
329 G4double Run::ComputeMscHighland()
330 {
331 // compute the width of the Gaussian central part of the MultipleScattering
332 // projected angular distribution.
333 // Eur. Phys. Jour. C15 (2000) page 166, formule 23.9
334
335 G4double t =
336 (fDetector->GetAbsorberThickness()) / (fDetector->GetAbsorberMaterial()->GetRadlen());
337 if (t < DBL_MIN) return 0.;
338
339 G4double T = fEkin;
340 G4double M = fParticle->GetPDGMass();
341 G4double z = std::abs(fParticle->GetPDGCharge() / eplus);
342
343 G4double bpc = T * (T + 2 * M) / (T + M);
344 G4double teta0 = 13.6 * MeV * z * std::sqrt(t) * (1. + 0.038 * std::log(t)) / bpc;
345 return teta0;
346 }
347
348 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
349