KVIDLine.h

/***************************************************************************
                          KVIDLine.h  -  description
                             -------------------
    begin                : Nov 10 2004
    copyright            : (C) 2004 by J.D. Frankland
    email                : frankland@ganil.fr
 
$Id: KVIDLine.h,v 1.16 2009/03/03 14:27:15 franklan Exp $
***************************************************************************/
 
#ifndef KVIDLine_H
#define KVIDLine_H
 
#include "TVector2.h"
#include "TMath.h"
#include "TH2.h"
#include "KVIDentifier.h"
class KVIDLine : public KVIDentifier {
 
   class MyVector2 {
      double x, y;
   public:
      MyVector2(double X, double Y) : x(X), y(Y) {}
      MyVector2(const MyVector2& T)
      {
         *this = T;
      }
      MyVector2& operator=(const MyVector2& T)
      {
         if (this != &T) {
            x = T.x;
            y = T.y;
         }
         return *this;
      }
      friend MyVector2 operator-(const MyVector2& a, const MyVector2& b)
      {
         return {a.x - b.x, a.y - b.y};
      }
      friend double operator*(const MyVector2& a, const MyVector2& b)
      {
         return a.x * b.x + a.y * b.y;
      }
      friend MyVector2 operator*(const MyVector2& a, double s)
      {
         return {a.x * s, a.y * s};
      }
      friend MyVector2 operator*(double s, const MyVector2& a)
      {
         return a * s;
      }
      double Mod2() const
      {
         return (*this) * (*this);
      }
      double Mod() const
      {
         return sqrt(Mod2());
      }
      MyVector2 Proj(const MyVector2& v) const
      {
         double scalar = v * (*this);
         scalar /= v.Mod2();
         return v * scalar;
      }
      MyVector2 Norm(const MyVector2& v) const
      {
         return (*this) - Proj(v);
      }
      Double_t Phi() const
      {
         return TMath::Pi() + TMath::ATan2(-y, -x);
      }
   };
 
public:
 
   KVIDLine();
   KVIDLine(const TGraph& gr);
   KVIDLine(const KVIDLine&);
 
   virtual ~ KVIDLine();
 
   void WaitForPrimitive() override;
   void ExecuteEvent(Int_t event, Int_t px, Int_t py) override;
 
   inline Double_t DistanceToLine(Double_t px, Double_t py, Int_t&);
   inline Double_t DistanceToLine(Double_t px, Double_t py, Double_t xp1,
                                  Double_t yp1, Double_t xp2, Double_t yp2,
                                  Int_t&);
   inline Bool_t WhereAmI(Double_t px, Double_t py, Option_t* opt);
   inline Bool_t PosRelToLine(Option_t* opt, Double_t px, Double_t py,
                              Double_t xp1, Double_t yp1, Double_t xp2,
                              Double_t yp2);
 
   inline void GetStartPoint(Double_t& x, Double_t& y) const;
   inline void GetEndPoint(Double_t& x, Double_t& y) const;
   inline Bool_t IsBetweenEndPoints(Double_t x, Double_t y,
                                    const Char_t* axis = "") const;
 
   static KVIDLine* MakeIDLine(TObject* obj, Double_t xdeb = -1., Double_t xfin = -1., Double_t np = 1., Bool_t save = kFALSE);
 
   ClassDefOverride(KVIDLine, 2)       //Base class for lines/cuts used for particle identification
};
 
 
inline Double_t KVIDLine::DistanceToLine(Double_t px, Double_t py,
      Int_t& i_near)
{
 
   Double_t distance, dist2;
   Int_t i;
   Double_t d;
   distance = dist2 = 9999.;
   Int_t i_nearest_point = 0, inear1 = 0, inear2 = 0;
 
   for (i = 0; i < fNpoints - 1; i++) {
      d = DistanceToLine(px, py, fX[i], fY[i], fX[i + 1], fY[i + 1],
                         i_nearest_point);
      if (d >= 0.) {
         if (d < distance) {
            distance = d;
            inear1 = i;
         }
      }
      else {
         if (-d < dist2) {
            dist2 = -d;
            inear2 = i + i_nearest_point;
         }
      }
   }
 
   i_near = inear1;
   if (distance < 9999.)
      return distance;
   i_near = inear2;
   if (inear2 > 0 && inear2 < (fNpoints - 1))
      return dist2;
   return -dist2;
}
 
 
 
inline Double_t KVIDLine::DistanceToLine(Double_t px, Double_t py,
      Double_t xp1, Double_t yp1,
      Double_t xp2, Double_t yp2,
      Int_t& i_nearest_point)
{
 
   MyVector2 P1(xp1, yp1), P2(xp2, yp2), P(px, py);
   MyVector2 P1P2 = P2 - P1;
   MyVector2 P1P = P - P1;
   MyVector2 P2P = P - P2;
   MyVector2 MP = P1P.Norm((P1P2));
   i_nearest_point = 0;
   Double_t sum = (P1P - MP).Mod() + (P2P - MP).Mod();
   if (TMath::Abs(sum - P1P2.Mod()) > 1.e-06) {
      if (P1P.Mod() < P2P.Mod()) {
         i_nearest_point = 0;
         return -(P1P.Mod());
      }
      i_nearest_point = 1;
      return -(P2P.Mod());
   }
   return MP.Mod();
}
 
 
 
