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 RunAction.cc
0027 /// \brief Implementation of the RunAction class
0028
0029 #include "RunAction.hh"
0030
0031 #include "DetectorConstruction.hh"
0032 #include "HistoManager.hh"
0033 #include "PrimaryGeneratorAction.hh"
0034
0035 #include "G4EmCalculator.hh"
0036 #include "G4Run.hh"
0037 #include "G4RunManager.hh"
0038 #include "G4SystemOfUnits.hh"
0039 #include "G4UnitsTable.hh"
0040 #include "Randomize.hh"
0041
0042 #include <iomanip>
0043
0044 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0045
0046 RunAction::RunAction(DetectorConstruction* det, PrimaryGeneratorAction* prim)
0047   : G4UserRunAction(), fDetector(det), fPrimary(prim), fProcCounter(0), fHistoManager(0)
0048 {
0049   fHistoManager = new HistoManager();
0050 }
0051
0052 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0053
0054 RunAction::~RunAction()
0055 {
0056   delete fHistoManager;
0057 }
0058
0059 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0060
0061 void RunAction::BeginOfRunAction(const G4Run*)
0062 {
0063   // save Rndm status
0064   ////G4RunManager::GetRunManager()->SetRandomNumberStore(true);
0065   CLHEP::HepRandom::showEngineStatus();
0066
0067   fProcCounter = new ProcessesCount;
0068   fTotalCount = 0;
0069
0070   fTruePL = fTruePL2 = fGeomPL = fGeomPL2 = 0.;
0071   fLDispl = fLDispl2 = fPsiSpa = fPsiSpa2 = 0.;
0072   fTetPrj = fTetPrj2 = 0.;
0073   fPhiCor = fPhiCor2 = 0.;
0074
0075   // histograms
0076   //
0077   G4AnalysisManager* analysisManager = G4AnalysisManager::Instance();
0078   if (analysisManager->IsActive()) {
0079     analysisManager->OpenFile();
0080   }
0081 }
0082
0083 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0084
0085 void RunAction::CountProcesses(G4String procName)
0086 {
0087   // does the process  already encounted ?
0088   size_t nbProc = fProcCounter->size();
0089   size_t i = 0;
0090   while ((i < nbProc) && ((*fProcCounter)[i]->GetName() != procName))
0091     i++;
0092   if (i == nbProc) fProcCounter->push_back(new OneProcessCount(procName));
0093
0094   (*fProcCounter)[i]->Count();
0095 }
0096
0097 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0098
0099 void RunAction::EndOfRunAction(const G4Run* aRun)
0100 {
0101   G4int NbOfEvents = aRun->GetNumberOfEvent();
0102   if (NbOfEvents == 0) return;
0103
0104   G4int prec = G4cout.precision(5);
0105
0106   G4Material* material = fDetector->GetMaterial();
0107   G4double density = material->GetDensity();
0108
0109   G4ParticleDefinition* particle = fPrimary->GetParticleGun()->GetParticleDefinition();
0110   G4String Particle = particle->GetParticleName();
0111   G4double energy = fPrimary->GetParticleGun()->GetParticleEnergy();
0112   G4cout << "\n The run consists of " << NbOfEvents << " " << Particle << " of "
0113          << G4BestUnit(energy, "Energy") << " through "
0114          << G4BestUnit(fDetector->GetBoxSize(), "Length") << " of " << material->GetName()
0115          << " (density: " << G4BestUnit(density, "Volumic Mass") << ")" << G4endl;
0116
0117   // frequency of processes
0118   G4cout << "\n Process calls frequency --->";
0119   for (size_t i = 0; i < fProcCounter->size(); i++) {
0120     G4String procName = (*fProcCounter)[i]->GetName();
0121     G4int count = (*fProcCounter)[i]->GetCounter();
0122     G4cout << "\t" << procName << " = " << count;
0123   }
0124
0125   if (fTotalCount > 0) {
0126     // compute path length and related quantities
0127     //
0128     G4double MeanTPL = fTruePL / fTotalCount;
0129     G4double MeanTPL2 = fTruePL2 / fTotalCount;
0130     G4double rmsTPL = std::sqrt(std::fabs(MeanTPL2 - MeanTPL * MeanTPL));
0131
0132     G4double MeanGPL = fGeomPL / fTotalCount;
0133     G4double MeanGPL2 = fGeomPL2 / fTotalCount;
0134     G4double rmsGPL = std::sqrt(std::fabs(MeanGPL2 - MeanGPL * MeanGPL));
0135
0136     G4double MeanLaD = fLDispl / fTotalCount;
0137     G4double MeanLaD2 = fLDispl2 / fTotalCount;
0138     G4double rmsLaD = std::sqrt(std::fabs(MeanLaD2 - MeanLaD * MeanLaD));
0139
0140     G4double MeanPsi = fPsiSpa / (fTotalCount);
0141     G4double MeanPsi2 = fPsiSpa2 / (fTotalCount);
0142     G4double rmsPsi = std::sqrt(std::fabs(MeanPsi2 - MeanPsi * MeanPsi));
0143
0144     G4double MeanTeta = fTetPrj / (2 * fTotalCount);
0145     G4double MeanTeta2 = fTetPrj2 / (2 * fTotalCount);
0146     G4double rmsTeta = std::sqrt(std::fabs(MeanTeta2 - MeanTeta * MeanTeta));
0147
0148     G4double MeanCorrel = fPhiCor / (fTotalCount);
0149     G4double MeanCorrel2 = fPhiCor2 / (fTotalCount);
0150     G4double rmsCorrel = std::sqrt(std::fabs(MeanCorrel2 - MeanCorrel * MeanCorrel));
0151
0152     G4cout << "\n\n truePathLength :\t" << G4BestUnit(MeanTPL, "Length") << " +- "
0153            << G4BestUnit(rmsTPL, "Length") << "\n geomPathLength :\t"
0154            << G4BestUnit(MeanGPL, "Length") << " +- " << G4BestUnit(rmsGPL, "Length")
0155            << "\n lateralDisplac :\t" << G4BestUnit(MeanLaD, "Length") << " +- "
0156            << G4BestUnit(rmsLaD, "Length") << "\n Psi            :\t" << MeanPsi / mrad << " mrad"
0157            << " +- " << rmsPsi / mrad << " mrad"
0158            << "  (" << MeanPsi / deg << " deg"
0159            << " +- " << rmsPsi / deg << " deg)" << G4endl;
0160
0161     G4cout << "\n Theta_plane    :\t" << rmsTeta / mrad << " mrad"
0162            << "  (" << rmsTeta / deg << " deg)"
0163            << "\n phi correlation:\t" << MeanCorrel << " +- " << rmsCorrel
0164            << "  (std::cos(phi_pos - phi_dir))" << G4endl;
0165
0166     // cross check from G4EmCalculator
0167     //
0168     G4cout << "\n Verification from G4EmCalculator. \n";
0169
0170     G4EmCalculator emCal;
0171
0172     // get transport mean free path (for multiple scattering)
0173     G4double MSmfp = emCal.GetMeanFreePath(energy, particle, "msc", material);
0174
0175     // get range from restricted dedx
0176     G4double range = emCal.GetRangeFromRestricteDEDX(energy, particle, material);
0177
0178     // effective facRange
0179     G4double efFacrange = MeanTPL / std::max(MSmfp, range);
0180     if (MeanTPL / range >= 0.99) efFacrange = 1.;
0181
0182     G4cout << "\n transport mean free path :\t" << G4BestUnit(MSmfp, "Length")
0183            << "\n range from restrict dE/dx:\t" << G4BestUnit(range, "Length")
0184            << "\n ---> effective facRange  :\t" << efFacrange << G4endl;
0185
0186     G4cout << "\n compute theta0 from Highland :\t" << ComputeMscHighland(MeanTPL) / mrad << " mrad"
0187            << "  (" << ComputeMscHighland(MeanTPL) / deg << " deg)" << G4endl;
0188   }
0189   else
0190     G4cout << G4endl;
0191
0192   // restore default format
0193   G4cout.precision(prec);
0194
0195   // delete and remove all contents in fProcCounter
0196   while (fProcCounter->size() > 0) {
0197     OneProcessCount* aProcCount = fProcCounter->back();
0198     fProcCounter->pop_back();
0199     delete aProcCount;
0200   }
0201   delete fProcCounter;
0202
0203   // save histograms
0204   G4AnalysisManager* analysisManager = G4AnalysisManager::Instance();
0205   if (analysisManager->IsActive()) {
0206     analysisManager->Write();
0207     analysisManager->CloseFile();
0208   }
0209
0210   // show Rndm status
0211   CLHEP::HepRandom::showEngineStatus();
0212 }
0213
0214 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0215
0216 G4double RunAction::ComputeMscHighland(G4double pathLength)
0217 {
0218   // compute the width of the Gaussian central part of the MultipleScattering
0219   // projected angular distribution.
0220   // Eur. Phys. Jour. C15 (2000) page 166, formule 23.9
0221
0222   G4double t = pathLength / (fDetector->GetMaterial()->GetRadlen());
0223   if (t < DBL_MIN) return 0.;
0224
0225   G4ParticleGun* particle = fPrimary->GetParticleGun();
0226   G4double T = particle->GetParticleEnergy();
0227   G4double M = particle->GetParticleDefinition()->GetPDGMass();
0228   G4double z = std::abs(particle->GetParticleDefinition()->GetPDGCharge() / eplus);
0229
0230   G4double bpc = T * (T + 2 * M) / (T + M);
0231   G4double teta0 = 13.6 * MeV * z * std::sqrt(t) * (1. + 0.038 * std::log(t)) / bpc;
0232   return teta0;
0233 }
0234
0235 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......