Fitting problem with TMinuit

ROOT
1

Dear Root Users,
I like to fit my data with user defined function, call the minimization using TMinuit. I like to estimate the error in the fitting parameter. From result it appears minimization routine is not working properly. Results are as follows

MATRIX FORCED POS-DEF BY ADDING -nan TO DIAGONAL
FCN=inf FROM MIGRAD
ERR MATRIX NOT POS-DEF
PARAMETER CORRELATION COEFFICIENTS all nan

Macro and .root file attached.
Another thing is how to extract chi2 and parameter from the minimization loop, so that one can plot chi2 vs parameter.
Thanks in advance,
Regards,

kaushik
error.txt (2.85 KB)
spectrum.root (3.88 KB)
minimise.C (4.45 KB)

2

Hi,

part of the function you use to perform your fit is not correctly written. For example, what did you mean to return here? (I indented the code)

Double_t func(Double_t *Ep,Double_t *par)
{
  Double_t arg=0;
  Double_t E=*Ep;
  if(par[1] !=0 && par[2] !=0){
    arg = (E - par[1])/par[2];
    Double_t fitval = par[0]*TMath::Exp(-0.5*((E-par[1])/par[2])*((E-par[1])/par[2])); 
   }
  return fitval;
}

It looks like the code could benefit from some reworking.

3

I think it is okay in this part, problem is in the minimization portion. Any way I welcome some more suggestion.

4

I think it is not ok.
The fitbal variable is defined in the scope after the if clause. This code does not even comoile as it is. I would start sorting out compilation errors, for example using aclic: root -b -q mymacro.C++.

5

Yes there are few errors when compile using .L minimise.C++, I solved them, now there are only few warnings for unused variable. corrected macro attached. But the main problem still persists.
minimise.C (4.47 KB)

6

Hi,

there was still a missing return statement (see indented and corrected code attached).
What one can see is that the calculation of chi2 has still issues. Indeed the first value of the “delta” variable in the fcn function is infinite. Either the calculation or the initial value of the parameters is to be sanitised.

//Content-Disposition: inline; filename="Vitt.C"

//
// 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(11);
// 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"
#include "TRandom.h"
#include "TMath.h"
#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>
#include <TGraphErrors.h>
#include <TFile.h>
#include "TH1F.h"

Double_t bincent[300], content[300], err[300];
extern Double_t func(Double_t *Ep, Double_t *par);
extern void fcn(Int_t &npar, Double_t *gin, Double_t &f, Double_t *par, Int_t iflag);
Int_t ncount = 0;

//______________________________________________________________________________
Int_t minimise()
{
   TFile *f1 = new TFile("spectrum.root");
   f1->ls();
   TH1F *hmass = (TH1F *)f1->Get("hx2");
   Int_t size = hmass->GetXaxis()->GetNbins();
   cout << size << endl;

   for (int i = 1; i <=  size; i++) {
      bincent[i] = hmass->GetBinCenter(i);
      content[i] = hmass->GetBinContent(i);
      err[i] = hmass->GetBinError(i);
   }



   const Int_t npars = 3;
   TMinuit *gMinuit = new TMinuit(npars); //initialize TMinuit with a maximum of 5 params
   gMinuit->SetFCN(fcn);

   Double_t arglist[300];
   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[3] = {220., 122., 28.6};
//static
   Double_t step[3] = {0.1, 0.1, 0.1};
   gMinuit->mnparm(0, "Constant", vstart[0], step[0], 0, 0, ierflg);
   gMinuit->mnparm(1, "Mean", vstart[1], step[1], 0, 0, ierflg);
   gMinuit->mnparm(2, "Sigma", vstart[2], step[2], 0, 0, ierflg);

// Now ready for minimization step
   arglist[0] = 500;
   arglist[1] = 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);

   Int_t iuext;
   TString chnam;   // The name of the parameter
   Double_t val;    // The current (external) value of the parameter
   Double_t errl;   // The current estimate of the parameter uncertainty
   Double_t xlolim; // The lower bound (or zero if no limits)
   Double_t xuplim; // The upper bound (or zero if no limits)
   Int_t iuint;     // The internal parameter number


   Double_t currentPar[npars] = {0};
   for (Int_t i = 0; i < npars; i++) {
      gMinuit->mnpout(i, chnam, currentPar[i], errl, xlolim, xuplim, iuint);
   }


   Int_t gcount = 300;
   TGraphErrors *gr = new TGraphErrors(gcount, bincent, content, 0, err);
   gr->Draw("AP");

   TF1 *fun_1 = new TF1("fun_1", func, 55, 195, 3);
   fun_1->SetParameters(currentPar);
   fun_1->SetLineColor(kBlue);
   fun_1->SetLineStyle(1);
   fun_1->SetLineWidth(2);
   fun_1->Draw("same");

   return 0;
}

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


//calculate chisquare
   Double_t chisq = 0;
   Double_t delta;
   for (Int_t i = 0; i < nbins; i++) {
      delta = (content[i] - func(&bincent[i], par)) / err[i];
      cout << "delta " << delta << endl;
      chisq += delta * delta;
   }
   f = chisq;
   ncount++;
}


Double_t func(Double_t *Ep, Double_t *par)
{
   Double_t arg = 0;
   Double_t E = *Ep;
   Double_t fitval;

   if (par[1] != 0 && par[2] != 0) {
      arg = (E - par[1]) / par[2];
      fitval =
         par[0] * TMath::Exp(-0.5 * ((E - par[1]) / par[2]) * ((E - par[1]) / par[2]));
   }

   return fitval;
}
7

Hi
There are few data point, where BinContent ane BinError are zero, These data points were creating the problem. I solve it by using following statement in chisquare calculation.
if(err[i]!=0){delta = (content[i]-func(&bincent[i],par))/err[i];}
Any suggestion for the second part of my question that is plotting chi2 vs parameter.

8

You have full control on your Chi2 as you calculate it in fcn. It is just a matter of saving (some of) its values along with the ones of the parameter(s) in some data structure and plot them.

9

chi2 I can extract from the function “fcn”, but how to extract the parameters for that particular chi2.

10

You have the double* variable called par in your code.

11

I know this, but simply
cout<<chi2<<par[i]<<endl; is not giving correct result.
Moreover I want to see all three parameters.
GetParameter(currentPar) is also not working.

12

You can access all the parameters with an integer index as that array contains the parameters (you do that in the func function - you use the parameters to calculate the value of the model you are using).

13

I got the parameters using
*par , *(par+1), *(par+2)
Thanks

14

Hi Kaushik,

Can you post the corrected minimise.c file.
I am trying to minimize chi squared of a function depends on 4 parameters and print the 4 parameters. So I thought I could get a help from your code.

Thanks.