#ifndef __JSIRENE__JPYTHIA__
#define __JSIRENE__JPYTHIA__

#include <cmath>

#include "JAAnet/JAAnetToolkit.hh"


/**
 * \author mdejong
 */

namespace JSIRENE {}
namespace JPP { using namespace JSIRENE; }

namespace JSIRENE {

  /**
   * Auxiliary class to determine EM-equivalent energy as a function of PDG particle code and energy.
   */
  struct JPythia {
    /**
     * Default constructor.
     */
    JPythia() 
    {}


    /**
     * Get EM-equivalent energy.
     *
     * \param   type    particle type   [PDG]
     * \param   Es      particle energy [GeV]
     * \return          EM-equivalent energy [GeV]
     */
    double operator()(const int type, const double Es) const
    {
      using namespace JAANET;

      switch (type) {
	
      case TRACK_TYPE_PHOTON:
      case TRACK_TYPE_NEUTRAL_PION:
      case TRACK_TYPE_NEUTRAL_ANTIPION:
      case TRACK_TYPE_ELECTRON:
      case TRACK_TYPE_ANTIELECTRON:
	return Es;
	
      case TRACK_TYPE_CHARGED_PION_PLUS:
      case TRACK_TYPE_CHARGED_PION_MINUS:
	return getEnergy(Es);

      default:
	return 0.0;
      }
    }


    /**
     * Get equivalent EM-energy for given pion energy.
     *
     * Reference:
     * Mona Dentler, "Investigation of the One-Particle Approximation in the ANTARES simulation package KM3",\n
     * Bachelorarbeit Erlangen Centre for Astroparticle Physics,\n
     * Friedrich-Alexander-Universit&auml;t, Erlangen-N&uuml;rnberg.
     *
     * \param   Es      particle energy [GeV]
     * \return          EM-equivalent energy [GeV]
     */
    static inline double getEnergy(const double Es) 
    {
      static const double a     =  72.425;
      static const double b     = -49.417;
      static const double c     =   5.858;
      static const double d     = 207.252;
      static const double e     = 132.784;
      static const double f     = -10.277;
      static const double g     = -19.441;
      static const double h     =  58.598;
      static const double i     =  53.161;
      static const double kref  =   2.698;
      
      static const double u     = (a - f) / kref;

      const double x = log10(Es);
      const double y = (e + x*(d + x*(c + x*(b + x*(a + x*u)))))  /  (i + x*(h + x*(g + x*(f + x*u))));

      return pow(10.0, y - kref);
    }
  };


  /**
   * Function object for relative light yield as a function of GEANT particle code.
   */
  static const JPythia pythia;
}

#endif