#include "TF1.h"
#include "TH1.h"
#include "TVirtualFitter.h"
#include "TCanvas.h"
#include "TFitResult.h"
#include "Math/DistFunc.h"
#include "TFile.h"
#include "TList.h"
#include "TMath.h"

TH1D *cQHPHO;
TH1D *cGHPHO;
TH1D *cDHPHO;
TH1D *tQh;
TH1D *tGh;

const int nBins = 15;

// testing stuff
TH1D *cDHPHODDP;
int nfitpoints = 0;

//___________________________________________________________________________
Double_t TF1FromHistos(Double_t *x, Double_t *par) 
{

  tQh = (TH1D*) cQHPHO->Clone("tQh");
  tGh = (TH1D*) cGHPHO->Clone("tGh");
  Double_t alphap = par[0];

  //  printf("par:%f \n",par[0]);

  Double_t betap = 1.0-alphap;
  tQh->Scale(alphap);
  tGh->Scale(betap);

  // printf("Integral of tQh:%f, integral of tGh:%f \n",tQh->Integral(),tGh->Integral());

  Int_t qcbin = tQh->GetBin(x[0]);
  Double_t qcvalue = tQh->GetBinContent(qcbin);
  Int_t gcbin = tGh->GetBin(x[0]);
  Double_t gcvalue = tGh->GetBinContent(gcbin);
  Double_t pcvalue = qcvalue+gcvalue;

  tQh->Reset();
  tGh->Reset();

  return pcvalue;
}
//______________________________________________________________________________
void ModChiSqPhojet(Int_t &, Double_t *, Double_t &chisq, Double_t * p, Int_t)
{

  // get pointer to function and histogram
  assert(TVirtualFitter::GetFitter() );
  TF1 * func = (TF1*) (TVirtualFitter::GetFitter() )->GetUserFunc();
  TH1 * hist = (TH1*) (TVirtualFitter::GetFitter() )->GetObjectFit();
  assert(func && hist); 
  
  const Int_t nbins = nBins;  
  
  Double_t delta= 0.;
  chisq=0.0;
  nfitpoints = 0;
  double x[1];
  for (Int_t i=1;i<nbins; i++) 
    {
      x[0] = hist->GetBinCenter(i);
      double data = hist->GetBinContent(i); // data at current bin
      double dataerror = hist->GetBinError(i); // error on the data at the current bin
      double theory = func->EvalPar(x,p);
      if(dataerror==0)
	continue;
      delta = (data-theory)/dataerror;
      chisq += delta*delta;
      nfitpoints++;
    }
  return;
}
//___________________________________________________
void TestFitsChiSquare_2()
{

  TFile file("Fit_Histos.root");
  TList *list = (TList*)file.Get("list");
      
  cQHPHO = dynamic_cast<TH1D*> (list->FindObject("cQHPHO"));
  if(!cQHPHO) {printf("ERROR: cQHPHO not available"); return;}
  // normalization of q
  Double_t QNorm = cQHPHO->Integral();
  cQHPHO->Scale(1.0/QNorm);

  cGHPHO = dynamic_cast<TH1D*> (list->FindObject("cGHPHO"));
  if(!cGHPHO) {printf("ERROR: cGHPHO not available"); return;}
  // normalization of g
  Double_t GNorm = cGHPHO->Integral();
  cGHPHO->Scale(1.0/GNorm);

  cDHPHO = dynamic_cast<TH1D*> (list->FindObject("cDHPHO"));
  if(!cDHPHO) {printf("ERROR: cDHPHO not available"); return;}

  cDHPHODDP = (TH1D*) cDHPHO->Clone("cDHPHODPP");

  Double_t fitRange[2] = {1,15}; 
  TF1 *MyTF1 = new TF1("MyTF1",TF1FromHistos,fitRange[0],fitRange[1],1);
  MyTF1->SetParName(0,"Alpha");
  MyTF1->SetParLimits(0,0.01,1.0); // from 0.01 to 1
  MyTF1->SetParameter(0,0.10); 

  Int_t FitResultP = cDHPHO->Fit(MyTF1,"ENO","",fitRange[0],fitRange[1]);
  Double_t vpar0_P = MyTF1->GetParameter(0);
  Double_t verrpar0_P = MyTF1->GetParError(0);
  Double_t vchi2_P = MyTF1->GetChisquare();
  Int_t vNDF_P = MyTF1->GetNDF();

  printf("___________________________________\n");
  printf("Fit result:%i \n",FitResultP);
  printf("Alpha:%f \n",vpar0_P);
  printf("Error on alpha:%f \n",verrpar0_P);
  printf("Chisquare:%f, NDF:%i \n",vchi2_P,vNDF_P);
  printf("____________________________________\n");

  // Clone the q/g histos and scale with the parameter output
  // From my Chi2
  TH1D *hQTF1 = (TH1D*) cQHPHO->Clone("hQTF1");
  hQTF1->Scale(vpar0_P);
  hQTF1->SetLineColor(4);
  hQTF1->SetMarkerColor(4);
  TH1D *hGTF1 = (TH1D*) cGHPHO->Clone("hGTF1");
  hGTF1->Scale(1.-vpar0_P);
  hGTF1->SetLineColor(2);
  hGTF1->SetMarkerColor(2);
  TH1D *hPTF1 = (TH1D*) hQTF1->Clone("hPTF1");
  hPTF1->Add(hGTF1);
  hPTF1->SetLineColor(13);
  hPTF1->SetMarkerColor(13);

  TCanvas *cdisplay = new TCanvas("cdisplay","Comparison of fit methods",10,10,810,610);
  cdisplay->SetFillColor(0);
  cdisplay->cd(1);
  cDHPHO->DrawCopy();
  hPTF1->DrawCopy("SAME");
  // hQTF1->DrawCopy("SAME");
  // hGTF1->DrawCopy("SAME");

  printf("With my own chisquare... \n");
  //CREATE A NEW TF1
  TF1 *DDP = new TF1("DDP",TF1FromHistos,fitRange[0],fitRange[1],1);
  DDP->SetParName(0,"Alpha2");
  DDP->SetParLimits(0,0.0,1.0); // from 0 to 1

  //PERFORM chi2 FIT USING MY DEFINED FUNCTION
  TVirtualFitter::Fitter(cDHPHODDP)->SetFCN(ModChiSqPhojet); 
  TFitResultPtr r2 = cDHPHODDP->Fit(DDP,"USNO+","",fitRange[0],fitRange[1]); 
  // need to set chi2 value and ndf in fitted function
  DDP->SetChisquare(r2->MinFcnValue() );
  DDP->SetNDF(nfitpoints - r2->NFreeParameters() );
  r2->Print();
  std::cout << "NDF = " << DDP->GetNDF() << std::endl;
  std::cout << "Fit chi2 p value = " << ROOT::Math::chisquared_cdf_c(DDP->GetChisquare(), DDP->GetNDF() ) << std::endl;
  
  Double_t vpar0_P2 = DDP->GetParameter(0);
  Double_t verrpar0_P2 = DDP->GetParError(0);
  printf("Alpha2:%f \n",vpar0_P2);
  printf("Error2:%f \n",verrpar0_P2);

}