// Author: Gero Flucke <mailto:gero.flucke@desy.de>
//______________________________________________
// class for fitting D0 distributions
//
//

#include "D0Fitter.h"

#include <TH1.h>
#include <TH2.h>
#include <TH3.h>
#include <TMath.h>
#include <TDatabasePDG.h>
#include <TROOT.h>
#include <TF1.h>
#include <TError.h>

#include <iostream>

using namespace std;

 Double_t D0Fitter::D0ExpBackgr(Double_t *x, Double_t *par)
{
  // background function:
  // par[0]: normalisation
  // par[1]: decay constant
  // x[0]:   deltaM
  const Double_t norm = par[1] ?
    (TMath::Exp(fgUpFitLimit * par[1]) - TMath::Exp(fgLowFitLimit * par[1])) / par[1] :
    (fgUpFitLimit - fgLowFitLimit);
  return (norm ? (par[0] * TMath::Exp(x[0] * par[1]) * fgBinSize / norm) : 0.);
}

//_____________________________________
 Double_t D0Fitter::D0LinearBackgr(Double_t *x, Double_t *par)
{
  // FIXME: Not used, to be adjusted to give normalisation etc.
  //
  // background function:
  // par[0]: constant term
  // par[1]: linear term
  // x[0]:   deltaM
  return (par[0] + x[0] * par[1]) * fgBinSize;

}

// //********************************************************

 Double_t D0Fitter::D0FitFct(Double_t* x, Double_t* par)
{
  // fit function:
  // par[0]: background normalisation
  // par[1]: background decay constant
  // par[2]: gauss' mean
  // par[3]: gauss' sigma
  // par[4]: N(D0)
  // x[0]:   mass(D0)
  if(par[3] == 0.) {
    ::Warning("D0Fitter::D0FitFct", "width = 0, reject point");
    TF1::RejectPoint();
    return 0.; // no width ==> senseless!
  }

  // fgBinSize is a static member of D0Fitter 
  // and should be determined before every fit (&Integral)
  // to ensure that it has the value of the histogram's
  // binSize (and not that of another histogram that is 
  // fitted by another instance of D0Fitter)!

  const Double_t dx = (x[0]-par[2]) / par[3];
  const Double_t gaus = 1./(D0Fitter::kSqrtOf2Pi*par[3]) * fgBinSize * par[4]
    * TMath::Exp(-dx*dx / 2.);

  return gaus + D0ExpBackgr(x, par);
}

ClassImp(D0Fitter)

const Double_t D0Fitter::kSqrtOf2Pi = TMath::Sqrt(2.*TMath::Pi());
const Double_t D0Fitter::fgBackgrSlopeDefault = -1.2;
const Double_t D0Fitter::fgSigmaDefault = .03;
const Double_t D0Fitter::fgMeanDefault = 
                         TDatabasePDG::Instance()->GetParticle(421)->Mass();
const Double_t D0Fitter::fgUpFitLimitDefault = D0Fitter::fgMeanDefault+ .25; 
// const Double_t D0Fitter::fgLowFitLimitDefault = D0Fitter::fgMeanDefault-.25;
const Double_t D0Fitter::fgLowFitLimitDefault = 1.7;
Double_t D0Fitter::fgUpFitLimit  = D0Fitter::fgUpFitLimitDefault;
Double_t D0Fitter::fgLowFitLimit = D0Fitter::fgLowFitLimitDefault;
Double_t D0Fitter::fgBinSize = 0.01;//-1.0;

