Logo ROOT   6.14/05
Reference Guide
HistoToWorkspaceFactory.cxx
Go to the documentation of this file.
1 // @(#)root/roostats:$Id: cranmer $
2 // Author: Kyle Cranmer, Akira Shibata
3 /*************************************************************************
4  * Copyright (C) 1995-2008, Rene Brun and Fons Rademakers. *
5  * All rights reserved. *
6  * *
7  * For the licensing terms see $ROOTSYS/LICENSE. *
8  * For the list of contributors see $ROOTSYS/README/CREDITS. *
9  *************************************************************************/
10 
11 ////////////////////////////////////////////////////////////////////////////////
12 
13 /*
14 BEGIN_HTML
15 <p>
16 </p>
17 END_HTML
18 */
19 //
20 
21 
22 #ifndef __CINT__
23 #include "RooGlobalFunc.h"
24 #endif
25 
26 // Roofit/Roostat include
27 #include "RooDataSet.h"
28 #include "RooRealVar.h"
29 #include "RooConstVar.h"
30 #include "RooAddition.h"
31 #include "RooProduct.h"
32 #include "RooProdPdf.h"
33 #include "RooAddPdf.h"
34 #include "RooGaussian.h"
35 #include "RooExponential.h"
36 #include "RooRandom.h"
37 #include "RooCategory.h"
38 #include "RooSimultaneous.h"
39 #include "RooMultiVarGaussian.h"
40 #include "RooNumIntConfig.h"
41 #include "RooMinuit.h"
42 #include "RooNLLVar.h"
43 #include "RooProfileLL.h"
44 #include "RooFitResult.h"
45 #include "RooDataHist.h"
46 #include "RooHistPdf.h"
47 #include "RooProduct.h"
48 #include "RooWorkspace.h"
49 #include "RooCustomizer.h"
50 #include "RooPlot.h"
51 #include "RooMsgService.h"
52 #include "RooStats/RooStatsUtils.h"
53 #include "RooStats/ModelConfig.h"
54 
55 #include "TH2F.h"
56 #include "TH3F.h"
57 #include "TFile.h"
58 #include "TCanvas.h"
59 #include "TH1.h"
60 #include "TLine.h"
61 #include "TTree.h"
62 #include "TMarker.h"
63 #include "TStopwatch.h"
64 #include "TROOT.h"
65 #include "TStyle.h"
66 #include "TVectorD.h"
67 #include "TMatrixDSym.h"
68 
69 // specific to this package
70 //#include "RooStats/HistFactory/Helper.h"
73 #include "Helper.h"
74 
75 #define VERBOSE
76 
77 #define alpha_Low "-5"
78 #define alpha_High "5"
79 #define NoHistConst_Low "0"
80 #define NoHistConst_High "2000"
81 
82 // use this order for safety on library loading
83 using namespace RooFit ;
84 using namespace RooStats ;
85 using namespace std ;
86 //using namespace RooMsgService ;
87 
89 
90 namespace RooStats{
91 namespace HistFactory{
92 
93  HistoToWorkspaceFactory::HistoToWorkspaceFactory() :
94  fNomLumi(0),
95  fLumiError(0),
96  fLowBin(0),
97  fHighBin(0),
98  fOut_f(0),
99  pFile(0)
100  {
101  }
102 
104  fclose(pFile);
105  }
106 
107  HistoToWorkspaceFactory::HistoToWorkspaceFactory(string filePrefix, string row, vector<string> syst, double nomL, double lumiE, int low, int high, TFile* f):
108  fFileNamePrefix(filePrefix),
109  fRowTitle(row),
110  fSystToFix(syst),
111  fNomLumi(nomL),
112  fLumiError(lumiE),
113  fLowBin(low),
114  fHighBin(high),
115  fOut_f(f) {
116 
117  // fResultsPrefixStr<<"results" << "_" << fNomLumi<< "_" << fLumiError<< "_" << fLowBin<< "_" << fHighBin;
118  fResultsPrefixStr<< "_" << fRowTitle;
119  while(fRowTitle.find("\\ ")!=string::npos){
120  int pos=fRowTitle.find("\\ ");
121  fRowTitle.replace(pos, 1, "");
122  }
123  pFile = fopen ((filePrefix+"_results.table").c_str(),"a");
124  //RooMsgService::instance().setGlobalKillBelow(RooFit::ERROR) ;
125 
126  }
127 
129 
130  stringstream ss;
131  ss << prefix << "_" << fNomLumi<< "_" << fLumiError<< "_" << fLowBin<< "_" << fHighBin<< "_"<<fRowTitle;
132 
133  return ss.str();
134  }
135 
136  void HistoToWorkspaceFactory::ProcessExpectedHisto(TH1* hist,RooWorkspace* proto, string prefix, string productPrefix, string systTerm, double low, double high, int lowBin, int highBin){
137  if(hist)
138  cout << "processing hist " << hist->GetName() << endl;
139  else
140  cout << "hist is empty" << endl;
141  RooArgSet argset(prefix.c_str());
142  string highStr = "inf";
143  for(Int_t i=lowBin; i<highBin; ++i){
144  std::stringstream str;
145  std::stringstream range;
146  str<<"_"<<i;
147  if(hist)
148  range<<"["<<hist->GetBinContent(i+1) << "," << low << "," << highStr << "]";
149  else
150  range<<"["<< low << "," << high << "]";
151  cout << "for bin N"+str.str() << " var " << prefix+str.str()+" with range " << range.str() << endl;
152  RooRealVar* var = (RooRealVar*) proto->factory((prefix+str.str()+range.str()).c_str());
153 
154  // now create the product of the overall efficiency times the sigma(params) for this estimate
155  if(! (productPrefix.empty() || systTerm.empty()) )
156  proto->factory(("prod:"+productPrefix+str.str()+"("+prefix+str.str()+","+systTerm+")").c_str() );
157 
158  var->setConstant();
159  argset.add(* var );
160  }
161  proto->defineSet(prefix.c_str(),argset);
162  // proto->Print();
163  }
164 
165  void HistoToWorkspaceFactory::AddMultiVarGaussConstraint(RooWorkspace* proto, string prefix,int lowBin, int highBin, vector<string>& likelihoodTermNames){
166  // these are the nominal predictions: eg. the mean of some space of variations
167  // later fill these in a loop over histogram bins
168  TVectorD mean(highBin-lowBin);
169  cout << "a" << endl;
170  for(Int_t i=lowBin; i<highBin; ++i){
171  std::stringstream str;
172  str<<"_"<<i;
173  RooRealVar* temp = proto->var((prefix+str.str()).c_str());
174  mean(i) = temp->getVal();
175  }
176 
177  TMatrixDSym Cov(highBin-lowBin);
178  for(int i=lowBin; i<highBin; ++i){
179  for(int j=0; j<highBin-lowBin; ++j){
180  if(i==j)
181  Cov(i,j) = sqrt(mean(i));
182  else
183  Cov(i,j) = 0;
184  }
185  }
186  // can't make MultiVarGaussian with factory yet, do it by hand
187  RooArgList floating( *(proto->set(prefix.c_str() ) ) );
188  RooMultiVarGaussian constraint((prefix+"Constraint").c_str(),"",
189  floating, mean, Cov);
190 
191  proto->import(constraint);
192 
193  likelihoodTermNames.push_back(constraint.GetName());
194 
195  }
196 
197 
198  void HistoToWorkspaceFactory::LinInterpWithConstraint(RooWorkspace* proto, TH1* nominal, vector<TH1*> lowHist, vector<TH1*> highHist,
199  vector<string> sourceName, string prefix, string productPrefix, string systTerm,
200  int lowBin, int highBin, vector<string>& likelihoodTermNames){
201  // these are the nominal predictions: eg. the mean of some space of variations
202  // later fill these in a loop over histogram bins
203 
204  // make list of abstract parameters that interpolate in space of variations
205  RooArgList params( ("alpha_Hist") );
206  // range is set using defined macro (see top of the page)
207  string range=string("[")+alpha_Low+","+alpha_High+"]";
208  for(unsigned int j=0; j<lowHist.size(); ++j){
209  std::stringstream str;
210  str<<"_"<<j;
211 
212  RooRealVar* temp = (RooRealVar*) proto->var(("alpha_"+sourceName.at(j)).c_str());
213  if(!temp){
214  temp = (RooRealVar*) proto->factory(("alpha_"+sourceName.at(j)+range).c_str());
215 
216  // now add a constraint term for these parameters
217  string command=("Gaussian::alpha_"+sourceName.at(j)+"Constraint(alpha_"+sourceName.at(j)+",nom_"+sourceName.at(j)+"[0.,-10,10],1.)");
218  cout << command << endl;
219  likelihoodTermNames.push_back( proto->factory( command.c_str() )->GetName() );
220  proto->var(("nom_"+sourceName.at(j)).c_str())->setConstant();
221  const_cast<RooArgSet*>(proto->set("globalObservables"))->add(*proto->var(("nom_"+sourceName.at(j)).c_str()));
222 
223  }
224 
225  params.add(* temp );
226 
227  }
228 
229  // now make function that linearly interpolates expectation between variations
230  for(Int_t i=lowBin; i<highBin; ++i){
231  std::stringstream str;
232  str<<"_"<<i;
233 
234  // get low/high variations to interpolate between
235  vector<double> low, high;
236  for(unsigned int j=0; j<lowHist.size(); ++j){
237  low.push_back( lowHist.at(j)->GetBinContent(i+1) );
238  high.push_back( highHist.at(j)->GetBinContent(i+1) );
239  cout << "for "+prefix+" bin "+str.str()+" creating linear interp of nominal " << nominal->GetBinContent(i+1)
240  << " in parameter " << sourceName.at(j)
241  << " between " << low.back() << " - " << high.back()
242  << " about " << 100.*fabs(low.back() - high.back() )/nominal->GetBinContent(i+1) << " % error"
243  << endl;
244  }
245 
246  // this is sigma(params), a piece-wise linear interpolation
247  LinInterpVar interp( (prefix+str.str()).c_str(), "", params, nominal->GetBinContent(i+1), low, high);
248 
249  // cout << "check: " << interp.getVal() << endl;
250  proto->import(interp); // individual params have already been imported in first loop of this function
251 
252  // now create the product of the overall efficiency times the sigma(params) for this estimate
253  proto->factory(("prod:"+productPrefix+str.str()+"("+prefix+str.str()+","+systTerm+")").c_str() );
254 
255  }
256 
257  }
258 
259  string HistoToWorkspaceFactory::AddNormFactor(RooWorkspace * proto, string & channel, string & sigmaEpsilon, EstimateSummary & es, bool doRatio){
260  string overallNorm_times_sigmaEpsilon ;
261  string prodNames;
262  vector<EstimateSummary::NormFactor> norm=es.normFactor;
263  if(norm.size()){
264  for(vector<EstimateSummary::NormFactor>::iterator itr=norm.begin(); itr!=norm.end(); ++itr){
265  cout << "making normFactor: " << itr->name << endl;
266  // remove "doRatio" and name can be changed when ws gets imported to the combined model.
267  std::stringstream range;
268  range<<"["<<itr->val<<","<<itr->low<<","<<itr->high<<"]";
269  //RooRealVar* var = 0;
270 
271  string varname;
272  if(!prodNames.empty()) prodNames+=",";
273  if(doRatio) {
274  varname=itr->name+"_"+channel;
275  }
276  else {
277  varname=itr->name;
278  }
279  proto->factory((varname+range.str()).c_str());
280  if(itr->constant){
281  // proto->var(varname.c_str())->setConstant();
282  // cout <<"setting " << varname << " constant"<<endl;
283  cout <<"WARNING: Const attribute to <NormFactor> tag is deprecated, will ignore."<<
284  " Instead, add \n\t<ParamSetting Const=\"True\">"<<varname<<"</ParamSetting>\n"<<
285  " to your top-level XML's <Measurment> entry"<< endl;
286  }
287  prodNames+=varname;
288  }
289  overallNorm_times_sigmaEpsilon = es.name+"_"+channel+"_overallNorm_x_sigma_epsilon";
290  proto->factory(("prod::"+overallNorm_times_sigmaEpsilon+"("+prodNames+","+sigmaEpsilon+")").c_str());
291  }
292 
293  if(!overallNorm_times_sigmaEpsilon.empty())
294  return overallNorm_times_sigmaEpsilon;
295  else
296  return sigmaEpsilon;
297  }
298 
299 
300  void HistoToWorkspaceFactory::AddEfficiencyTerms(RooWorkspace* proto, string prefix, string interpName,
301  map<string,pair<double,double> > systMap,
302  vector<string>& likelihoodTermNames, vector<string>& totSystTermNames){
303  // add variables for all the relative overall uncertainties we expect
304 
305  // range is set using defined macro (see top of the page)
306  string range=string("[0,")+alpha_Low+","+alpha_High+"]";
307  //string range="[0,-1,1]";
308  totSystTermNames.push_back(prefix);
309  //bool first=true;
310  RooArgSet params(prefix.c_str());
311  vector<double> lowVec, highVec;
312  for(map<string,pair<double,double> >::iterator it=systMap.begin(); it!=systMap.end(); ++it){
313  // add efficiency term
314  RooRealVar* temp = (RooRealVar*) proto->var((prefix+ it->first).c_str());
315  if(!temp){
316  temp = (RooRealVar*) proto->factory((prefix+ it->first +range).c_str());
317 
318  string command=("Gaussian::"+prefix+it->first+"Constraint("+prefix+it->first+",nom_"+prefix+it->first+"[0.,-10,10],1.)");
319  cout << command << endl;
320  likelihoodTermNames.push_back( proto->factory( command.c_str() )->GetName() );
321  proto->var(("nom_"+prefix+it->first).c_str())->setConstant();
322  const_cast<RooArgSet*>(proto->set("globalObservables"))->add(*proto->var(("nom_"+prefix+it->first).c_str()));
323 
324  }
325  params.add(*temp);
326 
327  // add constraint in terms of bifrucated gauss with low/high as sigmas
328  std::stringstream lowhigh;
329  double low = it->second.first;
330  double high = it->second.second;
331  lowVec.push_back(low);
332  highVec.push_back(high);
333 
334  }
335  if(systMap.size()>0){
336  // this is epsilon(alpha_j), a piece-wise linear interpolation
337  LinInterpVar interp( (interpName).c_str(), "", params, 1., lowVec, highVec);
338  proto->import(interp); // params have already been imported in first loop of this function
339  } else{
340  // some strange behavior if params,lowVec,highVec are empty.
341  //cout << "WARNING: No OverallSyst terms" << endl;
342  RooConstVar interp( (interpName).c_str(), "", 1.);
343  proto->import(interp); // params have already been imported in first loop of this function
344  }
345 
346  }
347 
348 
349  void HistoToWorkspaceFactory::MakeTotalExpected(RooWorkspace* proto, string totName, string /**/, string /**/,
350  int lowBin, int highBin, vector<string>& syst_x_expectedPrefixNames,
351  vector<string>& normByNames){
352 
353  // for ith bin calculate totN_i = lumi * sum_j expected_j * syst_j
354 
355  for(Int_t i=lowBin; i<highBin; ++i){
356  std::stringstream str;
357  str<<"_"<<i;
358  string command="sum::"+totName+str.str()+"(";
359  //vector<string>::iterator it=syst_x_expectedPrefixNames.begin();
360  string prepend="";
361  for(unsigned int j=0; j<syst_x_expectedPrefixNames.size();++j){
362  command+=prepend+normByNames.at(j)+"*"+syst_x_expectedPrefixNames.at(j)+str.str();
363  prepend=",";
364  }
365  command+=")";
366  cout << "function to calculate total: " << command << endl;
367  proto->factory(command.c_str());
368  }
369  }
370 
371  void HistoToWorkspaceFactory::AddPoissonTerms(RooWorkspace* proto, string prefix, string obsPrefix, string expPrefix, int lowBin, int highBin,
372  vector<string>& likelihoodTermNames){
373  /////////////////////////////////
374  // Relate observables to expected for each bin
375  // later modify variable named expPrefix_i to be product of terms
376  RooArgSet Pois(prefix.c_str());
377  for(Int_t i=lowBin; i<highBin; ++i){
378  std::stringstream str;
379  str<<"_"<<i;
380  //string command("Poisson::"+prefix+str.str()+"("+obsPrefix+str.str()+","+expPrefix+str.str()+")");
381  string command("Poisson::"+prefix+str.str()+"("+obsPrefix+str.str()+","+expPrefix+str.str()+",1)");//for no rounding
382  RooAbsArg* temp = (proto->factory( command.c_str() ) );
383 
384  // output
385  cout << "Poisson Term " << command << endl;
386  ((RooAbsPdf*) temp)->setEvalErrorLoggingMode(RooAbsReal::PrintErrors);
387  //cout << temp << endl;
388 
389  likelihoodTermNames.push_back( temp->GetName() );
390  Pois.add(* temp );
391  }
392  proto->defineSet(prefix.c_str(),Pois); // add argset to workspace
393  }
394 
395  void HistoToWorkspaceFactory::SetObsToExpected(RooWorkspace* proto, string obsPrefix, string expPrefix, int lowBin, int highBin){
396  /////////////////////////////////
397  // set observed to expected
398  TTree* tree = new TTree();
399  Double_t* obsForTree = new Double_t[highBin-lowBin];
400  RooArgList obsList("obsList");
401 
402  for(Int_t i=lowBin; i<highBin; ++i){
403  std::stringstream str;
404  str<<"_"<<i;
405  RooRealVar* obs = (RooRealVar*) proto->var((obsPrefix+str.str()).c_str());
406  cout << "expected number of events called: " << expPrefix << endl;
407  RooAbsReal* exp = proto->function((expPrefix+str.str()).c_str());
408  if(obs && exp){
409 
410  //proto->Print();
411  obs->setVal( exp->getVal() );
412  cout << "setting obs"+str.str()+" to expected = " << exp->getVal() << " check: " << obs->getVal() << endl;
413 
414  // add entry to array and attach to tree
415  obsForTree[i] = exp->getVal();
416  tree->Branch((obsPrefix+str.str()).c_str(), obsForTree+i ,(obsPrefix+str.str()+"/D").c_str());
417  obsList.add(*obs);
418  }else{
419  cout << "problem retrieving obs or exp " << obsPrefix+str.str() << obs << " " << expPrefix+str.str() << exp << endl;
420  }
421  }
422  tree->Fill();
423  RooDataSet* data = new RooDataSet("expData","", tree, obsList); // one experiment
424 
425  proto->import(*data);
426  delete[] obsForTree;
427  obsForTree = nullptr;
428  }
429 
430  void HistoToWorkspaceFactory::Customize(RooWorkspace* proto, const char* pdfNameChar, map<string,string> renameMap) {
431  cout << "in customizations" << endl;
432  string pdfName(pdfNameChar);
433  map<string,string>::iterator it;
434  string edit="EDIT::customized("+pdfName+",";
435  string precede="";
436  for(it=renameMap.begin(); it!=renameMap.end(); ++it) {
437  cout << it->first + "=" + it->second << endl;
438  edit+=precede + it->first + "=" + it->second;
439  precede=",";
440  }
441  edit+=")";
442  cout << edit<< endl;
443  proto->factory( edit.c_str() );
444  }
445 
446  //////////////////////////////////////////////////////////////////////////////
447  /// cout << "in edit, gammamap.size = " << gammaSyst.size() << ", unimap.size = " << uniformSyst.size() << endl;
448 
449  void HistoToWorkspaceFactory::EditSyst(RooWorkspace* proto, const char* pdfNameChar, map<string,double> gammaSyst, map<string,double> uniformSyst,map<string,double> logNormSyst) {
450  string pdfName(pdfNameChar);
451 
452  ModelConfig * combined_config = (ModelConfig *) proto->obj("ModelConfig");
453  // const RooArgSet * constrainedParams=combined_config->GetNuisanceParameters();
454  // RooArgSet temp(*constrainedParams);
455  string edit="EDIT::newSimPdf("+pdfName+",";
456  string editList;
457  string lastPdf=pdfName;
458  string precede="";
459  unsigned int numReplacements = 0;
460  unsigned int nskipped = 0;
461  map<string,double>::iterator it;
462 
463  // add gamma terms and their constraints
464  for(it=gammaSyst.begin(); it!=gammaSyst.end(); ++it) {
465  //cout << "edit for " << it->first << "with rel uncert = " << it->second << endl;
466  if(! proto->var(("alpha_"+it->first).c_str())){
467  //cout << "systematic not there" << endl;
468  nskipped++;
469  continue;
470  }
471  numReplacements++;
472 
473  double relativeUncertainty = it->second;
474  double scale = 1/sqrt((1+1/pow(relativeUncertainty,2)));
475 
476  // this is the Gamma PDF and in a form that doesn't have roundoff problems like the Poisson does
477  proto->factory(Form("beta_%s[1,0,10]",it->first.c_str()));
478  proto->factory(Form("y_%s[%f]",it->first.c_str(),1./pow(relativeUncertainty,2))) ;
479  proto->factory(Form("theta_%s[%f]",it->first.c_str(),pow(relativeUncertainty,2))) ;
480  proto->factory(Form("Gamma::beta_%sConstraint(beta_%s,sum::k_%s(y_%s,one[1]),theta_%s,zero[0])",
481  it->first.c_str(),
482  it->first.c_str(),
483  it->first.c_str(),
484  it->first.c_str(),
485  it->first.c_str())) ;
486 
487  /*
488  // this has some problems because N in poisson is rounded to nearest integer
489  proto->factory(Form("Poisson::beta_%sConstraint(y_%s[%f],prod::taub_%s(taus_%s[%f],beta_%s[1,0,5]))",
490  it->first.c_str(),
491  it->first.c_str(),
492  1./pow(relativeUncertainty,2),
493  it->first.c_str(),
494  it->first.c_str(),
495  1./pow(relativeUncertainty,2),
496  it->first.c_str()
497  ) ) ;
498  */
499  // combined->factory(Form("expr::alphaOfBeta('(beta-1)/%f',beta)",scale));
500  // combined->factory(Form("expr::alphaOfBeta_%s('(beta_%s-1)/%f',beta_%s)",it->first.c_str(),it->first.c_str(),scale,it->first.c_str()));
501  proto->factory(Form("PolyVar::alphaOfBeta_%s(beta_%s,{%f,%f})",it->first.c_str(),it->first.c_str(),-1./scale,1./scale));
502 
503  // set beta const status to be same as alpha
504  if(proto->var(Form("alpha_%s",it->first.c_str()))->isConstant())
505  proto->var(Form("beta_%s",it->first.c_str()))->setConstant(true);
506  else
507  proto->var(Form("beta_%s",it->first.c_str()))->setConstant(false);
508  // set alpha const status to true
509  // proto->var(Form("alpha_%s",it->first.c_str()))->setConstant(true);
510 
511  // replace alphas with alphaOfBeta and replace constraints
512  //cout << "alpha_"+it->first+"Constraint=beta_" + it->first+ "Constraint" << endl;
513  editList+=precede + "alpha_"+it->first+"Constraint=beta_" + it->first+ "Constraint";
514  precede=",";
515  // cout << "alpha_"+it->first+"=alphaOfBeta_"+ it->first << endl;
516  editList+=precede + "alpha_"+it->first+"=alphaOfBeta_"+ it->first;
517 
518  /*
519  if( proto->pdf(("alpha_"+it->first+"Constraint").c_str()) && proto->var(("alpha_"+it->first).c_str()) )
520  cout << " checked they are there" << proto->pdf(("alpha_"+it->first+"Constraint").c_str()) << " " << proto->var(("alpha_"+it->first).c_str()) << endl;
521  else
522  cout << "NOT THERE" << endl;
523  */
524 
525  // EDIT seems to die if the list of edits is too long. So chunck them up.
526  if(numReplacements%10 == 0 && numReplacements+nskipped!=gammaSyst.size()){
527  edit="EDIT::"+lastPdf+"_("+lastPdf+","+editList+")";
528  lastPdf+="_"; // append an underscore for the edit
529  editList=""; // reset edit list
530  precede="";
531  cout << "Going to issue this edit command\n" << edit<< endl;
532  proto->factory( edit.c_str() );
533  RooAbsPdf* newOne = proto->pdf(lastPdf.c_str());
534  if(!newOne)
535  cout << "\n\n ---------------------\n WARNING: failed to make EDIT\n\n" << endl;
536 
537  }
538  }
539 
540  // add uniform terms and their constraints
541  for(it=uniformSyst.begin(); it!=uniformSyst.end(); ++it) {
542  cout << "edit for " << it->first << "with rel uncert = " << it->second << endl;
543  if(! proto->var(("alpha_"+it->first).c_str())){
544  cout << "systematic not there" << endl;
545  nskipped++;
546  continue;
547  }
548  numReplacements++;
549 
550  // this is the Uniform PDF
551  proto->factory(Form("beta_%s[1,0,10]",it->first.c_str()));
552  proto->factory(Form("Uniform::beta_%sConstraint(beta_%s)",it->first.c_str(),it->first.c_str()));
553  proto->factory(Form("PolyVar::alphaOfBeta_%s(beta_%s,{-1,1})",it->first.c_str(),it->first.c_str()));
554 
555  // set beta const status to be same as alpha
556  if(proto->var(Form("alpha_%s",it->first.c_str()))->isConstant())
557  proto->var(Form("beta_%s",it->first.c_str()))->setConstant(true);
558  else
559  proto->var(Form("beta_%s",it->first.c_str()))->setConstant(false);
560  // set alpha const status to true
561  // proto->var(Form("alpha_%s",it->first.c_str()))->setConstant(true);
562 
563  // replace alphas with alphaOfBeta and replace constraints
564  cout << "alpha_"+it->first+"Constraint=beta_" + it->first+ "Constraint" << endl;
565  editList+=precede + "alpha_"+it->first+"Constraint=beta_" + it->first+ "Constraint";
566  precede=",";
567  cout << "alpha_"+it->first+"=alphaOfBeta_"+ it->first << endl;
568  editList+=precede + "alpha_"+it->first+"=alphaOfBeta_"+ it->first;
569 
570  if( proto->pdf(("alpha_"+it->first+"Constraint").c_str()) && proto->var(("alpha_"+it->first).c_str()) )
571  cout << " checked they are there" << proto->pdf(("alpha_"+it->first+"Constraint").c_str()) << " " << proto->var(("alpha_"+it->first).c_str()) << endl;
572  else
573  cout << "NOT THERE" << endl;
574 
575  // EDIT seems to die if the list of edits is too long. So chunck them up.
576  if(numReplacements%10 == 0 && numReplacements+nskipped!=gammaSyst.size()){
577  edit="EDIT::"+lastPdf+"_("+lastPdf+","+editList+")";
578  lastPdf+="_"; // append an underscore for the edit
579  editList=""; // reset edit list
580  precede="";
581  cout << edit<< endl;
582  proto->factory( edit.c_str() );
583  RooAbsPdf* newOne = proto->pdf(lastPdf.c_str());
584  if(!newOne)
585  cout << "\n\n ---------------------\n WARNING: failed to make EDIT\n\n" << endl;
586 
587  }
588  }
589 
590  /////////////////////////////////////////
591  ////////////////////////////////////
592 
593 
594  // add lognormal terms and their constraints
595  for(it=logNormSyst.begin(); it!=logNormSyst.end(); ++it) {
596  cout << "edit for " << it->first << "with rel uncert = " << it->second << endl;
597  if(! proto->var(("alpha_"+it->first).c_str())){
598  cout << "systematic not there" << endl;
599  nskipped++;
600  continue;
601  }
602  numReplacements++;
603 
604  double relativeUncertainty = it->second;
605  double kappa = 1+relativeUncertainty;
606  // when transforming beta -> alpha, need alpha=1 to be +1sigma value.
607  // the P(beta>kappa*\hat(beta)) = 16%
608  // and \hat(beta) is 1, thus
609  double scale = relativeUncertainty;
610  //double scale = kappa;
611 
612  // this is the LogNormal
613  proto->factory(Form("beta_%s[1,0,10]",it->first.c_str()));
614  proto->factory(Form("kappa_%s[%f]",it->first.c_str(),kappa));
615  proto->factory(Form("Lognormal::beta_%sConstraint(beta_%s,one[1],kappa_%s)",
616  it->first.c_str(),
617  it->first.c_str(),
618  it->first.c_str())) ;
619  proto->factory(Form("PolyVar::alphaOfBeta_%s(beta_%s,{%f,%f})",it->first.c_str(),it->first.c_str(),-1./scale,1./scale));
620  // proto->factory(Form("PolyVar::alphaOfBeta_%s(beta_%s,{%f,%f})",it->first.c_str(),it->first.c_str(),-1.,1./scale));
621 
622  // set beta const status to be same as alpha
623  if(proto->var(Form("alpha_%s",it->first.c_str()))->isConstant())
624  proto->var(Form("beta_%s",it->first.c_str()))->setConstant(true);
625  else
626  proto->var(Form("beta_%s",it->first.c_str()))->setConstant(false);
627  // set alpha const status to true
628  // proto->var(Form("alpha_%s",it->first.c_str()))->setConstant(true);
629 
630  // replace alphas with alphaOfBeta and replace constraints
631  cout << "alpha_"+it->first+"Constraint=beta_" + it->first+ "Constraint" << endl;
632  editList+=precede + "alpha_"+it->first+"Constraint=beta_" + it->first+ "Constraint";
633  precede=",";
634  cout << "alpha_"+it->first+"=alphaOfBeta_"+ it->first << endl;
635  editList+=precede + "alpha_"+it->first+"=alphaOfBeta_"+ it->first;
636 
637  if( proto->pdf(("alpha_"+it->first+"Constraint").c_str()) && proto->var(("alpha_"+it->first).c_str()) )
638  cout << " checked they are there" << proto->pdf(("alpha_"+it->first+"Constraint").c_str()) << " " << proto->var(("alpha_"+it->first).c_str()) << endl;
639  else
640  cout << "NOT THERE" << endl;
641 
642  // EDIT seems to die if the list of edits is too long. So chunck them up.
643  if(numReplacements%10 == 0 && numReplacements+nskipped!=gammaSyst.size()){
644  edit="EDIT::"+lastPdf+"_("+lastPdf+","+editList+")";
645  lastPdf+="_"; // append an underscore for the edit
646  editList=""; // reset edit list
647  precede="";
648  cout << edit<< endl;
649  proto->factory( edit.c_str() );
650  RooAbsPdf* newOne = proto->pdf(lastPdf.c_str());
651  if(!newOne)
652  cout << "\n\n ---------------------\n WARNING: failed to make EDIT\n\n" << endl;
653 
654  }
655  }
656 
657  /////////////////////////////////////////
658  ////////////////////////////////////
659 
660  // commit last bunch of edits
661  edit="EDIT::newSimPdf("+lastPdf+","+editList+")";
662  cout << edit<< endl;
663  proto->factory( edit.c_str() );
664  // proto->writeToFile(("results/model_"+fRowTitle+"_edited.root").c_str());
665  RooAbsPdf* newOne = proto->pdf("newSimPdf");
666  if(newOne){
667  // newOne->graphVizTree(("results/"+pdfName+"_"+fRowTitle+"newSimPdf.dot").c_str());
668  combined_config->SetPdf(*newOne);
669  }
670  else{
671  cout << "\n\n ---------------------\n WARNING: failed to make EDIT\n\n" << endl;
672  }
673  }
674 
676  // FILE * pFile;
677  pFile = fopen ((filename).c_str(),"w");
678 
679 
680  TIter iti = params->createIterator();
681  TIter itj = params->createIterator();
682  RooRealVar *myargi, *myargj;
683  fprintf(pFile," ") ;
684  while ((myargi = (RooRealVar *)iti.Next())) {
685  if(myargi->isConstant()) continue;
686  fprintf(pFile," & %s", myargi->GetName());
687  }
688  fprintf(pFile,"\\\\ \\hline \n" );
689  iti.Reset();
690  while ((myargi = (RooRealVar *)iti.Next())) {
691  if(myargi->isConstant()) continue;
692  fprintf(pFile,"%s", myargi->GetName());
693  itj.Reset();
694  while ((myargj = (RooRealVar *)itj.Next())) {
695  if(myargj->isConstant()) continue;
696  cout << myargi->GetName() << "," << myargj->GetName();
697  fprintf(pFile, " & %.2f", result->correlation(*myargi, *myargj));
698  }
699  cout << endl;
700  fprintf(pFile, " \\\\\n");
701  }
702  fclose(pFile);
703 
704  }
705 
706 
707  ///////////////////////////////////////////////
708  RooWorkspace* HistoToWorkspaceFactory::MakeSingleChannelModel(vector<EstimateSummary> summary, vector<string> systToFix, bool doRatio)
709  {
710 
711  if (summary.empty() ) {
712  Error("MakeSingleChannelModel","vector of EstimateSummry is empty - return a nullptr");
713  return 0;
714  }
715 
716  // to time the macro
717  TStopwatch t;
718  t.Start();
719  string channel=summary[0].channel;
720  cout << "\n\n-------------------\nStarting to process " << channel << " channel" << endl;
721 
722  //
723  // our main workspace that we are using to construct the model
724  //
725  RooWorkspace* proto = new RooWorkspace("proto","proto workspace");
726  ModelConfig * proto_config = new ModelConfig("ModelConfig", proto);
727  proto_config->SetWorkspace(*proto);
728 
729  RooArgSet likelihoodTerms("likelihoodTerms");
730  vector<string> likelihoodTermNames, totSystTermNames,syst_x_expectedPrefixNames, normalizationNames;
731 
732  string prefix, range;
733 
734  /////////////////////////////
735  // Make observables, set values to observed data if data is specified,
736  // otherwise use expected "Asimov" data
737  if (summary.at(0).name=="Data") {
738  ProcessExpectedHisto(summary.at(0).nominal,proto,"obsN","","",0,100000,fLowBin,fHighBin);
739  } else {
740  cout << "Will use expected (\"Asimov\") data set" << endl;
741  ProcessExpectedHisto(NULL,proto,"obsN","","",0,100000,fLowBin,fHighBin);
742  }
743 
744 
745 
746  /////////////////////////////
747  // shared parameters
748  // this is ratio of lumi to nominal lumi. We will include relative uncertainty in model
749  std::stringstream lumiStr;
750  // lumi range
751  lumiStr<<"["<<fNomLumi<<",0,"<<10.*fNomLumi<<"]";
752  proto->factory(("Lumi"+lumiStr.str()).c_str());
753  cout << "lumi str = " << lumiStr.str() << endl;
754 
755  std::stringstream lumiErrorStr;
756  // lumiErrorStr << "nominalLumi["<<fNomLumi << "]," << fLumiError ;
757  lumiErrorStr << "nominalLumi["<<fNomLumi << ",0,"<<fNomLumi+10*fLumiError<<"]," << fLumiError ;
758  proto->factory(("Gaussian::lumiConstraint(Lumi,"+lumiErrorStr.str()+")").c_str());
759  proto->var("nominalLumi")->setConstant();
760  proto->defineSet("globalObservables","nominalLumi");
761  likelihoodTermNames.push_back("lumiConstraint");
762  cout << "lumi Error str = " << lumiErrorStr.str() << endl;
763 
764  //proto->factory((string("SigXsecOverSM[1.,0.5,1..8]").c_str()));
765  ///////////////////////////////////
766  // loop through estimates, add expectation, floating bin predictions,
767  // and terms that constrain floating to expectation via uncertainties
768  vector<EstimateSummary>::iterator it = summary.begin();
769  for(; it!=summary.end(); ++it){
770  if(it->name=="Data") continue;
771 
772  string overallSystName = it->name+"_"+it->channel+"_epsilon";
773  string systSourcePrefix = "alpha_";
774  AddEfficiencyTerms(proto,systSourcePrefix, overallSystName,
775  it->overallSyst,
776  likelihoodTermNames, totSystTermNames);
777 
778  overallSystName=AddNormFactor(proto, channel, overallSystName, *it, doRatio);
779  // get histogram
780  TH1* nominal = it->nominal;
781  if(it->lowHists.size() == 0){
782  cout << it->name+"_"+it->channel+" has no variation histograms " <<endl;
783  string expPrefix=it->name+"_"+it->channel+"_expN";
784  string syst_x_expectedPrefix=it->name+"_"+it->channel+"_overallSyst_x_Exp";
785  ProcessExpectedHisto(nominal,proto,expPrefix,syst_x_expectedPrefix,overallSystName,atoi(NoHistConst_Low),atoi(NoHistConst_High),fLowBin,fHighBin);
786  syst_x_expectedPrefixNames.push_back(syst_x_expectedPrefix);
787  } else if(it->lowHists.size() != it->highHists.size()){
788  cout << "problem in "+it->name+"_"+it->channel
789  << " number of low & high variation histograms don't match" << endl;
790  return 0;
791  } else {
792  string constraintPrefix = it->name+"_"+it->channel+"_Hist_alpha"; // name of source for variation
793  string syst_x_expectedPrefix = it->name+"_"+it->channel+"_overallSyst_x_HistSyst";
794  LinInterpWithConstraint(proto, nominal, it->lowHists, it->highHists, it->systSourceForHist,
795  constraintPrefix, syst_x_expectedPrefix, overallSystName,
796  fLowBin, fHighBin, likelihoodTermNames);
797  syst_x_expectedPrefixNames.push_back(syst_x_expectedPrefix);
798  }
799 
800  // AddMultiVarGaussConstraint(proto, "exp"+it->first+"N", fLowBin, fHighBin, likelihoodTermNames);
801 
802  if(it->normName=="")
803  normalizationNames.push_back( "Lumi" );
804  else
805  normalizationNames.push_back( it->normName);
806  }
807  //proto->Print();
808 
809  ///////////////////////////////////
810  // for ith bin calculate totN_i = lumi * sum_j expected_j * syst_j
811  MakeTotalExpected(proto,channel+"_totN",channel+"_expN","Lumi",fLowBin,fHighBin,
812  syst_x_expectedPrefixNames, normalizationNames);
813 
814  /////////////////////////////////
815  // Relate observables to expected for each bin
816  AddPoissonTerms(proto, "Pois_"+channel, "obsN", channel+"_totN", fLowBin, fHighBin, likelihoodTermNames);
817 
818  /////////////////////////////////
819  // if no data histogram provided, make asimov data
820  if(summary.at(0).name!="Data"){
821  SetObsToExpected(proto, "obsN",channel+"_totN", fLowBin, fHighBin);
822  cout << " using asimov data" << endl;
823  } else{
824  SetObsToExpected(proto, "obsN","obsN", fLowBin, fHighBin);
825  cout << " using input data histogram" << endl;
826  }
827 
828  //////////////////////////////////////
829  // fix specified parameters
830  for(unsigned int i=0; i<systToFix.size(); ++i){
831  RooRealVar* temp = proto->var((systToFix.at(i)).c_str());
832  if(temp) temp->setConstant();
833  else cout << "could not find variable " << systToFix.at(i) << " could not set it to constant" << endl;
834  }
835 
836  //////////////////////////////////////
837  // final proto model
838  for(unsigned int i=0; i<likelihoodTermNames.size(); ++i){
839  // cout << likelihoodTermNames[i] << endl;
840  likelihoodTerms.add(* (proto->arg(likelihoodTermNames[i].c_str())) );
841  }
842  // likelihoodTerms.Print();
843 
844  proto->defineSet("likelihoodTerms",likelihoodTerms);
845  // proto->Print();
846 
847  cout <<"-----------------------------------------"<<endl;
848  cout <<"import model into workspace" << endl;
849  RooProdPdf* model = new RooProdPdf(("model_"+channel).c_str(),
850  "product of Poissons accross bins for a single channel",
851  likelihoodTerms);
852  proto->import(*model,RecycleConflictNodes());
853 
854  proto_config->SetPdf(*model);
855  proto_config->SetGlobalObservables(*proto->set("globalObservables"));
856 
857  proto->import(*proto_config,proto_config->GetName());
858  proto->importClassCode();
859  // proto->writeToFile(("results/model_"+channel+".root").c_str());
860 
861  return proto;
862  }
863 
864  RooWorkspace* HistoToWorkspaceFactory::MakeCombinedModel(vector<string> ch_names, vector<RooWorkspace*> chs)
865  {
866 
867  //
868  /// These things were used for debugging. Maybe useful in the future
869  //
870  // RooAbsReal::defaultIntegratorConfig()->setEpsAbs(1e-8) ;
871  // RooAbsReal::defaultIntegratorConfig()->setEpsRel(1e-8) ;
872  // RooMsgService::instance().setGlobalKillBelow(RooMsgService::WARNING);
873  // RooMsgService::instance().setGlobalKillBelow(RooMsgService::WARNING) ;
874  // cout << "MsgSvc: " << RooMsgService::instance().globalKillBelow() << " INFO "
875  // << RooMsgService::INFO << " WARNING " << RooMsgService::WARNING << endl;
876 
877  // RooArgSet* constrainedParams= new RooArgSet("constrainedParams");
878 
879  // check inputs (see JIRA-6890)
880  if (ch_names.empty() || chs.empty() ) {
881  Error("MakeCombinedModel","Input vectors are empty - return a nullptr");
882  return 0;
883  }
884  if (chs.size() < ch_names.size() ) {
885  Error("MakeCombinedModel","Input vector of workspace has an invalid size - return a nullptr");
886  return 0;
887  }
888 
889  map<string, RooAbsPdf*> pdfMap;
890  vector<RooAbsPdf*> models;
891  stringstream ss;
892 
893  RooArgSet globalObs;
894  for(unsigned int i = 0; i< ch_names.size(); ++i){
895  string channel_name=ch_names[i];
896 
897  if (ss.str().empty()) ss << channel_name ;
898  else ss << ',' << channel_name ;
899  RooWorkspace * ch=chs[i];
900 
901  RooAbsPdf* model = ch->pdf(("model_"+channel_name).c_str());
902  models.push_back(model);
903  globalObs.add(*ch->set("globalObservables"));
904 
905  // constrainedParams->add( * ch->set("constrainedParams") );
906  pdfMap[channel_name]=model;
907  }
908  //constrainedParams->Print();
909 
910  cout << "\n\n------------------\n Entering combination" << endl;
911  RooWorkspace* combined = new RooWorkspace("combined");
912 
913  RooCategory* channelCat = (RooCategory*) combined->factory(("channelCat["+ss.str()+"]").c_str());
914  RooSimultaneous * simPdf= new RooSimultaneous("simPdf","",pdfMap, *channelCat);
915  ModelConfig * combined_config = new ModelConfig("ModelConfig", combined);
916  combined_config->SetWorkspace(*combined);
917  // combined_config->SetNuisanceParameters(*constrainedParams);
918  combined->import(globalObs);
919  combined->defineSet("globalObservables",globalObs);
920  combined_config->SetGlobalObservables(*combined->set("globalObservables"));
921 
922  ////////////////////////////////////////////
923  // Make toy simultaneous dataset
924  cout <<"-----------------------------------------"<<endl;
925  cout << "create toy data for " << ss.str() << endl;
926 
927  const RooArgSet* obsN = chs[0]->set("obsN");
928 
929  RooDataSet * simData=new RooDataSet("simData","master dataset", *obsN,
930  Index(*channelCat), Import(ch_names[0].c_str(),*((RooDataSet*)chs[0]->data("expData"))));
931  for(unsigned int i = 1; i< ch_names.size(); ++i){
932  RooDataSet * simData_ch=new RooDataSet("simData","master dataset", *obsN,
933  Index(*channelCat), Import(ch_names[i].c_str(),*((RooDataSet*)chs[i]->data("expData"))));
934  simData->append(*simData_ch);
935  }
936  //for(int i=0; i<simData->numEntries(); ++i)
937  // simData->get(i)->Print("v");
938 
939  combined->import(*simData,RecycleConflictNodes());
940 
941  cout << "\n\n----------------\n Importing combined model" << endl;
942  combined->import(*simPdf,RecycleConflictNodes());
943  //combined->import(*simPdf, RenameVariable("SigXsecOverSM","SigXsecOverSM_comb"));
944  cout << "check pointer " << simPdf << endl;
945 
946  for(unsigned int i=0; i<fSystToFix.size(); ++i){
947  // make sure they are fixed
948  RooRealVar* temp = combined->var((fSystToFix.at(i)).c_str());
949  if(temp) {
950  temp->setConstant();
951  cout <<"setting " << fSystToFix.at(i) << " constant" << endl;
952  }
953  else cout << "could not find variable " << fSystToFix.at(i) << " could not set it to constant" << endl;
954  }
955 
956  ///
957  /// writing out the model in graphViz
958  ///
959  // RooAbsPdf* customized=combined->pdf("simPdf");
960  //combined_config->SetPdf(*customized);
961  combined_config->SetPdf(*simPdf);
962  // customized->graphVizTree(("results/"+fResultsPrefixStr.str()+"_simul.dot").c_str());
963  combined->import(*combined_config,combined_config->GetName());
964  combined->importClassCode();
965  // combined->writeToFile("results/model_combined.root");
966 
967  return combined;
968  }
969 
970  ///////////////////////////////////////////////
971  void HistoToWorkspaceFactory::FitModel(RooWorkspace * combined, string channel, string /*model_name*/, string data_name, bool /*doParamInspect*/)
972  {
973 
974  ModelConfig * combined_config = (ModelConfig *) combined->obj("ModelConfig");
975  RooDataSet * simData = (RooDataSet *) combined->obj(data_name.c_str());
976  // const RooArgSet * constrainedParams=combined_config->GetNuisanceParameters();
977  const RooArgSet * POIs=combined_config->GetParametersOfInterest();
978 
979  /*
980  RooRealVar* poi = (RooRealVar*) combined->var("SigXsecOverSM");
981  RooArgSet * params= new RooArgSet;
982  params->add(*poi);
983  combined_config->SetParameters(*params);
984 
985  RooAbsData* expData = combined->data("expData");
986  RooArgSet* temp = (RooArgSet*) combined->set("obsN")->Clone("temp");
987  temp->add(*poi);
988  RooAbsPdf* model=combined_config->GetPdf();
989  RooArgSet* constrainedParams = model->getParameters(temp);
990  combined->defineSet("constrainedParams", *constrainedParams);
991  */
992 
993  //RooAbsPdf* model=combined->pdf(model_name.c_str());
994  RooAbsPdf* model=combined_config->GetPdf();
995  // RooArgSet* allParams = model->getParameters(*simData);
996 
997  ///////////////////////////////////////
998  //Do combined fit
999  //RooMsgService::instance().setGlobalKillBelow(RooMsgService::INFO) ;
1000  cout << "\n\n---------------" << endl;
1001  cout << "---------------- Doing "<< channel << " Fit" << endl;
1002  cout << "---------------\n\n" << endl;
1003  // RooFitResult* result = model->fitTo(*simData, Minos(kTRUE), Save(kTRUE), PrintLevel(1));
1004  model->fitTo(*simData, Minos(kTRUE), PrintLevel(1));
1005  // PrintCovarianceMatrix(result, allParams, "results/"+FilePrefixStr(channel)+"_corrMatrix.table" );
1006 
1007  //
1008  // assuming there is only on poi
1009  //
1010  RooRealVar* poi = 0; // (RooRealVar*) POIs->first();
1011  // for results tables
1012  TIterator* params_itr=POIs->createIterator();
1013  TObject* params_obj=0;
1014  while((params_obj=params_itr->Next())){
1015  poi = (RooRealVar*) params_obj;
1016  cout << "printing results for " << poi->GetName() << " at " << poi->getVal()<< " high " << poi->getErrorLo() << " low " << poi->getErrorHi()<<endl;
1017  }
1018  fprintf(pFile, " %.4f / %.4f ", poi->getErrorLo(), poi->getErrorHi());
1019 
1020  RooAbsReal* nll = model->createNLL(*simData);
1021  RooAbsReal* profile = nll->createProfile(*poi);
1022  RooPlot* frame = poi->frame();
1023  FormatFrameForLikelihood(frame);
1024  TCanvas* c1 = new TCanvas( channel.c_str(), "",800,600);
1025  nll->plotOn(frame, ShiftToZero(), LineColor(kRed), LineStyle(kDashed));
1026  profile->plotOn(frame);
1027  frame->SetMinimum(0);
1028  frame->SetMaximum(2.);
1029  frame->Draw();
1030  // c1->SaveAs( ("results/"+FilePrefixStr(channel)+"_profileLR.eps").c_str() );
1031  c1->SaveAs( (fFileNamePrefix+"_"+channel+"_"+fRowTitle+"_profileLR.eps").c_str() );
1032 
1033  fOut_f->mkdir(channel.c_str())->mkdir("Summary")->cd();
1034 
1035  // an example of calculating profile for a nuisance parameter not poi
1036  /*
1037  RooRealVar* alpha_isrfsr = (RooRealVar*) combined->var("alpha_isrfsr");
1038  RooAbsReal* profile_isrfsr = nll->createProfile(*alpha_isrfsr);
1039  poi->setVal(0.55);
1040  poi->setConstant();
1041 
1042  RooPlot* frame_isrfsr = alpha_isrfsr->frame();
1043  profile_isrfsr->plotOn(frame_isrfsr, Precision(0.1));
1044  TCanvas c_isrfsr = new TCanvas( "combined", "",800,600);
1045  FormatFrameForLikelihood(frame_isrfsr, "alpha_{isrfsr}");
1046  frame_isrfsr->Draw();
1047  fOut_f->cd("Summary");
1048  c1->Write((FilePrefixStr(channel).str()+"_profileLR_alpha_isrfsr").c_str() );
1049  delete frame; delete c1;
1050  poi->setConstant(kFALSE);
1051  */
1052 
1053  RooCurve* curve=frame->getCurve();
1054  Int_t curve_N=curve->GetN();
1055  Double_t* curve_x=curve->GetX();
1056  delete frame; delete c1;
1057 
1058  //
1059  // Verbose output from MINUIT
1060  //
1061  /*
1062  RooMsgService::instance().setGlobalKillBelow(RooFit::DEBUG) ;
1063  profile->getVal();
1064  RooMinuit* minuit = ((RooProfileLL*) profile)->minuit();
1065  minuit->setPrintLevel(5) ; // Print MINUIT messages
1066  minuit->setVerbose(5) ; // Print RooMinuit messages with parameter
1067  // changes (corresponds to the Verbose() option of fitTo()
1068  */
1069 
1070  Double_t * x_arr = new Double_t[curve_N];
1071  Double_t * y_arr_nll = new Double_t[curve_N];
1072 // Double_t y_arr_prof_nll[curve_N];
1073 // Double_t y_arr_prof[curve_N];
1074 
1075  for(int i=0; i<curve_N; i++){
1076  double f=curve_x[i];
1077  poi->setVal(f);
1078  x_arr[i]=f;
1079  y_arr_nll[i]=nll->getVal();
1080  }
1081  TGraph * g = new TGraph(curve_N, x_arr, y_arr_nll);
1082  g->SetName((FilePrefixStr(channel)+"_nll").c_str());
1083  g->Write();
1084  delete g;
1085  delete [] x_arr;
1086  delete [] y_arr_nll;
1087 
1088  /** find out what's inside the workspace **/
1089  //combined->Print();
1090 
1091  }
1092 
1093 
1094 void HistoToWorkspaceFactory::FormatFrameForLikelihood(RooPlot* frame, string /*XTitle*/, string YTitle){
1095 
1098  gStyle->SetPadColor(0);
1099  gStyle->SetCanvasColor(255);
1100  gStyle->SetTitleFillColor(255);
1102  gStyle->SetStatColor(255);
1103 
1104  RooAbsRealLValue* var = frame->getPlotVar();
1105  double xmin = var->getMin();
1106  double xmax = var->getMax();
1107 
1108  frame->SetTitle("");
1109  // frame->GetXaxis()->SetTitle(XTitle.c_str());
1110  frame->GetXaxis()->SetTitle(var->GetTitle());
1111  frame->GetYaxis()->SetTitle(YTitle.c_str());
1112  frame->SetMaximum(2.);
1113  frame->SetMinimum(0.);
1114  TLine * line = new TLine(xmin,.5,xmax,.5);
1115  line->SetLineColor(kGreen);
1116  TLine * line90 = new TLine(xmin,2.71/2.,xmax,2.71/2.);
1117  line90->SetLineColor(kGreen);
1118  TLine * line95 = new TLine(xmin,3.84/2.,xmax,3.84/2.);
1119  line95->SetLineColor(kGreen);
1120  frame->addObject(line);
1121  frame->addObject(line90);
1122  frame->addObject(line95);
1123  }
1124 
1126  if(! file) return file;
1127  string path="";
1128  TDirectory* ptr=0;
1129  for(vector<string>::iterator itr=names.begin(); itr != names.end(); ++itr){
1130  if( ! path.empty() ) path+="/";
1131  path+=(*itr);
1132  ptr=file->GetDirectory(path.c_str());
1133  if( ! ptr ) ptr=file->mkdir((*itr).c_str());
1134  file=file->GetDirectory(path.c_str());
1135  }
1136  return ptr;
1137  }
1139  if(! file) return file;
1140  TDirectory* ptr=0;
1141  ptr=file->GetDirectory(name.c_str());
1142  if( ! ptr ) ptr=file->mkdir(name.c_str());
1143  return ptr;
1144  }
1145 
1146 }
1147 }
1148 
virtual RooAbsReal * createNLL(RooAbsData &data, const RooLinkedList &cmdList)
Construct representation of -log(L) of PDFwith given dataset.
Definition: RooAbsPdf.cxx:778
virtual Double_t getMin(const char *name=0) const
virtual const char * GetName() const
Returns name of object.
Definition: TNamed.h:47
virtual Int_t Write(const char *name=0, Int_t option=0, Int_t bufsize=0)
Write this object to the current directory.
Definition: TObject.cxx:785
TIterator * createIterator(Bool_t dir=kIterForward) const
virtual Bool_t add(const RooAbsArg &var, Bool_t silent=kFALSE)
Add the specified argument to list.
virtual RooPlot * plotOn(RooPlot *frame, const RooCmdArg &arg1=RooCmdArg(), const RooCmdArg &arg2=RooCmdArg(), const RooCmdArg &arg3=RooCmdArg(), const RooCmdArg &arg4=RooCmdArg(), const RooCmdArg &arg5=RooCmdArg(), const RooCmdArg &arg6=RooCmdArg(), const RooCmdArg &arg7=RooCmdArg(), const RooCmdArg &arg8=RooCmdArg(), const RooCmdArg &arg9=RooCmdArg(), const RooCmdArg &arg10=RooCmdArg()) const
Plot (project) PDF on specified frame.
float xmin
Definition: THbookFile.cxx:93
void MakeTotalExpected(RooWorkspace *proto, std::string totName, std::string, std::string, int lowBin, int highBin, std::vector< std::string > &syst_x_expectedPrefixNames, std::vector< std::string > &normByNames)
ModelConfig is a simple class that holds configuration information specifying how a model should be u...
Definition: ModelConfig.h:30
A RooCurve is a one-dimensional graphical representation of a real-valued function.
Definition: RooCurve.h:32
virtual Double_t getMax(const char *name=0) const
virtual Bool_t add(const RooAbsCollection &col, Bool_t silent=kFALSE)
Add a collection of arguments to this collection by calling add() for each element in the source coll...
Definition: RooArgSet.h:86
TAxis * GetYaxis() const
Definition: RooPlot.cxx:1117
void Start(Bool_t reset=kTRUE)
Start the stopwatch.
Definition: TStopwatch.cxx:58
void EditSyst(RooWorkspace *proto, const char *pdfNameChar, std::map< std::string, double > gammaSyst, std::map< std::string, double > uniformSyst, std::map< std::string, double > logNormSyst)
cout << "in edit, gammamap.size = " << gammaSyst.size() << ", unimap.size = " << uniformSyst.size() << endl;
void SetStatColor(Color_t color=19)
Definition: TStyle.h:367
void SetObsToExpected(RooWorkspace *proto, std::string obsPrefix, std::string expPrefix, int lowBin, int highBin)
TAxis * GetXaxis() const
Definition: RooPlot.cxx:1116
RooCmdArg LineColor(Color_t color)
TLine * line
virtual void SetWorkspace(RooWorkspace &ws)
Definition: ModelConfig.h:66
void addObject(TObject *obj, Option_t *drawOptions="", Bool_t invisible=kFALSE)
Add a generic object to this plot.
Definition: RooPlot.cxx:391
return c1
Definition: legend1.C:41
Definition: Rtypes.h:59
#define g(i)
Definition: RSha256.hxx:105
RooCmdArg PrintLevel(Int_t code)
RooProdPdf is an efficient implementation of a product of PDFs of the form.
Definition: RooProdPdf.h:31
Double_t getVal(const RooArgSet *set=0) const
Definition: RooAbsReal.h:64
RooCmdArg Minos(Bool_t flag=kTRUE)
R__EXTERN TStyle * gStyle
Definition: TStyle.h:406
virtual Double_t GetBinContent(Int_t bin) const
Return content of bin number bin.
Definition: TH1.cxx:4770
std::string AddNormFactor(RooWorkspace *, std::string &, std::string &, EstimateSummary &, bool)
RooAbsReal * function(const char *name) const
Retrieve function (RooAbsReal) with given name. Note that all RooAbsPdfs are also RooAbsReals...
#define f(i)
Definition: RSha256.hxx:104
int Int_t
Definition: RtypesCore.h:41
virtual TDirectory * mkdir(const char *name, const char *title="")
Create a sub-directory "a" or a hierarchy of sub-directories "a/b/c/...".
void SetTitle(const char *name)
Set the title of the RooPlot to &#39;title&#39;.
Definition: RooPlot.cxx:1098
RooAbsArg * arg(const char *name) const
Return RooAbsArg with given name. A null pointer is returned if none is found.
Definition: Rtypes.h:59
STL namespace.
RooCmdArg RecycleConflictNodes(Bool_t flag=kTRUE)
void PrintCovarianceMatrix(RooFitResult *result, RooArgSet *params, std::string filename)
#define NoHistConst_Low
void ProcessExpectedHisto(TH1 *hist, RooWorkspace *proto, std::string prefix, std::string productPrefix, std::string systTerm, double low, double high, int lowBin, int highBin)
Iterator abstract base class.
Definition: TIterator.h:30
#define alpha_Low
void Reset()
Definition: TCollection.h:252
void SetCanvasColor(Color_t color=19)
Definition: TStyle.h:321
Bool_t importClassCode(const char *pat="*", Bool_t doReplace=kFALSE)
Inport code of all classes in the workspace that have a class name that matches pattern &#39;pat&#39; and whi...
double sqrt(double)
virtual void SetMinimum(Double_t minimum=-1111)
Set minimum value of Y axis.
Definition: RooPlot.cxx:958
Double_t getErrorLo() const
Definition: RooRealVar.h:62
void SetFrameFillColor(Color_t color=1)
Definition: TStyle.h:349
virtual void SetName(const char *name="")
Set graph name.
Definition: TGraph.cxx:2207
void FormatFrameForLikelihood(RooPlot *frame, std::string XTitle=std::string("#sigma / #sigma_{SM}"), std::string YTitle=std::string("-log likelihood"))
RooCmdArg LineStyle(Style_t style)
double pow(double, double)
virtual void SetPdf(const RooAbsPdf &pdf)
Set the Pdf, add to the the workspace if not already there.
Definition: ModelConfig.h:75
void LinInterpWithConstraint(RooWorkspace *proto, TH1 *nominal, std::vector< TH1 *> lowHist, std::vector< TH1 *> highHist, std::vector< std::string > sourceName, std::string prefix, std::string productPrefix, std::string systTerm, int lowBin, int highBin, std::vector< std::string > &likelihoodTermNames)
RooConstVar represent a constant real-valued object.
Definition: RooConstVar.h:25
RooCmdArg ShiftToZero()
void AddEfficiencyTerms(RooWorkspace *proto, std::string prefix, std::string interpName, std::map< std::string, std::pair< double, double > > systMap, std::vector< std::string > &likelihoodTermNames, std::vector< std::string > &totSystTermNames)
RooRealVar represents a fundamental (non-derived) real valued object.
Definition: RooRealVar.h:36
virtual void setVal(Double_t value)
Set value of variable to &#39;value&#39;.
Definition: RooRealVar.cxx:204
virtual void SetLineColor(Color_t lcolor)
Set the line color.
Definition: TAttLine.h:40
VecExpr< UnaryOp< Fabs< T >, VecExpr< A, T, D >, T >, T, D > fabs(const VecExpr< A, T, D > &rhs)
void SetPadBorderMode(Int_t mode=1)
Definition: TStyle.h:334
void SetCanvasBorderMode(Int_t mode=1)
Definition: TStyle.h:323
void SetPadColor(Color_t color=19)
Definition: TStyle.h:332
void setConstant(Bool_t value=kTRUE)
virtual void SetMaximum(Double_t maximum=-1111)
Set maximum value of Y axis.
Definition: RooPlot.cxx:948
void Customize(RooWorkspace *proto, const char *pdfNameChar, std::map< std::string, std::string > renameMap)
Bool_t defineSet(const char *name, const RooArgSet &aset, Bool_t importMissing=kFALSE)
Define a named RooArgSet with given constituents.
TObject * Next()
Definition: TCollection.h:249
void AddMultiVarGaussConstraint(RooWorkspace *proto, std::string prefix, int lowBin, int highBin, std::vector< std::string > &likelihoodTermNames)
virtual void Error(const char *method, const char *msgfmt,...) const
Issue error message.
Definition: TObject.cxx:880
char * Form(const char *fmt,...)
A simple line.
Definition: TLine.h:23
Int_t GetN() const
Definition: TGraph.h:122
float xmax
Definition: THbookFile.cxx:93
RooDataSet is a container class to hold unbinned data.
Definition: RooDataSet.h:29
RooCategory represents a fundamental (non-derived) discrete value object.
Definition: RooCategory.h:24
RooPlot * frame(const RooCmdArg &arg1, const RooCmdArg &arg2=RooCmdArg::none(), const RooCmdArg &arg3=RooCmdArg::none(), const RooCmdArg &arg4=RooCmdArg::none(), const RooCmdArg &arg5=RooCmdArg::none(), const RooCmdArg &arg6=RooCmdArg::none(), const RooCmdArg &arg7=RooCmdArg::none(), const RooCmdArg &arg8=RooCmdArg::none()) const
TObject * obj(const char *name) const
Return any type of object (RooAbsArg, RooAbsData or generic object) with given name) ...
A RooPlot is a plot frame and a container for graphics objects within that frame. ...
Definition: RooPlot.h:41
Double_t * GetX() const
Definition: TGraph.h:129
RooCmdArg Import(const char *state, TH1 &histo)
Multivariate Gaussian p.d.f.
The Canvas class.
Definition: TCanvas.h:31
RooCmdArg Index(RooCategory &icat)
#define alpha_High
Namespace for the RooStats classes.
Definition: Asimov.h:20
RooAbsPdf * GetPdf() const
get model PDF (return NULL if pdf has not been specified or does not exist)
Definition: ModelConfig.h:222
const RooArgSet * GetParametersOfInterest() const
get RooArgSet containing the parameter of interest (return NULL if not existing)
Definition: ModelConfig.h:225
TDirectory * Mkdir(TDirectory *file, std::string name)
TDirectory * Makedirs(TDirectory *file, std::vector< std::string > names)
#define ClassImp(name)
Definition: Rtypes.h:359
RooCurve * getCurve(const char *name=0) const
Return a RooCurve pointer of the named object in this plot, or zero if the named object does not exis...
Definition: RooPlot.cxx:774
RooAbsPdf * pdf(const char *name) const
Retrieve p.d.f (RooAbsPdf) with given name. A null pointer is returned if not found.
double Double_t
Definition: RtypesCore.h:55
RooAbsReal is the common abstract base class for objects that represent a real value and implements f...
Definition: RooAbsReal.h:53
Describe directory structure in memory.
Definition: TDirectory.h:34
void append(RooDataSet &data)
Add all data points of given data set to this data set.
void AddPoissonTerms(RooWorkspace *proto, std::string prefix, std::string obsPrefix, std::string expPrefix, int lowBin, int highBin, std::vector< std::string > &likelihoodTermNames)
RooWorkspace * MakeCombinedModel(std::vector< std::string >, std::vector< RooWorkspace *>)
The TH1 histogram class.
Definition: TH1.h:56
std::vector< NormFactor > normFactor
void FitModel(RooWorkspace *, std::string, std::string, std::string, bool=false)
RooRealVar * var(const char *name) const
Retrieve real-valued variable (RooRealVar) with given name. A null pointer is returned if not found...
RooWorkspace * MakeSingleChannelModel(std::vector< RooStats::HistFactory::EstimateSummary > summary, std::vector< std::string > systToFix, bool doRatio=false)
Mother of all ROOT objects.
Definition: TObject.h:37
RooFactoryWSTool & factory()
Return instance to factory tool.
RooAbsRealLValue is the common abstract base class for objects that represent a real value that may a...
RooAbsPdf is the abstract interface for all probability density functions The class provides hybrid a...
Definition: RooAbsPdf.h:41
Bool_t import(const RooAbsArg &arg, const RooCmdArg &arg1=RooCmdArg(), const RooCmdArg &arg2=RooCmdArg(), const RooCmdArg &arg3=RooCmdArg(), const RooCmdArg &arg4=RooCmdArg(), const RooCmdArg &arg5=RooCmdArg(), const RooCmdArg &arg6=RooCmdArg(), const RooCmdArg &arg7=RooCmdArg(), const RooCmdArg &arg8=RooCmdArg(), const RooCmdArg &arg9=RooCmdArg())
Import a RooAbsArg object, e.g.
RooAbsRealLValue * getPlotVar() const
Definition: RooPlot.h:132
void SetTitleFillColor(Color_t color=1)
Definition: TStyle.h:381
Definition: file.py:1
virtual RooAbsReal * createProfile(const RooArgSet &paramsOfInterest)
Create a RooProfileLL object that eliminates all nuisance parameters in the present function...
Definition: RooAbsReal.cxx:461
virtual TObject * Next()=0
A Graph is a graphics object made of two arrays X and Y with npoints each.
Definition: TGraph.h:41
const RooArgSet * set(const char *name)
Return pointer to previously defined named set with given nmame If no such set is found a null pointe...
virtual TDirectory * GetDirectory(const char *namecycle, Bool_t printError=false, const char *funcname="GetDirectory")
Find a directory using apath.
Definition: TDirectory.cxx:400
const char * proto
Definition: civetweb.c:15049
virtual RooFitResult * fitTo(RooAbsData &data, const RooCmdArg &arg1=RooCmdArg::none(), const RooCmdArg &arg2=RooCmdArg::none(), const RooCmdArg &arg3=RooCmdArg::none(), const RooCmdArg &arg4=RooCmdArg::none(), const RooCmdArg &arg5=RooCmdArg::none(), const RooCmdArg &arg6=RooCmdArg::none(), const RooCmdArg &arg7=RooCmdArg::none(), const RooCmdArg &arg8=RooCmdArg::none())
Fit PDF to given dataset.
Definition: RooAbsPdf.cxx:1079
Definition: tree.py:1
#define NoHistConst_High
virtual void SaveAs(const char *filename="", Option_t *option="") const
Save Pad contents in a file in one of various formats.
Definition: TPad.cxx:5504
RooAbsArg is the common abstract base class for objects that represent a value (of arbitrary type) an...
Definition: RooAbsArg.h:66
virtual void SetGlobalObservables(const RooArgSet &set)
specify the global observables
Definition: ModelConfig.h:160
virtual const char * GetName() const
Returns name of object.
Definition: TObject.cxx:357
double exp(double)
virtual void SetTitle(const char *title="")
Set the title of the TNamed.
Definition: TNamed.cxx:164
RooSimultaneous facilitates simultaneous fitting of multiple PDFs to subsets of a given dataset...
const Bool_t kTRUE
Definition: RtypesCore.h:87
Double_t getErrorHi() const
Definition: RooRealVar.h:63
Double_t correlation(const RooAbsArg &par1, const RooAbsArg &par2) const
Definition: RooFitResult.h:117
char name[80]
Definition: TGX11.cxx:109
Bool_t isConstant() const
Definition: RooAbsArg.h:266
The RooWorkspace is a persistable container for RooFit projects.
Definition: RooWorkspace.h:43
virtual void Draw(Option_t *options=0)
Draw this plot and all of the elements it contains.
Definition: RooPlot.cxx:558
virtual const char * GetTitle() const
Returns title of object.
Definition: TNamed.h:48
Stopwatch class.
Definition: TStopwatch.h:28