| 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 medical/fanoCavity/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 "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........oooOO0OOooo........oooOO0OOooo......
50
51 Run::Run(DetectorConstruction* det, PrimaryGeneratorAction* kin)
52 : fDetector(det), fKinematic(kin), fProcCounter(0), fMateWall(0), fMateCavity(0)
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->GetCavityThickness();
63 fCavityRadius = fDetector->GetCavityRadius();
64 fSurfaceCavity = CLHEP::pi * fCavityRadius * fCavityRadius;
65 fVolumeCavity = fSurfaceCavity * fCavityThickness;
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 creation
75 //
76 fEsecondary = fEsecondary2 = 0.;
77 fNbSec = 0;
78
79 // charged particles and energy flow in cavity
80 //
81 fPartFlowCavity[0] = fPartFlowCavity[1] = 0;
82 fEnerFlowCavity[0] = fEnerFlowCavity[1] = 0.;
83
84 // total energy deposit and charged track segment in cavity
85 //
86 fEdepCavity = fEdepCavity2 = fTrkSegmCavity = 0.;
87 fNbEventCavity = 0;
88
89 // survey convergence
90 //
91 fOldEmean = fOldDose = 0.;
92
93 // stepLenth of charged particles
94 //
95 fStepWall = fStepWall2 = fStepCavity = fStepCavity2 = 0.;
96 fNbStepWall = fNbStepCavity = 0;
97
98 // histograms
99 //
100 G4AnalysisManager* analysisManager = G4AnalysisManager::Instance();
101 if (analysisManager->IsActive()) {
102 analysisManager->OpenFile();
103 }
104 }
105
106 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
107
108 Run::~Run() {}
109
110 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
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::floatfield);
116
117 if (numberOfEvent == 0) return;
118
119 // run conditions
120 //
121 G4ParticleDefinition* particle = fKinematic->GetParticleGun()->GetParticleDefinition();
122 G4String partName = particle->GetParticleName();
123 G4double energy = fKinematic->GetParticleGun()->GetParticleEnergy();
124
125 G4cout << "\n ======================== run summary ======================\n";
126
127 G4int prec = G4cout.precision(3);
128
129 G4cout << "\n The run consists of " << numberOfEvent << " " << partName << " of "
130 << G4BestUnit(energy, "Energy") << " through 2*" << G4BestUnit(fWallThickness, "Length")
131 << " of " << fMateWall->GetName()
132 << " (density: " << G4BestUnit(fDensityWall, "Volumic Mass") << ")" << G4endl;
133
134 G4cout << "\n the cavity is " << G4BestUnit(fCavityThickness, "Length") << " of "
135 << fMateCavity->GetName() << " (density: " << G4BestUnit(fDensityCavity, "Volumic Mass")
136 << "); Mass = " << G4BestUnit(fMassCavity, "Mass") << G4endl;
137
138 G4cout << "\n ============================================================\n";
139
140 // frequency of processes
141 //
142 G4cout << "\n Process calls frequency --->";
143 for (size_t i = 0; i < fProcCounter->size(); i++) {
144 G4String procName = (*fProcCounter)[i]->GetName();
145 G4int count = (*fProcCounter)[i]->GetCounter();
146 G4cout << " " << procName << "= " << count;
147 }
148 G4cout << G4endl;
149
150 // extract cross sections with G4EmCalculator
151 //
152 G4EmCalculator emCalculator;
153 G4cout << "\n Gamma crossSections in wall material :";
154 G4double sumc = 0.0;
155 for (size_t i = 0; i < fProcCounter->size(); i++) {
156 G4String procName = (*fProcCounter)[i]->GetName();
157 G4double massSigma =
158 emCalculator.ComputeCrossSectionPerVolume(energy, particle, procName, fMateWall)
159 / fDensityWall;
160 if (massSigma > 0.) {
161 sumc += massSigma;
162 G4cout << " " << procName << "= " << G4BestUnit(massSigma, "Surface/Mass");
163 }
164 }
165 G4cout << " --> total= " << G4BestUnit(sumc, "Surface/Mass") << G4endl;
166
167 // mean kinetic energy of secondary electrons
168 //
169 if (fNbSec == 0) return;
170 G4double meanEsecond = fEsecondary / fNbSec, meanEsecond2 = fEsecondary2 / fNbSec;
171 G4double varianceEsec = meanEsecond2 - meanEsecond * meanEsecond;
172 G4double dToverT = 0.;
173 if (varianceEsec > 0.) dToverT = std::sqrt(varianceEsec / fNbSec) / meanEsecond;
174 G4double csdaRange = emCalculator.GetCSDARange(meanEsecond, G4Electron::Electron(), fMateWall);
175
176 G4cout.precision(4);
177 G4cout << "\n Mean energy of secondary e- = " << G4BestUnit(meanEsecond, "Energy") << " +- "
178 << 100 * dToverT << " %"
179 << " (--> range in wall material = " << G4BestUnit(csdaRange, "Length") << ")" << G4endl;
180
181 // compute mass energy transfer coefficient
182 //
183 G4double massTransfCoef = sumc * meanEsecond / energy;
184
185 G4cout << " Mass_energy_transfer coef: " << G4BestUnit(massTransfCoef, "Surface/Mass") << G4endl;
186
187 // stopping power from EmCalculator
188 //
189 G4double dedxWall = emCalculator.GetDEDX(meanEsecond, G4Electron::Electron(), fMateWall);
190 dedxWall /= fDensityWall;
191 G4double dedxCavity = emCalculator.GetDEDX(meanEsecond, G4Electron::Electron(), fMateCavity);
192 dedxCavity /= fDensityCavity;
193
194 G4cout << "\n StoppingPower in wall = " << G4BestUnit(dedxWall, "Energy*Surface/Mass")
195 << "\n in cavity = " << G4BestUnit(dedxCavity, "Energy*Surface/Mass")
196 << G4endl;
197
198 // charged particle flow in cavity
199 //
200 G4cout << "\n Charged particle flow in cavity :"
201 << "\n Enter --> nbParticles = " << fPartFlowCavity[0]
202 << "\t Energy = " << G4BestUnit(fEnerFlowCavity[0], "Energy")
203 << "\n Exit --> nbParticles = " << fPartFlowCavity[1]
204 << "\t Energy = " << G4BestUnit(fEnerFlowCavity[1], "Energy") << G4endl;
205
206 if (fPartFlowCavity[0] == 0) return;
207
208 // beam energy fluence
209 //
210 G4double rBeam = fWallRadius * (fKinematic->GetBeamRadius());
211 G4double surfaceBeam = CLHEP::pi * rBeam * rBeam;
212
213 // error on Edep in cavity
214 //
215 if (fNbEventCavity == 0) return;
216 G4double meanEdep = fEdepCavity / fNbEventCavity;
217 G4double meanEdep2 = fEdepCavity2 / fNbEventCavity;
218 G4double varianceEdep = meanEdep2 - meanEdep * meanEdep;
219 G4double dEoverE = 0.;
220 if (varianceEdep > 0.) dEoverE = std::sqrt(varianceEdep / fNbEventCavity) / meanEdep;
221
222 // total dose in cavity
223 //
224 G4double doseCavity = fEdepCavity / fMassCavity;
225 G4double doseOverBeam = doseCavity * surfaceBeam / (numberOfEvent * energy);
226
227 // track length in cavity
228 G4double meantrack = fTrkSegmCavity / fPartFlowCavity[0];
229
230 G4cout.precision(4);
231 G4cout << "\n Total edep in cavity = " << G4BestUnit(fEdepCavity, "Energy") << " +- "
232 << 100 * dEoverE << " %"
233 << "\t Total charged trackLength = " << G4BestUnit(fTrkSegmCavity, "Length")
234 << " (mean value = " << G4BestUnit(meantrack, "Length") << ")"
235 << "\n Total dose in cavity = " << doseCavity / (MeV / mg) << " MeV/mg"
236 << "\n Dose/EnergyFluence = " << G4BestUnit(doseOverBeam, "Surface/Mass") << G4endl;
237
238 // ratio simulation/theory
239 //
240 G4double ratio = doseOverBeam / massTransfCoef;
241 G4double error = ratio * std::sqrt(dEoverE * dEoverE + dToverT * dToverT);
242
243 G4cout.precision(5);
244 G4cout << "\n (Dose/EnergyFluence)/Mass_energy_transfer = " << ratio << " +- " << error << G4endl;
245
246 // compute mean step size of charged particles
247 //
248 fStepWall /= fNbStepWall;
249 fStepWall2 /= fNbStepWall;
250 G4double rms = fStepWall2 - fStepWall * fStepWall;
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 = " << G4BestUnit(fStepWall, "Length") << " +- "
258 << G4BestUnit(rms, "Length") << "\t (nbSteps/track = " << double(fNbStepWall) / fNbSec
259 << ")";
260
261 fStepCavity /= fNbStepCavity;
262 fStepCavity2 /= fNbStepCavity;
263 rms = fStepCavity2 - fStepCavity * fStepCavity;
264 if (rms > 0.)
265 rms = std::sqrt(rms);
266 else
267 rms = 0.;
268
269 G4cout << "\n StepSize of ch. tracks in cavity = " << G4BestUnit(fStepCavity, "Length") << " +- "
270 << G4BestUnit(rms, "Length")
271 << "\t (nbSteps/track = " << double(fNbStepCavity) / fPartFlowCavity[0] << ")";
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 fProcCounter
280 while (fProcCounter->size() > 0) {
281 OneProcessCount* aProcCount = fProcCounter->back();
282 fProcCounter->pop_back();
283 delete aProcCount;
284 }
285 delete fProcCounter;
286 }
287
288 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
289
290 void Run::SurveyConvergence(G4int NbofEvents)
291 {
292 if (NbofEvents == 0) return;
293
294 // mean kinetic energy of secondary electrons
295 //
296 G4double meanEsecond = 0.;
297 if (fNbSec > 0) meanEsecond = fEsecondary / fNbSec;
298 G4double rateEmean = 0.;
299 // compute variation rate (%), iteration to iteration
300 if (fOldEmean > 0.) rateEmean = 100 * (meanEsecond / fOldEmean - 1.);
301 fOldEmean = meanEsecond;
302
303 // beam energy fluence
304 //
305 G4double rBeam = fWallRadius * (fKinematic->GetBeamRadius());
306 G4double surfaceBeam = CLHEP::pi * rBeam * rBeam;
307 G4double beamEnergy = fKinematic->GetParticleGun()->GetParticleEnergy();
308
309 // total dose in cavity
310 //
311 G4double doseCavity = fEdepCavity / fMassCavity;
312 G4double doseOverBeam = doseCavity * surfaceBeam / (NbofEvents * beamEnergy);
313 G4double rateDose = 0.;
314 // compute variation rate (%), iteration to iteration
315 if (fOldDose > 0.) rateDose = 100 * (doseOverBeam / fOldDose - 1.);
316 fOldDose = doseOverBeam;
317
318 std::ios::fmtflags mode = G4cout.flags();
319 G4cout.setf(std::ios::fixed, std::ios::floatfield);
320 G4int prec = G4cout.precision(3);
321
322 G4cout << " ---> NbofEvents= " << NbofEvents << " NbOfelectr= " << fNbSec
323 << " Tkin= " << G4BestUnit(meanEsecond, "Energy") << " (" << rateEmean << " %)"
324 << " NbOfelec in cav= " << fPartFlowCavity[0]
325 << " Dose/EnFluence= " << G4BestUnit(doseOverBeam, "Surface/Mass") << " (" << rateDose
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........oooOO0OOooo........oooOO0OOooo......
335
336 void Run::Merge(const G4Run* run)
337 {
338 const Run* localRun = static_cast<const Run*>(run);
339
340 fPartFlowCavity[0] += localRun->fPartFlowCavity[0];
341 fPartFlowCavity[1] += localRun->fPartFlowCavity[1];
342 fEnerFlowCavity[0] += localRun->fEnerFlowCavity[0];
343 fEnerFlowCavity[1] += localRun->fEnerFlowCavity[1];
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(); it != localRun->fProcCounter->end(); it++) {
367 OneProcessCount* process = *it;
368 for (G4int i = 0; i < process->GetCounter(); i++)
369 this->CountProcesses(process->GetName());
370 }
371
372 G4Run::Merge(run);
373 }
374
375 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
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]->GetName() != procName))
383 i++;
384 if (i == nbProc) fProcCounter->push_back(new OneProcessCount(procName));
385
386 (*fProcCounter)[i]->Count();
387 }
388
389 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
390