#include "TBinomialEfficiencyFitter.h"
#include "TMath.h"
#include "TPluginManager.h"
#include "TROOT.h"
#include "TH1.h"
#include "TF1.h"
#include "TVirtualFitter.h"
#include "TEnv.h"
#include <limits>
TVirtualFitter *TBinomialEfficiencyFitter::fgFitter = 0;
ClassImp(TBinomialEfficiencyFitter)
TBinomialEfficiencyFitter::TBinomialEfficiencyFitter() {
fNumerator = 0;
fDenominator = 0;
fFunction = 0;
fFitDone = kFALSE;
fAverage = kFALSE;
fRange = kFALSE;
}
TBinomialEfficiencyFitter::TBinomialEfficiencyFitter(const TH1 *numerator, const TH1 *denominator)
{
Set(numerator,denominator);
}
TBinomialEfficiencyFitter::~TBinomialEfficiencyFitter() {
delete fgFitter; fgFitter = 0;
}
void TBinomialEfficiencyFitter::Set(const TH1 *numerator, const TH1 *denominator)
{
fNumerator = (TH1*)numerator;
fDenominator = (TH1*)denominator;
fFitDone = kFALSE;
fAverage = kFALSE;
fRange = kFALSE;
}
TVirtualFitter* TBinomialEfficiencyFitter::GetFitter() {
return fgFitter;
}
Int_t TBinomialEfficiencyFitter::Fit(TF1 *f1, Option_t* option)
{
TString opt = option;
opt.ToUpper();
fAverage = opt.Contains("I");
fRange = opt.Contains("R");
fFunction = (TF1*)f1;
Int_t i, npar;
npar = f1->GetNpar();
if (npar <= 0) {
Error("Fit", "function %s has illegal number of parameters = %d", f1->GetName(), npar);
return -3;
}
if (!f1) return -1;
if (!fNumerator || !fDenominator) {
Error("Fit","No numerator or denominator histograms set");
return -5;
}
if (f1->GetNdim() != fNumerator->GetDimension()) {
Error("Fit","function %s dimension, %d, does not match histogram dimension, %d",
f1->GetName(), f1->GetNdim(), fNumerator->GetDimension());
return -4;
}
Int_t maxpar = npar;
if (!fgFitter) {
TPluginHandler *h;
TString fitterName = TVirtualFitter::GetDefaultFitter();
if (fitterName == "")
fitterName = gEnv->GetValue("Root.Fitter","Minuit");
if ((h = gROOT->GetPluginManager()->FindHandler("TVirtualFitter", fitterName ))) {
if (h->LoadPlugin() == -1)
return 0;
fgFitter = (TVirtualFitter*) h->ExecPlugin(1, maxpar);
}
if (!fgFitter) printf("ERROR fgFitter is NULL\n");
}
fgFitter->SetObjectFit(this);
fgFitter->Clear();
fgFitter->SetFCN(BinomialEfficiencyFitterFCN);
Int_t nfixed = 0;
Double_t al,bl,arglist[100];
for (i = 0; i < npar; i++) {
f1->GetParLimits(i,al,bl);
if (al*bl != 0 && al >= bl) {
al = bl = 0;
arglist[nfixed] = i+1;
nfixed++;
}
if (f1->GetParError(i) <= 0) f1->SetParError(i, 0.01);
fgFitter->SetParameter(i, f1->GetParName(i),
f1->GetParameter(i),
f1->GetParError(i), al,bl);
}
if(nfixed > 0)fgFitter->ExecuteCommand("FIX",arglist,nfixed);
Double_t plist[1];
plist[0] = 0.5;
fgFitter->ExecuteCommand("SET ERRDEF",plist,1);
fFitDone = kTRUE;
Int_t result = fgFitter->ExecuteCommand("MINIMIZE",0,0);
char parName[50];
Double_t par, we;
Double_t eplus,eminus,eparab,globcc,werr;
for (i=0;i<npar;i++) {
fgFitter->GetParameter(i,parName, par,we,al,bl);
fgFitter->GetErrors(i,eplus,eminus,eparab,globcc);
if (eplus > 0 && eminus < 0) werr = 0.5*(eplus-eminus);
else werr = we;
f1->SetParameter(i,par);
f1->SetParError(i,werr);
}
f1->SetNDF(f1->GetNumberFitPoints()-npar+nfixed);
return result;
}
void TBinomialEfficiencyFitter::ComputeFCN(Int_t& , Double_t* ,
Double_t& f, Double_t* par, Int_t ) {
int lowbin = fDenominator->GetXaxis()->GetFirst();
int highbin = fDenominator->GetXaxis()->GetLast();
fFunction->SetParameters(par);
if (fRange) {
double xmin, xmax;
fFunction->GetRange(xmin, xmax);
lowbin = fDenominator->GetXaxis()->FindBin(xmin);
highbin = fDenominator->GetXaxis()->FindBin(xmax);
}
f = 0.;
Int_t npoints = 0;
Double_t nmax = 0;
for (int bin = lowbin; bin <= highbin; ++bin) {
double xlow = fDenominator->GetBinLowEdge(bin);
double xup = fDenominator->GetBinLowEdge(bin+1);
double nDen = fDenominator->GetBinContent(bin);
double nNum = fNumerator->GetBinContent(bin);
if (nDen> nmax) nmax = nDen;
if (nDen <= 0.) continue;
npoints++;
double mu = (fAverage) ?
fFunction->Integral(xlow, xup) / (xup-xlow) :
fFunction->Eval(fDenominator->GetBinCenter(bin));
if (nNum != 0.)
if (mu > 0.)
f -= nNum * TMath::Log(mu);
else
f -= nmax * -1E30;
if (nDen - nNum != 0.)
if (1. - mu > 0.)
f -= (nDen - nNum) * TMath::Log(1. - mu);
else
f -= nmax * -1E30;
}
fFunction->SetNumberFitPoints(npoints);
fFunction->SetChisquare(f);
}
void BinomialEfficiencyFitterFCN(Int_t& npar, Double_t* gin, Double_t& f, Double_t* par, Int_t flag) {
TBinomialEfficiencyFitter* fitter = dynamic_cast<TBinomialEfficiencyFitter*>(TBinomialEfficiencyFitter::GetFitter()->GetObjectFit());
if (!fitter) {
fitter->Error("binomialFCN","Invalid fit object encountered!");
return;
}
fitter->ComputeFCN(npar, gin, f, par, flag);
}
This page has been automatically generated. If you have any comments or suggestions about the page layout send a mail to ROOT support, or contact the developers with any questions or problems regarding ROOT.