//*********************************************************

 D0Fitter::D0Fitter(TH1* hist) :fFitOption("ERI") // "M" leads to endless loop in ROOT < 4.X
{
  // default ctor

/*
fitoption
= "W"  Set all errors to 1
= "I" Use integral of function in bin instead of value at bin
                    center
= "L" Use Loglikelihood method (default is chisquare method)
= "U" Use a User specified fitting algorithm (via SetFCN)
= "Q" Quiet mode (minimum printing)
= "V" Verbose mode (default is between Q and V)
= "E" Perform better Errors estimation using Minos technique
= "M" More. Improve fit results
= "B" Take care of initial values of parameter, even for
                    predefined fcts as gaus, pol1,...
= "R" Use the Range specified in the function range
= "N" Do not store the graphics function, do not draw
= "0" Do not plot the result of the fit. By default the
             fitted function is drawn unless the
             option"N" above is specified.
= "+" Add this new fitted function to the list of fitted
             functions (by default, any previous function
             is deleted)
*/

  fLowFitLimit = fgLowFitLimitDefault;
  fUpFitLimit = fgUpFitLimitDefault;

  fFixMean = fFixSigma = fFixBackgrSlope = kFALSE;

  // FIXME
//   fD0FitFormu = new TF1("fGleichUmbenennen!", "expo(0) + gaus(2)",
// 			fLowFitLimit, fUpFitLimit);
  fD0FitFormu = new TF1("D0FitFormu", D0FitFct, fLowFitLimit, fUpFitLimit,5);
  gROOT->GetListOfFunctions()->Remove(fD0FitFormu);
  fD0FitFormu->SetBit(TFormula::kNotGlobal);

  this->SetHist(hist, kTRUE); // default fit values are set!
}

//___________________________________________
 D0Fitter::~D0Fitter()
{
  delete fD0FitFormu;
}

//___________________________________________
 void D0Fitter::Fit(Int_t mode)
{
  
  if(!fHist) {
    this->Error("Fit","No histogram! Fit canceled!");
    return;
  }

  Double_t backgrNormStart = this->DetermBackgrStartParam();
  Double_t nD0Start = this->DetermNDstarStart(backgrNormStart);
  this->DetermFitMode(mode); // which parameter should be fixed?
  this->DoTheFit(backgrNormStart, nD0Start);
  // add background function to hist:
  this->AddBackgrToHist();
}

//_______________________________________________
 void D0Fitter::DoTheFit(Double_t backgrNormStart, 
				    Double_t nD0Start)
{
  // does the fit having found startvalues

  TF1* f = new TF1("fGleichUmbenennen!", D0FitFct,
		   fLowFitLimit, fUpFitLimit, 5);
//   TF1* f = new TF1("fGleichUmbenennen!", "expo(0) + gaus(2)",
// 		   fLowFitLimit, fUpFitLimit);
  gROOT->GetListOfFunctions()->Remove(f);
  f->SetBit(TFormula::kNotGlobal);
  f->SetParameters(backgrNormStart, this->GetBackgrSlope(),
 		   this->GetMean(), this->GetSigma(), nD0Start);
  f->SetParNames("U_{n}","U_{s}", "#mu","#sigma","N(D^{0})");
  delete fD0FitFormu;
  fD0FitFormu = f;
  fD0FitFormu->SetName("D0FitFormu");

  if(fFixBackgrSlope) fD0FitFormu->SetParLimits(1,1,1);
  if(fFixMean)        fD0FitFormu->SetParLimits(2,1,1);
  if(fFixSigma)       fD0FitFormu->SetParLimits(3,1,1);

  fgBinSize = GetHistBinSize(); // important!!!!!!!!!!
  fgLowFitLimit = fLowFitLimit; // ...used in 
  fgUpFitLimit  = fUpFitLimit;  // ... backgrfct

  fHist->Fit(fD0FitFormu, 
	     this->UseLikelihoodOption(fHist) ? 
	     (fFitOption+"0l") : (fFitOption+'0'));
  if(!fFitOption.Contains("0")) {
    fHist->GetFunction("D0FitFormu")->ResetBit(1<<9);
  }
}

//*********************************************************

 const char* D0Fitter::GetHistName() const
{
  // returns name of the histogram that should be fitted
  if(fHist) return fHist->GetName();
  else {
    this->Warning("GetHistName","There is no histogram!");
    return "There is no histogram!";
  }
  
}

