| 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/TestEm2/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 "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........oo
45
46 Run::Run(DetectorConstruction* det, PrimaryGen
47 {
48 Reset();
49 }
50
51 //....oooOO0OOooo........oooOO0OOooo........oo
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 de
75 //
76 for (G4int i = 0; i < f_nLbin; ++i) {
77 fSumELongit[i] = fSumE2Longit[i] = fSumELo
78 }
79 for (G4int j = 0; j < f_nRbin; ++j) {
80 fSumERadial[j] = fSumE2Radial[j] = fSumERa
81 }
82 // initialize track length
83 fSumChargTrLength = fSum2ChargTrLength = fSu
84 }
85
86 //....oooOO0OOooo........oooOO0OOooo........oo
87
88 void Run::InitializePerEvent()
89 {
90 // initialize arrays of energy deposit per b
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........oo
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 * fChar
129 fSumNeutrTrLength += fNeutrTrLength;
130 fSum2NeutrTrLength += fNeutrTrLength * fNeut
131
132 // fill histograms
133 //
134
135 G4double Ekin = fKin->GetParticleGun()->GetP
136 G4double mass = fKin->GetParticleGun()->GetP
137 G4double radl = fDet->GetMaterial()->GetRadl
138
139 G4AnalysisManager* analysisManager = G4Analy
140 analysisManager->FillH1(1, 100. * dLCumul /
141 analysisManager->FillH1(2, fChargTrLength /
142 analysisManager->FillH1(3, fNeutrTrLength /
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_d
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_d
155 }
156 }
157
158 //....oooOO0OOooo........oooOO0OOooo........oo
159
160 void Run::Merge(const G4Run* run)
161 {
162 const Run* localRun = static_cast<const 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[
170 fSumELongitCumul[i] += localRun->fSumELong
171 fSumE2LongitCumul[i] += localRun->fSumE2Lo
172 }
173
174 for (G4int j = 0; j < f_nRbin; ++j) {
175 fSumERadial[j] += localRun->fSumERadial[j]
176 fSumE2Radial[j] += localRun->fSumE2Radial[
177 fSumERadialCumul[j] += localRun->fSumERadi
178 fSumE2RadialCumul[j] += localRun->fSumE2Ra
179 }
180
181 fSumChargTrLength += localRun->fSumChargTrLe
182 fSum2ChargTrLength += localRun->fSum2ChargTr
183 fSumNeutrTrLength += localRun->fSumNeutrTrLe
184 fSum2NeutrTrLength += localRun->fSum2NeutrTr
185
186 G4Run::Merge(run);
187 }
188
189 //....oooOO0OOooo........oooOO0OOooo........oo
190
191 void Run::EndOfRun(G4double edep, G4double rms
192 {
193 G4int NbOfEvents = GetNumberOfEvent();
194
195 G4double kinEnergy = fKin->GetParticleGun()-
196 assert(NbOfEvents * kinEnergy > 0);
197
198 fChargedStep /= G4double(NbOfEvents);
199 fNeutralStep /= G4double(NbOfEvents);
200
201 G4double mass = fKin->GetParticleGun()->GetP
202 G4double norme = 100. / (NbOfEvents * (kinEn
203
204 // longitudinal
205 //
206 G4double dLradl = fDet->GetdLradl();
207
208 MyVector MeanELongit(f_nLbin), rmsELongit(f_
209 MyVector MeanELongitCumul(f_nLbin), rmsELong
210
211 G4AnalysisManager* analysisManager = G4Analy
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 *
218
219 MeanELongitCumul[i] = norme * fSumELongitC
220 rmsELongitCumul[i] = norme
221 * std::sqrt(std::abs(
222
223 G4double bin = (i + 0.5) * dLradl;
224 analysisManager->FillH1(4, bin, MeanELongi
225 analysisManager->FillH1(5, bin, rmsELongit
226 bin = (i + 1) * dLradl;
227 analysisManager->FillH1(6, bin, MeanELongi
228 analysisManager->FillH1(7, bin, rmsELongit
229 }
230
231 // radial
232 //
233 G4double dRradl = fDet->GetdRradl();
234
235 MyVector MeanERadial(f_nRbin), rmsERadial(f_
236 MyVector MeanERadialCumul(f_nRbin), rmsERadi
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 *
242
243 MeanERadialCumul[i] = norme * fSumERadialC
244 rmsERadialCumul[i] = norme
245 * std::sqrt(std::abs(
246
247
248 G4double bin = (i + 0.5) * dRradl;
249 analysisManager->FillH1(8, bin, MeanERadia
250 analysisManager->FillH1(9, bin, rmsERadial
251 bin = (i + 1) * dRradl;
252 analysisManager->FillH1(10, bin, MeanERadi
253 analysisManager->FillH1(11, bin, rmsERadia
254 }
255
256 // find Moliere confinement
257 //
258 const G4double EMoliere = 90.;
259 G4double iMoliere = 0.;
260 if ((MeanERadialCumul[0] <= EMoliere) && (Me
261 G4int imin = 0;
262 while ((imin < f_nRbin - 1) && (MeanERadia
263 ++imin;
264 }
265
266 G4double del = MeanERadialCumul[imin + 1]
267 G4double ratio = (del > 0.0) ? (EMoliere -
268 iMoliere = 1. + imin + ratio;
269 }
270
271 // track length
272 //
273 norme = 1. / (NbOfEvents * (fDet->GetMateria
274 G4double MeanChargTrLength = norme * fSumCha
275 G4double rmsChargTrLength =
276 norme
277 * std::sqrt(std::abs(NbOfEvents * fSum2Cha
278
279 G4double MeanNeutrTrLength = norme * fSumNeu
280 G4double rmsNeutrTrLength =
281 norme
282 * std::sqrt(std::abs(NbOfEvents * fSum2Neu
283
284 // print
285 std::ios::fmtflags mode = G4cout.flags();
286 G4cout.setf(std::ios::fixed, std::ios::float
287 G4int prec = G4cout.precision(2);
288
289 if (fVerbose) {
290 G4cout << " LOGITUDINAL PR
291 << " CUMULATIVE LOGITUDINAL PR
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 + d
301
302 G4cout << std::setw(8) << inf << "->" <<
303 << MeanELongit[i] << "% " << std
304 << " 0->" << std::setw(5) <<
305 << MeanELongitCumul[i] << "% " <
306 << G4endl;
307 }
308
309 G4cout << G4endl << G4endl << G4endl;
310
311 G4cout << " RADIAL PROFIL
312 << " CUMULATIVE RADIAL PROFIL
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 + d
322
323 G4cout << std::setw(8) << inf << "->" <<
324 << MeanERadial[i] << "% " << std
325 << " 0->" << std::setw(5) <<
326 << MeanERadialCumul[i] << "% " <
327 << G4endl;
328 }
329 }
330
331 G4cout << "\n ===== SUMMARY ===== \n" << G4e
332
333 G4cout << " Total number of events: "
334 << " Mean number of charged steps: "
335 G4cout << " Mean number of neutral steps: "
336
337 G4cout << " energy deposit : " << std::setw(
338 << std::setw(7) << rmsELongitCumul[f_
339 G4cout << " charged traklen: " << std::setw(
340 << rmsChargTrLength << " radl" << G4e
341 G4cout << " neutral traklen: " << std::setw(
342 << rmsNeutrTrLength << " radl" << G4e
343
344 if (iMoliere > 0.) {
345 G4double RMoliere1 = iMoliere * fDet->Getd
346 G4double RMoliere2 = iMoliere * fDet->Getd
347 G4cout << "\n " << EMoliere << " % confine
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 Abs
362 G4double e = MeanELongitCumul[nLbin - 1] /
363 G4double r = rmsELongitCumul[nLbin - 1] /
364 acc.EmAcceptanceGauss("Edep", NbOfEvents,
365 acc.EmAcceptanceGauss("Erms", NbOfEvents,
366 acc.EndOfAcceptance();
367 }
368 limit = DBL_MAX;
369 }
370
371 //....oooOO0OOooo........oooOO0OOooo........oo
372