gambit is hosted by Hepforge, IPPP Durham
GAMBIT  v1.5.0-2191-ga4742ac
a Global And Modular Bsm Inference Tool
SimpleHist.cpp
Go to the documentation of this file.
1 // GAMBIT: Global and Modular BSM Inference Tool
2 // *********************************************
17 
21 
22 //#define DARKBIT_DEBUG
23 
24 namespace Gambit {
25  namespace DarkBit {
26 
27  SimpleHist::SimpleHist(int nBins, double Emin, double Emax, bool logscale):
28  nBins(nBins)
29  {
30 #ifdef DARKBIT_DEBUG
31  std::cout << "Emin: " << Emin << std::endl;
32  std::cout << "Emax: " << Emax << std::endl;
33  std::cout << "nBins: " << nBins << std::endl;
34 #endif
35  if(logscale)
36  {
37  if(Emin<=0)
38  {
39  DarkBit_error().raise(LOCAL_INFO,
40  "Lower histogram bin limit must be greater than 0 when using\n"
41  "logarithmic binning.");
42  }
43  double factor = pow(Emax/Emin,(1.0/nBins));
44  double binL=Emin;
45  for(int i=0; i<nBins; i++)
46  {
47 #ifdef DARKBIT_DEBUG
48  std::cout << binL << std::endl;
49 #endif
50  binLower.push_back(binL);
51  binVals.push_back(0.0);
52  wtSq.push_back(0.0);
53  binL*=factor;
54  }
55  binLower.push_back(binL);
56  }
57  else
58  {
59  double dE = (Emax-Emin)/nBins;
60  for(int i=0; i<nBins; i++)
61  {
62  double binL = Emin+i*dE;
63  binLower.push_back(binL);
64  binVals.push_back(0.0);
65  wtSq.push_back(0.0);
66  }
67  binLower.push_back(Emin+nBins*dE);
68  }
69  }
70 
71  SimpleHist::SimpleHist(std::vector<double> binLower) : binLower(binLower)
72  {
73  nBins=int(binLower.size())-1;
74  binVals=std::vector<double>(nBins,0.0);
75  wtSq =std::vector<double>(nBins,0.0);
76  }
77 
78  void SimpleHist::addEvent(double E, double weight)
79  {
80  int bin = findIndex(E);
81  if(bin>=0 and bin<nBins )
82  {
83  binVals[bin]+=weight;
84  wtSq[bin]+=weight*weight;
85  }
86  }
87 
88  void SimpleHist::addToBin(int bin, double weight)
89  {
90  binVals[bin]+=weight;
91  wtSq[bin]+=weight*weight;
92  }
93 
94  void SimpleHist::addBox(double Emin, double Emax, double weight)
95  {
96  int imin = findIndex(Emin);
97  int imax = findIndex(Emax);
98  double dE = Emax-Emin;
99  double norm = weight/dE;
100  if(imax<0 or imin>=nBins)
101  {
102  // Do nothing
103  }
104  else if(imin==imax)
105  {
106  addToBin(imin,weight);
107  }
108  else if(imin<0 and imax==0)
109  {
110  addToBin(0,(Emax-binLower[0])*norm);
111  }
112  else if(imin==(nBins-1) and imax>=nBins)
113  {
114  addToBin(nBins-1,(binLower[nBins]-Emin)*norm);
115  }
116  else
117  {
118  double binSize_low,binSize_high;
119  // Calculate part of lower bin covered by box
120  if(imin>=0)
121  binSize_low = binLower[imin+1]-Emin;
122  else
123  {
124  imin = 0;
125  binSize_low=binSize(imin);
126  }
127  // Calculate part of upper bin covered by box
128  if(imax<nBins)
129  binSize_high = Emax-binLower[imax];
130  else
131  {
132  imax = nBins-1;
133  binSize_high=binSize(imax);
134  }
135  // Add contribution to lower bin
136  addToBin(imin,binSize_low*norm);
137  // Add contribution to upper bin
138  addToBin(imax,binSize_high*norm);
139  // Add contributions to remaining bins
140  for(int i=imin+1;i<imax;i++)
141  {
142  addToBin(i,binSize(i)*norm);
143  }
144  }
145  }
146 
148  {
149  // Check that the number of bins is equal to avoid segfaults.
150  // It is up to the user to make sure the actual binning is identical.
151  if(in.nBins != nBins)
152  {
153  DarkBit_error().raise(LOCAL_INFO,
154  "SimpleHist::addHistAsWeights_sameBin requires identically binned\n"
155  "histograms.");
156  }
157  for(int i=0; i<nBins;i++)
158  {
159  addToBin(i,in.binVals[i]);
160  }
161  }
162 
163  double SimpleHist::getError(int bin) const
164  {
165  return sqrt(wtSq[bin]);
166  }
167 
168  double SimpleHist::getRelError(int bin) const
169  {
170  return ((binVals[bin]!=0) ? sqrt(wtSq[bin])/binVals[bin] : 0);
171  }
172 
174  {
175  for(int i=0;i<nBins;i++)
176  {
177  binVals[i]/=binSize(i);
178  wtSq[i] /=(binSize(i)*binSize(i));
179  }
180  }
181 
182  void SimpleHist::multiply(double x)
183  {
184  for(int i=0;i<nBins;i++)
185  {
186  binVals[i]*=x;
187  wtSq[i] *=x*x;
188  }
189  }
190 
191  int SimpleHist::findIndex(double val) const
192  {
193  if(val < binLower[0]) return -1;
194  std::vector<double>::const_iterator pos =
195  upper_bound(binLower.begin(),binLower.end(),val);
196  return pos - binLower.begin() -1;
197  }
198 
199  double SimpleHist::binSize(int bin) const
200  {
201  return binLower[bin+1]-binLower[bin];
202  }
203 
204  double SimpleHist::binCenter(int bin) const
205  {
206  return 0.5*(binLower[bin+1]+binLower[bin]);
207  }
208 
209  std::vector<double> SimpleHist::getBinCenters() const
210  {
211  std::vector<double> centers;
212  for(int i=0;i<nBins;i++)
213  {
214  centers.push_back(binCenter(i));
215  }
216  return centers;
217  }
218 
219  const std::vector<double>& SimpleHist::getBinValues() const
220  {
221  return binVals;
222  }
223 
224  void SimpleHist::getEdges(double& lower, double& upper) const
225  {
226  lower=binLower[0];
227  upper=binLower[nBins];
228  }
229 
230  } // namespace DarkBit
231 
232 } // namespace Gambit
233 #undef DARKBIT_DEBUG
error & DarkBit_error()
void divideByBinSize()
Divide all histogram bins by the respective bin size.
Definition: SimpleHist.cpp:173
Type definition header for DarkBit SimpleHist types.
void addHistAsWeights_sameBin(SimpleHist &in)
Add content of input histogram as weights.
Definition: SimpleHist.cpp:147
double getRelError(int bin) const
Get relative error for a specified bin.
Definition: SimpleHist.cpp:168
std::vector< double > binVals
Definition: SimpleHist.hpp:93
int findIndex(double val) const
Find bin index for given value.
Definition: SimpleHist.cpp:191
Histogram class for cascade decays.
Definition: SimpleHist.hpp:37
#define LOCAL_INFO
Definition: local_info.hpp:34
double binCenter(int bin) const
Get central value of bin.
Definition: SimpleHist.cpp:204
std::vector< double > wtSq
Sum of the squares of all weights.
Definition: SimpleHist.hpp:95
double binSize(int bin) const
Retrieve size of given bins.
Definition: SimpleHist.cpp:199
double getError(int bin) const
Get error for a specified bin.
Definition: SimpleHist.cpp:163
const std::vector< double > & getBinValues() const
Definition: SimpleHist.cpp:219
void addEvent(double E, double weight=1.0)
Add an entry to histogram.
Definition: SimpleHist.cpp:78
void addBox(double Emin, double Emax, double weight=1.0)
Add a box spectrum to the histogram.
Definition: SimpleHist.cpp:94
Header file that includes all GAMBIT headers required for a module source file.
void multiply(double x)
Multiply all bin contents by x.
Definition: SimpleHist.cpp:182
void addToBin(int bin, double weight=1.0)
Add an entry to a specified bin.
Definition: SimpleHist.cpp:88
DS5_MSPCTM DS_INTDOF int
std::vector< double > getBinCenters() const
Double get central values of all bins.
Definition: SimpleHist.cpp:209
std::vector< double > binLower
Definition: SimpleHist.hpp:92
Rollcall header for module DarkBit.
double pow(const double &a)
Outputs a^i.
void getEdges(double &lower, double &upper) const
Definition: SimpleHist.cpp:224
TODO: see if we can use this one:
Definition: Analysis.hpp:33
double E(const double x, const double y)