source: trunk/examples/extended/electromagnetic/TestEm1/src/RunAction.cc@ 1036

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

update
File size: 7.7 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.19 2006/06/29 16:37:25 gunter Exp $
27// GEANT4 tag $Name: $
28//
29//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
30//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
31
32#include "RunAction.hh"
33#include "DetectorConstruction.hh"
34#include "PrimaryGeneratorAction.hh"
35#include "HistoManager.hh"
36
37#include "G4Run.hh"
38#include "G4RunManager.hh"
39#include "G4UnitsTable.hh"
40#include "G4EmCalculator.hh"
41
42#include "Randomize.hh"
43#include <iomanip>
44
45//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
46
47RunAction::RunAction(DetectorConstruction* det, PrimaryGeneratorAction* kin,
48 HistoManager* histo)
49:detector(det), primary(kin), ProcCounter(0), histoManager(histo)
50{ }
51
52//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
53
54RunAction::~RunAction()
55{ }
56
57//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
58
59void RunAction::BeginOfRunAction(const G4Run* aRun)
60{
61 G4cout << "### Run " << aRun->GetRunID() << " start." << G4endl;
62
63 // save Rndm status
64 G4RunManager::GetRunManager()->SetRandomNumberStore(true);
65 CLHEP::HepRandom::showEngineStatus();
66
67 NbOfTraks0 = NbOfTraks1 = NbOfSteps0 = NbOfSteps1 = 0;
68 edep = 0.;
69 trueRange = trueRange2 = 0.;
70 projRange = projRange2 = 0.;
71 transvDev = transvDev2 = 0.;
72 ProcCounter = new ProcessesCount;
73
74 //histograms
75 //
76 histoManager->book();
77}
78
79//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
80
81void RunAction::CountProcesses(G4String procName)
82{
83 //does the process already encounted ?
84 size_t nbProc = ProcCounter->size();
85 size_t i = 0;
86 while ((i<nbProc)&&((*ProcCounter)[i]->GetName()!=procName)) i++;
87 if (i == nbProc) ProcCounter->push_back( new OneProcessCount(procName));
88
89 (*ProcCounter)[i]->Count();
90}
91
92//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
93
94void RunAction::EndOfRunAction(const G4Run* aRun)
95{
96 G4int NbOfEvents = aRun->GetNumberOfEvent();
97 if (NbOfEvents == 0) return;
98 G4double dNbOfEvents = double(NbOfEvents);
99
100 G4ParticleDefinition* particle = primary->GetParticleGun()
101 ->GetParticleDefinition();
102 G4String partName = particle->GetParticleName();
103 G4double energy = primary->GetParticleGun()->GetParticleEnergy();
104
105 G4double length = detector->GetSize();
106 G4Material* material = detector->GetMaterial();
107 G4double density = material->GetDensity();
108
109 G4cout << "\n ======================== run summary ======================\n";
110
111 G4int prec = G4cout.precision(5);
112
113 G4cout << "\n The run was: " << NbOfEvents << " " << partName << " of "
114 << G4BestUnit(energy,"Energy") << " through "
115 << G4BestUnit(length,"Length") << " of "
116 << material->GetName() << " (density: "
117 << G4BestUnit(density,"Volumic Mass") << ")" << G4endl;
118
119 G4cout << "\n ============================================================\n";
120
121 G4cout << "\n total energy deposit: "
122 << G4BestUnit(edep/dNbOfEvents, "Energy") << G4endl;
123
124 //nb of tracks and steps per event
125 //
126 G4cout << "\n nb tracks/event"
127 << " neutral: " << std::setw(10) << NbOfTraks0/dNbOfEvents
128 << " charged: " << std::setw(10) << NbOfTraks1/dNbOfEvents
129 << "\n nb steps/event"
130 << " neutral: " << std::setw(10) << NbOfSteps0/dNbOfEvents
131 << " charged: " << std::setw(10) << NbOfSteps1/dNbOfEvents
132 << G4endl;
133
134 //frequency of processes call
135 G4cout << "\n nb of process calls per event: \n ";
136 for (size_t i=0; i< ProcCounter->size();i++)
137 G4cout << std::setw(12) << (*ProcCounter)[i]->GetName();
138
139 G4cout << "\n ";
140 for (size_t j=0; j< ProcCounter->size();j++)
141 G4cout << std::setw(12) << ((*ProcCounter)[j]->GetCounter())/dNbOfEvents;
142 G4cout << G4endl;
143
144 //compute true and projected ranges, and transverse dispersion
145 //
146 trueRange /= NbOfEvents; trueRange2 /= NbOfEvents;
147 G4double trueRms = trueRange2 - trueRange*trueRange;
148 if (trueRms>0.) trueRms = std::sqrt(trueRms); else trueRms = 0.;
149
150 projRange /= NbOfEvents; projRange2 /= NbOfEvents;
151 G4double projRms = projRange2 - projRange*projRange;
152 if (projRms>0.) projRms = std::sqrt(projRms); else projRms = 0.;
153
154 transvDev /= 2*NbOfEvents; transvDev2 /= 2*NbOfEvents;
155 G4double trvsRms = transvDev2 - transvDev*transvDev;
156 if (trvsRms>0.) trvsRms = std::sqrt(trvsRms); else trvsRms = 0.;
157
158 //compare true range with csda range from PhysicsTables
159 //
160 G4EmCalculator emCalculator;
161 G4double rangeTable = 0.;
162 if (particle->GetPDGCharge() != 0.)
163 rangeTable = emCalculator.GetCSDARange(energy,particle,material);
164
165 G4cout << "\n---------------------------------------------------------\n";
166 G4cout << " Primary particle : " ;
167 G4cout << "\n true Range = " << G4BestUnit(trueRange,"Length")
168 << " rms = " << G4BestUnit(trueRms, "Length");
169
170 G4cout << "\n proj Range = " << G4BestUnit(projRange,"Length")
171 << " rms = " << G4BestUnit(projRms, "Length");
172
173 G4cout << "\n proj/true = " << projRange/trueRange;
174
175 G4cout << "\n transverse dispersion at end = "
176 << G4BestUnit(trvsRms,"Length");
177
178 G4cout << "\n mass true Range from simulation = "
179 << G4BestUnit(trueRange*density, "Mass/Surface")
180 << "\n from PhysicsTable (csda range) = "
181 << G4BestUnit(rangeTable*density, "Mass/Surface");
182 G4cout << "\n---------------------------------------------------------\n";
183 G4cout << G4endl;
184
185 // reset default precision
186 G4cout.precision(prec);
187
188 // delete and remove all contents in ProcCounter
189 while (ProcCounter->size()>0){
190 OneProcessCount* aProcCount=ProcCounter->back();
191 ProcCounter->pop_back();
192 delete aProcCount;
193 }
194 delete ProcCounter;
195
196 //save histograms
197 histoManager->save();
198
199 // show Rndm status
200 CLHEP::HepRandom::showEngineStatus();
201}
202
203//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
Note: See TracBrowser for help on using the repository browser.