Getting number of events (or fraction of events) in a range

[sarahwferg]

Hello,

I am trying to get the number (or fraction) of signal events vs background events in a range. I have tried a few different ways.

First off here is my code and plot:

#ifndef __CINT__
#include "RooGlobalFunc.h"
#endif
#include "RooRealVar.h"
#include "RooDataSet.h"
#include "RooDataHist.h"
#include "RooGaussian.h"
#include "TCanvas.h"
#include "RooPlot.h"
#include "TTree.h"
#include "TH1D.h"
#include "TRandom.h"

using namespace RooFit ;


void roofitmacro()
{
  //Import histo from file
  TFile f("analysis_hist.root","READ");

  TH1D *kmass = new TH1D("kmass","kaon mass with fit",100 , 450, 600);
  kmass = (TH1D*)f.Get("hkmass");

  //create roofit variable
        RooRealVar x("x","x",450,600) ;
  //create roofit gaussian
        RooRealVar mean("mean","mean",498,490,510) ;
        RooRealVar sigma("sigma","sigma",9,7,11) ;
        RooRealVar sig_yield("sig_yield","signal yield",1000000,0,3000000);
        RooGaussian gauss("gauss","signal p.d.f.",x,mean,sigma) ;
  //create roofit quadratic
        RooRealVar c1("c1","coefficient of x term",400, 300, 600.) ;
        RooRealVar c2("c2","coefficient of x^2 term",-0.4, -1, 0.1) ;
        RooRealVar bkg_yield("bkg_yield","background yield",3000000,0,4000000);
        RooPolynomial bkg("bkg","background", x, RooArgList(c1,c2)) ;
  //create composite model model(x) = sig_yield*gauss(x) + bkg_yield*bkg(x)
        RooAddPdf model("model","model",RooArgList(gauss,bkg),RooArgList(sig_yield,bkg_yield)) ;
  //imposrt data from hist
        RooDataHist data("data","dataset with K0 mass",x,kmass) ;
  //fit model to data
        model.fitTo(data,Extended(),Save());
  //roofit's version of a canvas:
        RooPlot* xframe = x.frame(Title("Kaon Mass with Gaussian Fit and Quadradic Background")) ;
  //plot
        data.plotOn(xframe) ;
        model.paramOn(xframe, Layout(0.5, 0.9, 0.9) );
        model.plotOn(xframe) ;
        model.plotOn(xframe,Components(bkg),LineColor(kGreen));
        xframe->Draw();

  //create canvas to save it on
        TCanvas *c = new TCanvas("c","c",1);
        xframe->Draw();
        c->Print("PLOTS/kmass_roofit.pdf");
}

I know I can do sig_yield.getVal() and it will give me the total num of events (233141) but I want the number of events in a 2sigma range around the peak.

I have tried using createIntegral(x) and createIntegral(x,Range("signal")) :

cout << "nevents sig total = " << sig_yield.getVal() << endl;
        double xmin = mean.getVal() - 2*sigma.getVal();
        double xmax = mean.getVal() + 2*sigma.getVal();
        x.setRange("signal",xmin,xmax);
        RooAbsReal *Imodel = model.createIntegral(x);
        cout << "Imodel = " << Imodel->getVal() << endl;
        RooAbsReal *Imodel_sig = model.createIntegral(x,Range("signal"));
        cout << "Imodel signal = " << Imodel_sig->getVal() << endl;
        RooAbsReal *Ibkg = bkg.createIntegral(x);
        cout << "Ibkg = " << Ibkg->getVal() << endl;
        RooAbsReal *Ibkg_sig = bkg.createIntegral(x,Range("signal"));
        cout << "Ibkg signal = " << Ibkg_sig->getVal() << endl;
        RooAbsReal *Isig = gauss.createIntegral(x);
        cout << "Isig = " << Isig->getVal() << endl;
        RooAbsReal *Isig_sig = gauss.createIntegral(x,Range("signal"));
        cout << "Isig signal = " << Isig_sig->getVal() << endl;
        cout << "Ibkg + Isig = " << (Ibkg->getVal() + Isig->getVal() ) << endl;
        cout << "Ibkg_sig + Isig_sig = " << (Ibkg_sig->getVal() + Isig_sig->getVal() ) << endl;

which outputs:

nevents sig total = 233141
[#1] INFO:Eval -- RooRealVar::setRange(x) new range named 'signal' created with bounds [480.094,514.946]
Imodel = 1.2572e+07
Imodel signal = 2.96325e+06
Ibkg = 1.44151e+07
Ibkg signal = 3.39767e+06
Isig = 21.8405
Isig signal = 20.8467
Ibkg + Isig = 1.44151e+07
Ibkg_sig + Isig_sig = 3.39769e+06

but these outputs do not make sense.

1. Background should not be larger than the total
2. Number of signal events is 233141 so integrating over gaussian signal should get me something close to that, not roughly 21. Maybe the ~21 number is a percentage? but if so why is the background integral not a percentage too?
3. Background yield is about 1.5 X 10^6 so why is background integral 1.4 X 10^7, an order of magnitude higher?
4. Why do I get more events if I add signal and background integrals than if I just integrate over the model? Shouldn’t I get roughly the same number?

So I do no know what I am misunderstanding or doing wrong.

Thanks in advance for your help.

Sarah

[moneta]

Hi,

I have experienced in the past fe issues with createIntegral in RooFit. I would try tp convert the
RooAbsPdf to a TF1 using the asTF function and then compute the integral in the range using the TF1::Integral function and eventually also its error using TF1::IntegralError, by passing the covariance matrix of the fitted function parameters

Best Regards

Lorenzo

[sarahwferg]

Hi Thanks for your reply.

I could not get your suggestion to work, however I figured out how to use the createIntegral() function:

        x.setRange("signal",mean.getVal()-2*sigma.getVal(),mean.getVal()+2*sigma.getVal());
        RooAbsReal* fsigregion_model = model.createIntegral(x,NormSet(x),Range("signal")); //The "NormSet(x)" normalizes it to the total number of events to give you the fraction n_signal_region_events/n_total_events
        RooAbsReal* fsigregion_bkg = bkg.createIntegral(x,NormSet(x),Range("signal")); 
        Double_t nsigevents = fsigregion_model->getVal()*(sig_yield.getVal()+bkg_yield.getVal())-fsigregion_bkg->getVal()*bkg_yield.getVal(); 
        Double_t fsig = nsigevents/(fsigregion_model->getVal()*(sig_yield.getVal()+bkg_yield.getVal()));

which outputs:

fsigregion_model = 0.262938 //This is the fraction of total events in the signal region
fsigregion_bkg = 0.182236 //This is the fraction of background events in the signal region
nsigevents = 232211//This is the number of signal events in the signal region
fsig = 0.379349 //This is the fraction of signal event in the signal region.

Which seems reasonable/correct.