Activation/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 oooOO0OOooo::Run class
0028
0029 #include "Run.hh"
0030
0031 #include "DetectorConstruction.hh"
0032 #include "HistoManager.hh"
0033 #include "PrimaryGeneratorAction.hh"
0034
0035 #include "G4AutoLock.hh"
0036 #include "G4SystemOfUnits.hh"
0037 #include "G4Threading.hh"
0038 #include "G4UnitsTable.hh"
0039
0040 // mutex in a file scope
0041
0042 namespace
0043 {
0044 // Mutex to lock updating the global ion map
0045 G4Mutex ionIdMapMutex = G4MUTEX_INITIALIZER;
0046 }  // namespace
0047
0048 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0049
0050 std::map<G4String, G4int> Run::fgIonMap;
0051 G4int Run::fgIonId = kMaxHisto1;
0052
0053 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0054
0055 Run::Run(DetectorConstruction* det) : fDetector(det) {}
0056
0057 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0058
0059 void Run::Merge(std::map<G4String, ParticleData>& destinationMap,
0060                 const std::map<G4String, ParticleData>& sourceMap) const
0061 {
0062   for (const auto& particleData : sourceMap) {
0063     G4String name = particleData.first;
0064     const ParticleData& localData = particleData.second;
0065     if (destinationMap.find(name) == destinationMap.end()) {
0066       destinationMap[name] = ParticleData(localData.fCount, localData.fEmean, localData.fEmin,
0067                                           localData.fEmax, localData.fTmean);
0068     }
0069     else {
0070       ParticleData& data = destinationMap[name];
0071       data.fCount += localData.fCount;
0072       data.fEmean += localData.fEmean;
0073       G4double emin = localData.fEmin;
0074       if (emin < data.fEmin) data.fEmin = emin;
0075       G4double emax = localData.fEmax;
0076       if (emax > data.fEmax) data.fEmax = emax;
0077       data.fTmean = localData.fTmean;
0078     }
0079   }
0080 }
0081
0082 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0083
0084 void Run::SetPrimary(G4ParticleDefinition* particle, G4double energy)
0085 {
0086   fParticle = particle;
0087   fEkin = energy;
0088 }
0089
0090 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0091
0092 void Run::CountProcesses(const G4VProcess* process)
0093 {
0094   if (process == nullptr) return;
0095   G4String procName = process->GetProcessName();
0096   std::map<G4String, G4int>::iterator it = fProcCounter.find(procName);
0097   if (it == fProcCounter.end()) {
0098     fProcCounter[procName] = 1;
0099   }
0100   else {
0101     fProcCounter[procName]++;
0102   }
0103 }
0104
0105 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0106
0107 void Run::ParticleCount(G4String name, G4double Ekin, G4double meanLife)
0108 {
0109   std::map<G4String, ParticleData>::iterator it = fParticleDataMap1.find(name);
0110   if (it == fParticleDataMap1.end()) {
0111     fParticleDataMap1[name] = ParticleData(1, Ekin, Ekin, Ekin, meanLife);
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 = meanLife;
0123   }
0124 }
0125
0126 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0127
0128 void Run::AddEdep(G4double edep)
0129 {
0130   fEnergyDeposit += edep;
0131   fEnergyDeposit2 += edep * edep;
0132 }
0133
0134 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0135
0136 void Run::AddEflow(G4double eflow)
0137 {
0138   fEnergyFlow += eflow;
0139   fEnergyFlow2 += eflow * eflow;
0140 }
0141 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0142
0143 void Run::ParticleFlux(G4String name, G4double Ekin)
0144 {
0145   std::map<G4String, ParticleData>::iterator it = fParticleDataMap2.find(name);
0146   if (it == fParticleDataMap2.end()) {
0147     fParticleDataMap2[name] = ParticleData(1, Ekin, Ekin, Ekin, -1 * ns);
0148   }
0149   else {
0150     ParticleData& data = it->second;
0151     data.fCount++;
0152     data.fEmean += Ekin;
0153     // update min max
0154     G4double emin = data.fEmin;
0155     if (Ekin < emin) data.fEmin = Ekin;
0156     G4double emax = data.fEmax;
0157     if (Ekin > emax) data.fEmax = Ekin;
0158     data.fTmean = -1 * ns;
0159   }
0160 }
0161
0162 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0163
0164 G4int Run::GetIonId(G4String ionName)
0165 {
0166   G4AutoLock lock(&ionIdMapMutex);
0167   // updating the global ion map needs to be locked
0168
0169   std::map<G4String, G4int>::const_iterator it = fgIonMap.find(ionName);
0170   if (it == fgIonMap.end()) {
0171     fgIonMap[ionName] = fgIonId;
0172     if (fgIonId < (kMaxHisto2 - 1)) fgIonId++;
0173   }
0174   return fgIonMap[ionName];
0175 }
0176
0177 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0178
0179 void Run::Merge(const G4Run* run)
0180 {
0181   const Run* localRun = static_cast<const Run*>(run);
0182
0183   // primary particle info
0184   //
0185   fParticle = localRun->fParticle;
0186   fEkin = localRun->fEkin;
0187
0188   // accumulate sums
0189   //
0190   fEnergyDeposit += localRun->fEnergyDeposit;
0191   fEnergyDeposit2 += localRun->fEnergyDeposit2;
0192   fEnergyFlow += localRun->fEnergyFlow;
0193   fEnergyFlow2 += localRun->fEnergyFlow2;
0194
0195   // map: processes count
0196   for (const auto& procCounter : localRun->fProcCounter) {
0197     G4String procName = procCounter.first;
0198     G4int localCount = procCounter.second;
0199     if (fProcCounter.find(procName) == fProcCounter.end()) {
0200       fProcCounter[procName] = localCount;
0201     }
0202     else {
0203       fProcCounter[procName] += localCount;
0204     }
0205   }
0206
0207   // map: created particles count
0208   Merge(fParticleDataMap1, localRun->fParticleDataMap1);
0209
0210   // map: particles flux count
0211   Merge(fParticleDataMap2, localRun->fParticleDataMap2);
0212
0213   G4Run::Merge(run);
0214 }
0215
0216 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0217
0218 void Run::EndOfRun()
0219 {
0220   G4int prec = 5, wid = prec + 2;
0221   G4int dfprec = G4cout.precision(prec);
0222
0223   // run condition
0224   //
0225   G4Material* material = fDetector->GetAbsorMaterial();
0226   G4double density = material->GetDensity();
0227
0228   G4String Particle = fParticle->GetParticleName();
0229   G4cout << "\n The run is " << numberOfEvent << " " << Particle << " of "
0230          << G4BestUnit(fEkin, "Energy") << " through "
0231          << G4BestUnit(fDetector->GetAbsorThickness(), "Length") << " of " << material->GetName()
0232          << " (density: " << G4BestUnit(density, "Volumic Mass") << ")" << G4endl;
0233
0234   if (numberOfEvent == 0) {
0235     G4cout.precision(dfprec);
0236     return;
0237   }
0238
0239   // frequency of processes
0240   //
0241   G4cout << "\n Process calls frequency :" << G4endl;
0242   G4int index = 0;
0243   for (const auto& procCounter : fProcCounter) {
0244     G4String procName = procCounter.first;
0245     G4int count = procCounter.second;
0246     G4String space = " ";
0247     if (++index % 3 == 0) space = "\n";
0248     G4cout << " " << std::setw(20) << procName << "=" << std::setw(7) << count << space;
0249   }
0250   G4cout << G4endl;
0251
0252   // particles count
0253   //
0254   G4cout << "\n List of generated particles (with meanLife != 0):" << G4endl;
0255
0256   for (const auto& particleData : fParticleDataMap1) {
0257     G4String name = particleData.first;
0258     ParticleData data = particleData.second;
0259     G4int count = data.fCount;
0260     G4double eMean = data.fEmean / count;
0261     G4double eMin = data.fEmin;
0262     G4double eMax = data.fEmax;
0263     G4double meanLife = data.fTmean;
0264
0265     G4cout << "  " << std::setw(13) << name << ": " << std::setw(7) << count
0266            << "  Emean = " << std::setw(wid) << G4BestUnit(eMean, "Energy") << "\t( "
0267            << G4BestUnit(eMin, "Energy") << " --> " << G4BestUnit(eMax, "Energy") << ")";
0268     if (meanLife >= 0.)
0269       G4cout << "\tmean life = " << G4BestUnit(meanLife, "Time") << G4endl;
0270     else
0271       G4cout << "\tstable" << G4endl;
0272   }
0273
0274   // compute mean Energy deposited and rms
0275   //
0276   G4int TotNbofEvents = numberOfEvent;
0277   fEnergyDeposit /= TotNbofEvents;
0278   fEnergyDeposit2 /= TotNbofEvents;
0279   G4double rmsEdep = fEnergyDeposit2 - fEnergyDeposit * fEnergyDeposit;
0280   if (rmsEdep > 0.)
0281     rmsEdep = std::sqrt(rmsEdep);
0282   else
0283     rmsEdep = 0.;
0284
0285   G4cout << "\n Mean energy deposit per event = " << G4BestUnit(fEnergyDeposit, "Energy")
0286          << ";  rms = " << G4BestUnit(rmsEdep, "Energy") << G4endl;
0287
0288   // compute mean Energy flow and rms
0289   //
0290   fEnergyFlow /= TotNbofEvents;
0291   fEnergyFlow2 /= TotNbofEvents;
0292   G4double rmsEflow = fEnergyFlow2 - fEnergyFlow * fEnergyFlow;
0293   if (rmsEflow > 0.)
0294     rmsEflow = std::sqrt(rmsEflow);
0295   else
0296     rmsEflow = 0.;
0297
0298   G4cout << " Mean energy flow per event    = " << G4BestUnit(fEnergyFlow, "Energy")
0299          << ";  rms = " << G4BestUnit(rmsEflow, "Energy") << G4endl;
0300
0301   // particles flux
0302   //
0303   G4cout << "\n List of particles emerging from the target :" << G4endl;
0304
0305   for (const auto& particleData : fParticleDataMap2) {
0306     G4String name = particleData.first;
0307     ParticleData data = particleData.second;
0308     G4int count = data.fCount;
0309     G4double eMean = data.fEmean / count;
0310     G4double eMin = data.fEmin;
0311     G4double eMax = data.fEmax;
0312     G4double Eflow = data.fEmean / TotNbofEvents;
0313
0314     G4cout << "  " << std::setw(13) << name << ": " << std::setw(7) << count
0315            << "  Emean = " << std::setw(wid) << G4BestUnit(eMean, "Energy") << "\t( "
0316            << G4BestUnit(eMin, "Energy") << " --> " << G4BestUnit(eMax, "Energy")
0317            << ") \tEflow/event = " << G4BestUnit(Eflow, "Energy") << G4endl;
0318   }
0319
0320   // histogram Id for populations
0321   //
0322   G4cout << "\n histo Id for populations :" << G4endl;
0323
0324   // Update the histogram titles according to the ion map
0325   // and print new titles
0326   G4AnalysisManager* analysisManager = G4AnalysisManager::Instance();
0327   for (const auto& ionMapElement : fgIonMap) {
0328     G4String ionName = ionMapElement.first;
0329     G4int h1Id = ionMapElement.second;
0330     // print new titles
0331     G4cout << " " << std::setw(20) << ionName << "  id = " << std::setw(3) << h1Id << G4endl;
0332
0333     // update histogram ids
0334     if (!analysisManager->GetH1(h1Id)) continue;
0335     // Skip inactive histograms, this is not necessary
0336     // but it  makes the code safe wrt modifications in future
0337     G4String title = analysisManager->GetH1Title(h1Id);
0338     title = ionName + title;
0339     analysisManager->SetH1Title(h1Id, title);
0340   }
0341   G4cout << G4endl;
0342
0343   // normalize histograms
0344   G4int ih = 2;
0345   G4double binWidth = analysisManager->GetH1Width(ih);
0346   G4double fac = (1. / (numberOfEvent * binWidth)) * (mm / MeV);
0347   analysisManager->ScaleH1(ih, fac);
0348
0349   for (ih = 14; ih < 24; ih++) {
0350     binWidth = analysisManager->GetH1Width(ih);
0351     G4double unit = analysisManager->GetH1Unit(ih);
0352     fac = (second / (binWidth * unit));
0353     analysisManager->ScaleH1(ih, fac);
0354   }
0355
0356   // remove all contents in fProcCounter, fCount
0357   fProcCounter.clear();
0358   fParticleDataMap1.clear();
0359   fParticleDataMap2.clear();
0360   fgIonMap.clear();
0361
0362   // restore default format
0363   G4cout.precision(dfprec);
0364 }
0365
0366 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......