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