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