// @(#)root/hist:$Id$
// Author: Frank Filthaut F.Filthaut@science.ru.nl  20/05/2002
// with additions by Bram Wijngaarden <dwijngaa@hef.kun.nl>

///////////////////////////////////////////////////////////////////////////////
//
// Fits MC fractions to data histogram (a la HMCMLL, see R. Barlow and C. Beeston,
// Comp. Phys. Comm. 77 (1993) 219-228, and http://www.hep.man.ac.uk/~roger/hfrac.f).
//
// The virtue of this fit is that it takes into account both data and Monte Carlo
// statistical uncertainties. The way in which this is done is through a standard
// likelihood fit using Poisson statistics; however, the template (MC) predictions
// are also varied within statistics, leading to additional contributions to the
// overall likelihood. This leads to many more fit parameters (one per bin per
// template), but the minimisation with respect to these additional parameters is
// done analytically rather than introducing them as formal fit parameters. Some
// special care needs to be taken in the case of bins with zero content. For more
// details please see the original publication cited above.
//
// An example application of this fit is given below. For a TH1* histogram
// ("data") fitted as the sum of three Monte Carlo sources ("mc"):
//
// {
//   TH1F *data;                              //data histogram
//   TH1F *mc0;                               // first MC histogram
//   TH1F *mc1;                               // second MC histogram
//   TH1F *mc2;                               // third MC histogram
//   ....                                     // retrieve histograms
//   TObjArray *mc = new TObjArray(3);        // MC histograms are put in this array
//   TFractionFitter* fit = new TFractionFitter(data, mc); // initialise
//   fit->Constrain(1,0.0,1.0);               // constrain fraction 1 to be between 0 and 1
//   fit->SetRangeX(1,15);                    // use only the first 15 bins in the fit
//   Int_t status = fit->Fit();               // perform the fit
//   std::cout << "fit status: " << status << std::endl;
//   if (status == 0) {                       // check on fit status
//     TH1F* result = (TH1F*) fit->GetPlot();
//     data->Draw("Ep");
//     result->Draw("same");
//   }
// }
//
//
// Assumptions
// ===========
// A few assumptions need to be made for the fit procedure to be carried out:
//
// (1) The total number of events in each template is not too small
//     (so that its Poisson uncertainty can be neglected).
// (2) The number of events in each bin is much smaller than the total
//     number of events in each template (so that multinomial
//     uncertainties can be replaced with Poisson uncertainties).
//
// Biased fit uncertainties may result if these conditions are not fulfilled
// (see e.g. arXiv:0803.2711).
//
// Instantiation
// =============
// A fit object is instantiated through
//     TFractionFitter* fit = new TFractionFitter(data, mc);
// A number of basic checks (intended to ensure that the template histograms
// represent the same "kind" of distribution as the data one) are carried out.
// The TVirtualFitter object is then addressed and all fit parameters (the
// template fractions) declared (initially unbounded).
//
// Applying constraints
// ====================
// Fit parameters can be constrained through
//     fit->Constrain(parameter #, lower bound, upper bound);
// Setting lower bound = upper bound = 0 removes the constraint (a la Minuit);
// however, a function
//     fit->Unconstrain(parameter #)
// is also provided to simplify this.
//
// Setting parameter values
// ========================
// The function
//     TVirtualFitter* vFit = fit->GetFitter();
// is provided for direct access to the TVirtualFitter object. This allows to
// set and fix parameter values, and set step sizes directly.
//
// Restricting the fit range
// =========================
// The fit range can be restricted through
//     fit->SetRangeX(first bin #, last bin #);
// and freed using
//     fit->ReleaseRangeX();
// For 2D histograms the Y range can be similarly restricted using
//     fit->SetRangeY(first bin #, last bin #);
//     fit->ReleaseRangeY();
// and for 3D histograms also
//     fit->SetRangeZ(first bin #, last bin #);
//     fit->ReleaseRangeZ();
// It is also possible to exclude individual bins from the fit through
//     fit->ExcludeBin(bin #);
// where the given bin number is assumed to follow the TH1::GetBin() numbering.
// Any bins excluded in this way can be included again using the corresponding
//     fit->IncludeBin(bin #);
//
// Weights histograms
// ==================
// Weights histograms (for a motivation see the above publication) can be specified
// for the individual MC sources through
//     fit->SetWeight(parameter #, pointer to weights histogram);
// and unset by specifying a null pointer.
//
// Obtaining fit results
// =====================
// The fit is carried out through
//     Int_t status = fit->Fit();
// where  status  is the code returned from the "MINIMIZE" command. For fits
// that converged, parameter values and errors can be obtained through
//     fit->GetResult(parameter #, value, error);
// and the histogram corresponding to the total Monte Carlo prediction (which
// is not the same as a simple weighted sum of the input Monte Carlo distributions)
// can be obtained by
//     TH1* result = fit->GetPlot();
//
// Using different histograms
// ==========================
// It is possible to change the histogram being fitted through
//     fit->SetData(TH1* data);
// and to change the template histogram for a given parameter number through
//     fit->SetMC(parameter #, TH1* MC);
// This can speed up code in case of multiple data or template histograms;
// however, it should be done with care as any settings are taken over from
// the previous fit. In addition, neither the dimensionality nor the numbers of
// bins of the histograms should change (in that case it is better to instantiate
// a new TFractionFitter object).
//
// Errors
// ======
// Any serious inconsistency results in an error.
//
///////////////////////////////////////////////////////////////////////////////

#include "Riostream.h"
#include "TH1.h"
#include "TMath.h"
#include "TClass.h"

#include "Fit/Fitter.h"
#include "TFitResult.h"
#include "Math/Functor.h"
#include "TFractionFitter.h"

ClassImp(TFractionFitter)

//______________________________________________________________________________
TFractionFitter::TFractionFitter() :
fFitDone(kFALSE),
fLowLimitX(0), fHighLimitX(0),
fLowLimitY(0), fHighLimitY(0),
fLowLimitZ(0), fHighLimitZ(0),
fData(0), fIntegralData(0),
fPlot(0)
{
// TFractionFitter default constructor.

fFractionFitter = 0;
fIntegralMCs   = 0;
fFractions     = 0;

fNpfits        = 0;
fNDF           = 0;
fChisquare     = 0;
fNpar          = 0;
}

//______________________________________________________________________________
TFractionFitter::TFractionFitter(TH1* data, TObjArray  *MCs, Option_t *option) :
fFitDone(kFALSE), fChisquare(0), fPlot(0)  {
// TFractionFitter constructor. Does a complete initialisation (including
// consistency checks, default fit range as the whole histogram but without
// under- and overflows, and declaration of the fit parameters). Note that
// the histograms are not copied, only references are used.
// Arguments:
//     data: histogram to be fitted
//     MCs:  array of TH1* corresponding template distributions
//    Option:  can be used to control the print level of the minimization algorithm
//             option = "Q"  : quite - no message is printed
//             option = "V"  : verbose - max print out
//             option = ""   : default: print initial fraction values and result

fData = data;
// Default: include all of the histogram (but without under- and overflows)
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) {
// Histogram containing template prediction
TString s = Form("Prediction for MC sample %i",par);
TH1* pred = (TH1*) ((TH1*)MCs->At(par))->Clone(s);
pred->SetTitle(s);
}
fIntegralMCs = new Double_t[fNpar];
fFractions = new Double_t[fNpar];

CheckConsistency();
fWeights.Expand(fNpar);

fFractionFitter = new ROOT::Fit::Fitter();

