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