Re: ROOFIT and TMinuit

Hello Anil,

  Although I don't have the data file you are fitting I presume the difference is due to the treatment of the empty bins. In bare ROOT and also in your case you are skipping the empty bins, while in Roofit they are considered with an error given by the 68% CL of the Poisson distribution with mu=0. The fit will be biased in both cases, since when the number of entry/bins is less than ~5, the normal approximation for the
bin error is not anymore valid.
If you have low statistics the reccomeded way to fit an histogram is by using a binned likelihood method using Poisson statistics for the bin content.
This can be done in ROOT by fitting an histogram using the TH1::Fit method with option "L".

In you case also you don;t need to use the TMinuit class, but when you can just use TH1::Fit also for a least square fit.

  Best Regards

  Lorenzo
On Aug 19, 2009, at 7:36 AM, Anil Singh wrote:

> Dear Rooters,
>
> I have certain confusions (may be very stupid )regarding the working
> of TMinuit
> and RooFit packages. I am vaguley aware that RooFit and TMinuit use
> Minuit2
> engine for minimization.
>
> I am trying to fit a resonance using Breit Wigner shape. The fit
> looks very nice
> with the code that I created using example given in minuit manual.
> But the
> chi square shows very high (~21) even if I limit myself to the range
> where fit is especially
> good. Suspecting some hidden errors I moved to use TMinuit class but
> the results
> were same. Now just when I was ready to accept the numbers I made
> last
> validation effort using the RooFit package. The chisquare of the fit
> now turned out
> to be "~1.8", while the visual features of plot and the values of
> fiited parameters
> looked same as is other two attempts. I dont know what to believe
> in. Here are
> the TMinuit and RooFit Program that I used.
>
>
> TMINUIT
> ====================================================================
>
> #include <iostream>
> #include <fstream>
> #include <cstdlib>
> #include <cmath>
> #include <string>
> #include <vector>
>
> #include <TMinuit.h>
> #include <TCanvas.h>
> #include <TStyle.h>
> #include <TROOT.h>
> #include <TF1.h>
>
>
>
> TFile* f1 = new TFile("MyAnalysis.root");
> TH1D* myHist =(TH1D*) gDirectory->Get("h_z_mass1");
> std::vector<double>* the_pos = new std::vector<double>() ;
> std::vector<double>* the_mes = new std::vector<double>();
> std::vector<double>* the_err = new std::vector<double>();
>
> double yield(){
> double sum=0;
> for(int i=0; i!=(*the_mes).size(); i++){
> sum = sum+(*the_mes)[i];
> }
> return sum;
> }
>
>
> double BreitWigner(double* xptr, double* par){
> double mean = par[0];
> double gamma=par[1];
> double cons=par[2];
> double m2 = mean*mean;
> double g2 = gamma*gamma;
> double g2overm2 = g2/m2;
> double x = *xptr;
> double s = x*x;
> double deltaS = s-m2;
> double lineshape=0;
> if (fabs(deltaS/m2)<16){
> double prop = deltaS*deltaS + s*s*g2overm2;
> lineshape = (2/3.14)*g2*m2/(deltaS*deltaS + s*s*(g2overm2));
> }
> return cons*lineshape;
> }
>
> void ReadData(){
> myHist->Sumw2();
> int n = myHist->GetNbinsX();
> for(int i=0; i<n; i++){
>
> double content = myHist->GetBinContent(i+1);
> double x = myHist->GetBinCenter(i+1);
> double er = myHist->GetBinError(i+1);
> if((x<40)||(x>200))
> continue;
> the_pos->push_back(x);
> the_mes->push_back(content);
> the_err->push_back(er);
> }
>
> }
>
>
>
> void fcn(Int_t &npar, Double_t *gin, Double_t &f, Double_t *par,
> Int_t iflag)
> {
> Int_t nbins = the_mes->size();
> // cout<<"-------------->>"<<nbins<<endl;
> //calculate chisquare
> Double_t chisq = 0;
> Double_t delta;
> for (int i=0;i<nbins; i++) {
> if((*the_err)[i]){
> double thepoint = (*the_pos)[i];
> delta = ((*the_mes)[i]-BreitWigner(&thepoint,par))/(*the_err)
> [i];
> chisq += delta*delta;}
> // cout<<"delta: "<<delta<<endl;
> }
> cout<<"chisq: "<<chisq/nbins<<endl;
> f = chisq;
>
> }
>
>
> void main(int argc, char **argv){
>
> TFile* file=new TFile("TMinuitFit.root","recreate");
> TF1* func = new TF1("funcplot", BreitWigner,40,195,3);
>
>
>
> ReadData();
> TMinuit *gMinuit = new TMinuit(3);
> gMinuit->SetFCN(fcn);
> Double_t arglist[3];
> Int_t ierflg = 0;
>
> double yld = yield();
> arglist[0] = 1;
>
>
>
> gMinuit->mnexcm("SET ERR", arglist ,1,ierflg);
> Double_t vstart[3] = {91,2.5,12000};
> Double_t step[3] = {10 , 0.01 , 1200};
>
> gMinuit->mnparm(0, "mean", vstart[0], step[0], 0,0,ierflg);
> gMinuit->mnparm(1, "gamma", vstart[1], step[1], 0,0,ierflg);
> gMinuit->mnparm(2, "cons", vstart[2], step[2], 0,0,ierflg);
>
> arglist[0]=5000;
> arglist[1]=0.1;
> gMinuit->mnexcm("SET STR", arglist,2,ierflg);
> gMinuit->mnexcm("MIGRAD", arglist, 2, ierflg);
> double err[3];
> double outpar[3];
> for(int i=0; i<3; i++){
> gMinuit->GetParameter(i,outpar[i],err[i]);
> }
> func->SetParameters(outpar);
> func->SetLineWidth(2);
> file->cd();
> myHist->Write();
> func->Write();
> file->Write();
> file->Close();
> cout<<"=================================="<<endl;
>
> Double_t amin,edm,errdef;
> Int_t nvpar,nparx,icstat;
> gMinuit->mnstat(amin,edm,errdef,nvpar,nparx,icstat);
> gMinuit->mnprin(3,amin);
> cout<<"AMIN: "<<amin<<std::endl;
>
> //
> =
> =
> =
> =
> =
> ======================================================================
>
>
> ROOFIT
>
>
>
>
> //
> =
> ======================================================================
>
>
>
> //Access Data from histograms
>
> TFile* f1 = new TFile("MyAnalysis.root");
> TH1D* myHist =(TH1D*) gDirectory->Get("h_z_mass1");
> myHist->Sumw2();
> //Declare variable for quantity under study.
> RooRealVar mass("mass","mass",40,200);
>
> //Declare a data set using histogram.
> RooDataHist data("mass","mass",mass,myHist);
>
> //Make a Briet-Wigner Model.
> RooRealVar mean("mean","mean",70,110);
> RooRealVar width("gamma","gamma",0,10);
> RooBreitWigner bw("bw","bw",mass,mean,width);
>
> RooRealVar n1("n1","peak",0,20000);
> RooAddPdf model("model","model",RooArgList(bw),n1);
>
> //Fit the Model to data
> // model.fitTo(data, RooFit::Range(40,200));
>
> RooChi2Var Chi2("Chi2","Chi2",model,data,true);
> RooMinuit m2(Chi2);
> m2.migrad();
> // m2.hesse();
> RooFitResult* r2 = m2.save() ;
>
> //Print the Parameter Info
> mean.Print();
> width.Print();
> n1.Print();
>
> //Plotting Bussiness
> RooPlot* xframe = mass.frame();
> data.plotOn(xframe, RooFit::Name("data"));
> model.plotOn(xframe,RooFit::Name("model"));
> xframe->Draw();
>
> //Goodness of Fit
> double chi2 = xframe->chiSquare("model","data",3);
> double ndoff = xframe->GetNbinsX();
> cout<<"chi2/ndoff: "<<chi2<<"/"<<ndoff<<'\t'<<
> chi2/ndoff<<endl;
>
> //
> =
> =
> =
> =
> ======================================================================
>
>
> Any advice is greatly appreciated.
> Anil Singh
>
>
>
>
>
>
Received on Wed Aug 19 2009 - 11:30:08 CEST