RE: [ROOT] Using TH1 Fit method within TMinuit

From: Rene Brun (Rene.Brun@cern.ch)
Date: Wed Mar 05 2003 - 19:59:16 MET


Hi Greg,

I see that the call to TMinuit recursively is not supported by
the interpreter. You must run with the compiler.
I have modified your script in attachement. Do
 root > .x fittest.C+(0)

I have modified the TFitter class to remove all dependencies on gMinuit
and use an internal member fMinuit instead. This is required in case
you would like to use more than two fitters in the inner loop (not your 
case).
TMinuit should be improved to accept recursivity with the interpreter.

Rene Brun

On 
Tue, 4 
Mar 2003, Greg Landsberg wrote:

> 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