Re: different fit behaviour in ROOT and ROOFIT

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);
>>>> }
>>>>
>>>>
>>>
>>>
>
>