On Aug 29, 2011, at 2:50 PM, Roberta Arnaldi wrote: > Dear All, > > sorry to bother again...but I still have problems in the comparison of root and roofit fits and I have really no idea on how to understand this difference. Fit ranges and procedure are the same. Furthermore I always find this discrepancy while fitting these spectra... > > Thanks in advance for your help! > > Roberta > > ---------- Forwarded message ---------- > Date: Mon, 22 Aug 2011 15:30:55 +0200 > From: Roberta Arnaldi > To: Lorenzo Moneta > Cc: "" > 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() ; >>>>>> } >>>>>> 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); >>>>>> } >>>>>> >>> >