// set print level
TString opt(option);
opt.ToUpper();
if (opt.Contains("Q") ) {
fFractionFitter->Config().MinimizerOptions().SetPrintLevel(0);
}
else if (opt.Contains("V") ) {
fFractionFitter->Config().MinimizerOptions().SetPrintLevel(2);
}
else
fFractionFitter->Config().MinimizerOptions().SetPrintLevel(1);

Double_t defaultFraction = 1.0/((Double_t)fNpar);
Double_t defaultStep = 0.01;
// set the parameters
std::vector<ROOT::Fit::ParameterSettings> & parameters = fFractionFitter->Config().ParamsSettings();
parameters.reserve(fNpar);
for (par = 0; par < fNpar; ++par) {
TString name("frac"); name += par;
parameters.push_back(ROOT::Fit::ParameterSettings(name.Data(), defaultFraction, defaultStep) );
}

if (fFractionFitter->Config().MinimizerOptions().ErrorDef() == 1.0 )
fFractionFitter->Config().MinimizerOptions().SetErrorDef(0.5);

}

//______________________________________________________________________________
TFractionFitter::~TFractionFitter() {
// TFractionFitter default destructor

if (fFractionFitter) delete fFractionFitter;
delete[] fIntegralMCs;
delete[] fFractions;
if (fPlot) delete fPlot;
fAji.Delete();
}

//______________________________________________________________________________
void TFractionFitter::SetData(TH1* data) {
// Change the histogram to be fitted to. Notes:
// - Parameter constraints and settings are retained from a possible previous fit.
// - Modifying the dimension or number of bins results in an error (in this case
//   rather instantiate a new TFractionFitter object)

fData = data;
fFitDone = kFALSE;
CheckConsistency();
}

//______________________________________________________________________________
void TFractionFitter::SetMC(Int_t parm, TH1* MC) {
// Change the histogram for template number <parm>. Notes:
// - Parameter constraints and settings are retained from a possible previous fit.
// - Modifying the dimension or number of bins results in an error (in this case
//   rather instantiate a new TFractionFitter object)

CheckParNo(parm);
fMCs.RemoveAt(parm);
fFitDone = kFALSE;
CheckConsistency();
}

//______________________________________________________________________________
void TFractionFitter::SetWeight(Int_t parm, TH1* weight) {
// Set bin by bin weights for template number <parm> (the parameter numbering
// follows that of the input template vector).
// Weights can be "unset" by passing a null pointer.
// Consistency of the weights histogram with the data histogram is checked at
// this point, and an error in case of problems.

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();
}
}

//______________________________________________________________________________
ROOT::Fit::Fitter* TFractionFitter::GetFitter() const {
// Give direct access to the underlying fitter class. This can be
// used e.g. to modify parameter values or step sizes.

return fFractionFitter;
}

//______________________________________________________________________________
void TFractionFitter::CheckParNo(Int_t parm) const {
// Function for internal use, checking parameter validity
// An invalid parameter results in an error.

if (parm < 0 || parm > fNpar) {
Error("CheckParNo","Invalid parameter number %d",parm);
}
}

//______________________________________________________________________________
void TFractionFitter::SetRangeX(Int_t low, Int_t high) {
// Set the X range of the histogram to be used in the fit.
// Use ReleaseRangeX() to go back to fitting the full histogram.
// The consistency check ensures that no empty fit range occurs (and also
// recomputes the bin content integrals).
// Arguments:
//     low:  lower X bin number
//     high: upper X bin number

fLowLimitX = (low > 0 ) ? low : 1;
fHighLimitX = ( high > 0 && high <= fData->GetNbinsX()) ? high : fData->GetNbinsX();
CheckConsistency();
}

//______________________________________________________________________________
void TFractionFitter::ReleaseRangeX() {
// Release restrictions on the X range of the histogram to be used in the fit.

fLowLimitX = 1;
fHighLimitX = fData->GetNbinsX();
CheckConsistency();
}

//______________________________________________________________________________
void TFractionFitter::SetRangeY(Int_t low, Int_t high) {
// Set the Y range of the histogram to be used in the fit (2D or 3D histograms only).
// Use ReleaseRangeY() to go back to fitting the full histogram.
// The consistency check ensures that no empty fit range occurs (and also
// recomputes the bin content integrals).
// Arguments:
//     low:  lower Y bin number
//     high: upper Y bin number

if (fData->GetDimension() < 2) {
Error("SetRangeY","Y range cannot be set for 1D histogram");
return;
}

fLowLimitY = (low > 0) ? low : 1;
fHighLimitY = (high > 0 && high <= fData->GetNbinsY()) ? high : fData->GetNbinsY();
CheckConsistency();
}

//______________________________________________________________________________
void TFractionFitter::ReleaseRangeY() {
// Release restrictions on the Y range of the histogram to be used in the fit.

fLowLimitY = 1;
fHighLimitY = fData->GetNbinsY();
CheckConsistency();
}

//______________________________________________________________________________
void TFractionFitter::SetRangeZ(Int_t low, Int_t high) {
// Set the Z range of the histogram to be used in the fit (3D histograms only).
// Use ReleaseRangeY() to go back to fitting the full histogram.
// The consistency check ensures that no empty fit range occurs (and also
// recomputes the bin content integrals).
// Arguments:
//     low:  lower Y bin number
//     high: upper Y bin number

if (fData->GetDimension() < 3) {
Error("SetRangeZ","Z range cannot be set for 1D or 2D histogram");
return;
}

fLowLimitZ = (low > 0) ? low : 1;
fHighLimitZ = (high > 0 && high <= fData->GetNbinsZ()) ? high : fData->GetNbinsZ();
CheckConsistency();
}

//______________________________________________________________________________
void TFractionFitter::ReleaseRangeZ() {
// Release restrictions on the Z range of the histogram to be used in the fit.

fLowLimitZ = 1;
fHighLimitZ = fData->GetNbinsZ();
CheckConsistency();
}

//______________________________________________________________________________
void TFractionFitter::ExcludeBin(Int_t bin) {
// Exclude the given bin from the fit. The bin numbering to be used is that
// of TH1::GetBin().

int excluded = fExcludedBins.size();
for (int b = 0; b < excluded; ++b) {
if (fExcludedBins[b] == bin) {
Error("ExcludeBin", "bin %d already excluded", bin);
return;
}
}
fExcludedBins.push_back(bin);
// This call serves to properly (re)determine the number of degrees of freeom
CheckConsistency();
}

//______________________________________________________________________________
void TFractionFitter::IncludeBin(Int_t bin) {
// Include the given bin in the fit, if it was excluded before using ExcludeBin().
// The bin numbering to be used is that of TH1::GetBin().

for (std::vector<Int_t>::iterator it = fExcludedBins.begin();
it != fExcludedBins.end(); ++it) {
if (*it == bin) {
fExcludedBins.erase(it);
// This call serves to properly (re)determine the number of degrees of freeom
CheckConsistency();
return;
}
}
Error("IncludeBin", "bin %d was not excluded", bin);
}

//______________________________________________________________________________
bool TFractionFitter::IsExcluded(Int_t bin) const {
// Function for internal use, checking whether the given bin is
// excluded from the fit or not.

for (unsigned int b = 0; b < fExcludedBins.size(); ++b)
if (fExcludedBins[b] == bin) return true;
return false;
}

//______________________________________________________________________________
void TFractionFitter::Constrain(Int_t parm, Double_t low, Double_t high) {
// Constrain the values of parameter number <parm> (the parameter numbering
// follows that of the input template vector).
// Use UnConstrain() to remove this constraint.

CheckParNo(parm);
assert( parm >= 0 && parm <    (int) fFractionFitter->Config().ParamsSettings().size() );
fFractionFitter->Config().ParSettings(parm).SetLimits(low,high);
}

