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 "RooProfileLL.h"
42#include "RooFitResult.h"
43#include "RooDataHist.h"
44#include "RooHistFunc.h"
45#include "RooHistPdf.h"
46#include "RooRealSumPdf.h"
47#include "RooWorkspace.h"
48#include "RooCustomizer.h"
49#include "RooPlot.h"
50#include "RooHelpers.h"
51#include "RooBinning.h"
52#include "RooBinWidthFunction.h"
58#include "HFMsgService.h"
59
60#include "TH1.h"
61#include "TStopwatch.h"
62#include "TVectorD.h"
63#include "TMatrixDSym.h"
64
65// specific to this package
71
72#include <algorithm>
73#include <memory>
74#include <utility>
75
76#define VERBOSE
77
78#define alpha_Low "-5"
79#define alpha_High "5"
80#define NoHistConst_Low "0"
81#define NoHistConst_High "2000"
82
83// use this order for safety on library loading
84using namespace RooFit ;
85using namespace RooStats ;
86using namespace std ;
87
89
90namespace RooStats{
91namespace HistFactory{
92
95
97 Configuration const& cfg) :
98 fSystToFix( measurement.GetConstantParams() ),
99 fParamValues( measurement.GetParamValues() ),
100 fNomLumi( measurement.GetLumi() ),
101 fLumiError( measurement.GetLumi()*measurement.GetLumiRelErr() ),
102 fLowBin( measurement.GetBinLow() ),
103 fHighBin( measurement.GetBinHigh() ),
104 fCfg{cfg} {
105
106 // Set Preprocess functions
108
109 }
110
111 void HistoToWorkspaceFactoryFast::ConfigureWorkspaceForMeasurement( const std::string& ModelName, RooWorkspace* ws_single, Measurement& measurement ) {
112
113 // Configure a workspace by doing any
114 // necessary post-processing and by
115 // creating a ModelConfig
116
117 // Make a ModelConfig and configure it
118 ModelConfig * proto_config = (ModelConfig *) ws_single->obj("ModelConfig");
119 if( proto_config == nullptr ) {
120 std::cout << "Error: Did not find 'ModelConfig' object in file: " << ws_single->GetName()
121 << std::endl;
122 throw hf_exc();
123 }
124
125 std::vector<std::string> poi_list = measurement.GetPOIList();
126 if( poi_list.empty() ) {
127 cxcoutWHF << "No Parametetrs of interest are set" << std::endl;
128 }
129
130
131 std::stringstream sstream;
132 sstream << "Setting Parameter(s) of Interest as: ";
133 for(unsigned int i = 0; i < poi_list.size(); ++i) {
134 sstream << poi_list.at(i) << " ";
135 }
136 cxcoutIHF << sstream.str() << endl;
137
138 RooArgSet params;
139 for( unsigned int i = 0; i < poi_list.size(); ++i ) {
140 std::string poi_name = poi_list.at(i);
141 RooRealVar* poi = (RooRealVar*) ws_single->var(poi_name);
142 if(poi){
143 params.add(*poi);
144 }
145 else {
146 std::cout << "WARNING: Can't find parameter of interest: " << poi_name
147 << " in Workspace. Not setting in ModelConfig." << std::endl;
148 //throw hf_exc();
149 }
150 }
151 proto_config->SetParametersOfInterest(params);
152
153 // Name of an 'edited' model, if necessary
154 std::string NewModelName = "newSimPdf"; // <- This name is hard-coded in HistoToWorkspaceFactoryFast::EditSyt. Probably should be changed to : std::string("new") + ModelName;
155
156 // Set the ModelConfig's Params of Interest
157 RooAbsData* expData = ws_single->data("asimovData");
158 if( !expData ) {
159 std::cout << "Error: Failed to find dataset: " << expData
160 << " in workspace" << std::endl;
161 throw hf_exc();
162 }
163 if(poi_list.size()!=0){
164 proto_config->GuessObsAndNuisance(*expData, RooMsgService::instance().isActive(static_cast<TObject*>(nullptr), RooFit::HistFactory, RooFit::INFO));
165 }
166
167 // Now, let's loop over any additional asimov datasets
168 // that we need to make
169
170 // Get the pdf
171 // Notice that we get the "new" pdf, this is the one that is
172 // used in the creation of these asimov datasets since they
173 // are fitted (or may be, at least).
174 RooAbsPdf* pdf = ws_single->pdf(NewModelName);
175 if( !pdf ) pdf = ws_single->pdf( ModelName );
176 const RooArgSet* observables = ws_single->set("observables");
177
178 // Create a SnapShot of the nominal values
179 std::string SnapShotName = "NominalParamValues";
180 ws_single->saveSnapshot(SnapShotName.c_str(), ws_single->allVars());
181
182 for( unsigned int i=0; i<measurement.GetAsimovDatasets().size(); ++i) {
183
184 // Set the variable values and "const" ness with the workspace
185 RooStats::HistFactory::Asimov& asimov = measurement.GetAsimovDatasets().at(i);
186 std::string AsimovName = asimov.GetName();
187
188 cxcoutPHF << "Generating additional Asimov Dataset: " << AsimovName << std::endl;
189 asimov.ConfigureWorkspace(ws_single);
190 std::unique_ptr<RooDataSet> asimov_dataset{static_cast<RooDataSet*>(AsymptoticCalculator::GenerateAsimovData(*pdf, *observables))};
191
192 cxcoutPHF << "Importing Asimov dataset" << std::endl;
193 bool failure = ws_single->import(*asimov_dataset, Rename(AsimovName.c_str()));
194 if( failure ) {
195 std::cout << "Error: Failed to import Asimov dataset: " << AsimovName
196 << std::endl;
197 throw hf_exc();
198 }
199
200 // Load the snapshot at the end of every loop iteration
201 // so we start each loop with a "clean" snapshot
202 ws_single->loadSnapshot(SnapShotName.c_str());
203 }
204
205 // Cool, we're done
206 return; // ws_single;
207 }
208
209
210 // We want to eliminate this interface and use the measurment directly
212
213 // This is a pretty light-weight wrapper function
214 //
215 // Take a fully configured measurement as well as
216 // one of its channels
217 //
218 // Return a workspace representing that channel
219 // Do this by first creating a vector of EstimateSummary's
220 // and this by configuring the workspace with any post-processing
221
222 // Get the channel's name
223 string ch_name = channel.GetName();
224
225 // Create a workspace for a SingleChannel from the Measurement Object
226 RooWorkspace* ws_single = this->MakeSingleChannelWorkspace(measurement, channel);
227 if( ws_single == nullptr ) {
228 cxcoutF(HistFactory) << "Error: Failed to make Single-Channel workspace for channel: " << ch_name
229 << " and measurement: " << measurement.GetName() << std::endl;
230 throw hf_exc();
231 }
232
233 // Finally, configure that workspace based on
234 // properties of the measurement
235 HistoToWorkspaceFactoryFast::ConfigureWorkspaceForMeasurement( "model_"+ch_name, ws_single, measurement );
236
237 return ws_single;
238
239 }
240
242
243 // This function takes a fully configured measurement
244 // which may contain several channels and returns
245 // a workspace holding the combined model
246 //
247 // This can be used, for example, within a script to produce
248 // a combined workspace on-the-fly
249 //
250 // This is a static function (for now) to make
251 // it a one-liner
252
254
255 // First, we create an instance of a HistFactory
256 HistoToWorkspaceFactoryFast factory( measurement );
257
258 // Loop over the channels and create the individual workspaces
259 vector<std::unique_ptr<RooWorkspace>> channel_workspaces;
260 vector<string> channel_names;
261
262 for( unsigned int chanItr = 0; chanItr < measurement.GetChannels().size(); ++chanItr ) {
263
264 HistFactory::Channel& channel = measurement.GetChannels().at( chanItr );
265
266 if( ! channel.CheckHistograms() ) {
267 cxcoutFHF << "MakeModelAndMeasurementsFast: Channel: " << channel.GetName()
268 << " has uninitialized histogram pointers" << std::endl;
269 throw hf_exc();
270 }
271
272 string ch_name = channel.GetName();
273 channel_names.push_back(ch_name);
274
275 // GHL: Renaming to 'MakeSingleChannelWorkspace'
276 RooWorkspace* ws_single = factory.MakeSingleChannelModel( measurement, channel );
277
278 channel_workspaces.emplace_back(ws_single);
279
280 }
281
282
283 // Now, combine the individual channel workspaces to
284 // form the combined workspace
285 RooWorkspace* ws = factory.MakeCombinedModel( channel_names, channel_workspaces );
286
287
288 // Configure the workspace
290
291 // Done. Return the pointer
292 return ws;
293
294 }
295
296/// Create observables of type RooRealVar. Creates 1 to 3 observables, depending on the type of the histogram.
298 RooArgList observables;
299
300 for (unsigned int idx=0; idx < fObsNameVec.size(); ++idx) {
301 if (!proto->var(fObsNameVec[idx])) {
302 const TAxis *axis = (idx == 0) ? hist->GetXaxis() : (idx == 1 ? hist->GetYaxis() : hist->GetZaxis());
303 Int_t nbins = axis->GetNbins();
304 double xmin = axis->GetXmin();
305 double xmax = axis->GetXmax();
306 // create observable
307 auto obs = static_cast<RooRealVar*>(proto->factory(
308 Form("%s[%f,%f]", fObsNameVec[idx].c_str(), xmin, xmax)));
309 if(strlen(axis->GetTitle())>0) obs->SetTitle(axis->GetTitle());
310 obs->setBins(nbins);
311 if (axis->IsVariableBinSize()) {
312 RooBinning binning(nbins, axis->GetXbins()->GetArray());
313 obs->setBinning(binning);
314 }
315 }
316
317 observables.add(*proto->var(fObsNameVec[idx]));
318 }
319
320 return observables;
321}
322
323 /// Create the nominal hist function from `hist`, and register it in the workspace.
325 const RooArgList& observables) const {
326 if(hist) {
327 cxcoutI(HistFactory) << "processing hist " << hist->GetName() << endl;
328 } else {
329 cxcoutF(HistFactory) << "hist is empty" << endl;
330 R__ASSERT(hist != 0);
331 return nullptr;
332 }
333
334 // determine histogram dimensionality
335 unsigned int histndim(1);
336 std::string classname = hist->ClassName();
337 if (classname.find("TH1")==0) { histndim=1; }
338 else if (classname.find("TH2")==0) { histndim=2; }
339 else if (classname.find("TH3")==0) { histndim=3; }
340 R__ASSERT( histndim==fObsNameVec.size() );
341
342 prefix += "_Hist_alphanominal";
343
344 RooDataHist histDHist((prefix + "DHist").c_str(),"",observables,hist);
345 RooHistFunc histFunc(prefix.c_str(),"",observables,histDHist,0);
346
347 proto->import(histFunc, RecycleConflictNodes());
348 auto histFuncInWS = static_cast<RooHistFunc*>(proto->arg(prefix.c_str()));
349
350 return histFuncInWS;
351 }
352
353 namespace {
354
355 void makeGaussianConstraint(RooAbsArg& param, RooWorkspace& proto, bool isUniform,
356 std::vector<std::string> & constraintTermNames) {
357 std::string paramName = param.GetName();
358 std::string constraintName = paramName + "Constraint";
359
360 // do nothing if the constraint term already exists
361 if(proto.pdf(constraintName)) return;
362
363 // case systematic is uniform (asssume they are like a Gaussian but with
364 // a large width (100 instead of 1)
365 const double gaussSigma = isUniform ? 100. : 1.0;
366 if (isUniform) {
367 cxcoutIHF << "Added a uniform constraint for " << paramName << " as a Gaussian constraint with a very large sigma " << std::endl;
368 }
369
370 std::stringstream command;
371 command << "Gaussian::" << constraintName << "(" << paramName << ",nom_" << paramName << "[0.,-10,10],"
372 << gaussSigma << ")";
373 constraintTermNames.emplace_back(proto.factory(command.str())->GetName());
374 auto * normParam = proto.var(std::string("nom_") + paramName);
375 normParam->setConstant();
376 const_cast<RooArgSet*>(proto.set("globalObservables"))->add(*normParam);
377 }
378
379 /// Make list of abstract parameters that interpolate in space of variations.
380 RooArgList makeInterpolationParameters(std::vector<HistoSys> const& histoSysList, RooWorkspace& proto) {
381 RooArgList params( ("alpha_Hist") );
382
383 // range is set using defined macro (see top of the page)
384 string range=string("[")+alpha_Low+","+alpha_High+"]";
385
386 for(auto const& histoSys : histoSysList) {
387 const std::string histoSysName = histoSys.GetName();
388 RooRealVar* temp = proto.var("alpha_" + histoSysName);
389 if(!temp){
390 temp = (RooRealVar*) proto.factory("alpha_" + histoSysName + range);
391 }
392 params.add(* temp );
393 }
394
395 return params;
396 }
397
398 /// Create a linear interpolation object that holds nominal and systematics, import it into the workspace,
399 /// and return a pointer to it.
400 RooAbsArg* makeLinInterp(RooArgList const& interpolationParams,
401 RooHistFunc* nominalFunc,
402 RooWorkspace* proto, const std::vector<HistoSys>& histoSysList,
403 const string& prefix,
404 const RooArgList& observables) {
405
406 // now make function that linearly interpolates expectation between variations
407 // get low/high variations to interpolate between
408 vector<double> low, high;
409 RooArgSet lowSet, highSet;
410 //ES// for(unsigned int j=0; j<lowHist.size(); ++j){
411 for(unsigned int j=0; j<histoSysList.size(); ++j){
412 std::stringstream str;
413 str<<"_"<<j;
414
415 const HistoSys& histoSys = histoSysList.at(j);
416 RooDataHist* lowDHist = new RooDataHist((prefix+str.str()+"lowDHist").c_str(),"",observables, histoSys.GetHistoLow());
417 RooDataHist* highDHist = new RooDataHist((prefix+str.str()+"highDHist").c_str(),"",observables, histoSys.GetHistoHigh());
418 RooHistFunc* lowFunc = new RooHistFunc((prefix+str.str()+"low").c_str(),"",observables,*lowDHist,0) ;
419 RooHistFunc* highFunc = new RooHistFunc((prefix+str.str()+"high").c_str(),"",observables,*highDHist,0) ;
420 lowSet.add(*lowFunc);
421 highSet.add(*highFunc);
422 }
423
424 // this is sigma(params), a piece-wise linear interpolation
425 PiecewiseInterpolation interp(prefix.c_str(),"",*nominalFunc,lowSet,highSet,interpolationParams);
426 interp.setPositiveDefinite();
427 interp.setAllInterpCodes(4); // LM: change to 4 (piece-wise linear to 6th order polynomial interpolation + linear extrapolation )
428 // KC: interpo codes 1 etc. don't have proper analytic integral.
429 RooArgSet observableSet(observables);
430 interp.setBinIntegrator(observableSet);
431 interp.forceNumInt();
432
433 proto->import(interp, RecycleConflictNodes()); // individual params have already been imported in first loop of this function
434
435 auto interpInWS = proto->arg(prefix.c_str());
436 assert(interpInWS);
437
438 return interpInWS;
439 }
440
441 }
442
443 // GHL: Consider passing the NormFactor list instead of the entire sample
444 std::unique_ptr<RooProduct> HistoToWorkspaceFactoryFast::CreateNormFactor(RooWorkspace* proto, string& channel, string& sigmaEpsilon, Sample& sample, bool doRatio){
445
446 std::vector<string> prodNames;
447
448 vector<NormFactor> normList = sample.GetNormFactorList();
449 vector<string> normFactorNames, rangeNames;
450
451
452 string overallNorm_times_sigmaEpsilon = sample.GetName() + "_" + channel + "_scaleFactors";
453 auto sigEps = proto->arg(sigmaEpsilon.c_str());
454 assert(sigEps);
455 auto normFactor = std::make_unique<RooProduct>(overallNorm_times_sigmaEpsilon.c_str(), overallNorm_times_sigmaEpsilon.c_str(), RooArgList(*sigEps));
456
457 if(normList.size() > 0){
458
459 for(vector<NormFactor>::iterator itr = normList.begin(); itr != normList.end(); ++itr){
460
461 NormFactor& norm = *itr;
462
463 string varname;
464 if(doRatio) {
465 varname = norm.GetName() + "_" + channel;
466 }
467 else {
468 varname=norm.GetName();
469 }
470
471 // GHL: Check that the NormFactor doesn't already exist
472 // (it may have been created as a function expression
473 // during preprocessing)
474 std::stringstream range;
475 range << "[" << norm.GetVal() << "," << norm.GetLow() << "," << norm.GetHigh() << "]";
476
477 if( proto->obj(varname) == nullptr) {
478 cxcoutI(HistFactory) << "making normFactor: " << norm.GetName() << endl;
479 // remove "doRatio" and name can be changed when ws gets imported to the combined model.
480 proto->factory(varname + range.str());
481 }
482
483 prodNames.push_back(varname);
484 rangeNames.push_back(range.str());
485 normFactorNames.push_back(varname);
486 }
487
488
489 for (const auto& name : prodNames) {
490 auto arg = proto->arg(name.c_str());
491 assert(arg);
492 normFactor->addTerm(arg);
493 }
494
495 }
496
497 unsigned int rangeIndex=0;
498 for( vector<string>::iterator nit = normFactorNames.begin(); nit!=normFactorNames.end(); ++nit){
499 if( count (normFactorNames.begin(), normFactorNames.end(), *nit) > 1 ){
500 cxcoutI(HistFactory) <<"<NormFactor Name =\""<<*nit<<"\"> is duplicated for <Sample Name=\""
501 << sample.GetName() << "\">, but only one factor will be included. \n Instead, define something like"
502 << "\n\t<Function Name=\""<<*nit<<"Squared\" Expression=\""<<*nit<<"*"<<*nit<<"\" Var=\""<<*nit<<rangeNames.at(rangeIndex)
503 << "\"> \nin your top-level XML's <Measurment> entry and use <NormFactor Name=\""<<*nit<<"Squared\" in your channel XML file."<< endl;
504 }
505 ++rangeIndex;
506 }
507
508 return normFactor;
509 }
510
512 string interpName,
513 std::vector<OverallSys>& systList,
514 vector<string>& constraintTermNames,
515 vector<string>& totSystTermNames) {
516
517 // add variables for all the relative overall uncertainties we expect
518 // range is set using defined macro (see top of the page)
519
520 string range=string("[0,")+alpha_Low+","+alpha_High+"]";
521 totSystTermNames.push_back(prefix);
522
523 RooArgSet params(prefix.c_str());
524 vector<double> lowVec, highVec;
525
526 std::map<std::string, double>::iterator itconstr;
527 for(unsigned int i = 0; i < systList.size(); ++i) {
528
529 OverallSys& sys = systList.at(i);
530 std::string strname = sys.GetName();
531 const char * name = strname.c_str();
532
533 // case of no systematic (is it possible)
534 if (meas.GetNoSyst().count(sys.GetName()) > 0 ) {
535 cxcoutI(HistFactory) << "HistoToWorkspaceFast::AddConstraintTerm - skip systematic " << sys.GetName() << std::endl;
536 continue;
537 }
538 // case systematic is a gamma constraint
539 if (meas.GetGammaSyst().count(sys.GetName()) > 0 ) {
540 double relerr = meas.GetGammaSyst().find(sys.GetName() )->second;
541 if (relerr <= 0) {
542 cxcoutI(HistFactory) << "HistoToWorkspaceFast::AddConstraintTerm - zero uncertainty assigned - skip systematic " << sys.GetName() << std::endl;
543 continue;
544 }
545 double tauVal = 1./(relerr*relerr);
546 double sqtau = 1./relerr;
547 RooAbsArg * beta = proto->factory(TString::Format("beta_%s[1,0,10]",name) );
548 // the global observable (y_s)
549 RooAbsArg * yvar = proto->factory(TString::Format("nom_%s[%f,0,10]",beta->GetName(),tauVal)) ;
550 // the rate of the gamma distribution (theta)
551 RooAbsArg * theta = proto->factory(TString::Format("theta_%s[%f]",name,1./tauVal));
552 // find alpha as function of beta
553 RooAbsArg* alphaOfBeta = proto->factory(TString::Format("PolyVar::alphaOfBeta_%s(beta_%s,{%f,%f})",name,name,-sqtau,sqtau));
554
555 // add now the constraint itself Gamma_beta_constraint(beta, y+1, tau, 0 )
556 // build the gamma parameter k = as y_s + 1
557 RooAbsArg * kappa = proto->factory(TString::Format("sum::k_%s(%s,1.)",name,yvar->GetName()) );
558 RooAbsArg * gamma = proto->factory(TString::Format("Gamma::%sConstraint(%s, %s, %s, 0.0)",beta->GetName(),beta->GetName(), kappa->GetName(), theta->GetName() ) );
560 alphaOfBeta->Print("t");
561 gamma->Print("t");
562 }
563 constraintTermNames.push_back(gamma->GetName());
564 // set global observables
565 RooRealVar * gobs = dynamic_cast<RooRealVar*>(yvar); assert(gobs);
566 gobs->setConstant(true);
567 const_cast<RooArgSet*>(proto->set("globalObservables"))->add(*yvar);
568
569 // add alphaOfBeta in the list of params to interpolate
570 params.add(*alphaOfBeta);
571 cxcoutI(HistFactory) << "Added a gamma constraint for " << name << std::endl;
572
573 }
574 else {
575 RooRealVar* alpha = (RooRealVar*) proto->var(prefix + sys.GetName());
576 if(!alpha) {
577 alpha = (RooRealVar*) proto->factory(prefix + sys.GetName() + range);
578 }
579 // add the Gaussian constraint part
580 const bool isUniform = meas.GetUniformSyst().count(sys.GetName()) > 0;
581 makeGaussianConstraint(*alpha, *proto, isUniform, constraintTermNames);
582
583 // check if exists a log-normal constraint
584 if (meas.GetLogNormSyst().count(sys.GetName()) == 0 && meas.GetGammaSyst().count(sys.GetName()) == 0 ) {
585 // just add the alpha for the parameters of the FlexibleInterpVar function
586 params.add(*alpha);
587 }
588 // case systematic is a log-normal constraint
589 if (meas.GetLogNormSyst().count(sys.GetName()) > 0 ) {
590 // log normal constraint for parameter
591 double relerr = meas.GetLogNormSyst().find(sys.GetName() )->second;
592 double tauVal = 1./relerr;
593 std::string tauName = "tau_" + sys.GetName();
594 proto->factory(TString::Format("%s[%f]",tauName.c_str(),tauVal ) );
595 double kappaVal = 1. + relerr;
596 std::string kappaName = "kappa_" + sys.GetName();
597 proto->factory(TString::Format("%s[%f]",kappaName.c_str(),kappaVal ) );
598 const char * alphaName = alpha->GetName();
599
600 std::string alphaOfBetaName = "alphaOfBeta_" + sys.GetName();
601 RooAbsArg * alphaOfBeta = proto->factory(TString::Format("expr::%s('%s*(pow(%s,%s)-1.)',%s,%s,%s)",alphaOfBetaName.c_str(),
602 tauName.c_str(),kappaName.c_str(),alphaName,
603 tauName.c_str(),kappaName.c_str(),alphaName ) );
604
605 cxcoutI(HistFactory) << "Added a log-normal constraint for " << name << std::endl;
606 if (RooMsgService::instance().isActive(static_cast<TObject*>(nullptr), RooFit::HistFactory, RooFit::DEBUG))
607 alphaOfBeta->Print("t");
608 params.add(*alphaOfBeta);
609 }
610
611 }
612 // add low/high vectors
613 double low = sys.GetLow();
614 double high = sys.GetHigh();
615 lowVec.push_back(low);
616 highVec.push_back(high);
617
618 } // end sys loop
619
620 if(systList.size() > 0){
621 // this is epsilon(alpha_j), a piece-wise linear interpolation
622 // LinInterpVar interp( (interpName).c_str(), "", params, 1., lowVec, highVec);
623
624 assert(!params.empty());
625 assert(int(lowVec.size()) == params.getSize() );
626
627 FlexibleInterpVar interp( (interpName).c_str(), "", params, 1., lowVec, highVec);
628 interp.setAllInterpCodes(4); // LM: change to 4 (piece-wise linear to 6th order polynomial interpolation + linear extrapolation )
629 //interp.setAllInterpCodes(0); // simple linear interpolation
630 proto->import(interp); // params have already been imported in first loop of this function
631 } else{
632 // some strange behavior if params,lowVec,highVec are empty.
633 //cout << "WARNING: No OverallSyst terms" << endl;
634 RooConstVar interp( (interpName).c_str(), "", 1.);
635 proto->import(interp); // params have already been imported in first loop of this function
636 }
637
638 // std::cout << "after creating FlexibleInterpVar " << std::endl;
639 // proto->Print();
640
641 }
642
643
645 const vector<RooProduct*>& sampleScaleFactors, std::vector<vector<RooAbsArg*>>& sampleHistFuncs) const {
646 assert(sampleScaleFactors.size() == sampleHistFuncs.size());
647
648 // for ith bin calculate totN_i = lumi * sum_j expected_j * syst_j
649
650 if (fObsNameVec.empty() && !fObsName.empty())
651 throw std::logic_error("HistFactory didn't process the observables correctly. Please file a bug report.");
652
653 auto firstHistFunc = dynamic_cast<const RooHistFunc*>(sampleHistFuncs.front().front());
654 if (!firstHistFunc) {
655 auto piecewiseInt = dynamic_cast<const PiecewiseInterpolation*>(sampleHistFuncs.front().front());
656 firstHistFunc = dynamic_cast<const RooHistFunc*>(piecewiseInt->nominalHist());
657 }
658 assert(firstHistFunc);
659
660 // Prepare a function to divide all bin contents by bin width to get a density:
661 const std::string binWidthFunctionName = totName + "_binWidth";
662 RooBinWidthFunction binWidth(binWidthFunctionName.c_str(), "Divide by bin width to obtain probability density", *firstHistFunc, true);
663 proto->import(binWidth);
664 auto binWidthWS = proto->function(binWidthFunctionName.c_str());
665 assert(binWidthWS);
666
667 // Loop through samples and create products of their functions:
668 RooArgSet coefList;
669 RooArgSet shapeList;
670 for (unsigned int i=0; i < sampleHistFuncs.size(); ++i) {
671 assert(!sampleHistFuncs[i].empty());
672 coefList.add(*sampleScaleFactors[i]);
673
674 std::vector<RooAbsArg*>& thisSampleHistFuncs = sampleHistFuncs[i];
675 thisSampleHistFuncs.push_back(binWidthWS);
676
677 if (thisSampleHistFuncs.size() == 1) {
678 // Just one function. Book it.
679 shapeList.add(*thisSampleHistFuncs.front());
680 } else {
681 // Have multiple functions. We need to multiply them.
682 std::string name = thisSampleHistFuncs.front()->GetName();
683 auto pos = name.find("Hist_alpha");
684 if (pos != std::string::npos) {
685 name = name.substr(0, pos) + "shapes";
686 } else if ( (pos = name.find("nominal")) != std::string::npos) {
687 name = name.substr(0, pos) + "shapes";
688 }
689
690 RooProduct shapeProduct(name.c_str(), thisSampleHistFuncs.front()->GetTitle(), RooArgSet(thisSampleHistFuncs.begin(), thisSampleHistFuncs.end()));
691 proto->import(shapeProduct, RecycleConflictNodes());
692 shapeList.add(*proto->function(name.c_str()));
693 }
694 }
695
696 // Sum all samples
697 RooRealSumPdf tot(totName.c_str(), totName.c_str(), shapeList, coefList, true);
698 tot.specialIntegratorConfig(true)->method1D().setLabel("RooBinIntegrator") ;
699 tot.specialIntegratorConfig(true)->method2D().setLabel("RooBinIntegrator") ;
700 tot.specialIntegratorConfig(true)->methodND().setLabel("RooBinIntegrator") ;
701 tot.forceNumInt();
702
703 // for mixed generation in RooSimultaneous
704 tot.setAttribute("GenerateBinned"); // for use with RooSimultaneous::generate in mixed mode
705
706 // Enable the binned likelihood optimization
708 tot.setAttribute("BinnedLikelihood");
709 }
710
711 proto->import(tot, RecycleConflictNodes());
712 }
713
714 //////////////////////////////////////////////////////////////////////////////
715
717
718 FILE* covFile = fopen ((filename).c_str(),"w");
719 fprintf(covFile," ") ;
720 for (auto const *myargi : static_range_cast<RooRealVar *>(*params)) {
721 if(myargi->isConstant()) continue;
722 fprintf(covFile," & %s", myargi->GetName());
723 }
724 fprintf(covFile,"\\\\ \\hline \n" );
725 for (auto const *myargi : static_range_cast<RooRealVar *>(*params)) {
726 if(myargi->isConstant()) continue;
727 fprintf(covFile,"%s", myargi->GetName());
728 for (auto const *myargj : static_range_cast<RooRealVar *>(*params)) {
729 if(myargj->isConstant()) continue;
730 cout << myargi->GetName() << "," << myargj->GetName();
731 fprintf(covFile, " & %.2f", result->correlation(*myargi, *myargj));
732 }
733 cout << endl;
734 fprintf(covFile, " \\\\\n");
735 }
736 fclose(covFile);
737
738 }
739
740
741 ///////////////////////////////////////////////
743
744 // check inputs (see JIRA-6890 )
745
746 if (channel.GetSamples().empty()) {
747 Error("MakeSingleChannelWorkspace",
748 "The input Channel does not contain any sample - return a nullptr");
749 return 0;
750 }
751
752 const TH1* channel_hist_template = channel.GetSamples().front().GetHisto();
753 if (channel_hist_template == nullptr) {
754 channel.CollectHistograms();
755 channel_hist_template = channel.GetSamples().front().GetHisto();
756 }
757 if (channel_hist_template == nullptr) {
758 std::ostringstream stream;
759 stream << "The sample " << channel.GetSamples().front().GetName()
760 << " in channel " << channel.GetName() << " does not contain a histogram. This is the channel:\n";
761 channel.Print(stream);
762 Error("MakeSingleChannelWorkspace", "%s", stream.str().c_str());
763 return 0;
764 }
765
766 if( ! channel.CheckHistograms() ) {
767 std::cout << "MakeSingleChannelWorkspace: Channel: " << channel.GetName()
768 << " has uninitialized histogram pointers" << std::endl;
769 throw hf_exc();
770 }
771
772
773
774 // Set these by hand inside the function
775 vector<string> systToFix = measurement.GetConstantParams();
776 bool doRatio=false;
777
778 // to time the macro
779 TStopwatch t;
780 t.Start();
781 //ES// string channel_name=summary[0].channel;
782 string channel_name = channel.GetName();
783
784 /// MB: reset observable names for each new channel.
785 fObsNameVec.clear();
786
787 /// MB: label observables x,y,z, depending on histogram dimensionality
788 /// GHL: Give it the first sample's nominal histogram as a template
789 /// since the data histogram may not be present
790 if (fObsNameVec.empty()) { GuessObsNameVec(channel_hist_template); }
791
792 for ( unsigned int idx=0; idx<fObsNameVec.size(); ++idx ) {
793 fObsNameVec[idx] = "obs_" + fObsNameVec[idx] + "_" + channel_name ;
794 }
795
796 if (fObsNameVec.empty()) {
797 fObsName= "obs_" + channel_name; // set name ov observable
798 fObsNameVec.push_back( fObsName );
799 }
800
801 if (fObsNameVec.empty() || fObsNameVec.size() >= 3) {
802 throw hf_exc("HistFactory is limited to 1- to 3-dimensional histograms.");
803 }
804
805 cxcoutP(HistFactory) << "\n-----------------------------------------\n"
806 << "\tStarting to process '"
807 << channel_name << "' channel with " << fObsNameVec.size() << " observables"
808 << "\n-----------------------------------------\n" << endl;
809
810 //
811 // our main workspace that we are using to construct the model
812 //
813 RooWorkspace* proto = new RooWorkspace(channel_name.c_str(), (channel_name+" workspace").c_str());
814 auto proto_config = make_unique<ModelConfig>("ModelConfig", proto);
815 proto_config->SetWorkspace(*proto);
816
817 // preprocess functions
818 vector<string>::iterator funcIter = fPreprocessFunctions.begin();
819 for(;funcIter!= fPreprocessFunctions.end(); ++funcIter){
820 cxcoutI(HistFactory) << "will preprocess this line: " << *funcIter <<endl;
821 proto->factory(*funcIter);
822 proto->Print();
823 }
824
825 RooArgSet likelihoodTerms("likelihoodTerms"), constraintTerms("constraintTerms");
826 vector<string> likelihoodTermNames, constraintTermNames, totSystTermNames;
827 // All histogram functions to be multiplied in each sample
828 std::vector<std::vector<RooAbsArg*>> allSampleHistFuncs;
829 std::vector<RooProduct*> sampleScaleFactors;
830
831 std::vector< pair<string,string> > statNamePairs;
832 std::vector< pair<const TH1*, std::unique_ptr<TH1>> > statHistPairs; // <nominal, error>
833 const std::string statFuncName = "mc_stat_" + channel_name;
834
835 string prefix, range;
836
837 /////////////////////////////
838 // shared parameters
839 // this is ratio of lumi to nominal lumi. We will include relative uncertainty in model
840 std::stringstream lumiStr;
841 // lumi range
842 lumiStr << "Lumi[" << fNomLumi << ",0," << 10.*fNomLumi << "]";
843 proto->factory(lumiStr.str());
844 cxcoutI(HistFactory) << "lumi str = " << lumiStr.str() << endl;
845
846 std::stringstream lumiErrorStr;
847 lumiErrorStr << "nominalLumi["<<fNomLumi << ",0,"<<fNomLumi+10*fLumiError<<"]," << fLumiError ;
848 proto->factory("Gaussian::lumiConstraint(Lumi,"+lumiErrorStr.str()+")");
849 proto->var("nominalLumi")->setConstant();
850 proto->defineSet("globalObservables","nominalLumi");
851 //likelihoodTermNames.push_back("lumiConstraint");
852 constraintTermNames.push_back("lumiConstraint");
853 cxcoutI(HistFactory) << "lumi Error str = " << lumiErrorStr.str() << endl;
854
855 //proto->factory("SigXsecOverSM[1.,0.5,1..8]");
856 ///////////////////////////////////
857 // loop through estimates, add expectation, floating bin predictions,
858 // and terms that constrain floating to expectation via uncertainties
859 // GHL: Loop over samples instead, which doesn't contain the data
860 for (Sample& sample : channel.GetSamples()) {
861 string overallSystName = sample.GetName() + "_" + channel_name + "_epsilon";
862
863 string systSourcePrefix = "alpha_";
864
865 // constraintTermNames and totSystTermNames are vectors that are passed
866 // by reference and filled by this method
867 AddConstraintTerms(proto,measurement, systSourcePrefix, overallSystName,
868 sample.GetOverallSysList(), constraintTermNames , totSystTermNames);
869
870 allSampleHistFuncs.emplace_back();
871 std::vector<RooAbsArg*>& sampleHistFuncs = allSampleHistFuncs.back();
872
873 // GHL: Consider passing the NormFactor list instead of the entire sample
874 auto normFactors = CreateNormFactor(proto, channel_name, overallSystName, sample, doRatio);
875 assert(normFactors);
876
877 // Create the string for the object
878 // that is added to the RooRealSumPdf
879 // for this channel
880// string syst_x_expectedPrefix = "";
881
882 // get histogram
883 //ES// TH1* nominal = it->nominal;
884 const TH1* nominal = sample.GetHisto();
885
886 // MB : HACK no option to have both non-hist variations and hist variations ?
887 // get histogram
888 // GHL: Okay, this is going to be non-trivial.
889 // We will loop over histosys's, which contain both
890 // the low hist and the high hist together.
891
892 // Logic:
893 // - If we have no HistoSys's, do part A
894 // - else, if the histo syst's don't match, return (we ignore this case)
895 // - finally, we take the syst's and apply the linear interpolation w/ constraint
896 string expPrefix = sample.GetName() + "_" + channel_name;
897 // create roorealvar observables
898 RooArgList observables = createObservables(sample.GetHisto(), proto);
899 RooHistFunc* nominalHistFunc = MakeExpectedHistFunc(sample.GetHisto(), proto, expPrefix, observables);
900 assert(nominalHistFunc);
901
902 if(sample.GetHistoSysList().empty()) {
903 // If no HistoSys
904 cxcoutI(HistFactory) << sample.GetName() + "_" + channel_name + " has no variation histograms " << endl;
905
906 sampleHistFuncs.push_back(nominalHistFunc);
907 } else {
908 // If there ARE HistoSys(s)
909 // name of source for variation
910 string constraintPrefix = sample.GetName() + "_" + channel_name + "_Hist_alpha";
911
912 // make list of abstract parameters that interpolate in space of variations
913 RooArgList interpParams = makeInterpolationParameters(sample.GetHistoSysList(), *proto);
914
915 // next, cerate the constraint terms
916 for(std::size_t i = 0; i < interpParams.size(); ++i) {
917 bool isUniform = measurement.GetUniformSyst().count(sample.GetHistoSysList()[i].GetName()) > 0;
918 makeGaussianConstraint(interpParams[i], *proto, isUniform, constraintTermNames);
919 }
920
921 // finally, create the interpolated function
922 sampleHistFuncs.push_back( makeLinInterp(interpParams, nominalHistFunc, proto,
923 sample.GetHistoSysList(), constraintPrefix, observables) );
924 }
925
926 sampleHistFuncs.front()->SetTitle( (nominal && strlen(nominal->GetTitle())>0) ? nominal->GetTitle() : sample.GetName().c_str() );
927
928 ////////////////////////////////////
929 // Add StatErrors to this Channel //
930 ////////////////////////////////////
931
932 if( sample.GetStatError().GetActivate() ) {
933
934 if( fObsNameVec.size() > 3 ) {
935 cxcoutF(HistFactory) << "Cannot include Stat Error for histograms of more than 3 dimensions."
936 << std::endl;
937 throw hf_exc();
938 } else {
939
940 // If we are using StatUncertainties, we multiply this object
941 // by the ParamHistFunc and then pass that to the
942 // RooRealSumPdf by appending it's name to the list
943
944 cxcoutI(HistFactory) << "Sample: " << sample.GetName() << " to be included in Stat Error "
945 << "for channel " << channel_name
946 << std::endl;
947
948 string UncertName = sample.GetName() + "_" + channel_name + "_StatAbsolUncert";
949 std::unique_ptr<TH1> statErrorHist;
950
951 if( sample.GetStatError().GetErrorHist() == nullptr ) {
952 // Make the absolute stat error
953 cxcoutI(HistFactory) << "Making Statistical Uncertainty Hist for "
954 << " Channel: " << channel_name
955 << " Sample: " << sample.GetName()
956 << std::endl;
957 statErrorHist.reset(MakeAbsolUncertaintyHist( UncertName, nominal));
958 } else {
959 // clone the error histograms because in case the sample has not error hist
960 // it is created in MakeAbsolUncertainty
961 // we need later to clean statErrorHist
962 statErrorHist.reset(static_cast<TH1*>(sample.GetStatError().GetErrorHist()->Clone()));
963 // We assume the (relative) error is provided.
964 // We must turn it into an absolute error
965 // using the nominal histogram
966 cxcoutI(HistFactory) << "Using external histogram for Stat Errors for "
967 << "\tChannel: " << channel_name
968 << "\tSample: " << sample.GetName()
969 << "\tError Histogram: " << statErrorHist->GetName() << std::endl;
970 // Multiply the relative stat uncertainty by the
971 // nominal to get the overall stat uncertainty
972 statErrorHist->Multiply( nominal );
973 statErrorHist->SetName( UncertName.c_str() );
974 }
975
976 // Save the nominal and error hists
977 // for the building of constraint terms
978 statHistPairs.emplace_back(nominal, std::move(statErrorHist));
979
980 // To do the 'conservative' version, we would need to do some
981 // intervention here. We would probably need to create a different
982 // ParamHistFunc for each sample in the channel. The would nominally
983 // use the same gamma's, so we haven't increased the number of parameters
984 // However, if a bin in the 'nominal' histogram is 0, we simply need to
985 // change the parameter in that bin in the ParamHistFunc for this sample.
986 // We also need to add a constraint term.
987 // Actually, we'd probably not use the ParamHistFunc...?
988 // we could remove the dependence in this ParamHistFunc on the ith gamma
989 // and then create the poisson term: Pois(tau | n_exp)Pois(data | n_exp)
990
991
992 // Next, try to get the common ParamHistFunc (it may have been
993 // created by another sample in this channel)
994 // or create it if it doesn't yet exist:
995 ParamHistFunc* paramHist = dynamic_cast<ParamHistFunc*>( proto->function(statFuncName.c_str()) );
996 if( paramHist == nullptr ) {
997
998 // Get a RooArgSet of the observables:
999 // Names in the list fObsNameVec:
1000 RooArgList theObservables;
1001 std::vector<std::string>::iterator itr = fObsNameVec.begin();
1002 for (int idx=0; itr!=fObsNameVec.end(); ++itr, ++idx ) {
1003 theObservables.add( *proto->var(*itr) );
1004 }
1005
1006 // Create the list of terms to
1007 // control the bin heights:
1008 std::string ParamSetPrefix = "gamma_stat_" + channel_name;
1009 double gammaMin = 0.0;
1010 double gammaMax = 10.0;
1012 ParamSetPrefix.c_str(),
1013 theObservables,
1014 gammaMin, gammaMax);
1015
1016 ParamHistFunc statUncertFunc(statFuncName.c_str(), statFuncName.c_str(),
1017 theObservables, statFactorParams );
1018
1019 proto->import( statUncertFunc, RecycleConflictNodes() );
1020
1021 paramHist = (ParamHistFunc*) proto->function( statFuncName.c_str() );
1022 }
1023
1024 // apply stat function to sample
1025 sampleHistFuncs.push_back(paramHist);
1026 }
1027 } // END: if DoMcStat
1028
1029
1030 ///////////////////////////////////////////
1031 // Create a ShapeFactor for this channel //
1032 ///////////////////////////////////////////
1033
1034 if( sample.GetShapeFactorList().size() > 0 ) {
1035
1036 if( fObsNameVec.size() > 3 ) {
1037 cxcoutF(HistFactory) << "Cannot include Stat Error for histograms of more than 3 dimensions."
1038 << std::endl;
1039 throw hf_exc();
1040 } else {
1041
1042 cxcoutI(HistFactory) << "Sample: " << sample.GetName() << " in channel: " << channel_name
1043 << " to be include a ShapeFactor."
1044 << std::endl;
1045
1046 for(unsigned int i=0; i < sample.GetShapeFactorList().size(); ++i) {
1047
1048 ShapeFactor& shapeFactor = sample.GetShapeFactorList().at(i);
1049
1050 std::string funcName = channel_name + "_" + shapeFactor.GetName() + "_shapeFactor";
1051 ParamHistFunc* paramHist = (ParamHistFunc*) proto->function( funcName.c_str() );
1052 if( paramHist == nullptr ) {
1053
1054 RooArgList theObservables;
1055 std::vector<std::string>::iterator itr = fObsNameVec.begin();
1056 for (int idx=0; itr!=fObsNameVec.end(); ++itr, ++idx ) {
1057 theObservables.add( *proto->var(*itr) );
1058 }
1059
1060 // Create the Parameters
1061 std::string funcParams = "gamma_" + shapeFactor.GetName();
1062
1063 // GHL: Again, we are putting hard ranges on the gamma's
1064 // We should change this to range from 0 to /inf
1065 RooArgList shapeFactorParams = ParamHistFunc::createParamSet(*proto,
1066 funcParams.c_str(),
1067 theObservables, 0, 1000);
1068
1069 // Create the Function
1070 ParamHistFunc shapeFactorFunc( funcName.c_str(), funcName.c_str(),
1071 theObservables, shapeFactorParams );
1072
1073 // Set an initial shape, if requested
1074 if( shapeFactor.GetInitialShape() != nullptr ) {
1075 TH1* initialShape = static_cast<TH1*>(shapeFactor.GetInitialShape()->Clone());
1076 cxcoutI(HistFactory) << "Setting Shape Factor: " << shapeFactor.GetName()
1077 << " to have initial shape from hist: "
1078 << initialShape->GetName()
1079 << std::endl;
1080 shapeFactorFunc.setShape( initialShape );
1081 }
1082
1083 // Set the variables constant, if requested
1084 if( shapeFactor.IsConstant() ) {
1085 cxcoutI(HistFactory) << "Setting Shape Factor: " << shapeFactor.GetName()
1086 << " to be constant" << std::endl;
1087 shapeFactorFunc.setConstant(true);
1088 }
1089
1090 proto->import( shapeFactorFunc, RecycleConflictNodes() );
1091 paramHist = (ParamHistFunc*) proto->function( funcName.c_str() );
1092
1093 } // End: Create ShapeFactor ParamHistFunc
1094
1095 sampleHistFuncs.push_back(paramHist);
1096 } // End loop over ShapeFactor Systematics
1097 }
1098 } // End: if ShapeFactorName!=""
1099
1100
1101 ////////////////////////////////////////
1102 // Create a ShapeSys for this channel //
1103 ////////////////////////////////////////
1104
1105 if( !sample.GetShapeSysList().empty() ) {
1106
1107 if( fObsNameVec.size() > 3 ) {
1108 cxcoutF(HistFactory) << "Cannot include Stat Error for histograms of more than 3 dimensions."
1109 << std::endl;
1110 throw hf_exc();
1111 } else {
1112
1113 // List of ShapeSys ParamHistFuncs
1114 std::vector<string> ShapeSysNames;
1115
1116 for( unsigned int i = 0; i < sample.GetShapeSysList().size(); ++i) {
1117
1118 // Create the ParamHistFunc's
1119 // Create their constraint terms and add them
1120 // to the list of constraint terms
1121
1122 // Create a single RooProduct over all of these
1123 // paramHistFunc's
1124
1125 // Send the name of that product to the RooRealSumPdf
1126
1127 RooStats::HistFactory::ShapeSys& shapeSys = sample.GetShapeSysList().at(i);
1128
1129 cxcoutI(HistFactory) << "Sample: " << sample.GetName() << " in channel: " << channel_name
1130 << " to include a ShapeSys." << std::endl;
1131
1132 std::string funcName = channel_name + "_" + shapeSys.GetName() + "_ShapeSys";
1133 ShapeSysNames.push_back( funcName );
1134 ParamHistFunc* paramHist = (ParamHistFunc*) proto->function( funcName.c_str() );
1135 if( paramHist == nullptr ) {
1136
1137 //std::string funcParams = "gamma_" + it->shapeFactorName;
1138 //paramHist = CreateParamHistFunc( proto, fObsNameVec, funcParams, funcName );
1139
1140 RooArgList theObservables;
1141 std::vector<std::string>::iterator itr = fObsNameVec.begin();
1142 for(; itr!=fObsNameVec.end(); ++itr ) {
1143 theObservables.add( *proto->var(*itr) );
1144 }
1145
1146 // Create the Parameters
1147 std::string funcParams = "gamma_" + shapeSys.GetName();
1148 RooArgList shapeFactorParams = ParamHistFunc::createParamSet(*proto,
1149 funcParams.c_str(),
1150 theObservables, 0, 10);
1151
1152 // Create the Function
1153 ParamHistFunc shapeFactorFunc( funcName.c_str(), funcName.c_str(),
1154 theObservables, shapeFactorParams );
1155
1156 proto->import( shapeFactorFunc, RecycleConflictNodes() );
1157 paramHist = (ParamHistFunc*) proto->function( funcName.c_str() );
1158
1159 } // End: Create ShapeFactor ParamHistFunc
1160
1161 // Create the constraint terms and add
1162 // them to the workspace (proto)
1163 // as well as the list of constraint terms (constraintTermNames)
1164
1165 // The syst should be a fractional error
1166 const TH1* shapeErrorHist = shapeSys.GetErrorHist();
1167
1168 // Constraint::Type shapeConstraintType = Constraint::Gaussian;
1169 Constraint::Type systype = shapeSys.GetConstraintType();
1170 if( systype == Constraint::Gaussian) {
1171 systype = Constraint::Gaussian;
1172 }
1173 if( systype == Constraint::Poisson ) {
1174 systype = Constraint::Poisson;
1175 }
1176
1177 double minShapeUncertainty = 0.0;
1178 RooArgList shapeConstraints = createStatConstraintTerms(proto, constraintTermNames,
1179 *paramHist, shapeErrorHist,
1180 systype,
1181 minShapeUncertainty);
1182
1183 } // End: Loop over ShapeSys vector in this EstimateSummary
1184
1185 // Now that we have the list of ShapeSys ParamHistFunc names,
1186 // we create the total RooProduct
1187 // we multiply the expected functio
1188
1189
1190 for( unsigned int i = 0; i < ShapeSysNames.size(); ++i ) {
1191 auto func = proto->function(ShapeSysNames.at(i).c_str());
1192 assert(func);
1193 sampleHistFuncs.push_back(func);
1194 }
1195
1196 } // End: NumObsVar == 1
1197
1198 } // End: !GetShapeSysList.empty()
1199
1200
1201 // GHL: This was pretty confusing before,
1202 // hopefully using the measurement directly
1203 // will improve it
1204 auto lumi = proto->arg("Lumi");
1205 if( !sample.GetNormalizeByTheory() ) {
1206 if (!lumi) {
1207 TString lumiParamString;
1208 lumiParamString += measurement.GetLumi();
1209 lumiParamString.ReplaceAll(' ', TString());
1210 lumi = proto->factory(("Lumi[" + lumiParamString + "]").Data());
1211 } else {
1212 static_cast<RooAbsRealLValue*>(lumi)->setVal(measurement.GetLumi());
1213 }
1214 }
1215 assert(lumi);
1216 normFactors->addTerm(lumi);
1217
1218 // Append the name of the "node"
1219 // that is to be summed with the
1220 // RooRealSumPdf
1221 proto->import(*normFactors, RecycleConflictNodes());
1222 auto normFactorsInWS = dynamic_cast<RooProduct*>(proto->arg(normFactors->GetName()));
1223 assert(normFactorsInWS);
1224
1225 sampleScaleFactors.push_back(normFactorsInWS);
1226 } // END: Loop over EstimateSummaries
1227
1228 // If a non-zero number of samples call for
1229 // Stat Uncertainties, create the statFactor functions
1230 if(!statHistPairs.empty()) {
1231
1232 // Create the histogram of (binwise)
1233 // stat uncertainties:
1234 unique_ptr<TH1> fracStatError( MakeScaledUncertaintyHist( channel_name + "_StatUncert" + "_RelErr", statHistPairs) );
1235 if( fracStatError == nullptr ) {
1236 cxcoutE(HistFactory) << "Error: Failed to make ScaledUncertaintyHist for: "
1237 << channel_name + "_StatUncert" + "_RelErr" << std::endl;
1238 throw hf_exc();
1239 }
1240
1241 // Using this TH1* of fractinal stat errors,
1242 // create a set of constraint terms:
1243 ParamHistFunc* chanStatUncertFunc = (ParamHistFunc*) proto->function( statFuncName.c_str() );
1244 cxcoutI(HistFactory) << "About to create Constraint Terms from: "
1245 << chanStatUncertFunc->GetName()
1246 << " params: " << chanStatUncertFunc->paramList()
1247 << std::endl;
1248
1249 // Get the constraint type and the
1250 // rel error threshold from the (last)
1251 // EstimateSummary looped over (but all
1252 // should be the same)
1253
1254 // Get the type of StatError constraint from the channel
1255 Constraint::Type statConstraintType = channel.GetStatErrorConfig().GetConstraintType();
1256 if( statConstraintType == Constraint::Gaussian) {
1257 cxcoutI(HistFactory) << "Using Gaussian StatErrors in channel: " << channel.GetName() << std::endl;
1258 }
1259 if( statConstraintType == Constraint::Poisson ) {
1260 cxcoutI(HistFactory) << "Using Poisson StatErrors in channel: " << channel.GetName() << std::endl;
1261 }
1262
1263 double statRelErrorThreshold = channel.GetStatErrorConfig().GetRelErrorThreshold();
1264 RooArgList statConstraints = createStatConstraintTerms(proto, constraintTermNames,
1265 *chanStatUncertFunc, fracStatError.get(),
1266 statConstraintType,
1267 statRelErrorThreshold);
1268
1269 } // END: Loop over stat Hist Pairs
1270
1271
1272 ///////////////////////////////////
1273 // for ith bin calculate totN_i = lumi * sum_j expected_j * syst_j
1274 MakeTotalExpected(proto, channel_name+"_model",
1275 sampleScaleFactors, allSampleHistFuncs);
1276 likelihoodTermNames.push_back(channel_name+"_model");
1277
1278 //////////////////////////////////////
1279 // fix specified parameters
1280 for(unsigned int i=0; i<systToFix.size(); ++i){
1281 RooRealVar* temp = proto->var(systToFix.at(i));
1282 if(temp) {
1283 // set the parameter constant
1284 temp->setConstant();
1285
1286 // remove the corresponding auxiliary observable from the global observables
1287 RooRealVar* auxMeas = nullptr;
1288 if(systToFix.at(i)=="Lumi"){
1289 auxMeas = proto->var("nominalLumi");
1290 } else {
1291 auxMeas = proto->var(std::string("nom_") + temp->GetName());
1292 }
1293
1294 if(auxMeas){
1295 const_cast<RooArgSet*>(proto->set("globalObservables"))->remove(*auxMeas);
1296 } else{
1297 cxcoutE(HistFactory) << "could not corresponding auxiliary measurement "
1298 << TString::Format("nom_%s",temp->GetName()) << endl;
1299 }
1300 } else {
1301 cxcoutE(HistFactory) << "could not find variable " << systToFix.at(i)
1302 << " could not set it to constant" << endl;
1303 }
1304 }
1305
1306 //////////////////////////////////////
1307 // final proto model
1308 for(unsigned int i=0; i<constraintTermNames.size(); ++i){
1309 RooAbsArg* proto_arg = (proto->arg(constraintTermNames[i].c_str()));
1310 if( proto_arg==nullptr ) {
1311 cxcoutF(HistFactory) << "Error: Cannot find arg set: " << constraintTermNames.at(i)
1312 << " in workspace: " << proto->GetName() << std::endl;
1313 throw hf_exc();
1314 }
1315 constraintTerms.add( *proto_arg );
1316 // constraintTerms.add(* proto_arg(proto->arg(constraintTermNames[i].c_str())) );
1317 }
1318 for(unsigned int i=0; i<likelihoodTermNames.size(); ++i){
1319 RooAbsArg* proto_arg = (proto->arg(likelihoodTermNames[i].c_str()));
1320 if( proto_arg==nullptr ) {
1321 cxcoutF(HistFactory) << "Error: Cannot find arg set: " << likelihoodTermNames.at(i)
1322 << " in workspace: " << proto->GetName() << std::endl;
1323 throw hf_exc();
1324 }
1325 likelihoodTerms.add( *proto_arg );
1326 }
1327 proto->defineSet("constraintTerms",constraintTerms);
1328 proto->defineSet("likelihoodTerms",likelihoodTerms);
1329
1330 // list of observables
1331 RooArgList observables;
1332 std::string observablesStr;
1333
1334 std::vector<std::string>::iterator itr = fObsNameVec.begin();
1335 for(; itr!=fObsNameVec.end(); ++itr ) {
1336 observables.add( *proto->var(*itr) );
1337 if (!observablesStr.empty()) { observablesStr += ","; }
1338 observablesStr += *itr;
1339 }
1340
1341 // We create two sets, one for backwards compatability
1342 // The other to make a consistent naming convention
1343 // between individual channels and the combined workspace
1344 proto->defineSet("observables", TString::Format("%s",observablesStr.c_str()));
1345 proto->defineSet("observablesSet", TString::Format("%s",observablesStr.c_str()));
1346
1347 // Create the ParamHistFunc
1348 // after observables have been made
1349 cxcoutP(HistFactory) << "\n-----------------------------------------\n"
1350 << "\timport model into workspace"
1351 << "\n-----------------------------------------\n" << endl;
1352
1353 auto model = make_unique<RooProdPdf>(
1354 ("model_"+channel_name).c_str(), // MB : have changed this into conditional pdf. Much faster for toys!
1355 "product of Poissons accross bins for a single channel",
1356 constraintTerms, Conditional(likelihoodTerms,observables));
1357 // can give channel a title by setting title of corresponding data histogram
1358 if (channel.GetData().GetHisto() && strlen(channel.GetData().GetHisto()->GetTitle())>0) {
1359 model->SetTitle(channel.GetData().GetHisto()->GetTitle());
1360 }
1361 proto->import(*model,RecycleConflictNodes());
1362
1363 proto_config->SetPdf(*model);
1364 proto_config->SetObservables(observables);
1365 proto_config->SetGlobalObservables(*proto->set("globalObservables"));
1366 // proto->writeToFile(("results/model_"+channel+".root").c_str());
1367 // fill out nuisance parameters in model config
1368 // proto_config->GuessObsAndNuisance(*proto->data("asimovData"));
1369 proto->import(*proto_config,proto_config->GetName());
1370 proto->importClassCode();
1371
1372 ///////////////////////////
1373 // make data sets
1374 // THis works and is natural, but the memory size of the simultaneous dataset grows exponentially with channels
1375 // New Asimov Generation: Use the code in the Asymptotic calculator
1376 // Need to get the ModelConfig...
1377 int asymcalcPrintLevel = 0;
1378 if (RooMsgService::instance().isActive(static_cast<TObject*>(nullptr), RooFit::HistFactory, RooFit::INFO)) asymcalcPrintLevel = 1;
1379 if (RooMsgService::instance().isActive(static_cast<TObject*>(nullptr), RooFit::HistFactory, RooFit::DEBUG)) asymcalcPrintLevel = 2;
1380 AsymptoticCalculator::SetPrintLevel(asymcalcPrintLevel);
1381 unique_ptr<RooAbsData> asimov_dataset(AsymptoticCalculator::GenerateAsimovData(*model, observables));
1382 proto->import(dynamic_cast<RooDataSet&>(*asimov_dataset), Rename("asimovData"));
1383
1384 // GHL: Determine to use data if the hist isn't 'nullptr'
1385 if(TH1 const* mnominal = channel.GetData().GetHisto()) {
1386 // This works and is natural, but the memory size of the simultaneous
1387 // dataset grows exponentially with channels.
1388 RooDataSet dataset{"obsData","",*proto->set("observables"), RooFit::WeightVar("weightVar")};
1389 ConfigureHistFactoryDataset( dataset, *mnominal, *proto, fObsNameVec );
1390 proto->import(dataset);
1391 } // End: Has non-null 'data' entry
1392
1393
1394 for(auto const& data : channel.GetAdditionalData()) {
1395 if(data.GetName().empty()) {
1396 cxcoutF(HistFactory) << "Error: Additional Data histogram for channel: " << channel.GetName()
1397 << " has no name! The name always needs to be set for additional datasets, "
1398 << "either via the \"Name\" tag in the XML or via RooStats::HistFactory::Data::SetName()." << std::endl;
1399 throw hf_exc();
1400 }
1401 std::string const& dataName = data.GetName();
1402 TH1 const* mnominal = data.GetHisto();
1403 if( !mnominal ) {
1404 cxcoutF(HistFactory) << "Error: Additional Data histogram for channel: " << channel.GetName()
1405 << " with name: " << dataName << " is nullptr" << std::endl;
1406 throw hf_exc();
1407 }
1408
1409 // THis works and is natural, but the memory size of the simultaneous dataset grows exponentially with channels
1410 RooDataSet dataset{dataName.c_str(), "", *proto->set("observables"), RooFit::WeightVar("weightVar")};
1411 ConfigureHistFactoryDataset( dataset, *mnominal, *proto, fObsNameVec );
1412 proto->import(dataset);
1413
1414 }
1415
1416 if (RooMsgService::instance().isActive(static_cast<TObject*>(nullptr), RooFit::HistFactory, RooFit::INFO))
1417 proto->Print();
1418
1419 return proto;
1420 }
1421
1422
1424 TH1 const& mnominal,
1426 std::vector<std::string> const& obsNameVec) {
1427
1428 // Take a RooDataSet and fill it with the entries
1429 // from a TH1*, using the observable names to
1430 // determine the columns
1431
1432 if (obsNameVec.empty() ) {
1433 Error("ConfigureHistFactoryDataset","Invalid input - return");
1434 return;
1435 }
1436
1437 TAxis const* ax = mnominal.GetXaxis();
1438 TAxis const* ay = mnominal.GetYaxis();
1439 TAxis const* az = mnominal.GetZaxis();
1440
1441 for (int i=1; i<=ax->GetNbins(); ++i) { // 1 or more dimension
1442
1443 double xval = ax->GetBinCenter(i);
1444 proto.var( obsNameVec[0] )->setVal( xval );
1445
1446 if(obsNameVec.size()==1) {
1447 double fval = mnominal.GetBinContent(i);
1448 obsDataUnbinned.add( *proto.set("observables"), fval );
1449 } else { // 2 or more dimensions
1450
1451 for(int j=1; j<=ay->GetNbins(); ++j) {
1452 double yval = ay->GetBinCenter(j);
1453 proto.var( obsNameVec[1] )->setVal( yval );
1454
1455 if(obsNameVec.size()==2) {
1456 double fval = mnominal.GetBinContent(i,j);
1457 obsDataUnbinned.add( *proto.set("observables"), fval );
1458 } else { // 3 dimensions
1459
1460 for(int k=1; k<=az->GetNbins(); ++k) {
1461 double zval = az->GetBinCenter(k);
1462 proto.var( obsNameVec[2] )->setVal( zval );
1463 double fval = mnominal.GetBinContent(i,j,k);
1464 obsDataUnbinned.add( *proto.set("observables"), fval );
1465 }
1466 }
1467 }
1468 }
1469 }
1470 }
1471
1473 {
1474 fObsNameVec.clear();
1475
1476 // determine histogram dimensionality
1477 unsigned int histndim(1);
1478 std::string classname = hist->ClassName();
1479 if (classname.find("TH1")==0) { histndim=1; }
1480 else if (classname.find("TH2")==0) { histndim=2; }
1481 else if (classname.find("TH3")==0) { histndim=3; }
1482
1483 for ( unsigned int idx=0; idx<histndim; ++idx ) {
1484 if (idx==0) { fObsNameVec.push_back("x"); }
1485 if (idx==1) { fObsNameVec.push_back("y"); }
1486 if (idx==2) { fObsNameVec.push_back("z"); }
1487 }
1488 }
1489
1490
1491 RooWorkspace* HistoToWorkspaceFactoryFast::MakeCombinedModel(vector<string> ch_names, vector<std::unique_ptr<RooWorkspace>>& chs)
1492 {
1494
1495 // check first the inputs (see JIRA-6890)
1496 if (ch_names.empty() || chs.empty() ) {
1497 Error("MakeCombinedModel","Input vectors are empty - return a nullptr");
1498 return 0;
1499 }
1500 if (chs.size() < ch_names.size() ) {
1501 Error("MakeCombinedModel","Input vector of workspace has an invalid size - return a nullptr");
1502 return 0;
1503 }
1504
1505 //
1506 /// These things were used for debugging. Maybe useful in the future
1507 //
1508
1509 map<string, RooAbsPdf*> pdfMap;
1510 vector<RooAbsPdf*> models;
1511
1512 RooArgList obsList;
1513 for(unsigned int i = 0; i< ch_names.size(); ++i){
1514 ModelConfig * config = (ModelConfig *) chs[i]->obj("ModelConfig");
1515 obsList.add(*config->GetObservables());
1516 }
1517 cxcoutI(HistFactory) <<"full list of observables:\n" << obsList << std::endl;
1518
1519 RooArgSet globalObs;
1520 stringstream channelString;
1521 channelString << "channelCat[";
1522 for(unsigned int i = 0; i< ch_names.size(); ++i){
1523 string channel_name=ch_names[i];
1524 if (i == 0 && isdigit(channel_name[0])) {
1525 throw std::invalid_argument("The first channel name for HistFactory cannot start with a digit. Got " + channel_name);
1526 }
1527 if (channel_name.find(',') != std::string::npos) {
1528 throw std::invalid_argument("Channel names for HistFactory cannot contain ','. Got " + channel_name);
1529 }
1530
1531 if (i == 0) channelString << channel_name ;
1532 else channelString << ',' << channel_name ;
1533 RooWorkspace * ch=chs[i].get();
1534
1535 RooAbsPdf* model = ch->pdf("model_"+channel_name);
1536 if(!model) cout <<"failed to find model for channel"<<endl;
1537 // cout << "int = " << model->createIntegral(*obsN)->getVal() << endl;;
1538 models.push_back(model);
1539 globalObs.add(*ch->set("globalObservables"), /*silent=*/true); // silent because observables might exist in other channel.
1540
1541 // constrainedParams->add( * ch->set("constrainedParams") );
1542 pdfMap[channel_name]=model;
1543 }
1544 channelString << "]";
1545
1546 cxcoutP(HistFactory) << "\n-----------------------------------------\n"
1547 << "\tEntering combination"
1548 << "\n-----------------------------------------\n" << endl;
1549 RooWorkspace* combined = new RooWorkspace("combined");
1550 // RooWorkspace* combined = chs[0];
1551
1552
1553 RooCategory* channelCat = dynamic_cast<RooCategory*>( combined->factory(channelString.str()) );
1554 if (!channelCat) throw std::runtime_error("Unable to construct a category from string " + channelString.str());
1555
1556 auto simPdf= std::make_unique<RooSimultaneous>("simPdf","",pdfMap, *channelCat);
1557 auto combined_config = std::make_unique<ModelConfig>("ModelConfig", combined);
1558 combined_config->SetWorkspace(*combined);
1559 // combined_config->SetNuisanceParameters(*constrainedParams);
1560
1561 combined->import(globalObs);
1562 combined->defineSet("globalObservables",globalObs);
1563 combined_config->SetGlobalObservables(*combined->set("globalObservables"));
1564
1565 combined->factory("weightVar[0,-1e10,1e10]");
1566 obsList.add(*combined->var("weightVar"));
1567 combined->defineSet("observables",{obsList, *channelCat}, /*importMissing=*/true);
1568 combined_config->SetObservables(*combined->set("observables"));
1569
1570
1571 // Now merge the observable datasets across the channels
1572 for(RooAbsData * data : chs[0]->allData()) {
1573 // We are excluding the Asimov data, because it needs to be regenerated
1574 // later after the parameter values are set.
1575 if(std::string("asimovData") != data->GetName()) {
1576 MergeDataSets(combined, chs, ch_names, data->GetName(), obsList, channelCat);
1577 }
1578 }
1579
1580
1581 if (RooMsgService::instance().isActive(static_cast<TObject*>(nullptr), RooFit::HistFactory, RooFit::INFO))
1582 combined->Print();
1583
1584 cxcoutP(HistFactory) << "\n-----------------------------------------\n"
1585 << "\tImporting combined model"
1586 << "\n-----------------------------------------\n" << endl;
1587 combined->import(*simPdf,RecycleConflictNodes());
1588
1589 std::map< std::string, double>::iterator param_itr = fParamValues.begin();
1590 for( ; param_itr != fParamValues.end(); ++param_itr ){
1591 // make sure they are fixed
1592 std::string paramName = param_itr->first;
1593 double paramVal = param_itr->second;
1594
1595 if(RooRealVar* temp = combined->var( paramName )) {
1596 temp->setVal( paramVal );
1597 cxcoutI(HistFactory) <<"setting " << paramName << " to the value: " << paramVal << endl;
1598 } else
1599 cxcoutE(HistFactory) << "could not find variable " << paramName << " could not set its value" << endl;
1600 }
1601
1602
1603 for(unsigned int i=0; i<fSystToFix.size(); ++i){
1604 // make sure they are fixed
1605 if(RooRealVar* temp = combined->var(fSystToFix[i])) {
1606 temp->setConstant();
1607 cxcoutI(HistFactory) <<"setting " << fSystToFix.at(i) << " constant" << endl;
1608 } else
1609 cxcoutE(HistFactory) << "could not find variable " << fSystToFix.at(i) << " could not set it to constant" << endl;
1610 }
1611
1612 ///
1613 /// writing out the model in graphViz
1614 ///
1615 // RooAbsPdf* customized=combined->pdf("simPdf");
1616 //combined_config->SetPdf(*customized);
1617 combined_config->SetPdf(*simPdf);
1618 // combined_config->GuessObsAndNuisance(*simData);
1619 // customized->graphVizTree(("results/"+fResultsPrefixStr.str()+"_simul.dot").c_str());
1620 combined->import(*combined_config,combined_config->GetName());
1621 combined->importClassCode();
1622 // combined->writeToFile("results/model_combined.root");
1623
1624
1625 ////////////////////////////////////////////
1626 // Make toy simultaneous dataset
1627 cxcoutP(HistFactory) << "\n-----------------------------------------\n"
1628 << "\tcreate toy data for " << channelString.str()
1629 << "\n-----------------------------------------\n" << endl;
1630
1631
1632 // now with weighted datasets
1633 // First Asimov
1634
1635 // Create Asimov data for the combined dataset
1636 std::unique_ptr<RooDataSet> asimov_combined{static_cast<RooDataSet*>(AsymptoticCalculator::GenerateAsimovData(
1637 *combined->pdf("simPdf"),
1638 obsList))};
1639 if( asimov_combined ) {
1640 combined->import( *asimov_combined, Rename("asimovData"));
1641 }
1642 else {
1643 std::cout << "Error: Failed to create combined asimov dataset" << std::endl;
1644 throw hf_exc();
1645 }
1646
1647 return combined;
1648 }
1649
1650
1652 std::vector<std::unique_ptr<RooWorkspace>>& wspace_vec,
1653 std::vector<std::string> const& channel_names,
1654 std::string const& dataSetName,
1655 RooArgList const& obsList,
1656 RooCategory* channelCat) {
1657
1658 // Create the total dataset
1659 std::unique_ptr<RooDataSet> simData;
1660
1661 // Loop through channels, get their individual datasets,
1662 // and add them to the combined dataset
1663 for(unsigned int i = 0; i< channel_names.size(); ++i){
1664
1665 // Grab the dataset for the existing channel
1666 cxcoutPHF << "Merging data for channel " << channel_names[i].c_str() << std::endl;
1667 RooDataSet* obsDataInChannel = (RooDataSet*) wspace_vec[i]->data(dataSetName.c_str());
1668 if( !obsDataInChannel ) {
1669 std::cout << "Error: Can't find DataSet: " << dataSetName
1670 << " in channel: " << channel_names.at(i)
1671 << std::endl;
1672 throw hf_exc();
1673 }
1674
1675 // Create the new Dataset
1676 auto tempData = std::make_unique<RooDataSet>(channel_names[i].c_str(),"",
1677 obsList, Index(*channelCat),
1678 WeightVar("weightVar"),
1679 Import(channel_names[i].c_str(),*obsDataInChannel));
1680 if(simData) {
1681 simData->append(*tempData);
1682 }
1683 else {
1684 simData = std::move(tempData);
1685 }
1686 } // End Loop Over Channels
1687
1688 // Check that we successfully created the dataset
1689 // and import it into the workspace
1690 if(simData) {
1691 combined->import(*simData, Rename(dataSetName.c_str()));
1692 return static_cast<RooDataSet*>(combined->data(dataSetName));
1693 }
1694 else {
1695 std::cout << "Error: Unable to merge observable datasets" << std::endl;
1696 throw hf_exc();
1697 return nullptr;
1698 }
1699 }
1700
1701
1703
1704 // Take a nominal TH1* and create
1705 // a TH1 representing the binwise
1706 // errors (taken from the nominal TH1)
1707
1708 TH1* ErrorHist = (TH1*) Nominal->Clone( Name.c_str() );
1709 ErrorHist->Reset();
1710
1711 Int_t numBins = Nominal->GetNbinsX()*Nominal->GetNbinsY()*Nominal->GetNbinsZ();
1712 Int_t binNumber = 0;
1713
1714 // Loop over bins
1715 for( Int_t i_bin = 0; i_bin < numBins; ++i_bin) {
1716
1717 binNumber++;
1718 // Ignore underflow / overflow
1719 while( Nominal->IsBinUnderflow(binNumber) || Nominal->IsBinOverflow(binNumber) ){
1720 binNumber++;
1721 }
1722
1723 double histError = Nominal->GetBinError( binNumber );
1724
1725 // Check that histError != NAN
1726 if( histError != histError ) {
1727 std::cout << "Warning: In histogram " << Nominal->GetName()
1728 << " bin error for bin " << i_bin
1729 << " is NAN. Not using Error!!!"
1730 << std::endl;
1731 throw hf_exc();
1732 //histError = sqrt( histContent );
1733 //histError = 0;
1734 }
1735
1736 // Check that histError ! < 0
1737 if( histError < 0 ) {
1738 std::cout << "Warning: In histogram " << Nominal->GetName()
1739 << " bin error for bin " << binNumber
1740 << " is < 0. Setting Error to 0"
1741 << std::endl;
1742 //histError = sqrt( histContent );
1743 histError = 0;
1744 }
1745
1746 ErrorHist->SetBinContent( binNumber, histError );
1747
1748 }
1749
1750 return ErrorHist;
1751
1752 }
1753
1754 // Take a list of < nominal, absolError > TH1* pairs
1755 // and construct a single histogram representing the
1756 // total fractional error as:
1757
1758 // UncertInQuad(bin i) = Sum: absolUncert*absolUncert
1759 // Total(bin i) = Sum: Value
1760 //
1761 // TotalFracError(bin i) = Sqrt( UncertInQuad(i) ) / TotalBin(i)
1762 std::unique_ptr<TH1> HistoToWorkspaceFactoryFast::MakeScaledUncertaintyHist( const std::string& Name, std::vector< std::pair<const TH1*, std::unique_ptr<TH1>> > const& HistVec ) const {
1763
1764
1765 unsigned int numHists = HistVec.size();
1766
1767 if( numHists == 0 ) {
1768 cxcoutE(HistFactory) << "Warning: Empty Hist Vector, cannot create total uncertainty" << std::endl;
1769 return nullptr;
1770 }
1771
1772 const TH1* HistTemplate = HistVec.at(0).first;
1773 Int_t numBins = HistTemplate->GetNbinsX()*HistTemplate->GetNbinsY()*HistTemplate->GetNbinsZ();
1774
1775 // Check that all histograms
1776 // have the same bins
1777 for( unsigned int i = 0; i < HistVec.size(); ++i ) {
1778
1779 const TH1* nominal = HistVec.at(i).first;
1780 const TH1* error = HistVec.at(i).second.get();
1781
1782 if( nominal->GetNbinsX()*nominal->GetNbinsY()*nominal->GetNbinsZ() != numBins ) {
1783 cxcoutE(HistFactory) << "Error: Provided hists have unequal bins" << std::endl;
1784 return nullptr;
1785 }
1786 if( error->GetNbinsX()*error->GetNbinsY()*error->GetNbinsZ() != numBins ) {
1787 cxcoutE(HistFactory) << "Error: Provided hists have unequal bins" << std::endl;
1788 return nullptr;
1789 }
1790 }
1791
1792 std::vector<double> TotalBinContent( numBins, 0.0);
1793 std::vector<double> HistErrorsSqr( numBins, 0.0);
1794
1795 Int_t binNumber = 0;
1796
1797 // Loop over bins
1798 for( Int_t i_bins = 0; i_bins < numBins; ++i_bins) {
1799
1800 binNumber++;
1801 while( HistTemplate->IsBinUnderflow(binNumber) || HistTemplate->IsBinOverflow(binNumber) ){
1802 binNumber++;
1803 }
1804
1805 for( unsigned int i_hist = 0; i_hist < numHists; ++i_hist ) {
1806
1807 const TH1* nominal = HistVec.at(i_hist).first;
1808 const TH1* error = HistVec.at(i_hist).second.get();
1809
1810 //Int_t binNumber = i_bins + 1;
1811
1812 double histValue = nominal->GetBinContent( binNumber );
1813 double histError = error->GetBinContent( binNumber );
1814 /*
1815 std::cout << " Getting Bin content for Stat Uncertainty"
1816 << " Nom name: " << nominal->GetName()
1817 << " Err name: " << error->GetName()
1818 << " HistNumber: " << i_hist << " bin: " << binNumber
1819 << " Value: " << histValue << " Error: " << histError
1820 << std::endl;
1821 */
1822
1823 if( histError != histError ) {
1824 cxcoutE(HistFactory) << "In histogram " << error->GetName()
1825 << " bin error for bin " << binNumber
1826 << " is NAN. Not using error!!";
1827 throw hf_exc();
1828 }
1829
1830 TotalBinContent.at(i_bins) += histValue;
1831 HistErrorsSqr.at(i_bins) += histError*histError; // Add in quadrature
1832
1833 }
1834 }
1835
1836 binNumber = 0;
1837
1838 // Creat the output histogram
1839 TH1* ErrorHist = (TH1*) HistTemplate->Clone( Name.c_str() );
1840 ErrorHist->Reset();
1841
1842 // Fill the output histogram
1843 for( Int_t i = 0; i < numBins; ++i) {
1844
1845 // Int_t binNumber = i + 1;
1846 binNumber++;
1847 while( ErrorHist->IsBinUnderflow(binNumber) || ErrorHist->IsBinOverflow(binNumber) ){
1848 binNumber++;
1849 }
1850
1851 double ErrorsSqr = HistErrorsSqr.at(i);
1852 double TotalVal = TotalBinContent.at(i);
1853
1854 if( TotalVal <= 0 ) {
1855 cxcoutW(HistFactory) << "Warning: Sum of histograms for bin: " << binNumber
1856 << " is <= 0. Setting error to 0"
1857 << std::endl;
1858
1859 ErrorHist->SetBinContent( binNumber, 0.0 );
1860 continue;
1861 }
1862
1863 double RelativeError = sqrt(ErrorsSqr) / TotalVal;
1864
1865 // If we otherwise get a NAN
1866 // it's an error
1867 if( RelativeError != RelativeError ) {
1868 cxcoutE(HistFactory) << "Error: bin " << i << " error is NAN\n"
1869 << " HistErrorsSqr: " << ErrorsSqr
1870 << " TotalVal: " << TotalVal;
1871 throw hf_exc();
1872 }
1873
1874 // 0th entry in vector is
1875 // the 1st bin in TH1
1876 // (we ignore underflow)
1877
1878 // Error and bin content are interchanged because for some reason, the other functions
1879 // use the bin content to convey the error ...
1880 ErrorHist->SetBinError(binNumber, TotalVal);
1881 ErrorHist->SetBinContent(binNumber, RelativeError);
1882
1883 cxcoutI(HistFactory) << "Making Total Uncertainty for bin " << binNumber
1884 << " Error = " << sqrt(ErrorsSqr)
1885 << " CentralVal = " << TotalVal
1886 << " RelativeError = " << RelativeError << "\n";
1887
1888 }
1889
1890 return std::unique_ptr<TH1>(ErrorHist);
1891}
1892
1893
1894
1896 createStatConstraintTerms( RooWorkspace* proto, vector<string>& constraintTermNames,
1897 ParamHistFunc& paramHist, const TH1* uncertHist,
1898 Constraint::Type type, double minSigma ) {
1899
1900
1901 // Take a RooArgList of RooAbsReal's and
1902 // create N constraint terms (one for
1903 // each gamma) whose relative uncertainty
1904 // is the value of the ith RooAbsReal
1905 //
1906 // The integer "type" controls the type
1907 // of constraint term:
1908 //
1909 // type == 0 : NONE
1910 // type == 1 : Gaussian
1911 // type == 2 : Poisson
1912 // type == 3 : LogNormal
1913
1914 RooArgList ConstraintTerms;
1915
1916 RooArgList paramSet = paramHist.paramList();
1917
1918 // Must get the full size of the TH1
1919 // (No direct method to do this...)
1920 Int_t numBins = uncertHist->GetNbinsX()*uncertHist->GetNbinsY()*uncertHist->GetNbinsZ();
1921 Int_t numParams = paramSet.getSize();
1922 // Int_t numBins = uncertHist->GetNbinsX()*uncertHist->GetNbinsY()*uncertHist->GetNbinsZ();
1923
1924 // Check that there are N elements
1925 // in the RooArgList
1926 if( numBins != numParams ) {
1927 std::cout << "Error: In createStatConstraintTerms, encountered bad number of bins" << std::endl;
1928 std::cout << "Given histogram with " << numBins << " bins,"
1929 << " but require exactly " << numParams << std::endl;
1930 throw hf_exc();
1931 }
1932
1933 Int_t TH1BinNumber = 0;
1934 for( Int_t i = 0; i < paramSet.getSize(); ++i) {
1935
1936 TH1BinNumber++;
1937
1938 while( uncertHist->IsBinUnderflow(TH1BinNumber) || uncertHist->IsBinOverflow(TH1BinNumber) ){
1939 TH1BinNumber++;
1940 }
1941
1942 RooRealVar& gamma = (RooRealVar&) (paramSet[i]);
1943
1944 cxcoutI(HistFactory) << "Creating constraint for: " << gamma.GetName()
1945 << ". Type of constraint: " << type << std::endl;
1946
1947 // Get the sigma from the hist
1948 // (the relative uncertainty)
1949 const double sigmaRel = uncertHist->GetBinContent(TH1BinNumber);
1950
1951 // If the sigma is <= 0,
1952 // do cont create the term
1953 if( sigmaRel <= 0 ){
1954 cxcoutI(HistFactory) << "Not creating constraint term for "
1955 << gamma.GetName()
1956 << " because sigma = " << sigmaRel
1957 << " (sigma<=0)"
1958 << " (TH1 bin number = " << TH1BinNumber << ")"
1959 << std::endl;
1960 gamma.setConstant(true);
1961 continue;
1962 }
1963
1964 // set reasonable ranges for gamma parameters
1965 gamma.setMax( 1 + 5*sigmaRel );
1966 gamma.setMin( 0. );
1967
1968 // Make Constraint Term
1969 std::string constrName = string(gamma.GetName()) + "_constraint";
1970 std::string nomName = string("nom_") + gamma.GetName();
1971 std::string sigmaName = string(gamma.GetName()) + "_sigma";
1972 std::string poisMeanName = string(gamma.GetName()) + "_poisMean";
1973
1974 if( type == Constraint::Gaussian ) {
1975
1976 // Type 1 : RooGaussian
1977
1978 // Make sigma
1979
1980 RooConstVar constrSigma( sigmaName.c_str(), sigmaName.c_str(), sigmaRel );
1981
1982 // Make "observed" value
1983 RooRealVar constrNom(nomName.c_str(), nomName.c_str(), 1.0,0,10);
1984 constrNom.setConstant( true );
1985
1986 // Make the constraint:
1987 RooGaussian gauss( constrName.c_str(), constrName.c_str(),
1988 constrNom, gamma, constrSigma );
1989
1990 proto->import( gauss, RecycleConflictNodes() );
1991
1992 // Give reasonable starting point for pre-fit errors by setting it to the absolute sigma
1993 // Mostly useful for pre-fit plotting.
1994 gamma.setError(sigmaRel);
1995 } else if( type == Constraint::Poisson ) {
1996
1997 double tau = 1/sigmaRel/sigmaRel; // this is correct Poisson equivalent to a Gaussian with mean 1 and stdev sigma
1998
1999 // Make nominal "observed" value
2000 RooRealVar constrNom(nomName.c_str(), nomName.c_str(), tau);
2001 constrNom.setMin(0);
2002 constrNom.setConstant( true );
2003
2004 // Make the scaling term
2005 std::string scalingName = string(gamma.GetName()) + "_tau";
2006 RooConstVar poissonScaling( scalingName.c_str(), scalingName.c_str(), tau);
2007
2008 // Make mean for scaled Poisson
2009 RooProduct constrMean( poisMeanName.c_str(), poisMeanName.c_str(), RooArgSet(gamma, poissonScaling) );
2010 //proto->import( constrSigma, RecycleConflictNodes() );
2011 //proto->import( constrSigma );
2012
2013 // Type 2 : RooPoisson
2014 RooPoisson pois(constrName.c_str(), constrName.c_str(), constrNom, constrMean);
2015 pois.setNoRounding(true);
2016 proto->import( pois, RecycleConflictNodes() );
2017
2018 if (std::string(gamma.GetName()).find("gamma_stat") != std::string::npos) {
2019 // Give reasonable starting point for pre-fit errors.
2020 // Mostly useful for pre-fit plotting.
2021 gamma.setError(sigmaRel);
2022 }
2023
2024 } else {
2025
2026 std::cout << "Error: Did not recognize Stat Error constraint term type: "
2027 << type << " for : " << paramHist.GetName() << std::endl;
2028 throw hf_exc();
2029 }
2030
2031 // If the sigma value is less
2032 // than a supplied threshold,
2033 // set the variable to constant
2034 if( sigmaRel < minSigma ) {
2035 cxcoutW(HistFactory) << "Warning: Bin " << i << " = " << sigmaRel
2036 << " and is < " << minSigma
2037 << ". Setting: " << gamma.GetName() << " to constant"
2038 << std::endl;
2039 gamma.setConstant(true);
2040 }
2041
2042 constraintTermNames.push_back( constrName );
2043 ConstraintTerms.add( *proto->pdf(constrName) );
2044
2045 // Add the "observed" value to the
2046 // list of global observables:
2047 RooArgSet* globalSet = const_cast<RooArgSet*>(proto->set("globalObservables"));
2048
2049 RooRealVar* nomVarInWorkspace = proto->var(nomName);
2050 if( ! globalSet->contains(*nomVarInWorkspace) ) {
2051 globalSet->add( *nomVarInWorkspace );
2052 }
2053
2054 } // end loop over parameters
2055
2056 return ConstraintTerms;
2057
2058}
2059
2060} // namespace RooStats
2061} // namespace HistFactory
2062
#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
size_t size(const MatrixT &matrix)
retrieve the size of a square matrix
#define cxcoutI(a)
Definition: RooMsgService.h:89
#define cxcoutW(a)
Definition: RooMsgService.h:97
#define cxcoutF(a)
#define cxcoutE(a)
#define cxcoutP(a)
Definition: RooMsgService.h:93
#define ClassImp(name)
Definition: Rtypes.h:375
#define R__ASSERT(e)
Definition: TError.h:118
Option_t Option_t TPoint TPoint const char GetTextMagnitude GetFillStyle GetLineColor GetLineWidth GetMarkerStyle GetTextAlign GetTextColor GetTextSize void data
Option_t Option_t TPoint TPoint const char GetTextMagnitude GetFillStyle GetLineColor GetLineWidth GetMarkerStyle GetTextAlign GetTextColor GetTextSize void char Point_t Rectangle_t WindowAttributes_t Float_t Float_t Float_t Int_t Int_t UInt_t UInt_t Rectangle_t Int_t Int_t Window_t TString Int_t GCValues_t GetPrimarySelectionOwner GetDisplay GetScreen GetColormap GetNativeEvent const char const char dpyName wid window const char font_name cursor keysym reg const char only_if_exist regb h Point_t winding char text const char depth char const char Int_t count const char ColorStruct_t color const char filename
Option_t Option_t TPoint TPoint const char GetTextMagnitude GetFillStyle GetLineColor GetLineWidth GetMarkerStyle GetTextAlign GetTextColor GetTextSize void char Point_t Rectangle_t WindowAttributes_t Float_t Float_t Float_t Int_t Int_t UInt_t UInt_t Rectangle_t result
Option_t Option_t TPoint TPoint const char GetTextMagnitude GetFillStyle GetLineColor GetLineWidth GetMarkerStyle GetTextAlign GetTextColor GetTextSize void char Point_t Rectangle_t WindowAttributes_t Float_t Float_t Float_t Int_t Int_t UInt_t UInt_t Rectangle_t Int_t Int_t Window_t TString Int_t GCValues_t GetPrimarySelectionOwner GetDisplay GetScreen GetColormap GetNativeEvent const char const char dpyName wid window const char font_name cursor keysym reg const char only_if_exist regb h Point_t winding char text const char depth char const char Int_t count const char ColorStruct_t color const char Pixmap_t Pixmap_t PictureAttributes_t attr const char char ret_data h unsigned char height h Atom_t Int_t ULong_t ULong_t unsigned char prop_list Atom_t Atom_t Atom_t Time_t type
char name[80]
Definition: TGX11.cxx:110
float xmin
Definition: THbookFile.cxx:95
float xmax
Definition: THbookFile.cxx:95
char * Form(const char *fmt,...)
Formats a string in a circular formatting buffer.
Definition: TString.cxx:2468
const char * proto
Definition: civetweb.c:17502
A class which maps the current values of a RooRealVar (or a set of RooRealVars) to one of a number of...
Definition: ParamHistFunc.h:24
const RooArgSet * get(Int_t masterIdx) const
Definition: ParamHistFunc.h:46
void setConstant(bool constant)
const RooArgList & paramList() const
Definition: ParamHistFunc.h:34
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)
The PiecewiseInterpolation is a class that can morph distributions into each other,...
RooAbsArg is the common abstract base class for objects that represent a value and a "shape" in RooFi...
Definition: RooAbsArg.h:72
void Print(Option_t *options=nullptr) const override
Print the object to the defaultPrintStream().
Definition: RooAbsArg.h:324
void setAttribute(const Text_t *name, bool value=true)
Set (default) or clear a named boolean attribute of this object.
Definition: RooAbsArg.cxx:245
bool contains(const RooAbsArg &var) const
Check if collection contains an argument with the same name as var.
bool empty() const
Int_t getSize() const
Return the number of elements in the collection.
virtual bool add(const RooAbsArg &var, bool silent=false)
Add the specified argument to list.
Storage_t::size_type size() const
RooAbsData is the common abstract base class for binned and unbinned datasets.
Definition: RooAbsData.h:61
RooAbsRealLValue is the common abstract base class for objects that represent a real value that may a...
void setConstant(bool value=true)
RooNumIntConfig * specialIntegratorConfig() const
Returns the specialized integrator configuration for this RooAbsReal.
virtual void forceNumInt(bool flag=true)
Definition: RooAbsReal.h:186
RooArgList is a container object that can hold multiple RooAbsArg objects.
Definition: RooArgList.h:22
RooArgSet is a container object that can hold multiple RooAbsArg objects.
Definition: RooArgSet.h:56
RooBinWidthFunction is a class that returns the bin width (or volume) given a RooHistFunc.
Class RooBinning is an implements RooAbsBinning in terms of an array of boundary values,...
Definition: RooBinning.h:27
RooCategory is an object to represent discrete states.
Definition: RooCategory.h:28
bool 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:26
The RooDataHist is a container class to hold N-dimensional binned data.
Definition: RooDataHist.h:39
RooDataSet is a container class to hold unbinned data.
Definition: RooDataSet.h:57
void add(const RooArgSet &row, double weight, double weightError)
Add one ore more rows of data.
RooFitResult is a container class to hold the input and output of a PDF fit to a dataset.
Definition: RooFitResult.h:40
Plain Gaussian p.d.f.
Definition: RooGaussian.h:24
Switches the message service to a different level while the instance is alive.
Definition: RooHelpers.h:42
RooHistFunc implements a real-valued function sampled from a multidimensional histogram.
Definition: RooHistFunc.h:31
static RooMsgService & instance()
Return reference to singleton instance.
bool isActive(const RooAbsArg *self, RooFit::MsgTopic facility, RooFit::MsgLevel level)
Check if logging is active for given object/topic/RooFit::MsgLevel combination.
RooCategory & method2D()
RooCategory & methodND()
RooCategory & method1D()
Poisson pdf.
Definition: RooPoisson.h:19
void setNoRounding(bool flag=true)
Switch off/on rounding of x to the nearest integer.
Definition: RooPoisson.h:37
A RooProduct represents the product of a given set of RooAbsReal objects.
Definition: RooProduct.h:29
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:40
void setVal(double value) override
Set value of variable to 'value'.
Definition: RooRealVar.cxx:254
void setMin(const char *name, double value)
Set minimum of name range to given value.
Definition: RooRealVar.cxx:460
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:30
std::vector< RooStats::HistFactory::Data > & GetAdditionalData()
retrieve vector of additional data objects
Definition: Channel.h:65
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:72
RooStats::HistFactory::Data & GetData()
get data object
Definition: Channel.h:59
std::vector< RooStats::HistFactory::Sample > & GetSamples()
get vector of samples for this channel
Definition: Channel.h:77
std::string GetName() const
get name of channel
Definition: Channel.h:43
This class provides helper functions for creating likelihood models from histograms.
std::unique_ptr< RooProduct > CreateNormFactor(RooWorkspace *proto, std::string &channel, std::string &sigmaEpsilon, Sample &sample, bool doRatio)
RooWorkspace * MakeSingleChannelModel(Measurement &measurement, Channel &channel)
RooWorkspace * MakeSingleChannelWorkspace(Measurement &measurement, Channel &channel)
RooArgList createStatConstraintTerms(RooWorkspace *proto, std::vector< std::string > &constraintTerms, ParamHistFunc &paramHist, const TH1 *uncertHist, Constraint::Type type, double minSigma)
std::unique_ptr< TH1 > MakeScaledUncertaintyHist(const std::string &Name, std::vector< std::pair< const TH1 *, std::unique_ptr< TH1 > > > const &HistVec) const
void SetFunctionsToPreprocess(std::vector< std::string > lines)
RooHistFunc * MakeExpectedHistFunc(const TH1 *hist, RooWorkspace *proto, std::string prefix, const RooArgList &observables) const
Create the nominal hist function from hist, and register it in the workspace.
TH1 * MakeAbsolUncertaintyHist(const std::string &Name, const TH1 *Hist)
RooDataSet * MergeDataSets(RooWorkspace *combined, std::vector< std::unique_ptr< RooWorkspace > > &wspace_vec, std::vector< std::string > const &channel_names, std::string const &dataSetName, RooArgList const &obsList, RooCategory *channelCat)
static void ConfigureWorkspaceForMeasurement(const std::string &ModelName, RooWorkspace *ws_single, Measurement &measurement)
void MakeTotalExpected(RooWorkspace *proto, const std::string &totName, const std::vector< RooProduct * > &sampleScaleFactors, std::vector< std::vector< RooAbsArg * > > &sampleHistFuncs) const
void ConfigureHistFactoryDataset(RooDataSet &obsData, TH1 const &nominal, RooWorkspace &proto, std::vector< std::string > const &obsNameVec)
static void PrintCovarianceMatrix(RooFitResult *result, RooArgSet *params, std::string filename)
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)
RooArgList createObservables(const TH1 *hist, RooWorkspace *proto) const
Create observables of type RooRealVar. Creates 1 to 3 observables, depending on the type of the histo...
RooWorkspace * MakeCombinedModel(std::vector< std::string >, std::vector< std::unique_ptr< RooWorkspace > > &)
The RooStats::HistFactory::Measurement class can be used to construct a model by combining multiple R...
Definition: Measurement.h:31
std::map< std::string, double > & GetGammaSyst()
Definition: Measurement.h:122
std::map< std::string, double > & GetLogNormSyst()
Definition: Measurement.h:124
std::map< std::string, double > & GetNoSyst()
Definition: Measurement.h:125
std::vector< std::string > & GetPOIList()
get vector of PoI names
Definition: Measurement.h:51
std::map< std::string, double > & GetUniformSyst()
Definition: Measurement.h:123
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:80
std::vector< std::string > GetPreprocessFunctions() const
Returns a list of defined preprocess function expressions.
double GetLumi()
retrieve integrated luminosity
Definition: Measurement.h:89
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
const std::string & GetName() const
Definition: Systematics.h:56
std::string GetName() const
get name of sample
Definition: Sample.h:82
std::vector< RooStats::HistFactory::NormFactor > & GetNormFactorList()
Definition: Sample.h:109
*Un*constrained bin-by-bin variation of affected histogram.
Definition: Systematics.h:270
const TH1 * GetInitialShape() const
Definition: Systematics.h:286
Constrained bin-by-bin variation of affected histogram.
Definition: Systematics.h:221
Constraint::Type GetConstraintType() const
Definition: Systematics.h:260
const TH1 * GetErrorHist() const
Definition: Systematics.h:252
Constraint::Type GetConstraintType() const
Definition: Systematics.h:384
ModelConfig is a simple class that holds configuration information specifying how a model should be u...
Definition: ModelConfig.h:30
void GuessObsAndNuisance(const RooAbsData &data, bool printModelConfig=true)
Makes sensible guesses of observables, parameters of interest and nuisance parameters if one or multi...
Definition: ModelConfig.cxx:66
virtual void SetParametersOfInterest(const RooArgSet &set)
Specify parameters of interest.
Definition: ModelConfig.h:100
const RooArgSet * GetObservables() const
get RooArgSet for observables (return nullptr if not existing)
Definition: ModelConfig.h:246
The RooWorkspace is a persistable container for RooFit projects.
Definition: RooWorkspace.h:43
TObject * obj(RooStringView name) const
Return any type of object (RooAbsArg, RooAbsData or generic object) with given name)
void Print(Option_t *opts=nullptr) const override
Print contents of the workspace.
RooAbsPdf * pdf(RooStringView name) const
Retrieve p.d.f (RooAbsPdf) with given name. A null pointer is returned if not found.
RooArgSet allVars() const
Return set with all variable objects.
bool 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.
bool saveSnapshot(RooStringView, const char *paramNames)
Save snapshot of values and attributes (including "Constant") of given parameters.
bool importClassCode(const char *pat="*", bool doReplace=false)
Inport code of all classes in the workspace that have a class name that matches pattern 'pat' and whi...
RooFactoryWSTool & factory()
Return instance to factory tool.
RooRealVar * var(RooStringView name) const
Retrieve real-valued variable (RooRealVar) with given name. A null pointer is returned if not found.
RooAbsData * data(RooStringView name) const
Retrieve dataset (binned or unbinned) with given name. A null pointer is returned if not found.
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...
bool loadSnapshot(const char *name)
Load the values and attributes of the parameters in the snapshot saved with the given name.
bool defineSet(const char *name, const RooArgSet &aset, bool importMissing=false)
Define a named RooArgSet with given constituents.
const Double_t * GetArray() const
Definition: TArrayD.h:43
Class to manage histogram axis.
Definition: TAxis.h:30
Bool_t IsVariableBinSize() const
Definition: TAxis.h:137
const char * GetTitle() const override
Returns title of object.
Definition: TAxis.h:130
virtual Double_t GetBinCenter(Int_t bin) const
Return center of bin.
Definition: TAxis.cxx:478
const TArrayD * GetXbins() const
Definition: TAxis.h:131
Double_t GetXmax() const
Definition: TAxis.h:135
Double_t GetXmin() const
Definition: TAxis.h:134
Int_t GetNbins() const
Definition: TAxis.h:121
TH1 is the base class of all histogram classes in ROOT.
Definition: TH1.h:58
TAxis * GetZaxis()
Definition: TH1.h:324
virtual Int_t GetNbinsY() const
Definition: TH1.h:296
virtual Double_t GetBinError(Int_t bin) const
Return value of error associated to bin number bin.
Definition: TH1.cxx:8929
virtual Int_t GetNbinsZ() const
Definition: TH1.h:297
virtual void Reset(Option_t *option="")
Reset this histogram: contents, errors, etc.
Definition: TH1.cxx:7091
TAxis * GetXaxis()
Definition: TH1.h:322
virtual Int_t GetNbinsX() const
Definition: TH1.h:295
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:9072
TAxis * GetYaxis()
Definition: TH1.h:323
Bool_t IsBinUnderflow(Int_t bin, Int_t axis=0) const
Return true if the bin is underflow.
Definition: TH1.cxx:5178
Bool_t IsBinOverflow(Int_t bin, Int_t axis=0) const
Return true if the bin is overflow.
Definition: TH1.cxx:5146
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:9088
virtual Double_t GetBinContent(Int_t bin) const
Return content of bin number bin.
Definition: TH1.cxx:5025
TObject * Clone(const char *newname="") const override
Make a complete copy of the underlying object.
Definition: TH1.cxx:2727
virtual void SetTitle(const char *title="")
Set the title of the TNamed.
Definition: TNamed.cxx:164
const char * GetName() const override
Returns name of object.
Definition: TNamed.h:47
const char * GetTitle() const override
Returns title of object.
Definition: TNamed.h:48
Mother of all ROOT objects.
Definition: TObject.h:41
virtual const char * ClassName() const
Returns name of class to which the object belongs.
Definition: TObject.cxx:207
virtual void Error(const char *method, const char *msgfmt,...) const
Issue error message.
Definition: TObject.cxx:970
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:136
TString & ReplaceAll(const TString &s1, const TString &s2)
Definition: TString.h:693
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:2357
RooCmdArg RecycleConflictNodes(bool flag=true)
RooCmdArg Rename(const char *suffix)
RooCmdArg Conditional(const RooArgSet &pdfSet, const RooArgSet &depSet, bool depsAreCond=false)
RooCmdArg Index(RooCategory &icat)
RooCmdArg WeightVar(const char *name="weight", bool reinterpretAsWeight=false)
RooCmdArg Import(const char *state, TH1 &histo)
double beta(double x, double y)
Calculates the beta function.
double gamma(double x)
VecExpr< UnaryOp< Sqrt< T >, VecExpr< A, T, D >, T >, T, D > sqrt(const VecExpr< A, T, D > &rhs)
The namespace RooFit contains mostly switches that change the behaviour of functions of PDFs (or othe...
Definition: Common.h:18
@ HistFactory
Definition: RooGlobalFunc.h:64
@ ObjectHandling
Definition: RooGlobalFunc.h:63
Namespace for the RooStats classes.
Definition: Asimov.h:19
static constexpr double gauss
const char * Name
Definition: TXMLSetup.cxx:67