CBAlgoCenterOfMass.cxx

Go to the documentation of this file.

#include "larreco/RecoAlg/CMTool/CMTAlgMerge/CBAlgoCenterOfMass.h"

#include "TString.h"

namespace cmtool {

  //----------------------------------------
  CBAlgoCenterOfMass::CBAlgoCenterOfMass() : CBoolAlgoBase()
  //----------------------------------------
  {

    SetDebug(false);
    SetMaxHitsSmallClus(10);
    SetMinHitsBigClus(40);
    SetMaxDistance(20.);
    SetLengthReach(3.0);
    UseCOMInPoly(true);
    UseCOMInCone(true);
    UseCOMNearClus(true);
  }

  //---------------------------------------------------------------------------
  bool CBAlgoCenterOfMass::Bool(const ::cluster::ClusterParamsAlg& cluster1,
                                const ::cluster::ClusterParamsAlg& cluster2)
  //---------------------------------------------------------------------------
  {

    int small = 0;
    int hitssmall = 0;

    //determine if clusters ok and if so which is big and which is small
    if ((cluster1.GetHitVector().size() > _minHits) and (cluster2.GetHitVector().size() < _maxHits))
      small = 2;
    else if ((cluster2.GetHitVector().size() > _minHits) and
             (cluster1.GetHitVector().size() < _maxHits))
      small = 1;
    else
      return false;

    //Define COM values on w & t
    double COM_t_s = 0;
    double COM_w_s = 0;
    double Q_s = 0;
    //Define cone parameters
    double opening_angle;
    double direc_angle;
    double length;
    double start_w;
    double start_t;
    double end_w;
    double end_t;
    Polygon2D poly;

    //Get Hit vector for small cluster
    //std::vector<util::PxHit> hitss;
    std::vector<util::PxHit> hitss;
    if (small == 1) {
      hitss = cluster1.GetHitVector();
      hitssmall = hitss.size();
      direc_angle = cluster2.GetParams().angle_2d;
      opening_angle = cluster2.GetParams().opening_angle * (180. / 3.14);
      length = cluster2.GetParams().length;
      start_w = cluster2.GetParams().start_point.w;
      start_t = cluster2.GetParams().start_point.t;
      end_w = cluster2.GetParams().end_point.w;
      end_t = cluster2.GetParams().end_point.t;
      poly = cluster2.GetParams().PolyObject;
    }
    if (small == 2) {
      hitss = cluster2.GetHitVector();
      hitssmall = hitss.size();
      direc_angle = cluster1.GetParams().angle_2d;
      opening_angle = cluster1.GetParams().opening_angle * (180. / 3.14);
      length = cluster1.GetParams().length;
      start_w = cluster1.GetParams().start_point.w;
      start_t = cluster1.GetParams().start_point.t;
      end_w = cluster1.GetParams().end_point.w;
      end_t = cluster1.GetParams().end_point.t;
      poly = cluster1.GetParams().PolyObject;
    }

    if (_debug) {
      std::cout << "Big Cluster:" << std::endl;
      std::cout << "\tOpening Angle: " << opening_angle << std::endl;
      std::cout << "\tLength: " << length << std::endl;
      std::cout << "\tStart Point: (" << start_w << ", " << end_w << ")" << std::endl;
      std::cout << "\tDirection Angle: " << direc_angle << std::endl;
      std::cout << std::endl;
    }

    for (auto& hit : hitss) {
      COM_t_s += hit.t * hit.charge;
      COM_w_s += hit.w * hit.charge;
      Q_s += hit.charge;
    }
    COM_t_s /= Q_s;
    COM_w_s /= Q_s;

    if (_debug) {
      std::cout << "Small Cluster: " << std::endl;
      std::cout << "N Hits: " << hitssmall << std::endl;
      std::cout << "COM: (w,t) -> (" << COM_w_s << ", " << COM_t_s << ")" << std::endl;
      std::cout << std::endl;
    }

    //Get COM
    std::pair<float, float> COM_s;
    COM_s = std::make_pair(COM_w_s, COM_t_s);

    //Get rotate COM to see if in cone
    double COM_w_rot = (COM_w_s - start_w) * cos(direc_angle * 3.14 / 180.) +
                       (COM_t_s - start_t) * sin(direc_angle * 3.14 / 180.);
    double COM_t_rot = -(COM_w_s - start_w) * sin(direc_angle * 3.14 / 180.) +
                       (COM_t_s - start_t) * cos(direc_angle * 3.14 / 180.);

    //look for polygon overlap
    if ((poly.PointInside(COM_s)) and _COMinPolyAlg) {
      if (_verbose) { std::cout << "Polygon Overlap -> Merge!" << std::endl << std::endl; }
      return true;
    }
    //look for COM in cone
    if (_COMinConeAlg and (COM_w_rot < length * _lengthReach) and (COM_w_rot > 0) and
        (abs(COM_t_rot) < abs(COM_w_rot * sin(opening_angle * 3.14 / 180.)))) {
      if (_verbose) { std::cout << "COM in Cone -> Merge! " << std::endl; }
      return true;
    }
    //look for COM close to start-end of other cluster
    if (_COMNearClus and
        ShortestDistanceSquared(COM_w_s, COM_t_s, start_w, start_t, end_w, end_t) < _MaxDist) {
      if (_verbose) { std::cout << "COM close to start-end -> Merge!" << std::endl; }
      return true;
    }

    return false;
  }

