#include<vector>
#include<iostream>

#include "TFile.h"
#include "TTree.h"
#include "TBranch.h"
#include "TFitResultPtr.h"
#include "Fit/UnBinData.h"
#include "Math/MinimizerOptions.h"
#include "HFitInterface.h"
#include "Fit/DataRange.h"
#include "Fit/FitData.h"
#include "Fit/Fitter.h"
#include "Math/WrappedMultiTF1.h"
#include "Math/Minimizer.h"

#include "TF1.h"
#include "TH1.h"
#include "TGraphAsymmErrors.h"

#define Blind

int Binning    = 30;
Double_t xMin  = 4766;
Double_t xMax  = 5966;
Double_t xbMin = 5166;
Double_t xbMax = 5526;

class ScaledTF1Wrapper
{
public:
  int _color;
  ScaledTF1Wrapper    ()                        {_myTF1 = NULL  ; _color = 0;};
  ScaledTF1Wrapper    (TF1 *theTF1)             {_myTF1 = theTF1; _color = _myTF1->GetLineColor();};
  void SetTF1         (TF1 *theTF1)             {_myTF1 = theTF1; _color = _myTF1->GetLineColor();};
  int  GetColor       ()                        {return _color;};
  Double_t theFunction(Double_t *x, Double_t *p){return p[_myTF1->GetNpar()]*_myTF1->EvalPar(x,p);};

private:
  TF1 *_myTF1;
};

void Plotter(std::vector<Double_t> DataPoints, TF1 *Function)
{
  TH1D *histo = new TH1D("histo","Histogram filled with the data",Binning,Function->GetXmin(),Function->GetXmax());
  histo->SetBinErrorOption(TH1::kPoisson);
  for(uint i = 0; i<DataPoints.size(); i++) histo->Fill(DataPoints[i]);
  std::vector<Double_t> X,Y,EYLow,EYHigh,EXLow,EXHigh;
  int NumberOfPoints = 0;
  for(int i = 1; i<=histo->GetNbinsX(); i++)
    {
      X     .push_back(histo->GetBinCenter(i));
      Y     .push_back(histo->GetBinContent(i));
      EYLow .push_back(histo->GetBinErrorLow(i));
      EYHigh.push_back(histo->GetBinErrorUp(i));
      EXLow .push_back(histo->GetXaxis()->GetBinWidth(i)/2.);
      EXHigh.push_back(histo->GetXaxis()->GetBinWidth(i)/2.);
      NumberOfPoints++;
    }
  TGraphAsymmErrors *AsymmPlot = new TGraphAsymmErrors(NumberOfPoints,&X[0],&Y[0],&EXLow[0],&EXHigh[0],&EYLow[0],&EYHigh[0]);
  AsymmPlot->SetMarkerStyle(20);
  AsymmPlot->SetMarkerSize(1.1);
  AsymmPlot->SetTitle("A ROOT plot of \"mass\"");
  AsymmPlot->GetXaxis()->SetRangeUser(histo->GetXaxis()->GetXmin(),histo->GetXaxis()->GetXmax());
  AsymmPlot->GetYaxis()->SetRangeUser(0,histo->GetMaximum()+(histo->GetMaximum()/10.));
  AsymmPlot->GetXaxis()->SetTitle("mass");
  char YTitle[1000];
  sprintf(YTitle,"Events / ( %0.1f )",(Function->GetXmax()-Function->GetXmin())/Binning);
  AsymmPlot->GetYaxis()->SetTitle(YTitle);
  ScaledTF1Wrapper *FunctionScaler = new ScaledTF1Wrapper(Function);
  TF1 *ScaledFunction = new TF1("ScaledFunction",FunctionScaler,&ScaledTF1Wrapper::theFunction,Function->GetXmin(),Function->GetXmax(),Function->GetNpar()+1,"ScaledTF1Wrapper","theFunction");
  ScaledFunction->SetParameter(Function->GetNpar(),histo->GetBinWidth(1));
  ScaledFunction->SetParName  (Function->GetNpar(),"BinningScale");
  for(int i = 0; i<Function->GetNpar(); i++)
    {
      ScaledFunction->SetParameter(i,Function->GetParameter(i));
      ScaledFunction->SetParName  (i,Function->GetParName(i));
    }
  ScaledFunction->SetLineColor(4);
  ScaledFunction->SetLineWidth(3);
  AsymmPlot     ->Draw("AP");
  ScaledFunction->Draw("same");
}

