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