| Geant4 Cross Reference |
1 // 1
2 // *******************************************
3 // * License and Disclaimer
4 // *
5 // * The Geant4 software is copyright of th
6 // * the Geant4 Collaboration. It is provided
7 // * conditions of the Geant4 Software License
8 // * LICENSE and available at http://cern.ch/
9 // * include a list of copyright holders.
10 // *
11 // * Neither the authors of this software syst
12 // * institutes,nor the agencies providing fin
13 // * work make any representation or warran
14 // * regarding this software system or assum
15 // * use. Please see the license in the file
16 // * for the full disclaimer and the limitatio
17 // *
18 // * This code implementation is the result
19 // * technical work of the GEANT4 collaboratio
20 // * By using, copying, modifying or distri
21 // * any work based on the software) you ag
22 // * use in resulting scientific publicati
23 // * acceptance of all terms of the Geant4 Sof
24 // *******************************************
25 //
26 /// \file XSHistoManager.cc
27 /// \brief Create a set of profiles for XS st
28 //
29 // Adapted from hadronic/Hadr00/src/HistoMana
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 G4VAnalys
40 ///
41 /// The profiles can be dumped to all usual f
42 /// (via G4VAnalysisManager).
43 /// They are also dumped in a format compatib
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........oo
64
65 XSHistoManager::XSHistoManager()
66 : fMessenger(new XSHistoManagerMessenger(thi
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::In
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(
90 }
91
92 // *******************************************
93 // Set output files names: 2 formats supported
94 // *******************************************
95 void XSHistoManager::SetOutputFileName(const G
96 {
97 fOutputFileName = outputFileName;
98 fRootOutputFileName = outputFileName + ".roo
99 fFlairOutputFileName = outputFileName + ".hi
100 }
101
102 // *******************************************
103 // Set the particle considered for XS study.
104 // *******************************************
105 void XSHistoManager::SetParticle(const G4Strin
106 {
107 fParticle = G4ParticleTable::GetParticleTabl
108 }
109
110 // *******************************************
111 // Set the target element considered for XS st
112 // *******************************************
113 void XSHistoManager::SetElement(const G4String
114 {
115 fElement = G4NistManager::Instance()->FindOr
116 // Also needs to set material!
117 SetMaterial(elementName);
118 }
119
120 // *******************************************
121 // Set the target material considered for XS s
122 // *******************************************
123 void XSHistoManager::SetMaterial(const G4Strin
124 {
125 // Check that material is not set already.
126 if (fMaterial) {
127 G4ExceptionDescription msg;
128 msg << "Please use UI command /allXS/eleme
129 << " OR UI command /allXS/nonElementar
130 << " BUT NOT BOTH!" << G4endl;
131 G4Exception("XSHistoManager::SetMaterial",
132 FatalException, msg);
133 }
134
135 fMaterial = G4NistManager::Instance()->FindO
136 }
137
138 // *******************************************
139 // Open output file + create all profiles cons
140 // All profiles are G4H1, created via G4VAnaly
141 // *******************************************
142 void XSHistoManager::Book()
143 {
144 // Check all XSHistoManager data is set prop
145 CheckInput();
146
147 // Open file.
148 if (!fAnalysisManager->OpenFile(fRootOutputF
149 G4ExceptionDescription msg;
150 msg << "Booking profiles: cannot open file
151 G4Exception("XSHistoManager::Book", "Canno
152 }
153 G4cout << "### XSHistoManager::Book: Success
154 << " for dumping profiles." << G4endl
155
156 // Create all G4H1, and keep track of each h
157 const G4int elasticXSProfileIndex =
158 fAnalysisManager->CreateH1("ElasticXS", "E
159 fMaxKineticEner
160 fXSProfileIndex.insert(std::make_pair(fElast
161
162 const G4int inelasticXSProfileIndex =
163 fAnalysisManager->CreateH1("InelasticXS",
164 fMaxKineticEner
165 fXSProfileIndex.insert(std::make_pair(fInela
166
167 const G4int captureXSProfileIndex =
168 fAnalysisManager->CreateH1("CaptureXS", "C
169 fMaxKineticEner
170 fXSProfileIndex.insert(std::make_pair(fCaptu
171
172 const G4int fissionXSProfileIndex =
173 fAnalysisManager->CreateH1("FissionXS", "F
174 fMaxKineticEner
175 fXSProfileIndex.insert(std::make_pair(fFissi
176
177 const G4int chargeExchangeXSProfileIndex = f
178 "ChargeExchangeXS", "Charge exchange XS",
179 fRootEnergyUnit, fFunctionName, fBinScheme
180 fXSProfileIndex.insert(std::make_pair(fCharg
181
182 const G4int totalXSProfileIndex =
183 fAnalysisManager->CreateH1("TotalXS", "Tot
184 fMaxKineticEner
185 fXSProfileIndex.insert(std::make_pair(fTotal
186
187 const G4int elasticPerVolumeXSProfileIndex =
188 "ElasticPerVolumeXS", "Elastic XS per volu
189 fRootEnergyUnit, fFunctionName, fBinScheme
190 fXSProfileIndex.insert(std::make_pair(fElast
191
192 const G4int inelasticPerVolumeXSProfileIndex
193 "InelasticPerVolumeXS", "Inelastic XS per
194 fMaxKineticEnergy, fRootEnergyUnit, fFunct
195 fXSProfileIndex.insert(
196 std::make_pair(fInelasticPerVolumeXSIndex,
197 }
198
199 // *******************************************
200 // Fill all plots, then dump them into relevan
201 // *******************************************
202 void XSHistoManager::EndOfRun()
203 {
204 G4cout << "### XSHistoManager::EndOfRun: Com
205 << " in " << (fElement ? fElement->Ge
206
207 G4HadronicProcessStore* const store = G4Hadr
208
209 // Fill XS profiles.
210 const G4double logMinKineticEnergy = std::lo
211 const G4double logMaxKineticEnergy = std::lo
212 const G4double deltaLogKineticEnergy = (logM
213
214 G4double logKineticEnergy = logMinKineticEne
215
216 // Loop on all kinetic energies of interest.
217 for (G4int binIndex = 0; binIndex < fNumBins
218 logKineticEnergy += deltaLogKineticEnergy;
219 const G4double kineticEnergy = std::pow(10
220
221 G4double totalXS = 0.;
222 if (fElement) {
223 // ELASTIC (ELEMENTARY MATERIAL)
224 const G4double elasticXS =
225 store->GetElasticCrossSectionPerAtom(f
226 fAnalysisManager->FillH1(fXSProfileIndex
227 totalXS += elasticXS;
228
229 // INELASTIC (ELEMENTARY MATERIAL)
230 const G4double inelasticXS =
231 store->GetInelasticCrossSectionPerAtom
232 fAnalysisManager->FillH1(fXSProfileIndex
233 inelasticXS / b
234 totalXS += inelasticXS;
235
236 if (fParticle == G4Neutron::Definition()
237 // NEUTRON CAPTURE (ELEMENTARY MATERIA
238 const G4double captureXS =
239 store->GetCaptureCrossSectionPerAtom
240 fAnalysisManager->FillH1(fXSProfileInd
241 totalXS += captureXS;
242
243 // FISSION (ELEMENTARY MATERIAL)
244 const G4double fissionXS =
245 store->GetFissionCrossSectionPerAtom
246 totalXS += fissionXS;
247 fAnalysisManager->FillH1(fXSProfileInd
248 }
249
250 // CHARGE EXCHANGE (ELEMENTARY MATERIAL)
251 const G4double chargeExchangeXS =
252 store->GetChargeExchangeCrossSectionPe
253 fAnalysisManager->FillH1(fXSProfileIndex
254 chargeExchangeX
255 totalXS += chargeExchangeXS;
256
257 // TOTAL (ELEMENTARY MATERIAL)
258 fAnalysisManager->FillH1(fXSProfileIndex
259 }
260
261 if (fMaterial) {
262 const G4double materialSurfacicDensity =
263 (fMaterial ? fMaterial->GetDensity() /
264
265 // ELASTIC
266 const G4double elasticPerVolumeXS =
267 store->GetElasticCrossSectionPerVolume
268 fAnalysisManager->FillH1(fXSProfileIndex
269 elasticPerVolum
270
271 // INELASTIC
272 const G4double inelasticPerVolumeXS =
273 store->GetInelasticCrossSectionPerVolu
274 fAnalysisManager->FillH1(fXSProfileIndex
275 inelasticPerVol
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:
299 G4Exception("XSHistoManager::CheckInput()"
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
307 << " or UI command /allXS/nonElementar
308 << G4endl;
309 G4Exception("XSHistoManager::CheckInput()"
310 FatalException, msg);
311 }
312 }
313
314 // *******************************************
315 // DUMP G4H1 PLOTS INTO ROOT FILE (via G4VAnal
316 // *******************************************
317 void XSHistoManager::DumpAllG4H1IntoRootFile()
318 {
319 if (!fAnalysisManager->Write()) {
320 G4ExceptionDescription message;
321 message << "Could not write ROOT file.";
322 G4Exception("XSHistoManager::EndOfRun()",
323 }
324 G4cout << "### All profiles saved to " << fR
325 }
326
327 // *******************************************
328 // DUMP G4H1 PLOTS INTO FLAIR FILE (via tools:
329 // *******************************************
330 void XSHistoManager::DumpAllG4H1IntoFlairFile(
331 {
332 std::ofstream output;
333 output.open(fFlairOutputFileName, std::ios_b
334 auto const rootAnalysisManager = dynamic_cas
335
336 G4int indexInOutputFile = 1;
337 for (G4int xsIndex = fElasticXSIndex; xsInde
338 const G4int histoIndex = fXSProfileIndex.a
339 const G4String& histoName = fAnalysisManag
340 const auto& histo = rootAnalysisManager->G
341
342 tools::histo::flair::dumpG4H1ProfileInFlai
343
344
345 ++indexInOutputFile;
346 }
347 output.close();
348 G4cout << "### All profiles saved to " << fF
349 }
350
351 //....oooOO0OOooo........oooOO0OOooo........oo
352