Re: again Ifit.C

From: Rene Brun (Rene.Brun@cern.ch)
Date: Wed Sep 17 1997 - 15:02:49 MEST


Wolfgang Korsch wrote:
> 
>    Hi again,
>      I'd like to use your Ifit.C function and I would like to extract the
>      final values for par[0] to par[3] (fitted values).
>      How do I do that? Can I use GetParameter ?
> 

I have modified the example Ifit.C that I posted a few months ago to illustrate
how to get the final values of parameters and the errors.
Note that you can also get access to the "assymetric errors" by calling mnerrs.
For this, look an example in TH1::Fit.

Rene Brun

//
//   Example of a program to fit non-equidistant data points
//   =======================================================
//
//   The fitting function fcn is a simple chisquare function
//   The data consists of 5 data points (arrays x,y,z) + the errors in errorsz
//   More details on the various functions or parameters for these functions
//   can be obtained in an interactive ROOT session with:
//    Root > TMinuit *minuit = new TMinuit(10);
//    Root > minuit->mnhelp("*",0)  to see the list of possible keywords
//    Root > minuit->mnhelp("SET",0) explains most parameters
//

Float_t z[5],x[5],y[5],errorz[5];

//______________________________________________________________________________
void fcn(Int_t &npar, Double_t *gin, Double_t &f, Double_t *par, Int_t iflag)
{
   const Int_t nbins = 5;
   Int_t i;

//calculate chisquare
   Double_t chisq = 0;
   Double_t delta;
   for (i=0;i<nbins; i++) {
     delta  = (z[i]-func(x[i],y[i],par))/errorz[i];
     chisq += delta*delta;
   }
   f = chisq;
}

//______________________________________________________________________________
Double_t func(float x,float y,Double_t *par)
{
 Double_t value=( (par[0]*par[0])/(x*x)-1)/ ( par[1]+par[2]*y-par[3]*y*y);
 return value;
}

//______________________________________________________________________________
void Ifit()
{
// The z values
	z[0]=1;
	z[1]=0.96;
	z[2]=0.89;
	z[3]=0.85;
	z[4]=0.78;
// The errors on z values
        Float_t error = 0.01;
	errorz[0]=error;
	errorz[1]=error;
	errorz[2]=error;
	errorz[3]=error;
	errorz[4]=error;
// the x values
	x[0]=1.5751;
	x[1]=1.5825;
	x[2]=1.6069;
	x[3]=1.6339;
	x[4]=1.6706;
// the y values
	y[0]=1.0642;
	y[1]=0.97685;
	y[2]=1.13168;
	y[3]=1.128654;
	y[4]=1.44016;

   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 parameters values and errors
   Double_t pvalue, perror, plow,phigh;
   char pname[3];
   for (Int_t i=0;i<4;i++) {
      sprintf(pname,"a%d",i+1);
      gMinuit->mnpout(i,pname,pvalue,perror,plow,phigh,ierflg);
      printf("Par %d, name=%s = %g +- %g\n",i,pname,pvalue,perror);
   }

// Print results using Minuit printing functions
   Double_t amin,edm,errdef;
   Int_t nvpar,nparx,icstat;
   gMinuit->mnstat(amin,edm,errdef,nvpar,nparx,icstat);
   gMinuit->mnprin(3,amin);
   
}



This archive was generated by hypermail 2b29 : Tue Jan 04 2000 - 00:26:20 MET