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

Last change on this file since 812 was 807, checked in by garnier, 17 years ago

update
File size: 9.6 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.2 2006/06/29 16:49:17 gunter Exp $
27// GEANT4 tag $Name: geant4-09-01-patch-02 $
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 totalCrossSection = 0.;
142 for (size_t i=0; i< ProcCounter->size();i++) {
143 G4String procName = (*ProcCounter)[i]->GetName();
144 if (procName != "Transportation")
145 totalCrossSection += ComputeTheory(procName, NbOfEvents);
146 }
147
148 MeanFreePath = 1./totalCrossSection;
149 massCrossSection = totalCrossSection/density;
150
151 G4cout << " Theory: "
152 << "total CrossSection = " << totalCrossSection*cm << " /cm"
153 << "\t MeanFreePath = " << G4BestUnit(MeanFreePath,"Length")
154 << "\t massicCrossSection = " << massCrossSection*g/cm2 << " cm2/g"
155 << G4endl;
156
157
158 G4cout.setf(mode,std::ios::floatfield);
159 G4cout.precision(prec);
160
161 // delete and remove all contents in ProcCounter
162 while (ProcCounter->size()>0){
163 OneProcessCount* aProcCount=ProcCounter->back();
164 ProcCounter->pop_back();
165 delete aProcCount;
166 }
167 delete ProcCounter;
168
169 histoManager->save();
170
171 // show Rndm status
172 CLHEP::HepRandom::showEngineStatus();
173}
174
175//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
176
177G4double RunAction::ComputeTheory(G4String process, G4int NbOfMu)
178{
179 G4Material* material = detector->GetMaterial();
180 G4double length = detector->GetSize();
181 G4double ekin = primary->GetParticleGun()->GetParticleEnergy();
182 MuCrossSections crossSections;
183
184 G4int id = 0; G4double cut = 0.;
185 if (process == "muIoni") {id = 1; cut = GetEnergyCut(material,1); }
186 if (process == "muPairProd") {id = 2; cut = 2*(GetEnergyCut(material,1)
187 + electron_mass_c2); }
188 if (process == "muBrems") {id = 3; cut = GetEnergyCut(material,0); }
189 if (process == "muNucl") id = 4;
190 if (id == 0) return 0.;
191
192 G4int nbOfBins = 100;
193 G4double binMin = -10., binMax = 0., binWidth = (binMax-binMin)/nbOfBins;
194
195#ifdef G4ANALYSIS_USE
196 //create histo for theoritical crossSections, with same bining as simulation
197 //
198 const G4String label[] = { "0", "1", "2", "3", "4", "5", "6", "7", "8", "9",
199 "10", "11", "12", "13", "14", "15", "16", "17", "18", "19"};
200
201 AIDA::IHistogram1D* histoMC = 0; AIDA::IHistogram1D* histoTh = 0;
202 if (histoManager->HistoExist(id)) {
203 histoMC = histoManager->GetHisto(id);
204 nbOfBins = histoManager->GetNbins(id);
205 binMin = histoManager->GetVmin (id);
206 binMax = histoManager->GetVmax (id);
207 binWidth = histoManager->GetBinWidth(id);
208
209 G4String labelTh = label[MaxHisto + id];
210 G4String titleTh = histoManager->GetTitle(id) + " (Th)";
211 histoTh = histoManager->GetHistogramFactory()
212 ->createHistogram1D(labelTh,titleTh,nbOfBins,binMin,binMax);
213 }
214#endif
215
216 //compute and plot differential crossSection, as function of energy transfert.
217 //compute and return integrated crossSection for a given process.
218 //(note: to compare with simulation, the integrated crossSection is function
219 // of the energy cut.)
220 //
221 G4double lgeps, etransf, sigmaE, dsigma, NbProcess;
222 G4double sigmaTot = 0.;
223 const G4double ln10 = std::log(10.);
224
225 for (G4int ibin=0; ibin<nbOfBins; ibin++) {
226 lgeps = binMin + (ibin+0.5)*binWidth;
227 etransf = ekin*std::pow(10.,lgeps);
228 sigmaE = crossSections.CR_Macroscopic(process,material,ekin,etransf);
229 dsigma = sigmaE*etransf*binWidth*ln10;
230 if (etransf > cut) sigmaTot += dsigma;
231 NbProcess = NbOfMu*length*dsigma;
232#ifdef G4ANALYSIS_USE
233 if (histoTh) histoTh->fill(lgeps,NbProcess);
234#endif
235 }
236
237#ifdef G4ANALYSIS_USE
238 //compare simulation and theory
239 //
240 if (histoMC && histoTh) histoManager->GetHistogramFactory()
241 ->divide(label[2*MaxHisto+id], *histoMC, *histoTh);
242#endif
243
244 //return integrated crossSection
245 //
246 return sigmaTot;
247}
248
249//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
250
251#include "G4ProductionCutsTable.hh"
252
253G4double RunAction::GetEnergyCut(G4Material* material, G4int idParticle)
254{
255 G4ProductionCutsTable* table = G4ProductionCutsTable::GetProductionCutsTable();
256
257 size_t index = 0;
258 while ( (table->GetMaterialCutsCouple(index)->GetMaterial() != material) &&
259 (index < table->GetTableSize())) index++;
260
261 return (*(table->GetEnergyCutsVector(idParticle)))[index];
262}
263
264//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
265
Note: See TracBrowser for help on using the repository browser.