Cluster.cxx

Go to the documentation of this file.

#include "DDTMonopole/Cluster.h"

#include "SlowMonopoleTrigger/Constants.h"

#include <iostream>



mono::Cluster::Cluster(fhicl::ParameterSet const& p) :
  max_cell_gap_cut_ (p.get<int>("max_cell_gap_cut")),
  max_plane_gap_cut_(p.get<int>("max_plane_gap_cut")),
  track_min_beta_   (p.get<double>("track_min_beta")),
  track_max_beta_   (p.get<double>("track_max_beta"))
{
}



rb::Cluster mono::Cluster::get() const
{
  return monopole_hits_;
}



bool mono::Cluster::make(mono::HitList const& hits)
{
  monopole_hits_.Clear();

  bool result = false;

  std::map<std::string, mono::HitList> hitmap;
  split_by_view(hits, hitmap);

  std::cout << "MF: (Total Hits, X Hits, Y Hits) = ("
            << hits.size()
            << ", " << hitmap["x"].size()
            << ", " << hitmap["y"].size()
            << ")" << std::endl;

  if (hitmap["x"].size() < novaddt::smt::Constants::TRACK_MIN_HITS ||
      hitmap["y"].size() < novaddt::smt::Constants::TRACK_MIN_HITS)
    return false;

  hitmap["x_surface"];
  hitmap["x_contained"];
  
  find_view_matched_hits(hitmap);

  if(hitmap["x_contained"].size() < novaddt::smt::Constants::TRACK_MIN_HITS)
    return false;

  for (auto entry_hit = hitmap["x_surface"].cbegin();
       entry_hit     != hitmap["x_surface"].cend() && !result; 
       ++entry_hit)
  {
    // we add 1, so that we avoid combinations between the same hit
    mono::HitList::const_reverse_iterator until(entry_hit + 1);

    auto contained_first = find(hitmap["x_contained"], (*entry_hit)->TDC());
    
    std::find_if(hitmap["x_surface"].crbegin(), until, 
                 [&](art::Ptr<rb::CellHit> const& exit_hit)
    {
      auto contained_last =
        find(contained_first, hitmap["x_contained"].cend(), exit_hit->TDC());

      if (std::distance(contained_first, contained_last) <
          novaddt::smt::Constants::TRACK_MIN_HITS)
        return false;
  
      Track track(*entry_hit, exit_hit);

      if (!good_track(track))
        return false;

      std::set<int> plane_set = { (*entry_hit)->Plane() };
      std::set<int> cell_set =  { (*entry_hit)->Cell() };

      int32_t t_max_plane_gap = track.start_rb()->TDC();
      if (track.slope_plane_is_valid())
        t_max_plane_gap += abs(track.slope_time_plane()) * max_plane_gap_cut_;

      int32_t t_max_cell_gap  = track.start_rb()->TDC();
      if (track.slope_cell_is_valid())
        t_max_cell_gap += abs(track.slope_time_cell()) * max_cell_gap_cut_;
      int32_t t_max_gap = std::min(t_max_plane_gap, t_max_cell_gap);
      
      unsigned n_contained_hit = 0;
      bool contains_gap_exceeding_cut = std::any_of
        (contained_first, contained_last,
         [&](art::Ptr<rb::CellHit> const& hit)
         {
           if (hit->TDC() > t_max_gap)
           {
             int64_t time_difference =
               static_cast<int64_t>(hit->TDC()) -
               static_cast<int64_t>(track.start_rb()->TDC());

             if (hit->TDC() > t_max_cell_gap)
             {
               double expected_cell =
               static_cast<double>(track.start_rb()->Cell());
               if (track.slope_time_is_valid())
                 expected_cell += track.slope_cell_time() * time_difference;
               
               if (gap_exceeds_cut(cell_set, expected_cell, max_cell_gap_cut_))
                 return true;
             }

             if (hit->TDC() > t_max_plane_gap)
             {
               double expected_plane =
                 static_cast<double>(track.start_rb()->Plane());
               if (track.slope_time_is_valid())
                 expected_plane += track.slope_plane_time() * time_difference;
               
               if (gap_exceeds_cut
                   (plane_set, expected_plane, max_plane_gap_cut_))
                 return true;
             }
           }

           if (hit_is_in_time_with_road(hit, track))
           {
             cell_set.insert(hit->Cell());
             plane_set.insert(hit->Plane());
           }   

           return false;
         });

      if (contains_gap_exceeding_cut)
        return false;

      plane_set.insert(exit_hit->Plane());
      if (plane_set.size() < novaddt::smt::Constants::TRACK_MIN_UNIQUE_PLANES)
        return false;
      
      if (find_max_gap(plane_set) > max_plane_gap_cut_)
        return false;

      cell_set.insert(exit_hit->Cell());
      if (find_max_gap(cell_set) > max_cell_gap_cut_)
        return false;

      result = true;
      monopole_hits_.Add(*entry_hit);
      monopole_hits_.Add(exit_hit);
      std::for_each
        (contained_first, contained_last,
         [&](art::Ptr<rb::CellHit> const& hit)
         {
           if (hit_is_in_time_with_road(hit, track))
             monopole_hits_.Add(hit);
         });

      return true;
    });
  }

  return result;
}



