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