////////////////////////////////////////////////////////////
// Author: Gero Flucke <mailto:flucke@mail.desy.de>
//____________________________________
// GFFitHistFrac
//   Author:      Gero Flucke
//   Date:        Sept 19th, 2002
//   last update: $Date: 2002/10/29 14:59:09 $  
//   by:          $Author: flucke $
//
// This class performs fits to get the best fractions of different 
// MC-contributions.
//
// Limitations: Error propagation is not 100% valid in case of more than
//              two MC-contributions
//
////////////////////////////////////////////////////////////


#include "TH1.h"
#include "TF1.h"
#include "TObjArray.h"

#include "GFFitHistFrac.h"

ClassImp(GFFitHistFrac)

using namespace std;

TObjArray* GFFitHistFrac::fgMcHists = NULL; // must contain TH1*, with same binning!
 Double_t GFFitHistFrac::FitFunc1D(Double_t* x, Double_t* par)
{
  // p[0] factor in front of fgMcHists->At(1),
  // p[1] factor in front of fgMcHists->At(2) etc.
  // (1. - p[1] - p[2] - ...) factor in front of fgMcHists->At(0)

  TH1* firstHist = static_cast<TH1*>(fgMcHists->First());
  Int_t bin = firstHist->GetXaxis()->FindBin(x[0]);
  Int_t numHists = fgMcHists->GetEntriesFast();
  Double_t result = 0.;
  Double_t fractionFirst = 1.;
  
  for(Int_t i = 1; i < numHists; ++i){ // 0'th hist (=firstHist) to be treated specially!
    result += 
      par[i-1] * static_cast<TH1*>(fgMcHists->At(i))->GetBinContent(bin);
    fractionFirst -= par[i-1];
  }
  result += fractionFirst * firstHist->GetBinContent(bin);

  return result;
}

//____________________________________
 GFFitHistFrac::GFFitHistFrac() : 
  fHistDataNorm(NULL), fHistsMcNorm(NULL), 
  fScalesMc(0), fFractions(0), fFractionErrors(0)
{

}

//____________________________________
 GFFitHistFrac::GFFitHistFrac(const TH1* histData, const TObjArray* histsMc) : 
  fHistDataNorm(NULL), fHistsMcNorm(NULL), 
  fScalesMc(0), fFractions(0), fFractionErrors(0)
{
  this->SetHistDataMc(histData, histsMc);
}

//____________________________________
 GFFitHistFrac::~GFFitHistFrac()
{
  delete fHistDataNorm;
  if(fHistsMcNorm) fHistsMcNorm->Delete();
  delete fHistsMcNorm;
}

//____________________________________
 Bool_t GFFitHistFrac::SetHistData(const TH1* hist)
{
  if(fHistDataNorm) {
    delete fHistDataNorm;
    fScaleData = 1.;
  }
  if(hist){
    fHistDataNorm = static_cast<TH1*>(hist->Clone());
    fScaleData = fHistDataNorm->Integral();//"width"); ??   //FIXME
    if(fScaleData) {
      fHistDataNorm->Scale(1./fScaleData);
      fHistDataNorm->SetMaximum(fHistDataNorm->GetMaximum()/fScaleData);
    }
  }
  return (Bool_t) hist;
}

//____________________________________
 Bool_t GFFitHistFrac::SetHistsMc(const TObjArray* hists)
{
  if(fHistsMcNorm) {
    fHistsMcNorm->Delete();
  }
  else fHistsMcNorm = new TObjArray;

  Int_t numMc = hists->GetEntries();
  fScalesMc.Set(numMc);
  TIter next(hists);
  Int_t count = 0;
  while(TObject* hOrig = next()){
    TH1* h = static_cast<TH1*>(hOrig->Clone());
    fScalesMc[count] = h->Integral();//"width"); ??   //FIXME: 
    if(fScalesMc[count]) {
      h->Scale(1./fScalesMc[count]);
      h->SetMaximum(h->GetMaximum()/fScalesMc[count]);
    }
    fHistsMcNorm->Add(h);
    ++count;
  }
  return (Bool_t) hists;
}

//____________________________________
 Bool_t GFFitHistFrac::SetHistDataMc(const TH1* histData, const TObjArray* histsMc)
{
  return(this->SetHistData(histData) && this->SetHistsMc(histsMc));
}