void mono::Cluster::split_by_view
(mono::HitList hits, std::map<std::string, mono::HitList> & hitmap) const
{
  auto bound = std::stable_partition
    (hits.begin(), hits.end(),
     [](art::Ptr<rb::CellHit> const& h)
     { return h->View() == geo::View_t::kX; });

  hitmap["x"] = mono::HitList(std::distance(hits.begin(), bound));
  std::move(hits.begin(), bound, hitmap["x"].begin());

  hitmap["y"] = mono::HitList(std::distance(bound, hits.end()));
  std::move(bound, hits.end(), hitmap["y"].begin());
}



void mono::Cluster::find_view_matched_hits
(std::map<std::string, mono::HitList> & hitmap) const
{
  for (auto const& x_hit : hitmap["x"])
  {
    auto y_first =
      find(hitmap["y"], x_hit->TDC() -
           novaddt::smt::Constants::VIEW_MATCH_TDC_CUT);
    auto y_last  =
      find(hitmap["y"], x_hit->TDC() +
           novaddt::smt::Constants::VIEW_MATCH_TDC_CUT);

    auto y_match = std::find_if
      (y_first, y_last, [&](art::Ptr<rb::CellHit> const& y_hit)
      { return hits_are_view_matched(x_hit, y_hit); });

    if (y_match != y_last)
    {
      // we found a matched hit
      hitmap["x_contained"].push_back(x_hit);

      // is this hit on the surface?
      if (hit_is_on_surface(x_hit))
      {
        hitmap["x_surface"].push_back(x_hit);
      } else {
        // now we need to look for additional matched y-hits and check whether
        // they are on the surface
        auto surface_y_match = std::find_if
          (y_match, y_last, [&](art::Ptr<rb::CellHit> const& y_hit)
          {
            if (!hits_are_view_matched(x_hit, y_hit))
              return false;

            if (!hit_is_on_surface(y_hit))
              return false;

            return true;
          });

        if (surface_y_match != y_last)
          hitmap["x_surface"].push_back(x_hit);
      }
    }
  }
}



bool mono::Cluster::hits_are_view_matched
(art::Ptr<rb::CellHit> const& x_hit, art::Ptr<rb::CellHit> const& y_hit) const
{
  int plane_difference =
    static_cast<int>(x_hit->Plane()) - static_cast<int>(y_hit->Plane());
  
  if (std::abs(plane_difference) <
      novaddt::smt::Constants::VIEW_MATCH_PLANE_CUT)
    return true;

  return false;
}



bool mono::Cluster::hit_is_on_surface
(art::Ptr<rb::CellHit> const& hit) const
{
  if (hit->Cell() < novaddt::smt::Constants::SURFACE_MIN_CELL)
    return true;

  if (hit->Cell() > novaddt::smt::Constants::SURFACE_MAX_CELL)
    return true;

  if (hit->Plane() < novaddt::smt::Constants::SURFACE_MIN_PLANE)
    return true;

  if (hit->Plane() > novaddt::smt::Constants::SURFACE_MAX_PLANE)
    return true;

  return false;
}