//******************************************************

 Double_t D0Fitter::GetHistBinSize() const
{
  // evaluates the bin size of the histogram
  // and prints warnings if it's not constant
  if(!fHist){
    this->Error("GetHistBinSize","There is no histogram! Return 100.");
    return 100.;
  }
  Int_t ind = fHist->FindBin(fgMeanDefault);
  Double_t binSize = fHist->GetBinLowEdge(ind+1) - fHist->GetBinLowEdge(ind);
  
  // I don't check binSize not to be 0,
  // but I do check that it's more or less constant:
  Int_t nBins = fHist->GetNbinsX();
  for(Int_t i = 1; i <= nBins; i++){
    Double_t actBinSize = 
      Double_t(fHist->GetBinLowEdge(i+1) - fHist->GetBinLowEdge(i));
    if(TMath::Abs(binSize - actBinSize)/binSize >0.0001){//0.01% deviat.
      this->Error("GetHistBinSize",
		  "histogram does NOT have constant binSize!");
      break;
    }
  }
  return binSize;
}

//_____________________________________
 void D0Fitter::SetHist(TH1* hist,
		       Bool_t defaultSigmaMeanBackgr)
{
  // interface to set a TH1* hist (may not be a TH2* or TH3*!)
  // if defaultSigmaMeanBack is kTRUE (its default is kFALSE),
  //     default values are set 
  //     for Sigma, Mean, amd background parameters
  // else they are taken as they are (e.g. from the last fit)
  // ??????
  // by default defaultSigmaMeanBackgr is kFALSE, if kTRUE, 
  // Sigma, Mean and background parameters are set to the default.

  fHist = hist;

  if(defaultSigmaMeanBackgr){ // sets defaults, 
    fD0FitFormu->SetParameters(10000., this->GetBackgrSlopeDefault(),
			       this->GetMeanDefault(), 
			       this->GetSigmaDefault(), 100.);
  } // else last values are used for next fit

  if(!hist) {
    return;
    //     this->Warning("SetHist", "No histogram set!");
  }  else if(hist->InheritsFrom(TH2::Class())||hist->InheritsFrom(TH3::Class())){
    this->Error("SetHist","Histogram given to fit is a TH2 or a TH3!");
  }

  cout << "Fitted histogram is " << fHist->GetName()  << endl;
}

//********************************************************

 Double_t D0Fitter::DetermBackgrStartParam()
{
  // determines start parameter for background fit
  Double_t startBackgr = 
    fHist->Integral() - 
    0.25 * fHist->Integral(fHist->FindBin(this->GetMean() - 2.5 * this->GetSigma()),// S/U = 1:3
			   fHist->FindBin(this->GetMean() + 2.5 * this->GetSigma()));
  cout << "Start parameter for backgr norm: " << startBackgr
       << endl;
  return startBackgr;
}

//________________________________
 Double_t D0Fitter::DetermNDstarStart(Double_t backgrNormStart)
{
  // determines start value of N(D0)
  Double_t signalLow = this->GetMean() - 2. * this->GetSigma();
  Double_t signalUp  = signalLow       + 4. * this->GetSigma();
  Int_t binD0Low = fHist->FindBin(signalLow);
  Int_t binD0Up = fHist->FindBin(signalUp);
  Double_t integral = fHist->Integral(binD0Low, binD0Up);

  TF1 backGr("backGrD0tmp", D0ExpBackgr, fLowFitLimit, fUpFitLimit, 2);
  backGr.SetParameters(backgrNormStart, this->GetBackgrSlope());

  fgBinSize = this->GetHistBinSize(); // fgBinSize is used in BackgrFct!
  fgLowFitLimit = fLowFitLimit; // dito used in 
  fgUpFitLimit  = fUpFitLimit;  // ... backgrfct

  Double_t backgrIntegr = 
    backGr.Integral(signalLow, signalUp) / fgBinSize;
  
  Double_t nD0Start = integral-backgrIntegr;
  if(nD0Start < 0.) nD0Start = 30.;
  cout << "Start parameter for N(D0): " << nD0Start
       << endl;
  return nD0Start;
}

