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