| 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/TestEm2/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 "EmAcceptance.hh"
36 #include "PrimaryGeneratorAction.hh"
37
38 #include "G4Run.hh"
39 #include "G4SystemOfUnits.hh"
40 #include "G4UnitsTable.hh"
41
42 #include <iomanip>
43
44 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
45
46 Run::Run(DetectorConstruction* det, PrimaryGeneratorAction* kin) : fDet(det), fKin(kin)
47 {
48 Reset();
49 }
50
51 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
52
53 void Run::Reset()
54 {
55 f_nLbin = fDet->GetnLtot();
56 f_dEdL.resize(f_nLbin);
57 fSumELongit.resize(f_nLbin);
58 fSumELongitCumul.resize(f_nLbin);
59 fSumE2Longit.resize(f_nLbin);
60 fSumE2LongitCumul.resize(f_nLbin);
61
62 f_nRbin = fDet->GetnRtot();
63 f_dEdR.resize(f_nRbin);
64 fSumERadial.resize(f_nRbin);
65 fSumERadialCumul.resize(f_nRbin);
66 fSumE2Radial.resize(f_nRbin);
67 fSumE2RadialCumul.resize(f_nRbin);
68
69 fChargedStep = 0.0;
70 fNeutralStep = 0.0;
71
72 fVerbose = 0;
73
74 // initialize arrays of cumulative energy deposition
75 //
76 for (G4int i = 0; i < f_nLbin; ++i) {
77 fSumELongit[i] = fSumE2Longit[i] = fSumELongitCumul[i] = fSumE2LongitCumul[i] = 0.;
78 }
79 for (G4int j = 0; j < f_nRbin; ++j) {
80 fSumERadial[j] = fSumE2Radial[j] = fSumERadialCumul[j] = fSumE2RadialCumul[j] = 0.;
81 }
82 // initialize track length
83 fSumChargTrLength = fSum2ChargTrLength = fSumNeutrTrLength = fSum2NeutrTrLength = 0.;
84 }
85
86 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
87
88 void Run::InitializePerEvent()
89 {
90 // initialize arrays of energy deposit per bin
91 for (G4int i = 0; i < f_nLbin; ++i) {
92 f_dEdL[i] = 0.;
93 }
94
95 for (G4int j = 0; j < f_nRbin; ++j) {
96 f_dEdR[j] = 0.;
97 }
98
99 // initialize tracklength
100 fChargTrLength = fNeutrTrLength = 0.;
101 }
102
103 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
104
105 void Run::FillPerEvent()
106 {
107 // accumulate statistic
108 //
109 G4double dLCumul = 0.;
110 for (G4int i = 0; i < f_nLbin; ++i) {
111 fSumELongit[i] += f_dEdL[i];
112 fSumE2Longit[i] += f_dEdL[i] * f_dEdL[i];
113 dLCumul += f_dEdL[i];
114 fSumELongitCumul[i] += dLCumul;
115 fSumE2LongitCumul[i] += dLCumul * dLCumul;
116 }
117
118 G4double dRCumul = 0.;
119 for (G4int j = 0; j < f_nRbin; j++) {
120 fSumERadial[j] += f_dEdR[j];
121 fSumE2Radial[j] += f_dEdR[j] * f_dEdR[j];
122 dRCumul += f_dEdR[j];
123 fSumERadialCumul[j] += dRCumul;
124 fSumE2RadialCumul[j] += dRCumul * dRCumul;
125 }
126
127 fSumChargTrLength += fChargTrLength;
128 fSum2ChargTrLength += fChargTrLength * fChargTrLength;
129 fSumNeutrTrLength += fNeutrTrLength;
130 fSum2NeutrTrLength += fNeutrTrLength * fNeutrTrLength;
131
132 // fill histograms
133 //
134
135 G4double Ekin = fKin->GetParticleGun()->GetParticleEnergy();
136 G4double mass = fKin->GetParticleGun()->GetParticleDefinition()->GetPDGMass();
137 G4double radl = fDet->GetMaterial()->GetRadlen();
138
139 G4AnalysisManager* analysisManager = G4AnalysisManager::Instance();
140 analysisManager->FillH1(1, 100. * dLCumul / (Ekin + mass));
141 analysisManager->FillH1(2, fChargTrLength / radl);
142 analysisManager->FillH1(3, fNeutrTrLength / radl);
143
144 // profiles
145 G4double norm = 100. / (Ekin + mass);
146 G4double dLradl = fDet->GetdLradl();
147 for (G4int i = 0; i < f_nLbin; ++i) {
148 G4double bin = (i + 0.5) * dLradl;
149 analysisManager->FillP1(0, bin, norm * f_dEdL[i] / dLradl);
150 }
151 G4double dRradl = fDet->GetdRradl();
152 for (G4int j = 0; j < f_nRbin; ++j) {
153 G4double bin = (j + 0.5) * dRradl;
154 analysisManager->FillP1(1, bin, norm * f_dEdR[j] / dRradl);
155 }
156 }
157
158 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
159
160 void Run::Merge(const G4Run* run)
161 {
162 const Run* localRun = static_cast<const Run*>(run);
163
164 fChargedStep += localRun->fChargedStep;
165 fNeutralStep += localRun->fNeutralStep;
166
167 for (G4int i = 0; i < f_nLbin; ++i) {
168 fSumELongit[i] += localRun->fSumELongit[i];
169 fSumE2Longit[i] += localRun->fSumE2Longit[i];
170 fSumELongitCumul[i] += localRun->fSumELongitCumul[i];
171 fSumE2LongitCumul[i] += localRun->fSumE2LongitCumul[i];
172 }
173
174 for (G4int j = 0; j < f_nRbin; ++j) {
175 fSumERadial[j] += localRun->fSumERadial[j];
176 fSumE2Radial[j] += localRun->fSumE2Radial[j];
177 fSumERadialCumul[j] += localRun->fSumERadialCumul[j];
178 fSumE2RadialCumul[j] += localRun->fSumE2RadialCumul[j];
179 }
180
181 fSumChargTrLength += localRun->fSumChargTrLength;
182 fSum2ChargTrLength += localRun->fSum2ChargTrLength;
183 fSumNeutrTrLength += localRun->fSumNeutrTrLength;
184 fSum2NeutrTrLength += localRun->fSum2NeutrTrLength;
185
186 G4Run::Merge(run);
187 }
188
189 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
190
191 void Run::EndOfRun(G4double edep, G4double rms, G4double& limit)
192 {
193 G4int NbOfEvents = GetNumberOfEvent();
194
195 G4double kinEnergy = fKin->GetParticleGun()->GetParticleEnergy();
196 assert(NbOfEvents * kinEnergy > 0);
197
198 fChargedStep /= G4double(NbOfEvents);
199 fNeutralStep /= G4double(NbOfEvents);
200
201 G4double mass = fKin->GetParticleGun()->GetParticleDefinition()->GetPDGMass();
202 G4double norme = 100. / (NbOfEvents * (kinEnergy + mass));
203
204 // longitudinal
205 //
206 G4double dLradl = fDet->GetdLradl();
207
208 MyVector MeanELongit(f_nLbin), rmsELongit(f_nLbin);
209 MyVector MeanELongitCumul(f_nLbin), rmsELongitCumul(f_nLbin);
210
211 G4AnalysisManager* analysisManager = G4AnalysisManager::Instance();
212
213 G4int i;
214 for (i = 0; i < f_nLbin; ++i) {
215 MeanELongit[i] = norme * fSumELongit[i];
216 rmsELongit[i] =
217 norme * std::sqrt(std::abs(NbOfEvents * fSumE2Longit[i] - fSumELongit[i] * fSumELongit[i]));
218
219 MeanELongitCumul[i] = norme * fSumELongitCumul[i];
220 rmsELongitCumul[i] = norme
221 * std::sqrt(std::abs(NbOfEvents * fSumE2LongitCumul[i]
222 - fSumELongitCumul[i] * fSumELongitCumul[i]));
223 G4double bin = (i + 0.5) * dLradl;
224 analysisManager->FillH1(4, bin, MeanELongit[i] / dLradl);
225 analysisManager->FillH1(5, bin, rmsELongit[i] / dLradl);
226 bin = (i + 1) * dLradl;
227 analysisManager->FillH1(6, bin, MeanELongitCumul[i]);
228 analysisManager->FillH1(7, bin, rmsELongitCumul[i]);
229 }
230
231 // radial
232 //
233 G4double dRradl = fDet->GetdRradl();
234
235 MyVector MeanERadial(f_nRbin), rmsERadial(f_nRbin);
236 MyVector MeanERadialCumul(f_nRbin), rmsERadialCumul(f_nRbin);
237
238 for (i = 0; i < f_nRbin; ++i) {
239 MeanERadial[i] = norme * fSumERadial[i];
240 rmsERadial[i] =
241 norme * std::sqrt(std::abs(NbOfEvents * fSumE2Radial[i] - fSumERadial[i] * fSumERadial[i]));
242
243 MeanERadialCumul[i] = norme * fSumERadialCumul[i];
244 rmsERadialCumul[i] = norme
245 * std::sqrt(std::abs(NbOfEvents * fSumE2RadialCumul[i]
246 - fSumERadialCumul[i] * fSumERadialCumul[i]));
247
248 G4double bin = (i + 0.5) * dRradl;
249 analysisManager->FillH1(8, bin, MeanERadial[i] / dRradl);
250 analysisManager->FillH1(9, bin, rmsERadial[i] / dRradl);
251 bin = (i + 1) * dRradl;
252 analysisManager->FillH1(10, bin, MeanERadialCumul[i]);
253 analysisManager->FillH1(11, bin, rmsERadialCumul[i]);
254 }
255
256 // find Moliere confinement
257 //
258 const G4double EMoliere = 90.;
259 G4double iMoliere = 0.;
260 if ((MeanERadialCumul[0] <= EMoliere) && (MeanERadialCumul[f_nRbin - 1] >= EMoliere)) {
261 G4int imin = 0;
262 while ((imin < f_nRbin - 1) && (MeanERadialCumul[imin] < EMoliere)) {
263 ++imin;
264 }
265
266 G4double del = MeanERadialCumul[imin + 1] - MeanERadialCumul[imin];
267 G4double ratio = (del > 0.0) ? (EMoliere - MeanERadialCumul[imin]) / del : 0.0;
268 iMoliere = 1. + imin + ratio;
269 }
270
271 // track length
272 //
273 norme = 1. / (NbOfEvents * (fDet->GetMaterial()->GetRadlen()));
274 G4double MeanChargTrLength = norme * fSumChargTrLength;
275 G4double rmsChargTrLength =
276 norme
277 * std::sqrt(std::abs(NbOfEvents * fSum2ChargTrLength - fSumChargTrLength * fSumChargTrLength));
278
279 G4double MeanNeutrTrLength = norme * fSumNeutrTrLength;
280 G4double rmsNeutrTrLength =
281 norme
282 * std::sqrt(std::abs(NbOfEvents * fSum2NeutrTrLength - fSumNeutrTrLength * fSumNeutrTrLength));
283
284 // print
285 std::ios::fmtflags mode = G4cout.flags();
286 G4cout.setf(std::ios::fixed, std::ios::floatfield);
287 G4int prec = G4cout.precision(2);
288
289 if (fVerbose) {
290 G4cout << " LOGITUDINAL PROFILE "
291 << " CUMULATIVE LOGITUDINAL PROFILE" << G4endl << G4endl;
292
293 G4cout << " bin "
294 << " Mean rms "
295 << " bin "
296 << " Mean rms \n"
297 << G4endl;
298
299 for (i = 0; i < f_nLbin; ++i) {
300 G4double inf = i * dLradl, sup = inf + dLradl;
301
302 G4cout << std::setw(8) << inf << "->" << std::setw(5) << sup << " radl: " << std::setw(7)
303 << MeanELongit[i] << "% " << std::setw(9) << rmsELongit[i] << "% "
304 << " 0->" << std::setw(5) << sup << " radl: " << std::setw(7)
305 << MeanELongitCumul[i] << "% " << std::setw(7) << rmsELongitCumul[i] << "% "
306 << G4endl;
307 }
308
309 G4cout << G4endl << G4endl << G4endl;
310
311 G4cout << " RADIAL PROFILE "
312 << " CUMULATIVE RADIAL PROFILE" << G4endl << G4endl;
313
314 G4cout << " bin "
315 << " Mean rms "
316 << " bin "
317 << " Mean rms \n"
318 << G4endl;
319
320 for (i = 0; i < f_nRbin; ++i) {
321 G4double inf = i * dRradl, sup = inf + dRradl;
322
323 G4cout << std::setw(8) << inf << "->" << std::setw(5) << sup << " radl: " << std::setw(7)
324 << MeanERadial[i] << "% " << std::setw(9) << rmsERadial[i] << "% "
325 << " 0->" << std::setw(5) << sup << " radl: " << std::setw(7)
326 << MeanERadialCumul[i] << "% " << std::setw(7) << rmsERadialCumul[i] << "% "
327 << G4endl;
328 }
329 }
330
331 G4cout << "\n ===== SUMMARY ===== \n" << G4endl;
332
333 G4cout << " Total number of events: " << NbOfEvents << "\n"
334 << " Mean number of charged steps: " << fChargedStep << G4endl;
335 G4cout << " Mean number of neutral steps: " << fNeutralStep << "\n" << G4endl;
336
337 G4cout << " energy deposit : " << std::setw(7) << MeanELongitCumul[f_nLbin - 1] << " % E0 +- "
338 << std::setw(7) << rmsELongitCumul[f_nLbin - 1] << " % E0" << G4endl;
339 G4cout << " charged traklen: " << std::setw(7) << MeanChargTrLength << " radl +- " << std::setw(7)
340 << rmsChargTrLength << " radl" << G4endl;
341 G4cout << " neutral traklen: " << std::setw(7) << MeanNeutrTrLength << " radl +- " << std::setw(7)
342 << rmsNeutrTrLength << " radl" << G4endl;
343
344 if (iMoliere > 0.) {
345 G4double RMoliere1 = iMoliere * fDet->GetdRradl();
346 G4double RMoliere2 = iMoliere * fDet->GetdRlength();
347 G4cout << "\n " << EMoliere << " % confinement: radius = " << RMoliere1 << " radl ("
348 << G4BestUnit(RMoliere2, "Length") << ")"
349 << "\n"
350 << G4endl;
351 }
352
353 G4cout.setf(mode, std::ios::floatfield);
354 G4cout.precision(prec);
355
356 // Acceptance
357
358 G4int nLbin = fDet->GetnLtot();
359 if (limit < DBL_MAX) {
360 EmAcceptance acc;
361 acc.BeginOfAcceptance("Total Energy in Absorber", NbOfEvents);
362 G4double e = MeanELongitCumul[nLbin - 1] / 100.;
363 G4double r = rmsELongitCumul[nLbin - 1] / 100.;
364 acc.EmAcceptanceGauss("Edep", NbOfEvents, e, edep, rms, limit);
365 acc.EmAcceptanceGauss("Erms", NbOfEvents, r, rms, rms, 2.0 * limit);
366 acc.EndOfAcceptance();
367 }
368 limit = DBL_MAX;
369 }
370
371 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
372