// Author: Gero Flucke <mailto:gero.flucke@desy.de>
//______________________________________________
// last change: $Date: 2005/11/26 20:04:56 $
// by         : $author$
//
// ABC for special fitting ('helper') classes 
// to fit the Delta-m distribution for D* analysis
// Use DstarDmFitter to do these fits!

#include "DstarDmFitterSpecial_I.h"

#include <TVirtualFitter.h>
#include <TROOT.h>
#include <TH1.h>
#include <TF1.h>
#include <TMath.h>
#include <TDatabasePDG.h>

#include <iostream>

using namespace std;

ClassImp(DstarDmFitterSpecial_I)

Double_t DstarDmFitterSpecial_I::fgSigmaDefault            = 0.0010;
Double_t DstarDmFitterSpecial_I::fgMeanDefault             = 0.14544;
Double_t DstarDmFitterSpecial_I::fgBackgrExpDefault        = 0.45;
Double_t DstarDmFitterSpecial_I::fgBackgrSqrDefault        =  18.; 
Double_t DstarDmFitterSpecial_I::fgUpperBackgrSqrDefault   =  36.; 
Double_t DstarDmFitterSpecial_I::fgLowerBackgrSqrDefault   =   0.; 
Double_t DstarDmFitterSpecial_I::fgUpperBackgrExpDefault   = 0.75; 
Double_t DstarDmFitterSpecial_I::fgLowerBackgrExpDefault   = 0.15; 
const Double_t DstarDmFitterSpecial_I::fgUpFitLimitDefault = 0.1675;
Double_t DstarDmFitterSpecial_I::fgUpFitLimit = 
DstarDmFitterSpecial_I::fgUpFitLimitDefault;
const Double_t DstarDmFitterSpecial_I::fgMassPion =
                            TDatabasePDG::Instance()->GetParticle(211)->Mass();
const Double_t DstarDmFitterSpecial_I::kSqrtOf2Pi = TMath::Sqrt(2.*TMath::Pi());
Double_t DstarDmFitterSpecial_I::fgBinSize = -1.0;     // garbage init
Double_t DstarDmFitterSpecial_I::fgShiftWcHist = -1.0; // garbage init
Double_t DstarDmFitterSpecial_I::fgEndRightCharge = -1.; // garbage init
Bool_t DstarDmFitterSpecial_I::fgDebug = kFALSE;


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