//______________________________________________________________________________
void TFractionFitter::UnConstrain(Int_t parm) {
// Remove the constraints on the possible values of parameter <parm>.

CheckParNo(parm);
fFractionFitter->Config().ParSettings(parm).RemoveLimits();
}

//______________________________________________________________________________
void TFractionFitter::CheckConsistency() {
// Function used internally to check the consistency between the
// various histograms. Checks are performed on nonexistent or empty
// histograms, the precise histogram class, and the number of bins.
// In addition, integrals over the "allowed" bin ranges are computed.
// Any inconsistency results in a error.

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) {
if (IsExcluded(fData->GetBin(x, y, z))) continue;
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) {
Int_t bin = fData->GetBin(x, y, z);
if (IsExcluded(bin)) continue;
Double_t MCEvents = h->GetBinContent(bin);
if (MCEvents < 0) {
Error("CheckConsistency", "Number of MC events (bin = %d, par = %d) cannot be negative: "
" their distribution is binomial (see paper)", bin, par);
}
fIntegralMCs[par] += MCEvents;
}
}
}
if (fIntegralMCs[par] <= 0) {
Error("CheckConsistency","Empty MC histogram #%d",par);
}
}
}

//______________________________________________________________________________
TFitResultPtr TFractionFitter::Fit() {
// Perform the fit with the default UP value.
// The value returned is the minimisation status.

// remove any existing output histogram
if (fPlot) {
delete fPlot; fPlot = 0;
}

// Make sure the correct likelihood computation is used
ROOT::Math::Functor fcnFunction(this, &TFractionFitter::EvaluateFCN, fNpar);
fFractionFitter->SetFCN(static_cast<ROOT::Math::IMultiGenFunction&>(fcnFunction));

// fit
Bool_t status = fFractionFitter->FitFCN();
if (!status)   Warning("Fit","Abnormal termination of minimization.");

fFitDone = kTRUE;

// determine goodness of fit
ComputeChisquareLambda();

// create a new result class
TFitResult* fr = new TFitResult(fFractionFitter->Result());
TString name = TString::Format("TFractionFitter_result_of_%s",fData->GetName() );
fr->SetName(name); fr->SetTitle(name);
return TFitResultPtr(fr);
}

//______________________________________________________________________________
void TFractionFitter::ErrorAnalysis(Double_t UP) {
// Set UP to the given value (see class TMinuit), and perform a MINOS minimisation.

if (! fFitDone) {
Error("ErrorAnalysis","Fit not yet performed");
return;
}

Double_t up = UP > 0 ? UP : 0.5;
fFractionFitter->Config().MinimizerOptions().SetErrorDef(up);
Bool_t status = fFractionFitter->CalculateMinosErrors();
if (!status) {
Error("ErrorAnalysis","Error return from MINOS: %d",fFractionFitter->Result().Status());
}
}

//______________________________________________________________________________
void TFractionFitter::GetResult(Int_t parm, Double_t& value, Double_t& error) const {
// Obtain the fit result for parameter <parm> (the parameter numbering
// follows that of the input template vector).

CheckParNo(parm);
if (! fFitDone) {
Error("GetResult","Fit not yet performed");
return;
}
value = fFractionFitter->Result().Parameter(parm);
error = fFractionFitter->Result().Error(parm);
}

//______________________________________________________________________________
TH1* TFractionFitter::GetPlot() {
// Return the "template prediction" corresponding to the fit result (this is not
// the same as the weighted sum of template distributions, as template statistical
// uncertainties are taken into account).
// Note that the name of this histogram will simply be the same as that of the
// "data" histogram, prefixed with the string "Fraction fit to hist: ".
// Note also that the histogram is managed by the TFractionFitter class, so the returned pointer will be invalid if
// the class is deleted

if (! fFitDone) {
Error("GetPlot","Fit not yet performed");
return 0;
}
if (! fPlot) {
Double_t f = 0;
const Double_t * par = fFractionFitter->Result().GetParams();
assert(par);
ComputeFCN(f, par, 3);
}
return fPlot;
}

//______________________________________________________________________________
void TFractionFitter::GetRanges(Int_t& minX, Int_t& maxX, Int_t& minY, Int_t& maxY,
Int_t& minZ, Int_t& maxZ) const {
// Used internally to obtain the bin ranges according to the dimensionality of
// the histogram and the limits set by hand.

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(Double_t& f, const Double_t* xx, Int_t flag)
{
// Used internally to compute the likelihood value.

// normalise the fit parameters
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) {
if (IsExcluded(fData->GetBin(x, y, z))) continue;
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();
}
// likelihood computation
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);
if (IsExcluded(bin)) continue;

// Solve for the "predictions"
int k0 = 0;
Double_t ti = 0.0; Double_t aki = 0.0;
FindPrediction(bin, 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 &t_i, int& k_0, Double_t &A_ki) const {
// Function used internally to obtain the template prediction in the individual bins
// 'bin' <=> 'i' (paper)
// 'par' <=> 'j' (paper)

std::vector<Double_t> wgtFrac(fNpar); // weighted fractions (strengths of the sources)
std::vector<Double_t> a_ji(fNpar); // number of observed MC events for bin 'i' and par (source) 'j'
Double_t d_i = fData->GetBinContent(bin); // number of events in the real data for bin 'i'

// Cache the weighted fractions and the number of observed MC events
// Sanity check: none of the fractions should be == 0
for (Int_t par = 0; par < fNpar; ++par) {
a_ji[par] = ((TH1*)fMCs.At(par))->GetBinContent(bin);
TH1* hw = (TH1*)fWeights.At(par);
wgtFrac[par] = hw ? hw->GetBinContent(bin) * fFractions[par] : fFractions[par];
if (wgtFrac[par] == 0) {
Error("FindPrediction", "Fraction[%d] = 0!", par);
return;
}
}

// Case data = 0
if (TMath::Nint(d_i) == 0) {
t_i = 1;
k_0 = -1;
A_ki = 0;
return;
}

// Case one or more of the MC bin contents == 0 -> find largest fraction
// k_0 stores the source index of the largest fraction
k_0 = 0;
Double_t maxWgtFrac = wgtFrac[0];
for (Int_t par = 1; par < fNpar; ++par) {
if (wgtFrac[par] > maxWgtFrac) {
k_0 = par;
maxWgtFrac = wgtFrac[par];
}
}
Double_t t_min = -1 / maxWgtFrac; // t_i cannot be smaller than this value (see paper, par 5)

// Determine if there are more sources which have the same maximum contribution (fraction)
Int_t nMax = 1; Double_t contentsMax = a_ji[k_0];
for (Int_t par = 0; par < fNpar; ++par) {
if (par == k_0) continue;
if (wgtFrac[par] == maxWgtFrac) {
nMax++;
contentsMax += a_ji[par];
}
}

// special action if there is a zero in the number of entries for the MC source with
// the largest strength (fraction) -> See Paper, Paragraph 5
if (contentsMax == 0) {
A_ki = d_i / (1.0 + maxWgtFrac);
for (Int_t par = 0; par < fNpar; ++par) {
if (par == k_0 || wgtFrac[par] == maxWgtFrac) continue;
A_ki -= a_ji[par] * wgtFrac[par] / (maxWgtFrac - wgtFrac[par]);
}
if (A_ki > 0) {
A_ki /= nMax;
t_i = t_min;
return;
}
}
k_0 = -1;

// Case of nonzero histogram contents: solve for t_i using Newton's method
// The equation that needs to be solved:
//    func(t_i) = \sum\limits_j{\frac{ p_j a_{ji} }{1 + p_j t_i}} - \frac{d_i}{1 - t_i} = 0
t_i = 0; Double_t step = 0.2;
Int_t maxIter = 100000; // maximum number of iterations
for(Int_t i = 0; i < maxIter; ++i) {
if (t_i >= 1 || t_i < t_min) {
step /= 10;
t_i = 0;
}
Double_t func   = - d_i / (1.0 - t_i);
Double_t deriv = func / (1.0 - t_i);
for (Int_t par = 0; par < fNpar; ++par) {
Double_t r = 1.0 / (t_i + 1.0 / wgtFrac[par]);
func   += a_ji[par] * r;
deriv -= a_ji[par] * r * r;
}
if (TMath::Abs(func) < 1e-12) return; // solution found
Double_t delta = - func / deriv; // update delta
if (TMath::Abs(delta) > step)
delta = (delta > 0) ? step : -step; // correct delta if it becomes too large
t_i += delta;
if (TMath::Abs(delta) < 1e-13) return; // solution found
} // the loop breaks when the solution is found, or when the maximum number of iterations is exhausted

Warning("FindPrediction", "Did not find solution for t_i in %d iterations", maxIter);

return;
}

