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 #include "MuCrossSections.hh"
39 #include "G4ProductionCutsTable.hh"
40 
41 #include "G4Run.hh"
42 #include "G4RunManager.hh"
43 #include "G4UnitsTable.hh"
44 
45 #include "G4PhysicalConstants.hh"
46 #include "G4SystemOfUnits.hh"
47 #include "Randomize.hh"
48 
49 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
50 
51 RunAction::RunAction(DetectorConstruction* det, PrimaryGeneratorAction* prim,
52  HistoManager* HistM)
53  : G4UserRunAction(),
54  fDetector(det), fPrimary(prim), fProcCounter(0), fHistoManager(HistM)
55 {
56  fMucs = new MuCrossSections();
57 }
58 
59 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
60 
61 RunAction::~RunAction()
62 {
63  delete fMucs;
64 }
65 
66 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
67 
68 void RunAction::BeginOfRunAction(const G4Run* aRun)
69 {
70  G4cout << "### Run " << aRun->GetRunID() << " start." << G4endl;
71 
72  // save Rndm status
73  CLHEP::HepRandom::showEngineStatus();
74 
75  fProcCounter = new ProcessesCount();
76  fHistoManager->Book();
77 }
78 
79 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
80 
81 void RunAction::CountProcesses(const G4String& procName)
82 {
83  //does the process already encounted ?
84  size_t n = fProcCounter->size();
85  for(size_t i = 0; i<n; ++i) {
86  if((*fProcCounter)[i]->GetName()==procName) {
87  (*fProcCounter)[i]->Count();
88  return;
89  }
90  }
91  OneProcessCount* count = new OneProcessCount(procName);
92  count->Count();
93  fProcCounter->push_back(count);
94 }
95 
96 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
97 
98 void RunAction::EndOfRunAction(const G4Run* aRun)
99 {
100  G4int NbOfEvents = aRun->GetNumberOfEvent();
101  if (NbOfEvents == 0) return;
102 
103  // std::ios::fmtflags mode = G4cout.flags();
104  G4int prec = G4cout.precision(2);
105 
106  const G4Material* material = fDetector->GetMaterial();
107  G4double length = fDetector->GetSize();
108  G4double density = material->GetDensity();
109 
110  G4String particle = fPrimary->GetParticleGun()->GetParticleDefinition()
111  ->GetParticleName();
112  G4double energy = fPrimary->GetParticleGun()->GetParticleEnergy();
113 
114  G4cout << "\n The run consists of " << NbOfEvents << " "<< particle << " of "
115  << G4BestUnit(energy,"Energy") << " through "
116  << G4BestUnit(length,"Length") << " of "
117  << material->GetName() << " (density: "
118  << G4BestUnit(density,"Volumic Mass") << ")" << G4endl;
119 
120  //total number of process calls
121  G4double countTot = 0.;
122  G4cout << "\n Number of process calls --->";
123  for (size_t i=0; i< fProcCounter->size();++i) {
124  G4String procName = (*fProcCounter)[i]->GetName();
125  if (procName != "Transportation") {
126  G4int count = (*fProcCounter)[i]->GetCounter();
127  G4cout << "\t" << procName << " : " << count;
128  countTot += count;
129  }
130  }
131  G4cout << G4endl;
132 
133  //compute totalCrossSection, meanFreePath and massicCrossSection
134  //
135  G4double totalCrossSection = countTot/(NbOfEvents*length);
136  G4double MeanFreePath = 1./totalCrossSection;
137  G4double massCrossSection = totalCrossSection/density;
138 
139  G4cout.precision(5);
140  G4cout << "\n Simulation: "
141  << "total CrossSection = " << totalCrossSection*cm << " /cm"
142  << "\t MeanFreePath = " << G4BestUnit(MeanFreePath,"Length")
143  << "\t massicCrossSection = " << massCrossSection*g/cm2 << " cm2/g"
144  << G4endl;
145 
146  //compute theoretical predictions
147  //
148  if(particle == "mu+" || particle == "mu-") {
149  totalCrossSection = 0.;
150  for (size_t i=0; i< fProcCounter->size();++i) {
151  G4String procName = (*fProcCounter)[i]->GetName();
152  if (procName != "Transportation") {
153  totalCrossSection += ComputeTheory(procName, NbOfEvents);
154  }
155  }
156 
157  MeanFreePath = 1./totalCrossSection;
158  massCrossSection = totalCrossSection/density;
159 
160  G4cout << " Theory: "
161  << "total CrossSection = " << totalCrossSection*cm << " /cm"
162  << "\t MeanFreePath = " << G4BestUnit(MeanFreePath,"Length")
163  << "\t massicCrossSection = " << massCrossSection*g/cm2 << " cm2/g"
164  << G4endl;
165  }
166 
167  // G4cout.setf(mode,std::ios::floatfield);
168  G4cout.precision(prec);
169 
170  // delete and remove all contents in fProcCounter
171  size_t n = fProcCounter->size();
172  for(size_t i = 0; i<n; ++i) { delete (*fProcCounter)[i]; }
173  delete fProcCounter;
174 
175  fHistoManager->Save();
176 
177  // show Rndm status
178  //CLHEP::HepRandom::showEngineStatus();
179 }
180 
181 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
182 
183 G4double RunAction::ComputeTheory(const G4String& process, G4int NbOfMu)
184 {
185  const G4Material* material = fDetector->GetMaterial();
186  G4double ekin = fPrimary->GetParticleGun()->GetParticleEnergy();
187 
188  G4int id = 0; G4double cut = 1.e-10*ekin;
189  if (process == "muIoni") {id = 11; cut = GetEnergyCut(material,1);}
190  else if (process == "muPairProd") {id = 12; cut = 2*(GetEnergyCut(material,1)
191  + electron_mass_c2); }
192  else if (process == "muBrems") {id = 13; cut = GetEnergyCut(material,0);}
193  else if (process == "muonNuclear"){id = 14; cut = 100*MeV;}
194  if (id == 0) { return 0.; }
195 
196  G4int nbOfBins = 100;
197  //G4double binMin = -10.;
198  G4double binMin = std::log10(cut/ekin);
199  G4double binMax = 0.;
200  G4double binWidth = (binMax-binMin)/G4double(nbOfBins);
201 
202  //create histo for theoretical crossSections, with same bining as simulation
203  //
204  G4AnalysisManager* analysisManager = G4AnalysisManager::Instance();
205 
206  G4H1* histoTh = 0;
207  if (fHistoManager->HistoExist(id)) {
208  histoTh = analysisManager->GetH1(fHistoManager->GetHistoID(id));
209  nbOfBins = fHistoManager->GetNbins(id);
210  binMin = fHistoManager->GetVmin (id);
211  binMax = fHistoManager->GetVmax (id);
212  binWidth = fHistoManager->GetBinWidth(id);
213  }
214 
215  //compute and plot differential crossSection, as function of energy transfert.
216  //compute and return integrated crossSection for a given process.
217  //(note: to compare with simulation, the integrated crossSection is function
218  // of the energy cut.)
219  //
220  G4double lgeps, etransf, sigmaE, dsigma;
221  G4double sigmaTot = 0.;
222  const G4double ln10 = std::log(10.);
223  G4double length = fDetector->GetSize();
224 
225  //G4cout << "MU: " << process << " E= " << ekin
226  // <<" binMin= " << binMin << " binW= " << binWidth << G4endl;
227 
228  for (G4int ibin=0; ibin<nbOfBins; ibin++) {
229  lgeps = binMin + (ibin+0.5)*binWidth;
230  etransf = ekin*std::pow(10.,lgeps);
231  sigmaE = fMucs->CR_Macroscopic(process,material,ekin,etransf);
232  dsigma = sigmaE*etransf*binWidth*ln10;
233  if (etransf > cut) sigmaTot += dsigma;
234  if (histoTh) {
235  G4double NbProcess = NbOfMu*length*dsigma;
236  histoTh->fill(lgeps, NbProcess);
237  }
238  }
239 
240  //return integrated crossSection
241  //
242  return sigmaTot;
243 }
244 
245 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
246 
247 G4double RunAction::GetEnergyCut(const G4Material* material, G4int idParticle)
248 {
249  G4ProductionCutsTable* table = G4ProductionCutsTable::GetProductionCutsTable();
250 
251  size_t index = 0;
252  while ( (table->GetMaterialCutsCouple(index)->GetMaterial() != material) &&
253  (index < table->GetTableSize())) index++;
254 
255  return (*(table->GetEnergyCutsVector(idParticle)))[index];
256 }
257 
258 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
259