SlowMonopoleTrigger.cxx

Go to the documentation of this file.

#include "DDTTools/DDTExtensions.h"
#include "DDTTools/DDTHelpers.h"
#include "DDTMonopole/SlowMonopoleTrigger.h"
#include "DDTMonopole/Track.h"

#include <DDTBaseDataProducts/CompareDAQHit.h>
#include <SlowMonopoleTrigger/Constants.h>

#include <boost/format.hpp>

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

using namespace std::chrono;



mono::SlowMonopoleTrigger::SlowMonopoleTrigger
(novaddt::HitList const& h, unsigned long long const& event_time,
 Configuration const& c) :
  hits_(h),
  event_time_(event_time),
  max_plane_gap_cut_(c.max_plane_gap_cut),
  sparsification_factor_(c.sparsification_factor)
{
  max_cell_gap_cut_    = 20;
  track_max_beta_      = 5e-3;
  track_min_beta_      = 4e-5;
  trigger_time_buffer_ = 256;

  auto t0 = high_resolution_clock::now();
  algorithm();
  auto t1 = high_resolution_clock::now();

  parameter_["algorithm_time"] = duration_cast<milliseconds>(t1-t0).count();
  
  std::cout << "MF: Algorithm Time: "
            << duration_cast<milliseconds>(t1-t0).count() << " ms"
            << std::endl;
}



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



unsigned mono::SlowMonopoleTrigger::trigger_decision_time_buffer() const
{
  return trigger_time_buffer_;
}




std::map<std::string, novaddt::HitList>
mono::SlowMonopoleTrigger::essential_hit_lists() const
{
  return essential_hit_lists_;
}



bool mono::SlowMonopoleTrigger::parameter_exists(std::string const name) const
{
  if (parameter_.find(name) == parameter_.end())
    return false;

  return true;
}



double mono::SlowMonopoleTrigger::parameter(std::string const name)
{
  if (parameter_.find(name) == parameter_.end())
  {
    std::cerr << "\n\tParameter " << name << " is not available.\n"
              << std::endl;
    assert(false);
  }

  return parameter_[name];
}



