#ifndef __JMATCH3B__
#define __JMATCH3B__

#include "JTools/JConstants.hh"

#include "JTrigger/JMatch.hh"


namespace JTRIGGER {

  namespace {
    using JTOOLS::getInverseSpeedOfLight;
    using JTOOLS::getIndexOfRefraction;
    using JTOOLS::getTanThetaC;
    using JTOOLS::getSinThetaC;
  }


  /** 
   * 3D match criterion with road width.
   * This match algorithm is intented for muon signals.
   *
   * B. Bakker, "Trigger studies for the Antares and KM3NeT detector.",
   * Master thesis, University of Amsterdam.
   */
  template<class JHit_t>
  class JMatch3B  :
    public JMatch<JHit_t>
  {
  public:
    /**
     * Constructor.
     *
     * \param  road_width_m     maximal road width [m]
     * \param  Tmax_ns          maximal extra time [ns]
     */
    JMatch3B(const double road_width_m,
	     const double Tmax_ns = 0.0) :
      roadWidth_m (road_width_m),
      TMaxExtra_ns(Tmax_ns)
    {
      //
      // calculation D2 in thesis is wrong, here correct
      //
      const double tt2 = getTanThetaC() * getTanThetaC();

      D0 = roadWidth_m;
      D1 = roadWidth_m * 2.0;
      D2 = roadWidth_m * 0.5 * sqrt(tt2 + 10.0 + 9.0/tt2);

      D02 = D0 * D0;
      D12 = D1 * D1;
      D22 = D2 * D2;

      const double R  = roadWidth_m;
      const double Rs = R * getSinThetaC();

      R2  = R  * R;
      Rs2 = Rs * Rs;
      Rst = Rs * getTanThetaC();
      Rt  = R  * getTanThetaC();
    }


    /**
     * Clone object.
     *
     * \return          match result
     */
    JMatch<JHit_t>* clone() const
    {
      return new JMatch3B<JHit_t>(*this);
    }


    /**
     * Match operator.
     *
     * \param  first    hit
     * \param  second   hit
     * \return          match result
     */
    virtual bool operator()(const JHit_t& first, const JHit_t& second) const
    {
      //
      x  = first.getX() - second.getX();
      y  = first.getY() - second.getY();
      z  = first.getZ() - second.getZ();
      d2 = x*x + y*y + z*z;
      t  = fabs(first.getT() - second.getT());

      if (d2 < D02)
	dmax = sqrt(d2) * getIndexOfRefraction();
      else 
	dmax = sqrt(d2 - Rs2) + Rst;

      if (t > dmax * getInverseSpeedOfLight()  +  TMaxExtra_ns)
	return false;

      if      (d2 > D22)
	dmin = sqrt(d2 - R2) - Rt;
      else if (d2 > D12)
	dmin = sqrt(d2 - D12);
      else
	return true;

      return t >= dmin * getInverseSpeedOfLight()  -  TMaxExtra_ns;
      //
      /*
      x = first.getX() - second.getX();
      y = first.getY() - second.getY();
      z = first.getZ() - second.getZ();
      d = sqrt(x*x + y*y + z*z);
      t = fabs(first.getT() - second.getT());

      if (d < D0)
	dmax = d * getIndexOfRefraction();
      else 
	dmax = sqrt((d + roadWidth_m*getSinThetaC()) * (d - roadWidth_m*getSinThetaC()))  +  roadWidth_m*getSinThetaC() * getTanThetaC();

      if (t > dmax * getInverseSpeedOfLight()  +  TMaxExtra_ns)
	return false;

      if      (d > D2)
	dmin = sqrt((d + roadWidth_m) * (d - roadWidth_m))  -  roadWidth_m * getTanThetaC();
      else if (d > D1) 
	dmin = sqrt((d + D1) * (d - D1));
      else
        return true;

      return t >= dmin * getInverseSpeedOfLight()  -  TMaxExtra_ns;
      */
    }


    double roadWidth_m;
    double TMaxExtra_ns;

  private:
    double D0;
    double D1;
    double D2;
    mutable double x;
    mutable double y;
    mutable double z;
    mutable double d;
    mutable double t;
    mutable double dmin;
    mutable double dmax;

    double D02;
    double D12;
    double D22;
    double Rs2;
    double Rst;
    double Rt;
    double R2;
    mutable double d2;
  };
}

#endif