doc/master/rf709__BarlowBeeston_8C.html

#include "RooRealVar.h"

#include "RooGaussian.h"

#include "RooUniform.h"

#include "RooDataSet.h"

#include "RooDataHist.h"

#include "RooHistFunc.h"

#include "RooRealSumPdf.h"

#include "RooProduct.h"

#include "RooConstVar.h"

#include "RooPoisson.h"

#include "RooProdPdf.h"

#include "RooFitResult.h"

#include "RooPlot.h"

#include "RooMsgService.h"

#include "RooStats/HistFactory/ParamHistFunc.h"


#include "TCanvas.h"

#include "TPaveText.h"


#include <cmath>

#include <iostream>

#include <memory>


RooArgSet createConstraints(RooArgList const &gammas, RooArgList const &samples, RooArgSet &globalObservables)

{

   RooArgSet constraints;

   for (std::size_t i = 0; i < gammas.size(); ++i) {

      auto *gamma = static_cast<RooRealVar *>(&gammas[i]);


      auto nomName = "nom_" + std::string(gamma->GetName());

      auto tauName = std::string(gamma->GetName()) + "_tau";

      auto meanName = std::string(gamma->GetName()) + "_mean";

      auto poissonName = std::string(gamma->GetName()) + "_constraint";


      double nomVal = 0.0;

      for (auto *func : static_range_cast<RooHistFunc *>(samples)) {

         nomVal += func->dataHist().weight(i);

      }


      // Like in HistFactory, the gamma of an empty bin is fixed and not

      // constrained. Otherwise, it would be a flat direction of the likelihood.

      // The other gammas get the relative Poisson uncertainty as initial error.

      if (nomVal == 0.0) {

         gamma->setConstant(true);

         continue;

      }

      gamma->setError(1.0 / std::sqrt(nomVal));


      auto nominal = std::make_unique<RooRealVar>(nomName.c_str(), nomName.c_str(), nomVal, 0.0, RooNumber::infinity());

      nominal->setConstant(true);

      globalObservables.add(*nominal);


      auto tau = std::make_unique<RooConstVar>(tauName.c_str(), tauName.c_str(), nominal->getVal());


      auto mean = std::make_unique<RooProduct>(meanName.c_str(), meanName.c_str(), *gamma, *tau);

      mean->addOwnedComponents(std::move(tau));


      auto poisson = std::make_unique<RooPoisson>(poissonName.c_str(), poissonName.c_str(), *nominal, *mean);

      poisson->addOwnedComponents(std::move(nominal));

      poisson->addOwnedComponents(std::move(mean));


      constraints.addOwned(std::move(poisson));

   }

   return constraints;

}


void rf709_BarlowBeeston()