bool mono::Cluster::hit_is_between_track_end_points
(art::Ptr<rb::CellHit> const& hit, Track const& track) const
{
  if (hit->Plane() < track.min_plane())
    return false;

  if (hit->Plane() > track.max_plane())
    return false;

  if (hit->Cell() < track.min_cell())
    return false;

  if (hit->Cell() > track.max_cell())
    return false;

  return true;
}



bool mono::Cluster::good_track(Track const& track) const
{
  if (track.distance() < novaddt::smt::Constants::TRACK_MIN_DISTANCE)
    return false;

  if (track.beta() < track_min_beta_)
    return false;

  if (track.beta() > track_max_beta_)
    return false;

  if (abs(track.plane_difference()) <= max_plane_gap_cut_)
    return false;
  
  return true; 
}



bool mono::Cluster::hit_is_on_road
(art::Ptr<rb::CellHit> const& hit, Track const& track) const
{
  if (!hit_is_between_track_end_points(hit, track))
    return false;

  // If the plane of the hit and track end points are all equal,
  // the hit is on the road.
  if (!track.slope_plane_is_valid())
    return true;
  
  double road_center_cell = track.slope_cell_plane() *
    (static_cast<double>(hit->Plane()) -
     static_cast<double>(track.start_rb()->Plane())) +
    static_cast<double>(track.start_rb()->Cell());

  if (hit->Cell() > road_center_cell + novaddt::smt::Constants::ROAD_HALF_WIDTH)
    return false;

  if (hit->Cell() < road_center_cell - novaddt::smt::Constants::ROAD_HALF_WIDTH)
    return false;

  return true;
}



bool mono::Cluster::hit_is_in_time_with_road
(art::Ptr<rb::CellHit> const& hit, Track const& track) const
{
  if (!hit_is_on_road(hit, track))
    return false;

  // If the plane of the hit and track end points are all equal,
  // we cannot calculate the time_slope and make no statement on the time.
  if (!track.slope_plane_is_valid())
    return false;


  int32_t relative_hit_time = hit->TDC() - track.start_rb()->TDC();
  
  double expected_time = track.slope_time_plane() *
    (static_cast<double>(hit->Plane()) -
     static_cast<double>(track.start_rb()->Plane()));

  if (relative_hit_time >
      expected_time + novaddt::smt::Constants::ROAD_TIME_VARIATION)
    return false;

  if (relative_hit_time <
      expected_time - novaddt::smt::Constants::ROAD_TIME_VARIATION)
    return false;

  return true;
}



bool mono::Cluster::gap_exceeds_cut
(std::set<int> collection, int const& object, int const& max_gap_cut) const
{
  // Note that we copy the collection in here on purpose, so that we do not
  // accidentally alter it.
  collection.insert(object);

  // We need to encode the max_plane cut differently here, so as to keep
  // this method as general as possible.
  if (find_max_gap(collection) > max_gap_cut)
    return true;

  return false;
}



int mono::Cluster::find_max_gap
(std::set<int> const& collection) const
{
  int result = -9999;

  std::adjacent_find
  (collection.begin(), collection.end(), [&](int object, int next_object)
  {
    int gap = next_object - object;
    result = std::max(gap, result);
    
    // This return is meaningless, but a boolean return is required by
    // adjacent find.  We choose false, so the loop goes until the end.
    return false;
  });

  return result;
}



mono::HitList::const_iterator mono::Cluster::find
(mono::HitList::const_iterator begin,
 mono::HitList::const_iterator end, int32_t tdc) const
{
  return std::lower_bound
    (begin, end, tdc,
     [](art::Ptr<rb::CellHit> const& h, int32_t const& tdc_val)
     { return h->TDC() < tdc_val; });
}



mono::HitList::const_iterator mono::Cluster::find
(mono::HitList const& hits, int32_t tdc) const
{
  return find(hits.cbegin(), hits.cend(), tdc);
}