source: trunk/examples/extended/electromagnetic/TestEm18/src/RunAction.cc @ 1337

Last change on this file since 1337 was 1337, checked in by garnier, 14 years ago

tag geant4.9.4 beta 1 + modifs locales
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-04-beta-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.