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 
35 #include "DetectorConstruction.hh"
36 #include "PrimaryGeneratorAction.hh"
37 #include "EmAcceptance.hh"
38 
39 #include "G4Run.hh"
40 #include "G4UnitsTable.hh"
41 #include "G4SystemOfUnits.hh"
42 #include <iomanip>
43 
44 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
45 
46 Run::Run(DetectorConstruction* det, PrimaryGeneratorAction* kin)
47  :G4Run(),fDet(det),fKin(kin),
48  f_nLbin(kMaxBin),f_nRbin(kMaxBin)
49 {
50  Reset();
51 }
52 
53 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
54 
55 void Run::Reset()
56 {
57  f_nLbin = fDet->GetnLtot();
58  f_dEdL.resize(f_nLbin);
59  fSumELongit.resize(f_nLbin);
60  fSumELongitCumul.resize(f_nLbin);
61  fSumE2Longit.resize(f_nLbin);
62  fSumE2LongitCumul.resize(f_nLbin);
63 
64  f_nRbin = fDet->GetnRtot();
65  f_dEdR.resize(f_nRbin);
66  fSumERadial.resize(f_nRbin);
67  fSumERadialCumul.resize(f_nRbin);
68  fSumE2Radial.resize(f_nRbin);
69  fSumE2RadialCumul.resize(f_nRbin);
70 
71  fChargedStep = 0.0;
72  fNeutralStep = 0.0;
73 
74  fVerbose = 0;
75 
76  //initialize arrays of cumulative energy deposition
77  //
78  for (G4int i=0; i<f_nLbin; ++i) {
79  fSumELongit[i]=fSumE2Longit[i]=fSumELongitCumul[i]=fSumE2LongitCumul[i]=0.;
80  }
81  for (G4int j=0; j<f_nRbin; ++j) {
82  fSumERadial[j]=fSumE2Radial[j]=fSumERadialCumul[j]=fSumE2RadialCumul[j]=0.;
83  }
84  //initialize track length
85  fSumChargTrLength=fSum2ChargTrLength=fSumNeutrTrLength=fSum2NeutrTrLength=0.;
86 
87 }
88 
89 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
90 
91 Run::~Run()
92 {}
93 
94 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
95 
96 void Run::InitializePerEvent()
97 {
98  //initialize arrays of energy deposit per bin
99  for (G4int i=0; i<f_nLbin; ++i)
100  { f_dEdL[i] = 0.; }
101 
102  for (G4int j=0; j<f_nRbin; ++j)
103  { f_dEdR[j] = 0.; }
104 
105  //initialize tracklength
106  fChargTrLength = fNeutrTrLength = 0.;
107 }
108 
109 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
110 
111 void Run::FillPerEvent()
112 {
113  //accumulate statistic
114  //
115  G4double dLCumul = 0.;
116  for (G4int i=0; i<f_nLbin; ++i)
117  {
118  fSumELongit[i] += f_dEdL[i];
119  fSumE2Longit[i] += f_dEdL[i]*f_dEdL[i];
120  dLCumul += f_dEdL[i];
121  fSumELongitCumul[i] += dLCumul;
122  fSumE2LongitCumul[i] += dLCumul*dLCumul;
123  }
124 
125  G4double dRCumul = 0.;
126  for (G4int j=0; j<f_nRbin; j++)
127  {
128  fSumERadial[j] += f_dEdR[j];
129  fSumE2Radial[j] += f_dEdR[j]*f_dEdR[j];
130  dRCumul += f_dEdR[j];
131  fSumERadialCumul[j] += dRCumul;
132  fSumE2RadialCumul[j] += dRCumul*dRCumul;
133  }
134 
135  fSumChargTrLength += fChargTrLength;
136  fSum2ChargTrLength += fChargTrLength*fChargTrLength;
137  fSumNeutrTrLength += fNeutrTrLength;
138  fSum2NeutrTrLength += fNeutrTrLength*fNeutrTrLength;
139 
140  //fill histograms
141  //
142 
143  G4double Ekin=fKin->GetParticleGun()->GetParticleEnergy();
144  G4double mass=fKin->GetParticleGun()->GetParticleDefinition()->GetPDGMass();
145  G4double radl=fDet->GetMaterial()->GetRadlen();
146 
147  G4AnalysisManager* analysisManager = G4AnalysisManager::Instance();
148  analysisManager->FillH1(1, 100.*dLCumul/(Ekin+mass));
149  analysisManager->FillH1(2, fChargTrLength/radl);
150  analysisManager->FillH1(3, fNeutrTrLength/radl);
151 
152  //profiles
153  G4double norm = 100./(Ekin+mass);
154  G4double dLradl = fDet->GetdLradl();
155  for (G4int i=0; i<f_nLbin; ++i) {
156  G4double bin = (i+0.5)*dLradl;
157  analysisManager->FillP1(0, bin, norm*f_dEdL[i]/dLradl);
158  }
159  G4double dRradl = fDet->GetdRradl();
160  for (G4int j=0; j<f_nRbin; ++j) {
161  G4double bin = (j+0.5)*dRradl;
162  analysisManager->FillP1(1, bin, norm*f_dEdR[j]/dRradl);
163  }
164 }
165 
166 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
167 
168 void Run::Merge(const G4Run* run)
169 {
170  const Run* localRun = static_cast<const Run*>(run);
171 
172  fChargedStep += localRun->fChargedStep;
173  fNeutralStep += localRun->fNeutralStep;
174 
175  for (G4int i=0; i<f_nLbin; ++i) {
176  fSumELongit[i] += localRun->fSumELongit[i];
177  fSumE2Longit[i] += localRun->fSumE2Longit[i];
178  fSumELongitCumul[i] += localRun->fSumELongitCumul[i];
179  fSumE2LongitCumul[i] += localRun->fSumE2LongitCumul[i];
180  }
181 
182  for (G4int j=0; j<f_nRbin; ++j) {
183  fSumERadial[j] += localRun->fSumERadial[j];
184  fSumE2Radial[j] += localRun->fSumE2Radial[j];
185  fSumERadialCumul[j] += localRun->fSumERadialCumul[j];
186  fSumE2RadialCumul[j] += localRun->fSumE2RadialCumul[j];
187  }
188 
189  fSumChargTrLength += localRun->fSumChargTrLength;
190  fSum2ChargTrLength += localRun->fSum2ChargTrLength;
191  fSumNeutrTrLength += localRun->fSumNeutrTrLength;
192  fSum2NeutrTrLength += localRun->fSum2NeutrTrLength;
193 
194  G4Run::Merge(run);
195 }
196 
197 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
198 
199 void Run::EndOfRun(G4double edep, G4double rms, G4double& limit)
200 {
201  G4int NbOfEvents = GetNumberOfEvent();
202 
203  G4double kinEnergy = fKin->GetParticleGun()->GetParticleEnergy();
204  assert(NbOfEvents*kinEnergy > 0);
205 
206  fChargedStep /= G4double(NbOfEvents);
207  fNeutralStep /= G4double(NbOfEvents);
208 
209  G4double mass=fKin->GetParticleGun()->GetParticleDefinition()->GetPDGMass();
210  G4double norme = 100./(NbOfEvents*(kinEnergy+mass));
211 
212  //longitudinal
213  //
214  G4double dLradl = fDet->GetdLradl();
215 
216  MyVector MeanELongit(f_nLbin), rmsELongit(f_nLbin);
217  MyVector MeanELongitCumul(f_nLbin), rmsELongitCumul(f_nLbin);
218 
219  G4AnalysisManager* analysisManager = G4AnalysisManager::Instance();
220 
221  G4int i;
222  for (i=0; i<f_nLbin; ++i) {
223  MeanELongit[i] = norme*fSumELongit[i];
224  rmsELongit[i] =
225  norme*std::sqrt(std::abs(NbOfEvents*fSumE2Longit[i]
226  - fSumELongit[i]*fSumELongit[i]));
227 
228  MeanELongitCumul[i] = norme*fSumELongitCumul[i];
229  rmsELongitCumul[i] = norme*std::sqrt(std::abs(NbOfEvents*
230  fSumE2LongitCumul[i] - fSumELongitCumul[i]*fSumELongitCumul[i]));
231  G4double bin = (i+0.5)*dLradl;
232  analysisManager->FillH1(4, bin,MeanELongit[i]/dLradl);
233  analysisManager->FillH1(5, bin, rmsELongit[i]/dLradl);
234  bin = (i+1)*dLradl;
235  analysisManager->FillH1(6, bin,MeanELongitCumul[i]);
236  analysisManager->FillH1(7, bin, rmsELongitCumul[i]);
237  }
238 
239  //radial
240  //
241  G4double dRradl = fDet->GetdRradl();
242 
243  MyVector MeanERadial(f_nRbin), rmsERadial(f_nRbin);
244  MyVector MeanERadialCumul(f_nRbin), rmsERadialCumul(f_nRbin);
245 
246  for (i=0; i<f_nRbin; ++i) {
247  MeanERadial[i] = norme*fSumERadial[i];
248  rmsERadial[i] = norme*std::sqrt(std::abs(NbOfEvents*fSumE2Radial[i]
249  - fSumERadial[i]*fSumERadial[i]));
250 
251  MeanERadialCumul[i] = norme*fSumERadialCumul[i];
252  rmsERadialCumul[i] =
253  norme*std::sqrt(std::abs(NbOfEvents*fSumE2RadialCumul[i]
254  - fSumERadialCumul[i]*fSumERadialCumul[i]));
255 
256  G4double bin = (i+0.5)*dRradl;
257  analysisManager->FillH1(8, bin,MeanERadial[i]/dRradl);
258  analysisManager->FillH1(9, bin, rmsERadial[i]/dRradl);
259  bin = (i+1)*dRradl;
260  analysisManager->FillH1(10, bin,MeanERadialCumul[i]);
261  analysisManager->FillH1(11, bin, rmsERadialCumul[i]);
262  }
263 
264  //find Moliere confinement
265  //
266  const G4double EMoliere = 90.;
267  G4double iMoliere = 0.;
268  if ((MeanERadialCumul[0] <= EMoliere) &&
269  (MeanERadialCumul[f_nRbin-1] >= EMoliere)) {
270  G4int imin = 0;
271  while( (imin < f_nRbin-1) && (MeanERadialCumul[imin] < EMoliere) )
272  { ++imin; }
273 
274  G4double del = MeanERadialCumul[imin+1] - MeanERadialCumul[imin];
275  G4double ratio =
276  (del > 0.0) ? (EMoliere - MeanERadialCumul[imin])/del : 0.0;
277  iMoliere = 1. + imin + ratio;
278  }
279 
280  //track length
281  //
282  norme = 1./(NbOfEvents*(fDet->GetMaterial()->GetRadlen()));
283  G4double MeanChargTrLength = norme*fSumChargTrLength;
284  G4double rmsChargTrLength =
285  norme*std::sqrt(std::abs(NbOfEvents*fSum2ChargTrLength
286  - fSumChargTrLength*fSumChargTrLength));
287 
288  G4double MeanNeutrTrLength = norme*fSumNeutrTrLength;
289  G4double rmsNeutrTrLength =
290  norme*std::sqrt(std::abs(NbOfEvents*fSum2NeutrTrLength
291  - fSumNeutrTrLength*fSumNeutrTrLength));
292 
293  //print
294  std::ios::fmtflags mode = G4cout.flags();
295  G4cout.setf(std::ios::fixed,std::ios::floatfield);
296  G4int prec = G4cout.precision(2);
297 
298  if (fVerbose) {
299 
300  G4cout << " LOGITUDINAL PROFILE "
301  << " CUMULATIVE LOGITUDINAL PROFILE" << G4endl << G4endl;
302 
303  G4cout << " bin " << " Mean rms "
304  << " bin " << " Mean rms \n" << G4endl;
305 
306  for (i=0; i<f_nLbin; ++i) {
307  G4double inf=i*dLradl, sup=inf+dLradl;
308 
309  G4cout << std::setw(8) << inf << "->"
310  << std::setw(5) << sup << " radl: "
311  << std::setw(7) << MeanELongit[i] << "% "
312  << std::setw(9) << rmsELongit[i] << "% "
313  << " 0->" << std::setw(5) << sup << " radl: "
314  << std::setw(7) << MeanELongitCumul[i] << "% "
315  << std::setw(7) << rmsELongitCumul[i] << "% "
316  <<G4endl;
317  }
318 
319  G4cout << G4endl << G4endl << G4endl;
320 
321  G4cout << " RADIAL PROFILE "
322  << " CUMULATIVE RADIAL PROFILE" << G4endl << G4endl;
323 
324  G4cout << " bin " << " Mean rms "
325  << " bin " << " Mean rms \n" << G4endl;
326 
327  for (i=0; i<f_nRbin; ++i) {
328  G4double inf=i*dRradl, sup=inf+dRradl;
329 
330  G4cout << std::setw(8) << inf << "->"
331  << std::setw(5) << sup << " radl: "
332  << std::setw(7) << MeanERadial[i] << "% "
333  << std::setw(9) << rmsERadial[i] << "% "
334  << " 0->" << std::setw(5) << sup << " radl: "
335  << std::setw(7) << MeanERadialCumul[i] << "% "
336  << std::setw(7) << rmsERadialCumul[i] << "% "
337  <<G4endl;
338  }
339  }
340 
341  G4cout << "\n ===== SUMMARY ===== \n" << G4endl;
342 
343  G4cout << " Total number of events: " << NbOfEvents << "\n"
344  << " Mean number of charged steps: " << fChargedStep << G4endl;
345  G4cout << " Mean number of neutral steps: " << fNeutralStep
346  << "\n" << G4endl;
347 
348  G4cout << " energy deposit : "
349  << std::setw(7) << MeanELongitCumul[f_nLbin-1] << " % E0 +- "
350  << std::setw(7) << rmsELongitCumul[f_nLbin-1] << " % E0" << G4endl;
351  G4cout << " charged traklen: "
352  << std::setw(7) << MeanChargTrLength << " radl +- "
353  << std::setw(7) << rmsChargTrLength << " radl" << G4endl;
354  G4cout << " neutral traklen: "
355  << std::setw(7) << MeanNeutrTrLength << " radl +- "
356  << std::setw(7) << rmsNeutrTrLength << " radl" << G4endl;
357 
358  if (iMoliere > 0. ) {
359  G4double RMoliere1 = iMoliere*fDet->GetdRradl();
360  G4double RMoliere2 = iMoliere*fDet->GetdRlength();
361  G4cout << "\n " << EMoliere << " % confinement: radius = "
362  << RMoliere1 << " radl ("
363  << G4BestUnit( RMoliere2, "Length") << ")" << "\n" << G4endl;
364  }
365 
366  G4cout.setf(mode,std::ios::floatfield);
367  G4cout.precision(prec);
368 
369  // Acceptance
370 
371  G4int nLbin = fDet->GetnLtot();
372  if (limit < DBL_MAX) {
373  EmAcceptance acc;
374  acc.BeginOfAcceptance("Total Energy in Absorber",NbOfEvents);
375  G4double e = MeanELongitCumul[nLbin-1]/100.;
376  G4double r = rmsELongitCumul[nLbin-1]/100.;
377  acc.EmAcceptanceGauss("Edep",NbOfEvents,e,edep,rms,limit);
378  acc.EmAcceptanceGauss("Erms",NbOfEvents,r,rms,rms,2.0*limit);
379  acc.EndOfAcceptance();
380  }
381  limit = DBL_MAX;
382 }
383 
384 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......