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