#ifndef FLUXHH
#define FLUXHH

/*! author: R. Coniglione, T. Heid */

#include <map>
#include <stdlib.h>
#include <functional>
#include "TH2.h"
#include "TClass.h"
#include "TROOT.h"
#include "TFunction.h"
#include "TInterpreter.h"
#include "km3net-dataformat/offline/Trk.hh"
#include "nufromgamma.hh"


/*! \brief Flux defines the interface to an (atmospheric or astrophysical) neutrino flux. */

struct Flux : public TObject 
{
    /*! \brief return the total flux (Gev-1, m-2, s-1,) */
    virtual double dNdE( int pdg_type, double loge ) =0;
    /*! \brief printing the flux info */
    virtual void Print() =0;

    virtual void prnt(ostream& out ) const {
      out << " Flux___ ";
    }

    ClassDef( Flux , 1 )
};



/*! \brief Diffuse cosmic or atmospheric flux, with a dependence on zenith angle */

struct DiffuseFlux : public Flux
{

  int nsteps; //number of integration steps

  /*! \brief return number of neutrinis per unit energy, time, and area (Gev-1, m-2, s-1, sr-1) */
  virtual double dNdEdOmega( int pdg_type, double loge, double costheta ) const =0;

  /*! \brief return the flux integrated over the full sphere (Gev-1, m-2, s-1,) */
  double dNdE( int pdg_type, double loge)
  {
    const double dct = 2.0/nsteps;
    double int_flux = 0;

    for (int i = 0; i < nsteps ; i++ )
    {
      double ct = -1 + dct * (i +0.5);
      int_flux += 2 * TMath::Pi() * dct * this->dNdEdOmega( pdg_type, loge , ct );
    } 

    return int_flux;    
  }
  
  /*! \brief return flux (Gev-1, m-2, s-1, sr-1) */
  double dNdEdOmega( const Trk& t ) const
  {
    return this->dNdEdOmega( t.type, log10( t.E), -t.dir.z ); 
  }

  DiffuseFlux(): nsteps(100) {}


  ClassDef( DiffuseFlux , 1 )
};



/*! ScalableFlux is a flux with a (changable) normalisation */


struct ScalableFlux : public Flux
{
    double norm;

    ClassDef ( ScalableFlux, 1 )
};


/*! PowerLaw flux */

struct PowerLawFlux : public ScalableFlux
{
  double gamma;

  PowerLawFlux() {
      norm = 1;
      gamma = 2.;
  }
  PowerLawFlux( double _norm, double _gamma ) {
      norm = _norm;
      gamma = _gamma;
  }

  PowerLawFlux(const PowerLawFlux& master) {
      norm = master.norm;
      gamma = master.gamma;
  }

  double dNdE( int pdg_type, double loge )
  {
    return norm * pow(10, -gamma * loge);
  }
  
  void Print(){
      printf(" gamma                = %f", gamma );
      printf(" norma                = %f", norm );
  }

  ClassDef ( PowerLawFlux, 1 )

};


/*! Flux taken from an expression that is evaluated by root TF1 */

struct NuExprFlux : public ScalableFlux
{
  TF1 nuflux;

  NuExprFlux() {
      norm = 1;
  }

  NuExprFlux( double _norm, string expr ) {
      norm = _norm;
      nuflux = TF1("nuflux",expr.c_str(), 2, 8 ); //emin emax hardcoded which is not great
  }

  NuExprFlux(const NuExprFlux& master) {
      norm = master.norm;
      nuflux = TF1("nuflux",master.nuflux.GetExpFormula().Data(),2,8) ;
  }

  double dNdE( int pdg_type, double loge )
  {
    return norm * nuflux.Eval(pow(10,loge));
  }
  
  void Print(){
      printf(" norma                = %f", norm );
			printf(" expression           = %s", nuflux.GetExpFormula().Data());
			for (int i = 0; i < nuflux.GetNpar(); i++) {
				 printf(" par #%d = %f", i, nuflux.GetParameter(i)); 
			}
  }

