CbmRoot
Loading...
Searching...
No Matches
CbmL1RichENNRingFinderParallel.cxx
Go to the documentation of this file.
1/* Copyright (C) 2010-2021 Frankfurt Institute for Advanced Studies, Goethe-Universitaet Frankfurt, Frankfurt
2 SPDX-License-Identifier: GPL-3.0-only
3 Authors: Ivan Kisel, Sergey Gorbunov, Igor Kulakov [committer] */
4
5/*
6 *====================================================================
7 *
8 * CBM Level 1 Reconstruction
9 *
10 * Authors: I.Kisel, S.Gorbunov
11 *
12 * e-mail : ikisel@kip.uni-heidelberg.de
13 *
14 *====================================================================
15 *
16 * Standalone RICH ring finder based on the Elastic Neural Net
17 *
18 *====================================================================
19 */
20
21#define PRINT_TIMING
22
23// CBM includes
24#include "CbmL1RichENNRingFinderParallel.h"
25
26#include "CbmEvent.h"
27#include "CbmRichHit.h"
28#include "CbmRichRing.h"
29#include "KfUtils.h"
30#include "TClonesArray.h"
31#include "TStopwatch.h"
32#include "assert.h"
33
34#include <algorithm>
35#include <cmath>
36#include <iostream>
37#include <vector>
38
39using std::cout;
40using std::endl;
41using std::fabs;
42using std::ios;
43using std::sqrt;
44using std::vector;
45
46namespace ca = cbm::algo::ca;
47namespace kf = cbm::algo::kf;
48
49CbmL1RichENNRingFinderParallel::CbmL1RichENNRingFinderParallel(Int_t /*verbose*/) : fRecoTime(0), fNEvents(0)
50{
51#ifdef PRINT_TIMING
52 TString name_tmp[NTimers] = {
53
54 "All",
55
56 "Ring finding",
57 "Ring finding: Prepare data",
58 "Ring finding: Init of param",
59 "Ring finding: Hit selection",
60 "Ring finding: Final fit",
61 "Ring finding: Store ring",
62
63 "Selection",
64
65 "Sort",
66 "Repack",
67 "Copy data"};
68 for (int i = 0; i < NTimers; i++) {
69 fTimersNames[i] = name_tmp[i];
70 }
71#endif // PRINT_TIMING
72}
73
74CbmL1RichENNRingFinderParallel::~CbmL1RichENNRingFinderParallel() {}
75
76
77void CbmL1RichENNRingFinderParallel::Init() {}
78
79Int_t CbmL1RichENNRingFinderParallel::DoFind(CbmEvent* event, TClonesArray* HitArray, TClonesArray* /*ProjArray*/,
80 TClonesArray* RingArray)
81{
82 if (!RingArray || !HitArray) return 0;
83
84 RingArray->Clear();
85#ifdef PRINT_TIMING
86 for (int i = 0; i < NTimers; i++) {
87 fTimers[i].Reset();
88 }
89#endif // PRINT_TIMING
90
91#ifdef PRINT_TIMING
92 GetTimer("Copy data").Start(0);
93#endif // PRINT_TIMING
94 vector<ENNRingHit> Up;
95 vector<ENNRingHit> Down;
96
97 const Int_t nhits = event ? event->GetNofData(ECbmDataType::kRichHit) : HitArray->GetEntriesFast();
98 for (Int_t i0 = 0; i0 < nhits; ++i0) {
99 Int_t i = event ? event->GetIndex(ECbmDataType::kRichHit, i0) : i0;
100 CbmRichHit* hit = dynamic_cast<CbmRichHit*>(HitArray->At(i));
101 if (!hit) continue;
102 ENNRingHit tmp;
103 tmp.x = hit->GetX();
104 tmp.y = hit->GetY();
105 tmp.outIndex = i;
106 tmp.quality = 0;
107 if (tmp.y > 0.) {
108 Up.push_back(tmp);
109 }
110 else {
111 Down.push_back(tmp);
112 }
113 }
114
115
116#ifdef PRINT_TIMING
117 GetTimer("Copy data").Stop();
118#endif // PRINT_TIMING
119
120#ifdef PRINT_TIMING
121 GetTimer("Sort").Start(0);
122#endif // PRINT_TIMING
123 sort(Up.begin(), Up.end(), ENNHit::Compare);
124 sort(Down.begin(), Down.end(), ENNHit::Compare);
125
126 // save local-out indices correspondece
127 cbm::algo::ca::Vector<Int_t> outIndicesUp; // indices in HitArray indexed by index in Up
128 cbm::algo::ca::Vector<Int_t> outIndicesDown; // indices in HitArray indexed by index in Down
129 outIndicesUp.reset(Up.size());
130 outIndicesDown.reset(Down.size());
131 for (THitIndex k = 0; k < Up.size(); k++) {
132 Up[k].localIndex = k;
133 outIndicesUp[k] = Up[k].outIndex;
134 }
135 for (THitIndex k = 0; k < Down.size(); k++) {
136 Down[k].localIndex = k;
137 outIndicesDown[k] = Down[k].outIndex;
138 }
139
140#ifdef PRINT_TIMING
141 GetTimer("Sort").Stop();
142#endif // PRINT_TIMING
143
144#ifdef PRINT_TIMING
145 GetTimer("Repack").Start(0);
146#endif // PRINT_TIMING
147
148 // save local-out indices correspondece
149 cbm::algo::ca::Vector<ENNHitV> UpV;
150 cbm::algo::ca::Vector<ENNHitV> DownV;
151 UpV.reset((Up.size() + fvec::size() - 1) / fvec::size());
152 DownV.reset((Down.size() + fvec::size() - 1) / fvec::size());
153 for (THitIndex k = 0; k < Up.size(); k++) {
154 int k_4 = k % fvec::size(), k_V = k / fvec::size();
155 ENNHitV& hits = UpV[k_V]; // TODO change on ENNHitV
156 hits.CopyHit(Up[k], k_4);
157 }
158 for (THitIndex k = 0; k < Down.size(); k++) {
159 int k_4 = k % fvec::size(), k_V = k / fvec::size();
160 ENNHitV& hits = DownV[k_V];
161 hits.CopyHit(Down[k], k_4);
162 }
163
164#ifdef PRINT_TIMING
165 GetTimer("Repack").Stop();
166#endif // PRINT_TIMING
167
168 vector<ENNRing> R;
169
170 float HitSize = .5;
171 THitIndex MinRingHits = 5;
172 float RMin = 2.;
173 float RMax = 5.;
174
175 TStopwatch timer;
176
177 ENNRingFinder(Up.size(), UpV, R, HitSize, MinRingHits, RMin, RMax);
178 ENNRingFinder(Down.size(), DownV, R, HitSize, MinRingHits, RMin, RMax);
179
180 timer.Stop();
181
182#ifdef PRINT_TIMING
183 GetTimer("Copy data").Start(0);
184#endif // PRINT_TIMING
185 Int_t NRings = 0;
186 for (vector<ENNRing>::iterator i = R.begin(); i != R.end(); ++i) {
187 if (!i->on) continue;
188 new ((*RingArray)[NRings]) CbmRichRing();
189 if (event != nullptr) event->AddData(ECbmDataType::kRichRing, NRings);
190 CbmRichRing* ring = dynamic_cast<CbmRichRing*>(RingArray->At(NRings));
191 NRings++;
192 ring->SetCenterX(i->x);
193 ring->SetCenterY(i->y);
194 ring->SetRadius(i->r);
195 ring->SetChi2(i->chi2);
196 const THitIndex NHits = i->localIHits.size();
197 for (THitIndex j = 0; j < NHits; j++) {
198 if (i->y > 0.)
199 ring->AddHit(outIndicesUp[i->localIHits[j]]);
200 else
201 ring->AddHit(outIndicesDown[i->localIHits[j]]);
202 }
203 }
204#ifdef PRINT_TIMING
205 GetTimer("Copy data").Stop();
206#endif // PRINT_TIMING
207
208 fNEvents++;
209 fRecoTime += timer.RealTime();
210 cout.setf(ios::fixed);
211 cout.setf(ios::showpoint);
212 cout.precision(4);
213 cout << "L1ENN Reco time stat/this [ms] : " << fRecoTime * 1000. / fNEvents << "/" << timer.RealTime() * 1000.
214 << endl;
215
216#ifdef PRINT_TIMING
217 static float timeStat[NTimers];
218 static bool firstCall = 1;
219 if (firstCall) {
220 for (int i = 0; i < NTimers; i++) {
221 timeStat[i] = 0;
222 }
223 firstCall = 0;
224 }
225
226 cout.width(30);
227 cout << "Timing statEvents / curEvent[ms] | statEvents[%] : " << endl;
228 cout.flags(ios::left | ios::oct | ios::showpoint | ios::fixed);
229 cout.precision(3);
230 for (int i = 0; i < NTimers; i++) {
231 timeStat[i] += fTimers[i].RealTime();
232 cout.width(30);
233 cout << fTimersNames[i] << " timer = " << timeStat[i] * 1000. / fNEvents << " / " << fTimers[i].RealTime() * 1000.
234 << " | " << timeStat[i] / timeStat[1] * 100 << endl;
235 }
236#endif // PRINT_TIMING
237 return NRings;
238}
239
240
241void CbmL1RichENNRingFinderParallel::ENNRingFinder(const int NHits, cbm::algo::ca::Vector<ENNHitV>& HitsV,
242 vector<ENNRing>& Rings, float HitSize, THitIndex MinRingHits,
243 fvec /*RMin*/, fvec RMax)
244{
245#ifdef PRINT_TIMING
246 GetTimer("All").Start(0);
247#endif // PRINT_TIMING
248 // INITIALIZATION
249
250 const fvec MinRingHits_v = MinRingHits;
251 const fvec Rejection = .5;
252 const float ShadowSize = HitSize / 4;
253 const int StartHitMaxQuality = 15;
254 const int SearchHitMaxQuality = 100; // TODO DELME
255
256 const fvec R2MinCut = 3. * 3., R2MaxCut = 7. * 7.;
257 // const fvec R2Min = RMin*RMin, R2Max = RMax*RMax;
258 const fvec HitSize2 = 2 * HitSize;
259 // const fvec HitSize4 = 4 * HitSize;
260 const fvec HitSizeSq_v = HitSize * HitSize;
261 const float HitSizeSq = HitSizeSq_v[0];
262 // const fvec HitSizeSq4 = HitSize2 * HitSize2;
263 const float AreaSize = 2 * RMax[0] + HitSize;
264 const float AreaSize2 = AreaSize * AreaSize;
265
266 // typedef vector<ENNRingHit>::iterator iH;
267 // typedef vector<ENNRingHit*>::iterator iP;
268
269 THitIndex NInRings = Rings.size();
270
271 // MAIN LOOP OVER HITS
272#ifdef PRINT_TIMING
273 GetTimer("Ring finding").Start(0);
274#endif // PRINT_TIMING
275 cbm::algo::ca::Vector<ENNSearchHitV> SearchArea;
276 cbm::algo::ca::Vector<ENNHitV> PickUpArea;
277 const int MaxAreaHits = 200; // TODO take into account NHits
278
279 SearchArea.reset(MaxAreaHits, ENNSearchHitV());
280 PickUpArea.reset(MaxAreaHits);
281
282 // TODO 1
283#if 0
284 cbm::algo::ca::Vector<ENNRingV> rings_tmp; // simd hits for tmp store
285 rings_tmp.reset(100); // TODO use NRings
286 int nRings_tmp = 0;
287#endif
288
289 // ENNRingHit* ileft = &(Hits[0]), *iright = ileft;//, i_main = ileft;
290 int ileft = 0, iright = ileft;
291 // int ileft[fvec::size()] = {0, 0, 0, 0};
292 // int iright[fvec::size()] = {0, 0, 0, 0};
293
294 THitIndex i_mains[fvec::size()] = {0};
295
296 THitIndex i_main_array[NHits]; // need for proceed in paralled almost independent areas
297 for (THitIndex i = 0; i < NHits; i++) {
298 i_main_array[i] = i; // TODO
299 }
300
301
302 for (THitIndex ii_main = 0; ii_main < NHits;) {
303
304#ifdef PRINT_TIMING
305 GetTimer("Ring finding: Prepare data").Start(0);
306#endif // PRINT_TIMING
307
308
309 fvec S0, S1, S2, S3, S4, S5, S6, S7;
310 fvec X = 0, Y = 0, R = 0, R2 = 0;
311
312 fmask validRing(fmask::One()); // mask of the valid rings
313
314 fvec SearchAreaSize = 0; // number of hits to fit and search ring
315 fvec PickUpAreaSize = 0;
316
317 for (size_t i_4 = 0; (i_4 < fvec::size()) && (ii_main < NHits); ii_main++) {
318 const THitIndex i_main = i_main_array[ii_main];
319 const int i_main_4 = i_main % fvec::size(), i_main_V = i_main / fvec::size();
320 ENNHitV* i = &HitsV[i_main_V];
321 if (i->quality[i_main_4] >= StartHitMaxQuality) continue; // already found hit
322
323 i_mains[i_4] = i_main;
324
325 float left = i->x[i_main_4] - AreaSize;
326 float right = i->x[i_main_4] + AreaSize;
327
328 // while( (HitsV[ileft[i_4]/fvec::size()] .x[ileft[i_4]%fvec::size()] < left ) ) ++ileft[i_4];
329 // while( (iright[i_4] < NHits) && (HitsV[iright[i_4]/fvec::size()].x[ileft[i_4]%fvec::size()] < right) ) ++iright[i_4];
330 // for( int j = ileft[i_4]; j < iright[i_4]; ++j ){
331 while ((HitsV[ileft / fvec::size()].x[ileft % fvec::size()] < left))
332 ++ileft; // TODO SIMDize
333 while ((iright < NHits) && (HitsV[iright / fvec::size()].x[ileft % fvec::size()] < right))
334 ++iright;
335
336 for (int j = ileft; j < iright; ++j) {
337 const int j_4 = j % fvec::size(), j_V = j / fvec::size();
338 ENNSearchHitV& sHit = SearchArea[int(SearchAreaSize[i_4])];
339 sHit.CopyHit(HitsV[j_V], j_4, i_4);
340
341 sHit.ly[i_4] = sHit.y[i_4] - i->y[i_main_4];
342 // if( fabs(sHit.ly) > AreaSize ) continue;
343 sHit.lx[i_4] = sHit.x[i_4] - i->x[i_main_4];
344 sHit.S0 = sHit.lx * sHit.lx;
345 sHit.S1 = sHit.ly * sHit.ly;
346 sHit.lr2 = sHit.S0 + sHit.S1;
347 // if(( sHit.lr2 > AreaSize2 ) || ( j == i )) continue; // no difference in speed
348 if (sHit.lr2[i_4] > AreaSize2) continue;
349 if (kf::utils::IsUnlikely(j == i_main)) continue;
350
351 if (sHit.quality[i_4] >= SearchHitMaxQuality) { // CHECKME do we really need PickUpArea
352 PickUpArea[static_cast<int>(PickUpAreaSize[i_4]++)].CopyHit(HitsV[j_V], j_4, i_4);
353 continue;
354 }
355
356 SearchAreaSize[i_4]++;
357 } // hits
358
359 // if( NAreaHits+1 < MinRingHits ) continue;
360 if (SearchAreaSize[i_4] < MinRingHits - 1) continue;
361 i_4++;
362 } // i_4
363
364 ENNHitV iHit;
365 int MaxSearchAreaSize = 0;
366 int MaxPickUpAreaSize = 0;
367 for (size_t i_4 = 0; i_4 < fvec::size(); i_4++) {
368 iHit.CopyHit(HitsV[i_mains[i_4] / fvec::size()], i_mains[i_4] % fvec::size(), i_4);
369 MaxSearchAreaSize = (MaxSearchAreaSize < SearchAreaSize[i_4]) ? int(SearchAreaSize[i_4]) : MaxSearchAreaSize;
370 MaxPickUpAreaSize = (MaxPickUpAreaSize < PickUpAreaSize[i_4]) ? int(PickUpAreaSize[i_4]) : MaxPickUpAreaSize;
371 }
372#ifdef PRINT_TIMING
373 GetTimer("Ring finding: Prepare data").Stop();
374 // assert( MaxSearchAreaSize[i_4] <= MaxAreaHits );
375#endif // PRINT_TIMING
376
377#ifdef PRINT_TIMING
378 GetTimer("Ring finding: Init of param").Start(0);
379#endif
380 for (int isa = 0; isa < MaxSearchAreaSize; isa++) { // TODO don't work w\o this because of nan in wights
381 ENNSearchHitV& sHit = SearchArea[isa];
382 const fmask validHit = (fvec(isa) < SearchAreaSize) & validRing;
383 sHit.lx = iif(validHit, sHit.lx, fvec::Zero());
384 sHit.ly = iif(validHit, sHit.ly, fvec::Zero());
385 sHit.lr2 = iif(validHit, sHit.lr2, fvec::Zero());
386 }
387
388 // initialize hits in the search area
389 fvec Dmax = 0.;
390 S0 = S1 = S2 = S3 = S4 = 0.;
391 for (int ih = 0; ih < MaxSearchAreaSize; ih++) {
392 ENNSearchHitV& sHit = SearchArea[ih];
393 const fmask validHit = (fvec(ih) < SearchAreaSize) & validRing;
394
395 fvec& lr2 = sHit.lr2;
396 fvec lr = sqrt(lr2);
397 Dmax = iif((lr > Dmax) & validHit, lr, Dmax);
398
399 sHit.S2 = sHit.lx * sHit.ly;
400 sHit.S3 = sHit.lx * lr2;
401 sHit.S4 = sHit.ly * lr2;
402 sHit.C = -lr * fvec(.5);
403
404 const fvec w = iif(validHit, iif(lr > fvec(1.E-4), fvec(1.) / lr, fvec::One()), fvec::Zero());
405 const fvec w2 = w * w;
406 sHit.Cx = w * sHit.lx;
407 sHit.Cy = w * sHit.ly;
408 S0 += w2 * sHit.S0;
409 S1 += w2 * sHit.S1;
410 S2 += w2 * sHit.S2;
411 S3 += w2 * sHit.S3;
412 S4 += w2 * sHit.S4;
413 }
414
415 // end of initialization of the search area
416#ifdef PRINT_TIMING
417 GetTimer("Ring finding: Init of param").Stop();
418 // assert( MaxSearchAreaSize[i_4] <= MaxAreaHits );
419#endif // PRINT_TIMING
420
421#ifdef PRINT_TIMING
422 // loop for minimization of E and noise rejection
423 GetTimer("Ring finding: Hit selection").Start(0);
424#endif // PRINT_TIMING
425 do {
426 // calculate parameters of the current ring
427 fvec tmp = S0 * S1 - S2 * S2;
428
429 // if( fabs(tmp) < 1.E-10 ) break; // CHECKME
430 tmp = fvec(0.5) / tmp;
431 X = (S3 * S1 - S4 * S2) * tmp;
432 Y = (S4 * S0 - S3 * S2) * tmp;
433 R2 = X * X + Y * Y;
434 R = sqrt(R2);
435
436 const fvec Dcut = Dmax * Rejection; // cut for noise hits
437 // float RingCut = HitSize4 * R ; // closeness
438 S0 = S1 = S2 = S3 = S4 = fvec::Zero();
439 Dmax = -1.;
440
441 for (THitIndex ih = 0; ih < MaxSearchAreaSize; ih++) {
442 const fmask validHit = (fvec(ih) < SearchAreaSize) & validRing;
443 // ENNHit *j = &(SearchArea[ih]);
444 ENNSearchHitV& sHit = SearchArea[ih];
445 const fvec dx = sHit.lx - X;
446 const fvec dy = sHit.ly - Y;
447 const fvec d = abs(sqrt(dx * dx + dy * dy) - R);
448 const fvec dSq = d * d;
449 sHit.on_ring = (d <= HitSize) & validHit;
450 const fvec dp = iif(sHit.on_ring, fvec(-1.), abs(sHit.C + sHit.Cx * X + sHit.Cy * Y));
451 Dmax = iif(((dp <= Dcut) & (dp > Dmax)), dp, Dmax);
452
453 fvec w = iif((sHit.on_ring), fvec(1.) / (HitSizeSq_v + abs(dSq)), fvec(1.) / (fvec(1.e-5) + abs(dSq)));
454 w = iif((dp <= Dcut) & validHit, w, fvec::Zero());
455 S0 += w * sHit.S0;
456 S1 += w * sHit.S1;
457 S2 += w * sHit.S2;
458 S3 += w * sHit.S3;
459 S4 += w * sHit.S4;
460 }
461
462 } while ((Dmax > fvec(0.)).isNotEmpty());
463
464
465#ifdef PRINT_TIMING
466 GetTimer("Ring finding: Hit selection").Stop();
467 // store the ring if it is found
468
469 GetTimer("Ring finding: Final fit").Start(0);
470#endif // PRINT_TIMING
471
472 fvec NRingHits = 1;
473 { // final fit of 3 parameters (X,Y,R)
474 S0 = S1 = S2 = S3 = S4 = S5 = S6 = S7 = 0.0;
475 for (THitIndex ih = 0; ih < MaxSearchAreaSize; ih++) {
476 ENNSearchHitV& sHit = SearchArea[ih];
477 const fvec w = iif(SearchArea[ih].on_ring, fvec::One(), fvec::Zero());
478 S0 += w * sHit.S0;
479 S1 += w * sHit.S1;
480 S2 += w * sHit.S2;
481 S3 += w * sHit.S3;
482 S4 += w * sHit.S4;
483 S5 += w * sHit.lx;
484 S6 += w * sHit.ly;
485 S7 += w * sHit.lr2;
486 NRingHits += w;
487 }
488 const fvec s0 = S6 * S0 - S2 * S5;
489 const fvec s1 = S0 * S1 - S2 * S2;
490 const fvec s2 = S0 * S4 - S2 * S3;
491 // if( fabs(s0) < 1.E-6 || fabs(s1) < 1.E-6 ) continue; // CHECKME
492 const fvec tmp = s1 * (S5 * S5 - NRingHits * S0) + s0 * s0;
493 const fvec A = ((S0 * S7 - S3 * S5) * s1 - s2 * s0) / tmp;
494 Y = (s2 + s0 * A) / s1 * fvec(0.5);
495 X = (S3 + S5 * A - S2 * Y * 2) / S0 * fvec(0.5);
496 R2 = X * X + Y * Y - A;
497 // if( R2 < 0 ) continue; // will be rejected letter by R2 > R2Min
498 R = sqrt(R2);
499 } // end of the final fit
500
501 validRing =
502 !((NRingHits < MinRingHits_v) | (R2 > R2MaxCut) | (R2 < R2MinCut)); // ghost suppresion // TODO constants
503// cout << validRing << endl;
504#ifdef PRINT_TIMING
505 GetTimer("Ring finding: Final fit").Stop();
506#endif // PRINT_TIMING
507
508
509#ifdef PRINT_TIMING
510 GetTimer("Ring finding: Store ring").Start(0);
511#endif // PRINT_TIMING
512 if (kf::utils::IsUnlikely(validRing.isEmpty())) continue;
513
515#if 0 // TODO 1
516 {
517 ENNRingV &ringV = rings_tmp[nRings_tmp++];
518
519 ringV.localIHits.reserve(25);
520 ringV.x = iHit.x + X;
521 ringV.y = iHit.y + Y;
522 ringV.r = R;
523 ringV.localIHits.push_back(iHit.localIndex);
524 ringV.NHits = 1;
525 ringV.chi2 = 0;
526 cbm::algo::ca::Vector<fvec> Shadow; // save hit indices of hits, which's quality will be changed
527 Shadow.reserve(25);
528 Shadow.push_back(iHit.localIndex);
529 for( THitIndex ih = 0; ih < MaxSearchAreaSize; ih++){
530 fvec validHit = ih < SearchAreaSize;
531
532 ENNSearchHitV &sHit = SearchArea[ih];
533 const fvec dx = sHit.x - ringV.x;
534 const fvec dy = sHit.y - ringV.y;
535 const fvec d = abs( sqrt(dx*dx+dy*dy) - ringV.r );
536 validHit = validHit & ( d <= HitSize );
537 ringV.chi2 += d*d;
538 ringV.localIHits.push_back( iif( validHit, sHit.localIndex, fvec(-1.) ) );
539 ringV.NHits += iif(validHit, fvec::One(), fvec::Zero());
540 validHit = validHit & ( d <= ShadowSize ); // TODO check *4
541 if ( validHit.isEmpty() ) continue; // CHECKME
542 Shadow.push_back( iif( validHit, sHit.localIndex, fvec(-1.) ) );
543 }
544 for( int ipu = 0; ipu < MaxPickUpAreaSize; ipu++ ) {
545 fvec validHit = ipu < PickUpAreaSize;
546
547 ENNHitV &puHit = PickUpArea[ipu];
548 const fvec dx = puHit.x - ringV.x;
549 const fvec dy = puHit.y - ringV.y;
550 const fvec d = abs( sqrt(dx*dx+dy*dy) - ringV.r );
551 validHit = validHit & ( d <= HitSize );
552 if ( validHit.isEmpty() ) continue;
553 ringV.chi2 += d*d;
554 ringV.localIHits.push_back( iif( validHit, puHit.localIndex, fvec(-1.) ) );
555 ringV.NHits += iif(validHit, fvec::One(), fvec::Zero());
556 validHit = validHit & ( d <= ShadowSize ); // TODO check *4
557 if ( validHit.isEmpty() ) continue; // CHECKME
558 Shadow.push_back( iif( validHit, puHit.localIndex, fvec(-1.) ) );
559 }
560
561 ringV.chi2 = ringV.chi2 / (( ringV.NHits - 3)*HitSizeSq);
562 const fvec quality = ringV.NHits;
563
564
565 // for( iP j = ringV.Hits.begin(); j != ringV.Hits.end(); ++j ){
566 // if( (*j)->quality<quality ) (*j)->quality = quality;
567 // }
568
569 //quality *= ShadowOpacity;
570 for( int i_4 = 0; (i_4 < fvec::size()); i_4++) {
571 const int NShadow = Shadow.size();
572 for( int is = 0; is < NShadow; is++ ) { // CHECKME change loops to speed up?
573 cout << i_4 << Shadow[is] << endl;
574 float ih_f = (Shadow[is])[i_4]-200;
575 if (ih_f == -1) continue;
576 int ih = static_cast<int>(ih_f); // TODO ! problem in conversion...
577 float ih_f2 = static_cast<float>(ih);
578 cout << ih_f << " " << ih << " " << ih_f2 << " " << ih%fvec::size() << " " << ih/fvec::size() << endl;
579
580 const THitIndex ih_4 = ih%fvec::size();
581 ENNHitV & hitV = HitsV[ih/fvec::size()];
582
583 // hitV.quality[ih_4] = ( hitV.quality[ih_4] < quality[i_4] ) ? quality[i_4] : hitV.quality[ih_4];
584 // shHit->quality = iif( shHit->quality < quality, quality, shHit->quality );
585 }
586 } // i_4
587 }
588#endif // 0
590
591 for (size_t i_4 = 0; (i_4 < fvec::size()); i_4++) {
592 // if( NRingHits < MinRingHits || R2 > R2Max || R2 < R2Min ) continue;
593
594 if (/*kf::utils::IsUnlikely*/ (!validRing[i_4])) continue;
595
596 ENNRing voidRing;
597 Rings.push_back(voidRing);
598 ENNRing& ring = Rings.back();
599 ring.localIHits.reserve(15);
600 ring.x = iHit.x[i_4] + X[i_4];
601 ring.y = iHit.y[i_4] + Y[i_4];
602 ring.r = R[i_4];
603 ring.localIHits.push_back(static_cast<THitIndex>(iHit.localIndex[i_4]));
604 ring.NHits = 1;
605 ring.chi2 = 0;
606 vector<THitIndex> Shadow; // save hit indices of hits, which's quality will be changed
607 Shadow.reserve(15);
608 Shadow.push_back(static_cast<THitIndex>(iHit.localIndex[i_4]));
609 for (THitIndex ih = 0; ih < SearchAreaSize[i_4]; ih++) {
610 ENNSearchHitV& sHit = SearchArea[ih];
611 const float dx = sHit.x[i_4] - ring.x;
612 const float dy = sHit.y[i_4] - ring.y;
613 const float d = fabs(sqrt(dx * dx + dy * dy) - R[i_4]);
614 if (kf::utils::IsUnlikely(d <= HitSize)) {
615 ring.chi2 += d * d;
616 ring.localIHits.push_back(int(sHit.localIndex[i_4]));
617 ring.NHits++;
618 if (d <= ShadowSize) Shadow.push_back(static_cast<THitIndex>(sHit.localIndex[i_4]));
619 }
620 }
621 for (int ipu = 0; ipu < PickUpAreaSize[i_4]; ipu++) {
622 ENNHitV& puHit = PickUpArea[ipu];
623 const float dx = puHit.x[i_4] - ring.x;
624 const float dy = puHit.y[i_4] - ring.y;
625 const float d = fabs(sqrt(dx * dx + dy * dy) - ring.r);
626 if (kf::utils::IsUnlikely(d <= HitSize)) {
627 ring.chi2 += d * d;
628 ring.localIHits.push_back(static_cast<THitIndex>(puHit.localIndex[i_4]));
629 ring.NHits++;
630 if (d <= ShadowSize) Shadow.push_back(static_cast<THitIndex>(puHit.localIndex[i_4]));
631 }
632 }
633 if (kf::utils::IsUnlikely(ring.NHits < MinRingHits)) {
634 Rings.pop_back();
635 continue;
636 }
637 ring.chi2 = ring.chi2 / ((ring.NHits - 3) * HitSizeSq);
638 int quality = ring.NHits;
639
640 // for( iP j = ring.Hits.begin(); j != ring.Hits.end(); ++j ){
641 // if( (*j)->quality<quality ) (*j)->quality = quality;
642 // }
643
644 //quality *= ShadowOpacity;
645 const int NShadow = Shadow.size();
646 for (int is = 0; is < NShadow; is++) {
647 const THitIndex ih = Shadow[is];
648 const THitIndex ih_4 = ih % fvec::size();
649 ENNHitV& hitV = HitsV[ih / fvec::size()];
650
651 hitV.quality[ih_4] = (hitV.quality[ih_4] < quality) ? quality : hitV.quality[ih_4];
652 // shHit->quality = iif( shHit->quality < quality, quality, shHit->quality );
653 }
654 } // i_4
655
656
657#ifdef PRINT_TIMING
658 GetTimer("Ring finding: Store ring").Stop();
659#endif // PRINT_TIMING
660 } // i_main
661#ifdef PRINT_TIMING
662 GetTimer("Ring finding").Stop();
663
664 // SELECTION OF RINGS
665 GetTimer("Selection").Start(0);
666#endif // PRINT_TIMING
667 typedef vector<ENNRing>::iterator iR;
668 iR Rbeg = Rings.begin() + NInRings;
669 iR Rend = Rings.end();
670
671//#define NEW_SELECTION
672#ifdef NEW_SELECTION // TODO optimize. at the moment just creates additional ghosts
673
674 sort(Rings.begin() + NInRings, Rings.end(), ENNRing::CompareENNHRings);
675
676 const int NHitsV = HitsV.size();
677 for (int ih = 0; ih < NHitsV; ih++)
678 HitsV[ih].quality = 0.;
679
680 for (iR ir = Rbeg; ir != Rend; ++ir) {
681 ir->on = 0;
682 ir->NOwn = ir->NHits;
683 }
684
685 for (iR i = Rbeg; i != Rend; ++i) {
686 if (i->skip) continue;
687 for (iR j = i + 1; j != Rend; ++j) {
688 if (j->skip) continue;
689
690 float dist = j->r + i->r + HitSize2[0];
691 float distCentr = sqrt((j->x - i->x) * (j->x - i->x) + (j->y - i->y) * (j->y - i->y));
692 if (distCentr > dist) continue;
693 Int_t NOverlaped = 0;
694
695 const THitIndex maxI = i->localIHits.size();
696 const THitIndex maxJ = j->localIHits.size();
697 for (THitIndex n = 0; n < maxI; n++)
698 for (THitIndex m = 0; m < maxJ; ++m)
699 if (i->localIHits[n] == j->localIHits[m]) NOverlaped++;
700 ENNRing* BigRing = 0;
701 if (NOverlaped > 0.7 * maxI) BigRing = &(*i);
702 if (NOverlaped > 0.7 * maxJ) BigRing = &(*j);
703 if (BigRing != 0) {
704 std::vector<THitIndex> newIndices;
705 for (THitIndex n = 0; n < maxI; n++) {
706 bool IsNew = 1;
707 for (THitIndex m = 0; m < maxJ; ++m)
708 if (i->localIHits[n] == j->localIHits[m]) IsNew = 0;
709 if (IsNew) newIndices.push_back(i->localIHits[n]);
710 }
711 if (maxI > maxJ) {
712 j->x = i->x;
713 j->y = i->y;
714 j->r = i->r;
715 }
716 const THitIndex newISize = newIndices.size();
717 for (THitIndex in = 0; in < newISize; in++)
718 j->localIHits.push_back(newIndices[in]);
719 i->skip = 1;
720 i->on = 0;
721 break;
722 }
723
724 } // j
725 } // i
726
727
728 for (iR i = Rbeg; i != Rend; ++i) {
729 if (!(i->skip)) {
730 i->on = 1;
731 float S0, S1, S2, S3, S4, S5, S6, S7, X, Y, R;
732 S0 = S1 = S2 = S3 = S4 = S5 = S6 = S7 = 0.0;
733
734
735 const THitIndex firstIh = i->localIHits[0];
736 const ENNHitV& firstHit = HitsV[firstIh / fvec::size()];
737 const int firstIh_4 = firstIh % fvec::size();
738 const THitIndex maxI = i->localIHits.size();
739
740 vector<ENNSearchHitV> shits;
741 shits.reset(maxI);
742 for (THitIndex iih = 0; iih < maxI; iih++) {
743 const THitIndex ih = i->localIHits[iih];
744 const ENNHitV& hit = HitsV[ih / fvec::size()];
745 const int ih_4 = ih % fvec::size();
746 ENNSearchHitV& shit = shits[iih];
747
748 shit.ly[0] = hit.y[ih_4] - firstHit.y[firstIh_4];
749 shit.lx[0] = hit.x[ih_4] - firstHit.x[firstIh_4];
750 shit.S0[0] = shit.lx[0] * shit.lx[0];
751 shit.S1[0] = shit.ly[0] * shit.ly[0];
752 shit.lr2[0] = shit.S0[0] + shit.S1[0];
753 float lr2 = shit.lr2[0];
754 float lr = sqrt(lr2);
755 shit.S2[0] = shit.lx[0] * shit.ly[0];
756 shit.S3[0] = shit.lx[0] * lr2;
757 shit.S4[0] = shit.ly[0] * lr2;
758 shit.C[0] = -lr * 0.5;
759 float w;
760 if (lr > 1.E-4)
761 w = 1. / lr;
762 else
763 w = 1.;
764 shit.Cx[0] = w * shit.lx[0];
765 shit.Cy[0] = w * shit.ly[0];
766 }
767 float Dmax = -1.;
768
769 X = i->x - firstHit.x[firstIh_4];
770 Y = i->y - firstHit.y[firstIh_4];
771 R = i->r;
772 int search_stop = 0;
773 do { // final fit of 3 parameters (X,Y,R)
774 float Dcut = Dmax * Rejection[0];
775 int n = 0;
776 S0 = S1 = S2 = S3 = S4 = S5 = S6 = S7 = 0.0;
777
778 for (THitIndex ih = 0; ih < maxI; ih++) {
779 ENNSearchHitV& shit = shits[ih];
780
781 float dx = shit.lx[0] - X;
782 float dy = shit.ly[0] - Y;
783 float d = fabs(sqrt(dx * dx + dy * dy) - R);
784 shit.on_ring[0] = (d <= HitSize2[0]);
785 float w;
786 if (shit.on_ring[0]) {
787 n++;
788 w = 1;
789 }
790 else {
791 float dp = fabs(shit.C[0] + shit.Cx[0] * X + shit.Cy[0] * Y);
792 if (dp > Dcut) continue;
793 if (dp > Dmax) Dmax = dp;
794 n++;
795 w = 10. / (1.e-5 + fabs(d * d));
796 }
797
798 S0 += w * shit.S0[0];
799 S1 += w * shit.S1[0];
800 S2 += w * shit.S2[0];
801 S3 += w * shit.S3[0];
802 S4 += w * shit.S4[0];
803 S5 += w * shit.lx[0];
804 S6 += w * shit.ly[0];
805 S7 += w * shit.lr2[0];
806 } // ih
807
808 float s0 = S6 * S0 - S2 * S5;
809 float s1 = S0 * S1 - S2 * S2;
810 float s2 = S0 * S4 - S2 * S3;
811
812 if (fabs(s0) < 1.E-6 || fabs(s1) < 1.E-6) continue;
813 float tmp = s1 * (S5 * S5 - n * S0) + s0 * s0;
814 float A = ((S0 * S7 - S3 * S5) * s1 - s2 * s0) / tmp;
815 Y = (s2 + s0 * A) / s1 / 2;
816 X = (S3 + S5 * A - S2 * Y * 2) / S0 / 2;
817 float R2 = X * X + Y * Y - A;
818
819 if (R2 < 0) continue;
820 R = sqrt(R2);
821 if (Dmax <= 0) search_stop++;
822 } while (search_stop < 2.);
823
824 i->r = R;
825 i->x = X + firstHit.x[firstIh_4];
826 i->y = Y + firstHit.y[firstIh_4];
827
828 if (R < 2.5 || R > 7.5) {
829 i->on = 0;
830 i->skip = 1;
831 continue;
832 }
833
834 std::vector<THitIndex> newHits;
835 i->chi2 = 0;
836
837 for (THitIndex iih = 0; iih < maxI; iih++) {
838 const THitIndex ih = i->localIHits[iih];
839 const ENNHitV& hit = HitsV[ih / fvec::size()];
840 ENNSearchHitV& shit = shits[iih];
841 const int ih_4 = ih % fvec::size();
842
843 float dx = hit.x[ih_4] - i->x;
844 float dy = hit.y[ih_4] - i->y;
845 float d = fabs(sqrt(dx * dx + dy * dy) - i->r);
846 shit.on_ring[ih_4] = (d <= HitSize);
847 if (shit.on_ring[ih_4]) {
848 newHits.push_back(ih);
849 i->chi2 += d * d;
850 }
851 }
852
853 i->localIHits.clear();
854 i->localIHits = newHits;
855 i->NHits = i->localIHits.size();
856 i->NOwn = i->NHits;
857
858 if (i->localIHits.size() < MinRingHits) {
859 i->on = 0;
860 i->skip = 1;
861 continue;
862 }
863 i->chi2 = i->chi2 / ((i->localIHits.size() - 3) * HitSize * HitSize); //.3/.3;
864 }
865 }
866
867 sort(Rings.begin() + NInRings, Rings.end(), ENNRing::CompareENNHRings);
868 iR best;
869 for (iR i = Rbeg; i != Rend; ++i) {
870 i->on = 0;
871 if (i->skip) continue;
872 if (i->NOwn > 5) {
873 best = i;
874 const THitIndex maxI = i->localIHits.size();
875 for (THitIndex n = 0; n < maxI; n++) {
876 const THitIndex ih = i->localIHits[n];
877 ENNHitV& hit = HitsV[ih / fvec::size()];
878 const int ih_4 = ih % fvec::size();
879 hit.quality[ih_4] = 1;
880 }
881 for (iR j = i + 1; j != Rend; ++j) {
882 if (i->skip) continue;
883 float dist = j->r + best->r + HitSize2[0];
884 if (fabs(j->x - best->x) > dist || fabs(j->y - best->y) > dist) continue;
885 j->NOwn = 0;
886 const THitIndex maxJ = j->localIHits.size();
887 for (THitIndex m = 0; m < maxJ; m++) {
888 const THitIndex ihm = j->localIHits[m];
889 const ENNHitV& hitm = HitsV[ihm / fvec::size()];
890 if (hitm.quality[ihm % fvec::size()] == 0) j->NOwn++;
891 }
892 }
893 i->on = 1;
894 i->skip = 1;
895 }
896 else
897 i->skip = 1;
898 }
899#else // NEW_SELECTION
900
901 const int NHitsV = HitsV.size();
902 for (int ih = 0; ih < NHitsV; ih++)
903 HitsV[ih].quality = 0.;
904 for (iR ir = Rbeg; ir != Rend; ++ir) {
905 ir->on = 0;
906 ir->NOwn = ir->NHits;
907 ir->skip = ((ir->NHits < MinRingHits) || (ir->NHits <= 6 && ir->chi2 > .3));
908 }
909
910 do {
911 iR best = Rend;
912 THitIndex bestOwn = 0;
913 float bestChi2 = 1.E20;
914 for (iR ir = Rbeg; ir != Rend; ++ir) {
915 // const bool skip = ir->skip || ( ( ir->NOwn < 1.0*MinRingHits ) ||
916 // ( ir->NHits < 10 && ir->NOwn < 1.00*ir->NHits ) ||
917 // ( ir->NOwn < .60*ir->NHits ) );
918 if (ir->skip) continue; // faster with if
919 const bool skip = ((ir->NOwn < 1.0 * MinRingHits) || (ir->NHits < 10 && ir->NOwn < 1.00 * ir->NHits)
920 || (ir->NOwn < .60 * ir->NHits));
921 ir->skip = skip;
922 const bool isBetter = !skip & ((ir->NOwn > 1.2 * bestOwn) || (ir->NOwn >= 1. * bestOwn && ir->chi2 < bestChi2));
923 bestOwn = (isBetter) ? ir->NOwn : bestOwn; //Hits;
924 bestChi2 = (isBetter) ? ir->chi2 : bestChi2;
925 best = (isBetter) ? ir : best;
926 }
927 if (best == Rend) break;
928 best->skip = 1;
929 best->on = 1;
930 const THitIndex NHitsBest = best->localIHits.size();
931 for (THitIndex iih = 0; iih < NHitsBest; iih++) {
932 const THitIndex ih = best->localIHits[iih];
933 HitsV[ih / fvec::size()].quality[ih % fvec::size()] = 1;
934 }
935 for (iR ir = Rbeg; ir != Rend; ++ir) {
936 if (ir->skip) continue;
937 float dist = ir->r + best->r + HitSize2[0];
938 if (fabs(ir->x - best->x) > dist || fabs(ir->y - best->y) > dist) continue;
939 ir->NOwn = 0;
940 const THitIndex NHitsCur = ir->localIHits.size();
941 for (THitIndex iih = 0; iih < NHitsCur; iih++) {
942 const THitIndex ih = ir->localIHits[iih];
943 ir->NOwn += (HitsV[ih / fvec::size()].quality[ih % fvec::size()] == 0);
944 }
945 }
946 } while (1);
947#endif // else NEW_SELECTION
948
949#ifdef PRINT_TIMING
950 GetTimer("Selection").Stop();
951 GetTimer("All").Stop();
952#endif // PRINT_TIMING
953}
954
955TStopwatch& CbmL1RichENNRingFinderParallel::GetTimer(TString t)
956{
957 int i = 0;
958 for (; (i < NTimers) && (fTimersNames[i] != t); i++)
959 ;
960 assert(i < NTimers);
961 return fTimers[i];
962}
hits
static vector< vector< QAHit > > hits
Definition CbmTofHitFinderTBQa.cxx:108
sqrt
friend fvec sqrt(const fvec &a)
Definition KfSimdPseudo.h:118
iif
friend fvec iif(fmask a, fvec b, fvec c)
Definition KfSimdPseudo.h:100
fvec
fvec()
Definition KfSimdPseudo.h:6
CbmEvent
Class characterising one event by a collection of links (indices) to data objects,...
Definition CbmEvent.h:34
CbmL1RichENNRingFinderParallel::ENNRingFinder
void ENNRingFinder(const int NHits, cbm::algo::ca::Vector< ENNHitV > &HitsV, std::vector< ENNRing > &Rings, float HitSize=1., THitIndex MinRingHits=5, fvec RMin=2., fvec RMax=6.)
Definition CbmL1RichENNRingFinderParallel.cxx:241
CbmL1RichENNRingFinderParallel::DoFind
Int_t DoFind(CbmEvent *event, TClonesArray *hitArray, TClonesArray *projArray, TClonesArray *ringArray)
Definition CbmL1RichENNRingFinderParallel.cxx:79
CbmPixelHit::GetY
double GetY() const
Definition CbmPixelHit.h:74
CbmPixelHit::GetX
double GetX() const
Definition CbmPixelHit.h:73
CbmRichRing::SetChi2
void SetChi2(double chi2)
Definition CbmRichRing.h:61
CbmRichRing::AddHit
void AddHit(uint32_t pHit)
Definition CbmRichRing.h:31
CbmRichRing::SetCenterY
void SetCenterY(float y)
Definition CbmRichRing.h:53
CbmRichRing::SetCenterX
void SetCenterX(float x)
Definition CbmRichRing.h:52
CbmRichRing::SetRadius
void SetRadius(float r)
Definition CbmRichRing.h:54
cbm::algo::ca::Vector::push_back
void push_back(Tinput value)
Pushes back a value to the vector.
Definition CaVector.h:176
cbm::algo::ca::Vector::reserve
void reserve(std::size_t count)
Reserves a new size for the vector.
Definition CaVector.h:162
cbm::algo::ca::Vector::reset
void reset(std::size_t count, Tinput... value)
Clears vector and resizes it to the selected size with selected values.
Definition CaVector.h:121
cbm::algo::ca
TODO: SZh 8.11.2022: add selection of parameterisation.
Definition CaBranch.h:14
CbmL1RichENNRingFinderParallel::ENNHit::Compare
static bool Compare(const ENNHit &h1, const ENNHit &h2)
Definition CbmL1RichENNRingFinderParallel.h:52
CbmL1RichENNRingFinderParallel::ENNRing::CompareENNHRings
static bool CompareENNHRings(const ENNRing &r1, const ENNRing &r2)
Definition CbmL1RichENNRingFinderParallel.h:155