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 #include "PrimaryGeneratorAction.hh"
36 #include "HistoManager.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 
46 Run::Run(DetectorConstruction* det)
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 
75 void Run::SetPrimary(G4ParticleDefinition* particle, G4double energy)
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  fMscProjecTheta += localRun->fMscProjecTheta;
109  fMscProjecTheta2 += localRun->fMscProjecTheta2;
110 
111 
112  fNbGamma += localRun->fNbGamma;
113  fNbElect += localRun->fNbElect;
114  fNbPosit += localRun->fNbPosit;
115 
116  fTransmit[0] += localRun->fTransmit[0];
117  fTransmit[1] += localRun->fTransmit[1];
118  fReflect[0] += localRun->fReflect[0];
119  fReflect[1] += localRun->fReflect[1];
120 
121  fMscEntryCentral += localRun->fMscEntryCentral;
122 
123  fEnergyLeak[0] += localRun->fEnergyLeak[0];
124  fEnergyLeak[1] += localRun->fEnergyLeak[1];
125  fEnergyLeak2[0] += localRun->fEnergyLeak2[0];
126  fEnergyLeak2[1] += localRun->fEnergyLeak2[1];
127 
128  G4Run::Merge(run);
129 }
130 
131 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
132 
133 void Run::EndOfRun()
134 {
135  // compute mean and rms
136  //
137  G4int TotNbofEvents = numberOfEvent;
138  if (TotNbofEvents == 0) return;
139 
140  G4double EnergyBalance = fEnergyDeposit + fEnergyLeak[0] + fEnergyLeak[1];
141  EnergyBalance /= TotNbofEvents;
142 
143  fEnergyDeposit /= TotNbofEvents; fEnergyDeposit2 /= TotNbofEvents;
144  G4double rmsEdep = fEnergyDeposit2 - fEnergyDeposit*fEnergyDeposit;
145  if (rmsEdep>0.) rmsEdep = std::sqrt(rmsEdep/TotNbofEvents);
146  else rmsEdep = 0.;
147 
148  fTrakLenCharged /= TotNbofEvents; fTrakLenCharged2 /= TotNbofEvents;
149  G4double rmsTLCh = fTrakLenCharged2 - fTrakLenCharged*fTrakLenCharged;
150  if (rmsTLCh>0.) rmsTLCh = std::sqrt(rmsTLCh/TotNbofEvents);
151  else rmsTLCh = 0.;
152 
153  fTrakLenNeutral /= TotNbofEvents; fTrakLenNeutral2 /= TotNbofEvents;
154  G4double rmsTLNe = fTrakLenNeutral2 - fTrakLenNeutral*fTrakLenNeutral;
155  if (rmsTLNe>0.) rmsTLNe = std::sqrt(rmsTLNe/TotNbofEvents);
156  else rmsTLNe = 0.;
157 
158  fNbStepsCharged /= TotNbofEvents; fNbStepsCharged2 /= TotNbofEvents;
159  G4double rmsStCh = fNbStepsCharged2 - fNbStepsCharged*fNbStepsCharged;
160  if (rmsStCh>0.) rmsStCh = std::sqrt(rmsStCh/TotNbofEvents);
161  else rmsStCh = 0.;
162 
163  fNbStepsNeutral /= TotNbofEvents; fNbStepsNeutral2 /= TotNbofEvents;
164  G4double rmsStNe = fNbStepsNeutral2 - fNbStepsNeutral*fNbStepsNeutral;
165  if (rmsStNe>0.) rmsStNe = std::sqrt(rmsStNe/TotNbofEvents);
166  else rmsStNe = 0.;
167 
168  G4double Gamma = (G4double)fNbGamma/TotNbofEvents;
169  G4double Elect = (G4double)fNbElect/TotNbofEvents;
170  G4double Posit = (G4double)fNbPosit/TotNbofEvents;
171 
172  G4double transmit[2];
173  transmit[0] = 100.*fTransmit[0]/TotNbofEvents;
174  transmit[1] = 100.*fTransmit[1]/TotNbofEvents;
175 
176  G4double reflect[2];
177  reflect[0] = 100.*fReflect[0]/TotNbofEvents;
178  reflect[1] = 100.*fReflect[1]/TotNbofEvents;
179 
180  G4double rmsMsc = 0., tailMsc = 0.;
181  if (fMscEntryCentral > 0) {
182  fMscProjecTheta /= fMscEntryCentral; fMscProjecTheta2 /= fMscEntryCentral;
183  rmsMsc = fMscProjecTheta2 - fMscProjecTheta*fMscProjecTheta;
184  if (rmsMsc > 0.) { rmsMsc = std::sqrt(rmsMsc); }
185  if(fTransmit[1] > 0.0) {
186  tailMsc = 100.- (100.*fMscEntryCentral)/(2*fTransmit[1]);
187  }
188  }
189 
190  fEnergyLeak[0] /= TotNbofEvents; fEnergyLeak2[0] /= TotNbofEvents;
191  G4double rmsEl0 = fEnergyLeak2[0] - fEnergyLeak[0]*fEnergyLeak[0];
192  if (rmsEl0>0.) rmsEl0 = std::sqrt(rmsEl0/TotNbofEvents);
193  else rmsEl0 = 0.;
194 
195  fEnergyLeak[1] /= TotNbofEvents; fEnergyLeak2[1] /= TotNbofEvents;
196  G4double rmsEl1 = fEnergyLeak2[1] - fEnergyLeak[1]*fEnergyLeak[1];
197  if (rmsEl1>0.) rmsEl1 = std::sqrt(rmsEl1/TotNbofEvents);
198  else rmsEl1 = 0.;
199 
200 
201  //Stopping Power from input Table.
202  //
203  const G4Material* material = fDetector->GetAbsorberMaterial();
204  G4double length = fDetector->GetAbsorberThickness();
205  G4double density = material->GetDensity();
206  G4String partName = fParticle->GetParticleName();
207 
208  G4EmCalculator emCalculator;
209  G4double dEdxTable = 0., dEdxFull = 0.;
210  if (fParticle->GetPDGCharge()!= 0.) {
211  dEdxTable = emCalculator.GetDEDX(fEkin,fParticle,material);
212  dEdxFull = emCalculator.ComputeTotalDEDX(fEkin,fParticle,material);
213  }
214  G4double stopTable = dEdxTable/density;
215  G4double stopFull = dEdxFull /density;
216 
217  //Stopping Power from simulation.
218  //
219  G4double meandEdx = fEnergyDeposit/length;
220  G4double stopPower = meandEdx/density;
221 
222  G4cout << "\n ======================== run summary ======================\n";
223 
224  G4int prec = G4cout.precision(3);
225 
226  G4cout << "\n The run was " << TotNbofEvents << " " << partName << " of "
227  << G4BestUnit(fEkin,"Energy") << " through "
228  << G4BestUnit(length,"Length") << " of "
229  << material->GetName() << " (density: "
230  << G4BestUnit(density,"Volumic Mass") << ")" << G4endl;
231 
232  G4cout.precision(4);
233 
234  G4cout << "\n Total energy deposit in absorber per event = "
235  << G4BestUnit(fEnergyDeposit,"Energy") << " +- "
236  << G4BestUnit(rmsEdep, "Energy")
237  << G4endl;
238 
239  G4cout << "\n -----> Mean dE/dx = " << meandEdx/(MeV/cm) << " MeV/cm"
240  << "\t(" << stopPower/(MeV*cm2/g) << " MeV*cm2/g)"
241  << G4endl;
242 
243  G4cout << "\n From formulas :" << G4endl;
244  G4cout << " restricted dEdx = " << dEdxTable/(MeV/cm) << " MeV/cm"
245  << "\t(" << stopTable/(MeV*cm2/g) << " MeV*cm2/g)"
246  << G4endl;
247 
248  G4cout << " full dEdx = " << dEdxFull/(MeV/cm) << " MeV/cm"
249  << "\t(" << stopFull/(MeV*cm2/g) << " MeV*cm2/g)"
250  << G4endl;
251 
252  G4cout << "\n Leakage : primary = "
253  << G4BestUnit(fEnergyLeak[0],"Energy") << " +- "
254  << G4BestUnit(rmsEl0, "Energy")
255  << " secondaries = "
256  << G4BestUnit(fEnergyLeak[1],"Energy") << " +- "
257  << G4BestUnit(rmsEl1, "Energy")
258  << G4endl;
259 
260  G4cout << " Energy balance : edep + eleak = "
261  << G4BestUnit(EnergyBalance,"Energy")
262  << G4endl;
263 
264  G4cout << "\n Total track length (charged) in absorber per event = "
265  << G4BestUnit(fTrakLenCharged,"Length") << " +- "
266  << G4BestUnit(rmsTLCh, "Length") << G4endl;
267 
268  G4cout << " Total track length (neutral) in absorber per event = "
269  << G4BestUnit(fTrakLenNeutral,"Length") << " +- "
270  << G4BestUnit(rmsTLNe, "Length") << G4endl;
271 
272  G4cout << "\n Number of steps (charged) in absorber per event = "
273  << fNbStepsCharged << " +- " << rmsStCh << G4endl;
274 
275  G4cout << " Number of steps (neutral) in absorber per event = "
276  << fNbStepsNeutral << " +- " << rmsStNe << G4endl;
277 
278  G4cout << "\n Number of secondaries per event : Gammas = " << Gamma
279  << "; electrons = " << Elect
280  << "; positrons = " << Posit << G4endl;
281 
282  G4cout << "\n Number of events with the primary particle transmitted = "
283  << transmit[1] << " %" << G4endl;
284 
285  G4cout << " Number of events with at least 1 particle transmitted "
286  << "(same charge as primary) = " << transmit[0] << " %" << G4endl;
287 
288  G4cout << "\n Number of events with the primary particle reflected = "
289  << reflect[1] << " %" << G4endl;
290 
291  G4cout << " Number of events with at least 1 particle reflected "
292  << "(same charge as primary) = " << reflect[0] << " %" << G4endl;
293 
294  // compute width of the Gaussian central part of the MultipleScattering
295  //
296  G4cout << "\n MultipleScattering:"
297  << "\n rms proj angle of transmit primary particle = "
298  << rmsMsc/mrad << " mrad (central part only)" << G4endl;
299 
300  G4cout << " computed theta0 (Highland formula) = "
301  << ComputeMscHighland()/mrad << " mrad" << G4endl;
302 
303  G4cout << " central part defined as +- "
304  << fMscThetaCentral/mrad << " mrad; "
305  << " Tail ratio = " << tailMsc << " %" << G4endl;
306 
307  // normalize histograms
308  //
309  G4AnalysisManager* analysisManager = G4AnalysisManager::Instance();
310 
311  G4int ih = 1;
312  G4double binWidth = analysisManager->GetH1Width(ih);
313  G4double fac = 1./(TotNbofEvents*binWidth);
314  analysisManager->ScaleH1(ih,fac);
315 
316  ih = 10;
317  binWidth = analysisManager->GetH1Width(ih);
318  fac = 1./(TotNbofEvents*binWidth);
319  analysisManager->ScaleH1(ih,fac);
320 
321  ih = 12;
322  analysisManager->ScaleH1(ih,1./TotNbofEvents);
323 
324  // reset default precision
325  G4cout.precision(prec);
326 }
327 
328 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
329 
330 G4double Run::ComputeMscHighland()
331 {
332  //compute the width of the Gaussian central part of the MultipleScattering
333  //projected angular distribution.
334  //Eur. Phys. Jour. C15 (2000) page 166, formule 23.9
335 
336  G4double t = (fDetector->GetAbsorberThickness())
337  /(fDetector->GetAbsorberMaterial()->GetRadlen());
338  if (t < DBL_MIN) return 0.;
339 
340  G4double T = fEkin;
341  G4double M = fParticle->GetPDGMass();
342  G4double z = std::abs(fParticle->GetPDGCharge()/eplus);
343 
344  G4double bpc = T*(T+2*M)/(T+M);
345  G4double teta0 = 13.6*MeV*z*std::sqrt(t)*(1.+0.038*std::log(t))/bpc;
346  return teta0;
347 }
348 
349 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......