Re: Code Example

From: Rene Brun (Rene.Brun@cern.ch)
Date: Wed Mar 19 1997 - 12:37:44 MET


Matthias Vitt wrote:
> 
> Rene,
> Thank you very much for your help, it works very well! (too bad my
> measurements aren't better.. :-)  ) One thing: I looked up the
> documentation for mnexcm, but I couldn't find the options like "SET
> ERR","FIX" etc.  It would be very helpful if you add some comments about
> the flags you set in the code and the options you set. Maybe some words
> about how to set a parameter fixed during an optimization would be
> helpful , too. Is there a function in the TMinuit class that returns the
> values of the fit-parameters ? (I bet there is one, but I couldn't find
> it)

As promised yesterday (thanks Matthias), I post an example
of fitting procedure to be executed in batch. In the comments, 
I indicate how to get more help on the various Minuit functions
or parameters.
I remind you that the current implementation of the class TMinuit
is a translation of the Fortran Minuit package. We are now in the
process of redesigning this class with a better and simpler
user interface. This new version will not be part of 1.00, but
the next version.

Rene Brun


//
//   Example of a fitting program
//   ============================
//
//   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
//
#include "TROOT.h"
#include "TMinuit.h"

Int_t Error;
Float_t z[5],x[5],y[5],errorz[5];
extern Double_t func(float x,float y,Double_t *par);
extern void fcn(Int_t &npar, Double_t *gin, Double_t &f, Double_t *x, Int_t iflag);
Int_t ncount = 0;

//______________________________________________________________________________
main()
{
// 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;

   TROOT minexam("TestMinuit","Test of Minuit");

   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);
   
}
//______________________________________________________________________________
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;
   ncount++;
}

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;
}



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