Hadr05/src/Run.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 Run.cc
0027 /// \brief Implementation of the Run class
0028 //
0029 //
0030 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0031 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0032
0033 #include "Run.hh"
0034
0035 #include "DetectorConstruction.hh"
0036 #include "HistoManager.hh"
0037 #include "PrimaryGeneratorAction.hh"
0038
0039 #include "G4ParticleDefinition.hh"
0040 #include "G4SystemOfUnits.hh"
0041 #include "G4Track.hh"
0042 #include "G4UnitsTable.hh"
0043
0044 #include <iomanip>
0045
0046 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0047
0048 Run::Run(DetectorConstruction* det) : fDetector(det)
0049 {
0050   // initialize energy deposited per absorber
0051   //
0052   for (G4int k = 0; k < kMaxAbsor; k++) {
0053     fSumEAbs[k] = fSum2EAbs[k] = fSumLAbs[k] = fSum2LAbs[k] = 0.;
0054   }
0055
0056   // initialize total energy deposited
0057   //
0058   fEdepTot = fEdepTot2 = 0.;
0059
0060   // initialize leakage
0061   //
0062   fEnergyLeak[0] = fEnergyLeak[1] = 0.;
0063   fEleakTot = fEleakTot2 = 0.;
0064
0065   // initialize total energy released
0066   //
0067   fEtotal = fEtotal2 = 0.;
0068
0069   // initialize Eflow
0070   //
0071   G4int nbPlanes = (fDetector->GetNbOfLayers()) * (fDetector->GetNbOfAbsor()) + 2;
0072   fEnergyFlow.resize(nbPlanes);
0073   for (G4int k = 0; k < nbPlanes; k++) {
0074     fEnergyFlow[k] = 0.;
0075   }
0076 }
0077
0078 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0079
0080 void Run::SetPrimary(G4ParticleDefinition* particle, G4double energy)
0081 {
0082   fParticle = particle;
0083   fEkin = energy;
0084 }
0085
0086 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0087
0088 void Run::CountProcesses(const G4VProcess* process)
0089 {
0090   if (process == nullptr) return;
0091   G4String procName = process->GetProcessName();
0092   std::map<G4String, G4int>::iterator it = fProcCounter.find(procName);
0093   if (it == fProcCounter.end()) {
0094     fProcCounter[procName] = 1;
0095   }
0096   else {
0097     fProcCounter[procName]++;
0098   }
0099 }
0100
0101 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0102
0103 void Run::SumEdepPerAbsorber(G4int kAbs, G4double EAbs, G4double LAbs)
0104 {
0105   // accumulate statistic with restriction
0106   //
0107   fSumEAbs[kAbs] += EAbs;
0108   fSum2EAbs[kAbs] += EAbs * EAbs;
0109   fSumLAbs[kAbs] += LAbs;
0110   fSum2LAbs[kAbs] += LAbs * LAbs;
0111 }
0112
0113 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0114
0115 void Run::SumEnergies(G4double edeptot, G4double eleak0, G4double eleak1)
0116 {
0117   fEdepTot += edeptot;
0118   fEdepTot2 += edeptot * edeptot;
0119
0120   fEnergyLeak[0] += eleak0;
0121   fEnergyLeak[1] += eleak1;
0122   G4double eleaktot = eleak0 + eleak1;
0123   fEleakTot += eleaktot;
0124   fEleakTot2 += eleaktot * eleaktot;
0125
0126   G4double etotal = edeptot + eleaktot;
0127   fEtotal += etotal;
0128   fEtotal2 += etotal * etotal;
0129 }
0130
0131 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0132
0133 void Run::SumEnergyFlow(G4int plane, G4double Eflow)
0134 {
0135   fEnergyFlow[plane] += Eflow;
0136 }
0137
0138 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0139
0140 void Run::Merge(const G4Run* run)
0141 {
0142   const Run* localRun = static_cast<const Run*>(run);
0143
0144   // pass information about primary particle
0145   fParticle = localRun->fParticle;
0146   fEkin = localRun->fEkin;
0147
0148   // accumulate sums
0149   //
0150   for (G4int k = 0; k < kMaxAbsor; k++) {
0151     fSumEAbs[k] += localRun->fSumEAbs[k];
0152     fSum2EAbs[k] += localRun->fSum2EAbs[k];
0153     fSumLAbs[k] += localRun->fSumLAbs[k];
0154     fSum2LAbs[k] += localRun->fSum2LAbs[k];
0155   }
0156
0157   fEdepTot += localRun->fEdepTot;
0158   fEdepTot2 += localRun->fEdepTot2;
0159
0160   fEnergyLeak[0] += localRun->fEnergyLeak[0];
0161   fEnergyLeak[1] += localRun->fEnergyLeak[1];
0162
0163   fEleakTot += localRun->fEleakTot;
0164   fEleakTot2 += localRun->fEleakTot2;
0165
0166   fEtotal += localRun->fEtotal;
0167   fEtotal2 += localRun->fEtotal2;
0168
0169   G4int nbPlanes = (fDetector->GetNbOfLayers()) * (fDetector->GetNbOfAbsor()) + 2;
0170   for (G4int k = 0; k < nbPlanes; k++) {
0171     fEnergyFlow[k] += localRun->fEnergyFlow[k];
0172   }
0173
0174   // map: processes count
0175   std::map<G4String, G4int>::const_iterator itp;
0176   for (itp = localRun->fProcCounter.begin(); itp != localRun->fProcCounter.end(); ++itp) {
0177     G4String procName = itp->first;
0178     G4int localCount = itp->second;
0179     if (fProcCounter.find(procName) == fProcCounter.end()) {
0180       fProcCounter[procName] = localCount;
0181     }
0182     else {
0183       fProcCounter[procName] += localCount;
0184     }
0185   }
0186
0187   G4Run::Merge(run);
0188 }
0189
0190 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0191
0192 void Run::EndOfRun()
0193 {
0194   // run condition
0195   //
0196   G4String Particle = fParticle->GetParticleName();
0197   G4cout << "\n ---> The run is " << numberOfEvent << " " << Particle << " of "
0198          << G4BestUnit(fEkin, "Energy") << " through calorimeter" << G4endl;
0199
0200   // frequency of processes
0201   //
0202   G4cout << "\n Process calls frequency :" << G4endl;
0203   G4int index = 0;
0204   std::map<G4String, G4int>::iterator it;
0205   for (it = fProcCounter.begin(); it != fProcCounter.end(); it++) {
0206     G4String procName = it->first;
0207     G4int count = it->second;
0208     G4String space = " ";
0209     if (++index % 3 == 0) space = "\n";
0210     G4cout << " " << std::setw(22) << procName << "=" << std::setw(10) << count << space;
0211   }
0212
0213   G4cout << G4endl;
0214   G4int nEvt = numberOfEvent;
0215   G4double norm = G4double(nEvt);
0216   if (norm > 0) norm = 1. / norm;
0217   G4double qnorm = std::sqrt(norm);
0218
0219   // energy deposit per absorber
0220   //
0221   G4double beamEnergy = fEkin;
0222   G4double sqbeam = std::sqrt(beamEnergy / GeV);
0223
0224   G4double MeanEAbs, MeanEAbs2, rmsEAbs, resolution, rmsres;
0225   G4double MeanLAbs, MeanLAbs2, rmsLAbs;
0226
0227   std::ios::fmtflags mode = G4cout.flags();
0228   G4int prec = G4cout.precision(2);
0229   G4cout << "\n------------------------------------------------------------\n";
0230   G4cout << std::setw(16) << "material" << std::setw(22) << "Edep        rmsE" << std::setw(31)
0231          << "sqrt(E0(GeV))*rmsE/Edep" << std::setw(23) << "total tracklen \n \n";
0232
0233   for (G4int k = 1; k <= fDetector->GetNbOfAbsor(); k++) {
0234     MeanEAbs = fSumEAbs[k] * norm;
0235     MeanEAbs2 = fSum2EAbs[k] * norm;
0236     rmsEAbs = std::sqrt(std::abs(MeanEAbs2 - MeanEAbs * MeanEAbs));
0237
0238     resolution = 100. * sqbeam * rmsEAbs / MeanEAbs;
0239     rmsres = resolution * qnorm;
0240
0241     MeanLAbs = fSumLAbs[k] * norm;
0242     MeanLAbs2 = fSum2LAbs[k] * norm;
0243     rmsLAbs = std::sqrt(std::abs(MeanLAbs2 - MeanLAbs * MeanLAbs));
0244
0245     // print
0246     //
0247     G4cout << std::setw(2) << k << std::setw(14) << fDetector->GetAbsorMaterial(k)->GetName()
0248            << std::setprecision(5) << std::setw(10) << G4BestUnit(MeanEAbs, "Energy")
0249            << std::setprecision(4) << std::setw(8) << G4BestUnit(rmsEAbs, "Energy") << std::setw(10)
0250            << resolution << " +- " << std::setprecision(3) << std::setw(5) << rmsres << " %"
0251            << std::setprecision(4) << std::setw(12) << G4BestUnit(MeanLAbs, "Length") << " +- "
0252            << std::setprecision(3) << std::setw(5) << G4BestUnit(rmsLAbs, "Length") << G4endl;
0253   }
0254
0255   // total energy deposited
0256   //
0257   fEdepTot /= nEvt;
0258   fEdepTot2 /= nEvt;
0259   G4double rmsEdep = std::sqrt(std::abs(fEdepTot2 - fEdepTot * fEdepTot));
0260
0261   G4cout << "\n Total energy deposited = " << std::setprecision(4) << G4BestUnit(fEdepTot, "Energy")
0262          << " +- " << G4BestUnit(rmsEdep, "Energy") << G4endl;
0263
0264   // Energy leakage
0265   //
0266   fEnergyLeak[0] /= nEvt;
0267   fEnergyLeak[1] /= nEvt;
0268   fEleakTot /= nEvt;
0269   fEleakTot2 /= nEvt;
0270   G4double rmsEleak = std::sqrt(std::abs(fEleakTot2 - fEleakTot * fEleakTot));
0271
0272   G4cout << " Leakage :  primary = " << G4BestUnit(fEnergyLeak[0], "Energy")
0273          << "   secondaries = " << G4BestUnit(fEnergyLeak[1], "Energy")
0274          << "  ---> total = " << G4BestUnit(fEleakTot, "Energy") << " +- "
0275          << G4BestUnit(rmsEleak, "Energy") << G4endl;
0276
0277   // total energy released
0278   //
0279   fEtotal /= nEvt;
0280   fEtotal2 /= nEvt;
0281   G4double rmsEtotal = std::sqrt(std::abs(fEtotal2 - fEtotal * fEtotal));
0282
0283   G4cout << " Total energy released :  Edep + Eleak = " << G4BestUnit(fEtotal, "Energy") << " +- "
0284          << G4BestUnit(rmsEtotal, "Energy") << G4endl;
0285   G4cout << "------------------------------------------------------------\n";
0286
0287   // Energy flow
0288   //
0289   G4AnalysisManager* analysis = G4AnalysisManager::Instance();
0290   G4int Idmax = (fDetector->GetNbOfLayers()) * (fDetector->GetNbOfAbsor());
0291   for (G4int Id = 1; Id <= Idmax + 1; Id++) {
0292     analysis->FillH1(2 * kMaxAbsor + 1, (G4double)Id, fEnergyFlow[Id] / nEvt);
0293   }
0294
0295   // normalize histograms
0296   //
0297   for (G4int ih = kMaxAbsor + 1; ih < 2 * kMaxAbsor + 1; ih++) {
0298     analysis->ScaleH1(ih, norm / MeV);
0299   }
0300
0301   // remove all contents in fProcCounter
0302   fProcCounter.clear();
0303
0304   G4cout.setf(mode, std::ios::floatfield);
0305   G4cout.precision(prec);
0306 }
0307
0308 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......