fit within a fit trouble

From: Stephen Bailey (bailey@physics.harvard.edu)
Date: Wed Mar 22 2000 - 21:09:56 MET 

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