| 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 tools_histo_flair.cc
27 /// \brief Tools to dump any G4H1 into Flair-compatible format.
28 //
29 // Author: G.Hugo, 08 December 2022
30 //
31 // ***************************************************************************
32 //
33 // tools_histo_flair
34 //
35 /// Tools to dump any G4H1 into Flair-compatible format.
36 ///
37 /// These tools are fully application-agnostic,
38 /// and could also be placed outside of the G4 examples,
39 /// as a core G4 Analysis Manager extension.
40 //
41 // ***************************************************************************
42
43 #include "tools_histo_flair.hh"
44
45 #include "G4SystemOfUnits.hh"
46 #include "G4UnitsTable.hh"
47 #include "G4ios.hh"
48
49 #include <cstring>
50 #include <fstream>
51 #include <iomanip>
52 #include <memory>
53 #include <stdarg.h>
54 #include <stdlib.h>
55
56 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
57
58 namespace tools
59 {
60 namespace histo
61 {
62 namespace flair
63 {
64
65 // ***************************************************************************
66 // Returns abscissa name.
67 // ***************************************************************************
68 G4String getAbscissaName(const Abscissa abscissaKind)
69 {
70 G4String abscissaName = "";
71 if (abscissaKind == Abscissa::KineticEnergy) {
72 abscissaName = "kE";
73 }
74 else if (abscissaKind == Abscissa::Z) {
75 abscissaName = "Z";
76 }
77 else if (abscissaKind == Abscissa::A) {
78 abscissaName = "A";
79 }
80 return abscissaName;
81 }
82
83 // ***************************************************************************
84 // Returns abscissa unit.
85 // ***************************************************************************
86 G4String getAbscissaUnit(const Abscissa abscissaKind)
87 {
88 G4String abscissaUnit = "";
89 if (abscissaKind == Abscissa::KineticEnergy) {
90 abscissaUnit = "GeV";
91 }
92 return abscissaUnit;
93 }
94
95 // ***************************************************************************
96 // Returns abscissa unit value.
97 // ***************************************************************************
98 G4double getAbscissaUnitValue(const Abscissa abscissaKind)
99 {
100 G4double abscissaUnitValue = 1.;
101 if (abscissaKind == Abscissa::KineticEnergy) {
102 abscissaUnitValue = G4UnitDefinition::GetValueOf(getAbscissaUnit(Abscissa::KineticEnergy));
103 }
104 return abscissaUnitValue;
105 }
106
107 // ***************************************************************************
108 // String format like printf() for strings.
109 // This sformat function is extracted as-is from geoviewer,
110 // and was written by V. Vlachoudis in 2010.
111 // ***************************************************************************
112 G4String sformat(const char* fmt, ...)
113 {
114 std::unique_ptr<char[]> formatted;
115 G4int n = (G4int)strlen(fmt) * 2; // Reserve two times as much as the length of the fmt_str.
116 while (true) {
117 // Wrap the plain char array into the unique_ptr.
118 formatted.reset(new char[n]);
119 strcpy(&formatted[0], fmt);
120 va_list ap;
121 va_start(ap, fmt);
122 G4int final_n = vsnprintf(&formatted[0], n, fmt, ap);
123 va_end(ap);
124 if (final_n < 0 || final_n >= n) {
125 n += std::abs(final_n - n + 1);
126 }
127 else {
128 break;
129 }
130 }
131 return G4String(formatted.get());
132 }
133
134 // ***************************************************************************
135 // Returns MC error.
136 // ***************************************************************************
137 G4double computeError(const G4double mean, const G4double sumSquares, const G4int numEvents)
138 {
139 const G4double var = std::max(0., sumSquares / numEvents - mean * mean);
140 const G4double err =
141 ((mean != 0. && numEvents > 1) ? std::sqrt(var / (numEvents - 1)) / mean : 1.);
142
143 return err * 100.;
144 }
145
146 // ***************************************************************************
147 // Dump G4H1 (histo mode) to a Flair-compatible format.
148 // ***************************************************************************
149 void dumpG4H1HistoInFlairFormat(std::ofstream& output, const G4int indexInOutputFile,
150 const G4String& histoName, G4H1* const histo,
151 const Abscissa abscissaKind, const G4String& binSchemeName,
152 const G4int numEvents, const G4double sumSquaredEventTotals,
153 const G4double sumSquaredEventInRangeTotals)
154 {
155 const G4bool isProfile = false;
156 dumpG4H1InFlairFormat(output, indexInOutputFile, histoName, histo, abscissaKind, binSchemeName,
157 numEvents, isProfile, sumSquaredEventTotals, sumSquaredEventInRangeTotals);
158 }
159
160 // ***************************************************************************
161 // Dump G4H1 (profile mode, ie no stats) to a Flair-compatible format.
162 // ***************************************************************************
163 void dumpG4H1ProfileInFlairFormat(std::ofstream& output, const G4int indexInOutputFile,
164 const G4String& histoName, G4H1* const histo,
165 const Abscissa abscissaKind, const G4String& binSchemeName)
166 {
167 const G4bool isProfile = true;
168 const G4int numEvents = 1;
169 dumpG4H1InFlairFormat(output, indexInOutputFile, histoName, histo, abscissaKind, binSchemeName,
170 numEvents, isProfile);
171 }
172
173 // ***************************************************************************
174 // Dump G4H1 (profile + histo modes) to a Flair-compatible format.
175 // ***************************************************************************
176 void dumpG4H1InFlairFormat(std::ofstream& output, const G4int indexInOutputFile,
177 const G4String& histoName, G4H1* const histo,
178 const Abscissa abscissaKind, const G4String& binSchemeName,
179 const G4int numEvents, const G4bool isProfile,
180 const G4double sumSquaredEventTotals,
181 const G4double sumSquaredEventInRangeTotals)
182 {
183 const G4int numBins = histo->axis().bins();
184 const G4double minAbscissa = histo->axis().lower_edge();
185 const G4double maxAbscissa = histo->axis().upper_edge();
186
187 const G4String abscissaName = getAbscissaName(abscissaKind);
188 const G4String abscissaUnit = getAbscissaUnit(abscissaKind);
189 const G4double abscissaUnitValue = getAbscissaUnitValue(abscissaKind);
190
191 // START HEADER.
192 if (indexInOutputFile != 1) {
193 output << G4endl << G4endl;
194 }
195 output << std::left << std::setw(13) << "# Detector: " << indexInOutputFile << " " << histoName
196 << G4endl;
197 output << std::left << std::setw(13) << "# Dim: " << histo->get_dimension() << G4endl;
198 output << std::left << std::setw(13) << "# Entries: " << numEvents << G4endl;
199
200 if (!isProfile) {
201 // Total integral
202 const G4double total = static_cast<G4double>(histo->sum_all_bin_heights()) / numEvents;
203 const G4double totalError = computeError(total, sumSquaredEventTotals, numEvents);
204 output << std::left << std::setw(13)
205 << "# Total: " << sformat("%-15g [1/pr] +/- %6.2f [%] \n", total, totalError);
206
207 // In range integral
208 const G4double inRange = static_cast<G4double>(histo->sum_bin_heights()) / numEvents;
209 const G4double inRangeError = computeError(inRange, sumSquaredEventInRangeTotals, numEvents);
210 output << std::left << std::setw(13)
211 << "# InRange: " << sformat("%-15g [1/pr] +/- %6.2f [%] \n", inRange, inRangeError);
212
213 // Underflow
214 const G4double underflow = static_cast<G4double>(histo->bins_sum_w().at(0)) / numEvents;
215 const G4double underflowError = computeError(underflow, histo->bins_sum_w2().at(0), numEvents);
216 output << std::left << std::setw(13)
217 << "# Under: " << sformat("%-15g [1/pr] +/- %6.2f [%] \n", underflow, underflowError);
218
219 // Overflow
220 const G4double overflow =
221 static_cast<G4double>(histo->bins_sum_w().at(numBins + 1)) / numEvents;
222 const G4double overflowError =
223 computeError(overflow, histo->bins_sum_w2().at(numBins + 1), numEvents);
224 output << std::left << std::setw(13)
225 << "# Over: " << sformat("%-15g [1/pr] +/- %6.2f [%] \n", overflow, overflowError);
226 }
227
228 // END HEADER.
229 output << std::left << std::setw(13) << "# Log: " << binSchemeName << G4endl;
230 output << std::left << std::setw(13) << "# Histogram: " << numBins << " "
231 << minAbscissa / abscissaUnitValue << " " << maxAbscissa / abscissaUnitValue << " ["
232 << abscissaUnit << "]" << G4endl;
233
234 output << std::left << std::setw(13) << "# Column: "
235 << "1 " << abscissaName << "_min "
236 << "[" << abscissaUnit << "]" << G4endl;
237 output << std::left << std::setw(13) << "# Column: "
238 << "2 " << abscissaName << "_max "
239 << "[" << abscissaUnit << "]" << G4endl;
240 output << std::left << std::setw(13) << "# Column: "
241 << "3 dN/d(" << abscissaName << ") "
242 << "[1" << (abscissaUnit.empty() ? "" : ("/" + abscissaUnit)) << "/pr]" << G4endl;
243 if (!isProfile) {
244 output << std::left << std::setw(13) << "# Column: "
245 << "4 error "
246 << "[%]" << G4endl;
247 }
248
249 // HISTO / PROFILE DATA.
250 for (G4int binIndex = 0; binIndex < numBins; ++binIndex) {
251 const G4double mean = histo->bin_height(binIndex) / numEvents;
252 const G4double err = computeError(mean, histo->bins_sum_w2().at(binIndex + 1), numEvents);
253
254 // histo mode
255 if (!isProfile) {
256 output << sformat("%15g %15g %15g %6.2f\n",
257 histo->axis().bin_lower_edge(binIndex) / abscissaUnitValue,
258 histo->axis().bin_upper_edge(binIndex) / abscissaUnitValue,
259 mean / histo->axis().bin_width(binIndex) * abscissaUnitValue, err);
260 }
261 // profile mode
262 else {
263 output << sformat("%15g %15g %15g\n",
264 histo->axis().bin_lower_edge(binIndex) / abscissaUnitValue,
265 histo->axis().bin_upper_edge(binIndex) / abscissaUnitValue,
266 mean / histo->axis().bin_width(binIndex) * abscissaUnitValue);
267 }
268 } // loop on all bins
269 } // dumpG4H1InFlairFormat
270
271 } // namespace flair
272 } // namespace histo
273 } // namespace tools
274
275 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
276