Decay time fit with errors that are a function of decay time

Dear experts,
I am trying to do an extension of the decay time fit example in Roofit rf307 example, where I want to build a PDF that has the bias and sigma,each as a function of decay time.
(also want to include an acceptance function, but first just getting this part to work without acceptance.) The basic example works if bias and sigma are constants, but when I make them RooForumlaVars with a simple linear function, then the fit clearly doesn’t describe the toy data.

The code is a simple adaptation as shown below.

Any suggestions on what I am doing wrong here, or is there a different way to go about this kind of model?

Thanks so much!

#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 ;
void BDECAY3()
   // 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
   RooRealVar dt("dt","dt",0,10) ;
   RooRealVar dterr("dterr","per-event error on dt",0.01,10) ;
   // Build a gaussian resolution model scaled by the per-event error = gauss(dt,bias,sigma*dterr)
   RooRealVar bias("bias","bias",0);//,-10,10) ;
   RooConstVar sigma1_0("sigma1_0","sigma1_0",0.1);
   RooConstVar sigma1_1("sigma1_1","sigma1_1",0.1);
   RooFormulaVar sigma1("sigma1","sigma1_0+sigma1_1*dt",RooArgList(dt,sigma1_0, sigma1_1));
   RooRealVar sigma("sigma","per-event error scale factor",1.0,0.1,10) ;
   RooGaussModel gm("gm1","gauss model scaled bt per-event error",dt,bias,dterr,sigma1) ;
   // Construct decay(dt) (x) gauss1(dt|dterr)
   RooRealVar tau("tau","tau",1.60,1.40,1.80);
   RooDecay decay_gm("decay_gm","decay",dt,tau,gm,RooDecay::SingleSided) ;
   // 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 = pdfDtErr.generate(dterr,10000) ;
   // 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(decay_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 = model.generate(RooArgSet(dt,dterr),50000) ;
   // 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,ConditionalObservables(dterr)) ;
   // 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() ;


Hi @sblusk1,

First of all, welcome to the ROOT forum. Maybe @jonas or @moneta can help with that.


Thanks @jalopezg, I will wait to see if @moneta or @jonas can point me in the right direction. Key point here is that the average decay time error is a function of the decay time.