| Geant4 Cross Reference |
1 //
2 // ********************************************************************
3 // * License and Disclaimer *
4 // * *
5 // * The Geant4 software is copyright of the Copyright Holders of *
6 // * the Geant4 Collaboration. It is provided under the terms and *
7 // * conditions of the Geant4 Software License, included in the file *
8 // * LICENSE and available at http://cern.ch/geant4/license . These *
9 // * include a list of copyright holders. *
10 // * *
11 // * Neither the authors of this software system, nor their employing *
12 // * institutes,nor the agencies providing financial support for this *
13 // * work make any representation or warranty, express or implied, *
14 // * regarding this software system or assume any liability for its *
15 // * use. Please see the license in the file LICENSE and URL above *
16 // * for the full disclaimer and the limitation of liability. *
17 // * *
18 // * This code implementation is the result of the scientific and *
19 // * technical work of the GEANT4 collaboration. *
20 // * By using, copying, modifying or distributing the software (or *
21 // * any work based on the software) you agree to acknowledge its *
22 // * use in resulting scientific publications, and indicate your *
23 // * acceptance of all terms of the Geant4 Software license. *
24 // ********************************************************************
25 //
26 /// \file XSHistoManager.cc
27 /// \brief Create a set of profiles for XS study.
28 //
29 // Adapted from hadronic/Hadr00/src/HistoManager.cc
30 // Author: G.Hugo, 06 January 2023
31 //
32 // ***************************************************************************
33 //
34 // XSHistoManager
35 //
36 /// Create a set of profiles for XS study.
37 ///
38 /// All profiles are G4H1.
39 /// They are created and filled via G4VAnalysisManager.
40 ///
41 /// The profiles can be dumped to all usual formats, including ROOT
42 /// (via G4VAnalysisManager).
43 /// They are also dumped in a format compatible with Flair
44 /// (via tools::histo::flair).
45 //
46 // ***************************************************************************
47
48 #include "XSHistoManager.hh"
49
50 #include "G4Element.hh"
51 #include "G4HadronicProcessStore.hh"
52 #include "G4Material.hh"
53 #include "G4NistManager.hh"
54 #include "G4ParticleDefinition.hh"
55 #include "G4ParticleTable.hh"
56 #include "G4ios.hh"
57
58 // #include "G4AnalysisManager.hh"
59 #include "tools_histo_flair.hh"
60
61 #include "G4RootAnalysisManager.hh"
62
63 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
64
65 XSHistoManager::XSHistoManager()
66 : fMessenger(new XSHistoManagerMessenger(this)),
67 fOutputFileName("all_XS"),
68 fRootOutputFileName("all_XS.root"),
69 fFlairOutputFileName("all_XS.hist"),
70 fParticle(nullptr),
71 fElement(nullptr),
72 fMaterial(nullptr),
73 fNumBins(10000),
74 fMinKineticEnergy(1. * keV),
75 fMaxKineticEnergy(10. * TeV),
76 fFunctionName("none"),
77 fBinSchemeName("log"),
78 fRootEnergyUnit("MeV"),
79 fAnalysisManager(G4RootAnalysisManager::Instance()),
80 fElasticXSIndex(0),
81 fInelasticXSIndex(1),
82 fCaptureXSIndex(2),
83 fFissionXSIndex(3),
84 fChargeExchangeXSIndex(4),
85 fTotalXSIndex(5),
86 fElasticPerVolumeXSIndex(6),
87 fInelasticPerVolumeXSIndex(7)
88 {
89 // G4NistManager::Instance()->ListMaterials("all");
90 }
91
92 // ***************************************************************************
93 // Set output files names: 2 formats supported, ROOT and Flair.
94 // ***************************************************************************
95 void XSHistoManager::SetOutputFileName(const G4String& outputFileName)
96 {
97 fOutputFileName = outputFileName;
98 fRootOutputFileName = outputFileName + ".root";
99 fFlairOutputFileName = outputFileName + ".hist";
100 }
101
102 // ***************************************************************************
103 // Set the particle considered for XS study.
104 // ***************************************************************************
105 void XSHistoManager::SetParticle(const G4String& particleName)
106 {
107 fParticle = G4ParticleTable::GetParticleTable()->FindParticle(particleName);
108 }
109
110 // ***************************************************************************
111 // Set the target element considered for XS study.
112 // ***************************************************************************
113 void XSHistoManager::SetElement(const G4String& elementName)
114 {
115 fElement = G4NistManager::Instance()->FindOrBuildElement(elementName);
116 // Also needs to set material!
117 SetMaterial(elementName);
118 }
119
120 // ***************************************************************************
121 // Set the target material considered for XS study.
122 // ***************************************************************************
123 void XSHistoManager::SetMaterial(const G4String& materialName)
124 {
125 // Check that material is not set already.
126 if (fMaterial) {
127 G4ExceptionDescription msg;
128 msg << "Please use UI command /allXS/elementName"
129 << " OR UI command /allXS/nonElementaryMaterialName,"
130 << " BUT NOT BOTH!" << G4endl;
131 G4Exception("XSHistoManager::SetMaterial", "A target material is already defined.",
132 FatalException, msg);
133 }
134
135 fMaterial = G4NistManager::Instance()->FindOrBuildMaterial("G4_" + materialName);
136 }
137
138 // ***************************************************************************
139 // Open output file + create all profiles considered for XS study.
140 // All profiles are G4H1, created via G4VAnalysisManager.
141 // ***************************************************************************
142 void XSHistoManager::Book()
143 {
144 // Check all XSHistoManager data is set properly.
145 CheckInput();
146
147 // Open file.
148 if (!fAnalysisManager->OpenFile(fRootOutputFileName)) {
149 G4ExceptionDescription msg;
150 msg << "Booking profiles: cannot open file " << fRootOutputFileName << G4endl;
151 G4Exception("XSHistoManager::Book", "Cannot open file", FatalException, msg);
152 }
153 G4cout << "### XSHistoManager::Book: Successfully opened file " << fRootOutputFileName
154 << " for dumping profiles." << G4endl;
155
156 // Create all G4H1, and keep track of each histo index in fXSProfileIndex.
157 const G4int elasticXSProfileIndex =
158 fAnalysisManager->CreateH1("ElasticXS", "Elastic XS", fNumBins, fMinKineticEnergy,
159 fMaxKineticEnergy, fRootEnergyUnit, fFunctionName, fBinSchemeName);
160 fXSProfileIndex.insert(std::make_pair(fElasticXSIndex, elasticXSProfileIndex));
161
162 const G4int inelasticXSProfileIndex =
163 fAnalysisManager->CreateH1("InelasticXS", "Inelastic XS", fNumBins, fMinKineticEnergy,
164 fMaxKineticEnergy, fRootEnergyUnit, fFunctionName, fBinSchemeName);
165 fXSProfileIndex.insert(std::make_pair(fInelasticXSIndex, inelasticXSProfileIndex));
166
167 const G4int captureXSProfileIndex =
168 fAnalysisManager->CreateH1("CaptureXS", "Capture XS", fNumBins, fMinKineticEnergy,
169 fMaxKineticEnergy, fRootEnergyUnit, fFunctionName, fBinSchemeName);
170 fXSProfileIndex.insert(std::make_pair(fCaptureXSIndex, captureXSProfileIndex));
171
172 const G4int fissionXSProfileIndex =
173 fAnalysisManager->CreateH1("FissionXS", "Fission XS", fNumBins, fMinKineticEnergy,
174 fMaxKineticEnergy, fRootEnergyUnit, fFunctionName, fBinSchemeName);
175 fXSProfileIndex.insert(std::make_pair(fFissionXSIndex, fissionXSProfileIndex));
176
177 const G4int chargeExchangeXSProfileIndex = fAnalysisManager->CreateH1(
178 "ChargeExchangeXS", "Charge exchange XS", fNumBins, fMinKineticEnergy, fMaxKineticEnergy,
179 fRootEnergyUnit, fFunctionName, fBinSchemeName);
180 fXSProfileIndex.insert(std::make_pair(fChargeExchangeXSIndex, chargeExchangeXSProfileIndex));
181
182 const G4int totalXSProfileIndex =
183 fAnalysisManager->CreateH1("TotalXS", "Total XS", fNumBins, fMinKineticEnergy,
184 fMaxKineticEnergy, fRootEnergyUnit, fFunctionName, fBinSchemeName);
185 fXSProfileIndex.insert(std::make_pair(fTotalXSIndex, totalXSProfileIndex));
186
187 const G4int elasticPerVolumeXSProfileIndex = fAnalysisManager->CreateH1(
188 "ElasticPerVolumeXS", "Elastic XS per volume", fNumBins, fMinKineticEnergy, fMaxKineticEnergy,
189 fRootEnergyUnit, fFunctionName, fBinSchemeName);
190 fXSProfileIndex.insert(std::make_pair(fElasticPerVolumeXSIndex, elasticPerVolumeXSProfileIndex));
191
192 const G4int inelasticPerVolumeXSProfileIndex = fAnalysisManager->CreateH1(
193 "InelasticPerVolumeXS", "Inelastic XS per volume", fNumBins, fMinKineticEnergy,
194 fMaxKineticEnergy, fRootEnergyUnit, fFunctionName, fBinSchemeName);
195 fXSProfileIndex.insert(
196 std::make_pair(fInelasticPerVolumeXSIndex, inelasticPerVolumeXSProfileIndex));
197 }
198
199 // ***************************************************************************
200 // Fill all plots, then dump them into relevant formats.
201 // ***************************************************************************
202 void XSHistoManager::EndOfRun()
203 {
204 G4cout << "### XSHistoManager::EndOfRun: Compute & fill XS for " << fParticle->GetParticleName()
205 << " in " << (fElement ? fElement->GetName() : fMaterial->GetName()) << G4endl;
206
207 G4HadronicProcessStore* const store = G4HadronicProcessStore::Instance();
208
209 // Fill XS profiles.
210 const G4double logMinKineticEnergy = std::log10(fMinKineticEnergy);
211 const G4double logMaxKineticEnergy = std::log10(fMaxKineticEnergy);
212 const G4double deltaLogKineticEnergy = (logMaxKineticEnergy - logMinKineticEnergy) / fNumBins;
213
214 G4double logKineticEnergy = logMinKineticEnergy - deltaLogKineticEnergy / 2.;
215
216 // Loop on all kinetic energies of interest.
217 for (G4int binIndex = 0; binIndex < fNumBins; ++binIndex) {
218 logKineticEnergy += deltaLogKineticEnergy;
219 const G4double kineticEnergy = std::pow(10., logKineticEnergy) * MeV;
220
221 G4double totalXS = 0.;
222 if (fElement) {
223 // ELASTIC (ELEMENTARY MATERIAL)
224 const G4double elasticXS =
225 store->GetElasticCrossSectionPerAtom(fParticle, kineticEnergy, fElement, fMaterial);
226 fAnalysisManager->FillH1(fXSProfileIndex[fElasticXSIndex], kineticEnergy, elasticXS / barn);
227 totalXS += elasticXS;
228
229 // INELASTIC (ELEMENTARY MATERIAL)
230 const G4double inelasticXS =
231 store->GetInelasticCrossSectionPerAtom(fParticle, kineticEnergy, fElement, fMaterial);
232 fAnalysisManager->FillH1(fXSProfileIndex[fInelasticXSIndex], kineticEnergy,
233 inelasticXS / barn);
234 totalXS += inelasticXS;
235
236 if (fParticle == G4Neutron::Definition()) {
237 // NEUTRON CAPTURE (ELEMENTARY MATERIAL)
238 const G4double captureXS =
239 store->GetCaptureCrossSectionPerAtom(fParticle, kineticEnergy, fElement, fMaterial);
240 fAnalysisManager->FillH1(fXSProfileIndex[fCaptureXSIndex], kineticEnergy, captureXS / barn);
241 totalXS += captureXS;
242
243 // FISSION (ELEMENTARY MATERIAL)
244 const G4double fissionXS =
245 store->GetFissionCrossSectionPerAtom(fParticle, kineticEnergy, fElement, fMaterial);
246 totalXS += fissionXS;
247 fAnalysisManager->FillH1(fXSProfileIndex[fFissionXSIndex], kineticEnergy, fissionXS / barn);
248 }
249
250 // CHARGE EXCHANGE (ELEMENTARY MATERIAL)
251 const G4double chargeExchangeXS =
252 store->GetChargeExchangeCrossSectionPerAtom(fParticle, kineticEnergy, fElement, fMaterial);
253 fAnalysisManager->FillH1(fXSProfileIndex[fChargeExchangeXSIndex], kineticEnergy,
254 chargeExchangeXS / barn);
255 totalXS += chargeExchangeXS;
256
257 // TOTAL (ELEMENTARY MATERIAL)
258 fAnalysisManager->FillH1(fXSProfileIndex[fTotalXSIndex], kineticEnergy, totalXS / barn);
259 }
260
261 if (fMaterial) {
262 const G4double materialSurfacicDensity =
263 (fMaterial ? fMaterial->GetDensity() / (g / cm2) : 1.);
264
265 // ELASTIC
266 const G4double elasticPerVolumeXS =
267 store->GetElasticCrossSectionPerVolume(fParticle, kineticEnergy, fMaterial);
268 fAnalysisManager->FillH1(fXSProfileIndex[fElasticPerVolumeXSIndex], kineticEnergy,
269 elasticPerVolumeXS / materialSurfacicDensity);
270
271 // INELASTIC
272 const G4double inelasticPerVolumeXS =
273 store->GetInelasticCrossSectionPerVolume(fParticle, kineticEnergy, fMaterial);
274 fAnalysisManager->FillH1(fXSProfileIndex[fInelasticPerVolumeXSIndex], kineticEnergy,
275 inelasticPerVolumeXS / materialSurfacicDensity);
276 }
277 }
278
279 // DUMP G4H1 PLOTS INTO ROOT FILE
280 DumpAllG4H1IntoRootFile();
281
282 // DUMP G4H1 PLOTS INTO FLAIR FILE
283 DumpAllG4H1IntoFlairFile();
284
285 // Close and clear fAnalysisManager.
286 fAnalysisManager->CloseFile();
287 fAnalysisManager->Clear();
288 }
289
290 // ***************************************************************************
291 // Checks that particle and material are set
292 // (all others have relevant default values).
293 // ***************************************************************************
294 void XSHistoManager::CheckInput()
295 {
296 if (!fParticle) {
297 G4ExceptionDescription msg;
298 msg << "Please add a particle to study XS: UI command /allXS/particleName" << G4endl;
299 G4Exception("XSHistoManager::CheckInput()", "Print XS: no input particle defined.",
300 FatalException, msg);
301 }
302
303 if (!fMaterial) {
304 G4ExceptionDescription msg;
305 msg << "Please add a material to study XS:"
306 << " UI command /allXS/elementName for an elementary material,"
307 << " or UI command /allXS/nonElementaryMaterialName for a compound/mixture material."
308 << G4endl;
309 G4Exception("XSHistoManager::CheckInput()", "Print XS: no target material defined.",
310 FatalException, msg);
311 }
312 }
313
314 // ***************************************************************************
315 // DUMP G4H1 PLOTS INTO ROOT FILE (via G4VAnalysisManager).
316 // ***************************************************************************
317 void XSHistoManager::DumpAllG4H1IntoRootFile() const
318 {
319 if (!fAnalysisManager->Write()) {
320 G4ExceptionDescription message;
321 message << "Could not write ROOT file.";
322 G4Exception("XSHistoManager::EndOfRun()", "I/O Error", FatalException, message);
323 }
324 G4cout << "### All profiles saved to " << fRootOutputFileName << G4endl;
325 }
326
327 // ***************************************************************************
328 // DUMP G4H1 PLOTS INTO FLAIR FILE (via tools::histo::flair).
329 // ***************************************************************************
330 void XSHistoManager::DumpAllG4H1IntoFlairFile() const
331 {
332 std::ofstream output;
333 output.open(fFlairOutputFileName, std::ios_base::out);
334 auto const rootAnalysisManager = dynamic_cast<G4RootAnalysisManager*>(fAnalysisManager);
335
336 G4int indexInOutputFile = 1;
337 for (G4int xsIndex = fElasticXSIndex; xsIndex <= fInelasticPerVolumeXSIndex; ++xsIndex) {
338 const G4int histoIndex = fXSProfileIndex.at(xsIndex);
339 const G4String& histoName = fAnalysisManager->GetH1Name(histoIndex);
340 const auto& histo = rootAnalysisManager->GetH1(histoIndex);
341
342 tools::histo::flair::dumpG4H1ProfileInFlairFormat(output, indexInOutputFile, histoName, histo,
343 tools::histo::flair::Abscissa::KineticEnergy,
344 fBinSchemeName);
345 ++indexInOutputFile;
346 }
347 output.close();
348 G4cout << "### All profiles saved to " << fFlairOutputFileName << G4endl;
349 }
350
351 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
352