Dear all,
I am seriously struggling to get my head around some problems I am having with the most basic of the tasks, fitting a straight line.
I have generated a minimal working example, I just want to fit a straight line in the sidebands. I want to do it with and extended p.d.f…
1- Why if I fit it with or without an extended p.d.f. I get a different answer for the slope of the line?
2- Why if I do not reduce the observable to the minimum I get a different answer from the fit?
3- Ultimately what is the best way to compute the expected number of events in a region between the sidebands? I was trying with integration but the propagation of uncertainties is not trivial and extending the p.d.f. is not working.
Here is the example:
[code]#include “RooRealVar.h”
#include “RooDataSet.h”
#include “RooChebychev.h”
#include “TFile.h”
#include “TTree.h”
#include “TAxis.h”
#include “RooPlot.h”
#include “RooFitResult.h”
#include “RooExtendPdf.h”
#include
#include
using namespace RooFit;
int minimal(double sigL, double sigR, double minL, double MAXL, double minR, double MAXR)
{
// — Define variable and its ranges
//RooRealVar *mHNL = new RooRealVar(“Lambda0_M”,“m_{HNL}”,1500,6000,“MeV/c^{2}”);
RooRealVar *mHNL = new RooRealVar(“Lambda0_M”,“m_{HNL}”,minL,MAXR,“MeV/c^{2}”);
mHNL->setRange(“R”,minL,MAXR) ;
mHNL->setRange(“R1”,minL,MAXL) ;
mHNL->setRange(“R2”,minR,MAXR) ;
mHNL->setRange(“Sig”,sigL,sigR) ;
// — Define fitting function
RooRealVar *a0 = new RooRealVar(“a0”, “coefficient of x^1 term”,-1.,1.);
RooChebychev *p0 = new RooChebychev(“p0”,“p0”,*mHNL,RooArgSet(*a0));
// — Extend
RooRealVar *nbkg = new RooRealVar(“nbkg”,“number of bkg in whole region”,0,1e8);
RooExtendPdf *ep0 = new RooExtendPdf(“ep0”,“extended cheb PDF”,*p0,*nbkg,“R”);
// — Load file to dataset
RooDataSet* mHNL_dataset = p0->generate(*mHNL,2000) ;
// — Fit
RooFitResult* fitep0 = ep0->fitTo(*mHNL_dataset,RooFit::Range(“R1,R2”),RooFit::Save(),Extended(kTRUE));
cout<<" +++ Fit results +++ "<<endl;
fitep0->Print();
cout<<"Fit status= “<status()<<”; covQual= "<covQual()<<endl;
// — Plot
RooPlot* frame = mHNL->frame(RooFit::Bins(100)) ;
mHNL_dataset->plotOn(frame, RooFit::DrawOption(“Zp”)) ;
ep0->plotOn(frame);
frame->Draw();
return 0;
}[/code]