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