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.
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);
>> }
>>
>
>
Received on Wed Jul 27 2011 - 17:48:26 CEST
This archive was generated by hypermail 2.2.0 : Wed Jul 27 2011 - 23:50:01 CEST