CMS 3D CMS Logo

BestTrackSelection.cc
Go to the documentation of this file.
2 
3 #include "helper.h" // to_hex, to_binary
4 
6  int endcap,
7  int sector,
8  int bx,
9  int bxWindow,
10  int maxRoadsPerZone,
11  int maxTracks,
12  bool useSecondEarliest,
13  bool bugSameSectorPt0) {
14  verbose_ = verbose;
15  endcap_ = endcap;
16  sector_ = sector;
17  bx_ = bx;
18 
19  bxWindow_ = bxWindow;
20  maxRoadsPerZone_ = maxRoadsPerZone;
22  useSecondEarliest_ = useSecondEarliest;
23  bugSameSectorPt0_ = bugSameSectorPt0;
24 }
25 
26 void BestTrackSelection::process(const std::deque<EMTFTrackCollection>& extended_best_track_cands,
27  EMTFTrackCollection& best_tracks) const {
28  int num_cands = 0;
29  for (const auto& cands : extended_best_track_cands) {
30  for (const auto& cand : cands) {
31  if (cand.Rank() > 0) {
32  num_cands += 1;
33  }
34  }
35  }
36  bool early_exit = (num_cands == 0);
37 
38  if (early_exit)
39  return;
40 
41  if (!useSecondEarliest_) {
42  cancel_one_bx(extended_best_track_cands, best_tracks);
43  } else {
44  cancel_multi_bx(extended_best_track_cands, best_tracks);
45  }
46 
47  if (verbose_ > 0) { // debug
48  for (const auto& track : best_tracks) {
49  std::cout << "track: " << track.Winner() << " rank: " << to_hex(track.Rank())
50  << " ph_deltas: " << array_as_string(track.PtLUT().delta_ph)
51  << " th_deltas: " << array_as_string(track.PtLUT().delta_th) << " phi: " << track.Phi_fp()
52  << " theta: " << track.Theta_fp() << " cpat: " << array_as_string(track.PtLUT().cpattern)
53  << " bx: " << track.BX() << std::endl;
54  for (int i = 0; i < emtf::NUM_STATIONS + 1; ++i) { // stations 0-4
55  if (track.PtLUT().bt_vi[i] != 0)
56  std::cout << ".. track segments: st: " << i << " v: " << track.PtLUT().bt_vi[i]
57  << " h: " << track.PtLUT().bt_hi[i] << " c: " << track.PtLUT().bt_ci[i]
58  << " s: " << track.PtLUT().bt_si[i] << std::endl;
59  }
60  }
61  }
62 }
63 
64 void BestTrackSelection::cancel_one_bx(const std::deque<EMTFTrackCollection>& extended_best_track_cands,
65  EMTFTrackCollection& best_tracks) const {
66  const int max_z = emtf::NUM_ZONES; // = 4 zones
67  const int max_n = maxRoadsPerZone_; // = 3 candidates per zone
68  const int max_zn = max_z * max_n; // = 12 total candidates
69  emtf_assert(maxTracks_ <= max_zn);
70 
71  // Emulate the arrays used in firmware
72  typedef std::array<int, 3> segment_ref_t;
73  std::vector<std::vector<segment_ref_t> > segments(max_zn,
74  std::vector<segment_ref_t>()); // 2D array [zn][num segments]
75 
76  std::vector<std::vector<bool> > larger(max_zn, std::vector<bool>(max_zn, false)); // 2D array [zn][zn]
77  std::vector<std::vector<bool> > winner(max_zn, std::vector<bool>(max_zn, false));
78 
79  std::vector<bool> exists(max_zn, false); // 1D array [zn]
80  std::vector<bool> killed(max_zn, false);
81  std::vector<int> rank(max_zn, 0);
82  //std::vector<int> good_bx(max_zn, 0);
83 
84  // Initialize arrays: rank, segments
85  for (int z = 0; z < max_z; ++z) {
86  const EMTFTrackCollection& tracks = extended_best_track_cands.at(z);
87  const int ntracks = tracks.size();
88  emtf_assert(ntracks <= max_n);
89 
90  for (int n = 0; n < ntracks; ++n) {
91  const int zn = (n * max_z) + z; // for (i = 0; i < 12; i = i+1) rank[i%4][i/4]
92  const EMTFTrack& track = tracks.at(n);
93 
94  rank.at(zn) = track.Rank();
95 
96  for (const auto& conv_hit : track.Hits()) {
97  emtf_assert(conv_hit.Valid());
98 
99  // A segment identifier (chamber, strip, bx)
100  const segment_ref_t segment = {{conv_hit.PC_station() * 9 + conv_hit.PC_chamber(),
101  (conv_hit.PC_segment() % 2),
102  0}}; // FW doesn't check whether a segment is CSC or RPC
103  //int chamber_index = int(conv_hit.Subsystem() == TriggerPrimitive::kRPC)*9*5;
104  //chamber_index += conv_hit.PC_station()*9;
105  //chamber_index += conv_hit.PC_chamber();
106  //const segment_ref_t segment = {{chamber_index, conv_hit.Strip(), 0}};
107  segments.at(zn).push_back(segment);
108  }
109  } // end loop over n
110  } // end loop over z
111 
112  // Simultaneously compare each rank with each other
113  int i = 0, j = 0, ri = 0, rj = 0, gt = 0, eq = 0, sum = 0;
114 
115  for (i = 0; i < max_zn; ++i) {
116  for (j = 0; j < max_zn; ++j) {
117  larger[i][j] = false;
118  }
119  larger[i][i] = true; // result of comparison with itself
120  //ri = rank[i%4][i/4]; // first index loops zone, second loops candidate. Zone loops faster, so we give equal priority to zones
121  ri = rank[i];
122 
123  for (j = 0; j < max_zn; ++j) {
124  // i&j bits show which rank is larger
125  // the comparison scheme below avoids problems
126  // when there are two or more tracks with the same rank
127  //rj = rank[j%4][j/4];
128  rj = rank[j];
129  gt = ri > rj;
130  eq = ri == rj;
131  if ((i < j && (gt || eq)) || (i > j && gt))
132  larger[i][j] = true;
133  }
134  // "larger" array shows the result of comparison for each rank
135 
136  // track exists if quality != 0
137  exists[i] = (ri != 0);
138  }
139 
140  // ghost cancellation, only in the current BX so far
141  for (i = 0; i < max_zn - 1; ++i) { // candidate loop
142  for (j = i + 1; j < max_zn; ++j) { // comparison candidate loop
143  int shared_segs = 0;
144 
145  // count shared segments
146  for (const auto& isegment : segments.at(i)) { // loop over all pairs of hits
147  for (const auto& jsegment : segments.at(j)) {
148  if (isegment == jsegment) { // same chamber and same segment
149  shared_segs += 1;
150  }
151  }
152  }
153 
154  if (shared_segs > 0) { // a single shared segment means it's ghost
155  // kill candidate that has lower rank
156  if (larger[i][j])
157  killed[j] = true;
158  else
159  killed[i] = true;
160  }
161  }
162  }
163 
164  // remove ghosts according to kill mask
165  //exists = exists & (~kill1);
166  for (i = 0; i < max_zn; ++i) {
167  exists[i] = exists[i] & (!killed[i]);
168  }
169 
170  bool anything_exists = (std::find(exists.begin(), exists.end(), 1) != exists.end());
171  if (!anything_exists)
172  return;
173 
174  // update "larger" array
175  for (i = 0; i < max_zn; ++i) {
176  for (j = 0; j < max_zn; ++j) {
177  //if (exists[i]) larger[i] = larger[i] | (~exists); // if this track exists make it larger than all non-existing tracks
178  //else larger[i] = 0; // else make it smaller than anything
179  if (exists[i])
180  larger[i][j] = larger[i][j] | (!exists[j]);
181  else
182  larger[i][j] = false;
183  }
184  }
185 
186  if (verbose_ > 0) { // debug
187  std::cout << "exists: ";
188  for (i = max_zn - 1; i >= 0; --i) {
189  std::cout << exists[i];
190  }
191  std::cout << std::endl;
192  std::cout << "killed: ";
193  for (i = max_zn - 1; i >= 0; --i) {
194  std::cout << killed[i];
195  }
196  std::cout << std::endl;
197  for (j = 0; j < max_zn; ++j) {
198  std::cout << "larger: ";
199  for (i = max_zn - 1; i >= 0; --i) {
200  std::cout << larger[j][i];
201  }
202  std::cout << std::endl;
203  }
204  }
205 
206  // count zeros in the comparison results. The best track will have none, the next will have one, the third will have two
207  // skip the bits corresponding to the comparison of the track with itself
208  for (i = 0; i < max_zn; ++i) {
209  sum = 0;
210  for (j = 0; j < max_zn; ++j) {
211  if (larger[i][j] == 0)
212  sum += 1;
213  }
214 
215  if (sum < maxTracks_) {
216  winner[sum][i] = true; // assign positional winner codes
217  }
218 
219  if (sum < maxTracks_ && bugSameSectorPt0_ && sum > 0) {
220  // just keep the best track and kill the rest of them
221  winner[sum][i] = false;
222  }
223  }
224 
225  // Output best tracks according to winner signals
226  best_tracks.clear();
227 
228  for (int o = 0; o < maxTracks_; ++o) { // output candidate loop
229  int z = 0, n = 0;
230  for (i = 0; i < max_zn; ++i) { // winner bit loop
231  if (winner[o][i]) {
232  n = i / max_z;
233  z = i % max_z;
234 
235  const EMTFTrackCollection& tracks = extended_best_track_cands.at(z);
236  const EMTFTrack& track = tracks.at(n);
237  best_tracks.push_back(track);
238 
239  // Update winner, BX
240  best_tracks.back().set_track_num(best_tracks.size() - 1);
241  best_tracks.back().set_winner(o);
242  best_tracks.back().set_bx(best_tracks.back().First_BX());
243  }
244  }
245  }
246 }
247 
248 void BestTrackSelection::cancel_multi_bx(const std::deque<EMTFTrackCollection>& extended_best_track_cands,
249  EMTFTrackCollection& best_tracks) const {
250  const int max_h = bxWindow_; // = 3 bx history
251  const int max_z = emtf::NUM_ZONES; // = 4 zones
252  const int max_n = maxRoadsPerZone_; // = 3 candidates per zone
253  const int max_zn = max_z * max_n; // = 12 total candidates
254  const int max_hzn = max_h * max_zn; // = 36 total candidates
255  emtf_assert(maxTracks_ <= max_hzn);
256 
257  const int delayBX = bxWindow_ - 1;
258  const int num_h = extended_best_track_cands.size() / max_z; // num of bx history so far
259 
260  // Emulate the arrays used in firmware
261  typedef std::array<int, 3> segment_ref_t;
262  std::vector<std::vector<segment_ref_t> > segments(max_hzn,
263  std::vector<segment_ref_t>()); // 2D array [hzn][num segments]
264 
265  std::vector<std::vector<bool> > larger(max_hzn, std::vector<bool>(max_hzn, false)); // 2D array [hzn][hzn]
266  std::vector<std::vector<bool> > winner(max_hzn, std::vector<bool>(max_hzn, false));
267 
268  std::vector<bool> exists(max_hzn, false); // 1D array [hzn]
269  std::vector<bool> killed(max_hzn, false);
270  std::vector<int> rank(max_hzn, 0);
271  std::vector<int> good_bx(max_hzn, 0);
272 
273  // Initialize arrays: rank, good_bx, segments
274  for (int h = 0; h < num_h; ++h) {
275  // 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]
276  // extended_best_track_cands[4..7] has 4 zones for road/pattern BX-1 with possible tracks from [BX-3, BX-2, BX-1]
277  // extended_best_track_cands[8..11] has 4 zones for road/pattern BX-2 with possible tracks from [BX-4, BX-3, BX-2]
278 
279  for (int z = 0; z < max_z; ++z) {
280  const EMTFTrackCollection& tracks = extended_best_track_cands.at(h * max_z + z);
281  const int ntracks = tracks.size();
282  emtf_assert(ntracks <= max_n);
283 
284  for (int n = 0; n < ntracks; ++n) {
285  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]
286  const EMTFTrack& track = tracks.at(n);
287  int cand_bx = track.Second_BX();
288  cand_bx -= (bx_ - delayBX); // convert track.second_bx=[BX-2, BX-1, BX-0] --> cand_bx=[0,1,2]
289 
290  rank.at(hzn) = track.Rank();
291  if (cand_bx == 0)
292  good_bx.at(hzn) = 1; // kill this rank if it's not the right BX
293 
294  for (const auto& conv_hit : track.Hits()) {
295  emtf_assert(conv_hit.Valid());
296 
297  // Notes from Alex (2017-03-16):
298  //
299  // What happens currently is that RPC hits are inserted instead of
300  // CSC hits, if CSC hits are missing. So the track address will point
301  // at a CSC chamber, but the actual hit may have come from
302  // corresponding RPC located behind it. If the same substitution
303  // happened in the neighboring sector, then the cancellation in uGMT
304  // will work correctly. If the substitution happened in only one
305  // sector, then RPC hit may "trump" a CSC hit, or vice versa.
306 
307  // A segment identifier (chamber, strip, bx)
308  const segment_ref_t segment = {{conv_hit.PC_station() * 9 + conv_hit.PC_chamber(),
309  (conv_hit.PC_segment() % 2),
310  conv_hit.BX()}}; // FW doesn't check whether a segment is CSC or RPC
311  //int chamber_index = int(conv_hit.Subsystem() == TriggerPrimitive::kRPC)*9*5;
312  //chamber_index += conv_hit.PC_station()*9;
313  //chamber_index += conv_hit.PC_chamber();
314  //const segment_ref_t segment = {{chamber_index, conv_hit.Strip(), conv_hit.BX()}};
315  segments.at(hzn).push_back(segment);
316  }
317  } // end loop over n
318  } // end loop over z
319  } // end loop over h
320 
321  // Simultaneously compare each rank with each other
322  int i = 0, j = 0, ri = 0, rj = 0, sum = 0;
323 
324  for (i = 0; i < max_hzn; ++i) {
325  //for (j = 0; j < max_hzn; ++j) {
326  // larger[i][j] = 0;
327  //}
328  larger[i][i] = true; // result of comparison with itself
329  //ri = rank[i%4][i/4]; // first index loops zone, second loops candidate. Zone loops faster, so we give equal priority to zones
330  ri = rank[i];
331 
332  for (j = i + 1; j < max_hzn; ++j) {
333  // i&j bits show which rank is larger
334  //rj = rank[j%4][j/4];
335  rj = rank[j];
336  if (ri >= rj)
337  larger[i][j] = true;
338  else
339  larger[j][i] = true;
340  }
341  // "larger" array shows the result of comparison for each rank
342 
343  // track exists if quality != 0
344  exists[i] = (ri != 0);
345  }
346 
347  // ghost cancellation, over 3 BXs
348  for (i = 0; i < max_hzn - 1; ++i) { // candidate loop
349  for (j = i + 1; j < max_hzn; ++j) { // comparison candidate loop
350  int shared_segs = 0;
351 
352  // count shared segments
353  for (const auto& isegment : segments.at(i)) { // loop over all pairs of hits
354  for (const auto& jsegment : segments.at(j)) {
355  if (isegment == jsegment) { // same chamber and same segment
356  shared_segs += 1;
357  }
358  }
359  }
360 
361  if (shared_segs > 0) { // a single shared segment means it's ghost
362  // kill candidate that has lower rank
363  if (larger[i][j])
364  killed[j] = true;
365  else
366  killed[i] = true;
367  }
368  }
369  }
370 
371  // remove ghosts according to kill mask
372  //exists = exists & (~kill1);
373  for (i = 0; i < max_hzn; ++i) {
374  exists[i] = exists[i] & (!killed[i]);
375  }
376 
377  // remove tracks that are not at correct BX number
378  //exists = exists & good_bx;
379  for (i = 0; i < max_hzn; ++i) {
380  exists[i] = exists[i] & good_bx[i];
381  }
382 
383  bool anything_exists = (std::find(exists.begin(), exists.end(), 1) != exists.end());
384  if (!anything_exists)
385  return;
386 
387  // update "larger" array
388  for (i = 0; i < max_hzn; ++i) {
389  for (j = 0; j < max_hzn; ++j) {
390  //if (exists[i]) larger[i] = larger[i] | (~exists); // if this track exists make it larger than all non-existing tracks
391  //else larger[i] = 0; // else make it smaller than anything
392  if (exists[i])
393  larger[i][j] = larger[i][j] | (!exists[j]);
394  else
395  larger[i][j] = false;
396  }
397  }
398 
399  if (verbose_ > 0) { // debug
400  std::cout << "exists: ";
401  for (i = max_hzn - 1; i >= 0; --i) {
402  std::cout << exists[i];
403  if ((i % max_zn) == 0 && i != 0)
404  std::cout << "_";
405  }
406  std::cout << std::endl;
407  std::cout << "killed: ";
408  for (i = max_hzn - 1; i >= 0; --i) {
409  std::cout << killed[i];
410  if ((i % max_zn) == 0 && i != 0)
411  std::cout << "_";
412  }
413  std::cout << std::endl;
414  for (j = 0; j < max_hzn; ++j) {
415  std::cout << "larger: ";
416  for (i = max_hzn - 1; i >= 0; --i) {
417  std::cout << larger[j][i];
418  if ((i % max_zn) == 0 && i != 0)
419  std::cout << "_";
420  }
421  std::cout << std::endl;
422  }
423  }
424 
425  // count zeros in the comparison results. The best track will have none, the next will have one, the third will have two
426  for (i = 0; i < max_hzn; ++i) {
427  sum = 0;
428  for (j = 0; j < max_hzn; ++j) {
429  if (larger[i][j] == 0)
430  sum += 1;
431  }
432 
433  if (sum < maxTracks_) {
434  winner[sum][i] = true; // assign positional winner codes
435  }
436 
437  if (sum < maxTracks_ && bugSameSectorPt0_ && sum > 0) {
438  // just keep the best track and kill the rest of them
439  winner[sum][i] = false;
440  }
441  }
442 
443  // Output best tracks according to winner signals
444  best_tracks.clear();
445 
446  for (int o = 0; o < maxTracks_; ++o) { // output candidate loop
447  int h = 0, n = 0, z = 0;
448  for (i = 0; i < max_hzn; ++i) { // winner bit loop
449  if (winner[o][i]) {
450  h = (i / max_z / max_n);
451  n = (i / max_z) % max_n;
452  z = i % max_z;
453 
454  const EMTFTrackCollection& tracks = extended_best_track_cands.at(h * max_z + z);
455  const EMTFTrack& track = tracks.at(n);
456  best_tracks.push_back(track);
457 
458  // Update winner, BX
459  best_tracks.back().set_track_num(best_tracks.size() - 1);
460  best_tracks.back().set_winner(o);
461  best_tracks.back().set_bx(best_tracks.back().Second_BX());
462  }
463  }
464  }
465 }
static char to_hex(unsigned int i)
Definition: types.cc:27
void process(const std::deque< EMTFTrackCollection > &extended_best_track_cands, EMTFTrackCollection &best_tracks) const
l1t::EMTFTrackCollection EMTFTrackCollection
Definition: Common.h:27
constexpr int NUM_ZONES
Definition: Common.h:54
bool verbose
int killed
Definition: mps_check.py:40
void find(edm::Handle< EcalRecHitCollection > &hits, DetId thisDet, std::vector< EcalRecHitCollection::const_iterator > &hit, bool debug=false)
Definition: FindCaloHit.cc:19
void configure(int verbose, int endcap, int sector, int bx, int bxWindow, int maxRoadsPerZone, int maxTracks, bool useSecondEarliest, bool bugSameSectorPt0)
auto const & tracks
cannot be loose
#define emtf_assert(expr)
Definition: DebugTools.h:18
constexpr int NUM_STATIONS
Definition: Common.h:58
void cancel_one_bx(const std::deque< EMTFTrackCollection > &extended_best_track_cands, EMTFTrackCollection &best_tracks) const
The Signals That Services Can Subscribe To This is based on ActivityRegistry h
Helper function to determine trigger accepts.
Definition: Activities.doc:4
void cancel_multi_bx(const std::deque< EMTFTrackCollection > &extended_best_track_cands, EMTFTrackCollection &best_tracks) const