Pol01/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 "PrimaryGeneratorAction.hh"
0033 #include "ProcessesCount.hh"
0034
0035 #include "G4AccumulableManager.hh"
0036 #include "G4Electron.hh"
0037 #include "G4EmCalculator.hh"
0038 #include "G4Gamma.hh"
0039 #include "G4ParticleDefinition.hh"
0040 #include "G4PhysicalConstants.hh"
0041 #include "G4Positron.hh"
0042 #include "G4Run.hh"
0043 #include "G4RunManager.hh"
0044 #include "G4SystemOfUnits.hh"
0045 #include "G4UnitsTable.hh"
0046
0047 #include <iomanip>
0048
0049 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0050
0051 RunAction::RunAction(DetectorConstruction* det, PrimaryGeneratorAction* prim)
0052   : fDetector(det), fPrimary(prim), fAnalysisManager(0), fTotalEventCount(0)
0053 {
0054   fGamma = G4Gamma::Gamma();
0055   fElectron = G4Electron::Electron();
0056   fPositron = G4Positron::Positron();
0057
0058   auto accumulableManager = G4AccumulableManager::Instance();
0059   auto fPhotonStats = new ParticleStatistics("PhotonStats");
0060   auto fElectronStats = new ParticleStatistics("ElectronStats");
0061   auto fPositronStats = new ParticleStatistics("PositronStats");
0062   auto fProcCounter = new ProcessesCount("ProcCounter");
0063
0064   accumulableManager->Register(fPhotonStats);
0065   accumulableManager->Register(fElectronStats);
0066   accumulableManager->Register(fPositronStats);
0067
0068   accumulableManager->Register(fProcCounter);
0069
0070   BookHisto();
0071 }
0072
0073 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0074
0075 RunAction::~RunAction() {}
0076
0077 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0078
0079 void RunAction::BeginOfRunAction(const G4Run* aRun)
0080 {
0081   G4cout << "### Run " << aRun->GetRunID() << " start." << G4endl;
0082
0083   auto accumulableManager = G4AccumulableManager::Instance();
0084   accumulableManager->Reset();
0085
0086   // save Rndm status
0087   //  G4RunManager::GetRunManager()->SetRandomNumberStore(false);
0088   //  CLHEP::HepRandom::showEngineStatus();
0089
0090   fTotalEventCount = 0;
0091
0092   // Open file histogram file
0093   fAnalysisManager->OpenFile();
0094 }
0095
0096 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0097
0098 void RunAction::FillData(const G4ParticleDefinition* particle, G4double kinEnergy,
0099                          G4double costheta, G4double phi, G4double longitudinalPolarization)
0100 {
0101   auto accManager = G4AccumulableManager::Instance();
0102   G4int id = -1;
0103   if (particle == fGamma) {
0104     dynamic_cast<ParticleStatistics*>(accManager->GetAccumulable("PhotonStats"))
0105       ->FillData(kinEnergy, costheta, longitudinalPolarization);
0106     if (fAnalysisManager) {
0107       id = 1;
0108     }
0109   }
0110   else if (particle == fElectron) {
0111     dynamic_cast<ParticleStatistics*>(accManager->GetAccumulable("ElectronStats"))
0112       ->FillData(kinEnergy, costheta, longitudinalPolarization);
0113     if (fAnalysisManager) {
0114       id = 5;
0115     }
0116   }
0117   else if (particle == fPositron) {
0118     dynamic_cast<ParticleStatistics*>(accManager->GetAccumulable("PositronStats"))
0119       ->FillData(kinEnergy, costheta, longitudinalPolarization);
0120     if (fAnalysisManager) {
0121       id = 9;
0122     }
0123   }
0124   if (id > 0) {
0125     fAnalysisManager->FillH1(id, kinEnergy, 1.0);
0126     fAnalysisManager->FillH1(id + 1, costheta, 1.0);
0127     fAnalysisManager->FillH1(id + 2, phi, 1.0);
0128     fAnalysisManager->FillH1(id + 3, longitudinalPolarization, 1.0);
0129   }
0130 }
0131
0132 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0133
0134 void RunAction::BookHisto()
0135 {
0136   // Always creating analysis manager
0137   fAnalysisManager = G4AnalysisManager::Instance();
0138   fAnalysisManager->SetDefaultFileType("root");
0139   fAnalysisManager->SetActivation(true);
0140   fAnalysisManager->SetVerboseLevel(1);
0141
0142   // Open file histogram file
0143   fAnalysisManager->SetFileName("pol01");
0144
0145   fAnalysisManager->SetFirstHistoId(1);
0146
0147   // Creating an 1-dimensional histograms in the root directory of the tree
0148   const G4String id[] = {"h1", "h2", "h3", "h4", "h5", "h6", "h7", "h8", "h9", "h10", "h11", "h12"};
0149   const G4String title[] = {
0150     "Gamma Energy distribution",  // 1
0151     "Gamma Cos(Theta) distribution",  // 2
0152     "Gamma Phi angular distribution",  // 3
0153     "Gamma longitudinal Polarization",  // 4
0154     "Electron Energy distribution",  // 5
0155     "Electron Cos(Theta) distribution",  // 6
0156     "Electron Phi angular distribution",  // 7
0157     "Electron longitudinal Polarization",  // 8
0158     "Positron Energy distribution",  // 9
0159     "Positron Cos(Theta) distribution",  // 10
0160     "Positron Phi angular distribution",  // 11
0161     "Positron longitudinal Polarization"  // 12
0162   };
0163   G4double vmin, vmax;
0164   G4int nbins = 120;
0165   for (int i = 0; i < 12; ++i) {
0166     G4int j = i - i / 4 * 4;
0167     if (0 == j) {
0168       vmin = 0.;
0169       vmax = 12. * MeV;
0170     }
0171     else if (1 == j) {
0172       vmin = -1.;
0173       vmax = 1.;
0174     }
0175     else if (2 == j) {
0176       vmin = 0.;
0177       vmax = pi;
0178     }
0179     else {
0180       vmin = -1.5;
0181       vmax = 1.5;
0182     }
0183     G4int ih = fAnalysisManager->CreateH1(id[i], title[i], nbins, vmin, vmax);
0184     fAnalysisManager->SetH1Activation(ih, false);
0185   }
0186 }
0187
0188 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0189
0190 void RunAction::SaveHisto(G4int nevents)
0191 {
0192   if (fAnalysisManager) {
0193     if (IsMaster()) {
0194       G4double norm = 1.0 / G4double(nevents);
0195       for (int i = 0; i < 12; ++i) {
0196         fAnalysisManager->ScaleH1(i, norm);
0197       }
0198     }
0199     fAnalysisManager->Write();
0200     fAnalysisManager->CloseFile();
0201   }
0202 }
0203
0204 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0205
0206 void RunAction::CountProcesses(G4String& procName)
0207 {
0208   auto accManager = G4AccumulableManager::Instance();
0209   dynamic_cast<ProcessesCount*>(accManager->GetAccumulable("ProcCounter"))->Count(procName);
0210 }
0211
0212 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0213
0214 void RunAction::EndOfRunAction(const G4Run* aRun)
0215 {
0216   //  Total number of events in run (all threads)
0217   G4int NbOfEvents = aRun->GetNumberOfEventToBeProcessed();
0218   //  G4int NbOfEvents = aRun->GetNumberOfEvent();
0219
0220   if (NbOfEvents == 0) return;
0221
0222   G4int prec = G4cout.precision(5);
0223
0224   G4Material* material = fDetector->GetMaterial();
0225   G4double density = material->GetDensity();
0226
0227   if (fPrimary != nullptr) {
0228     G4ParticleDefinition* particle = fPrimary->GetParticleGun()->GetParticleDefinition();
0229     G4String Particle = particle->GetParticleName();
0230     G4double energy = fPrimary->GetParticleGun()->GetParticleEnergy();
0231     G4cout << "\n The run consists of " << fTotalEventCount << " " << Particle << " of "
0232            << G4BestUnit(energy, "Energy") << " through "
0233            << G4BestUnit(fDetector->GetBoxSizeZ(), "Length") << " of " << material->GetName()
0234            << " (density: " << G4BestUnit(density, "Volumic Mass") << ")" << G4endl;
0235   }
0236   // cross check from G4EmCalculator
0237   //   G4cout << "\n Verification from G4EmCalculator. \n";
0238   //   G4EmCalculator emCal;
0239
0240   auto accManager = G4AccumulableManager::Instance();
0241   accManager->Merge();
0242
0243   if (IsMaster()) {
0244     // frequency of processes
0245     G4cout << "\n Process calls frequency --->\n";
0246     dynamic_cast<ProcessesCount*>(accManager->GetAccumulable("ProcCounter"))->Print();
0247     G4cout << " Gamma: \n";
0248     dynamic_cast<ParticleStatistics*>(accManager->GetAccumulable("PhotonStats"))
0249       ->PrintResults(NbOfEvents);
0250     G4cout << " Electron: \n";
0251     dynamic_cast<ParticleStatistics*>(accManager->GetAccumulable("ElectronStats"))
0252       ->PrintResults(NbOfEvents);
0253     G4cout << " Positron: \n";
0254     dynamic_cast<ParticleStatistics*>(accManager->GetAccumulable("PositronStats"))
0255       ->PrintResults(NbOfEvents);
0256     G4cout << G4endl;
0257   }
0258
0259   // restore default format
0260   G4cout.precision(prec);
0261
0262   // write out histograms
0263   SaveHisto(NbOfEvents);
0264
0265   if (IsMaster()) {
0266     // show Rndm status
0267     CLHEP::HepRandom::showEngineStatus();
0268   }
0269 }
0270
0271 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0272
0273 void RunAction::EventFinished()
0274 {
0275   auto accManager = G4AccumulableManager::Instance();
0276
0277   ++fTotalEventCount;
0278
0279   dynamic_cast<ParticleStatistics*>(accManager->GetAccumulable("PhotonStats"))->EventFinished();
0280   dynamic_cast<ParticleStatistics*>(accManager->GetAccumulable("ElectronStats"))->EventFinished();
0281   dynamic_cast<ParticleStatistics*>(accManager->GetAccumulable("PositronStats"))->EventFinished();
0282 }
0283
0284 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0285
0286 RunAction::ParticleStatistics::ParticleStatistics(const G4String& name)
0287   : G4VAccumulable(name),
0288     fCurrentNumber(0),
0289     fTotalNumber(0),
0290     fTotalNumber2(0),
0291     fSumEnergy(0),
0292     fSumEnergy2(0),
0293     fSumPolarization(0),
0294     fSumPolarization2(0),
0295     fSumCosTheta(0),
0296     fSumCosTheta2(0)
0297 {}
0298
0299 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0300
0301 RunAction::ParticleStatistics::~ParticleStatistics() {}
0302
0303 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0304
0305 void RunAction::ParticleStatistics::EventFinished()
0306 {
0307   fTotalNumber += fCurrentNumber;
0308   fTotalNumber2 += fCurrentNumber * fCurrentNumber;
0309   fCurrentNumber = 0;
0310 }
0311 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0312
0313 void RunAction::ParticleStatistics::FillData(G4double kinEnergy, G4double costheta,
0314                                              G4double longitudinalPolarization)
0315 {
0316   ++fCurrentNumber;
0317   fSumEnergy += kinEnergy;
0318   fSumEnergy2 += kinEnergy * kinEnergy;
0319   fSumPolarization += longitudinalPolarization;
0320   fSumPolarization2 += longitudinalPolarization * longitudinalPolarization;
0321   fSumCosTheta += costheta;
0322   fSumCosTheta2 += costheta * costheta;
0323 }
0324
0325 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0326
0327 void RunAction::ParticleStatistics::PrintResults(G4int totalNumberOfEvents)
0328 {
0329   G4cout << "Mean Number per Event :" << G4double(fTotalNumber) / G4double(totalNumberOfEvents)
0330          << "\n";
0331   if (fTotalNumber == 0) fTotalNumber = 1;
0332   G4double energyMean = fSumEnergy / fTotalNumber;
0333   G4double energyRms = std::sqrt(fSumEnergy2 / fTotalNumber - energyMean * energyMean);
0334   G4cout << "Mean Energy :" << G4BestUnit(energyMean, "Energy") << " +- "
0335          << G4BestUnit(energyRms, "Energy") << "\n";
0336   G4double polarizationMean = fSumPolarization / fTotalNumber;
0337   G4double polarizationRms =
0338     std::sqrt(fSumPolarization2 / fTotalNumber - polarizationMean * polarizationMean);
0339   G4cout << "Mean Polarization :" << polarizationMean << " +- " << polarizationRms << "\n";
0340 }
0341
0342 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0343
0344 void RunAction::ParticleStatistics::Reset()
0345 {
0346   fCurrentNumber = 0;
0347   fTotalNumber = fTotalNumber2 = 0;
0348   fSumEnergy = fSumEnergy2 = 0;
0349   fSumPolarization = fSumPolarization2 = 0;
0350   fSumCosTheta = fSumCosTheta2 = 0;
0351 }
0352
0353 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0354
0355 void RunAction::ParticleStatistics::Merge(const G4VAccumulable& other)
0356 {
0357   auto rstat = dynamic_cast<const RunAction::ParticleStatistics&>(other);
0358
0359   fCurrentNumber += rstat.fCurrentNumber;
0360   fTotalNumber += rstat.fTotalNumber;
0361   fTotalNumber2 += rstat.fTotalNumber2;
0362   fSumEnergy += rstat.fSumEnergy;
0363   fSumEnergy2 += rstat.fSumEnergy2;
0364   fSumPolarization += rstat.fSumPolarization;
0365   fSumPolarization2 += rstat.fSumPolarization2;
0366   fSumCosTheta += rstat.fSumCosTheta;
0367   fSumCosTheta2 += rstat.fSumCosTheta2;
0368 }
0369
0370 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......