Re: fit within a fit trouble

From: Rene Brun (Rene.Brun@cern.ch)
Date: Thu Mar 23 2000 - 10:24:15 MET


Hi Stephen,
You are making a good remark about a fit within a fit.
I have modified your program to make it possible.
I still have to investigate if it is possible to automatically
swap the fitters without any intervention in the user's code.

Rene Brun


// Toy example of fit within a fit.
// The TNtuple nt contains a Gaussian mass (m) distribution with an
// exponential lifetime (t).  It is on top of a background with a
// random mass and a shorter exponential lifetime distribution.

// The function otto() trys to maximize the statistical power S^2/(S+B)
// with a lifetime cut.  The minimization function projects the ntuple
// into a histogram and fits the histogram, extracting S and B and
// returns -S*S/(S+B);

// But apparently fitting the histogram within the running fit causes
// the fitter to crash with a segmentation violation

TNtuple* nt;
TFitter *fitter, *hfitter; //<============
TMinuit *fit, *hfit;  //<============
//extern TMinuit *gMinuit;
void initOtto(void) {
  nt = new TNtuple("nt", "Test Ntuple", "m:t");
  TRandom rand;

  // Fill with the signal and background
  for(int i = 0; i < 1000; i++) {
    nt.Fill( rand.Gaus(), rand.Exp(2) );
    nt.Fill( rand.Rndm() * 10 - 5, rand.Exp(1) );
  }

}


void otto(void) {

  // Create and set up the histogram fitter
  hfitter = new TFitter(25);  //<==============
  hfit = gMinuit;
  // Create and set up the fitter
  fitter  = new TFitter(1); //<===========
  fitter->SetFCN(fcn);
  fit = gMinuit;
  
  Double_t arg[5];
 
  // Quiet it down a bit
  arg[0] = -1;  fitter->ExecuteCommand("SET PRI", arg, 1);

  fitter->SetParameter(0, "t cut", 0.1, 0.01, 0, 0);
 
  // Now ready for minimization step
  arg[0] = 500;
  arg[1] = 1.0;
  fitter->ExecuteCommand("MIGRAD", arg, 2);
 
  // Print results
  Double_t amin;
  fitter->PrintResults(1, amin);

}  // end otto()
 

void fcn(Int_t &npar, Double_t *gin, Double_t &r, Double_t *par, Int_t
flag)
{
  TH1S h1("h1", "Test Histogram", 50, -5, 5);
  TF1 f1("f1", "gaus(0)+pol1(3)", -5, 5);
  f1.SetParameters(40, 0, 1, 2, 0);
  
  char cut[50];
  sprintf(cut, "t>%g", par[0]);

  h1.Clear();
  nt->Project("h1", "m", cut);
  TVirtualFitter::SetFitter(hfitter,25); //<==============
  gMinuit = hfit;
  h1.Fit("f1","0");
  TVirtualFitter::SetFitter(fitter,1);  //<============
  gMinuit = fit;

  // Normalization here is probably off, but that isn't the point.
  Double_t s = f1->GetParameter(0) * 2.507 * f1->GetParameter(2);
  Double_t b = f1->GetParameter(3);


  r = -(s*s)/(s+b);
  printf("s=%g, b=%g, r=%g\n",s,b,r);
}  // end fcn()


Stephen Bailey wrote:
> 
> Hi.
> 
> I'm trying to use a fit within a fit: the minimization function
> of my fit involves another fit to a histogram, but I get a segmentation
> violation using ROOT 2.23/11 on Linux.
> 
> Here is an example session followed by my code.  Am I doing something
> wrong or is this a bug?  Is there a solution or workaround?  Thanks.
> 
> Stephen
> 
> root [0] .L otto.cc
> root [1] initOtto();
> root [2] otto();
>  FCN=34.3176 FROM MIGRAD    STATUS=CONVERGED     165 CALLS         166
> TOTAL
>                      EDM=1.31629e-05    STRATEGY= 1      ERROR MATRIX
> ACCURATE
>   EXT PARAMETER                                   STEP         FIRST
>   NO.   NAME      VALUE            ERROR          SIZE      DERIVATIVE
>    1  p0           7.29523e+01   3.36079e+00   8.08107e-03  -3.91228e-04
>    2  p1          -6.99217e-02   4.98261e-02   1.51244e-04  -3.12920e-02
>    3  p2           1.05729e+00   3.89826e-02   9.41497e-05  -5.54527e-02
>    4  p3           1.72997e+01   8.81791e-01   2.06990e-03   2.13463e-03
>    5  p4          -3.70281e-01   2.21413e-01   6.13668e-04  -1.58787e-02
> <TCanvas::MakeDefCanvas>: created default TCanvas with name c1
> 
>  *** Break *** segmentation violation
> Root > Function otto() busy flag cleared
> 
> // Toy example of fit within a fit.
> // The TNtuple nt contains a Gaussian mass (m) distribution with an
> // exponential lifetime (t).  It is on top of a background with a
> // random mass and a shorter exponential lifetime distribution.
> 
> // The function otto() trys to maximize the statistical power S^2/(S+B)
> // with a lifetime cut.  The minimization function projects the ntuple
> // into a histogram and fits the histogram, extracting S and B and
> // returns -S*S/(S+B);
> 
> // But apparently fitting the histogram within the running fit causes
> // the fitter to crash with a segmentation violation
> 
> TNtuple* nt;
> 
> void initOtto(void) {
>   nt = new TNtuple("nt", "Test Ntuple", "m:t");
>   TRandom rand;
> 
>   // Fill with the signal and background
>   for(int i = 0; i < 1000; i++) {
>     nt.Fill( rand.Gaus(), rand.Exp(2) );
>     nt.Fill( rand.Rndm() * 10 - 5, rand.Exp(1) );
>   }
> 
> }
> 
> void otto(void) {
> 
>   // Create and set up the fitter
>   TFitter *fitter = new TFitter(1);
>   fitter->SetFCN(fcn);
> 
>   Double_t arg[5];
> 
>   // Quiet it down a bit
>   arg[0] = -1;  fitter->ExecuteCommand("SET PRI", arg, 1);
> 
>   fitter->SetParameter(0, "t cut", 0.1, 0.01, 0, 0);
> 
>   // Now ready for minimization step
>   arg[0] = 500;
>   arg[1] = 1.0;
>   fitter->ExecuteCommand("MIGRAD", arg, 2);
> 
>   // Print results
>   Double_t amin;
>   fitter->PrintResults(1, amin);
> 
> }  // end otto()
> 
> 
> void fcn(Int_t &npar, Double_t *gin, Double_t &r, Double_t *par, Int_t flag)
> {
>   TH1S h1("h1", "Test Histogram", 50, -5, 5);
>   TF1 f1("f1", "gaus(0)+pol1(3)", -5, 5);
>   f1.SetParameters(40, 0, 1, 2, 0);
> 
>   char cut[50];
>   sprintf(cut, "t>%g", par[0]);
> 
>   h1.Clear();
>   nt->Project("h1", "m", cut);
> 
>   h1.Fit("f1");
> 
>   // Normalization here is probably off, but that isn't the point.
>   Double_t s = f1->GetParameter(0) * 2.507 * f1->GetParameter(2);
>   Double_t b = f1->GetParameter(3);
> 
>   r = -(s*s)/(s+b);
> }  // end fcn()



This archive was generated by hypermail 2b29 : Tue Jan 02 2001 - 11:50:21 MET