rf702_efficiencyfit_2D

Special pdf's: unbinned maximum likelihood fit of an efficiency eff(x) function to a dataset D(x,cut), cut is a category encoding a selection whose efficiency as function of x should be described by eff(x)

Author: Clemens Lange, Wouter Verkerke (C++ version)
This notebook tutorial was automatically generated with ROOTBOOK-izer from the macro found in the ROOT repository on Wednesday, November 30, 2022 at 11:24 AM.

In [1]:
%%cpp -d
#include "RooRealVar.h"
#include "RooDataSet.h"
#include "RooGaussian.h"
#include "RooConstVar.h"
#include "RooCategory.h"
#include "RooEfficiency.h"
#include "RooPolynomial.h"
#include "RooProdPdf.h"
#include "RooFormulaVar.h"
#include "TCanvas.h"
#include "TAxis.h"
#include "TH1.h"
#include "RooPlot.h"
using namespace RooFit;

Arguments are defined.

In [2]:
bool flat = false;

Construct efficiency function e(x,y)

Declare variables x,mean,sigma with associated name, title, initial value and allowed range

In [3]:
RooRealVar x("x", "x", -10, 10);
RooRealVar y("y", "y", -10, 10);

Efficiency function eff(x;a,b)

In [4]:
RooRealVar ax("ax", "ay", 0.6, 0, 1);
RooRealVar bx("bx", "by", 5);
RooRealVar cx("cx", "cy", -1, -10, 10);

RooRealVar ay("ay", "ay", 0.2, 0, 1);
RooRealVar by("by", "by", 5);
RooRealVar cy("cy", "cy", -1, -10, 10);

RooFormulaVar effFunc("effFunc", "((1-ax)+ax*cos((x-cx)/bx))*((1-ay)+ay*cos((y-cy)/by))",
                      RooArgList(ax, bx, cx, x, ay, by, cy, y));

Acceptance state cut (1 or 0)

In [5]:
RooCategory cut("cut", "cutr", { {"accept", 1}, {"reject", 0} });

Construct conditional efficiency pdf E(cut|x,y)

Construct efficiency pdf eff(cut|x)

In [6]:
RooEfficiency effPdf("effPdf", "effPdf", effFunc, cut, "accept");

Generate data (x,y,cut) from a toy model

Construct global shape pdf shape(x) and product model(x,cut) = eff(cut|x)*shape(x) (These are only needed to generate some toy MC here to be used later)

In [7]:
RooPolynomial shapePdfX("shapePdfX", "shapePdfX", x, RooConst(flat ? 0 : -0.095));
RooPolynomial shapePdfY("shapePdfY", "shapePdfY", y, RooConst(flat ? 0 : +0.095));
RooProdPdf shapePdf("shapePdf", "shapePdf", RooArgSet(shapePdfX, shapePdfY));
RooProdPdf model("model", "model", shapePdf, Conditional(effPdf, cut));

Generate some toy data from model

