Geant4 Cross Reference

Cross-Referencing   Geant4
Geant4/examples/extended/electromagnetic/TestEm15/src/RunAction.cc

Version: [ ReleaseNotes ] [ 1.0 ] [ 1.1 ] [ 2.0 ] [ 3.0 ] [ 3.1 ] [ 3.2 ] [ 4.0 ] [ 4.0.p1 ] [ 4.0.p2 ] [ 4.1 ] [ 4.1.p1 ] [ 5.0 ] [ 5.0.p1 ] [ 5.1 ] [ 5.1.p1 ] [ 5.2 ] [ 5.2.p1 ] [ 5.2.p2 ] [ 6.0 ] [ 6.0.p1 ] [ 6.1 ] [ 6.2 ] [ 6.2.p1 ] [ 6.2.p2 ] [ 7.0 ] [ 7.0.p1 ] [ 7.1 ] [ 7.1.p1 ] [ 8.0 ] [ 8.0.p1 ] [ 8.1 ] [ 8.1.p1 ] [ 8.1.p2 ] [ 8.2 ] [ 8.2.p1 ] [ 8.3 ] [ 8.3.p1 ] [ 8.3.p2 ] [ 9.0 ] [ 9.0.p1 ] [ 9.0.p2 ] [ 9.1 ] [ 9.1.p1 ] [ 9.1.p2 ] [ 9.1.p3 ] [ 9.2 ] [ 9.2.p1 ] [ 9.2.p2 ] [ 9.2.p3 ] [ 9.2.p4 ] [ 9.3 ] [ 9.3.p1 ] [ 9.3.p2 ] [ 9.4 ] [ 9.4.p1 ] [ 9.4.p2 ] [ 9.4.p3 ] [ 9.4.p4 ] [ 9.5 ] [ 9.5.p1 ] [ 9.5.p2 ] [ 9.6 ] [ 9.6.p1 ] [ 9.6.p2 ] [ 9.6.p3 ] [ 9.6.p4 ] [ 10.0 ] [ 10.0.p1 ] [ 10.0.p2 ] [ 10.0.p3 ] [ 10.0.p4 ] [ 10.1 ] [ 10.1.p1 ] [ 10.1.p2 ] [ 10.1.p3 ] [ 10.2 ] [ 10.2.p1 ] [ 10.2.p2 ] [ 10.2.p3 ] [ 10.3 ] [ 10.3.p1 ] [ 10.3.p2 ] [ 10.3.p3 ] [ 10.4 ] [ 10.4.p1 ] [ 10.4.p2 ] [ 10.4.p3 ] [ 10.5 ] [ 10.5.p1 ] [ 10.6 ] [ 10.6.p1 ] [ 10.6.p2 ] [ 10.6.p3 ] [ 10.7 ] [ 10.7.p1 ] [ 10.7.p2 ] [ 10.7.p3 ] [ 10.7.p4 ] [ 11.0 ] [ 11.0.p1 ] [ 11.0.p2 ] [ 11.0.p3, ] [ 11.0.p4 ] [ 11.1 ] [ 11.1.1 ] [ 11.1.2 ] [ 11.1.3 ] [ 11.2 ] [ 11.2.1 ] [ 11.2.2 ] [ 11.3.0 ]

  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 /// \file electromagnetic/TestEm15/src/RunAction.cc
 27 /// \brief Implementation of the RunAction class
 28 //
 29 //
 30 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
 31 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
 32 
 33 #include "RunAction.hh"
 34 
 35 #include "DetectorConstruction.hh"
 36 #include "HistoManager.hh"
 37 #include "PrimaryGeneratorAction.hh"
 38 
 39 #include "G4EmCalculator.hh"
 40 #include "G4Run.hh"
 41 #include "G4RunManager.hh"
 42 #include "G4SystemOfUnits.hh"
 43 #include "G4UnitsTable.hh"
 44 #include "Randomize.hh"
 45 
 46 #include <iomanip>
 47 
 48 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
 49 
 50 RunAction::RunAction(DetectorConstruction* det, PrimaryGeneratorAction* prim)
 51   : G4UserRunAction(), fDetector(det), fPrimary(prim), fProcCounter(0), fHistoManager(0)
 52 {
 53   fHistoManager = new HistoManager();
 54 }
 55 
 56 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
 57 
 58 RunAction::~RunAction()
 59 {
 60   delete fHistoManager;
 61 }
 62 
 63 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
 64 
 65 void RunAction::BeginOfRunAction(const G4Run*)
 66 {
 67   // save Rndm status
 68   ////G4RunManager::GetRunManager()->SetRandomNumberStore(true);
 69   CLHEP::HepRandom::showEngineStatus();
 70 
 71   fProcCounter = new ProcessesCount;
 72   fTotalCount = 0;
 73 
 74   fTruePL = fTruePL2 = fGeomPL = fGeomPL2 = 0.;
 75   fLDispl = fLDispl2 = fPsiSpa = fPsiSpa2 = 0.;
 76   fTetPrj = fTetPrj2 = 0.;
 77   fPhiCor = fPhiCor2 = 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::CountProcesses(G4String procName)
 90 {
 91   // does the process  already encounted ?
 92   size_t nbProc = fProcCounter->size();
 93   size_t i = 0;
 94   while ((i < nbProc) && ((*fProcCounter)[i]->GetName() != procName))
 95     i++;
 96   if (i == nbProc) fProcCounter->push_back(new OneProcessCount(procName));
 97 
 98   (*fProcCounter)[i]->Count();
 99 }
100 
101 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
102 
103 void RunAction::EndOfRunAction(const G4Run* aRun)
104 {
105   G4int NbOfEvents = aRun->GetNumberOfEvent();
106   if (NbOfEvents == 0) return;
107 
108   G4int prec = G4cout.precision(5);
109 
110   G4Material* material = fDetector->GetMaterial();
111   G4double density = material->GetDensity();
112 
113   G4ParticleDefinition* particle = fPrimary->GetParticleGun()->GetParticleDefinition();
114   G4String Particle = particle->GetParticleName();
115   G4double energy = fPrimary->GetParticleGun()->GetParticleEnergy();
116   G4cout << "\n The run consists of " << NbOfEvents << " " << Particle << " of "
117          << G4BestUnit(energy, "Energy") << " through "
118          << G4BestUnit(fDetector->GetBoxSize(), "Length") << " of " << material->GetName()
119          << " (density: " << G4BestUnit(density, "Volumic Mass") << ")" << G4endl;
120 
121   // frequency of processes
122   G4cout << "\n Process calls frequency --->";
123   for (size_t i = 0; i < fProcCounter->size(); i++) {
124     G4String procName = (*fProcCounter)[i]->GetName();
125     G4int count = (*fProcCounter)[i]->GetCounter();
126     G4cout << "\t" << procName << " = " << count;
127   }
128 
129   if (fTotalCount > 0) {
130     // compute path length and related quantities
131     //
132     G4double MeanTPL = fTruePL / fTotalCount;
133     G4double MeanTPL2 = fTruePL2 / fTotalCount;
134     G4double rmsTPL = std::sqrt(std::fabs(MeanTPL2 - MeanTPL * MeanTPL));
135 
136     G4double MeanGPL = fGeomPL / fTotalCount;
137     G4double MeanGPL2 = fGeomPL2 / fTotalCount;
138     G4double rmsGPL = std::sqrt(std::fabs(MeanGPL2 - MeanGPL * MeanGPL));
139 
140     G4double MeanLaD = fLDispl / fTotalCount;
141     G4double MeanLaD2 = fLDispl2 / fTotalCount;
142     G4double rmsLaD = std::sqrt(std::fabs(MeanLaD2 - MeanLaD * MeanLaD));
143 
144     G4double MeanPsi = fPsiSpa / (fTotalCount);
145     G4double MeanPsi2 = fPsiSpa2 / (fTotalCount);
146     G4double rmsPsi = std::sqrt(std::fabs(MeanPsi2 - MeanPsi * MeanPsi));
147 
148     G4double MeanTeta = fTetPrj / (2 * fTotalCount);
149     G4double MeanTeta2 = fTetPrj2 / (2 * fTotalCount);
150     G4double rmsTeta = std::sqrt(std::fabs(MeanTeta2 - MeanTeta * MeanTeta));
151 
152     G4double MeanCorrel = fPhiCor / (fTotalCount);
153     G4double MeanCorrel2 = fPhiCor2 / (fTotalCount);
154     G4double rmsCorrel = std::sqrt(std::fabs(MeanCorrel2 - MeanCorrel * MeanCorrel));
155 
156     G4cout << "\n\n truePathLength :\t" << G4BestUnit(MeanTPL, "Length") << " +- "
157            << G4BestUnit(rmsTPL, "Length") << "\n geomPathLength :\t"
158            << G4BestUnit(MeanGPL, "Length") << " +- " << G4BestUnit(rmsGPL, "Length")
159            << "\n lateralDisplac :\t" << G4BestUnit(MeanLaD, "Length") << " +- "
160            << G4BestUnit(rmsLaD, "Length") << "\n Psi            :\t" << MeanPsi / mrad << " mrad"
161            << " +- " << rmsPsi / mrad << " mrad"
162            << "  (" << MeanPsi / deg << " deg"
163            << " +- " << rmsPsi / deg << " deg)" << G4endl;
164 
165     G4cout << "\n Theta_plane    :\t" << rmsTeta / mrad << " mrad"
166            << "  (" << rmsTeta / deg << " deg)"
167            << "\n phi correlation:\t" << MeanCorrel << " +- " << rmsCorrel
168            << "  (std::cos(phi_pos - phi_dir))" << G4endl;
169 
170     // cross check from G4EmCalculator
171     //
172     G4cout << "\n Verification from G4EmCalculator. \n";
173 
174     G4EmCalculator emCal;
175 
176     // get transport mean free path (for multiple scattering)
177     G4double MSmfp = emCal.GetMeanFreePath(energy, particle, "msc", material);
178 
179     // get range from restricted dedx
180     G4double range = emCal.GetRangeFromRestricteDEDX(energy, particle, material);
181 
182     // effective facRange
183     G4double efFacrange = MeanTPL / std::max(MSmfp, range);
184     if (MeanTPL / range >= 0.99) efFacrange = 1.;
185 
186     G4cout << "\n transport mean free path :\t" << G4BestUnit(MSmfp, "Length")
187            << "\n range from restrict dE/dx:\t" << G4BestUnit(range, "Length")
188            << "\n ---> effective facRange  :\t" << efFacrange << G4endl;
189 
190     G4cout << "\n compute theta0 from Highland :\t" << ComputeMscHighland(MeanTPL) / mrad << " mrad"
191            << "  (" << ComputeMscHighland(MeanTPL) / deg << " deg)" << G4endl;
192   }
193   else
194     G4cout << G4endl;
195 
196   // restore default format
197   G4cout.precision(prec);
198 
199   // delete and remove all contents in fProcCounter
200   while (fProcCounter->size() > 0) {
201     OneProcessCount* aProcCount = fProcCounter->back();
202     fProcCounter->pop_back();
203     delete aProcCount;
204   }
205   delete fProcCounter;
206 
207   // save histograms
208   G4AnalysisManager* analysisManager = G4AnalysisManager::Instance();
209   if (analysisManager->IsActive()) {
210     analysisManager->Write();
211     analysisManager->CloseFile();
212   }
213 
214   // show Rndm status
215   CLHEP::HepRandom::showEngineStatus();
216 }
217 
218 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
219 
220 G4double RunAction::ComputeMscHighland(G4double pathLength)
221 {
222   // compute the width of the Gaussian central part of the MultipleScattering
223   // projected angular distribution.
224   // Eur. Phys. Jour. C15 (2000) page 166, formule 23.9
225 
226   G4double t = pathLength / (fDetector->GetMaterial()->GetRadlen());
227   if (t < DBL_MIN) return 0.;
228 
229   G4ParticleGun* particle = fPrimary->GetParticleGun();
230   G4double T = particle->GetParticleEnergy();
231   G4double M = particle->GetParticleDefinition()->GetPDGMass();
232   G4double z = std::abs(particle->GetParticleDefinition()->GetPDGCharge() / eplus);
233 
234   G4double bpc = T * (T + 2 * M) / (T + M);
235   G4double teta0 = 13.6 * MeV * z * std::sqrt(t) * (1. + 0.038 * std::log(t)) / bpc;
236   return teta0;
237 }
238 
239 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
240