gambit is hosted by Hepforge, IPPP Durham
GAMBIT  v1.5.0-2191-ga4742ac
a Global And Modular Bsm Inference Tool
AnalysisData.hpp
Go to the documentation of this file.
1 // GAMBIT: Global and Modular BSM Inference Tool
2 // *********************************************
26 
27 #pragma once
28 
30 #include "Eigen/Core"
32 
33 #include <string>
34 #include <map>
35 #include <iostream>
36 #include <sstream>
37 #include <vector>
38 #include <cmath>
39 #include <cfloat>
40 #include <limits>
41 #include <memory>
42 #include <iomanip>
43 #include <algorithm>
45 
46 // #define ANALYSISDATA_DEBUG
47 
48 #ifdef ANALYSISDATA_DEBUG
49 #include <iostream>
50 #endif
51 
52 namespace Gambit {
53  namespace ColliderBit {
54 
55 
58  {
59 
61  SignalRegionData(const EventCounter& scounter,
62  double nobs, const std::pair<double,double>& nbkg,
63  double nsigscaled=0)
64  : SignalRegionData(scounter.name(), nobs, scounter.weight_sum(), nbkg.first, scounter.weight_sum_err(), nbkg.second, nsigscaled)
65  {}
66 
68  SignalRegionData(const std::string& sr,
69  double nobs, const EventCounter& scounter, const std::pair<double,double>& nbkg,
70  double nsigscaled=0)
71  : SignalRegionData(sr, nobs, scounter.weight_sum(), nbkg.first, scounter.weight_sum_err(), nbkg.second, nsigscaled)
72  {}
73 
75  SignalRegionData(const std::string& sr,
76  double nobs, const std::pair<double,double>& nsigMC, const std::pair<double,double>& nbkg,
77  double nsigscaled=0)
78  : SignalRegionData(sr, nobs, nsigMC.first, nbkg.first, nsigMC.second, nbkg.second, nsigscaled)
79  {}
80 
82  SignalRegionData(const std::string& sr,
83  double nobs, double nsigMC, double nbkg,
84  double nsigMCsys, double nbkgerr, double nsigscaled=0)
85  : sr_label(sr),
86  n_obs(nobs), n_sig_MC(nsigMC), n_sig_scaled(nsigscaled), n_bkg(nbkg),
87  n_sig_MC_sys(nsigMCsys), n_bkg_err(nbkgerr)
88  {}
89 
92 
94  bool check() const {
95  bool consistent = true;
97  return consistent;
98  }
99 
101  double scalefactor() const { return n_sig_MC == 0 ? 1 : n_sig_scaled / n_sig_MC; }
102 
103  double calc_n_sig_MC_stat() const { return sqrt(n_sig_MC); }
104 
105  double calc_n_sig_MC_err() const
106  {
107  double n_sig_MC_stat = calc_n_sig_MC_stat();
108  return sqrt( n_sig_MC_stat * n_sig_MC_stat + n_sig_MC_sys * n_sig_MC_sys );
109  }
110 
111  double calc_n_sig_scaled_stat() const { return scalefactor() * calc_n_sig_MC_stat(); }
112 
113  double calc_n_sig_scaled_sys() const { return scalefactor() * n_sig_MC_sys; }
114 
115  double calc_n_sig_scaled_err() const { return scalefactor() * calc_n_sig_MC_err(); }
116 
117  double calc_n_sigbkg_err() const
118  {
119  double n_sig_scaled_err = calc_n_sig_scaled_err();
120  return sqrt( n_sig_scaled_err * n_sig_scaled_err + n_bkg_err * n_bkg_err );
121  }
122 
124 
126 
127  std::string sr_label;
128 
129 
131 
132  double n_obs = 0;
133  double n_sig_MC = 0;
134  double n_sig_scaled = 0;
135  double n_bkg = 0;
136  double n_sig_MC_sys = 0;
137  double n_bkg_err = 0;
138 
139 
140  };
141 
142 
149  {
150 
153  {
154  #ifdef ANALYSISDATA_DEBUG
155  std::cerr << "DEBUG: AnalysisData: " << this << " - Constructed (default ctor)" << std::endl;
156  #endif
157  clear();
158  }
159 
161  AnalysisData(const std::string& name) :
162  analysis_name(name)
163  {
164  #ifdef ANALYSISDATA_DEBUG
165  std::cerr << "DEBUG: AnalysisData: " << this << " - Constructed (ctor with analysis name)" << std::endl;
166  #endif
167  clear();
168  }
169 
170  // A copy constructor only used for debugging
171  #ifdef ANALYSISDATA_DEBUG
172  AnalysisData(const AnalysisData& copy) :
173  analysis_name(copy.analysis_name),
174  srdata(copy.srdata),
175  srdata_identifiers(copy.srdata_identifiers),
176  srcov(copy.srcov)
177  {
178  std::cerr << "DEBUG: AnalysisData: " << this << " - Copy-constructed from " << &copy << std::endl;
179  }
180  #endif
181 
182  // A destructor only used for debugging
183  #ifdef ANALYSISDATA_DEBUG
184  ~AnalysisData()
185  {
186  std::cerr << "DEBUG: AnalysisData: " << this << " - Destructed" << std::endl;
187  }
188  #endif
189 
194  AnalysisData(const std::vector<SignalRegionData>& srds, const Eigen::MatrixXd& cov=Eigen::MatrixXd())
195  : srdata(srds), srcov(cov)
196  {
197  #ifdef ANALYSISDATA_DEBUG
198  std::cerr << "DEBUG: AnalysisData: " << this << " - Constructed (special ctor)" << std::endl;
199  #endif
200  check();
201  }
202 
205  void clear()
206  {
207  for (auto& sr : srdata)
208  {
209  sr.n_sig_MC = 0;
210  sr.n_sig_scaled = 0;
211  sr.n_sig_MC_sys = 0;
212  }
213  srcov = Eigen::MatrixXd();
214  #ifdef ANALYSISDATA_DEBUG
215  std::cerr << "DEBUG: AnalysisData: " << this << " - Cleared" << std::endl;
216  #endif
217  }
218 
220  size_t size() const
221  {
222  // check();
223  return srdata.size();
224  }
225 
227  bool empty() const { return size() == 0; }
228 
230  bool hasCorrs() const
231  {
232  // check(); // bjf> This was wrong! Needs to be !=, not ==
233  return srcov.rows() != 0;
234  }
235 
238  void add(const SignalRegionData& srd)
239  {
240  std::string key = analysis_name + srd.sr_label;
241  auto loc = srdata_identifiers.find(key);
242  if (loc == srdata_identifiers.end())
243  {
244  // If the signal region doesn't exist in this object yet, add it
245  srdata.push_back(srd);
246  srdata_identifiers[key] = srdata.size() - 1;
247  }
248  else
249  {
250  // If it does, just update the signal count in the existing SignalRegionData object
251  srdata[loc->second].n_sig_MC = srd.n_sig_MC;
252  }
253  check();
254  }
255 
257  bool check() const
258  {
259  for (const SignalRegionData& srd : srdata) srd.check();
260  assert(srcov.rows() == 0 || srcov.rows() == (int) srdata.size());
261  // for (int isr = 0; isr < srcov.rows(); ++isr) {
262  // const double& srbg = srdata[isr].n_bkg_err;
263  // #ifdef ANALYSISDATA_DEBUG
264  // std::cerr << "DEBUG: AnalysisData: isr:" << isr << ", srbg:" << srbg << ", srcov(isr,isr):" << srcov(isr,isr) << std::endl;
265  // #endif
266  // assert(fabs(srcov(isr,isr) - srbg*srbg) < 1e-2);
267  // }
268  return true;
269  }
270 
273  void pythonize_me() const;
274 
276  std::string analysis_name;
277 
279  SignalRegionData& operator[] (size_t i) { return srdata[i]; }
281  const SignalRegionData& operator[] (size_t i) const { return srdata[i]; }
282 
284  std::vector<SignalRegionData>::iterator begin() { return srdata.begin(); }
285  std::vector<SignalRegionData>::const_iterator begin() const { return srdata.begin(); }
286  std::vector<SignalRegionData>::iterator end() { return srdata.end(); }
287  std::vector<SignalRegionData>::const_iterator end() const { return srdata.end(); }
288 
290  std::vector<SignalRegionData> srdata;
291 
293  std::map<std::string, int> srdata_identifiers;
294 
296  Eigen::MatrixXd srcov;
297 
298  };
299 
300 
301  }
302 }
SignalRegionData(const std::string &sr, double nobs, double nsigMC, double nbkg, double nsigMCsys, double nbkgerr, double nsigscaled=0)
Constructor with separate n & nsys args.
double n_sig_MC
The number of simulated model events passing selection for this signal region.
bool check() const
Consistency check.
Eigen::MatrixXd srcov
Optional covariance matrix between SRs (0x0 null matrix = no correlation info)
EventCounter class.
AnalysisData(const std::string &name)
Constructor with analysis name.
SignalRegionData(const std::string &sr, double nobs, const std::pair< double, double > &nsigMC, const std::pair< double, double > &nbkg, double nsigscaled=0)
Constructor with {n,nsys} pair args.
std::vector< SignalRegionData >::iterator begin()
Iterators (sugar for direct access to this->srdata)
double n_bkg_err
The absolute error of n_bkg.
void add(const SignalRegionData &srd)
Add a SignalRegionData.
double n_bkg
The number of standard model events expected to pass the selection for this signal region...
SignalRegionData(const EventCounter &scounter, double nobs, const std::pair< double, double > &nbkg, double nsigscaled=0)
Constructor with EventCounter arg for the signal count and SR name.
A container for the result of an analysis, potentially with many signal regions and correlations...
double n_obs
The number of events passing selection for this signal region as reported by the experiment.
std::vector< SignalRegionData >::const_iterator end() const
SignalRegionData(const std::string &sr, double nobs, const EventCounter &scounter, const std::pair< double, double > &nbkg, double nsigscaled=0)
Constructor with EventCounter arg for the signal count, but separate name.
double scalefactor() const
Uncertainty calculators.
std::string sr_label
A label for the particular signal region of the analysis.
A simple class for counting events of type HEPUtils::Event.
AnalysisData()
Default constructor.
std::map< std::string, int > srdata_identifiers
Map of names and indices of all entries in srdata, for easy lookup.
std::string analysis_name
Analysis name.
double n_sig_MC_sys
The absolute systematic error of n_sig_MC.
Simple header file for turning compiler warnings back on after having included one of the begin_ignor...
std::vector< SignalRegionData > srdata
List of signal regions&#39; data summaries.
A simple container for the result of one signal region from one analysis.
std::vector< SignalRegionData >::iterator end()
DS5_MSPCTM DS_INTDOF int
Pragma directives to suppress compiler warnings coming from including Eigen library headers...
void clear()
Clear the list of SignalRegionData, and nullify the covariance matrix.
SignalRegionData()
Default constructor.
bool hasCorrs() const
Is there non-null correlation data?
double n_sig_scaled
n_sig_MC, scaled to luminosity * cross-section
bool check() const
Check that the SRData list and the covariance matrix are consistent.
bool empty() const
Is this container empty of signal regions?
std::vector< SignalRegionData >::const_iterator begin() const
AnalysisData(const std::vector< SignalRegionData > &srds, const Eigen::MatrixXd &cov=Eigen::MatrixXd())
Constructor from a list of SignalRegionData and an optional correlation (or covariance?) matrix.
TODO: see if we can use this one:
Definition: Analysis.hpp:33
size_t size() const
Number of analyses.