BestTrackSelection.cc

Go to the documentation of this file.

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

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

void BestTrackSelection::configure(int verbose,
                                   int endcap,
                                   int sector,
                                   int bx,
                                   int bxWindow,
                                   int maxRoadsPerZone,
                                   int maxTracks,
                                   bool useSecondEarliest,
                                   bool bugSameSectorPt0) {
  verbose_ = verbose;
  endcap_ = endcap;
  sector_ = sector;
  bx_ = bx;

  bxWindow_ = bxWindow;
  maxRoadsPerZone_ = maxRoadsPerZone;
  maxTracks_ = maxTracks;
  useSecondEarliest_ = useSecondEarliest;
  bugSameSectorPt0_ = bugSameSectorPt0;
}

void BestTrackSelection::process(const std::deque<EMTFTrackCollection>& extended_best_track_cands,
                                 EMTFTrackCollection& best_tracks) const {
  int num_cands = 0;
  for (const auto& cands : extended_best_track_cands) {
    for (const auto& cand : cands) {
      if (cand.Rank() > 0) {
        num_cands += 1;
      }
    }
  }
  bool early_exit = (num_cands == 0);

  if (early_exit)
    return;

  if (!useSecondEarliest_) {
    cancel_one_bx(extended_best_track_cands, best_tracks);
  } else {
    cancel_multi_bx(extended_best_track_cands, best_tracks);
  }

  if (verbose_ > 0) {  // debug
    for (const auto& track : best_tracks) {
      std::cout << "track: " << track.Winner() << " rank: " << to_hex(track.Rank())
                << " ph_deltas: " << array_as_string(track.PtLUT().delta_ph)
                << " th_deltas: " << array_as_string(track.PtLUT().delta_th) << " phi: " << track.Phi_fp()
                << " theta: " << track.Theta_fp() << " cpat: " << array_as_string(track.PtLUT().cpattern)
                << " bx: " << track.BX() << std::endl;
      for (int i = 0; i < emtf::NUM_STATIONS + 1; ++i) {  // stations 0-4
        if (track.PtLUT().bt_vi[i] != 0)
          std::cout << ".. track segments: st: " << i << " v: " << track.PtLUT().bt_vi[i]
                    << " h: " << track.PtLUT().bt_hi[i] << " c: " << track.PtLUT().bt_ci[i]
                    << " s: " << track.PtLUT().bt_si[i] << std::endl;
      }
    }
  }
}