  std::string expression() const {
    return nuflux.GetExpFormula().Data();
  }

  virtual void prnt( ostream& out ) const {
    cout <<" NuExprFlux " << expression();
  }

  ClassDef ( NuExprFlux, 1 )

};



inline std::ostream& operator<<( std::ostream& out, const Flux& f )
{
  f.prnt( out );
  return out;
}

/*! A flux that computes the neutrino spectra from a gamma-ray flux */

struct NuFromGammaFlux : public ScalableFlux
{
  nufromgamma nuflux;
  bool use_spline; //use spline to speed up calculation

  NuFromGammaFlux() {
      norm = 1;
  }

  NuFromGammaFlux(string _gammafluxeq, double _norm=1, bool _use_spline=true )  {
    nuflux= nufromgamma(_gammafluxeq);
    norm=_norm;
    use_spline=_use_spline;
  }

  NuFromGammaFlux(const NuFromGammaFlux& master) {
      norm = master.norm;
      use_spline = master.use_spline;
      nuflux= nufromgamma(master.nuflux.f_gamma.GetExpFormula().Data());
  }

  void SetFluxFromGamma(string _gammafluxeq) { 
    nuflux= nufromgamma(_gammafluxeq);
  }


  virtual double dNdE( int pdg_type, double loge )
  {
    double val;
    /*if (use_spline) val = nuflux.F_nu_totS(loge);
      else val = nuflux.F_nu_tot(pow(10,loge));*/

    if(pdg_type>0){ //numu, nue and nutau same flux?
      if (use_spline) val = nuflux.F_nuS(loge);
      else val = nuflux.F_nu(pow(10,loge));
    }  
    else if(pdg_type<0){//anumu, anue and anutau same flux?
      if (use_spline) val = nuflux.F_antinuS(loge);
      else val = nuflux.F_antinu(pow(10,loge));
    }

    return norm*val;
  }

  void Print(){
      printf(" norm                = %f", norm );
      printf(" gamma equiation     = %s", nuflux.f_gamma.GetExpFormula().Data() );
  }


  ClassDef ( NuFromGammaFlux, 1 )

};

namespace KM3NET_FLUX {
	/*! \brief simple hashing. */
	inline unsigned int simple_hash( const char* s )
	{
		// code googled somewhere
		unsigned hash = 0;
		while (*s)
			{
				hash = hash * 101  +  *s++;
			}
		return hash;
	}

	/*! \brief simple hashing. */
	inline unsigned int simple_hash( const string s )
	{
		return simple_hash( s.c_str() );
	}
}


/*! \brief A isotropic, cosmic neutrino flux that you can set with an expression.

    example:
\code 
    IsotropicFlux f("1e-4/(E*E)"); // E-2 flux : 1e-4 (E/GeV)^-2 (GeV-1 m-2 s-1 sr-1 )
*/

struct IsotropicFlux : public DiffuseFlux
{
  typedef double (*_fluxfunc_ptr_t) ( double E );
  _fluxfunc_ptr_t _fluxfunc;
  
  IsotropicFlux( string expression )
  {
    ostringstream funcname("__getweight");
		funcname << KM3NET_FLUX::simple_hash( expression.c_str() );
		//funcname = "__getweight"+str(simple_hash( expression.c_str() ));
    
    string s =
    " double "+funcname.str()+"(double E) { \n"
    "    return " + expression + ";        \n"
    " } ";

    gInterpreter->Declare( s.c_str() );
    TFunction* pp = gROOT->GetGlobalFunction( funcname.str().c_str() );
    if (!pp) { fprintf( stderr, "internal error"); }
    
    void* p = gInterpreter->FindSym( pp->GetMangledName() );
    _fluxfunc = (_fluxfunc_ptr_t) p;
  }

