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