#ifndef __pythia__
#define __pythia__

#include "JAAnet/JPDB.hh"

/**
 * \file
 * This file contains converted Fortran code from km3.
 * \author mdejong
 */

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

namespace JSIRENE {

  using JAANET::JPDB;
   
  // *** high-energy fit to weights for a high-energy pion (same as v4r5 version)
  inline double opa_weight_high_e(const double ekin)
  {
    static const double Ma    =  72.425;
    static const double Mb    = -49.417;
    static const double Mc    =   5.858;
    static const double Md    = 207.252;
    static const double Me    = 132.784;
    static const double Mf    = -10.277;
    static const double Mg    = -19.441;
    static const double Mh    =  58.598;
    static const double Mi    =  53.161;
    static const double Mkref =   2.698;
    //	 leading_coef Mlc=(Ma-Mf)/Mkref
    static const double Mlc   = (Ma-Mf)/Mkref;

    double num, denom, lognp, Ee, logE;
      
    //     	 implement the energy-fractions according to M. Dentler: ANTARES-SOFT-2012-009. E must be the kinetic energy in GeV. Here, 'everything else' is treated as a pion. Other than charged pions, mostly kaons of all varieties and protons fall into this category. These fits were made mostly at high energies - at low energies, other fits are used.
      	 
    if (ekin > 0.2) 
      logE=log10(ekin);
    else 
      // nphotons constant below 0.2 GeV (lowest fit point)
      return 0.292;

    num   = Me + Md*logE + Mc*pow(logE,2) + Mb*pow(logE,3) + Ma *pow(logE,4) + Mlc*pow(logE,5);
    denom = Mi + Mh*logE + Mg*pow(logE,2) + Mf*pow(logE,3) + Mlc*pow(logE,4);
    
    if (denom <= 0) {
      return 0.;
    }

    lognp = num/denom;
    Ee    = pow(10.0, lognp-Mkref);
    
    return Ee/ekin;
  }
    
  inline double ngamma_elec(const double ekin)
  {
    static const double eleca =  1.33356e5;
    static const double elecb =  1.66113e2;
    static const double elecc = 16.4949;
    static const double elecd =  1.5385e5;
    static const double elece =  6.04871e5;
    
    return (eleca*ekin+elecb) * exp(-ekin/elecc) + (elecd*ekin+elece)* (1.-exp(-ekin/elecc));
  }
	    
  inline double weight_pion(const double ekin)
  {
    static const double pia =  0.538346;
    static const double pib =  1.32041;
    static const double pic =  0.737415;
    static const double pid = -0.813861;
    static const double pie = -2.22444;
    static const double pif = -2.15795;
    static const double pig = -6.47242;
    static const double pih = -2.7567;
    static const double pix =  8.83498;
      
    if      (ekin < 0.06) 
      return 1e4*pia/ngamma_elec(ekin);
    else if (ekin < 0.15) 
      return pix*1e4*ekin/ngamma_elec(ekin);
    else {
      const double lek=log(ekin);
      return (pib*1e5*ekin + (pow(ekin,(1.-1./pic))) * (pid*1e4 + 1e4*pie*lek + 1e4*pif*pow(lek,2) + 1e3*pig*pow(lek,3) + 1e2*pih*pow(lek,4)))/ngamma_elec(ekin);
    }
  }

  // **** calculates the expected number of photons for pions as a function of energy
  inline double weight_kaon(const double ekin)
  {
    static const double ka  = 12.7537;
    static const double kb  = -1.052;
    static const double kc  =  3.49559;
    static const double kd  = 16.7914;
    static const double ke  = -0.355066;
    static const double kf  =  2.77116;

    if (ekin > 0.26) 
      return (1e4*ka*(ekin+kb)*(1.-exp(-ekin/kc)) + 1e4*(kd*ekin+ke)*exp(-ekin/kc))/ngamma_elec(ekin);
    else
      return kf*1e4/ngamma_elec(ekin);
  }      

  // **** calculates the expected number of photons for K0-short as a function of energy
  inline double weight_kshort(const double ekin)
  {
    static const double k0sa  = 0.351242;
    static const double k0sb  = 0.613076;
    static const double k0sc  = 6.24682;
    static const double k0sd  = 2.8858;

    return (k0sa*1e5 + k0sb*1e5*ekin + ekin*k0sc*1e4*log(k0sd + 1./ekin))/ngamma_elec(ekin);
  }

  // **** calculates the expected number of photons for K0-short as a function of energy
  inline double weight_klong(const double ekin)
  {
    static const double k0la  = 2.18152;
    static const double k0lb  = -0.632798;
    static const double k0lc  = 0.999514;
    static const double k0ld  = 1.36052;
    static const double k0le  = 4.22484;
    static const double k0lf  = -0.307859;
      
    if (ekin < 1.5)
      return (1e4*k0la + ekin*1e5*k0lb*log(ekin)  + 1e5*k0lc*pow(ekin,1.5))/ngamma_elec(ekin);
    else
      return (ekin*k0ld*1e5 + pow(ekin,1.-1./k0le)*k0lf*1e5)/ngamma_elec(ekin);
  }

