Geant4  9.6.p02
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
RunAction.cc
Go to the documentation of this file.
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 //
28 //
29 // $Id$
30 //
31 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
32 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
33 
34 #include "RunAction.hh"
35 
36 #include "DetectorConstruction.hh"
37 #include "PrimaryGeneratorAction.hh"
38 #include "HistoManager.hh"
39 
40 #include "G4Run.hh"
41 #include "G4RunManager.hh"
42 #include "G4HadronicProcessStore.hh"
43 #include "G4UnitsTable.hh"
44 #include "G4SystemOfUnits.hh"
45 
46 #include "Randomize.hh"
47 #include <iomanip>
48 
49 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
50 
51 RunAction::RunAction(DetectorConstruction* det, PrimaryGeneratorAction* prim)
52  : fDetector(det), fPrimary(prim)
53 {
54  fHistoManager = new HistoManager();
55 }
56 
57 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
58 
59 RunAction::~RunAction()
60 {
61  delete fHistoManager;
62 }
63 
64 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
65 
66 void RunAction::BeginOfRunAction(const G4Run* aRun)
67 {
68  G4cout << "### Run " << aRun->GetRunID() << " start." << G4endl;
69 
70  // save Rndm status
71  G4RunManager::GetRunManager()->SetRandomNumberStore(false);
72  CLHEP::HepRandom::showEngineStatus();
73 
74  fTotalCount = fGammaCount = 0;
75  fSumTrack = fSumTrack2 = 0.;
76  for (G4int i=0; i<3; i++) { fPbalance[i] = 0. ; }
77  for (G4int i=0; i<3; i++) { fNbGamma[i] = 0 ; }
78 
79  //histograms
80  //
81  G4AnalysisManager* analysisManager = G4AnalysisManager::Instance();
82  if ( analysisManager->IsActive() ) {
83  analysisManager->OpenFile();
84  }
85 }
86 
87 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
88 
89 void RunAction::ParticleCount(G4String name, G4double Ekin)
90 {
91  fParticleCount[name]++;
92  fEmean[name] += Ekin;
93  //update min max
94  if (fParticleCount[name] == 1) fEmin[name] = fEmax[name] = Ekin;
95  if (Ekin < fEmin[name]) fEmin[name] = Ekin;
96  if (Ekin > fEmax[name]) fEmax[name] = Ekin;
97 }
98 
99 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
100 
101 void RunAction::CountNuclearChannel(G4String name, G4double Q)
102 {
103  fNuclChannelCount[name]++;
104  fNuclChannelQ[name] += Q;
105 }
106 
107 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
108 
109 void RunAction::Balance(G4double Pbal)
110 {
111  fPbalance[0] += Pbal;
112  //update min max
113  if (fTotalCount == 1) fPbalance[1] = fPbalance[2] = Pbal;
114  if (Pbal < fPbalance[1]) fPbalance[1] = Pbal;
115  if (Pbal > fPbalance[2]) fPbalance[2] = Pbal;
116 }
117 
118 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
119 
120 void RunAction::CountGamma(G4int nGamma)
121 {
122  fGammaCount++;
123  fNbGamma[0] += nGamma;
124  //update min max
125  if (fGammaCount == 1) fNbGamma[1] = fNbGamma[2] = nGamma;
126  if (nGamma < fNbGamma[1]) fNbGamma[1] = nGamma;
127  if (nGamma > fNbGamma[2]) fNbGamma[2] = nGamma;
128 }
129 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
130 
131 void RunAction::EndOfRunAction(const G4Run* aRun)
132 {
133  G4int NbOfEvents = aRun->GetNumberOfEvent();
134  if (NbOfEvents == 0) return;
135 
136  G4int prec = 5, wid = prec + 2;
137  G4int dfprec = G4cout.precision(prec);
138 
139  G4Material* material = fDetector->GetMaterial();
140  G4double density = material->GetDensity();
141  G4int survive = 0;
142 
143  G4ParticleDefinition* particle =
144  fPrimary->GetParticleGun()->GetParticleDefinition();
145  G4String Particle = particle->GetParticleName();
146  G4double energy = fPrimary->GetParticleGun()->GetParticleEnergy();
147  G4cout << "\n The run consists of " << NbOfEvents << " "<< Particle << " of "
148  << G4BestUnit(energy,"Energy") << " through "
149  << G4BestUnit(fDetector->GetSize(),"Length") << " of "
150  << material->GetName() << " (density: "
151  << G4BestUnit(density,"Volumic Mass") << ")" << G4endl;
152 
153  //frequency of processes
154  G4cout << "\n Process calls frequency --->";
155  std::map<const G4VProcess*,G4int>::iterator it;
156  for (it = fProcCounter.begin(); it != fProcCounter.end(); it++) {
157  G4String procName = it->first->GetProcessName();
158  G4int count = it->second;
159  G4cout << "\t" << procName << "= " << count;
160  if (procName == "Transportation") survive = count;
161  }
162 
163  if (survive > 0) {
164  G4cout << "\n\n Nb of incident particles surviving after "
165  << G4BestUnit(fDetector->GetSize(),"Length") << " of "
166  << material->GetName() << " : " << survive << G4endl;
167  }
168 
169  if (fTotalCount == 0) fTotalCount = 1; //force printing anyway
170 
171  //compute mean free path and related quantities
172  //
173  G4double MeanFreePath = fSumTrack /fTotalCount;
174  G4double MeanTrack2 = fSumTrack2/fTotalCount;
175  G4double rms = std::sqrt(std::fabs(MeanTrack2 - MeanFreePath*MeanFreePath));
176  G4double CrossSection = 0.0;
177  if(MeanFreePath > 0.0) { CrossSection = 1./MeanFreePath; }
178  G4double massicMFP = MeanFreePath*density;
179  G4double massicCS = 0.0;
180  if(massicMFP > 0.0) { massicCS = 1./massicMFP; }
181 
182  G4cout << "\n\n MeanFreePath:\t" << G4BestUnit(MeanFreePath,"Length")
183  << " +- " << G4BestUnit( rms,"Length")
184  << "\tmassic: " << G4BestUnit(massicMFP, "Mass/Surface")
185  << "\n CrossSection:\t" << CrossSection*cm << " cm^-1 "
186  << "\t\tmassic: " << G4BestUnit(massicCS, "Surface/Mass")
187  << G4endl;
188 
189  //check cross section from G4HadronicProcessStore
190  //
191  G4cout << "\n Verification : "
192  << "crossSections from G4HadronicProcessStore:";
193 
194  G4HadronicProcessStore* store = G4HadronicProcessStore::Instance();
195  G4double sumc = 0.0;
196  for (it = fProcCounter.begin(); it != fProcCounter.end(); it++) {
197  const G4VProcess* process = it->first;
198  G4double xs =
199  store->GetCrossSectionPerVolume(particle,energy,process,material);
200  G4double massSigma = xs/density;
201  sumc += massSigma;
202  G4String procName = process->GetProcessName();
203  G4cout << "\n " << procName << "= "
204  << G4BestUnit(massSigma, "Surface/Mass");
205  }
206  G4cout << "\n total = "
207  << G4BestUnit(sumc, "Surface/Mass") << G4endl;
208 
209  //nuclear channel count
210  //
211  G4cout << "\n List of nuclear reactions: \n" << G4endl;
212 
213  std::map<G4String,G4int>::iterator ic;
214  for (ic = fNuclChannelCount.begin(); ic != fNuclChannelCount.end(); ic++) {
215  G4String name = ic->first;
216  G4int count = ic->second;
217  G4double Q = fNuclChannelQ[name]/count;
218 
219  G4cout << " " << std::setw(50) << name << ": " << std::setw(7) << count
220  << " Q = " << std::setw(wid) << G4BestUnit(Q, "Energy")
221  << G4endl;
222  }
223 
224  //Gamma count
225  //
226  if (fGammaCount > 0) {
227  G4cout << "\n" << std::setw(58) << "Number of gamma: N = "
228  << fNbGamma[1] << " --> " << fNbGamma[2] << G4endl;
229  }
230 
231  //particles count
232  //
233  G4cout << "\n List of generated particles: \n" << G4endl;
234 
235  std::map<G4String,G4int>::iterator ip;
236  for (ip = fParticleCount.begin(); ip != fParticleCount.end(); ip++) {
237  G4String name = ip->first;
238  G4int count = ip->second;
239  G4double eMean = fEmean[name]/count;
240  G4double eMin = fEmin[name], eMax = fEmax[name];
241 
242  G4cout << " " << std::setw(13) << name << ": " << std::setw(7) << count
243  << " Emean = " << std::setw(wid) << G4BestUnit(eMean, "Energy")
244  << "\t( " << G4BestUnit(eMin, "Energy")
245  << " --> " << G4BestUnit(eMax, "Energy")
246  << ")" << G4endl;
247  }
248 
249  //energy momentum balance
250  //
251  if (fTotalCount > 1) {
252  G4double Pbmean = fPbalance[0]/fTotalCount;
253  G4cout << "\n Momentum balance: Pmean = "
254  << std::setw(wid) << G4BestUnit(Pbmean, "Energy")
255  << "\t( " << G4BestUnit(fPbalance[1], "Energy")
256  << " --> " << G4BestUnit(fPbalance[2], "Energy")
257  << ") \n" << G4endl;
258  }
259 
260  //restore default format
261  G4cout.precision(dfprec);
262 
263  // remove all contents in fProcCounter
264  fProcCounter.clear();
265  // remove all contents in fNuclChannel
266  fNuclChannelCount.clear();
267  fNuclChannelQ.clear();
268  // remove all contents in fParticleCount
269  fParticleCount.clear();
270  fEmean.clear(); fEmin.clear(); fEmax.clear();
271 
272  //save histograms
273  G4AnalysisManager* analysisManager = G4AnalysisManager::Instance();
274  if ( analysisManager->IsActive() ) {
275  analysisManager->Write();
276  analysisManager->CloseFile();
277  }
278 
279  // show Rndm status
280  CLHEP::HepRandom::showEngineStatus();
281 }
282 
283 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
Generated on Sat May 25 2013 14:31:59 for Geant4 by   doxygen 1.8.4