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