RE: [ROOT] Using TH1 Fit method within TMinuit

From: Greg Landsberg (landsberg@hep.brown.edu)
Date: Thu Mar 06 2003 - 14:11:26 MET


Dear Rene,

Many thanks for the fix and explanation. In the meantime I found an
interim solution, which is to save the gMinuit pointer, then create a
new TMinuit object just before calling the h->Fit(), destroy it
afterwards, and assign the saved value back to gMinuit. I am not sure
that it is a complete fix, as MIGRAD acts a bit strange after that and
barely converges. (Perhaps the error matrix might be still destroyed
after such a switch.) Perhaps someone could look into that in more
detail as a remedy for the non-recursivity of TMinuit in the CINT mode.

Please, keep me posted on future TMinuit fixes related to this problem.

Best,

Greg

> -----Original Message-----
> From: Rene Brun [mailto:brun@pcbrun.cern.ch]
> Sent: Wednesday, March 05, 2003 1:59 PM
> To: Greg Landsberg
> Cc: roottalk@pcroot.cern.ch
> Subject: RE: [ROOT] Using TH1 Fit method within TMinuit
> 
> 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