FinalState/src/FinalStateHistoManager.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 FinalStateHistoManager.cc
0027 /// \brief Implementation of the FinalStateHistoManager class
0028 ///
0029 ///  Create a set of histos for final state study.
0030
0031 //
0032 //  Author: G.Hugo, 08 December 2022
0033 //
0034 // ***************************************************************************
0035 //
0036 //      FinalStateHistoManager
0037 //
0038 ///  Create a set of histos for final state study.
0039 ///  In practice, the interactions studied here are hadron nuclear inelastic interactions
0040 ///  (though the code is fully generic).
0041 ///
0042 ///  Energy spectra are plotted for all encountered secondaries
0043 ///  (one histo per secondary).
0044 ///  In addition, the residual nuclei Z and A distributions are plotted.
0045 ///
0046 ///  All histograms are G4H1.
0047 ///  They are created and filled solely via G4VAnalysisManager.
0048 ///
0049 ///  The histograms can be dumped to all usual formats, including ROOT
0050 ///  (via G4VAnalysisManager).
0051 ///  An interesting added feature here, is that the plots, while being allocated
0052 ///  and filled via G4VAnalysisManager, are also dumped
0053 ///  in a Flair-compatible format (via tools::histo::flair).
0054 ///
0055 ///  NB 1: Note that instead of a hardcoded number associated to a hardcoded set of particles,
0056 ///  particle PDG IDs are used to index the histos.
0057 ///  This allows a dynamic storage of all particles encountered in the final states.
0058 ///
0059 ///  NB 2: tools::histo::flair code, which allows the dump of any G4H1
0060 ///  into Flair-compatible format, is fully application-agnostic,
0061 ///  and is placed in FlukaCern/utils.
0062 ///  It could also be added as an extension of core G4 Analysis Manager.
0063 //
0064 // ***************************************************************************
0065
0066 #include "FinalStateHistoManager.hh"
0067
0068 #include "G4RootAnalysisManager.hh"
0069 // #include "G4AnalysisManager.hh"
0070
0071 #include "tools_histo_flair.hh"
0072
0073 #include "G4DynamicParticle.hh"
0074 #include "G4Exception.hh"
0075 #include "G4ParticleTable.hh"
0076 #include "G4ios.hh"
0077 #include "g4hntools_defs.hh"
0078
0079 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
0080
0081 FinalStateHistoManager::FinalStateHistoManager()
0082   : fOutputFileName("all_secondaries"),
0083     fRootOutputFileName(fOutputFileName + ".root"),
0084     fFlairOutputFileName(fOutputFileName + ".hist"),
0085     fNumBins(90),
0086     fMinKineticEnergy(10. * keV),
0087     fMaxKineticEnergy(10. * TeV),
0088     fFunctionName("none"),
0089     fBinSchemeName("log"),
0090     fRootEnergyUnit("MeV"),
0091     fNucleiZMax(25),
0092     fNucleiAMax(50),
0093     fNumEvents(0),
0094     fAnalysisManager(G4RootAnalysisManager::Instance()),
0095     fNucleiZScoreIndex(0),
0096     fNucleiAScoreIndex(1)
0097 {
0098   // fAnalysisManager = G4AnalysisManager::Instance();
0099   // fAnalysisManager->SetDefaultFileType("root");
0100   // fAnalysisManager->SetVerboseLevel(0);
0101   // fOutputFileName += fAnalysisManager->GetFileType();
0102 }
0103
0104 // ***************************************************************************
0105 // Open output file + create residual nuclei histograms considered for final state study.
0106 // The histograms are G4H1, created via G4VAnalysisManager.
0107 // ***************************************************************************
0108 void FinalStateHistoManager::Book()
0109 {
0110   // Open file.
0111   if (!fAnalysisManager->OpenFile(fRootOutputFileName)) {
0112     G4ExceptionDescription msg;
0113     msg << "Booking histograms: cannot open file " << fRootOutputFileName << G4endl;
0114     G4Exception("FinalStateHistoManager::Book", "Cannot open file", FatalException, msg);
0115   }
0116   G4cout << "### FinalStateHistoManager::Book: Successfully opended file " << fRootOutputFileName
0117          << " for dumping histograms." << G4endl;
0118
0119   // Create the residual nuclei distributions (in Z and A).
0120   const G4int nucleiZHistoIndex = fAnalysisManager->CreateH1(
0121     "nucleiZ", "Residual nuclei distribution in Z", fNucleiZMax, 0.5, fNucleiZMax + 0.5);
0122   auto nucleiZHistoWrapper = std::make_unique<G4H1Wrapper>(fAnalysisManager, nucleiZHistoIndex);
0123   fNucleiData.insert(std::make_pair(fNucleiZScoreIndex, std::move(nucleiZHistoWrapper)));
0124
0125   const G4int nucleiAHistoIndex = fAnalysisManager->CreateH1(
0126     "nucleiA", "Residual nuclei distribution in A", fNucleiAMax, 0.5, fNucleiAMax + 0.5);
0127   auto nucleiAHistoWrapper = std::make_unique<G4H1Wrapper>(fAnalysisManager, nucleiAHistoIndex);
0128   fNucleiData.insert(std::make_pair(fNucleiAScoreIndex, std::move(nucleiAHistoWrapper)));
0129 }
0130
0131 // ***************************************************************************
0132 // Keep track of the total number of events (used later on for normalization).
0133 // ***************************************************************************
0134 void FinalStateHistoManager::BeginOfEvent()
0135 {
0136   fNumEvents++;
0137 }
0138
0139 // ***************************************************************************
0140 // Fill all plots (WITHIN event, ie the interaction).
0141 // ***************************************************************************
0142 void FinalStateHistoManager::ScoreSecondary(const G4DynamicParticle* const secondary)
0143 {
0144   // SELECT SPECIFIC SECONDARIES ONLY
0145   // Select by angle with beam direction
0146   /* if ( (std::pow(secondary->GetMomentumDirection().x(), 2.)
0147      + std::pow(secondary->GetMomentumDirection().y(), 2.))
0148      <= 0.0001 ) {*/
0149
0150   // Select by production tag
0151   /* if (secondary->GetProductionTag() == 6) {*/
0152
0153   // Primary track
0154   /* if(track->GetParentID() == 0) {*/
0155
0156   const auto& particle = secondary->GetDefinition();
0157
0158   // SECONDARIES ENERGY SPECTRA
0159   // Dynamic creation of histos, so that all encountered particles have their own histos.
0160
0161   // Check whether a particle has already been encountered.
0162   const auto found = fParticleData.find(secondary->GetPDGcode());
0163
0164   G4H1Wrapper* particleHistoWrapper = nullptr;
0165   // If the particle has already been encountered, use the corresponding histos.
0166   if (found != fParticleData.end()) {
0167     particleHistoWrapper = found->second.get();
0168   }
0169   // Otherwise, create histos for that particle.
0170   else {
0171     const G4String& particleName = particle->GetParticleName();
0172     const G4int particlePDG = secondary->GetPDGcode();
0173
0174     const G4String histoTitle = (particlePDG == 0 ? G4String("Particle pdg==0 spectrum")
0175                                                   : particleName + G4String(" spectrum"));
0176     const G4int histoIndex =
0177       fAnalysisManager->CreateH1(particleName, histoTitle, fNumBins, fMinKineticEnergy,
0178                                  fMaxKineticEnergy, fRootEnergyUnit, fFunctionName, fBinSchemeName);
0179     auto histoWrapper = std::make_unique<G4H1Wrapper>(fAnalysisManager, histoIndex);
0180     particleHistoWrapper = histoWrapper.get();
0181     fParticleData.insert(std::make_pair(particlePDG, std::move(histoWrapper)));
0182   }
0183
0184   // Fill the G4H1Wrapper.
0185   const G4double kineticEnergy = secondary->GetKineticEnergy();
0186   particleHistoWrapper->Fill(kineticEnergy, 1.);
0187
0188   // NUCLEI DISTRIBUTIONS IN Z AND A
0189   if (particle->GetParticleType() == "nucleus") {
0190     // Fill the G4H1Wrapper.
0191     const G4double Z = particle->GetPDGCharge() / eplus;
0192     fNucleiData[fNucleiZScoreIndex]->Fill(Z, 1.);
0193
0194     // Fill the G4H1Wrapper.
0195     const G4double A = particle->GetBaryonNumber();
0196     fNucleiData[fNucleiAScoreIndex]->Fill(A, 1.);
0197   }
0198
0199   //} // select secondaries
0200 }
0201
0202 // ***************************************************************************
0203 // End of event: all event-level G4H1 are flushed into the Analysis Manager G4H1.
0204 // ***************************************************************************
0205 void FinalStateHistoManager::EndOfEvent()
0206 {
0207   for (const auto& particleIt : fParticleData) {
0208     particleIt.second->EndOfEvent();
0209   }
0210   for (const auto& nucleiScoreIt : fNucleiData) {
0211     nucleiScoreIt.second->EndOfEvent();
0212   }
0213 }
0214
0215 // ***************************************************************************
0216 // Printout secondary counts + dump all plots into relevant formats.
0217 // ***************************************************************************
0218 void FinalStateHistoManager::EndOfRun() const
0219 {
0220   // PRINTOUT SECONDARYS COUNTS (FULL ENERGY RANGE).
0221
0222   // Order the histos by particles names.
0223   std::map<G4String, const G4H1Wrapper*> particlesHistos;
0224
0225   for (const auto& particleIt : fParticleData) {
0226     const G4int particlePdg = particleIt.first;
0227     const G4String particleName =
0228       G4ParticleTable::GetParticleTable()->FindParticle(particlePdg)->GetParticleName();
0229
0230     const G4H1Wrapper* const particleHisto = particleIt.second.get();
0231     particlesHistos.insert(std::make_pair(particleName, particleHisto));
0232   }
0233
0234   // Printout secondarys counts (full energy range)
0235   // Values are averaged over the number of events.
0236   G4cout << "========================================================" << G4endl;
0237   G4cout << "Number of events                     " << fNumEvents << G4endl << G4endl;
0238   for (const auto& particleIt : particlesHistos) {
0239     // Note that the info is directly obtained from the histogram:
0240     // it is the integral over the full energy range.
0241     const G4int count = particleIt.second->GetG4H1()->sum_all_bin_heights();
0242
0243     const G4double averageCount = static_cast<G4double>(count) / fNumEvents;
0244
0245     G4cout << "Average (per event) number of " << particleIt.first << "              "
0246            << averageCount << G4endl;
0247   }
0248   G4cout << "========================================================" << G4endl;
0249   G4cout << G4endl;
0250
0251   // DUMP G4H1 PLOTS INTO ROOT FILE
0252   DumpAllG4H1IntoRootFile();
0253
0254   // DUMP G4H1 PLOTS INTO FLAIR FILE
0255   DumpAllG4H1IntoFlairFile(particlesHistos);
0256
0257   // Close and clear fAnalysisManager.
0258   fAnalysisManager->CloseFile();
0259   fAnalysisManager->Clear();
0260 }
0261
0262 // ***************************************************************************
0263 // DUMP G4H1 PLOTS INTO ROOT FILE (via G4VAnalysisManager).
0264 // ***************************************************************************
0265 void FinalStateHistoManager::DumpAllG4H1IntoRootFile() const
0266 {
0267   if (!fAnalysisManager->Write()) {
0268     G4ExceptionDescription message;
0269     message << "Could not write ROOT file.";
0270     G4Exception("FinalStateHistoManager::EndOfRun()", "I/O Error", FatalException, message);
0271   }
0272   G4cout << "### All histograms saved to " << fRootOutputFileName << G4endl;
0273 }
0274
0275 // ***************************************************************************
0276 // DUMP G4H1 PLOTS INTO FLAIR FILE (via tools::histo::flair).
0277 // ***************************************************************************
0278 void FinalStateHistoManager::DumpAllG4H1IntoFlairFile(
0279   const std::map<G4String, const G4H1Wrapper*>& particlesHistos) const
0280 {
0281   std::ofstream output;
0282   output.open(fFlairOutputFileName, std::ios_base::out);
0283   G4int indexInOutputFile = 1;
0284
0285   // SECONDARIES ENERGY SPECTRA
0286   for (const auto& particleIt : particlesHistos) {
0287     const G4String& histoName = particleIt.first;
0288     const auto& histo = particleIt.second->GetG4H1();
0289
0290     tools::histo::flair::dumpG4H1HistoInFlairFormat(
0291       output, indexInOutputFile, histoName, histo, tools::histo::flair::Abscissa::KineticEnergy,
0292       fBinSchemeName, fNumEvents, particleIt.second->GetSumSquaredEventTotals(),
0293       particleIt.second->GetSumSquaredEventInRangeTotals());
0294     ++indexInOutputFile;
0295   }
0296
0297   // RESIDUAL NUCLEI DISTRIBUTIONS
0298   for (const auto& plotIt : fNucleiData) {
0299     const auto& histo = plotIt.second->GetG4H1();
0300     const G4String& histoName = (plotIt.first == fNucleiZScoreIndex ? "nucleiZ" : "nucleiA");
0301     const auto& abscissaKind =
0302       (plotIt.first == fNucleiZScoreIndex ? tools::histo::flair::Abscissa::Z
0303                                           : tools::histo::flair::Abscissa::A);
0304
0305     tools::histo::flair::dumpG4H1HistoInFlairFormat(
0306       output, indexInOutputFile, histoName, histo, abscissaKind, fBinSchemeName, fNumEvents,
0307       plotIt.second->GetSumSquaredEventTotals(), plotIt.second->GetSumSquaredEventInRangeTotals());
0308     ++indexInOutputFile;
0309   }
0310
0311   output.close();
0312   G4cout << "### All histograms saved to " << fFlairOutputFileName << G4endl;
0313 }
0314
0315 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......