#include <iostream>
#include <string>
#include <vector>

#include "PLFit.h"

#include <TH1D.h>
#include <TFile.h>
#include <TCanvas.h>
#include <TPad.h>
#include <TROOT.h>
#include <TStyle.h>
#include <TEnv.h>
#include <TLegend.h>
#include <TPostScript.h>
#include <TLine.h>
#include <TMath.h>
#include <TGraph.h>

using namespace std;

PLFit::PLFit(double xMin, double xMax, int nBin, vector<double> truePdfPar,
	     double nuTot, vector<double> trueScalePar, int fitOption){

  m_xMin = xMin;
  m_xMax = xMax;
  m_nBin = nBin;
  m_fitOption = fitOption;
  m_truePdfPar = truePdfPar;
  m_trueScalePar = trueScalePar;
  m_nuTot = nuTot;

  m_dataHist = new TH1D("m_dataHist", "data histogram", 
                        m_nBin, m_xMin, m_xMax);
  m_trueHist = new TH1D("m_trueHist", "true histogram", 
                       m_nBin, m_xMin, m_xMax);
  m_modelHist = new TH1D("m_modelHist", "model histogram", 
                       m_nBin, m_xMin, m_xMax);

// set up true distribution (unknown in a real analysis)

  for (int i=1; i<=nBin; i++){
    double xMid = m_trueHist->GetBinCenter(i);
    double deltaX = m_trueHist->GetBinWidth(i);
    double s = this->scaleFactor(xMid, trueScalePar);
    double pdf = this->modelPdf(xMid, truePdfPar);
    double nu = nuTot*pdf*deltaX*s;
    m_trueHist->SetBinContent(i, nu);
  }

  //  this->checkHist(m_trueHist);

}

void PLFit::generateData(TRandom* ran){
  m_dataHist->Reset();
  for (int i=1; i<=m_nBin; i++){
    double nu = m_trueHist->GetBinContent(i);
    int n = ran->Poisson(nu);
    double nd = static_cast<double>(n);
    m_dataHist->SetBinContent(i, nd);
  }
}

void PLFit::setModelHist(vector<double> pdfPar, vector<double> scalePar){
  m_modelHist->Reset();
  for (int i=1; i<=m_nBin; i++){
    double xMid = m_modelHist->GetBinCenter(i);
    double deltaX = m_modelHist->GetBinWidth(i);
    double s = this->scaleFactor(xMid, scalePar);
    double pdf = this->modelPdf(xMid, pdfPar);
    double nu = m_nuTot*pdf*deltaX*s;
    m_modelHist->SetBinContent(i, nu);
  }
}

double PLFit::modelPdf(double x, vector<double> par){
  double f = 2.*(1.-x);
  return f;
}

double PLFit::scaleFactor(double x, vector<double> beta){
  double s;
  if ( beta.size() == 0 ) {
    s = 1.;
  }
  else {
    s = 0.;
    int nBern = beta.size() - 1;
    for (int k=0; k<=nBern; k++){
      s += beta[k] * PLUtil::BernsteinPoly(x, k, nBern);
    }
  }
  return s;
}

void PLFit::integrateModel(int a, int b, double& y, double& sigma, 
		     vector<double> thetaHat, vector<vector<double> > V){

// integrates the model from bins a to b (inclusive) and computes the error
// using error propagation with numerically estimated derivatives.

  vector<double> pdfPar;          // for now not used
  this->setModelHist(pdfPar, thetaHat);
  y = 0.;
  for(int i=a; i<=b; i++){
    y += modelHist()->GetBinContent(i);
  }

  const double reasonableFactor = 0.2;
  int npar = thetaHat.size();
  vector<double> delta(npar);
  vector<double> deriv(npar);
  for (int i=0; i<npar; i++){
    delta[i] = reasonableFactor * sqrt(V[i][i]);  // step size for derivative
  }

  for (int i=0; i<npar; i++){
    vector<double> theta = thetaHat;
    theta[i] += delta[i];
    this->setModelHist(pdfPar, theta);
    double yplus = 0.;
    for(int j=a; j<=b; j++){
      yplus += modelHist()->GetBinContent(j);
    }
    theta[i] -= 2.*delta[i];
    this->setModelHist(pdfPar, theta);
    double yminus = 0.;
    for(int j=a; j<=b; j++){
      yminus += modelHist()->GetBinContent(j);
    }
    deriv[i] = (yplus - yminus)/(2.*delta[i]);
  }

  double sigma2 = 0.;
  for (int i=0; i<npar; i++){
    for (int j=0; j<npar; j++){
      sigma2 += deriv[i]*deriv[j]*V[i][j];
    }
  }
  sigma = sqrt(sigma2);

}

void PLFit::checkHist(TH1D* h){
  
  int n = h->GetNbinsX();
  double xMin = h->GetXaxis()->GetXmin();
  double xMax = h->GetXaxis()->GetXmax();
  cout << "n    = " << n << endl;
  cout << "xMin = " << xMin << endl;
  cout << "xMax = " << xMax << endl;
  
  // recall for TH1 objects, bin 0 is underflow, 1 to n give the histogram,
  // n+1 is overflow;
  
  for (int i=0; i<n; i++){
    cout << i+1 << "  " << h->GetBinContent(i+1) << endl;
  }
  
}

void PLFit::plotSetup(TCanvas* canvas){
  
  canvas->SetFillColor(0);
  canvas->SetBorderMode(0);
  
  gROOT->SetStyle("Plain");
  canvas->UseCurrentStyle();
  gROOT->ForceStyle();
  
}