//  ProtonESgen.cc
// Contact person: Christopher Jones <christopher.jones@physics.ox.ac.uk>
// See  ProtonESgen.hh for more details
//———————————————————————//
///
/// Generates an neutrino-proton elastic scattering event, based on the
/// cross-section as function of neutrino energy and the electron's
/// recoil energy. This is the main vertex generator used for the supernova generator and is set to default as a supernova generator
/// This generator is very similar to the electron neutrino elastic scattering vertex generator.
/// Any neutrino spectrum can be input into this vertex as a RATDB file and these spectra are found in the pes_spectrum.RATDB files
/// C. Jones (Oxford) -30/04/2013
///
////////////////////////////////////////////////////////////////////

// - RAT includes
#include <RAT/ProtonESgen.hh>
#include <RAT/DB.hh>
#include <RAT/Log.hh>

#include <RAT/ProtonESgenMessenger.hh>

// - CLHEP includes
#include <CLHEP/Random/Randomize.h>

// - ROOT includes
#include <TGraph.h>
#include <CLHEP/Units/SystemOfUnits.h>
#include <CLHEP/Units/PhysicalConstants.h>
#include <TMath.h>

using namespace CLHEP;

namespace RAT {

const double  RAT::ProtonESgen::fGf    = 1.166371e-11;        	// Fermi constant (MeV^-2)


// This class is a helper to take care of the type of spectrum that is going to be used and
// reads the RATDB entries accordingly.

// In this case arg_typedbname = "sn_spectrum for the proton-neutrino elastic scattering events from a supernova

ProtonESgen::ProtonESgen() :fFluxMax(0.),fGenLoaded(false),fDBName("pes_spectrum")
{
  fMessenger = new ProtonESgenMessenger(this);    	// create a messenger to allow user to change some parameters
  //Default values for Supernova spectrum
  fGenType = "SN_all";
  fNumFP = 59;
  LoadGenerator();

  // Initialize pointers
#ifdef G4DEBUG
  debug << "[ProtonESgen]::In constructor with gentype [" << fGenType << "]." << newline;
#endif

  fMassProton = proton_mass_c2;

#ifdef G4DEBUG
  debug << "[ProtonESgen]::ProtonESgen :Leaving constructor." << newline;
#endif

}

void ProtonESgen::LoadGenerator()
{
  if (fGenLoaded) return;

  /* All possible spectra should be saved in the same format as the spectra in the pes_spectrum.ratdb file.
     The user defined input neutrino spectrum is given by the user in this file.
     The user defined spectra should also be normalised */

  DBLinkPtr linkdb;
  linkdb = DB::Get()->GetLink(fDBName, fGenType);
  fFlux.Set(linkdb->GetDArray("spec_energy"), linkdb->GetDArray("spec_flux"));
  fEnuTbl = linkdb->GetDArray("spec_energy");
  fFluxTbl = linkdb->GetDArray("spec_flux");
  fEnuMin = linkdb->GetD("emin");
  fEnuMax = linkdb->GetD("emax");
  fProbmax = linkdb->GetD("probmax");
  fGenLoaded = true;
}


ProtonESgen::~ProtonESgen()
{
  //Default destructor
  if(fMessenger != NULL)
    {
      delete fMessenger;
      fMessenger = NULL;
    }
}


void ProtonESgen::GenerateEvent(const G4ThreeVector &nu_dir,G4LorentzVector &neutrino,G4LorentzVector &proton) const
{

  // Check if the generator has been loaded successfully
  if (!fGenLoaded)
    {
      G4Exception("[ProtonESgen]::GenerateEvent","ArgError",FatalErrorInArgument,"Vertex generation called but it seems that it is not ready yet.");
    }

  double eNu;
  GenInteraction(eNu);

  // Compute nu 4-momentum
  neutrino.setVect(nu_dir*eNu); // MeV (divide by c if need real units)
  neutrino.setE(eNu);

  double Te = SampleRecoilEnergy((double)eNu);

  G4double theta_e = acos(sqrt((Te*(fMassProton+eNu)*(fMassProton+eNu))/(2*fMassProton*eNu*eNu + eNu*eNu*Te)));

  double totEproton = Te + fMassProton;
  double p_proton = sqrt(totEproton*totEproton - fMassProton*fMassProton);

  G4ThreeVector p_mom(p_proton*nu_dir);
  G4double phi_e = twopi * HepUniformRand();
  G4ThreeVector rotation_axis = nu_dir.orthogonal();
  rotation_axis.rotate(phi_e,nu_dir);

  p_mom.rotate(theta_e,rotation_axis);

  proton.setVect(p_mom);
  proton.setE(totEproton);

}

//Reset the generator
void ProtonESgen::Reset()
{
  LoadGenerator();
}

//current method
void ProtonESgen::GenInteraction(double &energy) const
{
  bool passed=false;
  while(!passed)
    {
      // generate a random energy
      double prob = 0.0;
      energy = fEnuMin + (fEnuMax-fEnuMin)*HepUniformRand();
      prob = fProbmax*HepUniformRand();
      passed = fFlux(energy) >= prob;
    }
}

// This controls the spectrum which is chosen from the ratdb file
void ProtonESgen::SetGenType(const G4String &newgentype)
{
  if (fGenType != newgentype )
    {
#ifdef G4DEBUG
      detail << "[ProtonESgen]::SetGenType : Switching to chosen neutrino spectrum [ " << newgentype << " ]." << newline;
#endif
      fGenType = newgentype;
      fGenLoaded = false;
      LoadGenerator();
    }
}

void ProtonESgen::SetNumFreeProtons(double newNumFP)
{
  if (fNumFP != newNumFP)
    {
      fNumFP = newNumFP;
      info << "[ProtonESgen::SetNumFreeProtons]- Number of free protons set to :" << fNumFP << "1e30\n";
      fGenLoaded = false;
      LoadGenerator();
    }
}

/* This function calculates the Proton-Neutrino Cross section,
	for a given neutrino energy and proton recoil energy */
double ProtonESgen::CrossSection(double nuEnergy,double tEnergy,double numFP) const
{
  double numFreeProtonsSnoplus = numFP*1e30;        	// number of free protons in the SNOPLUS detector (default to 59)
  double crossSection=  0.0;
  double constantCrossA = 4.83;
  double C_a = 1.27/2.0;
  double C_v = 0.04;

  if (nuEnergy != 0.0)
    {
      crossSection = (fGf*fGf*fMassProton/pi)*( ((1.0-fMassProton*tEnergy)/2.0*nuEnergy*nuEnergy)*C_v*C_v +((1.0+fMassProton*tEnergy)/2.0*nuEnergy*nuEnergy)*C_a*C_a );
    }
  else
    {
      crossSection = constantCrossA*(1.0e-42);
    }

  crossSection = crossSection*numFreeProtonsSnoplus;
  return crossSection;
}

/* This function chooses a Proton kinetic energy,
	using the proton-neutrino ES cross section at a given neutrino energy */
double ProtonESgen::SampleRecoilEnergy(double Enu) const{

  G4double Te = 0.0;

  // Get the shape of the differential cross section.
  TGraph *dsigmadt = DrawdSigmadT(Enu);

  // Interpolate between the discrete points
  CLHEP::RandGeneral *rndm = new CLHEP::RandGeneral(dsigmadt->GetY(),dsigmadt->GetN(),0);

  Te = rndm->shoot()*(dsigmadt->GetX())[dsigmadt->GetN()-1];
#ifdef G4DEBUG
  debug << "[ProtonESgen]::SampleRecoilEnergy : Random \'shoot\' for Enu = " << Enu << " returned a recoil of " << Te << newline;
#endif
  // Clean up memory
  delete rndm;
  delete dsigmadt;

  return Te;
}

/* This function creates a temporary TGraph to create the Proton-nu ES cross section.*/
TGraph *ProtonESgen::DrawdSigmadT(double Enu) const {
#ifdef RATDEBUG
  debug << "[ProtonESgen]::DrawdSigmadT : Sampling Enu=" << Enu << " MeV." << newline;
#endif

  TGraph *g = new TGraph();
  const int npoints = 1000;
  const double emin = 0.0, emax = 2.0*Enu*Enu/fMassProton;
  const double xwid = emax - emin;
  const double xstep = xwid/(double)npoints;
#ifdef RATDEBUG
  debug << "[ProtonESgen]::DrawdSigmadT : Sampling data: emin=" << emin << " MeV."
	<< " emax=" << emax << " MeV xwid=" << xwid << " MeV  xstep=" << xstep
	<< newline;
#endif

  double Te,dsigma;
  for (int ip = 0; ip < npoints; ip++)
    {
    Te = emin + (double)ip*xstep;
    dsigma = CrossSection(Enu,Te,fNumFP);
    g->SetPoint(ip,Te,dsigma);
    }
  return g;
}


} // namespace RAT