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