RooFit convolution problem

[ripiccin]

Hello,

I’m trying to the define a fit function for my data that should be the convolution of two pdf, each defined as the sum of a smeared RooDecay pdf and a double gaussian. Here is the code

//variable definition
 RooRealVar t_val("t","t",20,25);
RooRealVar sigma_high0("#sigma_high","#sigma core",0.005,0.,0.02);
  RooRealVar sigma_high1("#sigma_high1","#sigma core",0.1,0.,0.2);
  RooRealVar sigma_res("#sigma_res","#sigma tail",0.005,0.,0.02);
  RooRealVar sigma_PD("#sigma_PD","#sigma photodiode",0.0005,0.,0.002);
  RooRealVar mean("#mu","#mu",22.37,22.3,22.5);
  RooRealVar tau("#tau","#alpha",0.03,0.,0.1);
  RooRealVar frac1("fraction1","fraction",0.8,0,1);
  RooRealVar frac0("fraction0","fraction",0.9,0,1);
  RooFormulaVar bias("bias","@0-@1",RooArgList(t_val,mean));

  //pdf def
  RooGaussian pdf_highgain0("pdf_highgain0","pdf_highgain0",t_val,mean,sigma_high0);
  RooGaussian pdf_highgain1("pdf_highgain1","pdf_highgain1",t_val,mean,sigma_high1);
  RooAddPdf pdf_highgain("pdf_highgain","pdf_highgain",RooArgList(pdf_highgain0,pdf_highgain1),RooArgList(frac0));
  RooGaussian res_PD("res_PD","photodiode resolution",t_val,mean,sigma_PD);
  RooGaussModel res("res","ph distribution",t_val,mean,sigma_res);
  RooDecay pdf_lowgain("pdf_lowgain","tail",t_val,tau,res,RooDecay::SingleSided);
  RooAddPdf calib_func("calib_func","calibration function",RooArgList(pdf_lowgain,pdf_highgain),RooArgList(frac1));
  RooFFTConvPdf calib_func1("calib_func1","convolution",t_val,calib_func,calib_func);
  RooDataHist data_calib("data_calib","dataset",t_val,hdata_calib);
  t_val.setBins(1000,"cache");

 sigma_high0.setVal(0.005);
 sigma_high1.setVal(0.005);
 sigma_res.setVal(0.005);
 sigma_PD.setVal(0.0005);
 mean.setVal(22.37);
 tau.setVal(0.03);
 frac1.setVal(0.8);
 frac0.setVal(0.9);
 mean.setConstant();

if I try to plot it, I get a flat constant line and a loop of errors like this

