#include <TRandom.h>
#include <TCanvas.h>
#include <TMath.h>
#include <TH1D.h>
#include <TF1.h>
#include <TFile.h>
#include <iostream>

using namespace std;



// Global variables
Double_t fXMax;
Double_t fXMin;
Double_t fBinWidth;



// Normalized pdf
Double_t PdfExponential(Double_t *x, Double_t *par)
{
  Double_t slope         = par[0];
  Double_t normalization = (TMath::Exp(fXMax * slope) - TMath::Exp(fXMin * slope)) / slope;

  return TMath::Exp(x[0] * slope) / normalization;
}



// Project pdf on histo
Double_t pdf_proj(Double_t *x, Double_t *par)
{
  return fBinWidth  * par[0] * PdfExponential(x, &par[1]);
}



// This function
void CompareFitsSimple()
{
  // Set seed with clock (different values each time)
  gRandom->SetSeed(0);

  // Set parameters
  fXMin       = 1.2;
  fXMax       = 2.4;

  Int_t    iT = 100;
  Int_t    B  = 4000;
  Double_t m;


  
  // Set output
  TFile *outFile = new TFile("CompareFitsSimple.root", "RECREATE");

  TH1D  *BChi2   = new TH1D("BChi2",   "Background, Chi2 fit",       100, (Double_t)B - 7. * TMath::Sqrt((Double_t)B) - 1., (Double_t)B + 7. * TMath::Sqrt((Double_t)B) + 1.);
  TH1D  *mChi2   = new TH1D("mChi2",   "Slope, Chi2 fit",            100, -1.65, -1.35);
  TH1D  *BChi2e  = new TH1D("BChi2e",  "Background err, Chi2 fit",   100, TMath::Sqrt((Double_t)B) - TMath::Sqrt(TMath::Sqrt((Double_t)B)), TMath::Sqrt((Double_t)B) + TMath::Sqrt(TMath::Sqrt((Double_t)B)));
  TH1D  *mChi2e  = new TH1D("mChi2e",  "Slope err, Chi2 fit",        100,  0.,    0.1);
  
  TH1D  *BChi2L  = new TH1D("BChi2L",  "Background, Chi2 L fit",     100, (Double_t)B - 1. * TMath::Sqrt((Double_t)B) - 1., (Double_t)B + 1. * TMath::Sqrt((Double_t)B) + 1.);
  TH1D  *mChi2L  = new TH1D("mChi2L",  "Slope, Chi2 L fit",          100, -1.65, -1.35);
  TH1D  *BChi2Le = new TH1D("BChi2Le", "Background err, Chi2 L fit", 100, TMath::Sqrt((Double_t)B) - TMath::Sqrt(TMath::Sqrt((Double_t)B)), TMath::Sqrt((Double_t)B) + TMath::Sqrt(TMath::Sqrt((Double_t)B)));
  TH1D  *mChi2Le = new TH1D("mChi2Le", "Slope err, Chi2 L fit",      100,  0.,    0.1);
  TH1D  *BChi2D  = new TH1D("BChi2D",  "Background Chi2 - Chi2L",    100,             - 3. * TMath::Sqrt((Double_t)B) - 1.,               3. * TMath::Sqrt((Double_t)B) + 1.);
  TH1D  *mChi2D  = new TH1D("mChi2D",  "Slope Chi2 - Chi2L",         100, -0.1,   0.1);
  
  
  
  // Functions used for generation
  TH1D       *genEv         = new TH1D("GenEv", "Generated events;m_{K^{-}#pi^{+}};Events per 4 MeV", 300, fXMin, fXMax);
  fBinWidth                 = genEv->GetBinWidth(1);

  TF1        *genBackground = new TF1("GenBkg", PdfExponential, fXMin, fXMax, 1);
  genBackground->SetParameter(0, -1.5);
  genBackground->SetNpx(10000);


  
  TF1        *dummyF        = new TF1("DummyF", pdf_proj, fXMin, fXMax, 2);
  dummyF->SetParameters(2., -1.5);
  dummyF->SetNpx(10000);

  

  TCanvas    *checkCanvas   = new TCanvas("CheckCanvas", "Check canvas", 1000, 0, 900, 900);

  Double_t   chi2B;
  Double_t   chi2m;
  cout << endl;
  for (Int_t iTime = 0; iTime < iT; ++iTime)
    {
      // Generate the current random sample
      genEv->Reset("M");
      
      // B
      for (Int_t i = 0; i < B; ++i)
      	{
      	  m = genBackground->GetRandom(fXMin, fXMax);
      	  genEv->Fill(m);
      	}
      
      dummyF->SetRange(1.2, 2.4);
      
      
      
      ///////
      // Chi2
      ///////
      dummyF->SetParameters((Double_t)B * 0.9, - 1.5);
      genEv->Fit(dummyF, "R", "PE", 1.2, 2.4);
      
      BChi2 ->Fill(dummyF->GetParameter(0));
      mChi2 ->Fill(dummyF->GetParameter(1));
      BChi2e->Fill(dummyF->GetParError( 0));
      mChi2e->Fill(dummyF->GetParError( 1));
      
      chi2B  = dummyF->GetParameter(0);
      chi2m  = dummyF->GetParameter(1);
      
      
      
      /////////
      // Chi2 L
      /////////
      dummyF->SetParameters((Double_t)B * 0.9, - 1.5);
      genEv->Fit(dummyF, "RL", "PE", 1.2, 2.4);
      
      BChi2L ->Fill(dummyF->GetParameter(0));
      mChi2L ->Fill(dummyF->GetParameter(1));
      BChi2Le->Fill(dummyF->GetParError( 0));
      mChi2Le->Fill(dummyF->GetParError( 1));
      BChi2D ->Fill(dummyF->GetParameter(0) - chi2B);
      mChi2D ->Fill(dummyF->GetParameter(1) - chi2m);
      
      
      
      checkCanvas->Update();
    }
  delete genEv;
  outFile->cd();
  outFile->Write();

  return;
}