#ifndef __JAANET__JPOWERLAWFLUX__
#define __JAANET__JPOWERLAWFLUX__

#include "JLang/JEquals.hh"

#include "JAAnet/JEvtWeightFactorFunction.hh"


/**
 * \author bjung
 */

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

namespace JAANET {

  /**
   * Example function object for computing power-law flux.
   */
  struct JPowerLawFlux :
    public JLANG::JEquals<JPowerLawFlux>
  {
    /**
     * Default constructor.
     */
    JPowerLawFlux() :
      normalisation(1.0),
      spectralIndex(0.0)
    {}
    
    
    /**
     * Constructor.
     *
     * \param  normalisation         normalisation   [GeV * m^-2 * sr^-1 * s^-1]
     * \param  spectralIndex         spectral index
     */
    JPowerLawFlux(const double normalisation,
		  const double spectralIndex) :
      normalisation(normalisation),
      spectralIndex(spectralIndex)
    {}
    
    
    /**
     * Get flux of given event.
     *
     * \param  evt             event
     * \return                 flux                  [GeV * m^-2 * sr^-1 * s^-1]
     */
    double operator()(const Evt& evt) const
    {
      using namespace std;
      using namespace JPP;

      if (has_neutrino(evt)) {

	const Trk& neutrino = get_neutrino(evt);
	
	return normalisation * pow(neutrino.E, -spectralIndex);

      } else {
	
	THROW(JNullPointerException, "JPowerLawFlux::operator(): No neutrino for event " << evt.id << '.' << endl);	
      }
    }


    /**
     * Check if this flux is identical to the given flux.
     *
     * \param  object                power-law flux object
     * \return                       true if this flux is identical to given flux; else flase
     */
    bool equals(const JPowerLawFlux& object) const
    {
      return (object.normalisation == this->normalisation &&
	      object.spectralIndex == this->spectralIndex);
    }
    
    
    /**
     * Stream input.
     *
     * \param  in                     input stream
     * \param  object                 power-law flux
     * \return                        input stream
     */
    inline friend std::istream& operator>>(std::istream&  in,
					   JPowerLawFlux& object)
    {
      return in >> object.normalisation
	        >> object.spectralIndex;
    }
    
    
    /**
     * Write power-law parameters to output stream.
     *
     * \param  out                    output stream
     * \param  object                 power-law flux
     * \return                        output stream
     */
    inline friend std::ostream& operator<<(std::ostream&        out,
					   const JPowerLawFlux& object)
    {
      const JFormat format(out);
      
      out << FIXED(5,3) << object.normalisation << ' '
	  << FIXED(5,3) << object.spectralIndex;

      return out;
    }
    
    
    double normalisation;   //!< normalisation      [GeV * m^-2 * sr^-1 * s^-1]
    double spectralIndex;   //!< spectral index     >= 0
  };


  /**
   * Auxiliary method for creating an interface to a power-law flux function.
   *
   * \param  normalisation         normalisation    [GeV * m^-2 * sr^-1 * s^-1]
   * \param  spectralIndex         spectral index
   */
  inline JEvtWeightFactorFunction<JPowerLawFlux, JFlux> make_powerLawFluxFunction(const double normalisation,
										  const double spectralIndex) {

    const JPowerLawFlux flux(normalisation, spectralIndex);
    
    return make_fluxFunction(flux);
  }
}

#endif