DstarDmFitterSpecial_I::DstarDmFitterSpecial_I() : fHist(NULL), fHistWc(NULL), 
  fIsMC(kFALSE), fDetectIsMC(kFALSE), fFitOption("EMRIL")
{
  // default ctor., determines only default fit option and IsMc flag

/*
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
= "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)
*/

}

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

 const char* DstarDmFitterSpecial_I::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 DstarDmFitterSpecial_I::GetWcBackgrNormScale() const
{
  // should be implemented in derived classes: 
  // if fitted with WC hist gives the ratio between the background normalisations: U_n,wc/U_n
  this->AbstractMethod("GetWcBackgrNormScale");

  return 0.;
}


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

 Double_t DstarDmFitterSpecial_I::GetHistBinSize() const
{
  // evaluates the bin size of the histogram (and -if existing- of WC hist)
  // 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);
  const 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!");
      return -100.;
    }
  }

  if(fHistWc){
    nBins = fHistWc->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",
		    "wrong charge histogram does NOT have constant binSize!");
	return -100.;
      }
    }
  }

  return binSize;
}


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

 TH1* DstarDmFitterSpecial_I::CreateMergedHist(const TH1* h, const TH1* hWc, 
					      Double_t& shiftWc) const
{
  // creates a hist with the bins of h, followed by an empty bin and then come the bins
  // of hWc. shiftWc is filled with the value that the bins egdes of hWc are shifted

  if(!h) {
    if(!hWc) return NULL; 
    else return NULL; // FIXME: we could do sth. here...
  } else if(!hWc){
    return static_cast<TH1*>(h->Clone());
  }

  TArrayD bins(h->GetNbinsX() + hWc->GetNbinsX() + 2);
  for(Int_t i = 0; i < h->GetNbinsX(); ++i){
    bins[i] = h->GetBinLowEdge(i+1);
  }
  bins[h->GetNbinsX()] = h->GetXaxis()->GetXmax();
  const Double_t startWcHist = bins[h->GetNbinsX()] + h->GetBinWidth(h->GetNbinsX());
  const Double_t offset = startWcHist - hWc->GetBinLowEdge(0);
  for(Int_t i = h->GetNbinsX() + 1; i < bins.GetSize()-1; ++i){
    bins[i] = hWc->GetBinLowEdge(i - h->GetNbinsX()) + offset;
  }
  bins[bins.GetSize()-1] = hWc->GetXaxis()->GetXmax() + offset;

  TH1* result = new TH1F(Form("%sAnd%s", h->GetName(), hWc->GetName()), 
			 Form("%s and %s", h->GetTitle(), hWc->GetTitle()),
			 bins.GetSize()-1, bins.GetArray());
  for(Int_t i = 1; i <= h->GetNbinsX(); ++i){
    result->SetBinContent(i, h->GetBinContent(i));
    if(h->GetSumw2N() || hWc->GetSumw2N()) result->SetBinError(i, h->GetBinError(i));
  }
  for(Int_t i = 1; i <= hWc->GetNbinsX(); ++i){
    result->SetBinContent(h->GetNbinsX() + i + 1, hWc->GetBinContent(i));
    if(h->GetSumw2N() || hWc->GetSumw2N()) 
      result->SetBinError(h->GetNbinsX() + i + 1, hWc->GetBinError(i));
  }
  
  shiftWc = offset;
  return result;
}


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

 void DstarDmFitterSpecial_I::SetHist(TH1* hist, Bool_t defaultSigmaMeanExpSq)
{
// for backward comp.
  this->SetHist(hist, NULL, defaultSigmaMeanExpSq);
}

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

 void DstarDmFitterSpecial_I::SetHist(TH1* hist, TH1* hWc, Bool_t defaultSigmaMeanExpSq)
{
  // interface to set a TH1* hist (may not be a TH2* or TH3*!) and WC background hist
  // if defaultSigmaMeanExpSq is kTRUE (it's default is kFALSE),
  //     default values are set 
  //     for Sigma, Mean, BackgrExp, BackgrSqr
  // else they are taken as they are (e.g. from the last fit)

  if(hist && hist->GetDimension() > 1) {
    this->Error("SetHist","Histogram given to fit has dimension %d!", hist->GetDimension());
  }
  if(hWc && hWc->GetDimension() > 1){
    this->Error("SetHist","Wrong charge histogram given to fit has dimension %d!", 
		hWc->GetDimension());
  }
  this->SetHistReal(hist, hWc, defaultSigmaMeanExpSq);
}

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

 Double_t DstarDmFitterSpecial_I::DetermBackgrStartParam()
{
  // determines start parameter for background fit
  // and checks whether fHist is backgroundless MC (but fNoIsMC == true): no MC assumption!
  // (for data determines only normalisation)
  
  // get integral of background region:
  Float_t backgrRangeLow = fgMeanDefault + 2.5*fgSigmaDefault;
  Int_t binLow = fHist->FindBin(backgrRangeLow);
  Int_t binHigh = fHist->FindBin(fUpFitLimit);
  Double_t backgrIntegr = fHist->Integral(binLow, binHigh);

  // get integral of signal region:
  Int_t binSigLow = fHist->FindBin(fgMeanDefault - 2.*fgSigmaDefault);
  Int_t binSigHigh = fHist->FindBin(fgMeanDefault + 2.*fgSigmaDefault);
  Double_t signalRegIntegr = fHist->Integral(binSigLow, binSigHigh);

  Double_t backgrNormStart = 0.;
  if(fDetectIsMC && (backgrIntegr / signalRegIntegr) < 0.2){ // MC ( ~0.04)
    cout << "assume MC: no background!" << endl;
    backgrNormStart= 0.;
    fIsMC = kTRUE;
  } else { // normal data 
    backgrNormStart =  fHist->Integral(1, binHigh) - 0.25 * signalRegIntegr;//assume S/U = 1:3
//     backgrNormStart = backgrIntegr / (binHigh-binLow+1);
//     backgrNormStart /= TMath::Power((fUpFitLimit - backgrRangeLow),fgBackgrExpDefault)
//       * GetHistBinSize();
    cout << "Start parameter of background normalisation: " 
	 << backgrNormStart << endl;
    fIsMC = kFALSE;
  }
  return backgrNormStart;
}

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

 Int_t DstarDmFitterSpecial_I::UseLikelihoodOption(TH1* hist) const
{
  // determine use of likelihood option in fit, dependent on given option and
  // hists with weights:
  // 2 if fitoption contains LL 
  // 0 if fitoption does not contain any L 
  //   or if exact one L, but filled with weights
  // 1 else
  //

  if(!fFitOption.Contains("L")) return 0;
  if(fFitOption.Contains("LL")) return 2;

  // so we have no 'L' but not 'LL': if we have weights, no likelihood!
  if(hist->GetSumw2N()){ // hist has some error values stored...
    for(Int_t i = 1; i <= hist->GetNbinsX(); ++i){
      if(!hist->GetBinContent(i)) continue; // skip zero bins - we might have set there error=1!
      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!), option 'LL' to enforce");
#else
	this->Warning("UseLikelihoodOption","No likelihood (weights!), option 'LL' to enforce");
#endif
	return 0;
      }
    }
  }

  return 1;

