Hadr03/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 "G4HadronicProcess.hh"
0036 #include "G4HadronicProcessStore.hh"
0037 #include "G4Neutron.hh"
0038 #include "G4ProcessTable.hh"
0039 #include "G4SystemOfUnits.hh"
0040 #include "G4UnitsTable.hh"
0041
0042 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0043
0044 Run::Run(DetectorConstruction* det) : fDetector(det)
0045 {
0046   for (G4int i = 0; i < 3; i++) {
0047     fPbalance[i] = 0.;
0048   }
0049   for (G4int i = 0; i < 3; i++) {
0050     fNbGamma[i] = 0;
0051   }
0052   fPbalance[1] = DBL_MAX;
0053   fNbGamma[1] = 10000;
0054 }
0055
0056 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0057
0058 void Run::SetPrimary(G4ParticleDefinition* particle, G4double energy)
0059 {
0060   fParticle = particle;
0061   fEkin = energy;
0062 }
0063
0064 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0065
0066 void Run::SetTargetXXX(G4bool flag)
0067 {
0068   fTargetXXX = flag;
0069 }
0070
0071 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0072
0073 void Run::CountProcesses(G4VProcess* process)
0074 {
0075   if (process == nullptr) return;
0076   G4String procName = process->GetProcessName();
0077   std::map<G4String, G4int>::iterator it = fProcCounter.find(procName);
0078   if (it == fProcCounter.end()) {
0079     fProcCounter[procName] = 1;
0080   }
0081   else {
0082     fProcCounter[procName]++;
0083   }
0084 }
0085
0086 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0087
0088 void Run::SumTrack(G4double trackl)
0089 {
0090   fTotalCount++;
0091   fSumTrack += trackl;
0092   fSumTrack2 += trackl * trackl;
0093 }
0094
0095 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0096
0097 void Run::CountNuclearChannel(G4String name, G4double Q)
0098 {
0099   std::map<G4String, NuclChannel>::iterator it = fNuclChannelMap.find(name);
0100   if (it == fNuclChannelMap.end()) {
0101     fNuclChannelMap[name] = NuclChannel(1, Q);
0102   }
0103   else {
0104     NuclChannel& data = it->second;
0105     data.fCount++;
0106     data.fQ += Q;
0107   }
0108 }
0109
0110 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0111
0112 void Run::ParticleCount(G4String name, G4double Ekin)
0113 {
0114   std::map<G4String, ParticleData>::iterator it = fParticleDataMap.find(name);
0115   if (it == fParticleDataMap.end()) {
0116     fParticleDataMap[name] = ParticleData(1, Ekin, Ekin, Ekin);
0117   }
0118   else {
0119     ParticleData& data = it->second;
0120     data.fCount++;
0121     data.fEmean += Ekin;
0122     // update min max
0123     G4double emin = data.fEmin;
0124     if (Ekin < emin) data.fEmin = Ekin;
0125     G4double emax = data.fEmax;
0126     if (Ekin > emax) data.fEmax = Ekin;
0127   }
0128 }
0129 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0130
0131 void Run::Balance(G4double Pbal)
0132 {
0133   fPbalance[0] += Pbal;
0134   // update min max
0135   if (fTotalCount == 1) fPbalance[1] = fPbalance[2] = Pbal;
0136   if (Pbal < fPbalance[1]) fPbalance[1] = Pbal;
0137   if (Pbal > fPbalance[2]) fPbalance[2] = Pbal;
0138 }
0139
0140 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0141
0142 void Run::CountGamma(G4int nGamma)
0143 {
0144   fGammaCount++;
0145   fNbGamma[0] += nGamma;
0146   // update min max
0147   if (fGammaCount == 1) fNbGamma[1] = fNbGamma[2] = nGamma;
0148   if (nGamma < fNbGamma[1]) fNbGamma[1] = nGamma;
0149   if (nGamma > fNbGamma[2]) fNbGamma[2] = nGamma;
0150 }
0151
0152 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0153
0154 void Run::Merge(const G4Run* run)
0155 {
0156   const Run* localRun = static_cast<const Run*>(run);
0157
0158   // primary particle info
0159   //
0160   fParticle = localRun->fParticle;
0161   fEkin = localRun->fEkin;
0162
0163   // accumulate sums
0164   //
0165   fTotalCount += localRun->fTotalCount;
0166   fGammaCount += localRun->fGammaCount;
0167   fSumTrack += localRun->fSumTrack;
0168   fSumTrack2 += localRun->fSumTrack2;
0169
0170   fPbalance[0] += localRun->fPbalance[0];
0171   G4double min, max;
0172   min = localRun->fPbalance[1];
0173   max = localRun->fPbalance[2];
0174   if (fPbalance[1] > min) fPbalance[1] = min;
0175   if (fPbalance[2] < max) fPbalance[2] = max;
0176
0177   fNbGamma[0] += localRun->fNbGamma[0];
0178   G4int nbmin, nbmax;
0179   nbmin = localRun->fNbGamma[1];
0180   nbmax = localRun->fNbGamma[2];
0181   if (fNbGamma[1] > nbmin) fNbGamma[1] = nbmin;
0182   if (fNbGamma[2] < nbmax) fNbGamma[2] = nbmax;
0183
0184   // map: processes count
0185   std::map<G4String, G4int>::const_iterator itp;
0186   for (itp = localRun->fProcCounter.begin(); itp != localRun->fProcCounter.end(); ++itp) {
0187     G4String procName = itp->first;
0188     G4int localCount = itp->second;
0189     if (fProcCounter.find(procName) == fProcCounter.end()) {
0190       fProcCounter[procName] = localCount;
0191     }
0192     else {
0193       fProcCounter[procName] += localCount;
0194     }
0195   }
0196
0197   // map: nuclear channels
0198   std::map<G4String, NuclChannel>::const_iterator itc;
0199   for (itc = localRun->fNuclChannelMap.begin(); itc != localRun->fNuclChannelMap.end(); ++itc) {
0200     G4String name = itc->first;
0201     const NuclChannel& localData = itc->second;
0202     if (fNuclChannelMap.find(name) == fNuclChannelMap.end()) {
0203       fNuclChannelMap[name] = NuclChannel(localData.fCount, localData.fQ);
0204     }
0205     else {
0206       NuclChannel& data = fNuclChannelMap[name];
0207       data.fCount += localData.fCount;
0208       data.fQ += localData.fQ;
0209     }
0210   }
0211
0212   // map: particles count
0213   std::map<G4String, ParticleData>::const_iterator itn;
0214   for (itn = localRun->fParticleDataMap.begin(); itn != localRun->fParticleDataMap.end(); ++itn) {
0215     G4String name = itn->first;
0216     const ParticleData& localData = itn->second;
0217     if (fParticleDataMap.find(name) == fParticleDataMap.end()) {
0218       fParticleDataMap[name] =
0219         ParticleData(localData.fCount, localData.fEmean, localData.fEmin, localData.fEmax);
0220     }
0221     else {
0222       ParticleData& data = fParticleDataMap[name];
0223       data.fCount += localData.fCount;
0224       data.fEmean += localData.fEmean;
0225       G4double emin = localData.fEmin;
0226       if (emin < data.fEmin) data.fEmin = emin;
0227       G4double emax = localData.fEmax;
0228       if (emax > data.fEmax) data.fEmax = emax;
0229     }
0230   }
0231
0232   G4Run::Merge(run);
0233 }
0234
0235 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0236
0237 void Run::EndOfRun(G4bool print)
0238 {
0239   G4int prec = 5, wid = prec + 2;
0240   G4int dfprec = G4cout.precision(prec);
0241
0242   // run condition
0243   //
0244   const G4Material* material = fDetector->GetMaterial();
0245   G4double density = material->GetDensity();
0246
0247   G4String Particle = fParticle->GetParticleName();
0248   G4cout << "\n The run is " << numberOfEvent << " " << Particle << " of "
0249          << G4BestUnit(fEkin, "Energy") << " through " << G4BestUnit(fDetector->GetSize(), "Length")
0250          << " of " << material->GetName() << " (density: " << G4BestUnit(density, "Volumic Mass")
0251          << ")" << G4endl;
0252
0253   if (numberOfEvent == 0) {
0254     G4cout.precision(dfprec);
0255     return;
0256   }
0257
0258   // frequency of processes
0259   //
0260   G4cout << "\n Process calls frequency:" << G4endl;
0261   G4int survive = 0;
0262   std::map<G4String, G4int>::iterator it;
0263   for (it = fProcCounter.begin(); it != fProcCounter.end(); it++) {
0264     G4String procName = it->first;
0265     G4int count = it->second;
0266     G4cout << "\t" << procName << "= " << count;
0267     if (procName == "Transportation") survive = count;
0268   }
0269   G4cout << G4endl;
0270
0271   if (survive > 0) {
0272     G4cout << "\n Nb of incident particles surviving after "
0273            << G4BestUnit(fDetector->GetSize(), "Length") << " of " << material->GetName() << " : "
0274            << survive << G4endl;
0275   }
0276
0277   if (fTotalCount == 0) fTotalCount = 1;  // force printing anyway
0278
0279   // compute mean free path and related quantities
0280   //
0281   G4double MeanFreePath = fSumTrack / fTotalCount;
0282   G4double MeanTrack2 = fSumTrack2 / fTotalCount;
0283   G4double rms = std::sqrt(std::fabs(MeanTrack2 - MeanFreePath * MeanFreePath));
0284   G4double CrossSection = 0.0;
0285   if (MeanFreePath > 0.0) {
0286     CrossSection = 1. / MeanFreePath;
0287   }
0288   G4double massicMFP = MeanFreePath * density;
0289   G4double massicCS = 0.0;
0290   if (massicMFP > 0.0) {
0291     massicCS = 1. / massicMFP;
0292   }
0293
0294   G4cout << "\n\n MeanFreePath:\t" << G4BestUnit(MeanFreePath, "Length") << " +- "
0295          << G4BestUnit(rms, "Length") << "\tmassic: " << G4BestUnit(massicMFP, "Mass/Surface")
0296          << "\n CrossSection:\t" << CrossSection * cm << " cm^-1 "
0297          << "\t\tmassic: " << G4BestUnit(massicCS, "Surface/Mass") << G4endl;
0298
0299   // cross section per atom (only for single material)
0300   //
0301   if (material->GetNumberOfElements() == 1) {
0302     G4double nbAtoms = material->GetTotNbOfAtomsPerVolume();
0303     G4double crossSection = CrossSection / nbAtoms;
0304     G4cout << " crossSection per atom:\t" << G4BestUnit(crossSection, "Surface") << G4endl;
0305   }
0306   // check cross section from G4HadronicProcessStore
0307   //
0308   G4cout << "\n Verification: "
0309          << "crossSections from G4HadronicProcessStore" << G4endl;
0310
0311   G4ProcessTable* processTable = G4ProcessTable::GetProcessTable();
0312   G4HadronicProcessStore* store = G4HadronicProcessStore::Instance();
0313   G4double sumc1 = 0.0, sumc2 = 0.0;
0314   const G4Element* element =
0315     (material->GetNumberOfElements() == 1) ? material->GetElement(0) : nullptr;
0316   for (it = fProcCounter.begin(); it != fProcCounter.end(); ++it) {
0317     G4String procName = it->first;
0318     const G4VProcess* process = processTable->FindProcess(procName, fParticle);
0319     PrintXS(process, material, element, store, density, sumc1, sumc2);
0320   }
0321   if (sumc1 > 0.0) {
0322     G4cout << "\n"
0323            << std::setw(20) << "total"
0324            << " = " << G4BestUnit(sumc1, "Surface/Mass") << "\t";
0325     if (sumc2 > 0.0) {
0326       G4cout << G4BestUnit(sumc2, "Surface");
0327     }
0328     G4cout << G4endl;
0329   }
0330   else {
0331     G4cout << " not available" << G4endl;
0332   }
0333
0334   // nuclear channel count
0335   //
0336   G4cout << "\n List of nuclear reactions: \n" << G4endl;
0337   std::map<G4String, NuclChannel>::iterator ic;
0338   for (ic = fNuclChannelMap.begin(); ic != fNuclChannelMap.end(); ic++) {
0339     G4String name = ic->first;
0340     NuclChannel data = ic->second;
0341     G4int count = data.fCount;
0342     G4double Q = data.fQ / count;
0343     if (print)
0344       G4cout << "  " << std::setw(60) << name << ": " << std::setw(7) << count
0345              << "   Q = " << std::setw(wid) << G4BestUnit(Q, "Energy") << G4endl;
0346   }
0347
0348   // Gamma count
0349   //
0350   if (print && (fGammaCount > 0)) {
0351     G4cout << "\n"
0352            << std::setw(58) << "number of gamma or e- (ic): N = " << fNbGamma[1] << " --> "
0353            << fNbGamma[2] << G4endl;
0354   }
0355
0356   if (print && fTargetXXX) {
0357     G4cout << "\n   --> NOTE: XXXX because neutronHP is unable to return target nucleus" << G4endl;
0358   }
0359
0360   // particles count
0361   //
0362   G4cout << "\n List of generated particles:" << G4endl;
0363
0364   std::map<G4String, ParticleData>::iterator itn;
0365   for (itn = fParticleDataMap.begin(); itn != fParticleDataMap.end(); itn++) {
0366     G4String name = itn->first;
0367     ParticleData data = itn->second;
0368     G4int count = data.fCount;
0369     G4double eMean = data.fEmean / count;
0370     G4double eMin = data.fEmin;
0371     G4double eMax = data.fEmax;
0372     if (print)
0373       G4cout << "  " << std::setw(13) << name << ": " << std::setw(7) << count
0374              << "  Emean = " << std::setw(wid) << G4BestUnit(eMean, "Energy") << "\t( "
0375              << G4BestUnit(eMin, "Energy") << " --> " << G4BestUnit(eMax, "Energy") << ")"
0376              << G4endl;
0377   }
0378
0379   // energy momentum balance
0380   //
0381   if (fTotalCount > 1) {
0382     G4double Pbmean = fPbalance[0] / fTotalCount;
0383     G4cout << "\n   Momentum balance: Pmean = " << std::setw(wid) << G4BestUnit(Pbmean, "Energy")
0384            << "\t( " << G4BestUnit(fPbalance[1], "Energy") << " --> "
0385            << G4BestUnit(fPbalance[2], "Energy") << ") \n"
0386            << G4endl;
0387   }
0388
0389   // normalize histograms
0390   ////G4AnalysisManager* analysisManager = G4AnalysisManager::Instance();
0391   ////G4double factor = 1./numberOfEvent;
0392   ////analysisManager->ScaleH1(3,factor);
0393
0394   // remove all contents in fProcCounter, fCount
0395   fProcCounter.clear();
0396   fNuclChannelMap.clear();
0397   fParticleDataMap.clear();
0398
0399   // restore default format
0400   G4cout.precision(dfprec);
0401 }
0402
0403 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0404
0405 void Run::PrintXS(const G4VProcess* proc, const G4Material* mat, const G4Element* elm,
0406                   G4HadronicProcessStore* store, G4double density, G4double& sum1, G4double& sum2)
0407 {
0408   if (nullptr == proc) {
0409     return;
0410   }
0411   G4double xs1 = store->GetCrossSectionPerVolume(fParticle, fEkin, proc, mat);
0412   G4double massSigma = xs1 / density;
0413   sum1 += massSigma;
0414   if (nullptr != elm) {
0415     G4double xs2 = store->GetCrossSectionPerAtom(fParticle, fEkin, proc, elm, mat);
0416     sum2 += xs2;
0417     G4cout << "\n"
0418            << std::setw(20) << proc->GetProcessName() << " = "
0419            << G4BestUnit(massSigma, "Surface/Mass") << "\t" << G4BestUnit(xs2, "Surface");
0420   }
0421   else {
0422     G4cout << "\n"
0423            << std::setw(20) << proc->GetProcessName() << " = "
0424            << G4BestUnit(massSigma, "Surface/Mass");
0425   }
0426 }
0427
0428 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......