/*
 * IBDCrossSec.cc
 *
 *  Created on: Aug 23, 2015
 *      Author: nbarros
 */

#include <RAT/IBDCrossSec.hh>
#include <RAT/IBDCrossSecMessenger.hh>
#include <RAT/Log.hh>

// -- CLHEP includes
#include <CLHEP/Units/PhysicalConstants.h>

//-- STL includes
#include <sstream>
#include <cstring>
#include <cstdlib>
#include <cmath>

// -- ROOT includes
#include <TFile.h>
#include <TH2F.h>

using CLHEP::MeV;
using CLHEP::proton_mass_c2;
using CLHEP::neutron_mass_c2;
using CLHEP::electron_mass_c2;
using CLHEP::pi;
using CLHEP::hbarc;

namespace RAT {

const map<IBDCrossSec::Strategy,string> IBDCrossSec::_str_map =  IBDCrossSec::init_map();

IBDCrossSec::IBDCrossSec() : _strategy(Vogel),_input_file(0), _xs_table(0) {

  fMessenger = new IBDCrossSecMessenger(this);
}

IBDCrossSec::~IBDCrossSec() {
  if (_input_file) _input_file->Close();
  delete _input_file;
  delete fMessenger;
}

void IBDCrossSec::SetStrategy(Strategy strat) {
  switch (strat) {
  case VogelTable:
    LoadTable(strat);
    _strategy = strat;
    break;
  case Vogel:
    _strategy = strat;
    break;
  default:
    Log::Die("Unknown strategy. Be sure that a valid value is selected.",1);
  }
}

double IBDCrossSec::dSigmadTheta(const double Enu, const double cosThetaLab) {
  // This gives room to implement different calculation strategies to improve performance
  // and/or precision
  switch(_strategy) {
  // Strategy 0 is the default
  case Vogel:
    return dSigmadThetaVogel(Enu,cosThetaLab); break;
  case VogelTable:
  default:
    if (!_xs_table) LoadTable(_strategy);
    return dSigmadThetaVogelTable(Enu,cosThetaLab);
  }
}

double IBDCrossSec::dSigmadThetaVogel(const double Enu, const double cosThetaLab) {
  double sigma = 0.0;
  // Cross section constants.  Some for overall scale are just
  // to allow absolute comparison to published article
  // (Vogel & Beacom, 1999)
  // Physical Review D, 60(5), 053003.
  // http://doi.org/10.1103/PhysRevD.60.053003
  static const double kGFERMI = 1.16639e-11 / MeV / MeV;
  static const double cosThetaC = (0.9741+0.9756)/2;
  static const double gDELTA = neutron_mass_c2 - proton_mass_c2;

  // Radiative correction constant
  static const double RadCor = 0.024;

  // check for threshold
  static const double EminBeta = (
         (neutron_mass_c2+electron_mass_c2)*
         (neutron_mass_c2+electron_mass_c2)
         -proton_mass_c2*proton_mass_c2
         )  /2.0/proton_mass_c2;

  // If energy is below threshold just return 0
  if(Enu<EminBeta) return 0.0;

  // overall scale
  static const double Sigma0 = kGFERMI*kGFERMI*cosThetaC*cosThetaC/pi*(1+RadCor);

  // couplings
  static const double f = 1.00;
  static const double f2 = 3.706;
  static const double g = 1.26;

  /// -- Calculation of the cross section proper


  // -- order 0 terms
  double E0 = Enu - gDELTA;
  if(E0<electron_mass_c2) E0=electron_mass_c2;
  double p0 = sqrt(E0*E0-electron_mass_c2*electron_mass_c2);
  double v0 = p0/E0;

  // -- order 1 terms
  const double Ysquared = (gDELTA*gDELTA-electron_mass_c2*electron_mass_c2)/2;
  double E1 = E0*(1.0-Enu/proton_mass_c2*(1-v0*cosThetaLab))-Ysquared/proton_mass_c2;
  if(E1<electron_mass_c2) E1=electron_mass_c2;
  double p1 = sqrt(E1*E1-electron_mass_c2*electron_mass_c2);
  double v1 = p1/E1;

  double Gamma =
    2*(f+f2)*g*((2*E0+gDELTA)*(1-v0*cosThetaLab)-
    electron_mass_c2*electron_mass_c2/E0)
    +(f*f+g*g)*(gDELTA*(1+v0*cosThetaLab)+electron_mass_c2*electron_mass_c2/E0)
    +(f*f+3*g*g)*((E0+gDELTA)*(1-cosThetaLab/v0)-gDELTA)
    +(f*f-g*g)*((E0+gDELTA)*(1-cosThetaLab/v0)-gDELTA)*v0*cosThetaLab;

  sigma =
    ((f*f+3*g*g)+(f*f-g*g)*v1*cosThetaLab)*E1*p1
    -Gamma/proton_mass_c2*E0*p0;
  sigma *= Sigma0/2;
  sigma *= hbarc * hbarc; // Convert from MeV^-2 to mm^2 (native units for GEANT4)

  return sigma;
}

double IBDCrossSec::dSigmadThetaVogelTable(const double Enu, const double cosThetaLab) {

  return _xs_table->Interpolate(Enu,cosThetaLab);
}

void IBDCrossSec::LoadTable(Strategy strat) {
  info << "IBDCrossSec::LoadTable : Loading pre-computed cross section table. Strategy : " << _str_map.at(strat) << newline;
  // Check if the strategy is of table type
  if (strat != VogelTable) {
    warn << "IBDCrossSec::LoadTable : Strategy " << _str_map.at(strat) << " does not require loading a table. Aborting." << newline;
    return;
  }

  // See if the $RATDecayDataDir environment variable exists.
  std::stringstream dirName;
  char* ibdDataDir = getenv("IBDDATADIR");
  if( ibdDataDir == NULL )
    {
      //No specific directory set, so use the rat default
      dirName << getenv("RATROOT") << "/data";
    }
  else
    {
      dirName << ibdDataDir;
    }
  dirName << "/cross_section_tables.root";

  _input_file = TFile::Open(dirName.str().c_str());
  if (!_input_file) {
    Log::Die("IBDCrossSec::LoadTable : Input file not found " + dirName.str(),1);
  }

  switch (strat) {
  case VogelTable:
    _xs_table = (TH2F*)_input_file->Get("ibd_dSigmadCosTheta_strat0");
    break;
  default:
    Log::Die("IBDCrossSec::LoadTable : Unknown strategy.");
  }
}

double IBDCrossSec::Sigma(const double Enu) {
  // This method calculates an integral over all angles
  // First do the most basic check
  if (Enu == 0) return 0.0;
  if (Enu < 0) {
    std::stringstream ss;
    ss << "[ESCrossSec]::Sigma : Invalid neutrino Energy ( Enu = "
        << Enu << " ).";
    RAT::Log::Die(ss.str(),1);
  }

  double   sigma = 0.0;
  double cos_theta;
  const int num_steps = 200;
  const double dstep = 2./(double)num_steps;// changing between 200 and 2000 bins changes result by 0.05%
  for (int i = 0; i < num_steps; ++i ) {
      cos_theta = -1.0 + i*dstep;
      sigma += dSigmadTheta(Enu,cos_theta)*dstep;
  }
  return sigma;

}
} /* namespace RAT */