#ifdef OLD
//______________________________________________________________________________
void TFractionFitFCN(Int_t& npar, Double_t* gin, Double_t& f, Double_t* par, Int_t flag) {
// Function called by the minimisation package. The actual functionality is passed
// on to the TFractionFitter::ComputeFCN member function.

TFractionFitter* fitter = dynamic_cast<TFractionFitter*>(fFractionFitter->GetObjectFit());
if (!fitter) {
Error("TFractionFitFCN","Invalid fit object encountered!");
return;
}
fitter->ComputeFCN(npar, gin, f, par, flag);
}
#endif

//______________________________________________________________________________
Double_t TFractionFitter::GetChisquare() const
{
// Return the likelihood ratio Chi-squared (chi2) for the fit.
// The value is computed when the fit is executed successfully.
// Chi2 calculation is based on the "likelihood ratio" lambda,
// lambda = L(y;n) / L(m;n),
// where L(y;n) is the likelihood of the fit result <y> describing the data <n>
// and L(m;n) is the likelihood of an unknown "true" underlying distribution
// <m> describing the data <n>. Since <m> is unknown, the data distribution is
// lambda = L(y;n) / L(n;n).
// Note that this ratio is 1 if the fit is perfect. The chi2 value is then
// computed according to
// chi2 = -2*ln(lambda).
// This parameter can be shown to follow a Chi-square distribution. See for
// example S. Baker and R. Cousins, "Clarification of the use of chi-square
// and likelihood functions in fits to histograms", Nucl. Instr. Meth. A221,
// pp. 437-442 (1984)

return fChisquare;
}

//______________________________________________________________________________
Int_t TFractionFitter::GetNDF() const
{
// return the number of degrees of freedom in the fit
// the fNDF parameter has been previously computed during a fit.
// The number of degrees of freedom corresponds to the number of points
// used in the fit minus the number of templates.

if (fNDF == 0) return fNpfits-fNpar;
return fNDF;
}

//______________________________________________________________________________
Double_t TFractionFitter::GetProb() const
{
// return the fit probability

Int_t ndf = fNpfits - fNpar;
if (ndf <= 0) return 0;
return TMath::Prob(fChisquare,ndf);
}

