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 //
28 //
29 //
30 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
31 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
32 
33 #include "Run.hh"
34 #include "DetectorConstruction.hh"
35 
36 #include "EventAction.hh"
37 #include "HistoManager.hh"
38 #include "PrimaryGeneratorAction.hh"
39 
40 #include "G4Material.hh"
41 #include "G4Event.hh"
42 #include "G4SystemOfUnits.hh"
43 #include "G4UnitsTable.hh"
44 #include <iomanip>
45 
46 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
47 
48 Run::Run(DetectorConstruction* detector)
49 : G4Run(),
50  fDetector(detector),
51  fParticle(0), fEkin(0.),
52  fTrackLen(0.), fTrackLen2(0.),
53  fProjRange(0.), fProjRange2(0.),
54  fNbOfSteps(0), fNbOfSteps2(0),
55  fStepSize(0.), fStepSize2(0.)
56 {
57  for (G4int i=0; i<3; ++i) { fStatus[i] = 0; fTotEdep[i] = 0.; }
58  fTotEdep[1] = joule;
59  for (G4int i=0; i<kMaxAbsor; ++i) {
60  fEdeposit[i] = 0.; fEmin[i] = joule; fEmax[i] = 0.;
61  fCsdaRange[i] = 0.; fXfrontNorm[i] = 0.;
62  }
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::AddEdep (G4int i, G4double e)
81 {
82  if (e > 0.) {
83  fEdeposit[i] += e;
84  if (e < fEmin[i]) fEmin[i] = e;
85  if (e > fEmax[i]) fEmax[i] = e;
86  }
87 }
88 
89 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
90 
91 void Run::AddTotEdep (G4double e)
92 {
93  if (e > 0.) {
94  fTotEdep[0] += e;
95  if (e < fTotEdep[1]) fTotEdep[1] = e;
96  if (e > fTotEdep[2]) fTotEdep[2] = e;
97  }
98 }
99 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
100 
101 void Run::AddTrackLength (G4double t)
102 {
103  fTrackLen += t;
104  fTrackLen2 += t*t;
105 }
106 
107 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
108 
109 void Run::AddProjRange (G4double x)
110 {
111  fProjRange += x;
112  fProjRange2 += x*x;
113 }
114 
115 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
116 
117 void Run::AddStepSize (G4int nb, G4double st)
118 {
119  fNbOfSteps += nb;
120  fNbOfSteps2 += nb*nb;
121  fStepSize += st ;
122  fStepSize2 += st*st;
123 }
124 
125 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
126 
127 void Run::AddTrackStatus (G4int i)
128 {
129  fStatus[i]++ ;
130 }
131 
132 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
133 
134 void Run::SetCsdaRange (G4int i, G4double value)
135 {
136  fCsdaRange[i] = value;
137 }
138 
139 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
140 
141 void Run::SetXfrontNorm (G4int i, G4double value)
142 {
143  fXfrontNorm[i] = value;
144 }
145 
146 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
147 
148 G4double Run::GetCsdaRange (G4int i)
149 {
150  return fCsdaRange[i];
151 }
152 
153 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
154 
155 G4double Run::GetXfrontNorm (G4int i)
156 {
157  return fXfrontNorm[i];
158 }
159 
160 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
161 
162 void Run::Merge(const G4Run* run)
163 {
164  const Run* localRun = static_cast<const Run*>(run);
165 
166  // pass information about primary particle
167  fParticle = localRun->fParticle;
168  fEkin = localRun->fEkin;
169 
170  // accumulate sums
171  fTrackLen += localRun->fTrackLen;
172  fTrackLen2 += localRun->fTrackLen2;
173  fProjRange += localRun->fProjRange;
174  fProjRange2 += localRun->fProjRange2;
175  fNbOfSteps += localRun->fNbOfSteps ;
176  fNbOfSteps2 += localRun->fNbOfSteps2;
177  fStepSize += localRun->fStepSize;
178  fStepSize2 += localRun->fStepSize2;
179 
180  G4int nbOfAbsor = fDetector->GetNbOfAbsor();
181  for (G4int i=1; i<=nbOfAbsor; ++i) {
182  fEdeposit[i] += localRun->fEdeposit[i];
183  fCsdaRange[i] = localRun->fCsdaRange[i];
184  fXfrontNorm[i] = localRun->fXfrontNorm[i];
185  // min, max
186  G4double min,max;
187  min = localRun->fEmin[i]; max = localRun->fEmax[i];
188  if (fEmin[i] > min) fEmin[i] = min;
189  if (fEmax[i] < max) fEmax[i] = max;
190  }
191 
192  for (G4int i=0; i<3; ++i) fStatus[i] += localRun->fStatus[i];
193 
194  // total Edep
195  fTotEdep[0] += localRun->fTotEdep[0];
196  G4double min,max;
197  min = localRun->fTotEdep[1]; max = localRun->fTotEdep[2];
198 if (fTotEdep[1] > min) fTotEdep[1] = min;
199 if (fTotEdep[2] < max) fTotEdep[2] = max;
200 
201  G4Run::Merge(run);
202 }
203 
204 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
205 
206 void Run::EndOfRun()
207 {
208  std::ios::fmtflags mode = G4cout.flags();
209  G4cout.setf(std::ios::fixed,std::ios::floatfield);
210  G4int prec = G4cout.precision(2);
211 
212  //run conditions
213  //
214  G4String partName = fParticle->GetParticleName();
215  G4int nbOfAbsor = fDetector->GetNbOfAbsor();
216 
217  G4cout << "\n ======================== run summary =====================\n";
218  G4cout
219  << "\n The run is " << numberOfEvent << " "<< partName << " of "
220  << G4BestUnit(fEkin,"Energy")
221  << " through " << nbOfAbsor << " absorbers: \n";
222  for (G4int i=1; i<= nbOfAbsor; i++) {
223  G4Material* material = fDetector->GetAbsorMaterial(i);
224  G4double thickness = fDetector->GetAbsorThickness(i);
225  G4double density = material->GetDensity();
226  G4cout << std::setw(5) << i
227  << std::setw(10) << G4BestUnit(thickness,"Length") << " of "
228  << material->GetName() << " (density: "
229  << G4BestUnit(density,"Volumic Mass") << ")" << G4endl;
230  }
231 
232  if (numberOfEvent == 0) {
233  G4cout.setf(mode,std::ios::floatfield);
234  G4cout.precision(prec);
235  return;
236  }
237 
238  G4cout.precision(3);
239  G4double rms (0);
240 
241  for (G4int i=1; i<= nbOfAbsor; i++) {
242  fEdeposit[i] /= numberOfEvent;
243 
244  G4cout
245  << "\n Edep in absorber " << i << " = "
246  << G4BestUnit(fEdeposit[i],"Energy")
247  << "\t(" << G4BestUnit(fEmin[i], "Energy")
248  << "-->" << G4BestUnit(fEmax[i], "Energy")
249  << ")";
250  }
251  G4cout << G4endl;
252 
253  if (nbOfAbsor > 1) {
254  fTotEdep[0] /= numberOfEvent;
255  G4cout
256  << "\n Edep in all absorbers = " << G4BestUnit(fTotEdep[0],"Energy")
257  << "\t(" << G4BestUnit(fTotEdep[1], "Energy")
258  << "-->" << G4BestUnit(fTotEdep[2], "Energy")
259  << ")" << G4endl;
260  }
261 
262  //compute track length of primary track
263  //
264  fTrackLen /= numberOfEvent; fTrackLen2 /= numberOfEvent;
265  rms = fTrackLen2 - fTrackLen*fTrackLen;
266  if (rms>0.) rms = std::sqrt(rms); else rms = 0.;
267 
268  G4cout.precision(3);
269  G4cout
270  << "\n Track length of primary track = " << G4BestUnit(fTrackLen,"Length")
271  << " +- " << G4BestUnit( rms,"Length");
272 
273  //compare with csda range
274  //
275  G4int NbOfAbsor = fDetector->GetNbOfAbsor();
276  if (NbOfAbsor == 1) {
277  G4cout
278  << "\n Range from EmCalculator = " << G4BestUnit(fCsdaRange[1],"Length")
279  << " (from full dE/dx)" << G4endl;
280  }
281 
282  //compute projected range of primary track
283  //
284  fProjRange /= numberOfEvent; fProjRange2 /= numberOfEvent;
285  rms = fProjRange2 - fProjRange*fProjRange;
286  if (rms>0.) rms = std::sqrt(rms); else rms = 0.;
287 
288  G4cout
289  << "\n Projected range = " << G4BestUnit(fProjRange,"Length")
290  << " +- " << G4BestUnit( rms,"Length")
291  << G4endl;
292 
293  //nb of steps and step size of primary track
294  //
295  G4double dNofEvents = double(numberOfEvent);
296  G4double fNbSteps = fNbOfSteps/dNofEvents,
297  fNbSteps2 = fNbOfSteps2/dNofEvents;
298  rms = fNbSteps2 - fNbSteps*fNbSteps;
299  if (rms>0.) rms = std::sqrt(rms); else rms = 0.;
300 
301  G4cout.precision(2);
302  G4cout << "\n Nb of steps of primary track = " << fNbSteps << " +- " << rms;
303 
304  fStepSize /= numberOfEvent; fStepSize2 /= numberOfEvent;
305  rms = fStepSize2 - fStepSize*fStepSize;
306  if (rms>0.) rms = std::sqrt(rms); else rms = 0.;
307 
308  G4cout.precision(3);
309  G4cout
310  << "\t Step size= " << G4BestUnit(fStepSize,"Length")
311  << " +- " << G4BestUnit( rms,"Length")
312  << G4endl;
313 
314  //transmission coefficients
315  //
316  G4double absorbed = 100.*fStatus[0]/dNofEvents;
317  G4double transmit = 100.*fStatus[1]/dNofEvents;
318  G4double reflected = 100.*fStatus[2]/dNofEvents;
319 
320  G4cout.precision(2);
321  G4cout
322  << "\n absorbed = " << absorbed << " %"
323  << " transmit = " << transmit << " %"
324  << " reflected = " << reflected << " %" << G4endl;
325 
326  // normalize histograms of longitudinal energy profile
327  //
328  G4AnalysisManager* analysisManager = G4AnalysisManager::Instance();
329  G4int ih = 1;
330  G4double binWidth = analysisManager->GetH1Width(ih)
331  *analysisManager->GetH1Unit(ih);
332  G4double fac = (1./(numberOfEvent*binWidth))*(mm/MeV);
333  analysisManager->ScaleH1(ih,fac);
334 
335  ih = 8;
336  binWidth = analysisManager->GetH1Width(ih);
337  fac = (1./(numberOfEvent*binWidth))*(g/(MeV*cm2));
338  analysisManager->ScaleH1(ih,fac);
339 
340  // reset default formats
341  G4cout.setf(mode,std::ios::floatfield);
342  G4cout.precision(prec);
343 }
344 
345 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......