// Author: Gero Flucke <mailto:gero.flucke@desy.de>
//_______________________________
//
//   last update: $Date: 2005/11/26 18:47:45 $
//   by:          $Author: flucke $
//
// fitting class for D*-Analysis:
//                  fits delta-mass-plots (m_{Kpipis} - m_{Kpi})
//  (for use in ROOT 3.00_06 or higher only!) 
//
// performs the 'smart' fits with fixed mean, 
// sigma and U_sqr for DstarDmFitter,
// so please use that class! 

#include "DstarDmFitterMeSiSq.h"
#include "DstarDmFitterSpecial_I.h"

#include <TH1.h>
#include <TF1.h>
#include <TMath.h>

#include <iostream>

using namespace std;

ClassImp(DstarDmFitterMeSiSq)
  
Double_t DstarDmFitterMeSiSq::fgSigma     = DstarDmFitterMeSiSq::GetSigmaDefault();
Double_t DstarDmFitterMeSiSq::fgMean      = DstarDmFitterMeSiSq::GetMeanDefault();
Double_t DstarDmFitterMeSiSq::fgBackgrSqr = DstarDmFitterMeSiSq::GetBackgrSqrDefault();


//********************************************************
//******* constructor, destructor, set and get ***********
//********************************************************

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

 DstarDmFitterMeSiSq::DstarDmFitterMeSiSq(TH1* hist)
{
  // contructor: hist is the histogram to fit (may not be a TH2*!)

  // initialise some data members:
  fIsMC = kFALSE;
  this->SetUpFitLimit();

  fDstarFitFormu = new TF1("dstarFitFormuMss", FitFct, 0.135, fUpFitLimit,3);

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

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

 DstarDmFitterMeSiSq::DstarDmFitterMeSiSq(const DstarDmFitterMeSiSq& source)
{
  fHist = source.GetHistPtr();
  this->SetUpFitLimit(source.GetUpFitLimit());
  fIsMC = source.HistIsMC();
  this->SetDetectIsMC(source.IsDetectIsMC());
  // fgBinSize???
  fFitOption = source.fFitOption;

  fDstarFitFormu = new TF1("dstarFitFormuMss", FitFct, 0.135, fUpFitLimit, 3);
  fDstarFitFormu->SetParameters(source.GetBackgrNorm(),	source.GetBackgrExp(),
				source.GetNDstar());

  fgBackgrSqr = source.GetBackgrSqr();
  fgSigma     = source.GetSigma();
  fgMean      = source.GetMean();
  
  fDstarFitFormu->SetParError(0, source.GetBackgrNormErr());
  fDstarFitFormu->SetParError(2, source.GetNDstarErr());
}

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

 DstarDmFitterMeSiSq::DstarDmFitterMeSiSq(const DstarDmFitterSpecial_I& source)
{
  fHist = source.GetHistPtr();
  this->SetUpFitLimit(source.GetUpFitLimit());

  fIsMC = source.HistIsMC();
  this->SetDetectIsMC(source.IsDetectIsMC());
  // fgBinSize???
  fFitOption = source.GetFitOption();

  fDstarFitFormu = new TF1("dstarFitFormuMss", FitFct, 0.135, fUpFitLimit, 3);
  fDstarFitFormu->SetParameters(source.GetBackgrNorm(),	source.GetBackgrExp(),
				source.GetNDstar());

  fgBackgrSqr = source.GetBackgrSqr();
  fgSigma     = source.GetSigma();
  fgMean      = source.GetMean();
  
  fDstarFitFormu->SetParError(0, source.GetBackgrNormErr());
  fDstarFitFormu->SetParError(2, source.GetNDstarErr());
}

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

 DstarDmFitterMeSiSq::~DstarDmFitterMeSiSq()
{
  //deletes all possibly created objects that members point at
  if(fDstarFitFormu)  {delete fDstarFitFormu; fDstarFitFormu = 0;}
}

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

DstarDmFitterMeSiSq& DstarDmFitterMeSiSq::operator=
(const DstarDmFitterMeSiSq& rhs)
{
  // get the same/aequivalent data member values from rhs
  if (this != &rhs) {
    fHist = rhs.GetHistPtr();
    this->SetUpFitLimit(rhs.GetUpFitLimit());
    fIsMC = rhs.HistIsMC();
    this->SetDetectIsMC(rhs.IsDetectIsMC());

    fDstarFitFormu->SetParameter(0, rhs.GetBackgrNorm());
    fDstarFitFormu->SetParError(0, rhs.GetBackgrNormErr());

    this->SetBackgrExp(rhs.GetBackgrExp());
    this->SetBackgrSqr(rhs.GetBackgrSqr());
    this->SetMean(rhs.GetMean());
    this->SetSigma(rhs.GetSigma());

    fDstarFitFormu->SetParameter(2, rhs.GetNDstar());
    fDstarFitFormu->SetParError(2, rhs.GetNDstarErr());
  }
  return (*this);
}

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

DstarDmFitterMeSiSq& DstarDmFitterMeSiSq::operator=
(const DstarDmFitterSpecial_I& rhs)
{
  // get the same/aequivalent data member values from rhs
  if (this != &rhs) {
    fHist = rhs.GetHistPtr();
    this->SetUpFitLimit(rhs.GetUpFitLimit());
    fIsMC = rhs.HistIsMC();
    this->SetDetectIsMC(rhs.IsDetectIsMC());

    fDstarFitFormu->SetParameter(0, rhs.GetBackgrNorm());
    fDstarFitFormu->SetParError(0, rhs.GetBackgrNormErr());

    this->SetBackgrExp(rhs.GetBackgrExp());
    this->SetBackgrSqr(rhs.GetBackgrSqr());
    this->SetMean(rhs.GetMean());
    this->SetSigma(rhs.GetSigma());

    fDstarFitFormu->SetParameter(1, rhs.GetNDstar());
    fDstarFitFormu->SetParError(1, rhs.GetNDstarErr());
  }
  return (*this);
}

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

 void DstarDmFitterMeSiSq::SetHistReal(TH1* histo, TH1* hDummy,
				Bool_t defaultSigmaMeanExpSqr)
{
  // by default defaultSigmaMeanExpSqr is kFALSE, 
  // if kTRUE,
  // Sigma, Mean, BackgrExp and BackgrSqr are set to the default.

  if((fHist = histo)){
    cout << "Fitted histogram is " << fHist->GetName()  << endl;
  }

  if(hDummy) {
    this->Error("SetHistReal", "WC-background hist not yet supported, ignore %s",
		hDummy->GetName());
  }

  if(defaultSigmaMeanExpSqr){ // sets defaults, 
    fgSigma     = fgSigmaDefault;
    fgMean      = fgMeanDefault;
    fDstarFitFormu->SetParameter(1, fgBackgrExpDefault);
    fgBackgrSqr = fgBackgrSqrDefault;
  } // else last values are used for next fit
  // take care that in that case no other instance 
  // has changed the values of fgSigma, fgMean, fgBackgrSqr
}

//******************************************************
  
 void DstarDmFitterMeSiSq::SetBackgrExp(Double_t e)
{
  // set exponent of background
  fDstarFitFormu->SetParameter(1, e);
}

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

 Double_t DstarDmFitterMeSiSq::GetBackgrNorm() const
{
  // get normalisation factor of background
  return fDstarFitFormu->GetParameter(0);
}

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

 Double_t DstarDmFitterMeSiSq::GetBackgrNormErr() const
{
  // get error of normalisation factor of background
  return fDstarFitFormu->GetParError(0);
}

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

 Double_t DstarDmFitterMeSiSq::GetBackgrExp() const
{
  // get exponent of background
  return fDstarFitFormu->GetParameter(1);
}

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

 Double_t DstarDmFitterMeSiSq::GetNDstar() const
{
  // get N(D*)
  return fDstarFitFormu->GetParameter(2);
}

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

 Double_t DstarDmFitterMeSiSq::GetNDstarErr() const
{
  // get sigma(N(D*))
  return fDstarFitFormu->GetParError(2);
}

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

 Double_t DstarDmFitterMeSiSq::GetSignalBack() const
{
  // Get ratio of signal over background of last fit
  // in (fgMean +- 2*fgSigma)-region.
  // Gives 0.0 if MC or any other strange result,
  // but doesn't check whether you have not fitted yet.
  // may give wrong results if another instance of this class 
  // has set it's mean etc.!
  if(fIsMC || !fDstarFitFormu) return 0.;

  Double_t signalRegionLow = this->GetMean() - 2.*this->GetSigma();
  Double_t signalRegionHigh = this->GetMean() + 2.*this->GetSigma();

  TF1 backGr("backGrMss", BackgrFct, 0.135, fUpFitLimit, 2);
  backGr.SetParameter(0, this->GetBackgrNorm());
  backGr.SetParameter(1, this->GetBackgrExp());

  fgBinSize = this->GetHistBinSize(); // needed to ensure 
  fgUpFitLimit = fUpFitLimit;         // dito

  Double_t backgrIntegr = backGr.Integral(signalRegionLow, signalRegionHigh);
  Double_t signalIntegr = fDstarFitFormu->Integral(signalRegionLow,
						   signalRegionHigh);

  return 
    backgrIntegr ? (signalIntegr - backgrIntegr) / backgrIntegr : 0.0;
}
//******************************************************

 TF1* DstarDmFitterMeSiSq::GetFittedFunc() const 
{
  // pointer to the fitted function associated with the histogram
  TH1* h = this->GetHistPtr();
  return (h ? h->GetFunction("dstarFitFormuMss") : NULL);
}

//******************************************************
 TF1* DstarDmFitterMeSiSq::GetFittedBackgrFunc() const
{
  // pointer to the background part of the fitted function associated with the histogram
  TH1* h = this->GetHistPtr();
  return (h ? h->GetFunction("backGrMss") : NULL);
}


/*********************************************************/
/****************** the fit: *****************************/
/*********************************************************/

 Bool_t DstarDmFitterMeSiSq::Fit(Int_t mode)
{
  // performs the Delta-m fit.
  // mode is ignored: depends on class itself!
  if(!fHist) {
    this->Error("Fit", "No histogram! Fit canceled!");
    return kFALSE;
  }
  if(fHistWc){
    this->Error("Fit", "smart fit does not support wrong charge hist!");
  }

  const Double_t backgrNormStart = this->DetermBackgrStartParam();//also fIsMC set!
  const Double_t nDstarStart = this->DetermNDstarStart(backgrNormStart);

  const Int_t fitResult = this->DoTheFit(backgrNormStart, nDstarStart);
  if(fitResult) return kFALSE;
  this->AddBackgrToHist();
  return kTRUE;
}

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

 void DstarDmFitterMeSiSq::AddBackgrToHist()
{
  TF1* backGr = new TF1("backGrMss", BackgrFct, 
			fgMassPion, fUpFitLimit, 2);
  backGr->SetParameters(this->GetBackgrNorm(), this->GetBackgrExp());
  backGr->SetParent(fHist);
  backGr->SetBit(TFormula::kNotGlobal);
  backGr->SetLineStyle(2);
  backGr->SetLineWidth(2);
  backGr->SetFillStyle(0); // no filling: Why do we need to set it here...?
  TList* hFuncList = fHist->GetListOfFunctions();
  hFuncList->Add(backGr); // now fHist is responsible for deleting backGr
}


/********************************************************/
/******** private member functions: *********************/
/********************************************************/

 Int_t DstarDmFitterMeSiSq::DoTheFit(Double_t backgrNormStart, 
			     Double_t nDstarStart)
{
  // does the fit having found startvalues
  if(fDstarFitFormu) delete fDstarFitFormu; // to get rid of setlimits!
  fDstarFitFormu = new TF1("dstarFitFormuMss", FitFct,
			  fgMassPion, fUpFitLimit, 3);
  fDstarFitFormu->SetParNames("U_{n}","U_{exp}","N(D*)");
  fDstarFitFormu->SetParameters(backgrNormStart, this->GetBackgrExp(), nDstarStart);

  if(fIsMC) {
    fDstarFitFormu->SetParLimits(0,1,1);
    fDstarFitFormu->SetParLimits(1,1,1);
  }

  fgBinSize = this->GetHistBinSize(); // important!!!!!!!!
  fgUpFitLimit = fUpFitLimit;         // dito
  if(fgBinSize <= 0.) {
    return -1;
  }

  Double_t oldParams[10]; // we have max 7 parameters...
  fDstarFitFormu->GetParameters(oldParams);
  const Int_t fitResult = this->RealFit(fHist, fDstarFitFormu);
  if(fitResult != 0)  fDstarFitFormu->SetParameters(oldParams);

  return fitResult;
}

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

 Double_t DstarDmFitterMeSiSq::DetermNDstarStart(Double_t backgrNormStart)
{
  // determines start value of N(D*)
  Int_t binDsLow = fHist->FindBin(fgMeanDefault - 2.*fgSigmaDefault);
  Int_t binDsUp = fHist->FindBin(fgMeanDefault + 2.*fgSigmaDefault);
  Double_t integral = fHist->Integral(binDsLow, binDsUp);

  TF1 backGr("backGrMss", BackgrFct, fgMassPion, fUpFitLimit, 2);
  backGr.SetParNames("U_{n}", "U_{exp}");
  backGr.SetParameters(backgrNormStart, this->GetBackgrExp());

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

  Double_t backgrIntegr = 
    backGr.Integral(fgMeanDefault - 2.*fgSigmaDefault, 
		      fgMeanDefault + 2.*fgSigmaDefault) / fgBinSize;

  return integral-backgrIntegr > 0. ? integral-backgrIntegr : 30.;
}

/********************************************************/
/************ functions used for fits: ******************/
/********************************************************/

 Double_t DstarDmFitterMeSiSq::BackgrFct(Double_t *x, Double_t *par)
{
  // background fitting fct:
  // par[0]: normalisation
  // par[1]: exponent
  // x[0]:   deltaM

  if(x[0] <= fgMassPion) {
    TF1::RejectPoint();
    return 0.;
  }

  // calculate normalisation factor (cf. log book 2003-04-05)
  // (integral from fgMassPion to upfitlimit == par[0])
  Double_t relLimit = fgUpFitLimit - fgMassPion;
  Double_t relPowP1Plus1 = TMath::Power(relLimit, par[1]+1.); 
  Double_t norm = (par[1] == 1. || par[1] == 2. || par[1] == 3.) ? 0. :
    (1.-fgBackgrSqr*fgMassPion*fgMassPion)    *relPowP1Plus1 / (par[1]+1.)
    - fgBackgrSqr * (relPowP1Plus1        *relLimit*relLimit / (par[1]+3.)
		+ 2.*fgMassPion* relPowP1Plus1*relLimit / (par[1]+2.)
		);
  
  if(norm != 0.) {
    return (par[0] * fgBinSize / norm) * TMath::Power((x[0] - fgMassPion), par[1])
      * (1. - fgBackgrSqr * x[0]*x[0]);
  } else {
    TF1::RejectPoint();
    return 0.;
  }

//   return (x[0] - fgMassPion) < 0. ? 0. :
//     par[0] * fgBinSize
//     * TMath::Power((x[0] - fgMassPion), par[1]) 
//     * (1. - fgBackgrSqr * x[0]*x[0]);
}

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

 Double_t DstarDmFitterMeSiSq::FitFct(Double_t* x, Double_t* par)
{
  // fitting function:
  //par[0]: normalisation of background
  //par[1]: exponent
  //par[2]: N(D*)
  //x[0]:   deltaM

  if(x[0] <= fgMassPion) {// lower than mPion not possible for deltaM!
    TF1::RejectPoint();
    return 0.;
  }
  // fgBinSize is a static member of DstarDmFitterMeSiSq 
  // 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 DstarDmFitterMeSiSq)!
  // [It has to be static because it is used 
  // in this static function DstarDmFitterMeSiSq::FitFct(...)!]

  // DstarDmFitterSpecial_I::kSqrtOf2Pi = sqrt(2*pi) in Double_t-precision
  Double_t gaus = 1./(DstarDmFitterSpecial_I::kSqrtOf2Pi*fgSigma) * fgBinSize * par[2]
    * TMath::Exp(-((x[0]-fgMean)*(x[0]-fgMean)) / (2.*fgSigma*fgSigma));
  return gaus + BackgrFct(x, par);
}