ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
RunAction.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file RunAction.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 "RunAction.hh"
34 
35 #include "DetectorConstruction.hh"
36 #include "PrimaryGeneratorAction.hh"
37 #include "HistoManager.hh"
38 
39 #include "G4Run.hh"
40 #include "G4RunManager.hh"
41 #include "G4UnitsTable.hh"
42 #include "G4EmCalculator.hh"
43 
44 #include "Randomize.hh"
45 #include "G4SystemOfUnits.hh"
46 #include <iomanip>
47 
48 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
49 
50 RunAction::RunAction(DetectorConstruction* det, PrimaryGeneratorAction* prim)
51  : G4UserRunAction(),fDetector(det), fPrimary(prim), fProcCounter(0),
52  fHistoManager(0)
53 {
54  fHistoManager = new HistoManager();
55 }
56 
57 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
58 
59 RunAction::~RunAction()
60 {
61  delete fHistoManager;
62 }
63 
64 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
65 
66 void RunAction::BeginOfRunAction(const G4Run*)
67 {
68  // save Rndm status
70  CLHEP::HepRandom::showEngineStatus();
71 
72  fProcCounter = new ProcessesCount;
73  fTotalCount = 0;
74 
75  fTruePL = fTruePL2 = fGeomPL = fGeomPL2 = 0.;
76  fLDispl = fLDispl2 = fPsiSpa = fPsiSpa2 = 0.;
77  fTetPrj = fTetPrj2 = 0.;
78  fPhiCor = fPhiCor2 = 0.;
79 
80  //histograms
81  //
82  G4AnalysisManager* analysisManager = G4AnalysisManager::Instance();
83  if ( analysisManager->IsActive() ) {
84  analysisManager->OpenFile();
85  }
86 }
87 
88 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
89 
90 void RunAction::CountProcesses(G4String procName)
91 {
92  //does the process already encounted ?
93  size_t nbProc = fProcCounter->size();
94  size_t i = 0;
95  while ((i<nbProc)&&((*fProcCounter)[i]->GetName()!=procName)) i++;
96  if (i == nbProc) fProcCounter->push_back( new OneProcessCount(procName));
97 
98  (*fProcCounter)[i]->Count();
99 }
100 
101 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
102 
103 void RunAction::EndOfRunAction(const G4Run* aRun)
104 {
105  G4int NbOfEvents = aRun->GetNumberOfEvent();
106  if (NbOfEvents == 0) return;
107 
108  G4int prec = G4cout.precision(5);
109 
110  G4Material* material = fDetector->GetMaterial();
111  G4double density = material->GetDensity();
112 
113  G4ParticleDefinition* particle =
114  fPrimary->GetParticleGun()->GetParticleDefinition();
115  G4String Particle = particle->GetParticleName();
116  G4double energy = fPrimary->GetParticleGun()->GetParticleEnergy();
117  G4cout << "\n The run consists of " << NbOfEvents << " "<< Particle << " of "
118  << G4BestUnit(energy,"Energy") << " through "
119  << G4BestUnit(fDetector->GetBoxSize(),"Length") << " of "
120  << material->GetName() << " (density: "
121  << G4BestUnit(density,"Volumic Mass") << ")" << G4endl;
122 
123  //frequency of processes
124  G4cout << "\n Process calls frequency --->";
125  for (size_t i=0; i< fProcCounter->size();i++) {
126  G4String procName = (*fProcCounter)[i]->GetName();
127  G4int count = (*fProcCounter)[i]->GetCounter();
128  G4cout << "\t" << procName << " = " << count;
129  }
130 
131  if (fTotalCount > 0) {
132 
133  //compute path length and related quantities
134  //
135  G4double MeanTPL = fTruePL /fTotalCount;
136  G4double MeanTPL2 = fTruePL2/fTotalCount;
137  G4double rmsTPL = std::sqrt(std::fabs(MeanTPL2 - MeanTPL*MeanTPL));
138 
139  G4double MeanGPL = fGeomPL /fTotalCount;
140  G4double MeanGPL2 = fGeomPL2/fTotalCount;
141  G4double rmsGPL = std::sqrt(std::fabs(MeanGPL2 - MeanGPL*MeanGPL));
142 
143  G4double MeanLaD = fLDispl /fTotalCount;
144  G4double MeanLaD2 = fLDispl2/fTotalCount;
145  G4double rmsLaD = std::sqrt(std::fabs(MeanLaD2 - MeanLaD*MeanLaD));
146 
147  G4double MeanPsi = fPsiSpa /(fTotalCount);
148  G4double MeanPsi2 = fPsiSpa2/(fTotalCount);
149  G4double rmsPsi = std::sqrt(std::fabs(MeanPsi2 - MeanPsi*MeanPsi));
150 
151  G4double MeanTeta = fTetPrj /(2*fTotalCount);
152  G4double MeanTeta2 = fTetPrj2/(2*fTotalCount);
153  G4double rmsTeta = std::sqrt(std::fabs(MeanTeta2 - MeanTeta*MeanTeta));
154 
155  G4double MeanCorrel = fPhiCor /(fTotalCount);
156  G4double MeanCorrel2 = fPhiCor2/(fTotalCount);
157  G4double rmsCorrel =
158  std::sqrt(std::fabs(MeanCorrel2-MeanCorrel*MeanCorrel));
159 
160  G4cout << "\n\n truePathLength :\t" << G4BestUnit(MeanTPL,"Length")
161  << " +- " << G4BestUnit( rmsTPL,"Length")
162  << "\n geomPathLength :\t" << G4BestUnit(MeanGPL,"Length")
163  << " +- " << G4BestUnit( rmsGPL,"Length")
164  << "\n lateralDisplac :\t" << G4BestUnit(MeanLaD,"Length")
165  << " +- " << G4BestUnit( rmsLaD,"Length")
166  << "\n Psi :\t" << MeanPsi/mrad << " mrad"
167  << " +- " << rmsPsi /mrad << " mrad"
168  << " (" << MeanPsi/deg << " deg"
169  << " +- " << rmsPsi /deg << " deg)"
170  << G4endl;
171 
172  G4cout << "\n Theta_plane :\t" << rmsTeta/mrad << " mrad"
173  << " (" << rmsTeta/deg << " deg)"
174  << "\n phi correlation:\t" << MeanCorrel
175  << " +- " << rmsCorrel
176  << " (std::cos(phi_pos - phi_dir))"
177  << G4endl;
178 
179 
180  //cross check from G4EmCalculator
181  //
182  G4cout << "\n Verification from G4EmCalculator. \n";
183 
184  G4EmCalculator emCal;
185 
186  //get transport mean free path (for multiple scattering)
187  G4double MSmfp = emCal.GetMeanFreePath(energy,particle,"msc",material);
188 
189  //get range from restricted dedx
190  G4double range = emCal.GetRangeFromRestricteDEDX(energy,particle,material);
191 
192  //effective facRange
193  G4double efFacrange = MeanTPL/std::max(MSmfp, range);
194  if (MeanTPL/range >= 0.99) efFacrange = 1.;
195 
196  G4cout << "\n transport mean free path :\t" << G4BestUnit(MSmfp,"Length")
197  << "\n range from restrict dE/dx:\t" << G4BestUnit(range,"Length")
198  << "\n ---> effective facRange :\t" << efFacrange
199  << G4endl;
200 
201  G4cout << "\n compute theta0 from Highland :\t"
202  << ComputeMscHighland(MeanTPL)/mrad << " mrad"
203  << " (" << ComputeMscHighland(MeanTPL)/deg << " deg)"
204  << G4endl;
205 
206  } else
207  G4cout<< G4endl;
208 
209  //restore default format
210  G4cout.precision(prec);
211 
212  // delete and remove all contents in fProcCounter
213  while (fProcCounter->size()>0){
214  OneProcessCount* aProcCount=fProcCounter->back();
215  fProcCounter->pop_back();
216  delete aProcCount;
217  }
218  delete fProcCounter;
219 
220  //save histograms
221  G4AnalysisManager* analysisManager = G4AnalysisManager::Instance();
222  if ( analysisManager->IsActive() ) {
223  analysisManager->Write();
224  analysisManager->CloseFile();
225  }
226 
227  // show Rndm status
228  CLHEP::HepRandom::showEngineStatus();
229 }
230 
231 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
232 
233 G4double RunAction::ComputeMscHighland(G4double pathLength)
234 {
235  //compute the width of the Gaussian central part of the MultipleScattering
236  //projected angular distribution.
237  //Eur. Phys. Jour. C15 (2000) page 166, formule 23.9
238 
239  G4double t = pathLength/(fDetector->GetMaterial()->GetRadlen());
240  if (t < DBL_MIN) return 0.;
241 
242  G4ParticleGun* particle = fPrimary->GetParticleGun();
243  G4double T = particle->GetParticleEnergy();
244  G4double M = particle->GetParticleDefinition()->GetPDGMass();
245  G4double z = std::abs(particle->GetParticleDefinition()->GetPDGCharge()/eplus);
246 
247  G4double bpc = T*(T+2*M)/(T+M);
248  G4double teta0 = 13.6*MeV*z*std::sqrt(t)*(1.+0.038*std::log(t))/bpc;
249  return teta0;
250 }
251 
252 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......