void BestTrackSelection::cancel_one_bx(const std::deque<EMTFTrackCollection>& extended_best_track_cands,
                                       EMTFTrackCollection& best_tracks) const {
  const int max_z = emtf::NUM_ZONES;   // = 4 zones
  const int max_n = maxRoadsPerZone_;  // = 3 candidates per zone
  const int max_zn = max_z * max_n;    // = 12 total candidates
  if (not(maxTracks_ <= max_zn)) {
    edm::LogError("L1T") << "maxTracks_ = " << maxTracks_ << ", max_zn = " << max_zn;
    return;
  }

  // Emulate the arrays used in firmware
  typedef std::array<int, 3> segment_ref_t;
  std::vector<std::vector<segment_ref_t> > segments(max_zn,
                                                    std::vector<segment_ref_t>());  // 2D array [zn][num segments]

  std::vector<std::vector<bool> > larger(max_zn, std::vector<bool>(max_zn, false));  // 2D array [zn][zn]
  std::vector<std::vector<bool> > winner(max_zn, std::vector<bool>(max_zn, false));

  std::vector<bool> exists(max_zn, false);  // 1D array [zn]
  std::vector<bool> killed(max_zn, false);
  std::vector<int> rank(max_zn, 0);
  //std::vector<int>  good_bx(max_zn, 0);

  // Initialize arrays: rank, segments
  for (int z = 0; z < max_z; ++z) {
    const EMTFTrackCollection& tracks = extended_best_track_cands.at(z);
    const int ntracks = tracks.size();
    if (not(ntracks <= max_n)) {
      edm::LogError("L1T") << "ntracks = " << ntracks << ", max_n = " << max_n;
      return;
    }

    for (int n = 0; n < ntracks; ++n) {
      const int zn = (n * max_z) + z;  // for (i = 0; i < 12; i = i+1) rank[i%4][i/4]
      const EMTFTrack& track = tracks.at(n);

      rank.at(zn) = track.Rank();

      for (const auto& conv_hit : track.Hits()) {
        if (not(conv_hit.Valid())) {
          edm::LogError("L1T") << "conv_hit.Valid() = " << conv_hit.Valid();
          return;
        }

        // A segment identifier (chamber, strip, bx)
        const segment_ref_t segment = {{conv_hit.PC_station() * 9 + conv_hit.PC_chamber(),
                                        (conv_hit.PC_segment() % 2),
                                        0}};  // FW doesn't check whether a segment is CSC or RPC
        //int chamber_index  = int(conv_hit.Subsystem() == TriggerPrimitive::kRPC)*9*5;
        //chamber_index     += conv_hit.PC_station()*9;
        //chamber_index     += conv_hit.PC_chamber();
        //const segment_ref_t segment = {{chamber_index, conv_hit.Strip(), 0}};
        segments.at(zn).push_back(segment);
      }
    }  // end loop over n
  }    // end loop over z

  // Simultaneously compare each rank with each other
  int i = 0, j = 0, ri = 0, rj = 0, gt = 0, eq = 0, sum = 0;

  for (i = 0; i < max_zn; ++i) {
    for (j = 0; j < max_zn; ++j) {
      larger[i][j] = false;
    }
    larger[i][i] = true;  // result of comparison with itself
    //ri = rank[i%4][i/4]; // first index loops zone, second loops candidate. Zone loops faster, so we give equal priority to zones
    ri = rank[i];

    for (j = 0; j < max_zn; ++j) {
      // i&j bits show which rank is larger
      // the comparison scheme below avoids problems
      // when there are two or more tracks with the same rank
      //rj = rank[j%4][j/4];
      rj = rank[j];
      gt = ri > rj;
      eq = ri == rj;
      if ((i < j && (gt || eq)) || (i > j && gt))
        larger[i][j] = true;
    }
    // "larger" array shows the result of comparison for each rank

    // track exists if quality != 0
    exists[i] = (ri != 0);
  }

  // ghost cancellation, only in the current BX so far
  for (i = 0; i < max_zn - 1; ++i) {    // candidate loop
    for (j = i + 1; j < max_zn; ++j) {  // comparison candidate loop
      int shared_segs = 0;

      // count shared segments
      for (const auto& isegment : segments.at(i)) {  // loop over all pairs of hits
        for (const auto& jsegment : segments.at(j)) {
          if (isegment == jsegment) {  // same chamber and same segment
            shared_segs += 1;
          }
        }
      }

      if (shared_segs > 0) {  // a single shared segment means it's ghost
        // kill candidate that has lower rank
        if (larger[i][j])
          killed[j] = true;
        else
          killed[i] = true;
      }
    }
  }

  // remove ghosts according to kill mask
  //exists = exists & (~kill1);
  for (i = 0; i < max_zn; ++i) {
    exists[i] = exists[i] & (!killed[i]);
  }

  bool anything_exists = (std::find(exists.begin(), exists.end(), 1) != exists.end());
  if (!anything_exists)
    return;

  // update "larger" array
  for (i = 0; i < max_zn; ++i) {
    for (j = 0; j < max_zn; ++j) {
      //if  (exists[i]) larger[i] = larger[i] | (~exists); // if this track exists make it larger than all non-existing tracks
      //else  larger[i] = 0; // else make it smaller than anything
      if (exists[i])
        larger[i][j] = larger[i][j] | (!exists[j]);
      else
        larger[i][j] = false;
    }
  }

  if (verbose_ > 0) {  // debug
    std::cout << "exists: ";
    for (i = max_zn - 1; i >= 0; --i) {
      std::cout << exists[i];
    }
    std::cout << std::endl;
    std::cout << "killed: ";
    for (i = max_zn - 1; i >= 0; --i) {
      std::cout << killed[i];
    }
    std::cout << std::endl;
    for (j = 0; j < max_zn; ++j) {
      std::cout << "larger: ";
      for (i = max_zn - 1; i >= 0; --i) {
        std::cout << larger[j][i];
      }
      std::cout << std::endl;
    }
  }

  // count zeros in the comparison results. The best track will have none, the next will have one, the third will have two
  // skip the bits corresponding to the comparison of the track with itself
  for (i = 0; i < max_zn; ++i) {
    sum = 0;
    for (j = 0; j < max_zn; ++j) {
      if (larger[i][j] == 0)
        sum += 1;
    }

    if (sum < maxTracks_) {
      winner[sum][i] = true;  // assign positional winner codes
    }

    if (sum < maxTracks_ && bugSameSectorPt0_ && sum > 0) {
      // just keep the best track and kill the rest of them
      winner[sum][i] = false;
    }
  }

  // Output best tracks according to winner signals
  best_tracks.clear();

  for (int o = 0; o < maxTracks_; ++o) {  // output candidate loop
    int z = 0, n = 0;
    for (i = 0; i < max_zn; ++i) {  // winner bit loop
      if (winner[o][i]) {
        n = i / max_z;
        z = i % max_z;

        const EMTFTrackCollection& tracks = extended_best_track_cands.at(z);
        const EMTFTrack& track = tracks.at(n);
        best_tracks.push_back(track);

        // Update winner, BX
        best_tracks.back().set_track_num(best_tracks.size() - 1);
        best_tracks.back().set_winner(o);
        best_tracks.back().set_bx(best_tracks.back().First_BX());
      }
    }
  }
}

