#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"


/**
 * \file
 * Toolkit for JSirene.
 *
 * \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;


  /**
   * Auxiliary data structure for determination of number of photo-electrons.
   */
  const struct number_of_photo_electrons_type {
    /**
     * Get numbers of photo-electrons for PMTs given the expectation 
     * values of the number of photo-electrons on a module and PMTs.
     *
     * The number of photo-electrons on the PMTs are evaluated by pick-and-drop statistics 
     * from the generated number of photo-electrons when the option is set to <tt>true</tt>.
     * Otherwise, the numbers are directly evaluated using Poisson statistics
     * from the expectation values of the number of photo-electrons on the PMTs.
     *
     * It is recommended to set the option to <tt>true</tt> when the expectation value
     * of number of photo-electrons in the module is less than 25 or so,
     * which corresponds to a 5-sigma probability for zero photo-electrons in the module.
     *
     * \param  NPE       expectation value of number of photo-electrons in module
     * \param  N         generated number of photo-electrons in module
     * \param  npe       list of expectation values of number of photo-electrons in PMTs
     * \param  option    option
     * \return           list of number of photo-electrons in PMTs
     */
    const std::vector<size_t>& operator()(const double NPE, const int N, const std::vector<double>& npe, const bool option) const
    {
      using namespace std;

      ns.resize(npe.size());
      vs.resize(npe.size());

      if (option) {

	ns[0] = 0;
	vs[0] = npe[0];

	for (size_t i = 1; i != npe.size(); ++i) {
	  ns[i] = 0;
	  vs[i] = vs[i-1] + npe[i];
	}

	for (int i = N; i != 0; --i) {

	  const double x = gRandom->Rndm() * NPE;

	  if (x < *vs.rbegin()) {
	    ns[distance(vs.begin(), lower_bound(vs.begin(), vs.end(), x))] += 1;
	  }
	}

      } else {

	for (size_t i = 0; i != npe.size(); ++i) {
	  ns[i] = gRandom->Poisson(npe[i]);
	}
      }

      return ns;
    }


    /**
     * 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 size_t operator()(const size_t npe) const
    {
      static const size_t NPE = 20;
    
      if (npe > NPE) {
      
	const size_t n = (size_t) (NPE * log10((double) npe / (double) NPE));
    
	if      (n == 0)
	  return 1;
	else
	  return n;
      }

      return 1;
    }

  private:
    mutable std::vector<size_t> ns;
    mutable std::vector<double> vs;

  } getNumberOfPhotoElectrons;


  /**
   * 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