#include<iostream>
#include<vector>

#include "TFile.h"
#include "TTree.h"
#include "TBranch.h"
#include "TMath.h"
#include "TFitResult.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 Blinding

#ifdef Blinding
bool BlindingEnable = true;
Double_t Region1Min = 3;
Double_t Region1Max = 6;
#else
bool BlindingEnable = false;
Double_t Region1Min;
Double_t Region1Max;
#endif

Int_t Binning = 30;
Double_t RegionMin = 0;
Double_t RegionMax = 10;

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, std::vector<Double_t> RegionMin, std::vector<Double_t> RegionMax, bool ExtendedFit, std::vector<Double_t> FitData)
{
  ROOT::Fit::DataOptions DataOptions;
  ROOT::Fit::DataRange RangeOfData;
  for(uint iRange = 0; iRange<RegionMin.size(); iRange++) RangeOfData.AddRange(iRange,RegionMin[iRange],RegionMax[iRange]);
  ROOT::Fit::UnBinData *UnBinnedData = new ROOT::Fit::UnBinData(DataOptions,RangeOfData,FitData.size(),1,false);
  for(uint iData = 0; iData<FitData.size(); iData++) UnBinnedData->Add(FitData[iData]);
  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 iParam = 0; iParam<FitFunction->GetNpar(); iParam++)
    {
      FitParametersFromTF1[iParam]       = FitFunction->GetParameter(iParam);
      FitParametersErrorsFromTF1[iParam] = FitFunction->GetParError(iParam);
    }
  myFitConfig.SetParamsSettings(FitFunction->GetNpar(),FitParametersFromTF1,FitParametersErrorsFromTF1);

  for(int iParam = 0; iParam<FitFunction->GetNpar(); iParam++)
    {
      ROOT::Fit::ParameterSettings &FitParameterSettings = myFitConfig.ParSettings(iParam);
      Double_t FitParameterLowLimit, FitParameterUpLimit;
      FitFunction->GetParLimits(iParam,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(1000);
  myMinimizerOptions.Print();
  myFitConfig.SetMinimizerOptions(myMinimizerOptions);

  myFitter.Fit(*UnBinnedData,ExtendedFit);
  myFitter.GetMinimizer()->PrintResults();
  myFitter.CalculateHessErrors();
  const ROOT::Fit::FitResult &fitResult = myFitter.Result();
  for(uint iParam = 0; iParam<fitResult.NPar(); iParam++)
    {
      FitFunction->SetParameter(iParam,fitResult.Parameter(iParam));
      FitFunction->SetParError (iParam,fitResult.ParError(iParam));
    }
  std::cout << "Fit Status: " << fitResult.Status() << std::endl;
  fitResult.Print(std::cout);
  fitResult.PrintCovMatrix(std::cout);
  return fitResult;
}

Double_t NormExpo(Double_t aExpo, Double_t xx)
{
  if(BlindingEnable)
    {
      return (exp(aExpo*RegionMax)-exp(aExpo*RegionMin)-exp(aExpo*Region1Max)+exp(aExpo*Region1Min))/aExpo;
    }
  else return (exp(aExpo*RegionMax)-exp(aExpo*RegionMin))/aExpo;
}

Double_t Exponential(Double_t *x, Double_t *par)
{
  Double_t xx    = x[0];
  Double_t aExpo = par[0];
  Double_t nExpo = par[1];
  if(BlindingEnable)
    {
      if(xx > Region1Min && xx < Region1Max) return 0;
      else return nExpo*(1/NormExpo(aExpo,xx))*exp(aExpo*xx);
    }
  else return nExpo*(1/NormExpo(aExpo,xx))*exp(aExpo*xx);
}

int RootFitFit()
{
  TFile *inFile = new TFile("ToyExponentialTree.root","OPEN");
  TTree *myTree = (TTree*)inFile->Get("ExpoData");
  Double_t mass;
  myTree->SetBranchAddress("mass",&mass);
  Int_t nEntries = myTree->GetEntries();

  std::vector<Double_t> MassPoints;
  for(int iEntry = 0; iEntry<nEntries; iEntry++)
    {
      myTree->GetEntry(iEntry);
      if(BlindingEnable)
	{
	  if(mass < Region1Min || mass > Region1Max) MassPoints.push_back(mass);
	}
      else MassPoints.push_back(mass);
    }

  std::cout << "I got #" << MassPoints.size() << " data points" << std::endl;

  TF1 *Expo = new TF1("Expo",Exponential,RegionMin,RegionMax,2);
  Expo->SetParName(0,"aExpo");
  Expo->SetParName(1,"nExpo");
  Expo->SetParameter(0,-.5);   Expo->SetParLimits(0,-1.00000e+02,1.00000e+02);
  Expo->SetParameter(1,10000); Expo->SetParLimits(1,0,100000);
  Expo->SetParError(0,20.0);
  Expo->SetParError(1,5000);

  std::vector<Double_t> RegionMinVec, RegionMaxVec;
  //if(BlindingEnable)
  //  {
  //    RegionMinVec.push_back(RegionMin);
  //    RegionMaxVec.push_back(Region1Min);
  //    RegionMinVec.push_back(Region1Max);
  //    RegionMaxVec.push_back(RegionMax);
  //  }
  //else
  //  {
  RegionMinVec.push_back(RegionMin);
  RegionMaxVec.push_back(RegionMax);
      //  }
  
  ROOT::Fit::FitResult ResultObtained = RootFitter(Expo,RegionMinVec,RegionMaxVec,true,MassPoints);
  Plotter(MassPoints,Expo); 
  
  return 0;
}