BDSIM
BDSIM is a Geant4 extension toolkit for simulation of particle transport in accelerator beamlines.
HistogramAccumulator.cc
1/*
2Beam Delivery Simulation (BDSIM) Copyright (C) Royal Holloway,
3University of London 2001 - 2022.
4
5This file is part of BDSIM.
6
7BDSIM is free software: you can redistribute it and/or modify
8it under the terms of the GNU General Public License as published
9by the Free Software Foundation version 3 of the License.
10
11BDSIM is distributed in the hope that it will be useful, but
12WITHOUT ANY WARRANTY; without even the implied warranty of
13MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
14GNU General Public License for more details.
15
16You should have received a copy of the GNU General Public License
17along with BDSIM. If not, see <http://www.gnu.org/licenses/>.
18*/
19#include "HistogramAccumulator.hh"
20
21#include "TH1.h"
22#include "TH1D.h"
23#include "TH2D.h"
24#include "TH3D.h"
25#include "BDSBH4DBase.hh"
26
27#include <cmath>
28#include <stdexcept>
29#include <string>
30
32
34 nDimensions(1),
35 n(0),
36 terminated(false),
37 mean(nullptr),
38 variance(nullptr),
39 result(nullptr)
40{;}
41
43 const std::string& resultHistNameIn,
44 const std::string& resultHistTitleIn):
45 HistogramAccumulator(baseHistogram, (int)baseHistogram->GetDimension(), resultHistNameIn, resultHistTitleIn)
46{;}
47
49 int nDimensionsIn,
50 const std::string& resultHistNameIn,
51 const std::string& resultHistTitleIn):
52 nDimensions(nDimensionsIn),
53 n(0),
54 terminated(false),
55 resultHistName(resultHistNameIn),
56 resultHistTitle(resultHistTitleIn),
57 mean(nullptr),
58 variance(nullptr),
59 result(nullptr)
60{
61 std::string meanName = resultHistName + "_Mean";
62 std::string variName = resultHistName + "_Vari";
63 switch (nDimensions)
64 {
65 case 1:
66 {
67 mean = dynamic_cast<TH1D*>(baseHistogram->Clone(meanName.c_str()));
68 variance = dynamic_cast<TH1D*>(baseHistogram->Clone(variName.c_str()));
69 result = dynamic_cast<TH1D*>(baseHistogram->Clone(resultHistName.c_str()));
70 break;
71 }
72 case 2:
73 {
74 mean = dynamic_cast<TH2D*>(baseHistogram->Clone(meanName.c_str()));
75 variance = dynamic_cast<TH2D*>(baseHistogram->Clone(variName.c_str()));
76 result = dynamic_cast<TH2D*>(baseHistogram->Clone(resultHistName.c_str()));
77 break;
78 }
79 case 3:
80 {
81 mean = dynamic_cast<TH3D*>(baseHistogram->Clone(meanName.c_str()));
82 variance = dynamic_cast<TH3D*>(baseHistogram->Clone(variName.c_str()));
83 result = dynamic_cast<TH3D*>(baseHistogram->Clone(resultHistName.c_str()));
84 break;
85 }
86 case 4:
87 {
88#ifdef USE_BOOST
89 mean = dynamic_cast<BDSBH4DBase*>(baseHistogram)->Clone(meanName.c_str());
90 variance = dynamic_cast<BDSBH4DBase*>(baseHistogram)->Clone(variName.c_str());
91 result = dynamic_cast<BDSBH4DBase*>(baseHistogram)->Clone(resultHistName.c_str());
92 break;
93#endif
94 }
95 default:
96 {throw std::domain_error("Invalid number of dimensions"); break;}
97 }
98 if (mean && variance && result)
99 {// technically these could be nullptr
100 if (nDimensions != 4)
101 {
102 // empty contents
103 mean->Reset();
104 variance->Reset();
105 result->Reset();
106 // set title
107 result->SetTitle(resultHistTitle.c_str());
108 mean->SetTitle(meanName.c_str());
109 variance->SetTitle(variName.c_str());
110 }
111 else
112 {
113#ifdef USE_BOOST
114 // empty contents
115 static_cast<BDSBH4DBase*>(mean)->Reset_BDSBH4D();
116 static_cast<BDSBH4DBase*>(variance)->Reset_BDSBH4D();
117 static_cast<BDSBH4DBase*>(result)->Reset_BDSBH4D();
118 // set title
119 static_cast<BDSBH4DBase*>(result)->SetTitle(resultHistTitle.c_str());
120 static_cast<BDSBH4DBase*>(mean)->SetTitle(meanName.c_str());
121 static_cast<BDSBH4DBase*>(variance)->SetTitle(variName.c_str());
122#endif
123 }
124 }
125}
126
128{
129 // deleting histograms removes them from currently open output file
130 delete mean;
131 delete variance;
132 // leak result here as ROOT annoyingly requires this to be left
133}
134
136{
137 // temporary variables
138 double newMean = 0;
139 double newVari = 0;
140 const double error = 0; // needed to pass reference to unused parameter
141 const unsigned long nEntriesToAccumulate = 1;
142
143 n++;
144 switch (nDimensions)
145 {
146 case 1:
147 {
148 TH1D* h1 = dynamic_cast<TH1D*>(mean);
149 TH1D* h1e = dynamic_cast<TH1D*>(variance);
150 TH1D* ht = dynamic_cast<TH1D*>(newValue);
151 for (int j = 0; j <= h1->GetNbinsX() + 1; ++j)
152 {
153 AccumulateSingleValue(h1->GetBinContent(j),
154 h1e->GetBinContent(j),
155 ht->GetBinContent(j),
156 error, n, nEntriesToAccumulate,
157 newMean, newVari);
158 h1->SetBinContent(j, newMean);
159 h1e->SetBinContent(j, newVari);
160 }
161 break;
162 }
163 case 2:
164 {
165 TH2D* h1 = dynamic_cast<TH2D*>(mean);
166 TH2D* h1e = dynamic_cast<TH2D*>(variance);
167 TH2D* ht = dynamic_cast<TH2D*>(newValue);
168 for (int j = 0; j <= h1->GetNbinsX() + 1; ++j)
169 {
170 for (int k = 0; k <= h1->GetNbinsY() + 1; ++k)
171 {
172 AccumulateSingleValue(h1->GetBinContent(j,k),
173 h1e->GetBinContent(j,k),
174 ht->GetBinContent(j,k),
175 error, n, nEntriesToAccumulate,
176 newMean, newVari);
177 h1->SetBinContent(j, k, newMean);
178 h1e->SetBinContent(j, k, newVari);
179 }
180 }
181 break;
182 }
183 case 3:
184 {
185 TH3D* h1 = dynamic_cast<TH3D*>(mean);
186 TH3D* h1e = dynamic_cast<TH3D*>(variance);
187 TH3D* ht = dynamic_cast<TH3D*>(newValue);
188 for (int j = 0; j <= h1->GetNbinsX() + 1; ++j)
189 {
190 for (int k = 0; k <= h1->GetNbinsY() + 1; ++k)
191 {
192 for (int l = 0; l <= h1->GetNbinsZ() + 1; ++l)
193 {
194 AccumulateSingleValue(h1->GetBinContent(j,k,l),
195 h1e->GetBinContent(j,k,l),
196 ht->GetBinContent(j,k,l),
197 error, n, nEntriesToAccumulate,
198 newMean, newVari);
199 h1->SetBinContent(j, k, l, newMean);
200 h1e->SetBinContent(j, k, l, newVari);
201 }
202 }
203 }
204 break;
205 }
206 case 4:
207 {
208#ifdef USE_BOOST
209 BDSBH4DBase* h1 = dynamic_cast<BDSBH4DBase*>(mean);
210 BDSBH4DBase* h1e = dynamic_cast<BDSBH4DBase*>(variance);
211 BDSBH4DBase* ht = dynamic_cast<BDSBH4DBase*>(newValue);
212 for (int j = -1; j <= h1->GetNbinsX(); ++j)
213 {
214 for (int k = -1; k <= h1->GetNbinsY(); ++k)
215 {
216 for (int l = -1; l <= h1->GetNbinsZ(); ++l)
217 {
218 for (int e = -1; e <= h1->GetNbinsE(); ++e)
219 {
220 AccumulateSingleValue(h1->At(j,k,l,e),
221 h1e->At(j,k,l,e),
222 ht->At(j,k,l,e),
223 error, n, nEntriesToAccumulate,
224 newMean, newVari);
225 h1->Set_BDSBH4D(j,k,l,e, newMean);
226 h1e->Set_BDSBH4D(j,k,l,e, newVari);
227 }
228
229 }
230 }
231 }
232 break;
233#endif
234 }
235 default:
236 {break;}
237 }
238}
239
241{
242 // error on mean is sqrt(1/n) * std = sqrt(1/n) * sqrt(1/(n-1)) * sqrt(variance)
243 // the only variable is the variance, so take the rest out as a factor.
244 const double nD = (double)n; // cast only once
245 const double factor = std::sqrt(1./(nD * (nD - 1))); // nan if n = 1 -> won't be used
246 double mn = 0; // temporary variable for mean
247 double err = 0; // temporary variable for standard error on mean
248 double var = 0; // temporary variable for variance
249
250 // note here we set the std to 0 if there's only one entry (ie n = 1) to avoid
251 // division by zero and nans
252 switch (nDimensions)
253 {
254 case 1:
255 {
256 for (int j = 0; j <= result->GetNbinsX() + 1; ++j)
257 {
258 mn = mean->GetBinContent(j);
259 var = variance->GetBinContent(j);
260 err = n > 1 ? factor*std::sqrt(var) : 0;
261 result->SetBinContent(j, mn);
262 result->SetBinError(j, err);
263 }
264 break;
265 }
266 case 2:
267 {
268 for (int j = 0; j <= result->GetNbinsX() + 1; ++j)
269 {
270 for (int k = 0; k <= result->GetNbinsY() + 1; ++k)
271 {
272 mn = mean->GetBinContent(j,k);
273 var = variance->GetBinContent(j, k);
274 err = n > 1 ? factor*std::sqrt(var) : 0;
275 result->SetBinContent(j, k, mn);
276 result->SetBinError(j, k, err);
277 }
278 }
279 break;
280 }
281 case 3:
282 {
283 for (int j = 0; j <= result->GetNbinsX() + 1; ++j)
284 {
285 for (int k = 0; k <= result->GetNbinsY() + 1; ++k)
286 {
287 for (int l = 0; l <= result->GetNbinsZ() + 1; ++l)
288 {
289 mn = mean->GetBinContent(j,k,l);
290 var = variance->GetBinContent(j, k, l);
291 err = n > 1 ? factor*std::sqrt(var) : 0;
292 result->SetBinContent(j,k,l, mn);
293 result->SetBinError(j,k,l, err);
294 }
295 }
296 }
297 break;
298 }
299 case 4:
300 {
301#ifdef USE_BOOST
302 auto histCast = dynamic_cast<BDSBH4DBase*>(result);
303 auto mnCast = dynamic_cast<BDSBH4DBase*>(mean);
304 auto varCast = dynamic_cast<BDSBH4DBase*>(variance);
305 auto resCast = dynamic_cast<BDSBH4DBase*>(result);
306 int nBinsX = histCast->GetNbinsX();
307 int nBinsY = histCast->GetNbinsY();
308 int nBinsZ = histCast->GetNbinsZ();
309 int nBinsE = histCast->GetNbinsE();
310 for (int j = -1; j <= nBinsX; ++j)
311 {
312 for (int k = -1; k <= nBinsY; ++k)
313 {
314 for (int l = -1; l <= nBinsZ; ++l)
315 {
316 for (int e = -1; e <= nBinsE; ++e)
317 {
318 mn = mnCast->At(j, k, l, e);
319 var = varCast->At(j, k, l, e);
320 err = n > 1 ? factor*std::sqrt(var) : 0;
321 resCast->Set_BDSBH4D(j, k, l, e, mn);
322 resCast->SetError_BDSBH4D(j, k, l, e, err);
323 }
324 }
325 }
326 }
327 break;
328#endif
329 }
330 default:
331 {break;}
332 }
333 if(nDimensions==4)
334 {dynamic_cast<BDSBH4DBase*>(result)->SetEntries_BDSBH4D((double)n);}
335 else
336 {result->SetEntries((double)n);}
337
338 return result;
339}
340
342 double oldVari,
343 double x,
344 double /*xVari*/,
345 unsigned long nEntriesAccumulated,
346 unsigned long/*nEntriesToAccumulate*/,
347 double& newMean,
348 double& newVari) const
349{
350 newMean = oldMean + ((x - oldMean) / (double)nEntriesAccumulated);
351 newVari = oldVari + ((x - oldMean) * (x - newMean));
352}
Base class for the 4D histogram classes.
Definition: BDSBH4DBase.hh:37
Class to accumulate and merge histograms in different ways.
const std::string resultHistName
Name for resultant histogram.
int nDimensions
Number of dimensions.
virtual void AccumulateSingleValue(double oldMean, double oldVari, double x, double xVari, unsigned long nEntriesAccumulated, unsigned long nEntriesToAccumulate, double &newMean, double &newVari) const
const std::string resultHistTitle
Title for resultant histogram.
unsigned long n
Counter.
HistogramAccumulator()
Default constructor only for ROOT reflexivity - not intended for use.
virtual void Accumulate(TH1 *newValue)