#ifndef __JRECONSTRUCTION__JSTART__
#define __JRECONSTRUCTION__JSTART__

#include <istream>
#include <ostream>


/**
 * \file
 * Auxiliary method to locate start and end point of muon trajectory.
 * \author mdejong
 */
namespace JRECONSTRUCTION {}
namespace JPP { using namespace JRECONSTRUCTION; }

namespace JRECONSTRUCTION {

  /**
   * Auxiliary class for start or end point evaluation.
   */
  struct JStart {
    /**
     * Default constructor.
     */
    JStart() :
      Pmin1(0.0),
      Pmin2(0.0),
      Nmax2(0)
    {}

    
    /**
     * Constructor.
     *
     * \param  Pmin1              minimal probability single observation
     * \param  Pmin2              minimal probability for twofold observations
     * \param  Nmax2              maximal number for twofold observations
     */
    JStart(const double Pmin1,
	   const double Pmin2,
	   const int    Nmax2) :
      Pmin1(Pmin1),
      Pmin2(Pmin2),
      Nmax2(Nmax2)
    {}

    
    /**
     * Check validity of start or end point conditions.
     *
     * \return                     true if valid; else false
     */
    bool is_valid() const
    {
      return (Pmin1 >= 0.0 && Pmin1 <= 1.0 &&
	      Pmin2 >= 0.0 && Pmin2 <= 1.0);
    }

  
    /**
     * Get start point of muon trajectory.
     * 
     * The template parameter should correspond to a data type which has the member method: 
     * <pre>
     *     double getP() const;    // return probability
     * </pre>
     * The input data should have been sorted along the muon path beforehand.\n
     * A start point is triggered when an observed probability is less than Pmin1 or
     * when there are two probabilities less than Pmin2 within Nmax2 consecutive observations.\n
     * The start or end point can found by providing the data in ascending or descending order, respectively.
     *
     * \param  __begin            begin of data
     * \param  __end              end   of data
     * \return                    start point
     */
    template<class T>
    inline T find(T __begin, T __end) const
    {
      for (T p = __begin; p != __end; ++p) {

	if (p->getP() < this->Pmin1) {
	  return p;
	}

	if (p->getP() < this->Pmin2) {
	  for (T q = p; ++q != __end && distance(p,q) < Nmax2; ) {
	    if (q->getP() < this->Pmin2) {
	      return p;
	    }
	  }
	}
      }
      
      return __end;
    }

    
    /**
     * Read parameters from input.
     *
     * \param  in                  input stream
     * \param  parameters          parameters
     * \return                     input stream
     */
    friend inline std::istream& operator>>(std::istream& in, JStart& parameters)
    {
      in >> parameters.Pmin1 >> parameters.Pmin2 >> parameters.Nmax2;

      return in;
    }


    /**
     * Write parameters to output.
     *
     * \param  out                 output stream
     * \param  parameters          parameters
     * \return                     output stream
     */
    friend inline std::ostream& operator<<(std::ostream& out, const JStart& parameters)
    {
      out << parameters.Pmin1 << ' ' << parameters.Pmin2 << ' ' << parameters.Nmax2;

      return out;
    }
    

    double Pmin1;              //!< minimal probability single observation
    double Pmin2;              //!< minimal probability for twofold observations
    int    Nmax2;              //!< maximal number for twofold observations
  };
}

#endif