ROOFIT and TMinuit

From: Anil Singh <anil79_at_fnal.gov>
Date: Wed, 19 Aug 2009 10:36:51 +0500

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 - 07:36:58 CEST