ROOT::Fit::FitResult RootFitter(TF1 *FitFunction, Double_t RegionMin, Double_t RegionMax, bool ExtendedFit, std::vector<Double_t> FitData)
{
  //Unbinned Data set
  ROOT::Fit::DataOptions DataOptions;
  //ROOT::Fit::DataRange RangeOfData(RegionMin,RegionMax);
  ROOT::Fit::DataRange RangeOfData;
  RangeOfData.AddRange(0,RegionMin,RegionMax);
  ROOT::Fit::UnBinData *UnBinnedData = new ROOT::Fit::UnBinData(DataOptions,RangeOfData,FitData.size(),1,false);
  for(uint i = 0; i<FitData.size(); i++) UnBinnedData->Add(FitData[i]);
  ROOT::Math::WrappedMultiTF1 myFitFunction(*FitFunction,FitFunction->GetNdim());
  ROOT::Fit::Fitter myFitter;
  myFitter.SetFunction(myFitFunction,false);
  ROOT::Fit::FitConfig &myFitConfig               = myFitter.Config();
  ROOT::Math::MinimizerOptions myMinimizerOptions = myFitter.Config().MinimizerOptions();

  Double_t FitParametersFromTF1[FitFunction->GetNpar()], FitParametersErrorsFromTF1[FitFunction->GetNpar()];
  for(int i = 0; i<FitFunction->GetNpar(); i++)
    {
      FitParametersFromTF1[i]       = FitFunction->GetParameter(i);
      FitParametersErrorsFromTF1[i] = FitFunction->GetParError(i);
    }
  myFitConfig.SetParamsSettings(FitFunction->GetNpar(),FitParametersFromTF1,FitParametersErrorsFromTF1);
  
  for(int i = 0; i<FitFunction->GetNpar(); i++)
    {
      ROOT::Fit::ParameterSettings &FitParameterSettings = myFitConfig.ParSettings(i);
      Double_t FitParameterLowLimit, FitParameterUpLimit;
      FitFunction->GetParLimits(i,FitParameterLowLimit,FitParameterUpLimit);
      if(FitParameterLowLimit*FitParameterUpLimit != 0 && FitParameterLowLimit >= FitParameterUpLimit) FitParameterSettings.Fix();
      else if(FitParameterLowLimit < FitParameterUpLimit)
	{
	  if(!TMath::Finite(FitParameterUpLimit) && TMath::Finite(FitParameterLowLimit)) FitParameterSettings.SetLowerLimit(FitParameterLowLimit);
	  else if(!TMath::Finite(FitParameterLowLimit) && TMath::Finite(FitParameterUpLimit)) FitParameterSettings.SetUpperLimit(FitParameterUpLimit);
	  else FitParameterSettings.SetLimits(FitParameterLowLimit,FitParameterUpLimit);
	}
    }
  
  myMinimizerOptions.SetPrintLevel(4);
  myMinimizerOptions.SetTolerance(1.00);
  myMinimizerOptions.SetMaxFunctionCalls(2000);
  myMinimizerOptions.Print();
  myFitConfig.SetMinimizerOptions(myMinimizerOptions);
  
  myFitter.Fit(*UnBinnedData,ExtendedFit);
  myFitter.GetMinimizer()->PrintResults();
  myFitter.CalculateHessErrors();
  const ROOT::Fit::FitResult &fitResult = myFitter.Result();
  for(uint i = 0; i<fitResult.NPar(); i++)
    {
      FitFunction->SetParameter(i,fitResult.Parameter(i));
      FitFunction->SetParError(i,fitResult.ParError(i));
    }
  std::cout << "Fit status: " << fitResult.Status() << std::endl;
  fitResult.Print(std::cout);
  fitResult.PrintCovMatrix(std::cout);
  return fitResult;
}

