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