Belle II Software development
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
29using namespace std;
30
31namespace 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::rlrel
bool rlrel(vector< unsigned > a, vector< unsigned > b)
...rlrel...
Definition: PeakFinder.cc:157
Belle2::TRGCDCPeakFinder::~TRGCDCPeakFinder
virtual ~TRGCDCPeakFinder()
Destructor.
Definition: PeakFinder.cc:48
Belle2::TRGDebug::enterStage
static void enterStage(const std::string &stageName)
Declare that you enter new stage.
Definition: Debug.cc:24
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::udrel
bool udrel(vector< unsigned > a, vector< unsigned > b)
...udrel...
Definition: PeakFinder.cc:183
Belle2::TRGCDCPeakFinder::TRGCDCPeakFinder
TRGCDCPeakFinder(const std::string &name)
Contructor.
Definition: PeakFinder.cc:43
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::mirel
bool mirel(vector< unsigned > a, vector< unsigned > b)
...mirel...
Definition: PeakFinder.cc:212
Belle2::FindCP1
unsigned short FindCP1(vector< unsigned > a)
...Find center Pattern1 from Pattern 2
Definition: PeakFinder.cc:236
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::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::version
std::string version(void) const
returns version.
Definition: PeakFinder.cc:38
Belle2::FindP1C
unsigned FindP1C(vector< unsigned > a)
...Pattern1 Center...
Definition: PeakFinder.cc:303
Belle2
Abstract base class for different kinds of events.
Definition: MillepedeAlgorithm.h:17
std
STL namespace.