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