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 #include "EmAcceptance.hh"
38 
39 #include "G4ParticleTable.hh"
40 #include "G4ParticleDefinition.hh"
41 #include "G4Track.hh"
42 #include "G4Gamma.hh"
43 #include "G4Electron.hh"
44 #include "G4Positron.hh"
45 
46 #include "G4UnitsTable.hh"
47 #include "G4SystemOfUnits.hh"
48 
49 #include <iomanip>
50 
51 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
52 
53 Run::Run(DetectorConstruction* det)
54 : G4Run(),
55  fDetector(det),
56  fParticle(nullptr), fEkin(0.),
57  fChargedStep(0), fNeutralStep(0),
58  fN_gamma(0), fN_elec(0), fN_pos(0),
59  fApplyLimit(false)
60 {
61  //initialize cumulative quantities
62  //
63  for (G4int k=0; k<kMaxAbsor; k++) {
64  fSumEAbs[k] = fSum2EAbs[k] = fSumLAbs[k] = fSum2LAbs[k] = 0.;
65  fEnergyDeposit[k].clear();
66  fEdeptrue[k] = fRmstrue[k] = 1.;
67  fLimittrue[k] = DBL_MAX;
68  }
69 
70  //initialize Eflow
71  //
72  G4int nbPlanes = (fDetector->GetNbOfLayers())*(fDetector->GetNbOfAbsor()) + 2;
73  fEnergyFlow.resize(nbPlanes);
74  fLateralEleak.resize(nbPlanes);
75  for (G4int k=0; k<nbPlanes; k++) {fEnergyFlow[k] = fLateralEleak[k] = 0.; }
76 }
77 
78 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
79 
80 Run::~Run()
81 { }
82 
83 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
84 
85 void Run::SetPrimary(G4ParticleDefinition* particle, G4double energy)
86 {
87  fParticle = particle;
88  fEkin = energy;
89 }
90 
91 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
92 
93 void Run::FillPerEvent(G4int kAbs, G4double EAbs, G4double LAbs)
94 {
95  //accumulate statistic with restriction
96  //
97  if(fApplyLimit) fEnergyDeposit[kAbs].push_back(EAbs);
98  fSumEAbs[kAbs] += EAbs; fSum2EAbs[kAbs] += EAbs*EAbs;
99  fSumLAbs[kAbs] += LAbs; fSum2LAbs[kAbs] += LAbs*LAbs;
100 }
101 
102 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
103 
104 void Run::SumEnergyFlow(G4int plane, G4double Eflow)
105 {
106  fEnergyFlow[plane] += Eflow;
107 }
108 
109 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
110 
111 void Run::SumLateralEleak(G4int cell, G4double Eflow)
112 {
113  fLateralEleak[cell] += Eflow;
114 }
115 
116 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
117 
118 void Run::AddChargedStep()
119 {
120  fChargedStep += 1.0;
121 }
122 
123 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
124 
125 void Run::AddNeutralStep()
126 {
127  fNeutralStep += 1.0;
128 }
129 
130 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
131 
132 void Run::AddSecondaryTrack(const G4Track* track)
133 {
134  const G4ParticleDefinition* d = track->GetDefinition();
135  if(d == G4Gamma::Gamma()) { ++fN_gamma; }
136  else if (d == G4Electron::Electron()) { ++fN_elec; }
137  else if (d == G4Positron::Positron()) { ++fN_pos; }
138 }
139 
140 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
141 
142 void Run::Merge(const G4Run* run)
143 {
144  const Run* localRun = static_cast<const Run*>(run);
145 
146  // pass information about primary particle
147  fParticle = localRun->fParticle;
148  fEkin = localRun->fEkin;
149 
150  // accumulate sums
151  //
152  for (G4int k=0; k<kMaxAbsor; k++) {
153  fSumEAbs[k] += localRun->fSumEAbs[k];
154  fSum2EAbs[k] += localRun->fSum2EAbs[k];
155  fSumLAbs[k] += localRun->fSumLAbs[k];
156  fSum2LAbs[k] += localRun->fSum2LAbs[k];
157  }
158 
159  G4int nbPlanes = (fDetector->GetNbOfLayers())*(fDetector->GetNbOfAbsor()) + 2;
160  for (G4int k=0; k<nbPlanes; k++) {
161  fEnergyFlow[k] += localRun->fEnergyFlow[k];
162  fLateralEleak[k] += localRun->fLateralEleak[k];
163  }
164 
165 
166  fChargedStep += localRun->fChargedStep;
167  fNeutralStep += localRun->fNeutralStep;
168 
169  fN_gamma += localRun->fN_gamma;
170  fN_elec += localRun->fN_elec;
171  fN_pos += localRun->fN_pos;
172 
173  fApplyLimit = localRun->fApplyLimit;
174 
175  for (G4int k=0; k<kMaxAbsor; k++) {
176  fEdeptrue[k] = localRun->fEdeptrue[k];
177  fRmstrue[k] = localRun->fRmstrue[k];
178  fLimittrue[k] = localRun->fLimittrue[k];
179  }
180 
181  G4Run::Merge(run);
182 }
183 
184 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
185 
186 void Run::EndOfRun()
187 {
188  G4int nEvt = numberOfEvent;
189  G4double norm = G4double(nEvt);
190  if(norm > 0) norm = 1./norm;
191  G4double qnorm = std::sqrt(norm);
192 
193  fChargedStep *= norm;
194  fNeutralStep *= norm;
195 
196  //compute and print statistic
197  //
198  G4double beamEnergy = fEkin;
199  G4double sqbeam = std::sqrt(beamEnergy/GeV);
200 
201  G4double MeanEAbs,MeanEAbs2,rmsEAbs,resolution,rmsres;
202  G4double MeanLAbs,MeanLAbs2,rmsLAbs;
203 
204  std::ios::fmtflags mode = G4cout.flags();
205  G4int prec = G4cout.precision(2);
206  G4cout << "\n------------------------------------------------------------\n";
207  G4cout << std::setw(14) << "material"
208  << std::setw(17) << "Edep RMS"
209  << std::setw(33) << "sqrt(E0(GeV))*rmsE/Emean"
210  << std::setw(23) << "total tracklen \n \n";
211 
212  for (G4int k=1; k<=fDetector->GetNbOfAbsor(); k++)
213  {
214  MeanEAbs = fSumEAbs[k]*norm;
215  MeanEAbs2 = fSum2EAbs[k]*norm;
216  rmsEAbs = std::sqrt(std::abs(MeanEAbs2 - MeanEAbs*MeanEAbs));
217  //G4cout << "k= " << k << " RMS= " << rmsEAbs
218  // << " fApplyLimit: " << fApplyLimit << G4endl;
219  if(fApplyLimit) {
220  G4int nn = 0;
221  G4double sume = 0.0;
222  G4double sume2 = 0.0;
223  // compute trancated means
224  G4double lim = rmsEAbs * 2.5;
225  for(G4int i=0; i<nEvt; i++) {
226  G4double e = (fEnergyDeposit[k])[i];
227  if(std::abs(e - MeanEAbs) < lim) {
228  sume += e;
229  sume2 += e*e;
230  nn++;
231  }
232  }
233  G4double norm1 = G4double(nn);
234  if(norm1 > 0.0) norm1 = 1.0/norm1;
235  MeanEAbs = sume*norm1;
236  MeanEAbs2 = sume2*norm1;
237  rmsEAbs = std::sqrt(std::abs(MeanEAbs2 - MeanEAbs*MeanEAbs));
238  }
239 
240  resolution= 100.*sqbeam*rmsEAbs/MeanEAbs;
241  rmsres = resolution*qnorm;
242 
243  // Save mean and RMS
244  fSumEAbs[k] = MeanEAbs;
245  fSum2EAbs[k] = rmsEAbs;
246 
247  MeanLAbs = fSumLAbs[k]*norm;
248  MeanLAbs2 = fSum2LAbs[k]*norm;
249  rmsLAbs = std::sqrt(std::abs(MeanLAbs2 - MeanLAbs*MeanLAbs));
250 
251  //print
252  //
253  G4cout
254  << std::setw(14) << fDetector->GetAbsorMaterial(k)->GetName() << ": "
255  << std::setprecision(5)
256  << std::setw(6) << G4BestUnit(MeanEAbs,"Energy") << " : "
257  << std::setprecision(4)
258  << std::setw(5) << G4BestUnit( rmsEAbs,"Energy")
259  << std::setw(10) << resolution << " +- "
260  << std::setw(5) << rmsres << " %"
261  << std::setprecision(3)
262  << std::setw(10) << G4BestUnit(MeanLAbs,"Length") << " +- "
263  << std::setw(4) << G4BestUnit( rmsLAbs,"Length")
264  << G4endl;
265  }
266  G4cout << "\n------------------------------------------------------------\n";
267 
268  G4cout << " Beam particle "
269  << fParticle->GetParticleName()
270  << " E = " << G4BestUnit(beamEnergy,"Energy") << G4endl;
271  G4cout << " Mean number of gamma " << (G4double)fN_gamma*norm << G4endl;
272  G4cout << " Mean number of e- " << (G4double)fN_elec*norm << G4endl;
273  G4cout << " Mean number of e+ " << (G4double)fN_pos*norm << G4endl;
274  G4cout << std::setprecision(6)
275  << " Mean number of charged steps " << fChargedStep << G4endl;
276  G4cout << " Mean number of neutral steps " << fNeutralStep << G4endl;
277  G4cout << "------------------------------------------------------------\n";
278 
279  //Energy flow
280  //
281  G4AnalysisManager* analysis = G4AnalysisManager::Instance();
282  G4int Idmax = (fDetector->GetNbOfLayers())*(fDetector->GetNbOfAbsor());
283  for (G4int Id=1; Id<=Idmax+1; Id++) {
284  analysis->FillH1(2*kMaxAbsor+1, (G4double)Id, fEnergyFlow[Id]);
285  analysis->FillH1(2*kMaxAbsor+2, (G4double)Id, fLateralEleak[Id]);
286  }
287 
288  //Energy deposit from energy flow balance
289  //
290  G4double EdepTot[kMaxAbsor];
291  for (G4int k=0; k<kMaxAbsor; k++) EdepTot[k] = 0.;
292 
293  G4int nbOfAbsor = fDetector->GetNbOfAbsor();
294  for (G4int Id=1; Id<=Idmax; Id++) {
295  G4int iAbsor = Id%nbOfAbsor; if (iAbsor==0) iAbsor = nbOfAbsor;
296  EdepTot[iAbsor] += (fEnergyFlow[Id]-fEnergyFlow[Id+1]-fLateralEleak[Id]);
297  }
298 
299  G4cout << std::setprecision(3)
300  << "\n Energy deposition from Energy flow balance : \n"
301  << std::setw(10) << " material \t Total Edep \n \n";
302  G4cout.precision(6);
303 
304  for (G4int k=1; k<=nbOfAbsor; k++) {
305  EdepTot [k] *= norm;
306  G4cout << std::setw(10) << fDetector->GetAbsorMaterial(k)->GetName() << ":"
307  << "\t " << G4BestUnit(EdepTot [k],"Energy") << "\n";
308  }
309 
310  G4cout << "\n------------------------------------------------------------\n"
311  << G4endl;
312 
313  // Acceptance
314  EmAcceptance acc;
315  G4bool isStarted = false;
316  for (G4int j=1; j<=fDetector->GetNbOfAbsor(); j++) {
317  if (fLimittrue[j] < DBL_MAX) {
318  if (!isStarted) {
319  acc.BeginOfAcceptance("Sampling Calorimeter",nEvt);
320  isStarted = true;
321  }
322  MeanEAbs = fSumEAbs[j];
323  rmsEAbs = fSum2EAbs[j];
324  G4String mat = fDetector->GetAbsorMaterial(j)->GetName();
325  acc.EmAcceptanceGauss("Edep"+mat, nEvt, MeanEAbs,
326  fEdeptrue[j], fRmstrue[j], fLimittrue[j]);
327  acc.EmAcceptanceGauss("Erms"+mat, nEvt, rmsEAbs,
328  fRmstrue[j], fRmstrue[j], 2.0*fLimittrue[j]);
329  }
330  }
331  if(isStarted) acc.EndOfAcceptance();
332 
333  //normalize histograms
334  //
335  for (G4int ih = kMaxAbsor+1; ih < kMaxHisto; ih++) {
336  analysis->ScaleH1(ih,norm/MeV);
337  }
338 
339  G4cout.setf(mode,std::ios::floatfield);
340  G4cout.precision(prec);
341 }
342 
343 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
344 
345 void Run::SetEdepAndRMS(G4int i, G4double edep, G4double rms, G4double lim)
346 {
347  if (i>=0 && i<kMaxAbsor) {
348  fEdeptrue [i] = edep;
349  fRmstrue [i] = rms;
350  fLimittrue[i] = lim;
351  }
352 }
353 
354 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
355 
356 void Run::SetApplyLimit(G4bool val)
357 {
358  fApplyLimit = val;
359 }
360 
361 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......