  virtual double dNdEdOmega( int pdg_type, double loge, double costheta ) const
  {
    return dNdEdOmega( pow(10,loge) );
  }

  double dNdEdOmega( const Trk& t ) const
  {
    return dNdEdOmega( t.E );
  }

  double dNdEdOmega( double E ) const
  {
    return _fluxfunc( E );
  }

  void Print(){
      printf(" IsotropicFlux ");
  }


  ClassDef( IsotropicFlux, 1 );
};



// fit between 100 GeV and 1 PeV
inline double correctKneeFit1(double lE)
{
  double p0  = 5.78198    ;     
  double p1  = -9.3379    ;    
  double p2  = 7.15694    ;
  double p3  = -2.83412   ;   
  double p4  = 0.618725   ;     
  double p5  = -0.0703262 ;    
  double p6  = 0.00321813 ;     

  return  p0 + p1*pow(lE,1) + p2*pow(lE,2)+p3*pow(lE,3)+p4*pow(lE,4)+p5*pow(lE,5)+p6*pow(lE,6);
}

// fit between 100 GeV and 1 PeV
inline double correctKneeFit2(double lE)
{
  double  p0  = -0.63421   ;     
  double  p1  = 20.1904    ;    
  double  p2  = -11.415    ;
  double  p3  = 2.5948     ;   
  double  p4  = -0.294555  ;     
  double  p5  = 0.0166707  ;    
  double  p6  = -0.00037607;       
  return  p0 + p1*pow(lE,1) + p2*pow(lE,2)+p3*pow(lE,3)+p4*pow(lE,4)+p5*pow(lE,5)+p6*pow(lE,6);
}

/*! Correction function for atmospheric flux models */
inline double correctKnee(double logE)
{
  double lE = logE;
  double Phi_H3a;
  if (lE < 2.)
  {
    Phi_H3a = correctKneeFit1(2.);
  }
  else if(lE<6.)
  {
    Phi_H3a = correctKneeFit1(lE);
  }
  else if (lE<10.)
  {
    Phi_H3a = correctKneeFit2(lE);
  }
  else
  {
    Phi_H3a = correctKneeFit2(10.);
  }
  
  return Phi_H3a; 
}




/*! Honda 2006 atmospheric neutrino model */

struct Flux_Honda2006 : public DiffuseFlux
{
  bool apply_knee_correction = true;


