a Code for the Combination of Indirect and Direct Constraints on High Energy Physics Models Logo
Observable2D.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 "Observable2D.h"
9 #include "StandardModel.h"
10 #include <TFile.h>
11 #include <fstream>
12 #include <TROOT.h>
13 #include <limits>
14 
15 Observable2D::Observable2D(const std::string name_i,
16  const std::string thname_i,
17  const std::string thname2_i,
18  const std::string label_i,
19  const std::string label2_i,
20  const bool tMCMC_i,
21  const double min_i,
22  const double max_i,
23  const double min2_i,
24  const double max2_i,
25  ThObservable * tho_i,
26  ThObservable * tho2_i)
27 : Observable(name_i, thname_i, label_i, tMCMC_i, min_i, max_i, tho_i),
28  bin_min(2, 0.),
29  bin_max(2, 0.)
30 {
31  thname2 = thname2_i;
32  label2 = label2_i;
33  min2 = min2_i;
34  max2 = max2_i;
35  tho2 = tho2_i;
36  obsType2 = "";
37  filepath = "";
38  iterationNo2 = std::numeric_limits<int>::max();
39  IsEOF = false;
40 }
41 
42 Observable2D::Observable2D(const Observable& o1d)
43 : Observable(o1d),
44  bin_min(2, 0.),
45  bin_max(2, 0.)
46 {
47  thname2 = "";
48  label2 = "";
49  min2 = 0.;
50  max2 = 0.;
51  tho2 = NULL;
52  obsType2 = "";
53  filepath = "";
54  iterationNo2 = std::numeric_limits<int>::max();
55  IsEOF = false;
56 }
57 
58 Observable2D::Observable2D()
59 : Observable(),
60  bin_min(2, 0.),
61  bin_max(2, 0.)
62 {
63  thname2 = "";
64  label2 = "";
65  min2 = 0.;
66  max2 = 0.;
67  tho2 = NULL;
68  obsType2 = "";
69  filepath = "";
70  iterationNo2 = std::numeric_limits<int>::max();
71  IsEOF = false;
72 }
73 
74 Observable2D::Observable2D(const Observable2D& orig)
75 : Observable(orig.name, orig.thname, orig.label, orig.tMCMC, orig.min, orig.max, orig.tho),
76  bin_min(orig.bin_min),
77  bin_max(orig.bin_max)
78 {
79 
80  distr = orig.distr;
81  filename = orig.filename;
82  histoname = orig.histoname;
83  ave = orig.ave;
84  errg = orig.errg;
85  errg2 = orig.errg2;
86  errf = orig.errf;
87  errf2 = orig.errf2;
88  errgl = orig.errgl;
89  errgl2 = orig.errgl2;
90  errgr = orig.errgr;
91  errgr2 = orig.errgr2;
92 
93  thname2 = orig.thname2;
94  label2 = orig.label2;
95  min2 = orig.min2;
96  max2 = orig.max2;
97  tho2 = orig.tho2;
98  filepath = orig.filepath;
99  iterationNo2 = orig.iterationNo2;
100  IsEOF = orig.IsEOF;
101  inhisto2d = orig.inhisto2d;
102 }
103 
104 Observable2D::~Observable2D()
105 {}
106 
107 double Observable2D::computeTheoryValue2()
108 {
109  if (tho2->getModel().getIterationNo() == iterationNo2) {
110  return thValue2;
111  } else {
112  iterationNo2 = tho2->getModel().getIterationNo();
113  thValue2 = tho2->computeThValue();
114  return thValue2;
115  }
116 }
117 
118 void Observable2D::setLikelihoodFromHisto(std::string filename, std::string histoname)
119 {
120  this->filename = filename;
121  this->histoname = histoname;
122  TFile *lik2 = new TFile((filename + ".root").c_str(), "read");
123  TH2D *htmp2 = (TH2D*) (lik2->Get(histoname.c_str()));
124  if (htmp2 == NULL)
125  throw std::runtime_error("ERROR: nonexistent histogram called "
126  + histoname
127  + " in " + filename + ".root");
128  inhisto2d = (TH2D *) htmp2->Clone((filename + "/" + histoname).c_str());
129  inhisto2d->SetDirectory(gROOT);
130  std::cout << "added 2D input histogram " << name << std::endl;
131  setMin(inhisto2d->GetXaxis()->GetXmin());
132  setMax(inhisto2d->GetXaxis()->GetXmax());
133  setMin2(inhisto2d->GetYaxis()->GetXmin());
134  setMax2(inhisto2d->GetYaxis()->GetXmax());
135  lik2->Close();
136  delete lik2;
137 }
138 
139 double Observable2D::computeWeight(double th1, double th2)
140 {
141  double logprob;
142  if (distr.compare("file") == 0) {
143  int i = inhisto2d->FindBin(th1, th2);
144  if (inhisto2d->IsBinOverflow(i) || inhisto2d->IsBinUnderflow(i))
145  logprob = log(0.);
146  else
147  logprob = log(inhisto2d->GetBinContent(i));
148  //logprob = log(h->GetBinContent(h->GetXaxis()->FindBin(th1),h->GetYaxis()->FindBin(th2)));
149  } else if (distr.compare("weight") == 0) {
150  double wt2;
151  if (obsType2.compare("AsyGausObservable") == 0) wt2 = Observable::LogSplitGaussian(th2, ave2, errgl, errgr);
152  else wt2 = Observable::computeWeight(th2, ave2, errg2, errf2);
153  logprob = Observable::computeWeight(th1) + wt2;
154  } else
155  throw std::runtime_error("ERROR: 2D MonteCarloEngine::Weight() called without file for "
156  + name);
157  return (logprob);
158 }
159 
160 int Observable2D::ParseObservable2D(std::string& type,
161  boost::tokenizer<boost::char_separator<char> >* tok,
162  boost::tokenizer<boost::char_separator<char> >::iterator& beg,
163  std::string& infilename,
164  std::ifstream& ifile,
165  int lineNo,
166  int rank)
167 {
168  if (infilename.find("\\/") == std::string::npos) filepath = infilename.substr(0, infilename.find_last_of("\\/") + 1);
169  if (std::distance(tok->begin(), tok->end()) < 12) {
170  setName(*beg);
171  ++beg;
172  if (std::distance(tok->begin(), tok->end()) < 4) {
173  if(rank == 0) throw std::runtime_error("ERROR: lack of information on " + name + " in " + infilename + " at line number" + boost::lexical_cast<std::string>(lineNo));
174  else sleep (2);
175  }
176  std::string toMCMC = *beg;
177  if (toMCMC.compare("MCMC") == 0)
178  setTMCMC(true);
179  else if (toMCMC.compare("noMCMC") == 0)
180  setTMCMC(false);
181  else {
182  if (rank == 0) throw std::runtime_error("ERROR: wrong MCMC flag in Observable2D" + name + " at line number:" + boost::lexical_cast<std::string>(lineNo) + " of file " + infilename + ".\n");
183  else sleep(2);
184  }
185 
186  ++beg;
187  setDistr(*beg);
188  if (distr.compare("file") == 0) {
189  if (std::distance(tok->begin(), tok->end()) < 6) {
190  if(rank == 0) throw std::runtime_error("ERROR: lack of information on "+ *beg + " in " + infilename);
191  else sleep (2);
192  }
193  setFilename(filepath + *(++beg));
194  setHistoname(*(++beg));
195  }
196 
197  std::vector<double> min(2, 0.);
198  std::vector<double> max(2, 0.);
199  std::vector<double> ave(2, 0.);
200  std::vector<double> errg(2, 0.);
201  std::vector<double> errf(2, 0.);
202  std::vector<double> errgl(2, 0.);
203  std::vector<double> errgr(2, 0.);
204  std::vector<std::string> thname(2, "");
205  std::vector<std::string> label(2, "");
206  std::vector<std::string> type2D(2, "");
207  std::string line;
208  size_t pos = 0;
209  boost::char_separator<char> sep(" \t");
210  for (int i = 0; i < 2; i++) {
211  IsEOF = getline(ifile, line).eof();
212  if (line.empty() || line.at(0) == '#') {
213  if (rank == 0) throw std::runtime_error("ERROR: no comments or empty lines in Observable2D please! In file " + infilename + " at line number:" + boost::lexical_cast<std::string>(lineNo) + ".\n");
214  else sleep(2);
215  }
216  lineNo++;
217  boost::tokenizer<boost::char_separator<char> > mytok(line, sep);
218  beg = mytok.begin();
219  type2D[i] = *beg;
220  if (type2D[i].compare("Observable") != 0 && type2D[i].compare("BinnedObservable") != 0 && type2D[i].compare("FunctionObservable") != 0 && type2D[i].compare("AsyGausObservable") != 0) {
221  if (rank == 0) throw std::runtime_error("ERROR: in line no." + boost::lexical_cast<std::string>(lineNo) + " of file " + infilename + ", expecting an Observable or BinnedObservable or FunctionObservable type here...\n");
222  else sleep(2);
223  }
224  ++beg;
225  thname[i] = *beg;
226  ++beg;
227  label[i] = *beg;
228  while ((pos = label[i].find("~", pos)) != std::string::npos)
229  label[i].replace(pos++, 1, " ");
230  ++beg;
231  min[i] = atof((*beg).c_str());
232  ++beg;
233  max[i] = atof((*beg).c_str());
234  if (distr.compare("weight") == 0) {
235  ++beg;
236  ave[i] = atof((*beg).c_str());
237  ++beg;
238  if (type.compare("AsyGausObservable") == 0) errgl[i] = atof((*beg).c_str());
239  else errg[i] = atof((*beg).c_str());
240  ++beg;
241  if (type.compare("AsyGausObservable") == 0) errgr[i] = atof((*beg).c_str());
242  else errf[i] = atof((*beg).c_str());
243  if (type.compare("AsyGausObservable") != 0 && errg[i] == 0. && errf[i] == 0.) {
244  if (rank == 0) throw std::runtime_error("ERROR: The Gaussian and flat error in weight for " + name + " cannot both be 0. in the " + infilename + " file, line number:" + boost::lexical_cast<std::string>(lineNo) + ".\n");
245  else sleep(2);
246  }
247  if (type.compare("AsyGausObservable") == 0 && errgl[i] == 0. && errgr[i] == 0.) {
248  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 + " file, line number:" + boost::lexical_cast<std::string>(lineNo) + ".\n");
249  else sleep(2);
250  }
251  } else if (distr.compare("noweight") == 0 || distr.compare("file") == 0) {
252  if (type2D[i].compare("BinnedObservable") == 0 || type2D[i].compare("FunctionObservable") == 0) {
253  ++beg;
254  ++beg;
255  ++beg;
256  }
257  } else {
258  if (rank == 0) throw std::runtime_error("ERROR: wrong distribution flag in " + name + " in file " + infilename + ".\n");
259  else sleep(2);
260  }
261  if (type2D[i].compare("BinnedObservable") == 0) {
262  ++beg;
263  bin_min[i] = atof((*beg).c_str());
264  ++beg;
265  bin_max[i] = atof((*beg).c_str());
266  } else if (type2D[i].compare("FunctionObservable") == 0) {
267  ++beg;
268  bin_min[i] = atof((*beg).c_str());
269  ++beg;
270  }
271  }
272  setObsType(type2D[0]);
273  obsType2 = type2D[1];
274  setThname(thname[0]);
275  thname2 = thname[1];
276  setLabel(label[0]);
277  label2 = label[1];
278  setMin(min[0]);
279  min2 = min[1];
280  setMax(max[0]);
281  max2= max[1];
282  setAve(ave[0]);
283  ave2 = ave[1];
284  setErrg(errg[0]);
285  errg2 = errg[1];
286  setErrf(errf[0]);
287  errf2 = errf[1];
288  setErrgl(errgl[0]);
289  errgl2 = errgl[1];
290  setErrgr(errgr[0]);
291  errgr2 = errgr[1];
292  if (distr.compare("file") == 0) {
293  setLikelihoodFromHisto(filename, histoname);
294  if (rank == 0) std::cout << "added input histogram " << filename << "/" << histoname << std::endl;
295  }
296  return lineNo;
297  } else {
298  beg = ParseObservable(type, tok, beg, filepath, filename, rank);
299  ++beg;
300  std::string distr = *beg;
301  if (distr.compare("file") == 0) {
302  if (std::distance(tok->begin(), tok->end()) < 14) {
303  if (rank == 0) throw std::runtime_error("ERROR: lack of information on " + *beg + " in " + infilename + ".\n");
304  else sleep(2);
305  }
306  setFilename(filepath + *(++beg));
307  setHistoname(*(++beg));
308  setLikelihoodFromHisto(filename, histoname);
309  if (rank == 0) std::cout << "added input histogram " << filename << "/" << histoname << std::endl;
310  } else if (distr.compare("noweight") == 0) {
311  } else {
312  if (rank == 0) throw std::runtime_error("ERROR: wrong distribution flag in " + name);
313  else sleep(2);
314  }
315  setDistr(distr);
316  ++beg;
317  thname2 = *beg;
318  ++beg;
319  std::string label = *beg;
320  size_t pos = 0;
321  while ((pos = label.find("~", pos)) != std::string::npos)
322  label.replace(pos, 1, " ");
323  label2 = label;
324  ++beg;
325  min2 = atof((*beg).c_str());
326  ++beg;
327  max2 = atof((*beg).c_str());
328  ++beg;
329  if (beg != tok->end())
330  if (rank == 0) std::cout << "WARNING: unread information in observable2D " << name << std::endl;
331  return lineNo;
332  }
333 }
Observable2D::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: Observable2D.cpp:118
ThObservable::computeThValue
virtual double computeThValue()=0
A member to be overloaded by the respective theory observable. class that calculates the value of the...
Observable2D::ParseObservable2D
int ParseObservable2D(std::string &type, boost::tokenizer< boost::char_separator< char > > *tok, boost::tokenizer< boost::char_separator< char > >::iterator &beg, std::string &infilename, std::ifstream &ifile, int lineNo, int rank)
Definition: Observable2D.cpp:160
StandardModel::getIterationNo
int getIterationNo() const
Definition: StandardModel.h:582
Observable2D::computeWeight
virtual double computeWeight()
A method to compute the weight associated with the observable.
Definition: Observable2D.h:212
Observable2D::setMax2
void setMax2(double max2)
A set method to fix the maximum value for the second observable.
Definition: Observable2D.h:128
Observable::histoname
std::string histoname
The name of the histogram for the observable.
Definition: Observable.h:484
Observable2D::ave2
double ave2
The average value of the second observable.
Definition: Observable2D.h:295
Observable::setMin
void setMin(double min)
A set method to fix the minimum value for the observable.
Definition: Observable.h:314
Observable2D::tho2
ThObservable * tho2
A pointer to an object of the ThObservable class.
Definition: Observable2D.h:301
Observable::max
double max
The maximum valus of the observable.
Definition: Observable.h:491
Observable::setErrgl
void setErrgl(double errgl)
A set method to fix the left Gaussian error of the observable.
Definition: Observable.h:194
Observable::computeWeight
virtual double computeWeight()
A method to compute the weight associated with the observable.
Definition: Observable.h:113
Observable2D::thname2
std::string thname2
The name for the second observable as fixed in the ThObservable() class.
Definition: Observable2D.h:291
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
Observable::setErrf
void setErrf(double errf)
A set method to fix the flat error of the observable.
Definition: Observable.h:167
Observable::setErrg
void setErrg(double errg)
A set method to fix the gaussian error of the observable.
Definition: Observable.h:185
ThObservable::getModel
const StandardModel & getModel()
A get method to get the model.
Definition: ThObservable.h:100
StandardModel.h
Observable2D::bin_max
std::vector< double > bin_max
The maximum value of the bin.
Definition: Observable2D.h:304
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
Observable2D::label2
std::string label2
A label for the second observable.
Definition: Observable2D.h:292
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
Observable2D::errg2
double errg2
The Gaussian error of the second observable.
Definition: Observable2D.h:296
Observable2D::min2
double min2
The minimum value of the second observable.
Definition: Observable2D.h:293
Observable::setThname
void setThname(std::string thname)
A set method to fix the name of the observable as listed in ThFactory class.
Definition: Observable.h:377
Observable2D::errgr2
double errgr2
The right Gaussian error of the second observable.
Definition: Observable2D.h:299
Observable2D::IsEOF
bool IsEOF
A bolean that is true if the end of file is reached.
Definition: Observable2D.h:308
Observable2D::bin_min
std::vector< double > bin_min
The minimum value of the bin.
Definition: Observable2D.h:303
Observable2D::filepath
std::string filepath
The path to the file being parsed.
Definition: Observable2D.h:305
Observable::setName
void setName(std::string name)
A set method to fix the name for the observable.
Definition: Observable.h:332
Observable::LogSplitGaussian
double LogSplitGaussian(double x, double ave, double errl, double errr)
Definition: Observable.cpp:141
Observable::setAve
void setAve(double ave)
A set method to fix the average value of the observable.
Definition: Observable.h:131
Observable2D::errgl2
double errgl2
The left Gaussian error of the second observable.
Definition: Observable2D.h:298
Observable::errg
double errg
The gaussian error of the observable.
Definition: Observable.h:486
Observable2D::iterationNo2
int iterationNo2
Counts the iteration to help with caching.
Definition: Observable2D.h:306
Observable::errf
double errf
The flat error of the observable.
Definition: Observable.h:487
Observable::thname
std::string thname
The name for the observable as fixed in the ThObservable class.
Definition: Observable.h:480
Observable::setDistr
void setDistr(std::string distr)
A set method to fix the name of the distribution of the observable.
Definition: Observable.h:149
Observable
A class for observables.
Definition: Observable.h:28
Observable::setErrgr
void setErrgr(double errgr)
A set method to fix the right Gaussian error of the observable.
Definition: Observable.h:203
Observable2D::inhisto2d
TH2D * inhisto2d
2D Histogram containing the experimental likelihood for the observable.
Definition: Observable2D.h:302
Observable2D::errf2
double errf2
the flat error of the second observable.
Definition: Observable2D.h:297
Observable2D
A class for analyzing observables pairwise.
Definition: Observable2D.h:24
Observable::name
std::string name
A name for the observable.
Definition: Observable.h:479
Observable2D::computeTheoryValue2
double computeTheoryValue2()
A method to access the computed theory value of the second observable.
Definition: Observable2D.cpp:107
Observable2D::obsType2
std::string obsType2
Type of the second Observable. 0: Observable, 1: HiggsObservable, 2: BinnedObservable,...
Definition: Observable2D.h:300
Observable2D.h
Observable::label
std::string label
A label for the observable.
Definition: Observable.h:481
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
Observable2D::thValue2
double thValue2
The theory value of the second observable.
Definition: Observable2D.h:307
Observable::setFilename
void setFilename(std::string filename_i)
Definition: Observable.h:217
Observable2D::Observable2D
Observable2D()
The default constructor.
Definition: Observable2D.cpp:58
Observable::filename
std::string filename
The name of the file containing the experimental likelihood for the observable.
Definition: Observable.h:483
Observable::setLabel
void setLabel(std::string label)
A set method to fix the label for the observable.
Definition: Observable.h:278
Observable::setObsType
void setObsType(std::string &obsType_s)
A set method to set the Observable type.
Definition: Observable.h:395
Observable2D::~Observable2D
virtual ~Observable2D()
The default destructor.
Definition: Observable2D.cpp:104
Observable::distr
std::string distr
The name of the distribution of the the observable.
Definition: Observable.h:482
Observable::setTMCMC
void setTMCMC(bool tMCMC)
A set method to fix the observable's inclusion in the MCMC listing.
Definition: Observable.h:359
Observable2D::setMin2
void setMin2(double min2)
A set method to fix the minimum value for the second observable.
Definition: Observable2D.h:146
Observable::setMax
void setMax(double max)
A set method to fix the maximum value for the observable.
Definition: Observable.h:296
Observable2D::max2
double max2
The maximum valus of the second observable.
Definition: Observable2D.h:294
Observable::setHistoname
void setHistoname(std::string histoname_i)
A set method to set the name of the histogram containing the likelihood.
Definition: Observable.h:260