PrimitiveMatching.cc

Go to the documentation of this file.

#include "L1Trigger/L1TMuonEndCap/interface/PrimitiveMatching.h"

#include "helper.h"  // to_hex, to_binary, merge_sort3

namespace {
  const int bw_fph = 13;  // bit width of ph, full precision
  const int bpow = 7;     // (1 << bpow) is count of input ranks
  const int invalid_ph_diff = 0x1ff;  // 511 (9-bit)
}


void PrimitiveMatching::configure(
    int verbose, int endcap, int sector, int bx,
    bool fixZonePhi, bool useNewZones,
    bool bugSt2PhDiff, bool bugME11Dupes
) {
  verbose_ = verbose;
  endcap_  = endcap;
  sector_  = sector;
  bx_      = bx;

  fixZonePhi_      = fixZonePhi;
  useNewZones_     = useNewZones;
  bugSt2PhDiff_    = bugSt2PhDiff;
  bugME11Dupes_    = bugME11Dupes;
}

void PrimitiveMatching::process(
    const std::deque<EMTFHitCollection>& extended_conv_hits,
    const emtf::zone_array<EMTFRoadCollection>& zone_roads,
    emtf::zone_array<EMTFTrackCollection>& zone_tracks
) const {

  // Function to update fs_history encoded in fs_segment
  auto update_fs_history = [](int fs_segment, int this_bx, int hit_bx) {
    // 0 for current BX, 1 for previous BX, 2 for BX before that
    int fs_history = this_bx - hit_bx;
    fs_segment |= ((fs_history & 0x3)<<4);
    return fs_segment;
  };

  // Function to update bt_history encoded in bt_segment
  auto update_bt_history = [](int bt_segment, int this_bx, int hit_bx) {
    // 0 for current BX, 1 for previous BX, 2 for BX before that
    int bt_history = this_bx - hit_bx;
    bt_segment |= ((bt_history & 0x3)<<5);
    return bt_segment;
  };


  // Exit if no roads
  int num_roads = 0;
  for (const auto& roads : zone_roads)
    num_roads += roads.size();
  bool early_exit = (num_roads == 0);

  if (early_exit)
    return;


  if (verbose_ > 0) {  // debug
    for (const auto& roads : zone_roads) {
      for (const auto& road : roads) {
        std::cout << "pattern on match input: z: " << road.Zone()-1 << " r: " << road.Winner()
            << " ph_num: " << road.Key_zhit() << " ph_q: " << to_hex(road.Quality_code())
            << " ly: " << to_binary(road.Layer_code(), 3) << " str: " << to_binary(road.Straightness(), 3)
            << std::endl;
      }
    }
  }

  // Organize converted hits by (zone, station)
  std::array<EMTFHitCollection, emtf::NUM_ZONES*emtf::NUM_STATIONS> zs_conv_hits;

  bool use_fs_zone_code = true;  // use zone code as in firmware find_segment module

  std::deque<EMTFHitCollection>::const_iterator ext_conv_hits_it  = extended_conv_hits.begin();
  std::deque<EMTFHitCollection>::const_iterator ext_conv_hits_end = extended_conv_hits.end();

  for (; ext_conv_hits_it != ext_conv_hits_end; ++ext_conv_hits_it) {
    EMTFHitCollection::const_iterator conv_hits_it  = ext_conv_hits_it->begin();
    EMTFHitCollection::const_iterator conv_hits_end = ext_conv_hits_it->end();

    for (; conv_hits_it != conv_hits_end; ++conv_hits_it) {
      int istation = conv_hits_it->Station()-1;
      int zone_code = conv_hits_it->Zone_code();  // decide based on original zone code
      if (use_fs_zone_code)
        zone_code = conv_hits_it->FS_zone_code();  // decide based on new zone code

      // A hit can go into multiple zones
      for (int izone = 0; izone < emtf::NUM_ZONES; ++izone) {
        if (!zone_roads.at(izone).empty()) {

          if (zone_code & (1<<izone)) {
            const int zs = (izone*emtf::NUM_STATIONS) + istation;
            zs_conv_hits.at(zs).push_back(*conv_hits_it);

            // Update fs_history and bt_history depending on the processor BX
            // This update only goes into the hits associated to a track, it does not affect the original hit collection
            EMTFHit& conv_hit = zs_conv_hits.at(zs).back();   // pass by reference
            int old_fs_segment = conv_hit.FS_segment();
            int new_fs_segment = update_fs_history(old_fs_segment, bx_, conv_hit.BX());
            conv_hit.set_fs_segment( new_fs_segment );

            int old_bt_segment = conv_hit.BT_segment();
            int new_bt_segment = update_bt_history(old_bt_segment, bx_, conv_hit.BX());
            conv_hit.set_bt_segment( new_bt_segment );
          }
        }
      }

    }  // end loop over conv_hits
  }  // end loop over extended_conv_hits

  if (verbose_ > 1) {  // debug
    for (int izone = 0; izone < emtf::NUM_ZONES; ++izone) {
      for (int istation = 0; istation < emtf::NUM_STATIONS; ++istation) {
        const int zs = (izone*emtf::NUM_STATIONS) + istation;
        for (const auto& conv_hit : zs_conv_hits.at(zs)) {
          std::cout << "z: " << izone << " st: " << istation+1 << " cscid: " << conv_hit.CSC_ID()
              << " ph_zone_phi: " << conv_hit.Zone_hit() << " ph_low_prec: " << (conv_hit.Zone_hit()<<5)
              << " ph_high_prec: " << conv_hit.Phi_fp() << " ph_high_low_diff: " << (conv_hit.Phi_fp() - (conv_hit.Zone_hit()<<5))
              << std::endl;
        }
      }
    }
  }

  // Keep the best phi difference for every road by (zone, station)
  std::array<std::vector<hit_sort_pair_t>, emtf::NUM_ZONES*emtf::NUM_STATIONS> zs_phi_differences;

  // Get the best-matching hits by comparing phi difference between
  // pattern and segment
  for (int izone = 0; izone < emtf::NUM_ZONES; ++izone) {
    for (int istation = 0; istation < emtf::NUM_STATIONS; ++istation) {
      const int zs = (izone*emtf::NUM_STATIONS) + istation;

      // This leaves zone_roads.at(izone) and zs_conv_hits.at(zs) unchanged
      // zs_phi_differences.at(zs) gets filled with a pair of <phi_diff, conv_hit> for the
      // conv_hit with the lowest phi_diff from the pattern in this station and zone
      process_single_zone_station(
          izone+1, istation+1,
          zone_roads.at(izone),
          zs_conv_hits.at(zs),
          zs_phi_differences.at(zs)
      );

      if (not(zone_roads.at(izone).size() == zs_phi_differences.at(zs).size()))
    { edm::LogError("L1T") << "zone_roads.at(izone).size() = " << zone_roads.at(izone).size()
                   << ", zs_phi_differences.at(zs).size() = " << zs_phi_differences.at(zs).size(); return; }
    }  // end loop over stations
  }  // end loop over zones

  if (verbose_ > 1) {  // debug
    for (int izone = 0; izone < emtf::NUM_ZONES; ++izone) {
      const auto& roads = zone_roads.at(izone);
      for (unsigned iroad = 0; iroad < roads.size(); ++iroad) {
        const auto& road = roads.at(iroad);
        for (int istation = 0; istation < emtf::NUM_STATIONS; ++istation) {
          const int zs = (izone*emtf::NUM_STATIONS) + istation;
          int ph_diff = zs_phi_differences.at(zs).at(iroad).first;
          std::cout << "find seg: z: " << road.Zone()-1 << " r: " << road.Winner()
              << " st: " << istation << " ph_diff: " << ph_diff
              << std::endl;
        }
      }
    }
  }


  // Build all tracks in each zone
  for (int izone = 0; izone < emtf::NUM_ZONES; ++izone) {
    const EMTFRoadCollection& roads = zone_roads.at(izone);

    for (unsigned iroad = 0; iroad < roads.size(); ++iroad) {
      const EMTFRoad& road = roads.at(iroad);

      // Create a track
      EMTFTrack track;
      track.set_endcap   ( road.Endcap() );
      track.set_sector   ( road.Sector() );
      track.set_sector_idx ( road.Sector_idx() );
      track.set_bx       ( road.BX() );
      track.set_zone     ( road.Zone() );
      track.set_ph_num   ( road.Key_zhit() );
      track.set_ph_q     ( road.Quality_code() );
      track.set_rank     ( road.Quality_code() );
      track.set_winner   ( road.Winner() );

      track.clear_Hits();

      // Insert hits
      for (int istation = 0; istation < emtf::NUM_STATIONS; ++istation) {
        const int zs = (izone*emtf::NUM_STATIONS) + istation;

        const EMTFHitCollection& conv_hits = zs_conv_hits.at(zs);
        int       ph_diff      = zs_phi_differences.at(zs).at(iroad).first;
        hit_ptr_t conv_hit_ptr = zs_phi_differences.at(zs).at(iroad).second;

        if (ph_diff != invalid_ph_diff) {
          // Inserts the conv_hit with the lowest phi_diff, as well as its duplicate
          // (same strip and phi, different wire and theta), if a duplicate exists
          insert_hits(conv_hit_ptr, conv_hits, track);
        }
      }

      if (fixZonePhi_) {
        if (not(!track.Hits().empty()))
      { edm::LogError("L1T") << "track.Hits().empty() = " << track.Hits().empty(); return; }
      }

      // Output track
      zone_tracks.at(izone).push_back(track);

    }  // end loop over roads
  }  // end loop over zones

  if (verbose_ > 0) {  // debug
    for (const auto& tracks : zone_tracks) {
      for (const auto& track : tracks) {
        for (const auto& hit : track.Hits()) {
          std::cout << "match seg: z: " << track.Zone()-1 << " pat: " << track.Winner() <<  " st: " << hit.Station()
              << " vi: " << to_binary(0b1, 2) << " hi: " << ((hit.FS_segment()>>4) & 0x3)
              << " ci: " << ((hit.FS_segment()>>1) & 0x7) << " si: " << (hit.FS_segment() & 0x1)
              << " ph: " << hit.Phi_fp() << " th: " << hit.Theta_fp()
              << std::endl;
        }
      }
    }
  }

}

