CbmRoot
CbmRichRingFinderHoughSimd.cxx
Go to the documentation of this file.
1 
8 //#include "../L1/L1Algo/L1Types.h"
10 #include "../L1/L1Algo/vectors/P4_F32vec4.h"
11 #include <emmintrin.h>
12 
14 
15 void CbmRichRingFinderHoughSimd::HoughTransformGroup(unsigned short int indmin,
16  unsigned short int indmax,
17  Int_t iPart) {
18  // Float_t r12, r13, r23;
19  // Float_t rx0, rx1, rx2, ry0, ry1,ry2; //rx[3], ry[3];//, x[3], y[3];
20  //Float_t xc, yc, r;
21  //Float_t xcs, ycs; // xcs = xc - fCurMinX
22  Int_t *intX, *intY, *intR;
23  Int_t* indXY;
24 
25  unsigned short int iH11, iH12, iH13, iH14, iH2, iH3;
26  Int_t nofHitsNorm = fHitInd[0].size() + 1;
27  Int_t iPmulNofHits;
28 
29  //Float_t t5, t10, t19, det, t6, t7;
30  Float_t dx = 1.0f / fDx, dy = 1.0f / fDy, dr = 1.0f / fDr;
31  //Float_t iH1X, iH1Y, iH2X, iH2Y, iH3X, iH3Y;
32 
33  fvec xcs, ycs;
34  fvec fCurMinXV = fvec(fCurMinX), fCurMinYV = fvec(fCurMinY);
35  fvec xc, yc, r;
36  fvec r12, r13, r23;
37  fvec rx0, rx1, rx2, ry0, ry1, ry2; //rx[3], ry[3];//, x[3], y[3];
38  fvec t5, t10, t19, det, t6, t7;
39  fvec iH1X, iH1Y, iH2X, iH2Y, iH3X, iH3Y;
40  // fvec intXfv, intYfv, intRfv;
41 
42  __m128i intXv, intYv, intRv, indXYv;
43  __m128i fNofBinsXv = _mm_set1_epi32(fNofBinsX);
44  __m128i fNofBinsXYv = _mm_set1_epi32(fNofBinsXY + 1);
45  __m128i m128i_0 = _mm_set1_epi32(0);
46  fvec fMinDistanceSqv = fvec(fMinDistanceSq);
47  fvec fMaxDistanceSqv = fvec(fMaxDistanceSq);
48  __m128i fNofBinsRv = _mm_set1_epi32(fNofBinsR + 1);
49  __m128i m128i_10 = _mm_set1_epi32(10);
50 
51  fvec dxv = fvec(1.0f / fDx), dyv = fvec(1.0f / fDy), drv = fvec(1.0f / fDr);
52  fvec fvec05 = fvec(0.5f);
53  fvec fvec0 = fvec(0.f);
54  Int_t nofHits = fHitInd[iPart].size();
55  if (nofHits <= fMinNofHitsInArea) return;
56  iPmulNofHits = iPart * nofHitsNorm;
57 
58  vector<unsigned short> hitIndPart;
59  hitIndPart.assign(fHitInd[iPart].begin(), fHitInd[iPart].end());
60 
61  CbmRichHoughHitVec h2, h3;
62 
63  for (unsigned short int iHit1 = 0; iHit1 < (nofHits & ~3); iHit1 += 4) {
64  iH11 = hitIndPart[iHit1];
65  iH12 = hitIndPart[iHit1 + 1];
66  iH13 = hitIndPart[iHit1 + 2];
67  iH14 = hitIndPart[iHit1 + 3];
68 
69  iH1X =
70  fvec(fData[iH11]->fX, fData[iH12]->fX, fData[iH13]->fX, fData[iH14]->fX);
71  iH1Y =
72  fvec(fData[iH11]->fY, fData[iH12]->fY, fData[iH13]->fY, fData[iH14]->fY);
73 
74  for (unsigned short int iHit2 = iHit1 + 1; iHit2 < nofHits; iHit2++) {
75  iH2 = hitIndPart[iHit2];
76  h2 = fDataV[iH2];
77  iH2X = h2.fX;
78  iH2Y = h2.fY;
79 
80  rx0 = iH1X - iH2X; //rx12
81  ry0 = iH1Y - iH2Y; //ry12
82  r12 = rx0 * rx0 + ry0 * ry0;
83  if (_mm_comineq_ss(_mm_cmpgt_ss(r12, fMinDistanceSqv), fvec0) == 1)
84  continue;
85  if (_mm_comineq_ss(_mm_cmplt_ss(r12, fMaxDistanceSqv), fvec0) == 1)
86  continue;
87 
88 
89  t10 = fvec(fData[iH11]->fX2plusY2 - fData[iH2]->fX2plusY2,
90  fData[iH12]->fX2plusY2 - fData[iH2]->fX2plusY2,
91  fData[iH13]->fX2plusY2 - fData[iH2]->fX2plusY2,
92  fData[iH14]->fX2plusY2 - fData[iH2]->fX2plusY2);
93  for (unsigned short int iHit3 = iHit2 + 1; iHit3 < nofHits; iHit3++) {
94  //iH3 = hitIndPart[iHit3];
95  h3 = fDataV[hitIndPart[iHit3]];
96  // iH3X = h3.fX;
97  // iH3Y = h3.fY;
98  t5 = h2.fX2plusY2 - h3.fX2plusY2;
99 
100  rx1 = iH1X - iH3X; //rx13
101  ry1 = iH1Y - iH3Y; //ry13
102  r13 = rx1 * rx1 + ry1 * ry1;
103  if (_mm_comineq_ss(_mm_cmpgt_ss(r13, fMinDistanceSqv), fvec0) == 1)
104  continue;
105  if (_mm_comineq_ss(_mm_cmplt_ss(r13, fMaxDistanceSqv), fvec0) == 1)
106  continue;
107 
108  rx2 = iH2X - h3.fX; //rx23
109  ry2 = iH2Y - h3.fY; //ry23
110  r23 = rx2 * rx2 + ry2 * ry2;
111  if (_mm_comineq_ss(_mm_cmpgt_ss(r23, fMinDistanceSqv), fvec0) == 1)
112  continue;
113  if (_mm_comineq_ss(_mm_cmplt_ss(r23, fMaxDistanceSqv), fvec0) == 1)
114  continue;
115 
116  t19 = fvec05 / (rx2 * ry0 - rx0 * ry2);
117 
118  xc = (t5 * ry0 - t10 * ry2) * t19;
119  xcs = xc - fCurMinXV;
120  //intX = int( xcs *dx);
121  //if (intX < 0 || intX >= fNofBinsX ) continue;
122 
123  yc = (t10 * rx2 - t5 * rx0) * t19;
124  ycs = yc - fCurMinYV;
125  //intY = int( ycs *dy);
126  //if (intY < 0 || intY >= fNofBinsY ) continue;
127 
128  //radius calculation
129  t6 = iH1X - xc;
130  t7 = iH1Y - yc;
131  r = sqrt(t6 * t6 + t7 * t7);
132 
133  //intR = int(r *dr);
134  //if (intR < 0 || intR > fNofBinsR) continue;
135  //indXY = intX * fNofBinsX + intY;
136  // intXfv = xcs *dxv;
137  // intYfv = ycs *dyv;
138  // intRfv = r *drv;
139  intRv = _mm_cvtps_epi32(r * drv);
140  // if (_mm_movemask_epi8( _mm_cmpgt_epi32(intRv, m128i_10) ) == 0x0000)continue;
141  // if (_mm_movemask_epi8( _mm_cmplt_epi32(intRv, fNofBinsRv) ) == 0x0000) continue;
142 
143  intXv = _mm_cvtps_epi32(xcs * dxv);
144  intYv = _mm_cvtps_epi32(ycs * dyv);
145  indXYv = intXv * fNofBinsXv + intYv;
146 
147  // if (_mm_movemask_epi8( _mm_cmpgt_epi32(indXYv, m128i_0) ) == 0x0000) continue;
148  // if (_mm_movemask_epi8( _mm_cmplt_epi32(indXYv, fNofBinsXYv) ) == 0x0000) continue;
149 
150 
151  indXY = (int*) &indXYv;
152  intR = (int*) &intRv;
153 
154  if (intR[0] > 9 && intR[0] < fNofBinsR && indXY[0] > -1
155  && indXY[0] < fNofBinsXY) {
156  fHist[indXY[0]]++;
157  fHistR[intR[0]]++;
158  }
159  if (intR[1] > 9 && intR[1] < fNofBinsR && indXY[1] > -1
160  && indXY[1] < fNofBinsXY) {
161  fHist[indXY[1]]++;
162  fHistR[intR[1]]++;
163  }
164 
165  if (intR[2] > 9 && intR[2] < fNofBinsR && indXY[2] > -1
166  && indXY[2] < fNofBinsXY) {
167  fHist[indXY[2]]++;
168  fHistR[intR[2]]++;
169  }
170 
171  if (intR[3] > 9 && intR[3] < fNofBinsR && indXY[3] > -1
172  && indXY[3] < fNofBinsXY) {
173  fHist[indXY[3]]++;
174  fHistR[intR[3]]++;
175  }
176  } //iHit1
177  } //iHit2
178  } //iHit3
179 }
180 
181 
183  Int_t indmin, indmax;
184  Float_t x0, y0;
185 
186  fDataV.clear();
187  fDataV.reserve(fData.size());
188  for (UInt_t iHit = 0; iHit < fData.size(); iHit++) {
190  h.fX = fvec(fData[iHit].fX);
191  h.fY = fvec(fData[iHit].fY);
192  h.fX2plusY2 = fvec(fData[iHit].fX2plusY2);
193  fDataV.push_back(h);
194  }
195 
196  for (UInt_t iHit = 0; iHit < fData.size(); iHit++) {
197  if (fData[iHit].fIsUsed == true) continue;
198 
199  x0 = fData[iHit].fX;
200  y0 = fData[iHit].fY;
201 
202  fCurMinX = x0 - fMaxDistance;
203  fCurMinY = y0 - fMaxDistance;
204 
205  DefineLocalAreaAndHits(x0, y0, &indmin, &indmax);
206  HoughTransform(indmin, indmax);
207  FindPeak(indmin, indmax);
208 
209  } //main loop over hits
210 }
CbmRichHoughHitVec::fY
fvec fY
Definition: CbmRichRingFinderHoughSimd.h:20
CbmRichRingFinderHoughImpl::fHitInd
vector< vector< unsigned int > > fHitInd
Definition: CbmRichRingFinderHoughImpl.h:85
f
float f
Definition: L1/vectors/P4_F32vec4.h:24
F32vec4
Definition: L1/vectors/P4_F32vec4.h:47
CbmRichRingFinderHoughImpl::fCurMinY
float fCurMinY
Definition: CbmRichRingFinderHoughImpl.h:77
sqrt
friend F32vec4 sqrt(const F32vec4 &a)
Definition: L1/vectors/P4_F32vec4.h:41
fvec
F32vec4 fvec
Definition: L1/vectors/P4_F32vec4.h:249
CbmRichRingFinderHoughImpl::HoughTransform
virtual void HoughTransform(unsigned int indmin, unsigned int indmax)
Run HoughTransformGroup for each group of hits.
Definition: CbmRichRingFinderHoughImpl.cxx:228
CbmRichRingFinderHoughSimd.h
SIMDized ring finder based on Hough Transform method.
CbmRichRingFinderHoughImpl::fDy
float fDy
Definition: CbmRichRingFinderHoughImpl.h:51
CbmRichRingFinderHoughSimd::HoughTransformReconstruction
virtual void HoughTransformReconstruction()
Run HT for each hit.
Definition: CbmRichRingFinderHoughSimd.cxx:182
CbmRichRingFinderHoughImpl::fDr
float fDr
Definition: CbmRichRingFinderHoughImpl.h:52
CbmRichRingFinderHoughSimd::HoughTransformGroup
virtual void HoughTransformGroup(unsigned short int indmin, unsigned short int indmax, Int_t iPart)
Definition: CbmRichRingFinderHoughSimd.cxx:15
CbmRichRingFinderHoughImpl::fMaxDistance
float fMaxDistance
Definition: CbmRichRingFinderHoughImpl.h:42
CbmRichRingFinderHoughImpl::fMinNofHitsInArea
unsigned short fMinNofHitsInArea
Definition: CbmRichRingFinderHoughImpl.h:64
CbmRichRingFinderHoughImpl::fNofBinsXY
unsigned short fNofBinsXY
Definition: CbmRichRingFinderHoughImpl.h:55
CbmRichRingFinderHoughImpl::DefineLocalAreaAndHits
virtual void DefineLocalAreaAndHits(float x0, float y0, int *indmin, int *indmax)
Find hits in a local area.
Definition: CbmRichRingFinderHoughImpl.cxx:179
h
Data class with information on a STS local track.
CbmRichRingFinderHoughImpl::fHist
vector< unsigned short > fHist
Definition: CbmRichRingFinderHoughImpl.h:82
CbmRichRingFinderHoughImpl::FindPeak
void FindPeak(int indmin, int indmax)
Definition: CbmRichRingFinderHoughImpl.cxx:318
CbmRichRingFinderHoughImpl::fMinDistanceSq
float fMinDistanceSq
Definition: CbmRichRingFinderHoughImpl.h:44
CbmRichRingFinderHoughImpl::fHistR
vector< unsigned short > fHistR
Definition: CbmRichRingFinderHoughImpl.h:83
CbmRichRingFinderHoughImpl::fNofBinsR
unsigned short fNofBinsR
Definition: CbmRichRingFinderHoughImpl.h:60
CbmRichRingFinderHoughImpl::fCurMinX
float fCurMinX
Definition: CbmRichRingFinderHoughImpl.h:76
CbmRichHoughHitVec::fX2plusY2
fvec fX2plusY2
Definition: CbmRichRingFinderHoughSimd.h:21
CbmRichRingFinderHoughSimd::fDataV
std::vector< CbmRichHoughHitVec > fDataV
Definition: CbmRichRingFinderHoughSimd.h:47
CbmRichRingFinderHoughImpl::fDx
float fDx
Definition: CbmRichRingFinderHoughImpl.h:50
CbmRichRingFinderHoughSimd::CbmRichRingFinderHoughSimd
CbmRichRingFinderHoughSimd()
Definition: CbmRichRingFinderHoughSimd.cxx:13
CbmRichHoughHitVec::fX
fvec fX
Definition: CbmRichRingFinderHoughSimd.h:19
CbmRichRingFinderHoughImpl::fNofBinsX
unsigned short fNofBinsX
Definition: CbmRichRingFinderHoughImpl.h:53
CbmRichHoughHitVec
Definition: CbmRichRingFinderHoughSimd.h:17
CbmRichRingFinderHoughImpl::fData
vector< CbmRichHoughHit > fData
Definition: CbmRichRingFinderHoughImpl.h:81
CbmRichRingFinderHoughImpl::fMaxDistanceSq
float fMaxDistanceSq
Definition: CbmRichRingFinderHoughImpl.h:45