#include "Riostream.h"
#include "TH1.h"
#include "TMath.h"
#include "TClass.h"
#include "TFractionFitter.h"
TVirtualFitter *fractionFitter=0;
ClassImp(TFractionFitter)
TFractionFitter::TFractionFitter() :
fFitDone(kFALSE),fData(0), fPlot(0) {
fractionFitter = 0;
fIntegralMCs = 0;
fFractions = 0;
fNpfits = 0;
fNDF = 0;
fChisquare = 0;
}
TFractionFitter::TFractionFitter(TH1* data, TObjArray *MCs) :
fFitDone(kFALSE), fChisquare(0), fPlot(0) {
fData = data;
fLowLimitX = 1;
fHighLimitX = fData->GetNbinsX();
if (fData->GetDimension() > 1) {
fLowLimitY = 1;
fHighLimitY = fData->GetNbinsY();
if (fData->GetDimension() > 2) {
fLowLimitZ = 1;
fHighLimitZ = fData->GetNbinsZ();
}
}
fNpar = MCs->GetEntries();
Int_t par;
for (par = 0; par < fNpar; ++par) {
fMCs.Add(MCs->At(par));
TString s = Form("Prediction for MC sample %i",par);
TH1* pred = (TH1*) ((TH1*)MCs->At(par))->Clone(s);
pred->SetTitle(s);
fAji.Add(pred);
}
fIntegralMCs = new Double_t[fNpar];
fFractions = new Double_t[fNpar];
CheckConsistency();
fWeights.Expand(fNpar);
fractionFitter = TVirtualFitter::Fitter(this, fNpar);
fractionFitter->Clear();
fractionFitter->SetObjectFit(this);
fractionFitter->SetFCN(TFractionFitFCN);
Double_t defaultFraction = 1.0/((Double_t)fNpar);
Double_t defaultStep = 0.01;
for (par = 0; par < fNpar; ++par) {
TString name("frac"); name += par;
fractionFitter->SetParameter(par, name.Data(), defaultFraction, defaultStep, 0, 0);
}
}
TFractionFitter::~TFractionFitter() {
delete fractionFitter;
delete[] fIntegralMCs;
delete[] fFractions;
}
void TFractionFitter::SetData(TH1* data) {
fData = data;
fFitDone = kFALSE;
CheckConsistency();
}
void TFractionFitter::SetMC(Int_t parm, TH1* MC) {
CheckParNo(parm);
fMCs.RemoveAt(parm);
fMCs.AddAt(MC,parm);
fFitDone = kFALSE;
CheckConsistency();
}
void TFractionFitter::SetWeight(Int_t parm, TH1* weight) {
CheckParNo(parm);
if (fWeights[parm]) {
fWeights.RemoveAt(parm);
}
if (weight) {
if (weight->GetNbinsX() != fData->GetNbinsX() ||
(fData->GetDimension() > 1 && weight->GetNbinsY() != fData->GetNbinsY()) ||
(fData->GetDimension() > 2 && weight->GetNbinsZ() != fData->GetNbinsZ())) {
Error("SetWeight","Inconsistent weights histogram for source %d", parm);
return;
}
TString ts = "weight hist: "; ts += weight->GetName();
fWeights.AddAt(weight,parm);
}
}
TVirtualFitter* TFractionFitter::GetFitter() const {
return fractionFitter;
}
void TFractionFitter::CheckParNo(Int_t parm) const {
if (parm < 0 || parm > fNpar) {
Error("CheckParNo","Invalid parameter number %d",parm);
}
}
void TFractionFitter::SetRangeX(Int_t low, Int_t high) {
if (low <=0 || high <= 0) {
Error("SetRangeX","Invalid fit range");
return;
}
fLowLimitX = (low > 0) ? low : 1;
fHighLimitX = (high <= fData->GetNbinsX()) ? high : fData->GetNbinsX();
CheckConsistency();
}
void TFractionFitter::ReleaseRangeX() {
fLowLimitX = 1;
fHighLimitX = fData->GetNbinsX();
CheckConsistency();
}
void TFractionFitter::SetRangeY(Int_t low, Int_t high) {
if (fData->GetDimension() < 2) {
Error("SetRangeY","Y range cannot be set for 1D histogram");
return;
}
if (low <=0 || high <= 0) {
Error("SetRangeY","Invalid fit range");
return;
}
fLowLimitY = (low > 0) ? low : 1;
fHighLimitY = (high <= fData->GetNbinsY()) ? high : fData->GetNbinsY();
CheckConsistency();
}
void TFractionFitter::ReleaseRangeY() {
fLowLimitY = 1;
fHighLimitY = fData->GetNbinsY();
CheckConsistency();
}
void TFractionFitter::SetRangeZ(Int_t low, Int_t high) {
if (fData->GetDimension() < 3) {
Error("SetRangeZ","Z range cannot be set for 1D or 2D histogram");
return;
}
if (low <=0 || high <= 0) {
Error("SetRangeZ","Invalid fit range");
return;
}
fLowLimitZ = (low > 0) ? low : 1;
fHighLimitZ = (high <= fData->GetNbinsZ()) ? high : fData->GetNbinsZ();
CheckConsistency();
}
void TFractionFitter::ReleaseRangeZ() {
fLowLimitZ = 1;
fHighLimitZ = fData->GetNbinsZ();
CheckConsistency();
}
void TFractionFitter::Constrain(Int_t parm, Double_t low, Double_t high) {
CheckParNo(parm);
Double_t plist[3];
plist[0] = (Double_t) parm;
plist[1] = low;
plist[2] = high;
fractionFitter->ExecuteCommand("SET LIMIT", plist, 3);
}
void TFractionFitter::UnConstrain(Int_t parm) {
CheckParNo(parm);
Double_t plist[3];
plist[0] = (Double_t) parm;
plist[1] = 0;
plist[2] = 0;
fractionFitter->ExecuteCommand("SET LIMIT", plist, 3);
}
void TFractionFitter::CheckConsistency() {
if (! fData) {
Error("CheckConsistency","Nonexistent data histogram");
return;
}
Int_t minX, maxX, minY, maxY, minZ, maxZ;
Int_t x,y,z,par;
GetRanges(minX, maxX, minY, maxY, minZ, maxZ);
fIntegralData = 0;
fNpfits = 0;
for (z = minZ; z <= maxZ; ++z) {
for (y = minY; y <= maxY; ++y) {
for (x = minX; x <= maxX; ++x) {
fNpfits++;
fIntegralData += fData->GetBinContent(x, y, z);
}
}
}
if (fIntegralData <= 0) {
Error("CheckConsistency","Empty data histogram");
return;
}
TClass* cl = fData->Class();
fNDF = fNpfits - fNpar;
if (fNpar < 2) {
Error("CheckConsistency","Need at least two MC histograms");
return;
}
for (par = 0; par < fNpar; ++par) {
TH1 *h = (TH1*)fMCs.At(par);
if (! h) {
Error("CheckConsistency","Nonexistent MC histogram for source #%d",par);
return;
}
if ((! h->Class()->InheritsFrom(cl)) || h->GetNbinsX() != fData->GetNbinsX() ||
(fData->GetDimension() > 1 && h->GetNbinsY() != fData->GetNbinsY()) ||
(fData->GetDimension() > 2 && h->GetNbinsZ() != fData->GetNbinsZ())) {
Error("CheckConsistency","Histogram inconsistency for source #%d",par);
return;
}
fIntegralMCs[par] = 0;
for (z = minZ; z <= maxZ; ++z) {
for (y = minY; y <= maxY; ++y) {
for (x = minX; x <= maxX; ++x) {
fIntegralMCs[par] += h->GetBinContent(x, y, z);
}
}
}
if (fIntegralMCs[par] <= 0) {
Error("CheckConsistency","Empty MC histogram #%d",par);
}
}
}
Int_t TFractionFitter::Fit() {
Double_t plist[1];
plist[0] = 0.5;
fractionFitter->ExecuteCommand("SET ERRDEF",plist,1);
if (fPlot) {
delete fPlot; fPlot = 0;
}
fractionFitter->SetObjectFit(this);
Int_t status = fractionFitter->ExecuteCommand("MINIMIZE",0,0);
if (status == 0) fFitDone = kTRUE;
ComputeChisquareLambda();
return status;
}
void TFractionFitter::ErrorAnalysis(Double_t UP) {
if (! fFitDone) {
Error("ErrorAnalysis","Fit not yet performed");
return;
}
fractionFitter->SetObjectFit(this);
Double_t plist[1];
plist[0] = UP > 0 ? UP : 0.5;
fractionFitter->ExecuteCommand("SET ERRDEF",plist,1);
Int_t status = fractionFitter->ExecuteCommand("MINOS",0,0);
if (status != 0) {
Error("ErrorAnalysis","Error return from MINOS: %d",status);
}
}
void TFractionFitter::GetResult(Int_t parm, Double_t& value, Double_t& error) const {
CheckParNo(parm);
if (! fFitDone) {
Error("GetResult","Fit not yet performed");
return;
}
char parname[100];
Double_t vlow, vhigh;
fractionFitter->GetParameter(parm, parname, value, error, vlow, vhigh);
}
TH1* TFractionFitter::GetPlot() {
if (! fFitDone) {
Error("GetPlot","Fit not yet performed");
return 0;
}
if (! fPlot) {
Double_t plist[1];
plist[0] = 3;
fractionFitter->ExecuteCommand("CALL FCN", plist, 1);
}
return fPlot;
}
void TFractionFitter::GetRanges(Int_t& minX, Int_t& maxX, Int_t& minY, Int_t& maxY,
Int_t& minZ, Int_t& maxZ) const {
if (fData->GetDimension() < 2) {
minY = maxY = minZ = maxZ = 0;
minX = fLowLimitX;
maxX = fHighLimitX;
} else if (fData->GetDimension() < 3) {
minZ = maxZ = 0;
minX = fLowLimitX;
maxX = fHighLimitX;
minY = fLowLimitY;
maxY = fHighLimitY;
} else {
minX = fLowLimitX;
maxX = fHighLimitX;
minY = fLowLimitY;
maxY = fHighLimitY;
minZ = fLowLimitZ;
maxZ = fHighLimitZ;
}
}
void TFractionFitter::ComputeFCN(Int_t& , Double_t* ,
Double_t& f, Double_t* xx, Int_t flag)
{
Int_t bin, mc;
Int_t minX, maxX, minY, maxY, minZ, maxZ;
Int_t x,y,z;
GetRanges(minX, maxX, minY, maxY, minZ, maxZ);
for (mc = 0; mc < fNpar; ++mc) {
Double_t tot;
TH1 *h = (TH1*)fMCs[mc];
TH1 *hw = (TH1*)fWeights[mc];
if (hw) {
tot = 0;
for (z = minZ; z <= maxZ; ++z) {
for (y = minY; y <= maxY; ++y) {
for (x = minX; x <= maxX; ++x) {
Double_t weight = hw->GetBinContent(x, y, z);
if (weight <= 0) {
Error("ComputeFCN","Invalid weight encountered for MC source %d",mc);
return;
}
tot += weight * h->GetBinContent(x, y, z);
}
}
}
} else tot = fIntegralMCs[mc];
fFractions[mc] = xx[mc] * fIntegralData / tot;
}
if (flag == 3) {
TString ts = "Fraction fit to hist: "; ts += fData->GetName();
fPlot = (TH1*) fData->Clone(ts.Data());
fPlot->Reset();
}
Double_t result = 0;
for (z = minZ; z <= maxZ; ++z) {
for (y = minY; y <= maxY; ++y) {
for (x = minX; x <= maxX; ++x) {
bin = fData->GetBin(x, y, z);
int k0;
Double_t ti; Double_t aki;
FindPrediction(bin, fFractions, ti, k0, aki);
Double_t prediction = 0;
for (mc = 0; mc < fNpar; ++mc) {
TH1 *h = (TH1*)fMCs[mc];
TH1 *hw = (TH1*)fWeights[mc];
Double_t binPrediction;
Double_t binContent = h->GetBinContent(bin);
Double_t weight = hw ? hw->GetBinContent(bin) : 1;
if (k0 >= 0 && fFractions[mc] == fFractions[k0]) {
binPrediction = aki;
} else {
binPrediction = binContent > 0 ? binContent / (1+weight*fFractions[mc]*ti) : 0;
}
prediction += fFractions[mc]*weight*binPrediction;
result -= binPrediction;
if (binContent > 0 && binPrediction > 0)
result += binContent*TMath::Log(binPrediction);
if (flag == 3) {
((TH1*)fAji.At(mc))->SetBinContent(bin, binPrediction);
}
}
if (flag == 3) {
fPlot->SetBinContent(bin, prediction);
}
result -= prediction;
Double_t found = fData->GetBinContent(bin);
if (found > 0 && prediction > 0)
result += found*TMath::Log(prediction);
}
}
}
f = -result;
}
void TFractionFitter::FindPrediction(int bin, Double_t *fractions, Double_t &ti, int& k0, Double_t &aki) const {
Int_t par;
TH1 *hw;
for (par = 0; par < fNpar; ++par) {
hw = (TH1*)fWeights.At(par);
Double_t weightedFraction = hw ?
hw->GetBinContent(bin) * fractions[par] : fractions[par];
if (weightedFraction == 0) {
Error("FindPrediction","Fraction[%d] = 0!",par);
return;
}
}
if (TMath::Nint(fData->GetBinContent(bin)) == 0) {
ti = 1;
k0 = -1;
aki=0;
return;
}
k0 = 0;
TH1 *hw0 = (TH1*)fWeights.At(0);
Double_t refWeightedFraction = hw0 ?
hw0->GetBinContent(bin) * fractions[0] : fractions[0];
for (par = 1; par < fNpar; ++par) {
hw = (TH1*)fWeights.At(par);
Double_t weightedFraction = hw ?
hw->GetBinContent(bin) * fractions[par] : fractions[par];
if (weightedFraction > refWeightedFraction) {
k0 = par;
refWeightedFraction = weightedFraction;
}
}
int nMax = 1;
Double_t contentsMax = ((TH1*)fMCs.At(k0))->GetBinContent(bin);
for (par = 0; par < fNpar; ++par) {
if (par == k0) continue;
hw = (TH1*)fWeights.At(par);
Double_t weightedFraction = hw ?
hw->GetBinContent(bin) * fractions[par] : fractions[par];
if (weightedFraction == refWeightedFraction) {
nMax++;
contentsMax += ((TH1*)fMCs.At(par))->GetBinContent(bin);
}
}
Double_t tmin = -1/refWeightedFraction;
if (contentsMax == 0) {
aki = fData->GetBinContent(bin)/(1+refWeightedFraction);
for (par = 0; par < fNpar; ++par) {
hw = (TH1*)fWeights.At(par);
if (par != k0) {
Double_t weightedFraction = hw ?
hw->GetBinContent(bin) * fractions[par] : fractions[par];
if (weightedFraction != refWeightedFraction)
aki -= ((TH1*)fMCs.At(par))->GetBinContent(bin)*weightedFraction/ (refWeightedFraction - weightedFraction);
}
}
if (aki > 0) {
if (nMax) aki /= nMax;
ti = tmin;
return;
}
}
k0 = -1;
ti = 0;
for (Double_t step = 0.2;;) {
if (ti >= 1 || ti < tmin) {
step /= 10;
ti = 0;
}
Double_t aFunction = - fData->GetBinContent(bin)/(1-ti);
Double_t aDerivative = aFunction / (1-ti);
for (par = 0; par < fNpar; ++par) {
TH1 *h = (TH1*)fMCs.At(par);
hw = (TH1*)fWeights.At(par);
Double_t weightedFraction = hw ?
hw->GetBinContent(bin) * fractions[par] : fractions[par];
Double_t d = 1/(ti+1/weightedFraction);
aFunction += h->GetBinContent(bin)*d;
aDerivative -= h->GetBinContent(bin)*d*d;
}
if (TMath::Abs(aFunction) < 1e-12) break;
Double_t delta = -aFunction/aDerivative;
if (TMath::Abs(delta) > step)
delta = (delta > 0) ? step : -step;
ti += delta;
if (TMath::Abs(delta) < 1e-13) break;
}
return;
}
void TFractionFitFCN(Int_t& npar, Double_t* gin, Double_t& f, Double_t* par, Int_t flag) {
TFractionFitter* fitter = dynamic_cast<TFractionFitter*>(fractionFitter->GetObjectFit());
if (!fitter) {
fitter->Error("TFractionFitFCN","Invalid fit object encountered!");
return;
}
fitter->ComputeFCN(npar, gin, f, par, flag);
}
Double_t TFractionFitter::GetChisquare() const
{
return fChisquare;
}
Int_t TFractionFitter::GetNDF() const
{
if (fNDF == 0) return fNpfits-fNpar;
return fNDF;
}
Double_t TFractionFitter::GetProb() const
{
Int_t ndf = fNpfits - fNpar;
if (ndf <= 0) return 0;
return TMath::Prob(fChisquare,ndf);
}
void TFractionFitter::ComputeChisquareLambda()
{
if ( !fFitDone ) {
Error("ComputeChisquareLambda","Fit not yet (successfully) performed");
fChisquare = 0;
return;
}
if (! fPlot)
GetPlot();
Int_t minX, maxX, minY, maxY, minZ, maxZ;
GetRanges(minX, maxX, minY, maxY, minZ, maxZ);
Double_t logLyn = 0;
Double_t logLmn = 0;
for(Int_t x = minX; x <= maxX; x++) {
for(Int_t y = minY; y <= maxY; y++) {
for(Int_t z = minZ; z <= maxZ; z++) {
Double_t di = fData->GetBinContent(x, y, z);
Double_t fi = fPlot->GetBinContent(x, y, z);
if(fi != 0) logLyn += di * TMath::Log(fi) - fi;
if(di != 0) logLmn += di * TMath::Log(di) - di;
for(Int_t j = 0; j < fNpar; j++) {
Double_t aji = ((TH1*)fMCs.At(j))->GetBinContent(x, y, z);
Double_t bji = ((TH1*)fAji.At(j))->GetBinContent(x, y, z);
if(bji != 0) logLyn += aji * TMath::Log(bji) - bji;
if(aji != 0) logLmn += aji * TMath::Log(aji) - aji;
}
}
}
}
fChisquare = -2*logLyn + 2*logLmn;
return;
}
TH1* TFractionFitter::GetMCPrediction(Int_t parm) const
{
CheckParNo(parm);
if ( !fFitDone ) {
Error("GetMCPrediction","Fit not yet performed");
return 0;
}
return (TH1*) fAji.At(parm);
}
Last change: Wed Jun 25 08:39:17 2008
Last generated: 2008-06-25 08:39
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.