void PrimitiveMatching::process_single_zone_station(
    int zone, int station,
    const EMTFRoadCollection& roads,
    const EMTFHitCollection& conv_hits,
    std::vector<hit_sort_pair_t>& phi_differences
) const {
  // max phi difference between pattern and segment
  // This doesn't depend on the pattern straightness - any hit within the max phi difference may match
  int max_ph_diff = (station == 1) ? 15 : 7;
  //int bw_ph_diff = (station == 1) ? 5 : 4; // ph difference bit width
  //int invalid_ph_diff = (station == 1) ? 31 : 15;  // invalid difference

  if (fixZonePhi_) {
    if (station == 1) {
      max_ph_diff = 496;  // width of pattern in ME1 + rounding error 15*32+16
      //bw_ph_diff = 9;
      //invalid_ph_diff = 0x1ff;
    } else if (station == 2) {
      if (bugSt2PhDiff_)
        max_ph_diff = 16;   // just rounding error for ME2 (pattern must match ME2 hit phi if there was one)
      else
        max_ph_diff = 240;  // same as ME3,4
      //bw_ph_diff = 5;
      //invalid_ph_diff = 0x1f;
    } else {
      max_ph_diff = 240;  // width of pattern in ME3,4 + rounding error 7*32+16
      //bw_ph_diff = 8;
      //invalid_ph_diff = 0xff;
    }
  }

  auto abs_diff = [](int a, int b) { return std::abs(a-b); };

  // Simple sort by ph_diff
  struct {
    typedef hit_sort_pair_t value_type;
    bool operator()(const value_type& lhs, const value_type& rhs) const {
      return lhs.first <= rhs.first;
    }
  } less_ph_diff_cmp;

  // Emulation of FW sorting with 3-way comparator
  struct {
    typedef hit_sort_pair_t value_type;
    int operator()(const value_type& a, const value_type& b, const value_type& c) const {
      int r = 0;
      r |= bool(a.first <= b.first);
      r <<= 1;
      r |= bool(b.first <= c.first);
      r <<= 1;
      r |= bool(c.first <= a.first);

      int rr = 0;
      switch(r) {
      //case 0b000 : rr = 3; break;  // invalid
      case 0b001 : rr = 2; break;  // c
      case 0b010 : rr = 1; break;  // b
      case 0b011 : rr = 1; break;  // b
      case 0b100 : rr = 0; break;  // a
      case 0b101 : rr = 2; break;  // c
      case 0b110 : rr = 0; break;  // a
      //case 0b111 : rr = 0; break;  // invalid
      default    : rr = 0; break;
      }
      return rr;
    }
  } less_ph_diff_cmp3;


  // ___________________________________________________________________________
  // For each road, find the segment with min phi difference in every station

  EMTFRoadCollection::const_iterator roads_it  = roads.begin();
  EMTFRoadCollection::const_iterator roads_end = roads.end();

  for (; roads_it != roads_end; ++roads_it) {
    int ph_pat = roads_it->Key_zhit();     // pattern key phi value
    int ph_q   = roads_it->Quality_code(); // pattern quality code
    if (not(ph_pat >= 0 && ph_q > 0))
      { edm::LogError("L1T") << "ph_pat = " << ph_pat << ", ph_q = " << ph_q; return; }

    if (fixZonePhi_) {
      ph_pat <<= 5;  // add missing 5 lower bits to pattern phi
    }

    std::vector<hit_sort_pair_t> tmp_phi_differences;

    EMTFHitCollection::const_iterator conv_hits_it  = conv_hits.begin();
    EMTFHitCollection::const_iterator conv_hits_end = conv_hits.end();

    for (; conv_hits_it != conv_hits_end; ++conv_hits_it) {
      int ph_seg     = conv_hits_it->Phi_fp();  // ph from segments
      int ph_seg_red = ph_seg >> (bw_fph-bpow-1);  // remove unused low bits
      if (not(ph_seg >= 0))
    { edm::LogError("L1T") << "ph_seg = " << ph_seg; return; }

      if (fixZonePhi_) {
        ph_seg_red = ph_seg;  // use full-precision phi
      }

      // Get abs phi difference
      int ph_diff = abs_diff(ph_pat, ph_seg_red);
      if (ph_diff > max_ph_diff)
        ph_diff = invalid_ph_diff;  // difference is too high, cannot be the same pattern

      if (ph_diff != invalid_ph_diff)
        tmp_phi_differences.push_back(std::make_pair(ph_diff, conv_hits_it));  // make a key-value pair
    }

    // _________________________________________________________________________
    // Sort to find the segment with min phi difference

    if (!tmp_phi_differences.empty()) {
      // Because the sorting is sensitive to FW ordering, use the exact FW sorting.
      // This implementation still slightly differs from FW because I prefer to
      // use a sorting function that is as generic as possible.
      bool use_fw_sorting = true;

      if (useNewZones_)  use_fw_sorting = false;

      if (use_fw_sorting) {
        // zone_cham = 4 for [fs_01, fs_02, fs_03, fs_11], or 7 otherwise
        // tot_diff = 27 or 45 in FW; it is 27 or 54 in the C++ merge_sort3 impl
        const int max_drift = 3; // should use bxWindow from the config
        const int zone_cham = ((zone == 1 && (2 <= station && station <= 4)) || (zone == 2 && station == 2)) ? 4 : 7;
        const int seg_ch    = 2;
        const int tot_diff  = (max_drift*zone_cham*seg_ch) + ((zone_cham == 4) ? 3 : 12);  // provide padding for 3-input comparators

        std::vector<hit_sort_pair_t> fw_sort_array(tot_diff, std::make_pair(invalid_ph_diff, conv_hits_end));

        // FW doesn't check if the hit is CSC or RPC
        std::vector<hit_sort_pair_t>::const_iterator phdiffs_it  = tmp_phi_differences.begin();
        std::vector<hit_sort_pair_t>::const_iterator phdiffs_end = tmp_phi_differences.end();

        for (; phdiffs_it != phdiffs_end; ++phdiffs_it) {
          //int ph_diff    = phdiffs_it->first;
          int fs_segment = phdiffs_it->second->FS_segment();

          // Calculate the index to put into the fw_sort_array
          int fs_history = ((fs_segment>>4) & 0x3);
          int fs_chamber = ((fs_segment>>1) & 0x7);
          fs_segment = (fs_segment & 0x1);
          unsigned fw_sort_array_index = (fs_history * zone_cham * seg_ch) + (fs_chamber * seg_ch) + fs_segment;

          if (not(fs_history < max_drift && fs_chamber < zone_cham && fs_segment < seg_ch))
        { edm::LogError("L1T") << "fs_history = " << fs_history << ", max_drift = " << max_drift
                   << ", fs_chamber = " << fs_chamber << ", zone_cham = " << zone_cham
                   << ", fs_segment = " << fs_segment << ", seg_ch = " << seg_ch; return; }
          if (not(fw_sort_array_index < fw_sort_array.size()))
        { edm::LogError("L1T") << "fw_sort_array_index = " << fw_sort_array_index
                   << ", fw_sort_array.size() = " << fw_sort_array.size(); return; }
          fw_sort_array.at(fw_sort_array_index) = *phdiffs_it;
        }

        // Debug
        //std::cout << "phdiffs" << std::endl;
        //for (unsigned i = 0; i < fw_sort_array.size(); ++i)
        //  std::cout << fw_sort_array.at(i).first << " ";
        //std::cout << std::endl;

        // Debug
        //std::cout << "Before sort" << std::endl;
        //for (unsigned i = 0; i < fw_sort_array.size(); ++i)
        //  std::cout << fw_sort_array.at(i).second->FS_segment() << " ";
        //std::cout << std::endl;

        // Find the best phi difference according to FW sorting
        //merge_sort3(fw_sort_array.begin(), fw_sort_array.end(), less_ph_diff_cmp, less_ph_diff_cmp3);
        merge_sort3_with_hint(fw_sort_array.begin(), fw_sort_array.end(), less_ph_diff_cmp, less_ph_diff_cmp3, ((tot_diff == 54) ? tot_diff/2 : tot_diff/3));

        // Store the best phi difference
        phi_differences.push_back(fw_sort_array.front());

        // Debug
        //std::cout << "After sort" << std::endl;
        //for (unsigned i = 0; i < fw_sort_array.size(); ++i)
        //  std::cout << fw_sort_array.at(i).second->FS_segment() << " ";
        //std::cout << std::endl;

      } else {  // use C++ sorting
        struct {
          typedef hit_sort_pair_t value_type;
          bool operator()(const value_type& lhs, const value_type& rhs) const {
            // If different types, prefer CSC over RPC; else prefer the closer hit in dPhi
            if (lhs.second->Subsystem() != rhs.second->Subsystem())
              return (lhs.second->Subsystem() == TriggerPrimitive::kCSC);
            else
              return lhs.first <= rhs.first;
          }
        } tmp_less_ph_diff_cmp;

        // Find best phi difference
        std::stable_sort(tmp_phi_differences.begin(), tmp_phi_differences.end(), tmp_less_ph_diff_cmp);

        // Store the best phi difference
        phi_differences.push_back(tmp_phi_differences.front());
      }

    } else {
      // No segment found
      phi_differences.push_back(std::make_pair(invalid_ph_diff, conv_hits_end));  // make a key-value pair
    }

  }  // end loop over roads
}

