#ifndef __JSIRENE__
#define __JSIRENE__

#include <utility>

#include "TRandom3.h"

#include "antcc/Header.hh"
#include "antcc/Track.hh" 
#include "antcc/Event.hh" 

#include "JGeometry3D/JPoint3D.hh"
#include "JGeometry3D/JVersor3D.hh"
#include "JGeometry3D/JAxis3D.hh"
#include "JDetector/JDetector.hh"
#include "JSirene/JAntcc.hh"


namespace JSIRENE {

  //typedef  unsigned char    uchar;
  //typedef  unsigned short   ushort;
  //typedef  unsigned int     uint;

  namespace {
    using JGEOMETRY3D::JPoint3D;
    using JGEOMETRY3D::JVersor3D;
    using JGEOMETRY3D::JAxis3D;
  }


  inline JPoint3D         getPosition (const Vec3D& v)    { return JPoint3D    (v.x, v.y, v.z); }
  inline JVersor3D        getDirection(const Vec3D& v)    { return JVersor3D(v.x, v.y, v.z); }

  inline JPoint3D         getPosition (const BaseTrack& track)  { return getPosition (track.position ()); }
  inline JVersor3D        getDirection(const BaseTrack& track)  { return getDirection(track.direction()); }

  inline JAxis3D getAxis(const BaseTrack& track) 
  { 
    return JAxis3D(getPosition(track), getDirection(track));
  }


  // GEANT particle codes

  inline bool is_photon  (const McTrack& track) { return track.type() ==  1 || track.type() ==  7; }
  inline bool is_electron(const McTrack& track) { return track.type() ==  2 || track.type() ==  3; }
  inline bool is_neutrino(const McTrack& track) { return track.type() ==  4; }
  inline bool is_muon    (const McTrack& track) { return track.type() ==  5 || track.type() ==  6; }
  inline bool is_pion    (const McTrack& track) { return track.type() ==  8 || track.type() ==  9; }
  inline bool is_neutron (const McTrack& track) { return track.type() == 13; }
  inline bool is_proton  (const McTrack& track) { return track.type() == 14; }
  inline bool is_tau     (const McTrack& track) { return track.type() == 33 || track.type() == 34; }


  /**
   * 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 j = N; j != 0; --j) {
	if (gRandom->Rndm() * NPE < npe)
	  ++n;
      }
      
      return n;
    }
      
    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 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)
  {
    const int n = npe >> 4;
    
    if (n ==  0) return  1;
    if (n >= 32) return 32;

    return n;
  }


  /**
   * Get intersection points of straight track with cylinder.
   *
   * \param  track   Monte Carlo track
   * \param  can     Monte Carlo detector can
   * \return         up and down stream positions along track
   */
  inline std::pair<double, double> getIntersection(const BaseTrack&        track, 
						   const MONTE_CARLO::can& can) 
  {
    double path[] = { 0.0, 0.0 };

    const Vec3D& pos = track.position();
    const Vec3D& dir = track.direction();

    if (fabs(dir.z) != 0.0) {

      // intersection with bottom or top surface

      const double Z[] = { 
	dir.z > 0 ? can.zmin : can.zmax,
	dir.z > 0 ? can.zmax : can.zmin
      };

      for (int i = 0; i != sizeof(Z)/sizeof(Z[0]); ++i) {

	const double u = (Z[i] - pos.z) / dir.z;
	const double x =  pos.x + u * dir.x;
	const double y =  pos.y + u * dir.y;
      
	if (sqrt(x*x + y*y) <= can.r)
	  path[i] = u;
      }
    }
  
    if (fabs(dir.z) != 1.0) {
    
      // intersection with cylinder wall

      const double a =   (dir.x * dir.x  +  dir.y * dir.y);
      const double b = 2*(dir.x * pos.x  +  dir.y * pos.y);
      const double c =   (pos.x * pos.x  +  pos.y * pos.y)  -  can.r * can.r;
    
      const double q = b*b - 4*a*c;
    
      if (q >= 0) {
      
	const double u[] = {
	  (-b - sqrt(q)) / (2*a),
	  (-b + sqrt(q)) / (2*a)
	};
      
	for (int i = 0; i != sizeof(u) / sizeof(u[0]); ++i) {
	
	  const double z = pos.z + u[i] * dir.z;
	
	  if (z >= can.zmin && z <= can.zmax)
	    path[i] = u[i];
	}
      }
    }

    return std::make_pair(path[0], path[1]);
  }
}

#endif