| 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 medical/fanoCavity/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 "G4Electron.hh"
40 #include "G4EmCalculator.hh"
41 #include "G4Run.hh"
42 #include "G4RunManager.hh"
43 #include "G4SystemOfUnits.hh"
44 #include "G4UnitsTable.hh"
45 #include "Randomize.hh"
46
47 #include <iomanip>
48
49 //....oooOO0OOooo........oooOO0OOooo........oo
50
51 Run::Run(DetectorConstruction* det, PrimaryGen
52 : fDetector(det), fKinematic(kin), fProcCoun
53
54 {
55 // geometry
56 //
57 fWallThickness = fDetector->GetWallThickness
58 fWallRadius = fDetector->GetWallRadius();
59 fMateWall = fDetector->GetWallMaterial();
60 fDensityWall = fMateWall->GetDensity();
61
62 fCavityThickness = fDetector->GetCavityThick
63 fCavityRadius = fDetector->GetCavityRadius()
64 fSurfaceCavity = CLHEP::pi * fCavityRadius *
65 fVolumeCavity = fSurfaceCavity * fCavityThic
66 fMateCavity = fDetector->GetCavityMaterial()
67 fDensityCavity = fMateCavity->GetDensity();
68 fMassCavity = fVolumeCavity * fDensityCavity
69
70 // process counter
71 //
72 fProcCounter = new ProcessesCount;
73
74 // kinetic energy of charged secondary a cre
75 //
76 fEsecondary = fEsecondary2 = 0.;
77 fNbSec = 0;
78
79 // charged particles and energy flow in cavi
80 //
81 fPartFlowCavity[0] = fPartFlowCavity[1] = 0;
82 fEnerFlowCavity[0] = fEnerFlowCavity[1] = 0.
83
84 // total energy deposit and charged track se
85 //
86 fEdepCavity = fEdepCavity2 = fTrkSegmCavity
87 fNbEventCavity = 0;
88
89 // survey convergence
90 //
91 fOldEmean = fOldDose = 0.;
92
93 // stepLenth of charged particles
94 //
95 fStepWall = fStepWall2 = fStepCavity = fStep
96 fNbStepWall = fNbStepCavity = 0;
97
98 // histograms
99 //
100 G4AnalysisManager* analysisManager = G4Analy
101 if (analysisManager->IsActive()) {
102 analysisManager->OpenFile();
103 }
104 }
105
106 //....oooOO0OOooo........oooOO0OOooo........oo
107
108 Run::~Run() {}
109
110 //....oooOO0OOooo........oooOO0OOooo........oo
111
112 void Run::EndOfRun()
113 { // Only call by Master thread
114 std::ios::fmtflags mode = G4cout.flags();
115 G4cout.setf(std::ios::fixed, std::ios::float
116
117 if (numberOfEvent == 0) return;
118
119 // run conditions
120 //
121 G4ParticleDefinition* particle = fKinematic-
122 G4String partName = particle->GetParticleNam
123 G4double energy = fKinematic->GetParticleGun
124
125 G4cout << "\n ======================== run s
126
127 G4int prec = G4cout.precision(3);
128
129 G4cout << "\n The run consists of " << numbe
130 << G4BestUnit(energy, "Energy") << "
131 << " of " << fMateWall->GetName()
132 << " (density: " << G4BestUnit(fDensi
133
134 G4cout << "\n the cavity is " << G4BestUnit(
135 << fMateCavity->GetName() << " (densi
136 << "); Mass = " << G4BestUnit(fMassCa
137
138 G4cout << "\n ==============================
139
140 // frequency of processes
141 //
142 G4cout << "\n Process calls frequency --->";
143 for (size_t i = 0; i < fProcCounter->size();
144 G4String procName = (*fProcCounter)[i]->Ge
145 G4int count = (*fProcCounter)[i]->GetCount
146 G4cout << " " << procName << "= " << coun
147 }
148 G4cout << G4endl;
149
150 // extract cross sections with G4EmCalculato
151 //
152 G4EmCalculator emCalculator;
153 G4cout << "\n Gamma crossSections in wall ma
154 G4double sumc = 0.0;
155 for (size_t i = 0; i < fProcCounter->size();
156 G4String procName = (*fProcCounter)[i]->Ge
157 G4double massSigma =
158 emCalculator.ComputeCrossSectionPerVolum
159 / fDensityWall;
160 if (massSigma > 0.) {
161 sumc += massSigma;
162 G4cout << " " << procName << "= " << G4
163 }
164 }
165 G4cout << " --> total= " << G4BestUnit(sum
166
167 // mean kinetic energy of secondary electron
168 //
169 if (fNbSec == 0) return;
170 G4double meanEsecond = fEsecondary / fNbSec,
171 G4double varianceEsec = meanEsecond2 - meanE
172 G4double dToverT = 0.;
173 if (varianceEsec > 0.) dToverT = std::sqrt(v
174 G4double csdaRange = emCalculator.GetCSDARan
175
176 G4cout.precision(4);
177 G4cout << "\n Mean energy of secondary e- =
178 << 100 * dToverT << " %"
179 << " (--> range in wall material = "
180
181 // compute mass energy transfer coefficient
182 //
183 G4double massTransfCoef = sumc * meanEsecond
184
185 G4cout << " Mass_energy_transfer coef: " <<
186
187 // stopping power from EmCalculator
188 //
189 G4double dedxWall = emCalculator.GetDEDX(mea
190 dedxWall /= fDensityWall;
191 G4double dedxCavity = emCalculator.GetDEDX(m
192 dedxCavity /= fDensityCavity;
193
194 G4cout << "\n StoppingPower in wall = " <<
195 << "\n in cavity = " <<
196 << G4endl;
197
198 // charged particle flow in cavity
199 //
200 G4cout << "\n Charged particle flow in cavit
201 << "\n Enter --> nbParticles = "
202 << "\t Energy = " << G4BestUnit(fEner
203 << "\n Exit --> nbParticles = "
204 << "\t Energy = " << G4BestUnit(fEner
205
206 if (fPartFlowCavity[0] == 0) return;
207
208 // beam energy fluence
209 //
210 G4double rBeam = fWallRadius * (fKinematic->
211 G4double surfaceBeam = CLHEP::pi * rBeam * r
212
213 // error on Edep in cavity
214 //
215 if (fNbEventCavity == 0) return;
216 G4double meanEdep = fEdepCavity / fNbEventCa
217 G4double meanEdep2 = fEdepCavity2 / fNbEvent
218 G4double varianceEdep = meanEdep2 - meanEdep
219 G4double dEoverE = 0.;
220 if (varianceEdep > 0.) dEoverE = std::sqrt(v
221
222 // total dose in cavity
223 //
224 G4double doseCavity = fEdepCavity / fMassCav
225 G4double doseOverBeam = doseCavity * surface
226
227 // track length in cavity
228 G4double meantrack = fTrkSegmCavity / fPartF
229
230 G4cout.precision(4);
231 G4cout << "\n Total edep in cavity = " << G4
232 << 100 * dEoverE << " %"
233 << "\t Total charged trackLength = "
234 << " (mean value = " << G4BestUnit(
235 << "\n Total dose in cavity = " << do
236 << "\n Dose/EnergyFluence = " << G4
237
238 // ratio simulation/theory
239 //
240 G4double ratio = doseOverBeam / massTransfCo
241 G4double error = ratio * std::sqrt(dEoverE *
242
243 G4cout.precision(5);
244 G4cout << "\n (Dose/EnergyFluence)/Mass_ener
245
246 // compute mean step size of charged particl
247 //
248 fStepWall /= fNbStepWall;
249 fStepWall2 /= fNbStepWall;
250 G4double rms = fStepWall2 - fStepWall * fSte
251 if (rms > 0.)
252 rms = std::sqrt(rms);
253 else
254 rms = 0.;
255
256 G4cout.precision(4);
257 G4cout << "\n StepSize of ch. tracks in wall
258 << G4BestUnit(rms, "Length") << "\t (
259 << ")";
260
261 fStepCavity /= fNbStepCavity;
262 fStepCavity2 /= fNbStepCavity;
263 rms = fStepCavity2 - fStepCavity * fStepCavi
264 if (rms > 0.)
265 rms = std::sqrt(rms);
266 else
267 rms = 0.;
268
269 G4cout << "\n StepSize of ch. tracks in cavi
270 << G4BestUnit(rms, "Length")
271 << "\t (nbSteps/track = " << double(f
272
273 G4cout << G4endl;
274
275 // reset default formats
276 G4cout.setf(mode, std::ios::floatfield);
277 G4cout.precision(prec);
278
279 // delete and remove all contents in fProcCo
280 while (fProcCounter->size() > 0) {
281 OneProcessCount* aProcCount = fProcCounter
282 fProcCounter->pop_back();
283 delete aProcCount;
284 }
285 delete fProcCounter;
286 }
287
288 //....oooOO0OOooo........oooOO0OOooo........oo
289
290 void Run::SurveyConvergence(G4int NbofEvents)
291 {
292 if (NbofEvents == 0) return;
293
294 // mean kinetic energy of secondary electron
295 //
296 G4double meanEsecond = 0.;
297 if (fNbSec > 0) meanEsecond = fEsecondary /
298 G4double rateEmean = 0.;
299 // compute variation rate (%), iteration to
300 if (fOldEmean > 0.) rateEmean = 100 * (meanE
301 fOldEmean = meanEsecond;
302
303 // beam energy fluence
304 //
305 G4double rBeam = fWallRadius * (fKinematic->
306 G4double surfaceBeam = CLHEP::pi * rBeam * r
307 G4double beamEnergy = fKinematic->GetParticl
308
309 // total dose in cavity
310 //
311 G4double doseCavity = fEdepCavity / fMassCav
312 G4double doseOverBeam = doseCavity * surface
313 G4double rateDose = 0.;
314 // compute variation rate (%), iteration to
315 if (fOldDose > 0.) rateDose = 100 * (doseOve
316 fOldDose = doseOverBeam;
317
318 std::ios::fmtflags mode = G4cout.flags();
319 G4cout.setf(std::ios::fixed, std::ios::float
320 G4int prec = G4cout.precision(3);
321
322 G4cout << " ---> NbofEvents= " << NbofEvents
323 << " Tkin= " << G4BestUnit(meanEsec
324 << " NbOfelec in cav= " << fPartFlo
325 << " Dose/EnFluence= " << G4BestUni
326 << " %) \n"
327 << G4endl;
328
329 // reset default formats
330 G4cout.setf(mode, std::ios::floatfield);
331 G4cout.precision(prec);
332 }
333
334 //....oooOO0OOooo........oooOO0OOooo........oo
335
336 void Run::Merge(const G4Run* run)
337 {
338 const Run* localRun = static_cast<const Run*
339
340 fPartFlowCavity[0] += localRun->fPartFlowCav
341 fPartFlowCavity[1] += localRun->fPartFlowCav
342 fEnerFlowCavity[0] += localRun->fEnerFlowCav
343 fEnerFlowCavity[1] += localRun->fEnerFlowCav
344 fEdepCavity += localRun->fEdepCavity;
345 fEdepCavity2 += localRun->fEdepCavity2;
346 fTrkSegmCavity += localRun->fTrkSegmCavity;
347 fNbEventCavity += localRun->fNbEventCavity;
348
349 fStepWall += localRun->fStepWall;
350 fStepWall2 += localRun->fStepWall2;
351 fStepCavity += localRun->fStepCavity;
352 fStepCavity2 += localRun->fStepCavity2;
353 fNbStepWall += localRun->fNbStepWall;
354 fNbStepCavity += localRun->fNbStepCavity;
355
356 fEsecondary += localRun->fEsecondary;
357 fEsecondary2 += localRun->fEsecondary2;
358
359 fNbSec += localRun->fNbSec;
360
361 // ??? G4double fOldEmean
362 // ??? G4Double fOldDose;
363
364 // Merge ProcessCount varaibles
365 std::vector<OneProcessCount*>::iterator it;
366 for (it = localRun->fProcCounter->begin(); i
367 OneProcessCount* process = *it;
368 for (G4int i = 0; i < process->GetCounter(
369 this->CountProcesses(process->GetName())
370 }
371
372 G4Run::Merge(run);
373 }
374
375 //....oooOO0OOooo........oooOO0OOooo........oo
376
377 void Run::CountProcesses(G4String procName)
378 {
379 // does the process already encounted ?
380 size_t nbProc = fProcCounter->size();
381 size_t i = 0;
382 while ((i < nbProc) && ((*fProcCounter)[i]->
383 i++;
384 if (i == nbProc) fProcCounter->push_back(new
385
386 (*fProcCounter)[i]->Count();
387 }
388
389 //....oooOO0OOooo........oooOO0OOooo........oo
390