  virtual double dNdEdOmega( int pdg_type, double loge, double costheta ) const
  {

    static constexpr double honda2006_nue[] = { -1.56213, -2.20552, 0.526815, 3.13682, 2.81648, -1.05559,
                                                -0.00122178, 0.00200366, 0.0405091, -0.21125, -3.10896,
                                                -0.00506147, -0.0527292, -0.230067, -0.395967,
                                                500, 500, 500, 500, 1.197e-4, 1.73, 0.85, 1.528e-4, 1.58, 0.205,
                                                -1.7, 8.e-6, 1.02573e-01, -6.82870e-02, 9.58633e-01, 4.07253e-02,
                                                8.17285e-01, 1., .94, .92, .92, .9, .9, 4.22e-2, 1.09, 0.115, 3.01e-3,
                                                1.08, 0.85, 2.11e-4, 1.23, 0.205, 0.016, 0.077
                                              };

    static constexpr double honda2006_anue[] = { -1.38115e+00, -2.10431e+00, 5.39290e-01, 3.06680e+00, 2.77104e+00,
                                                 -1.19389e+00, -7.36101e-04, -8.93156e-04, 3.74032e-02, -1.76907e-01,
                                                 -3.13244e+00, -3.78853e-03, -3.54783e-02, -1.36228e-01, -1.97266e-01,
                                                 1500, 1500, 1500, 500, 6.24e-5, 1.52, 0.85, 1.8336e-4, 1.58, 0.205, -1.7,
                                                 8.e-6, 1.02573e-01, -6.82870e-02, 9.58633e-01, 4.07253e-02, 8.17285e-01,
                                                 1., 1., 1., .96, .95, .96, 3.38e-2, 1.09, 0.115, 9.83e-4, 0.95, 0.85, 2.11e-4,
                                                 1.23, 0.205, 0.016, 0.077
                                               };

    static constexpr double honda2006_numu[] = { -8.79638e-01, -1.17432e+00, -9.39143e-01, 1.34865e-01, 6.88957e-01,
                                                 -1.43770e+00, 7.25078e-04, -8.69003e-03, 3.39606e-02, -1.04246e-01,
                                                 -2.76920e+00, 2.65441e-03, -1.07521e-02, -1.26859e-01, -4.20873e-01,
                                                 700, 700, 700, 700, 8.90369e-08, 2.77000e+00, 1.15000e-01, 2.649640e-08,
                                                 1.18000e+00, 8.50000e-01, -1.70000e+00, 1.000000, 1.02573e-01,
                                                 -6.82870e-02, 9.58633e-01, 4.07253e-02, 8.17285e-01, 0.97, 0.985,
                                                 0.995, 1.005, 1.015, 1.005, 1., 1., 1., 1., 1., 1., 1., 1., 1., 0., 0.
                                               };

    static constexpr double honda2006_anumu[] = { -1.14682e+00, -1.68745e+00, -1.10493e+00, 4.10708e-01, 9.74675e-01,
                                                  -1.39129e+00, 7.54875e-04, -9.12622e-03, 3.56607e-02, -1.12264e-01,
                                                  -2.81577e+00, 1.37765e-03, -1.88374e-02, -1.52999e-01, -4.52546e-01,
                                                  700, 700, 700, 700, 7.85269e-08, 2.77000e+00, 1.15000e-01, 1.36649e-08,
                                                  1.18000e+00, 8.50000e-01, -1.70000e+00, 1.000, 1.02573e-01, -6.82870e-02,
                                                  9.58633e-01, 4.07253e-02, 8.17285e-01, 0.96, 0.97, 0.98, 1.005, 1.015, 1.02,
                                                  1., 1., 1., 1., 1., 1., 1., 1., 1., 0., 0.
                                                };



    double conventional;
    //double costheta = cos(zenith_rad);
    double lE = loge;

    // I have absolutely no idea what the hell is going on here
    if (costheta >= 0.) {costheta = -costheta;}

    switch (pdg_type)
    {
    case 12:
      conventional = fluxFunctionConventional(honda2006_nue, lE, costheta);
      break;
    case -12:
      conventional = fluxFunctionConventional(honda2006_anue, lE, costheta);
      break;
    case 14:
      conventional = fluxFunctionConventional(honda2006_numu, lE, costheta);
      break;
    case -14:
      conventional = fluxFunctionConventional(honda2006_anumu, lE, costheta);
      break;
    default:
      conventional = 0.;
      break;
    }

    if (apply_knee_correction) conventional *= correctKnee( lE );
    return conventional;
  }


