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