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 #include "DetectorConstruction.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 "G4SystemOfUnits.hh"
46 #include <iomanip>
47 
48 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
49 
50 RunAction::RunAction(DetectorConstruction* det, PrimaryGeneratorAction* kin)
51 :fDetector(det), fPrimary(kin), fHistoManager(0)
52 {
53  // Book predefined histograms
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  //initialisation
71  fEnergyDeposit = fEnergyDeposit2 = 0.;
72  fTrakLenCharged = fTrakLenCharged2 = 0.;
73  fTrakLenNeutral = fTrakLenNeutral2 = 0.;
74  fNbStepsCharged = fNbStepsCharged2 = 0.;
75  fNbStepsNeutral = fNbStepsNeutral2 = 0.;
76  fMscProjecTheta = fMscProjecTheta2 = 0.;
77  fMscThetaCentral = 3*ComputeMscHighland();
78 
79  fNbGamma = fNbElect = fNbPosit = 0;
80 
81  fTransmit[0] = fTransmit[1] = fReflect[0] = fReflect[1] = 0;
82 
83  fMscEntryCentral = 0;
84 
85  fEnergyLeak[0] = fEnergyLeak[1] = fEnergyLeak2[0] = fEnergyLeak2[1] = 0.;
86 
87  G4AnalysisManager* analysisManager = G4AnalysisManager::Instance();
88  if ( analysisManager->IsActive() ) {
89  analysisManager->OpenFile();
90  }
91 
92  // save Rndm status
94  CLHEP::HepRandom::showEngineStatus();
95 }
96 
97 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
98 
99 void RunAction::EndOfRunAction(const G4Run* aRun)
100 {
101  // compute mean and rms
102  //
103  G4int TotNbofEvents = aRun->GetNumberOfEvent();
104  if (TotNbofEvents == 0) return;
105 
106  G4double EnergyBalance = fEnergyDeposit + fEnergyLeak[0] + fEnergyLeak[1];
107  EnergyBalance /= TotNbofEvents;
108 
109  fEnergyDeposit /= TotNbofEvents; fEnergyDeposit2 /= TotNbofEvents;
110  G4double rmsEdep = fEnergyDeposit2 - fEnergyDeposit*fEnergyDeposit;
111  if (rmsEdep>0.) rmsEdep = std::sqrt(rmsEdep/TotNbofEvents);
112  else rmsEdep = 0.;
113 
114  fTrakLenCharged /= TotNbofEvents; fTrakLenCharged2 /= TotNbofEvents;
115  G4double rmsTLCh = fTrakLenCharged2 - fTrakLenCharged*fTrakLenCharged;
116  if (rmsTLCh>0.) rmsTLCh = std::sqrt(rmsTLCh/TotNbofEvents);
117  else rmsTLCh = 0.;
118 
119  fTrakLenNeutral /= TotNbofEvents; fTrakLenNeutral2 /= TotNbofEvents;
120  G4double rmsTLNe = fTrakLenNeutral2 - fTrakLenNeutral*fTrakLenNeutral;
121  if (rmsTLNe>0.) rmsTLNe = std::sqrt(rmsTLNe/TotNbofEvents);
122  else rmsTLNe = 0.;
123 
124  fNbStepsCharged /= TotNbofEvents; fNbStepsCharged2 /= TotNbofEvents;
125  G4double rmsStCh = fNbStepsCharged2 - fNbStepsCharged*fNbStepsCharged;
126  if (rmsStCh>0.) rmsStCh = std::sqrt(rmsTLCh/TotNbofEvents);
127  else rmsStCh = 0.;
128 
129  fNbStepsNeutral /= TotNbofEvents; fNbStepsNeutral2 /= TotNbofEvents;
130  G4double rmsStNe = fNbStepsNeutral2 - fNbStepsNeutral*fNbStepsNeutral;
131  if (rmsStNe>0.) rmsStNe = std::sqrt(rmsTLCh/TotNbofEvents);
132  else rmsStNe = 0.;
133 
134  G4double Gamma = (G4double)fNbGamma/TotNbofEvents;
135  G4double Elect = (G4double)fNbElect/TotNbofEvents;
136  G4double Posit = (G4double)fNbPosit/TotNbofEvents;
137 
138  G4double transmit[2];
139  transmit[0] = 100.*fTransmit[0]/TotNbofEvents;
140  transmit[1] = 100.*fTransmit[1]/TotNbofEvents;
141 
142  G4double reflect[2];
143  reflect[0] = 100.*fReflect[0]/TotNbofEvents;
144  reflect[1] = 100.*fReflect[1]/TotNbofEvents;
145 
146  G4double rmsMsc = 0., tailMsc = 0.;
147  if (fMscEntryCentral > 0) {
148  fMscProjecTheta /= fMscEntryCentral; fMscProjecTheta2 /= fMscEntryCentral;
149  rmsMsc = fMscProjecTheta2 - fMscProjecTheta*fMscProjecTheta;
150  if (rmsMsc > 0.) { rmsMsc = std::sqrt(rmsMsc); }
151  if(fTransmit[1] > 0.0) {
152  tailMsc = 100.- (100.*fMscEntryCentral)/(2*fTransmit[1]);
153  }
154  }
155 
156  fEnergyLeak[0] /= TotNbofEvents; fEnergyLeak2[0] /= TotNbofEvents;
157  G4double rmsEl0 = fEnergyLeak2[0] - fEnergyLeak[0]*fEnergyLeak[0];
158  if (rmsEl0>0.) rmsEl0 = std::sqrt(rmsEl0/TotNbofEvents);
159  else rmsEl0 = 0.;
160 
161  fEnergyLeak[1] /= TotNbofEvents; fEnergyLeak2[1] /= TotNbofEvents;
162  G4double rmsEl1 = fEnergyLeak2[1] - fEnergyLeak[1]*fEnergyLeak[1];
163  if (rmsEl1>0.) rmsEl1 = std::sqrt(rmsEl1/TotNbofEvents);
164  else rmsEl1 = 0.;
165 
166 
167  //Stopping Power from input Table.
168  //
169  G4Material* material = fDetector->GetAbsorberMaterial();
170  G4double length = fDetector->GetAbsorberThickness();
171  G4double density = material->GetDensity();
172 
173  G4ParticleDefinition* particle = fPrimary->GetParticleGun()
174  ->GetParticleDefinition();
175  G4String partName = particle->GetParticleName();
176  G4double energy = fPrimary->GetParticleGun()->GetParticleEnergy();
177 
178  G4EmCalculator emCalculator;
179  G4double dEdxTable = 0., dEdxFull = 0.;
180  if (particle->GetPDGCharge()!= 0.) {
181  dEdxTable = emCalculator.GetDEDX(energy,particle,material);
182  dEdxFull = emCalculator.ComputeTotalDEDX(energy,particle,material);
183  }
184  G4double stopTable = dEdxTable/density;
185  G4double stopFull = dEdxFull /density;
186 
187  //Stopping Power from simulation.
188  //
189  G4double meandEdx = fEnergyDeposit/length;
190  G4double stopPower = meandEdx/density;
191 
192  G4cout << "\n ======================== run summary ======================\n";
193 
194  G4int prec = G4cout.precision(3);
195 
196  G4cout << "\n The run was " << TotNbofEvents << " " << partName << " of "
197  << G4BestUnit(energy,"Energy") << " through "
198  << G4BestUnit(length,"Length") << " of "
199  << material->GetName() << " (density: "
200  << G4BestUnit(density,"Volumic Mass") << ")" << G4endl;
201 
202  G4cout.precision(4);
203 
204  G4cout << "\n Total energy deposit in absorber per event = "
205  << G4BestUnit(fEnergyDeposit,"Energy") << " +- "
206  << G4BestUnit(rmsEdep, "Energy")
207  << G4endl;
208 
209  G4cout << "\n -----> Mean dE/dx = " << meandEdx/(MeV/cm) << " MeV/cm"
210  << "\t(" << stopPower/(MeV*cm2/g) << " MeV*cm2/g)"
211  << G4endl;
212 
213  G4cout << "\n From formulas :" << G4endl;
214  G4cout << " restricted dEdx = " << dEdxTable/(MeV/cm) << " MeV/cm"
215  << "\t(" << stopTable/(MeV*cm2/g) << " MeV*cm2/g)"
216  << G4endl;
217 
218  G4cout << " full dEdx = " << dEdxFull/(MeV/cm) << " MeV/cm"
219  << "\t(" << stopFull/(MeV*cm2/g) << " MeV*cm2/g)"
220  << G4endl;
221 
222  G4cout << "\n Leakage : primary = "
223  << G4BestUnit(fEnergyLeak[0],"Energy") << " +- "
224  << G4BestUnit(rmsEl0, "Energy")
225  << " secondaries = "
226  << G4BestUnit(fEnergyLeak[1],"Energy") << " +- "
227  << G4BestUnit(rmsEl1, "Energy")
228  << G4endl;
229 
230  G4cout << " Energy balance : edep + eleak = "
231  << G4BestUnit(EnergyBalance,"Energy")
232  << G4endl;
233 
234  G4cout << "\n Total track length (charged) in absorber per event = "
235  << G4BestUnit(fTrakLenCharged,"Length") << " +- "
236  << G4BestUnit(rmsTLCh, "Length") << G4endl;
237 
238  G4cout << " Total track length (neutral) in absorber per event = "
239  << G4BestUnit(fTrakLenNeutral,"Length") << " +- "
240  << G4BestUnit(rmsTLNe, "Length") << G4endl;
241 
242  G4cout << "\n Number of steps (charged) in absorber per event = "
243  << fNbStepsCharged << " +- " << rmsStCh << G4endl;
244 
245  G4cout << " Number of steps (neutral) in absorber per event = "
246  << fNbStepsNeutral << " +- " << rmsStNe << G4endl;
247 
248  G4cout << "\n Number of secondaries per event : Gammas = " << Gamma
249  << "; electrons = " << Elect
250  << "; positrons = " << Posit << G4endl;
251 
252  G4cout << "\n Number of events with the primary particle transmitted = "
253  << transmit[1] << " %" << G4endl;
254 
255  G4cout << " Number of events with at least 1 particle transmitted "
256  << "(same charge as primary) = " << transmit[0] << " %" << G4endl;
257 
258  G4cout << "\n Number of events with the primary particle reflected = "
259  << reflect[1] << " %" << G4endl;
260 
261  G4cout << " Number of events with at least 1 particle reflected "
262  << "(same charge as primary) = " << reflect[0] << " %" << G4endl;
263 
264  // compute width of the Gaussian central part of the MultipleScattering
265  //
266  G4cout << "\n MultipleScattering:"
267  << "\n rms proj angle of transmit primary particle = "
268  << rmsMsc/mrad << " mrad (central part only)" << G4endl;
269 
270  G4cout << " computed theta0 (Highland formula) = "
271  << ComputeMscHighland()/mrad << " mrad" << G4endl;
272 
273  G4cout << " central part defined as +- "
274  << fMscThetaCentral/mrad << " mrad; "
275  << " Tail ratio = " << tailMsc << " %" << G4endl;
276 
277  G4cout.precision(prec);
278 
279  // normalize histograms
280  //
281  G4AnalysisManager* analysisManager = G4AnalysisManager::Instance();
282 
283  G4int ih = 1;
284  G4double binWidth = analysisManager->GetH1Width(ih);
285  G4double unit = analysisManager->GetH1Unit(ih);
286  G4double fac = unit/(TotNbofEvents*binWidth);
287  analysisManager->ScaleH1(ih,fac);
288 
289  ih = 10;
290  binWidth = analysisManager->GetH1Width(ih);
291  unit = analysisManager->GetH1Unit(ih);
292  fac = unit/(TotNbofEvents*binWidth);
293  analysisManager->ScaleH1(ih,fac);
294 
295  ih = 12;
296  analysisManager->ScaleH1(ih,1./TotNbofEvents);
297 
298  // save histograms
299  if ( analysisManager->IsActive() ) {
300  analysisManager->Write();
301  analysisManager->CloseFile();
302  }
303 
304  // show Rndm status
305  CLHEP::HepRandom::showEngineStatus();
306 }
307 
308 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
309 
310 G4double RunAction::ComputeMscHighland()
311 {
312  //compute the width of the Gaussian central part of the MultipleScattering
313  //projected angular distribution.
314  //Eur. Phys. Jour. C15 (2000) page 166, formule 23.9
315 
316  G4double t = (fDetector->GetAbsorberThickness())
317  /(fDetector->GetAbsorberMaterial()->GetRadlen());
318  if (t < DBL_MIN) return 0.;
319 
320  G4ParticleGun* particle = fPrimary->GetParticleGun();
321  G4double T = particle->GetParticleEnergy();
322  G4double M = particle->GetParticleDefinition()->GetPDGMass();
323  G4double z = std::abs(particle->GetParticleDefinition()->GetPDGCharge()/eplus);
324 
325  G4double bpc = T*(T+2*M)/(T+M);
326  G4double teta0 = 13.6*MeV*z*std::sqrt(t)*(1.+0.038*std::log(t))/bpc;
327  return teta0;
328 }
329 
330 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
Generated on Sat May 25 2013 14:31:59 for Geant4 by   doxygen 1.8.4