| 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 // This example is provided by the Geant4-DNA collaboration
27 // Any report or published results obtained using the Geant4-DNA software
28 // shall cite the following Geant4-DNA collaboration publication:
29 // Med. Phys. 37 (2010) 4692-4708
30 // and papers
31 // M. Batmunkh et al. J Radiat Res Appl Sci 8 (2015) 498-507
32 // O. Belov et al. Physica Medica 32 (2016) 1510-1520
33 // The Geant4-DNA web site is available at http://geant4-dna.org
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........oooOO0OOooo........oooOO0OOooo......
62
63 Run::Run(DetectorConstruction* det) : G4Run(), fDetector(det), fParticle(0)
64 {
65 fSoma3DEdep = new G4double[fDetector->GetnbSomacomp()];
66 for (G4int i = 0; i < fDetector->GetnbSomacomp(); i++) {
67 fSoma3DEdep[i] = 0.;
68 }
69 fDend3DEdep = new G4double[fDetector->GetnbDendritecomp()];
70 for (G4int i = 0; i < fDetector->GetnbDendritecomp(); i++) {
71 fDend3DEdep[i] = 0.;
72 }
73 fAxon3DEdep = new G4double[fDetector->GetnbAxoncomp()];
74 for (G4int i = 0; i < fDetector->GetnbAxoncomp(); i++) {
75 fAxon3DEdep[i] = 0.;
76 }
77
78 fEdepAll = fEdepAll_err = fEdepMedium = fEdepMedium_err = fEdepSlice = fEdepSlice_err =
79 fEdepSoma = fEdepSoma_err = 0.;
80 fEdepDend = fEdepDend_err = fEdepAxon = fEdepAxon_err = fEdepNeuron = fEdepNeuron_err = 0.;
81 }
82
83 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
84
85 Run::~Run()
86 {
87 delete[] fSoma3DEdep;
88 delete[] fDend3DEdep;
89 delete[] fAxon3DEdep;
90 }
91
92 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
93
94 void Run::SetPrimary(G4ParticleDefinition* particle, G4double energy)
95 {
96 fParticle = particle;
97 fEkin = energy;
98 }
99
100 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
101
102 void Run::AddPrimaryLET(G4double let)
103 {
104 fLET += let;
105 fLET2 += let * let;
106 }
107
108 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
109
110 void Run::SetTrackLength(G4double t)
111 {
112 ftrackLength = t;
113 fTrackLen += t;
114 fTrackLen2 += t * t;
115 }
116
117 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
118
119 void Run::CountProcesses(const G4VProcess* process)
120 {
121 G4String procName = process->GetProcessName();
122 std::map<G4String, G4int>::iterator it = fProcCounter.find(procName);
123 if (it == fProcCounter.end()) {
124 fProcCounter[procName] = 1;
125 }
126 else {
127 fProcCounter[procName]++;
128 }
129 }
130
131 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
132
133 void Run::ParticleCount(G4String name, G4double Ekin)
134 {
135 std::map<G4String, ParticleData>::iterator it = fParticleDataMap1.find(name);
136 if (it == fParticleDataMap1.end()) {
137 fParticleDataMap1[name] = ParticleData(1, Ekin, Ekin, Ekin);
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........oooOO0OOooo........oooOO0OOooo......
152
153 void Run::ParticleCountNeuron(G4String name, G4double Ekin)
154 {
155 std::map<G4String, ParticleData>::iterator it = fParticleDataMap2.find(name);
156 if (it == fParticleDataMap2.end()) {
157 fParticleDataMap2[name] = ParticleData(1, Ekin, Ekin, Ekin);
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........oooOO0OOooo........oooOO0OOooo......
172 /*
173 void Run::MoleculeCount(G4String, G4double)
174 {
175 //fMoleculeNumber = G4MoleculeCounter::Instance()
176 // ->GetNMoleculesAtTime(moleculename, Gtime);
177 }
178 */
179 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
180
181 void Run::MoleculeCountNeuron(G4Molecule* molecule)
182 {
183 G4String moleculename = molecule->GetName();
184 std::map<G4String, G4int>::iterator it = fMoleculeNumber.find(moleculename);
185 if (it == fMoleculeNumber.end()) {
186 fMoleculeNumber[moleculename] = 1;
187 }
188 else {
189 fMoleculeNumber[moleculename]++;
190 }
191 }
192
193 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
194
195 void Run::AddEflow(G4double eflow)
196 {
197 fEnergyFlow += eflow;
198 fEnergyFlow2 += eflow * eflow;
199 }
200
201 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
202
203 void Run::Merge(const G4Run* run)
204 {
205 const Run* localRun = static_cast<const Run*>(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 itp;
219 for (itp = localRun->fProcCounter.begin(); itp != localRun->fProcCounter.end(); ++itp) {
220 G4String procName = itp->first;
221 G4int localCount = itp->second;
222 if (fProcCounter.find(procName) == fProcCounter.end()) {
223 fProcCounter[procName] = localCount;
224 }
225 else {
226 fProcCounter[procName] += localCount;
227 }
228 }
229
230 // map: created particles count outside neuron structure
231 std::map<G4String, ParticleData>::const_iterator itc;
232 for (itc = localRun->fParticleDataMap1.begin(); itc != localRun->fParticleDataMap1.end(); ++itc) {
233 G4String name = itc->first;
234 const ParticleData& localData = itc->second;
235 if (fParticleDataMap1.find(name) == fParticleDataMap1.end()) {
236 fParticleDataMap1[name] =
237 ParticleData(localData.fCount, localData.fEmean, localData.fEmin, localData.fEmax);
238 }
239 else {
240 ParticleData& data = fParticleDataMap1[name];
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 neuron structure
251 std::map<G4String, ParticleData>::const_iterator itn;
252 for (itn = localRun->fParticleDataMap2.begin(); itn != localRun->fParticleDataMap2.end(); ++itn) {
253 G4String name = itn->first;
254 const ParticleData& localData = itn->second;
255 if (fParticleDataMap2.find(name) == fParticleDataMap2.end()) {
256 fParticleDataMap2[name] =
257 ParticleData(localData.fCount, localData.fEmean, localData.fEmin, localData.fEmax);
258 }
259 else {
260 ParticleData& data = fParticleDataMap2[name];
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 itm;
271 for (itm = localRun->fMoleculeNumber.begin(); itm != localRun->fMoleculeNumber.end(); ++itm) {
272 G4String moleculeName = itm->first;
273 G4int localCount = itm->second;
274 if (fMoleculeNumber.find(moleculeName) == fMoleculeNumber.end()) {
275 fMoleculeNumber[moleculeName] = localCount;
276 }
277 else {
278 fMoleculeNumber[moleculeName] += localCount;
279 }
280 }
281
282 // hits compartments in neuron compartments
283 //
284 for (G4int i = 0; i < fDetector->GetnbSomacomp(); i++) {
285 fSoma3DEdep[i] += localRun->fSoma3DEdep[i];
286 }
287 for (G4int i = 0; i < fDetector->GetnbDendritecomp(); i++) {
288 fDend3DEdep[i] += localRun->fDend3DEdep[i];
289 }
290 for (G4int i = 0; i < fDetector->GetnbAxoncomp(); i++) {
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........oooOO0OOooo........oooOO0OOooo......
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->GetParticleName();
326 G4double GunArea;
327
328 G4Material* material = fDetector->GetTargetMaterial();
329
330 // compute track length of primary track
331 //
332 fTrackLen /= numberOfEvent;
333 fTrackLen2 /= numberOfEvent;
334 G4double rms = fTrackLen2 - fTrackLen * fTrackLen;
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 * EdepTotal;
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(fEkin, fParticle, material);
356 }
357
358 // Stopping Power from simulation.
359 G4double meandEdx = (EdepTotal / fTrackLen) / (keV / um);
360 G4double meandEdxerr = (rmst / fTrackLen) / (keV / um);
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 %f", &tmp, &tmp, &tmp, &tmp, &En, &tmp, &tmp, &tmp);
367 if (ncols < 0) break;
368 if (En > 0) Ntrav++;
369 }
370 fclose(fp);
371 // The surface area is calculated as spherical medium
372 GunArea = fDetector->GetTotSurfMedium();
373 // Fluence dose of single paticle track
374 G4double FluenceDoseperBeam = 0.160218 * (dEdxFull) / (GunArea * std::pow(10, 18));
375
376 G4cout << "\n ======= The summary of simulation results 'neuron' ========\n";
377 G4cout << "\n Primary particle = " << Particle
378 << "\n Kinetic energy of beam = " << fEkin / MeV << " A*MeV"
379 << "\n Particle traversals the neuron = " << Ntrav << " of " << numberOfEvent
380 << "\n Full LET of beam as formulas = " << dEdxFull / (keV / um) << " keV/um"
381 << "\n Mean LET of beam as simulation = " << meandEdx << " +- " << meandEdxerr
382 << " keV/um"
383 << "\n Mean track length of beam = " << fTrackLen / um << " +- " << rms << " um"
384 << "\n Particle fluence = " << numberOfEvent / (GunArea / (cm * cm))
385 << " particles/cm^2"
386 << "\n Fluence dose (full) = "
387 << numberOfEvent * FluenceDoseperBeam / (joule / kg) << " Gy"
388 << "\n Fluence dose ber beam = " << FluenceDoseperBeam / (joule / kg) << " Gy"
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 process:" << G4endl;
399
400 G4int index = 0;
401 std::map<G4String, G4int>::iterator it;
402 for (it = fProcCounter.begin(); it != fProcCounter.end(); it++) {
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 << "=" << std::setw(7) << count << space;
408 }
409 G4cout << G4endl;
410
411 // particles count outside neuron structure
412 //
413 G4cout << "\n List of generated particles outside neuron structure:" << G4endl;
414
415 std::map<G4String, ParticleData>::iterator itc;
416 for (itc = fParticleDataMap1.begin(); itc != fParticleDataMap1.end(); 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 / TotNbofEvents;
425
426 G4cout << " " << std::setw(13) << name << " : " << std::setw(7) << count
427 << " Emean = " << std::setw(wid) << G4BestUnit(eMean, "Energy") << "\t( "
428 << G4BestUnit(eMin, "Energy") << " --> " << G4BestUnit(eMax, "Energy")
429 << ") \tEflow/event = " << G4BestUnit(Eflow, "Energy") << G4endl;
430 }
431
432 // particles count inside neuron structure
433 //
434 G4cout << "\n Number of secondary particles inside neuron structure:" << G4endl;
435
436 std::map<G4String, ParticleData>::iterator itn;
437 for (itn = fParticleDataMap2.begin(); itn != fParticleDataMap2.end(); itn++) {
438 G4String name = itn->first;
439 ParticleData data = itn->second;
440 G4int count = data.fCount;
441
442 G4cout << " " << std::setw(13) << name << " : " << std::setw(7) << count << G4endl;
443 }
444
445 // molecules count inside neuron
446 // Time cut (from 1 ps to 10 ps) in class ITTrackingAction.cc
447 G4cout << "\n Number of molecular products inside neuron structure:"
448 << "\n time: 1 ps - 10 ps " << G4endl;
449 // if (1 ns < time <= 10 ns ) MoleculeCount(molname, time) ;
450 G4int ind = 0;
451 std::map<G4String, G4int>::iterator itm;
452 for (itm = fMoleculeNumber.begin(); itm != fMoleculeNumber.end(); itm++) {
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) << moleculeName << " : " << std::setw(7) << count << G4endl;
459 }
460
461 // compute total Energy and Dose deposited for all events
462 G4cout << "\n Total energy (MeV) deposition :" << G4endl;
463
464 G4cout << " " << std::setw(13) << "All volume: " << std::setw(7) << fEdepAll / MeV << "\n "
465 << " " << std::setw(13) << "Bounding slice: " << std::setw(7)
466 << (fEdepSlice + fEdepNeuron) / MeV << "\n "
467 << " " << std::setw(13) << "Neuron: " << std::setw(7) << fEdepNeuron / MeV << "\n "
468 << " " << std::setw(13) << "Soma: " << std::setw(7) << fEdepSoma / MeV << "\n "
469 << " " << std::setw(13) << "Dendrites: " << std::setw(7) << fEdepDend / MeV << "\n "
470 << " " << std::setw(13) << "Axon: " << std::setw(7) << fEdepAxon / MeV << G4endl;
471
472 // compute mean Energy and Dose deposited in hit compartments
473 //
474 // G4AnalysisManager* analys = G4AnalysisManager::Instance();
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->GetnbSomacomp(); i++) {
483 if (fSoma3DEdep[i] > 0.0) {
484 somaCompHit++;
485 somaCompEdep += fSoma3DEdep[i];
486 somaCompEdep2 += fSoma3DEdep[i] * fSoma3DEdep[i];
487
488 std::ofstream WriteDataInSoma("Soma3DEdep.out", std::ios::app);
489 // Index of targeted compartments
490 WriteDataInSoma << fDetector->GetPosSomacomp(i).x() << '\t' << " "
491 << fDetector->GetPosSomacomp(i).y() << '\t' << " "
492 << fDetector->GetPosSomacomp(i).z() << '\t'
493 << " "
494 // Edep in compartments
495 << fSoma3DEdep[i] / keV << '\t' << " " << G4endl;
496 }
497 }
498 // compute mean Energy and Dose deposited in compartments;
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 * 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 * 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.out", std::ios::app);
526 for (G4int i = 0; i < fDetector->GetnbDendritecomp(); i++) {
527 if (fDend3DEdep[i] > 0.0) {
528 DendCompHit++;
529 DendCompEdep += fDend3DEdep[i];
530 DendCompEdep2 += fDend3DEdep[i] * fDend3DEdep[i];
531
532 WriteDataInDend //<< i+1 << '\t' << " "
533 // position of compartments
534 << fDetector->GetPosDendcomp(i).x() << '\t' << " " << fDetector->GetPosDendcomp(i).y()
535 << '\t' << " " << fDetector->GetPosDendcomp(i).z() << '\t' << " "
536 << fDetector->GetDistADendSoma(i) << '\t' << " " << fDetector->GetDistBDendSoma(i) << '\t'
537 << " " << fDend3DEdep[i] / keV << '\t' << " " << G4endl;
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 * 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 * 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.out", std::ios::app);
567 for (G4int i = 0; i < fDetector->GetnbAxoncomp(); i++) {
568 if (fAxon3DEdep[i] > 0.0) {
569 AxonCompHit++;
570 AxonCompEdep += fAxon3DEdep[i];
571 AxonCompEdep2 += fAxon3DEdep[i] * fAxon3DEdep[i];
572
573 WriteDataInAxon //<< i+1 << '\t' << " "
574 // position of compartments
575 << fDetector->GetPosAxoncomp(i).x() << '\t' << " " << fDetector->GetPosAxoncomp(i).y()
576 << '\t' << " " << fDetector->GetPosAxoncomp(i).z() << '\t' << " "
577 << fDetector->GetDistAxonsoma(i) << '\t' << " " << fAxon3DEdep[i] / keV << '\t' << " "
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 * 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 * AxonCompDose;
595 if (rmsDoseA > 0.)
596 rmsDoseA = std::sqrt(rmsDoseA);
597 else
598 rmsDoseA = 0.;
599 }
600
601 G4cout << "\n Number of compartments traversed by particle tracks :" << G4endl;
602 G4cout << " " << std::setw(13) << "Soma: " << std::setw(7) << somaCompHit
603 << " of total: " << fDetector->GetnbSomacomp() << "\n "
604 << " " << std::setw(13) << "Dendrites: " << std::setw(7) << DendCompHit
605 << " of total: " << fDetector->GetnbDendritecomp() << "\n "
606 << " " << std::setw(13) << "Axon: " << std::setw(7) << AxonCompHit
607 << " of total: " << fDetector->GetnbAxoncomp() << "\n " << G4endl;
608 G4cout << "\n Dendritic (or Axon) compartmental energy deposits \n"
609 << " at the distance from Soma have been written into *.out files:"
610 << "\n Dend3DEdep.out, Axon3DEdep.out, Soma3DEdep.out"
611 << "\n Outputs of energy deposition per event written in data file:"
612 << "\n OutputPerEvent.out" << G4endl;
613
614 // remove all contents in fProcCounter, fCount
615 fProcCounter.clear();
616 fParticleDataMap2.clear();
617
618 // restore default format
619 G4cout.precision(dfprec);
620 }
621
622 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
623