////////////////////////////////////////////////////////////////////
/// \class RAT::ESCrossSec
///
/// \brief Calculates neutrino-electron elastic scattering.
/// (based on original QSNO code by F. Duncan, M. Chen and Y. Takeuchi).
///
/// \author Nuno Barros <nfbarros@hep.upenn.edu> -- contact person
///
/// REVISION HISTORY:\n
/// 23-NOV-2010 - Nuno Barros 	- Imported and adapted into SNO+ RAT
///             				- Shamelessly taken from SNO QPhysics class PNuE
///
/// 22-JUN-2012 - Nuno Barros	- Cleaned up code and solved a small problem with unit conversion.
///
/// 07-JUL-2012 - Nuno Barros 	- Removed all Geant4 streamers replacing them by RAT log objects.
///								- Changed geant4 data types by C++ data types (suggestion by CIC).
///								- Revised constness of members.
///
///
/// \details Calculates the neutrino-electron ES cross section.
///			The calculation is performed either for the ES_e (W+Z) channel, or the ES_mu,tau (Z) channel.
/// 		Some remarks concerning the calculations.
///			\param E,Enu 	: Neutrino energy (MeV).
///			\param T,Te		: Recoil electron energy (MeV).
///			\return Cross section in units of $10^{-42}cm^{2}$
///
///			Available strategies for the ES cross section calculation (set by messenger):
///			- 1 : Original routine from QSNO::PNuE (Bahcall).
///			- 2 : Improved routine from QSNO::PNuE (without rad. corrections).
///			- 3 : Improved routine from QSNO::PNuE (with rad. corrections - analytical).
///			- 4 (default) : Improved routine from QSNO::PNuE (with rad. corrections - table).
///
///
////////////////////////////////////////////////////////////////////

#ifndef __RAT_ESCrossSec__
#define __RAT_ESCrossSec__

// G4 headers
#include <globals.hh>
// RAT headers
#include <RAT/LinearInterp.hh>
#include <RAT/GLG4StringUtil.hh>
// ROOT headers
#include <TGraph.h>

namespace RAT {

  // Forward declarations within the namespace
  class ESCrossSecMessenger;


    class ESCrossSec {

    public:
	enum NuEType {nue,nuebar,numu,numubar};

	ESCrossSec(const char* flavor = "nue");

	~ESCrossSec();

	// Set's the defaults for the calculation
	void Defaults();


	/**
	 * @brief Calculate the total cross section for the neutrino energy Enu.
	 *
	 * @param Enu
	 * @return total cross section in units of \f$ 10^{-42} cm^{2} \f$ .
	 */
	double Sigma(const double Enu) const;


	/**
	 * @brief Calculate the differential cross section for the neutrino energy Enu.
	 *
	 * @param Enu Incoming neutrino energy (MeV).
	 * @param Te  Recoil electron energy (MeV).
	 * @return Differencial cross section \f$ \frac{d\sigma}{dT} \f$ in units of \f$ 10^{-42} cm^{2} \f$ .
	 */
	double dSigmadT(const double Enu,const double Te) const;

	/**
	 * Integrate the differential cross section between recoil energies T1 and T2.
	 *
	 * @param Enu Incoming neutrino energy (MeV).
	 * @param T1  Lower limit of recoil energy interval (MeV).
	 * @param T2  Upper limit of recoil energy interval (MeV).
	 * @return \f$ \left.\frac{d\sigma}{dT}\right|_{\left[T1;T2\right]} \f$ in units of \f$ 10^{-42} cm^{2} \f$ .
	 */
	double IntegraldSigmadT(const double Enu,const double T1,const double T2) const ;

	/**
	 * @brief Calculate the differential cross section as a function of the recoil angle.
	 *
	 * @param Enu Incoming neutrino energy.
	 * @param CosTh Cosine of laboratory recoil angle of the electron.
	 * @return  \f$ \frac{d\sigma}{d\cos \theta} \f$ in units of \f$ 10^{-42} cm^{2} \f$ .
	 */
	double dSigmadCosTh(const double Enu,const double CosTh) const;

	/**
	 * @brief 3D equivalent of ESCrossSec::dSigmadCosTh
	 *
	 * @param Enu Incoming neutrino energy.
	 * @param theta Cosine (FIXME) of laboratory recoil angle of the electron.
	 * @param phi 	Azimuthal recoil angle.
	 * @return  \f$ \frac{d\sigma}{d\Omega} \f$ in units of \f$ 10^{-42} cm^{2} \f$ .
	 *
	 * @note Usually not used. Needs further verification.
	 *
	 */
	double dSigmadOmega(const double Enu,const double theta, const double phi) const;


	//-----------------------------------------------------------------------
	/**
	 * @brief Setter for which calculation to use.
	 *
	 * The available calculations are:
	 * 	-# 1 : No radiative corrections (Bahcall first calculation).
	 * 	-# 2 : No radiative corrections from table.
	 * 	-# 3 : Full analytical calculation with radiative corrections.
	 * 	-# 4 : Full calculations with radiative corrections from table (default).
	 *
	 * @param ii Option for calculation to use.
	 *
	 */
	void SetRadiativeCorrection(const int ii);

