Geant4 Cross Reference

Cross-Referencing   Geant4
Geant4/examples/extended/medical/dna/neuron/src/Run.cc

Version: [ ReleaseNotes ] [ 1.0 ] [ 1.1 ] [ 2.0 ] [ 3.0 ] [ 3.1 ] [ 3.2 ] [ 4.0 ] [ 4.0.p1 ] [ 4.0.p2 ] [ 4.1 ] [ 4.1.p1 ] [ 5.0 ] [ 5.0.p1 ] [ 5.1 ] [ 5.1.p1 ] [ 5.2 ] [ 5.2.p1 ] [ 5.2.p2 ] [ 6.0 ] [ 6.0.p1 ] [ 6.1 ] [ 6.2 ] [ 6.2.p1 ] [ 6.2.p2 ] [ 7.0 ] [ 7.0.p1 ] [ 7.1 ] [ 7.1.p1 ] [ 8.0 ] [ 8.0.p1 ] [ 8.1 ] [ 8.1.p1 ] [ 8.1.p2 ] [ 8.2 ] [ 8.2.p1 ] [ 8.3 ] [ 8.3.p1 ] [ 8.3.p2 ] [ 9.0 ] [ 9.0.p1 ] [ 9.0.p2 ] [ 9.1 ] [ 9.1.p1 ] [ 9.1.p2 ] [ 9.1.p3 ] [ 9.2 ] [ 9.2.p1 ] [ 9.2.p2 ] [ 9.2.p3 ] [ 9.2.p4 ] [ 9.3 ] [ 9.3.p1 ] [ 9.3.p2 ] [ 9.4 ] [ 9.4.p1 ] [ 9.4.p2 ] [ 9.4.p3 ] [ 9.4.p4 ] [ 9.5 ] [ 9.5.p1 ] [ 9.5.p2 ] [ 9.6 ] [ 9.6.p1 ] [ 9.6.p2 ] [ 9.6.p3 ] [ 9.6.p4 ] [ 10.0 ] [ 10.0.p1 ] [ 10.0.p2 ] [ 10.0.p3 ] [ 10.0.p4 ] [ 10.1 ] [ 10.1.p1 ] [ 10.1.p2 ] [ 10.1.p3 ] [ 10.2 ] [ 10.2.p1 ] [ 10.2.p2 ] [ 10.2.p3 ] [ 10.3 ] [ 10.3.p1 ] [ 10.3.p2 ] [ 10.3.p3 ] [ 10.4 ] [ 10.4.p1 ] [ 10.4.p2 ] [ 10.4.p3 ] [ 10.5 ] [ 10.5.p1 ] [ 10.6 ] [ 10.6.p1 ] [ 10.6.p2 ] [ 10.6.p3 ] [ 10.7 ] [ 10.7.p1 ] [ 10.7.p2 ] [ 10.7.p3 ] [ 10.7.p4 ] [ 11.0 ] [ 11.0.p1 ] [ 11.0.p2 ] [ 11.0.p3, ] [ 11.0.p4 ] [ 11.1 ] [ 11.1.1 ] [ 11.1.2 ] [ 11.1.3 ] [ 11.2 ] [ 11.2.1 ] [ 11.2.2 ] [ 11.3.0 ]

Diff markup

Differences between /examples/extended/medical/dna/neuron/src/Run.cc (Version 11.3.0) and /examples/extended/medical/dna/neuron/src/Run.cc (Version 10.7.p4)


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