CBAlgoCenterOfMassSmall.cxx

Go to the documentation of this file.

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

#include "TString.h"

namespace cmtool {

  //----------------------------------------
  CBAlgoCenterOfMassSmall::CBAlgoCenterOfMassSmall() : CBoolAlgoBase()
  //----------------------------------------
  {

    SetDebug(false);
    SetMaxHitsSmallClus(30);
    SetMaxDistance(20.);
    SetMaxCOMDistance(25.);
    UseCOMInPoly(true);
    UseCOMClose(true);
    UseCOMNearClus(true);
  }

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

    int Nhits1 = 0;
    int Nhits2 = 0;

    //Both clusters should have less hits than some threshold
    if ((cluster1.GetHitVector().size() > _maxHits) or (cluster2.GetHitVector().size() > _maxHits))
      return false;

    //Define COM values on w & t
    double COM_t_1 = 0;
    double COM_w_1 = 0;
    double Q_1 = 0;
    double start_w_1;
    double start_t_1;
    double end_w_1;
    double end_t_1;
    Polygon2D poly1;
    double COM_t_2 = 0;
    double COM_w_2 = 0;
    double Q_2 = 0;
    double start_w_2;
    double start_t_2;
    double end_w_2;
    double end_t_2;
    Polygon2D poly2;

    //Get Hit vector for cluster 1
    //std::vector<util::PxHit> hits1;
    std::vector<util::PxHit> hits1;
    hits1 = cluster1.GetHitVector();
    Nhits1 = hits1.size();
    poly1 = cluster1.GetParams().PolyObject;
    start_w_1 = cluster1.GetParams().start_point.w;
    start_t_1 = cluster1.GetParams().start_point.t;
    end_w_1 = cluster1.GetParams().end_point.w;
    end_t_1 = cluster1.GetParams().end_point.t;
    //Get Hit vector for cluster 2
    //std::vector<util::PxHit> hits2;
    std::vector<util::PxHit> hits2;
    hits2 = cluster2.GetHitVector();
    Nhits2 = hits2.size();
    poly2 = cluster2.GetParams().PolyObject;
    start_w_2 = cluster2.GetParams().start_point.w;
    start_t_2 = cluster2.GetParams().start_point.t;
    end_w_2 = cluster2.GetParams().end_point.w;
    end_t_2 = cluster2.GetParams().end_point.t;

    //Find COM for cluster 1
    for (auto& hit : hits1) {
      COM_t_1 += hit.t * hit.charge;
      COM_w_1 += hit.w * hit.charge;
      Q_1 += hit.charge;
    }
    COM_t_1 /= Q_1;
    COM_w_1 /= Q_1;

    //Find COM for cluster 2
    for (auto& hit : hits2) {
      COM_t_2 += hit.t * hit.charge;
      COM_w_2 += hit.w * hit.charge;
      Q_2 += hit.charge;
    }
    COM_t_2 /= Q_2;
    COM_w_2 /= Q_2;

    if (_debug) {
      std::cout << "Cluster 1: " << std::endl;
      std::cout << "N Hits: " << Nhits1 << std::endl;
      std::cout << "COM: (w,t) -> (" << COM_w_1 << ", " << COM_t_1 << ")" << std::endl;
      std::cout << "Cluster 2: " << std::endl;
      std::cout << "N Hits: " << Nhits2 << std::endl;
      std::cout << "COM: (w,t) -> (" << COM_w_2 << ", " << COM_t_2 << ")" << std::endl;
      std::cout << std::endl;
    }

    //Get COM
    std::pair<float, float> COM_1;
    COM_1 = std::make_pair(COM_w_1, COM_t_1);
    std::pair<float, float> COM_2;
    COM_2 = std::make_pair(COM_w_2, COM_t_2);

    //look for polygon overlap
    if (((poly2.PointInside(COM_1)) and _COMinPolyAlg) or
        ((poly1.PointInside(COM_2)) and _COMinPolyAlg)) {
      if (_verbose) { std::cout << "Polygon Overlap -> Merge!" << std::endl << std::endl; }
      return true;
    }

    //look for COM of 1 close to COM of 2
    double distCOMs =
      (COM_w_1 - COM_w_2) * (COM_w_1 - COM_w_2) + (COM_t_1 - COM_t_2) * (COM_t_1 - COM_t_2);
    if (_COMsClose and (distCOMs < _MaxCOMDistSquared)) {
      if (_verbose) { std::cout << "COMs close to each other -> Merge!" << std::endl << std::endl; }
      return true;
    }

    //look for COM close to start-end of other cluster
    if (_COMNearClus and
        ((ShortestDistanceSquared(COM_w_1, COM_t_1, start_w_2, start_t_2, end_w_2, end_t_2) <
          _MaxDist) or
         (ShortestDistanceSquared(COM_w_2, COM_t_2, start_w_1, start_t_1, end_w_1, end_t_1) <
          _MaxDist))) {
      if (_verbose) { std::cout << "COM close to start-end -> Merge!" << std::endl; }
      return true;
    }

    return false;
  }

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

  double CBAlgoCenterOfMassSmall::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

}