//********************************************************

 void D0Fitter::AddBackgrToHist()
{
  TF1* backGr = new TF1("backGrD0", D0ExpBackgr, 
			fLowFitLimit, fUpFitLimit, 2);
  backGr->SetParameters(this->GetBackgrNorm(), this->GetBackgrSlope());
  backGr->SetParent(fHist);
  gROOT->GetListOfFunctions()->Remove(backGr);
  backGr->SetBit(TFormula::kNotGlobal);
  backGr->SetLineStyle(2);
  backGr->SetLineWidth(2);
  backGr->SetFillStyle(0); // no filling: Why do we need to set it here...? Since 3.04_02!
  TList* hFuncList = fHist->GetListOfFunctions();
  hFuncList->Add(backGr); // now fHist is responsible for 
  //                    deleting backGr
}

//********************************************************

 Bool_t D0Fitter::UseLikelihoodOption(TH1* hist)
{
  // determine use of likelihood option in fit
  // false if hists filled with weights
  // true if more than 1 bin from lowfitlimit to upfitlimit <= 5 entries

  if(hist->GetSumw2N()){ // hist has some error values stored...
    for(Int_t i = 1; i <= hist->GetNbinsX(); ++i){
      const Double_t error = hist->GetBinError(i);
      if(TMath::Abs(TMath::Abs(hist->GetBinContent(i)) - error*error) > 1.e-10){
#if ROOT_VERSION_CODE > ROOT_VERSION(3,3,0)
        this->Info("UseLikelihoodOption", "No likelihood (weights!)");
#else
        this->Warning("UseLikelihoodOption", "No likelihood (weights!)");
#endif
        return kFALSE;
      }
    }
  }

  Int_t numLowStatBins = 0;
  const Int_t lastBin = TMath::Min(hist->FindBin(this->GetUpFitLimit()), hist->GetNbinsX());
  for(Int_t i = TMath::Max(1, hist->FindBin(this->GetLowFitLimit())); i < lastBin; ++i){
    if(hist->GetBinContent(i) <= 5) {
      if(++numLowStatBins > 1) {
	cout << "Switched to likelihood!" << endl;
	return kTRUE;
      }
    }
  }

  return kFALSE;
}

//******************************************************

 void D0Fitter::DetermFitMode(Int_t mode)
{
  //interprets mode and saves the result in the foreseen Bool_t's
  fFixMean        =  mode        % 10;
  fFixSigma       = (mode /= 10) % 10;
  fFixBackgrSlope = (mode /= 10) % 10;
}

//******************************************************

 Double_t D0Fitter::GetBackgrNorm() const
{
  // get 
  return fD0FitFormu->GetParameter(0);
}

//******************************************************

 Double_t D0Fitter::GetBackgrSlope() const
{
  // get 
  return fD0FitFormu->GetParameter(1);
}

//******************************************************

 Double_t D0Fitter::GetMean() const
{
  // get mean of D0-Peak
  return fD0FitFormu->GetParameter(2);
}

//******************************************************

 Double_t D0Fitter::GetSigma() const
{
  // get sigma of D0-peak
  return fD0FitFormu->GetParameter(3);
}

//**********************************************

 Double_t D0Fitter::GetND0() const
{
  // get N(D0)
  return fD0FitFormu->GetParameter(4);
}

//******************************************************

 Double_t D0Fitter::GetND0Err() const
{
  // get sigma(N(D*))
  return fD0FitFormu->GetParError(4);
}

//******************************************************

 void D0Fitter::SetBackgrSlope(Double_t m)
{
  // set 
  fD0FitFormu->SetParameter(1, m);
}

//******************************************************

 void D0Fitter::SetMean(Double_t m)
{
  // set mean of D0-Peak
  fD0FitFormu->SetParameter(2, m);
}

//******************************************************

 void D0Fitter::SetSigma(Double_t s)
{
  // set sigma of D0-peak
  fD0FitFormu->SetParameter(3, s);
}