#ifndef __JSIRENE__JSIRENETOOLKIT__
#define __JSIRENE__JSIRENETOOLKIT__

#include <vector>
#include <cmath>

#include "TRandom3.h"

#include "JPhysics/JConstants.hh"
#include "JPhysics/JPDFTypes.hh"
#include "JPhysics/JGeane.hh"
#include "JAAnet/JAAnetToolkit.hh"
#include "JGeometry3D/JPosition3D.hh"
#include "JGeometry3D/JTime.hh"
#include "JGeometry3D/JVersor3Z.hh"
#include "JGeometry2D/JVector2D.hh"
#include "JLang/JComparator.hh"


/**
 * \author mdejong
 */


/**
 * Detector simulations.
 */
namespace JSIRENE {}
namespace JPP { using namespace JSIRENE; }

namespace JSIRENE {

  using JPHYSICS::JGeane;
  using JPHYSICS::gWater;
  using JPHYSICS::MASS_MUON;
  using JPHYSICS::JPDFType_t;
  using JPHYSICS::getSinThetaC;
  using JPHYSICS::getSpeedOfLight;
  
  using JGEOMETRY3D::JPosition3D;
  using JGEOMETRY3D::JTime;
  using JGEOMETRY3D::JVersor3Z;
  using JGEOMETRY2D::JVector2D;
  
  using JAANET::JHitType_t;


  /**
   * Get number of photo-electrons of a hit given the expectation 
   * values of the number of photo-electrons on a module and PMT.
   *
   * The return value is evaluated by pick-and-drop statistics from
   * the generated number of photo-electrons when the expectation 
   * value of the number of photo-electrons on a module deviates 
   * less than 5 sigmas from 0 (i.e.\ when it is less than 25).
   * Otherwise, the return value is evaluated by Poisson statistics
   * from the expectation value of the number of photo-electrons on PMT.
   *
   * \param  NPE       expectation value of npe on module
   * \param  N         generated   value of npe on module
   * \param  npe       expectation value of npe on PMT
   * \return           number of photo-electrons on PMT
   */
  inline int getNumberOfPhotoElectrons(const double NPE,
				       const int    N,
				       const double npe)
  {
    if (NPE < 25.0) {
      
      int n = 0;
			    
      for (int i = N; i != 0; --i) {
	if (gRandom->Rndm() * NPE < npe) {
	  ++n;
	}
      }
      
      return n;

    } else {
      
      return gRandom->Poisson(npe);
    }
  }


  /**
   * Get number of photo-electrons of a hit given number of photo-electrons on PMT.
   *
   * The number of photo-electrons of a hit may be larger than unity
   * to limit the overall number of hits and consequently the memory coonsumption as well as 
   * the number of times the arrival time needs to be evaluated which is CPU intensive.
   *
   * \param  npe       number of photo-electrons on PMT
   * \return           number of photo-electrons of hit
   */
  inline int getNumberOfPhotoElectrons(const int npe)
  {
    static const int NPE = 20;
    
    if (npe > NPE) {
      
      const int n = (int) (NPE * log10((double) npe / (double) NPE));
    
      if      (n == 0)
	return 1;
      else
	return n;
    }

    return 1;
  }


  /**
   * Get hit type corresponding to given PDF type.
   *
   * \param  pdf       PDF type
   * \param  shower    force origin from neutrino interaction shower
   * \return           hit type
   */
  inline JHitType_t getHitType(const JPDFType_t pdf, const bool shower = false)
  {
    using namespace JAANET;
    using namespace JPHYSICS;
    
    switch (pdf) {

    case DIRECT_LIGHT_FROM_MUON:
      return HIT_TYPE_MUON_DIRECT;

    case SCATTERED_LIGHT_FROM_MUON:
      return HIT_TYPE_MUON_SCATTERED;

    case DIRECT_LIGHT_FROM_DELTARAYS:
      return HIT_TYPE_DELTARAYS_DIRECT;

    case SCATTERED_LIGHT_FROM_DELTARAYS:
      return HIT_TYPE_DELTARAYS_SCATTERED;

    case DIRECT_LIGHT_FROM_EMSHOWER:
      return (shower ? HIT_TYPE_SHOWER_DIRECT : HIT_TYPE_BREMSSTRAHLUNG_DIRECT );

    case SCATTERED_LIGHT_FROM_EMSHOWER:
      return (shower ? HIT_TYPE_SHOWER_SCATTERED : HIT_TYPE_BREMSSTRAHLUNG_SCATTERED);

    default:
      return HIT_TYPE_UNKNOWN;      
    }    
  }
  

  /**
   * Point along muon trajectory.
   */
  struct JPoint :
    JPosition3D,
    JTime
  {
    /**
     * Default constructor.
     */
    JPoint() :
      JPosition3D(),
      JTime(),
      __E (0.0),
      __Es(0.0)
    {}


    /**
     * Constructor.
     *
     * \param  z         position [m]
     * \param  t         time     [ns]
     * \param  E         energy   [GeV]
     */
    JPoint(const double z,
	   const double t,
	   const double E) :
      JPosition3D(0.0, 0.0, z),
      JTime(t),
      __E (E),
      __Es(0.0)
    {}


    /**
     * Get muon energy.
     *
     * \return           energy [GeV]
     */
    double getE() const
    {
      return __E;
    }


    /**
     * Get shower energy.
     *
     * \return           energy [GeV]
     */
    double getEs() const
    {
      return __Es;
    }

  protected:
    double __E;
    double __Es;
  };


