Hi all,
I was wondering if there’s a way to create a RooEfficiency from an efficiency histogram I was given. The purpose of which is to fit my model*efficiency as specified in the RooFit tutorial
Thanks,
Y.
Hi all,
I was wondering if there’s a way to create a RooEfficiency from an efficiency histogram I was given. The purpose of which is to fit my model*efficiency as specified in the RooFit tutorial
Thanks,
Y.
I found my solution:
TH1F* map = createEffHist();
RooHistPdf* hPdf = mapToPdf(h, "hPdf", x);
int interpolationOrder = 1;
bool impDens = kFALSE;
RooDataHist* h = new RooDataHist("tmp", "tmp", argSet, Import(*map, impDens));
RooHistPdf* pdf = new RooHistPdf("name", "desc", argSet, *h, interpolationOrder);
RooCategory cut("cut","cutr") ;
cut.defineType("accept",1) ;
cut.defineType("reject",0) ;
RooEfficiency hEffPdf("hEffPdf", "hEffPdf", *pdf, cut, "accept");
Posted too soon. Apparently there’s a difference between using a RooFormulaVar and a histogram to formulate your efficiency.
It looks like when I use RooFormulaVar, and the efficiency is formulated as a fraction, I do get the right fraction of events generated and passing my efficiency. When I use a histogram, the fraction I get is always 1.0 in areas where efficiency was >0, even if it was 0.1.
You can see the differences in the following images:
Here is the code which I used to produce them:
#include <RooRealVar.h>
#include <RooPlot.h>
#include <RooArgList.h>
#include <RooArgSet.h>
#include <RooDataSet.h>
#include <RooGaussian.h>
#include <RooExponential.h>
#include <RooPolynomial.h>
#include <RooProdPdf.h>
#include <RooFormulaVar.h>
#include <RooCategory.h>
#include <RooEfficiency.h>
#include <RooDataHist.h>
#include <RooHistPdf.h>
#include <TCanvas.h>
#include <TH1F.h>
#define USEHISTEFF false
using namespace RooFit ;
TH1F* createEffHist() {
TH1::SetDefaultSumw2();
double eff = 0.5;
TH1F* h = new TH1F("heff", "heff", 100, -10, 10);
for (int i=50-10; i <= 50+10; ++i) {
h->SetBinContent(i, eff);
h->SetBinError(i, eff*(1.0-eff)/100.0);
}
return h;
}
RooHistPdf* mapToPdf(TH1F* map, const char* pdfName, RooArgSet argSet) {
int interpolationOrder = 1;
bool impDens = kFALSE;
RooDataHist* h = new RooDataHist("tmp", "tmp", argSet, Import(*map, impDens));
RooHistPdf* pdf = new RooHistPdf(pdfName, pdfName, argSet, *h, interpolationOrder);
return pdf;
}
void RooFitEfficiency() {
// Observable
RooRealVar x("x","x",-10,10) ;
// Efficiency function eff(x;a,b)
RooRealVar a("a","a",2.0,0.0,10.0) ;
RooRealVar b("b","b",0.5, 0.0, 1.0) ;
RooRealVar c("c","c",-1,-10,10) ;
RooRealVar lambda("lambda", "#lambda", -0.1, -10, 1.0);
RooExponential shapePdf("shapePdf", "shapePdf", x, lambda);
// Acceptance state cut (1 or 0)
RooCategory cut("cut","cutr") ;
cut.defineType("accept",1) ;
cut.defineType("reject",0) ;
// Construct efficiency p.d.f eff(cut|x)
#if USEHISTEFF
TH1F* h = createEffHist();
RooHistPdf* hPdf = mapToPdf(h, "hPdf", x);
RooFormulaVar effFunc("effFunc", "effFunc", *hPdf);
RooEfficiency effPdf("effPdf", "effPdf", *hPdf, cut, "accept");
#else
// RooFormulaVar effFunc("effFunc","(1-a)+a*cos((x-c)/b)",RooArgList(a,b,c,x)) ;
RooFormulaVar effFunc("effFunc","0.0*(abs(x)>a) + b*(abs(x)<a)",RooArgList(a,b,x)) ;
RooEfficiency effPdf("effPdf","effPdf",effFunc,cut,"accept") ;
#endif
RooProdPdf model("model", "model", shapePdf, Conditional(effPdf, cut));
//**** DATA ****//
RooDataSet* data = model.generate(RooArgSet(x,cut),100000);
// Plot distribution of all events and accepted fraction of events on frame
RooPlot* frame1 = x.frame(Bins(20),Title("Data (all, accepted)")) ;
data->plotOn(frame1) ;
data->plotOn(frame1,Cut("cut==cut::accept"),MarkerColor(kRed),LineColor(kRed)) ;
model.plotOn(frame1,LineColor(kBlue));
// Plot accept/reject efficiency on data overlay fitted efficiency curve
RooPlot* frame2 = x.frame(Bins(20),Title("Fitted efficiency")) ;
data->plotOn(frame2,Efficiency(cut)) ; // needs ROOT version >= 5.21
effPdf.plotOn(frame2, LineColor(kRed)) ;
effFunc.plotOn(frame2, LineColor(kBlue)) ;
TCanvas* can1 = new TCanvas("can1","can1", 0, 0, 800, 400);
can1->Divide(2);
can1->cd(1); frame1->Draw();
can1->cd(2); frame2->Draw();
}