Trigger.cxx

Go to the documentation of this file.

#include "SlowMonopoleTrigger/Trigger.h"

#include "DDTBaseDataProducts/CompareDAQHit.h"
#include "SlowMonopoleTrigger/Constants.h"
#include "SlowMonopoleTrigger/Track.h"

#include <algorithm>
#include <bitset>
#include <chrono>
#include <iostream>

using namespace std::chrono;



novaddt::smt::Trigger::Trigger(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")),
  sparsification_factor_         (p.get<unsigned>("sparsification_factor")),
  track_max_beta_                (p.get<double>("track_max_beta")),
  track_min_beta_                (p.get<double>("track_min_beta")),
  trigger_time_buffer_           (p.get<unsigned>("trigger_time_buffer")),
  trigger_timeout_               (p.get<double>("trigger_timeout")),
  trigger_use_timeout_           (p.get<bool>("trigger_use_timeout")),
  trigger_write_all_timeout_data_(p.get<bool>("trigger_write_all_timeout_data"))
{
}



std::vector<novaddt::TriggerDecision> 
novaddt::smt::Trigger::trigger_decisions() const
{
  return trigger_decisions_;
}



bool novaddt::smt::Trigger::run_algorithm
(art::Handle<novaddt::HitList> const& hits)
{
  trigger_decisions_.clear();
  std::bitset<Extra_Info::BITSET_SIZE> extra_info;
  extra_info.set(Extra_Info::CHECK_BIT, true);

  auto t_start = high_resolution_clock::now();
  bool result = false;
  novaddt::DAQHit trigger_entry_hit, trigger_exit_hit;

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

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

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

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

  unsigned n_good_track = 0;
  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
    novaddt::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, 
    [&](novaddt::DAQHit const& exit_hit)
    {
      auto contained_last =
        find(contained_first, hitmap["x_contained"].cend(), exit_hit.TDC());

      if (std::distance(contained_first, contained_last) <
          Constants::TRACK_MIN_HITS)
        return false;

      Track track(*entry_hit, exit_hit);

      if (!good_track(track))
        return false;

      if (++n_good_track % sparsification_factor_)
        return false;

      auto time_spent = duration_cast<milliseconds>
        (high_resolution_clock::now() - t_start).count();

      if (trigger_use_timeout_ &&
          time_spent > trigger_timeout_)
      {
        result = true;
        trigger_entry_hit = *entry_hit;

        if (trigger_write_all_timeout_data_)
          trigger_exit_hit = hitmap["x_surface"].back();
        else
          trigger_exit_hit = *entry_hit;

        extra_info.set(Extra_Info::TIMEOUT_BIT, true);

        return true;
      }
      
      std::set<int> plane_set = {entry_hit->Plane().val};
      std::set<int> cell_set =  {entry_hit->Cell().val};

      uint64_t t_max_plane_gap = track.start().TDC().val;
      if (track.slope_plane_is_valid())
        t_max_plane_gap += abs(track.slope_time_plane()) * max_plane_gap_cut_;

      uint64_t t_max_cell_gap  = track.start().TDC().val;
      if (track.slope_cell_is_valid())
        t_max_cell_gap += abs(track.slope_time_cell()) * max_cell_gap_cut_;
      uint64_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,
         [&](novaddt::DAQHit const& hit)
         {
           if (hit.TDC().val > t_max_gap)
           {
             int64_t time_difference =
               static_cast<int64_t>(hit.TDC().val) -
               static_cast<int64_t>(track.start().TDC().val);

             if (hit.TDC().val > t_max_cell_gap)
             {
               double expected_cell =
               static_cast<double>(track.start().Cell().val);
               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().val > t_max_plane_gap)
             {
               double expected_plane =
                 static_cast<double>(track.start().Plane().val);
               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().val);
             plane_set.insert(hit.Plane().val);
           }   

           return false;
         });

      if (contains_gap_exceeding_cut)
        return false;

      plane_set.insert(exit_hit.Plane().val);
      if (plane_set.size() < 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().val);
      if (find_max_gap(cell_set) > max_cell_gap_cut_)
        return false;

      result = true;
      trigger_entry_hit = *entry_hit;
      trigger_exit_hit = exit_hit;

      return true;
    });
  }

  if (result)
  {
    auto t_begin = trigger_entry_hit.TDC().val - trigger_time_buffer_;
    auto t_end   = trigger_exit_hit.TDC().val  + trigger_time_buffer_;
    
    trigger_decisions_.
      emplace_back(t_begin, t_end - t_begin,
                   daqdataformats::trigID::TRIG_ID_DATA_SLOWMONO, 1,
                   extra_info.to_ulong());
  }
  
  return result;
}



