source: trunk/WuerzburgSoft/Thomas/mphys/MFit.cc@ 14912

#include "MFit.h"

#include <iomanip.h>

#include <TF1.h>
#include <TH1.h>

#include <TMinuit.h>

ClassImp(MFit);

TF1 *MFit::fgFunc=NULL;
TH1 *MFit::fgHist=NULL;

Int_t MFit::fgBinMin;
Int_t MFit::fgBinMax;

void MFit::Func(Int_t &npar, Double_t *gin, Double_t &f, Double_t *par, Int_t flag)
{
    Double_t chisq = 0;

    for (int i=fgBinMin+1; i<fgBinMax; i++)
    {
        Double_t x = fgHist->GetBinCenter(i);

        Double_t y = fgFunc->EvalPar(&x, par);
        Double_t d = (y - fgHist->GetBinContent(i))/fgHist->GetBinError(i);

        chisq += d*d;
    }

    f = chisq / fgFunc->GetNDF();
    fgFunc->SetChisquare(f);
}

void MFit::FuncLog(Int_t &npar, Double_t *gin, Double_t &f, Double_t *par, Int_t flag)
{
    Double_t chisq = 0;

    for (int i=fgBinMin+1; i<fgBinMax; i++)
    {
        Double_t x1 = fgHist->GetBinLowEdge(i+1);
        Double_t x0 = fgHist->GetBinLowEdge(i);

        Double_t x = sqrt(x1*x0);

        Double_t y = fgFunc->EvalPar(&x, par);
        Double_t d = (y - fgHist->GetBinContent(i))/fgHist->GetBinError(i);

        chisq += d*d;
    }

    f = chisq / fgFunc->GetNDF();
    fgFunc->SetChisquare(f);
}

void MFit::DeleteF()
{
    if (fFunc)
        if (fFunc->TestBit(kIsOwner))
            delete fFunc;

    if (fMinuit)
        delete fMinuit;
}

void MFit::DeleteH()
{
    if (fHist)
        if (fHist->TestBit(kIsOwner))
            delete fHist;
}

void MFit::SetFunc(TF1 *f)
{
    DeleteF();

    fFunc = f;

    fMinuit = new TMinuit(fFunc->GetNpar());
    fMinuit->SetPrintLevel(-1);
}

void MFit::SetFunc(TString formula, Axis_t min, Axis_t max)
{
    TF1 *f = new TF1("", formula, min, max);
    f->SetBit(kIsOwner);

    SetFunc(f);
}

void MFit::SetRange(Axis_t min, Axis_t max)
{
    if (fFunc)
        fFunc->SetRange(min, max);
}

void MFit::SetHist(TH1 *h, Bool_t candelete)
{
    DeleteH();

    fHist = h;

    if (candelete)
        fHist->SetBit(kIsOwner);
}

void MFit::Fit(Bool_t log)
{
    fgBinMin = fHist->FindBin(fFunc->GetXmin()); // excluded
    fgBinMax = fHist->FindBin(fFunc->GetXmax()); // excluded

    fFunc->SetNDF(fgBinMax-fgBinMin-1);

    fgHist = fHist;
    fgFunc = fFunc;

    fMinuit->SetFCN(log ? FuncLog : Func);
    fMinuit->Migrad();

    for (int i=0; i<fFunc->GetNpar(); i++)
    {
        Double_t par, err;
        fMinuit->GetParameter(i, par, err);
        fFunc->SetParameter(i, par);
        fFunc->SetParError(i, err);
    }
}

void MFit::Print(Option_t *opt) const
{
    for (int i=0; i<fFunc->GetNpar(); i++)
        cout << "Par["<<i<<"]: " << GetParameter(i) << " +- " << GetParError(i) << endl;
    cout << "ChiSq:   " << fFunc->GetChisquare() << endl;
    cout << "NDF:     " << fFunc->GetNDF() << endl;
}

void MFit::SetParameter(Int_t n, TString name, Axis_t start, Axis_t min, Axis_t max)
{
    if (n>=fFunc->GetNpar())
        return;

    fMinuit->DefineParameter(n, name, start, min/10,  min, max);
}

Double_t MFit::GetParameter(Int_t n) const
{
    return fFunc->GetParameter(n);
}

Double_t MFit::GetParError(Int_t n) const
{
    return fFunc->GetParError(n);
}

Double_t MFit::operator[](Int_t n) const
{
    return fFunc->GetParameter(n);
}

TF1 *MFit::DrawCopy(Option_t *o) const
{
    return fFunc->DrawCopy(o);
}