  // 1e4 for cm-2 to m-2
  double fluxFunctionConventional(const double *function, double lE, double costheta) const
  {
    double Elim;
    //bug in E connection Nov 08
    if ( costheta > 0.5 || costheta < -0.5) { Elim = function[15]; }
    else if ( (0.3 < costheta && costheta <= 0.5) || (-0.5 <= costheta && costheta < -0.3) ) { Elim = function[16]; }
    else if ( (0.1 < costheta && costheta <= 0.3) || (-0.3 <= costheta && costheta < -0.1) ) { Elim = function[17]; }
    else { Elim = function[18]; }

    if ( pow(10., lE) < Elim ) {
      return 1e4 * pow(10, (function[0] * pow(costheta, 5.)
                            + function[1] * pow(costheta, 4.)
                            + function[2] * pow(costheta, 3.)
                            + function[3] * pow(costheta, 2.)
                            + function[4] * costheta
                            + function[5]
                            + function[6] * pow(lE, 5.)
                            + function[7] * pow(lE, 4.)
                            + function[8] * pow(lE, 3.)
                            + function[9] * pow(lE, 2.)
                            + function[10] * lE
                            + function[11] * pow(lE, 4.) * costheta
                            + function[12] * pow(lE, 3.) * pow(costheta, 2.)
                            + function[13] * pow(lE, 2.) * pow(costheta, 3.)
                            + function[14] * lE * pow(costheta, 4.)));
    }

    else { //check costheta
      if (costheta < 0 ) { costheta = -costheta;}

      // Costheta-dependent normalization tweak to patch LE/HE regions cleanly
      // Normalization points in data file match angles below
      double normCosTheta[6] = {0.1, 0.2, 0.4, 0.6, 0.8, 1.0};
      int normIdx = 0;
      double normTweak = 1.0;
      // Tweak is not accurate very close to zero
      if (costheta > 0.0) {
        for (int i = 0; i < 6; i++)
          normIdx = (costheta > normCosTheta[i]) ? i : normIdx;
        normTweak = function[32 + normIdx] +

                    (function[32 + normIdx + 1] - function[32 + normIdx]) / (normCosTheta[normIdx + 1] - normCosTheta[normIdx]) *
                    (costheta - normCosTheta[normIdx]);
      }
      //numu and nue
      double lETeV = lE - 3;
      return 1e4 * (normTweak * function[26] * pow(pow(10., lETeV), function[25])
                    * (function[19] / (1. + function[20] * pow(10., lETeV)
                                       * sqrt((pow(costheta, 2) + pow(function[27], 2) + function[28] * pow(costheta, function[29]) + function[30] * pow(costheta, function[31])) / (1. + pow(function[27], 2) + function[28] + function[30])) / function[21])
                       + function[22] / (1. + function[23] * pow(10., lETeV)
                                         * sqrt((pow(costheta, 2) + pow(function[27], 2) + function[28] * pow(costheta, function[29]) + function[30] * pow(costheta, function[31])) / (1. + pow(function[27], 2) + function[28] + function[30])) / function[24])
                       + (function[38] / (1. + function[39] * pow(10., lETeV)
                                          * sqrt((pow(costheta, 2) + pow(function[27], 2) + function[28] * pow(costheta, function[29]) + function[30] * pow(costheta, function[31])) / (1. + pow(function[27], 2) + function[28] + function[30])) / function[40])
                          + function[41] / (1. + function[42] * pow(10., lETeV)
                                            * sqrt((pow(costheta, 2) + pow(function[27], 2) + function[28] * pow(costheta, function[29]) + function[30] * pow(costheta, function[31])) / (1. + pow(function[27], 2) + function[28] + function[30])) / function[43])
                          + function[44] / (1. + function[45] * pow(10., lETeV)
                                            * sqrt((pow(costheta, 2) + pow(function[27], 2) + function[28] * pow(costheta, function[29]) + function[30] * pow(costheta, function[31])) / (1. + pow(function[27], 2) + function[28] + function[30])) / function[46])
                         )
                       * function[48] * (1. - exp(-function[47] / pow(10., lETeV) / sqrt((pow(costheta, 2) + pow(function[27], 2) + function[28] * pow(costheta, function[29]) + function[30] * pow(costheta, function[31])) / (1. + pow(function[27], 2) + function[28] + function[30])))
                                        )) / pow(10., lE));

    }
  }

  void Print() {
    printf("Honda2006");
  }

  ClassDef( Flux_Honda2006, 1)
};



/*! Prompt neutrino flux according to Sarsevic */
struct Flux_prompt_Sarsevic_std : public DiffuseFlux
{
  bool apply_knee_correction = true;

