Geant4 Cross Reference

Cross-Referencing   Geant4
Geant4/examples/extended/electromagnetic/TestEm11/src/Run.cc

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

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