ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
Run.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file Run.cc
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 //
41 
42 #include "Run.hh"
43 #include "DetectorConstruction.hh"
44 #include "PrimaryGeneratorAction.hh"
45 #include "Analysis.hh"
46 //#include "NeuronLoadDataFile.hh"
47 #include "G4UnitsTable.hh"
48 #include "G4SystemOfUnits.hh"
49 #include "G4ios.hh"
50 #include "G4Molecule.hh"
51 #include "G4MoleculeCounter.hh"
52 #include "G4MoleculeGun.hh"
53 #include "G4H2O.hh"
54 #include <G4Scheduler.hh>
55 #include "G4MoleculeTable.hh"
56 #include <cmath>
57 #include "G4EmCalculator.hh"
58 #include <iomanip>
59 
60 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
61 
62 Run::Run(DetectorConstruction* det)
63 : G4Run(), fDetector(det),
64  fParticle(0), fEkin(0.),
65  fLET(0.), fLET2(0.),
66  fTrackLen(0.), fTrackLen2(0.)
67 {
68  //fNeuronLoadParamz = new NeuronLoadDataFile();
69 
70  fSoma3DEdep = new G4double[fDetector->GetnbSomacomp()];
71  for (G4int i=0; i<fDetector->GetnbSomacomp(); i++)
72  {
73  fSoma3DEdep[i]=0.;
74  }
75  fDend3DEdep = new G4double[fDetector->GetnbDendritecomp()];
76  for (G4int i=0; i<fDetector->GetnbDendritecomp(); i++)
77  {
78  fDend3DEdep[i]=0.;
79  }
80  fAxon3DEdep = new G4double[fDetector->GetnbAxoncomp()];
81  for (G4int i=0; i<fDetector->GetnbAxoncomp(); i++)
82  {
83  fAxon3DEdep[i]=0.;
84  }
85 
86  fEdepAll = fEdepAll_err=fEdepMedium= fEdepMedium_err= fEdepSlice=
87  fEdepSlice_err=fEdepSoma=fEdepSoma_err=0. ;
88  fEdepDend = fEdepDend_err=fEdepAxon= fEdepAxon_err=fEdepNeuron=
89  fEdepNeuron_err =0. ;
90  fEnergyFlow = fEnergyFlow2 = 0.;
91 
92 }
93 
94 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
95 
96 Run::~Run()
97 {
98  delete[] fSoma3DEdep;
99  delete[] fDend3DEdep;
100  delete[] fAxon3DEdep;
101  }
102 
103 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
104 
105 void Run::SetPrimary(G4ParticleDefinition* particle, G4double energy)
106 {
107  fParticle = particle;
108  fEkin = energy;
109 }
110 
111 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
112 
113 void Run::AddPrimaryLET(G4double let)
114 {
115  fLET += let;
116  fLET2 += let*let;
117 }
118 
119 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
120 
121 void Run::SetTrackLength(G4double t)
122 {
123  ftrackLength = t;
124  fTrackLen += t;
125  fTrackLen2 += t*t;
126 }
127 
128 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
129 
130 void Run::CountProcesses(const G4VProcess* process)
131 {
132  G4String procName = process->GetProcessName();
133  std::map<G4String,G4int>::iterator it = fProcCounter.find(procName);
134  if ( it == fProcCounter.end()) {
135  fProcCounter[procName] = 1;
136  }
137  else {
138  fProcCounter[procName]++;
139  }
140 }
141 
142 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
143 
144 void Run::ParticleCount(G4String name, G4double Ekin)
145 {
146  std::map<G4String, ParticleData>::iterator it = fParticleDataMap1.find(name);
147  if ( it == fParticleDataMap1.end()) {
148  fParticleDataMap1[name] = ParticleData(1, Ekin, Ekin, Ekin);
149  }
150  else {
151  ParticleData& data = it->second;
152  data.fCount++;
153  data.fEmean += Ekin;
154  //update min max
155  G4double emin = data.fEmin;
156  if (Ekin < emin) data.fEmin = Ekin;
157  G4double emax = data.fEmax;
158  if (Ekin > emax) data.fEmax = Ekin;
159  }
160 }
161 
162 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
163 
164 void Run::ParticleCountNeuron(G4String name, G4double Ekin)
165 {
166  std::map<G4String, ParticleData>::iterator it = fParticleDataMap2.find(name);
167  if ( it == fParticleDataMap2.end()) {
168  fParticleDataMap2[name] = ParticleData(1, Ekin, Ekin, Ekin);
169  }
170  else {
171  ParticleData& data = it->second;
172  data.fCount++;
173  data.fEmean += Ekin;
174  //update min max
175  G4double emin = data.fEmin;
176  if (Ekin < emin) data.fEmin = Ekin;
177  G4double emax = data.fEmax;
178  if (Ekin > emax) data.fEmax = Ekin;
179  }
180 }
181 
182 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
183 
184 void Run::MoleculeCount(G4String, G4double)
185 {
186 //fMoleculeNumber = G4MoleculeCounter::Instance()
187 // ->GetNMoleculesAtTime(moleculename, Gtime);
188 }
189 
190 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
191 
192 void Run::MoleculeCountNeuron(G4Molecule* molecule)
193 {
194  G4String moleculename = molecule->GetName();
195  std::map<G4String,G4int>::iterator it = fMoleculeNumber.find(moleculename);
196  if ( it == fMoleculeNumber.end()) {
197  fMoleculeNumber[moleculename] = 1;
198  }
199  else {
200  fMoleculeNumber[moleculename]++;
201  }
202 }
203 
204 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
205 
206 void Run::AddEflow(G4double eflow)
207 {
208  fEnergyFlow += eflow;
209  fEnergyFlow2 += eflow*eflow;
210 }
211 
212 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
213 
214 void Run::Merge(const G4Run* run)
215 {
216  const Run* localRun = static_cast<const Run*>(run);
217 
218  //primary particle info
219  //
220  fParticle = localRun->fParticle;
221  fEkin = localRun->fEkin;
222  ftrackLength = localRun->ftrackLength;
223  fTrackLen += localRun->fTrackLen;
224  fTrackLen2 += localRun->fTrackLen2;
225  fLET += localRun->fLET;
226  fLET2 += localRun->fLET2;
227 
228  //map: processes count in simulation medium
229  std::map<G4String,G4int>::const_iterator itp;
230  for ( itp = localRun->fProcCounter.begin();
231  itp != localRun->fProcCounter.end(); ++itp ) {
232 
233  G4String procName = itp->first;
234  G4int localCount = itp->second;
235  if ( fProcCounter.find(procName) == fProcCounter.end()) {
236  fProcCounter[procName] = localCount;
237  }
238  else {
239  fProcCounter[procName] += localCount;
240  }
241  }
242 
243  //map: created particles count outside neuron structure
244  std::map<G4String,ParticleData>::const_iterator itc;
245  for (itc = localRun->fParticleDataMap1.begin();
246  itc != localRun->fParticleDataMap1.end(); ++itc) {
247 
248  G4String name = itc->first;
249  const ParticleData& localData = itc->second;
250  if ( fParticleDataMap1.find(name) == fParticleDataMap1.end()) {
251  fParticleDataMap1[name]
252  = ParticleData(localData.fCount,
253  localData.fEmean,
254  localData.fEmin,
255  localData.fEmax);
256  }
257  else {
258  ParticleData& data = fParticleDataMap1[name];
259  data.fCount += localData.fCount;
260  data.fEmean += localData.fEmean;
261  G4double emin = localData.fEmin;
262  if (emin < data.fEmin) data.fEmin = emin;
263  G4double emax = localData.fEmax;
264  if (emax > data.fEmax) data.fEmax = emax;
265  }
266  }
267 
268  //map: created particles count inside neuron structure
269  std::map<G4String,ParticleData>::const_iterator itn;
270  for (itn = localRun->fParticleDataMap2.begin();
271  itn != localRun->fParticleDataMap2.end(); ++itn) {
272 
273  G4String name = itn->first;
274  const ParticleData& localData = itn->second;
275  if ( fParticleDataMap2.find(name) == fParticleDataMap2.end()) {
276  fParticleDataMap2[name]
277  = ParticleData(localData.fCount,
278  localData.fEmean,
279  localData.fEmin,
280  localData.fEmax);
281  }
282  else {
283  ParticleData& data = fParticleDataMap2[name];
284  data.fCount += localData.fCount;
285  data.fEmean += localData.fEmean;
286  G4double emin = localData.fEmin;
287  if (emin < data.fEmin) data.fEmin = emin;
288  G4double emax = localData.fEmax;
289  if (emax > data.fEmax) data.fEmax = emax;
290  }
291  }
292 
293  //map: molecule count
294  //fMoleculeNumber += localRun->fMoleculeNumber;
295 
296  std::map<G4String,G4int>::const_iterator itm;
297  for ( itm = localRun->fMoleculeNumber.begin();
298  itm != localRun->fMoleculeNumber.end(); ++itm ) {
299 
300  G4String moleculeName = itm->first;
301  G4int localCount = itm->second;
302  if ( fMoleculeNumber.find(moleculeName) == fMoleculeNumber.end()) {
303  fMoleculeNumber[moleculeName] = localCount;
304  }
305  else {
306  fMoleculeNumber[moleculeName] += localCount;
307  }
308  }
309 
310  // hits compartments in neuron compartments
311  //
312  for (G4int i=0; i<fDetector->GetnbSomacomp(); i++)
313  {
314  fSoma3DEdep[i] += localRun->fSoma3DEdep[i];
315  }
316  for (G4int i=0; i<fDetector->GetnbDendritecomp(); i++)
317  {
318  fDend3DEdep[i] +=localRun->fDend3DEdep[i];
319  }
320  for (G4int i=0; i<fDetector->GetnbAxoncomp(); i++)
321  {
322  fAxon3DEdep[i] +=localRun->fAxon3DEdep[i];
323  }
324 
325  // accumulate sums
326  //
327  fEdepAll += localRun->fEdepAll; fEdepAll_err += localRun->fEdepAll_err;
328  fEdepMedium += localRun->fEdepMedium;
329  fEdepMedium_err += localRun->fEdepMedium_err;
330  fEdepSlice += localRun->fEdepSlice;
331  fEdepSlice_err += localRun->fEdepSlice_err;
332  fEdepSoma += localRun->fEdepSoma; fEdepSoma_err += localRun->fEdepSoma_err;
333  fEdepDend += localRun->fEdepDend; fEdepDend_err += localRun->fEdepDend_err;
334  fEdepAxon += localRun->fEdepAxon; fEdepAxon_err+= localRun->fEdepAxon_err;
335  fEdepNeuron += localRun->fEdepNeuron;
336  fEdepNeuron_err += localRun->fEdepNeuron_err;
337 
338  fEnergyFlow += localRun->fEnergyFlow;
339  fEnergyFlow2 += localRun->fEnergyFlow2;
340 
341  G4Run::Merge(run);
342 }
343 
344 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
345 
346 void Run::EndOfRun()
347 {
348  G4int prec = 5, wid = prec + 2;
349  G4int dfprec = G4cout.precision(prec);
350 
351  // Characteristics of Primary particle
352  G4String Particle = fParticle->GetParticleName();
353  G4double GunArea ;
354  //GunArea = fParticle->GetGunArea();
355  G4Material* material = fDetector->GetTargetMaterial();
356  //G4double density = material->GetDensity();
357  //compute track length of primary track
358  //
359  fTrackLen /= numberOfEvent; fTrackLen2 /= numberOfEvent;
360  G4double rms = fTrackLen2 - fTrackLen*fTrackLen;
361  if (rms>0.) rms = std::sqrt(rms); else rms = 0.;
362 
363  G4int TotNbofEvents = numberOfEvent;
364  G4double EdepTotal = fEdepAll;
365  G4double EdepTotal2 = fEdepAll_err;
366  EdepTotal /= TotNbofEvents; EdepTotal2 /= TotNbofEvents;
367  G4double rmst = EdepTotal2 - EdepTotal*EdepTotal;
368  if (rmst>0.) rmst = std::sqrt(rmst);
369  else rmst = 0.;
370 
371  //Stopping Power from input Table.
372  G4EmCalculator emCalculator;
373  //G4double dEdxTable = 0.,
374  G4double dEdxFull = 0.;
375  if (fParticle->GetPDGCharge()!= 0.) {
376  // dEdxTable = emCalculator.GetDEDX(fEkin,fParticle,material);
377  dEdxFull = emCalculator.ComputeTotalDEDX(fEkin,fParticle,material);
378  }
379 
380  //Stopping Power from simulation.
381  G4double meandEdx = (EdepTotal/fTrackLen)/(keV/um);
382  G4double meandEdxerr = (rmst/fTrackLen)/(keV/um);
383  //G4double stopPower = (meandEdx/density)/(MeV*cm2/g);
384  G4int ncols=0;
385  G4int nlines = 0;
386  G4float tmp, En;
387  G4int Ntrav = 0;
388  FILE * fp = fopen("OutputPerEvent.out","r");
389  while (1) {
390  ncols = fscanf(fp," %f %f %f %f %f %f %f %f",
391  &tmp, &tmp, &tmp, &tmp, &En, &tmp, &tmp, &tmp);
392  if (ncols < 0) break;
393  if (En>0) Ntrav++;
394  nlines++;}
395  fclose(fp);
396  // The surface area is calculated as spherical medium
397  GunArea = fDetector->GetTotSurfMedium();
398  // Fluence dose of single paticle track
399  G4double FluenceDoseperBeam = 0.160218*(dEdxFull)/(GunArea*std::pow(10,18)) ;
400 
401  G4cout << "\n ======= The summary of simulation results 'neuron' ========\n";
402  G4cout
403  << "\n Primary particle = " << Particle
404  << "\n Kinetic energy of beam = " << fEkin/MeV<<" A*MeV"
405  << "\n Particle traversals the neuron = " << Ntrav<<" of "<<numberOfEvent
406  << "\n Full LET of beam as formulas = " <<dEdxFull/(keV/um) << " keV/um"
407  << "\n Mean LET of beam as simulation = "
408  << meandEdx << " +- " << meandEdxerr << " keV/um"
409  << "\n Mean track length of beam = "
410  << fTrackLen/um << " +- " << rms << " um"
411  << "\n Particle fluence = "
412  << numberOfEvent/(GunArea/(cm*cm))<<" particles/cm^2"
413  << "\n Fluence dose (full) = "
414  << numberOfEvent*FluenceDoseperBeam/(joule/kg)<<" Gy"
415  << "\n Fluence dose ber beam = "
416  << FluenceDoseperBeam/(joule/kg) <<" Gy" << G4endl;
417 
418  if (numberOfEvent == 0) { G4cout.precision(dfprec); return;}
419 
420  //frequency of processes in all volume
421  //
422  G4cout << "\n List of generated physical process:" << G4endl;
423 
424  G4int index = 0;
425  std::map<G4String,G4int>::iterator it;
426  for (it = fProcCounter.begin(); it != fProcCounter.end(); it++) {
427  G4String procName = it->first;
428  G4int count = it->second;
429  G4String space = " "; if (++index%1 == 0) space = "\n";
430  G4cout << " " << std::setw(20) << procName << "="<< std::setw(7) << count
431  << space;
432  }
433  G4cout << G4endl;
434 
435  //particles count outside neuron structure
436  //
437  G4cout << "\n List of generated particles outside neuron structure:"
438  << G4endl;
439 
440  std::map<G4String,ParticleData>::iterator itc;
441  for (itc = fParticleDataMap1.begin(); itc != fParticleDataMap1.end(); itc++) {
442  G4String name = itc->first;
443  ParticleData data = itc->second;
444  G4int count = data.fCount;
445  G4double eMean = data.fEmean/count;
446  G4double eMin = data.fEmin;
447  G4double eMax = data.fEmax;
448  //-----> secondary particles flux
449  G4double Eflow = data.fEmean/TotNbofEvents;
450 
451  G4cout << " " << std::setw(13) << name << " : " << std::setw(7) << count
452  << " Emean = " << std::setw(wid) << G4BestUnit(eMean, "Energy")
453  << "\t( " << G4BestUnit(eMin, "Energy")
454  << " --> " << G4BestUnit(eMax, "Energy")
455  << ") \tEflow/event = " << G4BestUnit(Eflow, "Energy")
456  << G4endl;
457  }
458 
459  //particles count inside neuron structure
460  //
461  G4cout << "\n Number of secondary particles inside neuron structure:"
462  << G4endl;
463 
464  std::map<G4String,ParticleData>::iterator itn;
465  for (itn = fParticleDataMap2.begin(); itn != fParticleDataMap2.end(); itn++) {
466  G4String name = itn->first;
467  ParticleData data = itn->second;
468  G4int count = data.fCount;
469  //G4double eMean = data.fEmean/count;
470  //G4double eMin = data.fEmin;
471  //G4double eMax = data.fEmax;
472  //-----> secondary particles flux
473  //G4double Eflow = data.fEmean/TotNbofEvents;
474 
475  G4cout << " " << std::setw(13) << name << " : " << std::setw(7) << count
476  //<< " Emean = " << std::setw(wid) << G4BestUnit(eMean, "Energy")
477  //<< "\t( " << G4BestUnit(eMin, "Energy")
478  //<< " --> " << G4BestUnit(eMax, "Energy")
479  //<< ") \tEflow/event = " << G4BestUnit(Eflow, "Energy")
480  << G4endl;
481  }
482 
483  //molecules count inside neuron
484  // Time cut (from 1 ps to 10 ps) in class ITTrackingAction.cc
485  G4cout << "\n Number of molecular products inside neuron structure:"
486  << "\n time: 1 ps - 10 ps "<< G4endl;
487  // if (1 ns < time <= 10 ns ) MoleculeCount(molname, time) ;
488  G4int ind = 0;
489  std::map<G4String,G4int>::iterator itm;
490  for (itm = fMoleculeNumber.begin(); itm != fMoleculeNumber.end(); itm++) {
491  G4String moleculeName = itm->first;
492  G4int count = itm->second;
493  G4String space = " "; if (++ind%3 == 0) space = "\n";
494 
495  G4cout << " " << std::setw(13) << moleculeName << " : " << std::setw(7)
496  << count << G4endl;
497  }
498 
499  // compute total Energy and Dose deposited for all events
500  //
501  // EdepMedum + EdepSlice + EdepNeuron --> EdepAll
502  // EdepSoma + EdepDend + EdepAxon + EdepSpines --> EdepNeuron
503  G4cout << "\n Total energy (MeV) deposition :" << G4endl;
504 
505  G4cout << " "
506  << std::setw(13) << "All volume: " << std::setw(7) << fEdepAll/MeV<< "\n "
507  << " " << std::setw(13) << "Bounding slice: "
508  << std::setw(7) << (fEdepSlice+fEdepNeuron)/MeV << "\n "
509  << " " << std::setw(13) << "Neuron: " << std::setw(7)
510  << fEdepNeuron/MeV << "\n "
511  << " " << std::setw(13) << "Soma: " << std::setw(7)
512  << fEdepSoma/MeV<< "\n "
513  << " " << std::setw(13) << "Dendrites: " << std::setw(7)
514  << fEdepDend/MeV << "\n "
515  << " " << std::setw(13) << "Axon: " << std::setw(7)
516  << fEdepAxon/MeV
517  << G4endl;
518 
519  // compute mean Energy and Dose deposited in hit compartments
520  //
521  //G4AnalysisManager* analys = G4AnalysisManager::Instance();
522  G4int somaCompHit = 0;
523  G4double somaCompEdep = 0.;
524  G4double somaCompDose = 0.;
525  G4double somaCompEdep2 = 0.;
526  G4double somaCompDose2 = 0.;
527  // Remove old outputs
528  remove ("Soma3DEdep.out");
529  for (G4int i=0; i<fDetector->GetnbSomacomp(); i++)
530  {
531  if (fSoma3DEdep[i] > 0.0)
532  {
533  somaCompHit ++;
534  somaCompEdep += fSoma3DEdep[i] ;
535  //somaCompDose += fSoma3DEdep[i]/fDetector->GetMassSomacomp(i) ;
536  somaCompEdep2 += fSoma3DEdep[i]*fSoma3DEdep[i] ;
537  //somaCompDose2 += (fSoma3DEdep[i]/fDetector->GetMassSomacomp(i))*
538  // (fSoma3DEdep[i]/fDetector->GetMassSomacomp(i));
539  /*G4double distSoma = 0.;
540  //Fill ntuple #1
541  analys->FillNtupleDColumn(1,0,i+1);
542  analys->FillNtupleDColumn(1,1,distSoma);
543  analys->FillNtupleDColumn(1,2,fSoma3DEdep[i]/keV);
544  analys->FillNtupleDColumn(1,3,(fSoma3DEdep[i]/joule)/
545  (fDetector->GetMassSomacomp(i)/kg));
546  analys->AddNtupleRow(1);
547  */
548 
549  std::ofstream WriteDataInSoma("Soma3DEdep.out", std::ios::app);
550  // Index of targeted compartments
551  WriteDataInSoma //<< i+1 << '\t' << " "
552  // position of compartments
553  << fDetector->GetPosSomacomp(i).x()<< '\t' << " "
554  << fDetector->GetPosSomacomp(i).y()<< '\t' << " "
555  << fDetector->GetPosSomacomp(i).z()<< '\t' << " "
556  // Edep in compartments
557  << fSoma3DEdep[i]/keV << '\t' << " "
558  // Dose in compartments
559  //<< (fSoma3DEdep[i]/joule)/(fDetector->GetMassSomacomp(i)/kg)
560  // Dose in whole structure of Soma
561  //<< (fSoma3DEdep[i]/joule)/(fDetector->GetMassTotSo1()/kg)
562  //<< '\t' << " "
563  << G4endl;
564  }
565  }
566  // compute mean Energy and Dose deposited in compartments;
567  // +- RMS : Root Mean Square
568  G4double rmsEdepS =0.;
569  G4double rmsDoseS =0.;
570  if (somaCompHit >0)
571  {
572  somaCompEdep /= somaCompHit; somaCompEdep2 /= somaCompHit;
573  rmsEdepS = somaCompEdep2 - somaCompEdep*somaCompEdep;
574  if (rmsEdepS>0.) rmsEdepS = std::sqrt(rmsEdepS);
575  else rmsEdepS = 0.;
576  somaCompDose /= somaCompHit; somaCompDose2 /= somaCompHit;
577  rmsDoseS = somaCompDose2 - somaCompDose*somaCompDose;
578  if (rmsDoseS>0.) rmsDoseS = std::sqrt(rmsDoseS);
579  else rmsDoseS = 0.;
580  }
581 
582  G4int DendCompHit = 0;
583  G4double DendCompEdep = 0.;
584  G4double DendCompDose = 0.;
585  G4double DendCompEdep2 = 0.;
586  G4double DendCompDose2 = 0.;
587  remove ("Dend3DEdep.out");
588  for (G4int i=0; i<fDetector->GetnbDendritecomp(); i++)
589  {
590  if (fDend3DEdep[i] > 0.0)
591  {
592  DendCompHit ++;
593  DendCompEdep += fDend3DEdep[i] ;
594  //DendCompDose += fDend3DEdep[i]/fDetector->GetMassDendcomp(i) ;
595  DendCompEdep2 += fDend3DEdep[i]*fDend3DEdep[i] ;
596  //DendCompDose2 += (fDend3DEdep[i]/fDetector->GetMassDendcomp(i))*
597  // (fDend3DEdep[i]/fDetector->GetMassDendcomp(i));
598  //Fill ntuple #2
599  /*analys->FillNtupleDColumn(2,0,i+1);
600  analys->FillNtupleDColumn(2,1,fDetector->GetDistDendsoma(i));
601  analys->FillNtupleDColumn(2,2,fDend3DEdep[i]/keV);
602  analys->FillNtupleDColumn(2,3,(fDend3DEdep[i]/joule)/
603  (fDetector->GetMassDendcomp(i)/kg));
604  analys->AddNtupleRow(2);
605  */
606 
607  std::ofstream WriteDataInDend("Dend3DEdep.out", std::ios::app);
608  WriteDataInDend //<< i+1 << '\t' << " "
609  // position of compartments
610  << fDetector->GetPosDendcomp(i).x()<< '\t' << " "
611  << fDetector->GetPosDendcomp(i).y()<< '\t' << " "
612  << fDetector->GetPosDendcomp(i).z()<< '\t' << " "
613  << fDetector->GetDistADendSoma(i)<< '\t' << " "
614  << fDetector->GetDistBDendSoma(i)<< '\t' << " "
615  << fDend3DEdep[i]/keV << '\t' << " "
616  //<< (fDend3DEdep[i]/joule)/(fDetector->GetMassDendcomp(i)/kg)
617  << G4endl;
618  }
619  }
620  // +- RMS : Root Mean Square
621  G4double rmsEdepD =0.;
622  G4double rmsDoseD =0.;
623  if (DendCompHit >0)
624  {
625  DendCompEdep /= DendCompHit; DendCompEdep2 /= DendCompHit;
626  rmsEdepD = DendCompEdep2 - DendCompEdep*DendCompEdep;
627  if (rmsEdepD>0.) rmsEdepD = std::sqrt(rmsEdepD);
628  else rmsEdepD = 0.;
629  DendCompDose /= DendCompHit; DendCompDose2 /= DendCompHit;
630  rmsDoseD = DendCompDose2 - DendCompDose*DendCompDose;
631  if (rmsDoseD>0.) rmsDoseD = std::sqrt(rmsDoseD);
632  else rmsDoseD = 0.;
633  }
634 
635  G4int AxonCompHit = 0;
636  G4double AxonCompEdep = 0.;
637  G4double AxonCompDose = 0.;
638  G4double AxonCompEdep2 = 0.;
639  G4double AxonCompDose2 = 0.;
640  remove ("Axon3DEdep.out");
641  for (G4int i=0; i<fDetector->GetnbAxoncomp(); i++)
642  {
643  if (fAxon3DEdep[i] > 0.0)
644  {
645  AxonCompHit ++;
646  AxonCompEdep += fAxon3DEdep[i] ;
647  //AxonCompDose += fAxon3DEdep[i]/fDetector->GetMassAxoncomp(i) ;
648  AxonCompEdep2 += fAxon3DEdep[i]*fAxon3DEdep[i] ;
649  //AxonCompDose2 += (fAxon3DEdep[i]/fDetector->GetMassAxoncomp(i))*
650  // (fAxon3DEdep[i]/fDetector->GetMassAxoncomp(i));
651  //Fill ntuple #3
652  /*analys->FillNtupleDColumn(3,0,i+1);
653  analys->FillNtupleDColumn(3,1,fDetector->GetDistAxonsoma(i));
654  analys->FillNtupleDColumn(3,2,fAxon3DEdep[i]/keV);
655  analys->FillNtupleDColumn(3,3,(fAxon3DEdep[i]/joule)/
656  (fDetector->GetMassAxoncomp(i)/kg));
657  analys->AddNtupleRow(3);
658  */
659 
660  std::ofstream WriteDataInAxon("Axon3DEdep.out", std::ios::app);
661  WriteDataInAxon //<< i+1 << '\t' << " "
662  // position of compartments
663  << fDetector->GetPosAxoncomp(i).x()<< '\t' << " "
664  << fDetector->GetPosAxoncomp(i).y()<< '\t' << " "
665  << fDetector->GetPosAxoncomp(i).z()<< '\t' << " "
666  << fDetector->GetDistAxonsoma(i) << '\t' << " "
667  << fAxon3DEdep[i]/keV << '\t' << " "
668  //<< (fAxon3DEdep[i]/joule)/(fDetector->GetMassAxoncomp(i)/kg)
669  << G4endl;
670  }
671  }
672  // +- RMS : Root Mean Square
673  G4double rmsEdepA =0.;
674  G4double rmsDoseA =0.;
675  if (AxonCompHit >0)
676  {
677  AxonCompEdep /= AxonCompHit; AxonCompEdep2 /= AxonCompHit;
678  rmsEdepA = AxonCompEdep2 - AxonCompEdep*AxonCompEdep;
679  if (rmsEdepA>0.) rmsEdepA = std::sqrt(rmsEdepA);
680  else rmsEdepA = 0.;
681  AxonCompDose /= AxonCompHit; AxonCompDose2 /= AxonCompHit;
682  rmsDoseA = AxonCompDose2 - AxonCompDose*AxonCompDose;
683  if (rmsDoseA>0.) rmsDoseA = std::sqrt(rmsDoseA);
684  else rmsDoseA = 0.;
685  }
686 
687  G4cout << "\n Number of compartments traversed by particle tracks :"
688  << G4endl;
689  G4cout << " " << std::setw(13) << "Soma: " << std::setw(7) << somaCompHit
690  << " of total: "<< fDetector->GetnbSomacomp() << "\n "
691  << " " << std::setw(13) << "Dendrites: " << std::setw(7) << DendCompHit
692  << " of total: "<< fDetector->GetnbDendritecomp() << "\n "
693  << " " << std::setw(13) << "Axon: " << std::setw(7) << AxonCompHit
694  << " of total: "<< fDetector->GetnbAxoncomp() << "\n "
695  << G4endl;
696 /*
697  G4cout << "\n Mean energy (keV) and dose (Gy) deposition in "
698  << "hitting compartments :" << G4endl;
699  G4cout
700  << " " << std::setw(13) << "Soma: " << std::setw(7) << somaCompEdep/keV
701  << " +- "<< rmsEdepS/keV << " and "
702  << somaCompDose/(joule/kg)<< " +- "<< rmsDoseS/(joule/kg)<< "\n "
703  << " " << std::setw(13) << "Dendrites: " << std::setw(7)
704  << DendCompEdep/keV
705  << " +- "<< rmsEdepD/keV << " and "
706  << DendCompDose/(joule/kg)<< " +- "<< rmsDoseD/(joule/kg)<< "\n "
707  << " " << std::setw(13) << "Axon: " << std::setw(7) << AxonCompEdep/keV
708  << " +- "<< rmsEdepA/keV << " and "
709  << AxonCompDose/(joule/kg)<< " +- "<< rmsDoseA/(joule/kg)<< "\n "
710  << G4endl; */
711  /*
712  // compute mean Energy and Dose deposited per event
713  //
714  fEdepNeuron /= TotNbofEvents; fEdepNeuron_err /= TotNbofEvents;
715  G4double rmsEdep = fEdepNeuron_err - fEdepNeuron*fEdepNeuron;
716  if (rmsEdep>0.) rmsEdep = std::sqrt(rmsEdep);
717  else rmsEdep = 0.;
718  G4double fDoseNeuron = ((fEdepNeuron/joule)/ (fDetector->
719  GetTotMassNeuron()/kg));
720  G4double fDoseNeuron_err = ((fEdepNeuron_err/joule)/ (fDetector->
721  GetTotMassNeuron()/kg));
722  G4double rmsDose = fDoseNeuron_err - fDoseNeuron*fDoseNeuron;
723  if (rmsDose>0.) rmsDose = std::sqrt(rmsDose);
724  else rmsDose = 0.;
725  G4cout << "\n Mean energy (keV) deposition per event:"
726  << G4endl;
727  G4cout
728  << " " << std::setw(13) << "Neuron: " << std::setw(7) << fEdepNeuron/keV
729  << " +- "<< rmsEdep/keV << " and "
730  << std::setw(wid) << fDoseNeuron
731  << " +- "<< rmsDose << "\n "
732  << G4endl; */
733  G4cout<< "\n Dendritic (or Axon) compartmental energy deposits \n"
734  << " at the distance from Soma have been written into *.out files:"
735  << "\n Dend3DEdep.out, Axon3DEdep.out, Soma3DEdep.out"
736  << "\n Outputs of energy deposition per event written in data file:"
737  << "\n OutputPerEvent.out"
738  << "\n " << G4endl;
739 
740  //remove all contents in fProcCounter, fCount
741  fProcCounter.clear();
742  fParticleDataMap2.clear();
743 
744  //restore default format
745  G4cout.precision(dfprec);
746 }
747 
748 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......