| 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/fanoCavity2/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 : G4Run(),
53 fDetector(det),
54 fKinematic(kin),
55 fProcCounter(0),
56 fEdepCavity(0.),
57 fEdepCavity2(0.),
58 fTrkSegmCavity(0.),
59 fNbEventCavity(0),
60 fStepWall(0.),
61 fStepWall2(0.),
62 fStepCavity(0.),
63 fStepCavity2(0.),
64 fNbStepWall(0),
65 fNbStepCavity(0),
66 fEnergyGun(0.),
67 fMassWall(0.),
68 fMassCavity(0.),
69 fIsMaster(isMaster)
70 {
71 // run conditions
72 //
73 G4ParticleDefinition* particleGun = fKinemat
74 G4String partName = particleGun->GetParticle
75 fEnergyGun = fKinematic->GetParticleGun()->G
76
77 // geometry : effective wall volume
78 //
79 G4double cavityThickness = fDetector->GetCav
80 G4Material* mateCavity = fDetector->GetCavit
81 G4double densityCavity = mateCavity->GetDens
82 fMassCavity = cavityThickness * densityCavit
83
84 G4double wallThickness = fDetector->GetWallT
85 G4Material* mateWall = fDetector->GetWallMat
86 G4double densityWall = mateWall->GetDensity(
87
88 G4EmCalculator emCal;
89 G4double RangeWall = emCal.GetCSDARange(fEne
90 G4double factor = 1.2;
91 G4double effWallThick = factor * RangeWall;
92 if ((effWallThick > wallThickness) || (effWa
93 fMassWall = 2 * effWallThick * densityWall;
94
95 G4double massTotal = fMassWall + fMassCavity
96 G4double fMassWallRatio = fMassWall / massTo
97 fKinematic->RunInitialisation(effWallThick,
98
99 G4double massRatio = fMassCavity / fMassWall
100
101 // check radius
102 //
103 G4double worldRadius = fDetector->GetWallRad
104 G4double RangeCavity = emCal.GetCSDARange(fE
105
106 // G4String partName = particleGun->GetPa
107
108 std::ios::fmtflags mode = G4cout.flags();
109 G4cout.setf(std::ios::fixed, std::ios::float
110 G4int prec = G4cout.precision(3);
111
112 G4cout << "\n ===================== run cond
113
114 G4cout << "\n The run will be " << numberOfE
115 << G4BestUnit(fEnergyGun, "Energy") <
116 << " of " << mateWall->GetName()
117 << " (density: " << G4BestUnit(densit
118 << "); Mass/cm2 = " << G4BestUnit(fMa
119 << "\n csdaRange: " << G4BestUnit(Ran
120
121 G4cout << "\n the cavity is " << G4BestUnit(
122 << mateCavity->GetName() << " (densit
123 << "); Mass/cm2 = " << G4BestUnit(fMa
124 << " --> massRatio = " << std::setpre
125
126 G4cout.precision(3);
127 G4cout << " Wall radius: " << G4BestUnit(wor
128 << "; range in cavity: " << G4BestUni
129
130 G4cout << "\n ==============================
131
132 // stopping power from EmCalculator
133 //
134 G4double dedxWall = emCal.GetDEDX(fEnergyGun
135 dedxWall /= densityWall;
136 G4double dedxCavity = emCal.GetDEDX(fEnergyG
137 dedxCavity /= densityCavity;
138
139 G4cout << std::setprecision(4)
140 << "\n StoppingPower in wall = " <<
141 << "\n in cavity = " <<
142 << G4endl;
143
144 // process counter
145 //
146 fProcCounter = new ProcessesCount;
147
148 // charged particles and energy flow in cavi
149 //
150 fPartFlowCavity[0] = fPartFlowCavity[1] = 0;
151 fEnerFlowCavity[0] = fEnerFlowCavity[1] = 0.
152
153 // total energy deposit and charged track se
154 //
155 fEdepCavity = fEdepCavity2 = fTrkSegmCavity
156 fNbEventCavity = 0;
157
158 // stepLenth of charged particles
159 //
160 fStepWall = fStepWall2 = fStepCavity = fStep
161 fNbStepWall = fNbStepCavity = 0;
162
163 // reset default formats
164 G4cout.setf(mode, std::ios::floatfield);
165 G4cout.precision(prec);
166
167 // histograms
168 //
169 G4AnalysisManager* analysisManager = G4Analy
170 if (analysisManager->IsActive()) {
171 analysisManager->OpenFile();
172 }
173 }
174
175 //....oooOO0OOooo........oooOO0OOooo........oo
176
177 Run::~Run() {}
178
179 //....oooOO0OOooo........oooOO0OOooo........oo
180
181 void Run::CountProcesses(G4String procName)
182 {
183 // does the process already encounted ?
184 size_t nbProc = fProcCounter->size();
185 size_t i = 0;
186 while ((i < nbProc) && ((*fProcCounter)[i]->
187 i++;
188 if (i == nbProc) fProcCounter->push_back(new
189
190 (*fProcCounter)[i]->Count();
191 }
192
193 //....oooOO0OOooo........oooOO0OOooo........oo
194
195 void Run::SurveyConvergence(G4int NbofEvents)
196 {
197 if (NbofEvents == 0) return;
198
199 // beam fluence
200 //
201 G4int Nwall = fKinematic->GetWallCount();
202 G4int Ncavity = fKinematic->GetCavityCount()
203 G4double Iwall = Nwall / fMassWall;
204 G4double Icavity = Ncavity / fMassCavity;
205 G4double Iratio = Icavity / Iwall;
206 G4double Itot = NbofEvents / (fMassWall + fM
207 G4double energyFluence = fEnergyGun * Itot;
208
209 // total dose in cavity
210 //
211 G4double doseCavity = fEdepCavity / fMassCav
212 G4double ratio = doseCavity / energyFluence;
213 G4double err = 100 * (ratio - 1.);
214
215 std::ios::fmtflags mode = G4cout.flags();
216 G4cout.setf(std::ios::fixed, std::ios::float
217 G4int prec = G4cout.precision(5);
218
219 G4cout << "--->evntNb= " << NbofEvents << "
220 << " Ic/Iw= " << Iratio << " Ne-_cav=
221 << " doseCavity/Ebeam= " << ratio <<
222 << G4endl;
223
224 // reset default formats
225 G4cout.setf(mode, std::ios::floatfield);
226 G4cout.precision(prec);
227 }
228
229 //....oooOO0OOooo........oooOO0OOooo........oo
230
231 void Run::EndOfRun()
232 { // Only call by Master thread
233
234 std::ios::fmtflags mode = G4cout.flags();
235 G4cout.setf(std::ios::fixed, std::ios::float
236 G4int prec = G4cout.precision(3);
237
238 if (numberOfEvent == 0) return;
239
240 // frequency of processes
241 //
242 G4cout << "\n Process calls frequency --->";
243 for (size_t i = 0; i < fProcCounter->size();
244 G4String procName = (*fProcCounter)[i]->Ge
245 G4int count = (*fProcCounter)[i]->GetCount
246 G4cout << " " << procName << "= " << coun
247 }
248 G4cout << G4endl;
249
250 // charged particle flow in cavity
251 //
252 G4cout << "\n Charged particle flow in cavit
253 << "\n Enter --> nbParticles = "
254 << "\t Energy = " << G4BestUnit(fEner
255 << "\n Exit --> nbParticles = "
256 << "\t Energy = " << G4BestUnit(fEner
257
258 if (fPartFlowCavity[0] == 0) return;
259
260 G4int Nwall = fKinematic->GetWallCount();
261 G4int Ncavity = fKinematic->GetCavityCount()
262
263 G4double Iwall = Nwall / fMassWall;
264 G4double Icavity = Ncavity / fMassCavity;
265 G4double Iratio = Icavity / Iwall;
266 G4double Itot = numberOfEvent / (fMassWall +
267 G4double energyFluence = fEnergyGun * Itot;
268
269 G4cout.precision(5);
270 G4cout << "\n beamFluence in wall = " << Nwa
271 << "\t Icav/Iwall = " << Iratio
272 << "\t energyFluence = " << energyFlu
273
274 // error on Edep in cavity
275 //
276 if (fNbEventCavity == 0) return;
277 G4double meanEdep = fEdepCavity / fNbEventCa
278 G4double meanEdep2 = fEdepCavity2 / fNbEvent
279 G4double varianceEdep = meanEdep2 - meanEdep
280 G4double dEoverE = 0.;
281 if (varianceEdep > 0.) dEoverE = std::sqrt(v
282
283 // total dose in cavity
284 //
285 G4double doseCavity = fEdepCavity / fMassCav
286
287 G4double ratio = doseCavity / energyFluence,
288
289 G4cout << "\n Total edep in cavity = " << G4
290 << 100 * dEoverE << " %"
291 << "\n Total dose in cavity = " << do
292 << " +- " << 100 * dEoverE << " %"
293 << "\n\n DoseCavity/EnergyFluence = "
294
295 // track length in cavity
296 G4double meantrack = fTrkSegmCavity / fPartF
297
298 G4cout.precision(4);
299 G4cout << "\n Total charged trackLength in c
300 << " (mean value = " << G4BestUnit(
301
302 // compute mean step size of charged particl
303 //
304 fStepWall /= fNbStepWall;
305 fStepWall2 /= fNbStepWall;
306 G4double rms = fStepWall2 - fStepWall * fSte
307 if (rms > 0.)
308 rms = std::sqrt(rms);
309 else
310 rms = 0.;
311 G4double nbTrackWall = fKinematic->GetWallCo
312
313 G4cout << "\n StepSize of ch. tracks in wall
314 << G4BestUnit(rms, "Length") << "\t (
315 << ")";
316
317 fStepCavity /= fNbStepCavity;
318 fStepCavity2 /= fNbStepCavity;
319 rms = fStepCavity2 - fStepCavity * fStepCavi
320 if (rms > 0.)
321 rms = std::sqrt(rms);
322 else
323 rms = 0.;
324
325 G4cout << "\n StepSize of ch. tracks in cavi
326 << G4BestUnit(rms, "Length")
327 << "\t (nbSteps/track = " << double(f
328
329 G4cout << G4endl;
330
331 // reset default formats
332 G4cout.setf(mode, std::ios::floatfield);
333 G4cout.precision(prec);
334
335 // delete and remove all contents in fProcCo
336 while (fProcCounter->size() > 0) {
337 OneProcessCount* aProcCount = fProcCounter
338 fProcCounter->pop_back();
339 delete aProcCount;
340 }
341 delete fProcCounter;
342
343 // show Rndm status
344 CLHEP::HepRandom::showEngineStatus();
345 }
346
347 //....oooOO0OOooo........oooOO0OOooo........oo
348
349 void Run::Merge(const G4Run* run)
350 {
351 const Run* localRun = static_cast<const Run*
352
353 // Merge Run variables
354 fPartFlowCavity[0] += localRun->fPartFlowCav
355 fPartFlowCavity[1] += localRun->fPartFlowCav
356 fEnerFlowCavity[0] += localRun->fEnerFlowCav
357 fEnerFlowCavity[1] += localRun->fEnerFlowCav
358 fEdepCavity += localRun->fEdepCavity;
359 fEdepCavity2 += localRun->fEdepCavity2;
360 fTrkSegmCavity += localRun->fTrkSegmCavity;
361 fNbEventCavity += localRun->fNbEventCavity;
362 fStepWall += localRun->fStepWall;
363 fStepWall2 += localRun->fStepWall2;
364 fStepCavity += localRun->fStepCavity;
365 fStepCavity2 += localRun->fStepCavity2;
366 fNbStepWall += localRun->fNbStepWall;
367 fNbStepCavity += localRun->fNbStepCavity;
368
369 // Merge PrimaryGenerator variables
370 fKinematic->AddWallCount(localRun->fKinemati
371 fKinematic->AddCavityCount(localRun->fKinema
372
373 // Merge ProcessCount varaibles
374 std::vector<OneProcessCount*>::iterator it;
375 for (it = localRun->fProcCounter->begin(); i
376 OneProcessCount* process = *it;
377 for (G4int i = 0; i < process->GetCounter(
378 this->CountProcesses(process->GetName())
379 }
380
381 G4Run::Merge(run);
382 }
383
384 //....oooOO0OOooo........oooOO0OOooo........oo
385