void novaddt::smt::Trigger::split_by_view
(novaddt::HitList hits, std::map<std::string, novaddt::HitList> & hitmap) const
{
  auto bound = std::stable_partition
    (hits.begin(), hits.end(),
     [](novaddt::DAQHit const& h)
     { return h.View().val == daqchannelmap::X_VIEW; });

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

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



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

    auto y_match = std::find_if
      (y_first, y_last, [&](novaddt::DAQHit 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, [&](novaddt::DAQHit 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 novaddt::smt::Trigger::hits_are_view_matched
(novaddt::DAQHit const& x_hit, novaddt::DAQHit const& y_hit) const
{
  int plane_difference =
    static_cast<int>(x_hit.Plane().val) - static_cast<int>(y_hit.Plane().val);
  
  if (std::abs(plane_difference) < Constants::VIEW_MATCH_PLANE_CUT)
    return true;

  return false;
}



bool novaddt::smt::Trigger::hit_is_on_surface
(novaddt::DAQHit const& hit) const
{
  if (hit.Cell().val < Constants::SURFACE_MIN_CELL)
    return true;

  if (hit.Cell().val > Constants::SURFACE_MAX_CELL)
    return true;

  if (hit.Plane().val < Constants::SURFACE_MIN_PLANE)
    return true;

  if (hit.Plane().val > Constants::SURFACE_MAX_PLANE)
    return true;

  return false;
}



bool novaddt::smt::Trigger::good_track(Track const& track) const
{
  if (track.distance() < 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 novaddt::smt::Trigger::hit_is_on_road
(novaddt::DAQHit 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().val) -
     static_cast<double>(track.start().Plane().val)) +
    static_cast<double>(track.start().Cell().val);

  if (hit.Cell().val > road_center_cell + Constants::ROAD_HALF_WIDTH)
    return false;

  if (hit.Cell().val < road_center_cell - Constants::ROAD_HALF_WIDTH)
    return false;

  return true;
}



bool novaddt::smt::Trigger::hit_is_in_time_with_road
(novaddt::DAQHit 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;


  int64_t relative_hit_time =
    static_cast<int64_t>(hit.TDC().val) -
    static_cast<int64_t>(track.start().TDC().val);
  
  double expected_time = track.slope_time_plane() *
    (static_cast<double>(hit.Plane().val) -
     static_cast<double>(track.start().Plane().val));

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

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

  return true;
}



bool novaddt::smt::Trigger::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 novaddt::smt::Trigger::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;
}



bool novaddt::smt::Trigger::hit_is_between_track_end_points
(novaddt::DAQHit const& hit, Track const& track) const
{
  if (hit.Plane().val < track.min_plane())
    return false;

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

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

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

  return true;
}



novaddt::HitList::const_iterator novaddt::smt::Trigger::find
(novaddt::HitList::const_iterator begin,
 novaddt::HitList::const_iterator end, novaddt::TDC tdc) const
{
  return
    std::lower_bound(begin, end, tdc, novaddt::CompareDAQHit<novaddt::TDC>());
}



novaddt::HitList::const_iterator novaddt::smt::Trigger::find
(novaddt::HitList const& hits, novaddt::TDC tdc) const
{
  return find(hits.cbegin(), hits.cend(), tdc);
}