source: trunk/examples/extended/electromagnetic/TestEm11/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: 8.4 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.7 2007/08/19 20:52:53 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 "PhysicsList.hh"
35#include "StepMax.hh"
36#include "PrimaryGeneratorAction.hh"
37#include "HistoManager.hh"
38
39#include "G4Run.hh"
40#include "G4RunManager.hh"
41#include "G4UnitsTable.hh"
42#include "G4EmCalculator.hh"
43
44#include "Randomize.hh"
45#include <iomanip>
46
47//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
48
49RunAction::RunAction(DetectorConstruction* det, PhysicsList* phys,
50                     PrimaryGeneratorAction* kin, HistoManager* histo)
51:detector(det),physics(phys),kinematic(kin),histoManager(histo)
52{ }
53
54//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
55
56RunAction::~RunAction()
57{ }
58
59//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
60
61void RunAction::BeginOfRunAction(const G4Run* aRun)
62{ 
63  G4cout << "### Run " << aRun->GetRunID() << " start." << G4endl;
64 
65  // save Rndm status
66  G4RunManager::GetRunManager()->SetRandomNumberStore(true);
67  CLHEP::HepRandom::showEngineStatus();
68 
69  //initialize total energy deposit
70  //
71  Edeposit = Edeposit2 = 0.; 
72   
73  //initialize track legth of primary
74  //
75  trackLen = trackLen2 = 0.;
76   
77  //initialize projected range
78  //
79  projRange = projRange2 = 0.;
80   
81  //initialize mean step size
82  //
83  nbOfSteps = nbOfSteps2 = 0;  stepSize = stepSize2 = 0.;
84 
85  //initialize track status
86  //
87  status[0] = status[1] = status[2] = 0;
88 
89  //get csdaRange from EmCalculator
90  //
91  G4EmCalculator emCalculator;
92  G4Material* material = detector->GetAbsorMaterial();
93  G4ParticleDefinition* particle = kinematic->GetParticleGun()
94                                          ->GetParticleDefinition();
95  G4double energy = kinematic->GetParticleGun()->GetParticleEnergy();
96  csdaRange = DBL_MAX;
97  if (particle->GetPDGCharge() != 0.)
98    csdaRange = emCalculator.GetCSDARange(energy,particle,material);
99       
100  //histograms
101  //
102  histoManager->book();
103   
104  //set StepMax from histos
105  //
106  G4double stepMax = histoManager->ComputeStepMax(csdaRange);
107  physics->GetStepMaxProcess()->SetMaxStep(stepMax);           
108}
109
110//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
111
112void RunAction::EndOfRunAction(const G4Run* aRun)
113{
114  std::ios::fmtflags mode = G4cout.flags();
115  G4cout.setf(std::ios::fixed,std::ios::floatfield);
116 
117  G4int NbofEvents = aRun->GetNumberOfEvent();
118  if (NbofEvents == 0) return;
119  G4double fNbofEvents = double(NbofEvents);
120
121  //run conditions
122  // 
123  G4Material* material = detector->GetAbsorMaterial();
124  G4double density = material->GetDensity();
125   
126  G4ParticleDefinition* particle = kinematic->GetParticleGun()
127                                          ->GetParticleDefinition();
128  G4String partName = particle->GetParticleName();                       
129  G4double energy = kinematic->GetParticleGun()->GetParticleEnergy();
130 
131  G4cout << "\n ======================== run summary ======================\n";
132 
133  G4int prec = G4cout.precision(2);
134 
135  G4cout << "\n The run consists of " << NbofEvents << " "<< partName << " of "
136         << G4BestUnit(energy,"Energy") << " through " 
137         << G4BestUnit(detector->GetAbsorSizeX(),"Length") << " of "
138         << material->GetName() << " (density: " 
139         << G4BestUnit(density,"Volumic Mass") << ")" << G4endl;
140         
141  G4cout << "\n ============================================================\n";
142 
143  //compute total energy deposit
144  //
145  Edeposit /= NbofEvents; Edeposit2 /= NbofEvents;
146  G4double rms = Edeposit2 - Edeposit*Edeposit;       
147  if (rms>0.) rms = std::sqrt(rms); else rms = 0.;
148
149  G4cout.precision(3);       
150  G4cout
151    << "\n Total Energy deposited        = " << G4BestUnit(Edeposit,"Energy")
152    << " +- "                                << G4BestUnit( rms,"Energy")
153    << G4endl;
154         
155  //compute track length of primary track
156  //
157  trackLen /= NbofEvents; trackLen2 /= NbofEvents;
158  rms = trackLen2 - trackLen*trackLen;       
159  if (rms>0.) rms = std::sqrt(rms); else rms = 0.;
160
161  G4cout.precision(3);       
162  G4cout
163    << "\n Track length of primary track = " << G4BestUnit(trackLen,"Length")
164    << " +- "                                << G4BestUnit( rms,"Length");
165   
166  //compare with csda range
167  //
168  //G4EmCalculator emCalculator;
169  //G4double csdaRange = 0.;
170  //if (particle->GetPDGCharge() != 0.)
171  //  csdaRange = emCalculator.GetCSDARange(energy,particle,material);
172  G4cout
173    << "\n Range from EmCalculator       = " << G4BestUnit(csdaRange,"Length")
174    << " (from full dE/dx)" << G4endl;
175                 
176  //compute projected range of primary track
177  //
178  projRange /= NbofEvents; projRange2 /= NbofEvents;
179  rms = projRange2 - projRange*projRange;       
180  if (rms>0.) rms = std::sqrt(rms); else rms = 0.;
181   
182  G4cout
183    << "\n Projected range               = " << G4BestUnit(projRange,"Length")
184    << " +- "                                << G4BestUnit( rms,"Length")   
185    << G4endl;
186   
187  //nb of steps and step size of primary track
188  //
189  G4double fNbSteps = nbOfSteps/fNbofEvents, fNbSteps2 = nbOfSteps2/fNbofEvents;
190  rms = fNbSteps2 - fNbSteps*fNbSteps;       
191  if (rms>0.) rms = std::sqrt(rms); else rms = 0.;
192
193  G4cout.precision(2);       
194  G4cout << "\n Nb of steps of primary track  = " << fNbSteps << " +- " << rms;
195   
196  stepSize /= NbofEvents; stepSize2 /= NbofEvents;
197  rms = stepSize2 - stepSize*stepSize;       
198  if (rms>0.) rms = std::sqrt(rms); else rms = 0.;
199
200  G4cout.precision(3);       
201  G4cout
202    << "\t Step size= " << G4BestUnit(stepSize,"Length")
203    << " +- "           << G4BestUnit( rms,"Length")
204    << G4endl;
205   
206  //transmission coefficients
207  //
208  G4double absorbed  = 100.*status[0]/fNbofEvents;
209  G4double transmit  = 100.*status[1]/fNbofEvents;
210  G4double reflected = 100.*status[2]/fNbofEvents; 
211
212  G4cout.precision(2);       
213  G4cout
214    << "\n absorbed = "  << absorbed  << " %"
215    << "   transmit = "  << transmit  << " %"
216    << "   reflected = " << reflected << " %" << G4endl;
217     
218  // normalize histograms of longitudinal energy profile
219  //
220  G4int ih = 1;
221  G4double binWidth = histoManager->GetBinWidth(ih);
222  G4double fac = (1./(NbofEvents*binWidth))*(mm/MeV);
223  histoManager->Scale(ih,fac);
224 
225  ih = 8;
226  binWidth = histoManager->GetBinWidth(ih);
227  G4double range = histoManager->GetcsdaRange();
228  fac = (1./(NbofEvents*binWidth*range*density))*(g/(MeV*cm2));
229  histoManager->Scale(ih,fac);
230   
231   // reset default formats
232  G4cout.setf(mode,std::ios::floatfield);
233  G4cout.precision(prec);
234 
235  // save histograms
236  histoManager->save();
237 
238  // show Rndm status
239  CLHEP::HepRandom::showEngineStatus();
240}
241
242//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
Note: See TracBrowser for help on using the repository browser.