//______________________________________________________________________________
void TFractionFitter::ComputeChisquareLambda()
{
// Method used internally to compute the likelihood ratio chi2
// See the function GetChisquare() for details

if ( !fFitDone ) {
Error("ComputeChisquareLambda","Fit not yet (successfully) performed");
fChisquare = 0;
return;
}

// fPlot must be initialized and filled. Leave this to the GetPlot() method.
if (! fPlot)
GetPlot();

Int_t minX, maxX, minY, maxY, minZ, maxZ;
GetRanges(minX, maxX, minY, maxY, minZ, maxZ);

Double_t logLyn = 0; // likelihood of prediction
Double_t logLmn = 0; // likelihood of data ("true" distribution)
for(Int_t x = minX; x <= maxX; x++) {
for(Int_t y = minY; y <= maxY; y++) {
for(Int_t z = minZ; z <= maxZ; z++) {
if (IsExcluded(fData->GetBin(x, y, z))) continue;
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
{
// Return the adjusted MC template (Aji) for template (parm).
// Note that the (Aji) times fractions only sum to the total prediction
// of the fit if all weights are 1.
// Note also that the histogram is managed by the TFractionFitter class, so the returned pointer will be invalid if
// the class is deleted

CheckParNo(parm);
if ( !fFitDone ) {
Error("GetMCPrediction","Fit not yet performed");
return 0;
}
return (TH1*) fAji.At(parm);
}

TFractionFitter.cxx:1
TFractionFitter.cxx:2
TFractionFitter.cxx:3
TFractionFitter.cxx:4
TFractionFitter.cxx:5
TFractionFitter.cxx:6
TFractionFitter.cxx:7
TFractionFitter.cxx:8
TFractionFitter.cxx:9
TFractionFitter.cxx:10
TFractionFitter.cxx:11
TFractionFitter.cxx:12
TFractionFitter.cxx:13
TFractionFitter.cxx:14
TFractionFitter.cxx:15
TFractionFitter.cxx:16
TFractionFitter.cxx:17
TFractionFitter.cxx:18
TFractionFitter.cxx:19
TFractionFitter.cxx:20
TFractionFitter.cxx:21
TFractionFitter.cxx:22
TFractionFitter.cxx:23
TFractionFitter.cxx:24
TFractionFitter.cxx:25
TFractionFitter.cxx:26
TFractionFitter.cxx:27
TFractionFitter.cxx:28
TFractionFitter.cxx:29
TFractionFitter.cxx:30
TFractionFitter.cxx:31
TFractionFitter.cxx:32
TFractionFitter.cxx:33
TFractionFitter.cxx:34
TFractionFitter.cxx:35
TFractionFitter.cxx:36
TFractionFitter.cxx:37
TFractionFitter.cxx:38
TFractionFitter.cxx:39
TFractionFitter.cxx:40
TFractionFitter.cxx:41
TFractionFitter.cxx:42
TFractionFitter.cxx:43
TFractionFitter.cxx:44
TFractionFitter.cxx:45
TFractionFitter.cxx:46
TFractionFitter.cxx:47
TFractionFitter.cxx:48
TFractionFitter.cxx:49
TFractionFitter.cxx:50
TFractionFitter.cxx:51
TFractionFitter.cxx:52
TFractionFitter.cxx:53
TFractionFitter.cxx:54
TFractionFitter.cxx:55
TFractionFitter.cxx:56
TFractionFitter.cxx:57
TFractionFitter.cxx:58
TFractionFitter.cxx:59
TFractionFitter.cxx:60
TFractionFitter.cxx:61
TFractionFitter.cxx:62
TFractionFitter.cxx:63
TFractionFitter.cxx:64
TFractionFitter.cxx:65
TFractionFitter.cxx:66
TFractionFitter.cxx:67
TFractionFitter.cxx:68
TFractionFitter.cxx:69
TFractionFitter.cxx:70
TFractionFitter.cxx:71
TFractionFitter.cxx:72
TFractionFitter.cxx:73
TFractionFitter.cxx:74
TFractionFitter.cxx:75
TFractionFitter.cxx:76
TFractionFitter.cxx:77
TFractionFitter.cxx:78
TFractionFitter.cxx:79
TFractionFitter.cxx:80
TFractionFitter.cxx:81
TFractionFitter.cxx:82
TFractionFitter.cxx:83
TFractionFitter.cxx:84
TFractionFitter.cxx:85
TFractionFitter.cxx:86
TFractionFitter.cxx:87
TFractionFitter.cxx:88
TFractionFitter.cxx:89
TFractionFitter.cxx:90
TFractionFitter.cxx:91
TFractionFitter.cxx:92
TFractionFitter.cxx:93
TFractionFitter.cxx:94
TFractionFitter.cxx:95
TFractionFitter.cxx:96
TFractionFitter.cxx:97
TFractionFitter.cxx:98
TFractionFitter.cxx:99
TFractionFitter.cxx:100
TFractionFitter.cxx:101
TFractionFitter.cxx:102
TFractionFitter.cxx:103
TFractionFitter.cxx:104
TFractionFitter.cxx:105
TFractionFitter.cxx:106
TFractionFitter.cxx:107
TFractionFitter.cxx:108
TFractionFitter.cxx:109
TFractionFitter.cxx:110
TFractionFitter.cxx:111
TFractionFitter.cxx:112
TFractionFitter.cxx:113
TFractionFitter.cxx:114
TFractionFitter.cxx:115
TFractionFitter.cxx:116
TFractionFitter.cxx:117
TFractionFitter.cxx:118
TFractionFitter.cxx:119
TFractionFitter.cxx:120
TFractionFitter.cxx:121
TFractionFitter.cxx:122
TFractionFitter.cxx:123
TFractionFitter.cxx:124
TFractionFitter.cxx:125
TFractionFitter.cxx:126
TFractionFitter.cxx:127
TFractionFitter.cxx:128
TFractionFitter.cxx:129
TFractionFitter.cxx:130
TFractionFitter.cxx:131
TFractionFitter.cxx:132
TFractionFitter.cxx:133
TFractionFitter.cxx:134
TFractionFitter.cxx:135
TFractionFitter.cxx:136
TFractionFitter.cxx:137
TFractionFitter.cxx:138
TFractionFitter.cxx:139
TFractionFitter.cxx:140
TFractionFitter.cxx:141
TFractionFitter.cxx:142
TFractionFitter.cxx:143
TFractionFitter.cxx:144
TFractionFitter.cxx:145
TFractionFitter.cxx:146
TFractionFitter.cxx:147
TFractionFitter.cxx:148
TFractionFitter.cxx:149
TFractionFitter.cxx:150
TFractionFitter.cxx:151
TFractionFitter.cxx:152
TFractionFitter.cxx:153
TFractionFitter.cxx:154
TFractionFitter.cxx:155
TFractionFitter.cxx:156
TFractionFitter.cxx:157
TFractionFitter.cxx:158
TFractionFitter.cxx:159
TFractionFitter.cxx:160
TFractionFitter.cxx:161
TFractionFitter.cxx:162
TFractionFitter.cxx:163
TFractionFitter.cxx:164
TFractionFitter.cxx:165
TFractionFitter.cxx:166
TFractionFitter.cxx:167
TFractionFitter.cxx:168
TFractionFitter.cxx:169
TFractionFitter.cxx:170
TFractionFitter.cxx:171
TFractionFitter.cxx:172
TFractionFitter.cxx:173
TFractionFitter.cxx:174
TFractionFitter.cxx:175
TFractionFitter.cxx:176
TFractionFitter.cxx:177
TFractionFitter.cxx:178
TFractionFitter.cxx:179
TFractionFitter.cxx:180
TFractionFitter.cxx:181
TFractionFitter.cxx:182
TFractionFitter.cxx:183
TFractionFitter.cxx:184
TFractionFitter.cxx:185
TFractionFitter.cxx:186
TFractionFitter.cxx:187
TFractionFitter.cxx:188
TFractionFitter.cxx:189
TFractionFitter.cxx:190
TFractionFitter.cxx:191
TFractionFitter.cxx:192
TFractionFitter.cxx:193
TFractionFitter.cxx:194
TFractionFitter.cxx:195
TFractionFitter.cxx:196
TFractionFitter.cxx:197
TFractionFitter.cxx:198
TFractionFitter.cxx:199
TFractionFitter.cxx:200
TFractionFitter.cxx:201
TFractionFitter.cxx:202
TFractionFitter.cxx:203
TFractionFitter.cxx:204
TFractionFitter.cxx:205
TFractionFitter.cxx:206
TFractionFitter.cxx:207
TFractionFitter.cxx:208
TFractionFitter.cxx:209
TFractionFitter.cxx:210
TFractionFitter.cxx:211
TFractionFitter.cxx:212
TFractionFitter.cxx:213
TFractionFitter.cxx:214
TFractionFitter.cxx:215
TFractionFitter.cxx:216
TFractionFitter.cxx:217
TFractionFitter.cxx:218
TFractionFitter.cxx:219
TFractionFitter.cxx:220
TFractionFitter.cxx:221
TFractionFitter.cxx:222
TFractionFitter.cxx:223
TFractionFitter.cxx:224
TFractionFitter.cxx:225
TFractionFitter.cxx:226
TFractionFitter.cxx:227
TFractionFitter.cxx:228
TFractionFitter.cxx:229
TFractionFitter.cxx:230
TFractionFitter.cxx:231
TFractionFitter.cxx:232
TFractionFitter.cxx:233
TFractionFitter.cxx:234
TFractionFitter.cxx:235
TFractionFitter.cxx:236
TFractionFitter.cxx:237
TFractionFitter.cxx:238
TFractionFitter.cxx:239
TFractionFitter.cxx:240
TFractionFitter.cxx:241
TFractionFitter.cxx:242
TFractionFitter.cxx:243
TFractionFitter.cxx:244
TFractionFitter.cxx:245
TFractionFitter.cxx:246
TFractionFitter.cxx:247
TFractionFitter.cxx:248
TFractionFitter.cxx:249
TFractionFitter.cxx:250
TFractionFitter.cxx:251
TFractionFitter.cxx:252
TFractionFitter.cxx:253
TFractionFitter.cxx:254
TFractionFitter.cxx:255
TFractionFitter.cxx:256
TFractionFitter.cxx:257
TFractionFitter.cxx:258
TFractionFitter.cxx:259
TFractionFitter.cxx:260
TFractionFitter.cxx:261
TFractionFitter.cxx:262
TFractionFitter.cxx:263
TFractionFitter.cxx:264
TFractionFitter.cxx:265
TFractionFitter.cxx:266
TFractionFitter.cxx:267
TFractionFitter.cxx:268
TFractionFitter.cxx:269
TFractionFitter.cxx:270
TFractionFitter.cxx:271
TFractionFitter.cxx:272
TFractionFitter.cxx:273
TFractionFitter.cxx:274
TFractionFitter.cxx:275
TFractionFitter.cxx:276
TFractionFitter.cxx:277
TFractionFitter.cxx:278
TFractionFitter.cxx:279
TFractionFitter.cxx:280
TFractionFitter.cxx:281
TFractionFitter.cxx:282
TFractionFitter.cxx:283
TFractionFitter.cxx:284
TFractionFitter.cxx:285
TFractionFitter.cxx:286
TFractionFitter.cxx:287
TFractionFitter.cxx:288
TFractionFitter.cxx:289
TFractionFitter.cxx:290
TFractionFitter.cxx:291
TFractionFitter.cxx:292
TFractionFitter.cxx:293
TFractionFitter.cxx:294
TFractionFitter.cxx:295
TFractionFitter.cxx:296
TFractionFitter.cxx:297
TFractionFitter.cxx:298
TFractionFitter.cxx:299
TFractionFitter.cxx:300
TFractionFitter.cxx:301
TFractionFitter.cxx:302
TFractionFitter.cxx:303
TFractionFitter.cxx:304
TFractionFitter.cxx:305
TFractionFitter.cxx:306
TFractionFitter.cxx:307
TFractionFitter.cxx:308
TFractionFitter.cxx:309
TFractionFitter.cxx:310
TFractionFitter.cxx:311
TFractionFitter.cxx:312
TFractionFitter.cxx:313
TFractionFitter.cxx:314
TFractionFitter.cxx:315
TFractionFitter.cxx:316
TFractionFitter.cxx:317
TFractionFitter.cxx:318
TFractionFitter.cxx:319
TFractionFitter.cxx:320
TFractionFitter.cxx:321
TFractionFitter.cxx:322
TFractionFitter.cxx:323
TFractionFitter.cxx:324
TFractionFitter.cxx:325
TFractionFitter.cxx:326
TFractionFitter.cxx:327
TFractionFitter.cxx:328
TFractionFitter.cxx:329
TFractionFitter.cxx:330
TFractionFitter.cxx:331
TFractionFitter.cxx:332
TFractionFitter.cxx:333
TFractionFitter.cxx:334
TFractionFitter.cxx:335
TFractionFitter.cxx:336
TFractionFitter.cxx:337
TFractionFitter.cxx:338
TFractionFitter.cxx:339
TFractionFitter.cxx:340
TFractionFitter.cxx:341
TFractionFitter.cxx:342
TFractionFitter.cxx:343
TFractionFitter.cxx:344
TFractionFitter.cxx:345
TFractionFitter.cxx:346
TFractionFitter.cxx:347
TFractionFitter.cxx:348
TFractionFitter.cxx:349
TFractionFitter.cxx:350
TFractionFitter.cxx:351
TFractionFitter.cxx:352
TFractionFitter.cxx:353
TFractionFitter.cxx:354
TFractionFitter.cxx:355
TFractionFitter.cxx:356
TFractionFitter.cxx:357
TFractionFitter.cxx:358
TFractionFitter.cxx:359
TFractionFitter.cxx:360
TFractionFitter.cxx:361
TFractionFitter.cxx:362
TFractionFitter.cxx:363
TFractionFitter.cxx:364
TFractionFitter.cxx:365
TFractionFitter.cxx:366
TFractionFitter.cxx:367
TFractionFitter.cxx:368
TFractionFitter.cxx:369
TFractionFitter.cxx:370
TFractionFitter.cxx:371
TFractionFitter.cxx:372
TFractionFitter.cxx:373
TFractionFitter.cxx:374
TFractionFitter.cxx:375
TFractionFitter.cxx:376
TFractionFitter.cxx:377
TFractionFitter.cxx:378
TFractionFitter.cxx:379
TFractionFitter.cxx:380
TFractionFitter.cxx:381
TFractionFitter.cxx:382
TFractionFitter.cxx:383
TFractionFitter.cxx:384
TFractionFitter.cxx:385
TFractionFitter.cxx:386
TFractionFitter.cxx:387
TFractionFitter.cxx:388
TFractionFitter.cxx:389
TFractionFitter.cxx:390
TFractionFitter.cxx:391
TFractionFitter.cxx:392
TFractionFitter.cxx:393
TFractionFitter.cxx:394
TFractionFitter.cxx:395
TFractionFitter.cxx:396
TFractionFitter.cxx:397
TFractionFitter.cxx:398
TFractionFitter.cxx:399
TFractionFitter.cxx:400
TFractionFitter.cxx:401
TFractionFitter.cxx:402
TFractionFitter.cxx:403
TFractionFitter.cxx:404
TFractionFitter.cxx:405
TFractionFitter.cxx:406
TFractionFitter.cxx:407
TFractionFitter.cxx:408
TFractionFitter.cxx:409
TFractionFitter.cxx:410
TFractionFitter.cxx:411
TFractionFitter.cxx:412
TFractionFitter.cxx:413
TFractionFitter.cxx:414
TFractionFitter.cxx:415
TFractionFitter.cxx:416
TFractionFitter.cxx:417
TFractionFitter.cxx:418
TFractionFitter.cxx:419
TFractionFitter.cxx:420
TFractionFitter.cxx:421
TFractionFitter.cxx:422
TFractionFitter.cxx:423
TFractionFitter.cxx:424
TFractionFitter.cxx:425
TFractionFitter.cxx:426
TFractionFitter.cxx:427
TFractionFitter.cxx:428
TFractionFitter.cxx:429
TFractionFitter.cxx:430
TFractionFitter.cxx:431
TFractionFitter.cxx:432
TFractionFitter.cxx:433
TFractionFitter.cxx:434
TFractionFitter.cxx:435
TFractionFitter.cxx:436
TFractionFitter.cxx:437
TFractionFitter.cxx:438
TFractionFitter.cxx:439
TFractionFitter.cxx:440
TFractionFitter.cxx:441
TFractionFitter.cxx:442
TFractionFitter.cxx:443
TFractionFitter.cxx:444
TFractionFitter.cxx:445
TFractionFitter.cxx:446
TFractionFitter.cxx:447
TFractionFitter.cxx:448
TFractionFitter.cxx:449
TFractionFitter.cxx:450
TFractionFitter.cxx:451
TFractionFitter.cxx:452
TFractionFitter.cxx:453
TFractionFitter.cxx:454
TFractionFitter.cxx:455
TFractionFitter.cxx:456
TFractionFitter.cxx:457
TFractionFitter.cxx:458
TFractionFitter.cxx:459
TFractionFitter.cxx:460
TFractionFitter.cxx:461
TFractionFitter.cxx:462
TFractionFitter.cxx:463
TFractionFitter.cxx:464
TFractionFitter.cxx:465
TFractionFitter.cxx:466
TFractionFitter.cxx:467
TFractionFitter.cxx:468
TFractionFitter.cxx:469
TFractionFitter.cxx:470
TFractionFitter.cxx:471
TFractionFitter.cxx:472
TFractionFitter.cxx:473
TFractionFitter.cxx:474
TFractionFitter.cxx:475
TFractionFitter.cxx:476
TFractionFitter.cxx:477
TFractionFitter.cxx:478
TFractionFitter.cxx:479
TFractionFitter.cxx:480
TFractionFitter.cxx:481
TFractionFitter.cxx:482
TFractionFitter.cxx:483
TFractionFitter.cxx:484
TFractionFitter.cxx:485
TFractionFitter.cxx:486
TFractionFitter.cxx:487
TFractionFitter.cxx:488
TFractionFitter.cxx:489
TFractionFitter.cxx:490
TFractionFitter.cxx:491
TFractionFitter.cxx:492
TFractionFitter.cxx:493
TFractionFitter.cxx:494
TFractionFitter.cxx:495
TFractionFitter.cxx:496
TFractionFitter.cxx:497
TFractionFitter.cxx:498
TFractionFitter.cxx:499
TFractionFitter.cxx:500
TFractionFitter.cxx:501
TFractionFitter.cxx:502
TFractionFitter.cxx:503
TFractionFitter.cxx:504
TFractionFitter.cxx:505
TFractionFitter.cxx:506
TFractionFitter.cxx:507
TFractionFitter.cxx:508
TFractionFitter.cxx:509
TFractionFitter.cxx:510
TFractionFitter.cxx:511
TFractionFitter.cxx:512
TFractionFitter.cxx:513
TFractionFitter.cxx:514
TFractionFitter.cxx:515
TFractionFitter.cxx:516
TFractionFitter.cxx:517
TFractionFitter.cxx:518
TFractionFitter.cxx:519
TFractionFitter.cxx:520
TFractionFitter.cxx:521
TFractionFitter.cxx:522
TFractionFitter.cxx:523
TFractionFitter.cxx:524
TFractionFitter.cxx:525
TFractionFitter.cxx:526
TFractionFitter.cxx:527
TFractionFitter.cxx:528
TFractionFitter.cxx:529
TFractionFitter.cxx:530
TFractionFitter.cxx:531
TFractionFitter.cxx:532
TFractionFitter.cxx:533
TFractionFitter.cxx:534
TFractionFitter.cxx:535
TFractionFitter.cxx:536
TFractionFitter.cxx:537
TFractionFitter.cxx:538
TFractionFitter.cxx:539
TFractionFitter.cxx:540
TFractionFitter.cxx:541
TFractionFitter.cxx:542
TFractionFitter.cxx:543
TFractionFitter.cxx:544
TFractionFitter.cxx:545
TFractionFitter.cxx:546
TFractionFitter.cxx:547
TFractionFitter.cxx:548
TFractionFitter.cxx:549
TFractionFitter.cxx:550
TFractionFitter.cxx:551
TFractionFitter.cxx:552
TFractionFitter.cxx:553
TFractionFitter.cxx:554
TFractionFitter.cxx:555
TFractionFitter.cxx:556
TFractionFitter.cxx:557
TFractionFitter.cxx:558
TFractionFitter.cxx:559
TFractionFitter.cxx:560
TFractionFitter.cxx:561
TFractionFitter.cxx:562
TFractionFitter.cxx:563
TFractionFitter.cxx:564
TFractionFitter.cxx:565
TFractionFitter.cxx:566
TFractionFitter.cxx:567
TFractionFitter.cxx:568
TFractionFitter.cxx:569
TFractionFitter.cxx:570
TFractionFitter.cxx:571
TFractionFitter.cxx:572
TFractionFitter.cxx:573
TFractionFitter.cxx:574
TFractionFitter.cxx:575
TFractionFitter.cxx:576
TFractionFitter.cxx:577
TFractionFitter.cxx:578
TFractionFitter.cxx:579
TFractionFitter.cxx:580
TFractionFitter.cxx:581
TFractionFitter.cxx:582
TFractionFitter.cxx:583
TFractionFitter.cxx:584
TFractionFitter.cxx:585
TFractionFitter.cxx:586
TFractionFitter.cxx:587
TFractionFitter.cxx:588
TFractionFitter.cxx:589
TFractionFitter.cxx:590
TFractionFitter.cxx:591
TFractionFitter.cxx:592
TFractionFitter.cxx:593
TFractionFitter.cxx:594
TFractionFitter.cxx:595
TFractionFitter.cxx:596
TFractionFitter.cxx:597
TFractionFitter.cxx:598
TFractionFitter.cxx:599
TFractionFitter.cxx:600
TFractionFitter.cxx:601
TFractionFitter.cxx:602
TFractionFitter.cxx:603
TFractionFitter.cxx:604
TFractionFitter.cxx:605
TFractionFitter.cxx:606
TFractionFitter.cxx:607
TFractionFitter.cxx:608
TFractionFitter.cxx:609
TFractionFitter.cxx:610
TFractionFitter.cxx:611
TFractionFitter.cxx:612
TFractionFitter.cxx:613
TFractionFitter.cxx:614
TFractionFitter.cxx:615
TFractionFitter.cxx:616
TFractionFitter.cxx:617
TFractionFitter.cxx:618
TFractionFitter.cxx:619
TFractionFitter.cxx:620
TFractionFitter.cxx:621
TFractionFitter.cxx:622
TFractionFitter.cxx:623
TFractionFitter.cxx:624
TFractionFitter.cxx:625
TFractionFitter.cxx:626
TFractionFitter.cxx:627
TFractionFitter.cxx:628
TFractionFitter.cxx:629
TFractionFitter.cxx:630
TFractionFitter.cxx:631
TFractionFitter.cxx:632
TFractionFitter.cxx:633
TFractionFitter.cxx:634
TFractionFitter.cxx:635
TFractionFitter.cxx:636
TFractionFitter.cxx:637
TFractionFitter.cxx:638
TFractionFitter.cxx:639
TFractionFitter.cxx:640
TFractionFitter.cxx:641
TFractionFitter.cxx:642
TFractionFitter.cxx:643
TFractionFitter.cxx:644
TFractionFitter.cxx:645
TFractionFitter.cxx:646
TFractionFitter.cxx:647
TFractionFitter.cxx:648
TFractionFitter.cxx:649
TFractionFitter.cxx:650
TFractionFitter.cxx:651
TFractionFitter.cxx:652
TFractionFitter.cxx:653
TFractionFitter.cxx:654
TFractionFitter.cxx:655
TFractionFitter.cxx:656
TFractionFitter.cxx:657
TFractionFitter.cxx:658
TFractionFitter.cxx:659
TFractionFitter.cxx:660
TFractionFitter.cxx:661
TFractionFitter.cxx:662
TFractionFitter.cxx:663
TFractionFitter.cxx:664
TFractionFitter.cxx:665
TFractionFitter.cxx:666
TFractionFitter.cxx:667
TFractionFitter.cxx:668
TFractionFitter.cxx:669
TFractionFitter.cxx:670
TFractionFitter.cxx:671
TFractionFitter.cxx:672
TFractionFitter.cxx:673
TFractionFitter.cxx:674
TFractionFitter.cxx:675
TFractionFitter.cxx:676
TFractionFitter.cxx:677
TFractionFitter.cxx:678
TFractionFitter.cxx:679
TFractionFitter.cxx:680
TFractionFitter.cxx:681
TFractionFitter.cxx:682
TFractionFitter.cxx:683
TFractionFitter.cxx:684
TFractionFitter.cxx:685
TFractionFitter.cxx:686
TFractionFitter.cxx:687
TFractionFitter.cxx:688
TFractionFitter.cxx:689
TFractionFitter.cxx:690
TFractionFitter.cxx:691
TFractionFitter.cxx:692
TFractionFitter.cxx:693
TFractionFitter.cxx:694
TFractionFitter.cxx:695
TFractionFitter.cxx:696
TFractionFitter.cxx:697
TFractionFitter.cxx:698
TFractionFitter.cxx:699
TFractionFitter.cxx:700
TFractionFitter.cxx:701
TFractionFitter.cxx:702
TFractionFitter.cxx:703
TFractionFitter.cxx:704
TFractionFitter.cxx:705
TFractionFitter.cxx:706
TFractionFitter.cxx:707
TFractionFitter.cxx:708
TFractionFitter.cxx:709
TFractionFitter.cxx:710
TFractionFitter.cxx:711
TFractionFitter.cxx:712
TFractionFitter.cxx:713
TFractionFitter.cxx:714
TFractionFitter.cxx:715
TFractionFitter.cxx:716
TFractionFitter.cxx:717
TFractionFitter.cxx:718
TFractionFitter.cxx:719
TFractionFitter.cxx:720
TFractionFitter.cxx:721
TFractionFitter.cxx:722
TFractionFitter.cxx:723
TFractionFitter.cxx:724
TFractionFitter.cxx:725
TFractionFitter.cxx:726
TFractionFitter.cxx:727
TFractionFitter.cxx:728
TFractionFitter.cxx:729
TFractionFitter.cxx:730
TFractionFitter.cxx:731
TFractionFitter.cxx:732
TFractionFitter.cxx:733
TFractionFitter.cxx:734
TFractionFitter.cxx:735
TFractionFitter.cxx:736
TFractionFitter.cxx:737
TFractionFitter.cxx:738
TFractionFitter.cxx:739
TFractionFitter.cxx:740
TFractionFitter.cxx:741
TFractionFitter.cxx:742
TFractionFitter.cxx:743
TFractionFitter.cxx:744
TFractionFitter.cxx:745
TFractionFitter.cxx:746
TFractionFitter.cxx:747
TFractionFitter.cxx:748
TFractionFitter.cxx:749
TFractionFitter.cxx:750
TFractionFitter.cxx:751
TFractionFitter.cxx:752
TFractionFitter.cxx:753
TFractionFitter.cxx:754
TFractionFitter.cxx:755
TFractionFitter.cxx:756
TFractionFitter.cxx:757
TFractionFitter.cxx:758
TFractionFitter.cxx:759
TFractionFitter.cxx:760
TFractionFitter.cxx:761
TFractionFitter.cxx:762
TFractionFitter.cxx:763
TFractionFitter.cxx:764
TFractionFitter.cxx:765
TFractionFitter.cxx:766
TFractionFitter.cxx:767
TFractionFitter.cxx:768
TFractionFitter.cxx:769
TFractionFitter.cxx:770
TFractionFitter.cxx:771
TFractionFitter.cxx:772
TFractionFitter.cxx:773
TFractionFitter.cxx:774
TFractionFitter.cxx:775
TFractionFitter.cxx:776
TFractionFitter.cxx:777
TFractionFitter.cxx:778
TFractionFitter.cxx:779
TFractionFitter.cxx:780
TFractionFitter.cxx:781
TFractionFitter.cxx:782
TFractionFitter.cxx:783
TFractionFitter.cxx:784
TFractionFitter.cxx:785
TFractionFitter.cxx:786
TFractionFitter.cxx:787
TFractionFitter.cxx:788
TFractionFitter.cxx:789
TFractionFitter.cxx:790
TFractionFitter.cxx:791
TFractionFitter.cxx:792
TFractionFitter.cxx:793
TFractionFitter.cxx:794
TFractionFitter.cxx:795
TFractionFitter.cxx:796
TFractionFitter.cxx:797
TFractionFitter.cxx:798
TFractionFitter.cxx:799
TFractionFitter.cxx:800
TFractionFitter.cxx:801
TFractionFitter.cxx:802
TFractionFitter.cxx:803
TFractionFitter.cxx:804
TFractionFitter.cxx:805
TFractionFitter.cxx:806
TFractionFitter.cxx:807
TFractionFitter.cxx:808
TFractionFitter.cxx:809
TFractionFitter.cxx:810
TFractionFitter.cxx:811
TFractionFitter.cxx:812
TFractionFitter.cxx:813
TFractionFitter.cxx:814
TFractionFitter.cxx:815
TFractionFitter.cxx:816
TFractionFitter.cxx:817
TFractionFitter.cxx:818
TFractionFitter.cxx:819
TFractionFitter.cxx:820
TFractionFitter.cxx:821
TFractionFitter.cxx:822
TFractionFitter.cxx:823
TFractionFitter.cxx:824
TFractionFitter.cxx:825
TFractionFitter.cxx:826
TFractionFitter.cxx:827
TFractionFitter.cxx:828
TFractionFitter.cxx:829
TFractionFitter.cxx:830
TFractionFitter.cxx:831
TFractionFitter.cxx:832
TFractionFitter.cxx:833
TFractionFitter.cxx:834
TFractionFitter.cxx:835
TFractionFitter.cxx:836
TFractionFitter.cxx:837
TFractionFitter.cxx:838
TFractionFitter.cxx:839
TFractionFitter.cxx:840
TFractionFitter.cxx:841
TFractionFitter.cxx:842
TFractionFitter.cxx:843
TFractionFitter.cxx:844
TFractionFitter.cxx:845
TFractionFitter.cxx:846
TFractionFitter.cxx:847
TFractionFitter.cxx:848
TFractionFitter.cxx:849
TFractionFitter.cxx:850
TFractionFitter.cxx:851
TFractionFitter.cxx:852
TFractionFitter.cxx:853
TFractionFitter.cxx:854
TFractionFitter.cxx:855
TFractionFitter.cxx:856
TFractionFitter.cxx:857
TFractionFitter.cxx:858
TFractionFitter.cxx:859
TFractionFitter.cxx:860
TFractionFitter.cxx:861
TFractionFitter.cxx:862
TFractionFitter.cxx:863
TFractionFitter.cxx:864
TFractionFitter.cxx:865
TFractionFitter.cxx:866
TFractionFitter.cxx:867
TFractionFitter.cxx:868
TFractionFitter.cxx:869
TFractionFitter.cxx:870
TFractionFitter.cxx:871
TFractionFitter.cxx:872
TFractionFitter.cxx:873
TFractionFitter.cxx:874
TFractionFitter.cxx:875
TFractionFitter.cxx:876
TFractionFitter.cxx:877
TFractionFitter.cxx:878
TFractionFitter.cxx:879
TFractionFitter.cxx:880
TFractionFitter.cxx:881
TFractionFitter.cxx:882
TFractionFitter.cxx:883
TFractionFitter.cxx:884
TFractionFitter.cxx:885
TFractionFitter.cxx:886
TFractionFitter.cxx:887
TFractionFitter.cxx:888
TFractionFitter.cxx:889
TFractionFitter.cxx:890
TFractionFitter.cxx:891
TFractionFitter.cxx:892
TFractionFitter.cxx:893
TFractionFitter.cxx:894
TFractionFitter.cxx:895
TFractionFitter.cxx:896
TFractionFitter.cxx:897
TFractionFitter.cxx:898
TFractionFitter.cxx:899
TFractionFitter.cxx:900
TFractionFitter.cxx:901
TFractionFitter.cxx:902
TFractionFitter.cxx:903
TFractionFitter.cxx:904
TFractionFitter.cxx:905
TFractionFitter.cxx:906
TFractionFitter.cxx:907
TFractionFitter.cxx:908
TFractionFitter.cxx:909
TFractionFitter.cxx:910
TFractionFitter.cxx:911
TFractionFitter.cxx:912
TFractionFitter.cxx:913
TFractionFitter.cxx:914
TFractionFitter.cxx:915
TFractionFitter.cxx:916
TFractionFitter.cxx:917
TFractionFitter.cxx:918
TFractionFitter.cxx:919
TFractionFitter.cxx:920
TFractionFitter.cxx:921
TFractionFitter.cxx:922
TFractionFitter.cxx:923
TFractionFitter.cxx:924
TFractionFitter.cxx:925
TFractionFitter.cxx:926
TFractionFitter.cxx:927
TFractionFitter.cxx:928
TFractionFitter.cxx:929
TFractionFitter.cxx:930
TFractionFitter.cxx:931
TFractionFitter.cxx:932
TFractionFitter.cxx:933
TFractionFitter.cxx:934
TFractionFitter.cxx:935
TFractionFitter.cxx:936
TFractionFitter.cxx:937
TFractionFitter.cxx:938
TFractionFitter.cxx:939
TFractionFitter.cxx:940
TFractionFitter.cxx:941
TFractionFitter.cxx:942
TFractionFitter.cxx:943
TFractionFitter.cxx:944
TFractionFitter.cxx:945
TFractionFitter.cxx:946
TFractionFitter.cxx:947
TFractionFitter.cxx:948
TFractionFitter.cxx:949
TFractionFitter.cxx:950
TFractionFitter.cxx:951
TFractionFitter.cxx:952
TFractionFitter.cxx:953
TFractionFitter.cxx:954
TFractionFitter.cxx:955
TFractionFitter.cxx:956
TFractionFitter.cxx:957
TFractionFitter.cxx:958
TFractionFitter.cxx:959
TFractionFitter.cxx:960
TFractionFitter.cxx:961