Belle II Software  release-05-01-25
CDCTriggerHoughtrafo.cc
1 /*************************************************************************
2 * BASF2 (Belle Analysis Framework 2) *
3 * Copyright(C) 2013-2014 Belle II Collaboration *
4 * *
5 * Author: The Belle II Collaboration *
6 * Contributors: Michael Schnell, Sara Neuhaus *
7 * *
8 **************************************************************************/
9 
10 #include <trg/cdc/modules/houghtracking/CDCTrigger2DFinderModule.h>
11 
12 #include <framework/logging/Logger.h>
13 
14 #include <cmath>
15 
16 #include <root/TMatrix.h>
17 
18 /* defines */
19 #define CDC_SUPER_LAYERS 9
20 
21 using namespace std;
22 using namespace Belle2;
23 
24 unsigned short
25 CDCTrigger2DFinderModule::countSL(bool* layer)
26 {
27  if (m_requireSL0 && !layer[0]) return 0;
28  unsigned short lcnt = 0;
29  for (int i = 0; i < CDC_SUPER_LAYERS; ++i) {
30  if (layer[i] == true) ++lcnt;
31  }
32  return lcnt;
33 }
34 
35 bool
36 CDCTrigger2DFinderModule::shortTrack(bool* layer)
37 {
38  unsigned short lcnt = 0;
39  // check axial super layers (even layer number),
40  // break at first layer without hit
41  for (int i = 0; i < CDC_SUPER_LAYERS; i += 2) {
42  if (layer[i] == true) ++lcnt;
43  else break;
44  }
45  return (lcnt >= m_minHitsShort);
46 }
47 
48 /*
49 * Find the intercept in hough space.
50 * In each iteration the Hough plane is divided in quarters.
51 * For each quarter containing enough hits, the process is repeated
52 * until the quarters correspond to single Hough cells.
53 * Return zero on success.
54 * x = m
55 * y = a
56 */
57 int
58 CDCTrigger2DFinderModule::fastInterceptFinder(cdcMap& hits,
59  double x1_s, double x2_s,
60  double y1_s, double y2_s,
61  unsigned iterations,
62  unsigned ix_s, unsigned iy_s)
63 {
64  string indent = "";
65  for (unsigned i = 0; i < iterations; ++i) indent += " ";
66  B2DEBUG(150, indent << "intercept finder iteration " << iterations
67  << " x1 " << x1_s * 180 / M_PI << " x2 " << x2_s * 180 / M_PI
68  << " y1 " << y1_s << " y2 " << y2_s
69  << " ix " << ix_s << " iy " << iy_s);
70 
71  int i, j, iHit;
72  double unitx, unity;
73  double y1, y2;
74  double m, a;
75  cdcPair hp;
76  unsigned short iSL;
77  // list of hit indices crossing the current rectangle
78  vector<unsigned> idx_list;
79 
80  // limits of divided rectangles
81  double x1_d, x2_d, y1_d, y2_d;
82 
83  // cell indices in full plane
84  unsigned ix, iy;
85 
86  unitx = ((x2_s - x1_s) / 2.0);
87  unity = ((y2_s - y1_s) / 2.0);
88 
89  // divide x axis in half
90  for (i = 0; i < 2 ; ++i) {
91  // divide y axis in half
92  for (j = 0; j < 2; ++j) {
93  x1_d = x1_s + i * unitx;
94  x2_d = x1_s + (i + 1) * unitx;
95  y1_d = y1_s + j * unity;
96  y2_d = y1_s + (j + 1) * unity;
97 
98  ix = 2 * ix_s + i;
99  iy = 2 * iy_s + j;
100 
101  // skip extra cells outside of the Hough plane
102  if (x2_d <= -M_PI || x1_d >= M_PI || y2_d <= -maxR + shiftR || y1_d >= maxR + shiftR) {
103  B2DEBUG(150, indent << "skip Hough cell outside of plane limits");
104  continue;
105  }
106 
107  // if the cell size in phi is too large, the hit calculation is not reliable
108  // -> continue to next iteration without hit check
109  if (iterations != maxIterations && unitx > M_PI / 2.) {
110  fastInterceptFinder(hits, x1_d, x2_d, y1_d, y2_d, iterations + 1, ix, iy);
111  continue;
112  }
113 
114  idx_list.clear();
115  bool layerHit[CDC_SUPER_LAYERS] = {false}; /* For layer filter */
116  for (auto it = hits.begin(); it != hits.end(); ++it) {
117  iHit = it->first;
118  hp = it->second;
119  iSL = hp.first;
120  m = hp.second.X();
121  a = hp.second.Y();
122  // calculate Hough curve with slightly enlarged limits to avoid errors due to rounding
123  y1 = m * sin(x1_d - 1e-10) - a * cos(x1_d - 1e-10);
124  y2 = m * sin(x2_d + 1e-10) - a * cos(x2_d + 1e-10);
125  // skip decreasing half of the sine (corresponds to curl back half of the track)
126  if (iterations == maxIterations && y1 > y2) continue;
127  if (!((y1 > y2_d && y2 > y2_d) || (y1 < y1_d && y2 < y1_d))) {
128  if (iterations == maxIterations) {
129  idx_list.push_back(iHit);
130  }
131  layerHit[iSL] = true; /* layer filter */
132  }
133  }
134  unsigned short nSL = countSL(layerHit);
135  B2DEBUG(150, indent << "i " << i << " j " << j
136  << " layerHit " << int(layerHit[0]) << int(layerHit[2])
137  << int(layerHit[4]) << int(layerHit[6]) << int(layerHit[8])
138  << " nSL " << nSL);
139  if (nSL >= m_minHits || shortTrack(layerHit)) {
140  if (iterations != maxIterations) {
141  fastInterceptFinder(hits, x1_d, x2_d, y1_d, y2_d, iterations + 1, ix, iy);
142  } else {
143  TVector2 v1(x1_d, y1_d);
144  TVector2 v2(x2_d, y2_d);
145  houghCand.push_back(CDCTriggerHoughCand(idx_list, make_pair(v1, v2),
146  nSL, houghCand.size()));
147  if (m_storePlane > 0) {
148  (*m_houghPlane)[ix - (nCells - m_nCellsPhi) / 2][iy - (nCells - m_nCellsR) / 2] = nSL;
149  }
150  }
151  } else if (m_storePlane > 1) {
152  // to store the full Hough plane, we need to continue the iteration
153  // to minimal cell size everywhere
154  for (unsigned min = m_minHits - 1; min > 0; --min) {
155  if (nSL >= min) {
156  if (iterations != maxIterations) {
157  fastInterceptFinder(hits, x1_d, x2_d, y1_d, y2_d, iterations + 1, ix, iy);
158  } else {
159  (*m_houghPlane)[ix - (nCells - m_nCellsPhi) / 2][iy - (nCells - m_nCellsR) / 2] = nSL;
160  }
161  break;
162  }
163  }
164  }
165  }
166  }
167 
168  return 0;
169 }
170 
171 /*
172 * Peak finding method: connected regions & center of gravity
173 */
174 void
175 CDCTrigger2DFinderModule::connectedRegions()
176 {
177  vector<vector<CDCTriggerHoughCand>> regions;
178  vector<CDCTriggerHoughCand> cpyCand = houghCand;
179 
180  // debug: print candidate list
181  B2DEBUG(50, "houghCand number " << cpyCand.size());
182  for (unsigned icand = 0; icand < houghCand.size(); ++icand) {
183  coord2dPair coord = houghCand[icand].getCoord();
184  B2DEBUG(100, houghCand[icand].getID()
185  << " nSL " << houghCand[icand].getSLcount()
186  << " x1 " << coord.first.X() << " x2 " << coord.second.X()
187  << " y1 " << coord.first.Y() << " y2 " << coord.second.Y());
188  }
189 
190  // combine connected cells to regions
191  while (cpyCand.size() > 0) {
192  B2DEBUG(100, "make new region");
193  vector<CDCTriggerHoughCand> region;
194  vector<CDCTriggerHoughCand> rejected;
195  CDCTriggerHoughCand start = cpyCand[0];
196  cpyCand.erase(cpyCand.begin());
197  region.push_back(start);
198  addNeighbors(start, cpyCand, region, rejected, start.getSLcount());
199  regions.push_back(region);
200  for (auto cand = cpyCand.begin(); cand != cpyCand.end();) {
201  if (inList(*cand, region) || inList(*cand, rejected))
202  cpyCand.erase(cand);
203  else
204  ++cand;
205  }
206  }
207 
208  // find center of gravity for each region
209  for (unsigned ir = 0; ir < regions.size(); ++ir) {
210  B2DEBUG(50, "region " << ir << " (" << regions[ir].size() << " cells).");
211  // skip regions with size below cut
212  if (regions[ir].size() < m_minCells) {
213  B2DEBUG(50, "Skipping region with " << regions[ir].size() << " cells.");
214  continue;
215  }
216  double xfirst = regions[ir][0].getCoord().first.X();
217  double x = 0;
218  double y = 0;
219  int n = 0;
220  vector<unsigned> mergedList;
221  int xmin = m_nCellsPhi;
222  int xmax = 0;
223  int ymin = m_nCellsR;
224  int ymax = 0;
225  vector<TVector2> cellIds = {};
226  for (unsigned ir2 = 0; ir2 < regions[ir].size(); ++ir2) {
227  coord2dPair hc = regions[ir][ir2].getCoord();
228  B2DEBUG(100, " " << regions[ir][ir2].getID()
229  << " nSL " << regions[ir][ir2].getSLcount()
230  << " x1 " << hc.first.X() << " x2 " << hc.second.X()
231  << " y1 " << hc.first.Y() << " y2 " << hc.second.Y());
232  int ix = floor((hc.first.X() + M_PI) / 2. / M_PI * m_nCellsPhi + 0.5);
233  int iy = floor((hc.first.Y() + maxR - shiftR) / 2. / maxR * m_nCellsR + 0.5);
234  x += (hc.first.X() + hc.second.X());
235  if (xfirst - hc.first.X() > M_PI) {
236  x += 4 * M_PI;
237  ix += m_nCellsPhi;
238  } else if (hc.first.X() - xfirst > M_PI) {
239  x -= 4 * M_PI;
240  ix -= m_nCellsPhi;
241  }
242  y += (hc.first.Y() + hc.second.Y());
243  n += 1;
244  vector<unsigned> idList = regions[ir][ir2].getIdList();
245  mergeIdList(mergedList, mergedList, idList);
246  xmin = min(xmin, ix);
247  xmax = max(xmax, ix);
248  ymin = min(ymin, iy);
249  ymax = max(ymax, iy);
250  cellIds.push_back(TVector2(ix, iy));
251  }
252  x *= 0.5 / n;
253  if (x > M_PI)
254  x -= 2 * M_PI;
255  else if (x <= -M_PI)
256  x += 2 * M_PI;
257  y *= 0.5 / n;
258  B2DEBUG(50, "x " << x << " y " << y);
259 
260  // select 1 hit per super layer
261  vector<unsigned> selectedList = {};
262  vector<unsigned> unselectedList = {};
263  selectHits(mergedList, selectedList, unselectedList);
264 
265  // save track
266  const CDCTriggerTrack* track =
267  m_tracks.appendNew(x, 2. * y, 0.);
268  // relations to selected hits
269  for (unsigned i = 0; i < selectedList.size(); ++i) {
270  unsigned its = selectedList[i];
271  track->addRelationTo(m_segmentHits[its]);
272  }
273  // relations to additional hits get a negative weight
274  for (unsigned i = 0; i < unselectedList.size(); ++i) {
275  unsigned its = unselectedList[i];
276  track->addRelationTo(m_segmentHits[its], -1.);
277  }
278  // save detail information about the cluster
279  const CDCTriggerHoughCluster* cluster =
280  m_clusters.appendNew(xmin, xmax, ymin, ymax, cellIds);
281  track->addRelationTo(cluster);
282  }
283 }
284 
285 void
286 CDCTrigger2DFinderModule::addNeighbors(const CDCTriggerHoughCand& center,
287  const vector<CDCTriggerHoughCand>& candidates,
288  vector<CDCTriggerHoughCand>& merged,
289  vector<CDCTriggerHoughCand>& rejected,
290  unsigned short nSLmax) const
291 {
292  for (unsigned icand = 0; icand < candidates.size(); ++icand) {
293  B2DEBUG(120, "compare center " << center.getID()
294  << " to " << candidates[icand].getID());
295  if (inList(candidates[icand], merged) || inList(candidates[icand], rejected)) {
296  B2DEBUG(120, " " << candidates[icand].getID() << " already in list");
297  continue;
298  }
299  bool reject = inList(center, rejected);
300  if (connected(center, candidates[icand])) {
301  if (m_onlyLocalMax && candidates[icand].getSLcount() < nSLmax) {
302  B2DEBUG(100, " lower than highest SLcount, rejected");
303  rejected.push_back(candidates[icand]);
304  } else if (m_onlyLocalMax && !reject && candidates[icand].getSLcount() > nSLmax) {
305  B2DEBUG(100, " new highest SLcount, clearing list");
306  nSLmax = candidates[icand].getSLcount();
307  for (unsigned imerged = 0; imerged < merged.size(); ++imerged) {
308  rejected.push_back(merged[imerged]);
309  }
310  merged.clear();
311  merged.push_back(candidates[icand]);
312  } else if (m_onlyLocalMax && candidates[icand].getSLcount() > center.getSLcount()) {
313  B2DEBUG(100, " connected to rejected cell, skip");
314  continue;
315  } else if (m_onlyLocalMax && reject) {
316  B2DEBUG(100, " connected to rejected cell, rejected");
317  rejected.push_back(candidates[icand]);
318  } else {
319  B2DEBUG(100, " connected");
320  merged.push_back(candidates[icand]);
321  }
322  vector<CDCTriggerHoughCand> cpyCand = candidates;
323  cpyCand.erase(cpyCand.begin() + icand);
324  addNeighbors(candidates[icand], cpyCand, merged, rejected, nSLmax);
325  }
326  }
327 }
328 
329 bool
330 CDCTrigger2DFinderModule::inList(const CDCTriggerHoughCand& a,
331  const vector<CDCTriggerHoughCand>& list) const
332 {
333  for (unsigned i = 0; i < list.size(); ++i) {
334  if (a == list[i]) return true;
335  }
336  return false;
337 }
338 
339 bool
340 CDCTrigger2DFinderModule::connected(const CDCTriggerHoughCand& a,
341  const CDCTriggerHoughCand& b) const
342 {
343  double ax1 = a.getCoord().first.X();
344  double ax2 = a.getCoord().second.X();
345  double ay1 = a.getCoord().first.Y();
346  double ay2 = a.getCoord().second.Y();
347  double bx1 = b.getCoord().first.X();
348  double bx2 = b.getCoord().second.X();
349  double by1 = b.getCoord().first.Y();
350  double by2 = b.getCoord().second.Y();
351  // direct neighbors
352  bool direct = ((ax2 == bx1 && ay1 == by1) || // right
353  (ax1 == bx2 && ay1 == by1) || // left
354  (ax1 == bx1 && ay2 == by1) || // above
355  (ax1 == bx1 && ay1 == by2) || // below
356  (ax1 + 2. * M_PI == bx2 && ay1 == by1) ||
357  (ax2 == bx1 + 2. * M_PI && ay1 == by1));
358  // diagonal connections
359  bool diagRise = ((ax2 == bx1 && ay2 == by1) || // right above
360  (ax1 == bx2 && ay1 == by2) || // left below
361  (ax1 + 2. * M_PI == bx2 && ay1 == by2) ||
362  (ax2 == bx1 + 2. * M_PI && ay2 == by1));
363  bool diagFall = ((ax1 == bx2 && ay2 == by1) || // left above
364  (ax2 == bx1 && ay1 == by2) || // right below
365  (ax2 == bx1 + 2. * M_PI && ay1 == by2) ||
366  (ax1 + 2. * M_PI == bx2 && ay2 == by1));
367  if (m_connect == 4) return direct;
368  else if (m_connect == 6) return (direct || diagRise);
369  else if (m_connect == 8) return (direct || diagRise || diagFall);
370  else B2WARNING("Unknown option for connect " << m_connect << ", using default.");
371  return (direct || diagRise);
372 }
373 
374 /*
375 * Merge Id lists.
376 */
377 void
378 CDCTrigger2DFinderModule::mergeIdList(std::vector<unsigned>& merged,
379  std::vector<unsigned>& a,
380  std::vector<unsigned>& b)
381 {
382  bool found;
383 
384  for (auto it = a.begin(); it != a.end(); ++it) {
385  found = false;
386  for (auto it_in = merged.begin(); it_in != merged.end(); ++it_in) {
387  if (*it_in == *it) {
388  found = true;
389  break;
390  }
391  }
392  if (!found) {
393  merged.push_back(*it);
394  }
395  }
396 
397  for (auto it = b.begin(); it != b.end(); ++it) {
398  found = false;
399  for (auto it_in = merged.begin(); it_in != merged.end(); ++it_in) {
400  if (*it_in == *it) {
401  found = true;
402  break;
403  }
404  }
405  if (!found) {
406  merged.push_back(*it);
407  }
408  }
409 }
410 
411 void
412 CDCTrigger2DFinderModule::findAllCrossingHits(std::vector<unsigned>& list,
413  double x1, double x2,
414  double y1, double y2,
415  const cdcMap& inputMap)
416 {
417  double m, a, y1_h, y2_h;
418  for (int iHit = 0; iHit < m_segmentHits.getEntries(); iHit++) {
419  unsigned short iSL = m_segmentHits[iHit]->getISuperLayer();
420  if (iSL % 2) continue;
421  // associate only the hits actually used to form the cluster
422  if (m_suppressClone && inputMap.find(iHit) == inputMap.end()) continue;
423  //TODO: add options: center cell / active priority cell / all priority cells
424  vector<double> phi = {0, 0, 0};
425  phi[0] = m_segmentHits[iHit]->getSegmentID() - TSoffset[iSL];
426  phi[1] = phi[0] + 0.5;
427  phi[2] = phi[0] - 0.5;
428  vector<double> r = {radius[iSL][0], radius[iSL][1], radius[iSL][1]};
429  for (unsigned i = 0; i < 3; ++i) {
430  phi[i] *= 2. * M_PI / (TSoffset[iSL + 1] - TSoffset[iSL]);
431  m = cos(phi[i]) / r[i];
432  a = sin(phi[i]) / r[i];
433  // calculate Hough curve with slightly enlarged limits to avoid errors due to rounding
434  y1_h = m * sin(x1 - 1e-10) - a * cos(x1 - 1e-10);
435  y2_h = m * sin(x2 + 1e-10) - a * cos(x2 + 1e-10);
436  // skip decreasing half of the sine (corresponds to curl back half of the track)
437  if (y1_h > y2_h) continue;
438  if (!((y1_h > y2 && y2_h > y2) || (y1_h < y1 && y2_h < y1))) {
439  list.push_back(iHit);
440  break;
441  }
442  }
443  }
444 }
445 
446 /*
447  * Select one hit per super layer
448  */
449 void
450 CDCTrigger2DFinderModule::selectHits(std::vector<unsigned>& list,
451  std::vector<unsigned>& selected,
452  std::vector<unsigned>& unselected)
453 {
454  std::vector<int> bestPerSL(5, -1);
455  for (unsigned i = 0; i < list.size(); ++i) {
456  unsigned iax = m_segmentHits[list[i]]->getISuperLayer() / 2;
457  bool firstPriority = (m_segmentHits[list[i]]->getPriorityPosition() == 3);
458  if (bestPerSL[iax] < 0) {
459  bestPerSL[iax] = i;
460  } else {
461  unsigned itsBest = list[bestPerSL[iax]];
462  bool firstBest = (m_segmentHits[itsBest]->getPriorityPosition() == 3);
463  // selection rules:
464  // first priority, higher ID
465  if ((firstPriority && !firstBest) ||
466  (firstPriority == firstBest &&
467  m_segmentHits[list[i]]->getSegmentID() > m_segmentHits[itsBest]->getSegmentID())) {
468  bestPerSL[iax] = i;
469  }
470  }
471  }
472 
473  for (unsigned i = 0; i < list.size(); ++i) {
474  unsigned iax = m_segmentHits[list[i]]->getISuperLayer() / 2;
475  if (int(i) == bestPerSL[iax]) selected.push_back(list[i]);
476  else unselected.push_back(list[i]);
477  }
478 }
479 
480 /*
481 * Alternative peak finding with nested patterns
482 */
483 void
484 CDCTrigger2DFinderModule::patternClustering(const cdcMap& inputMap)
485 {
486  // fill a matrix of 2 x 2 squares
487  TMatrix plane2(m_nCellsPhi / 2, m_nCellsR / 2);
488  TMatrix planeIcand(m_nCellsPhi, m_nCellsR);
489  for (unsigned icand = 0; icand < houghCand.size(); ++icand) {
490  double x = (houghCand[icand].getCoord().first.X() +
491  houghCand[icand].getCoord().second.X()) / 2.;
492  unsigned ix = floor((x + M_PI) / 2. / M_PI * m_nCellsPhi);
493  double y = (houghCand[icand].getCoord().first.Y() +
494  houghCand[icand].getCoord().second.Y()) / 2.;
495  unsigned iy = floor((y + maxR - shiftR) / 2. / maxR * m_nCellsR);
496  plane2[ix / 2][iy / 2] += 1 << ((ix % 2) + 2 * (iy % 2));
497  planeIcand[ix][iy] = icand;
498  B2DEBUG(100, "candidate " << icand << " at ix " << ix << " iy " << iy);
499  }
500  // look for clusters of 2 x 2 squares in a (rX x rY) region
501  unsigned nX = m_nCellsPhi / 2;
502  unsigned nY = m_nCellsR / 2;
503  unsigned rX = m_clusterSizeX;
504  unsigned rY = m_clusterSizeY;
505  for (unsigned ix = 0; ix < nX; ++ix) {
506  for (unsigned iy = 0; iy < nY; ++iy) {
507  if (!plane2[ix][iy]) continue;
508  // check if we are in a lower left corner
509  unsigned ileft = (ix - 1 + nX) % nX;
510  B2DEBUG(100, "ix " << ix << " iy " << iy);
511  if (connectedLR(plane2[ileft][iy], plane2[ix][iy]) ||
512  (iy > 0 && connectedUD(plane2[ix][iy - 1], plane2[ix][iy])) ||
513  (iy > 0 && connectedDiag(plane2[ileft][iy - 1], plane2[ix][iy]))) {
514  B2DEBUG(100, "skip connected square");
515  continue;
516  }
517  // form cluster
518  vector<unsigned> pattern(rX * rY, 0);
519  pattern[0] = plane2[ix][iy];
520  vector<TVector2> cellIds = {TVector2(2 * ix, 2 * iy)};
521  for (unsigned ix2 = 0; ix2 < rX; ++ix2) {
522  for (unsigned iy2 = 0; iy2 < rY; ++iy2) {
523  if (iy + iy2 >= nY) continue;
524  unsigned ip = ix2 + rX * iy2;
525  unsigned iright = (ix + ix2 + nX) % nX;
526  ileft = (iright - 1 + nX) % nX;
527  B2DEBUG(100, "ix2 " << ix2 << " ileft " << ileft << " iright " << iright);
528  if ((ix2 > 0 && // check left/right connection
529  pattern[ip - 1] &&
530  connectedLR(plane2[ileft][iy + iy2], plane2[iright][iy + iy2])) ||
531  (iy2 > 0 && // check up/down connection
532  pattern[ip - rX] &&
533  connectedUD(plane2[iright][iy + iy2 - 1], plane2[iright][iy + iy2])) ||
534  (ix2 > 0 && iy2 > 0 && // check diagonal connection
535  pattern[ip - rX - 1] &&
536  connectedDiag(plane2[ileft][iy + iy2 - 1], plane2[iright][iy + iy2]))) {
537  pattern[ip] = plane2[iright][iy + iy2];
538  B2DEBUG(100, "connect cell " << iright << " " << iy + iy2);
539  cellIds.push_back(TVector2(2 * (ix + ix2), 2 * (iy + iy2)));
540  }
541  }
542  }
543  B2DEBUG(100, "cluster starting at " << ix << " " << iy);
544  // check if cluster continues beyond defined area
545  bool overflowRight = false;
546  bool overflowTop = false;
547  for (unsigned iy2 = 0; iy2 < rY; ++iy2) {
548  unsigned ip = rX - 1 + rX * iy2;
549  if (!pattern[ip]) continue;
550  unsigned iright = (ix + rX + nX) % nX;
551  ileft = (iright - 1 + nX) % nX;
552  if (connectedLR(plane2[ileft][iy + iy2], plane2[iright][iy + iy2]) ||
553  ((iy + iy2 + 1 < nY) &&
554  connectedDiag(plane2[ileft][iy + iy2], plane2[iright][iy + iy2 + 1]))) {
555  setReturnValue(false);
556  overflowRight = true;
557  }
558  }
559  if (iy + rY < nY) {
560  for (unsigned ix2 = 0; ix2 < rX; ++ix2) {
561  unsigned ip = ix2 + rX * (rY - 1);
562  if (!pattern[ip]) continue;
563  unsigned iright = (ix + ix2 + 1 + nX) % nX;
564  ileft = (iright - 1 + nX) % nX;
565  if (connectedUD(plane2[ileft][iy + rY - 1], plane2[ileft][iy + rY]) ||
566  connectedDiag(plane2[ileft][iy + rY - 1], plane2[iright][iy + rY])) {
567  setReturnValue(false);
568  overflowTop = true;
569  }
570  }
571  }
572  if (overflowRight && !overflowTop) {
573  B2DEBUG(100, "cluster extends right of " << rX << " x " << rY << " area");
574  } else if (overflowTop && !overflowRight) {
575  B2DEBUG(100, "cluster extends above " << rX << " x " << rY << " area");
576  } else if (overflowRight && overflowTop) {
577  B2DEBUG(100, "cluster extends right and above " << rX << " x " << rY << " area");
578  }
579  // find corners of cluster
580  unsigned topRight2 = topRightSquare(pattern);
581  unsigned topRight = topRightCorner(pattern[topRight2]);
582  unsigned bottomLeft = bottomLeftCorner(pattern[0]);
583  B2DEBUG(100, "topRight2 " << topRight2 << " topRight " << topRight << " bottomLeft " << bottomLeft);
584  // average over corners to find cluster center
585  unsigned ixTR = 2 * (topRight2 % m_clusterSizeX) + (topRight % 2);
586  unsigned ixBL = bottomLeft % 2;
587  unsigned iyTR = 2 * (topRight2 / m_clusterSizeX) + (topRight / 2);
588  unsigned iyBL = bottomLeft / 2;
589  B2DEBUG(100, "ixTR " << ixTR << " ixBL " << ixBL << " iyTR " << iyTR << " iyBL " << iyBL);
590  // skip size 1 clusters
591  if (m_minCells > 1 && ixTR == ixBL && iyTR == iyBL) {
592  B2DEBUG(100, "skipping cluster of size 1");
593  continue;
594  }
595  float centerX = 2 * ix + (ixTR + ixBL) / 2.;
596  if (centerX >= m_nCellsPhi) centerX -= m_nCellsPhi;
597  float centerY = 2 * iy + (iyTR + iyBL) / 2.;
598  B2DEBUG(100, "center at cell (" << centerX << ", " << centerY << ")");
599  // convert to coordinates
600  double x = -M_PI + (centerX + 0.5) * 2. * M_PI / m_nCellsPhi;
601  double y = -maxR + shiftR + (centerY + 0.5) * 2. * maxR / m_nCellsR;
602  B2DEBUG(100, "center coordinates (" << x << ", " << y << ")");
603  // get list of related hits
604  vector <unsigned> idList = {};
605  if (m_hitRelationsFromCorners) {
606  unsigned icandTR = planeIcand[(ixTR + 2 * ix) % m_nCellsPhi][iyTR + 2 * iy];
607  unsigned icandBL = planeIcand[ixBL + 2 * ix][iyBL + 2 * iy];
608  vector<unsigned> candIdListTR = houghCand[icandTR].getIdList();
609  vector<unsigned> candIdListBL = houghCand[icandBL].getIdList();
610  mergeIdList(idList, candIdListTR, candIdListBL);
611  B2DEBUG(100, "merge id lists from candidates " << icandTR << " and " << icandBL);
612  } else {
613  double dx = M_PI / m_nCellsPhi;
614  double dy = maxR / m_nCellsR;
615  double x1 = (round(centerX) == centerX) ? x - dx : x - 2 * dx;
616  double x2 = (round(centerX) == centerX) ? x + dx : x + 2 * dx;
617  double y1 = (round(centerY) == centerY) ? y - dy : y - 2 * dy;
618  double y2 = (round(centerY) == centerY) ? y + dy : y + 2 * dy;
619  findAllCrossingHits(idList, x1, x2, y1, y2, inputMap);
620  }
621  if (idList.size() == 0) {
622  setReturnValue(false);
623  B2DEBUG(100, "id list empty");
624  }
625 
626  // select 1 hit per super layer
627  vector<unsigned> selectedList = {};
628  vector<unsigned> unselectedList = {};
629  selectHits(idList, selectedList, unselectedList);
630 
631  // save track
632  const CDCTriggerTrack* track =
633  m_tracks.appendNew(x, 2. * y, 0.);
634  // relations to selected hits
635  for (unsigned i = 0; i < selectedList.size(); ++i) {
636  unsigned its = selectedList[i];
637  track->addRelationTo(m_segmentHits[its]);
638  }
639  // relations to additional hits get a negative weight
640  for (unsigned i = 0; i < unselectedList.size(); ++i) {
641  unsigned its = unselectedList[i];
642  track->addRelationTo(m_segmentHits[its], -1.);
643  }
644  // save detail information about the cluster
645  const CDCTriggerHoughCluster* cluster =
646  m_clusters.appendNew(2 * ix, 2 * (ix + m_clusterSizeX) - 1,
647  2 * iy, 2 * (iy + m_clusterSizeY) - 1,
648  cellIds);
649  track->addRelationTo(cluster);
650  }
651  }
652 }
653 
654 /*
655 * connection definitions for 2 x 2 squares
656 */
657 bool
658 CDCTrigger2DFinderModule::connectedLR(unsigned patternL, unsigned patternR)
659 {
660  // connected if
661  // . x | x . or . . | . .
662  // . . | . . . x | x .
663  bool connectDirect = (((patternL >> 3) & 1) && ((patternR >> 2) & 1)) ||
664  (((patternL >> 1) & 1) && ((patternR >> 0) & 1));
665  // connected if
666  // . . | x .
667  // . x | . .
668  bool connectRise = ((patternL >> 1) & 1) && ((patternR >> 2) & 1);
669  // connected if
670  // . x | . .
671  // . . | x .
672  bool connectFall = ((patternL >> 3) & 1) && ((patternR >> 0) & 1);
673 
674  if (m_connect == 4) return connectDirect;
675  else if (m_connect == 6) return (connectDirect || connectRise);
676  else if (m_connect == 8) return (connectDirect || connectRise || connectFall);
677  else B2WARNING("Unknown option for connect " << m_connect << ", using default.");
678  return (connectDirect || connectRise);
679 }
680 
681 bool
682 CDCTrigger2DFinderModule::connectedUD(unsigned patternD, unsigned patternU)
683 {
684  // connected if
685  // . . . .
686  // x . . x
687  // --- or ---
688  // x . . x
689  // . . . .
690  bool connectDirect = (((patternU >> 0) & 1) && ((patternD >> 2) & 1)) ||
691  (((patternU >> 1) & 1) && ((patternD >> 3) & 1));
692  // connected if
693  // . .
694  // . x
695  // ---
696  // x .
697  // . .
698  bool connectRise = ((patternU >> 1) & 1) && ((patternD >> 2) & 1);
699  // connected if
700  // . .
701  // x .
702  // ---
703  // . x
704  // . .
705  bool connectFall = ((patternU >> 0) & 1) && ((patternD >> 3) & 1);
706 
707  if (m_connect == 4) return connectDirect;
708  else if (m_connect == 6) return (connectDirect || connectRise);
709  else if (m_connect == 8) return (connectDirect || connectRise || connectFall);
710  else B2WARNING("Unknown option for connect " << m_connect << ", using default.");
711  return (connectDirect || connectRise);
712 }
713 
714 bool
715 CDCTrigger2DFinderModule::connectedDiag(unsigned patternLD, unsigned patternRU)
716 {
717  if (m_connect == 4) return false;
718 
719  // connected if
720  // . .
721  // x .
722  // . x
723  // . .
724  return (((patternRU >> 0) & 1) && ((patternLD >> 3) & 1));
725 }
726 
727 unsigned
728 CDCTrigger2DFinderModule::topRightSquare(vector<unsigned>& pattern)
729 {
730  // scan from top right corner until an active square is found
731  for (unsigned index = pattern.size() - 1; index > 0; --index) {
732  if (!pattern[index]) continue;
733  // check for ambiguity
734  unsigned ix = index % m_clusterSizeX;
735  unsigned iy = index / m_clusterSizeX;
736  if (ix < m_clusterSizeX - 1 && iy > 0) {
737  bool unique = true;
738  for (unsigned index2 = index - 1; index2 > 0; --index2) {
739  if (!pattern[index2]) continue;
740  unsigned ix2 = index2 % m_clusterSizeX;
741  unsigned iy2 = index2 / m_clusterSizeX;
742  if (iy2 < iy && ix2 > ix) {
743  unique = false;
744  break;
745  }
746  }
747  if (!unique) {
748  setReturnValue(false);
749  B2DEBUG(100, "topRightSquare not unique");
750  }
751  }
752  return index;
753  }
754  return 0;
755 }
756 
757 unsigned
758 CDCTrigger2DFinderModule::topRightCorner(unsigned pattern)
759 {
760  // scan pattern from right to left:
761  // 2 3
762  // 0 1
763  if ((pattern >> 3) & 1) return 3;
764  if ((pattern >> 1) & 1) {
765  if ((pattern >> 2) & 1) {
766  setReturnValue(false);
767  B2DEBUG(100, "topRightCorner not unique");
768  }
769  return 1;
770  }
771  if ((pattern >> 2) & 1) return 2;
772  return 0;
773 }
774 
775 unsigned
776 CDCTrigger2DFinderModule::bottomLeftCorner(unsigned pattern)
777 {
778  // scan pattern from left to right:
779  // 2 3
780  // 0 1
781  if (pattern & 1) return 0;
782  if ((pattern >> 2) & 1) {
783  if ((pattern >> 1) & 1) {
784  setReturnValue(false);
785  B2DEBUG(100, "bottomLeftCorner not unique");
786  }
787  return 2;
788  }
789  if ((pattern >> 1) & 1) return 1;
790  return 3;
791 }
Belle2::CDCTriggerHoughCand::getID
unsigned getID() const
Get candidate number.
Definition: CDCTriggerHoughETFModule.h:70
Belle2::getID
int getID(const std::vector< double > &breaks, double t)
get id of the time point t
Definition: calibTools.h:71
Belle2::CDCTriggerTrack
Track created by the CDC trigger.
Definition: CDCTriggerTrack.h:15
Belle2::cdcMap
std::map< int, cdcPair > cdcMap
Map of <counter, cdcPair>, for hits with indices.
Definition: CDCTriggerHoughETFModule.h:44
Belle2
Abstract base class for different kinds of events.
Definition: MillepedeAlgorithm.h:19
Belle2::coord2dPair
std::pair< TVector2, TVector2 > coord2dPair
Typedef for the coord2dPair, which represents the DATCONHoughCand.
Definition: DATCONHoughCand.h:32
Belle2::CDCTriggerHoughCand
Hough Candidates class.
Definition: CDCTriggerHoughETFModule.h:54
Belle2::CDCTriggerHoughCand::getSLcount
unsigned short getSLcount() const
Get super layer count.
Definition: CDCTriggerHoughETFModule.h:68
Belle2::cdcPair
std::pair< unsigned short, TVector2 > cdcPair
Pair of <iSuperLayer, (x, y)>, for hits in conformal space.
Definition: CDCTriggerHoughETFModule.h:42
Belle2::CDCTriggerHoughCluster
Cluster created by the Hough finder of the CDC trigger.
Definition: CDCTriggerHoughCluster.h:13