#include <RAT/Log.hh>
#include <RAT/DU/ReactorNuOsc.hh>
#include <RAT/LLAUtil.hh>
#include <CLHEP/Random/Randomize.h>         // HepUniformRand()
#include <TMath.h>

using namespace RAT;
using namespace RAT::DU;

using namespace std;

#include <string>
#include <sstream>                          // istringstream
#include <cmath>
#include <climits>                          // INT_MAX
#include <locale>                           // std::locale, std::tolower

void
ReactorNuOsc::BeginOfRun()
{
  fOscProb = 0.0;
  fLink = DB::Get()->GetLink("OSCILLATIONS", "PDG2016");
  fDm21Sqr = fLink->GetD("deltamsqr21");
  fSSqrTheta12 = fLink->GetD("sinsqrtheta12");
  fSSqrTheta13 = fLink->GetD("sinsqrtheta13");

  info << "[ReactorNuOsc]::BeginOfRun: ";
    info << "Values of oscillation parameters:" << newline;
  info << "[ReactorNuOsc]  Delta m_12^2 [eV^2] : ";
    info << "dm21Sqr     = " << fDm21Sqr << newline;
  info << "[ReactorNuOsc]  sin^2(theta12)      : ";
    info << "sSqrTheta12 = " << fSSqrTheta12 << newline;
  info << "[ReactorNuOsc]  sin^2(theta13)      : ";
    info << "sSqrTheta13 = " << fSSqrTheta13 << newline;

  fSSqr2Theta12 = pow(sin(2.0 * asin(sqrt(fSSqrTheta12))), 2.0);
  fS4 = pow(fSSqrTheta13, 2.0);
  fC4 = pow(1.0-fSSqrTheta13, 2.0);
}


void ReactorNuOsc::OscProbability(const double &nuE,
                                      const double &baseline)
{
  // Intermediate steps in probability calculation
  const double scale = 1.267e3; // for nuE in [MeV] and baseline in [km]
  const double sSqrDmBE = pow(sin(scale * fDm21Sqr * baseline / nuE), 2.0);

  //Probability to oscillate (from Satoko Asahi [DocDB-3763]):
  fOscProb = 1.0 - (fC4 * (1.0 - fSSqr2Theta12 * sSqrDmBE) + fS4);
#ifdef RATDEBUG
  debug << "[ReactorNuOsc]::OscProbability: ";
    debug << "Oscillation probability = " << fOscProb << newline;
#endif
}

bool ReactorNuOsc::Sample() {
  const double rndm = CLHEP::HepUniformRand();

  // If the antineutrino has NOT oscillated, keeps the event
  if(fOscProb < rndm)
    return true;
  else
    return false;
}

bool ReactorNuOsc::ReactorOsc(const double nuKE, std::string coreData) {
  // Ensure correct string split when reactor names have spaces
  // (e.g. "ST. LAURENT B")
  std::string::size_type fnd = coreData.find_last_of(" ");
  unsigned int cNo=999;
  std::string rName="";
  if (fnd != std::string::npos) {
    rName = coreData.substr(0,fnd);
    std::istringstream(coreData.substr(fnd+1)) >> cNo;
#ifdef RATDEBUG
  debug << "[ReactorNuOsc]::DSEvent: Reactor information" << newline;
  debug << "[ReactorNuOsc]  Name = " << rName << newline;
  debug << "[ReactorNuOsc]  Core = " << cNo << newline;
#endif
  } else {
    Log::Die("[ReactorNuOsc]::DSEvent: Unknown metaInfo string format: "
      + coreData);
  }

  // Get reactor coordinates
  fLink = DB::Get()->GetLink("REACTOR", rName);
  fLatitude  = fLink->GetDArray("latitude");
  fLongitute = fLink->GetDArray("longitude");
  fAltitude  = fLink->GetDArray("altitude");
#ifdef RATDEBUG
  debug << "[ReactorNuOsc]::DSEvent: Reactor coordinates" << newline;
  debug << "[ReactorNuOsc]  latitude  = " << fLatitude[cNo] << newline;
  debug << "[ReactorNuOsc]  longitude = " << fLongitute[cNo] << newline;
  debug << "[ReactorNuOsc]  altitude  = " << fAltitude[cNo] << newline;
#endif

  // Calculate the distance [km] from the core to SNO+
  const double distance = GetReactorDistanceLLA(fLongitute[cNo], fLatitude[cNo], fAltitude[cNo]);
#ifdef RATDEBUG
  debug << "[ReactorNuOsc]::DSEvent: SNO+ to " << rName << newline;
  debug << "[ReactorNuOsc]  distance  = " << distance << newline;
#endif

  // If the antineutrino has NOT oscillated, keeps the event
  OscProbability(nuKE, distance);
  if(Sample() == true)
    return true;
  else
    return false;
}