ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
Run.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file Run.cc
1 //
2 // ********************************************************************
3 // * License and Disclaimer *
4 // * *
5 // * The Geant4 software is copyright of the Copyright Holders of *
6 // * the Geant4 Collaboration. It is provided under the terms and *
7 // * conditions of the Geant4 Software License, included in the file *
8 // * LICENSE and available at http://cern.ch/geant4/license . These *
9 // * include a list of copyright holders. *
10 // * *
11 // * Neither the authors of this software system, nor their employing *
12 // * institutes,nor the agencies providing financial support for this *
13 // * work make any representation or warranty, express or implied, *
14 // * regarding this software system or assume any liability for its *
15 // * use. Please see the license in the file LICENSE and URL above *
16 // * for the full disclaimer and the limitation of liability. *
17 // * *
18 // * This code implementation is the result of the scientific and *
19 // * technical work of the GEANT4 collaboration. *
20 // * By using, copying, modifying or distributing the software (or *
21 // * any work based on the software) you agree to acknowledge its *
22 // * use in resulting scientific publications, and indicate your *
23 // * acceptance of all terms of the Geant4 Software license. *
24 // ********************************************************************
25 //
26 // This example is provided by the Geant4-DNA collaboration
27 // Any report or published results obtained using the Geant4-DNA software
28 // shall cite the following Geant4-DNA collaboration publications:
29 // Med. Phys. 37 (2010) 4692-4708
30 // Phys. Med. 31 (2015) 861-874
31 // The Geant4-DNA web site is available at http://geant4-dna.org
32 //
35 
36 #include "Run.hh"
37 #include "DetectorConstruction.hh"
38 #include "HistoManager.hh"
39 #include "MyFile.hh"
40 
41 #ifdef MYFILE
42  #include "MyPrimaryGeneratorActionFromFile.hh"
43 #else
44  #include "PrimaryGeneratorAction.hh"
45 #endif
46 
47 #include "G4Material.hh"
48 #include "G4SystemOfUnits.hh"
49 #include "G4UnitsTable.hh"
50 
51 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
52 
53 Run::Run(const DetectorConstruction* detector)
54 : G4Run(),
55  fDetector(detector),
56  fParticle(0), fEkin(0.),
57  fCytoEdeposit(0.), fCytoEdeposit2(0.),
58  fNuclEdeposit(0.), fNuclEdeposit2(0.),
59  fTrackLen(0.), fTrackLen2(0.),
60  fProjRange(0.), fProjRange2(0.),
61  fNbOfSteps(0), fNbOfSteps2(0),
62  fStepSize(0.), fStepSize2(0.)
63 {}
64 
65 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
66 
67 Run::~Run()
68 {}
69 
70 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
71 
72 void Run::SetPrimary (G4ParticleDefinition* particle, G4double energy)
73 {
74  fParticle = particle;
75  fEkin = energy;
76 }
77 
78 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
79 
80 void Run::AddCytoEdep (G4double e)
81 {
82  fCytoEdeposit += e;
83  fCytoEdeposit2 += e*e;
84 }
85 
86 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
87 
88 void Run::AddNuclEdep (G4double e)
89 {
90  fNuclEdeposit += e;
91  fNuclEdeposit2 += e*e;
92 }
93 
94 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
95 
96 void Run::AddTrackLength (G4double t)
97 {
98  fTrackLen += t;
99  fTrackLen2 += t*t;
100 }
101 
102 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
103 
104 void Run::AddProjRange (G4double x)
105 {
106  fProjRange += x;
107  fProjRange2 += x*x;
108 }
109 
110 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
111 
112 void Run::AddStepSize (G4int nb, G4double st)
113 {
114  fNbOfSteps += nb;
115  fNbOfSteps2 += nb*nb;
116  fStepSize += st ;
117  fStepSize2 += st*st;
118 }
119 
120 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
121 
122 void Run::Merge(const G4Run* run)
123 {
124  const Run* localRun = static_cast<const Run*>(run);
125 
126  // pass information about primary particle
127  fParticle = localRun->fParticle;
128  fEkin = localRun->fEkin;
129 
130  // accumulate sums
131 
132  fCytoEdeposit += localRun->fCytoEdeposit;
133  fCytoEdeposit2 += localRun->fCytoEdeposit2;
134  fNuclEdeposit += localRun->fNuclEdeposit;
135  fNuclEdeposit2 += localRun->fNuclEdeposit2;
136 
137  fTrackLen += localRun->fTrackLen;
138  fTrackLen2 += localRun->fTrackLen2;
139  fProjRange += localRun->fProjRange;
140  fProjRange2 += localRun->fProjRange2;
141  fNbOfSteps += localRun->fNbOfSteps ;
142  fNbOfSteps2 += localRun->fNbOfSteps2;
143  fStepSize += localRun->fStepSize;
144  fStepSize2 += localRun->fStepSize2;
145 
146  G4Run::Merge(run);
147 }
148 
149 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
150 
151 void Run::EndOfRun()
152 {
153  std::ios::fmtflags mode = G4cout.flags();
154  G4cout.setf(std::ios::fixed,std::ios::floatfield);
155  G4int prec = G4cout.precision(2);
156 
157  //run conditions
158  //
159  G4String partName = fParticle->GetParticleName();
160 
161  G4cout << "\n ======================== run summary =====================\n";
162  G4cout
163  << "\n The run is " << numberOfEvent << " "<< partName << " of "
164  << G4BestUnit(fEkin,"Energy") << G4endl;
165 
166  if (numberOfEvent == 0) {
167  G4cout.setf(mode,std::ios::floatfield);
168  G4cout.precision(prec);
169  return;
170  }
171 
172  //
173 
174  fCytoEdeposit /= numberOfEvent; fCytoEdeposit2 /= numberOfEvent;
175  G4double rmsCyto = fCytoEdeposit2 - fCytoEdeposit*fCytoEdeposit;
176  if (rmsCyto>0.) rmsCyto = std::sqrt(rmsCyto); else rmsCyto = 0.;
177 
178  G4cout.precision(3);
179  G4cout
180  << "\n Total Energy deposited in cytoplasm = " << G4BestUnit(fCytoEdeposit,"Energy")
181  << " +- " << G4BestUnit( rmsCyto,"Energy")
182  << G4endl;
183 
184  G4double sValueCyto=fCytoEdeposit/fDetector->GetCytoMass();
185  G4double rmsSValueCyto=rmsCyto/fDetector->GetCytoMass();
186 
187  G4cout.precision(3);
188  G4cout
189  << "\n S value for cytoplasm (C<-C) = " << sValueCyto/gray << " Gy/Bq.s "
190  << " +- " << rmsSValueCyto/gray
191  << " Gy/Bq.s "
192  << G4endl;
193 
194  //
195 
196  fNuclEdeposit /= numberOfEvent; fNuclEdeposit2 /= numberOfEvent;
197  G4double rmsNucl = fNuclEdeposit2 - fNuclEdeposit*fNuclEdeposit;
198  if (rmsNucl>0.) rmsNucl = std::sqrt(rmsNucl); else rmsNucl = 0.;
199 
200  G4cout.precision(3);
201  G4cout
202  << "\n Total Energy deposited in nucleus = " << G4BestUnit(fNuclEdeposit,"Energy")
203  << " +- " << G4BestUnit( rmsNucl,"Energy")
204  << G4endl;
205 
206  G4double sValueNucl=fNuclEdeposit/fDetector->GetNuclMass();
207  G4double rmsSValueNucl=rmsNucl/fDetector->GetNuclMass();
208 
209  G4cout.precision(3);
210  G4cout
211  << "\n S value for nucleus (N<-C) = " << sValueNucl/gray << " Gy/Bq.s "
212  << " +- " << rmsSValueNucl/gray
213  << " Gy/Bq.s "
214  << G4endl;
215 
216  //compute track length of primary track
217  //
218  fTrackLen /= numberOfEvent; fTrackLen2 /= numberOfEvent;
219  G4double rms = fTrackLen2 - fTrackLen*fTrackLen;
220  if (rms>0.) rms = std::sqrt(rms); else rms = 0.;
221 
222  G4cout.precision(3);
223  G4cout
224  << "\n Track length of primary track = " << G4BestUnit(fTrackLen,"Length")
225  << " +- " << G4BestUnit( rms,"Length");
226 
227  //compute projected range of primary track
228  //
229  fProjRange /= numberOfEvent; fProjRange2 /= numberOfEvent;
230  rms = fProjRange2 - fProjRange*fProjRange;
231  if (rms>0.) rms = std::sqrt(rms); else rms = 0.;
232 
233  G4cout
234  << "\n Projected range = " << G4BestUnit(fProjRange,"Length")
235  << " +- " << G4BestUnit( rms,"Length")
236  << G4endl;
237 
238  //nb of steps and step size of primary track
239  //
240  G4double dNofEvents = double(numberOfEvent);
241  G4double fNbSteps = fNbOfSteps/dNofEvents,
242  fNbSteps2 = fNbOfSteps2/dNofEvents;
243  rms = fNbSteps2 - fNbSteps*fNbSteps;
244  if (rms>0.) rms = std::sqrt(rms); else rms = 0.;
245 
246  G4cout.precision(2);
247  G4cout << "\n Nb of steps of primary track = " << fNbSteps << " +- " << rms
248  << G4endl;
249 
250  fStepSize /= numberOfEvent; fStepSize2 /= numberOfEvent;
251  rms = fStepSize2 - fStepSize*fStepSize;
252  if (rms>0.) rms = std::sqrt(rms); else rms = 0.;
253 
254  G4cout.precision(3);
255  G4cout
256  << "\n Step size = " << G4BestUnit(fStepSize,"Length")
257  << " +- " << G4BestUnit( rms,"Length")
258  << G4endl;
259 
260  // normalize histograms of longitudinal energy profile
261  //
262  G4AnalysisManager* analysisManager = G4AnalysisManager::Instance();
263  G4int ih = 1;
264  G4double binWidth = analysisManager->GetH1Width(ih);
265  G4double fac = (1./(numberOfEvent*binWidth))*(mm/MeV);
266  analysisManager->ScaleH1(ih,fac);
267 
268  // reset default formats
269  G4cout.setf(mode,std::ios::floatfield);
270  G4cout.precision(prec);
271 
272  //output file
273  FILE *myFile;
274  myFile = fopen ("s.txt","a");
275  fprintf (myFile, "%e %e %e %e %e %e %e \n",
276  fDetector->GetNuclRadius()/nm,
277  fDetector->GetCytoThickness()/nm,
278  fEkin/eV,
279  sValueCyto/gray,
280  rmsSValueCyto/gray,
281  sValueNucl/gray,
282  rmsSValueNucl/gray
283  );
284  fclose (myFile);
285 
286 }