#include <TH1D.h>
#include <TMath.h>
using namespace ROOT;

#include <RAT/NearAVAngular.hh>
#include <RAT/DB.hh>
#include <RAT/DS/Entry.hh>
#include <RAT/DS/FitVertex.hh>
#include <RAT/DS/ClassifierResult.hh>
#include <RAT/DU/Utility.hh>
#include <RAT/DU/PMTInfo.hh>
#include <RAT/DU/GroupVelocity.hh>
using namespace RAT;
using namespace RAT::Methods;

#include <string>
using namespace std;


void NearAVAngular::Initialise(const std::string& param)
{
  fIndex = param;
}

void NearAVAngular::BeginOfRun(DS::Run& run)
{
  DB *db = DB::Get();
  string index = fIndex;
  if (fIndex.empty()) {index = db->GetLink("GEO", "inner_av")->GetS("material");}

  // Default parameters are those of "labppo_scintillator"
  DBLinkPtr dbNAVLink = db->GetLink("FIT_NEAR_AV_ANGULAR", "labppo_scintillator");
  DBLinkGroup grp = db->GetLinkGroup("FIT_NEAR_AV_ANGULAR");
  DBLinkGroup::iterator it;

  // Check for the material-specific entry and overwrite the default values only if found
  for (it = grp.begin(); it != grp.end(); ++it)
  {
      if (index == it->first)
      {
          dbNAVLink = db->GetLink("FIT_NEAR_AV_ANGULAR", index);
          break;
      }
  }

  fNhitsWindow = dbNAVLink->GetIArray("nhits_windows");
  fWindowWidth = dbNAVLink->GetI("nhits_window_width");
  fFitCoeffs0 = dbNAVLink->GetDArray("fit_coefficients0");
  fFitCoeffs1 = dbNAVLink->GetDArray("fit_coefficients1");
  fFitCoeffs2 = dbNAVLink->GetDArray("fit_coefficients2");
  fRatioCuts = dbNAVLink->GetDArray("ratio_cuts");
  fRatioBins = dbNAVLink->GetDArray("ratio_bins");
  fBinWidth = dbNAVLink->GetD("ratio_bin_width");
  fNegativeErrors = dbNAVLink->GetDArray("negative_errors");
  fPositiveErrors = dbNAVLink->GetDArray("positive_errors");
  fThetaError = dbNAVLink->GetD("theta_error");
  fPhiError = dbNAVLink->GetD("phi_error");

  fAVRadius = db->GetLink("SOLID", "acrylic_vessel_outer")->GetD("r_sphere");
}

void NearAVAngular::EndOfRun(DS::Run& run)
{
  fNhitsWindow.clear();
  fFitCoeffs0.clear();
  fFitCoeffs1.clear();
  fFitCoeffs2.clear();
  fRatioCuts.clear();
  fRatioBins.clear();
  fNegativeErrors.clear();
  fPositiveErrors.clear();
}