void PrimitiveMatching::insert_hits(
    hit_ptr_t conv_hit_ptr, const EMTFHitCollection& conv_hits,
    EMTFTrack& track
) const {
  EMTFHitCollection::const_iterator conv_hits_it  = conv_hits.begin();
  EMTFHitCollection::const_iterator conv_hits_end = conv_hits.end();

  const bool is_csc_me11 = (conv_hit_ptr->Subsystem() == TriggerPrimitive::kCSC) &&
      (conv_hit_ptr->Station() == 1) && (conv_hit_ptr->Ring() == 1 || conv_hit_ptr->Ring() == 4);

  // Find all possible duplicated hits, insert them
  for (; conv_hits_it != conv_hits_end; ++conv_hits_it) {
    const EMTFHit& conv_hit_i = *conv_hits_it;
    const EMTFHit& conv_hit_j = *conv_hit_ptr;

    // All these must match: [bx_history][station][chamber][segment]
    if (
      (conv_hit_i.Subsystem()  == conv_hit_j.Subsystem()) &&
      (conv_hit_i.PC_station() == conv_hit_j.PC_station()) &&
      (conv_hit_i.PC_chamber() == conv_hit_j.PC_chamber()) &&
      (conv_hit_i.Ring()       == conv_hit_j.Ring()) &&  // because of ME1/1
      (conv_hit_i.Strip()      == conv_hit_j.Strip()) &&
      //(conv_hit_i.Wire()       == conv_hit_j.Wire()) &&
      (conv_hit_i.Pattern()    == conv_hit_j.Pattern()) &&
      (conv_hit_i.BX()         == conv_hit_j.BX()) &&
      ( (conv_hit_i.Is_RPC()     == false) ||
    ( (conv_hit_i.Strip_low()  == conv_hit_j.Strip_low()) &&
      (conv_hit_i.Strip_hi()   == conv_hit_j.Strip_hi()) &&
      (conv_hit_i.Roll()       == conv_hit_j.Roll()) &&
      (conv_hit_i.Phi_fp()     == conv_hit_j.Phi_fp()) &&
      (conv_hit_i.Theta_fp()   == conv_hit_j.Theta_fp()) ) ) &&
      true
    ) {
      // All duplicates with the same strip but different wire must have same phi_fp
      if (not(conv_hit_i.Phi_fp() == conv_hit_j.Phi_fp()))
    { edm::LogError("L1T") << "conv_hit_i.Phi_fp() = " << conv_hit_i.Phi_fp() 
                   << ", conv_hit_j.Phi_fp() = " << conv_hit_j.Phi_fp(); return; }

      track.push_Hit( conv_hit_i );

    } else if (
      (bugME11Dupes_ && is_csc_me11) &&  // if reproduce ME1/1 theta duplication bug, do not check 'ring', 'strip' and 'pattern'
      (conv_hit_i.Subsystem()  == conv_hit_j.Subsystem()) &&
      (conv_hit_i.PC_station() == conv_hit_j.PC_station()) &&
      (conv_hit_i.PC_chamber() == conv_hit_j.PC_chamber()) &&
      //(conv_hit_i.Ring()       == conv_hit_j.Ring()) &&  // because of ME1/1
      //(conv_hit_i.Strip()      == conv_hit_j.Strip()) &&
      //(conv_hit_i.Wire()       == conv_hit_j.Wire()) &&
      //(conv_hit_i.Pattern()    == conv_hit_j.Pattern()) &&
      (conv_hit_i.BX()         == conv_hit_j.BX()) &&
      //(conv_hit_i.Strip_low()  == conv_hit_j.Strip_low()) && // For RPC clusters
      //(conv_hit_i.Strip_hi()   == conv_hit_j.Strip_hi()) &&  // For RPC clusters
      //(conv_hit_i.Roll()       == conv_hit_j.Roll()) &&      // For RPC clusters
      true
    ) {
      // Dirty hack
      EMTFHit tmp_hit = conv_hit_j;
      tmp_hit.set_theta_fp( conv_hit_i.Theta_fp() );
      track.push_Hit( tmp_hit );
    }
  }

  // Sort by station
  struct {
    typedef EMTFHit value_type;
    bool operator()(const value_type& lhs, const value_type& rhs) const {
      return lhs.Station() < rhs.Station();
    }
  } less_station_cmp;

  EMTFHitCollection tmp_hits = track.Hits();
  std::stable_sort(tmp_hits.begin(), tmp_hits.end(), less_station_cmp);
  track.set_Hits( tmp_hits );
}