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