Fitting a function in the form of an integral with variable limits to a data set with Roofit

[Rahul]

Hello Roofit experts,
I wanted to fit a function in the form of an integral with variable limits to a data set. Fo this, I made a code which first does the integration. The integration works if Tmb and Tbe are constants . But I need to make them as a paramenter for the fit. Can you please tell me as to how I do that with the present form of my RooGenericPdf fmb and fbe ?

#ifndef __CINT__
#include "RooGlobalFunc.h"
#endif
#include "RooRealVar.h"
#include "RooGenericPdf.h"
#include "RooDataSet.h"
#include "RooProdPdf.h"
#include "RooGaussian.h"
#include "RooConstVar.h"
#include "RooExponential.h"
#include "TCanvas.h"
#include "TAxis.h"
#include "RooPlot.h"
#include "RooFitResult.h"
#include "RooDataHist.h"
#include "RooDataSet.h"
#include "TH1F.h"
#include "TFile.h"
#include "TH1F.h"
#include "TList.h"
#include "TLegend.h"
#include "TStyle.h"
#include "TLine.h"
#include "TCanvas.h"
#include <cmath>

using namespace RooFit ;

void testIntegration()
{

//setup  variables
    RooRealVar p("p","p",0.0,0.0,22.4);
    RooRealVar c("c","c",0.0,0.0,1.0);
    RooConstVar mass("mass","mass",0.139);
   
    RooConstVar Tmb("Tmb","Tmb",0.086);
    RooConstVar Tbe("Tbe","Tbe",0.091);

//RooRealVar Tmb("Tmb","Tmb",0.0,0.6);// Needs to to Give as Paramenter for the fit
//RooRealVar Tbe("Tbe","Tbe",0.0,0.6);// Needs to to Give as Paramenter for the fit


//Create pdf
    RooGenericPdf fmb("fmbA","fmbA","pow(p,2)*(exp(-sqrt((pow(p,2) + pow(mass,2)))/Tmb)) ",RooArgList(p,mass,Tmb));
    RooGenericPdf fbe("fbeA","fbeA","pow(p,2)/(exp(sqrt((pow(p,2) + pow(mass,2)))/Tbe) - 1)",RooArgList(p,mass,Tbe));


//Create Limit Function
    RooFormulaVar max("max","max","((-TMath::Log(2.0)*22.4*c) + sqrt( pow((TMath::Log(2.0)*22.4),2) - (pow(mass,2)*(1-pow(c,2))))) /(1-pow(c,2))",RooArgList(c,mass));
    RooRealVar    min("min","min",0.);

//Set Range  For integration
    p.setRange("range",min,max);

//Do the integration
    RooAbsReal* pdfMB = fmb.createIntegral(p,p,"range");
    RooAbsReal* pdfBE = fbe.createIntegral(p,p,"range");

//Create a frame
    RooPlot* frame = c.frame(Name("frame"),Title("Result of Integration Integration"));



    TFile* File1 = TFile::Open("/home/rahul/Documents/result.root");
    if(!File1) printf("File1 Not Found");
    TList* list1 = dynamic_cast<TList*> (File1->Get("HistList"));
    if(!list1) printf("List1 Not Found");
    TH1F* fHist1 = (TH1F*)list1->FindObject("fCG[%d]");
    if(!fHist1) printf("Hist Not Found");
    TH1F* fHist2 = (TH1F*)list1->FindObject("fCL[%d]");
    if(!fHist2) printf("Hist2Not Found");



    RooDataHist dhG("dhG","dhG",c, fHist1) ;
    RooDataHist dhL("dhL","dhL",c, fHist2) ;



    pdfBE->plotOn(frameMB,LineColor(kRed));
    pdfMB->plotOn(frameMB,LineColor(kBlack));

    // pdfMB->fitTo(dhG);//	chi2FitTo() ??

    TCanvas* MB = new TCanvas("CMB","CMB");
    frameMB->Draw();

}

I pasted the code above.
Regards
Rahul

[Axel]

Hi Rahul,

Unless someone else jumps in, @StephanH will likely be able to help when he’s back at work in about two weeks!

Cheers, Axel

[StephanH]

Hello Rahul,

I’m not sure if the ranges can be adjusted dynamically. Please try as follows:

• Set a default range for p. I.e., use p.setRange(min, max).
• When creating the integral, don’t specify a range. It should fall back to the default range now, and use a parametric binning.
• Try to fit and integrate.

Let me know if it works out.

[Rahul]

Hello Hageboeck,

Thanks for the reply!

I tried to do p.setRange(min, max) followed by not specifying the range in the createIntegral() method.

 RooAbsReal* pdfMB = fmb.createIntegral(p,p);
 RooAbsReal* pdfBE = fbe.createIntegral(p,p);

But the integration gives curve with strange and unexpected shape.
The log is as follows

