Geant4  10.03.p03
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
Run.cc
Go to the documentation of this file.
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 // $Id: Run.cc 71376 2013-06-14 07:44:50Z maire $
30 //
31 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
32 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
33 
34 #include "Run.hh"
35 #include "DetectorConstruction.hh"
36 #include "PrimaryGeneratorAction.hh"
37 #include "HistoManager.hh"
38 
39 #include "G4EmCalculator.hh"
40 #include "G4SystemOfUnits.hh"
41 #include "G4UnitsTable.hh"
42 
43 #include <iomanip>
44 
45 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
46 
48 : G4Run(),
49  fDetector(det),
50  fParticle(0), fEkin(0.)
51 {
52  fEnergyDeposit = fEnergyDeposit2 = 0.;
53  fTrakLenCharged = fTrakLenCharged2 = 0.;
54  fTrakLenNeutral = fTrakLenNeutral2 = 0.;
55  fNbStepsCharged = fNbStepsCharged2 = 0.;
56  fNbStepsNeutral = fNbStepsNeutral2 = 0.;
57  fMscProjecTheta = fMscProjecTheta2 = 0.;
58  fMscThetaCentral = 0.;
59 
60  fNbGamma = fNbElect = fNbPosit = 0;
61 
62  fTransmit[0] = fTransmit[1] = fReflect[0] = fReflect[1] = 0;
63 
64  fMscEntryCentral = 0;
65 
66  fEnergyLeak[0] = fEnergyLeak[1] = fEnergyLeak2[0] = fEnergyLeak2[1] = 0.;
67  }
68 
69 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
70 
71 Run::~Run()
72 { }
73 
74 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
75 
77 {
78  fParticle = particle;
79  fEkin = energy;
80 
81  //compute theta0
82  fMscThetaCentral = 3*ComputeMscHighland();
83 }
84 
85 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
86 
87 void Run::Merge(const G4Run* run)
88 {
89  const Run* localRun = static_cast<const Run*>(run);
90 
91  // pass information about primary particle
92  fParticle = localRun->fParticle;
93  fEkin = localRun->fEkin;
94 
95  fMscThetaCentral = localRun->fMscThetaCentral;
96 
97  // accumulate sums
98  //
99  fEnergyDeposit += localRun->fEnergyDeposit;
100  fEnergyDeposit2 += localRun->fEnergyDeposit2;
101  fTrakLenCharged += localRun->fTrakLenCharged;
102  fTrakLenCharged2 += localRun->fTrakLenCharged2;
103  fTrakLenNeutral += localRun->fTrakLenNeutral;
104  fTrakLenNeutral2 += localRun->fTrakLenNeutral2;
105  fNbStepsCharged += localRun->fNbStepsCharged;
106  fNbStepsCharged2 += localRun->fNbStepsCharged2;
107  fNbStepsNeutral += localRun->fNbStepsNeutral;
108  fNbStepsNeutral2 += localRun->fNbStepsNeutral2;
109  fMscProjecTheta += localRun->fMscProjecTheta;
110  fMscProjecTheta2 += localRun->fMscProjecTheta2;
111 
112 
113  fNbGamma += localRun->fNbGamma;
114  fNbElect += localRun->fNbElect;
115  fNbPosit += localRun->fNbPosit;
116 
117  fTransmit[0] += localRun->fTransmit[0];
118  fTransmit[1] += localRun->fTransmit[1];
119  fReflect[0] += localRun->fReflect[0];
120  fReflect[1] += localRun->fReflect[1];
121 
122  fMscEntryCentral += localRun->fMscEntryCentral;
123 
124  fEnergyLeak[0] += localRun->fEnergyLeak[0];
125  fEnergyLeak[1] += localRun->fEnergyLeak[1];
126  fEnergyLeak2[0] += localRun->fEnergyLeak2[0];
127  fEnergyLeak2[1] += localRun->fEnergyLeak2[1];
128 
129  G4Run::Merge(run);
130 }
131 
132 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
133 
134 void Run::EndOfRun()
135 {
136  // compute mean and rms
137  //
138  G4int TotNbofEvents = numberOfEvent;
139  if (TotNbofEvents == 0) return;
140 
141  G4double EnergyBalance = fEnergyDeposit + fEnergyLeak[0] + fEnergyLeak[1];
142  EnergyBalance /= TotNbofEvents;
143 
144  fEnergyDeposit /= TotNbofEvents; fEnergyDeposit2 /= TotNbofEvents;
145  G4double rmsEdep = fEnergyDeposit2 - fEnergyDeposit*fEnergyDeposit;
146  if (rmsEdep>0.) rmsEdep = std::sqrt(rmsEdep/TotNbofEvents);
147  else rmsEdep = 0.;
148 
149  fTrakLenCharged /= TotNbofEvents; fTrakLenCharged2 /= TotNbofEvents;
150  G4double rmsTLCh = fTrakLenCharged2 - fTrakLenCharged*fTrakLenCharged;
151  if (rmsTLCh>0.) rmsTLCh = std::sqrt(rmsTLCh/TotNbofEvents);
152  else rmsTLCh = 0.;
153 
154  fTrakLenNeutral /= TotNbofEvents; fTrakLenNeutral2 /= TotNbofEvents;
155  G4double rmsTLNe = fTrakLenNeutral2 - fTrakLenNeutral*fTrakLenNeutral;
156  if (rmsTLNe>0.) rmsTLNe = std::sqrt(rmsTLNe/TotNbofEvents);
157  else rmsTLNe = 0.;
158 
159  fNbStepsCharged /= TotNbofEvents; fNbStepsCharged2 /= TotNbofEvents;
160  G4double rmsStCh = fNbStepsCharged2 - fNbStepsCharged*fNbStepsCharged;
161  if (rmsStCh>0.) rmsStCh = std::sqrt(rmsStCh/TotNbofEvents);
162  else rmsStCh = 0.;
163 
164  fNbStepsNeutral /= TotNbofEvents; fNbStepsNeutral2 /= TotNbofEvents;
165  G4double rmsStNe = fNbStepsNeutral2 - fNbStepsNeutral*fNbStepsNeutral;
166  if (rmsStNe>0.) rmsStNe = std::sqrt(rmsStNe/TotNbofEvents);
167  else rmsStNe = 0.;
168 
169  G4double Gamma = (G4double)fNbGamma/TotNbofEvents;
170  G4double Elect = (G4double)fNbElect/TotNbofEvents;
171  G4double Posit = (G4double)fNbPosit/TotNbofEvents;
172 
173  G4double transmit[2];
174  transmit[0] = 100.*fTransmit[0]/TotNbofEvents;
175  transmit[1] = 100.*fTransmit[1]/TotNbofEvents;
176 
177  G4double reflect[2];
178  reflect[0] = 100.*fReflect[0]/TotNbofEvents;
179  reflect[1] = 100.*fReflect[1]/TotNbofEvents;
180 
181  G4double rmsMsc = 0., tailMsc = 0.;
182  if (fMscEntryCentral > 0) {
183  fMscProjecTheta /= fMscEntryCentral; fMscProjecTheta2 /= fMscEntryCentral;
184  rmsMsc = fMscProjecTheta2 - fMscProjecTheta*fMscProjecTheta;
185  if (rmsMsc > 0.) { rmsMsc = std::sqrt(rmsMsc); }
186  if(fTransmit[1] > 0.0) {
187  tailMsc = 100.- (100.*fMscEntryCentral)/(2*fTransmit[1]);
188  }
189  }
190 
191  fEnergyLeak[0] /= TotNbofEvents; fEnergyLeak2[0] /= TotNbofEvents;
192  G4double rmsEl0 = fEnergyLeak2[0] - fEnergyLeak[0]*fEnergyLeak[0];
193  if (rmsEl0>0.) rmsEl0 = std::sqrt(rmsEl0/TotNbofEvents);
194  else rmsEl0 = 0.;
195 
196  fEnergyLeak[1] /= TotNbofEvents; fEnergyLeak2[1] /= TotNbofEvents;
197  G4double rmsEl1 = fEnergyLeak2[1] - fEnergyLeak[1]*fEnergyLeak[1];
198  if (rmsEl1>0.) rmsEl1 = std::sqrt(rmsEl1/TotNbofEvents);
199  else rmsEl1 = 0.;
200 
201 
202  //Stopping Power from input Table.
203  //
204  G4Material* material = fDetector->GetAbsorberMaterial();
205  G4double length = fDetector->GetAbsorberThickness();
206  G4double density = material->GetDensity();
207  G4String partName = fParticle->GetParticleName();
208 
209  G4EmCalculator emCalculator;
210  G4double dEdxTable = 0., dEdxFull = 0.;
211  if (fParticle->GetPDGCharge()!= 0.) {
212  dEdxTable = emCalculator.GetDEDX(fEkin,fParticle,material);
213  dEdxFull = emCalculator.ComputeTotalDEDX(fEkin,fParticle,material);
214  }
215  G4double stopTable = dEdxTable/density;
216  G4double stopFull = dEdxFull /density;
217 
218  //Stopping Power from simulation.
219  //
220  G4double meandEdx = fEnergyDeposit/length;
221  G4double stopPower = meandEdx/density;
222 
223  G4cout << "\n ======================== run summary ======================\n";
224 
225  G4int prec = G4cout.precision(3);
226 
227  G4cout << "\n The run was " << TotNbofEvents << " " << partName << " of "
228  << G4BestUnit(fEkin,"Energy") << " through "
229  << G4BestUnit(length,"Length") << " of "
230  << material->GetName() << " (density: "
231  << G4BestUnit(density,"Volumic Mass") << ")" << G4endl;
232 
233  G4cout.precision(4);
234 
235  G4cout << "\n Total energy deposit in absorber per event = "
236  << G4BestUnit(fEnergyDeposit,"Energy") << " +- "
237  << G4BestUnit(rmsEdep, "Energy")
238  << G4endl;
239 
240  G4cout << "\n -----> Mean dE/dx = " << meandEdx/(MeV/cm) << " MeV/cm"
241  << "\t(" << stopPower/(MeV*cm2/g) << " MeV*cm2/g)"
242  << G4endl;
243 
244  G4cout << "\n From formulas :" << G4endl;
245  G4cout << " restricted dEdx = " << dEdxTable/(MeV/cm) << " MeV/cm"
246  << "\t(" << stopTable/(MeV*cm2/g) << " MeV*cm2/g)"
247  << G4endl;
248 
249  G4cout << " full dEdx = " << dEdxFull/(MeV/cm) << " MeV/cm"
250  << "\t(" << stopFull/(MeV*cm2/g) << " MeV*cm2/g)"
251  << G4endl;
252 
253  G4cout << "\n Leakage : primary = "
254  << G4BestUnit(fEnergyLeak[0],"Energy") << " +- "
255  << G4BestUnit(rmsEl0, "Energy")
256  << " secondaries = "
257  << G4BestUnit(fEnergyLeak[1],"Energy") << " +- "
258  << G4BestUnit(rmsEl1, "Energy")
259  << G4endl;
260 
261  G4cout << " Energy balance : edep + eleak = "
262  << G4BestUnit(EnergyBalance,"Energy")
263  << G4endl;
264 
265  G4cout << "\n Total track length (charged) in absorber per event = "
266  << G4BestUnit(fTrakLenCharged,"Length") << " +- "
267  << G4BestUnit(rmsTLCh, "Length") << G4endl;
268 
269  G4cout << " Total track length (neutral) in absorber per event = "
270  << G4BestUnit(fTrakLenNeutral,"Length") << " +- "
271  << G4BestUnit(rmsTLNe, "Length") << G4endl;
272 
273  G4cout << "\n Number of steps (charged) in absorber per event = "
274  << fNbStepsCharged << " +- " << rmsStCh << G4endl;
275 
276  G4cout << " Number of steps (neutral) in absorber per event = "
277  << fNbStepsNeutral << " +- " << rmsStNe << G4endl;
278 
279  G4cout << "\n Number of secondaries per event : Gammas = " << Gamma
280  << "; electrons = " << Elect
281  << "; positrons = " << Posit << G4endl;
282 
283  G4cout << "\n Number of events with the primary particle transmitted = "
284  << transmit[1] << " %" << G4endl;
285 
286  G4cout << " Number of events with at least 1 particle transmitted "
287  << "(same charge as primary) = " << transmit[0] << " %" << G4endl;
288 
289  G4cout << "\n Number of events with the primary particle reflected = "
290  << reflect[1] << " %" << G4endl;
291 
292  G4cout << " Number of events with at least 1 particle reflected "
293  << "(same charge as primary) = " << reflect[0] << " %" << G4endl;
294 
295  // compute width of the Gaussian central part of the MultipleScattering
296  //
297  G4cout << "\n MultipleScattering:"
298  << "\n rms proj angle of transmit primary particle = "
299  << rmsMsc/mrad << " mrad (central part only)" << G4endl;
300 
301  G4cout << " computed theta0 (Highland formula) = "
302  << ComputeMscHighland()/mrad << " mrad" << G4endl;
303 
304  G4cout << " central part defined as +- "
305  << fMscThetaCentral/mrad << " mrad; "
306  << " Tail ratio = " << tailMsc << " %" << G4endl;
307 
308  // normalize histograms
309  //
310  G4AnalysisManager* analysisManager = G4AnalysisManager::Instance();
311 
312  G4int ih = 1;
313  G4double binWidth = analysisManager->GetH1Width(ih);
314  G4double unit = analysisManager->GetH1Unit(ih);
315  G4double fac = unit/(TotNbofEvents*binWidth);
316  analysisManager->ScaleH1(ih,fac);
317 
318  ih = 10;
319  binWidth = analysisManager->GetH1Width(ih);
320  unit = analysisManager->GetH1Unit(ih);
321  fac = unit/(TotNbofEvents*binWidth);
322  analysisManager->ScaleH1(ih,fac);
323 
324  ih = 12;
325  analysisManager->ScaleH1(ih,1./TotNbofEvents);
326 
327  // reset default precision
328  G4cout.precision(prec);
329 }
330 
331 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
332 
334 {
335  //compute the width of the Gaussian central part of the MultipleScattering
336  //projected angular distribution.
337  //Eur. Phys. Jour. C15 (2000) page 166, formule 23.9
338 
339  G4double t = (fDetector->GetAbsorberThickness())
340  /(fDetector->GetAbsorberMaterial()->GetRadlen());
341  if (t < DBL_MIN) return 0.;
342 
343  G4double T = fEkin;
344  G4double M = fParticle->GetPDGMass();
345  G4double z = std::abs(fParticle->GetPDGCharge()/eplus);
346 
347  G4double bpc = T*(T+2*M)/(T+M);
348  G4double teta0 = 13.6*MeV*z*std::sqrt(t)*(1.+0.038*std::log(t))/bpc;
349  return teta0;
350 }
351 
352 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
virtual void Merge(const G4Run *)
Definition: G4Run.cc:54
G4int numberOfEvent
Definition: G4Run.hh:59
static constexpr double cm2
Definition: G4SIunits.hh:120
Run()
Definition: Run.cc:43
const G4String & GetName() const
Definition: G4Material.hh:178
G4double ComputeMscHighland()
Definition: Run.cc:333
G4double GetDensity() const
Definition: G4Material.hh:180
G4double ComputeTotalDEDX(G4double kinEnergy, const G4ParticleDefinition *, const G4Material *, G4double cut=DBL_MAX)
#define G4BestUnit(a, b)
#define G4_USE_G4BESTUNIT_FOR_VERBOSE 1
int G4int
Definition: G4Types.hh:78
const G4String & GetParticleName() const
string material
Definition: eplot.py:19
function g(Y1, Y2, PT2)
Definition: hijing1.383.f:5205
static const double prec
Definition: RanecuEngine.cc:58
G4GLOB_DLL std::ostream G4cout
G4Material * GetAbsorberMaterial()
G4double GetDEDX(G4double kinEnergy, const G4ParticleDefinition *, const G4Material *, const G4Region *r=nullptr)
static constexpr double cm
Definition: G4SIunits.hh:119
static constexpr double eplus
Definition: G4SIunits.hh:199
Definition: G4Run.hh:46
static constexpr double mrad
Definition: G4SIunits.hh:150
G4double GetRadlen() const
Definition: G4Material.hh:220
G4double GetPDGMass() const
G4double energy(const ThreeVector &p, const G4double m)
#define DBL_MIN
Definition: templates.hh:75
tuple z
Definition: test.py:28
static const G4double fac
Detector construction class to define materials and geometry.
void EndOfRun()
Definition: Run.cc:147
#define G4endl
Definition: G4ios.hh:61
static constexpr double MeV
Definition: G4SIunits.hh:214
void SetPrimary(G4ParticleDefinition *particle, G4double energy)
Definition: Run.cc:77
double G4double
Definition: G4Types.hh:76
Definition: Run.hh:46
G4double GetPDGCharge() const
virtual void Merge(const G4Run *)
Definition: Run.cc:115
G4CsvAnalysisManager G4AnalysisManager
Definition: g4csv_defs.hh:77
~Run()
Definition: Run.cc:72