//____________________________________
 Bool_t GFFitHistFrac::CheckConsistency()
{
  // checks whether we have data and MC histograms 
  // and whether they have the same binning
 
  // check NULL pointers:
  if(!fHistDataNorm || !fHistsMcNorm){
    this->Error("CheckConsistency", "NULL pointer: dataHist = 0X%X, mcHists = 0X%X",
		fHistDataNorm, fHistsMcNorm);
    return kFALSE;
  }

  // enough entries?
  Int_t numMc = fHistsMcNorm->GetEntriesFast();
  if(numMc < 2) {
    this->Error("CheckConsistency", "Only %d histograms to fit fractions!", numMc);
    return kFALSE;
  }
  
  TIter histMcIter(fHistsMcNorm);
  while(TH1* mc = static_cast<TH1*>(histMcIter.Next())){
    for(Int_t i = 1; i < fHistDataNorm->GetNbinsX()+1; ++i){ // +1: test also upper edge!
      if(fHistDataNorm->GetBinLowEdge(i) != mc->GetBinLowEdge(i)){
	this->Error("CheckConsistency", "bin %i has no common lower edge!", i);
	return kFALSE;
      }
    }
  }
  return kTRUE;
}

//____________________________________
 TArrayD GFFitHistFrac::Fit()
{
  // performs the fit, returns an array with a fraction for each MC-hist
  // each fraction shows how much the MC hist contributes to the shape(!)
  // of the data histogram 

  // clear results:
  fFractions.Reset();
  fFractionErrors.Reset();

  // fit if hists are consistent:
  if(this->CheckConsistency()){
    // initialise fractions to be equal (summing up to 1.) 
    // and errors to be zero:
    Int_t numFractions = fHistsMcNorm->GetEntriesFast();
    fFractions.Set(numFractions);
    fFractions.Reset(1./(numFractions));
    fFractionErrors.Set(numFractions); // new places in the array are set to 0!
    
    fgMcHists = fHistsMcNorm; // global/static variable used in fitting func.!!!
    TF1 mcHistsFunc("mcHistsFunc", GFFitHistFrac::FitFunc1D,
		    fHistDataNorm->GetXaxis()->GetXmin(),
		    fHistDataNorm->GetXaxis()->GetXmax(),
		    numFractions-1);
    // initialise with equal 'fractions','+1' is tricky: 
    mcHistsFunc.SetParameters(fFractions.GetArray()+1); 
    fHistDataNorm->Fit(&mcHistsFunc, "EMRN");//"0+"
    // filling the output arrays:
    fFractions[0] = 1.;
    Double_t firstErrorSqr = 0;
    for(Int_t i = 1; i < numFractions; ++i){ // starting at '1'!
      fFractions[i] = mcHistsFunc.GetParameter(i-1);
      fFractions[0] -= fFractions[i];
      fFractionErrors[i] = mcHistsFunc.GetParError(i-1);
      // FIXME: assuming uncorrelated fitting parameters:
      // (OK for just 2 MC hists: there is only one parameter!)
      firstErrorSqr += fFractionErrors[i] * fFractionErrors[i];         
    }
    fFractionErrors[0] = TMath::Sqrt(firstErrorSqr);// not 100% correct, cf. up
    
//     fgMcHists = NULL; leads to crash: needed by ROOT to draw fitting function
  }

  return fFractions;
}

//____________________________________
 TArrayD GFFitHistFrac::Fit(const TH1* dataHist, const TObjArray* mcHists)
{
  // setting histograms and then Fit()
  this->SetHistDataMc(dataHist, mcHists);
  return this->Fit();
}

//____________________________________
 TArrayD GFFitHistFrac::GetFactors() const
{
  // after Fit() returning array of factors each MC-histogram has to be 
  // multiplied with to get the best agreement with data
  TArrayD factors(fFractions.GetSize());
  for(Int_t i = 0; i < factors.GetSize(); ++i){
    factors[i] = fScalesMc.At(i) ? 
      fFractions.At(i) * fScaleData / fScalesMc.At(i) : 0.;
  }
  return factors;
}

//____________________________________
 TArrayD GFFitHistFrac::GetFactorErrors() const
{
  // returns array of errors belonging to array from GetFractions()
  // FIXME: assuming fractions of the others were independend!! ???
  // (is fine in case of just fitting 2 MC hists to data!)

  TArrayD result(fFractionErrors.GetSize());
  this->Warning("GetFractionErrors", "Not yet implemented");
  return result;
}