#include <iostream>
#include <fstream>
#include <sstream>
#include <vector>
#include <cmath>
#include <cstdlib>

#include "PLFit.h"
#include "PLUtil.h"
#include "fitPar.h"

#include <TApplication.h>
#include <TCanvas.h>
#include <TPad.h>
#include <TMath.h>
#include <TH1D.h>
#include <TF1.h>
#include <TFile.h>
#include <TMinuit.h>
#include <TRandom3.h>
#include <TPostScript.h>
#include <sys/stat.h>
#include <sys/types.h>

using namespace std;

PLFit* plf;      // must be global to communicate with fcn (see below)

int main(int argc, char *argv[]) {

// Important parameter values

  double nuTot = 10000;       // mean total number of entries
  int nBin = 50;
  double xMin = 0.;
  double xMax = 1.;
  int nparFitMax = 6;

// Set up histograms, function objects, etc.

  TFile* file = new TFile("PLFit.root", "recreate");
  vector<TH1D*> chi2Hist;
  for (int i=0; i<=nparFitMax; i++){
    stringstream ss1, ss2;
    string istr1, istr2;
    ss1 << max(i-1,0);
    ss1 >> istr1;
    ss2 << i;
    ss2 >> istr2;
    string histName = "h_chi2_" + istr1 + istr2;
    string histTitle = "chi2_" + istr1 + istr2;
    TH1D* h = new TH1D(histName.c_str(), histTitle.c_str(), 400, 0., 400.);
    chi2Hist.push_back(h);
  }

  vector<TH1D*> deltaChi2Hist;
  for (int i=0; i<nparFitMax; i++){
    stringstream ss1, ss2;
    string istr1, istr2;
    ss1 << i;
    ss1 >> istr1;
    ss2 << i+1;
    ss2 >> istr2;
    string histName = "h_deltaChi2_" + istr1 + istr2;
    string histTitle = "deltaCchi2_" + istr1 + istr2;
    TH1D* h = new TH1D(histName.c_str(), histTitle.c_str(), 400, 0., 400.);
    deltaChi2Hist.push_back(h);
  }

  int nBern = 2;
  int npar = nBern+1;
  vector<double> trueBeta(npar);
  trueBeta[0] = 1.;
  trueBeta[1] = 0.5;
  trueBeta[2] = 1.5;

  int fitOption = 0;
  vector<double> truePdfPar;

  plf = new PLFit(xMin, xMax, nBin, truePdfPar, nuTot, trueBeta, fitOption);

  vector<double> modelBeta(npar);
  for (int i=0; i<npar; i++){
    modelBeta[i] = 1.;
  }

  int seed = 12345;
  TRandom3* ran = new TRandom3(seed);

  int numExp = 1;
  for (int j=0; j<numExp; j++){

    vector<double> pdfPar;
    pdfPar = truePdfPar;
    plf->setModelHist(pdfPar, modelBeta);
    plf->generateData(ran);

// test hypothetical distribution using the data set

    vector<double> chi2Array;
    for (int nparFit=0; nparFit<=nparFitMax; nparFit++){

// fit parameters and compute profile likelihood ratio 

      vector<double> parVec;
      vector<vector<double> > V;      // covariance matrix of fitted parameters
      for (int i=0; i<nparFit; i++){
        parVec.push_back(1.);
      }
      double chi2Min;
      fitPar(plf, nparFit, chi2Min, parVec, V);
      double pval = TMath::Prob(chi2Min, nBin - nparFit);
      chi2Array.push_back(chi2Min);
      // cout << nparFit << "  " << chi2Min << "  " << pval << endl;

      vector<vector<double> > rho;
      for (int k=0; k<nparFit; k++){
        vector<double> rowK;
        for (int l=0; l<nparFit; l++){
          double r = V[k][l] / sqrt(V[k][k]*V[l][l]);
          rowK.push_back(r);
        }
        rho.push_back(rowK);
      }

      //  if ( nparFit == 3 ) {
      //    for (int k=0; k<3; k++){
      //      double dpar = sqrt(V[k][k]);
      //    printf("%12d %12.6f %12.6f \n", k, parVec[k], dpar);
      //    }
      //    cout << endl;
      //    for (int k=0; k<3; k++){
      //      printf("%12.6f %12.6f %12.6f \n", V[k][0], V[k][1], V[k][2]);
      //    }
      //    cout << endl;
      //    for (int k=0; k<3; k++){
      //      printf("%12.6f %12.6f %12.6f \n", rho[k][0], rho[k][1], rho[k][2]);
      //    }
      //    cout << endl;
      //  }

      double y, sigma;
      plf->integrateModel(30, 39, y, sigma, parVec, V);
      // cout << nparFit << "  "  << y << "  " << sigma << endl;

      if ( j%max(numExp/10,1) == 0 && j > 0 && nparFit == nparFitMax ) {
        cout << "generated experiment " << j << endl;
      }

      // cout << nparFit << "  " << chi2Min << "  " << pval << endl;

    }          // end loop over number of fitted parameters


    for (int n=0; n<=nparFitMax; n++){
      chi2Hist[n]->Fill(chi2Array[n]);
      int ndof = plf->nBin() - n;
      double pval0 = TMath::Prob(chi2Array[n], ndof);
      if ( n < nparFitMax ) {
        double deltaChi2 = chi2Array[n] - chi2Array[n+1];
        deltaChi2Hist[n]->Fill(deltaChi2);
        double pval01 = TMath::Prob(deltaChi2, 1);
        cout << n << "  &  " << chi2Array[n] << "  &  " << pval0 << 
	  "  &  " << deltaChi2 << "  &  " << pval01 << "  \\\\ " << endl;
      }
      else {
        cout << n << "  " << chi2Array[n] << "  " << pval0 << endl;

      }
    }

  }           // end loop over experiments

// store histograms and close up

  file->Write();
  file->Close();

  return 0;

}