In [8]:
RooDataSet *data = model.generate(RooArgSet(x, y, cut), 10000);
[#0] WARNING:Generation -- RooAcceptReject::ctor(effPdf_Int[]_Norm[cut]) WARNING: performing accept/reject sampling on a p.d.f in 2 dimensions without prior knowledge on maximum value of p.d.f. Determining maximum value by taking 200000 trial samples. If p.d.f contains sharp peaks smaller than average distance between trial sampling points these may be missed and p.d.f. may be sampled incorrectly.
[#0] WARNING:Generation -- RooAcceptReject::ctor(effPdf_Int[]_Norm[cut]): WARNING: 200000 trial samples requested by p.d.f for 2-dimensional accept/reject sampling, this may take some time
input_line_59:2:2: warning: 'data' shadows a declaration with the same name in the 'std' namespace; use '::data' to reference this declaration
 RooDataSet *data = model.generate(RooArgSet(x, y, cut), 10000);
 ^

Fit conditional efficiency pdf to data

Fit conditional efficiency pdf to data

In [9]:
effPdf.fitTo(*data, ConditionalObservables(RooArgSet(x, y)));
input_line_60:2:16: error: reference to 'data' is ambiguous
 effPdf.fitTo(*data, ConditionalObservables(RooArgSet(x, y)));
               ^
input_line_59:2:14: note: candidate found by name lookup is '__cling_N527::data'
 RooDataSet *data = model.generate(RooArgSet(x, y, cut), 10000);
             ^
/usr/include/c++/9/bits/range_access.h:318:5: note: candidate found by name lookup is 'std::data'
    data(initializer_list<_Tp> __il) noexcept
    ^
/usr/include/c++/9/bits/range_access.h:289:5: note: candidate found by name lookup is 'std::data'
    data(_Container& __cont) noexcept(noexcept(__cont.data()))
    ^
/usr/include/c++/9/bits/range_access.h:299:5: note: candidate found by name lookup is 'std::data'
    data(const _Container& __cont) noexcept(noexcept(__cont.data()))
    ^
/usr/include/c++/9/bits/range_access.h:309:5: note: candidate found by name lookup is 'std::data'
    data(_Tp (&__array)[_Nm]) noexcept
    ^

Plot fitted, data efficiency

Make 2D histograms of all data, selected data and efficiency function

In [10]:
TH1 *hh_data_all = data->createHistogram("hh_data_all", x, Binning(8), YVar(y, Binning(8)));
TH1 *hh_data_sel = data->createHistogram("hh_data_sel", x, Binning(8), YVar(y, Binning(8)), Cut("cut==cut::accept"));
TH1 *hh_eff = effFunc.createHistogram("hh_eff", x, Binning(50), YVar(y, Binning(50)));
input_line_61:2:21: error: reference to 'data' is ambiguous
 TH1 *hh_data_all = data->createHistogram("hh_data_all", x, Binning(8), YVar(y, Binning(8)));
                    ^
input_line_59:2:14: note: candidate found by name lookup is '__cling_N527::data'
 RooDataSet *data = model.generate(RooArgSet(x, y, cut), 10000);
             ^
/usr/include/c++/9/bits/range_access.h:318:5: note: candidate found by name lookup is 'std::data'
    data(initializer_list<_Tp> __il) noexcept
    ^
/usr/include/c++/9/bits/range_access.h:289:5: note: candidate found by name lookup is 'std::data'
    data(_Container& __cont) noexcept(noexcept(__cont.data()))
    ^
/usr/include/c++/9/bits/range_access.h:299:5: note: candidate found by name lookup is 'std::data'
    data(const _Container& __cont) noexcept(noexcept(__cont.data()))
    ^
/usr/include/c++/9/bits/range_access.h:309:5: note: candidate found by name lookup is 'std::data'
    data(_Tp (&__array)[_Nm]) noexcept
    ^
input_line_61:3:20: error: reference to 'data' is ambiguous
TH1 *hh_data_sel = data->createHistogram("hh_data_sel", x, Binning(8), YVar(y, Binning(8)), Cut("cut==cut::accept"));
                   ^
input_line_59:2:14: note: candidate found by name lookup is '__cling_N527::data'
 RooDataSet *data = model.generate(RooArgSet(x, y, cut), 10000);
             ^
/usr/include/c++/9/bits/range_access.h:318:5: note: candidate found by name lookup is 'std::data'
    data(initializer_list<_Tp> __il) noexcept
    ^
/usr/include/c++/9/bits/range_access.h:289:5: note: candidate found by name lookup is 'std::data'
    data(_Container& __cont) noexcept(noexcept(__cont.data()))
    ^
/usr/include/c++/9/bits/range_access.h:299:5: note: candidate found by name lookup is 'std::data'
    data(const _Container& __cont) noexcept(noexcept(__cont.data()))
    ^
/usr/include/c++/9/bits/range_access.h:309:5: note: candidate found by name lookup is 'std::data'
    data(_Tp (&__array)[_Nm]) noexcept
    ^

Some adjustment for good visualization

In [11]:
hh_data_all->SetMinimum(0);
hh_data_sel->SetMinimum(0);
hh_eff->SetMinimum(0);
hh_eff->SetLineColor(kBlue);
input_line_62:5:22: error: cannot take the address of an rvalue of type 'EColor'
hh_eff->SetLineColor(kBlue);
                     ^~~~~
Error while creating dynamic expression for:
  hh_eff->SetLineColor(kBlue)

Draw all frames on a canvas

In [12]:
TCanvas *ca = new TCanvas("rf702_efficiency_2D", "rf702_efficiency_2D", 1200, 400);
ca->Divide(3);
ca->cd(1);
gPad->SetLeftMargin(0.15);
hh_data_all->GetZaxis()->SetTitleOffset(1.8);
hh_data_all->Draw("lego");
ca->cd(2);
gPad->SetLeftMargin(0.15);
hh_data_sel->GetZaxis()->SetTitleOffset(1.8);
hh_data_sel->Draw("lego");
ca->cd(3);
gPad->SetLeftMargin(0.15);
hh_eff->GetZaxis()->SetTitleOffset(1.8);
hh_eff->Draw("surf");

return;
input_line_64:2:3: error: use of undeclared identifier 'hh_data_all'
 (hh_data_all->GetZaxis()->SetTitleOffset(1.8))
  ^
Error in <HandleInterpreterException>: Error evaluating expression (hh_data_all->GetZaxis()->SetTitleOffset(1.8))
Execution of your code was aborted.

Draw all canvases

In [13]:
gROOT->GetListOfCanvases()->Draw()