Belle II Software development
SingleHoughSpaceFastInterceptFinder.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#include <tracking/vxdHoughTracking/findlets/SingleHoughSpaceFastInterceptFinder.h>
9#include <tracking/vxdHoughTracking/entities/VXDHoughState.h>
10#include <tracking/spacePointCreation/SpacePoint.h>
11#include <tracking/spacePointCreation/SpacePointTrackCand.h>
12#include <tracking/trackFindingCDC/utilities/StringManipulation.h>
13#include <tracking/trackFindingCDC/utilities/Algorithms.h>
14#include <vxd/dataobjects/VxdID.h>
15#include <framework/core/ModuleParamList.h>
16#include <framework/core/ModuleParamList.templateDetails.h>
17
18using namespace Belle2;
19using namespace TrackFindingCDC;
20using namespace vxdHoughTracking;
21
25
26void SingleHoughSpaceFastInterceptFinder::exposeParameters(ModuleParamList* moduleParamList, const std::string& prefix)
27{
28 Super::exposeParameters(moduleParamList, prefix);
29
30 moduleParamList->addParameter(TrackFindingCDC::prefixed(prefix, "maximumRecursionLevel"), m_maxRecursionLevel,
31 "Maximum recursion level for the fast Hough trafo algorithm.", m_maxRecursionLevel);
32
33 moduleParamList->addParameter(TrackFindingCDC::prefixed(prefix, "nAngleSectors"), m_nAngleSectors,
34 "Number of angle sectors (= x-axis) dividing the Hough space.", m_nAngleSectors);
35
36 moduleParamList->addParameter(TrackFindingCDC::prefixed(prefix, "nVerticalSectors"), m_nVerticalSectors,
37 "Number of vertical sectors (= y-axis) dividing the Hough space.", m_nVerticalSectors);
38
39 moduleParamList->addParameter(TrackFindingCDC::prefixed(prefix, "verticalHoughSpaceSize"), m_verticalHoughSpaceSize,
40 "Vertical size of the Hough space.", m_verticalHoughSpaceSize);
41
42 moduleParamList->addParameter(TrackFindingCDC::prefixed(prefix, "HoughSpaceMinimumX"), m_minimumX,
43 "Minimum x value of the Hough space.", m_minimumX);
44
45 moduleParamList->addParameter(TrackFindingCDC::prefixed(prefix, "HoughSpaceMaximumX"), m_maximumX,
46 "Maximum x value of the Hough space.", m_maximumX);
47
48 moduleParamList->addParameter(TrackFindingCDC::prefixed(prefix, "minimumHSClusterSize"), m_MinimumHSClusterSize,
49 "Maximum x value of the Hough space.", m_MinimumHSClusterSize);
50
51 moduleParamList->addParameter(TrackFindingCDC::prefixed(prefix, "maximumHSClusterSize"), m_MaximumHSClusterSize,
52 "Maximum x value of the Hough space.", m_MaximumHSClusterSize);
53
54 moduleParamList->addParameter(TrackFindingCDC::prefixed(prefix, "maximumHSClusterSizeX"), m_MaximumHSClusterSizeX,
55 "Maximum x value of the Hough space.", m_MaximumHSClusterSizeX);
56
57 moduleParamList->addParameter(TrackFindingCDC::prefixed(prefix, "maximumHSClusterSizeY"), m_MaximumHSClusterSizeY,
58 "Maximum x value of the Hough space.", m_MaximumHSClusterSizeY);
59
60}
61
63{
65
66 const uint maxRecursionLevelFromSectors = ceil(log2(std::max(m_nAngleSectors, m_nVerticalSectors))) - 1;
67 m_maxRecursionLevel = std::max(maxRecursionLevelFromSectors, m_maxRecursionLevel);
68 if (m_maxRecursionLevel > 14) {
69 B2ERROR("The maximum number of recursions (maximumRecursionLevel) must not be larger than 14, but it is " <<
71 ", please choose a smaller value for maximumRecursionLevel, and / or for nAngleSectors and / or nVerticalSectors.");
72 }
74 for (uint i = 0; i < m_nAngleSectors; i++) {
75 double x = m_minimumX + m_unitX * (double)i;
76 double xc = x + 0.5 * m_unitX;
77
78 m_HSXLUT[i] = x;
79 m_HSSinValuesLUT[i] = sin(x);
80 m_HSCosValuesLUT[i] = cos(x);
81 m_HSCenterSinValuesLUT[i] = sin(xc);
82 m_HSCenterCosValuesLUT[i] = cos(xc);
83 m_HSXCenterLUT[i] = xc;
84 }
88
90 for (uint i = 0; i <= m_nVerticalSectors; i++) {
93 }
94 B2DEBUG(29, "HS size x: " << (m_maximumX - m_minimumX) << " HS size y: " << m_verticalHoughSpaceSize <<
95 " unitX: " << m_unitX << " unitY: " << m_unitY);
96}
97
98
99void SingleHoughSpaceFastInterceptFinder::apply(std::vector<VXDHoughState>& hits,
100 std::vector<std::vector<VXDHoughState*>>& rawTrackCandidates)
101{
102 m_trackCandidates.clear();
103 m_activeSectors.clear();
104
105 const std::vector<VXDHoughState*> currentEventHitList = TrackFindingCDC::as_pointers<VXDHoughState>(hits);
106
107 fastInterceptFinder2d(currentEventHitList, 0, m_nAngleSectors, 0, m_nVerticalSectors, 0);
108
110
111 for (auto& trackCand : m_trackCandidates) {
112 // sort for layer, and 2D radius in case of same layer before storing as SpacePointTrackCand
113 // outer hit goes first, as later on tracks are build from outside to inside
114 std::sort(trackCand.begin(), trackCand.end(),
115 [](const VXDHoughState * a, const VXDHoughState * b) {
116 return
117 (a->getDataCache().layer > b->getDataCache().layer) or
118 (a->getDataCache().layer == b->getDataCache().layer
119 and a->getHit()->getPosition().Perp() > b->getHit()->getPosition().Perp());
120 });
121
122 rawTrackCandidates.emplace_back(trackCand);
123 }
124
125 B2DEBUG(29, "m_trackCandidates.size: " << m_trackCandidates.size());
126
127}
128
129
130void SingleHoughSpaceFastInterceptFinder::fastInterceptFinder2d(const std::vector<VXDHoughState*>& hits, uint xmin, uint xmax,
131 uint ymin,
132 uint ymax, uint currentRecursion)
133{
134 std::vector<VXDHoughState*> containedHits;
135 containedHits.reserve(hits.size());
136 std::bitset<8> layerHits; /* For layer filter */
137
138 if (currentRecursion == m_maxRecursionLevel + 1) return;
139
140 // these int-divisions can cause {min, center} or {center, max} to be the same, which is a desired behaviour
141 const uint centerx = xmin + (uint)((xmax - xmin) / 2);
142 const uint centery = ymin + (uint)((ymax - ymin) / 2);
143 const uint xIndexCache[3] = {xmin, centerx, xmax};
144 const uint yIndexCache[3] = {ymin, centery, ymax};
145
146 for (int i = 0; i < 2 ; ++i) {
147 const uint left = xIndexCache[i];
148 const uint right = xIndexCache[i + 1];
149 const uint localIndexX = left;
150
151 if (left == right) continue;
152
153 const double& sinLeft = m_HSSinValuesLUT[left];
154 const double& cosLeft = m_HSCosValuesLUT[left];
155 const double& sinRight = m_HSSinValuesLUT[right];
156 const double& cosRight = m_HSCosValuesLUT[right];
157
158 // the sin and cos of the current center can't be stored in a LUT, as the number of possible centers
159 // is quite large and the logic would become rather complex
160 const double sinCenter = m_HSCenterSinValuesLUT[(left + right) / 2];
161 const double cosCenter = m_HSCenterCosValuesLUT[(left + right) / 2];
162
163 for (int j = 0; j < 2; ++j) {
164 const uint lowerIndex = yIndexCache[j];
165 const uint upperIndex = yIndexCache[j + 1];
166
167 if (lowerIndex == upperIndex) continue;
168
169 const uint localIndexY = lowerIndex;
170 const double& localUpperCoordinate = m_HSYLUT[lowerIndex];
171 const double& localLowerCoordinate = m_HSYLUT[upperIndex];
172
173 // reset layerHits and containedHits
174 layerHits = 0;
175 containedHits.clear();
176 for (VXDHoughState* hit : hits) {
177
178 const VXDHoughState::DataCache& hitData = hit->getDataCache();
179 const double& m = hitData.xConformal;
180 const double& a = hitData.yConformal;
181
182 const double derivativeyLeft = m * -sinLeft + a * cosLeft;
183 const double derivativeyRight = m * -sinRight + a * cosRight;
184 const double derivativeyCenter = m * -sinCenter + a * cosCenter;
185
186 // Only interested in the rising arm of the sinosoidal curves.
187 // Thus if derivative on both sides of the cell is negative, ignore and continue.
188 if (derivativeyLeft < 0 and derivativeyRight < 0 and derivativeyCenter < 0) continue;
189
190 const double yLeft = m * cosLeft + a * sinLeft;
191 const double yRight = m * cosRight + a * sinRight;
192 const double yCenter = m * cosCenter + a * sinCenter;
193
194 /* Check if HS-parameter curve is inside (or outside) actual sub-HS */
195 if ((yLeft <= localUpperCoordinate and yRight >= localLowerCoordinate) or
196 (yCenter <= localUpperCoordinate and yLeft >= localLowerCoordinate and yRight >= localLowerCoordinate) or
197 (yCenter >= localLowerCoordinate and yLeft <= localUpperCoordinate and yRight <= localUpperCoordinate)) {
198 layerHits[hitData.layer] = true;
199 containedHits.emplace_back(hit);
200 }
201 }
202
203 if (layerFilter(layerHits) > 0) {
204 // recursive call of fastInterceptFinder2d, until currentRecursion == m_maxRecursionLevel
205 if (currentRecursion < m_maxRecursionLevel) {
206 fastInterceptFinder2d(containedHits, left, right, lowerIndex, upperIndex, currentRecursion + 1);
207 } else {
208 m_activeSectors.insert({std::make_pair(localIndexX, localIndexY), std::make_pair(-layerFilter(layerHits), containedHits) });
209 }
210 }
211 }
212 }
213}
214
215
217{
218 m_clusterCount = 1;
219
220 for (auto& currentCell : m_activeSectors) {
221
222 // cell content meanings:
223 // -3, -4 : active sector, not yet visited
224 // 0 : non-active sector (will never be visited, only checked)
225 // 1,2,3...: index of the clusters
226 if (currentCell.second.first > -1) continue;
227
228 m_clusterInitialPosition = std::make_pair(currentCell.first.first, currentCell.first.second);
229 m_clusterSize = 1;
230 currentCell.second.first = m_clusterCount;
231
233 for (VXDHoughState* hit : currentCell.second.second) {
234 m_currentTrackCandidate.emplace_back(hit);
235 }
236
237 // Check for HS sectors connected to each other which could form a cluster
238 DepthFirstSearch(currentCell.first.first, currentCell.first.second);
239 // if cluster valid (i.e. not too small and not too big): finalize!
241
244 }
246 }
247}
248
249void SingleHoughSpaceFastInterceptFinder::DepthFirstSearch(uint lastIndexX, uint lastIndexY)
250{
252
253 for (uint currentIndexY = lastIndexY; currentIndexY >= lastIndexY - 1; currentIndexY--) {
254 if (abs((int)m_clusterInitialPosition.second - (int)currentIndexY) >= m_MaximumHSClusterSizeY or
255 m_clusterSize >= m_MaximumHSClusterSize or currentIndexY > m_nVerticalSectors) return;
256 for (uint currentIndexX = lastIndexX; currentIndexX <= lastIndexX + 1; currentIndexX++) {
257 if (abs((int)m_clusterInitialPosition.first - (int)currentIndexX) >= m_MaximumHSClusterSizeX or
258 m_clusterSize >= m_MaximumHSClusterSize or currentIndexX > m_nAngleSectors) return;
259
260 // The cell (currentIndexX, currentIndexY) is the current one has already been checked, so continue
261 if (lastIndexX == currentIndexX && lastIndexY == currentIndexY) continue;
262
263 // first check bounds to avoid out-of-bound array access
264 // as they are uints, they are always >= 0, and in case of an overflow they would be too large
265 if (currentIndexX < m_nAngleSectors and currentIndexY < m_nVerticalSectors) {
266
267 auto activeSector = m_activeSectors.find({currentIndexX, currentIndexY});
268 // Only continue searching if the current cluster is smaller than the maximum cluster size
269 if (activeSector != m_activeSectors.end() and activeSector->second.first < 0 /*and m_clusterSize < m_MaximumHSClusterSize*/) {
270 activeSector->second.first = m_clusterCount;
272
273 // No need to check whether currentIndex exists as a key in m_activeSectors as they were created at the same time
274 // so it's certain the key exists.
275 for (VXDHoughState* hit : activeSector->second.second) {
276 if (not TrackFindingCDC::is_in(hit, m_currentTrackCandidate)) {
277 m_currentTrackCandidate.emplace_back(hit);
278 }
279 }
280
281 // search in the next Hough Space cells...
282 DepthFirstSearch(currentIndexX, currentIndexY);
283 }
284 }
285 }
286 }
287}
The Module parameter list class.
virtual void exposeParameters(ModuleParamList *moduleParamList, const std::string &prefix)
std::array< double, 16384 > m_HSCenterSinValuesLUT
sine values of the Hough Space sector center coordinates
uint m_nVerticalSectors
number of sectors of the Hough Space on the vertical axis
uint m_MinimumHSClusterSize
minimum cluster size of sectors belonging to intercepts in the Hough Space
std::array< double, 16384 > m_HSYCenterLUT
y values of the Hough Space sector centers
std::array< double, 16385 > m_HSYLUT
y values of the Hough Space sector boarders
std::array< double, 16384 > m_HSCenterCosValuesLUT
cosine values of the Hough Space sector center coordinates
std::map< std::pair< uint, uint >, std::pair< int, std::vector< VXDHoughState * > >, paircompare > m_activeSectors
Map containing only active HS sectors, i.e.
uint m_nAngleSectors
number of sectors of the Hough Space on the horizontal axis
std::array< double, 16385 > m_HSXLUT
x values of the Hough Space sector boarders
double m_minimumX
minimum x value of the Hough Space, defaults to the value for u-side
std::vector< std::vector< VXDHoughState * > > m_trackCandidates
vector containing track candidates, consisting of the found intersection values in the Hough Space
double m_verticalHoughSpaceSize
vertical size of the Hough Space, defaults to the value for u-side
unsigned short layerFilter(const std::bitset< 8 > &layer)
layer filter, checks if at least hits from 3 layers are in a set of hits
TrackFindingCDC::Findlet< VXDHoughState, std::vector< VXDHoughState * > > Super
Parent class.
void DepthFirstSearch(uint lastIndexX, uint lastIndexY)
Perform depth first search recursive algorithm to find clusters in the Hough Space.
std::pair< uint, uint > m_clusterInitialPosition
start cell of the recursive cluster finding in the Hough Space
uint m_MaximumHSClusterSizeX
maximum cluster size in x of sectors belonging to intercepts in the Hough Space
void fastInterceptFinder2d(const std::vector< VXDHoughState * > &hits, uint xmin, uint xmax, uint ymin, uint ymax, uint currentRecursion)
find intercepts in the 2D Hough Space by recursively calling itself until no hits are assigned to a g...
std::array< double, 16385 > m_HSSinValuesLUT
Look-Up-Tables for values as cache to speed up calculation sine values of the Hough Space sector boar...
std::vector< VXDHoughState * > m_currentTrackCandidate
the current track candidate
double m_maximumX
maximum x value of the Hough Space, defaults to the value for u-side
uint m_MaximumHSClusterSizeY
maximum cluster size in y of sectors belonging to intercepts in the Hough Space
void exposeParameters(ModuleParamList *moduleParamList, const std::string &prefix) override
Expose the parameters of the sub findlets.
uint m_MaximumHSClusterSize
maximum cluster size of sectors belonging to intercepts in the Hough Space
void apply(std::vector< VXDHoughState > &hits, std::vector< std::vector< VXDHoughState * > > &rawTrackCandidates) override
Load in the prepared hits and create track candidates for further processing like hit filtering and f...
std::array< double, 16385 > m_HSCosValuesLUT
cosine values of the Hough Space sector boarder coordinates
std::array< double, 16384 > m_HSXCenterLUT
x values of the Hough Space sector centers
uint m_maxRecursionLevel
maximum number of recursive calls of FastInterceptFinder2d
Simple container for hit information to be used during intercept finding.
void addParameter(const std::string &name, T &paramVariable, const std::string &description, const T &defaultValue)
Adds a new parameter to the module list.
Abstract base class for different kinds of events.
Cache containing the most important information of this state which will often be needed.
float yConformal
conformal transformed y coordinate of this hit
unsigned short layer
Geometrical Layer this state is based on.
float xConformal
conformal transformed x coordinate of this hit