[#0] ERROR:Eval -- RooAbsReal::logEvalError(pdf_highgain_fft) evaluation error, 
 origin       : RooAddPdf::pdf_highgain_fft[ fraction0 * pdf_highgain0_fft + [%] * pdf_highgain1_fft ]
 message      : p.d.f value is Not-a-Number
 server values: !refCoefNorm=(t = 20.4713), !pdfs=(pdf_highgain0_fft = 0/0,pdf_highgain1_fft = 0/0), !coefficients=(fraction0 = 0.9)

what am I missing? or maybe have I reached a limit of what roofit con do?

thanks
Emanuele

[jonas]

Hi, welcome to the ROOT forum!

I don’t quite understand what you want to do in this line:

RooFFTConvPdf calib_func1("calib_func1","convolution",t_val,calib_func,calib_func);

You really want to compute the convolution of a PDF with itself?

I think you need to rethink mathematically what you want to do. But in any case, the problem is this:

RooRealVar t_val("t","t",20,25);

The definition range of your observable is from 20 to 25, and outside of it all your PDFs before the convolution will be zero. That’s how the ranges in RooFit work.

That means, mathematically, the convolution of your two PDFs will be non-zero in the range 40 to 50, which is outside the variable definition range. Inside the range, the convolution is always zero then, which results in the error and the flat line you see.

But I guess in reality, you want to do the convolution of some signal shape that is defined in the 20 to 25 range, and a resolution function, right? A resolution function is usually defined around zero, like a Gaussian smearing always has mean zero. To account for this, you will also have to extend the definition range of t_val, maybe even to negative values.

You can see something like that in one of our tutorials:
https://root.cern/doc/master/rf208__convolution_8C.html

I hope this helps!

Cheers,
Jonas

[ripiccin]

Hi Jonas,

many thanks for your reply. Yes actually that was an error I have corrected just after sending the post

this RooFFTConvPdf calib_func1(“calib_func1”,“convolution”,t_val,calib_func,calib_func);
is actually
RooFFTConvPdf calib_func1(“calib_func1”,“convolution”,t_val,calib_func,calib_func);

RooFFTConvPdf calib_func1(“calib_func1”,“convolution”,t_val,calib_func,pdf_PD);
where
RooDecay pdf_PD(“pdf_PD”,“tail”,t_val,tau_PD,gaus_PD,RooDecay::SingleSided);

The point is that I am not doing the convolution with a gaussian but with a RooDecay

matematically we have

calib_func1= calib_func X pdf_PD

I have set the range from -2,25 to include negative values, but still errors and flat 0 convoluted pdf

RooAbsMinimizerFcn: Minimized function has error status.
Returning maximum FCN so far (-inf) to force MIGRAD to back out of this region. Error log follows.
Parameter values: #sigma_PD=0.000146014 #sigma_high=0.005 #sigma_res=0.005 #tau=0.03 fraction1=0.8
RooNLLVar::nll_calib_func1_data_calib[ paramSet=(#mu,#sigma_PD,#sigma_high,#sigma_res,#tau,fraction1) ]
function value is NAN @ paramSet=(#mu = 22.37,#sigma_PD = 0.000146014 +/- 0.000728735,#sigma_high = 0.005 +/- 0.00420735,#sigma_res = 0.005 +/- 0.00728735,#tau = 0.03 +/- 0.038561,fraction1 = 0.8 +/- 0.336588)
RooFFTConvPdf::calib_func1[ calib_func(t) (*) pdf_PD(t) ]
getLogVal() top-level p.d.f evaluates to zero @ !x=t=22.1498, !xprime=NULL, !pdf1=calib_func=0, !pdf2=pdf_PD=0, !params=(#mu = 22.37,#sigma_high = 0.005 +/- 0.00420735,#sigma_res = 0.005 +/- 0.00728735,#tau = 0.03 +/- 0.038561,fraction1 = 0.8 +/- 0.336588,#sigma_PD = 0.000146014 +/- 0.000728735), !cacheObs=()
getLogVal() top-level p.d.f evaluates to zero @ !x=t=22.1503, !xprime=NULL, !pdf1=calib_func=0, !pdf2=pdf_PD=0, !params=(#mu = 22.37,#sigma_high = 0.005 +/- 0.00420735,#sigma_res = 0.005 +/- 0.00728735,#tau = 0.03 +/- 0.038561,fraction1 = 0.8 +/- 0.336588,#sigma_PD = 0.000146014 +/- 0.000728735), !cacheObs=()
getLogVal() top-level p.d.f evaluates to zero @ !x=t=22.1508, !xprime=NULL, !pdf1=calib_func=0, !pdf2=pdf_PD=0, !params=(#mu = 22.37,#sigma_high = 0.005 +/- 0.00420735,#sigma_res = 0.005 +/- 0.00728735,#tau = 0.03 +/- 0.038561,fraction1 = 0.8 +/- 0.336588,#sigma_PD = 0.000146014 +/- 0.000728735), !cacheObs=()
getLogVal() top-level p.d.f evaluates to zero @ !x=t=22.1513, !xprime=NULL, !pdf1=calib_func=0, !pdf2=pdf_PD=0, !params=(#mu = 22.37,#sigma_high = 0.005 +/- 0.00420735,#sigma_res = 0.005 +/- 0.00728735,#tau = 0.03 +/- 0.038561,fraction1 = 0.8 +/- 0.336588,#sigma_PD = 0.000146014 +/- 0.000728735), !cacheObs=()
getLogVal() top-level p.d.f evaluates to zero @ !x=t=22.1518, !xprime=NULL, !pdf1=calib_func=0, !pdf2=pdf_PD=0, !params=(#mu = 22.37,#sigma_high = 0.005 +/- 0.00420735,#sigma_res = 0.005 +/- 0.00728735,#tau = 0.03 +/- 0.038561,fraction1 = 0.8 +/- 0.336588,#sigma_PD = 0.000146014 +/- 0.000728735), !cacheObs=()
getLogVal() top-level p.d.f evaluates to zero @ !x=t=22.1523, !xprime=NULL, !pdf1=calib_func=0, !pdf2=pdf_PD=0, !params=(#mu = 22.37,#sigma_high = 0.005 +/- 0.00420735,#sigma_res = 0.005 +/- 0.00728735,#tau = 0.03 +/- 0.038561,fraction1 = 0.8 +/- 0.336588,#sigma_PD = 0.000146014 +/- 0.000728735), !cacheObs=()
getLogVal() top-level p.d.f evaluates to zero @ !x=t=22.1528, !xprime=NULL, !pdf1=calib_func=0, !pdf2=pdf_PD=0, !params=(#mu = 22.37,#sigma_high = 0.005 +/- 0.00420735,#sigma_res = 0.005 +/- 0.00728735,#tau = 0.03 +/- 0.038561,fraction1 = 0.8 +/- 0.336588,#sigma_PD = 0.000146014 +/- 0.000728735), !cacheObs=()
getLogVal() top-level p.d.f evaluates to zero @ !x=t=22.1533, !xprime=NULL, !pdf1=calib_func=0, !pdf2=pdf_PD=0, !params=(#mu = 22.37,#sigma_high = 0.005 +/- 0.00420735,#sigma_res = 0.005 +/- 0.00728735,#tau = 0.03 +/- 0.038561,fraction1 = 0.8 +/- 0.336588,#sigma_PD = 0.000146014 +/- 0.000728735), !cacheObs=()
getLogVal() top-level p.d.f evaluates to zero @ !x=t=22.1538, !xprime=NULL, !pdf1=calib_func=0, !pdf2=pdf_PD=0, !params=(#mu = 22.37,#sigma_high = 0.005 +/- 0.00420735,#sigma_res = 0.005 +/- 0.00728735,#tau = 0.03 +/- 0.038561,fraction1 = 0.8 +/- 0.336588,#sigma_PD = 0.000146014 +/- 0.000728735), !cacheObs=()
getLogVal() top-level p.d.f evaluates to zero @ !x=t=22.1543, !xprime=NULL, !pdf1=calib_func=0, !pdf2=pdf_PD=0, !params=(#mu = 22.37,#sigma_high = 0.005 +/- 0.00420735,#sigma_res = 0.005 +/- 0.00728735,#tau = 0.03 +/- 0.038561,fraction1 = 0.8 +/- 0.336588,#sigma_PD = 0.000146014 +/- 0.000728735), !cacheObs=()
getLogVal() top-level p.d.f evaluates to zero @ !x=t=22.1548, !xprime=NULL, !pdf1=calib_func=0, !pdf2=pdf_PD=0, !params=(#mu = 22.37,#sigma_high = 0.005 +/- 0.00420735,#sigma_res = 0.005 +/- 0.00728735,#tau = 0.03 +/- 0.038561,fraction1 = 0.8 +/- 0.336588,#sigma_PD = 0.000146014 +/- 0.000728735), !cacheObs=()
getLogVal() top-level p.d.f evaluates to zero @ !x=t=22.1553, !xprime=NULL, !pdf1=calib_func=0, !pdf2=pdf_PD=0, !params=(#mu = 22.37,#sigma_high = 0.005 +/- 0.00420735,#sigma_res = 0.005 +/- 0.00728735,#tau = 0.03 +/- 0.038561,fraction1 = 0.8 +/- 0.336588,#sigma_PD = 0.000146014 +/- 0.000728735), !cacheObs=()
… (remaining 990 messages suppressed)

bests

Emanuele