  virtual double dNdEdOmega( int pdg_type, double loge, double costheta ) const
  {

    static constexpr double e_joint_nu_prompt = 100.;
// 100
// with all of these, we have "e joint" as first param, and the second column is *I think* the error
    static constexpr double sarcevic_std_anue[] = {
      0.0938106731958, /*0.00968084556496*/
      0.127447604594, /*0.0104518800441*/
      0.0690640747739, /*0.00928723182021*/
      -0.00202883752944, /*0.00733072332339*/
      -0.00740634544259, /*0.00464412823313*/
      -6.44177730199, /*0.00191354750982*/
      -2.1934999156e-05, /*2.35828231602e-07*/
      -0.00060750673143, /*2.12405966392e-06*/
      0.00516577885096, /*1.54299379371e-05*/
      -0.0442845804904, /*9.88758722982e-05*/
      -2.25082715085, /*0.000603342281596*/
      1.71594420306e-05, /*2.73369422978e-06*/
      0.000466436022016, /*3.10876388484e-05*/
      0.00487549671107, /*0.000249015033389*/
      0.0210161773176
    }; /*0.00178410485142};*/

    //0.0676506993903*/

    //100.
    static constexpr double sarcevic_std_nue[] = {
      0.0938106731958, /*0.00968084556496*/
      0.127447604594, /*0.0104518800441*/
      0.0690640747739, /*0.00928723182021*/
      -0.00202883752944, /*0.00733072332339*/
      -0.00740634544259, /*0.00464412823313*/
      -6.44177730199, /*0.00191354750982*/
      -2.1934999156e-05, /*2.35828231602e-07*/
      -0.00060750673143, /*2.12405966392e-06*/
      0.00516577885096, /*1.54299379371e-05*/
      -0.0442845804904, /*9.88758722982e-05*/
      -2.25082715085, /*0.000603342281596*/
      1.71594420306e-05, /*2.73369422978e-06*/
      0.000466436022016, /*3.10876388484e-05*/
      0.00487549671107, /*0.000249015033389*/
      0.0210161773176 /*0.00178410485142};*/
    };
    //0.0676506993903*/

    //100
    static constexpr double sarcevic_std_anumu[] = {
      0.0938106731958, /*0.00968084556496*/
      0.127447604594, /*0.0104518800441*/
      0.0690640747739, /*0.00928723182021*/
      -0.00202883752944, /*0.00733072332339*/
      -0.00740634544259, /*0.00464412823313*/
      -6.44177730199, /*0.00191354750982*/
      -2.1934999156e-05, /*2.35828231602e-07*/
      -0.00060750673143, /*2.12405966392e-06*/
      0.00516577885096, /*1.54299379371e-05*/
      -0.0442845804904, /*9.88758722982e-05*/
      -2.25082715085, /*0.000603342281596*/
      1.71594420306e-05, /*2.73369422978e-06*/
      0.000466436022016, /*3.10876388484e-05*/
      0.00487549671107, /*0.000249015033389*/
      0.0210161773176
    }; /*0.00178410485142}; */

    //0.0676506993903

    //100
    static constexpr double sarcevic_std_numu[] = {
      0.0938106731958, /*0.00968084556496*/
      0.127447604594, /*0.0104518800441*/
      0.0690640747739, /*0.00928723182021*/
      -0.00202883752944, /*0.00733072332339*/
      -0.00740634544259, /*0.00464412823313*/
      -6.44177730199, /*0.00191354750982*/
      -2.1934999156e-05, /*2.35828231602e-07*/
      -0.00060750673143, /*2.12405966392e-06*/
      0.00516577885096, /*1.54299379371e-05*/
      -0.0442845804904, /*9.88758722982e-05*/
      -2.25082715085, /*0.000603342281596*/
      1.71594420306e-05, /*2.73369422978e-06*/
      0.000466436022016, /*3.10876388484e-05*/
      0.00487549671107, /*0.000249015033389*/
      0.0210161773176
    }; /*0.00178410485142}; */

//0.0676506993903


    //double costheta = cos(zenith_rad);
    double lE = loge;
    double prompt;

    double energy_gev = pow(10, loge );

    if (energy_gev <= e_joint_nu_prompt) {return 0.;}

    if (costheta >= 0.) {costheta = -costheta;}

    switch (pdg_type)
    {
    case 12:
      prompt = fluxFunctionPrompt(sarcevic_std_nue, lE, costheta);
      break;
    case -12:
      prompt = fluxFunctionPrompt(sarcevic_std_anue, lE, costheta);
      break;
    case 14:
      prompt = fluxFunctionPrompt(sarcevic_std_numu, lE, costheta);
      break;
    case -14:
      prompt = fluxFunctionPrompt(sarcevic_std_anumu, lE, costheta);
      break;
    default:
      prompt = 0.;
      break;
    }

    if (apply_knee_correction) prompt *= correctKnee( lE );
    return prompt;
  }


