CrossSection/src/XSHistoManager.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 XSHistoManager.cc
0027 /// \brief Implementation of the XSHistoManager class
0028
0029 //
0030 //  Adapted from hadronic/Hadr00/src/HistoManager.cc
0031 //  Author: G.Hugo, 06 January 2023
0032 //
0033 // ***************************************************************************
0034 //
0035 //      XSHistoManager
0036 //
0037 ///  Create a set of profiles for XS study.
0038 ///
0039 ///  All profiles are G4H1.
0040 ///  They are created and filled via G4VAnalysisManager.
0041 ///
0042 ///  The profiles can be dumped to all usual formats, including ROOT
0043 ///  (via G4VAnalysisManager).
0044 ///  They are also dumped in a format compatible with Flair
0045 ///  (via tools::histo::flair).
0046 //
0047 // ***************************************************************************
0048
0049 #include "XSHistoManager.hh"
0050
0051 #include "G4Element.hh"
0052 #include "G4HadronicProcessStore.hh"
0053 #include "G4Material.hh"
0054 #include "G4NistManager.hh"
0055 #include "G4ParticleDefinition.hh"
0056 #include "G4ParticleTable.hh"
0057 #include "G4ios.hh"
0058
0059 // #include "G4AnalysisManager.hh"
0060 #include "tools_histo_flair.hh"
0061
0062 #include "G4RootAnalysisManager.hh"
0063
0064 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0065
0066 XSHistoManager::XSHistoManager()
0067   : fMessenger(new XSHistoManagerMessenger(this)),
0068     fOutputFileName("all_XS"),
0069     fRootOutputFileName("all_XS.root"),
0070     fFlairOutputFileName("all_XS.hist"),
0071     fParticle(nullptr),
0072     fElement(nullptr),
0073     fMaterial(nullptr),
0074     fNumBins(10000),
0075     fMinKineticEnergy(1. * keV),
0076     fMaxKineticEnergy(10. * TeV),
0077     fFunctionName("none"),
0078     fBinSchemeName("log"),
0079     fRootEnergyUnit("MeV"),
0080     fAnalysisManager(G4RootAnalysisManager::Instance()),
0081     fElasticXSIndex(0),
0082     fInelasticXSIndex(1),
0083     fCaptureXSIndex(2),
0084     fFissionXSIndex(3),
0085     fChargeExchangeXSIndex(4),
0086     fTotalXSIndex(5),
0087     fElasticPerVolumeXSIndex(6),
0088     fInelasticPerVolumeXSIndex(7)
0089 {
0090   // G4NistManager::Instance()->ListMaterials("all");
0091 }
0092
0093 // ***************************************************************************
0094 // Set output files names: 2 formats supported, ROOT and Flair.
0095 // ***************************************************************************
0096 void XSHistoManager::SetOutputFileName(const G4String& outputFileName)
0097 {
0098   fOutputFileName = outputFileName;
0099   fRootOutputFileName = outputFileName + ".root";
0100   fFlairOutputFileName = outputFileName + ".hist";
0101 }
0102
0103 // ***************************************************************************
0104 // Set the particle considered for XS study.
0105 // ***************************************************************************
0106 void XSHistoManager::SetParticle(const G4String& particleName)
0107 {
0108   fParticle = G4ParticleTable::GetParticleTable()->FindParticle(particleName);
0109 }
0110
0111 // ***************************************************************************
0112 // Set the target element considered for XS study.
0113 // ***************************************************************************
0114 void XSHistoManager::SetElement(const G4String& elementName)
0115 {
0116   fElement = G4NistManager::Instance()->FindOrBuildElement(elementName);
0117   // Also needs to set material!
0118   SetMaterial(elementName);
0119 }
0120
0121 // ***************************************************************************
0122 // Set the target material considered for XS study.
0123 // ***************************************************************************
0124 void XSHistoManager::SetMaterial(const G4String& materialName)
0125 {
0126   // Check that material is not set already.
0127   if (fMaterial) {
0128     G4ExceptionDescription msg;
0129     msg << "Please use UI command /allXS/elementName"
0130         << " OR UI command /allXS/nonElementaryMaterialName,"
0131         << " BUT NOT BOTH!" << G4endl;
0132     G4Exception("XSHistoManager::SetMaterial", "A target material is already defined.",
0133                 FatalException, msg);
0134   }
0135
0136   fMaterial = G4NistManager::Instance()->FindOrBuildMaterial("G4_" + materialName);
0137 }
0138
0139 // ***************************************************************************
0140 // Open output file + create all profiles considered for XS study.
0141 // All profiles are G4H1, created via G4VAnalysisManager.
0142 // ***************************************************************************
0143 void XSHistoManager::Book()
0144 {
0145   // Check all XSHistoManager data is set properly.
0146   CheckInput();
0147
0148   // Open file.
0149   if (!fAnalysisManager->OpenFile(fRootOutputFileName)) {
0150     G4ExceptionDescription msg;
0151     msg << "Booking profiles: cannot open file " << fRootOutputFileName << G4endl;
0152     G4Exception("XSHistoManager::Book", "Cannot open file", FatalException, msg);
0153   }
0154   G4cout << "### XSHistoManager::Book: Successfully opened file " << fRootOutputFileName
0155          << " for dumping profiles." << G4endl;
0156
0157   // Create all G4H1, and keep track of each histo index in fXSProfileIndex.
0158   const G4int elasticXSProfileIndex =
0159     fAnalysisManager->CreateH1("ElasticXS", "Elastic XS", fNumBins, fMinKineticEnergy,
0160                                fMaxKineticEnergy, fRootEnergyUnit, fFunctionName, fBinSchemeName);
0161   fXSProfileIndex.insert(std::make_pair(fElasticXSIndex, elasticXSProfileIndex));
0162
0163   const G4int inelasticXSProfileIndex =
0164     fAnalysisManager->CreateH1("InelasticXS", "Inelastic XS", fNumBins, fMinKineticEnergy,
0165                                fMaxKineticEnergy, fRootEnergyUnit, fFunctionName, fBinSchemeName);
0166   fXSProfileIndex.insert(std::make_pair(fInelasticXSIndex, inelasticXSProfileIndex));
0167
0168   const G4int captureXSProfileIndex =
0169     fAnalysisManager->CreateH1("CaptureXS", "Capture XS", fNumBins, fMinKineticEnergy,
0170                                fMaxKineticEnergy, fRootEnergyUnit, fFunctionName, fBinSchemeName);
0171   fXSProfileIndex.insert(std::make_pair(fCaptureXSIndex, captureXSProfileIndex));
0172
0173   const G4int fissionXSProfileIndex =
0174     fAnalysisManager->CreateH1("FissionXS", "Fission XS", fNumBins, fMinKineticEnergy,
0175                                fMaxKineticEnergy, fRootEnergyUnit, fFunctionName, fBinSchemeName);
0176   fXSProfileIndex.insert(std::make_pair(fFissionXSIndex, fissionXSProfileIndex));
0177
0178   const G4int chargeExchangeXSProfileIndex = fAnalysisManager->CreateH1(
0179     "ChargeExchangeXS", "Charge exchange XS", fNumBins, fMinKineticEnergy, fMaxKineticEnergy,
0180     fRootEnergyUnit, fFunctionName, fBinSchemeName);
0181   fXSProfileIndex.insert(std::make_pair(fChargeExchangeXSIndex, chargeExchangeXSProfileIndex));
0182
0183   const G4int totalXSProfileIndex =
0184     fAnalysisManager->CreateH1("TotalXS", "Total XS", fNumBins, fMinKineticEnergy,
0185                                fMaxKineticEnergy, fRootEnergyUnit, fFunctionName, fBinSchemeName);
0186   fXSProfileIndex.insert(std::make_pair(fTotalXSIndex, totalXSProfileIndex));
0187
0188   const G4int elasticPerVolumeXSProfileIndex = fAnalysisManager->CreateH1(
0189     "ElasticPerVolumeXS", "Elastic XS per volume", fNumBins, fMinKineticEnergy, fMaxKineticEnergy,
0190     fRootEnergyUnit, fFunctionName, fBinSchemeName);
0191   fXSProfileIndex.insert(std::make_pair(fElasticPerVolumeXSIndex, elasticPerVolumeXSProfileIndex));
0192
0193   const G4int inelasticPerVolumeXSProfileIndex = fAnalysisManager->CreateH1(
0194     "InelasticPerVolumeXS", "Inelastic XS per volume", fNumBins, fMinKineticEnergy,
0195     fMaxKineticEnergy, fRootEnergyUnit, fFunctionName, fBinSchemeName);
0196   fXSProfileIndex.insert(
0197     std::make_pair(fInelasticPerVolumeXSIndex, inelasticPerVolumeXSProfileIndex));
0198 }
0199
0200 // ***************************************************************************
0201 // Fill all plots, then dump them into relevant formats.
0202 // ***************************************************************************
0203 void XSHistoManager::EndOfRun()
0204 {
0205   G4cout << "### XSHistoManager::EndOfRun: Compute & fill XS for " << fParticle->GetParticleName()
0206          << " in " << (fElement ? fElement->GetName() : fMaterial->GetName()) << G4endl;
0207
0208   G4HadronicProcessStore* const store = G4HadronicProcessStore::Instance();
0209
0210   // Fill XS profiles.
0211   const G4double logMinKineticEnergy = std::log10(fMinKineticEnergy);
0212   const G4double logMaxKineticEnergy = std::log10(fMaxKineticEnergy);
0213   const G4double deltaLogKineticEnergy = (logMaxKineticEnergy - logMinKineticEnergy) / fNumBins;
0214
0215   G4double logKineticEnergy = logMinKineticEnergy - deltaLogKineticEnergy / 2.;
0216
0217   // Loop on all kinetic energies of interest.
0218   for (G4int binIndex = 0; binIndex < fNumBins; ++binIndex) {
0219     logKineticEnergy += deltaLogKineticEnergy;
0220     const G4double kineticEnergy = std::pow(10., logKineticEnergy) * MeV;
0221
0222     G4double totalXS = 0.;
0223     if (fElement) {
0224       // ELASTIC (ELEMENTARY MATERIAL)
0225       const G4double elasticXS =
0226         store->GetElasticCrossSectionPerAtom(fParticle, kineticEnergy, fElement, fMaterial);
0227       fAnalysisManager->FillH1(fXSProfileIndex[fElasticXSIndex], kineticEnergy, elasticXS / barn);
0228       totalXS += elasticXS;
0229
0230       // INELASTIC (ELEMENTARY MATERIAL)
0231       const G4double inelasticXS =
0232         store->GetInelasticCrossSectionPerAtom(fParticle, kineticEnergy, fElement, fMaterial);
0233       fAnalysisManager->FillH1(fXSProfileIndex[fInelasticXSIndex], kineticEnergy,
0234                                inelasticXS / barn);
0235       totalXS += inelasticXS;
0236
0237       if (fParticle == G4Neutron::Definition()) {
0238         // NEUTRON CAPTURE (ELEMENTARY MATERIAL)
0239         const G4double captureXS =
0240           store->GetCaptureCrossSectionPerAtom(fParticle, kineticEnergy, fElement, fMaterial);
0241         fAnalysisManager->FillH1(fXSProfileIndex[fCaptureXSIndex], kineticEnergy, captureXS / barn);
0242         totalXS += captureXS;
0243
0244         // FISSION (ELEMENTARY MATERIAL)
0245         const G4double fissionXS =
0246           store->GetFissionCrossSectionPerAtom(fParticle, kineticEnergy, fElement, fMaterial);
0247         totalXS += fissionXS;
0248         fAnalysisManager->FillH1(fXSProfileIndex[fFissionXSIndex], kineticEnergy, fissionXS / barn);
0249       }
0250
0251       // CHARGE EXCHANGE (ELEMENTARY MATERIAL)
0252       const G4double chargeExchangeXS =
0253         store->GetChargeExchangeCrossSectionPerAtom(fParticle, kineticEnergy, fElement, fMaterial);
0254       fAnalysisManager->FillH1(fXSProfileIndex[fChargeExchangeXSIndex], kineticEnergy,
0255                                chargeExchangeXS / barn);
0256       totalXS += chargeExchangeXS;
0257
0258       // TOTAL (ELEMENTARY MATERIAL)
0259       fAnalysisManager->FillH1(fXSProfileIndex[fTotalXSIndex], kineticEnergy, totalXS / barn);
0260     }
0261
0262     if (fMaterial) {
0263       const G4double materialSurfacicDensity =
0264         (fMaterial ? fMaterial->GetDensity() / (g / cm2) : 1.);
0265
0266       // ELASTIC
0267       const G4double elasticPerVolumeXS =
0268         store->GetElasticCrossSectionPerVolume(fParticle, kineticEnergy, fMaterial);
0269       fAnalysisManager->FillH1(fXSProfileIndex[fElasticPerVolumeXSIndex], kineticEnergy,
0270                                elasticPerVolumeXS / materialSurfacicDensity);
0271
0272       // INELASTIC
0273       const G4double inelasticPerVolumeXS =
0274         store->GetInelasticCrossSectionPerVolume(fParticle, kineticEnergy, fMaterial);
0275       fAnalysisManager->FillH1(fXSProfileIndex[fInelasticPerVolumeXSIndex], kineticEnergy,
0276                                inelasticPerVolumeXS / materialSurfacicDensity);
0277     }
0278   }
0279
0280   // DUMP G4H1 PLOTS INTO ROOT FILE
0281   DumpAllG4H1IntoRootFile();
0282
0283   // DUMP G4H1 PLOTS INTO FLAIR FILE
0284   DumpAllG4H1IntoFlairFile();
0285
0286   // Close and clear fAnalysisManager.
0287   fAnalysisManager->CloseFile();
0288   fAnalysisManager->Clear();
0289 }
0290
0291 // ***************************************************************************
0292 // Checks that particle and material are set
0293 // (all others have relevant default values).
0294 // ***************************************************************************
0295 void XSHistoManager::CheckInput()
0296 {
0297   if (!fParticle) {
0298     G4ExceptionDescription msg;
0299     msg << "Please add a particle to study XS: UI command /allXS/particleName" << G4endl;
0300     G4Exception("XSHistoManager::CheckInput()", "Print XS: no input particle defined.",
0301                 FatalException, msg);
0302   }
0303
0304   if (!fMaterial) {
0305     G4ExceptionDescription msg;
0306     msg << "Please add a material to study XS:"
0307         << " UI command /allXS/elementName for an elementary material,"
0308         << " or UI command /allXS/nonElementaryMaterialName for a compound/mixture material."
0309         << G4endl;
0310     G4Exception("XSHistoManager::CheckInput()", "Print XS: no target material defined.",
0311                 FatalException, msg);
0312   }
0313 }
0314
0315 // ***************************************************************************
0316 // DUMP G4H1 PLOTS INTO ROOT FILE (via G4VAnalysisManager).
0317 // ***************************************************************************
0318 void XSHistoManager::DumpAllG4H1IntoRootFile() const
0319 {
0320   if (!fAnalysisManager->Write()) {
0321     G4ExceptionDescription message;
0322     message << "Could not write ROOT file.";
0323     G4Exception("XSHistoManager::EndOfRun()", "I/O Error", FatalException, message);
0324   }
0325   G4cout << "### All profiles saved to " << fRootOutputFileName << G4endl;
0326 }
0327
0328 // ***************************************************************************
0329 // DUMP G4H1 PLOTS INTO FLAIR FILE (via tools::histo::flair).
0330 // ***************************************************************************
0331 void XSHistoManager::DumpAllG4H1IntoFlairFile() const
0332 {
0333   std::ofstream output;
0334   output.open(fFlairOutputFileName, std::ios_base::out);
0335   auto const rootAnalysisManager = dynamic_cast<G4RootAnalysisManager*>(fAnalysisManager);
0336
0337   G4int indexInOutputFile = 1;
0338   for (G4int xsIndex = fElasticXSIndex; xsIndex <= fInelasticPerVolumeXSIndex; ++xsIndex) {
0339     const G4int histoIndex = fXSProfileIndex.at(xsIndex);
0340     const G4String& histoName = fAnalysisManager->GetH1Name(histoIndex);
0341     const auto& histo = rootAnalysisManager->GetH1(histoIndex);
0342
0343     tools::histo::flair::dumpG4H1ProfileInFlairFormat(output, indexInOutputFile, histoName, histo,
0344                                                       tools::histo::flair::Abscissa::KineticEnergy,
0345                                                       fBinSchemeName);
0346     ++indexInOutputFile;
0347   }
0348   output.close();
0349   G4cout << "### All profiles saved to " << fFlairOutputFileName << G4endl;
0350 }
0351
0352 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......