#ifndef __JFIT__JSHOWERNPE__
#define __JFIT__JSHOWERNPE__

#include <cmath>

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

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


/**
 * \author mdejong, adomi
 */

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

namespace JFIT {
  
  /**
   * Auxiliary class for handling EM shower light yield.
   */
  struct JShowerNPE :
    public JK40
  {
    /**
     * Default constructor.
     */
    JShowerNPE() :
      JK40(),
      y1(0.0)
    {}


    /**
     * Constructor.
     *
     * \param  y0                 light yield due to random background         [npe]
     * \param  y1                 light yield due to shower                    [npe/GeV]
     */
    JShowerNPE(const double y0,
	       const double y1) :
      JK40(y0),
      y1(y1)
    {}


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

    
    /**
     * Expected number of photo-electrons for shower hypothesis as a function of shower energy.
     *
     * \param  E_GeV              energy [GeV]
     * \return                    light yield [npe]
     */
    inline double getH1(const double E_GeV) const
    {
      return E_GeV * this->getY1();
    }


    /**
     * Get probability for observing a hit or not as a function of shower 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 shower 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 shower                    [npe/GeV]
  };
}

#endif