Re: different fit behaviour in ROOT and ROOFIT

Dear Smbat,

thanks for the hint...indeed in the TestROOFIT.C version there was a convergence problem! However I modified the initial parameters, enlarging the limits and now I checked that the fit is indeed converging. In attachment you can find the updated TestROOFIT version.

However, even in this case, the number of JPsi I extract from ROOFIT is smaller than the one I obtain from ROOT. Now I have: ROOFIT: 2086+-159
ROOT: 2226 +- 189
So the difference should be in something else!

 Thanks again,

   Roberta   

Smbat Grigoryan wrote:
> Hi Roberta,
> It seems the difference is due to a bad fit convergence. How you have defined the "cbmean"
> limits 3 - 3.13 in TestROOFIT.C ? It gives CB integral: 2090 +- 158, CB width: 72 MeV,
> while 3 - 3.14 gives CB integral: 2076 +- 159, CB width: 71 MeV, and 3 - 3.15 gives:
> CB integral: 2710 +- 172, CB width: 110 MeV.
> Smbat
> ________________________________________
> From: owner-roottalk_at_root.cern.ch [owner-roottalk_at_root.cern.ch] on behalf of Roberta Arnaldi [arnaldi_at_to.infn.it]
> Sent: Monday, August 22, 2011 15:30
> To: Lorenzo Moneta
> Cc: <roottalk_at_lxroot01.cern.ch>
> Subject: Re: [ROOT] different fit behaviour in ROOT and ROOFIT
>
> Dear Lorenzo,
>
> I still have some differences performing fits to an histogram using ROOT
> or ROOFIT (see the attached mail thread). The fit is based on a
> CrystalBall function for the signal and a double exponential to describe
> the background.
> If you are back from holidays, can you please have a look to my code. I
> attach the .root file and the the two macros to perform the fit with
> ROOT or ROOFIT.
>
> Thanks in advance for your help!
>
> Roberta
>
> Roberta Arnaldi wrote:
>
>> Hi Lorenzo,
>>
>> thanks for your help and for looking into this problem when you'll be
>> back!
>> I checked, but the range of the function is the same in both cases.
>> The root file with the histo was already attached to my first mail,
>> but I'm resending it.
>>
>> Have nice holidays,
>>
>> Roberta
>>
>>
>>
>> Lorenzo Moneta wrote:
>>
>>> Hi Roberta,
>>> This was a possible difference looking at your code, maybe check also
>>> the range of the function in both cases. I can have a look at it, but
>>> I am right now in vacation and I could look at it only when I am back
>>> in the office in about 10 days. You should send me also the root file
>>> with the histograms data so I can run your code,
>>> Cheers, Lorenzo
>>> On Jul 27, 2011, at 5:52 PM, Roberta Arnaldi wrote:
>>>
>>>
>>>
>>>> Hi Lorenzo,
>>>>
>>>> thanks a lot for the very fast answer!
>>>> However, even if I remove from ROOT the option "I", I get even
>>>> higher values wrt ROOFIT.
>>>>
>>>> - ROOT without "I" : CB integral: 2226 +- 189, CB width: 78 MeV
>>>> - ROOT with "I": CB integral: 2220 +- 183, CB width: 76 MeV
>>>> - ROOFIT: CB integral: 2090 +- 158, CB width: 72 MeV
>>>>
>>>> So this does not seem to explain the difference...
>>>>
>>>> Thanks,
>>>>
>>>> Roberta
>>>>
>>>>
>>>> Lorenzo Moneta wrote:
>>>>
>>>>
>>>>> Hi Roberta, You are using in ROOT the option "I", integral of
>>>>> function in the bin. I don't think this is supported in RooFit and
>>>>> I guess this makes a difference in your case Best Regards
>>>>>
>>>>> Lorenzo
>>>>> On Jul 27, 2011, at 5:28 PM, Roberta Arnaldi wrote:
>>>>>
>>>>>
>>>>>
>>>>>
>>>>>> Dear All,
>>>>>>
>>>>>> I'm fitting an invariant mass spectrum histogram with a sum of
>>>>>> functions (Crystal Ball for the peak and two exponentials for the
>>>>>> background), using ROOT or ROOTFIT.
>>>>>> Initial parameters, fit limits...are the same in both cases.
>>>>>>
>>>>>> The value I want to extract is the number of signal events under
>>>>>> the Crystal Ball (CB) and the width of the Crystal Ball itself.
>>>>>> Even if the fit quality is good in both cases, the number of
>>>>>> events under the Crystal Ball is systematically different (~6-10%
>>>>>> depending on the fit), as well as the CB width. In particular,
>>>>>> ROOFIT gives always lower values. As an example:
>>>>>>
>>>>>> - ROOFIT: CB integral: 2090 +- 158, CB width: 72 MeV
>>>>>> - ROOT: CB integral: 2220 +- 183, CB width: 76 MeV
>>>>>>
>>>>>> Do you have any suggestion? I'm using root v5.28.00d and in ROOT I
>>>>>> use the options "LI".
>>>>>> I attach to the mail the root file and two macro I use to perform
>>>>>> the fit with root or roofit.
>>>>>>
>>>>>> Best regards and thanks in advance for your help,
>>>>>>
>>>>>> Roberta
>>>>>> #ifndef __CINT__
>>>>>> #include "RooGlobalFunc.h"
>>>>>> #endif
>>>>>> #include "RooRealVar.h"
>>>>>> #include "RooDataSet.h"
>>>>>> #include "RooGaussian.h"
>>>>>> #include "RooExponential.h"
>>>>>> #include "RooCBShape.h"
>>>>>> #include "RooConstVar.h"
>>>>>> #include "RooDataHist.h"
>>>>>> #include "RooAddPdf.h"
>>>>>> #include "RooWorkspace.h"
>>>>>> #include "RooFitResult.h"
>>>>>> #include "RooExtendPdf.h"
>>>>>> #include "RooPlot.h"
>>>>>> #include "TCanvas.h"
>>>>>> #include "TAxis.h"
>>>>>> #include "TStyle.h"
>>>>>> #include "TFile.h"
>>>>>> #include "TH1.h"
>>>>>> using namespace RooFit ;
>>>>>>
>>>>>> void TestROOFIT(){
>>>>>> gStyle->SetOptTitle(0) ;
>>>>>> gStyle->SetOptStat(0) ;
>>>>>> gStyle->SetPalette(1) ;
>>>>>> gStyle->SetCanvasColor(10) ;
>>>>>> gStyle->SetFrameFillColor(10) ;
>>>>>>
>>>>>> TFile *f = new TFile("Histo.root");
>>>>>> TH1D *histo = (TH1D*) f->Get("histo");
>>>>>> //------------------------------------
>>>>>> // define x range and frame
>>>>>> //------------------------------------
>>>>>>
>>>>>> RooRealVar x("DimuonMass","Mass",1.5,5.,"GeV/c^{2}") ;
>>>>>> RooPlot* frame = x.frame(Title("Test Roofit")) ;
>>>>>> TCanvas *c1 = new TCanvas("c1","c1",20,20,600,600);
>>>>>> gPad->SetLogy(1);
>>>>>>
>>>>>> //------------------------------------
>>>>>> // import histo //------------------------------------
>>>>>>
>>>>>> RooDataSet *datat; RooDataHist *datah;
>>>>>> datah = new RooDataHist("datah","datah",x,Import(*histo));
>>>>>> datah->plotOn(frame,Name("datah"),DataError(RooAbsData::SumW2),MarkerColor(kRed),MarkerSize(0.8),XErrorSize(0.));
>>>>>> datah->statOn(frame);
>>>>>>
>>>>>> //------------------------------------
>>>>>> // define Crystal Ball
>>>>>> //------------------------------------
>>>>>> RooRealVar cbmean("cb mean","cb_mean",3.096,3.,3.13);
>>>>>> RooRealVar cbsigma("cb sigma","cb_sigma",0.08,0.065,0.11);
>>>>>> RooRealVar cbn("cb n","cb_n",3.6,0.1,100.);
>>>>>> RooRealVar cbalpha("cb alpha","cb_alpha",1,0.1,10.);
>>>>>> RooCBShape cb("cb","cb",x,cbmean,cbsigma,cbalpha,cbn) ;
>>>>>>
>>>>>> cbalpha.setVal(0.98); cbalpha.setConstant(kTRUE);
>>>>>> cbn.setVal(5.2); cbn.setConstant(kTRUE);
>>>>>> //----------------------------------------------------
>>>>>> // fit expo1 //-----------------------------------------------------
>>>>>>
>>>>>> RooRealVar alpha1("bck1 alpha","bck1_alpha",-2.282635,-5.,0.) ;
>>>>>> RooRealVar nbck1("nbck1","nbck1",500,0.,1e10);
>>>>>> RooExponential exp1("exp1","exp1",x,alpha1) ;
>>>>>> x.setRange("rangeBck1",1.5,2.7);
>>>>>> RooExtendPdf ebck1_step1("ebck1_step1","ebck1_step1",exp1,nbck1);
>>>>>> RooAddPdf bck1_step1("bck1_step1","exp
>>>>>> bck1",RooArgList(ebck1_step1));
>>>>>>
>>>>>> //-----------------------------------------------------
>>>>>> // fit expo2 //
>>>>>> ----------------------------------------------------- RooRealVar
>>>>>> alpha2("bck2 alpha","bck2_alpha",-0.559789,-10.,0.) ;
>>>>>> RooRealVar nbck2("nbck2","nbck2",500,0.,1e10);
>>>>>> RooExponential exp2("exp2","exp2",x,alpha2) ;
>>>>>> RooExtendPdf ebck2_step1("ebck2_step1","ebck2_step1",exp2,nbck2);
>>>>>> RooAddPdf bck2_step1("bck2_step1","exp
>>>>>> bck2",RooArgList(ebck2_step1));
>>>>>> x.setRange("rangeBck2",3.5,5.);
>>>>>>
>>>>>> //------------------------------------ ----------------
>>>>>> // define expo1+expo2
>>>>>> //----------------------------------------------------- RooAddPdf
>>>>>> bck("bck","exp1+exp2",RooArgList(ebck1_step1,ebck2_step1));
>>>>>> //------------------------------------
>>>>>> // fit with signal+bck
>>>>>> //------------------------------------
>>>>>>
>>>>>> RooRealVar nsig("N.JPsi","nsignal",500,0.,1000000);
>>>>>> RooExtendPdf esig("esig","esig",cb,nsig);
>>>>>> RooAddPdf
>>>>>> sum("sum","bck1+bck2+cb",RooArgList(esig,ebck1_step1,ebck2_step1));
>>>>>> RooFitResult* r;
>>>>>> x.setRange("fitrange",2.,5.);
>>>>>> r = sum.fitTo(*datah,Range("fitrange"),Extended(kTRUE),Save()) ;
>>>>>>
>>>>>> //------------------------------------
>>>>>> // plot
>>>>>> //------------------------------------
>>>>>> sum.plotOn(frame,Components("cb"),LineColor(kGreen)) ;
>>>>>> sum.plotOn(frame,Components("exp1,exp2"),LineColor(kOrange)) ;
>>>>>> sum.plotOn(frame,Name("sum"),LineColor(kBlue)) ;
>>>>>> sum.paramOn(frame,Parameters(RooArgSet(cbmean,cbsigma,nsig)),Layout(0.6),Format("NELU",AutoPrecision(2)));
>>>>>>
>>>>>> frame->Draw() ;
>>>>>>
>>>>>> }
>>>>>> <Histo.root>void TestROOT(){
>>>>>>
>>>>>> gStyle->SetOptStat(0);
>>>>>> gStyle->SetOptTitle(0);
>>>>>> gStyle->SetCanvasColor(10);
>>>>>> gStyle->SetFrameFillColor(10);
>>>>>>
>>>>>> //----------------------------------------------------
>>>>>> // Open file
>>>>>> //----------------------------------------------------
>>>>>> TFile *f = new TFile("Histo.root");
>>>>>> TH1D *histo = (TH1D*) f->Get("histo"); TCanvas *c = new
>>>>>> TCanvas("c","c",20,20,600,600);
>>>>>> //----------------------------------------------------
>>>>>> // Fit the invariant mass spectrum // with a Crystal Ball function
>>>>>> + double exponential
>>>>>> //----------------------------------------------------
>>>>>> Double_t FitMin=2., FitMax=5.; TF1 *functot = new
>>>>>> TF1("functot",FuncTot,FitMin,FitMax,9);
>>>>>> functot->SetParameter(0,14.425808); // background parameters
>>>>>> functot->SetParameter(1,-2.282635); // background parameters
>>>>>> functot->SetParameter(2,6.466449); // background parameters
>>>>>> functot->SetParameter(3,-0.559789); // background parameters
>>>>>> functot->SetParameter(4,2859.295924); // JPsi normalization
>>>>>> functot->SetParameter(5,3.130); // JPsi mass position
>>>>>> functot->SetParameter(6,0.08); // JPsi width
>>>>>> functot->FixParameter(7,0.98); // Crystal Ball tails
>>>>>> functot->FixParameter(8,5.2); // Crystal Ball tails
>>>>>> functot->SetLineColor(kBlue);
>>>>>> functot->SetLineWidth(2);
>>>>>>
>>>>>> TFitResultPtr r = histo->Fit(functot,"RIlS");
>>>>>>
>>>>>> TMatrixDSym cov = r->GetCovarianceMatrix();
>>>>>> Double_t *fullmat;
>>>>>> fullmat = cov.GetMatrixArray();
>>>>>> Double_t psimat[25];
>>>>>> for(Int_t i=0;i<5;i++){
>>>>>> for(Int_t j=0;j<5;j++) psimat[5*i+j]=fullmat[40+j+9*i];
>>>>>> }
>>>>>>
>>>>>> histo->GetXaxis()->SetRangeUser(2.,5.);
>>>>>> histo->SetMinimum(1.);
>>>>>> gPad->SetLogy(1);
>>>>>> histo->Draw("e");
>>>>>> histo->SetMarkerStyle(20);
>>>>>> histo->SetMarkerColor(2);
>>>>>> histo->SetMarkerSize(0.7);
>>>>>>
>>>>>> TF1 *psifix = new TF1("psifix",FuncJpsi,0.,5.,9); for(int
>>>>>> i=0;i<9;i++) psifix->SetParameter(i,functot->GetParameter(i));
>>>>>> psifix->SetLineColor(kGreen);
>>>>>> psifix->Draw("same");
>>>>>> Double_t binWidth= histo->GetBinWidth(1);
>>>>>> Double_t NPsi=psifix->Integral(0.,5.)/binWidth;
>>>>>>
>>>>>> TF1 *psifix2 = new TF1("psifix2",FuncJpsi2,0.,5.,5); for(int
>>>>>> i=0;i<5;i++) psifix2->SetParameter(i,functot->GetParameter(i+4));
>>>>>> Double_t psipar[5];
>>>>>> for(int i=0;i<5;i++) psipar[i]=functot->GetParameter(4+i);
>>>>>> Double_t ErrPsiCorrParam =
>>>>>> psifix2->IntegralError(0.,5.,psipar,psimat)/binWidth;
>>>>>>
>>>>>> char text3[100];
>>>>>> sprintf(text3,"N_{J/#psi}= %3.0f #pm %2.0f",NPsi,ErrPsiCorrParam);
>>>>>> TLatex *l1 = new
>>>>>> TLatex(3.35,psifix2->Integral(0.,5.)/binWidth*8/5.,text3);
>>>>>> l1->Draw();
>>>>>>
>>>>>> char text5[100];
>>>>>> sprintf(text5,"#sigma_{J/#psi}= %5.3f #pm %5.3f
>>>>>> GeV/c^{2}",functot->GetParameter(6),functot->GetParError(6));
>>>>>> TLatex *l3 = new
>>>>>> TLatex(3.35,psifix2->Integral(0.,5.)/binWidth*0.3*5./1.5,text5);
>>>>>> l3->SetTextSize(0.035);
>>>>>> l3->Draw();
>>>>>> }
>>>>>>
>>>>>> Double_t FuncJpsi(Double_t *x, Double_t *par){ //Crystal Ball
>>>>>> Double_t FitJPsi;
>>>>>> Double_t t = (x[0]-par[5])/par[6];
>>>>>> if (t > (-par[7])){
>>>>>> FitJPsi = par[4]*TMath::Exp(-t*t/2.);
>>>>>> } else if (t <= (-par[7])) {
>>>>>> Double_t AA =
>>>>>> TMath::Power(par[8]/TMath::Abs(par[7]),par[8])*TMath::Exp(-TMath::Abs(par[7])*TMath::Abs(par[7])/2.);
>>>>>>
>>>>>> Double_t BB = par[8]/TMath::Abs(par[7])-TMath::Abs(par[7]);
>>>>>> if(TMath::Power((BB-t),par[8])!=0){
>>>>>> FitJPsi = par[4]*AA/TMath::Power((BB-t),par[8]);
>>>>>> } else FitJPsi = 0;
>>>>>> } return FitJPsi;
>>>>>> }
>>>>>>
>>>>>> Double_t FuncJpsi2(Double_t *x, Double_t *par){ //Crystal Ball
>>>>>> Double_t FitJPsi;
>>>>>> Double_t t = (x[0]-par[1])/par[2];
>>>>>> if (t > (-par[3])){
>>>>>> FitJPsi = par[0]*TMath::Exp(-t*t/2.);
>>>>>> } else if (t <= (-par[3])) {
>>>>>> Double_t AA =
>>>>>> TMath::Power(par[4]/TMath::Abs(par[3]),par[4])*TMath::Exp(-TMath::Abs(par[3])*TMath::Abs(par[3])/2.);
>>>>>>
>>>>>> Double_t BB = par[4]/TMath::Abs(par[3])-TMath::Abs(par[3]);
>>>>>> if(TMath::Power((BB-t),par[4])!=0){
>>>>>> FitJPsi = par[0]*AA/TMath::Power((BB-t),par[4]);
>>>>>> } else FitJPsi = 0;
>>>>>> } return FitJPsi;
>>>>>> }
>>>>>>
>>>>>> Double_t FuncBck1(Double_t *x, Double_t *par){ //exponential
>>>>>> Double_t FitBck1 = exp(par[0]+par[1]*x[0]);
>>>>>> return FitBck1;
>>>>>> }
>>>>>> Double_t FuncBck2(Double_t *x, Double_t *par){ //exponential
>>>>>> Double_t FitBck2 = exp(par[2]+par[3]*x[0]);
>>>>>> return FitBck2;
>>>>>> }
>>>>>>
>>>>>> Double_t FuncBck(Double_t *x, Double_t *par){ //exponential
>>>>>> return FuncBck1(x,par)+FuncBck2(x,par);
>>>>>> }
>>>>>>
>>>>>> Double_t FuncTot(Double_t *x, Double_t *par){ // total fit
>>>>>> return FuncBck(x,par)+FuncJpsi(x,par);
>>>>>> }
>>>>>>
>>>>>>
>>>>>
>>>
>
>