Gaussian pdf with per event error

Roofit and RooStats
1

Hi,
I have an observable x and its per-event error estimate x_err. The pdf for x is assumed to be a gaussian, its width being a rescale of x_err by a factor S. In other words, I need to define
pdf(x) = gauss(x|mean, S * x_err)
so I try via

[code]// observables
x = new RooRealVar(“x”, “x”, m_minX, m_maxX);
x_err = new RooRealVar(“x_err”, “x_err”, m_minXerr, m_maxXerr);

// ranges
x->setRange(“integration_range”, m_minX, m_maxX);

// parameters
scale_factor = new RooRealVar(“scale_factor”, “per-event error scale factor”, 4, 0.1, 10);
mean = new RooRealVar(“mean”, “mean”, m_initialGuess, m_minMean, m_maxMean);
sigma = new RooProduct(“sigma”,“sigma”, RooArgList(*scale_factor, *x_err));

// signal pdf
gauss = new RooGaussian(“gauss”, “gauss(invMass,mean,sigma)”, *invMass, *mean, *sigma);

// background parameter and pdf
slope = new RooRealVar(“slope”, “slope”, -1, -4., +4.0);
expo = new RooExponential(“expo”, “expo”, *x, *slope);

// sig+bkg mixture parameter
f_gaus_sig = new RooRealVar(“f_gaus_sig”, “gauss signal fraction”, 0.09, 0., 1.);

// composite model and fit
total_pdf = new RooAddPdf(“pdf”, “pdf”, RooArgList(*gauss, *expo), RooArgList(*f_gauss_sig));
total_pdf->fitTo(dataset, RooFit::ConditionalObservables(*x_err), RooFit::Range(“integration_range”));[/code]
but in this way scale_factor seems not to be treated as a fit parameter (and actually the fit returns a lot of errors)…

What is my mistake?

2

Hi,

I see no mistake in the code you posted here. If you provide me with a complete example,
I might be able to help you more.

Wouter

3

Hi,

how did you solve this problem?
I tried to rewrite roofit tutorial 307 https://root.cern.ch/root/html/tutorials/roofit/rf307_fullpereventerrors.C.html without decay gm, but I failed.

#ifndef CINT
#include “RooGlobalFunc.h”
#endif
#include “RooRealVar.h”
#include “RooDataSet.h”
#include “RooGaussian.h”
#include “RooGaussModel.h”
#include “RooConstVar.h”
#include “RooDecay.h”
#include “RooLandau.h”
#include “RooProdPdf.h”
#include “RooHistPdf.h”
#include “RooPlot.h”
#include “TCanvas.h”
#include “TAxis.h”
#include “TH1.h”
using namespace RooFit ;
#include “TFile.h”
#include “TTree.h”

void bfit()
{
// B - p h y s i c s p d f w i t h p e r - e v e n t G a u s s i a n r e s o l u t i o n
// ----------------------------------------------------------------------------------------------

// Observables
TTree *ts;
TFile *f1 = new TFile(TString(in));
ts = (TTree*)f1->Get("Bd");
TTree * t2 = ts-> CloneTree(10000)  ;   

RooRealVar dt("Bd_VTX_mass","Bd_VTX_mass",5000,5700) ;
RooRealVar dterr("Bd_VTX_massError","Bd_VTX_massError",0,360);
RooDataSet* data= new RooDataSet("alldata","alldata",RooArgSet(dt,dterr),Import(*t2)) ;   


// Build a gaussian resolution model scaled by the per-event error = gauss(dt,bias,sigma*dterr)
RooRealVar bias("bias","bias",0,5000,5600) ;
RooRealVar sigma("sigma","per-event error scale factor",1,0.1,10) ;
RooGaussModel gm("gm1","gauss model scaled bt per-event error",dt,bias,sigma,dterr) ;


// C o n s t r u c t   e m p i r i c a l   p d f   f o r   p e r - e v e n t   e r r o r
// -----------------------------------------------------------------

// Use landau p.d.f to get empirical distribution with long tail
//RooLandau pdfDtErr("pdfDtErr","pdfDtErr",dterr,RooConst(1),RooConst(0.25)) ;
RooDataSet* expDataDterr = new RooDataSet("erdata","erdata",dterr,Import(*t2)) ;

// Construct a histogram pdf to describe the shape of the dtErr distribution
RooDataHist* expHistDterr = expDataDterr->binnedClone() ;
RooHistPdf pdfErr("pdfErr","pdfErr",dterr,*expHistDterr) ;


// C o n s t r u c t   c o n d i t i o n a l   p r o d u c t   d e c a y _ d m ( d t | d t e r r ) * p d f ( d t e r r )
// ----------------------------------------------------------------------------------------------------------------------

// Construct production of conditional decay_dm(dt|dterr) with empirical pdfErr(dterr)
RooProdPdf model("model","model",pdfErr,Conditional(gm,dt)) ;

// (Alternatively you could also use the landau shape pdfDtErr)
//RooProdPdf model("model","model",pdfDtErr,Conditional(decay_gm,dt)) ;



// S a m p l e,   f i t   a n d   p l o t   p r o d u c t   m o d e l 
// ------------------------------------------------------------------

// Specify external dataset with dterr values to use model_dm as conditional p.d.f.
//RooDataSet* data = gm.generate(RooArgSet(dt,dterr),10000) ;



// F i t   c o n d i t i o n a l   d e c a y _ d m ( d t | d t e r r )
// ---------------------------------------------------------------------

// Specify dterr as conditional observable
 model.fitTo(*data) ;



// P l o t   c o n d i t i o n a l   d e c a y _ d m ( d t | d t e r r )
// ---------------------------------------------------------------------


// Make two-dimensional plot of conditional p.d.f in (dt,dterr)
TH1* hh_model = model.createHistogram("hh_model",dt,Binning(50),YVar(dterr,Binning(50))) ;
hh_model->SetLineColor(kBlue) ;


// Make projection of data an dt
RooPlot* frame = dt.frame(Title("Projection of model(dt|dterr) on dt")) ;
data->plotOn(frame) ;
model.plotOn(frame) ;


// Draw all frames on canvas
TCanvas* c = new TCanvas("rf307_fullpereventerrors","rf307_fullperventerrors",800, 400);
c->Divide(2) ;
c->cd(1) ; gPad->SetLeftMargin(0.20) ; hh_model->GetZaxis()->SetTitleOffset(2.5) ; hh_model->Draw("surf") ;
c->cd(2) ; gPad->SetLeftMargin(0.15) ; frame->GetYaxis()->SetTitleOffset(1.6) ; frame->Draw() ;

}