	/**
	 * @brief Getter of the cross section calculation type being used.
	 *
	 * @return index of calculation type.
	 */
	inline int GetRadiativeCorrection( ) const {return fRadiativeCorrection;};

	/**
	 * @brief Setter for the Weak mixing angle.
	 *
	 * @param sintw Weak mixing angle : \f$ \sin \theta_{W} \f$
	 *
	 */
	void SetSinThetaW(const double &sintw ) ; //{fsinthetaW2 = sintw;};

	/**
	 * @brief Getter for the Weak mixing angle.
	 *
	 * @return Weak mixing angle : \f$ \sin \theta_{W} \f$
	 *
	 */
	double GetSinThetaW() const {return fsinthetaW2;};

	/**
	 * @brief Setter for the neutrino type to be used.
	 *
	 * @param reaction neutrino type (nue,numu,nuebar,numubar)
	 */
	void SetReaction(const std::string &reaction);

	/**
	 * @brief Getter for the neutrino type to be used.
	 *
	 * @return reaction neutrino type (nue,numu,nuebar,numubar)
	 */
    const std::string& GetReaction() const {return fReactionStr;};

	/**
	 * @brief Return a TGraph with the differential cross section for an incoming neutrino with energy Enu.
	 * @param Enu Incoming neutrino energy (MeV)
	 * @return TGraph with the shape of \f$ \frac{d\sigma}{dT} \f$ in units of \f$ 10^{-42} cm^{2} \f$ .
	 */
	TGraph *DrawdSigmadT(const double Enu) const;

	/**
	 * Returns the global normalization of the cross section calculation.
	 * For precision reasons, the cross-section is performed on a different scale, and therefore
	 * any result returned by the calculation is missing this scale, which has to be applied separately.
	 *
	 * @return cross section scaling factor (1e-42)
	 */
	double CrossSecNorm() const {return 1e-42;};

    protected:

	/**
	 * @brief Internal function to load the data from the DB.
	 */
	void LoadTablesDB();

	void CalcG();


    private:

	/**
	 * Rename of ESCrossSec::Sigma. To be discontinued.
	 *
	 * @see ESCrossSec::Sigma
	 * @param Enu Incoming neutrino energy.
	 * @return total cross section in units of \f$ 10^{-42} cm^{2} \f$ .
	 */
    double SigmaLab(double Enu) const;


	NuEType fReaction;   	/// Reaction type
	std::string fReactionStr;	/// String characterizing the reaction type

	double fEmin,fEmax;	/// Auxiliary variables to deal with the energies

	// Some constants
	static const double  fGf; 		/// Fermi constant (GeV^-2)
	static const double  fhbarc; 	/// hbar*c (MeV*fm)
	static const double  fhbarc2;	/// hbar*c^2(GeV^2 mb)
	static const double  falpha;	/// radiative correction term


	/**
	 * This variable is defined as static (and not const) because it can be changed
	 * in the macro file. However, the change should propagate to all instances of the class
	 */
	static double  fsinthetaW2;

	double  fMe;			/// electron mass
	double  sigmaoverme; 	/// \f$ \frac{\sigma}{m_{e^{-}}}\f$

	double  fgL;
	double  fgR;

	//-----------------------------------------------------------------------
	// add by y.t. 16-JAN-2003
	int     fRadiativeCorrection;  // flag to use radiative correction or not
	//double  fL(double x);        // internal routine


	/// Vars to manipulate the tables. Since the tables are the same for all instances, these can be static
	static const double  fTableTeMin;
	static const double  fTableTeMax;

	// Total cross section table
	int    fNDataTot;
	float  fEnuStepTot;
	LinearInterp<double> fTableTot_e;

	// Differential cross section table
	int                  fNDataDif;
	float                fEnuStepDif;
	LinearInterp<double> fTableDif_e;
	std::vector<double>  fTableDif;

	ESCrossSecMessenger *fMessenger;

    };

//
// Some inline definitions to make the code a bit faster
//

    inline void ESCrossSec::CalcG()
    {
	// calculate the gL and gR for
	// the reaction type

	if ( fReaction == nue)           {

	    fgL =  0.5 + fsinthetaW2;
	    fgR =  fsinthetaW2;

	} else if (fReaction == nuebar)  {

	    fgL =   fsinthetaW2;
	    fgR =   0.5 + fsinthetaW2;

	} else if (fReaction == numu)    {

	    fgL =   -0.5 + fsinthetaW2;
	    fgR =    fsinthetaW2;

	} else if (fReaction == numubar) {

	    fgL =   fsinthetaW2;
	    fgR =   -0.5 + fsinthetaW2;

	}
	else throw std::invalid_argument("Unknown reaction " + fReactionStr);
    }
}



#endif