source: trunk/examples/extended/electromagnetic/TestEm18/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.2 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 2007/02/16 11:59:47 maire 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#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), 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* run)
60{
61 G4cout << "### Run " << run->GetRunID() << " start." << G4endl;
62
63 //initialisation
64 //
65 energyDeposit = 0.;
66
67 nbCharged = nbNeutral = 0;
68 energyCharged = energyNeutral = 0.;
69 emin[0] = emin[1] = DBL_MAX;
70 emax[0] = emax[1] = 0.;
71
72 nbSteps = 0;
73 trackLength = 0.;
74
75 histoManager->book();
76
77 // do not save Rndm status
78 G4RunManager::GetRunManager()->SetRandomNumberStore(false);
79 CLHEP::HepRandom::showEngineStatus();
80}
81
82//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
83
84void RunAction::EndOfRunAction(const G4Run* aRun)
85{
86 G4int nbEvents = aRun->GetNumberOfEvent();
87 if (nbEvents == 0) return;
88
89 G4Material* material = detector->GetMaterial();
90 G4double length = detector->GetSize();
91 G4double density = material->GetDensity();
92
93 G4ParticleDefinition* particle = primary->GetParticleGun()
94 ->GetParticleDefinition();
95 G4String partName = particle->GetParticleName();
96 G4double eprimary = primary->GetParticleGun()->GetParticleEnergy();
97
98 G4int prec = G4cout.precision(3);
99 G4cout << "\n ======================== run summary ======================\n";
100 G4cout << "\n The run was " << nbEvents << " " << partName << " of "
101 << G4BestUnit(eprimary,"Energy") << " through "
102 << G4BestUnit(length,"Length") << " of "
103 << material->GetName() << " (density: "
104 << G4BestUnit(density,"Volumic Mass") << ")";
105 G4cout << "\n ===========================================================\n";
106 G4cout << G4endl;
107
108 if (particle->GetPDGCharge() == 0.) return;
109
110 G4cout.precision(5);
111
112 //track length
113 //
114 G4double trackLPerEvent = trackLength/nbEvents;
115 G4double nbStepPerEvent = double(nbSteps)/nbEvents;
116 G4double stepSize = trackLength/nbSteps;
117
118 G4cout
119 << "\n trackLength= "
120 << G4BestUnit(trackLPerEvent, "Length")
121 << "\t nb of steps= " << nbStepPerEvent
122 << " stepSize= " << G4BestUnit(stepSize, "Length")
123 << G4endl;
124
125 //charged secondaries (ionization, direct pair production)
126 //
127 G4double energyPerEvent = energyCharged/nbEvents;
128 G4double nbPerEvent = double(nbCharged)/nbEvents;
129 G4double meanEkin = 0.;
130 if (nbCharged) meanEkin = energyCharged/nbCharged;
131
132 G4cout
133 << "\n d-rays : eLoss/primary= "
134 << G4BestUnit(energyPerEvent, "Energy")
135 << "\t nb of d-rays= " << nbPerEvent
136 << " <Tkin>= " << G4BestUnit(meanEkin, "Energy")
137 << " Tmin= " << G4BestUnit(emin[0], "Energy")
138 << " Tmax= " << G4BestUnit(emax[0], "Energy")
139 << G4endl;
140
141 //neutral secondaries (bremsstrahlung)
142 //
143 energyPerEvent = energyNeutral/nbEvents;
144 nbPerEvent = double(nbNeutral)/nbEvents;
145 meanEkin = 0.;
146 if (nbNeutral) meanEkin = energyNeutral/nbNeutral;
147
148 G4cout
149 << "\n brems : eLoss/primary= "
150 << G4BestUnit(energyPerEvent, "Energy")
151 << "\t nb of gammas= " << nbPerEvent
152 << " <Tkin>= " << G4BestUnit(meanEkin, "Energy")
153 << " Tmin= " << G4BestUnit(emin[1], "Energy")
154 << " Tmax= " << G4BestUnit(emax[1], "Energy")
155 << G4endl;
156
157
158 G4EmCalculator emCal;
159
160 //local energy deposit
161 //
162 energyPerEvent = energyDeposit/nbEvents;
163 //
164 G4double r0 = emCal.GetRangeFromRestricteDEDX(eprimary,particle,material);
165 G4double r1 = r0 - trackLPerEvent;
166 G4double etry = eprimary - energyPerEvent;
167 G4double efinal = 0.;
168 if (r1 > 0.) efinal = GetEnergyFromRestrictedRange(r1,particle,material,etry);
169 G4double dEtable = eprimary - efinal;
170 G4double ratio = 0.;
171 if (dEtable > 0.) ratio = energyPerEvent/dEtable;
172
173 G4cout
174 << "\n deposit : eLoss/primary= "
175 << G4BestUnit(energyPerEvent, "Energy")
176 << "\t <dEcut > table= "
177 << G4BestUnit(dEtable, "Energy")
178 << " ---> simul/reference= " << ratio
179 << G4endl;
180
181 //total energy transferred
182 //
183 G4double energyTotal = energyDeposit + energyCharged + energyNeutral;
184 energyPerEvent = energyTotal/nbEvents;
185 //
186 r0 = emCal.GetCSDARange(eprimary,particle,material);
187 r1 = r0 - trackLPerEvent;
188 etry = eprimary - energyPerEvent;
189 efinal = 0.;
190 if (r1 > 0.) efinal = GetEnergyFromCSDARange(r1,particle,material,etry);
191 dEtable = eprimary - efinal;
192 ratio = 0.;
193 if (dEtable > 0.) ratio = energyPerEvent/dEtable;
194
195 G4cout
196 << "\n total : eLoss/primary= "
197 << G4BestUnit(energyPerEvent, "Energy")
198 << "\t <dEfull> table= "
199 << G4BestUnit(dEtable, "Energy")
200 << " ---> simul/reference= " << ratio
201 << G4endl;
202
203 G4cout.precision(prec);
204
205 histoManager->save();
206
207 // show Rndm status
208 CLHEP::HepRandom::showEngineStatus();
209}
210
211//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
212
213G4double RunAction::GetEnergyFromRestrictedRange(G4double range,
214 G4ParticleDefinition* particle, G4Material* material, G4double Etry)
215{
216 G4EmCalculator emCal;
217
218 G4double Energy = Etry, dE = 0., dEdx;
219 G4double r, dr;
220 G4double err = 1., errmax = 0.00001;
221 G4int iter = 0 , itermax = 10;
222 while (err > errmax && iter < itermax) {
223 iter++;
224 Energy -= dE;
225 r = emCal.GetRangeFromRestricteDEDX(Energy,particle,material);
226 dr = r - range;
227 dEdx = emCal.GetDEDX(Energy,particle,material);
228 dE = dEdx*dr;
229 err = std::abs(dE)/Energy;
230 }
231 if (iter == itermax) {
232 G4cout
233 << "\n ---> warning: RunAction::GetEnergyFromRestRange() did not converge"
234 << " Etry = " << G4BestUnit(Etry,"Energy")
235 << " Energy = " << G4BestUnit(Energy,"Energy")
236 << " err = " << err
237 << " iter = " << iter << G4endl;
238 }
239
240 return Energy;
241}
242
243//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
244
245G4double RunAction::GetEnergyFromCSDARange(G4double range,
246 G4ParticleDefinition* particle, G4Material* material, G4double Etry)
247{
248 G4EmCalculator emCal;
249
250 G4double Energy = Etry, dE = 0., dEdx;
251 G4double r, dr;
252 G4double err = 1., errmax = 0.00001;
253 G4int iter = 0 , itermax = 10;
254 while (err > errmax && iter < itermax) {
255 iter++;
256 Energy -= dE;
257 r = emCal.GetCSDARange(Energy,particle,material);
258 dr = r - range;
259 dEdx = emCal.ComputeTotalDEDX(Energy,particle,material);
260 dE = dEdx*dr;
261 err = std::abs(dE)/Energy;
262 }
263 if (iter == itermax) {
264 G4cout
265 << "\n ---> warning: RunAction::GetEnergyFromCSDARange() did not converge"
266 << " Etry = " << G4BestUnit(Etry,"Energy")
267 << " Energy = " << G4BestUnit(Energy,"Energy")
268 << " err = " << err
269 << " iter = " << iter << G4endl;
270 }
271
272 return Energy;
273}
274
275//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
Note: See TracBrowser for help on using the repository browser.