| Geant4 Cross Reference |
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