NeutronSource/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 "G4SystemOfUnits.hh"
0040 #include "G4UnitsTable.hh"
0041
0042 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0043
0044 Run::Run(DetectorConstruction* det) : fDetector(det) {}
0045
0046 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0047
0048 void Run::SetPrimary(G4ParticleDefinition* particle, G4double energy)
0049 {
0050   fParticle = particle;
0051   fEkin = energy;
0052 }
0053
0054 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0055
0056 void Run::CountProcesses(const G4VProcess* process)
0057 {
0058   if (process == nullptr) return;
0059   G4String procName = process->GetProcessName();
0060   std::map<G4String, G4int>::iterator it = fProcCounter.find(procName);
0061   if (it == fProcCounter.end()) {
0062     fProcCounter[procName] = 1;
0063   }
0064   else {
0065     fProcCounter[procName]++;
0066   }
0067 }
0068
0069 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0070
0071 void Run::ParticleCount(G4String name, G4double Ekin, G4double meanLife)
0072 {
0073   std::map<G4String, ParticleData>::iterator it = fParticleDataMap1.find(name);
0074   if (it == fParticleDataMap1.end()) {
0075     fParticleDataMap1[name] = ParticleData(1, Ekin, Ekin, Ekin, meanLife);
0076   }
0077   else {
0078     ParticleData& data = it->second;
0079     data.fCount++;
0080     data.fEmean += Ekin;
0081     // update min max
0082     G4double emin = data.fEmin;
0083     if (Ekin < emin) data.fEmin = Ekin;
0084     G4double emax = data.fEmax;
0085     if (Ekin > emax) data.fEmax = Ekin;
0086     data.fTmean = meanLife;
0087   }
0088 }
0089
0090 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0091
0092 void Run::AddEdep(G4double edep)
0093 {
0094   fEnergyDeposit += edep;
0095   fEnergyDeposit2 += edep * edep;
0096 }
0097
0098 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0099
0100 void Run::AddEflow(G4double eflow)
0101 {
0102   fEnergyFlow += eflow;
0103   fEnergyFlow2 += eflow * eflow;
0104 }
0105 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0106
0107 void Run::ParticleFlux(G4String name, G4double Ekin)
0108 {
0109   std::map<G4String, ParticleData>::iterator it = fParticleDataMap2.find(name);
0110   if (it == fParticleDataMap2.end()) {
0111     fParticleDataMap2[name] = ParticleData(1, Ekin, Ekin, Ekin, -1 * ns);
0112   }
0113   else {
0114     ParticleData& data = it->second;
0115     data.fCount++;
0116     data.fEmean += Ekin;
0117     // update min max
0118     G4double emin = data.fEmin;
0119     if (Ekin < emin) data.fEmin = Ekin;
0120     G4double emax = data.fEmax;
0121     if (Ekin > emax) data.fEmax = Ekin;
0122     data.fTmean = -1 * ns;
0123   }
0124 }
0125
0126 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0127
0128 void Run::Merge(const G4Run* run)
0129 {
0130   const Run* localRun = static_cast<const Run*>(run);
0131
0132   // primary particle info
0133   //
0134   fParticle = localRun->fParticle;
0135   fEkin = localRun->fEkin;
0136
0137   // accumulate sums
0138   //
0139   fEnergyDeposit += localRun->fEnergyDeposit;
0140   fEnergyDeposit2 += localRun->fEnergyDeposit2;
0141   fEnergyFlow += localRun->fEnergyFlow;
0142   fEnergyFlow2 += localRun->fEnergyFlow2;
0143
0144   // map: processes count
0145   std::map<G4String, G4int>::const_iterator itp;
0146   for (itp = localRun->fProcCounter.begin(); itp != localRun->fProcCounter.end(); ++itp) {
0147     G4String procName = itp->first;
0148     G4int localCount = itp->second;
0149     if (fProcCounter.find(procName) == fProcCounter.end()) {
0150       fProcCounter[procName] = localCount;
0151     }
0152     else {
0153       fProcCounter[procName] += localCount;
0154     }
0155   }
0156
0157   // map: created particles count
0158   std::map<G4String, ParticleData>::const_iterator itc;
0159   for (itc = localRun->fParticleDataMap1.begin(); itc != localRun->fParticleDataMap1.end(); ++itc) {
0160     G4String name = itc->first;
0161     const ParticleData& localData = itc->second;
0162     if (fParticleDataMap1.find(name) == fParticleDataMap1.end()) {
0163       fParticleDataMap1[name] = ParticleData(localData.fCount, localData.fEmean, localData.fEmin,
0164                                              localData.fEmax, localData.fTmean);
0165     }
0166     else {
0167       ParticleData& data = fParticleDataMap1[name];
0168       data.fCount += localData.fCount;
0169       data.fEmean += localData.fEmean;
0170       G4double emin = localData.fEmin;
0171       if (emin < data.fEmin) data.fEmin = emin;
0172       G4double emax = localData.fEmax;
0173       if (emax > data.fEmax) data.fEmax = emax;
0174       data.fTmean = localData.fTmean;
0175     }
0176   }
0177
0178   // map: particles flux count
0179   std::map<G4String, ParticleData>::const_iterator itn;
0180   for (itn = localRun->fParticleDataMap2.begin(); itn != localRun->fParticleDataMap2.end(); ++itn) {
0181     G4String name = itn->first;
0182     const ParticleData& localData = itn->second;
0183     if (fParticleDataMap2.find(name) == fParticleDataMap2.end()) {
0184       fParticleDataMap2[name] = ParticleData(localData.fCount, localData.fEmean, localData.fEmin,
0185                                              localData.fEmax, localData.fTmean);
0186     }
0187     else {
0188       ParticleData& data = fParticleDataMap2[name];
0189       data.fCount += localData.fCount;
0190       data.fEmean += localData.fEmean;
0191       G4double emin = localData.fEmin;
0192       if (emin < data.fEmin) data.fEmin = emin;
0193       G4double emax = localData.fEmax;
0194       if (emax > data.fEmax) data.fEmax = emax;
0195       data.fTmean = localData.fTmean;
0196     }
0197   }
0198
0199   G4Run::Merge(run);
0200 }
0201
0202 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0203
0204 void Run::EndOfRun()
0205 {
0206   G4int prec = 5, wid = prec + 2;
0207   G4int dfprec = G4cout.precision(prec);
0208
0209   // run condition
0210   //
0211   G4Material* material = fDetector->GetAbsorMaterial();
0212   G4String Particle = fParticle->GetParticleName();
0213   G4cout << "\n The run is " << numberOfEvent << " " << Particle << " of "
0214          << G4BestUnit(fEkin, "Energy") << " within " << material->GetName()
0215          << " (D =  " << G4BestUnit(2 * (fDetector->GetAbsorRadius()), "Length")
0216          << " L = " << G4BestUnit(fDetector->GetAbsorLength(), "Length") << ")" << G4endl;
0217
0218   if (numberOfEvent == 0) {
0219     G4cout.precision(dfprec);
0220     return;
0221   }
0222
0223   // frequency of processes
0224   //
0225   G4cout << "\n Process calls frequency :" << G4endl;
0226   G4int index = 0;
0227   std::map<G4String, G4int>::iterator it;
0228   for (it = fProcCounter.begin(); it != fProcCounter.end(); it++) {
0229     G4String procName = it->first;
0230     G4int count = it->second;
0231     G4String space = " ";
0232     if (++index % 3 == 0) space = "\n";
0233     G4cout << " " << std::setw(20) << procName << "=" << std::setw(7) << count << space;
0234   }
0235   G4cout << G4endl;
0236
0237   // particles count
0238   //
0239   G4cout << "\n List of generated particles (with meanLife != 0) :" << G4endl;
0240
0241   std::map<G4String, ParticleData>::iterator itc;
0242   for (itc = fParticleDataMap1.begin(); itc != fParticleDataMap1.end(); itc++) {
0243     G4String name = itc->first;
0244     ParticleData data = itc->second;
0245     G4int count = data.fCount;
0246     G4double eMean = data.fEmean / count;
0247     G4double eMin = data.fEmin;
0248     G4double eMax = data.fEmax;
0249     G4double meanLife = data.fTmean;
0250
0251     G4cout << "  " << std::setw(13) << name << ": " << std::setw(7) << count
0252            << "  Emean = " << std::setw(wid) << G4BestUnit(eMean, "Energy") << "\t( "
0253            << G4BestUnit(eMin, "Energy") << " --> " << G4BestUnit(eMax, "Energy") << ")";
0254     if (meanLife >= 0.)
0255       G4cout << "\tmean life = " << G4BestUnit(meanLife, "Time") << G4endl;
0256     else
0257       G4cout << "\tstable" << G4endl;
0258   }
0259
0260   // compute mean Energy deposited and rms
0261   //
0262   G4int TotNbofEvents = numberOfEvent;
0263   fEnergyDeposit /= TotNbofEvents;
0264   fEnergyDeposit2 /= TotNbofEvents;
0265   G4double rmsEdep = fEnergyDeposit2 - fEnergyDeposit * fEnergyDeposit;
0266   if (rmsEdep > 0.)
0267     rmsEdep = std::sqrt(rmsEdep);
0268   else
0269     rmsEdep = 0.;
0270
0271   G4cout << "\n Mean energy deposit per event = " << G4BestUnit(fEnergyDeposit, "Energy")
0272          << ";  rms = " << G4BestUnit(rmsEdep, "Energy") << G4endl;
0273
0274   // compute mean Energy flow and rms
0275   //
0276   fEnergyFlow /= TotNbofEvents;
0277   fEnergyFlow2 /= TotNbofEvents;
0278   G4double rmsEflow = fEnergyFlow2 - fEnergyFlow * fEnergyFlow;
0279   if (rmsEflow > 0.)
0280     rmsEflow = std::sqrt(rmsEflow);
0281   else
0282     rmsEflow = 0.;
0283
0284   G4cout << " Mean energy flow per event    = " << G4BestUnit(fEnergyFlow, "Energy")
0285          << ";  rms = " << G4BestUnit(rmsEflow, "Energy") << G4endl;
0286
0287   // particles flux
0288   //
0289   G4cout << "\n List of particles emerging from the container :" << G4endl;
0290
0291   std::map<G4String, ParticleData>::iterator itn;
0292   for (itn = fParticleDataMap2.begin(); itn != fParticleDataMap2.end(); itn++) {
0293     G4String name = itn->first;
0294     ParticleData data = itn->second;
0295     G4int count = data.fCount;
0296     G4double eMean = data.fEmean / count;
0297     G4double eMin = data.fEmin;
0298     G4double eMax = data.fEmax;
0299     G4double Eflow = data.fEmean / TotNbofEvents;
0300
0301     G4cout << "  " << std::setw(13) << name << ": " << std::setw(7) << count
0302            << "  Emean = " << std::setw(wid) << G4BestUnit(eMean, "Energy") << "\t( "
0303            << G4BestUnit(eMin, "Energy") << " --> " << G4BestUnit(eMax, "Energy")
0304            << ") \tEflow/event = " << G4BestUnit(Eflow, "Energy") << G4endl;
0305   }
0306
0307   // remove all contents in fProcCounter, fCount
0308   fProcCounter.clear();
0309   fParticleDataMap2.clear();
0310
0311   // restore default format
0312   G4cout.precision(dfprec);
0313 }
0314
0315 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......