#include "RooSNLikelihood.hh"

// Class constructor for RooSNLikelihood.
RooSNLikelihood::RooSNLikelihood(const char *name, const char *title, RooAbsReal& _t0) :

  RooAbsReal(name, title),
  t0("t0", "Dependent", this, _t0)
{

  // Here, I want to set some of the ROOT stuff for later. I was originally thinking of using
  // a dedicated function like SetROOTPointers, but we shall just do it in the constructor. 

  std::cout << "                                                      " << std::endl;
  std::cout << "//////////////////////////////////////////////////////" << std::endl;
  std::cout << "/////RooSNLikelihood::RooSNLikelihood(Constructor)////" << std::endl;
  std::cout << "//////////////////////////////////////////////////////" << std::endl;
  std::cout << "                                                      " << std::endl;
}

// Class copy constructor of RooSNLikelihood.
RooSNLikelihood::RooSNLikelihood(const RooSNLikelihood& other, const char* name) :

  RooAbsReal(other, name),
  t0("t0",this, other.t0)
{
  std::cout << "                                                      " << std::endl;
  std::cout << "//////////////////////////////////////////////////////" << std::endl;
  std::cout << "//RooSNLikelihood::RooSNLikelihood(Copy constructor)//" << std::endl; 
  std::cout << "//////////////////////////////////////////////////////" << std::endl;
  std::cout << "                                                      " << std::endl;

  // Make sure to setup the files as well. This is necessary for functions like plotOn
  // in the RooAbsReal class since it creates a copy of the RooSNLikelihood object.
    
  std::cout << "                                                                    " << std::endl;
  std::cout << " Copying over varibales declarations to new RooSNLikelihood object! " << std::endl;
  std::cout << "                                                                    " << std::endl;

  SetDataFilename(other.fDataFilename);
  SetMeanFilename(other.fMeanFilename);
  ReadMeanEventEvolution();
  ReadSimulationData();

  std::cout << "                                                                    " << std::endl;
  std::cout << "Finished copying over varibales declarations to new RooSNLikelihood object! " << std::endl;
  std::cout << "                                                                    " << std::endl;

}


void RooSNLikelihood::SetDataFilename(const char* name){

  std::cout << "                                                      " << std::endl;
  std::cout << "//////////////////////////////////////////////////////" << std::endl;
  std::cout << "///////////RooSNLikelihood::SetDataFilename///////////" << std::endl;
  std::cout << "//////////////////////////////////////////////////////" << std::endl;
  std::cout << "                                                      " << std::endl;

  fDataFilename = name;

  std::cout << "For the data file we are using: " << name << std::endl;

}

void RooSNLikelihood::SetMeanFilename(const char* name){

  std::cout << "                                                      " << std::endl;
  std::cout << "//////////////////////////////////////////////////////" << std::endl;
  std::cout << "///////////RooSNLikelihood::SetMeanFilename///////////" << std::endl;
  std::cout << "//////////////////////////////////////////////////////" << std::endl;
  std::cout << "                                                      " << std::endl;

  std::cout << "For the mean file name we are using: " << name << std::endl;

  fMeanFilename = name;

}

void RooSNLikelihood::ReadMeanEventEvolution(){

  std::cout << "                                                      " << std::endl;
  std::cout << "//////////////////////////////////////////////////////" << std::endl;
  std::cout << "///////RooSNLikelihood::ReadMeanEventEvolution////////" << std::endl;
  std::cout << "//////////////////////////////////////////////////////" << std::endl;
  std::cout << "                                                      " << std::endl;


  // Okay, we need to open the file (with the TGraph) and
  // then get the TGraph and detach it. 
  fMeanROOTTFile = new TFile(fMeanFilename, "OPEN");

  TGraph* temp  = (TGraph*)fMeanROOTTFile->Get("mean-time");

  fMeanEventEvolution = (TGraph*)temp->Clone();

  std::cout << "The mean event evolution TGraph has been opened successfully!" << std::endl;

  // Close the file.
  fMeanROOTTFile->Close();

}

