TestEm15/src/RunAction.cc

0001 //
0002 // ********************************************************************
0003 // * License and Disclaimer                                           *
0004 // *                                                                  *
0005 // * The  Geant4 software  is  copyright of the Copyright Holders  of *
0006 // * the Geant4 Collaboration.  It is provided  under  the terms  and *
0007 // * conditions of the Geant4 Software License,  included in the file *
0008 // * LICENSE and available at  http://cern.ch/geant4/license .  These *
0009 // * include a list of copyright holders.                             *
0010 // *                                                                  *
0011 // * Neither the authors of this software system, nor their employing *
0012 // * institutes,nor the agencies providing financial support for this *
0013 // * work  make  any representation or  warranty, express or implied, *
0014 // * regarding  this  software system or assume any liability for its *
0015 // * use.  Please see the license in the file  LICENSE  and URL above *
0016 // * for the full disclaimer and the limitation of liability.         *
0017 // *                                                                  *
0018 // * This  code  implementation is the result of  the  scientific and *
0019 // * technical work of the GEANT4 collaboration.                      *
0020 // * By using,  copying,  modifying or  distributing the software (or *
0021 // * any work based  on the software)  you  agree  to acknowledge its *
0022 // * use  in  resulting  scientific  publications,  and indicate your *
0023 // * acceptance of all terms of the Geant4 Software license.          *
0024 // ********************************************************************
0025 //
0026 /// \file electromagnetic/TestEm15/src/RunAction.cc
0027 /// \brief Implementation of the RunAction class
0028 //
0029 //
0030 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0031 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0032
0033 #include "RunAction.hh"
0034
0035 #include "DetectorConstruction.hh"
0036 #include "HistoManager.hh"
0037 #include "PrimaryGeneratorAction.hh"
0038
0039 #include "G4EmCalculator.hh"
0040 #include "G4Run.hh"
0041 #include "G4RunManager.hh"
0042 #include "G4SystemOfUnits.hh"
0043 #include "G4UnitsTable.hh"
0044 #include "Randomize.hh"
0045
0046 #include <iomanip>
0047
0048 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0049
0050 RunAction::RunAction(DetectorConstruction* det, PrimaryGeneratorAction* prim)
0051   : G4UserRunAction(), fDetector(det), fPrimary(prim), fProcCounter(0), fHistoManager(0)
0052 {
0053   fHistoManager = new HistoManager();
0054 }
0055
0056 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0057
0058 RunAction::~RunAction()
0059 {
0060   delete fHistoManager;
0061 }
0062
0063 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0064
0065 void RunAction::BeginOfRunAction(const G4Run*)
0066 {
0067   // save Rndm status
0068   ////G4RunManager::GetRunManager()->SetRandomNumberStore(true);
0069   CLHEP::HepRandom::showEngineStatus();
0070
0071   fProcCounter = new ProcessesCount;
0072   fTotalCount = 0;
0073
0074   fTruePL = fTruePL2 = fGeomPL = fGeomPL2 = 0.;
0075   fLDispl = fLDispl2 = fPsiSpa = fPsiSpa2 = 0.;
0076   fTetPrj = fTetPrj2 = 0.;
0077   fPhiCor = fPhiCor2 = 0.;
0078
0079   // histograms
0080   //
0081   G4AnalysisManager* analysisManager = G4AnalysisManager::Instance();
0082   if (analysisManager->IsActive()) {
0083     analysisManager->OpenFile();
0084   }
0085 }
0086
0087 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0088
0089 void RunAction::CountProcesses(G4String procName)
0090 {
0091   // does the process  already encounted ?
0092   size_t nbProc = fProcCounter->size();
0093   size_t i = 0;
0094   while ((i < nbProc) && ((*fProcCounter)[i]->GetName() != procName))
0095     i++;
0096   if (i == nbProc) fProcCounter->push_back(new OneProcessCount(procName));
0097
0098   (*fProcCounter)[i]->Count();
0099 }
0100
0101 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0102
0103 void RunAction::EndOfRunAction(const G4Run* aRun)
0104 {
0105   G4int NbOfEvents = aRun->GetNumberOfEvent();
0106   if (NbOfEvents == 0) return;
0107
0108   G4int prec = G4cout.precision(5);
0109
0110   G4Material* material = fDetector->GetMaterial();
0111   G4double density = material->GetDensity();
0112
0113   G4ParticleDefinition* particle = fPrimary->GetParticleGun()->GetParticleDefinition();
0114   G4String Particle = particle->GetParticleName();
0115   G4double energy = fPrimary->GetParticleGun()->GetParticleEnergy();
0116   G4cout << "\n The run consists of " << NbOfEvents << " " << Particle << " of "
0117          << G4BestUnit(energy, "Energy") << " through "
0118          << G4BestUnit(fDetector->GetBoxSize(), "Length") << " of " << material->GetName()
0119          << " (density: " << G4BestUnit(density, "Volumic Mass") << ")" << G4endl;
0120
0121   // frequency of processes
0122   G4cout << "\n Process calls frequency --->";
0123   for (size_t i = 0; i < fProcCounter->size(); i++) {
0124     G4String procName = (*fProcCounter)[i]->GetName();
0125     G4int count = (*fProcCounter)[i]->GetCounter();
0126     G4cout << "\t" << procName << " = " << count;
0127   }
0128
0129   if (fTotalCount > 0) {
0130     // compute path length and related quantities
0131     //
0132     G4double MeanTPL = fTruePL / fTotalCount;
0133     G4double MeanTPL2 = fTruePL2 / fTotalCount;
0134     G4double rmsTPL = std::sqrt(std::fabs(MeanTPL2 - MeanTPL * MeanTPL));
0135
0136     G4double MeanGPL = fGeomPL / fTotalCount;
0137     G4double MeanGPL2 = fGeomPL2 / fTotalCount;
0138     G4double rmsGPL = std::sqrt(std::fabs(MeanGPL2 - MeanGPL * MeanGPL));
0139
0140     G4double MeanLaD = fLDispl / fTotalCount;
0141     G4double MeanLaD2 = fLDispl2 / fTotalCount;
0142     G4double rmsLaD = std::sqrt(std::fabs(MeanLaD2 - MeanLaD * MeanLaD));
0143
0144     G4double MeanPsi = fPsiSpa / (fTotalCount);
0145     G4double MeanPsi2 = fPsiSpa2 / (fTotalCount);
0146     G4double rmsPsi = std::sqrt(std::fabs(MeanPsi2 - MeanPsi * MeanPsi));
0147
0148     G4double MeanTeta = fTetPrj / (2 * fTotalCount);
0149     G4double MeanTeta2 = fTetPrj2 / (2 * fTotalCount);
0150     G4double rmsTeta = std::sqrt(std::fabs(MeanTeta2 - MeanTeta * MeanTeta));
0151
0152     G4double MeanCorrel = fPhiCor / (fTotalCount);
0153     G4double MeanCorrel2 = fPhiCor2 / (fTotalCount);
0154     G4double rmsCorrel = std::sqrt(std::fabs(MeanCorrel2 - MeanCorrel * MeanCorrel));
0155
0156     G4cout << "\n\n truePathLength :\t" << G4BestUnit(MeanTPL, "Length") << " +- "
0157            << G4BestUnit(rmsTPL, "Length") << "\n geomPathLength :\t"
0158            << G4BestUnit(MeanGPL, "Length") << " +- " << G4BestUnit(rmsGPL, "Length")
0159            << "\n lateralDisplac :\t" << G4BestUnit(MeanLaD, "Length") << " +- "
0160            << G4BestUnit(rmsLaD, "Length") << "\n Psi            :\t" << MeanPsi / mrad << " mrad"
0161            << " +- " << rmsPsi / mrad << " mrad"
0162            << "  (" << MeanPsi / deg << " deg"
0163            << " +- " << rmsPsi / deg << " deg)" << G4endl;
0164
0165     G4cout << "\n Theta_plane    :\t" << rmsTeta / mrad << " mrad"
0166            << "  (" << rmsTeta / deg << " deg)"
0167            << "\n phi correlation:\t" << MeanCorrel << " +- " << rmsCorrel
0168            << "  (std::cos(phi_pos - phi_dir))" << G4endl;
0169
0170     // cross check from G4EmCalculator
0171     //
0172     G4cout << "\n Verification from G4EmCalculator. \n";
0173
0174     G4EmCalculator emCal;
0175
0176     // get transport mean free path (for multiple scattering)
0177     G4double MSmfp = emCal.GetMeanFreePath(energy, particle, "msc", material);
0178
0179     // get range from restricted dedx
0180     G4double range = emCal.GetRangeFromRestricteDEDX(energy, particle, material);
0181
0182     // effective facRange
0183     G4double efFacrange = MeanTPL / std::max(MSmfp, range);
0184     if (MeanTPL / range >= 0.99) efFacrange = 1.;
0185
0186     G4cout << "\n transport mean free path :\t" << G4BestUnit(MSmfp, "Length")
0187            << "\n range from restrict dE/dx:\t" << G4BestUnit(range, "Length")
0188            << "\n ---> effective facRange  :\t" << efFacrange << G4endl;
0189
0190     G4cout << "\n compute theta0 from Highland :\t" << ComputeMscHighland(MeanTPL) / mrad << " mrad"
0191            << "  (" << ComputeMscHighland(MeanTPL) / deg << " deg)" << G4endl;
0192   }
0193   else
0194     G4cout << G4endl;
0195
0196   // restore default format
0197   G4cout.precision(prec);
0198
0199   // delete and remove all contents in fProcCounter
0200   while (fProcCounter->size() > 0) {
0201     OneProcessCount* aProcCount = fProcCounter->back();
0202     fProcCounter->pop_back();
0203     delete aProcCount;
0204   }
0205   delete fProcCounter;
0206
0207   // save histograms
0208   G4AnalysisManager* analysisManager = G4AnalysisManager::Instance();
0209   if (analysisManager->IsActive()) {
0210     analysisManager->Write();
0211     analysisManager->CloseFile();
0212   }
0213
0214   // show Rndm status
0215   CLHEP::HepRandom::showEngineStatus();
0216 }
0217
0218 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0219
0220 G4double RunAction::ComputeMscHighland(G4double pathLength)
0221 {
0222   // compute the width of the Gaussian central part of the MultipleScattering
0223   // projected angular distribution.
0224   // Eur. Phys. Jour. C15 (2000) page 166, formule 23.9
0225
0226   G4double t = pathLength / (fDetector->GetMaterial()->GetRadlen());
0227   if (t < DBL_MIN) return 0.;
0228
0229   G4ParticleGun* particle = fPrimary->GetParticleGun();
0230   G4double T = particle->GetParticleEnergy();
0231   G4double M = particle->GetParticleDefinition()->GetPDGMass();
0232   G4double z = std::abs(particle->GetParticleDefinition()->GetPDGCharge() / eplus);
0233
0234   G4double bpc = T * (T + 2 * M) / (T + M);
0235   G4double teta0 = 13.6 * MeV * z * std::sqrt(t) * (1. + 0.038 * std::log(t)) / bpc;
0236   return teta0;
0237 }
0238
0239 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......