// true if more than 1 bin above fgMeanDefault+ 2.*fgSigmaDefault has <= 5 entries
//   Int_t numLowStatBins = 0;
//   const Int_t lastBin = TMath::Min(hist->FindBin(this->GetUpFitLimit()),  hist->GetNbinsX());
//   for(Int_t i = TMath::Max(hist->FindBin(fgMeanDefault+ 2.*fgSigmaDefault),1); i < lastBin;++i){
//     if(hist->GetBinContent(i) <= 5) {
//       if(++numLowStatBins > 1) {
// 	cout << "Switched to likelihood!" << endl;
// 	return 1;
//       }
//     }
//   }

//   return 0;
}

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

 Bool_t DstarDmFitterSpecial_I::TreatZeroErrorBins(TH1* hist) const
{
  // Assign an error to all bins with error = 0. to keep them, in the fit also for chi^2-fit.
  // The error is set to 1. if no squares of weights are stored with the hist, otherwise
  // to the mean of the errors of the first bin on the left and on the right that have 
  // error != 0. Only bins in the fitrange (m(pi) -> upfitlimit) are used for that.
  // False if no hist given or all bins in range have error == 0.

  if(!hist) return kFALSE;

  const Int_t firstBin = hist->FindBin(fgMassPion);
  const Int_t lastBin = TMath::Min(hist->FindBin(this->GetUpFitLimit()),  hist->GetNbinsX());

  const Bool_t storedWeights = hist->GetSumw2N();
  Int_t numProbBins = 0;
  TArrayD newError(lastBin + 1); // little longer than needed: no mix up of indices...
  for(Int_t i = firstBin; i <= lastBin; ++i){
    if(hist->GetBinError(i) == 0.){
      ++numProbBins;
      if(!storedWeights){
	newError[i] = 1.;
      } else {
	Int_t binLeft = TMath::Max(i - 1, firstBin);
	Double_t errorLeft = 0.;
	while(binLeft >= firstBin){
	  errorLeft = hist->GetBinError(binLeft--);// postfix decrement
	  if(errorLeft) break;
	}
	Int_t binRight = TMath::Min(i + 1, lastBin); 
	Double_t errorRight = 0.;
	while(binRight <= lastBin){
	  errorRight = hist->GetBinError(binRight++); // postfix increment
	  if(errorRight) break;
	}
	// newError[i] is stored and not directly set: binLeft could not have 0 any more...
	if(errorLeft && errorRight){
	  newError[i] = (errorLeft + errorRight)/2.;
	} else if(errorLeft){
	  newError[i] = errorLeft;
	} else if(errorRight){
	  newError[i] = errorRight;
	} else {
	  this->Warning("TreatZeroErrorBins",
			"No bin has any error, set err(%d) = max(1,sqrt(|bincontent|)", i);
	  newError[i] = TMath::Max(1., TMath::Sqrt(TMath::Abs(hist->GetBinContent(i))));
	}
      }
    }
  }

  if(numProbBins == lastBin - firstBin + 1){ // all bins had no error!
    return kFALSE;
  } else {
    Int_t manip = 0;
    for(Int_t i = firstBin; i <= lastBin; ++i){
      if(newError[i]){ // == 0 means: not manipulated
	if(fgDebug) this->Info("TreatZeroErrorBins", "manipulate error of bin %d to %f", 
			       i, newError[i]);
	hist->SetBinError(i, newError[i]);
	++manip;
      }
    }
    if(manip) this->Info("TreatZeroErrorBins", "manipulate error of %d bins", manip);
    return kTRUE;
  }
}

 TString DstarDmFitterSpecial_I::GetSimpleFitOption(const char* opt) const
{
  TString option(opt);
  option.ToUpper();
  option.ReplaceAll("I", NULL);
  option.ReplaceAll("M", NULL);
  option.ReplaceAll("E", NULL);
  option.ReplaceAll("L", NULL);
  option += "NQ"; // add 'Q'uiet
  return option;
}

 Bool_t DstarDmFitterSpecial_I::IsEquivalentOpt(const TString& optOrig, 
					       const TString& opt2) const
{
  TString origToManip(optOrig);
  origToManip.ToUpper();
  TString opt2ToManip(opt2);
  opt2ToManip.ToUpper();

  origToManip.ReplaceAll("N", NULL);
  opt2ToManip.ReplaceAll("N", NULL);
  origToManip.ReplaceAll("Q", NULL);
  opt2ToManip.ReplaceAll("Q", NULL);
  origToManip.ReplaceAll("0", NULL);
  opt2ToManip.ReplaceAll("0", NULL);


  if(origToManip.Contains("R") && opt2ToManip.Contains("R")){
    origToManip.ReplaceAll("R", NULL);
    opt2ToManip.ReplaceAll("R", NULL);
  }
  if(origToManip.Contains("I") && opt2ToManip.Contains("I")){
    origToManip.ReplaceAll("I", NULL);
    opt2ToManip.ReplaceAll("I", NULL);
  }
  if(origToManip.Contains("M") && opt2ToManip.Contains("M")){
    origToManip.ReplaceAll("M", NULL);
    opt2ToManip.ReplaceAll("M", NULL);
  }
  if(origToManip.Contains("E") && opt2ToManip.Contains("E")){
    origToManip.ReplaceAll("E", NULL);
    opt2ToManip.ReplaceAll("E", NULL);
  }
  if(origToManip.Contains("LL") && opt2ToManip.Contains("LL")){
    origToManip.ReplaceAll("LL", NULL);
    opt2ToManip.ReplaceAll("LL", NULL);
  } else if(origToManip.CountChar('L') == 1 && opt2ToManip.CountChar('L') == 1){//not OK for LXL
    origToManip.ReplaceAll("L", NULL);
    opt2ToManip.ReplaceAll("L", NULL);
  }
  
  return (origToManip.Length() == 0 && opt2ToManip.Length() == 0);
}

 Int_t DstarDmFitterSpecial_I::RealFit(TH1* histToFit, TF1* func) const
{
  // calling histToFit->Fit(func, option) where option first is a simple version
  // of fFitOption and then increasing (L, I, E, M) when given/needed
  // return result of histToFit->Fit(func, fFitOption) or the first failing of prior tries.
  // If likelihood method is NOT used, call TreatZeroErrorBins for fHist and fHistWc,
  // if there is a problem, return -10 and don't fit.

  //  TVirtualFitter::SetMaxIterations(500);
  //  TVirtualFitter::SetPrecision();//default
  //  TVirtualFitter::SetPrecision(5.*TVirtualFitter::GetPrecision());

  TString fitOpt = this->GetSimpleFitOption(fFitOption);
  if(fHist != histToFit) this->Warning("RealFit", "test only right charge hist for liklihood");
  const Int_t useLikelihood = this->UseLikelihoodOption(fHist);
  if(useLikelihood == 0){
    if(!this->TreatZeroErrorBins(fHist)){
      this->Error("RealFit", "Cannot fit hist with all errors = 0");
      return -10;
    }
    if(fHistWc && !this->TreatZeroErrorBins(fHistWc)){
      this->Error("RealFit", "Cannot fit histWc with all errors = 0");
      return -10;
    }
  }


  Int_t fitResult = this->RealFit(histToFit, func, fitOpt); 
  if(fitResult == 0 && useLikelihood) {
    fitOpt += 'L';
    if(useLikelihood == 2) fitOpt += 'L';
    if(this->IsEquivalentOpt(fitOpt, fFitOption)) fitOpt = fFitOption+'0';
    fitResult = this->RealFit(histToFit, func, fitOpt); 
  }
  if(fitResult == 0 && fFitOption.Contains("I", TString::kIgnoreCase) 
     && fitOpt != (fFitOption+'0')){
    fitOpt += 'I';
    if(this->IsEquivalentOpt(fitOpt, fFitOption)) fitOpt = fFitOption+'0';
    fitResult = this->RealFit(histToFit, func, fitOpt); 
  }
  if(fitResult == 0 && fFitOption.Contains("M", TString::kIgnoreCase) 
     && fitOpt != (fFitOption+'0')){
    fitOpt += 'M';
    if(this->IsEquivalentOpt(fitOpt, fFitOption)) fitOpt = fFitOption+'0';
    fitResult = this->RealFit(histToFit, func, fitOpt); 
  }
  if(fitResult == 0 && fFitOption.Contains("E", TString::kIgnoreCase) 
     && fitOpt != (fFitOption+'0')){
    fitOpt += 'E';
    if(this->IsEquivalentOpt(fitOpt, fFitOption)) fitOpt = fFitOption+'0';
    fitResult = this->RealFit(histToFit, func, fitOpt); 
  }
  if(fitResult == 0 && fitOpt != (fFitOption+'0')){
    fitOpt = fFitOption+'0';
    if(useLikelihood == 0) fitOpt.ReplaceAll("L", NULL);
    fitResult = this->RealFit(histToFit, func, fitOpt); 
  }
  
  if(fitResult != 0){
    this->Error("RealFit", "fit finally unsuccessful, last used option %s", fitOpt.Data());
    delete histToFit->GetListOfFunctions()->Remove(histToFit->FindObject(func->GetName()));
  } else if(!fFitOption.Contains("0")){
    histToFit->GetFunction(func->GetName())->ResetBit(TF1::kNotDraw);
  }

  return fitResult;
}

 Int_t DstarDmFitterSpecial_I::RealFit(TH1* histToFit, TF1* func, const char* fitOpt) const
{
  // fit hist with func and option
  // first time displaces start parameter by factor 0.98, if no sucess tries several times
  // if no success, func's params are reset to old values
  // return fit result of last attempt

  if(fgDebug) this->Info("RealFit","fit option %s", fitOpt);

  Double_t params[50];// 7 is maximum...
  func->GetParameters(params);
  this->MultiplyFreeParams(func, 0.98);
  Int_t fitResult = histToFit->Fit(func, fitOpt);
  if(fitResult != 0){
    TString nonQuietOpt(fitOpt);
    nonQuietOpt.ToUpper();
    nonQuietOpt.ReplaceAll("Q", NULL);
    Long_t tries = 1;// Long_t for TMath::Odd
    do {
      this->Warning("RealFit","%s failed, try again (%d) with modified parameters (non quiet)", 
		    fitOpt, tries);
      ++tries;
      func->SetParameters(params);
      Double_t factor = 1.;
      if(TMath::Odd(tries)) factor += (tries * 0.01);// first time: tries == 2, so 1.02, 1.04,
      else factor -= (tries * 0.01);// .97, .95, .93
      this->MultiplyFreeParams(func, factor);
      fitResult = histToFit->Fit(func, nonQuietOpt);
      if(fitResult == 0) break;
    } while (tries < 10);
  }
  
  if(fitResult != 0) func->SetParameters(params); // reset to old values

  return fitResult;
}

 Int_t DstarDmFitterSpecial_I::MultiplyFreeParams(TF1* func, Double_t factor) const
{
  // multiply all non-fixed parameters of func by factor, returns number of non-free params
  Int_t changed = 0;
  for(Int_t i = 0; i < func->GetNpar(); ++i){
    Double_t parMin = 0., parMax = 0.;
    func->GetParLimits(i, parMin, parMax);
    Bool_t fixed = kTRUE;
    if(!(parMin*parMax != 0. && parMin >= parMax)){ // cf.  TF1::GetNumberFreeParameters()
      func->SetParameter(i, func->GetParameter(i) * factor);// set if not fixed
      ++changed;
      fixed = kFALSE;
    }
    if(!fixed && func->GetParameter(i) <= 0. && strcmp(func->GetParName(i), "#sigma") == 0){
      this->Info("MultiplyFreeParams", "#sigma manipulated from %f to default %f", 
		 func->GetParameter(i), fgSigmaDefault);
      func->SetParameter(i, fgSigmaDefault);
    }
    if(!fixed && strcmp(func->GetParName(i), "#mu") == 0 && 
       (func->GetParameter(i) < 0.14 || func->GetParameter(i) > 0.15) ){
      this->Info("MultiplyFreeParams", "#mu manipulated from %f to default %f", 
		 func->GetParameter(i), fgMeanDefault);
      func->SetParameter(i, fgMeanDefault);
    }

    if(fgDebug) {
      cout << func->GetParName(i) << " = " << func->GetParameter(i) << " ";
      if(fixed) cout << "(fix) ";
    }
  }
  if(fgDebug) cout << " (from MultiplyFreeParams)" << endl;

  if(changed != func->GetNumberFreeParameters()){
    this->Error("MultiplyFreeParams", "Changed %d, not %d parameters.", 
		changed, func->GetNumberFreeParameters());
  }

  return changed;
}