#ifndef __JFIT__JNPE__
#define __JFIT__JNPE__

#include <cmath>

#include "JPhysics/JConstants.hh"
#include "JPhysics/JGeane.hh"
#include "JPhysics/JPDFToolkit.hh"

#include "JFit/JK40.hh"
#include "JFit/JFitToolkit.hh"


/**
 * \author mdejong
 */

namespace JFIT {}
namespace JPP { using namespace JFIT; }

namespace JFIT {
  
  /**
   * Auxiliary class for handling various light yields.
   *
   * Note that the effective light yield due to bremsstrahlung is proportional to the muon energy.\n
   * The position along the muon path is used to correct for the energy loss of the muon 
   * between the reference position (z = 0) and the point of emission of the light. 
   */
  struct JNPE :
    public JK40
  {
    /**
     * Default constructor.
     */
    JNPE() :
      JK40(),
      y1(0.0),
      yA(0.0),
      yB(0.0),
      z (0.0)
    {}


    /**
     * Constructor.
     *
     * \param  y0                 light yield due to random background         [npe]
     * \param  y1                 light yield due to minimum ionizing particle [npe]
     * \param  yA                 light yield due to delta-rays                [npe*m/GeV]
     * \param  yB                 light yield due to bremsstrahlung            [npe/GeV]
     * \param  z                  position along muon path [m]
     */
    JNPE(const double y0,
	 const double y1,
	 const double yA,
	 const double yB,
	 const double z) :
      JK40(y0),
      y1(y1),
      yA(yA),
      yB(yB),
      z (z)
    {}


    /**
     * Get light yield due to minimum ionizing particle.
     *
     * \return                    light yield [npe]
     */
    double getY1() const
    {
      return y1;
    }


    /**
     * Get light yield due to delta-rays
     *
     * \return                    light yield [npe*m/GeV]
     */
    double getYA() const
    {
      return yA;
    }


    /**
     * Get light yield due to bremsstrahlung.
     *
     * \return                    light yield [npe/GeV]
     */
    double getYB() const
    {
      return yB;
    }


    /**
     * Get position along muon path.
     *
     * \return                    position along muon path [m]
     */
    double getZ() const
    {
      return z;
    }

    
    /**
     * Expected number of photo-electrons for muon hypothesis as a function of muon energy.
     *
     * \param  E_GeV              energy [GeV]
     * \return                    light yield [npe]
     */
    inline double getH1(const double E_GeV) const
    {
      using namespace JPP;

      const double E = gWater.getE(E_GeV, this->getZ());

      if (E >= MASS_MUON * INDEX_OF_REFRACTION_WATER)
	return (this->getY1()                            + 
		this->getYA() * getDeltaRaysFromMuon(E)  +
		this->getYB() * E);
      else
	return 0.0;
    }


    /**
     * Get probability for observing a hit or not as a function of muon energy.
     *
     * \param  E_GeV              energy [GeV]
     * \param  hit                hit
     * \return                    probability
     */
    double getP(const double E_GeV, const bool hit) const
    {
      return JFIT::getP(this->getH1(E_GeV) + this->getH0(), hit);
    }


    /**
     * Get chi2 for observing a hit or not as a function of muon energy.
     *
     * \param  E_GeV              energy [GeV]
     * \param  hit                hit
     * \return                    probability
     */
    double getChi2(const double E_GeV, const bool hit) const
    {
      return JFIT::getChi2(this->getH1(E_GeV) + this->getH0(), hit);
    }

  protected:
    double y1;      //!< light yield due to minimum ionizing particle [npe]
    double yA;      //!< light yield due to delta-rays                [npe*m/GeV]
    double yB;      //!< light yield due to bremsstrahlung            [npe/GeV]
    double z;       //!< position along muon path [m]
  };
}

#endif