void BestTrackSelection::cancel_multi_bx(const std::deque<EMTFTrackCollection>& extended_best_track_cands,
                                         EMTFTrackCollection& best_tracks) const {
  const int max_h = bxWindow_;         // = 3 bx history
  const int max_z = emtf::NUM_ZONES;   // = 4 zones
  const int max_n = maxRoadsPerZone_;  // = 3 candidates per zone
  const int max_zn = max_z * max_n;    // = 12 total candidates
  const int max_hzn = max_h * max_zn;  // = 36 total candidates
  if (not(maxTracks_ <= max_hzn)) {
    edm::LogError("L1T") << "maxTracks_ = " << maxTracks_ << ", max_hzn = " << max_hzn;
    return;
  }

  const int delayBX = bxWindow_ - 1;
  const int num_h = extended_best_track_cands.size() / max_z;  // num of bx history so far

  // Emulate the arrays used in firmware
  typedef std::array<int, 3> segment_ref_t;
  std::vector<std::vector<segment_ref_t> > segments(max_hzn,
                                                    std::vector<segment_ref_t>());  // 2D array [hzn][num segments]

  std::vector<std::vector<bool> > larger(max_hzn, std::vector<bool>(max_hzn, false));  // 2D array [hzn][hzn]
  std::vector<std::vector<bool> > winner(max_hzn, std::vector<bool>(max_hzn, false));

  std::vector<bool> exists(max_hzn, false);  // 1D array [hzn]
  std::vector<bool> killed(max_hzn, false);
  std::vector<int> rank(max_hzn, 0);
  std::vector<int> good_bx(max_hzn, 0);

  // Initialize arrays: rank, good_bx, segments
  for (int h = 0; h < num_h; ++h) {
    // extended_best_track_cands[0..3] has 4 zones for road/pattern BX-0 (i.e. current) with possible tracks from [BX-2, BX-1, BX-0]
    // extended_best_track_cands[4..7] has 4 zones for road/pattern BX-1 with possible tracks from [BX-3, BX-2, BX-1]
    // extended_best_track_cands[8..11] has 4 zones for road/pattern BX-2 with possible tracks from [BX-4, BX-3, BX-2]

    for (int z = 0; z < max_z; ++z) {
      const EMTFTrackCollection& tracks = extended_best_track_cands.at(h * max_z + z);
      const int ntracks = tracks.size();
      if (not(ntracks <= max_n)) {
        edm::LogError("L1T") << "ntracks = " << ntracks << ", max_n = " << max_n;
        return;
      }

      for (int n = 0; n < ntracks; ++n) {
        const int hzn = (h * max_z * max_n) + (n * max_z) + z;  // for (i = 0; i < 12; i = i+1) rank[i%4][i/4]
        const EMTFTrack& track = tracks.at(n);
        int cand_bx = track.Second_BX();
        cand_bx -= (bx_ - delayBX);  // convert track.second_bx=[BX-2, BX-1, BX-0] --> cand_bx=[0,1,2]

        rank.at(hzn) = track.Rank();
        if (cand_bx == 0)
          good_bx.at(hzn) = 1;  // kill this rank if it's not the right BX

        for (const auto& conv_hit : track.Hits()) {
          if (not(conv_hit.Valid())) {
            edm::LogError("L1T") << "conv_hit.Valid() = " << conv_hit.Valid();
            return;
          }

          // Notes from Alex (2017-03-16):
          //
          //     What happens currently is that RPC hits are inserted instead of
          // CSC hits, if CSC hits are missing. So the track address will point
          // at a CSC chamber, but the actual hit may have come from
          // corresponding RPC located behind it. If the same substitution
          // happened in the neighboring sector, then the cancellation in uGMT
          // will work correctly. If the substitution happened in only one
          // sector, then RPC hit may "trump" a CSC hit, or vice versa.

          // A segment identifier (chamber, strip, bx)
          const segment_ref_t segment = {{conv_hit.PC_station() * 9 + conv_hit.PC_chamber(),
                                          (conv_hit.PC_segment() % 2),
                                          conv_hit.BX()}};  // FW doesn't check whether a segment is CSC or RPC
          //int chamber_index  = int(conv_hit.Subsystem() == TriggerPrimitive::kRPC)*9*5;
          //chamber_index     += conv_hit.PC_station()*9;
          //chamber_index     += conv_hit.PC_chamber();
          //const segment_ref_t segment = {{chamber_index, conv_hit.Strip(), conv_hit.BX()}};
          segments.at(hzn).push_back(segment);
        }
      }  // end loop over n
    }    // end loop over z
  }      // end loop over h

  // Simultaneously compare each rank with each other
  int i = 0, j = 0, ri = 0, rj = 0, sum = 0;

  for (i = 0; i < max_hzn; ++i) {
    //for (j = 0; j < max_hzn; ++j) {
    //  larger[i][j] = 0;
    //}
    larger[i][i] = true;  // result of comparison with itself
    //ri = rank[i%4][i/4]; // first index loops zone, second loops candidate. Zone loops faster, so we give equal priority to zones
    ri = rank[i];

    for (j = i + 1; j < max_hzn; ++j) {
      // i&j bits show which rank is larger
      //rj = rank[j%4][j/4];
      rj = rank[j];
      if (ri >= rj)
        larger[i][j] = true;
      else
        larger[j][i] = true;
    }
    // "larger" array shows the result of comparison for each rank

    // track exists if quality != 0
    exists[i] = (ri != 0);
  }

  // ghost cancellation, over 3 BXs
  for (i = 0; i < max_hzn - 1; ++i) {    // candidate loop
    for (j = i + 1; j < max_hzn; ++j) {  // comparison candidate loop
      int shared_segs = 0;

      // count shared segments
      for (const auto& isegment : segments.at(i)) {  // loop over all pairs of hits
        for (const auto& jsegment : segments.at(j)) {
          if (isegment == jsegment) {  // same chamber and same segment
            shared_segs += 1;
          }
        }
      }

      if (shared_segs > 0) {  // a single shared segment means it's ghost
        // kill candidate that has lower rank
        if (larger[i][j])
          killed[j] = true;
        else
          killed[i] = true;
      }
    }
  }

  // remove ghosts according to kill mask
  //exists = exists & (~kill1);
  for (i = 0; i < max_hzn; ++i) {
    exists[i] = exists[i] & (!killed[i]);
  }

  // remove tracks that are not at correct BX number
  //exists = exists & good_bx;
  for (i = 0; i < max_hzn; ++i) {
    exists[i] = exists[i] & good_bx[i];
  }

  bool anything_exists = (std::find(exists.begin(), exists.end(), 1) != exists.end());
  if (!anything_exists)
    return;

  // update "larger" array
  for (i = 0; i < max_hzn; ++i) {
    for (j = 0; j < max_hzn; ++j) {
      //if  (exists[i]) larger[i] = larger[i] | (~exists); // if this track exists make it larger than all non-existing tracks
      //else  larger[i] = 0; // else make it smaller than anything
      if (exists[i])
        larger[i][j] = larger[i][j] | (!exists[j]);
      else
        larger[i][j] = false;
    }
  }

  if (verbose_ > 0) {  // debug
    std::cout << "exists: ";
    for (i = max_hzn - 1; i >= 0; --i) {
      std::cout << exists[i];
      if ((i % max_zn) == 0 && i != 0)
        std::cout << "_";
    }
    std::cout << std::endl;
    std::cout << "killed: ";
    for (i = max_hzn - 1; i >= 0; --i) {
      std::cout << killed[i];
      if ((i % max_zn) == 0 && i != 0)
        std::cout << "_";
    }
    std::cout << std::endl;
    for (j = 0; j < max_hzn; ++j) {
      std::cout << "larger: ";
      for (i = max_hzn - 1; i >= 0; --i) {
        std::cout << larger[j][i];
        if ((i % max_zn) == 0 && i != 0)
          std::cout << "_";
      }
      std::cout << std::endl;
    }
  }

  // count zeros in the comparison results. The best track will have none, the next will have one, the third will have two
  for (i = 0; i < max_hzn; ++i) {
    sum = 0;
    for (j = 0; j < max_hzn; ++j) {
      if (larger[i][j] == 0)
        sum += 1;
    }

    if (sum < maxTracks_) {
      winner[sum][i] = true;  // assign positional winner codes
    }

    if (sum < maxTracks_ && bugSameSectorPt0_ && sum > 0) {
      // just keep the best track and kill the rest of them
      winner[sum][i] = false;
    }
  }

  // Output best tracks according to winner signals
  best_tracks.clear();

  for (int o = 0; o < maxTracks_; ++o) {  // output candidate loop
    int h = 0, n = 0, z = 0;
    for (i = 0; i < max_hzn; ++i) {  // winner bit loop
      if (winner[o][i]) {
        h = (i / max_z / max_n);
        n = (i / max_z) % max_n;
        z = i % max_z;

        const EMTFTrackCollection& tracks = extended_best_track_cands.at(h * max_z + z);
        const EMTFTrack& track = tracks.at(n);
        best_tracks.push_back(track);

        // Update winner, BX
        best_tracks.back().set_track_num(best_tracks.size() - 1);
        best_tracks.back().set_winner(o);
        best_tracks.back().set_bx(best_tracks.back().Second_BX());
      }
    }
  }
}