void RooSNLikelihood::ReadSimulationData(){

  std::cout << "                                                      " << std::endl;
  std::cout << "//////////////////////////////////////////////////////" << std::endl;
  std::cout << "/////////RooSNLikelihood::ReadSimulationData//////////" << std::endl;
  std::cout << "//////////////////////////////////////////////////////" << std::endl;
  std::cout << "                                                      " << std::endl;
  
  // Open the ROOT TFile where the data is located. 
  fDataROOTTFile = new TFile(fDataFilename, "OPEN");

  // Retrieve the TTree where the data is located. This will be 
  // under the background tree (where the underlying 1 Hz signal has
  // been added in).
  fDataROOTTTree = (TTree*)fDataROOTTFile->Get("background");

  // Now we will use the Draw() function to retrieve the event time
  // stamps and store them in the event vector.

  fDataROOTTTree->Draw("time", "", "goff");
  Int_t length = fDataROOTTTree->GetSelectedRows();

  for (Int_t v = 0; v < length; v++){
    fEvents.push_back( fDataROOTTTree->GetV1()[v] );
  }

  // Okay now we need to sort the vector.
  std::sort( fEvents.begin(), fEvents.end() );

  std::cout << "Data has been read correctly. There are " << fEvents.size() << " events in the burst" << std::endl;

  // Set the number of events and normalization.
  // These might be called by the user at a later date ( hint hint :) )
  SetNumberOfEvents();
  SetNormalization();

  // We can close the file now that we have the data.
  fDataROOTTFile->Close();

}

void RooSNLikelihood::SetNumberOfEvents(){
  fNumEvents = fEvents.size();
}

// Simple getter to pass the number of events. 
Int_t RooSNLikelihood::GetNumberOfEvents(){
  return fNumEvents;
}

// Set the normalization here. I would do it in evaluate()
// but I get errors when trying to do so. 
void RooSNLikelihood::SetNormalization(){
  fNormalization = fNumEvents / ( fMeanEventEvolution->Eval( fEvents.at(fNumEvents-1),0,0 ) );
}

// Simple getter to pass the normalization.
Double_t RooSNLikelihood::GetNormalization(){
  return fNormalization;
}

// Here we make use of the KSTest to measure the maximum deviation of the data set
// from the cummulative pdf. More information can be found at the following URL:
// https:://itl.nist.gov/div898/handbook/eda/section3/eda35g.htm
Double_t RooSNLikelihood::GetKSTest(){

  // Set the max deviation here.
  Double_t maxD = 0.0;

  // This will require that the dataset and the mean TGraph have 
  // been read in by the other functions in this class. In principle
  // this should always be the case but I will add a check later on to 
  // make sure the user doesn't forget.
  /*
  Int_t numEvtInBurst = fEvents.size();
  Double_t finalEvtTime = fEvents.at(numEvtInBurst-1);
  Double_t mean_at_tn = fMeanEventEvolution->Eval(finalEvtTime, 0, 0);
  Double_t norm = fNumEvents / ( fMeanEventEvolution->Eval( fEvents.at(fNumEvents-1),0,0 ) );  
  */

  // Return the statistic.
  return maxD;

}

Double_t RooSNLikelihood::evaluate() const
{

  // Okay, so here we will construct the likelihood function.
  // To do so we will make use of the mean event evolution (TGraph theMean).
  Double_t likelihood = 0.0;

  // First thing first let us get the number of events in the burst. 
  // Use this to sum over (as per the Likelihood section in the paper
  // I haven't finished...)  
  Int_t numEvtInBurst = fEvents.size();

  // Get the normalization constant. Ratio of events in burst compared to 
  // mean events at a time t_N.
  Double_t norm = (Double_t)numEvtInBurst / ( fMeanEventEvolution->Eval( fEvents.at(numEvtInBurst-1),0,0 ) );

  // Use this to sum over the events.
  for (Int_t j = 0; j < numEvtInBurst; j++){

    // Get the current event time.
    Double_t t = fEvents.at(j);

    // Evaluate the mean at this time minus t0.
    Double_t evalMean = fMeanEventEvolution->Eval(t-t0, 0, 0);

    // First likelihood term. 
    Double_t firstTerm = norm*evalMean;
    // Second likelihood term.
    Double_t secondTerm = -1*(j+1)*TMath::Log(evalMean*norm);
    // There is a third term that is constant for all points so we remove it. 
    
    likelihood += firstTerm + secondTerm;
  
  }
  

  return 2*likelihood;

}