Hi Chris,
See documentation of TH1::Fit to specify your own minimisation model.
There is an example. In particular, you must specify the option "U"
See your script modified in the attachement.
Note that your script has problems. I did not fix them, but your
function is called.
Rene Brun
On
Thu, 17 Jul 2003, Chris Roat wrote:
> 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