  /**
   * Muon trajectory.
   */
  struct JTrack :
    public std::vector<JPoint>
  {
    /**
     * Constructor.
     *
     * \param  point     muon starting point
     */
    JTrack(const JPoint& point) :
      std::vector<JPoint>(1, point)
    {}


    /**
     * Get muon energy at given position along trajectory.
     *
     * \param  z         position [m]
     * \return           energy [GeV]
     */
    double getE(const double z) const
    {
       if (!this->empty()) {
	
	const_iterator p = lower_bound(this->begin(), this->end(), z, JLANG::make_comparator(&JPoint::getZ));

	if (p != this->end() && p != this->begin()) {
	  
	  --p;
	
	  return p->getE() -  (z - p->getZ()) * gWater.getA();
	}
      }

      return MASS_MUON;
    }


    /**
     * Get muon position at given position along trajectory.
     *
     * \param  z         position [m]
     * \return           position
     */
    JVector2D getPosition(const double z) const
    {
      using namespace JPP;

      const double precision = 1.0e-2;

      if (!this->empty()) {
	
	const_iterator p = lower_bound(this->begin(), this->end(), z, JLANG::make_comparator(&JPoint::getZ));

	if (p == this->end()) { 
	  --p;
	}

	if (p == this->begin()) {

	  return JVector2D(p->getX(),
			   p->getY());

	} else {

	  JVector3D pos;

	  pos  = p->getPosition();

	  --p;

	  pos -= p->getPosition();

	  const double u = (pos.getZ() > precision ? (z - p->getZ()) / pos.getZ() : 0.0);

	  return JVector2D(p->getX() + u * pos.getX(),
			   p->getY() + u * pos.getY());
	}
      }

      return JVector2D(0.0, 0.0);
    }
  };


  /**
   * Vertex of energy loss of muon.
   */
  struct JVertex :
    JPoint,
    JVersor3Z
  {
    /**
     * Default constructor.
     */
    JVertex() :
      JPoint(),
      JVersor3Z()
    {}


    /**
     * Constructor.
     *
     * \param  z         position [m]
     * \param  t         time     [ns]
     * \param  E         energy   [GeV]
     */
    JVertex(const double z,
	    const double t,
	    const double E) :
      JPoint(z, t, E),
      JVersor3Z()
    {}


    /**
     * Get visible range of muon using default ionisation energy loss.
     * This method applies only ionisation energy loss.
     *
     * \return           range [m]
     */
    double getRange() const
    {
      return getRange(gWater.getA());
    }


    /**
     * Get visible range of muon using given ionisation energy loss.
     * This method applies only ionisation energy loss.
     *
     * \param  a         ionization parameter [GeV/m]
     * \return           range [m]
     */
    double getRange(double a) const
    {
      static const double Ethreshold = MASS_MUON / getSinThetaC();
      
      if (__E > Ethreshold)
        return (__E - Ethreshold) / a;
      else
        return 0.0;
    }
    

    /**
     * Get range of muon using given energy loss function.
     *
     * \param  geane     energy loss function
     * \return           range [m]
     */
    double getRange(const JGeane& geane) const
    {
      return geane(__E);
    }


    /**
     * Apply shower energy loss.
     *
     * \param  Ts        scattering angles
     * \param  Es        shower energy [GeV]
     */
    void applyEloss(const JVersor3Z& Ts, const double Es) 
    {
      static_cast<JVersor3Z&>(*this).add(Ts);

      this->__E  -= Es;
      this->__Es  = Es;
    }


    /**
     * Step using default ionisation energy loss.
     * This method applies only ionisation energy loss.
     *
     * \param  ds        step [m]
     * \return           this vertex
     */
    JVertex& step(const double ds)
    {
      return step(gWater.getA(), ds);
    }


    /**
     * Step using given ionisation energy loss.
     * This method applies only given ionisation energy loss.
     *
     * \param  a         ionization parameter [GeV/m]
     * \param  ds        step [m]
     * \return           this vertex
     */
    JVertex& step(const double a, const double ds)
    {
      __x += this->getDX() * ds;
      __y += this->getDY() * ds;
      __z += this->getDZ() * ds;
      __t += ds / getSpeedOfLight();
      __E -= ds * a;

      return *this;
    }


    /**
     * Step using given energy loss function.
     *
     * \param  geane     energy loss function
     * \param  ds        step [m]
     * \return           this vertex
     */
    JVertex& step(const JGeane& geane, const double ds)
    {
      __x += this->getDX() * ds;
      __y += this->getDY() * ds;
      __z += this->getDZ() * ds;
      __t += ds / getSpeedOfLight();
      __E  = geane(__E, ds);

      return *this;
    }
  };
}

#endif