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