DS::FitResult NearAVAngular::GetBestFit()
{
  fFitResult.Reset();

  // Fitter cannot fit events with Nhits outside the coordinated range
  if (fPMTData.empty())
  {
    throw MethodFailError("NearAVAngular::GetBestFit() - no Nhits in event - cannot fit position");
  }
  if (fPMTData.size() >= (fNhitsWindow.back() + fWindowWidth))
  {
    throw MethodFailError("NearAVAngular::GetBestFit() - Nhits is beyond the upper limit of the fitter");
  }

  // Get the Classifier Result, i.e. the ratio value
  double ratio = 0.0;
  try
  {
    DS::ClassifierResult nearAVClassification = fEvent->GetClassifierResult("nearAVAngular");
    if(nearAVClassification.GetValid() == false)
    {
      throw MethodFailError("NearAVAngular::GetBestFit() - nearAVAngular classification is invalid");
    }
    ratio = nearAVClassification.GetClassification("ratio");
  }
  catch(std::runtime_error& e)
  {
    throw MethodFailError("NearAVAngular::GetBestFit() - no nearAVAngular Classification - please run the NearAVAngular Classifier");
  }

  // Set all coordinated variables based on the Nhits window and ratio bin that the event is located in
  int eventNhitsWindow = -1;
  for (unsigned int i = 0; i < fNhitsWindow.size(); i++)
  {
    if ((fPMTData.size() >= fNhitsWindow[i]) && (fPMTData.size() < (fNhitsWindow[i] + fWindowWidth)))
    {
      eventNhitsWindow = i;
      break;
    }
  }
  int eventRatioBin = -1;
  for (unsigned int j = 0; j < fRatioBins.size(); j++)
  {
    if ((ratio >= fRatioBins[j]) && (ratio < (fRatioBins[j] + fBinWidth)))
    {
      eventRatioBin = j;
      break;
    }
  }
  if (eventRatioBin == -1)
  {
    throw MethodFailError("NearAVAngular::GetBestFit() - event ratio is beyond the upper limit of the fitter");
  }

  double a = fFitCoeffs2[eventNhitsWindow];
  double b = fFitCoeffs1[eventNhitsWindow];
  double c = fFitCoeffs0[eventNhitsWindow];
  double ratioCut = fRatioCuts[eventNhitsWindow];
  double negativeRadialError = fNegativeErrors[(fRatioBins.size() * eventNhitsWindow) + eventRatioBin];
  double positiveRadialError = fPositiveErrors[(fRatioBins.size() * eventNhitsWindow) + eventRatioBin];

  // If the event has ratio value not consistent with the near-AV region, do not fit and end fitter algorithm here
  if (ratio > ratioCut)
  {
    throw MethodFailError("NearAVAngular::GetBestFit() - NearAVAngular Classifier result is inconsistent with the near-AV region - use a different position fitter");
  }

  // Calculate the event radius from the classifier's ratio value using a pol2 function and the coefficients
  // Only the negative solution of the quadratic function is used
  double fitDistanceFromCentre = ((-1.0 * b) - sqrt((b * b) - (4 * a * (c - ratio)))) / (2 * a);

  // If reconstructed radius is outside the AV, set it to be on the AV shell
  if (fitDistanceFromCentre > fAVRadius)
  {
    fitDistanceFromCentre = fAVRadius;
  }

  // Calculate the reconstructed event axis w.r.t. to AV centre, and plot PMT hit times for later time reconstruction
  TVector3 eventAxis (0.0, 0.0, 0.0);
  TH1D eventTimes ("Temp", "Temp", 500, 0.0, 500.0);
  eventTimes.SetDirectory(0);

  const DU::PMTInfo& pmtInfo = DU::Utility::Get()->GetPMTInfo();
  for (vector<FitterPMT>::const_iterator iPMT = fPMTData.begin(); iPMT != fPMTData.end(); ++iPMT)
  {
    const TVector3 pmtVector = pmtInfo.GetPosition(iPMT->GetID());
    eventAxis += pmtVector;
    eventTimes.Fill(iPMT->GetTime());
  }

  // Fill in the full vertex information
  eventAxis = eventAxis.Unit();
  const TVector3 fitPosition = eventAxis * fitDistanceFromCentre;
  const double firstTime = RepositionHistogram(&eventTimes);
  const double transitTime = DU::Utility::Get()->GetGroupVelocity().CalcByDistance((fAVRadius - fitPosition.Mag()), 55.0, 2340.0);

  DS::FitVertex fitVertex;
  fitVertex.SetTime(firstTime - transitTime);
  fitVertex.SetTimeErrors(10.0);
  fitVertex.SetPosition(fitPosition);
  fitVertex.SetPositivePositionError(ROOT::Math::Polar3DVectorF(positiveRadialError, fThetaError, fPhiError));
  fitVertex.SetNegativePositionError(ROOT::Math::Polar3DVectorF(negativeRadialError, fThetaError, fPhiError));

  fFitResult.SetVertex(0, fitVertex);
  return fFitResult;
}

double NearAVAngular::RepositionHistogram(TH1D* histo)
{
  // Expect ~3 noise hits per event, thus first bin with more than 3 signals the event start
  const int maxBin = histo->FindFirstBinAbove(3.0);
  const int offsetBin = maxBin - histo->GetXaxis()->FindBin(250.0);

  TH1D reposHisto ("reposHisto", "reposHisto", 500, 0.0, 500.0);
  reposHisto.SetDirectory(0);

  for(int i = 1; i < 501; i++)  // Ignore the under and overflow bins
  {
    const int oldBin = i + offsetBin;
    if ((oldBin < 1) || (oldBin > 500))
    {
      reposHisto.SetBinContent(i, 0.0);
    }
    else
    {
      reposHisto.SetBinContent(i, histo->GetBinContent(oldBin));
    }
  }
  *histo = reposHisto;
  return histo->GetBinCenter(maxBin);
}