source: trunk/examples/extended/electromagnetic/TestEm17/src/RunAction.cc @ 1333

Last change on this file since 1333 was 1230, checked in by garnier, 15 years ago

update to geant4.9.3
File size: 9.9 KB
Line 
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//
26// $Id: RunAction.cc,v 1.4 2008/03/31 10:22:59 vnivanch Exp $
27// GEANT4 tag $Name: geant4-09-03-cand-01 $
28//
29//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
30//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
31
32#include "RunAction.hh"
33
34#include "DetectorConstruction.hh"
35#include "PrimaryGeneratorAction.hh"
36#include "HistoManager.hh"
37#include "MuCrossSections.hh"
38
39#include "G4Run.hh"
40#include "G4RunManager.hh"
41#include "G4UnitsTable.hh"
42
43#include "Randomize.hh"
44
45#ifdef G4ANALYSIS_USE
46 #include "AIDA/AIDA.h"
47#endif
48
49//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
50
51RunAction::RunAction(DetectorConstruction* det, PrimaryGeneratorAction* prim,
52                     HistoManager* HistM)
53  : detector(det), primary(prim), ProcCounter(0), histoManager(HistM)
54{}
55
56//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
57
58RunAction::~RunAction()
59{}
60
61//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
62
63void RunAction::BeginOfRunAction(const G4Run* aRun)
64{ 
65  G4cout << "### Run " << aRun->GetRunID() << " start." << G4endl;
66 
67  // save Rndm status
68  G4RunManager::GetRunManager()->SetRandomNumberStore(true);
69  CLHEP::HepRandom::showEngineStatus();
70
71  ProcCounter = new ProcessesCount;
72 
73  histoManager->book();
74}
75
76//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
77
78void RunAction::CountProcesses(G4String procName)
79{
80   //does the process  already encounted ?
81   size_t nbProc = ProcCounter->size();
82   size_t i = 0;
83   while ((i<nbProc)&&((*ProcCounter)[i]->GetName()!=procName)) i++;
84   if (i == nbProc) ProcCounter->push_back( new OneProcessCount(procName));
85
86   (*ProcCounter)[i]->Count();
87}
88
89//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
90
91void 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 = detector->GetMaterial();
100  G4double length  = detector->GetSize();
101  G4double density = material->GetDensity();
102   
103  G4String particle = primary->GetParticleGun()->GetParticleDefinition()
104                      ->GetParticleName();   
105  G4double energy = primary->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< ProcCounter->size();i++) {
117     G4String procName = (*ProcCounter)[i]->GetName();
118     if (procName != "Transportation") {
119       G4int count    = (*ProcCounter)[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< ProcCounter->size();i++) {
144      G4String procName = (*ProcCounter)[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 ProcCounter
163  while (ProcCounter->size()>0){
164    OneProcessCount* aProcCount=ProcCounter->back();
165    ProcCounter->pop_back();
166    delete aProcCount;
167  }
168  delete ProcCounter;
169 
170  histoManager->save();
171 
172  // show Rndm status
173  CLHEP::HepRandom::showEngineStatus();
174}
175
176//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
177
178G4double RunAction::ComputeTheory(G4String process, G4int NbOfMu)
179{   
180  G4Material* material = detector->GetMaterial();
181  G4double length = detector->GetSize();
182  G4double ekin = primary->GetParticleGun()->GetParticleEnergy();
183  MuCrossSections crossSections;
184
185  G4int id = 0; G4double cut = 0.;
186  if (process == "muIoni")          {id = 1; cut =    GetEnergyCut(material,1);}
187  else if (process == "muPairProd") {id = 2; cut = 2*(GetEnergyCut(material,1) 
188                                                      + electron_mass_c2); }
189  else if (process == "muBrems")    {id = 3; cut =    GetEnergyCut(material,0);}
190  else if (process == "muNucl")      id = 4;
191  else if (process == "hIoni")      {id = 5; cut =    GetEnergyCut(material,1);}
192  else if (process == "hPairProd")  {id = 6; cut = 2*(GetEnergyCut(material,1) 
193                                                      + electron_mass_c2); }
194  else if (process == "hBrems")     {id = 7; cut =    GetEnergyCut(material,0);}
195  if (id == 0) return 0.;
196 
197  G4int nbOfBins = 100;
198  G4double binMin = -10., binMax = 0., binWidth = (binMax-binMin)/nbOfBins;
199   
200#ifdef G4ANALYSIS_USE
201  //create histo for theoritical crossSections, with same bining as simulation
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  AIDA::IHistogram1D* histoMC = 0; AIDA::IHistogram1D* histoTh = 0; 
207  if (histoManager->HistoExist(id)) {
208    histoMC  = histoManager->GetHisto(id); 
209    nbOfBins = histoManager->GetNbins(id);
210    binMin   = histoManager->GetVmin (id);
211    binMax   = histoManager->GetVmax (id);
212    binWidth = histoManager->GetBinWidth(id);
213   
214    G4String labelTh = label[MaxHisto + id];
215    G4String titleTh = histoManager->GetTitle(id) + " (Th)";
216    histoTh = histoManager->GetHistogramFactory()
217          ->createHistogram1D(labelTh,titleTh,nbOfBins,binMin,binMax);   
218  }
219#endif
220 
221  //compute and plot differential crossSection, as function of energy transfert.
222  //compute and return integrated crossSection for a given process.
223  //(note: to compare with simulation, the integrated crossSection is function
224  //       of the energy cut.)
225  //
226  G4double lgeps, etransf, sigmaE, dsigma, NbProcess;
227  G4double sigmaTot = 0.;
228  const G4double ln10 = std::log(10.);
229   
230  for (G4int ibin=0; ibin<nbOfBins; ibin++) {
231    lgeps = binMin + (ibin+0.5)*binWidth;
232    etransf = ekin*std::pow(10.,lgeps);
233    sigmaE = crossSections.CR_Macroscopic(process,material,ekin,etransf);
234    dsigma = sigmaE*etransf*binWidth*ln10;
235    if (etransf > cut) sigmaTot += dsigma;   
236    NbProcess = NbOfMu*length*dsigma;
237#ifdef G4ANALYSIS_USE
238    if (histoTh) histoTh->fill(lgeps,NbProcess);
239#endif     
240  }
241 
242#ifdef G4ANALYSIS_USE
243  //compare simulation and theory
244  //
245  if (histoMC && histoTh) histoManager->GetHistogramFactory()
246                     ->divide(label[2*MaxHisto+id], *histoMC, *histoTh);
247#endif
248   
249  //return integrated crossSection
250  //
251  return sigmaTot;   
252}
253
254//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
255
256#include "G4ProductionCutsTable.hh"
257
258G4double RunAction::GetEnergyCut(G4Material* material, G4int idParticle)
259{ 
260 G4ProductionCutsTable* table = G4ProductionCutsTable::GetProductionCutsTable();
261 
262 size_t index = 0;
263 while ( (table->GetMaterialCutsCouple(index)->GetMaterial() != material) &&
264        (index < table->GetTableSize())) index++;
265
266 return (*(table->GetEnergyCutsVector(idParticle)))[index];
267} 
268
269//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
270               
Note: See TracBrowser for help on using the repository browser.