Hi,
I want to fit two histograms simultaneouslz with minuit. Since I don’t want to use RooFit but program a fit macro by myself I tried to use from the tutorial the “Ifit.C” example.
My idea how to do it is to fit both histograms in “fcn” – so I defined my histograms “*histo1” and “*histo2” outside of the functions as well as my function which I want to fit.
In the first step I tried to fit just one histogram which doesn’t worked…
ROOT always returns and division by zero for the “fcn” function.
Maybe somebody can gave my a hint how to use minuit correctly (especially how to choose the starting value and step sizes since I assume at the moment that it is the crux of the matter…)
Cheers & Thanks,
Thomas
Well my code looks like this:
#include "TMinuit.h"
TH1D* histo1 = (TH1D*) hHisto1->Clone(); //the histograms exist outside of my script
TH1D* histo2 = (TH1D*) hHisto2->Clone();
//My function I want to fit with one variable x and (at the moment) one parameter par[0]
Double_t HgFit(Double_t x, Double_t *par)
{ ...}
//function fcn for minuit
void fcn(Int_t &npar, Double_t *gin, Double_t &f, Double_t *par, Int_t iflag)
{
//variable definition
Int_t nbins1 = histo1->GetNbinsX();
Int_t nbins2 = histo2->GetNbinsX();
Double_t bincontent1=0., bincontent2=0., binerror1 = 0., binerror2 = 0.;
Double_t chisq1 = 0, chisq2 = 0;
Double_t delta1 = 0., delta2 = 0.;
for (Int_t i=1;i<nbins1>GetBinContent(i);
binerror1 = histo1->GetBinError(i);
delta1 = HgFit(bincontent1,par)/binerror1;
chisq1 += delta1*delta1;
}
f = chisq1;
}
//main function
void MyFit()
{
//both histograms should have the same binning ;)
if (!histo1->GetNbinsX()!=histo2->GetNbinsX())
{
cout<<"ERROR Histograms with unequal binning EXITING NOW";
return;
}
//from here I have not changed anything since I am not sure what it means in detail
TMinuit *gMinuit = new TMinuit(5); //initialize TMinuit with a maximum of 5 params
gMinuit->SetFCN(fcn);
Double_t arglist[10];
Int_t ierflg = 0;
arglist[0] = 1;
gMinuit->mnexcm("SET ERR", arglist ,1,ierflg);
// Set starting values and step sizes for parameters
static Double_t vstart[4] = {3, 1 , 0.1 , 0.01};
static Double_t step[4] = {0.1 , 0.1 , 0.01 , 0.001};
gMinuit->mnparm(0, "a1", vstart[0], step[0], 0,0,ierflg);
gMinuit->mnparm(1, "a2", vstart[1], step[1], 0,0,ierflg);
gMinuit->mnparm(2, "a3", vstart[2], step[2], 0,0,ierflg);
gMinuit->mnparm(3, "a4", vstart[3], step[3], 0,0,ierflg);
// Now ready for minimization step
arglist[0] = 500;
arglist[1] = 1.;
gMinuit->mnexcm("MIGRAD", arglist ,2,ierflg);
// Print results
Double_t amin,edm,errdef;
Int_t nvpar,nparx,icstat;
gMinuit->mnstat(amin,edm,errdef,nvpar,nparx,icstat);
gMinuit->mnprin(3,amin);
}