BDSIM
BDSIM is a Geant4 extension toolkit for simulation of particle transport in accelerator beamlines.
HistogramAccumulatorMerge.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 "HistogramAccumulatorMerge.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 <string>
29
31
33 HistogramAccumulator(nullptr, 0, "", "")
34{;}
35
37 int nDimensionsIn,
38 const std::string& resultHistNameIn,
39 const std::string& resultHistTitleIn):
40 HistogramAccumulator(baseHistogram,
41 nDimensionsIn,
42 resultHistNameIn,
43 resultHistTitleIn)
44{;}
45
47{
48 // temporary variables
49 double newMean = 0;
50 double newVari = 0;
51 double var = 0;
52
53 // Want the number of events accumulated so far. We purposively set the entries
54 // in the mean histogram as the number of events accumulated, not the number of
55 // histograms (ie files here).
56 unsigned long oldEntries;
57 unsigned long newEntries;
58
59 if (nDimensions == 4)
60 {
61 oldEntries = (unsigned long)static_cast<BDSBH4DBase*>(mean)->GetEntries_BDSBH4D();
62 newEntries = (unsigned long)static_cast<BDSBH4DBase*>(newValue)->GetEntries_BDSBH4D();
63 }
64 else
65 {
66 oldEntries = (unsigned long)mean->GetEntries(); // works for base class*
67 newEntries = (unsigned long)newValue->GetEntries(); // works for base class*
68 }
69
70 unsigned long newTotalEntries = oldEntries + newEntries;
71 const double nD = (double)newEntries;
72 const double factor = nD * (nD - 1);
73
74 switch (nDimensions)
75 {
76 case 1:
77 {
78 TH1D* h1 = dynamic_cast<TH1D*>(mean);
79 TH1D* h1e = dynamic_cast<TH1D*>(variance);
80 TH1D* ht = dynamic_cast<TH1D*>(newValue);
81 for (int j = 0; j <= h1->GetNbinsX() + 1; ++j)
82 {
83 var = std::pow(ht->GetBinError(j), 2) * factor;
84 AccumulateSingleValue(h1->GetBinContent(j),
85 h1e->GetBinContent(j),
86 ht->GetBinContent(j),
87 var,
88 oldEntries, newEntries,
89 newMean, newVari);
90 h1->SetBinContent(j, newMean);
91 h1e->SetBinContent(j, newVari);
92 }
93 break;
94 }
95 case 2:
96 {
97 TH2D* h1 = dynamic_cast<TH2D*>(mean);
98 TH2D* h1e = dynamic_cast<TH2D*>(variance);
99 TH2D* ht = dynamic_cast<TH2D*>(newValue);
100 for (int j = 0; j <= h1->GetNbinsX() + 1; ++j)
101 {
102 for (int k = 0; k <= h1->GetNbinsY() + 1; ++k)
103 {
104 var = std::pow(ht->GetBinError(j,k), 2) * factor;
105 AccumulateSingleValue(h1->GetBinContent(j,k),
106 h1e->GetBinContent(j,k),
107 ht->GetBinContent(j,k),
108 var,
109 oldEntries, newEntries,
110 newMean, newVari);
111 h1->SetBinContent(j, k, newMean);
112 h1e->SetBinContent(j, k, newVari);
113 }
114 }
115 break;
116 }
117 case 3:
118 {
119 TH3D* h1 = dynamic_cast<TH3D*>(mean);
120 TH3D* h1e = dynamic_cast<TH3D*>(variance);
121 TH3D* ht = dynamic_cast<TH3D*>(newValue);
122 for (int j = 0; j <= h1->GetNbinsX() + 1; ++j)
123 {
124 for (int k = 0; k <= h1->GetNbinsY() + 1; ++k)
125 {
126 for (int l = 0; l <= h1->GetNbinsZ() + 1; ++l)
127 {
128 var = std::pow(ht->GetBinError(j,k,l), 2) * factor;
129 AccumulateSingleValue(h1->GetBinContent(j,k,l),
130 h1e->GetBinContent(j,k,l),
131 ht->GetBinContent(j,k,l),
132 var,
133 oldEntries, newEntries,
134 newMean, newVari);
135 h1->SetBinContent(j, k, l, newMean);
136 h1e->SetBinContent(j, k, l, newVari);
137 }
138 }
139 }
140 break;
141 }
142 case 4:
143 {
144#ifdef USE_BOOST
145 BDSBH4DBase* h1 = dynamic_cast<BDSBH4DBase*>(mean);
146 BDSBH4DBase* h1e = dynamic_cast<BDSBH4DBase*>(variance);
147 BDSBH4DBase* ht = dynamic_cast<BDSBH4DBase*>(newValue);
148 for (int j = -1; j <= h1->GetNbinsX(); ++j)
149 {
150 for (int k = -1; k <= h1->GetNbinsY(); ++k)
151 {
152 for (int l = -1; l <= h1->GetNbinsZ(); ++l)
153 {
154 for (int e = -1; e <= h1->GetNbinsE(); ++e)
155 {
156 var = std::pow(ht->AtError(j,k,l,e), 2) * factor;
157 AccumulateSingleValue(h1->At(j,k,l,e),
158 h1e->At(j,k,l,e),
159 ht->At(j,k,l,e),
160 var,
161 oldEntries, newEntries,
162 newMean, newVari);
163 h1->Set_BDSBH4D(j, k, l, e, newMean);
164 h1e->Set_BDSBH4D(j, k, l, e, newVari);
165 }
166 }
167 }
168 }
169 break;
170#endif
171 }
172 default:
173 {break;}
174 }
175 if(nDimensions==4)
176 {
177 dynamic_cast<BDSBH4DBase*>(mean)->SetEntries_BDSBH4D(newTotalEntries);
178 dynamic_cast<BDSBH4DBase*>(variance)->SetEntries_BDSBH4D(newTotalEntries);
179 }
180 else
181 {
182 mean->SetEntries(newTotalEntries);
183 variance->SetEntries(newTotalEntries);
184 }
185 n = newTotalEntries; // updated to Terminate() works correctly
186}
187
189 double oldVari,
190 double x,
191 double xVari,
192 unsigned long nEntriesAccumulated,
193 unsigned long nEntriesToAccumulate,
194 double& newMean,
195 double& newVari) const
196{
197 double dMean = x - oldMean;
198 double dMean2 = std::pow(dMean, 2);
199 double nA = (double)nEntriesAccumulated;
200 double nB = (double)nEntriesToAccumulate;
201 double nT = nA + nB;
202
203 newMean = oldMean + nB * (dMean / nT);
204
205 double q = (nA * nB) * (dMean2 / nT);
206 newVari = oldVari + xVari + q;
207}
Base class for the 4D histogram classes.
Definition: BDSBH4DBase.hh:37
Class to combine histograms of mean with error on mean.
virtual void AccumulateSingleValue(double oldMean, double oldVari, double x, double xVari, unsigned long nEntriesAccumulated, unsigned long nEntriesToAccumulate, double &newMean, double &newVari) const
virtual void Accumulate(TH1 *newValue)
HistogramAccumulatorMerge()
Default constructor only for ROOT reflexivity - not intended for use.
Class to accumulate and merge histograms in different ways.
int nDimensions
Number of dimensions.
unsigned long n
Counter.