#ifndef __JAANET__JDIFFUSEFLUX__
#define __JAANET__JDIFFUSEFLUX__

#include <iostream>

#include "km3net-dataformat/offline/Evt.hh"
#include "km3net-dataformat/offline/Trk.hh"

#include "JAAnet/JFlux.hh"
#include "JAAnet/JAAnetToolkit.hh"

#include "flux/Flux.hh"


/**
 * \author bjung
 */

namespace JAANET {}
namespace JPP { using namespace JAANET; }

namespace JAANET {

  /**
   * Low-level interface for diffuse fluxes.
   */
  struct JDiffuseFlux :
    public JFlux,
    public DiffuseFlux
  {
    /**
     * Virtual destructor.
     */
    virtual ~JDiffuseFlux()
    {}
    

    /**
     * Get event-weight factor for given event.
     *
     * \param  evt                 event 
     * \return                     event-weight factor [GeV^-1 * m^-2 * sr^-1 * s^-1]
     */
    double getFactor(const Evt& evt) const override
    {
      const Trk& neutrino = get_neutrino(evt);
      const double log10E = log10(neutrino.E);
      const double costh  = neutrino.dir.z / neutrino.dir.len();
      
      return dNdEdOmega(neutrino.type, log10E, costh);
    }


    /**
     * Get event-weight factor for given particle PDG-identifier, energy and zenith-angle.
     *
     * \param  type              PDG particle type
     * \param  log10E            logarithmic neutrino energy [GeV]
     * \param  costh             cosine zenith angle
     * \return                   diffuse flux [GeV^-1 * m^-2 * sr^-1 * s^-1]
     */
    double getFactor(const int    type,
		     const double log10E,
		     const double costh) const
    {
      return dNdEdOmega(type, log10E, costh);
    }


    /**
     * Get event-weight factor for given particle PDG-identifier, energy and zenith-angle.
     *
     * \param  type              PDG particle type
     * \param  log10E            logarithmic neutrino energy [GeV]
     * \param  costh             cosine zenith angle
     * \return                   diffuse flux [GeV^-1 * m^-2 * sr^-1 * s^-1]
     */
    double getFlux(const int    type,
		   const double log10E,
		   const double costh) const
    {
      return dNdEdOmega(type, log10E, costh);
    }


    /**
     * Print flux information.
     */
    virtual void Print() const override
    {
      prnt(std::cout);
    }
  };
}

#endif