#ifndef __JAANET__JVOLUME__
#define __JAANET__JVOLUME__

#include <cmath>

#include "TAxis.h"

#include "km3net-dataformat/offline/Head.hh"

#include "JTools/JRange.hh"

#include "JAAnet/JHead.hh"


/**
 * \author mdejong
 */

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

namespace JAANET {

  /**
   * Auxiliary class for histogramming of effective volume.
   * In this, it is assumed that events are generated according \f$ \frac{dN}{dE} \propto E^{\alpha} \f$.
   * The weight is expressed in \f$ \mathrm{km}^{3} \f$.
   */
  struct JVolume {
    /**
     * Constructor.
     *
     * \param  head          Monte Carlo run header
     * \param  Elog10        application of log10 to energy
     */
    JVolume(const Head& head, 
	    const bool  Elog10 = false) :
      elog  (Elog10),
      E(0.0, 0.0),
      alpha (0.0),
      Wall  (1.0)
    {
      const JHead buffer(head);

      if        (buffer.is_valid(&JHead::spectrum) && 
		 buffer.is_valid(&JHead::cut_nu)   && 
		 buffer.is_valid(&JHead::genvol)) {

	alpha  = buffer.spectrum.alpha;
	E      = buffer.cut_nu.E;
	Wall   = (1.0e-9 * (getW(E.getUpperLimit()) -
			    getW(E.getLowerLimit())) *
		  buffer.genvol.volume / buffer.genvol.numberOfEvents);    // km^3

      } else if (buffer.is_valid(&JHead::cut_in) &&
		 buffer.is_valid(&JHead::livetime)) {

	E      = buffer.cut_in.E;
	Wall   = 1.0 / buffer.livetime.numberOfSeconds;                     // Hz
      }
    }


    /**
     * Get spectral index of energy distribution.
     *
     * \return               spectral index
     */
    inline double getAlpha() const
    {
      return alpha;
    }


    /**
     * Get generation dependent weight.
     *
     * \return               weight
     */
    inline double getWall() const
    {
      return Wall;
    }


    /**
     * Get minimal abscissa value.
     *
     * \return               abscissa value
     */
    inline Double_t getXmin() const
    {
      return getX(E.getLowerLimit());
    }


    /**
     * Get maximal abscissa value.
     *
     * \return               abscissa value
     */
    inline Double_t getXmax() const
    {
      return getX(E.getUpperLimit());
    }


    /**
     * Get abscissa value.
     *
     * \param  E             energy
     * \param  constrain     constrain
     * \return               abscissa value
     */
    inline Double_t getX(const Double_t E, double constrain = false) const
    {
      double x = (elog ? log10(E) : E);

      if (constrain) {
        if        (x < getXmin()) {
          return getXmin();
        } else if (x > getXmax()) {
          return getXmax();
        }
      }

      return x;
    }


    /**
     * Get energy.
     *
     * \param  x             abscissa value
     * \param  constrain     constrain
     * \return               energy
     */
    inline Double_t getE(const Double_t x, double constrain = false) const
    {
      const double Ex = (elog ? pow(10.0, x) : x);

      if (constrain)
	return E.constrain(Ex);
      else
	return Ex;
    }


    /**
     * Get bin width corrected energy spectrum dependent weight.
     *
     * \param  axis          axis
     * \param  E             energy
     * \return               weight [km^3]
     */
    inline Double_t getW(TAxis* axis, const Double_t E) const
    {
      const Int_t    index = axis->FindBin(getX(E));

      const Double_t xmin  = axis->GetBinLowEdge(index);
      const Double_t xmax  = axis->GetBinUpEdge (index);
      
      const Double_t Wmin  = getW(getE(xmin));
      const Double_t Wmax  = getW(getE(xmax));

      return Wall / (Wmax - Wmin);
    }


    /**
     * Get generation dependent integral value of given energy.
     *
     * \param  E             energy
     * \return               integral value
     */
    inline double getW(const double E) const
    {
      if (alpha != -1.0)
	return pow(E, 1.0 + alpha) / (1.0 + alpha);
      else
	return log(E);
    }


    /**
     * Multiply weight.
     *
     * \param  factor        factor
     * \return               this volume
     */
    JVolume& mul(const double factor) 
    {
      Wall *= factor;

      return *this;
    }


    /**
     * Divide weight.
     *
     * \param  factor        factor
     * \return               this volume
     */
    JVolume& div(const double factor) 
    {
      Wall /= factor;

      return *this;
    }


  protected:
    bool                   elog;     //!< histogram option
    JTOOLS::JRange<double> E;        //!< Energy range [GeV]
    double                 alpha;    //!< spectral index
    double                 Wall;     //!< generation volume
  };
}

#endif