  //-----------------------
  void CBAlgoCenterOfMass::Report()
  //-----------------------
  {}

  double CBAlgoCenterOfMass::ShortestDistanceSquared(double point_x,
                                                     double point_y,
                                                     double start_x,
                                                     double start_y,
                                                     double end_x,
                                                     double end_y) const
  {

    //This code finds the shortest distance between a point and a line segment.
    //code based off sample from
    //http://stackoverflow.com/questions/849211/shortest-distance-between-a-point-and-a-line-segment
    //note to self: rewrite this with TVector2 and compare time differences...
    //TVector2 code might be more understandable

    double distance_squared = -1;

    // Line segment: from ("V") = (start_x, start_y) to ("W")=(end_x, end_y)
    double length_squared = pow((end_x - start_x), 2) + pow((end_y - start_y), 2);

    // Treat the case start & end point overlaps
    if (length_squared < 0.1) {
      if (_verbose) {
        std::cout << std::endl;
        std::cout << Form(" Provided very short line segment: (%g,%g) => (%g,%g)",
                          start_x,
                          start_y,
                          end_x,
                          end_y)
                  << std::endl;
        std::cout << " Likely this means one of two clusters have start & end point identical."
                  << std::endl;
        std::cout << " Check the cluster output!" << std::endl;
        std::cout << std::endl;
        std::cout << Form(" At this time, the algorithm uses a point (%g,%g)", start_x, start_y)
                  << std::endl;
        std::cout << " to represent this cluster's location." << std::endl;
        std::cout << std::endl;
      }

      return (pow((point_x - start_x), 2) + pow((point_y - start_y), 2));
    }

    //Find shortest distance between point ("P")=(point_x,point_y) to this line segment
    double t = ((point_x - start_x) * (end_x - start_x) + (point_y - start_y) * (end_y - start_y)) /
               length_squared;

    if (t < 0.0)
      distance_squared = pow((point_x - start_x), 2) + pow((point_y - start_y), 2);

    else if (t > 1.0)
      distance_squared = pow((point_x - end_x), 2) + pow(point_y - end_y, 2);

    else
      distance_squared = pow((point_x - (start_x + t * (end_x - start_x))), 2) +
                         pow((point_y - (start_y + t * (end_y - start_y))), 2);

    return distance_squared;

  } //end ShortestDistanceSquared function

}