a Code for the Combination of Indirect and Direct Constraints on High Energy Physics Models Logo
Observable.cpp
Go to the documentation of this file.
1 /*
2  * Copyright (C) 2012 HEPfit Collaboration
3  *
4  *
5  * For the licensing terms see doc/COPYING.
6  */
7 
8 #include "Observable.h"
9 #include "StandardModel.h"
10 #include <TNamed.h>
11 #include <TFile.h>
12 #include <TROOT.h>
13 #include <TMath.h>
14 #include <TH1D.h>
15 #include <limits>
16 
17 
18 Observable::Observable (const std::string name_i,
19  const std::string thname_i,
20  const std::string label_i,
21  const bool tMCMC_i,
22  const double min_i,
23  const double max_i,
24  ThObservable * tho_i)
25 {
26  name = name_i;
27  thname = thname_i;
28  label = label_i;
29  min = min_i;
30  max = max_i;
31  tMCMC = tMCMC_i;
32  tho = tho_i;
33  distr = "";
34  filename = "";
35  histoname = "";
36  ave = 0.;
37  errg = 0.;
38  errf = 0.;
39  errgl = 0.;
40  errgr = 0.;
41  obsType = "";
42  bin_min = 0.;
43  bin_max = 0.;
44  iterationNo = std::numeric_limits<int>::max();
45  writeChain = false;
46  hasInverseCovariance = false;
47 }
48 
49 Observable::Observable(const Observable& orig)
50 {
51  name = orig.name;
52  thname = orig.thname;
53  label = orig.label;
54  min = orig.min;
55  max = orig.max;
56  tMCMC = orig.tMCMC;
57  tho = orig.tho;
58  distr = orig.distr;
59  filename = orig.filename;
60  histoname = orig.histoname;
61  ave = orig.ave;
62  errg = orig.errg;
63  errf = orig.errf;
64  errgl = orig.errgl;
65  errgr = orig.errgr;
66  obsType = orig.obsType;
67  bin_min = orig.bin_min;
68  bin_min = orig.bin_max;
69  iterationNo = orig.iterationNo;
70  inhisto = orig.inhisto;
71  writeChain = orig.writeChain;
72  hasInverseCovariance = orig.hasInverseCovariance;
73 }
74 
75 Observable::Observable()
76 {
77  name = "";
78  thname = "";
79  label = "";
80  min = 0.;
81  max = 0.;
82  tMCMC = false;
83  tho = NULL;
84  distr = "";
85  filename = "";
86  histoname = "";
87  ave = 0.;
88  errg = 0.;
89  errf = 0.;
90  errgl = 0.;
91  errgr = 0.;
92  obsType = "";
93  bin_min = 0.;
94  bin_max = 0.;
95  iterationNo = std::numeric_limits<int>::max();
96  writeChain = false;
97  hasInverseCovariance = false;
98 }
99 
100 Observable::~Observable() {}
101 
102 std::ostream& operator<<(std::ostream& output, const Observable& o)
103 {
104  output << "Observable name, tMCMC, min, max, distribution, distribution parameters" << std::endl;
105  output << o.name << " " << o.tMCMC << " " << o.min << " " << o.max << " "
106  << o.distr << " " << o.filename << " " << o.histoname << " " << o.ave
107  << " " << o.errg << " " << o.errf << std::endl;
108  return output;
109 }
110 
111 void Observable::setLikelihoodFromHisto(std::string filename, std::string histoname)
112  {
113  this->filename = filename;
114  this->histoname = histoname;
115  TFile *lik = new TFile((filename + ".root").c_str(), "read");
116  TH1D *htmp = (TH1D*) (lik->Get(histoname.c_str()));
117  if (htmp == NULL)
118  throw std::runtime_error("ERROR: nonexistent histogram called "
119  + histoname + " in "
120  + filename + ".root");
121  inhisto = (TH1D *) htmp->Clone((filename + "/" + histoname).c_str());
122  inhisto->SetDirectory(gROOT);
123  setMin(inhisto->GetXaxis()->GetXmin());
124  setMax(inhisto->GetXaxis()->GetXmax());
125  lik->Close();
126  delete lik;
127  }
128 
129 
130 double Observable::computeTheoryValue()
131 {
132  if (tho->getModel().getIterationNo() == iterationNo) {
133  return thValue;
134  } else {
135  iterationNo = tho->getModel().getIterationNo();
136  thValue = tho->computeThValue();
137  return thValue;
138  }
139 }
140 
141 double Observable::LogSplitGaussian(double x, double ave, double errl, double errr)
142 {
143  double sigma = (x > ave ? errr : errl);
144  return -0.5 * (x-ave) * (x-ave) / sigma / sigma;
145 }
146 
147 double Observable::LogGaussian(double x, double ave, double sigma)
148 {
149  return -0.5 * (x-ave) * (x-ave) / sigma / sigma;
150 }
151 
152 double Observable::computeWeight(double th)
153 {
154  double logprob;
155  if (distr.compare("weight") == 0) {
156  if (getObsType().compare("AsyGausObservable") == 0) logprob = LogSplitGaussian(th, ave, errgl, errgr);
157  else {
158  if (errf == 0.)
159  logprob = LogGaussian(th, ave, errg);
160  else if (errg == 0.) {
161  if (th < ave + errf && th > ave - errf)
162  logprob = 1.;
163  else
164  logprob = log(0.);
165  } else
166  logprob = log(TMath::Erf((th - ave + errf) / sqrt(2.) / errg)
167  - TMath::Erf((th - ave - errf) / sqrt(2.) / errg));
168  }
169  } else if (distr.compare("file") == 0) {
170  int i = inhisto->FindBin(th);
171  if (inhisto->IsBinOverflow(i) || inhisto->IsBinUnderflow(i))
172  logprob = log(0.);
173  else
174  logprob = log(inhisto->GetBinContent(i));
175  //logprob = log(h->GetBinContent(h->FindBin(th)));
176  } else
177  throw std::runtime_error("ERROR: MonteCarloEngine::Weight() called without weight for "
178  + name);
179  return (logprob);
180 }
181 
182 double Observable::computeWeight(double th, double ave_i, double errg_i, double errf_i)
183 {
184  double logprob;
185  if (distr.compare("weight") == 0) {
186  if (errf_i == 0.)
187  logprob = LogGaussian(th, ave_i, errg_i);
188  else if (errg_i == 0.) {
189  if (th < ave_i + errf_i && th > ave_i - errf_i)
190  logprob = 1.;
191  else
192  logprob = log(0.);
193  } else
194  logprob = log(TMath::Erf((th - ave_i + errf_i) / sqrt(2.) / errg_i)
195  - TMath::Erf((th - ave_i - errf_i) / sqrt(2.) / errg_i));
196  } else
197  throw std::runtime_error("ERROR: MonteCarloEngine::Weight() called without weight for "
198  + name);
199  return (logprob);
200 }
201 
202 boost::tokenizer<boost::char_separator<char> >::iterator & Observable::ParseObservable(std::string& type,
203  boost::tokenizer<boost::char_separator<char> >* tok,
204  boost::tokenizer<boost::char_separator<char> >::iterator & beg,
205  std::string& filepath,
206  std::string& infilename,
207  int rank)
208 {
209  if ((type.compare("Observable") == 0 || type.compare("HiggsObservable")) && std::distance(tok->begin(), tok->end()) < 8) {
210  if(rank == 0) throw std::runtime_error("ERROR: lack of information on " + *beg + " in " + infilename);
211  else sleep (2);
212  } else if (type.compare("BinnedObservable") == 0 && std::distance(tok->begin(), tok->end()) < 12) {
213  if(rank == 0) throw std::runtime_error("ERROR: lack of information on " + *beg + " in " + infilename);
214  else sleep (2);
215  } else if (type.compare("FunctionObservable") == 0 && std::distance(tok->begin(), tok->end()) < 11) {
216  if(rank == 0) throw std::runtime_error("ERROR: lack of information on " + *beg + " in " + infilename);
217  else sleep (2);
218  } else {
219  obsType = type;
220  name = *beg;
221  ++beg;
222  thname = *beg;
223  ++beg;
224  label = *beg;
225  size_t pos = 0;
226  while ((pos = label.find("~", pos)) != std::string::npos)
227  label.replace(pos++, 1, " ");
228  ++beg;
229  min = atof((*beg).c_str());
230  ++beg;
231  max = atof((*beg).c_str());
232  ++beg;
233  std::string toMCMC = *beg;
234  if (toMCMC.compare("MCMC") == 0)
235  tMCMC = true;
236  else if (toMCMC.compare("noMCMC") == 0)
237  tMCMC = false;
238  else {
239  if (rank == 0) throw std::runtime_error("ERROR: wrong MCMC flag in " + name + ".\n");
240  else sleep (2);
241  }
242 
243  if (obsType.compare("Observable") == 0 || obsType.compare("BinnedObservable") == 0 || obsType.compare("FunctionObservable") == 0 || type.compare("AsyGausObservable") == 0) {
244  ++beg;
245  distr = *beg;
246  if (distr.compare("file") == 0) {
247  if (std::distance(tok->begin(), tok->end()) < 10) {
248  if (rank == 0) throw std::runtime_error("ERROR: lack of information on " + *beg + " in " + infilename + ".\n");
249  else sleep(2);
250  } else if (type.compare("AsyGausObservable") == 0) {
251  if (rank == 0) throw std::runtime_error("ERROR: AsyGausObservable cannot be specified with a input histogram, use Observable keyword instead.\n");
252  else sleep(2);
253  } else {
254  filename = filepath + *(++beg);
255  histoname = *(++beg);
256  setLikelihoodFromHisto(filename, histoname);
257  if (rank == 0) std::cout << "added input histogram " << filename << "/" << histoname << std::endl;
258  }
259  } else if (distr.compare("weight") == 0) {
260  if (std::distance(tok->begin(), tok->end()) < 11) {
261  if(rank == 0) throw std::runtime_error("ERROR: lack of information on " + *beg + " in " + infilename + ".\n");
262  else sleep (2);
263  }
264  ++beg;
265  ave = atof((*beg).c_str());
266  ++beg;
267  if (type.compare("AsyGausObservable") == 0) errgl = atof((*beg).c_str());
268  else errg = atof((*beg).c_str());
269  ++beg;
270  if (type.compare("AsyGausObservable") == 0) errgr = atof((*beg).c_str());
271  else errf = atof((*beg).c_str());
272  if (type.compare("AsyGausObservable") != 0 && errf == 0. && errg == 0. && !hasInverseCovariance) {
273  if (rank == 0) throw std::runtime_error("ERROR: The Gaussian and flat error in weight for " + name + " cannot both be 0. in the " + infilename + " .\n");
274  else sleep(2);
275  }
276  if (type.compare("AsyGausObservable") == 0 && errgl == 0. && errgr == 0. && !hasInverseCovariance) {
277  if (rank == 0) throw std::runtime_error("ERROR: The left Gaussian and right Gaussian error in weight for " + name + " cannot both be 0. in the " + infilename + " .\n");
278  else sleep(2);
279  }
280  } else if (distr.compare("noweight") == 0 || distr.compare("writeChain") == 0) {
281  if (!tMCMC) writeChain = (distr.compare("writeChain") == 0);
282  else {
283  if (rank == 0) throw std::runtime_error("ERROR: writeChains cannot be requested since " + name + " is set to MCMC in " + infilename + " .\n");
284  else sleep(2);
285  }
286  if (obsType.compare("BinnedObservable") == 0 || obsType.compare("FunctionObservable") == 0) {
287  ++beg;
288  ++beg;
289  ++beg;
290  }
291  } else {
292  if (rank == 0) throw std::runtime_error("ERROR: wrong distribution flag in " + name + " in file " + infilename + ": can be only weight, noweight or writeChain.\n");
293  else sleep(2);
294  }
295  ++beg;
296  if (obsType.compare("BinnedObservable") == 0) {
297  bin_min = atof((*beg).c_str());
298  ++beg;
299  bin_max = atof((*beg).c_str());
300  ++beg;
301  } else if (obsType.compare("FunctionObservable") == 0) {
302  bin_min = atof((*beg).c_str());
303  ++beg;
304  }
305  if (beg != tok->end() && rank == 0) std::cout << "WARNING: unread information in observable " << name << std::endl;
306  }
307  }
308  return beg;
309 }
Observable::tho
ThObservable * tho
A pointer of to the object of the ThObservables class.
Definition: Observable.h:478
operator<<
std::ostream & operator<<(std::ostream &output, const Observable &o)
Definition: Observable.cpp:102
Observable::computeTheoryValue
double computeTheoryValue()
A method to access the computed theory value of the observable.
Definition: Observable.cpp:130
ThObservable::computeThValue
virtual double computeThValue()=0
A member to be overloaded by the respective theory observable. class that calculates the value of the...
StandardModel::getIterationNo
int getIterationNo() const
Definition: StandardModel.h:582
Observable::histoname
std::string histoname
The name of the histogram for the observable.
Definition: Observable.h:484
Observable::setMin
void setMin(double min)
A set method to fix the minimum value for the observable.
Definition: Observable.h:314
Observable::max
double max
The maximum valus of the observable.
Definition: Observable.h:491
Observable::tMCMC
bool tMCMC
The flag to include or exclude the observable from the MCMC run.
Definition: Observable.h:492
Observable::computeWeight
virtual double computeWeight()
A method to compute the weight associated with the observable.
Definition: Observable.h:113
Observable::obsType
std::string obsType
Type of the Observable. 0: Observable, 1: HiggsObservable, 2: BinnedObservable, 3: FunctionObservable...
Definition: Observable.h:495
Observable::setLikelihoodFromHisto
virtual void setLikelihoodFromHisto(std::string filename, std::string histoname)
A set method to set the likelihood from which the experimental likelihood of the observable will be r...
Definition: Observable.cpp:111
Observable::ParseObservable
boost::tokenizer< boost::char_separator< char > >::iterator & ParseObservable(std::string &type, boost::tokenizer< boost::char_separator< char > > *tok, boost::tokenizer< boost::char_separator< char > >::iterator &beg, std::string &filepath, std::string &infilename, int rank)
The parser for Observables.
Definition: Observable.cpp:202
ThObservable::getModel
const StandardModel & getModel()
A get method to get the model.
Definition: ThObservable.h:100
StandardModel.h
Observable::min
double min
The minimum value of the observable.
Definition: Observable.h:490
Observable::ave
double ave
The average value of the observable.
Definition: Observable.h:485
gslpp::log
complex log(const complex &z)
Definition: gslpp_complex.cpp:342
Observable::errgr
double errgr
The upper gaussian error of the observable.
Definition: Observable.h:489
Observable::thValue
double thValue
The theory value of the first observable.
Definition: Observable.h:499
Observable::writeChain
bool writeChain
The flag to write the chain for the observable from the MCMC run.
Definition: Observable.h:493
Observable::getObsType
std::string getObsType() const
A get method to get the Observable type.
Definition: Observable.h:404
Observable::LogSplitGaussian
double LogSplitGaussian(double x, double ave, double errl, double errr)
Definition: Observable.cpp:141
gslpp::sqrt
complex sqrt(const complex &z)
Definition: gslpp_complex.cpp:385
Observable::errg
double errg
The gaussian error of the observable.
Definition: Observable.h:486
Observable::errf
double errf
The flat error of the observable.
Definition: Observable.h:487
Observable.h
Observable::thname
std::string thname
The name for the observable as fixed in the ThObservable class.
Definition: Observable.h:480
Observable
A class for observables.
Definition: Observable.h:28
Observable::~Observable
virtual ~Observable()
The default destructor.
Definition: Observable.cpp:100
Observable::hasInverseCovariance
bool hasInverseCovariance
Definition: Observable.h:500
Observable::name
std::string name
A name for the observable.
Definition: Observable.h:479
Observable::label
std::string label
A label for the observable.
Definition: Observable.h:481
Observable::iterationNo
int iterationNo
A counter for the interation that helps with the observable caching.
Definition: Observable.h:498
Observable::Observable
Observable()
The default constructor.
Definition: Observable.cpp:75
ThObservable
A class for a model prediction of an observable.
Definition: ThObservable.h:25
Observable::errgl
double errgl
The lower gaussian error of the observable.
Definition: Observable.h:488
Observable::filename
std::string filename
The name of the file containing the experimental likelihood for the observable.
Definition: Observable.h:483
Observable::bin_max
double bin_max
The maximum valus of the observable bin.
Definition: Observable.h:497
Observable::distr
std::string distr
The name of the distribution of the the observable.
Definition: Observable.h:482
Observable::setMax
void setMax(double max)
A set method to fix the maximum value for the observable.
Definition: Observable.h:296
Observable::LogGaussian
double LogGaussian(double x, double ave, double sigma)
Definition: Observable.cpp:147
Observable::inhisto
TH1D * inhisto
1D Histogram containing the experimental likelihood for the observable
Definition: Observable.h:494
Observable::bin_min
double bin_min
The minimum value of the observable bin.
Definition: Observable.h:496