#ifndef __JDETECTOR__JDETECTORSUBSET__
#define __JDETECTOR__JDETECTORSUBSET__

#include <algorithm>
#include <limits>

#include "JGeometry3D/JPosition3D.hh"
#include "JGeometry3D/JDirection3D.hh"
#include "JGeometry3D/JAxis3D.hh"
#include "JGeometry3D/JTransformation3D.hh"
#include "JTools/JRange.hh"
#include "JDetector/JDetector.hh"
#include "JDetector/JModuleComparator.hh"


/**
 * \file
 *
 * Auxiliary class to extract a subset of optical modules from a detector.
 * \author mdejong
 */
namespace JDETECTOR {}
namespace JPP { using namespace JDETECTOR; }

namespace JDETECTOR {

  using JGEOMETRY3D::JPosition3D;
  using JGEOMETRY3D::JVector3D;
  using JGEOMETRY3D::JDirection3D;
  using JGEOMETRY3D::JAxis3D;
  using JGEOMETRY3D::JTransformation3D;
  using JTOOLS::JRange;
  

  /**
   * Detector subset without binary search functionality.
   */
  class JDetectorSubset_t :
    public std::vector<JModule>
  {
  public:
    /**
     * Constructor.
     *
     * The positions and orientations of the modules and the PMTs after the transformation
     * are relative to given particle track (see ANTARES-SOFT 2010-002).
     * Only modules within given road width and z-axis range are accepted.
     *
     * \param  detector       detector
     * \param  track          particle track
     * \param  Rmax           maximal distance of approach [m]
     * \param  Z              range of positions along z-axis [m]
     */
    JDetectorSubset_t(const JDetector&      detector,
		      const JAxis3D&        track,
		      const double          Rmax = std::numeric_limits<double>::max(),
		      const JRange<double>& Z    = JRange<double>()) 
    {
      const JTransformation3D transformation(track.getPosition(), track.getDirection());

      JModule module;

      for (JDetector::const_iterator i = detector.begin(); i != detector.end(); ++i) {

	JPosition3D pos(i->getPosition());

	pos.transform(transformation.getRotation(), transformation.getPosition());

	if (Z(pos.getZ()) && pos.getX() <= Rmax) {
	      
	  module = *i;
	    
	  module.transform(transformation);
	  
	  this->push_back(module);
	}
      }
    }


    /**
     * Constructor.
     *
     * The positions and orientations of the modules and the PMTs
     * are offset with respect to given position.
     * Only modules within given distance are accepted.
     *
     * \param  detector       detector
     * \param  position       position of vertex
     * \param  Dmax           maximal distance [m]
     */
    JDetectorSubset_t(const JDetector& detector,
		      const JVector3D& position,
		      const double     Dmax = std::numeric_limits<double>::max()) 
    {
      for (JDetector::const_iterator i = detector.begin(); i != detector.end(); ++i) {

	JPosition3D pos(i->getPosition());

	pos.sub(position);

	if (pos.getLength() <= Dmax) {
	  
	  JModule module(*i);

	  module.sub(position);
	  
	  this->push_back(module);
	}
      }
    }
  };


  /**
   * Detector subset with binary search functionality.
   */
  template<class JComparator_t>
  class JDetectorSubset :
    public JDetectorSubset_t
  {
  public:
    /**
     * Auxiliary class for range of iterators.
     */
    struct range_type {
      /**
       * Constructor.
       *
       * \param begin           begin of iteration
       * \param end             end   of iteration
       */
      range_type(const_iterator begin,
		 const_iterator end) :
	__begin(begin),
	__end  (end)
      {}

      inline const_iterator begin() const { return __begin; }   //!< begin of iteration
      inline const_iterator end()   const { return __end; }     //!< end   of iteration

      bool empty() const { return __begin == __end; }           //!< check emptyness
      
    protected:
      const_iterator __begin;
      const_iterator __end;
    };

    
    /**
     * Constructor.
     * 
     * See constructor JDetectorSubset_t::JDetectorSubset_t.
     * The modules are sorted according the specified comparator.
     *
     * \param  detector       detector
     * \param  track          particle track
     * \param  comparator     comparator
     * \param  Rmax           maximal distance of approach [m]
     * \param  Z              range of positions along z-axis [m]
     */
    JDetectorSubset(const JDetector&      detector,
		    const JAxis3D&        track,
		    const JComparator_t&  comparator,
		    const double          Rmax = std::numeric_limits<double>::max(),
		    const JRange<double>& Z    = JRange<double>()) :
      JDetectorSubset_t(detector, track, Rmax, Z),
      compare(comparator)
    {
      std::sort(this->begin(), this->end(), compare);
    }

    
    /**
     * Constructor.
     *
     * See constructor JDetectorSubset_t::JDetectorSubset_t.
     * The modules are sorted according the default comparator.
     *
     * \param  detector       detector
     * \param  track          particle track
     * \param  Rmax           maximal distance of approach [m]
     * \param  Z              range of positions along z-axis [m]
     */
    JDetectorSubset(const JDetector&      detector,
		    const JAxis3D&        track,
		    const double          Rmax = std::numeric_limits<double>::max(), 
		    const JRange<double>& Z    = JRange<double>()) :
      JDetectorSubset_t(detector, track, Rmax, Z),
      compare()
    {
      std::sort(this->begin(), this->end(), compare);
    }


    /**
     * Constructor.
     *
     * See constructor JDetectorSubset_t::JDetectorSubset_t.
     * The modules are sorted according the specified comparator.
     *
     * \param  detector       detector
     * \param  position       position of vertex
     * \param  comparator     comparator
     * \param  Dmax           maximal distance [m]
     */
    JDetectorSubset(const JDetector&     detector,
		    const JVector3D&     position,
		    const JComparator_t& comparator,
		    const double         Dmax = std::numeric_limits<double>::max()) :
      JDetectorSubset_t(detector, position, Dmax),
      compare(comparator)
    {
      std::sort(this->begin(), this->end(), compare);      
    }


    /**
     * Constructor.
     *
     * See constructor JDetectorSubset_t::JDetectorSubset_t.
     * The modules are sorted according the default comparator.
     *
     * \param  detector       detector
     * \param  position       position of vertex
     * \param  Dmax           maximal distance [m]
     */
    JDetectorSubset(const JDetector&     detector,
		    const JVector3D&     position,
		    const double         Dmax = std::numeric_limits<double>::max()) :
      JDetectorSubset_t(detector, position, Dmax),
      compare()
    {
      std::sort(this->begin(), this->end(), compare);      
    }


    /**
     * Get comparator.
     *
     * \return                compartor
     */
    const JComparator_t& getComparator() const
    {
      return compare;
    }


    /**
     * Find first module after given value according specified comparator.
     *
     * \param  value          value
     * \return                iterator to module
     */
    const_iterator lower_bound(const double value) const
    {
      return std::lower_bound(this->begin(), this->end(), value, compare);
    }


    /**
     * Get range of modules between given values according specified comparator.
     *
     * \param  xmin           minimal value
     * \param  xmax           maximal value
     * \return                range of iterators
     */
    range_type getRange(const double xmin,
			const double xmax) const
    {
      return range_type(this->lower_bound(xmin),
			this->lower_bound(xmax));
    }
    
  protected:
    const JComparator_t compare;
  };
}

#endif