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