Fitting to calculate chi squared

ROOT
1

I’m using a multi-parameter fit with TMinuit to calibrate a detector array. It was working well previously, but I wanted to add a condition for the location of centroids that must be fitted with a user-defined function inside the minimization function. This seems to interfere with Minuit.

[quote]MINUIT WARNING IN MIGRAD
============== Negative diagonal element 1 in Error Matrix
MINUIT WARNING IN MIGRAD
============== Negative diagonal element 2 in Error Matrix
MINUIT WARNING IN MIGRAD
============== Negative diagonal element 3 in Error Matrix
…
[/quote]
Then a lot of ‘nan’ values end up in the minimization, which then fails.

This happens just from running the fit, even without considering the fit results in the chi squared value.

[code]	detCent[nuc][line][ii]=SeGAfits->FitLine(sega_e_calDop[ii][nuc],lineE[nuc][line]);

// delta = (par[23+thisLine]-detCent[nuc][line][ii])/SeGAfits->GetParErr(2);
// chi2nuc[nuc] += 8deltadelta;
// chisq += 8deltadelta;
// cout << " " << delta*delta ;
[/code]

Is there some way to get around this problem?

Any suggestions or help will be much appreciated,
Matt

2

Hi Matt,

You must have a problem in your function . I see already in your snippet many indices that can go out of bounds .

You will have to show your code . Response time and likelihood
of an answer will go down exponentially with the size of the
code :slight_smile:

Eddy

3

The code is large and in several pieces, so I’ll try to spare you the details. The indices are all safe. If I replace the fitting with some do-nothing code that uses all of the same arrays, as below, it runs fine.

detCent[nuc][line][ii] = lineE[nuc][line]; sega_e_calDop[ii][nuc]->GetBinContent(200); //detCent[nuc][line][ii]=SeGAfits->FitLine(sega_e_calDop[ii][nuc],lineE[nuc][line]);
so it’s not an issue of the arrays being out of bounds.

And the fitting class is well tested. I’m fairly certain that the trouble comes from using my simple gaussian fitting class, which contains

TF1 *gausMod = new TF1("gausMod","[0]+[1]*exp(-(x-[2])^2/(2*[3]^2))",0,8000);
...
  gamSpectrum->Fit("gausMod","MQ+","",peakBinE-FIT_HW,peakBinE+FIT_HW);

inside the function

void fcn(Int_t &npar, Double_t *gin, Double_t &f, Double_t *par, Int_t iflag) { //calculate chisquare ... f = chisq; }
which defines the chi squared value for my minimization using TMinuit

TMinuit *gMinuit = new TMinuit(28); gMinuit->SetFCN(fcn);
So I suspect that the problem comes from fitting a TF1 inside the function fcn that is called by the TMinuit object. It seems to be corrupting the variables and parameters of the larger minimization problem. I wouldn’t expect this to happen since the TF1 function is fit inside another object, but I guess it has something to do with the way the minimization package works, which I know very little about.

Thanks again for any help,
Matt

4

Matt,

  • If you suspect the TF1 function, why not replace it with line
    of code .
  • In your Gauss you divide through parameter [3] . No protection
    is provided when [3] gets very small or even zero .

Eddy

5

I would replace it, but I still need to somehow accomplish a gaussian fit of a histogram. I’m not sure how to do that without using TF1. I will try to redefine the function to make use of the predefined fit functions ‘gaus’ and ‘pol0’ and see if this behaves better with TMinuit.

Parameter [3] isn’t allowed to go to zero. The parameters are all limited automatically inside my fitting class (it’s long so I didn’t want to paste it all here), which I use in many other programs to fit exactly these same histograms. I even tested it in this specific context, and it works fine until TMinuit tries to minimize after the TF1->Fit() has been run.

I’ll try to reproduce the problem in a small example code.

Thanks for the help…
gamFit.h (488 Bytes)
gamFit.C (2.42 KB)

6

You should check if your fit gaus fit of a TF1 fails. In that case you can get crazy parameter values which then will give you problems in your chi2 minimization.

Lorenzo

7

I did check the fits by turning off the quiet option and they are fitting very well. What they produce doesn’t matter now anyway since I’m not using their output for the chi2 minimization right now. Just running them crashes the minimization.

I tried to see if it would work better with a simpler fitting routine coded inline. So I used

ff = new TF1("ff","pol0(0)+gaus(1)");
for the definition and then in the fcn used by TMinuit I run

sega_e_calDop[ii][nuc]->Fit("ff","Q","",lineE[nuc][line]-80,lineE[nuc][line]+80);
This fit runs, but running this fit similarly corrupts the larger minimization.

I tried to create a smaller simpler code that crashes in the same way, but for my smaller code, it has no problem performing a fit inside the minimization function (I’ve attached this program below).

I have now discovered that if I rerun the minimization a second time, after it fails once in the same root session, then the minimization will not be corrupted by the fitting inside the objective function.

So for now I’ll just run it twice.

I also found a related topic in the forums:
http://root.cern.ch/phpBB2/viewtopic.php?t=4043&highlight=pol1+gaus
I noticed that you also responded to that post Lorenzo, do you think it would make any difference in my case if I were to use

TFitter * minuit = new TFitter(6);

instead of

TMinuit *gMinuit = new TMinuit(28);

Is there any difference between the two?

Thanks,
Matt
test.root (6.71 KB)
exFitInFit.C (2.46 KB)

8

Using the class TFitter or TMinuit, should be the same. TFitter is just a wrapper class for TMinuit.
IA difference would be if you use the class TVirtualFitter. This uses a global instance which is also used in TH1::Fit and there you could have a conflict between the inner and outer fit.
I would suggest you run all your code in compiled mode and fix first any compilation problem.

Lorenzo