////////////////////////////////////////////////////////////////////
/// \class RAT::ProtonESgen
///
/// \author Christopher Jones <christopher.jones@physics.ox.ac.uk> -- contact person
/// \date 09-04-2013
///
/// REVISION HISTORY:\n
///	09-04-2013 : Chris Jones -- first version
///
/// \brief Implements the generation of a proton neutrino-electron elastic scattering event.
/// This has been built from the electron neutrino elastic scattering generator
///
////////////////////////////////////////////////////////////////////////


#ifndef __RAT_ProtonESgen__
#define __RAT_ProtonESgen__

/// - RAT includes
#include <RAT/LinearInterp.hh>
#include <RAT/GLG4PosGen.hh>


/// - Geant4 includes
#include <G4ThreeVector.hh>
#include <G4LorentzVector.hh>

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

// ROOT headers
#include <TGraph.h>
#include <TMath.h>


namespace RAT {

/// Forward declarations in the namespace
class ProtonESgenMessenger;

class ProtonESgen{
public:

  /**
   *  Default constructor.
   */
  ProtonESgen();

  /** Default destructor. */
  ~ProtonESgen();


  /**
   * Generate random event vectors.
   *
   * \param[in] nu_dir Incoming neutrino direction (lab coordinates).
   * \param[out] neutrino Outgoing neutrino direction (lab coordinates).
   * \param[out] proton Outgoing proton direction (lab coordinates).
   */
  void GenerateEvent(const G4ThreeVector &nu_dir,G4LorentzVector &neutrino, G4LorentzVector &proton) const;
  void SetGenType(const G4String &newgentype);
  void SetNumFreeProtons(double newNumFP);

  /** Getter for the spectrum being used */
  inline const G4String & GetGenType() const { return fGenType;};

  /**
   * Getter of the DB entry to input the spectrum from.
   * @return name of the DB name.
   */
  const G4String GetDBName() const  {return fDBName;};

  void SetDBName(const G4String name) {fDBName = name;};

  double Flux(float energy) const { return fFlux(energy); };
  double CrossSection(double nuEnergy,double tEnergy,double numFP)const;
  double SampleRecoilEnergy(double Enu) const;
  TGraph *DrawdSigmadT(const double Enu) const;

private:

  static const double  fGf; 		/// Fermi constant (GeV^-2)

  /** Private member for load the database and cross-section data. */
  void LoadGenerator();

  /** Generate the interaction given the neutrino energy .*/
  void GenInteraction(double &energy) const;

  /** Resets internal vectors. To be removed. */
  void Reset();

  /** Generator type */
  G4String fGenType;

protected:
  /** ProtonESgen Messenger **/
  ProtonESgenMessenger *fMessenger;


  /**SPECTRUM PARAMETERS */

  /** Upper limit for the neutrino energy.*/
  double fEnuMax;
  /** Lower limit for the neutrino energy.*/
  double fEnuMin;

  double fFluxMax;

  /** The number of free protons for proton elastic scattering in units 10^30 **/
  double fNumFP;

  /** Upper limit on the probability of the normalised neutrino flux-cross section convolution.*/
  double fProbmax;

  /** Linear interpolation of the neutrino flux */
  LinearInterp<double> fFlux;

  /** vector of neutrino energy points in the neutrino spectrum shape. */
  std::vector<double> fEnuTbl;

  /** Normalized flux in the neutrino spectrum shape. */
  std::vector<double> fFluxTbl;

  /** Generator loaded flag. */
  G4bool  fGenLoaded;

  /** proton mass */
  G4double fMassProton;

  /** Name of the database entry to read the input spectrum from.
   * Defaults to sn_spectrum (ideally)
   */
  G4String fDBName;


};

} // namespace RAT

#endif