Geant4  9.6.p02
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
RunAction.cc
Go to the documentation of this file.
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 // $Id$
30 //
31 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
32 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
33 
34 #include "RunAction.hh"
35 
36 #include "DetectorConstruction.hh"
37 #include "PrimaryGeneratorAction.hh"
38 #include "HistoManager.hh"
39 #include "MuCrossSections.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  : fDetector(det), fPrimary(prim), fProcCounter(0), fHistoManager(HistM)
54 {}
55 
56 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
57 
58 RunAction::~RunAction()
59 {}
60 
61 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
62 
63 void RunAction::BeginOfRunAction(const G4Run* aRun)
64 {
65  G4cout << "### Run " << aRun->GetRunID() << " start." << G4endl;
66 
67  // save Rndm status
68  G4RunManager::GetRunManager()->SetRandomNumberStore(false);
69  CLHEP::HepRandom::showEngineStatus();
70 
71  fProcCounter = new ProcessesCount;
72 
73  fHistoManager->book();
74 }
75 
76 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
77 
78 void RunAction::CountProcesses(G4String procName)
79 {
80  //does the process already encounted ?
81  size_t nbProc = fProcCounter->size();
82  size_t i = 0;
83  while ((i<nbProc)&&((*fProcCounter)[i]->GetName()!=procName)) i++;
84  if (i == nbProc) fProcCounter->push_back( new OneProcessCount(procName));
85 
86  (*fProcCounter)[i]->Count();
87 }
88 
89 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
90 
91 void RunAction::EndOfRunAction(const G4Run* aRun)
92 {
93  G4int NbOfEvents = aRun->GetNumberOfEvent();
94  if (NbOfEvents == 0) return;
95 
96  std::ios::fmtflags mode = G4cout.flags();
97  G4int prec = G4cout.precision(2);
98 
99  G4Material* material = fDetector->GetMaterial();
100  G4double length = fDetector->GetSize();
101  G4double density = material->GetDensity();
102 
103  G4String particle = fPrimary->GetParticleGun()->GetParticleDefinition()
104  ->GetParticleName();
105  G4double energy = fPrimary->GetParticleGun()->GetParticleEnergy();
106 
107  G4cout << "\n The run consists of " << NbOfEvents << " "<< particle << " of "
108  << G4BestUnit(energy,"Energy") << " through "
109  << G4BestUnit(length,"Length") << " of "
110  << material->GetName() << " (density: "
111  << G4BestUnit(density,"Volumic Mass") << ")" << G4endl;
112 
113  //total number of process calls
114  G4double countTot = 0.;
115  G4cout << "\n Number of process calls --->";
116  for (size_t i=0; i< fProcCounter->size();i++) {
117  G4String procName = (*fProcCounter)[i]->GetName();
118  if (procName != "Transportation") {
119  G4int count = (*fProcCounter)[i]->GetCounter();
120  G4cout << "\t" << procName << " : " << count;
121  countTot += count;
122  }
123  }
124  G4cout << G4endl;
125 
126  //compute totalCrossSection, meanFreePath and massicCrossSection
127  //
128  G4double totalCrossSection = countTot/(NbOfEvents*length);
129  G4double MeanFreePath = 1./totalCrossSection;
130  G4double massCrossSection =totalCrossSection/density;
131 
132  G4cout.precision(5);
133  G4cout << "\n Simulation: "
134  << "total CrossSection = " << totalCrossSection*cm << " /cm"
135  << "\t MeanFreePath = " << G4BestUnit(MeanFreePath,"Length")
136  << "\t massicCrossSection = " << massCrossSection*g/cm2 << " cm2/g"
137  << G4endl;
138 
139  //compute theoritical predictions
140  //
141  if(particle == "mu+" || particle == "mu-") {
142  totalCrossSection = 0.;
143  for (size_t i=0; i< fProcCounter->size();i++) {
144  G4String procName = (*fProcCounter)[i]->GetName();
145  if (procName != "Transportation")
146  totalCrossSection += ComputeTheory(procName, NbOfEvents);
147  }
148 
149  MeanFreePath = 1./totalCrossSection;
150  massCrossSection = totalCrossSection/density;
151 
152  G4cout << " Theory: "
153  << "total CrossSection = " << totalCrossSection*cm << " /cm"
154  << "\t MeanFreePath = " << G4BestUnit(MeanFreePath,"Length")
155  << "\t massicCrossSection = " << massCrossSection*g/cm2 << " cm2/g"
156  << G4endl;
157  }
158 
159  G4cout.setf(mode,std::ios::floatfield);
160  G4cout.precision(prec);
161 
162  // delete and remove all contents in fProcCounter
163  while (fProcCounter->size()>0){
164  OneProcessCount* aProcCount=fProcCounter->back();
165  fProcCounter->pop_back();
166  delete aProcCount;
167  }
168  delete fProcCounter;
169 
170  fHistoManager->save();
171 
172  // show Rndm status
173  CLHEP::HepRandom::showEngineStatus();
174 }
175 
176 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
177 
178 G4double RunAction::ComputeTheory(G4String process, G4int NbOfMu)
179 {
180  G4Material* material = fDetector->GetMaterial();
181  G4double ekin = fPrimary->GetParticleGun()->GetParticleEnergy();
182  MuCrossSections crossSections;
183 
184  G4int id = 0; G4double cut = 0.;
185  if (process == "muIoni") {id = 1; cut = GetEnergyCut(material,1);}
186  else if (process == "muPairProd") {id = 2; cut = 2*(GetEnergyCut(material,1)
187  + electron_mass_c2); }
188  else if (process == "muBrems") {id = 3; cut = GetEnergyCut(material,0);}
189  else if (process == "muNucl") id = 4;
190  else if (process == "hIoni") {id = 5; cut = GetEnergyCut(material,1);}
191  else if (process == "hPairProd") {id = 6; cut = 2*(GetEnergyCut(material,1)
192  + electron_mass_c2); }
193  else if (process == "hBrems") {id = 7; cut = GetEnergyCut(material,0);}
194  if (id == 0) return 0.;
195 
196  G4int nbOfBins = 100;
197  G4double binMin = -10., binMax = 0., binWidth = (binMax-binMin)/nbOfBins;
198 
199  //create histo for theoritical crossSections, with same bining as simulation
200  //
201  G4AnalysisManager* analysisManager = G4AnalysisManager::Instance();
202 
203  const G4String label[] = { "0", "1", "2", "3", "4", "5", "6", "7", "8", "9",
204  "10", "11", "12", "13", "14", "15", "16", "17", "18", "19"};
205 
206  G4AnaH1* histoTh = 0;
208  if (fHistoManager->HistoExist(id)) {
210  nbOfBins = fHistoManager->GetNbins(id);
211  binMin = fHistoManager->GetVmin (id);
212  binMax = fHistoManager->GetVmax (id);
213  binWidth = fHistoManager->GetBinWidth(id);
214 
215  G4String labelTh = label[MaxHisto + id];
216  G4String titleTh = fHistoManager->GetTitle(id) + " (Th)";
217  G4int histThId = analysisManager
218  ->CreateH1(labelTh,titleTh,nbOfBins,binMin,binMax);
219  histoTh = analysisManager->GetH1(histThId);
220  }
221 
222  //compute and plot differential crossSection, as function of energy transfert.
223  //compute and return integrated crossSection for a given process.
224  //(note: to compare with simulation, the integrated crossSection is function
225  // of the energy cut.)
226  //
227  G4double lgeps, etransf, sigmaE, dsigma;
228  G4double sigmaTot = 0.;
229  const G4double ln10 = std::log(10.);
230  G4double length = fDetector->GetSize();
231 
232  for (G4int ibin=0; ibin<nbOfBins; ibin++) {
233  lgeps = binMin + (ibin+0.5)*binWidth;
234  etransf = ekin*std::pow(10.,lgeps);
235  sigmaE = crossSections.CR_Macroscopic(process,material,ekin,etransf);
236  dsigma = sigmaE*etransf*binWidth*ln10;
237  if (etransf > cut) sigmaTot += dsigma;
238  G4double NbProcess = NbOfMu*length*dsigma;
239  if (histoTh) histoTh->fill(lgeps, NbProcess);
240  }
241 
242  //compare simulation and theory
243  //
246 
247  //return integrated crossSection
248  //
249  return sigmaTot;
250 }
251 
252 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
253 
254 #include "G4ProductionCutsTable.hh"
255 
256 G4double RunAction::GetEnergyCut(G4Material* material, G4int idParticle)
257 {
258  G4ProductionCutsTable* table = G4ProductionCutsTable::GetProductionCutsTable();
259 
260  size_t index = 0;
261  while ( (table->GetMaterialCutsCouple(index)->GetMaterial() != material) &&
262  (index < table->GetTableSize())) index++;
263 
264  return (*(table->GetEnergyCutsVector(idParticle)))[index];
265 }
266 
267 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
268 
Generated on Sat May 25 2013 14:31:59 for Geant4 by   doxygen 1.8.4