inline Bool_t KVIDLine::PosRelToLine(Option_t* opt, Double_t px,
                                     Double_t py, Double_t xp1,
                                     Double_t yp1, Double_t xp2,
                                     Double_t yp2)
{
 
   MyVector2 P1(xp1, yp1), P2(xp2, yp2), P(px, py);
   MyVector2 P1P2 = P2 - P1;
   MyVector2 P1P = P - P1;
   MyVector2 MP = P1P.Norm((P1P2));
   Double_t phi = MP.Phi() * TMath::RadToDeg();
   Bool_t result = kFALSE;
   if (!strcmp(opt, "left")) {
      result = (phi > 90. && phi < 270.);
   }
   else if (!strcmp(opt, "right")) {
      result = (phi < 90. || phi > 270.);
   }
   else if (!strcmp(opt, "above")) {
      result = (phi > 0. && phi < 180.);
   }
   else if (!strcmp(opt, "below")) {
      result = (phi > 180. && phi < 360.);
   }
   return result;
}
 
 
inline Bool_t KVIDLine::WhereAmI(Double_t px, Double_t py, Option_t* opt)
{
 
   Double_t* XX, *YY;
   Double_t  xx,  yy;
   Int_t    sign  = 1;
 
   if (!strcmp(opt, "left")) {
      XX = fY;
      xx = py;
      YY = fX;
      yy = px;
   }
   else if (!strcmp(opt, "right")) {
      XX = fY;
      xx = py;
      YY = fX;
      yy = px;
      sign = -1;
   }
   else if (!strcmp(opt, "above")) {
      XX = fX;
      xx = px;
      YY = fY;
      yy = py;
      sign = -1;
   }
   else if (!strcmp(opt, "below")) {
      XX = fX;
      xx = px;
      YY = fY;
      yy = py;
   }
   else return kFALSE;
 
   Int_t i_start     = 0;
   Int_t i_stop      = fNpoints - 1;
   Int_t prev_i_stop = 0;
   Bool_t same_sign = TMath::Sign(1., XX[i_start] - xx) == TMath::Sign(1., XX[i_stop] - xx);
   while ((i_start < i_stop - 1) || same_sign) {
 
      if (same_sign && (prev_i_stop == 0))  break;
      else if (same_sign) {
         i_start = i_stop;
         i_stop  = prev_i_stop;
      }
      else {
         prev_i_stop = i_stop;
         i_stop = (Int_t)((i_start + i_stop) / 2 + 0.5);
      }
 
      same_sign = TMath::Sign(1., XX[i_start] - xx) == TMath::Sign(1., XX[i_stop] - xx);
   }
 
   Double_t a = (YY[i_stop] - YY[i_start]) / (XX[i_stop] - XX[i_start]);
   Double_t b = YY[i_start] - a * XX[i_start];
 
   Bool_t res = (sign * yy < sign * (a * xx + b));
 
   return res;
}
 
 
inline void KVIDLine::GetStartPoint(Double_t& x, Double_t& y) const
{
#if ROOT_VERSION_CODE >= ROOT_VERSION(4,0,3)
   GetPoint(0, x, y);
#else
   const_cast < KVIDLine* >(this)->GetPoint(0, x, y);
#endif
}
 
 
inline void KVIDLine::GetEndPoint(Double_t& x, Double_t& y) const
{
#if ROOT_VERSION_CODE >= ROOT_VERSION(4,0,3)
   GetPoint((GetN() - 1), x, y);
#else
   const_cast < KVIDLine* >(this)->GetPoint((GetN() - 1), x, y);
#endif
}
 
 
inline Bool_t KVIDLine::IsBetweenEndPoints(Double_t x, Double_t y,
      const Char_t* axis) const
{
 
   TString ax(axis);
   ax.ToUpper();
   Double_t x1, y1, x2, y2;
   GetStartPoint(x1, y1);
   GetEndPoint(x2, y2);
   Double_t xmin = TMath::Min(x1, x2);
   Double_t xmax = TMath::Max(x1, x2);
   Double_t ymin = TMath::Min(y1, y2);
   Double_t ymax = TMath::Max(y1, y2);
 
   Bool_t in_range_x = (x <= xmax && x >= xmin);
   Bool_t in_range_y = (y <= ymax && y >= ymin);
   if (ax == "X") {
      return in_range_x;
   }
   else if (ax == "Y") {
      return in_range_y;
   }
 
   return (in_range_x && in_range_y);
}
 
#endif