{

   using namespace RooFit;


   RooMsgService::instance().getStream(1).removeTopic(NumIntegration);


   // First, construct a likelihood model with a Gaussian signal on top of a uniform background

   RooRealVar x("x", "x", -20, 20);

   x.setBins(25);


   RooRealVar meanG("meanG", "meanG", 1, -10, 10);

   RooRealVar sigG("sigG", "sigG", 1.5, -10, 10);

   RooGaussian g("g", "Gauss", x, meanG, sigG);

   RooUniform u("u", "Uniform", x);


   // Generate the data to be fitted

   std::unique_ptr<RooDataSet> sigData(g.generate(x, 50));

   std::unique_ptr<RooDataSet> bkgData(u.generate(x, 1000));


   RooDataSet sumData("sumData", "Gauss + Uniform", x);

   sumData.append(*sigData);

   sumData.append(*bkgData);


   // Make histogram templates for signal and background.

   // Let's take a signal distribution with low statistics and a more accurate

   // background distribution.

   // Normally, these come from Monte Carlo simulations, but we will just generate them.

   std::unique_ptr<RooDataHist> dh_sig(g.generateBinned(x, 50));

   std::unique_ptr<RooDataHist> dh_bkg(u.generateBinned(x, 10000));


   // ***** Case 0 - 'Rigid templates' *****


   // Construct histogram shapes for signal and background

   RooHistFunc p_h_sig("p_h_sig", "p_h_sig", x, *dh_sig);

   RooHistFunc p_h_bkg("p_h_bkg", "p_h_bkg", x, *dh_bkg);


   // Construct scale factors for adding the two distributions

   RooRealVar Asig0("Asig", "Asig", 1, 0.01, 5000);

   RooRealVar Abkg0("Abkg", "Abkg", 1, 0.01, 5000);


   // Construct the sum model

   RooRealSumPdf model0("model0", "model0", {p_h_sig, p_h_bkg}, {Asig0, Abkg0}, true);


   // ***** Case 1 - 'Barlow Beeston' *****


   // Construct parameterized histogram shapes for signal and background

   RooArgList gammas_sig1{ParamHistFunc::createParamSet("p_ph_sig_gamma", x.numBins(), 0, 10)};

   RooArgList gammas_bkg1{ParamHistFunc::createParamSet("p_ph_bkg_gamma", x.numBins(), 0, 10)};


   ParamHistFunc p_phf_sig1{"p_phf_sig1", "p_phf_sig1", x, gammas_sig1};

   ParamHistFunc p_phf_bkg1{"p_phf_bkg1", "p_phf_bkg1", x, gammas_bkg1};


   RooProduct p_p_sig1{"p_p_sig1", "p_p_sig1", p_h_sig, p_phf_sig1};

   RooProduct p_p_bkg1{"p_p_bkg1", "p_p_bkg1", p_h_bkg, p_phf_bkg1};


   RooRealVar Asig1("Asig", "Asig", 1, 0.01, 5000);

   RooRealVar Abkg1("Abkg", "Abkg", 1, 0.01, 5000);


   // Construct the sum of these

   RooRealSumPdf model1_tmp("sp_ph", "sp_ph", {p_p_sig1, p_p_bkg1}, {Asig1, Abkg1}, true);


   // Construct the subsidiary poisson measurements constraining the histogram parameters

   // These ensure that the bin contents of the histograms are only allowed to vary within

   // the statistical uncertainty of the Monte Carlo.

   RooArgSet globs1;

   RooArgSet constraints_sig1{createConstraints(gammas_sig1, p_h_sig, globs1)};

   RooArgSet constraints_bkg1{createConstraints(gammas_bkg1, p_h_bkg, globs1)};


   // Construct the joint model with template PDFs and constraints

   RooProdPdf model1("model1", "model1", {model1_tmp, constraints_sig1, constraints_bkg1});


   // ***** Case 2 - 'Barlow Beeston' light (one parameter per bin for all samples) *****


   // Construct the histogram shapes, using the same parameters for signal and background.

   // There is one parameter per bin, such as `p_ph_gamma_bin_0`. This allows bin 0 to

   // fluctuate up and down. The SAME ParamHistFunc scales both the signal and the background

   // histogram, so the bins fluctuate synchronously. This reduces the number of parameters.

   RooArgList gammas2{ParamHistFunc::createParamSet("p_ph_gamma", x.numBins(), 0, 10)};


   ParamHistFunc p_phf2{"p_phf2", "p_phf2", x, gammas2};


   RooProduct p_p_sig2{"p_p_sig2", "p_p_sig2", p_h_sig, p_phf2};

   RooProduct p_p_bkg2{"p_p_bkg2", "p_p_bkg2", p_h_bkg, p_phf2};


   RooRealVar Asig2("Asig", "Asig", 1, 0.01, 5000);

   RooRealVar Abkg2("Abkg", "Abkg", 1, 0.01, 5000);


   // As before, construct the sum of signal2 and background2

   RooRealSumPdf model2_tmp("sp_ph", "sp_ph", {p_p_sig2, p_p_bkg2}, {Asig2, Abkg2}, true);


   // Construct the subsidiary poisson measurements constraining the statistical fluctuations

   RooArgSet globs2;

   RooArgSet constraints2{createConstraints(gammas2, {p_h_sig, p_h_bkg}, globs2)};


   // Construct the joint model

   RooProdPdf model2("model2", "model2", {model2_tmp, constraints2});


   // ************ Fit all models to data and plot *********************


   std::unique_ptr<RooFitResult> result0{model0.fitTo(sumData, PrintLevel(0), Save())};

   std::unique_ptr<RooFitResult> result1{model1.fitTo(sumData, PrintLevel(0), Save(), GlobalObservables(globs1))};

   std::unique_ptr<RooFitResult> result2{model2.fitTo(sumData, PrintLevel(0), Save(), GlobalObservables(globs2))};


   TCanvas *can = new TCanvas("can", "", 1500, 600);

   can->Divide(3, 1);


   TPaveText pt(-19.5, 1, -2, 25);

   pt.SetFillStyle(0);

   pt.SetBorderSize(0);


   can->cd(1);

   auto frame = x.frame(Title("No template uncertainties"));

   // Plot data to enable automatic determination of model0 normalisation:

   sumData.plotOn(frame);

   model0.plotOn(frame, LineColor(kBlue), VisualizeError(*result0));

   // Plot data again to show it on top of model0 error bands:

   sumData.plotOn(frame);

   // Plot model components

   model0.plotOn(frame, LineColor(kBlue));

   model0.plotOn(frame, Components(p_h_sig), LineColor(kAzure));

   model0.plotOn(frame, Components(p_h_bkg), LineColor(kRed));

   model0.paramOn(frame);


   sigData->plotOn(frame, MarkerColor(kBlue));

   frame->Draw();


   for (auto text : {"No template uncertainties", "are taken into account.", "This leads to low errors",

                     "for the parameters A, since", "the only source of errors", "are the data statistics."}) {

      pt.AddText(text);

   }

   pt.DrawClone();


   can->cd(2);

   frame = x.frame(Title("Barlow Beeston for Sig & Bkg separately"));

   sumData.plotOn(frame);

   model1.plotOn(frame, LineColor(kBlue), VisualizeError(*result1));

   // Plot data again to show it on top of error bands:

   sumData.plotOn(frame);

   model1.plotOn(frame, LineColor(kBlue));

   model1.plotOn(frame, Components(p_p_sig1), LineColor(kAzure));

   model1.plotOn(frame, Components(p_p_bkg1), LineColor(kRed));

   model1.paramOn(frame, Parameters({Asig1, Abkg1}));


   sigData->plotOn(frame, MarkerColor(kBlue));

   frame->Draw();


   pt.Clear();

   for (auto text : {"With gamma parameters, the", "signal & background templates", "can adapt to the data.",

                     "Note how the blue signal", "template changes its shape.", "This leads to higher errors",

                     "of the scale parameters A."}) {

      pt.AddText(text);

   }

   pt.DrawClone();


   can->cd(3);

   frame = x.frame(Title("Barlow Beeston light for (Sig+Bkg)"));

   sumData.plotOn(frame);

   model2.plotOn(frame, LineColor(kBlue), VisualizeError(*result2));

   // Plot data again to show it on top of model0 error bands:

   sumData.plotOn(frame);

   model2.plotOn(frame, LineColor(kBlue));

   model2.plotOn(frame, Components(p_p_sig2), LineColor(kAzure));

   model2.plotOn(frame, Components(p_p_bkg2), LineColor(kRed));

   model2.paramOn(frame, Parameters({Asig2, Abkg2}));


   sigData->plotOn(frame, MarkerColor(kBlue));

   frame->Draw();


   pt.Clear();

   for (auto text : {"When signal and background", "template share one gamma para-", "meter per bin, they adapt less.",

                     "The errors of the A parameters", "also shrink slightly."}) {

      pt.AddText(text);

   }

   pt.DrawClone();


   std::cout << "Asig [normal ] = " << Asig0.getVal() << " +/- " << Asig0.getError() << std::endl;

   std::cout << "Asig [BB     ] = " << Asig1.getVal() << " +/- " << Asig1.getError() << std::endl;

   std::cout << "Asig [BBlight] = " << Asig2.getVal() << " +/- " << Asig2.getError() << std::endl;

}

