source: trunk/examples/extended/electromagnetic/TestEm5/src/RunAction.cc @ 1333

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

update to geant4.9.3
File size: 11.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.29 2009/06/18 19:08:18 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#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* aRun)
60{
61  G4cout << "### Run " << aRun->GetRunID() << " start." << G4endl;
62
63  //initialisation
64  EnergyDeposit  = EnergyDeposit2  = 0.;
65  TrakLenCharged = TrakLenCharged2 = 0.;
66  TrakLenNeutral = TrakLenNeutral2 = 0.;
67  nbStepsCharged = nbStepsCharged2 = 0.;
68  nbStepsNeutral = nbStepsNeutral2 = 0.;
69  MscProjecTheta = MscProjecTheta2 = 0.;
70  MscThetaCentral = 3*ComputeMscHighland();
71
72  nbGamma = nbElect = nbPosit = 0;
73
74  Transmit[0] = Transmit[1] = Reflect[0] = Reflect[1] = 0;
75 
76  MscEntryCentral = 0;
77 
78  EnergyLeak[0] = EnergyLeak[1] = EnergyLeak2[0] = EnergyLeak2[1] = 0.;
79
80  histoManager->book();
81
82  // save Rndm status
83  G4RunManager::GetRunManager()->SetRandomNumberStore(true);
84  CLHEP::HepRandom::showEngineStatus();
85}
86
87//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
88
89void RunAction::EndOfRunAction(const G4Run* aRun)
90{
91  // compute mean and rms
92  //
93  G4int TotNbofEvents = aRun->GetNumberOfEvent();
94  if (TotNbofEvents == 0) return;
95 
96  G4double EnergyBalance = EnergyDeposit + EnergyLeak[0] + EnergyLeak[1];
97  EnergyBalance /= TotNbofEvents;
98
99  EnergyDeposit /= TotNbofEvents; EnergyDeposit2 /= TotNbofEvents;
100  G4double rmsEdep = EnergyDeposit2 - EnergyDeposit*EnergyDeposit;
101  if (rmsEdep>0.) rmsEdep = std::sqrt(rmsEdep/TotNbofEvents);
102  else            rmsEdep = 0.;
103
104  TrakLenCharged /= TotNbofEvents; TrakLenCharged2 /= TotNbofEvents;
105  G4double rmsTLCh = TrakLenCharged2 - TrakLenCharged*TrakLenCharged;
106  if (rmsTLCh>0.) rmsTLCh = std::sqrt(rmsTLCh/TotNbofEvents);
107  else            rmsTLCh = 0.;
108
109  TrakLenNeutral /= TotNbofEvents; TrakLenNeutral2 /= TotNbofEvents;
110  G4double rmsTLNe = TrakLenNeutral2 - TrakLenNeutral*TrakLenNeutral;
111  if (rmsTLNe>0.) rmsTLNe = std::sqrt(rmsTLNe/TotNbofEvents);
112  else            rmsTLNe = 0.;
113
114  nbStepsCharged /= TotNbofEvents; nbStepsCharged2 /= TotNbofEvents;
115  G4double rmsStCh = nbStepsCharged2 - nbStepsCharged*nbStepsCharged;
116  if (rmsStCh>0.) rmsStCh = std::sqrt(rmsTLCh/TotNbofEvents);
117  else            rmsStCh = 0.;
118
119  nbStepsNeutral /= TotNbofEvents; nbStepsNeutral2 /= TotNbofEvents;
120  G4double rmsStNe = nbStepsNeutral2 - nbStepsNeutral*nbStepsNeutral;
121  if (rmsStNe>0.) rmsStNe = std::sqrt(rmsTLCh/TotNbofEvents);
122  else            rmsStNe = 0.;
123
124  G4double Gamma = (G4double)nbGamma/TotNbofEvents;
125  G4double Elect = (G4double)nbElect/TotNbofEvents;
126  G4double Posit = (G4double)nbPosit/TotNbofEvents;
127
128  G4double transmit[2];
129  transmit[0] = 100.*Transmit[0]/TotNbofEvents;
130  transmit[1] = 100.*Transmit[1]/TotNbofEvents;
131
132  G4double reflect[2];
133  reflect[0] = 100.*Reflect[0]/TotNbofEvents;
134  reflect[1] = 100.*Reflect[1]/TotNbofEvents;
135
136  G4double rmsMsc = 0., tailMsc = 0.;
137  if (MscEntryCentral > 0) {
138    MscProjecTheta /= MscEntryCentral; MscProjecTheta2 /= MscEntryCentral;
139    rmsMsc = MscProjecTheta2 - MscProjecTheta*MscProjecTheta;
140    if (rmsMsc > 0.) rmsMsc = std::sqrt(rmsMsc);
141    tailMsc = 100.- (100.*MscEntryCentral)/(2*Transmit[1]);   
142  }
143 
144  EnergyLeak[0] /= TotNbofEvents; EnergyLeak2[0] /= TotNbofEvents;
145  G4double rmsEl0 = EnergyLeak2[0] - EnergyLeak[0]*EnergyLeak[0];
146  if (rmsEl0>0.) rmsEl0 = std::sqrt(rmsEl0/TotNbofEvents);
147  else           rmsEl0 = 0.;
148 
149  EnergyLeak[1] /= TotNbofEvents; EnergyLeak2[1] /= TotNbofEvents;
150  G4double rmsEl1 = EnergyLeak2[1] - EnergyLeak[1]*EnergyLeak[1];
151  if (rmsEl1>0.) rmsEl1 = std::sqrt(rmsEl1/TotNbofEvents);
152  else           rmsEl1 = 0.;   
153 
154     
155  //Stopping Power from input Table.
156  //
157  G4Material* material = detector->GetAbsorberMaterial();
158  G4double length  = detector->GetAbsorberThickness();
159  G4double density = material->GetDensity();
160   
161  G4ParticleDefinition* particle = primary->GetParticleGun()
162                                          ->GetParticleDefinition();
163  G4String partName = particle->GetParticleName();
164  G4double energy = primary->GetParticleGun()->GetParticleEnergy();
165
166  G4EmCalculator emCalculator;
167  G4double dEdxTable = 0., dEdxFull = 0.;
168  if (particle->GetPDGCharge()!= 0.) { 
169    dEdxTable = emCalculator.GetDEDX(energy,particle,material);
170    dEdxFull  = emCalculator.ComputeTotalDEDX(energy,particle,material);   
171  }
172  G4double stopTable = dEdxTable/density;
173  G4double stopFull  = dEdxFull /density; 
174   
175  //Stopping Power from simulation.
176  //   
177  G4double meandEdx  = EnergyDeposit/length;
178  G4double stopPower = meandEdx/density; 
179
180  G4cout << "\n ======================== run summary ======================\n";
181
182  G4int prec = G4cout.precision(3);
183 
184  G4cout << "\n The run was " << TotNbofEvents << " " << partName << " of "
185         << G4BestUnit(energy,"Energy") << " through " 
186         << G4BestUnit(length,"Length") << " of "
187         << material->GetName() << " (density: " 
188         << G4BestUnit(density,"Volumic Mass") << ")" << G4endl;
189 
190  G4cout.precision(4);
191 
192  G4cout << "\n Total energy deposit in absorber per event = "
193         << G4BestUnit(EnergyDeposit,"Energy") << " +- "
194         << G4BestUnit(rmsEdep,      "Energy") 
195         << G4endl;
196         
197  G4cout << "\n -----> Mean dE/dx = " << meandEdx/(MeV/cm) << " MeV/cm"
198         << "\t(" << stopPower/(MeV*cm2/g) << " MeV*cm2/g)"
199         << G4endl;
200         
201  G4cout << "\n From formulas :" << G4endl; 
202  G4cout << "   restricted dEdx = " << dEdxTable/(MeV/cm) << " MeV/cm"
203         << "\t(" << stopTable/(MeV*cm2/g) << " MeV*cm2/g)"
204         << G4endl;
205         
206  G4cout << "   full dEdx       = " << dEdxFull/(MeV/cm) << " MeV/cm"
207         << "\t(" << stopFull/(MeV*cm2/g) << " MeV*cm2/g)"
208         << G4endl;
209         
210  G4cout << "\n Leakage :  primary = "
211         << G4BestUnit(EnergyLeak[0],"Energy") << " +- "
212         << G4BestUnit(rmsEl0,       "Energy")
213         << "   secondaries = "
214         << G4BestUnit(EnergyLeak[1],"Energy") << " +- "
215         << G4BestUnit(rmsEl1,       "Energy")   
216         << G4endl;
217         
218  G4cout << " Energy balance :  edep + eleak = "
219         << G4BestUnit(EnergyBalance,"Energy")
220         << G4endl;     
221                         
222  G4cout << "\n Total track length (charged) in absorber per event = "
223         << G4BestUnit(TrakLenCharged,"Length") << " +- "
224         << G4BestUnit(rmsTLCh,       "Length") << G4endl;
225
226  G4cout << " Total track length (neutral) in absorber per event = "
227         << G4BestUnit(TrakLenNeutral,"Length") << " +- "
228         << G4BestUnit(rmsTLNe,       "Length") << G4endl;
229
230  G4cout << "\n Number of steps (charged) in absorber per event = "
231         << nbStepsCharged << " +- " << rmsStCh << G4endl;
232
233  G4cout << " Number of steps (neutral) in absorber per event = "
234         << nbStepsNeutral << " +- " << rmsStNe << G4endl;
235
236  G4cout << "\n Number of secondaries per event : Gammas = " << Gamma
237         << ";   electrons = " << Elect
238         << ";   positrons = " << Posit << G4endl;
239
240  G4cout << "\n Number of events with the primary particle transmitted = "
241         << transmit[1] << " %" << G4endl;
242
243  G4cout << " Number of events with at least  1 particle transmitted "
244         << "(same charge as primary) = " << transmit[0] << " %" << G4endl;
245
246  G4cout << "\n Number of events with the primary particle reflected = "
247         << reflect[1] << " %" << G4endl;
248
249  G4cout << " Number of events with at least  1 particle reflected "
250         << "(same charge as primary) = " << reflect[0] << " %" << G4endl;
251
252  // compute width of the Gaussian central part of the MultipleScattering
253  //
254  G4cout << "\n MultipleScattering:" 
255         << "\n  rms proj angle of transmit primary particle = "
256         << rmsMsc/mrad << " mrad (central part only)" << G4endl;
257
258  G4cout << "  computed theta0 (Highland formula)          = "
259         << ComputeMscHighland()/mrad << " mrad" << G4endl;
260           
261  G4cout << "  central part defined as +- "
262         << MscThetaCentral/mrad << " mrad; " 
263         << "  Tail ratio = " << tailMsc << " %" << G4endl;       
264
265  G4cout.precision(prec);
266 
267  // normalize histograms
268  //
269  G4int ih = 1;
270  G4double binWidth = histoManager->GetBinWidth(ih);
271  G4double unit     = histoManager->GetHistoUnit(ih); 
272  G4double fac = unit/(TotNbofEvents*binWidth);
273  histoManager->Scale(ih,fac);
274   
275  ih = 10;
276  binWidth = histoManager->GetBinWidth(ih);
277  unit     = histoManager->GetHistoUnit(ih); 
278  fac = unit/(TotNbofEvents*binWidth);
279  histoManager->Scale(ih,fac);
280   
281  // save histograms
282  histoManager->save();
283
284  // show Rndm status
285  CLHEP::HepRandom::showEngineStatus();
286}
287
288//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
289
290G4double RunAction::ComputeMscHighland()
291{
292  //compute the width of the Gaussian central part of the MultipleScattering
293  //projected angular distribution.
294  //Eur. Phys. Jour. C15 (2000) page 166, formule 23.9
295
296  G4double t = (detector->GetAbsorberThickness())
297              /(detector->GetAbsorberMaterial()->GetRadlen());
298  if (t < DBL_MIN) return 0.;
299
300  G4ParticleGun* particle = primary->GetParticleGun();
301  G4double T = particle->GetParticleEnergy();
302  G4double M = particle->GetParticleDefinition()->GetPDGMass();
303  G4double z = std::abs(particle->GetParticleDefinition()->GetPDGCharge()/eplus);
304
305  G4double bpc = T*(T+2*M)/(T+M);
306  G4double teta0 = 13.6*MeV*z*std::sqrt(t)*(1.+0.038*std::log(t))/bpc;
307  return teta0;
308}
309
310//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
Note: See TracBrowser for help on using the repository browser.