#ifndef __JSHOWEREH__
#define __JSHOWEREH__

#include <istream>
#include <ostream>

#include "JMath/JMath.hh"
#include "JFit/JShower3Z.hh"
#include "JFit/JEnergy.hh"


/**
 * \author adomi
 */

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

namespace JFIT {

  using JMATH::JMath;


  /**
   * Data structure for fit of straight line in positive z-direction with energy.
   * Note that the position coordinates are defined with respect to the given direction.
   */
  class JShowerEH : 
    public JShower3Z,
    public JEnergy,
    public JMath<JShowerEH>
  {
  public:

    /**
     * Default constructor.
     */
    JShowerEH() :
      JShower3Z(),
      __Eem(0.0),
      __Eh(0.0),
      __By(0.0)
    {}


    /**
     * Constructor.
     *
     * \param  line             line
     * \param  x                EM shower energy
     * \param  h                Hadronic shower energy  
     * \param  By               bjorken y
     */
    JShowerEH(const JShower3Z& line, const double x, const double h, const double By = 0.0) :
      JShower3Z(line),
      __Eem(x),
      __Eh(h),
      __By(By)
    {}

    /**
     * Constructor.
     *
     * \param  point          point
     * \param  dir            direction
     * \param  x              EM shower energy
     * \param  h              Hadronic shower energy   
     * \param  By             bjorken y
     */
    JShowerEH(const JPoint4D&  point,
	      const JVersor3Z& dir,
	      const double     x,
	      const double     h,
	      const double     By = 0.0):
      JShower3Z(JPoint4D(point), JVersor3Z(dir)),
      __Eem(x),
      __Eh(h),
      __By(By)
    {}


     /**
      * Get EM Energy.
      *
      * \return               EM Energy
      */
    double getEMEnergy() const
     {
       return __Eem;
     }

    /**
      * Get Hadronic Energy.
      *
      * \return               Hadronic Energy
      */
    double getHEnergy() const
    {
       return __Eh;
     }

    /**
     * Get bjorken y
     *
     * \return                  bjorken y
     */
    double getBy() const
    {
      return __By;
    }

    /**
     * Prefix unary minus.
     *
     * \return                  shower
     */
    JShowerEH& negate()
    {
      JShower3Z::negate();
      __Eem = -__Eem;
      __Eh  = -__Eh;
      __By  = -__By;

      return *this;
    }

    /**
     * Addition operator.
     *
     * \param  value            shower
     * \return                  shower
     */
    JShowerEH& add(const JShowerEH& value)
    {
      JShower3Z::add(value);
      __Eem += value.getEMEnergy();
      __Eh  += value.getHEnergy();
      __By  += value.getBy();

      return *this;
    }


    /**
     * Subtraction operator.
     *
     * \param  value            shower
     * \return                  shower
     */
    JShowerEH& sub(const JShowerEH& value)
    {
      JShower3Z::sub(value);
      __Eem -= value.getEMEnergy();
      __Eh  -= value.getHEnergy();
      __By  -= value.getBy();

      return *this;
    }


    /**
     * Multiplication operator.
     *
     * \param  value            multiplication factor
     * \return                  shower
     */
    JShowerEH& mul(const double value)
    {
      JShower3Z::mul(value);
      __Eem *= value;
      __Eh  *= value;
      __By  *= value;

      return *this;
    }


    /**
     * Division operator.
     *
     * \param  value            division      factor
     * \return                  shower
     */
    JShowerEH& div(const double value)
    {
      JShower3Z::div(value);
      __Eem /= value;
      __Eh  /= value;
      __By  /= value;

      return *this;
    }

    /**
     * Get EM energy.
     *
     * \return                  EM Energy [log(E/GeV)]
     */
    double getLogEem() const
    {
      return __Eem;
    }
    
    /**
     * Get EM energy.
     *
     * \return                  EM Energy [GeV]
     */
    double getEem() const
    {
      return pow(10.0, __Eem);
    }

    /**
     * Put EM energy.
     *
     * \param  E                Energy [GeV]
     */
    void putEem(const double E)
    {
      __Eem = log10(E);
    }

    /**
     * Get derivative of energy.
     *
     * \return                  dE/dx [GeV]
     */
    double getDEem() const
    {
      return getEem() * log(10.0);
    }

    /**
     * Get Hadronic energy.
     *
     * \return                  H Energy [log(E/GeV)]
     */
    double getLogEh() const
    {
      return __Eh;
    }
    
    /**
     * Get Hadronic energy.
     *
     * \return                  H Energy [GeV]
     */
    double getEh() const
    {
      return pow(10.0, __Eh);
    }

    /**
     * Put Hadronic energy.
     *
     * \param  E                Energy [GeV]
     */
    void putEh(const double E)
    {
      __Eh = log10(E);
    }

    /**
     * Get derivative of energy.
     *
     * \return                  dE/dx [GeV]
     */
    double getDEh() const
    {
      return getEh() * log(10.0);
    }

    
    /**
     * Read object from input.
     *
     * \param  in            input stream
     * \param  object        object
     * \return               input stream
     */
    friend inline std::istream& operator>>(std::istream& in, JShowerEH& object)
    {
      in >> static_cast<JShower3Z&>(object);
      in >> object.__Eem;
      in >> object.__Eh;
      in >> object.__By;

      return in;
    }


    /**
     * Write object to output.
     *
     * \param  out           output stream
     * \param  object        object
     * \return               output stream
     */
    friend inline std::ostream& operator<<(std::ostream& out, const JShowerEH& object)
    {
      out << static_cast<const JShower3Z&>(object);
      out << object.__Eem;
      out << object.__Eh;
      out << object.__By;

      return out;
    }

    typedef double JShowerEH::*parameter_type;

    static parameter_type pEem() { return &JShowerEH::__Eem; }
    static parameter_type pEh()  { return &JShowerEH::__Eh; }
    static parameter_type pBy()  { return &JShowerEH::__By; }

  protected:
    double __Eem;
    double __Eh;
    double __By;
  };
}

#endif