| 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 FinalStateHistoManager.hh
27 /// \brief Create a set of histos for final s
28 //
29 // Author: G.Hugo, 08 December 2022
30 //
31 // *******************************************
32 //
33 // FinalStateHistoManager
34 //
35 /// Create a set of histos for final state st
36 /// In practice, the interactions studied her
37 /// (though the code is fully generic).
38 ///
39 /// Energy spectra are plotted for all encoun
40 /// (one histo per secondary).
41 /// In addition, the residual nuclei Z and A
42 ///
43 /// All histograms are G4H1.
44 /// They are created and filled solely via G4
45 ///
46 /// The histograms can be dumped to all usual
47 /// (via G4VAnalysisManager).
48 /// An interesting added feature here, is tha
49 /// and filled via G4VAnalysisManager, are al
50 /// in a Flair-compatible format (via tools::
51 ///
52 /// NB 1: Note that instead of a hardcoded nu
53 /// particle PDG IDs are used to index the hi
54 /// This allows a dynamic storage of all part
55 ///
56 /// NB 2: tools::histo::flair code, which all
57 /// into Flair-compatible format, is fully ap
58 /// and is placed in FlukaCern/utils.
59 /// It could also be added as an extension of
60 //
61 // *******************************************
62
63 #include "FinalStateHistoManager.hh"
64
65 #include "G4RootAnalysisManager.hh"
66 // #include "G4AnalysisManager.hh"
67
68 #include "tools_histo_flair.hh"
69
70 #include "G4DynamicParticle.hh"
71 #include "G4Exception.hh"
72 #include "G4ParticleTable.hh"
73 #include "G4ios.hh"
74 #include "g4hntools_defs.hh"
75
76 //....oooOO0OOooo........oooOO0OOooo........oo
77
78 FinalStateHistoManager::FinalStateHistoManager
79 : fOutputFileName("all_secondaries"),
80 fRootOutputFileName(fOutputFileName + ".ro
81 fFlairOutputFileName(fOutputFileName + ".h
82 fNumBins(90),
83 fMinKineticEnergy(10. * keV),
84 fMaxKineticEnergy(10. * TeV),
85 fFunctionName("none"),
86 fBinSchemeName("log"),
87 fRootEnergyUnit("MeV"),
88 fNucleiZMax(25),
89 fNucleiAMax(50),
90 fNumEvents(0),
91 fAnalysisManager(G4RootAnalysisManager::In
92 fNucleiZScoreIndex(0),
93 fNucleiAScoreIndex(1)
94 {
95 // fAnalysisManager = G4AnalysisManager::Ins
96 // fAnalysisManager->SetDefaultFileType("roo
97 // fAnalysisManager->SetVerboseLevel(0);
98 // fOutputFileName += fAnalysisManager->GetF
99 }
100
101 // *******************************************
102 // Open output file + create residual nuclei h
103 // The histograms are G4H1, created via G4VAna
104 // *******************************************
105 void FinalStateHistoManager::Book()
106 {
107 // Open file.
108 if (!fAnalysisManager->OpenFile(fRootOutputF
109 G4ExceptionDescription msg;
110 msg << "Booking histograms: cannot open fi
111 G4Exception("FinalStateHistoManager::Book"
112 }
113 G4cout << "### FinalStateHistoManager::Book:
114 << " for dumping histograms." << G4en
115
116 // Create the residual nuclei distributions
117 const G4int nucleiZHistoIndex = fAnalysisMan
118 "nucleiZ", "Residual nuclei distribution i
119 auto nucleiZHistoWrapper = std::make_unique<
120 fNucleiData.insert(std::make_pair(fNucleiZSc
121
122 const G4int nucleiAHistoIndex = fAnalysisMan
123 "nucleiA", "Residual nuclei distribution i
124 auto nucleiAHistoWrapper = std::make_unique<
125 fNucleiData.insert(std::make_pair(fNucleiASc
126 }
127
128 // *******************************************
129 // Keep track of the total number of events (u
130 // *******************************************
131 void FinalStateHistoManager::BeginOfEvent()
132 {
133 fNumEvents++;
134 }
135
136 // *******************************************
137 // Fill all plots (WITHIN event, ie the intera
138 // *******************************************
139 void FinalStateHistoManager::ScoreSecondary(co
140 {
141 // SELECT SPECIFIC SECONDARIES ONLY
142 // Select by angle with beam direction
143 /* if ( (std::pow(secondary->GetMomentumDire
144 + std::pow(secondary->GetMomentumDirectio
145 <= 0.0001 ) {*/
146
147 // Select by production tag
148 /* if (secondary->GetProductionTag() == 6) {
149
150 // Primary track
151 /* if(track->GetParentID() == 0) {*/
152
153 const auto& particle = secondary->GetDefinit
154
155 // SECONDARIES ENERGY SPECTRA
156 // Dynamic creation of histos, so that all e
157
158 // Check whether a particle has already been
159 const auto found = fParticleData.find(second
160
161 G4H1Wrapper* particleHistoWrapper = nullptr;
162 // If the particle has already been encounte
163 if (found != fParticleData.end()) {
164 particleHistoWrapper = found->second.get()
165 }
166 // Otherwise, create histos for that particl
167 else {
168 const G4String& particleName = particle->G
169 const G4int particlePDG = secondary->GetPD
170
171 const G4String histoTitle = (particlePDG =
172
173 const G4int histoIndex =
174 fAnalysisManager->CreateH1(particleName,
175 fMaxKineticEn
176 auto histoWrapper = std::make_unique<G4H1W
177 particleHistoWrapper = histoWrapper.get();
178 fParticleData.insert(std::make_pair(partic
179 }
180
181 // Fill the G4H1Wrapper.
182 const G4double kineticEnergy = secondary->Ge
183 particleHistoWrapper->Fill(kineticEnergy, 1.
184
185 // NUCLEI DISTRIBUTIONS IN Z AND A
186 if (particle->GetParticleType() == "nucleus"
187 // Fill the G4H1Wrapper.
188 const G4double Z = particle->GetPDGCharge(
189 fNucleiData[fNucleiZScoreIndex]->Fill(Z, 1
190
191 // Fill the G4H1Wrapper.
192 const G4double A = particle->GetBaryonNumb
193 fNucleiData[fNucleiAScoreIndex]->Fill(A, 1
194 }
195
196 //} // select secondaries
197 }
198
199 // *******************************************
200 // End of event: all event-level G4H1 are flus
201 // *******************************************
202 void FinalStateHistoManager::EndOfEvent()
203 {
204 for (const auto& particleIt : fParticleData)
205 particleIt.second->EndOfEvent();
206 }
207 for (const auto& nucleiScoreIt : fNucleiData
208 nucleiScoreIt.second->EndOfEvent();
209 }
210 }
211
212 // *******************************************
213 // Printout secondary counts + dump all plots
214 // *******************************************
215 void FinalStateHistoManager::EndOfRun() const
216 {
217 // PRINTOUT SECONDARYS COUNTS (FULL ENERGY R
218
219 // Order the histos by particles names.
220 std::map<G4String, const G4H1Wrapper*> parti
221
222 for (const auto& particleIt : fParticleData)
223 const G4int particlePdg = particleIt.first
224 const G4String particleName =
225 G4ParticleTable::GetParticleTable()->Fin
226
227 const G4H1Wrapper* const particleHisto = p
228 particlesHistos.insert(std::make_pair(part
229 }
230
231 // Printout secondarys counts (full energy r
232 // Values are averaged over the number of ev
233 G4cout << "=================================
234 G4cout << "Number of events
235 for (const auto& particleIt : particlesHisto
236 // Note that the info is directly obtained
237 // it is the integral over the full energy
238 const G4int count = particleIt.second->Get
239
240 const G4double averageCount = static_cast<
241
242 G4cout << "Average (per event) number of "
243 << averageCount << G4endl;
244 }
245 G4cout << "=================================
246 G4cout << G4endl;
247
248 // DUMP G4H1 PLOTS INTO ROOT FILE
249 DumpAllG4H1IntoRootFile();
250
251 // DUMP G4H1 PLOTS INTO FLAIR FILE
252 DumpAllG4H1IntoFlairFile(particlesHistos);
253
254 // Close and clear fAnalysisManager.
255 fAnalysisManager->CloseFile();
256 fAnalysisManager->Clear();
257 }
258
259 // *******************************************
260 // DUMP G4H1 PLOTS INTO ROOT FILE (via G4VAnal
261 // *******************************************
262 void FinalStateHistoManager::DumpAllG4H1IntoRo
263 {
264 if (!fAnalysisManager->Write()) {
265 G4ExceptionDescription message;
266 message << "Could not write ROOT file.";
267 G4Exception("FinalStateHistoManager::EndOf
268 }
269 G4cout << "### All histograms saved to " <<
270 }
271
272 // *******************************************
273 // DUMP G4H1 PLOTS INTO FLAIR FILE (via tools:
274 // *******************************************
275 void FinalStateHistoManager::DumpAllG4H1IntoFl
276 const std::map<G4String, const G4H1Wrapper*>
277 {
278 std::ofstream output;
279 output.open(fFlairOutputFileName, std::ios_b
280 G4int indexInOutputFile = 1;
281
282 // SECONDARIES ENERGY SPECTRA
283 for (const auto& particleIt : particlesHisto
284 const G4String& histoName = particleIt.fir
285 const auto& histo = particleIt.second->Get
286
287 tools::histo::flair::dumpG4H1HistoInFlairF
288 output, indexInOutputFile, histoName, hi
289 fBinSchemeName, fNumEvents, particleIt.s
290 particleIt.second->GetSumSquaredEventInR
291 ++indexInOutputFile;
292 }
293
294 // RESIDUAL NUCLEI DISTRIBUTIONS
295 for (const auto& plotIt : fNucleiData) {
296 const auto& histo = plotIt.second->GetG4H1
297 const G4String& histoName = (plotIt.first
298 const auto& abscissaKind =
299 (plotIt.first == fNucleiZScoreIndex ? to
300 : to
301
302 tools::histo::flair::dumpG4H1HistoInFlairF
303 output, indexInOutputFile, histoName, hi
304 plotIt.second->GetSumSquaredEventTotals(
305 ++indexInOutputFile;
306 }
307
308 output.close();
309 G4cout << "### All histograms saved to " <<
310 }
311
312 //....oooOO0OOooo........oooOO0OOooo........oo
313