bool mono::SlowMonopoleTrigger::algorithm()
{
  // Timing Cheat Sheet
  //
  // std::map<std::string, double> dur =
  //   { {"setup", 0}, {"find", 0}, {"track", 0}, {"sets", 0} };
  // auto t_start = high_resolution_clock::now();
  // dur["setup"] +=
  // duration_cast<nanoseconds>(high_resolution_clock::now() - t_start).count();

  bool result = false;
  novaddt::DAQHit trigger_entry_hit, trigger_exit_hit;

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

  // std::cout << "MF: (x hits, y hits, min_hits) = ("
  //        << hitmap["x"].size()
  //        << ", " << hitmap["y"].size()
  //        << ", " << min_hits_
  //        << ")" << 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);

  // std::cout << "MF: (x surface, x contained, min_hits) = ("
  //        << hitmap["x_surface"].size()
  //        << ", " << hitmap["x_contained"].size()
  //        << ", " << min_hits_
  //        << ")" << std::endl;

  // for (auto const& hit : hitmap["x_surface"])
  //   ddthelpers::print(hit, event_time_);

  // std::cout << "\n\n\nContained Hits:" << std::endl;
  // for (auto const& hit : hitmap["x_contained"])
  //   ddthelpers::print(hit, event_time_);
  
  if(hitmap["x_contained"].size() < novaddt::smt::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) <
          novaddt::smt::Constants::TRACK_MIN_HITS)
        return false;

      Track track(*entry_hit, exit_hit);

      // std::cout << "\nMF (Track): (distance, beta, good) = ("
      //                << track.distance_in_cm() << ", "
      //                << track.beta() << ", "
      //                << std::boolalpha << good_track(track) << ")" << std::endl;
      // ddthelpers::print(*entry_hit, event_time_);
      // ddthelpers::print(exit_hit, event_time_);

      if (!good_track(track))
        return false;

      if (++n_good_track % sparsification_factor_)
        return false;
      
      std::set<int> plane_set = {entry_hit->Plane().val};
      std::set<int> cell_set =  {entry_hit->Cell().val};

      unsigned long long t_max_plane_gap = track.start_ddt().TDC().val;
      if (track.slope_plane_is_valid())
        t_max_plane_gap += abs(track.slope_time_plane()) * max_plane_gap_cut_;

      unsigned long long t_max_cell_gap  = track.start_ddt().TDC().val;
      if (track.slope_cell_is_valid())
        t_max_cell_gap += abs(track.slope_time_cell()) * max_cell_gap_cut_;
      unsigned long long 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_ddt().TDC().val);

             if (hit.TDC().val > t_max_cell_gap)
             {
               double expected_cell =
               static_cast<double>(track.start_ddt().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_ddt().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;
         });

      // std::cout << "MF: contains_gap_exceeding_cut_ = "
      //                << std::boolalpha << contains_gap_exceeding_cut
      //                << std::endl;
      
      if (contains_gap_exceeding_cut)
        return false;

      // std::cout << "MF: contains_gap_exceeding_cut_ = "
      //                << std::boolalpha << contains_gap_exceeding_cut
      //                << std::endl;
      
      plane_set.insert(exit_hit.Plane().val);
      // std::cout << "MF: (plane_set_size, plane_max_gap) = ("
      //                << plane_set.size()
      //                << ", " << find_max_gap(plane_set)
      //                << ")" << std::endl;

      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().val);
      // std::cout << "MF: (cell_set_size, cell_max_gap) = ("
      //                << cell_set.size()
      //                << ", " << find_max_gap(cell_set)
      //                << ")" << std::endl;
      
      if (find_max_gap(cell_set) > max_cell_gap_cut_)
        return false;

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

      std::cout << "MF: Trigger Decision Hits:" << std::endl;
      ddthelpers::print(trigger_entry_hit, event_time_);
      ddthelpers::print(trigger_exit_hit, event_time_);
      std::cout << "MF: Trigger Decision Track:"
                << "\nMF: (distance, velocity) = ("
                << track.distance() << ", "
                << track.beta() << ")" << std::endl;
      
      return true;
    });
  }

  // std::cout << "MF: n_good_track = " << n_good_track << std::endl;

  
  if (result)
  {
    std::cout << "MF: Trigger Decision!" << std::endl;

    auto t0 = trigger_entry_hit.TDC().val;
    auto dt = trigger_exit_hit.TDC().val - t0;
    
    trigger_decisions_.
      emplace_back(t0 - 2 * trigger_time_buffer_,
                   dt + 4 * trigger_time_buffer_,
                   daqdataformats::trigID::TRIG_ID_DATA_SLOWMONO, 1);

    essential_hit_lists_["mf_entry_exit"].push_back(trigger_entry_hit);
    essential_hit_lists_["mf_entry_exit"].push_back(trigger_exit_hit);
  }
  
  return result;
}



void mono::SlowMonopoleTrigger::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 mono::SlowMonopoleTrigger::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 - novaddt::smt::Constants::VIEW_MATCH_TDC_CUT);
    auto y_last  =
      find(hitmap["y"],
           x_hit.TDC().val + novaddt::smt::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 mono::SlowMonopoleTrigger::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) <
      novaddt::smt::Constants::VIEW_MATCH_PLANE_CUT)
    return true;

  return false;
}



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

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

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

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

  return false;
}



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

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

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

  // For simplicity, let's punt on all events that cross fewer planes
  // than the maximum gap cut:
  if (abs(track.plane_difference()) <= max_plane_gap_cut_)
    return false;
  
  return true; 
}



bool mono::SlowMonopoleTrigger::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_ddt().Plane().val)) +
    static_cast<double>(track.start_ddt().Cell().val);

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

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

  return true;
}



bool mono::SlowMonopoleTrigger::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_ddt().TDC().val);
  
  double expected_time = track.slope_time_plane() *
    (static_cast<double>(hit.Plane().val) -
     static_cast<double>(track.start_ddt().Plane().val));

  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::SlowMonopoleTrigger::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);

  // std::cout << "MF: Collection = ";
  // for (auto const& element : collection)
  //   std::cout << element << " ";
  // std::cout << std::endl;
  // std::cout << "MF: (max_gap, max_gap_cut) = ("
  //        << find_max_gap(collection)
  //        << ", " << max_gap_cut
  //        << ")" << std::endl;
  
  // 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::SlowMonopoleTrigger::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 mono::SlowMonopoleTrigger::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 mono::SlowMonopoleTrigger::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 mono::SlowMonopoleTrigger::find
(novaddt::HitList const& hits, novaddt::TDC tdc) const
{
  return find(hits.cbegin(), hits.cend(), tdc);
}