[#1] INFO:NumericIntegration -- RooRealIntegral::init(fmb_Int[p]) using numeric integrator RooIntegrator1D to calculate Int(p)
[#1] INFO:NumericIntegration -- RooRealIntegral::init(fbe_Int[p]) using numeric integrator RooIntegrator1D to calculate Int(p)
[#0] ERROR:Plotting -- RooRealIntegral::fbe_Int[p]_Norm[p]:createPlotProjection: "c" is not a dependent and will be ignored.
[#1] INFO:NumericIntegration -- RooRealIntegral::init(fbe_Int[p]_Norm[p]) using numeric integrator RooIntegrator1D to calculate Int(p)
[#1] INFO:NumericIntegration -- RooRealIntegral::init(fbe_Int[p]) using numeric integrator RooIntegrator1D to calculate Int(p)
[#0] ERROR:Plotting -- RooRealIntegral::fmb_Int[p]_Norm[p]:createPlotProjection: "c" is not a dependent and will be ignored.
[#1] INFO:NumericIntegration -- RooRealIntegral::init(fmb_Int[p]_Norm[p]) using numeric integrator RooIntegrator1D to calculate Int(p)
[#1] INFO:NumericIntegration -- RooRealIntegral::init(fmb_Int[p]) using numeric integrator RooIntegrator1D to calculate Int(p)

Before fitting, I wanted to confirm that the integration is done correctly. Basically, what I want to do is three integrals.

1. Integral [range: 0 to Pt_Square(max) ] (E f(E) dPt^2 ) with Pt_Square(max) = function ( x axis variable of the histogram which I want to fit it with, which ranges from 0, 4)

PS: The code I made for this is pasted below

void Integrate()
{

//setup our variables
    RooRealVar E("E","E",0.0,0.0,22.4);
    RooRealVar z("z","z",0.01,2.0,4.0);// Decides the Integration limits && The X axis variable of the histogram to fit the integrated results with
    RooConstVar mass("mass","mass",0.139);
    
    RooConstVar Tmb("Tmb","Tmb",0.119);
    RooConstVar Tbe("Tbe","Tbe",0.134);

    RooRealVar ptSq; // .. 

//Create pdf
    RooGenericPdf fmb("fmb","fmb","E*exp(-E/Tmb) ",RooArgList(E,Tmb));
    RooGenericPdf fbe("fbe","fbe","E/(exp(E/Tbe) - 1)",RooArgList(E,Tbe));

//Create Limit Function
    RooFormulaVar max("max","max","pow((exp(-z)*22.4),2)-pow(mass,2)",RooArgList(z,mass));
    RooRealVar    min("min","min",0.);
   
    //Set Range  For integration
    ptSq.setRange("range",min,max);// Integration variable

//Do the integration
    RooAbsReal* pdfMB = fmb.createIntegral(ptSq,"range");
    RooAbsReal* pdfBE = fbe.createIntegral(ptSq,"range");

//Create a frame and plot
    RooPlot* frame = z.frame(Name("frame"),Title("Mass = 0.139"));
    pdfBE->plotOn(frame,LineColor(kRed));
    pdfMB->plotOn(frame,LineColor(kBlack));
   
    TCanvas* MB = new TCanvas("CMB","CMB");
    frame->Draw();
}

2. Integral[range: [range: 0 to Pz_(max) ] ( f(E) dPz^2) with Pz_(max) = function (Pt^2), which is the x axis of the distribution which i want to fit with

3. Integral [range: 0 to p(max) ] (f(E)p^2dp) with pmax = a function of ‘c’, which is the x axis of the distribution I want to fit it with( c ranges from 0,1).

I made the first code in the thread for the third case.

I got this idea of setting the range in the createIntegral() from this post: https://root-forum.cern.ch/t/intergration-with-variable-limits/11741

Also I was referring to this old post in root forum: https://root-forum.cern.ch/t/fitting-an-integral-function/8517/6

I am a bit confused with why the integration is not working at all now. If possible could you please suggest me some clear and correct way of doing it with root or roofit?

Regards
Rahul

[StephanH]

Hi Rahul,

in the code you pasted, RooFit doesn’t know what to integrate over. Specifically, fmb does not depend on ptSq. It’s like asking it to integrate f(x) * dy. That’s obviously constant. To demonstrate that it works, I decided to integrate over E instead. It works, and I am also able to change the integration range using z.

using namespace RooFit;

void Integrate()
{

//setup our variables
    RooRealVar E("E","E",0.0,0.0,22.4);
    RooRealVar z("z","z",0.01,2.0,4.0);// Decides the Integration limits && The X axis variable of the histogram to fit the integrated results with
    RooConstVar mass("mass","mass",0.139);
    
    RooConstVar Tmb("Tmb","Tmb",0.119);
    RooConstVar Tbe("Tbe","Tbe",0.134);


//Create pdf
    RooGenericPdf fmb("fmb","fmb","E*exp(-E/Tmb) ",RooArgList(E,Tmb));
    RooGenericPdf fbe("fbe","fbe","E/(exp(E/Tbe) - 1)",RooArgList(E,Tbe));

//Create Limit Function
    RooFormulaVar max("max","max","pow((exp(-z)*22.4),2)-pow(mass,2)",RooArgList(z,mass));
    RooRealVar    min("min","min",0.);
   
    //Set Range  For integration
    RooRealVar ptSq("ptSq", "ptSq", 0.); // .. 
    ptSq.setRange("range", min,max);// Integration variable
    E.setRange("range", min, max);

//Do the integration
    RooAbsReal* pdfMB = fmb.createIntegral(E, E, "range");
    //RooAbsReal* pdfBE = fbe.createIntegral(ptSq, ptSq, "range");

    pdfMB->Print("V");
    
//Create a frame and plot
    RooPlot* frame = z.frame(Name("frame"),Title("Mass = 0.139"));
    //pdfBE->plotOn(frame,LineColor(kRed));
    pdfMB->plotOn(frame,LineColor(kBlack));
   
    TCanvas* MB = new TCanvas("CMB","CMB");
    frame->Draw();
}

I hope that’s enough to set you on a path towards implementing all three steps.

[Rahul]

Dear Stephen,
Thanks a lot for pointing out this foolish mistake of not changing the variables. I changed the form of equations appropriately and it gives the desired shape (at least) after the integration. The modified code and result is pasted below.

#ifndef __CINT__
#include "RooGlobalFunc.h"
#endif
#include "RooRealVar.h"
#include "RooAbsPdf.h"
#include "RooGenericPdf.h"
#include "RooDataSet.h"
#include "RooProdPdf.h"
#include "RooGaussian.h"
#include "RooConstVar.h"
#include "RooExponential.h"
#include "TCanvas.h"
#include "TAxis.h"
#include "RooPlot.h"
#include "RooFitResult.h"
#include "RooDataHist.h"
#include "RooDataSet.h"
#include "TH1F.h"
#include "TFile.h"
#include "TH1F.h"
#include "TList.h"
#include "TLegend.h"
#include "TStyle.h"
#include "TLine.h"
#include "TCanvas.h"
#include <cmath>
#include "TGraph.h"

using namespace RooFit ;

void ZetaIntegration()
{

//setup our variables
    RooRealVar PtSq("PtSq","PtSq",0.01,0.01,100);
    RooRealVar z("z","z",0.0,0.0,4.0);

    RooConstVar mass("mass","mass",0.139);

    RooConstVar sQs("sQs","sQs",6.66);

    RooConstVar Tmb("Tmb","Tmb",.119);
    RooConstVar Tbe("Tbe","Tbe",.134);

 //Create Limit Function
    RooFormulaVar max("max","max","fabs(pow(sQs*exp(-z),2) - pow(mass,2))",RooArgList(sQs,z,mass));
    RooConstVar    min("min","min",0.0);
    PtSq.setRange("range",min,max);





//Create pdf
    RooGenericPdf fmb("fmb","fmb","sqrt(pow((PtSq - pow(sQs*exp(-z),2) + pow(mass,2)),2) + PtSq + pow(mass,2))*exp(-sqrt(pow((PtSq - pow(sQs*exp(-z),2) + pow(mass,2)),2) + PtSq + pow(mass,2))/Tmb)",RooArgList(PtSq,sQs,z,mass,Tmb));
    RooGenericPdf fbe("fbe","fbe","sqrt(pow((PtSq - pow(sQs*exp(-z),2) + pow(mass,2)),2) + PtSq + pow(mass,2))/(exp(sqrt(pow((PtSq - pow(sQs*exp(-z),2) + pow(mass,2)),2) + PtSq + pow(mass,2))/Tbe) - 1)",RooArgList(PtSq,sQs,z,mass,Tbe));


//Do the integration
    RooAbsReal* pdfMB = fmb.RooAbsReal::createIntegral(RooArgList(PtSq),RooArgList(PtSq,z),"range");
    RooAbsReal* pdfBE = fbe.RooAbsReal::createIntegral(RooArgList(PtSq),RooArgList(PtSq,z),"range");

//Create a frame
   RooPlot* frame = z.frame(Name("frame"),Title("Mass = 0.139"));
   pdfBE->plotOn(frame,LineColor(kBlack));
   pdfMB->plotOn(frame,LineColor(kRed));
   
   TCanvas* MB = new TCanvas("CMB","CMB");
   frame->Draw("SAME");

}

However, when I do the integration using
ROOT::Math::Functor1D f1(& funcMB);
ROOT::Math::Integrator igmb(f1, ROOT::Math::IntegrationOneDim::kADAPTIVE);
technique for one of the integration, I mentioned above, I am getting a different than from roofit. I am cross checking it and I shall come up with some more details. I do hope that you will take a look at it.
Sincerely
Rahul

[StephanH]

Hi Rahul,

one thing I can tell you already now is that it’s almost sure that numerical integrations are being used. Therefore, depending on the method (and settings), you will get different results.
Can you have a look at the terminal when the RooFit integration starts? Like this we would at least know which method it is using.

In the end, it might just be a matter of precision (and therefore runtime) to get a decent result. If we don’t manage to have all methods converge to the same point, it’s a systematic uncertainty of your analysis.

[Rahul]

Hello Stephen,
Thanks a lot for the suggestions and comments. When I set the precision
as fmb->Integral(GetMin(x),Getmax(x),1e-6); I got the same result as in the case of Roofit. So your suggestion hit the bull’s eye. I also tried to use the following lines

// Activate debug-level messages for topic integration to be able to follow actions below
  RooMsgService::instance().addStream(DEBUG,Topic(Integration)) ;

  // for normalization integrals for MINUIT to succeed.
  RooAbsReal::defaultIntegratorConfig()->setEpsAbs(1e-8) ;
  RooAbsReal::defaultIntegratorConfig()->setEpsRel(1e-8) ;

// Construct a custom configuration which uses the adaptive Gauss-Kronrod technique
  // for closed 1D integrals
  RooNumIntConfig customConfig(*RooAbsReal::defaultIntegratorConfig()) ;
  customConfig.method1D().setLabel("RooAdaptiveGaussKronrodIntegrator1D");

  // Another possibility: Change global default for 1D numeric integration strategy on finite domains
  RooAbsReal::defaultIntegratorConfig()->method1D().setLabel("RooAdaptiveGaussKronrodIntegrator1D");

and all of them gave me the expected results. So It was indeed a successful mission with those three integrations I mentioned. I am using my corresponding ROOT code only for fitting now as the fitted curves are not normalized to data points in the case of Roofit. I am still working on it. But it gave consistent parameter values as I compared them with the ROOT.

The log of integration is as follows.

[#2] INFO:Integration -- RooRealIntegral::ctor(fmb_Int[po|range]_Norm[co,po]) Constructing integral of function fmb over observables(po) with normalization (po,co) with range identifier range
[#2] DEBUG:Integration -- fmb : Observable po has parameterized binning, add value dependence of boundary objects rather than shape of leaf
[#2] INFO:Integration -- fmb: Observable po is suitable for analytical integration (if supported by p.d.f)
[#2] INFO:Integration -- fmb: Observables (po) are numerically integrated
[#2] INFO:Integration -- RooRealIntegral::ctor(fbe_Int[po|range]_Norm[co,po]) Constructing integral of function fbe over observables(po) with normalization (po,co) with range identifier range
[#2] DEBUG:Integration -- fbe : Observable po has parameterized binning, add value dependence of boundary objects rather than shape of leaf
[#2] INFO:Integration -- fbe: Observable po is suitable for analytical integration (if supported by p.d.f)
[#2] INFO:Integration -- fbe: Observables (po) are numerically integrated
[#2] INFO:Integration -- RooRealIntegral::ctor(fbe_Int[po|range]_Norm[co,po]_Int[]_Norm[co]) Constructing integral of function fbe_Int[po|range]_Norm[co,po] over observables() with normalization (co) with range identifier <none>
[#1] INFO:NumericIntegration -- RooRealIntegral::init(fbe_Int[po|range]_Norm[co,po]) using numeric integrator RooAdaptiveGaussKronrodIntegrator1D to calculate Int(po)
[#2] INFO:Integration -- RooRealIntegral::ctor(fmb_Int[po|range]_Norm[co,po]_Int[]_Norm[co]) Constructing integral of function fmb_Int[po|range]_Norm[co,po] over observables() with normalization (co) with range identifier <none>
[#1] INFO:NumericIntegration -- RooRealIntegral::init(fmb_Int[po|range]_Norm[co,po]) using numeric integrator RooAdaptiveGaussKronrodIntegrator1D to calculate Int(po)
Info in <TCanvas::MakeDefCanvas>:  created default TCanvas with name c1
root [1] 

In addition there were few mistakes in parenthesis which I missed while entering the expression to integrate. So the whole problem now has converged to normalizing the fitted curve to the data point in case of RooFit. But I successfully did it with ROOT macro.
Thanks again for the valuable suggestions. It indeed helped me to solve some of my long running issues.
Sincerely
Rahul

[StephanH]

Maybe the tutorial rf110 can help here. It explains the normalisation integrations for PDFs.
It’s in $ROOTSYS/tutorials/roofit/rf110_normintegration.[C|py] (online).