ParamHistFunc.h

g
#define g(i)
Definition RSha256.hxx:105

static_range_cast
ROOT::RRangeCast< T, false, Range_t > static_range_cast(Range_t &&coll)
Definition RooAbsCollection.h:45

RooConstVar.h

RooDataHist.h

RooDataSet.h

RooFitResult.h

RooGaussian.h

RooHistFunc.h

RooMsgService.h

RooPlot.h

RooPoisson.h

RooProdPdf.h

RooProduct.h

RooRealSumPdf.h

RooRealVar.h

RooUniform.h

kRed
@ kRed
Definition Rtypes.h:66

kBlue
@ kBlue
Definition Rtypes.h:66

kAzure
@ kAzure
Definition Rtypes.h:67

TCanvas.h

text
Option_t Option_t TPoint TPoint const char text
Definition TGWin32VirtualXProxy.cxx:71

TPaveText.h

ParamHistFunc
A class which maps the current values of a RooRealVar (or a set of RooRealVars) to one of a number of...
Definition ParamHistFunc.h:24

ParamHistFunc::createParamSet
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.
Definition ParamHistFunc.cxx:286

RooAbsCollection::add
virtual bool add(const RooAbsArg &var, bool silent=false)
Add the specified argument to list.
Definition RooAbsCollection.cxx:440

RooAbsCollection::size
Storage_t::size_type size() const
Definition RooAbsCollection.h:274

RooAbsCollection::addOwned
virtual bool addOwned(RooAbsArg &var, bool silent=false)
Add an argument and transfer the ownership to the collection.
Definition RooAbsCollection.cxx:364

RooAbsPdf::fitTo
RooFit::OwningPtr< RooFitResult > fitTo(RooAbsData &data, CmdArgs_t const &... cmdArgs)
Fit PDF to given dataset.
Definition RooAbsPdf.h:149

RooArgList
RooArgList is a container object that can hold multiple RooAbsArg objects.
Definition RooArgList.h:22