  double fluxFunctionPrompt(const double *function, double lE, double costheta) const
  {
    return 1e4 * pow(10., function[0] * pow(costheta, 5.)
                     + function[1] * pow(costheta, 4.)
                     + function[2] * pow(costheta, 3.)
                     + function[3] * pow(costheta, 2.)
                     + function[4] * costheta
                     + function[5]
                     + function[6] * pow(lE, 5.)
                     + function[7] * pow(lE, 4.)
                     + function[8] * pow(lE, 3.)
                     + function[9] * pow(lE, 2.)
                     + function[10] * lE
                     + function[11] * pow(lE, 4.) * costheta
                     + function[12] * pow(lE, 3.) * pow(costheta, 2.)
                     + function[13] * pow(lE, 2.) * pow(costheta, 3.)
                     + function[14] * lE * pow(costheta, 4.));
  }

  void Print() {
    printf("prompt_Sarsevic");
  }

  ClassDef( Flux_prompt_Sarsevic_std , 1 )
};


/*! Official KM3NeT atmospheric flux model, consisting of Honda2006 for convential
    neutrinos + Endberg,Sarcevic et al. for the prompt component.
    A correction is applied that has to do with the cosmic ray knee.
    (implementation by R. Coniglione and T. Heid *)

    usage example (c++)

\code
    #include "Flux.hh"

    Flux_Atmospheric flux;

    # atm. nu_mu (14) flux at horizon (cos-zenith=0) at 1 TeV (loge=3)
    double f = flux.dNdEdOmega( 14, 3, 0 );

    # total flux, integrated over 4pi
    double ftot = flux.dNdEdOmega( 14, 3  );

\endcode
 */
struct Flux_Atmospheric : public DiffuseFlux
{
  virtual double dNdEdOmega( int pdg_type, double loge, double costheta ) const
  {
    static Flux_Honda2006 conv;
    static Flux_prompt_Sarsevic_std prompt;
    return conv.dNdEdOmega( pdg_type , loge , costheta ) + prompt.dNdEdOmega( pdg_type , loge , costheta );
  }
  void Print() {
    printf("Official KM3NeT atmospheric flux model");
  }

  ClassDef( Flux_Atmospheric, 1 );
};







struct Flux_Bartol2004 : public DiffuseFlux
{
  map< string, TH2F > fluxes;

  string get_type( int pdg_type )
  {
    switch ( pdg_type )
    {
    case 12  : return "nue";
    case 14  : return "num";
    case 16  : return "nut";
    case -12 : return "nbe";
    case -14 : return "nbm";
    case -16 : return "nbt";
    default :
      cout << " invalid neutrino type " << pdg_type << endl;
    }
    return "";
  }

