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