| 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 // This example is provided by the Geant4-DNA
27 // Any report or published results obtained us
28 // shall cite the following Geant4-DNA collabo
29 // Med. Phys. 37 (2010) 4692-4708
30 // and papers
31 // M. Batmunkh et al. J Radiat Res Appl Sci 8
32 // O. Belov et al. Physica Medica 32 (2016) 15
33 // The Geant4-DNA web site is available at htt
34 //
35 // -------------------------------------------
36 // November 2016
37 // -------------------------------------------
38 //
39 /// \file Run.cc
40 /// \brief Implementation of the Run class
41
42 #include "Run.hh"
43
44 #include "DetectorConstruction.hh"
45 #include "PrimaryGeneratorAction.hh"
46
47 #include "G4AnalysisManager.hh"
48 #include "G4EmCalculator.hh"
49 #include "G4H2O.hh"
50 #include "G4Molecule.hh"
51 #include "G4MoleculeCounter.hh"
52 #include "G4MoleculeGun.hh"
53 #include "G4MoleculeTable.hh"
54 #include "G4SystemOfUnits.hh"
55 #include "G4UnitsTable.hh"
56 #include "G4ios.hh"
57
58 #include <G4Scheduler.hh>
59 #include <iomanip>
60
61 //....oooOO0OOooo........oooOO0OOooo........oo
62
63 Run::Run(DetectorConstruction* det) : G4Run(),
64 {
65 fSoma3DEdep = new G4double[fDetector->GetnbS
66 for (G4int i = 0; i < fDetector->GetnbSomaco
67 fSoma3DEdep[i] = 0.;
68 }
69 fDend3DEdep = new G4double[fDetector->GetnbD
70 for (G4int i = 0; i < fDetector->GetnbDendri
71 fDend3DEdep[i] = 0.;
72 }
73 fAxon3DEdep = new G4double[fDetector->GetnbA
74 for (G4int i = 0; i < fDetector->GetnbAxonco
75 fAxon3DEdep[i] = 0.;
76 }
77
78 fEdepAll = fEdepAll_err = fEdepMedium = fEde
79 fEdepSoma = fEdepSoma_err = 0.;
80 fEdepDend = fEdepDend_err = fEdepAxon = fEde
81 }
82
83 //....oooOO0OOooo........oooOO0OOooo........oo
84
85 Run::~Run()
86 {
87 delete[] fSoma3DEdep;
88 delete[] fDend3DEdep;
89 delete[] fAxon3DEdep;
90 }
91
92 //....oooOO0OOooo........oooOO0OOooo........oo
93
94 void Run::SetPrimary(G4ParticleDefinition* par
95 {
96 fParticle = particle;
97 fEkin = energy;
98 }
99
100 //....oooOO0OOooo........oooOO0OOooo........oo
101
102 void Run::AddPrimaryLET(G4double let)
103 {
104 fLET += let;
105 fLET2 += let * let;
106 }
107
108 //....oooOO0OOooo........oooOO0OOooo........oo
109
110 void Run::SetTrackLength(G4double t)
111 {
112 ftrackLength = t;
113 fTrackLen += t;
114 fTrackLen2 += t * t;
115 }
116
117 //....oooOO0OOooo........oooOO0OOooo........oo
118
119 void Run::CountProcesses(const G4VProcess* pro
120 {
121 G4String procName = process->GetProcessName(
122 std::map<G4String, G4int>::iterator it = fPr
123 if (it == fProcCounter.end()) {
124 fProcCounter[procName] = 1;
125 }
126 else {
127 fProcCounter[procName]++;
128 }
129 }
130
131 //....oooOO0OOooo........oooOO0OOooo........oo
132
133 void Run::ParticleCount(G4String name, G4doubl
134 {
135 std::map<G4String, ParticleData>::iterator i
136 if (it == fParticleDataMap1.end()) {
137 fParticleDataMap1[name] = ParticleData(1,
138 }
139 else {
140 ParticleData& data = it->second;
141 data.fCount++;
142 data.fEmean += Ekin;
143 // update min max
144 G4double emin = data.fEmin;
145 if (Ekin < emin) data.fEmin = Ekin;
146 G4double emax = data.fEmax;
147 if (Ekin > emax) data.fEmax = Ekin;
148 }
149 }
150
151 //....oooOO0OOooo........oooOO0OOooo........oo
152
153 void Run::ParticleCountNeuron(G4String name, G
154 {
155 std::map<G4String, ParticleData>::iterator i
156 if (it == fParticleDataMap2.end()) {
157 fParticleDataMap2[name] = ParticleData(1,
158 }
159 else {
160 ParticleData& data = it->second;
161 data.fCount++;
162 data.fEmean += Ekin;
163 // update min max
164 G4double emin = data.fEmin;
165 if (Ekin < emin) data.fEmin = Ekin;
166 G4double emax = data.fEmax;
167 if (Ekin > emax) data.fEmax = Ekin;
168 }
169 }
170
171 //....oooOO0OOooo........oooOO0OOooo........oo
172 /*
173 void Run::MoleculeCount(G4String, G4double)
174 {
175 //fMoleculeNumber = G4MoleculeCounter::Instanc
176 // ->GetNMoleculesAtTime(mole
177 }
178 */
179 //....oooOO0OOooo........oooOO0OOooo........oo
180
181 void Run::MoleculeCountNeuron(G4Molecule* mole
182 {
183 G4String moleculename = molecule->GetName();
184 std::map<G4String, G4int>::iterator it = fMo
185 if (it == fMoleculeNumber.end()) {
186 fMoleculeNumber[moleculename] = 1;
187 }
188 else {
189 fMoleculeNumber[moleculename]++;
190 }
191 }
192
193 //....oooOO0OOooo........oooOO0OOooo........oo
194
195 void Run::AddEflow(G4double eflow)
196 {
197 fEnergyFlow += eflow;
198 fEnergyFlow2 += eflow * eflow;
199 }
200
201 //....oooOO0OOooo........oooOO0OOooo........oo
202
203 void Run::Merge(const G4Run* run)
204 {
205 const Run* localRun = static_cast<const Run*
206
207 // primary particle info
208 //
209 fParticle = localRun->fParticle;
210 fEkin = localRun->fEkin;
211 ftrackLength = localRun->ftrackLength;
212 fTrackLen += localRun->fTrackLen;
213 fTrackLen2 += localRun->fTrackLen2;
214 fLET += localRun->fLET;
215 fLET2 += localRun->fLET2;
216
217 // map: processes count in simulation medium
218 std::map<G4String, G4int>::const_iterator it
219 for (itp = localRun->fProcCounter.begin(); i
220 G4String procName = itp->first;
221 G4int localCount = itp->second;
222 if (fProcCounter.find(procName) == fProcCo
223 fProcCounter[procName] = localCount;
224 }
225 else {
226 fProcCounter[procName] += localCount;
227 }
228 }
229
230 // map: created particles count outside neur
231 std::map<G4String, ParticleData>::const_iter
232 for (itc = localRun->fParticleDataMap1.begin
233 G4String name = itc->first;
234 const ParticleData& localData = itc->secon
235 if (fParticleDataMap1.find(name) == fParti
236 fParticleDataMap1[name] =
237 ParticleData(localData.fCount, localDa
238 }
239 else {
240 ParticleData& data = fParticleDataMap1[n
241 data.fCount += localData.fCount;
242 data.fEmean += localData.fEmean;
243 G4double emin = localData.fEmin;
244 if (emin < data.fEmin) data.fEmin = emin
245 G4double emax = localData.fEmax;
246 if (emax > data.fEmax) data.fEmax = emax
247 }
248 }
249
250 // map: created particles count inside neuro
251 std::map<G4String, ParticleData>::const_iter
252 for (itn = localRun->fParticleDataMap2.begin
253 G4String name = itn->first;
254 const ParticleData& localData = itn->secon
255 if (fParticleDataMap2.find(name) == fParti
256 fParticleDataMap2[name] =
257 ParticleData(localData.fCount, localDa
258 }
259 else {
260 ParticleData& data = fParticleDataMap2[n
261 data.fCount += localData.fCount;
262 data.fEmean += localData.fEmean;
263 G4double emin = localData.fEmin;
264 if (emin < data.fEmin) data.fEmin = emin
265 G4double emax = localData.fEmax;
266 if (emax > data.fEmax) data.fEmax = emax
267 }
268 }
269
270 std::map<G4String, G4int>::const_iterator it
271 for (itm = localRun->fMoleculeNumber.begin()
272 G4String moleculeName = itm->first;
273 G4int localCount = itm->second;
274 if (fMoleculeNumber.find(moleculeName) ==
275 fMoleculeNumber[moleculeName] = localCou
276 }
277 else {
278 fMoleculeNumber[moleculeName] += localCo
279 }
280 }
281
282 // hits compartments in neuron compartments
283 //
284 for (G4int i = 0; i < fDetector->GetnbSomaco
285 fSoma3DEdep[i] += localRun->fSoma3DEdep[i]
286 }
287 for (G4int i = 0; i < fDetector->GetnbDendri
288 fDend3DEdep[i] += localRun->fDend3DEdep[i]
289 }
290 for (G4int i = 0; i < fDetector->GetnbAxonco
291 fAxon3DEdep[i] += localRun->fAxon3DEdep[i]
292 }
293
294 // accumulate sums
295 //
296 fEdepAll += localRun->fEdepAll;
297 fEdepAll_err += localRun->fEdepAll_err;
298 fEdepMedium += localRun->fEdepMedium;
299 fEdepMedium_err += localRun->fEdepMedium_err
300 fEdepSlice += localRun->fEdepSlice;
301 fEdepSlice_err += localRun->fEdepSlice_err;
302 fEdepSoma += localRun->fEdepSoma;
303 fEdepSoma_err += localRun->fEdepSoma_err;
304 fEdepDend += localRun->fEdepDend;
305 fEdepDend_err += localRun->fEdepDend_err;
306 fEdepAxon += localRun->fEdepAxon;
307 fEdepAxon_err += localRun->fEdepAxon_err;
308 fEdepNeuron += localRun->fEdepNeuron;
309 fEdepNeuron_err += localRun->fEdepNeuron_err
310
311 fEnergyFlow += localRun->fEnergyFlow;
312 fEnergyFlow2 += localRun->fEnergyFlow2;
313
314 G4Run::Merge(run);
315 }
316
317 //....oooOO0OOooo........oooOO0OOooo........oo
318
319 void Run::EndOfRun()
320 {
321 G4int prec = 5, wid = prec + 2;
322 G4int dfprec = G4cout.precision(prec);
323
324 // Characteristics of Primary particle
325 G4String Particle = fParticle->GetParticleNa
326 G4double GunArea;
327
328 G4Material* material = fDetector->GetTargetM
329
330 // compute track length of primary track
331 //
332 fTrackLen /= numberOfEvent;
333 fTrackLen2 /= numberOfEvent;
334 G4double rms = fTrackLen2 - fTrackLen * fTra
335 if (rms > 0.)
336 rms = std::sqrt(rms);
337 else
338 rms = 0.;
339
340 G4int TotNbofEvents = numberOfEvent;
341 G4double EdepTotal = fEdepAll;
342 G4double EdepTotal2 = fEdepAll_err;
343 EdepTotal /= TotNbofEvents;
344 EdepTotal2 /= TotNbofEvents;
345 G4double rmst = EdepTotal2 - EdepTotal * Ede
346 if (rmst > 0.)
347 rmst = std::sqrt(rmst);
348 else
349 rmst = 0.;
350
351 // Stopping Power from input Table.
352 G4EmCalculator emCalculator;
353 G4double dEdxFull = 0.;
354 if (fParticle->GetPDGCharge() != 0.) {
355 dEdxFull = emCalculator.ComputeTotalDEDX(f
356 }
357
358 // Stopping Power from simulation.
359 G4double meandEdx = (EdepTotal / fTrackLen)
360 G4double meandEdxerr = (rmst / fTrackLen) /
361 G4int ncols = 0;
362 G4float tmp, En;
363 G4int Ntrav = 0;
364 FILE* fp = fopen("OutputPerEvent.out", "r");
365 while (1) {
366 ncols = fscanf(fp, " %f %f %f %f %f %f %f
367 if (ncols < 0) break;
368 if (En > 0) Ntrav++;
369 }
370 fclose(fp);
371 // The surface area is calculated as spheric
372 GunArea = fDetector->GetTotSurfMedium();
373 // Fluence dose of single paticle track
374 G4double FluenceDoseperBeam = 0.160218 * (dE
375
376 G4cout << "\n ======= The summary of simulat
377 G4cout << "\n Primary particle
378 << "\n Kinetic energy of beam
379 << "\n Particle traversals the neuro
380 << "\n Full LET of beam as formulas
381 << "\n Mean LET of beam as simulatio
382 << " keV/um"
383 << "\n Mean track length of beam
384 << "\n Particle fluence
385 << " particles/cm^2"
386 << "\n Fluence dose (full)
387 << numberOfEvent * FluenceDoseperBeam
388 << "\n Fluence dose ber beam
389 << G4endl;
390
391 if (numberOfEvent == 0) {
392 G4cout.precision(dfprec);
393 return;
394 }
395
396 // frequency of processes in all volume
397 //
398 G4cout << "\n List of generated physical pro
399
400 G4int index = 0;
401 std::map<G4String, G4int>::iterator it;
402 for (it = fProcCounter.begin(); it != fProcC
403 G4String procName = it->first;
404 G4int count = it->second;
405 G4String space = " ";
406 if (++index % 1 == 0) space = "\n";
407 G4cout << " " << std::setw(20) << procName
408 }
409 G4cout << G4endl;
410
411 // particles count outside neuron structure
412 //
413 G4cout << "\n List of generated particles ou
414
415 std::map<G4String, ParticleData>::iterator i
416 for (itc = fParticleDataMap1.begin(); itc !=
417 G4String name = itc->first;
418 ParticleData data = itc->second;
419 G4int count = data.fCount;
420 G4double eMean = data.fEmean / count;
421 G4double eMin = data.fEmin;
422 G4double eMax = data.fEmax;
423 //-----> secondary particles flux
424 G4double Eflow = data.fEmean / TotNbofEven
425
426 G4cout << " " << std::setw(13) << name <<
427 << " Emean = " << std::setw(wid) <
428 << G4BestUnit(eMin, "Energy") << "
429 << ") \tEflow/event = " << G4BestUn
430 }
431
432 // particles count inside neuron structure
433 //
434 G4cout << "\n Number of secondary particles
435
436 std::map<G4String, ParticleData>::iterator i
437 for (itn = fParticleDataMap2.begin(); itn !=
438 G4String name = itn->first;
439 ParticleData data = itn->second;
440 G4int count = data.fCount;
441
442 G4cout << " " << std::setw(13) << name <<
443 }
444
445 // molecules count inside neuron
446 // Time cut (from 1 ps to 10 ps) in class I
447 G4cout << "\n Number of molecular products i
448 << "\n time: 1 ps - 10 ps " << G4e
449 // if (1 ns < time <= 10 ns ) MoleculeCount(
450 G4int ind = 0;
451 std::map<G4String, G4int>::iterator itm;
452 for (itm = fMoleculeNumber.begin(); itm != f
453 G4String moleculeName = itm->first;
454 G4int count = itm->second;
455 G4String space = " ";
456 if (++ind % 3 == 0) space = "\n";
457
458 G4cout << " " << std::setw(13) << molecul
459 }
460
461 // compute total Energy and Dose deposited f
462 G4cout << "\n Total energy (MeV) deposition
463
464 G4cout << " " << std::setw(13) << "All volu
465 << " " << std::setw(13) << "Bounding
466 << (fEdepSlice + fEdepNeuron) / MeV <
467 << " " << std::setw(13) << "Neuron:
468 << " " << std::setw(13) << "Soma:
469 << " " << std::setw(13) << "Dendrite
470 << " " << std::setw(13) << "Axon:
471
472 // compute mean Energy and Dose deposited in
473 //
474 // G4AnalysisManager* analys = G4AnalysisMan
475 G4int somaCompHit = 0;
476 G4double somaCompEdep = 0.;
477 G4double somaCompDose = 0.;
478 G4double somaCompEdep2 = 0.;
479 G4double somaCompDose2 = 0.;
480 // Remove old outputs
481 remove("Soma3DEdep.out");
482 for (G4int i = 0; i < fDetector->GetnbSomaco
483 if (fSoma3DEdep[i] > 0.0) {
484 somaCompHit++;
485 somaCompEdep += fSoma3DEdep[i];
486 somaCompEdep2 += fSoma3DEdep[i] * fSoma3
487
488 std::ofstream WriteDataInSoma("Soma3DEde
489 // Index of targeted compartments
490 WriteDataInSoma << fDetector->GetPosSoma
491 << fDetector->GetPosSoma
492 << fDetector->GetPosSoma
493 << " "
494 // Edep in compartments
495 << fSoma3DEdep[i] / keV
496 }
497 }
498 // compute mean Energy and Dose deposited in
499 // +- RMS : Root Mean Square
500 G4double rmsEdepS = 0.;
501 G4double rmsDoseS = 0.;
502 if (somaCompHit > 0) {
503 somaCompEdep /= somaCompHit;
504 somaCompEdep2 /= somaCompHit;
505 rmsEdepS = somaCompEdep2 - somaCompEdep *
506 if (rmsEdepS > 0.)
507 rmsEdepS = std::sqrt(rmsEdepS);
508 else
509 rmsEdepS = 0.;
510 somaCompDose /= somaCompHit;
511 somaCompDose2 /= somaCompHit;
512 rmsDoseS = somaCompDose2 - somaCompDose *
513 if (rmsDoseS > 0.)
514 rmsDoseS = std::sqrt(rmsDoseS);
515 else
516 rmsDoseS = 0.;
517 }
518
519 G4int DendCompHit = 0;
520 G4double DendCompEdep = 0.;
521 G4double DendCompDose = 0.;
522 G4double DendCompEdep2 = 0.;
523 G4double DendCompDose2 = 0.;
524 remove("Dend3DEdep.out");
525 std::ofstream WriteDataInDend("Dend3DEdep.ou
526 for (G4int i = 0; i < fDetector->GetnbDendri
527 if (fDend3DEdep[i] > 0.0) {
528 DendCompHit++;
529 DendCompEdep += fDend3DEdep[i];
530 DendCompEdep2 += fDend3DEdep[i] * fDend3
531
532 WriteDataInDend //<< i+1 <
533 // position of compartm
534 << fDetector->GetPosDendcomp(i).x() <<
535 << '\t' << " " << fDetector->GetPosD
536 << fDetector->GetDistADendSoma(i) << '
537 << " " << fDend3DEdep[i] / keV << '\
538 }
539 }
540 // +- RMS : Root Mean Square
541 G4double rmsEdepD = 0.;
542 G4double rmsDoseD = 0.;
543 if (DendCompHit > 0) {
544 DendCompEdep /= DendCompHit;
545 DendCompEdep2 /= DendCompHit;
546 rmsEdepD = DendCompEdep2 - DendCompEdep *
547 if (rmsEdepD > 0.)
548 rmsEdepD = std::sqrt(rmsEdepD);
549 else
550 rmsEdepD = 0.;
551 DendCompDose /= DendCompHit;
552 DendCompDose2 /= DendCompHit;
553 rmsDoseD = DendCompDose2 - DendCompDose *
554 if (rmsDoseD > 0.)
555 rmsDoseD = std::sqrt(rmsDoseD);
556 else
557 rmsDoseD = 0.;
558 }
559
560 G4int AxonCompHit = 0;
561 G4double AxonCompEdep = 0.;
562 G4double AxonCompDose = 0.;
563 G4double AxonCompEdep2 = 0.;
564 G4double AxonCompDose2 = 0.;
565 remove("Axon3DEdep.out");
566 std::ofstream WriteDataInAxon("Axon3DEdep.ou
567 for (G4int i = 0; i < fDetector->GetnbAxonco
568 if (fAxon3DEdep[i] > 0.0) {
569 AxonCompHit++;
570 AxonCompEdep += fAxon3DEdep[i];
571 AxonCompEdep2 += fAxon3DEdep[i] * fAxon3
572
573 WriteDataInAxon //<< i+1 <
574 // position of compartm
575 << fDetector->GetPosAxoncomp(i).x() <<
576 << '\t' << " " << fDetector->GetPosA
577 << fDetector->GetDistAxonsoma(i) << '\
578 << G4endl;
579 }
580 }
581 // +- RMS : Root Mean Square
582 G4double rmsEdepA = 0.;
583 G4double rmsDoseA = 0.;
584 if (AxonCompHit > 0) {
585 AxonCompEdep /= AxonCompHit;
586 AxonCompEdep2 /= AxonCompHit;
587 rmsEdepA = AxonCompEdep2 - AxonCompEdep *
588 if (rmsEdepA > 0.)
589 rmsEdepA = std::sqrt(rmsEdepA);
590 else
591 rmsEdepA = 0.;
592 AxonCompDose /= AxonCompHit;
593 AxonCompDose2 /= AxonCompHit;
594 rmsDoseA = AxonCompDose2 - AxonCompDose *
595 if (rmsDoseA > 0.)
596 rmsDoseA = std::sqrt(rmsDoseA);
597 else
598 rmsDoseA = 0.;
599 }
600
601 G4cout << "\n Number of compartments travers
602 G4cout << " " << std::setw(13) << "Soma: "
603 << " of total: " << fDetector->GetnbS
604 << " " << std::setw(13) << "Dendrite
605 << " of total: " << fDetector->GetnbD
606 << " " << std::setw(13) << "Axon: "
607 << " of total: " << fDetector->GetnbA
608 G4cout << "\n Dendritic (or Axon) compartmen
609 << " at the distance from Soma have b
610 << "\n Dend3DEdep.out, Axon3DEdep.out
611 << "\n Outputs of energy deposition p
612 << "\n OutputPerEvent.out" << G4endl;
613
614 // remove all contents in fProcCounter, fCou
615 fProcCounter.clear();
616 fParticleDataMap2.clear();
617
618 // restore default format
619 G4cout.precision(dfprec);
620 }
621
622 //....oooOO0OOooo........oooOO0OOooo........oo
623