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

Last change on this file since 1230 was 1230, checked in by garnier, 16 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.