[ROOT] User function for minimization of histogram fit

From: Chris Roat (croat@SLAC.stanford.edu)
Date: Thu Jul 17 2003 - 19:04:43 MEST 

Hi,

How can I pass TMinuit a user function for minimization when
fitting a histogram?  I'm having difficulties adapting the
tutorial Ifit.C (MINUIT is called directly, rather than
fitting a histogram).

Below is a macro with a modified chi-square function.  If I
run this program (RH 7.2, ROOT 3.05/05), I don't see a
difference if I comment out the SetFCN call.  Also, I added
the printf statement and it never gets executed.

Thanks for any help,
Chris



void diffMiniFunc() {
   TH1F *h = new TH1F("h","My histogram",100,-3,3);
   h->FillRandom("gaus",6000);
   h->Fit("gaus"); // fit first to get a gMinuit object
   gMinuit->SetFCN( MyH1FitChisquare );
   h->Fit("gaus");
}

void MyH1FitChisquare(Int_t &npar, Double_t *gin, Double_t &f, Double_t *u, Int_t flag)
{
//           Minimization function for H1s using a Chisquare method
//           ======================================================

  printf("I've been called!!\n");

   Double_t cu,eu,fu,fsum;
   Double_t dersum[100], grad[100];
   Double_t x[3];
   Int_t bin,binx,biny,binz,k;
   Axis_t binlow, binup, binsize;

   Int_t npfits = 0;


   TVirtualFitter *hFitter = TVirtualFitter::GetFitter();
   TH1 *hfit = (TH1*)hFitter->GetObjectFit();
   TF1 *f1   = (TF1*)hFitter->GetUserFunc();
   Foption_t Foption = hFitter->GetFitOption();

   f1->InitArgs(x,u);
   npar = f1->GetNpar();
   if (flag == 2) for (k=0;k<npar;k++) dersum[k] = gin[k] = 0;
   f = 0;
   Int_t hxfirst = hFitter->GetXfirst();
   Int_t hxlast  = hFitter->GetXlast();
   Int_t hyfirst = hFitter->GetYfirst();
   Int_t hylast  = hFitter->GetYlast();
   Int_t hzfirst = hFitter->GetZfirst();
   Int_t hzlast  = hFitter->GetZlast();
   TAxis *xaxis  = hfit->GetXaxis();
   TAxis *yaxis  = hfit->GetYaxis();
   TAxis *zaxis  = hfit->GetZaxis();

   for (binz=hzfirst;binz<=hzlast;binz++) {
      x[2]  = zaxis->GetBinCenter(binz);
      for (biny=hyfirst;biny<=hylast;biny++) {
         x[1]  = yaxis->GetBinCenter(biny);
         for (binx=hxfirst;binx<=hxlast;binx++) {
            x[0]  = xaxis->GetBinCenter(binx);
            if (!f1->IsInside(x)) continue;
            bin = hfit->GetBin(binx,biny,binz);
            cu  = hfit->GetBinContent(bin);
            TF1::RejectPoint(kFALSE);
            if (Foption.Integral) {
               binlow  = xaxis->GetBinLowEdge(binx);
               binsize = xaxis->GetBinWidth(binx);
               binup   = binlow + binsize;
               fu      = f1->Integral(binlow,binup,u)/binsize;
            } else {
               fu = f1->EvalPar(x,u);
            }
            if (TF1::RejectedPoint()) continue;
            if (Foption.W1) {
               eu = 1;
            } else {
               //eu  = hfit->GetBinError(bin);
	      // take error from value of fit
               eu  = TMath::Sqrt(fu);
               if (eu <= 0) continue;
            }
            if (flag == 2) {
               for (k=0;k<npar;k++) dersum[k] += 1; //should be the derivative
            }
            npfits++;
            if (flag == 2) {
               for (k=0;k<npar;k++) grad[k] += dersum[k]*(fu-cu)/eu; dersum[k] = 0;
            }
            fsum = (cu-fu)/eu;
            f += fsum*fsum;
         }
      }
   }
   f1->SetNumberFitPoints(npfits);
}

This archive was generated by hypermail 2b29 : Thu Jan 01 2004 - 17:50:13 MET