  Flux_Bartol2004 ( bool wget_files )
  {
    if (wget_files) retreive_files();

    // the following come from manual inspection of the files()
    const double logemin    = -1;
    const double logemax    = 4;
    const int nbins_loge    = int ( logemax - logemin ) * 10; // 10 bins/decade
    const double ctmin      = -1;
    const double ctmax      = 1;
    const int nbins_ct      = int (ctmax - ctmin ) * 10;

		vector<string> types{"nue", "nbe", "num", "nbm"};
		for( auto typeit = types.begin(); typeit != types.end(); typeit++) {
			string type = *typeit;
      fluxes[ type ] = TH2F( Form("flux_%s", type.c_str() ),
                             Form("flux of %s [ /m**2/sr/sec ]", type.c_str() ),
                             nbins_loge, logemin, logemax,
                             nbins_ct  , ctmin,   ctmax );

      read_file( fluxes[ type ], "fmin10_i0403z.kam_" + type );
      read_file( fluxes[ type ], "f210_3_z.kam_" + type      );
    }
  }


  TH2F* get_2dhist( string type )
  {
    return &fluxes[type];
  }

  double getflux( string type, double loge, double costheta )
  {
    return fluxes[type].Interpolate( loge, costheta );
  }



  virtual double dNdEdOmega( int pdg_type, double loge, double costheta )
  {
    if ( pdg_type == 16 || pdg_type == -16 ) return 0;
    return getflux_dnde( get_type( pdg_type ), loge, costheta );
  }

  double getflux_dnde( string type, double loge, double costheta )
  {
    const double e = pow(10, loge);
    return getflux( type, loge, costheta ) / e;
  }




  void retreive_files()
  {
    // see http://www-pnp.physics.ox.ac.uk/~barr/fluxfiles/0403i/index.html
    // use kamioka solar min and follow link to 'newer version of high-energy flux'

    system("wget -nc http://www-pnp.physics.ox.ac.uk/~barr/fluxfiles/0403i/fmin10_i0403z.kam_nue");
    system("wget -nc http://www-pnp.physics.ox.ac.uk/~barr/fluxfiles/0403i/fmin10_i0403z.kam_nbe");
    system("wget -nc http://www-pnp.physics.ox.ac.uk/~barr/fluxfiles/0403i/fmin10_i0403z.kam_num");
    system("wget -nc http://www-pnp.physics.ox.ac.uk/~barr/fluxfiles/0403i/fmin10_i0403z.kam_nbm");

    system("wget -nc http://www-pnp.physics.ox.ac.uk/~barr/fluxfiles/0408i/f210_3_z.kam_nue");
    system("wget -nc http://www-pnp.physics.ox.ac.uk/~barr/fluxfiles/0408i/f210_3_z.kam_nbe");
    system("wget -nc http://www-pnp.physics.ox.ac.uk/~barr/fluxfiles/0408i/f210_3_z.kam_num");
    system("wget -nc http://www-pnp.physics.ox.ac.uk/~barr/fluxfiles/0408i/f210_3_z.kam_nbm");
  }

  // from http://www-pnp.physics.ox.ac.uk/~barr/fluxfiles/
  // 1. Energy at centre of bin in GeV
  // 2. Zenith files: Cos(zenith) at the centre of the bin
  //    Azimuth files: Azimuth at centre of bin (degrees)
  // 3. Flux dN/dln E in /m**2/steradian/sec...  = E dN/de
  // 4. MC statistical error on the flux
  // 5. Number of unweighted events entering in the bin (not too useful)

  void read_file( TH2& h, string filename )
  {
    ifstream f( filename.c_str() );
    if (!f )
    {
      printf("failed to open file %s", filename.c_str() );
      printf("trying to wget files from the internet");
      exit(1);
    }

    printf("reading file %s", filename.c_str());
		string line;

    while ( f >> line )
    {
			if ( line[0] == '#'  ) continue;
			istringstream iss(line);
			double e, cos_theta, flux;
			iss >> e >> cos_theta >> flux;
			
      int bin = h.FindBin( log10(e), cos_theta );
      h.SetBinContent( bin , flux );
    }
  }
  void Print() {
    printf("Bartol2004");
  }

  ClassDef( Flux_Bartol2004 , 1 )
};

#endif