Logo ROOT   6.08/07
Reference Guide
HistoToWorkspaceFactoryFast.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 #include "RooDataSet.h"
27 #include "RooRealVar.h"
28 #include "RooConstVar.h"
29 #include "RooAddition.h"
30 #include "RooProduct.h"
31 #include "RooProdPdf.h"
32 #include "RooAddPdf.h"
33 #include "RooGaussian.h"
34 #include "RooPoisson.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 "RooHistFunc.h"
47 #include "RooHistPdf.h"
48 #include "RooRealSumPdf.h"
49 #include "RooProduct.h"
50 #include "RooWorkspace.h"
51 #include "RooCustomizer.h"
52 #include "RooPlot.h"
53 #include "RooMsgService.h"
54 #include "RooStats/RooStatsUtils.h"
55 #include "RooStats/ModelConfig.h"
59 
60 #include "TH2F.h"
61 #include "TH3F.h"
62 #include "TFile.h"
63 #include "TCanvas.h"
64 #include "TH1.h"
65 #include "TLine.h"
66 #include "TTree.h"
67 #include "TMarker.h"
68 #include "TStopwatch.h"
69 #include "TROOT.h"
70 #include "TStyle.h"
71 #include "TVectorD.h"
72 #include "TMatrixDSym.h"
73 
74 // specific to this package
79 #include "Helper.h"
80 
81 #include <algorithm>
82 
83 #define VERBOSE
84 
85 #define alpha_Low "-5"
86 #define alpha_High "5"
87 #define NoHistConst_Low "0"
88 #define NoHistConst_High "2000"
89 
90 // use this order for safety on library loading
91 using namespace RooFit ;
92 using namespace RooStats ;
93 using namespace std ;
94 
96 
97 namespace RooStats{
98 namespace HistFactory{
99 
100  HistoToWorkspaceFactoryFast::HistoToWorkspaceFactoryFast() :
101  fNomLumi(1.0), fLumiError(0),
102  fLowBin(0), fHighBin(0)
103  {}
104 
106  }
107 
109  fSystToFix( measurement.GetConstantParams() ),
110  fParamValues( measurement.GetParamValues() ),
111  fNomLumi( measurement.GetLumi() ),
112  fLumiError( measurement.GetLumi()*measurement.GetLumiRelErr() ),
113  fLowBin( measurement.GetBinLow() ),
114  fHighBin( measurement.GetBinHigh() ) {
115 
116  // Set Preprocess functions
118 
119  }
120 
121  void HistoToWorkspaceFactoryFast::ConfigureWorkspaceForMeasurement( const std::string& ModelName, RooWorkspace* ws_single, Measurement& measurement ) {
122 
123  // Configure a workspace by doing any
124  // necessary post-processing and by
125  // creating a ModelConfig
126 
127  // Make a ModelConfig and configure it
128  ModelConfig * proto_config = (ModelConfig *) ws_single->obj("ModelConfig");
129  if( proto_config == NULL ) {
130  std::cout << "Error: Did not find 'ModelConfig' object in file: " << ws_single->GetName()
131  << std::endl;
132  throw hf_exc();
133  }
134 
135  std::vector<std::string> poi_list = measurement.GetPOIList();
136  if( poi_list.size()==0 ) {
137  std::cout << "Warining: No Parametetrs of interest are set" << std::endl;
138  }
139 
140  cout << "Setting Parameter(s) of Interest as: ";
141  for(unsigned int i = 0; i < poi_list.size(); ++i) {
142  cout << poi_list.at(i) << " ";
143  }
144  cout << endl;
145 
146  RooArgSet * params= new RooArgSet;
147  for( unsigned int i = 0; i < poi_list.size(); ++i ) {
148  std::string poi_name = poi_list.at(i);
149  RooRealVar* poi = (RooRealVar*) ws_single->var( poi_name.c_str() );
150  if(poi){
151  params->add(*poi);
152  }
153  else {
154  std::cout << "WARNING: Can't find parameter of interest: " << poi_name
155  << " in Workspace. Not setting in ModelConfig." << std::endl;
156  //throw hf_exc();
157  }
158  }
159  proto_config->SetParametersOfInterest(*params);
160 
161  // Name of an 'edited' model, if necessary
162  std::string NewModelName = "newSimPdf"; // <- This name is hard-coded in HistoToWorkspaceFactoryFast::EditSyt. Probably should be changed to : std::string("new") + ModelName;
163 
164 #ifdef DO_EDIT_WS
165  // Activate Additional Constraint Terms
166  if( measurement.GetGammaSyst().size() > 0
167  || measurement.GetUniformSyst().size() > 0
168  || measurement.GetLogNormSyst().size() > 0
169  || measurement.GetNoSyst().size() > 0) {
170  HistoToWorkspaceFactoryFast::EditSyst( ws_single, (ModelName).c_str(),
171  measurement.GetGammaSyst(),
172  measurement.GetUniformSyst(),
173  measurement.GetLogNormSyst(),
174  measurement.GetNoSyst());
175 
176  proto_config->SetPdf( *ws_single->pdf( "newSimPdf" ) );
177  }
178 #endif
179 
180  // Set the ModelConfig's Params of Interest
181  RooAbsData* expData = ws_single->data("asimovData");
182  if( !expData ) {
183  std::cout << "Error: Failed to find dataset: " << expData
184  << " in workspace" << std::endl;
185  throw hf_exc();
186  }
187  if(poi_list.size()!=0){
188  proto_config->GuessObsAndNuisance(*expData);
189  }
190 
191  // Now, let's loop over any additional asimov datasets
192  // that we need to make
193 
194  // Get the pdf
195  // Notice that we get the "new" pdf, this is the one that is
196  // used in the creation of these asimov datasets since they
197  // are fitted (or may be, at least).
198  RooAbsPdf* pdf = ws_single->pdf(NewModelName.c_str());
199  if( !pdf ) pdf = ws_single->pdf( ModelName.c_str() );
200  const RooArgSet* observables = ws_single->set("observables");
201 
202  // Create a SnapShot of the nominal values
203  std::string SnapShotName = "NominalParamValues";
204  ws_single->saveSnapshot(SnapShotName.c_str(), ws_single->allVars());
205 
206  for( unsigned int i=0; i<measurement.GetAsimovDatasets().size(); ++i) {
207 
208  // Set the variable values and "const" ness with the workspace
209  RooStats::HistFactory::Asimov& asimov = measurement.GetAsimovDatasets().at(i);
210  std::string AsimovName = asimov.GetName();
211 
212  std::cout << "Generating additional Asimov Dataset: " << AsimovName << std::endl;
213  asimov.ConfigureWorkspace(ws_single);
214  RooDataSet* asimov_dataset =
216 
217  std::cout << "Importing Asimov dataset" << std::endl;
218  bool failure = ws_single->import(*asimov_dataset, Rename(AsimovName.c_str()));
219  if( failure ) {
220  std::cout << "Error: Failed to import Asimov dataset: " << AsimovName
221  << std::endl;
222  delete asimov_dataset;
223  throw hf_exc();
224  }
225 
226  // Load the snapshot at the end of every loop iteration
227  // so we start each loop with a "clean" snapshot
228  ws_single->loadSnapshot(SnapShotName.c_str());
229 
230  // we can now deleted the data set after having imported it
231  delete asimov_dataset;
232 
233  }
234 
235  // Cool, we're done
236  return; // ws_single;
237  }
238 
239 
240  // We want to eliminate this interface and use the measurment directly
242 
243  // This is a pretty light-weight wrapper function
244  //
245  // Take a fully configured measurement as well as
246  // one of its channels
247  //
248  // Return a workspace representing that channel
249  // Do this by first creating a vector of EstimateSummary's
250  // and this by configuring the workspace with any post-processing
251 
252  // Get the channel's name
253  string ch_name = channel.GetName();
254 
255  // Create a workspace for a SingleChannel from the Measurement Object
256  RooWorkspace* ws_single = this->MakeSingleChannelWorkspace(measurement, channel);
257  if( ws_single == NULL ) {
258  std::cout << "Error: Failed to make Single-Channel workspace for channel: " << ch_name
259  << " and measurement: " << measurement.GetName() << std::endl;
260  throw hf_exc();
261  }
262 
263  // Finally, configure that workspace based on
264  // properties of the measurement
265  HistoToWorkspaceFactoryFast::ConfigureWorkspaceForMeasurement( "model_"+ch_name, ws_single, measurement );
266 
267  return ws_single;
268 
269  }
270 
272 
273  // This function takes a fully configured measurement
274  // which may contain several channels and returns
275  // a workspace holding the combined model
276  //
277  // This can be used, for example, within a script to produce
278  // a combined workspace on-the-fly
279  //
280  // This is a static function (for now) to make
281  // it a one-liner
282 
283  // First, we create an instance of a HistFactory
284  HistoToWorkspaceFactoryFast factory( measurement );
285 
286  // Loop over the channels and create the individual workspaces
287  vector<RooWorkspace*> channel_workspaces;
288  vector<string> channel_names;
289 
290  for( unsigned int chanItr = 0; chanItr < measurement.GetChannels().size(); ++chanItr ) {
291 
292  HistFactory::Channel& channel = measurement.GetChannels().at( chanItr );
293 
294  if( ! channel.CheckHistograms() ) {
295  std::cout << "MakeModelAndMeasurementsFast: Channel: " << channel.GetName()
296  << " has uninitialized histogram pointers" << std::endl;
297  throw hf_exc();
298  }
299 
300  string ch_name = channel.GetName();
301  channel_names.push_back(ch_name);
302 
303  // GHL: Renaming to 'MakeSingleChannelWorkspace'
304  RooWorkspace* ws_single = factory.MakeSingleChannelModel( measurement, channel );
305 
306  channel_workspaces.push_back(ws_single);
307 
308  }
309 
310 
311  // Now, combine the individual channel workspaces to
312  // form the combined workspace
313  RooWorkspace* ws = factory.MakeCombinedModel( channel_names, channel_workspaces );
314 
315 
316  // Configure the workspace
318 
319  // Delete channel workspaces
320  for (vector<RooWorkspace*>::iterator iter = channel_workspaces.begin() ; iter != channel_workspaces.end() ; ++iter) {
321  delete *iter ;
322  }
323 
324  // Done. Return the pointer
325  return ws;
326 
327  }
328 
330  string prefix, string productPrefix,
331  string systTerm ) {
332  if(hist) {
333  cout << "processing hist " << hist->GetName() << endl;
334  } else {
335  cout << "hist is empty" << endl;
336  R__ASSERT(hist != 0);
337  return;
338  }
339 
340  /// require dimension >=1 or <=3
341  if (fObsNameVec.empty() && !fObsName.empty()) { fObsNameVec.push_back(fObsName); }
342  R__ASSERT( fObsNameVec.size()>=1 && fObsNameVec.size()<=3 );
343 
344  /// determine histogram dimensionality
345  unsigned int histndim(1);
346  std::string classname = hist->ClassName();
347  if (classname.find("TH1")==0) { histndim=1; }
348  else if (classname.find("TH2")==0) { histndim=2; }
349  else if (classname.find("TH3")==0) { histndim=3; }
350  R__ASSERT( histndim==fObsNameVec.size() );
351 
352  /// create roorealvar observables
353  RooArgList observables;
354  std::vector<std::string>::iterator itr = fObsNameVec.begin();
355  for (int idx=0; itr!=fObsNameVec.end(); ++itr, ++idx ) {
356  if ( !proto->var(itr->c_str()) ) {
357  TAxis* axis(0);
358  if (idx==0) { axis = hist->GetXaxis(); }
359  if (idx==1) { axis = hist->GetYaxis(); }
360  if (idx==2) { axis = hist->GetZaxis(); }
361  Int_t nbins = axis->GetNbins();
362  Double_t xmin = axis->GetXmin();
363  Double_t xmax = axis->GetXmax();
364  // create observable
365  proto->factory(Form("%s[%f,%f]",itr->c_str(),xmin,xmax));
366  proto->var(itr->c_str())->setBins(nbins);
367  }
368  observables.add( *proto->var(itr->c_str()) );
369  }
370 
371  RooDataHist* histDHist = new RooDataHist((prefix+"nominalDHist").c_str(),"",observables,hist);
372  RooHistFunc* histFunc = new RooHistFunc((prefix+"_nominal").c_str(),"",observables,*histDHist,0) ;
373 
374  proto->import(*histFunc);
375 
376  /// now create the product of the overall efficiency times the sigma(params) for this estimate
377  proto->factory(("prod:"+productPrefix+"("+prefix+"_nominal,"+systTerm+")").c_str() );
378 
379  delete histDHist;
380  delete histFunc;
381 
382  }
383 
384  void HistoToWorkspaceFactoryFast::AddMultiVarGaussConstraint(RooWorkspace* proto, string prefix,int lowBin, int highBin, vector<string>& constraintTermNames){
385  // these are the nominal predictions: eg. the mean of some space of variations
386  // later fill these in a loop over histogram bins
387 
388  TVectorD mean(highBin); //-lowBin); // MB: fix range
389  cout << "a" << endl;
390  for(Int_t i=lowBin; i<highBin; ++i){
391  std::stringstream str;
392  str<<"_"<<i;
393  RooRealVar* temp = proto->var((prefix+str.str()).c_str());
394  mean(i) = temp->getVal();
395  }
396 
397  TMatrixDSym Cov(highBin-lowBin);
398  for(int i=lowBin; i<highBin; ++i){
399  for(int j=0; j<highBin-lowBin; ++j){
400  if(i==j) { Cov(i,j) = sqrt(mean(i)); } // MB : this doesn't make sense to me if lowBin!=0 (?)
401  else { Cov(i,j) = 0; }
402  }
403  }
404 
405  // can't make MultiVarGaussian with factory yet, do it by hand
406  RooArgList floating( *(proto->set(prefix.c_str() ) ) );
407  RooMultiVarGaussian constraint((prefix+"Constraint").c_str(),"",
408  floating, mean, Cov);
409 
410  proto->import(constraint);
411 
412  constraintTermNames.push_back(constraint.GetName());
413  }
414 
416  std::vector<HistoSys> histoSysList,
417  string prefix, string productPrefix,
418  string systTerm,
419  vector<string>& constraintTermNames){
420 
421  // these are the nominal predictions: eg. the mean of some space of variations
422  // later fill these in a loop over histogram bins
423 
424  // require dimension >=1 or <=3
425  if (fObsNameVec.empty() && !fObsName.empty()) { fObsNameVec.push_back(fObsName); }
426  R__ASSERT( fObsNameVec.size()>=1 && fObsNameVec.size()<=3 );
427 
428  // determine histogram dimensionality
429  unsigned int histndim(1);
430  std::string classname = nominal->ClassName();
431  if (classname.find("TH1")==0) { histndim=1; }
432  else if (classname.find("TH2")==0) { histndim=2; }
433  else if (classname.find("TH3")==0) { histndim=3; }
434  R__ASSERT( histndim==fObsNameVec.size() );
435  // cout <<"In LinInterpWithConstriants and histndim = " << histndim <<endl;
436 
437  // create roorealvar observables
438  RooArgList observables;
439  std::vector<std::string>::iterator itr = fObsNameVec.begin();
440  for (int idx=0; itr!=fObsNameVec.end(); ++itr, ++idx ) {
441  if ( !proto->var(itr->c_str()) ) {
442  TAxis* axis(NULL);
443  if (idx==0) { axis = nominal->GetXaxis(); }
444  else if (idx==1) { axis = nominal->GetYaxis(); }
445  else if (idx==2) { axis = nominal->GetZaxis(); }
446  else {
447  std::cout << "Error: Too many observables. "
448  << "HistFactory only accepts up to 3 observables (3d) "
449  << std::endl;
450  throw hf_exc();
451  }
452  Int_t nbins = axis->GetNbins();
453  Double_t xmin = axis->GetXmin();
454  Double_t xmax = axis->GetXmax();
455  // create observable
456  proto->factory(Form("%s[%f,%f]",itr->c_str(),xmin,xmax));
457  proto->var(itr->c_str())->setBins(nbins);
458  }
459  observables.add( *proto->var(itr->c_str()) );
460  }
461 
462  RooDataHist* nominalDHist = new RooDataHist((prefix+"nominalDHist").c_str(),"",observables,nominal);
463  RooHistFunc* nominalFunc = new RooHistFunc((prefix+"nominal").c_str(),"",observables,*nominalDHist,0) ;
464 
465  // make list of abstract parameters that interpolate in space of variations
466  RooArgList params( ("alpha_Hist") );
467  // range is set using defined macro (see top of the page)
468  string range=string("[")+alpha_Low+","+alpha_High+"]";
469 
470  // Loop over the HistoSys list
471  for(unsigned int j=0; j<histoSysList.size(); ++j){
472  std::stringstream str;
473  str<<"_"<<j;
474 
475  HistoSys& histoSys = histoSysList.at(j);
476  string histoSysName = histoSys.GetName();
477 
478  RooRealVar* temp = (RooRealVar*) proto->var(("alpha_" + histoSysName).c_str());
479  if(!temp){
480 
481  temp = (RooRealVar*) proto->factory(("alpha_" + histoSysName + range).c_str());
482 
483  // now add a constraint term for these parameters
484  string command=("Gaussian::alpha_"+histoSysName+"Constraint(alpha_"+histoSysName+",nom_alpha_"+histoSysName+"[0.,-10,10],1.)");
485  cout << command << endl;
486  constraintTermNames.push_back( proto->factory( command.c_str() )->GetName() );
487  proto->var(("nom_alpha_"+histoSysName).c_str())->setConstant();
488  const_cast<RooArgSet*>(proto->set("globalObservables"))->add(*proto->var(("nom_alpha_"+histoSysName).c_str()));
489  }
490  params.add(* temp );
491  }
492 
493  // now make function that linearly interpolates expectation between variations
494  // get low/high variations to interpolate between
495  vector<double> low, high;
496  RooArgSet lowSet, highSet;
497  //ES// for(unsigned int j=0; j<lowHist.size(); ++j){
498  for(unsigned int j=0; j<histoSysList.size(); ++j){
499  std::stringstream str;
500  str<<"_"<<j;
501 
502  HistoSys& histoSys = histoSysList.at(j);
503  RooDataHist* lowDHist = new RooDataHist((prefix+str.str()+"lowDHist").c_str(),"",observables, histoSys.GetHistoLow());
504  RooDataHist* highDHist = new RooDataHist((prefix+str.str()+"highDHist").c_str(),"",observables, histoSys.GetHistoHigh());
505  RooHistFunc* lowFunc = new RooHistFunc((prefix+str.str()+"low").c_str(),"",observables,*lowDHist,0) ;
506  RooHistFunc* highFunc = new RooHistFunc((prefix+str.str()+"high").c_str(),"",observables,*highDHist,0) ;
507  lowSet.add(*lowFunc);
508  highSet.add(*highFunc);
509  }
510 
511  // this is sigma(params), a piece-wise linear interpolation
512  PiecewiseInterpolation interp(prefix.c_str(),"",*nominalFunc,lowSet,highSet,params);
513  interp.setPositiveDefinite();
514  interp.setAllInterpCodes(4); // LM: change to 4 (piece-wise linear to 6th order polynomial interpolation + linear extrapolation )
515  // KC: interpo codes 1 etc. don't have proper analytic integral.
516  RooArgSet observableSet(observables);
517  interp.setBinIntegrator(observableSet);
518  interp.forceNumInt();
519 
520  proto->import(interp); // individual params have already been imported in first loop of this function
521 
522  // now create the product of the overall efficiency times the sigma(params) for this estimate
523  proto->factory(("prod:"+productPrefix+"("+prefix+","+systTerm+")").c_str() );
524 
525  }
526 
527  // GHL: Consider passing the NormFactor list instead of the entire sample
528  string HistoToWorkspaceFactoryFast::AddNormFactor(RooWorkspace* proto, string& channel, string& sigmaEpsilon, Sample& sample, bool doRatio){
529  string overallNorm_times_sigmaEpsilon ;
530  string prodNames;
531 
532  vector<NormFactor> normList = sample.GetNormFactorList();
533  vector<string> normFactorNames, rangeNames;
534 
535  if(normList.size() > 0){
536 
537  for(vector<NormFactor>::iterator itr = normList.begin(); itr != normList.end(); ++itr){
538 
539  NormFactor& norm = *itr;
540 
541  string varname;
542  if(!prodNames.empty()) prodNames += ",";
543  if(doRatio) {
544  varname = norm.GetName() + "_" + channel;
545  }
546  else {
547  varname=norm.GetName();
548  }
549 
550  // GHL: Check that the NormFactor doesn't already exist
551  // (it may have been created as a function expression
552  // during preprocessing)
553  std::stringstream range;
554  range << "[" << norm.GetVal() << "," << norm.GetLow() << "," << norm.GetHigh() << "]";
555 
556  if( proto->obj(varname.c_str()) == NULL) {
557  cout << "making normFactor: " << norm.GetName() << endl;
558  // remove "doRatio" and name can be changed when ws gets imported to the combined model.
559  proto->factory((varname + range.str()).c_str());
560  }
561 
562  if(norm.GetConst()) {
563  // proto->var(varname.c_str())->setConstant();
564  // cout <<"setting " << varname << " constant"<<endl;
565  cout << "WARNING: Const attribute to <NormFactor> tag is deprecated, will ignore." <<
566  " Instead, add \n\t<ParamSetting Const=\"True\">" << varname << "</ParamSetting>\n" <<
567  " to your top-level XML's <Measurment> entry" << endl;
568  }
569  prodNames+=varname;
570  rangeNames.push_back(range.str());
571  normFactorNames.push_back(varname);
572  }
573 
574  overallNorm_times_sigmaEpsilon = sample.GetName() + "_" + channel + "_overallNorm_x_sigma_epsilon";
575  proto->factory(("prod::" + overallNorm_times_sigmaEpsilon + "(" + prodNames + "," + sigmaEpsilon + ")").c_str());
576  }
577 
578  unsigned int rangeIndex=0;
579  for( vector<string>::iterator nit = normFactorNames.begin(); nit!=normFactorNames.end(); ++nit){
580  if( count (normFactorNames.begin(), normFactorNames.end(), *nit) > 1 ){
581  cout <<"WARNING: <NormFactor Name =\""<<*nit<<"\"> is duplicated for <Sample Name=\""
582  << sample.GetName() << "\">, but only one factor will be included. \n Instead, define something like"
583  << "\n\t<Function Name=\""<<*nit<<"Squared\" Expresion=\""<<*nit<<"*"<<*nit<<"\" Var=\""<<*nit<<rangeNames.at(rangeIndex)
584  << "\"> \nin your top-level XML's <Measurment> entry and use <NormFactor Name=\""<<*nit<<"Squared\" in your channel XML file."<< endl;
585  }
586  ++rangeIndex;
587  }
588 
589  if(!overallNorm_times_sigmaEpsilon.empty())
590  return overallNorm_times_sigmaEpsilon;
591  else
592  return sigmaEpsilon;
593  }
594 
596  string interpName,
597  std::vector<OverallSys>& systList,
598  vector<string>& constraintTermNames,
599  vector<string>& totSystTermNames) {
600 
601  // add variables for all the relative overall uncertainties we expect
602  // range is set using defined macro (see top of the page)
603 
604  string range=string("[0,")+alpha_Low+","+alpha_High+"]";
605  totSystTermNames.push_back(prefix);
606 
607  RooArgSet params(prefix.c_str());
608  vector<double> lowVec, highVec;
609 
610  std::map<std::string, double>::iterator itconstr;
611  for(unsigned int i = 0; i < systList.size(); ++i) {
612 
613  OverallSys& sys = systList.at(i);
614  const char * name = sys.GetName().c_str();
615 
616  // case of no systematic (is it possible)
617  if (meas.GetNoSyst().count(sys.GetName()) > 0 ) {
618  std::cout << "HistoToWorkspaceFast::AddConstraintTerm - skip systematic " << sys.GetName() << std::endl;
619  continue;
620  }
621  // case systematic is a gamma constraint
622  if (meas.GetGammaSyst().count(sys.GetName()) > 0 ) {
623  double relerr = meas.GetGammaSyst().find(sys.GetName() )->second;
624  if (relerr <= 0) {
625  std::cout << "HistoToWorkspaceFast::AddConstraintTerm - zero uncertainty assigned - skip systematic " << sys.GetName() << std::endl;
626  continue;
627  }
628  double tauVal = 1./(relerr*relerr);
629  double sqtau = 1./relerr;
630  RooAbsArg * beta = proto->factory(TString::Format("beta_%s[1,0,10]",name) );
631  // the global observable (y_s)
632  RooAbsArg * yvar = proto->factory(TString::Format("nom_%s[%f,0,10]",beta->GetName(),tauVal)) ;
633  // the rate of the gamma distribution (theta)
634  RooAbsArg * theta = proto->factory(TString::Format("theta_%s[%f]",name,1./tauVal));
635  // find alpha as function of beta
636  RooAbsArg* alphaOfBeta = proto->factory(TString::Format("PolyVar::alphaOfBeta_%s(beta_%s,{%f,%f})",name,name,-sqtau,sqtau));
637 
638  // add now the constraint itself Gamma_beta_constraint(beta, y+1, tau, 0 )
639  // build the gamma parameter k = as y_s + 1
640  RooAbsArg * kappa = proto->factory(TString::Format("sum::k_%s(%s,1.)",name,yvar->GetName()) );
641  RooAbsArg * gamma = proto->factory(TString::Format("Gamma::%sConstraint(%s, %s, %s, 0.0)",beta->GetName(),beta->GetName(), kappa->GetName(), theta->GetName() ) );
642  alphaOfBeta->Print("t");
643  gamma->Print("t");
644  constraintTermNames.push_back(gamma->GetName());
645  // set global observables
646  RooRealVar * gobs = dynamic_cast<RooRealVar*>(yvar); assert(gobs);
647  gobs->setConstant(true);
648  const_cast<RooArgSet*>(proto->set("globalObservables"))->add(*yvar);
649 
650  // add alphaOfBeta in the list of params to interpolate
651  params.add(*alphaOfBeta);
652  std::cout << "Added a gamma constraint for " << name << std::endl;
653 
654  }
655  else {
656 
657  // add the Gaussian constraint part
658 
659  // case systematic is uniform (asssume they are like a gauaaian bbut with a large width
660  // (100 instead of 1)
661  double gaussSigma = 1;
662  if (meas.GetUniformSyst().count(sys.GetName()) > 0 ) {
663  gaussSigma = 100;
664  std::cout << "Added a uniform constraint for " << name << " as a gaussian constraint with a very large sigma " << std::endl;
665  }
666 
667  // add Gaussian constraint terms (normal + log-normal case)
668  RooRealVar* alpha = (RooRealVar*) proto->var((prefix + sys.GetName()).c_str());
669  if(!alpha) {
670 
671  alpha = (RooRealVar*) proto->factory((prefix + sys.GetName() + range).c_str());
672  RooAbsArg * nomAlpha = proto->factory(TString::Format("nom_%s[0.,-10,10]",alpha->GetName() ) );
673  RooAbsArg * gausConstraint = proto->factory(TString::Format("Gaussian::%sConstraint(%s,%s,%f)",alpha->GetName(),alpha->GetName(), nomAlpha->GetName(), gaussSigma) );
674  //cout << command << endl;
675  constraintTermNames.push_back( gausConstraint->GetName() );
676  proto->var(("nom_" + prefix + sys.GetName()).c_str())->setConstant();
677  const_cast<RooArgSet*>(proto->set("globalObservables"))->add(*nomAlpha);
678  }
679 
680 
681  // add constraint in terms of bifrucated gauss with low/high as sigmas
682  //std::stringstream lowhigh;
683 
684  // check if exists a log-normal constraint
685  if (meas.GetLogNormSyst().count(sys.GetName()) == 0 && meas.GetGammaSyst().count(sys.GetName()) == 0 ) {
686  // just add the alpha for the parameters of the FlexibleInterpVar function
687  params.add(*alpha);
688  }
689  // case systematic is a log-normal constraint
690  if (meas.GetLogNormSyst().count(sys.GetName()) > 0 ) {
691  // log normal constraint for parameter
692  double relerr = meas.GetLogNormSyst().find(sys.GetName() )->second;
693  double tauVal = 1./relerr;
694  std::string tauName = "tau_" + sys.GetName();
695  proto->factory(TString::Format("%s[%f]",tauName.c_str(),tauVal ) );
696  double kappaVal = 1. + relerr;
697  std::string kappaName = "kappa_" + sys.GetName();
698  proto->factory(TString::Format("%s[%f]",kappaName.c_str(),kappaVal ) );
699  const char * alphaName = alpha->GetName();
700 
701  std::string alphaOfBetaName = "alphaOfBeta_" + sys.GetName();
702  RooAbsArg * alphaOfBeta = proto->factory(TString::Format("expr::%s('%s*(pow(%s,%s)-1.)',%s,%s,%s)",alphaOfBetaName.c_str(),
703  tauName.c_str(),kappaName.c_str(),alphaName,
704  tauName.c_str(),kappaName.c_str(),alphaName ) );
705 
706  std::cout << "Added a log-normal constraint for " << name << std::endl;
707  alphaOfBeta->Print("t");
708  params.add(*alphaOfBeta);
709  }
710 
711  }
712  // add low/high vectors
713  double low = sys.GetLow();
714  double high = sys.GetHigh();
715  lowVec.push_back(low);
716  highVec.push_back(high);
717 
718  } // end sys loop
719 
720  if(systList.size() > 0){
721  // this is epsilon(alpha_j), a piece-wise linear interpolation
722  // LinInterpVar interp( (interpName).c_str(), "", params, 1., lowVec, highVec);
723 
724  assert( params.getSize() > 0);
725  assert(int(lowVec.size()) == params.getSize() );
726 
727  FlexibleInterpVar interp( (interpName).c_str(), "", params, 1., lowVec, highVec);
728  interp.setAllInterpCodes(4); // LM: change to 4 (piece-wise linear to 6th order polynomial interpolation + linear extrapolation )
729  //interp.setAllInterpCodes(0); // simple linear interpolation
730  proto->import(interp); // params have already been imported in first loop of this function
731  } else{
732  // some strange behavior if params,lowVec,highVec are empty.
733  //cout << "WARNING: No OverallSyst terms" << endl;
734  RooConstVar interp( (interpName).c_str(), "", 1.);
735  proto->import(interp); // params have already been imported in first loop of this function
736  }
737 
738  // std::cout << "after creating FlexibleInterpVar " << std::endl;
739  // proto->Print();
740 
741  }
742 
743 
745  vector<string>& syst_x_expectedPrefixNames,
746  vector<string>& normByNames){
747 
748  // for ith bin calculate totN_i = lumi * sum_j expected_j * syst_j
749  string command;
750  string coeffList="";
751  string shapeList="";
752  string prepend="";
753 
754  if (fObsNameVec.empty() && !fObsName.empty()) { fObsNameVec.push_back(fObsName); }
755 
756  double binWidth(1.0);
757  std::string obsNameVecStr;
758  std::vector<std::string>::iterator itr = fObsNameVec.begin();
759  for (; itr!=fObsNameVec.end(); ++itr) {
760  std::string obsName = *itr;
761  binWidth *= proto->var(obsName.c_str())->numBins()/(proto->var(obsName.c_str())->getMax() - proto->var(obsName.c_str())->getMin()) ; // MB: Note: requires fixed bin sizes
762  if (obsNameVecStr.size()>0) { obsNameVecStr += "_"; }
763  obsNameVecStr += obsName;
764  }
765 
766  //vector<string>::iterator it=syst_x_expectedPrefixNames.begin();
767  for(unsigned int j=0; j<syst_x_expectedPrefixNames.size();++j){
768  std::stringstream str;
769  str<<"_"<<j;
770  // repatative, but we need one coeff for each term in the sum
771  // maybe can be avoided if we don't use bin width as coefficient
772  command=string(Form("binWidth_%s_%d[%e]",obsNameVecStr.c_str(),j,binWidth));
773  proto->factory(command.c_str());
774  proto->var(Form("binWidth_%s_%d",obsNameVecStr.c_str(),j))->setConstant();
775  coeffList+=prepend+"binWidth_"+obsNameVecStr+str.str();
776 
777  command="prod::L_x_"+syst_x_expectedPrefixNames.at(j)+"("+normByNames.at(j)+","+syst_x_expectedPrefixNames.at(j)+")";
778  /*RooAbsReal* tempFunc =(RooAbsReal*) */
779  proto->factory(command.c_str());
780  shapeList+=prepend+"L_x_"+syst_x_expectedPrefixNames.at(j);
781  prepend=",";
782 
783  // add to num int to product
784  // tempFunc->specialIntegratorConfig(kTRUE)->method1D().setLabel("RooBinIntegrator") ;
785  // tempFunc->forceNumInt();
786 
787  }
788 
789  proto->defineSet("coefList",coeffList.c_str());
790  proto->defineSet("shapeList",shapeList.c_str());
791  // proto->factory(command.c_str());
792  RooRealSumPdf tot(totName.c_str(),totName.c_str(),*proto->set("shapeList"),*proto->set("coefList"),kTRUE);
793  tot.specialIntegratorConfig(kTRUE)->method1D().setLabel("RooBinIntegrator") ;
794  tot.specialIntegratorConfig(kTRUE)->method2D().setLabel("RooBinIntegrator") ;
795  tot.specialIntegratorConfig(kTRUE)->methodND().setLabel("RooBinIntegrator") ;
796  tot.forceNumInt();
797 
798  // for mixed generation in RooSimultaneous
799  tot.setAttribute("GenerateBinned"); // for use with RooSimultaneous::generate in mixed mode
800  // tot.setAttribute("GenerateUnbinned"); // we don't want that
801 
802  /*
803  // Use binned numeric integration
804  int nbins = 0;
805  if( fObsNameVec.size() == 1 ) {
806  nbins = proto->var(fObsNameVec.at(0).c_str())->numBins();
807 
808  cout <<"num bis for RooRealSumPdf = "<<nbins <<endl;
809  //int nbins = ((RooRealVar*) allVars.first())->numBins();
810  tot.specialIntegratorConfig(kTRUE)->getConfigSection("RooBinIntegrator").setRealValue("numBins",nbins);
811  tot.forceNumInt();
812 
813  } else {
814  cout << "Bin Integrator only supports 1-d. Will be slow." << std::endl;
815  }
816  */
817 
818 
819  proto->import(tot);
820 
821  }
822 
823  void HistoToWorkspaceFactoryFast::AddPoissonTerms(RooWorkspace* proto, string prefix, string obsPrefix, string expPrefix, int lowBin, int highBin,
824  vector<string>& likelihoodTermNames){
825  /////////////////////////////////
826  // Relate observables to expected for each bin
827  // later modify variable named expPrefix_i to be product of terms
828  RooArgSet Pois(prefix.c_str());
829  for(Int_t i=lowBin; i<highBin; ++i){
830  std::stringstream str;
831  str<<"_"<<i;
832  //string command("Poisson::"+prefix+str.str()+"("+obsPrefix+str.str()+","+expPrefix+str.str()+")");
833  string command("Poisson::"+prefix+str.str()+"("+obsPrefix+str.str()+","+expPrefix+str.str()+",1)");//for no rounding
834  RooAbsArg* temp = (proto->factory( command.c_str() ) );
835 
836  // output
837  cout << "Poisson Term " << command << endl;
838  ((RooAbsPdf*) temp)->setEvalErrorLoggingMode(RooAbsReal::PrintErrors);
839  //cout << temp << endl;
840 
841  likelihoodTermNames.push_back( temp->GetName() );
842  Pois.add(* temp );
843  }
844  proto->defineSet(prefix.c_str(),Pois); // add argset to workspace
845  }
846 
847  void HistoToWorkspaceFactoryFast::SetObsToExpected(RooWorkspace* proto, string obsPrefix, string expPrefix, int lowBin, int highBin){
848  /////////////////////////////////
849  // set observed to expected
850  TTree* tree = new TTree();
851  Double_t* obsForTree = new Double_t[highBin-lowBin];
852  RooArgList obsList("obsList");
853 
854  for(Int_t i=lowBin; i<highBin; ++i){
855  std::stringstream str;
856  str<<"_"<<i;
857  RooRealVar* obs = (RooRealVar*) proto->var((obsPrefix+str.str()).c_str());
858  cout << "expected number of events called: " << expPrefix << endl;
859  RooAbsReal* exp = proto->function((expPrefix+str.str()).c_str());
860  if(obs && exp){
861 
862  //proto->Print();
863  obs->setVal( exp->getVal() );
864  cout << "setting obs"+str.str()+" to expected = " << exp->getVal() << " check: " << obs->getVal() << endl;
865 
866  // add entry to array and attach to tree
867  obsForTree[i] = exp->getVal();
868  tree->Branch((obsPrefix+str.str()).c_str(), obsForTree+i ,(obsPrefix+str.str()+"/D").c_str());
869  obsList.add(*obs);
870  }else{
871  cout << "problem retrieving obs or exp " << obsPrefix+str.str() << obs << " " << expPrefix+str.str() << exp << endl;
872  }
873  }
874  tree->Fill();
875  RooDataSet* data = new RooDataSet("expData","", tree, obsList); // one experiment
876 
877  delete tree;
878  delete [] obsForTree;
879 
880  proto->import(*data);
881 
882  delete data;
883 
884  }
885 
886  //////////////////////////////////////////////////////////////////////////////
887 
889  map<string,double> gammaSyst,
890  map<string,double> uniformSyst,
891  map<string,double> logNormSyst,
892  map<string,double> noSyst) {
893  string pdfName(pdfNameChar);
894 
895  ModelConfig * combined_config = (ModelConfig *) proto->obj("ModelConfig");
896  if( combined_config==NULL ) {
897  std::cout << "Error: Failed to find object 'ModelConfig' in workspace: "
898  << proto->GetName() << std::endl;
899  throw hf_exc();
900  }
901  // const RooArgSet * constrainedParams=combined_config->GetNuisanceParameters();
902  // RooArgSet temp(*constrainedParams);
903  string edit="EDIT::newSimPdf("+pdfName+",";
904  string editList;
905  string lastPdf=pdfName;
906  string precede="";
907  unsigned int numReplacements = 0;
908  unsigned int nskipped = 0;
909  map<string,double>::iterator it;
910 
911 
912  // add gamma terms and their constraints
913  for(it=gammaSyst.begin(); it!=gammaSyst.end(); ++it) {
914  //cout << "edit for " << it->first << "with rel uncert = " << it->second << endl;
915  if(! proto->var(("alpha_"+it->first).c_str())){
916  //cout << "systematic not there" << endl;
917  nskipped++;
918  continue;
919  }
920  numReplacements++;
921 
922  double relativeUncertainty = it->second;
923  double scale = 1/sqrt((1+1/pow(relativeUncertainty,2)));
924 
925  // this is the Gamma PDF and in a form that doesn't have roundoff problems like the Poisson does
926  proto->factory(Form("beta_%s[1,0,10]",it->first.c_str()));
927  proto->factory(Form("y_%s[%f]",it->first.c_str(),1./pow(relativeUncertainty,2))) ;
928  proto->factory(Form("theta_%s[%f]",it->first.c_str(),pow(relativeUncertainty,2))) ;
929  proto->factory(Form("Gamma::beta_%sConstraint(beta_%s,sum::k_%s(y_%s,one[1]),theta_%s,zero[0])",
930  it->first.c_str(),
931  it->first.c_str(),
932  it->first.c_str(),
933  it->first.c_str(),
934  it->first.c_str())) ;
935 
936  /*
937  // this has some problems because N in poisson is rounded to nearest integer
938  proto->factory(Form("Poisson::beta_%sConstraint(y_%s[%f],prod::taub_%s(taus_%s[%f],beta_%s[1,0,5]))",
939  it->first.c_str(),
940  it->first.c_str(),
941  1./pow(relativeUncertainty,2),
942  it->first.c_str(),
943  it->first.c_str(),
944  1./pow(relativeUncertainty,2),
945  it->first.c_str()
946  ) ) ;
947  */
948  // combined->factory(Form("expr::alphaOfBeta('(beta-1)/%f',beta)",scale));
949  // 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()));
950  proto->factory(Form("PolyVar::alphaOfBeta_%s(beta_%s,{%f,%f})",it->first.c_str(),it->first.c_str(),-1./scale,1./scale));
951 
952  // set beta const status to be same as alpha
953  if(proto->var(Form("alpha_%s",it->first.c_str()))->isConstant()) {
954  proto->var(Form("beta_%s",it->first.c_str()))->setConstant(true);
955  }
956  else {
957  proto->var(Form("beta_%s",it->first.c_str()))->setConstant(false);
958  }
959  // set alpha const status to true
960  // proto->var(Form("alpha_%s",it->first.c_str()))->setConstant(true);
961 
962  // replace alphas with alphaOfBeta and replace constraints
963  editList+=precede + "alpha_"+it->first+"Constraint=beta_" + it->first+ "Constraint";
964  precede=",";
965  editList+=precede + "alpha_"+it->first+"=alphaOfBeta_"+ it->first;
966 
967  /*
968  if( proto->pdf(("alpha_"+it->first+"Constraint").c_str()) && proto->var(("alpha_"+it->first).c_str()) )
969  cout << " checked they are there" << proto->pdf(("alpha_"+it->first+"Constraint").c_str()) << " " << proto->var(("alpha_"+it->first).c_str()) << endl;
970  else
971  cout << "NOT THERE" << endl;
972  */
973 
974  // EDIT seems to die if the list of edits is too long. So chunck them up.
975  if(numReplacements%10 == 0 && numReplacements+nskipped!=gammaSyst.size()){
976  edit="EDIT::"+lastPdf+"_("+lastPdf+","+editList+")";
977  lastPdf+="_"; // append an underscore for the edit
978  editList=""; // reset edit list
979  precede="";
980  cout << "Going to issue this edit command\n" << edit<< endl;
981  proto->factory( edit.c_str() );
982  RooAbsPdf* newOne = proto->pdf(lastPdf.c_str());
983  if(!newOne)
984  cout << "\n\n ---------------------\n WARNING: failed to make EDIT\n\n" << endl;
985 
986  }
987  }
988 
989  // add uniform terms and their constraints
990  for(it=uniformSyst.begin(); it!=uniformSyst.end(); ++it) {
991  cout << "edit for " << it->first << "with rel uncert = " << it->second << endl;
992  if(! proto->var(("alpha_"+it->first).c_str())){
993  cout << "systematic not there" << endl;
994  nskipped++;
995  continue;
996  }
997  numReplacements++;
998 
999  // this is the Uniform PDF
1000  proto->factory(Form("beta_%s[1,0,10]",it->first.c_str()));
1001  proto->factory(Form("Uniform::beta_%sConstraint(beta_%s)",it->first.c_str(),it->first.c_str()));
1002  proto->factory(Form("PolyVar::alphaOfBeta_%s(beta_%s,{-1,1})",it->first.c_str(),it->first.c_str()));
1003 
1004  // set beta const status to be same as alpha
1005  if(proto->var(Form("alpha_%s",it->first.c_str()))->isConstant())
1006  proto->var(Form("beta_%s",it->first.c_str()))->setConstant(true);
1007  else
1008  proto->var(Form("beta_%s",it->first.c_str()))->setConstant(false);
1009  // set alpha const status to true
1010  // proto->var(Form("alpha_%s",it->first.c_str()))->setConstant(true);
1011 
1012  // replace alphas with alphaOfBeta and replace constraints
1013  cout << "alpha_"+it->first+"Constraint=beta_" + it->first+ "Constraint" << endl;
1014  editList+=precede + "alpha_"+it->first+"Constraint=beta_" + it->first+ "Constraint";
1015  precede=",";
1016  cout << "alpha_"+it->first+"=alphaOfBeta_"+ it->first << endl;
1017  editList+=precede + "alpha_"+it->first+"=alphaOfBeta_"+ it->first;
1018 
1019  if( proto->pdf(("alpha_"+it->first+"Constraint").c_str()) && proto->var(("alpha_"+it->first).c_str()) )
1020  cout << " checked they are there" << proto->pdf(("alpha_"+it->first+"Constraint").c_str()) << " " << proto->var(("alpha_"+it->first).c_str()) << endl;
1021  else
1022  cout << "NOT THERE" << endl;
1023 
1024  // EDIT seems to die if the list of edits is too long. So chunck them up.
1025  if(numReplacements%10 == 0 && numReplacements+nskipped!=gammaSyst.size()){
1026  edit="EDIT::"+lastPdf+"_("+lastPdf+","+editList+")";
1027  lastPdf+="_"; // append an underscore for the edit
1028  editList=""; // reset edit list
1029  precede="";
1030  cout << edit<< endl;
1031  proto->factory( edit.c_str() );
1032  RooAbsPdf* newOne = proto->pdf(lastPdf.c_str());
1033  if(!newOne)
1034  cout << "\n\n ---------------------\n WARNING: failed to make EDIT\n\n" << endl;
1035 
1036  }
1037  }
1038 
1039  /////////////////////////////////////////
1040  ////////////////////////////////////
1041 
1042 
1043  // add lognormal terms and their constraints
1044  for(it=logNormSyst.begin(); it!=logNormSyst.end(); ++it) {
1045  cout << "edit for " << it->first << "with rel uncert = " << it->second << endl;
1046  if(! proto->var(("alpha_"+it->first).c_str())){
1047  cout << "systematic not there" << endl;
1048  nskipped++;
1049  continue;
1050  }
1051  numReplacements++;
1052 
1053  double relativeUncertainty = it->second;
1054  double kappa = 1+relativeUncertainty;
1055  // when transforming beta -> alpha, need alpha=1 to be +1sigma value.
1056  // the P(beta>kappa*\hat(beta)) = 16%
1057  // and \hat(beta) is 1, thus
1058  double scale = relativeUncertainty;
1059  //double scale = kappa;
1060 
1061  const char * cname = it->first.c_str();
1062 
1063  // this is the LogNormal
1064  proto->factory(TString::Format("beta_%s[1,0,10]",cname));
1065  proto->factory(TString::Format("nom_beta_%s[1]",cname));
1066  proto->factory(TString::Format("kappa_%s[%f]",cname,kappa));
1067  proto->factory(TString::Format("Lognormal::beta_%sConstraint(beta_%s,nom_beta_%s,kappa_%s)",
1068  cname, cname, cname, cname)) ;
1069  proto->factory(TString::Format("PolyVar::alphaOfBeta_%s(beta_%s,{%f,%f})",cname,cname,-1./scale,1./scale));
1070 
1071 
1072  // set beta const status to be same as alpha
1073  if(proto->var(TString::Format("alpha_%s",cname))->isConstant())
1074  proto->var(TString::Format("beta_%s",cname))->setConstant(true);
1075  else
1076  proto->var(TString::Format("beta_%s",cname))->setConstant(false);
1077  // set alpha const status to true
1078  // proto->var(TString::Format("alpha_%s",cname))->setConstant(true);
1079 
1080  // replace alphas with alphaOfBeta and replace constraints
1081  cout << "alpha_"+it->first+"Constraint=beta_" + it->first+ "Constraint" << endl;
1082  editList+=precede + "alpha_"+it->first+"Constraint=beta_" + it->first+ "Constraint";
1083  precede=",";
1084  cout << "alpha_"+it->first+"=alphaOfBeta_"+ it->first << endl;
1085  editList+=precede + "alpha_"+it->first+"=alphaOfBeta_"+ it->first;
1086 
1087  if( proto->pdf(("alpha_"+it->first+"Constraint").c_str()) && proto->var(("alpha_"+it->first).c_str()) )
1088  cout << " checked they are there" << proto->pdf(("alpha_"+it->first+"Constraint").c_str()) << " " << proto->var(("alpha_"+it->first).c_str()) << endl;
1089  else
1090  cout << "NOT THERE" << endl;
1091 
1092  // EDIT seems to die if the list of edits is too long. So chunck them up.
1093  if(numReplacements%10 == 0 && numReplacements+nskipped!=gammaSyst.size()){
1094  edit="EDIT::"+lastPdf+"_("+lastPdf+","+editList+")";
1095  lastPdf+="_"; // append an underscore for the edit
1096  editList=""; // reset edit list
1097  precede="";
1098  cout << edit<< endl;
1099  proto->factory( edit.c_str() );
1100  RooAbsPdf* newOne = proto->pdf(lastPdf.c_str());
1101  if(!newOne)
1102  cout << "\n\n ---------------------\n WARNING: failed to make EDIT\n\n" << endl;
1103 
1104  }
1105  // add global observables
1106  const RooArgSet * gobs = proto->set("globalObservables");
1107  RooArgSet gobsNew(*gobs);
1108  gobsNew.add(*proto->var(TString::Format("nom_beta_%s",cname)) );
1109  proto->removeSet("globalObservables");
1110  proto->defineSet("globalObservables",gobsNew);
1111  gobsNew.Print();
1112 
1113  }
1114 
1115  /////////////////////////////////////////
1116 
1117  // MB: remove a systematic constraint
1118  for(it=noSyst.begin(); it!=noSyst.end(); ++it) {
1119 
1120  cout << "remove constraint for parameter" << it->first << endl;
1121  if(! proto->var(("alpha_"+it->first).c_str()) || ! proto->pdf(("alpha_"+it->first+"Constraint").c_str()) ) {
1122  cout << "systematic not there" << endl;
1123  nskipped++;
1124  continue;
1125  }
1126  numReplacements++;
1127 
1128  // dummy replacement pdf
1129  if ( !proto->var("one") ) { proto->factory("one[1.0]"); }
1130  proto->var("one")->setConstant();
1131 
1132  // replace constraints
1133  cout << "alpha_"+it->first+"Constraint=one" << endl;
1134  editList+=precede + "alpha_"+it->first+"Constraint=one";
1135  precede=",";
1136 
1137  // EDIT seems to die if the list of edits is too long. So chunck them up.
1138  if(numReplacements%10 == 0 && numReplacements+nskipped!=gammaSyst.size()){
1139  edit="EDIT::"+lastPdf+"_("+lastPdf+","+editList+")";
1140  lastPdf+="_"; // append an underscore for the edit
1141  editList=""; // reset edit list
1142  precede="";
1143  cout << edit << endl;
1144  proto->factory( edit.c_str() );
1145  RooAbsPdf* newOne = proto->pdf(lastPdf.c_str());
1146  if(!newOne) { cout << "\n\n ---------------------\n WARNING: failed to make EDIT\n\n" << endl; }
1147  }
1148  }
1149 
1150  /////////////////////////////////////////
1151 
1152  // commit last bunch of edits
1153  edit="EDIT::newSimPdf("+lastPdf+","+editList+")";
1154  cout << edit<< endl;
1155  proto->factory( edit.c_str() );
1156  // proto->writeToFile(("results/model_"+fRowTitle+"_edited.root").c_str());
1157  RooAbsPdf* newOne = proto->pdf("newSimPdf");
1158  if(newOne){
1159  // newOne->graphVizTree(("results/"+pdfName+"_"+fRowTitle+"newSimPdf.dot").c_str());
1160  combined_config->SetPdf(*newOne);
1161  }
1162  else{
1163  cout << "\n\n ---------------------\n WARNING: failed to make EDIT\n\n" << endl;
1164  }
1165  }
1166 
1168  // Change-> Now a static utility
1169 
1170  FILE* covFile = fopen ((filename).c_str(),"w");
1171 
1172  TIter iti = params->createIterator();
1173  TIter itj = params->createIterator();
1174  RooRealVar *myargi, *myargj;
1175  fprintf(covFile," ") ;
1176  while ((myargi = (RooRealVar *)iti.Next())) {
1177  if(myargi->isConstant()) continue;
1178  fprintf(covFile," & %s", myargi->GetName());
1179  }
1180  fprintf(covFile,"\\\\ \\hline \n" );
1181  iti.Reset();
1182  while ((myargi = (RooRealVar *)iti.Next())) {
1183  if(myargi->isConstant()) continue;
1184  fprintf(covFile,"%s", myargi->GetName());
1185  itj.Reset();
1186  while ((myargj = (RooRealVar *)itj.Next())) {
1187  if(myargj->isConstant()) continue;
1188  cout << myargi->GetName() << "," << myargj->GetName();
1189  fprintf(covFile, " & %.2f", result->correlation(*myargi, *myargj));
1190  }
1191  cout << endl;
1192  fprintf(covFile, " \\\\\n");
1193  }
1194  fclose(covFile);
1195 
1196  }
1197 
1198 
1199  ///////////////////////////////////////////////
1201 
1202  // check inputs (see JIRA-6890 )
1203 
1204  if (channel.GetSamples().empty()) {
1205  Error("MakeSingleChannelWorkspace","The input Channel does not contain any sample - return a nullptr");
1206  return 0;
1207  }
1208 
1209  // Set these by hand inside the function
1210  vector<string> systToFix = measurement.GetConstantParams();
1211  bool doRatio=false;
1212 
1213  // to time the macro
1214  TStopwatch t;
1215  t.Start();
1216  //ES// string channel_name=summary[0].channel;
1217  string channel_name = channel.GetName();
1218 
1219  /// MB: reset observable names for each new channel.
1220  fObsNameVec.clear();
1221 
1222  /// MB: label observables x,y,z, depending on histogram dimensionality
1223  /// GHL: Give it the first sample's nominal histogram as a template
1224  /// since the data histogram may not be present
1225  TH1* channel_hist_template = channel.GetSamples().at(0).GetHisto();
1226  if (fObsNameVec.empty()) { GuessObsNameVec(channel_hist_template); }
1227 
1228  for ( unsigned int idx=0; idx<fObsNameVec.size(); ++idx ) {
1229  fObsNameVec[idx] = "obs_" + fObsNameVec[idx] + "_" + channel_name ;
1230  }
1231 
1232  if (fObsNameVec.empty()) {
1233  fObsName= "obs_" + channel_name; // set name ov observable
1234  fObsNameVec.push_back( fObsName );
1235  }
1236 
1237  R__ASSERT( fObsNameVec.size()>=1 && fObsNameVec.size()<=3 );
1238 
1239  cout << "\n\n-------------------\nStarting to process " << channel_name << " channel with " << fObsNameVec.size() << " observables" << endl;
1240 
1241  //
1242  // our main workspace that we are using to construct the model
1243  //
1244  RooWorkspace* proto = new RooWorkspace(channel_name.c_str(), (channel_name+" workspace").c_str());
1245  ModelConfig * proto_config = new ModelConfig("ModelConfig", proto);
1246  proto_config->SetWorkspace(*proto);
1247 
1248  // preprocess functions
1249  vector<string>::iterator funcIter = fPreprocessFunctions.begin();
1250  for(;funcIter!= fPreprocessFunctions.end(); ++funcIter){
1251  cout <<"will preprocess this line: " << *funcIter <<endl;
1252  proto->factory(funcIter->c_str());
1253  proto->Print();
1254  }
1255 
1256  RooArgSet likelihoodTerms("likelihoodTerms"), constraintTerms("constraintTerms");
1257  vector<string> likelihoodTermNames, constraintTermNames, totSystTermNames, syst_x_expectedPrefixNames, normalizationNames;
1258 
1259  vector< pair<string,string> > statNamePairs;
1260  vector< pair<TH1*,TH1*> > statHistPairs; // <nominal, error>
1261  std::string statFuncName; // the name of the ParamHistFunc
1262  std::string statNodeName; // the name of the McStat Node
1263  // Constraint::Type statConstraintType=Constraint::Gaussian;
1264  // Double_t statRelErrorThreshold=0.0;
1265 
1266  string prefix, range;
1267 
1268  /////////////////////////////
1269  // shared parameters
1270  // this is ratio of lumi to nominal lumi. We will include relative uncertainty in model
1271  std::stringstream lumiStr;
1272  // lumi range
1273  lumiStr<<"["<<fNomLumi<<",0,"<<10.*fNomLumi<<"]";
1274  proto->factory(("Lumi"+lumiStr.str()).c_str());
1275  cout << "lumi str = " << lumiStr.str() << endl;
1276 
1277  std::stringstream lumiErrorStr;
1278  lumiErrorStr << "nominalLumi["<<fNomLumi << ",0,"<<fNomLumi+10*fLumiError<<"]," << fLumiError ;
1279  proto->factory(("Gaussian::lumiConstraint(Lumi,"+lumiErrorStr.str()+")").c_str());
1280  proto->var("nominalLumi")->setConstant();
1281  proto->defineSet("globalObservables","nominalLumi");
1282  //likelihoodTermNames.push_back("lumiConstraint");
1283  constraintTermNames.push_back("lumiConstraint");
1284  cout << "lumi Error str = " << lumiErrorStr.str() << endl;
1285 
1286  //proto->factory((string("SigXsecOverSM[1.,0.5,1..8]").c_str()));
1287  ///////////////////////////////////
1288  // loop through estimates, add expectation, floating bin predictions,
1289  // and terms that constrain floating to expectation via uncertainties
1290  // GHL: Loop over samples instead, which doesn't contain the data
1291  vector<Sample>::iterator it = channel.GetSamples().begin();
1292  for(; it!=channel.GetSamples().end(); ++it) {
1293 
1294  //ES// string overallSystName = it->name+"_"+it->channel+"_epsilon";
1295  Sample& sample = (*it);
1296  string overallSystName = sample.GetName() + "_" + channel_name + "_epsilon";
1297 
1298  string systSourcePrefix = "alpha_";
1299 
1300  // constraintTermNames and totSystTermNames are vectors that are passed
1301  // by reference and filled by this method
1302  AddConstraintTerms(proto,measurement, systSourcePrefix, overallSystName,
1303  sample.GetOverallSysList(), constraintTermNames , totSystTermNames);
1304 
1305  // GHL: Consider passing the NormFactor list instead of the entire sample
1306  overallSystName = AddNormFactor(proto, channel_name, overallSystName, sample, doRatio);
1307 
1308  // Create the string for the object
1309  // that is added to the RooRealSumPdf
1310  // for this channel
1311  string syst_x_expectedPrefix = "";
1312 
1313  // get histogram
1314  //ES// TH1* nominal = it->nominal;
1315  TH1* nominal = sample.GetHisto();
1316 
1317  // MB : HACK no option to have both non-hist variations and hist variations ?
1318  // get histogram
1319  // GHL: Okay, this is going to be non-trivial.
1320  // We will loop over histosys's, which contain both
1321  // the low hist and the high hist together.
1322 
1323  // Logic:
1324  // - If we have no HistoSys's, do part A
1325  // - else, if the histo syst's don't match, return (we ignore this case)
1326  // - finally, we take the syst's and apply the linear interpolation w/ constraint
1327 
1328  if(sample.GetHistoSysList().size() == 0) {
1329 
1330  // If no HistoSys
1331  cout << sample.GetName() + "_" + channel_name + " has no variation histograms " << endl;
1332  string expPrefix = sample.GetName() + "_" + channel_name; //+"_expN";
1333  syst_x_expectedPrefix = sample.GetName() + "_" + channel_name + "_overallSyst_x_Exp";
1334 
1335  ProcessExpectedHisto(sample.GetHisto(), proto, expPrefix, syst_x_expectedPrefix,
1336  overallSystName);
1337  }
1338  else {
1339  // If there ARE HistoSys(s)
1340  // name of source for variation
1341  string constraintPrefix = sample.GetName() + "_" + channel_name + "_Hist_alpha";
1342  syst_x_expectedPrefix = sample.GetName() + "_" + channel_name + "_overallSyst_x_HistSyst";
1343  // constraintTermNames are passed by reference and appended to,
1344  // overallSystName is a std::string for this sample
1345 
1346  LinInterpWithConstraint(proto, nominal, sample.GetHistoSysList(),
1347  constraintPrefix, syst_x_expectedPrefix, overallSystName,
1348  constraintTermNames);
1349  }
1350 
1351  ////////////////////////////////////
1352  // Add StatErrors to this Channel //
1353  ////////////////////////////////////
1354 
1355  if( sample.GetStatError().GetActivate() ) {
1356 
1357  if( fObsNameVec.size() > 3 ) {
1358  std::cout << "Cannot include Stat Error for histograms of more than 3 dimensions."
1359  << std::endl;
1360  throw hf_exc();
1361  } else {
1362 
1363  // If we are using StatUncertainties, we multiply this object
1364  // by the ParamHistFunc and then pass that to the
1365  // RooRealSumPdf by appending it's name to the list
1366 
1367  std::cout << "Sample: " << sample.GetName() << " to be included in Stat Error "
1368  << "for channel " << channel_name
1369  << std::endl;
1370 
1371  /*
1372  Constraint::Type type = channel.GetStatErrorConfig().GetConstraintType();
1373  statConstraintType = Constraint::Gaussian;
1374  if( type == Constraint::Gaussian) {
1375  std::cout << "Using Gaussian StatErrors" << std::endl;
1376  statConstraintType = Constraint::Gaussian;
1377  }
1378  if( type == Constraint::Poisson ) {
1379  std::cout << "Using Poisson StatErrors" << std::endl;
1380  statConstraintType = Constraint::Poisson;
1381  }
1382  */
1383 
1384  //statRelErrorThreshold = channel.GetStatErrorConfig().GetRelErrorThreshold();
1385 
1386  // First, get the uncertainty histogram
1387  // and push it back to our vectors
1388 
1389  //if( sample.GetStatError().GetErrorHist() ) {
1390  //statErrorHist = (TH1*) sample.GetStatError().GetErrorHist()->Clone();
1391  //}
1392  //if( statErrorHist == NULL ) {
1393 
1394  // We need to get the *ABSOLUTE* uncertainty for use in Stat Uncertainties
1395  // This can be done in one of two ways:
1396  // - Use the built-in Errors in the TH1 itself (they are aboslute)
1397  // - Take the supplied *RELATIVE* error and multiply by the nominal
1398  string UncertName = syst_x_expectedPrefix + "_StatAbsolUncert";
1399  TH1* statErrorHist = NULL;
1400 
1401  if( sample.GetStatError().GetErrorHist() == NULL ) {
1402  // Make the absolute stat error
1403  std::cout << "Making Statistical Uncertainty Hist for "
1404  << " Channel: " << channel_name
1405  << " Sample: " << sample.GetName()
1406  << std::endl;
1407  statErrorHist = MakeAbsolUncertaintyHist( UncertName, nominal );
1408  } else {
1409  // clone the error histograms because in case the sample has not error hist
1410  // it is created in MakeAbsolUncertainty
1411  // we need later to clean statErrorHist
1412  statErrorHist = (TH1*) sample.GetStatError().GetErrorHist()->Clone();
1413  // We assume the (relative) error is provided.
1414  // We must turn it into an absolute error
1415  // using the nominal histogram
1416  std::cout << "Using external histogram for Stat Errors for "
1417  << " Channel: " << channel_name
1418  << " Sample: " << sample.GetName()
1419  << std::endl;
1420  std::cout << "Error Histogram: " << statErrorHist->GetName() << std::endl;
1421  // Multiply the relative stat uncertainty by the
1422  // nominal to get the overall stat uncertainty
1423  statErrorHist->Multiply( nominal );
1424  statErrorHist->SetName( UncertName.c_str() );
1425  }
1426 
1427  // Save the nominal and error hists
1428  // for the building of constraint terms
1429  statHistPairs.push_back( pair<TH1*,TH1*>(nominal, statErrorHist) );
1430 
1431  // To do the 'conservative' version, we would need to do some
1432  // intervention here. We would probably need to create a different
1433  // ParamHistFunc for each sample in the channel. The would nominally
1434  // use the same gamma's, so we haven't increased the number of parameters
1435  // However, if a bin in the 'nominal' histogram is 0, we simply need to
1436  // change the parameter in that bin in the ParamHistFunc for this sample.
1437  // We also need to add a constraint term.
1438  // Actually, we'd probably not use the ParamHistFunc...?
1439  // we could remove the dependence in this ParamHistFunc on the ith gamma
1440  // and then create the poisson term: Pois(tau | n_exp)Pois(data | n_exp)
1441 
1442 
1443  // Next, try to get the ParamHistFunc (it may have been
1444  // created by another sample in this channel)
1445  // or create it if it doesn't yet exist:
1446  statFuncName = "mc_stat_" + channel_name;
1447  ParamHistFunc* paramHist = (ParamHistFunc*) proto->function( statFuncName.c_str() );
1448  if( paramHist == NULL ) {
1449 
1450  // Get a RooArgSet of the observables:
1451  // Names in the list fObsNameVec:
1452  RooArgList observables;
1453  std::vector<std::string>::iterator itr = fObsNameVec.begin();
1454  for (int idx=0; itr!=fObsNameVec.end(); ++itr, ++idx ) {
1455  observables.add( *proto->var(itr->c_str()) );
1456  }
1457 
1458  // Create the list of terms to
1459  // control the bin heights:
1460  std::string ParamSetPrefix = "gamma_stat_" + channel_name;
1461  Double_t gammaMin = 0.0;
1462  Double_t gammaMax = 10.0;
1463  RooArgList statFactorParams = ParamHistFunc::createParamSet(*proto,
1464  ParamSetPrefix.c_str(),
1465  observables,
1466  gammaMin, gammaMax);
1467 
1468  ParamHistFunc statUncertFunc(statFuncName.c_str(), statFuncName.c_str(),
1469  observables, statFactorParams );
1470 
1471  proto->import( statUncertFunc, RecycleConflictNodes() );
1472 
1473  paramHist = (ParamHistFunc*) proto->function( statFuncName.c_str() );
1474 
1475  } // END: If Statement: Create ParamHistFunc
1476 
1477  // Create the node as a product
1478  // of this function and the
1479  // expected value from MC
1480  statNodeName = sample.GetName() + "_" + channel_name + "_overallSyst_x_StatUncert";
1481 
1482  RooAbsReal* expFunc = (RooAbsReal*) proto->function( syst_x_expectedPrefix.c_str() );
1483  RooProduct nodeWithMcStat(statNodeName.c_str(), statNodeName.c_str(),
1484  RooArgSet(*paramHist, *expFunc) );
1485 
1486  proto->import( nodeWithMcStat, RecycleConflictNodes() );
1487 
1488  // Push back the final name of the node
1489  // to be used in the RooRealSumPdf
1490  // (node to be created later)
1491  syst_x_expectedPrefix = nodeWithMcStat.GetName();
1492 
1493  }
1494  } // END: if DoMcStat
1495 
1496 
1497  ///////////////////////////////////////////
1498  // Create a ShapeFactor for this channel //
1499  ///////////////////////////////////////////
1500 
1501  if( sample.GetShapeFactorList().size() > 0 ) {
1502 
1503  if( fObsNameVec.size() > 3 ) {
1504  std::cout << "Cannot include Stat Error for histograms of more than 3 dimensions."
1505  << std::endl;
1506  throw hf_exc();
1507  } else {
1508 
1509  std::cout << "Sample: " << sample.GetName() << " in channel: " << channel_name
1510  << " to be include a ShapeFactor."
1511  << std::endl;
1512 
1513  std::vector<ParamHistFunc*> paramHistFuncList;
1514  std::vector<std::string> shapeFactorNameList;
1515 
1516  for(unsigned int i=0; i < sample.GetShapeFactorList().size(); ++i) {
1517 
1518  ShapeFactor& shapeFactor = sample.GetShapeFactorList().at(i);
1519 
1520  std::string funcName = channel_name + "_" + shapeFactor.GetName() + "_shapeFactor";
1521  ParamHistFunc* paramHist = (ParamHistFunc*) proto->function( funcName.c_str() );
1522  if( paramHist == NULL ) {
1523 
1524  RooArgList observables;
1525  std::vector<std::string>::iterator itr = fObsNameVec.begin();
1526  for (int idx=0; itr!=fObsNameVec.end(); ++itr, ++idx ) {
1527  observables.add( *proto->var(itr->c_str()) );
1528  }
1529 
1530  // Create the Parameters
1531  std::string funcParams = "gamma_" + shapeFactor.GetName();
1532 
1533  // GHL: Again, we are putting hard ranges on the gamma's
1534  // We should change this to range from 0 to /inf
1535  RooArgList shapeFactorParams = ParamHistFunc::createParamSet(*proto,
1536  funcParams.c_str(),
1537  observables, 0, 1000);
1538 
1539  // Create the Function
1540  ParamHistFunc shapeFactorFunc( funcName.c_str(), funcName.c_str(),
1541  observables, shapeFactorParams );
1542 
1543  // Set an initial shape, if requested
1544  if( shapeFactor.GetInitialShape() != NULL ) {
1545  TH1* initialShape = shapeFactor.GetInitialShape();
1546  std::cout << "Setting Shape Factor: " << shapeFactor.GetName()
1547  << " to have initial shape from hist: "
1548  << initialShape->GetName()
1549  << std::endl;
1550  shapeFactorFunc.setShape( initialShape );
1551  }
1552 
1553  // Set the variables constant, if requested
1554  if( shapeFactor.IsConstant() ) {
1555  std::cout << "Setting Shape Factor: " << shapeFactor.GetName()
1556  << " to be constant" << std::endl;
1557  shapeFactorFunc.setConstant(true);
1558  }
1559 
1560  proto->import( shapeFactorFunc, RecycleConflictNodes() );
1561  paramHist = (ParamHistFunc*) proto->function( funcName.c_str() );
1562 
1563  } // End: Create ShapeFactor ParamHistFunc
1564 
1565  paramHistFuncList.push_back(paramHist);
1566  shapeFactorNameList.push_back(funcName);
1567 
1568  } // End loop over ShapeFactor Systematics
1569 
1570  // Now that we have the right ShapeFactor,
1571  // we multiply the expected function
1572 
1573  //std::string shapeFactorNodeName = syst_x_expectedPrefix + "_x_" + funcName;
1574  // Dynamically build the name as a long product
1575  std::string shapeFactorNodeName = syst_x_expectedPrefix;
1576  for( unsigned int i=0; i < shapeFactorNameList.size(); ++i) {
1577  shapeFactorNodeName += "_x_" + shapeFactorNameList.at(i);
1578  }
1579 
1580  RooAbsReal* expFunc = (RooAbsReal*) proto->function( syst_x_expectedPrefix.c_str() );
1581  RooArgSet nodesForProduct(*expFunc);
1582  for( unsigned int i=0; i < paramHistFuncList.size(); ++i) {
1583  nodesForProduct.add( *paramHistFuncList.at(i) );
1584  }
1585  //RooProduct nodeWithShapeFactor(shapeFactorNodeName.c_str(),
1586  // shapeFactorNodeName.c_str(),
1587  //RooArgSet(*paramHist, *expFunc) );
1588  RooProduct nodeWithShapeFactor(shapeFactorNodeName.c_str(),
1589  shapeFactorNodeName.c_str(),
1590  nodesForProduct );
1591 
1592  proto->import( nodeWithShapeFactor, RecycleConflictNodes() );
1593 
1594  // Push back the final name of the node
1595  // to be used in the RooRealSumPdf
1596  // (node to be created later)
1597  syst_x_expectedPrefix = nodeWithShapeFactor.GetName();
1598 
1599  }
1600  } // End: if ShapeFactorName!=""
1601 
1602 
1603  ////////////////////////////////////////
1604  // Create a ShapeSys for this channel //
1605  ////////////////////////////////////////
1606 
1607  if( sample.GetShapeSysList().size() != 0 ) {
1608 
1609  if( fObsNameVec.size() > 3 ) {
1610  std::cout << "Cannot include Stat Error for histograms of more than 3 dimensions."
1611  << std::endl;
1612  throw hf_exc();
1613  } else {
1614 
1615  // List of ShapeSys ParamHistFuncs
1616  std::vector<string> ShapeSysNames;
1617 
1618  for( unsigned int i = 0; i < sample.GetShapeSysList().size(); ++i) {
1619 
1620  // Create the ParamHistFunc's
1621  // Create their constraint terms and add them
1622  // to the list of constraint terms
1623 
1624  // Create a single RooProduct over all of these
1625  // paramHistFunc's
1626 
1627  // Send the name of that product to the RooRealSumPdf
1628 
1629  RooStats::HistFactory::ShapeSys& shapeSys = sample.GetShapeSysList().at(i);
1630 
1631  std::cout << "Sample: " << sample.GetName() << " in channel: " << channel_name
1632  << " to include a ShapeSys." << std::endl;
1633 
1634  std::string funcName = channel_name + "_" + shapeSys.GetName() + "_ShapeSys";
1635  ShapeSysNames.push_back( funcName );
1636  ParamHistFunc* paramHist = (ParamHistFunc*) proto->function( funcName.c_str() );
1637  if( paramHist == NULL ) {
1638 
1639  //std::string funcParams = "gamma_" + it->shapeFactorName;
1640  //paramHist = CreateParamHistFunc( proto, fObsNameVec, funcParams, funcName );
1641 
1642  RooArgList observables;
1643  std::vector<std::string>::iterator itr = fObsNameVec.begin();
1644  for(; itr!=fObsNameVec.end(); ++itr ) {
1645  observables.add( *proto->var(itr->c_str()) );
1646  }
1647 
1648  // Create the Parameters
1649  std::string funcParams = "gamma_" + shapeSys.GetName();
1650  RooArgList shapeFactorParams = ParamHistFunc::createParamSet(*proto,
1651  funcParams.c_str(),
1652  observables, 0, 10);
1653 
1654  // Create the Function
1655  ParamHistFunc shapeFactorFunc( funcName.c_str(), funcName.c_str(),
1656  observables, shapeFactorParams );
1657 
1658  proto->import( shapeFactorFunc, RecycleConflictNodes() );
1659  paramHist = (ParamHistFunc*) proto->function( funcName.c_str() );
1660 
1661  } // End: Create ShapeFactor ParamHistFunc
1662 
1663  // Create the constraint terms and add
1664  // them to the workspace (proto)
1665  // as well as the list of constraint terms (constraintTermNames)
1666 
1667  // The syst should be a fractional error
1668  TH1* shapeErrorHist = shapeSys.GetErrorHist();
1669 
1670  // Constraint::Type shapeConstraintType = Constraint::Gaussian;
1671  Constraint::Type systype = shapeSys.GetConstraintType();
1672  if( systype == Constraint::Gaussian) {
1673  systype = Constraint::Gaussian;
1674  }
1675  if( systype == Constraint::Poisson ) {
1676  systype = Constraint::Poisson;
1677  }
1678 
1679  Double_t minShapeUncertainty = 0.0;
1680  RooArgList shapeConstraints = createStatConstraintTerms(proto, constraintTermNames,
1681  *paramHist, shapeErrorHist,
1682  systype,
1683  minShapeUncertainty);
1684 
1685  } // End: Loop over ShapeSys vector in this EstimateSummary
1686 
1687  // Now that we have the list of ShapeSys ParamHistFunc names,
1688  // we create the total RooProduct
1689  // we multiply the expected functio
1690 
1691  std::string NodeName = syst_x_expectedPrefix;
1692  RooArgList ShapeSysForNode;
1693  RooAbsReal* expFunc = (RooAbsReal*) proto->function( syst_x_expectedPrefix.c_str() );
1694  ShapeSysForNode.add( *expFunc );
1695  for( unsigned int i = 0; i < ShapeSysNames.size(); ++i ) {
1696  std::string ShapeSysName = ShapeSysNames.at(i);
1697  ShapeSysForNode.add( *proto->function(ShapeSysName.c_str()) );
1698  NodeName = NodeName + "_x_" + ShapeSysName;
1699  }
1700 
1701  // Create the name for this NEW Node
1702  RooProduct nodeWithShapeFactor(NodeName.c_str(), NodeName.c_str(), ShapeSysForNode );
1703  proto->import( nodeWithShapeFactor, RecycleConflictNodes() );
1704 
1705  // Push back the final name of the node
1706  // to be used in the RooRealSumPdf
1707  // (node to be created later)
1708  syst_x_expectedPrefix = nodeWithShapeFactor.GetName();
1709 
1710  } // End: NumObsVar == 1
1711 
1712  } // End: GetShapeSysList.size() != 0
1713 
1714  // Append the name of the "node"
1715  // that is to be summed with the
1716  // RooRealSumPdf
1717  syst_x_expectedPrefixNames.push_back(syst_x_expectedPrefix);
1718 
1719  // GHL: This was pretty confusing before,
1720  // hopefully using the measurement directly
1721  // will improve it
1722  if( sample.GetNormalizeByTheory() ) {
1723  normalizationNames.push_back( "Lumi" );
1724  }
1725  else {
1726  TString lumiParamString;
1727  lumiParamString += measurement.GetLumi();
1728  lumiParamString.ReplaceAll(' ', TString());
1729  normalizationNames.push_back(lumiParamString.Data());
1730  }
1731 
1732  } // END: Loop over EstimateSummaries
1733  // proto->Print();
1734 
1735  // If a non-zero number of samples call for
1736  // Stat Uncertainties, create the statFactor functions
1737  if( statHistPairs.size() > 0 ) {
1738 
1739  // Create the histogram of (binwise)
1740  // stat uncertainties:
1741  TH1* fracStatError = MakeScaledUncertaintyHist( statNodeName + "_RelErr", statHistPairs );
1742  if( fracStatError == NULL ) {
1743  std::cout << "Error: Failed to make ScaledUncertaintyHist for: "
1744  << statNodeName << std::endl;
1745  throw hf_exc();
1746  }
1747 
1748  // Using this TH1* of fractinal stat errors,
1749  // create a set of constraint terms:
1750  ParamHistFunc* chanStatUncertFunc = (ParamHistFunc*) proto->function( statFuncName.c_str() );
1751  std::cout << "About to create Constraint Terms from: "
1752  << chanStatUncertFunc->GetName()
1753  << " params: " << chanStatUncertFunc->paramList()
1754  << std::endl;
1755 
1756  // Get the constraint type and the
1757  // rel error threshold from the (last)
1758  // EstimateSummary looped over (but all
1759  // should be the same)
1760 
1761  // Get the type of StatError constraint from the channel
1762  Constraint::Type statConstraintType = channel.GetStatErrorConfig().GetConstraintType();
1763  if( statConstraintType == Constraint::Gaussian) {
1764  std::cout << "Using Gaussian StatErrors in channel: " << channel.GetName() << std::endl;
1765  }
1766  if( statConstraintType == Constraint::Poisson ) {
1767  std::cout << "Using Poisson StatErrors in channel: " << channel.GetName() << std::endl;
1768  }
1769 
1770  double statRelErrorThreshold = channel.GetStatErrorConfig().GetRelErrorThreshold();
1771  RooArgList statConstraints = createStatConstraintTerms(proto, constraintTermNames,
1772  *chanStatUncertFunc, fracStatError,
1773  statConstraintType,
1774  statRelErrorThreshold);
1775 
1776 
1777  // clean stat hist pair (need to delete second histogram)
1778  for (unsigned int i = 0; i < statHistPairs.size() ; ++i )
1779  delete statHistPairs[i].second;
1780 
1781  statHistPairs.clear();
1782  //delete also histogram of stat uncertainties created in MakeScaledUncertaintyHist
1783  delete fracStatError;
1784 
1785  } // END: Loop over stat Hist Pairs
1786 
1787 
1788  ///////////////////////////////////
1789  // for ith bin calculate totN_i = lumi * sum_j expected_j * syst_j
1790  //MakeTotalExpected(proto,channel_name+"_model",channel_name,"Lumi",fLowBin,fHighBin,
1791  // syst_x_expectedPrefixNames, normalizationNames);
1792  MakeTotalExpected(proto, channel_name+"_model", //channel_name,"Lumi",fLowBin,fHighBin,
1793  syst_x_expectedPrefixNames, normalizationNames);
1794  likelihoodTermNames.push_back(channel_name+"_model");
1795 
1796  //////////////////////////////////////
1797  // fix specified parameters
1798  for(unsigned int i=0; i<systToFix.size(); ++i){
1799  RooRealVar* temp = proto->var((systToFix.at(i)).c_str());
1800  if(temp) {
1801  // set the parameter constant
1802  temp->setConstant();
1803 
1804  // remove the corresponding auxiliary observable from the global observables
1805  RooRealVar* auxMeas = NULL;
1806  if(systToFix.at(i)=="Lumi"){
1807  auxMeas = proto->var("nominalLumi");
1808  } else {
1809  auxMeas = proto->var(TString::Format("nom_%s",temp->GetName()));
1810  }
1811 
1812  if(auxMeas){
1813  const_cast<RooArgSet*>(proto->set("globalObservables"))->remove(*auxMeas);
1814  } else{
1815  cout << "could not corresponding auxiliary measurement "
1816  << TString::Format("nom_%s",temp->GetName()) << endl;
1817  }
1818  } else {
1819  cout << "could not find variable " << systToFix.at(i)
1820  << " could not set it to constant" << endl;
1821  }
1822  }
1823 
1824  //////////////////////////////////////
1825  // final proto model
1826  for(unsigned int i=0; i<constraintTermNames.size(); ++i){
1827  RooAbsArg* proto_arg = (proto->arg(constraintTermNames[i].c_str()));
1828  if( proto_arg==NULL ) {
1829  std::cout << "Error: Cannot find arg set: " << constraintTermNames.at(i)
1830  << " in workspace: " << proto->GetName() << std::endl;
1831  throw hf_exc();
1832  }
1833  constraintTerms.add( *proto_arg );
1834  // constraintTerms.add(* proto_arg(proto->arg(constraintTermNames[i].c_str())) );
1835  }
1836  for(unsigned int i=0; i<likelihoodTermNames.size(); ++i){
1837  RooAbsArg* proto_arg = (proto->arg(likelihoodTermNames[i].c_str()));
1838  if( proto_arg==NULL ) {
1839  std::cout << "Error: Cannot find arg set: " << likelihoodTermNames.at(i)
1840  << " in workspace: " << proto->GetName() << std::endl;
1841  throw hf_exc();
1842  }
1843  likelihoodTerms.add( *proto_arg );
1844  }
1845  proto->defineSet("constraintTerms",constraintTerms);
1846  proto->defineSet("likelihoodTerms",likelihoodTerms);
1847  // proto->Print();
1848 
1849  // list of observables
1850  RooArgList observables;
1851  std::string observablesStr;
1852 
1853  std::vector<std::string>::iterator itr = fObsNameVec.begin();
1854  for(; itr!=fObsNameVec.end(); ++itr ) {
1855  observables.add( *proto->var(itr->c_str()) );
1856  if (!observablesStr.empty()) { observablesStr += ","; }
1857  observablesStr += *itr;
1858  }
1859 
1860  // We create two sets, one for backwards compatability
1861  // The other to make a consistent naming convention
1862  // between individual channels and the combined workspace
1863  proto->defineSet("observables", TString::Format("%s",observablesStr.c_str()));
1864  proto->defineSet("observablesSet", TString::Format("%s",observablesStr.c_str()));
1865 
1866  // Create the ParamHistFunc
1867  // after observables have been made
1868  cout <<"-----------------------------------------"<<endl;
1869  cout <<"import model into workspace" << endl;
1870 
1871  RooProdPdf* model = new RooProdPdf(("model_"+channel_name).c_str(), // MB : have changed this into conditional pdf. Much faster for toys!
1872  "product of Poissons accross bins for a single channel",
1873  constraintTerms, Conditional(likelihoodTerms,observables)); //likelihoodTerms);
1874  proto->import(*model,RecycleConflictNodes());
1875 
1876  proto_config->SetPdf(*model);
1877  proto_config->SetObservables(observables);
1878  proto_config->SetGlobalObservables(*proto->set("globalObservables"));
1879  // proto->writeToFile(("results/model_"+channel+".root").c_str());
1880  // fill out nuisance parameters in model config
1881  // proto_config->GuessObsAndNuisance(*proto->data("asimovData"));
1882  proto->import(*proto_config,proto_config->GetName());
1883  proto->importClassCode();
1884 
1885  ///////////////////////////
1886  // make data sets
1887  // THis works and is natural, but the memory size of the simultaneous dataset grows exponentially with channels
1888  const char* weightName="weightVar";
1889  proto->factory(TString::Format("%s[0,-1e10,1e10]",weightName));
1890  proto->defineSet("obsAndWeight",TString::Format("%s,%s",weightName,observablesStr.c_str()));
1891 
1892  /* Old code for generating the asimov
1893  Asimov generation is now done later...
1894 
1895  RooAbsData* asimov_data = model->generateBinned(observables,ExpectedData());
1896 
1897  /// Asimov dataset
1898  RooDataSet* asimovDataUnbinned = new RooDataSet("asimovData","",*proto->set("obsAndWeight"),weightName);
1899  for(int i=0; i<asimov_data->numEntries(); ++i){
1900  asimov_data->get(i)->Print("v");
1901  //cout << "GREPME : " << i << " " << data->weight() <<endl;
1902  asimovDataUnbinned->add( *asimov_data->get(i), asimov_data->weight() );
1903  }
1904  proto->import(*asimovDataUnbinned);
1905  */
1906 
1907  // New Asimov Generation: Use the code in the Asymptotic calculator
1908  // Need to get the ModelConfig...
1909  RooDataSet* asimov_dataset = (RooDataSet*) AsymptoticCalculator::GenerateAsimovData(*model, observables);
1910  proto->import(*asimov_dataset, Rename("asimovData"));
1911 
1912  // GHL: Determine to use data if the hist isn't 'NULL'
1913  if(channel.GetData().GetHisto() != NULL) {
1914 
1915  Data& data = channel.GetData();
1916  TH1* mnominal = data.GetHisto();
1917  if( !mnominal ) {
1918  std::cout << "Error: Data histogram for channel: " << channel.GetName()
1919  << " is NULL" << std::endl;
1920  throw hf_exc();
1921  }
1922 
1923  // THis works and is natural, but the memory size of the simultaneous dataset grows exponentially with channels
1924  RooDataSet* obsDataUnbinned = new RooDataSet("obsData","",*proto->set("obsAndWeight"),weightName);
1925 
1926 
1927  ConfigureHistFactoryDataset( obsDataUnbinned, mnominal,
1928  proto, fObsNameVec );
1929 
1930  /*
1931  //ES// TH1* mnominal = summary.at(0).nominal;
1932  TH1* mnominal = data.GetHisto();
1933  TAxis* ax = mnominal->GetXaxis();
1934  TAxis* ay = mnominal->GetYaxis();
1935  TAxis* az = mnominal->GetZaxis();
1936 
1937  for (int i=1; i<=ax->GetNbins(); ++i) { // 1 or more dimension
1938  Double_t xval = ax->GetBinCenter(i);
1939  proto->var( fObsNameVec[0].c_str() )->setVal( xval );
1940  if (fObsNameVec.size()==1) {
1941  Double_t fval = mnominal->GetBinContent(i);
1942  obsDataUnbinned->add( *proto->set("obsAndWeight"), fval );
1943  } else { // 2 or more dimensions
1944  for (int j=1; j<=ay->GetNbins(); ++j) {
1945  Double_t yval = ay->GetBinCenter(j);
1946  proto->var( fObsNameVec[1].c_str() )->setVal( yval );
1947  if (fObsNameVec.size()==2) {
1948  Double_t fval = mnominal->GetBinContent(i,j);
1949  obsDataUnbinned->add( *proto->set("obsAndWeight"), fval );
1950  } else { // 3 dimensions
1951  for (int k=1; k<=az->GetNbins(); ++k) {
1952  Double_t zval = az->GetBinCenter(k);
1953  proto->var( fObsNameVec[2].c_str() )->setVal( zval );
1954  Double_t fval = mnominal->GetBinContent(i,j,k);
1955  obsDataUnbinned->add( *proto->set("obsAndWeight"), fval );
1956  }
1957  }
1958  }
1959  }
1960  }
1961  */
1962 
1963  proto->import(*obsDataUnbinned);
1964  delete obsDataUnbinned;
1965  } // End: Has non-null 'data' entry
1966 
1967 
1968  for(unsigned int i=0; i < channel.GetAdditionalData().size(); ++i) {
1969 
1970  Data& data = channel.GetAdditionalData().at(i);
1971  std::string dataName = data.GetName();
1972  TH1* mnominal = data.GetHisto();
1973  if( !mnominal ) {
1974  std::cout << "Error: Additional Data histogram for channel: " << channel.GetName()
1975  << " with name: " << dataName << " is NULL" << std::endl;
1976  throw hf_exc();
1977  }
1978 
1979  // THis works and is natural, but the memory size of the simultaneous dataset grows exponentially with channels
1980  RooDataSet* obsDataUnbinned = new RooDataSet(dataName.c_str(), dataName.c_str(),
1981  *proto->set("obsAndWeight"), weightName);
1982 
1983  ConfigureHistFactoryDataset( obsDataUnbinned, mnominal,
1984  proto, fObsNameVec );
1985 
1986  /*
1987  //ES// TH1* mnominal = summary.at(0).nominal;
1988  TH1* mnominal = data.GetHisto();
1989  TAxis* ax = mnominal->GetXaxis();
1990  TAxis* ay = mnominal->GetYaxis();
1991  TAxis* az = mnominal->GetZaxis();
1992 
1993  for (int i=1; i<=ax->GetNbins(); ++i) { // 1 or more dimension
1994  Double_t xval = ax->GetBinCenter(i);
1995  proto->var( fObsNameVec[0].c_str() )->setVal( xval );
1996  if (fObsNameVec.size()==1) {
1997  Double_t fval = mnominal->GetBinContent(i);
1998  obsDataUnbinned->add( *proto->set("obsAndWeight"), fval );
1999  } else { // 2 or more dimensions
2000  for (int j=1; j<=ay->GetNbins(); ++j) {
2001  Double_t yval = ay->GetBinCenter(j);
2002  proto->var( fObsNameVec[1].c_str() )->setVal( yval );
2003  if (fObsNameVec.size()==2) {
2004  Double_t fval = mnominal->GetBinContent(i,j);
2005  obsDataUnbinned->add( *proto->set("obsAndWeight"), fval );
2006  } else { // 3 dimensions
2007  for (int k=1; k<=az->GetNbins(); ++k) {
2008  Double_t zval = az->GetBinCenter(k);
2009  proto->var( fObsNameVec[2].c_str() )->setVal( zval );
2010  Double_t fval = mnominal->GetBinContent(i,j,k);
2011  obsDataUnbinned->add( *proto->set("obsAndWeight"), fval );
2012  }
2013  }
2014  }
2015  }
2016  }
2017  */
2018 
2019  proto->import(*obsDataUnbinned);
2020 
2021  delete obsDataUnbinned;
2022 
2023  } // End: Has non-null 'data' entry
2024 
2025  // clean up
2026  delete model;
2027  delete proto_config;
2028  delete asimov_dataset;
2029 
2030  proto->Print();
2031  return proto;
2032  }
2033 
2034 
2036  TH1* mnominal,
2038  std::vector<std::string> ObsNameVec) {
2039 
2040  // Take a RooDataSet and fill it with the entries
2041  // from a TH1*, using the observable names to
2042  // determine the columns
2043 
2044  if (ObsNameVec.empty() ) {
2045  Error("ConfigureHistFactoryDataset","Invalid input - return");
2046  return;
2047  }
2048 
2049  //ES// TH1* mnominal = summary.at(0).nominal;
2050  // TH1* mnominal = data.GetHisto();
2051  TAxis* ax = mnominal->GetXaxis();
2052  TAxis* ay = mnominal->GetYaxis();
2053  TAxis* az = mnominal->GetZaxis();
2054 
2055  for (int i=1; i<=ax->GetNbins(); ++i) { // 1 or more dimension
2056 
2057  Double_t xval = ax->GetBinCenter(i);
2058  proto->var( ObsNameVec[0].c_str() )->setVal( xval );
2059 
2060  if(ObsNameVec.size()==1) {
2061  Double_t fval = mnominal->GetBinContent(i);
2062  obsDataUnbinned->add( *proto->set("obsAndWeight"), fval );
2063  } else { // 2 or more dimensions
2064 
2065  for(int j=1; j<=ay->GetNbins(); ++j) {
2066  Double_t yval = ay->GetBinCenter(j);
2067  proto->var( ObsNameVec[1].c_str() )->setVal( yval );
2068 
2069  if(ObsNameVec.size()==2) {
2070  Double_t fval = mnominal->GetBinContent(i,j);
2071  obsDataUnbinned->add( *proto->set("obsAndWeight"), fval );
2072  } else { // 3 dimensions
2073 
2074  for(int k=1; k<=az->GetNbins(); ++k) {
2075  Double_t zval = az->GetBinCenter(k);
2076  proto->var( ObsNameVec[2].c_str() )->setVal( zval );
2077  Double_t fval = mnominal->GetBinContent(i,j,k);
2078  obsDataUnbinned->add( *proto->set("obsAndWeight"), fval );
2079  }
2080  }
2081  }
2082  }
2083  }
2084  }
2085 
2087  {
2088  fObsNameVec.clear();
2089 
2090  // determine histogram dimensionality
2091  unsigned int histndim(1);
2092  std::string classname = hist->ClassName();
2093  if (classname.find("TH1")==0) { histndim=1; }
2094  else if (classname.find("TH2")==0) { histndim=2; }
2095  else if (classname.find("TH3")==0) { histndim=3; }
2096 
2097  for ( unsigned int idx=0; idx<histndim; ++idx ) {
2098  if (idx==0) { fObsNameVec.push_back("x"); }
2099  if (idx==1) { fObsNameVec.push_back("y"); }
2100  if (idx==2) { fObsNameVec.push_back("z"); }
2101  }
2102  }
2103 
2104 
2105  RooWorkspace* HistoToWorkspaceFactoryFast::MakeCombinedModel(vector<string> ch_names, vector<RooWorkspace*> chs)
2106  {
2107 
2108 
2109  // check first the inputs (see JIRA-6890)
2110  if (ch_names.empty() || chs.empty() ) {
2111  Error("MakeCombinedModel","Input vectors are empty - return a nullptr");
2112  return 0;
2113  }
2114  if (chs.size() < ch_names.size() ) {
2115  Error("MakeCombinedModel","Input vector of workspace has an invalid size - return a nullptr");
2116  return 0;
2117  }
2118 
2119  //
2120  /// These things were used for debugging. Maybe useful in the future
2121  //
2122 
2123  map<string, RooAbsPdf*> pdfMap;
2124  vector<RooAbsPdf*> models;
2125  stringstream ss;
2126 
2127  RooArgList obsList;
2128  for(unsigned int i = 0; i< ch_names.size(); ++i){
2129  ModelConfig * config = (ModelConfig *) chs[i]->obj("ModelConfig");
2130  obsList.add(*config->GetObservables());
2131  }
2132  cout <<"full list of observables:"<<endl;
2133  obsList.Print();
2134 
2135  RooArgSet globalObs;
2136  for(unsigned int i = 0; i< ch_names.size(); ++i){
2137  string channel_name=ch_names[i];
2138 
2139  if (ss.str().empty()) ss << channel_name ;
2140  else ss << ',' << channel_name ;
2141  RooWorkspace * ch=chs[i];
2142 
2143  RooAbsPdf* model = ch->pdf(("model_"+channel_name).c_str());
2144  if(!model) cout <<"failed to find model for channel"<<endl;
2145  // cout << "int = " << model->createIntegral(*obsN)->getVal() << endl;;
2146  models.push_back(model);
2147  globalObs.add(*ch->set("globalObservables"));
2148 
2149  // constrainedParams->add( * ch->set("constrainedParams") );
2150  pdfMap[channel_name]=model;
2151  }
2152  //constrainedParams->Print();
2153 
2154  cout << "\n\n------------------\n Entering combination" << endl;
2155  RooWorkspace* combined = new RooWorkspace("combined");
2156  // RooWorkspace* combined = chs[0];
2157 
2158 
2159  RooCategory* channelCat = (RooCategory*) combined->factory(("channelCat["+ss.str()+"]").c_str());
2160  RooSimultaneous * simPdf= new RooSimultaneous("simPdf","",pdfMap, *channelCat);
2161  ModelConfig * combined_config = new ModelConfig("ModelConfig", combined);
2162  combined_config->SetWorkspace(*combined);
2163  // combined_config->SetNuisanceParameters(*constrainedParams);
2164 
2165  combined->import(globalObs);
2166  combined->defineSet("globalObservables",globalObs);
2167  combined_config->SetGlobalObservables(*combined->set("globalObservables"));
2168 
2169 
2170  ////////////////////////////////////////////
2171  // Make toy simultaneous dataset
2172  cout <<"-----------------------------------------"<<endl;
2173  cout << "create toy data for " << ss.str() << endl;
2174 
2175 
2176  // now with weighted datasets
2177  // First Asimov
2178  //RooDataSet * simData=NULL;
2179  combined->factory("weightVar[0,-1e10,1e10]");
2180  obsList.add(*combined->var("weightVar"));
2181 
2182  // Loop over channels and create the asimov
2183  /*
2184  for(unsigned int i = 0; i< ch_names.size(); ++i){
2185  cout << "merging data for channel " << ch_names[i].c_str() << endl;
2186  RooDataSet * tempData=new RooDataSet(ch_names[i].c_str(),"", obsList, Index(*channelCat),
2187  WeightVar("weightVar"),
2188  Import(ch_names[i].c_str(),*(RooDataSet*)chs[i]->data("asimovData")));
2189  if(simData){
2190  simData->append(*tempData);
2191  delete tempData;
2192  }else{
2193  simData = tempData;
2194  }
2195  }
2196 
2197  if (simData) combined->import(*simData,Rename("asimovData"));
2198  */
2199  RooDataSet* asimov_combined = (RooDataSet*) AsymptoticCalculator::GenerateAsimovData(*simPdf,
2200  obsList);
2201  if( asimov_combined ) {
2202  combined->import( *asimov_combined, Rename("asimovData"));
2203  }
2204  else {
2205  std::cout << "Error: Failed to create combined asimov dataset" << std::endl;
2206  throw hf_exc();
2207  }
2208  delete asimov_combined;
2209 
2210  // Now merge the observable datasets across the channels
2211  if(chs[0]->data("obsData") != NULL) {
2212  MergeDataSets(combined, chs, ch_names, "obsData", obsList, channelCat);
2213  }
2214 
2215  /*
2216  if(chs[0]->data("obsData") != NULL){
2217  RooDataSet * simData=NULL;
2218  //simData=NULL;
2219 
2220  // Loop through channels, get their individual datasets,
2221  // and add them to the combined dataset
2222  for(unsigned int i = 0; i< ch_names.size(); ++i){
2223  cout << "merging data for channel " << ch_names[i].c_str() << endl;
2224 
2225  RooDataSet* obsDataInChannel = (RooDataSet*) chs[i]->data("obsData");
2226  RooDataSet * tempData = new RooDataSet(ch_names[i].c_str(),"", obsList, Index(*channelCat),
2227  WeightVar("weightVar"),
2228  Import(ch_names[i].c_str(),*obsDataInChannel));
2229  // *(RooDataSet*) chs[i]->data("obsData")));
2230  if(simData) {
2231  simData->append(*tempData);
2232  delete tempData;
2233  }
2234  else {
2235  simData = tempData;
2236  }
2237  } // End Loop Over Channels
2238 
2239  // Check that we successfully created the dataset
2240  // and import it into the workspace
2241  if(simData) {
2242  combined->import(*simData, Rename("obsData"));
2243  }
2244  else {
2245  std::cout << "Error: Unable to merge observable datasets" << std::endl;
2246  throw hf_exc();
2247  }
2248 
2249  } // End 'if' on data != NULL
2250  */
2251 
2252  // Now, create any additional Asimov datasets that
2253  // are configured in the measurement
2254 
2255 
2256  // obsList.Print();
2257  // combined->import(obsList);
2258  // combined->Print();
2259 
2260  obsList.add(*channelCat);
2261  combined->defineSet("observables",obsList);
2262  combined_config->SetObservables(*combined->set("observables"));
2263 
2264  combined->Print();
2265 
2266  cout << "\n\n----------------\n Importing combined model" << endl;
2267  combined->import(*simPdf,RecycleConflictNodes());
2268  //combined->import(*simPdf, RenameVariable("SigXsecOverSM","SigXsecOverSM_comb"));
2269  // cout << "check pointer " << simPdf << endl;
2270  // cout << "check val " << simPdf->getVal() << endl;
2271 
2272  std::map< std::string, double>::iterator param_itr = fParamValues.begin();
2273  for( ; param_itr != fParamValues.end(); ++param_itr ){
2274  // make sure they are fixed
2275  std::string paramName = param_itr->first;
2276  double paramVal = param_itr->second;
2277 
2278  RooRealVar* temp = combined->var( paramName.c_str() );
2279  if(temp) {
2280  temp->setVal( paramVal );
2281  cout <<"setting " << paramName << " to the value: " << paramVal << endl;
2282  } else
2283  cout << "could not find variable " << paramName << " could not set its value" << endl;
2284  }
2285 
2286 
2287  for(unsigned int i=0; i<fSystToFix.size(); ++i){
2288  // make sure they are fixed
2289  RooRealVar* temp = combined->var((fSystToFix.at(i)).c_str());
2290  if(temp) {
2291  temp->setConstant();
2292  cout <<"setting " << fSystToFix.at(i) << " constant" << endl;
2293  } else
2294  cout << "could not find variable " << fSystToFix.at(i) << " could not set it to constant" << endl;
2295  }
2296 
2297  ///
2298  /// writing out the model in graphViz
2299  ///
2300  // RooAbsPdf* customized=combined->pdf("simPdf");
2301  //combined_config->SetPdf(*customized);
2302  combined_config->SetPdf(*simPdf);
2303  // combined_config->GuessObsAndNuisance(*simData);
2304  // customized->graphVizTree(("results/"+fResultsPrefixStr.str()+"_simul.dot").c_str());
2305  combined->import(*combined_config,combined_config->GetName());
2306  combined->importClassCode();
2307  // combined->writeToFile("results/model_combined.root");
2308 
2309  //clean up
2310  delete combined_config;
2311  delete simPdf;
2312 
2313  return combined;
2314  }
2315 
2316 
2318  std::vector<RooWorkspace*> wspace_vec,
2319  std::vector<std::string> channel_names,
2320  std::string dataSetName,
2321  RooArgList obsList,
2322  RooCategory* channelCat) {
2323 
2324  // Create the total dataset
2325  RooDataSet* simData=NULL;
2326 
2327  // Loop through channels, get their individual datasets,
2328  // and add them to the combined dataset
2329  for(unsigned int i = 0; i< channel_names.size(); ++i){
2330 
2331  // Grab the dataset for the existing channel
2332  std::cout << "Merging data for channel " << channel_names[i].c_str() << std::endl;
2333  RooDataSet* obsDataInChannel = (RooDataSet*) wspace_vec[i]->data(dataSetName.c_str());
2334  if( !obsDataInChannel ) {
2335  std::cout << "Error: Can't find DataSet: " << dataSetName
2336  << " in channel: " << channel_names.at(i)
2337  << std::endl;
2338  throw hf_exc();
2339  }
2340 
2341  // Create the new Dataset
2342  RooDataSet * tempData = new RooDataSet(channel_names[i].c_str(),"",
2343  obsList, Index(*channelCat),
2344  WeightVar("weightVar"),
2345  Import(channel_names[i].c_str(),*obsDataInChannel));
2346  if(simData) {
2347  simData->append(*tempData);
2348  delete tempData;
2349  }
2350  else {
2351  simData = tempData;
2352  }
2353  } // End Loop Over Channels
2354 
2355  // Check that we successfully created the dataset
2356  // and import it into the workspace
2357  if(simData) {
2358  combined->import(*simData, Rename(dataSetName.c_str()));
2359  delete simData;
2360  simData = (RooDataSet*) combined->data(dataSetName.c_str() );
2361  }
2362  else {
2363  std::cout << "Error: Unable to merge observable datasets" << std::endl;
2364  throw hf_exc();
2365  }
2366 
2367  return simData;
2368 
2369  }
2370 
2371 
2372  TH1* HistoToWorkspaceFactoryFast::MakeAbsolUncertaintyHist( const std::string& Name, const TH1* Nominal ) {
2373 
2374  // Take a nominal TH1* and create
2375  // a TH1 representing the binwise
2376  // errors (taken from the nominal TH1)
2377 
2378  TH1* ErrorHist = (TH1*) Nominal->Clone( Name.c_str() );
2379  ErrorHist->Reset();
2380 
2381  Int_t numBins = Nominal->GetNbinsX()*Nominal->GetNbinsY()*Nominal->GetNbinsZ();
2382  Int_t binNumber = 0;
2383 
2384  // Loop over bins
2385  for( Int_t i_bin = 0; i_bin < numBins; ++i_bin) {
2386 
2387  binNumber++;
2388  // Ignore underflow / overflow
2389  while( Nominal->IsBinUnderflow(binNumber) || Nominal->IsBinOverflow(binNumber) ){
2390  binNumber++;
2391  }
2392 
2393  Double_t histError = Nominal->GetBinError( binNumber );
2394 
2395  // Check that histError != NAN
2396  if( histError != histError ) {
2397  std::cout << "Warning: In histogram " << Nominal->GetName()
2398  << " bin error for bin " << i_bin
2399  << " is NAN. Not using Error!!!"
2400  << std::endl;
2401  throw hf_exc();
2402  //histError = sqrt( histContent );
2403  //histError = 0;
2404  }
2405 
2406  // Check that histError ! < 0
2407  if( histError < 0 ) {
2408  std::cout << "Warning: In histogram " << Nominal->GetName()
2409  << " bin error for bin " << binNumber
2410  << " is < 0. Setting Error to 0"
2411  << std::endl;
2412  //histError = sqrt( histContent );
2413  histError = 0;
2414  }
2415 
2416  ErrorHist->SetBinContent( binNumber, histError );
2417 
2418  }
2419 
2420  return ErrorHist;
2421 
2422  }
2423 
2424  TH1* HistoToWorkspaceFactoryFast::MakeScaledUncertaintyHist( const std::string& Name, std::vector< std::pair<TH1*, TH1*> > HistVec ) {
2425 
2426  // Take a list of < nominal, absolError > TH1* pairs
2427  // and construct a single histogram representing the
2428  // total fractional error as:
2429 
2430  // UncertInQuad(bin i) = Sum: absolUncert*absolUncert
2431  // Total(bin i) = Sum: Value
2432  //
2433  // TotalFracError(bin i) = Sqrt( UncertInQuad(i) ) / TotalBin(i)
2434 
2435 
2436  unsigned int numHists = HistVec.size();
2437 
2438  if( numHists == 0 ) {
2439  std::cout << "Warning: Empty Hist Vector, cannot create total uncertainty" << std::endl;
2440  return NULL;
2441  }
2442 
2443  TH1* HistTemplate = HistVec.at(0).first;
2444  Int_t numBins = HistTemplate->GetNbinsX()*HistTemplate->GetNbinsY()*HistTemplate->GetNbinsZ();
2445 
2446  // Check that all histograms
2447  // have the same bins
2448  for( unsigned int i = 0; i < HistVec.size(); ++i ) {
2449 
2450  TH1* nominal = HistVec.at(i).first;
2451  TH1* error = HistVec.at(i).second;
2452 
2453  if( nominal->GetNbinsX()*nominal->GetNbinsY()*nominal->GetNbinsZ() != numBins ) {
2454  std::cout << "Error: Provided hists have unequal bins" << std::endl;
2455  return NULL;
2456  }
2457  if( error->GetNbinsX()*error->GetNbinsY()*error->GetNbinsZ() != numBins ) {
2458  std::cout << "Error: Provided hists have unequal bins" << std::endl;
2459  return NULL;
2460  }
2461  }
2462 
2463  std::vector<double> TotalBinContent( numBins, 0.0);
2464  std::vector<double> HistErrorsSqr( numBins, 0.0);
2465 
2466  Int_t binNumber = 0;
2467 
2468  // Loop over bins
2469  for( Int_t i_bins = 0; i_bins < numBins; ++i_bins) {
2470 
2471  binNumber++;
2472  while( HistTemplate->IsBinUnderflow(binNumber) || HistTemplate->IsBinOverflow(binNumber) ){
2473  binNumber++;
2474  }
2475 
2476  for( unsigned int i_hist = 0; i_hist < numHists; ++i_hist ) {
2477 
2478  TH1* nominal = HistVec.at(i_hist).first;
2479  TH1* error = HistVec.at(i_hist).second;
2480 
2481  //Int_t binNumber = i_bins + 1;
2482 
2483  Double_t histValue = nominal->GetBinContent( binNumber );
2484  Double_t histError = error->GetBinContent( binNumber );
2485  /*
2486  std::cout << " Getting Bin content for Stat Uncertainty"
2487  << " Nom name: " << nominal->GetName()
2488  << " Err name: " << error->GetName()
2489  << " HistNumber: " << i_hist << " bin: " << binNumber
2490  << " Value: " << histValue << " Error: " << histError
2491  << std::endl;
2492  */
2493 
2494  if( histError != histError ) {
2495  std::cout << "Warning: In histogram " << error->GetName()
2496  << " bin error for bin " << binNumber
2497  << " is NAN. Not using error!!"
2498  << std::endl;
2499  throw hf_exc();
2500  //histError = 0;
2501  }
2502 
2503  TotalBinContent.at(i_bins) += histValue;
2504  HistErrorsSqr.at(i_bins) += histError*histError; // Add in quadrature
2505 
2506  }
2507  }
2508 
2509  binNumber = 0;
2510 
2511  // Creat the output histogram
2512  TH1* ErrorHist = (TH1*) HistTemplate->Clone( Name.c_str() );
2513  ErrorHist->Reset();
2514 
2515  // Fill the output histogram
2516  for( Int_t i = 0; i < numBins; ++i) {
2517 
2518  // Int_t binNumber = i + 1;
2519  binNumber++;
2520  while( ErrorHist->IsBinUnderflow(binNumber) || ErrorHist->IsBinOverflow(binNumber) ){
2521  binNumber++;
2522  }
2523 
2524  Double_t ErrorsSqr = HistErrorsSqr.at(i);
2525  Double_t TotalVal = TotalBinContent.at(i);
2526 
2527  if( TotalVal <= 0 ) {
2528  std::cout << "Warning: Sum of histograms for bin: " << binNumber
2529  << " is <= 0. Setting error to 0"
2530  << std::endl;
2531 
2532  ErrorHist->SetBinContent( binNumber, 0.0 );
2533  continue;
2534  }
2535 
2536  Double_t RelativeError = sqrt(ErrorsSqr) / TotalVal;
2537 
2538  // If we otherwise get a NAN
2539  // it's an error
2540  if( RelativeError != RelativeError ) {
2541  std::cout << "Error: bin " << i << " error is NAN" << std::endl;
2542  std::cout << " HistErrorsSqr: " << ErrorsSqr
2543  << " TotalVal: " << TotalVal
2544  << std::endl;
2545  throw hf_exc();
2546  }
2547 
2548  // 0th entry in vector is
2549  // the 1st bin in TH1
2550  // (we ignore underflow)
2551 
2552  ErrorHist->SetBinContent( binNumber, RelativeError );
2553 
2554  std::cout << "Making Total Uncertainty for bin " << binNumber
2555  << " Error = " << sqrt(ErrorsSqr)
2556  << " Val = " << TotalVal
2557  << " RelativeError = " << RelativeError
2558  << std::endl;
2559 
2560  }
2561 
2562  return ErrorHist;
2563 
2564 }
2565 
2566 
2567 
2569  createStatConstraintTerms( RooWorkspace* proto, vector<string>& constraintTermNames,
2570  ParamHistFunc& paramHist, TH1* uncertHist,
2571  Constraint::Type type, Double_t minSigma ) {
2572 
2573 
2574  // Take a RooArgList of RooAbsReal's and
2575  // create N constraint terms (one for
2576  // each gamma) whose relative uncertainty
2577  // is the value of the ith RooAbsReal
2578  //
2579  // The integer "type" controls the type
2580  // of constraint term:
2581  //
2582  // type == 0 : NONE
2583  // type == 1 : Gaussian
2584  // type == 2 : Poisson
2585  // type == 3 : LogNormal
2586 
2587  RooArgList ConstraintTerms;
2588 
2589  RooArgList paramSet = paramHist.paramList();
2590 
2591  // Must get the full size of the TH1
2592  // (No direct method to do this...)
2593  Int_t numBins = uncertHist->GetNbinsX()*uncertHist->GetNbinsY()*uncertHist->GetNbinsZ();
2594  Int_t numParams = paramSet.getSize();
2595  // Int_t numBins = uncertHist->GetNbinsX()*uncertHist->GetNbinsY()*uncertHist->GetNbinsZ();
2596 
2597  // Check that there are N elements
2598  // in the RooArgList
2599  if( numBins != numParams ) {
2600  std::cout << "Error: In createStatConstraintTerms, encountered bad number of bins" << std::endl;
2601  std::cout << "Given histogram with " << numBins << " bins,"
2602  << " but require exactly " << numParams << std::endl;
2603  throw hf_exc();
2604  }
2605 
2606  Int_t TH1BinNumber = 0;
2607  for( Int_t i = 0; i < paramSet.getSize(); ++i) {
2608 
2609  TH1BinNumber++;
2610 
2611  while( uncertHist->IsBinUnderflow(TH1BinNumber) || uncertHist->IsBinOverflow(TH1BinNumber) ){
2612  TH1BinNumber++;
2613  }
2614 
2615  RooRealVar& gamma = (RooRealVar&) (paramSet[i]);
2616 
2617  std::cout << "Creating constraint for: " << gamma.GetName()
2618  << ". Type of constraint: " << type << std::endl;
2619 
2620  // Get the sigma from the hist
2621  // (the relative uncertainty)
2622  Double_t sigma = uncertHist->GetBinContent( TH1BinNumber );
2623 
2624  // If the sigma is <= 0,
2625  // do cont create the term
2626  if( sigma <= 0 ){
2627  std::cout << "Not creating constraint term for "
2628  << gamma.GetName()
2629  << " because sigma = " << sigma
2630  << " (sigma<=0)"
2631  << " (TH1 bin number = " << TH1BinNumber << ")"
2632  << std::endl;
2633  gamma.setConstant(kTRUE);
2634  continue;
2635  }
2636 
2637  // set reasonable ranges for gamma parameters
2638  gamma.setMax( 1 + 5*sigma );
2639  // gamma.setMin( TMath::Max(1. - 5*sigma, 0.) );
2640  gamma.setMin( 0. );
2641 
2642  // Make Constraint Term
2643  std::string constrName = string(gamma.GetName()) + "_constraint";
2644  std::string nomName = string("nom_") + gamma.GetName();
2645  std::string sigmaName = string(gamma.GetName()) + "_sigma";
2646  std::string poisMeanName = string(gamma.GetName()) + "_poisMean";
2647 
2648  if( type == Constraint::Gaussian ) {
2649 
2650  // Type 1 : RooGaussian
2651 
2652  // Make sigma
2653 
2654  RooConstVar constrSigma( sigmaName.c_str(), sigmaName.c_str(), sigma );
2655  //proto->import( constrSigma, RecycleConflictNodes() );
2656  //proto->import( constrSigma );
2657 
2658  // Make "observed" value
2659  RooRealVar constrNom(nomName.c_str(), nomName.c_str(), 1.0,0,10);
2660  constrNom.setConstant( true );
2661 
2662  // Make the constraint:
2663  RooGaussian gauss( constrName.c_str(), constrName.c_str(),
2664  constrNom, gamma, constrSigma );
2665 
2666  proto->import( gauss, RecycleConflictNodes() );
2667  //proto->import( gauss );
2668 
2669  } else if( type == Constraint::Poisson ) {
2670 
2671  Double_t tau = 1/sigma/sigma; // this is correct Poisson equivalent to a Gaussian with mean 1 and stdev sigma
2672 
2673  // Make nominal "observed" value
2674  RooRealVar constrNom(nomName.c_str(), nomName.c_str(), tau);
2675  constrNom.setMin(0);
2676  constrNom.setConstant( true );
2677 
2678  // Make the scaling term
2679  std::string scalingName = string(gamma.GetName()) + "_tau";
2680  RooConstVar poissonScaling( scalingName.c_str(), scalingName.c_str(), tau);
2681 
2682  // Make mean for scaled Poisson
2683  RooProduct constrMean( poisMeanName.c_str(), poisMeanName.c_str(), RooArgSet(gamma, poissonScaling) );
2684  //proto->import( constrSigma, RecycleConflictNodes() );
2685  //proto->import( constrSigma );
2686 
2687  // Type 2 : RooPoisson
2688  RooPoisson pois(constrName.c_str(), constrName.c_str(), constrNom, constrMean);
2689  pois.setNoRounding(true);
2690  proto->import( pois, RecycleConflictNodes() );
2691 
2692  } else {
2693 
2694  std::cout << "Error: Did not recognize Stat Error constraint term type: "
2695  << type << " for : " << paramHist.GetName() << std::endl;
2696  throw hf_exc();
2697  }
2698 
2699  // If the sigma value is less
2700  // than a supplied threshold,
2701  // set the variable to constant
2702  if( sigma < minSigma ) {
2703  std::cout << "Warning: Bin " << i << " = " << sigma
2704  << " and is < " << minSigma
2705  << ". Setting: " << gamma.GetName() << " to constant"
2706  << std::endl;
2707  gamma.setConstant(kTRUE);
2708  }
2709 
2710  constraintTermNames.push_back( constrName );
2711  ConstraintTerms.add( *proto->pdf(constrName.c_str()) );
2712 
2713  // Add the "observed" value to the
2714  // list of global observables:
2715  RooArgSet* globalSet = const_cast<RooArgSet*>(proto->set("globalObservables"));
2716 
2717  RooRealVar* nomVarInWorkspace = proto->var(nomName.c_str());
2718  if( ! globalSet->contains(*nomVarInWorkspace) ) {
2719  globalSet->add( *nomVarInWorkspace );
2720  }
2721 
2722  } // end loop over parameters
2723 
2724  return ConstraintTerms;
2725 
2726 }
2727 
2728 } // namespace RooStats
2729 } // namespace HistFactory
2730 
std::vector< RooStats::HistFactory::Asimov > & GetAsimovDatasets()
Definition: Measurement.h:79
static RooAbsData * GenerateAsimovData(const RooAbsPdf &pdf, const RooArgSet &observables)
virtual const char * GetName() const
Returns name of object.
Definition: TNamed.h:51
Bool_t saveSnapshot(const char *name, const char *paramNames)
Save snapshot of values and attributes (including "Constant") of parameters &#39;params&#39; If importValues ...
Poisson pdf.
Definition: RooPoisson.h:19
RooStats::HistFactory::StatError & GetStatError()
Definition: Sample.h:117
TIterator * createIterator(Bool_t dir=kIterForward) const
RooWorkspace * MakeCombinedModel(std::vector< std::string >, std::vector< RooWorkspace *>)
virtual Bool_t add(const RooAbsArg &var, Bool_t silent=kFALSE)
Add the specified argument to list.
float xmin
Definition: THbookFile.cxx:93
std::map< std::string, double > & GetNoSyst()
Definition: Measurement.h:124
ModelConfig is a simple class that holds configuration information specifying how a model should be u...
Definition: ModelConfig.h:52
Bool_t IsBinUnderflow(Int_t bin, Int_t axis=0) const
Return true if the bin is underflow.
Definition: TH1.cxx:4772
const RooArgSet * GetObservables() const
get RooArgSet for observables (return NULL if not existing)
Definition: ModelConfig.h:259
void AddMultiVarGaussConstraint(RooWorkspace *proto, std::string prefix, int lowBin, int highBin, std::vector< std::string > &likelihoodTermNames)
const RooArgList & paramList() const
Definition: ParamHistFunc.h:39
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
void Start(Bool_t reset=kTRUE)
Start the stopwatch.
Definition: TStopwatch.cxx:58
std::vector< RooStats::HistFactory::NormFactor > & GetNormFactorList()
Definition: Sample.h:109
void AddPoissonTerms(RooWorkspace *proto, std::string prefix, std::string obsPrefix, std::string expPrefix, int lowBin, int highBin, std::vector< std::string > &likelihoodTermNames)
Bool_t IsBinOverflow(Int_t bin, Int_t axis=0) const
Return true if the bin is overflow.
Definition: TH1.cxx:4740
static 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, std::map< std::string, double > noSyst)
virtual void SetWorkspace(RooWorkspace &ws)
Definition: ModelConfig.h:88
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
void GuessObsAndNuisance(const RooAbsData &data)
guesses Observables and ParametersOfInterest if not already set
Definition: ModelConfig.cxx:32
RooWorkspace * MakeSingleChannelWorkspace(Measurement &measurement, Channel &channel)
virtual Double_t GetBinContent(Int_t bin) const
Return content of bin number bin.
Definition: TH1.cxx:4638
RooDataSet * MergeDataSets(RooWorkspace *combined, std::vector< RooWorkspace *> wspace_vec, std::vector< std::string > channel_names, std::string dataSetName, RooArgList obsList, RooCategory *channelCat)
std::string GetName()
Definition: Data.h:34
void ConfigureWorkspace(RooWorkspace *)
Definition: Asimov.cxx:26
virtual Int_t GetNbinsZ() const
Definition: TH1.h:303
#define R__ASSERT(e)
Definition: TError.h:98
RooCmdArg Conditional(const RooArgSet &pdfSet, const RooArgSet &depSet, Bool_t depsAreCond=kFALSE)
std::vector< std::string > & GetPOIList()
Definition: Measurement.h:51
RooAbsReal * function(const char *name) const
Retrieve function (RooAbsReal) with given name. Note that all RooAbsPdfs are also RooAbsReals...
int Int_t
Definition: RtypesCore.h:41
int nbins[3]
RooAbsArg * arg(const char *name) const
Return RooAbsArg with given name. A null pointer is returned if none is found.
STL namespace.
Class RooRealSumPdf implements a PDF constructed from a sum of functions:
Definition: RooRealSumPdf.h:24
RooCmdArg RecycleConflictNodes(Bool_t flag=kTRUE)
virtual void SetObservables(const RooArgSet &set)
specify the observables
Definition: ModelConfig.h:156
virtual Bool_t Multiply(TF1 *h1, Double_t c1=1)
Performs the operation: this = this*c1*f1 if errors are defined (see TH1::Sumw2), errors are also rec...
Definition: TH1.cxx:5363
std::vector< RooStats::HistFactory::HistoSys > & GetHistoSysList()
Definition: Sample.h:111
void setMax(const char *name, Double_t value)
Set maximum of name range to given value.
Definition: RooRealVar.cxx:418
void Reset()
Definition: TCollection.h:161
static void ConfigureWorkspaceForMeasurement(const std::string &ModelName, RooWorkspace *ws_single, Measurement &measurement)
std::string AddNormFactor(RooWorkspace *proto, std::string &channel, std::string &sigmaEpsilon, Sample &sample, bool doRatio)
virtual void Reset(Option_t *option="")
Reset this histogram: contents, errors, etc.
Definition: TH1.cxx:6418
double beta(double x, double y)
Calculates the beta function.
virtual Bool_t setLabel(const char *label, Bool_t printError=kTRUE)
Set value by specifying the name of the desired state If printError is set, a message will be printed...
const char * Name
Definition: TXMLSetup.cxx:67
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...
std::vector< std::string > GetPreprocessFunctions()
void SetFunctionsToPreprocess(std::vector< std::string > lines)
double sqrt(double)
Double_t GetXmin() const
Definition: TAxis.h:139
RooCategory & method1D()
virtual const char * ClassName() const
Returns name of class to which the object belongs.
Definition: TObject.cxx:188
RooDataSet is a container class to hold N-dimensional binned data.
Definition: RooDataHist.h:40
static TString Format(const char *fmt,...)
Static method which formats a string using a printf style format descriptor and return a TString...
Definition: TString.cxx:2335
you should not use this method at all Int_t Int_t Double_t Double_t Double_t Int_t Double_t Double_t Double_t tau
Definition: TRolke.cxx:630
RooHistFunc implements a real-valued function sampled from a multidimensional histogram.
Definition: RooHistFunc.h:29
void ProcessExpectedHisto(TH1 *hist, RooWorkspace *proto, std::string prefix, std::string productPrefix, std::string systTerm)
Plain Gaussian p.d.f.
Definition: RooGaussian.h:25
double pow(double, double)
std::vector< RooStats::HistFactory::Channel > & GetChannels()
Definition: Measurement.h:105
virtual Double_t GetBinCenter(Int_t bin) const
Return center of bin.
Definition: TAxis.cxx:464
virtual void Print(Option_t *options=0) const
Print TNamed name and title.
Definition: RooAbsArg.h:227
virtual void SetPdf(const RooAbsPdf &pdf)
Set the Pdf, add to the the workspace if not already there.
Definition: ModelConfig.h:97
RooConstVar represent a constant real-valued object.
Definition: RooConstVar.h:25
const Double_t sigma
void AddConstraintTerms(RooWorkspace *proto, Measurement &measurement, std::string prefix, std::string interpName, std::vector< OverallSys > &systList, std::vector< std::string > &likelihoodTermNames, std::vector< std::string > &totSystTermNames)
std::string GetName()
Definition: Sample.h:82
virtual void Print(Option_t *options=0) const
This method must be overridden when a class wants to print itself.
static RooArgList createParamSet(RooWorkspace &w, const std::string &, const RooArgList &Vars)
std::map< std::string, double > & GetUniformSyst()
Definition: Measurement.h:122
RooNumIntConfig * specialIntegratorConfig() const
Returns the specialized integrator configuration for this RooAbsReal.
std::string GetName()
Definition: Asimov.h:32
void LinInterpWithConstraint(RooWorkspace *proto, TH1 *nominal, std::vector< HistoSys >, std::string prefix, std::string productPrefix, std::string systTerm, std::vector< std::string > &likelihoodTermNames)
RooRealVar represents a fundamental (non-derived) real valued object.
Definition: RooRealVar.h:37
RooAbsData * data(const char *name) const
Retrieve dataset (binned or unbinned) with given name. A null pointer is returned if not found...
virtual void setVal(Double_t value)
Set value of variable to &#39;value&#39;.
Definition: RooRealVar.cxx:205
#define alpha_High
double gamma(double x)
Int_t getSize() const
Class to manage histogram axis.
Definition: TAxis.h:36
void setConstant(Bool_t value=kTRUE)
void setMin(const char *name, Double_t value)
Set minimum of name range to given value.
Definition: RooRealVar.cxx:388
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:158
virtual void SetBinContent(Int_t bin, Double_t content)
Set bin content see convention for numbering bins in TH1::GetBin In case the bin number is greater th...
Definition: TH1.cxx:8323
std::vector< RooStats::HistFactory::ShapeSys > & GetShapeSysList()
Definition: Sample.h:114
virtual void Error(const char *method, const char *msgfmt,...) const
Issue error message.
Definition: TObject.cxx:925
char * Form(const char *fmt,...)
TAxis * GetYaxis()
Definition: TH1.h:325
RooAbsData is the common abstract base class for binned and unbinned datasets.
Definition: RooAbsData.h:37
float xmax
Definition: THbookFile.cxx:93
Bool_t loadSnapshot(const char *name)
Load the values and attributes of the parameters in the snapshot saved with the given name...
void setNoRounding(bool flag=kTRUE)
Definition: RooPoisson.h:33
RooCmdArg Rename(const char *suffix)
RooArgSet allVars() const
Return set with all variable objects.
void SetObsToExpected(RooWorkspace *proto, std::string obsPrefix, std::string expPrefix, int lowBin, int highBin)
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:25
std::vector< RooStats::HistFactory::OverallSys > & GetOverallSysList()
Definition: Sample.h:108
RooProduct a RooAbsReal implementation that represent the product of a given set of other RooAbsReal ...
Definition: RooProduct.h:32
virtual void add(const RooArgSet &row, Double_t weight=1.0, Double_t weightError=0)
Add a data point, with its coordinates specified in the &#39;data&#39; argset, to the data set...
TObject * obj(const char *name) const
Return any type of object (RooAbsArg, RooAbsData or generic object) with given name) ...
std::vector< std::string > & GetConstantParams()
Definition: Measurement.h:60
RooCmdArg Import(const char *state, TH1 &histo)
virtual void SetName(const char *name)
Change the name of this histogram.
Definition: TH1.cxx:8071
RooWorkspace * MakeSingleChannelModel(Measurement &measurement, Channel &channel)
Multivariate Gaussian p.d.f.
RooCmdArg Index(RooCategory &icat)
Namespace for the RooStats classes.
Definition: Asimov.h:20
#define alpha_Low
#define ClassImp(name)
Definition: Rtypes.h:279
HistFactory::StatErrorConfig & GetStatErrorConfig()
Definition: Channel.h:67
Constraint::Type GetConstraintType()
Definition: Systematics.h:242
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
Bool_t removeSet(const char *name)
Remove a named set from the workspace.
void append(RooDataSet &data)
Add all data points of given data set to this data set.
int type
Definition: TGX11.cxx:120
RooCmdArg WeightVar(const char *name, Bool_t reinterpretAsWeight=kFALSE)
std::map< std::string, double > & GetLogNormSyst()
Definition: Measurement.h:123
The TH1 histogram class.
Definition: TH1.h:80
std::vector< RooStats::HistFactory::Sample > & GetSamples()
Definition: Channel.h:71
static void PrintCovarianceMatrix(RooFitResult *result, RooArgSet *params, std::string filename)
RooRealVar * var(const char *name) const
Retrieve real-valued variable (RooRealVar) with given name. A null pointer is returned if not found...
std::vector< RooStats::HistFactory::Data > & GetAdditionalData()
Definition: Channel.h:60
TAxis * GetZaxis()
Definition: TH1.h:326
RooFactoryWSTool & factory()
Return instance to factory tool.
TH1 * MakeScaledUncertaintyHist(const std::string &Name, std::vector< std::pair< TH1 *, TH1 *> > HistVec)
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.
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 * Clone(const char *newname=0) const
Make a complete copy of the underlying object.
Definition: TH1.cxx:2544
void setPositiveDefinite(bool flag=true)
#define NULL
Definition: Rtypes.h:82
const char * proto
Definition: civetweb.c:11652
Definition: tree.py:1
virtual void SetParametersOfInterest(const RooArgSet &set)
Definition: ModelConfig.h:115
double result[121]
Bool_t contains(const RooAbsArg &var) const
virtual Int_t GetNbinsX() const
Definition: TH1.h:301
Int_t GetNbins() const
Definition: TAxis.h:127
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:182
virtual const char * GetName() const
Returns name of object.
Definition: TObject.cxx:416
double exp(double)
TH1 * MakeAbsolUncertaintyHist(const std::string &Name, const TH1 *Hist)
const Bool_t kTRUE
Definition: Rtypes.h:91
void Print(Option_t *opts=0) const
Print contents of the workspace.
RooStats::HistFactory::Data & GetData()
Definition: Channel.h:55
RooSimultaneous facilitates simultaneous fitting of multiple PDFs to subsets of a given dataset...
double norm(double *x, double *p)
Definition: unuranDistr.cxx:40
Double_t GetXmax() const
Definition: TAxis.h:140
void ConfigureHistFactoryDataset(RooDataSet *obsData, TH1 *nominal, RooWorkspace *proto, std::vector< std::string > obsNameVec)
void MakeTotalExpected(RooWorkspace *proto, std::string totName, std::vector< std::string > &syst_x_expectedPrefixNames, std::vector< std::string > &normByNames)
RooArgList createStatConstraintTerms(RooWorkspace *proto, std::vector< std::string > &constraintTerms, ParamHistFunc &paramHist, TH1 *uncertHist, Constraint::Type type, Double_t minSigma)
Double_t correlation(const RooAbsArg &par1, const RooAbsArg &par2) const
Definition: RooFitResult.h:116
char name[80]
Definition: TGX11.cxx:109
std::vector< RooStats::HistFactory::ShapeFactor > & GetShapeFactorList()
Definition: Sample.h:115
Bool_t isConstant() const
Definition: RooAbsArg.h:266
TAxis * GetXaxis()
Definition: TH1.h:324
The RooWorkspace is a persistable container for RooFit projects.
Definition: RooWorkspace.h:42
virtual Int_t GetNbinsY() const
Definition: TH1.h:302
std::map< std::string, double > & GetGammaSyst()
Definition: Measurement.h:121
virtual Double_t GetBinError(Int_t bin) const
Return value of error associated to bin number bin.
Definition: TH1.cxx:8173
Stopwatch class.
Definition: TStopwatch.h:30