// Program to implement "Errors on Errors" ("EOE") model as in
// G. Cowan, "Statistical Models with Uncertain Error Parameters",
// Eur. Phys. J. C (2019) 79:133; e-print arXiv:1809.05778.

// Glen Cowan
// RHUL Physics
// May 2019

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

#include <TMinuit.h>
#include <TPad.h>
#include <TFile.h>
#include <TF1.h>
#include <TMatrixD.h>
#include <TH2D.h>
#include <TRandom3.h>
#include <TMath.h>
#include <Math/Random.h>
#include <Math/GSLRndmEngines.h>

#include "Eoe.h"
#include "fitPar.h"

using namespace std;

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

  // Read input steering file
  string inputFileName;
  if ( argc == 1 ) {
    cout << "Enter name of input file: ";
    cin >> inputFileName;
    cin.ignore();         // clear buffer
  }
  else { 
    inputFileName = argv[1];
  }

// Open output histogram file and book histograms
  TFile* histFile = new TFile("eoe.root", "recreate");
  TH1D* h_muHat = new TH1D("h_muHat", "muHat", 100, 0., 20.); 
  TH1D* h_q = new TH1D("h_q", "q", 100, 0., 40.);

// Create random number generators (ran for uniform, rgsl for gamma)
  int seed = 1234567;
  TRandom3* ran = new TRandom3(seed);
  ROOT::Math::Random<ROOT::Math::GSLRngMT>* rgsl = new 
              ROOT::Math::Random<ROOT::Math::GSLRngMT>();
  rgsl->SetSeed(431155);

// Create Eoe object (reads steering file w/ data) and do fit
  Eoe eoe(inputFileName, ran, rgsl);
  int L = eoe.L();     // number of measurements
  int M = eoe.M();     // num param of interest (for ave. M=1)
  int N = eoe.N();     // num nuis param (= num control meas.)
  eoe.fit();           // calls fitPar

// Extract parameters and std. dev.
  vector<double> parVec = eoe.parVec();
  double muHat = parVec[0];
  vector<double> sigParVec = eoe.sigParVec();
  double sigmaMuHat = sigParVec[0];
  vector<double> thetaVec(N);
  for (int i=0; i<N; i++){
    thetaVec[i] = parVec[i+1];
  }
  vector<double> sigma_u_hh_vec = eoe.sigma_u_HatHatVec(thetaVec);

// Get MINOS interval for parameter of interest mu
  double epl = eoe.getEplus();
  double emi = eoe.getEminus();
  double epa = eoe.getEparab();
  double minosa = 0.5 * (epl + abs(emi));
  double minosLo = muHat - abs(emi);
  double minosUp = muHat + epl;

 // Get goodness of fit
  double q = eoe.q();               // EOE Eq. (61)
  // double qBCMC = eoe.qBCMC();       // q', MC-based Bartlett correction

// Output
  
  cout << "Input data:  y     dy" << endl;
  for (int i=0; i<L; i++){
    cout << eoe.y(i) << "\t" << eoe.sigma_y(i) << endl;
  }
  cout << "# control measurements, sys err, error on error (r):" << endl;
  for (int i=0; i<N; i++){
    cout << eoe.u(i) << "\t" << eoe.s(i) << "\t" << eoe.r(i) << endl;
  }
  cout << "#  Fit results:" << endl;
  cout << "# muHat = " << muHat << endl;
  cout << "# MINOS interval low edge     = " << minosLo << endl;
  cout << "# MINOS interval upper edge   = " << minosUp << endl;
  cout << "# minos-, minos+, minos ave, parab = " << emi 
       << "  " << epl << "  " << minosa << "  " << epa << endl;
  cout << "# q = " << q << endl;
  // cout << "# qBCMC = " << qBCMC << endl;
  cout << "# p-value from q = " << TMath::Prob(q, L-1) << endl;
  // cout << "# p-value from qBCMC = " << TMath::Prob(qBCMC, L-1) << endl;
  // cout << "# p-value from toy MC = " << eoe.pVal() << endl;
  cout << "#" << endl;

  cout << "# covariance matrix of sys biases U_{ij} = cov[bHat_i, bHat_j]..."
       << endl;
  for (int j=0; j<L; j++){
    cout << "#";
    for (int i=0; i<L; i++){
      cout << "\t" << eoe.cov_bHat(i,j);
    }
    cout << endl;
  }
  cout << "#" << endl;  

  cout << "#  meas. no.  y   sigma_y (stat)   sigma_bHat (sys)" << endl;
  cout << "#" << endl;  
  for (int i=0; i<L; i++){
    cout << i+1 << "\t" << eoe.y(i) << "\t" << eoe.sigma_y(i) << "\t" 
	 << eoe.sigma_bHat(i) << "\t" << endl;
  }
  
// Store any histograms and close up
  histFile->Write();
  histFile->Close();

  return 0;

}