RooArgSet
RooArgSet is a container object that can hold multiple RooAbsArg objects.
Definition RooArgSet.h:24

RooDataSet
Container class to hold unbinned data.
Definition RooDataSet.h:32

RooGaussian
Plain Gaussian p.d.f.
Definition RooGaussian.h:24

RooHistFunc
A real-valued function sampled from a multidimensional histogram.
Definition RooHistFunc.h:29

RooMsgService::instance
static RooMsgService & instance()
Return reference to singleton instance.
Definition RooMsgService.cxx:346

RooMsgService::getStream
StreamConfig & getStream(Int_t id)
Definition RooMsgService.h:161

RooNumber::infinity
static constexpr double infinity()
Return internal infinity representation.
Definition RooNumber.h:25

RooProdPdf
Efficient implementation of a product of PDFs of the form.
Definition RooProdPdf.h:35

RooProduct
Represents the product of a given set of RooAbsReal objects.
Definition RooProduct.h:29

RooRealSumPdf
Implements a PDF constructed from a sum of functions:
Definition RooRealSumPdf.h:23

RooRealVar
Variable that can be changed from the outside.
Definition RooRealVar.h:37

RooUniform
Flat p.d.f.
Definition RooUniform.h:24

TAttFill::SetFillStyle
virtual void SetFillStyle(Style_t fstyle)
Set the fill area style.
Definition TAttFill.h:42

TCanvas
The Canvas class.
Definition TCanvas.h:23

TCanvas::cd
TVirtualPad * cd(Int_t subpadnumber=0) override
Set current canvas & pad.
Definition TCanvas.cxx:732

TObject::DrawClone
virtual TObject * DrawClone(Option_t *option="") const
Draw a clone of this object in the current selected pad with: gROOT->SetSelectedPad(c1).
Definition TObject.cxx:318

TPad::Divide
void Divide(Int_t nx=1, Int_t ny=1, Float_t xmargin=0.01, Float_t ymargin=0.01, Int_t color=0) override
Automatic pad generation by division.
Definition TPad.cxx:1361

TPaveText
A Pave (see TPave) with text, lines or/and boxes inside.
Definition TPaveText.h:21

TPaveText::AddText
virtual TText * AddText(Double_t x1, Double_t y1, const char *label)
Add a new Text line to this pavetext at given coordinates.
Definition TPaveText.cxx:189

TPaveText::Clear
void Clear(Option_t *option="") override
Clear all lines in this pavetext.
Definition TPaveText.cxx:215

TPave::SetBorderSize
virtual void SetBorderSize(Int_t bordersize=4)
Sets the border size of the TPave box and shadow.
Definition TPave.h:79

TVirtualPad::Draw
void Draw(Option_t *option="") override=0
Default Draw method for all objects.

pt
TPaveText * pt
Definition entrylist_figure1.C:7

RooFit::Parameters
RooCmdArg Parameters(const RooArgSet &params)
Definition RooGlobalFunc.cxx:809

RooFit::Save
RooCmdArg Save(bool flag=true)
Definition RooGlobalFunc.cxx:671

RooFit::GlobalObservables
RooCmdArg GlobalObservables(Args_t &&... argsOrArgSet)
Definition RooGlobalFunc.h:321

RooFit::PrintLevel
RooCmdArg PrintLevel(Int_t code)
Definition RooGlobalFunc.cxx:679

RooFit::Components
RooCmdArg Components(Args_t &&... argsOrArgSet)
Definition RooGlobalFunc.h:141

RooFit::LineColor
RooCmdArg LineColor(TColorNumber color)
Definition RooGlobalFunc.cxx:229

RooFit::VisualizeError
RooCmdArg VisualizeError(const RooDataSet &paramData, double Z=1)
Definition RooGlobalFunc.cxx:289

RooFit::MarkerColor
RooCmdArg MarkerColor(TColorNumber color)
Definition RooGlobalFunc.cxx:341

x
Double_t x[n]
Definition legend1.C:17

ROOT::Math::Cephes::gamma
double gamma(double x)
Definition SpecFuncCephes.cxx:339

RooFit::Detail::MathFuncs::poisson
double poisson(double x, double par)
Definition MathFuncs.h:239

RooFit
The namespace RooFit contains mostly switches that change the behaviour of functions of PDFs (or othe...
Definition CodegenImpl.h:73

rf709_BarlowBeeston
Definition rf709_BarlowBeeston.py:1

xmlio::Title
const char * Title
Definition TXMLSetup.cxx:67

RooMsgService::StreamConfig::removeTopic
void removeTopic(RooFit::MsgTopic oldTopic)
Definition RooMsgService.h:122