source: trunk/examples/extended/polarisation/Pol01/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.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.3 2006/11/17 11:44:46 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
34#include "DetectorConstruction.hh"
35#include "PrimaryGeneratorAction.hh"
36#include "HistoManager.hh"
37
38#include "G4Run.hh"
39#include "G4RunManager.hh"
40#include "G4UnitsTable.hh"
41#include "G4EmCalculator.hh"
42
43#include "Randomize.hh"
44#include <iomanip>
45
46//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
47
48RunAction::RunAction(DetectorConstruction* det, PrimaryGeneratorAction* prim,
49                     HistoManager* histo)
50  : detector(det), primary(prim), ProcCounter(0), histoManager(histo)
51{
52  totalEventCount=0;
53 }
54
55//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
56
57RunAction::~RunAction()
58{ }
59
60//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
61
62void RunAction::BeginOfRunAction(const G4Run* aRun)
63{ 
64  G4cout << "### Run " << aRun->GetRunID() << " start." << G4endl;
65 
66  // save Rndm status
67  //  G4RunManager::GetRunManager()->SetRandomNumberStore(false);
68  //  CLHEP::HepRandom::showEngineStatus();
69 
70  if (ProcCounter) delete ProcCounter;
71  ProcCounter = new ProcessesCount;
72  totalEventCount = 0;
73  photonStats.Clear();
74  electronStats.Clear();
75  positronStats.Clear();
76   
77  histoManager->book();
78}
79
80//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
81void RunAction::FillData(const G4String & particleName,
82                         G4double kinEnergy, G4double costheta, 
83                         G4double /* phi*/,
84                         G4double longitudinalPolarization)
85{
86  if (particleName=="gamma") 
87    photonStats.FillData(kinEnergy, costheta, longitudinalPolarization);
88  else if (particleName=="e-") 
89    electronStats.FillData(kinEnergy, costheta, longitudinalPolarization);
90  else if (particleName=="e+") 
91    positronStats.FillData(kinEnergy, costheta, longitudinalPolarization);
92}
93
94void RunAction::CountProcesses(G4String procName)
95{
96  // is the process already counted ?
97  // *AS* change to std::map?!
98  size_t nbProc = ProcCounter->size();
99  size_t i = 0;
100  while ((i<nbProc)&&((*ProcCounter)[i]->GetName()!=procName)) i++;
101  if (i == nbProc) ProcCounter->push_back( new OneProcessCount(procName));
102 
103  (*ProcCounter)[i]->Count();
104}
105
106//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
107
108void RunAction::EndOfRunAction(const G4Run* aRun)
109{
110  G4int NbOfEvents = aRun->GetNumberOfEvent();
111  if (NbOfEvents == 0) return;
112 
113  G4int  prec = G4cout.precision(5);
114   
115  G4Material* material = detector->GetMaterial();
116  G4double density = material->GetDensity();
117   
118  G4ParticleDefinition* particle = 
119                            primary->GetParticleGun()->GetParticleDefinition();
120  G4String Particle = particle->GetParticleName();   
121  G4double energy = primary->GetParticleGun()->GetParticleEnergy();
122  G4cout << "\n The run consists of " << NbOfEvents << " "<< Particle << " of "
123         << G4BestUnit(energy,"Energy") << " through " 
124         << G4BestUnit(detector->GetBoxSizeZ(),"Length") << " of "
125         << material->GetName() << " (density: " 
126         << G4BestUnit(density,"Volumic Mass") << ")" << G4endl;
127 
128  //frequency of processes
129  G4cout << "\n Process calls frequency --->\n";
130  for (size_t i=0; i< ProcCounter->size();i++) {
131     G4String procName = (*ProcCounter)[i]->GetName();
132     G4int    count    = (*ProcCounter)[i]->GetCounter(); 
133     G4cout << "\t" << procName << " = " << count<<"\n";
134  }
135 
136  if (totalEventCount == 0) return;
137 
138  G4cout<<" Gamma: \n";
139  photonStats.PrintResults(totalEventCount);
140  G4cout<<" Electron: \n";
141  electronStats.PrintResults(totalEventCount);
142  G4cout<<" Positron: \n";
143  positronStats.PrintResults(totalEventCount);
144
145  //cross check from G4EmCalculator
146  //  G4cout << "\n Verification from G4EmCalculator. \n"; 
147  //  G4EmCalculator emCal;
148
149  //restore default format       
150  G4cout.precision(prec);         
151
152  // write out histograms 
153  histoManager->save();
154
155  // show Rndm status
156  CLHEP::HepRandom::showEngineStatus();
157}
158
159//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
160
161void RunAction::EventFinished()
162{
163  ++totalEventCount;
164  photonStats.EventFinished();
165  electronStats.EventFinished();
166  positronStats.EventFinished();
167}
168
169//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
170
171RunAction::ParticleStatistics::ParticleStatistics()
172  : currentNumber(0),
173    totalNumber(0), totalNumber2(0),
174    sumEnergy(0), sumEnergy2(0),
175    sumPolarization(0), sumPolarization2(0),
176    sumCosTheta(0), sumCosTheta2(0)
177{}
178
179//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
180
181RunAction::ParticleStatistics::~ParticleStatistics()
182{}
183
184//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
185
186void RunAction::ParticleStatistics::EventFinished()
187{
188  totalNumber+=currentNumber;
189  totalNumber2+=currentNumber*currentNumber;
190  currentNumber=0;
191}
192//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
193
194void RunAction::ParticleStatistics:: FillData(G4double kinEnergy, 
195                                              G4double costheta,
196                                              G4double longitudinalPolarization)
197{
198  ++currentNumber;
199  sumEnergy+=kinEnergy;
200  sumEnergy2+=kinEnergy*kinEnergy;
201  sumPolarization+=longitudinalPolarization;
202  sumPolarization2+=longitudinalPolarization*longitudinalPolarization;
203  sumCosTheta+=costheta;
204  sumCosTheta2+=costheta*costheta;
205}
206
207//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
208
209void RunAction::ParticleStatistics::PrintResults(G4int totalNumberOfEvents)
210{
211  G4cout<<"Mean Number per Event :"
212        <<G4double(totalNumber)/G4double(totalNumberOfEvents)<<"\n";
213  if (totalNumber==0) totalNumber=1;
214  G4double energyMean=sumEnergy/totalNumber;
215  G4double energyRms=std::sqrt(sumEnergy2/totalNumber-energyMean*energyMean);
216  G4cout<<"Mean Energy :"<< G4BestUnit(energyMean,"Energy")
217        <<" +- "<<G4BestUnit(energyRms,"Energy")<<"\n";
218  G4double polarizationMean=sumPolarization/totalNumber;
219  G4double polarizationRms=
220    std::sqrt(sumPolarization2/totalNumber-polarizationMean*polarizationMean);
221  G4cout<<"Mean Polarization :"<< polarizationMean
222        <<" +- "<<polarizationRms<<"\n";
223}
224
225//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
226
227void RunAction::ParticleStatistics::Clear()
228{
229  currentNumber=0;
230  totalNumber=totalNumber2=0;
231  sumEnergy=sumEnergy2=0;
232  sumPolarization=sumPolarization2=0;
233  sumCosTheta=sumCosTheta2=0;
234}
235
236//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
Note: See TracBrowser for help on using the repository browser.