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

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()
{
  fLink = DB::Get()->GetLink("OSCILLATIONS");

  fDmSqr21 = fLink->GetD("deltamsqr21");
  fDmSqr32 = fLink->GetD("deltamsqr32");
  fSSqrTheta12 = fLink->GetD("sinsqrtheta12");
  fSSqrTheta13 = fLink->GetD("sinsqrtheta13");

#ifdef RATDEBUG
  debug << "[ReactorNuOsc]::BeginOfRun: ";
  debug << "Values of oscillation parameters:" << newline;
  debug << "[ReactorNuOsc]  Delta m^2_21 [eV^2] : ";
  debug << "DmSqr21     = " << fDmSqr21 << newline;
  debug << "[ReactorNuOsc]  Delta m^2_32 [eV^2] : ";
  debug << "DmSqr32     = " << fDmSqr32 << newline;
  debug << "[ReactorNuOsc]  sin^2(theta12)      : ";
  debug << "sSqrTheta12 = " << fSSqrTheta12 << newline;
  debug << "[ReactorNuOsc]  sin^2(theta13)      : ";
  debug << "sSqrTheta13 = " << fSSqrTheta13 << newline;
#endif

  fSSqr2Theta12 = 4.0*fSSqrTheta12*(1.0-fSSqrTheta12);
  fSSqr2Theta13 = 4.0*fSSqrTheta13*(1.0-fSSqrTheta13);
  fC4 = pow(1.0-fSSqrTheta13, 2.0);
  fOscProb = 0.0;
}


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

  // Electron (anti)neutrino oscillation probability (1 - survival probability)
  fOscProb = fC4 * fSSqr2Theta12 * sSqrDm21BE
           + fSSqr2Theta13 * ( (1.0-fSSqrTheta12)*sSqrDm31BE + fSSqrTheta12*sSqrDm32BE );

#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;
}