Double_t ChebNorm(Double_t aCheb)
{
#ifdef Blind
  return (-xMin*xMin+2*xMin*xMax-xMax*xMax+(1+aCheb)*xMin*(xbMin-xbMax)+(-1+aCheb)*xMax*(xbMin-xbMax)+aCheb*(-xbMin*xbMin+xbMax*xbMax))/(xMin-xMax);
#else
  return xMax - xMin;
#endif
}

Double_t ExpoNorm(Double_t aExpo)
{
#ifdef Blind
  return ((exp(aExpo * xMax) - exp(aExpo * xMin)) - (exp(aExpo * xbMax) - exp(aExpo * xbMin)))/aExpo;
#else
  return (exp(aExpo * xMax) - exp(aExpo * xMin))/aExpo;
#endif
}

Double_t CombinedModel(Double_t *x, Double_t *FitParameters)
{
  Double_t xx    = x[0];
  Double_t aCheb = FitParameters[0];
  Double_t aExpo = FitParameters[1];
  Double_t nCheb = FitParameters[2];
  Double_t nExpo = FitParameters[3];
  Double_t Ntot  = nExpo+nCheb;

  Double_t Exponential = (1/ExpoNorm(aExpo))*exp(aExpo*xx);
  Double_t Chebichev   = (1/ChebNorm(aCheb))*(1.0+aCheb*(2*xx-xMin-xMax)/(xMax-xMin));

#ifdef Blind
  if (xx < xbMin || xx > xbMax) return (nExpo+nCheb)*((nCheb/(nExpo+nCheb))*Chebichev+(1-(nCheb/(nCheb+nExpo)))*Exponential);
  else return 0;
#else
  return (nExpo+nCheb)*((nCheb/(nExpo+nCheb))*Chebichev+(1-(nCheb/(nCheb+nExpo)))*Exponential);
#endif
}

int BkgModelRoot()
{
  TFile *inFile = new TFile("BkgCombModelTree.root","OPEN");
  TTree *myTree = (TTree*)inFile->Get("CombModelData");
  Double_t massVal;
  myTree->SetBranchAddress("mass",&massVal);
  
  std::vector<Double_t> mass;

  for(int i = 0; i<myTree->GetEntries(); i++)
    {
      myTree->GetEntry(i);
#ifdef Blind
      if(massVal < xbMin || massVal > xbMax) mass.push_back(massVal);
#else
      mass.push_back(massVal);
#endif
    }

  std::cout << "I got #" << mass.size() << std::endl;

  TF1 *CombModel = new TF1("CombModel",CombinedModel,xMin,xMax,4);
  CombModel->SetParName(0,"aCheb");
  CombModel->SetParName(1,"aExpo");
  CombModel->SetParName(2,"nCheb");
  CombModel->SetParName(3,"nExpo");
  CombModel->SetParameter(0,-2.0164e-01); 
  CombModel->SetParError (0,2.0);
  CombModel->SetParLimits(0,-10.0,10.0);
  CombModel->SetParameter(1,-4.8184e-03); 
  CombModel->SetParError (1,2.0);    
  CombModel->SetParLimits(1,-10.0,10.0);
  CombModel->SetParameter(2,8000);        
  CombModel->SetParError (2,1000.0); 
  CombModel->SetParLimits(2,-10000.0,10000.0);
  CombModel->SetParameter(3,700); 
  CombModel->SetParError (3,2000.0); 
  CombModel->SetParLimits(3,-10000.0,10000.0);

  RootFitter(CombModel,xMin,xMax,true,mass);
  Plotter(mass,CombModel);
  
  return 0;
}