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