  // **** calculates the expected number of photons for protons as a function of energy
  inline double weight_proton(const double ekin)
  {     
    static const double pa   =  12.1281;
    static const double pb   = -17.1528;
    static const double pc   =   0.573158;
    static const double pd   =  34.1436;
    static const double pe   =  -0.28944;
    static const double pf   = -13.2619;
    static const double pg   =  24.1357;
      
    return (1e4*(pa*ekin+pb)*(1.-exp(-ekin/pc)) + 1e4*(pd*ekin +pe+ pf*pow(ekin,2) + pg*pow(ekin,3))*exp(-ekin/pc)) / ngamma_elec(ekin);
  }

  // **** calculates the expected number of photons for neutrons as a function of energy
  inline double weight_neutron(const double ekin)
  {
    static const double na   = 1.24605;
    static const double nb   = 0.63819;
    static const double nc   = -0.802822;
    static const double nd   = -0.935327;
    static const double ne   = -6.1126;
    static const double nf   = -1.96894;
    static const double ng   = 0.716954;
    static const double nh   = 2.68246;
    static const double ni   = -9.39464;
    static const double nj   = 15.0737;
      
    if (ekin > 0.5) {
      const double lek=log(ekin);
      return (na*1e5*ekin + 1e3*pow(ekin,1.-1./nb) * (100.*nc+100.*nd*lek + 10.*ne*pow(lek,2) + 11.*nf*pow(lek,3)))/ngamma_elec(ekin);
    } else
      return (1e3*ng + 1e4*nh*ekin + 1e4*ni*pow(ekin,2) + 1e4*nj*pow(ekin,3))/ngamma_elec(ekin);
  }


  // Routine to calculate the effective electron-shower energy using the
  // one-particle approximation.
  // includes low-energy fits per particle, and high-energy fits for pions
  inline double opa_efrac(const int ipart, const double ekin)
  {
    using namespace JPP;

    double etemp, weight_part;

    // neutrino or muon
    if (ipart == GEANT4_TYPE_NEUTRINO       || 
	ipart == GEANT4_TYPE_ANTIMUON       || 
	ipart == GEANT4_TYPE_MUON           ||
	ipart == GEANT4_TYPE_ANTITAU        || 
	ipart == GEANT4_TYPE_TAU) {
      return 0.0;
    }
      
    // 100% energy to EM showers if electron, positron, pi0, or gamma
    if (ipart == GEANT4_TYPE_PHOTON         || 
	ipart == GEANT4_TYPE_ANTIELECTRON   || 
	ipart == GEANT4_TYPE_ELECTRON       ||
	ipart == GEANT4_TYPE_NEUTRAL_PION   ||
	ipart == GEANT4_TYPE_ETA) {
      return 1.0;
    }

    // low-energy fits went up to 40 GeV      
    if (ekin <= 40) 
      etemp=ekin;
    else
      etemp=40.;
      
    // otherwise, calculate fractions the long way. First we get the
    // expected number of photons from the particle:
    if      (ipart == GEANT4_TYPE_PION_PLUS ||
	     ipart == GEANT4_TYPE_PION_MINUS)
      weight_part=weight_pion(etemp);
    else if (ipart == GEANT4_TYPE_KAON_LONG)
      weight_part=weight_klong(etemp);
    else if (ipart == GEANT4_TYPE_KAON_SHORT) 
      weight_part=weight_kshort(etemp);
    else if (ipart == GEANT4_TYPE_KAON_PLUS ||
	     ipart == GEANT4_TYPE_KAON_MINUS)
      weight_part=weight_kaon(etemp);
    else if (ipart == GEANT4_TYPE_NEUTRON) 
      weight_part=weight_neutron(etemp);
    else if (ipart == GEANT4_TYPE_PROTON    ||
	     ipart == GEANT4_TYPE_ANTIPROTON)
      weight_part=weight_proton(etemp);
    else
      //     	 when in doubt, treat it as a proton
      weight_part=weight_proton(etemp);
      	
    if (ekin <= 40)
      return weight_part;
    else if (ekin >= 1e7) 
      return opa_weight_high_e(ekin);
    else {
      // smooth transition between low-energy weights and high-energy pion weight
      const double wp40  = weight_part;
      const double whe40 = opa_weight_high_e(40.);
      const double whex  = opa_weight_high_e(ekin);	 
      return whex - (whe40-wp40)*(7.-log10(ekin))/5.398;
    }
  }

  
  /**
   * Get equivalent EM-energy for given pion energy.
   *
   * \param   type    particle type   [PDG]
   * \param   E       particle energy [GeV]
   * \return          EM-equivalent energy [GeV]
   */
  inline double pythia(const int type, const double E)
  {
    int   ipart = JPDB::getInstance().getPDG(type).geant;
    float ekin  = (float) E;

    return opa_efrac(ipart, ekin) * E;
  }
}

#endif