Belle II Software  release-08-01-10
PeakFinder.cc
1 /**************************************************************************
2  * basf2 (Belle II Analysis Software Framework) *
3  * Author: The Belle II Collaboration *
4  * *
5  * See git log for contributors and copyright holders. *
6  * This file is licensed under LGPL-3.0, see LICENSE.md. *
7  **************************************************************************/
8 
9 //-----------------------------------------------------------------------------
10 // Description : A class to find tracks usning Hough algorithm
11 //-----------------------------------------------------------------------------
12 
13 #define TRG_SHORT_NAMES
14 #define TRGCDC_SHORT_NAMES
15 
16 #include <cstdlib>
17 #include "trg/trg/Debug.h"
18 #include "trg/cdc/Circle.h"
19 #include "trg/cdc/PeakFinder.h"
20 #include "trg/cdc/HoughPlaneMulti2.h"
21 #include "trg/cdc/HoughTransformationCircle.h"
22 
24 
25 #include <iostream>
26 #include <math.h>
28 
29 using namespace std;
30 
31 namespace Belle2 {
37  string
38  TRGCDCPeakFinder::version(void) const
39  {
40  return string("TRGCDCPeakFinder 5.05");
41  }
42 
43  TRGCDCPeakFinder::TRGCDCPeakFinder(const string& name)
44  : _name(name)
45  {
46  }
47 
48  TRGCDCPeakFinder::~TRGCDCPeakFinder()
49  {
50  }
51 
52  void
53  TRGCDCPeakFinder::regions(TCHPlane& hp, const unsigned threshold) const
54  {
55 
56  TRGDebug::enterStage("Making regions");
57 
58  //...Search cells above threshold...
59  unsigned nCells = hp.nX() * hp.nY();
60  static unsigned* candidates =
61  (unsigned*) malloc(nCells * sizeof(unsigned));
62  unsigned nActive = 0;
63  for (unsigned j = 0; j < hp.nY(); j++) {
64  //minus x direction , plus -x direction
65  if ((hp.name()) == "circle hough minus")
66  for (unsigned i = 0; i < hp.nX(); i++) {
67  //...Threshold check...
68  const unsigned n = hp.entry(i, j);
69  if (n < threshold) continue;
70  candidates[nActive] = hp.serialId(i, j);
71  ++nActive;
72  } else
73  for (unsigned z = hp.nX(); z > 0 ; --z) {
74  //...Threshold check...
75  unsigned i = 0;
76  i = z - 1;
77  const unsigned n = hp.entry(i, j);
78  if (n < threshold) continue;
79  candidates[nActive] = hp.serialId(i, j);
80  ++nActive;
81  }
82  }
83 
84  if (TRGDebug::level())
85  cout << TRGDebug::tab() << "Active cells=" << nActive << endl;
86 
87  //...Loop over active cells...
88  const unsigned used = nCells;
89  for (unsigned i = 0; i < nActive; i++) {
90  if (candidates[i] == used) continue;
91  const unsigned id0 = candidates[i];
92  candidates[i] = used;
93 
94  //...Make a new region...
95  vector<unsigned>* region = new vector<unsigned>;
96  region->push_back(id0);
97  if (TRGDebug::level() > 2)
98  cout << TRGDebug::tab(4) << "new region made" << endl;
99 
100  //...Search neighbors...
101  for (unsigned j = 0; j < nActive; j++) {
102  if (candidates[j] == used) continue;
103  const unsigned id1 = candidates[j];
104 
105  unsigned x1 = 0;
106  unsigned y1 = 0;
107  hp.id(id1, x1, y1);
108  if (TRGDebug::level() > 2)
109  cout << TRGDebug::tab(8) << "cell:x=" << x1 << ",y=" << y1
110  << endl;
111 
112  for (unsigned k = 0; k < unsigned(region->size()); k++) {
113  unsigned id2 = (* region)[k];
114  unsigned x2 = 0;
115  unsigned y2 = 0;
116  hp.id(id2, x2, y2);
117  int difx = abs(int(x1) - int(x2));
118  int dify = abs(int(y1) - int(y2));
119  if (difx > (int) hp.nX() / 2) difx = hp.nX() - difx;
120  if (dify > (int) hp.nY() / 2) dify = hp.nY() - dify;
121 
122  if (TRGDebug::level() > 2) {
123  cout << TRGDebug::tab(12) << "x=" << x2 << ",y=" << y2
124  << ":difx=" << difx << ",dify=" << dify;
125  if ((difx < 2) && (dify < 3))
126  cout << " ... connected" << endl;
127  else
128  cout << endl;
129  }
130 
131  if ((difx < 2) && (dify < 3)) {
132  region->push_back(id1);
133  candidates[j] = used;
134  break;
135  }
136 
137  }
138  }
139  hp.setRegion(region);
140  }
141 
142  if (TRGDebug::level())
143  cout << TRGDebug::tab() << "Regions=" << hp.regions().size() << endl;
144 
145  TRGDebug::leaveStage("Making regions");
146  }
147 
149 
151  // right/left neighboring 2x2 squares
152  // a3 a4 | b3 b4
153  // a1 a2 | b1 b2
154  // unconnected if
155  // x x | . x or x . | x x
156  // x x | . x x . | x x
157  bool rlrel(vector<unsigned> a, vector<unsigned> b)
158  {
159  bool t = true;
160  //case 1
161  if (b[1] == 0 && b[3] == 0) {
162  t = false;
163  }
164  //case 2
165  if (a[2] == 0 && a[4] == 0) {
166  t = false;
167  }
168  return t;
169  }
171  // up/down neighboring 2x2 squares
172  // b3 b4
173  // b1 b2
174  // -----
175  // a3 a4
176  // a1 a2
177  // unconnected if
178  // x x x x . .
179  // x x . . . .
180  // --- or --- or ---
181  // . . . . x x
182  // . . x x x x
183  bool udrel(vector<unsigned> a, vector<unsigned> b)
184  {
185  bool t = true;
186  //case 1
187  if (a[1] == 0 && a[2] == 0 && a[3] == 0 && a[4] == 0) {
188  t = false;
189  }
190  //case 2
191  if (a[3] == 0 && a[4] == 0 && b[1] == 0 && b[2] == 0) {
192  t = false;
193  }
194  //case 3
195  if (b[1] == 0 && b[2] == 0 && b[3] == 0 && b[4] == 0) {
196  t = false;
197  }
198 
199  return t;
200  }
202  // diagonal neighboring 2x2 squares
203  // b3 b4
204  // b1 b2
205  // a3 a4
206  // a1 a2
207  // unconnected if
208  // x x x x . x
209  // . x x x . x
210  // x . or x . or x x
211  // x x x . x x
212  bool mirel(vector<unsigned> a, vector<unsigned> b)
213  {
214  bool t = true;
215  //case 1
216  if (a[4] == 0 && b[1] == 0) {
217  t = false;
218  }
219  //case 2
220  if (a[2] == 0 && a[4] == 0) {
221  t = false;
222  }
223  //case 3
224  if (b[1] == 0 && b[3] == 0) {
225  t = false;
226  }
227  return t;
228  }
229 
231  // Pattern 2 is a 3x2 rectangle of shape
232  // a4 a5
233  // a2 a3
234  // a0 a1
235  // function returns index of peak depending on pattern
236  unsigned short FindCP1(vector<unsigned> a)
237  {
238  unsigned center = 0;
239  //...1...
240  //(x on, . off, return O)
241  // . .
242  // . .
243  // O .
244  if (a[0] == 1 && a[1] == 0 && a[2] == 0 && a[3] == 0 && a[4] == 0 && a[5] == 0) {center = 0;}
245  //...2...
246  //(x on, . off, return O)
247  // . . . . . .
248  // . . x . . x
249  // O x O . O .
250  if (a[0] == 1 && a[1] == 1 && a[2] == 0 && a[3] == 0 && a[4] == 0 && a[5] == 0) {center = 0;}
251  if (a[0] == 1 && a[1] == 0 && a[2] == 1 && a[3] == 0 && a[4] == 0 && a[5] == 0) {center = 0;}
252  if (a[0] == 1 && a[1] == 0 && a[2] == 0 && a[3] == 1 && a[4] == 0 && a[5] == 0) {center = 0;}
253  //...3...
254  //(x on, . off, return O)
255  // . . . . . . x . . x x . . x
256  // x . . x O x O . O . . O . O
257  // O x x O x . x . x . x . x .
258  if (a[0] == 1 && a[1] == 1 && a[2] == 1 && a[3] == 0 && a[4] == 0 && a[5] == 0) {center = 0;}
259  if (a[0] == 1 && a[1] == 1 && a[2] == 0 && a[3] == 1 && a[4] == 0 && a[5] == 0) {center = 1;}
260  if (a[0] == 1 && a[1] == 0 && a[2] == 1 && a[3] == 1 && a[4] == 0 && a[5] == 0) {center = 2;}
261  if (a[0] == 1 && a[1] == 0 && a[2] == 1 && a[3] == 0 && a[4] == 1 && a[5] == 0) {center = 2;}
262  if (a[0] == 1 && a[1] == 0 && a[2] == 1 && a[3] == 0 && a[4] == 0 && a[5] == 1) {center = 2;}
263  if (a[0] == 1 && a[1] == 0 && a[2] == 0 && a[3] == 1 && a[4] == 1 && a[5] == 0) {center = 3;}
264  if (a[0] == 1 && a[1] == 0 && a[2] == 0 && a[3] == 1 && a[4] == 0 && a[5] == 1) {center = 3;}
265  //...4...
266  //(x on, . off, return O)
267  // . . x . . x x . . x x . . x x x x x
268  // x x O . O . . O . O O x O x O . . O
269  // O x x x x x x x x x x . x . x . x .
270  if (a[0] == 1 && a[1] == 1 && a[2] == 1 && a[3] == 1 && a[4] == 0 && a[5] == 0) {center = 0;}
271  if (a[0] == 1 && a[1] == 1 && a[2] == 1 && a[3] == 0 && a[4] == 1 && a[5] == 0) {center = 2;}
272  if (a[0] == 1 && a[1] == 1 && a[2] == 1 && a[3] == 0 && a[4] == 0 && a[5] == 1) {center = 2;}
273  if (a[0] == 1 && a[1] == 1 && a[2] == 0 && a[3] == 1 && a[4] == 1 && a[5] == 0) {center = 3;}
274  if (a[0] == 1 && a[1] == 1 && a[2] == 0 && a[3] == 1 && a[4] == 0 && a[5] == 1) {center = 3;}
275  if (a[0] == 1 && a[1] == 0 && a[2] == 1 && a[3] == 1 && a[4] == 1 && a[5] == 0) {center = 2;}
276  if (a[0] == 1 && a[1] == 0 && a[2] == 1 && a[3] == 1 && a[4] == 0 && a[5] == 1) {center = 2;}
277  if (a[0] == 1 && a[1] == 0 && a[2] == 1 && a[3] == 0 && a[4] == 1 && a[5] == 1) {center = 2;}
278  if (a[0] == 1 && a[1] == 0 && a[2] == 0 && a[3] == 1 && a[4] == 1 && a[5] == 1) {center = 3;}
279  //...5...
280  //(x on, . off, return O)
281  // x . . x x x x x x x
282  // O x O x O . . O O x
283  // x x x x x x x x x .
284  if (a[0] == 1 && a[1] == 1 && a[2] == 1 && a[3] == 1 && a[4] == 1 && a[5] == 0) {center = 2;}
285  if (a[0] == 1 && a[1] == 1 && a[2] == 1 && a[3] == 1 && a[4] == 0 && a[5] == 1) {center = 2;}
286  if (a[0] == 1 && a[1] == 1 && a[2] == 1 && a[3] == 0 && a[4] == 1 && a[5] == 1) {center = 2;}
287  if (a[0] == 1 && a[1] == 1 && a[2] == 0 && a[3] == 1 && a[4] == 1 && a[5] == 1) {center = 3;}
288  if (a[0] == 1 && a[1] == 0 && a[2] == 1 && a[3] == 1 && a[4] == 1 && a[5] == 1) {center = 2;}
289  //...6...
290  //(x on, . off, return O)
291  // x x
292  // O x
293  // x x
294  if (a[0] == 1 && a[1] == 1 && a[2] == 1 && a[3] == 1 && a[4] == 1 && a[5] == 1) {center = 2;}
295  return center;
296  }
297 
299  // Pattern 1 is a 2x2 square of shape
300  // a3 a4
301  // a1 a2
302  // function returns index of peak depending on pattern
303  unsigned FindP1C(vector<unsigned> a)
304  {
305  unsigned hits = 0;
306  unsigned short center = 0;
307  for (unsigned short k = 1; k < a.size(); k++) {
308  if (a[k] == 1) {hits++;}
309  }
310  //...1...
311  //(x on, . off, return O)
312  // . . . . O . . O
313  // O . . O . . . .
314  if (hits == 1) {
315  if (a[1] == 1 && a[2] == 0 && a[3] == 0 && a[4] == 0) {center = 1;}
316  if (a[1] == 0 && a[2] == 1 && a[3] == 0 && a[4] == 0) {center = 2;}
317  if (a[1] == 0 && a[2] == 0 && a[3] == 1 && a[4] == 0) {center = 3;}
318  if (a[1] == 0 && a[2] == 0 && a[3] == 0 && a[4] == 1) {center = 4;}
319 
320  }
321  //...2...
322  //(x on, . off, return O)
323  // . . x . . x x . . x O x
324  // O x O . O . . O . O . .
325  if (hits == 2) {
326  if (a[1] == 1 && a[2] == 1 && a[3] == 0 && a[4] == 0) {center = 1;}
327  if (a[1] == 1 && a[2] == 0 && a[3] == 1 && a[4] == 0) {center = 1;}
328  if (a[1] == 1 && a[2] == 0 && a[3] == 0 && a[4] == 1) {center = 1;}
329  if (a[1] == 0 && a[2] == 1 && a[3] == 1 && a[4] == 0) {center = 2;}
330  if (a[1] == 0 && a[2] == 1 && a[3] == 0 && a[4] == 1) {center = 2;}
331  if (a[1] == 0 && a[2] == 0 && a[3] == 1 && a[4] == 1) {center = 3;}
332 
333  }
334  //...3...
335  //(x on, . off, return O)
336  // x . . x O x x O
337  // O x x O x . . x
338  if (hits == 3) {
339  if (a[1] == 1 && a[2] == 1 && a[3] == 1 && a[4] == 0) {center = 1;}
340  if (a[1] == 1 && a[2] == 1 && a[3] == 0 && a[4] == 1) {center = 2;}
341  if (a[1] == 1 && a[2] == 0 && a[3] == 1 && a[4] == 1) {center = 3;}
342  if (a[1] == 0 && a[2] == 1 && a[3] == 1 && a[4] == 1) {center = 4;}
343 
344 
345  }
346  //...4...
347  //(x on, . off, return O)
348  // x x
349  // O x
350  if (hits == 4) {
351  if (a[1] == 1 && a[2] == 1 && a[3] == 1 && a[4] == 1) {center = 1;}
352 
353  }
354 
355  return center;
356  }
357 
358  void
359  TRGCDCPeakFinder::findPeaksTrasan(TCHPlaneMulti2& hp,
360  const unsigned threshold,
361  const bool centerIsPeak,
362  vector<unsigned>& peakSerialIds) const
363  {
364 
365  TRGDebug::enterStage("Peak Finding (trasan methode)");
366  if (TRGDebug::level())
367  cout << TRGDebug::tab() << "threshold=" << threshold
368  << ",plane name=[" << hp.name() << "]" << endl;
369 
370  //...Make ionnected regions (is this the best way???)...
371  regions(hp, threshold);
372 
373  //...Determine peaks...
374  const vector<vector<unsigned> *>& regions = hp.regions();
375  for (unsigned i = 0; i < (unsigned) regions.size(); i++) {
376 
377  if (TRGDebug::level() > 1)
378  cout << TRGDebug::tab() << "region " << i << " contents" << endl;
379 
380  //...Calculate size and center of a region...
381  const vector<unsigned>& r = * regions[i];
382  unsigned minX = hp.nX();
383  unsigned maxX = 0;
384  unsigned minY = hp.nY();
385  unsigned maxY = 0;
386  for (unsigned j = 0; j < (unsigned) r.size(); j++) {
387  const unsigned s = r[j];
388  unsigned x = 0;
389  unsigned y = 0;
390  hp.id(s, x, y);
391  if (x < minX) minX = x;
392  if (x > maxX) maxX = x;
393  if (y < minY) minY = y;
394  if (y > maxY) maxY = y;
395 
396  if (TRGDebug::level() > 1)
397  cout << TRGDebug::tab(4) << "x=" << x << ",y=" << y << endl;
398  }
399  const unsigned cX = minX + (maxX - minX) / 2;
400  const unsigned cY = minY + (maxY - minY) / 2;
401 
402  //...Determine a center of a region...
403  unsigned ncX = hp.nX() * hp.nY();
404  unsigned ncY = ncX;
405  if (TRGDebug::level() > 1)
406  cout << TRGDebug::tab() << "center of region:x=" << cX << ",y="
407  << cY << endl;
408 
409  if (! centerIsPeak) {
410 
411  if (TRGDebug::level() > 1)
412  cout << TRGDebug::tab() << "Searching a cell closest to the "
413  << "region center" << endl;
414 
415  //...Search for a cell which is the closest to the center...
416  float minDiff2 = float(hp.nX() * hp.nX() + hp.nY() * hp.nY());
417  for (unsigned j = 0; j < (unsigned) r.size(); j++) {
418  const unsigned s = r[j];
419  unsigned x = 0;
420  unsigned y = 0;
421  hp.id(s, x, y);
422 
423  const float diff2 = (float(x) - float(cX)) *
424  (float(x) - float(cX))
425  + (float(y) - float(cY)) *
426  (float(y) - float(cY));
427 
428  if (diff2 < minDiff2) {
429  minDiff2 = diff2;
430  ncX = x;
431  ncY = y;
432  }
433  if (TRGDebug::level() > 1)
434  cout << TRGDebug::tab(4) << "x=" << ncX << ",y=" << ncY
435  << ":diff2=" << diff2 << endl;
436  }
437  } else {
438 
439  //...Search for a peak...
440  float max = 0;
441  for (unsigned j = 0; j < (unsigned) r.size(); j++) {
442  const unsigned s = r[j];
443  const float entry = hp.entry(s);
444  if (max < entry) {
445  max = entry;
446  unsigned x = 0;
447  unsigned y = 0;
448  hp.id(s, x, y);
449  ncX = x;
450  ncY = y;
451  }
452  }
453  }
454 
455  //...Store the center cell...
456  const unsigned serialId = hp.serialId(ncX, ncY);
457  peakSerialIds.push_back(serialId);
458 
459  if (TRGDebug::level()) {
460  cout << TRGDebug::tab() << "region " << i << " final center:x="
461  << ncX << ",y=" << ncY << endl
462  << TRGDebug::tab(4) << "position in HP:x="
463  << hp.position(ncX, ncY).x() << ",y="
464  << hp.position(ncX, ncY).y() << endl;
465  }
466  }
467 
468  if (TRGDebug::level())
469  cout << TRGDebug::tab() << peakSerialIds.size() << " peak(s)"
470  << " found in total" << endl;
471 
472  TRGDebug::leaveStage("Peak Finding (trasan methode)");
473  return;
474  }
475 
476  void
477  TRGCDCPeakFinder::findPeaks(const TCHPlaneMulti2& hp,
478  const unsigned threshold,
479  vector<vector<unsigned>>& peaks) const
480  {
481 
482  const string sn = "Peak Finding";
483  TRGDebug::enterStage(sn);
484 
485  if (TRGDebug::level())
486  cout << TRGDebug::tab() << "threshold=" << threshold
487  << ",plane name=[" << hp.name() << "]" << endl;
488 
489  // p1p2(hp, threshold, peaks);
490  p1p2Methode(hp, threshold, peaks);
491 
492  TRGDebug::leaveStage(sn);
493 
494  return;
495  }
496 
497  void
498  TRGCDCPeakFinder::p1p2Methode(const TCHPlane& hp,
499  const unsigned threshold,
500  vector<vector<unsigned>>& peak_xy) const
501  {
502  const string sn = "p1p2";
503  TRGDebug::enterStage(sn);
504 
505  unsigned nCells = hp.nX() * hp.nY();
506  unsigned nX2 = hp.nX() / 2;
507  unsigned nY2 = hp.nY() / 2;
508 
509  //...Search cells above threshold...
510  static unsigned* candidates = (unsigned*) malloc(nCells * sizeof(unsigned));
511  unsigned nActive = 0;
512  for (unsigned j = 0; j < hp.nY(); ++j) {
513  //minus x direction, plus -x direction
514  if ((hp.name()) == "circle hough minus") {
515  for (unsigned i = 0; i < hp.nX(); ++i) {
516  //...Threshold check...
517  const unsigned n = hp.entry(i, j);
518  if (n < threshold) continue;
519  candidates[nActive] = hp.serialId(i, j);
520  ++nActive;
521  }
522  } else {
523  for (unsigned i = 0; i < hp.nX(); ++i) {
524  //...Threshold check...
525  const unsigned n = hp.entry(hp.nX() - i - 1, j);
526  if (n < threshold) continue;
527  candidates[nActive] = hp.serialId(hp.nX() - i - 1, j);
528  ++nActive;
529  }
530  }
531  }
532 
533  vector<unsigned> p;
534  vector<vector<unsigned>> p1m;
535  unsigned short no = 0;
536  //...create pattern1...begin
537  // divide the plane into squares of 2x2
538  // outer loop (n x m) goes over the squares,
539  // inner loop (j x k) goes over the cells in each square
540  for (unsigned n = 0; n < nY2; n++) {
541  for (unsigned m = 0; m < nX2; m++) {
542  unsigned a = m * 2;
543  unsigned b = n * 2;
544  bool ot = false;
545  ++no; // numbering starts at 1, not 0
546  p.push_back(no);
547 
548  //...find 4 cells...begin
549  for (unsigned j = 0; j < 2; j++) {
550  unsigned yy = b + j;
551  for (unsigned k = 0; k < 2; k++) {
552  unsigned xx = a + k;
553  // Plus plane transform (x axis mirrored compared to minus plane)
554  if ((hp.name()) == "circle hough plus") {
555  xx = hp.nX() - xx - 1;
556  }
557 
558  // go over the candidates and look for candidate in current cell
559  // if cell is peak candidate, add 1 to p, else add 0
560  unsigned short t = 0;
561  for (unsigned i = 0; i < nActive; i++) {
562  unsigned id1 = candidates[i];
563  unsigned x1 = 0;
564  unsigned y1 = 0;
565  hp.id(id1, x1, y1);
566  if (xx == x1 && yy == y1) {
567  t = 1;
568  ot = true;
569  break;
570  }
571  }
572  p.push_back(t);
573  }
574  }
575  //...find 4 cells...end
576 
577  // p = [n a b c d]
578  // n: number of 2x2 square
579  // a b c d: on/off flag for cells in 2x2 square
580  // shape:
581  // c d
582  // a b
583  if (ot == true) {
584  p1m.push_back(p);
585  }
586 
587  p.clear();
588  }
589  }
590  if (TRGDebug::level()) cout << TRGDebug::tab() << "size of p1m=" << p1m.size() << endl;
591  if (TRGDebug::level()) cout << TRGDebug::tab() << "~~~~~~~~~~~~~~~~~~~~~~~~~pattern1~~~~~~~~~~~~~~~~~~~" << endl;
592  //...create pattern1...end (output p1m)
593 
594 
595  //...Pattern2 & Find Peak...begin
596  if (TRGDebug::level()) cout << TRGDebug::tab() << ">>>>>>>>>>Pattern 2 & Find Peak Begin!!!>>>>>>>>>>" << endl;
597 
598  vector<unsigned> p0(5, 0);
599  vector<vector<unsigned>> op2;
600 
601  // loop over 2x2 peak candidates
602  for (unsigned short i = 0; i < p1m.size(); i++) {
603  unsigned short j = p1m[i][0]; // 2x2 square number (starting at 1)
604  unsigned short a = 0;
605  bool p1rel = false;
606  if (TRGDebug::level()) cout << TRGDebug::tab() << "no." << j << endl;
607 
608  // XYZ (begin)
609  // check for connections to neighboring 2x2 candidates
610  // if connection is found, continue to next candidate
611 
612  //X (horizontal connection to the left)
613  if ((j % nX2) == 1) {
614  a = j + nX2 - 1;
615  } else {
616  a = j - 1;
617  }
618  // loop over rest of candidates
619  for (unsigned k = 0; k < p1m.size(); k++) {
620  if (a == p1m[k][0]) {
621  // check connection to left neighbor
622  // by predefined subpattern in 2x2 square
623  if (!rlrel(p1m[k], p1m[i])) {
624  if (TRGDebug::level()) cout << TRGDebug::tab() << "no." << j << " & no." << a << " / X no rel" << endl;
625  p1rel = false;
626  } else {
627  if (TRGDebug::level()) cout << TRGDebug::tab() << "no." << j << " & no." << a << " / X rel" << endl;
628  p1rel = true;
629  }
630  break;
631  }
632  }
633  if (p1rel) {
634  continue;
635  }
636 
637  //Y (vertical connection to lower neighbor)
638  if (j > nX2) {
639  a = j - nX2;
640  // loop over rest of candidates
641  for (unsigned k = 0; k < p1m.size(); k++) {
642  if (a == p1m[k][0]) {
643  // check connection to lower neighbor
644  // by predefined subpattern in 2x2 square
645  if (!udrel(p1m[k], p1m[i])) {
646  if (TRGDebug::level()) cout << TRGDebug::tab() << "no." << j << " & no." << a << " / Y no rel" << endl;
647  p1rel = false;
648  } else {
649  if (TRGDebug::level()) cout << TRGDebug::tab() << "no." << j << " & no." << a << " / Y rel" << endl;
650  p1rel = true;
651  }
652  break;
653  }
654  }
655  if (p1rel) {
656  continue;
657  }
658 
659 
660  //Z (diagonal connection to lower left)
661  if ((j % nX2) == 1) {
662  a = j - 1;
663  } else {
664  a = j - nX2 - 1;
665  }
666  // loop over test of candidates
667  for (unsigned k = 0; k < p1m.size(); k++) {
668  if (a == p1m[k][0]) {
669  // check connection to lower left neighbor
670  // by predefined subpattern in 2x2 square
671  if (!mirel(p1m[k], p1m[i])) {
672  if (TRGDebug::level()) cout << TRGDebug::tab() << "no." << j << " & no." << a << " / Z no rel" << endl;
673  p1rel = false;
674  } else {
675  if (TRGDebug::level()) cout << TRGDebug::tab() << "no." << j << " & no." << a << " / Z rel" << endl;
676  p1rel = true;
677  }
678  break;
679  }
680  }
681  if (p1rel) {
682  continue;
683  }
684  }
685  // XYZ (end)
686 
687  // Pattern2 value (begin)
688  // make 3x2 rectangle of 2x2 squares
689  // p2v = [A B C D E F] (numbers of 2x2 squares in 3x2 rectangle)
690  // A: number of current 2x2 candidate
691  // shape of p2v:
692  // E F
693  // C D
694  // A B
695  vector<unsigned> p2v;
696  for (unsigned ip2 = 0; ip2 < 3; ++ip2) {
697  for (unsigned jp2 = 0; jp2 < 2; ++jp2) {
698  p2v.push_back(j + jp2 + ip2 * nX2);
699  }
700  if ((j % nX2) == 0) {
701  p2v[ip2 * 3 + 1] -= nX2;
702  }
703  }
704  // Pattern2 value(End)
705 
706  // make Pattern2(begin)
707  // get subpattern for each 2x2 square within 3x2 rectangle
708  // stored in op2
709 
710  // loop over 3x2 rectangle
711  for (unsigned short imp2 = 0; imp2 < p2v.size(); imp2++) {
712  unsigned short p2v_i = p2v[imp2];
713  // loop over 2x2 candidates
714  for (unsigned short jmp2 = 0; jmp2 < p1m.size(); jmp2++) {
715  unsigned short p1m_no = p1m[jmp2][0];
716  // if none of match number in p1m then pass this scan
717  if (p2v_i != p1m_no) {
718  if (jmp2 == (p1m.size() - 1)) {
719  // if no match is found use default (0 0 0 0 0)
720  op2.push_back(p0);
721  }
722  continue;
723  }
724  op2.push_back(p1m[jmp2]);
725 
726  break;
727  }
728  }
729  // make Pattern2(End)
730 
731  // Pattern2 relation(Begin)
732  // go over 3x2 rectangle and keep only cells connected to lower left 2x2 square
733  vector<vector<unsigned>> final_op2;
734  vector<unsigned> p2_state;
735  // A (start point)
736  final_op2.push_back(op2[0]);
737  p2_state.push_back(1);
738  // B (keep if connected to A)
739  if (rlrel(op2[0], op2[1])) {
740  final_op2.push_back(op2[1]);
741  p2_state.push_back(1);
742  } else {
743  final_op2.push_back(p0);
744  p2_state.push_back(0);
745  }
746  // C (keep if connected to A)
747  if (udrel(op2[0], op2[2])) {
748  final_op2.push_back(op2[2]);
749  p2_state.push_back(1);
750  } else {
751  final_op2.push_back(p0);
752  p2_state.push_back(0);
753  }
754  // D (keep connected to A, B or D)
755  if (mirel(op2[0], op2[3]) || udrel(op2[1], op2[3]) || rlrel(op2[2], op2[3])) {
756  final_op2.push_back(op2[3]);
757  p2_state.push_back(1);
758  } else {
759  final_op2.push_back(p0);
760  p2_state.push_back(0);
761  }
762  // E (keep if connected to C)
763  if (udrel(op2[2], op2[4])) {
764  final_op2.push_back(op2[4]);
765  p2_state.push_back(1);
766  } else {
767  final_op2.push_back(p0);
768  p2_state.push_back(0);
769  }
770  // F (keep if connected to C, D or E)
771  if (mirel(op2[2], op2[5]) || udrel(op2[3], op2[5]) || rlrel(op2[4], op2[5])) {
772  final_op2.push_back(op2[5]);
773  p2_state.push_back(1);
774  } else {
775  final_op2.push_back(p0);
776  p2_state.push_back(0);
777  }
778  // TODO: should compare connection to final_op2 instead of op2?
779  // otherwise e.g. op2 = A C-E with connection C-E
780  // would give final_op2 = A . E (should be A . .)
781 
782  // Pattern2 relation(End)
783 
784  // Find center peak(begin)
785  unsigned short fcpi = 0; // number of peak in 2x2 square (start: 1)
786  unsigned short fcpn = 0; // number of peak 2x2 square (start: 1)
787  unsigned short fcpx = 0; // x index in original hough plane (start: 0)
788  unsigned short fcpxs = 0; // x index in 2x2 square (0, 1)
789  unsigned short fcpy = 0; // y index in original hough plane (start: 0)
790  unsigned short fcpys = 0; // y index in 2x2 square (0, 1)
791 
792  // p2_state: vector of on/off flags for 3x3 square
793  // FindCP1(p2_state): peak within 3x3 square by predefined pattern
794  // op2[FindCP1(p2_state)]: 2x2 subpattern of this peak
795  // FindP1C(...): peak within 2x2 square by predefined pattern
796  fcpi = FindP1C(op2[FindCP1(p2_state)]);
797 
798  fcpn = op2[FindCP1(p2_state)][0];
799 
800  // get x index
801  if (fcpi >= 3) {
802  fcpxs = fcpi - 3;
803  } else {
804  fcpxs = fcpi - 1;
805  }
806  fcpx = ((fcpn - 1) % nX2) * 2 + fcpxs;
807  // Plus plane transform back to original numbering
808  if ((hp.name()) == "circle hough plus") {
809  fcpx = hp.nX() - fcpx - 1;
810  }
811 
812  // get y index
813  if (fcpi >= 3) {
814  fcpys = 1;
815  }
816  fcpy = fcpy + ((fcpn - 1) / nX2) * 2 + fcpys;
817 
818  if (TRGDebug::level()) cout << TRGDebug::tab() << "center of peak x=" << fcpx << " y=" << fcpy << endl;
819 
820  p.push_back(fcpx);
821  p.push_back(fcpy);
822  peak_xy.push_back(p);
823  p.clear();
824  // Find center peak(end)
825 
826  if (TRGDebug::level()) cout << TRGDebug::tab() << "~~~~~~~~~~Pattern 2 & Find Peak End!!!~~~~~~~~~~" << endl;
827  p2_state.clear();
828  final_op2.clear();
829  p2v.clear();
830  op2.clear();
831  } // end of loop over 2x2 candidates
832 
833  //... Pattern 2...end
834 
835  if (TRGDebug::level())
836  cout << TRGDebug::tab() << "total peaks=" << peak_xy.size() << endl;
837 
838  p1m.clear();
839 
840  TRGDebug::leaveStage(sn);
841  }
842 
844 } // namespace Belle2
Belle2::TRGDebug::tab
static std::string tab(void)
returns tab spaces.
Definition: Debug.cc:47
Belle2::TRGCDCPeakFinder::regions
void regions(TRGCDCHoughPlane &hp, const unsigned threshold) const
Makes regions.
Definition: PeakFinder.cc:53
Belle2::TRGCDCPeakFinder::~TRGCDCPeakFinder
virtual ~TRGCDCPeakFinder()
Destructor.
Definition: PeakFinder.cc:48
Belle2::rlrel
bool rlrel(vector< unsigned > a, vector< unsigned > b)
...rlrel...
Definition: PeakFinder.cc:157
Belle2::FindCP1
unsigned short FindCP1(vector< unsigned > a)
...Find center Pattern1 from Pattern 2
Definition: PeakFinder.cc:236
Belle2::udrel
bool udrel(vector< unsigned > a, vector< unsigned > b)
...udrel...
Definition: PeakFinder.cc:183
Belle2::TRGDebug::enterStage
static void enterStage(const std::string &stageName)
Declare that you enter new stage.
Definition: Debug.cc:24
Belle2::TRGCDCPeakFinder::findPeaksTrasan
void findPeaksTrasan(TRGCDCHoughPlaneMulti2 &hp, const unsigned threshold, const bool centerIsPeak, std::vector< unsigned > &peakSerialIds) const
do peak finding. This is a copy from "trasan".
Definition: PeakFinder.cc:359
Belle2::FindP1C
unsigned FindP1C(vector< unsigned > a)
...Pattern1 Center...
Definition: PeakFinder.cc:303
Belle2::TRGDebug::level
static int level(void)
returns the debug level.
Definition: Debug.cc:67
Belle2::TRGDebug::leaveStage
static void leaveStage(const std::string &stageName)
Declare that you leave a stage.
Definition: Debug.cc:34
Belle2::TRGCDCPeakFinder::findPeaks
void findPeaks(const TRGCDCHoughPlaneMulti2 &hp, const unsigned threshold, std::vector< std::vector< unsigned >> &peaks) const
do peak finding. Kaiyu's version using p1p2Methode.
Definition: PeakFinder.cc:477
Belle2::mirel
bool mirel(vector< unsigned > a, vector< unsigned > b)
...mirel...
Definition: PeakFinder.cc:212
Belle2::TRGCDCPeakFinder::p1p2Methode
void p1p2Methode(const TRGCDCHoughPlane &hp, const unsigned threshold, std::vector< std::vector< unsigned >> &peaks) const
Kaiyu's logic. Finds peaks from nested patterns.
Definition: PeakFinder.cc:498
Belle2
Abstract base class for different kinds of events.
Definition: MillepedeAlgorithm.h:17