Hi Rene,
I have isolated the problem with TMinuit and make a simple test program
that reproduces the error. The program fills a histogram h with a
Gaussian distribution, which has width and mean dependent on the values
of the parameters of the minimization. The function fcn is constructed
in such a way to have a minimum for the values of the first parameter of
0.9, and the second one of 1.1.
The main function fittest can be called either as fittest(0), in which
case the parameters of the Gaussian distribution are obtained from the
fit of the histogram via h->Fit("gaus"), or as fittest(1), where
parameters are returned via taking a mean and RMS of the histogram,
which yields similar results.
The code works just fine when fittest(1) is used and returns correct
value of the minimum. However, when fittest(0) is used, the TMinuit is
switched into minimizing the function "gaus" from the h->Fit("gaus") as
soon as the first instance of h->Fit appears. One can explicitly see
this via the printout loop put in fcn: the printout appears 10 times (as
expected) for fittest(1), and only once for fittest(0), which means that
fcn is never called on subsequent minimization steps.
I'd really appreciate if you can help to figure out what's wrong with
either my test code or the TMinuit itself.
Many thanks,
Greg
//----------- The code starts here -------------
#include "TH1.h"
#include <iostream>
#include <stdlib.h>
#include "TSystem.h"
#include "TMath.h"
#include "TMinuit.h"
using namespace std;
Bool_t fit = true;
Int_t ncall = 0;
void fittest(Int_t flag = 0)
{
Double_t x[2], err[2], tmp;
char name[10];
//
if (flag != 0) fit = false;
//
TMinuit *minuit = new TMinuit();
minuit->DefineParameter(0, "Mass Scale", 1., 0.05, 0.8, 1.2);
minuit->DefineParameter(1, "Width Scale", 1., 0.05, 0.8, 1.2);
//
minuit->SetFCN(fcn);
Double_t arglist[10];
Int_t ierflg = 0;
//
arglist[0] = 1;
minuit->mnexcm("SET ERR", arglist ,1, ierflg);
//
arglist[0] = 2000;
arglist[1] = 1.;
if (ierflg == 0) minuit->mnexcm("MINIMIZE", arglist ,1, ierflg);
//
if (ierflg != 0) cout << "Error in mnexcm: " << ierflg << endl;
//
for (Int_t i = 0; i < 2; i++)
{
minuit->GetParameter(i, x[i], err[i]);
cout << "Parameter " << i << " = " << x[i] << " +- " << err[i] <<
endl;
}
}
void fastfit(TH1F *h, Double_t *pars)
{
if (fit)
{
h->Fit("gaus","0","",0.,50.);
TF1 *f=h->GetFunction("gaus");
if (f == 0) cout << "Error: no function gaus exist\n";
else for (Int_t i = 0; i < 3; i++) pars[i] = f->GetParameter(i);
}
else
{
pars[0] = h->GetMaximum();
pars[1] = h->GetMean();
pars[2] = h->GetRMS();
}
return;
}
void fcn(Int_t &npar, Double_t *gin, Double_t &f, Double_t *par, Int_t
iflag)
{
Float_t mass, width;
getpar(par, mass, width);
ncall++;
//
f = 1000.*width + (mass - 22.5)*(mass - 22.5);
if (ncall < 10) cout << "iflag, mass, width, f = " << iflag << ", "
<< mass << ", " << width << ", " << f << endl;
//
return;
}
void getpar(Double_t *x, Float_t &mass, Float_t &width)
{
Double_t par[6];
//
TH1F *h = new TH1F("peak","peak",50,0.,50.);
TRandom *myrandom = new TRandom(654321);
for (Int_t i=0; i<10000;i++) // Fill the histogram with a Gaussian
{
Float_t M = myrandom->Gaus(25.*x[0],5.0+10.*(x[1]-1.1)*(x[1]-1.1));
h->Fill(M);
}
//
fastfit(h,par);
mass = par[1];
width = par[2];
delete h;
delete myrandom;
return;
}
//-------------------------- The code ends here ---------------------
> -----Original Message-----
> From: Rene Brun [mailto:brun@pcbrun.cern.ch]
> Sent: Tuesday, March 04, 2003 4:09 AM
> To: Greg Landsberg
> Cc: roottalk@pcroot.cern.ch
> Subject: Re: [ROOT] Using TH1 Fit method within TMinuit
>
> Hi Greg,
>
> You mention TMinuit and TVirtualFitter, but you do not indicate which
> one you are really using in your external loop and in your internal
loop.
> The following should work:
> - you create your own instance of TMinuit *myminuit = new
TMinuit(...);
> - you then use only the myminuit pointer (not gMinuit)
> - you let TH1::Fit calls TVirtualFitter
>
> Rene Brun
>
> On
> Mon, 3 Mar 2003, Greg Landsberg wrote:
>
> > Dear Rooters,
> >
> > I have faced the following problem when trying to use TMinuit or
> > TVirtualFitter class to do a complicated minimization. My
minimization
> > function, fcn, which is passed to TMinuit/TVirtualFitter via the
SetFCN
> > method fills a histogram that it further fits to a certain function
via
> > standard TH1F Fit() method. The results of this 'internal' fit are
used
> > in order to estimate the return value of fcn.
> >
> > However, as soon as the first instance of the TH1F Fit() is called,
the
> > Minuit is switched to optimize the histogram, not the function fcn!
This
> > is as if there were just a single copy of the TMinuit object that
exists
> > in memory, and is superceded with whatever the latest initialization
> > setting that took place (i.e. via TH1F Fit).
> >
> > This sounds like a bug to me, rather than a designed behavior. I
wonder
> > if anyone faced a similar problem and knows a workaround.
> >
> > Please, CC: your reply directly to me, as I don't always read the
> > roottalk.
> >
> > Many thanks,
> >
> > Greg Landsberg
> >
This archive was generated by hypermail 2b29 : Thu Jan 01 2004 - 17:50:09 MET