| 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 FinalStateHistoManager.hh
27 /// \brief Create a set of histos for final state study.
28 //
29 // Author: G.Hugo, 08 December 2022
30 //
31 // ***************************************************************************
32 //
33 // FinalStateHistoManager
34 //
35 /// Create a set of histos for final state study.
36 /// In practice, the interactions studied here are hadron nuclear inelastic interactions
37 /// (though the code is fully generic).
38 ///
39 /// Energy spectra are plotted for all encountered secondaries
40 /// (one histo per secondary).
41 /// In addition, the residual nuclei Z and A distributions are plotted.
42 ///
43 /// All histograms are G4H1.
44 /// They are created and filled solely via G4VAnalysisManager.
45 ///
46 /// The histograms can be dumped to all usual formats, including ROOT
47 /// (via G4VAnalysisManager).
48 /// An interesting added feature here, is that the plots, while being allocated
49 /// and filled via G4VAnalysisManager, are also dumped
50 /// in a Flair-compatible format (via tools::histo::flair).
51 ///
52 /// NB 1: Note that instead of a hardcoded number associated to a hardcoded set of particles,
53 /// particle PDG IDs are used to index the histos.
54 /// This allows a dynamic storage of all particles encountered in the final states.
55 ///
56 /// NB 2: tools::histo::flair code, which allows the dump of any G4H1
57 /// into Flair-compatible format, is fully application-agnostic,
58 /// and is placed in FlukaCern/utils.
59 /// It could also be added as an extension of core G4 Analysis Manager.
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........oooOO0OOooo........oooOO0OOooo....
77
78 FinalStateHistoManager::FinalStateHistoManager()
79 : fOutputFileName("all_secondaries"),
80 fRootOutputFileName(fOutputFileName + ".root"),
81 fFlairOutputFileName(fOutputFileName + ".hist"),
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::Instance()),
92 fNucleiZScoreIndex(0),
93 fNucleiAScoreIndex(1)
94 {
95 // fAnalysisManager = G4AnalysisManager::Instance();
96 // fAnalysisManager->SetDefaultFileType("root");
97 // fAnalysisManager->SetVerboseLevel(0);
98 // fOutputFileName += fAnalysisManager->GetFileType();
99 }
100
101 // ***************************************************************************
102 // Open output file + create residual nuclei histograms considered for final state study.
103 // The histograms are G4H1, created via G4VAnalysisManager.
104 // ***************************************************************************
105 void FinalStateHistoManager::Book()
106 {
107 // Open file.
108 if (!fAnalysisManager->OpenFile(fRootOutputFileName)) {
109 G4ExceptionDescription msg;
110 msg << "Booking histograms: cannot open file " << fRootOutputFileName << G4endl;
111 G4Exception("FinalStateHistoManager::Book", "Cannot open file", FatalException, msg);
112 }
113 G4cout << "### FinalStateHistoManager::Book: Successfully opended file " << fRootOutputFileName
114 << " for dumping histograms." << G4endl;
115
116 // Create the residual nuclei distributions (in Z and A).
117 const G4int nucleiZHistoIndex = fAnalysisManager->CreateH1(
118 "nucleiZ", "Residual nuclei distribution in Z", fNucleiZMax, 0.5, fNucleiZMax + 0.5);
119 auto nucleiZHistoWrapper = std::make_unique<G4H1Wrapper>(fAnalysisManager, nucleiZHistoIndex);
120 fNucleiData.insert(std::make_pair(fNucleiZScoreIndex, std::move(nucleiZHistoWrapper)));
121
122 const G4int nucleiAHistoIndex = fAnalysisManager->CreateH1(
123 "nucleiA", "Residual nuclei distribution in A", fNucleiAMax, 0.5, fNucleiAMax + 0.5);
124 auto nucleiAHistoWrapper = std::make_unique<G4H1Wrapper>(fAnalysisManager, nucleiAHistoIndex);
125 fNucleiData.insert(std::make_pair(fNucleiAScoreIndex, std::move(nucleiAHistoWrapper)));
126 }
127
128 // ***************************************************************************
129 // Keep track of the total number of events (used later on for normalization).
130 // ***************************************************************************
131 void FinalStateHistoManager::BeginOfEvent()
132 {
133 fNumEvents++;
134 }
135
136 // ***************************************************************************
137 // Fill all plots (WITHIN event, ie the interaction).
138 // ***************************************************************************
139 void FinalStateHistoManager::ScoreSecondary(const G4DynamicParticle* const secondary)
140 {
141 // SELECT SPECIFIC SECONDARIES ONLY
142 // Select by angle with beam direction
143 /* if ( (std::pow(secondary->GetMomentumDirection().x(), 2.)
144 + std::pow(secondary->GetMomentumDirection().y(), 2.))
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->GetDefinition();
154
155 // SECONDARIES ENERGY SPECTRA
156 // Dynamic creation of histos, so that all encountered particles have their own histos.
157
158 // Check whether a particle has already been encountered.
159 const auto found = fParticleData.find(secondary->GetPDGcode());
160
161 G4H1Wrapper* particleHistoWrapper = nullptr;
162 // If the particle has already been encountered, use the corresponding histos.
163 if (found != fParticleData.end()) {
164 particleHistoWrapper = found->second.get();
165 }
166 // Otherwise, create histos for that particle.
167 else {
168 const G4String& particleName = particle->GetParticleName();
169 const G4int particlePDG = secondary->GetPDGcode();
170
171 const G4String histoTitle = (particlePDG == 0 ? G4String("Particle pdg==0 spectrum")
172 : particleName + G4String(" spectrum"));
173 const G4int histoIndex =
174 fAnalysisManager->CreateH1(particleName, histoTitle, fNumBins, fMinKineticEnergy,
175 fMaxKineticEnergy, fRootEnergyUnit, fFunctionName, fBinSchemeName);
176 auto histoWrapper = std::make_unique<G4H1Wrapper>(fAnalysisManager, histoIndex);
177 particleHistoWrapper = histoWrapper.get();
178 fParticleData.insert(std::make_pair(particlePDG, std::move(histoWrapper)));
179 }
180
181 // Fill the G4H1Wrapper.
182 const G4double kineticEnergy = secondary->GetKineticEnergy();
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() / eplus;
189 fNucleiData[fNucleiZScoreIndex]->Fill(Z, 1.);
190
191 // Fill the G4H1Wrapper.
192 const G4double A = particle->GetBaryonNumber();
193 fNucleiData[fNucleiAScoreIndex]->Fill(A, 1.);
194 }
195
196 //} // select secondaries
197 }
198
199 // ***************************************************************************
200 // End of event: all event-level G4H1 are flushed into the Analysis Manager G4H1.
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 into relevant formats.
214 // ***************************************************************************
215 void FinalStateHistoManager::EndOfRun() const
216 {
217 // PRINTOUT SECONDARYS COUNTS (FULL ENERGY RANGE).
218
219 // Order the histos by particles names.
220 std::map<G4String, const G4H1Wrapper*> particlesHistos;
221
222 for (const auto& particleIt : fParticleData) {
223 const G4int particlePdg = particleIt.first;
224 const G4String particleName =
225 G4ParticleTable::GetParticleTable()->FindParticle(particlePdg)->GetParticleName();
226
227 const G4H1Wrapper* const particleHisto = particleIt.second.get();
228 particlesHistos.insert(std::make_pair(particleName, particleHisto));
229 }
230
231 // Printout secondarys counts (full energy range)
232 // Values are averaged over the number of events.
233 G4cout << "========================================================" << G4endl;
234 G4cout << "Number of events " << fNumEvents << G4endl << G4endl;
235 for (const auto& particleIt : particlesHistos) {
236 // Note that the info is directly obtained from the histogram:
237 // it is the integral over the full energy range.
238 const G4int count = particleIt.second->GetG4H1()->sum_all_bin_heights();
239
240 const G4double averageCount = static_cast<G4double>(count) / fNumEvents;
241
242 G4cout << "Average (per event) number of " << particleIt.first << " "
243 << averageCount << G4endl;
244 }
245 G4cout << "========================================================" << G4endl;
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 G4VAnalysisManager).
261 // ***************************************************************************
262 void FinalStateHistoManager::DumpAllG4H1IntoRootFile() const
263 {
264 if (!fAnalysisManager->Write()) {
265 G4ExceptionDescription message;
266 message << "Could not write ROOT file.";
267 G4Exception("FinalStateHistoManager::EndOfRun()", "I/O Error", FatalException, message);
268 }
269 G4cout << "### All histograms saved to " << fRootOutputFileName << G4endl;
270 }
271
272 // ***************************************************************************
273 // DUMP G4H1 PLOTS INTO FLAIR FILE (via tools::histo::flair).
274 // ***************************************************************************
275 void FinalStateHistoManager::DumpAllG4H1IntoFlairFile(
276 const std::map<G4String, const G4H1Wrapper*>& particlesHistos) const
277 {
278 std::ofstream output;
279 output.open(fFlairOutputFileName, std::ios_base::out);
280 G4int indexInOutputFile = 1;
281
282 // SECONDARIES ENERGY SPECTRA
283 for (const auto& particleIt : particlesHistos) {
284 const G4String& histoName = particleIt.first;
285 const auto& histo = particleIt.second->GetG4H1();
286
287 tools::histo::flair::dumpG4H1HistoInFlairFormat(
288 output, indexInOutputFile, histoName, histo, tools::histo::flair::Abscissa::KineticEnergy,
289 fBinSchemeName, fNumEvents, particleIt.second->GetSumSquaredEventTotals(),
290 particleIt.second->GetSumSquaredEventInRangeTotals());
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 == fNucleiZScoreIndex ? "nucleiZ" : "nucleiA");
298 const auto& abscissaKind =
299 (plotIt.first == fNucleiZScoreIndex ? tools::histo::flair::Abscissa::Z
300 : tools::histo::flair::Abscissa::A);
301
302 tools::histo::flair::dumpG4H1HistoInFlairFormat(
303 output, indexInOutputFile, histoName, histo, abscissaKind, fBinSchemeName, fNumEvents,
304 plotIt.second->GetSumSquaredEventTotals(), plotIt.second->GetSumSquaredEventInRangeTotals());
305 ++indexInOutputFile;
306 }
307
308 output.close();
309 G4cout << "### All histograms saved to " << fFlairOutputFileName << G4endl;
310 }
311
312 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
313