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
17using namespace Belle2;
18using namespace TrackFindingCDC;
19using namespace vxdHoughTracking;
20
24
25void SingleHoughSpaceFastInterceptFinder::exposeParameters(ModuleParamList* moduleParamList, const std::string& prefix)
26{
27 Super::exposeParameters(moduleParamList, prefix);
28
29 moduleParamList->addParameter(TrackFindingCDC::prefixed(prefix, "maximumRecursionLevel"), m_maxRecursionLevel,
30 "Maximum recursion level for the fast Hough trafo algorithm.", m_maxRecursionLevel);
31
32 moduleParamList->addParameter(TrackFindingCDC::prefixed(prefix, "nAngleSectors"), m_nAngleSectors,
33 "Number of angle sectors (= x-axis) dividing the Hough space.", m_nAngleSectors);
34
35 moduleParamList->addParameter(TrackFindingCDC::prefixed(prefix, "nVerticalSectors"), m_nVerticalSectors,
36 "Number of vertical sectors (= y-axis) dividing the Hough space.", m_nVerticalSectors);
37
38 moduleParamList->addParameter(TrackFindingCDC::prefixed(prefix, "verticalHoughSpaceSize"), m_verticalHoughSpaceSize,
39 "Vertical size of the Hough space.", m_verticalHoughSpaceSize);
40
41 moduleParamList->addParameter(TrackFindingCDC::prefixed(prefix, "HoughSpaceMinimumX"), m_minimumX,
42 "Minimum x value of the Hough space.", m_minimumX);
43
44 moduleParamList->addParameter(TrackFindingCDC::prefixed(prefix, "HoughSpaceMaximumX"), m_maximumX,
45 "Maximum x value of the Hough space.", m_maximumX);
46
47 moduleParamList->addParameter(TrackFindingCDC::prefixed(prefix, "minimumHSClusterSize"), m_MinimumHSClusterSize,
48 "Maximum x value of the Hough space.", m_MinimumHSClusterSize);
49
50 moduleParamList->addParameter(TrackFindingCDC::prefixed(prefix, "maximumHSClusterSize"), m_MaximumHSClusterSize,
51 "Maximum x value of the Hough space.", m_MaximumHSClusterSize);
52
53 moduleParamList->addParameter(TrackFindingCDC::prefixed(prefix, "maximumHSClusterSizeX"), m_MaximumHSClusterSizeX,
54 "Maximum x value of the Hough space.", m_MaximumHSClusterSizeX);
55
56 moduleParamList->addParameter(TrackFindingCDC::prefixed(prefix, "maximumHSClusterSizeY"), m_MaximumHSClusterSizeY,
57 "Maximum x value of the Hough space.", m_MaximumHSClusterSizeY);
58
59}
60
62{
64
65 const uint maxRecursionLevelFromSectors = ceil(log2(std::max(m_nAngleSectors, m_nVerticalSectors))) - 1;
66 m_maxRecursionLevel = std::max(maxRecursionLevelFromSectors, m_maxRecursionLevel);
67 if (m_maxRecursionLevel > 14) {
68 B2ERROR("The maximum number of recursions (maximumRecursionLevel) must not be larger than 14, but it is " <<
70 ", please choose a smaller value for maximumRecursionLevel, and / or for nAngleSectors and / or nVerticalSectors.");
71 }
73 for (uint i = 0; i < m_nAngleSectors; i++) {
74 double x = m_minimumX + m_unitX * (double)i;
75 double xc = x + 0.5 * m_unitX;
76
77 m_HSXLUT[i] = x;
78 m_HSSinValuesLUT[i] = sin(x);
79 m_HSCosValuesLUT[i] = cos(x);
80 m_HSCenterSinValuesLUT[i] = sin(xc);
81 m_HSCenterCosValuesLUT[i] = cos(xc);
82 m_HSXCenterLUT[i] = xc;
83 }
87
89 for (uint i = 0; i <= m_nVerticalSectors; i++) {
92 }
93 B2DEBUG(29, "HS size x: " << (m_maximumX - m_minimumX) << " HS size y: " << m_verticalHoughSpaceSize <<
94 " unitX: " << m_unitX << " unitY: " << m_unitY);
95}
96
97
98void SingleHoughSpaceFastInterceptFinder::apply(std::vector<VXDHoughState>& hits,
99 std::vector<std::vector<VXDHoughState*>>& rawTrackCandidates)
100{
101 m_trackCandidates.clear();
102 m_activeSectors.clear();
103
104 const std::vector<VXDHoughState*> currentEventHitList = TrackFindingCDC::as_pointers<VXDHoughState>(hits);
105
106 fastInterceptFinder2d(currentEventHitList, 0, m_nAngleSectors, 0, m_nVerticalSectors, 0);
107
109
110 for (auto& trackCand : m_trackCandidates) {
111 // sort for layer, and 2D radius in case of same layer before storing as SpacePointTrackCand
112 // outer hit goes first, as later on tracks are build from outside to inside
113 std::sort(trackCand.begin(), trackCand.end(),
114 [](const VXDHoughState * a, const VXDHoughState * b) {
115 return
116 (a->getDataCache().layer > b->getDataCache().layer) or
117 (a->getDataCache().layer == b->getDataCache().layer
118 and a->getHit()->getPosition().Perp() > b->getHit()->getPosition().Perp());
119 });
120
121 rawTrackCandidates.emplace_back(trackCand);
122 }
123
124 B2DEBUG(29, "m_trackCandidates.size: " << m_trackCandidates.size());
125
126}
127
128
129void SingleHoughSpaceFastInterceptFinder::fastInterceptFinder2d(const std::vector<VXDHoughState*>& hits, uint xmin, uint xmax,
130 uint ymin,
131 uint ymax, uint currentRecursion)
132{
133 std::vector<VXDHoughState*> containedHits;
134 containedHits.reserve(hits.size());
135 std::bitset<8> layerHits; /* For layer filter */
136
137 if (currentRecursion == m_maxRecursionLevel + 1) return;
138
139 // these int-divisions can cause {min, center} or {center, max} to be the same, which is a desired behaviour
140 const uint centerx = xmin + (uint)((xmax - xmin) / 2);
141 const uint centery = ymin + (uint)((ymax - ymin) / 2);
142 const uint xIndexCache[3] = {xmin, centerx, xmax};
143 const uint yIndexCache[3] = {ymin, centery, ymax};
144
145 for (int i = 0; i < 2 ; ++i) {
146 const uint left = xIndexCache[i];
147 const uint right = xIndexCache[i + 1];
148 const uint localIndexX = left;
149
150 if (left == right) continue;
151
152 const double& sinLeft = m_HSSinValuesLUT[left];
153 const double& cosLeft = m_HSCosValuesLUT[left];
154 const double& sinRight = m_HSSinValuesLUT[right];
155 const double& cosRight = m_HSCosValuesLUT[right];
156
157 // the sin and cos of the current center can't be stored in a LUT, as the number of possible centers
158 // is quite large and the logic would become rather complex
159 const double sinCenter = m_HSCenterSinValuesLUT[(left + right) / 2];
160 const double cosCenter = m_HSCenterCosValuesLUT[(left + right) / 2];
161
162 for (int j = 0; j < 2; ++j) {
163 const uint lowerIndex = yIndexCache[j];
164 const uint upperIndex = yIndexCache[j + 1];
165
166 if (lowerIndex == upperIndex) continue;
167
168 const uint localIndexY = lowerIndex;
169 const double& localUpperCoordinate = m_HSYLUT[lowerIndex];
170 const double& localLowerCoordinate = m_HSYLUT[upperIndex];
171
172 // reset layerHits and containedHits
173 layerHits = 0;
174 containedHits.clear();
175 for (VXDHoughState* hit : hits) {
176
177 const VXDHoughState::DataCache& hitData = hit->getDataCache();
178 const double& m = hitData.xConformal;
179 const double& a = hitData.yConformal;
180
181 const double derivativeyLeft = m * -sinLeft + a * cosLeft;
182 const double derivativeyRight = m * -sinRight + a * cosRight;
183 const double derivativeyCenter = m * -sinCenter + a * cosCenter;
184
185 // Only interested in the rising arm of the sinosoidal curves.
186 // Thus if derivative on both sides of the cell is negative, ignore and continue.
187 if (derivativeyLeft < 0 and derivativeyRight < 0 and derivativeyCenter < 0) continue;
188
189 const double yLeft = m * cosLeft + a * sinLeft;
190 const double yRight = m * cosRight + a * sinRight;
191 const double yCenter = m * cosCenter + a * sinCenter;
192
193 /* Check if HS-parameter curve is inside (or outside) actual sub-HS */
194 if ((yLeft <= localUpperCoordinate and yRight >= localLowerCoordinate) or
195 (yCenter <= localUpperCoordinate and yLeft >= localLowerCoordinate and yRight >= localLowerCoordinate) or
196 (yCenter >= localLowerCoordinate and yLeft <= localUpperCoordinate and yRight <= localUpperCoordinate)) {
197 layerHits[hitData.layer] = true;
198 containedHits.emplace_back(hit);
199 }
200 }
201
202 if (layerFilter(layerHits) > 0) {
203 // recursive call of fastInterceptFinder2d, until currentRecursion == m_maxRecursionLevel
204 if (currentRecursion < m_maxRecursionLevel) {
205 fastInterceptFinder2d(containedHits, left, right, lowerIndex, upperIndex, currentRecursion + 1);
206 } else {
207 m_activeSectors.insert({std::make_pair(localIndexX, localIndexY), std::make_pair(-layerFilter(layerHits), containedHits) });
208 }
209 }
210 }
211 }
212}
213
214
216{
217 m_clusterCount = 1;
218
219 for (auto& currentCell : m_activeSectors) {
220
221 // cell content meanings:
222 // -3, -4 : active sector, not yet visited
223 // 0 : non-active sector (will never be visited, only checked)
224 // 1,2,3...: index of the clusters
225 if (currentCell.second.first > -1) continue;
226
227 m_clusterInitialPosition = std::make_pair(currentCell.first.first, currentCell.first.second);
228 m_clusterSize = 1;
229 currentCell.second.first = m_clusterCount;
230
232 for (VXDHoughState* hit : currentCell.second.second) {
233 m_currentTrackCandidate.emplace_back(hit);
234 }
235
236 // Check for HS sectors connected to each other which could form a cluster
237 DepthFirstSearch(currentCell.first.first, currentCell.first.second);
238 // if cluster valid (i.e. not too small and not too big): finalize!
240
243 }
245 }
246}
247
248void SingleHoughSpaceFastInterceptFinder::DepthFirstSearch(uint lastIndexX, uint lastIndexY)
249{
251
252 for (uint currentIndexY = lastIndexY; currentIndexY >= lastIndexY - 1; currentIndexY--) {
253 if (abs((int)m_clusterInitialPosition.second - (int)currentIndexY) >= m_MaximumHSClusterSizeY or
254 m_clusterSize >= m_MaximumHSClusterSize or currentIndexY > m_nVerticalSectors) return;
255 for (uint currentIndexX = lastIndexX; currentIndexX <= lastIndexX + 1; currentIndexX++) {
256 if (abs((int)m_clusterInitialPosition.first - (int)currentIndexX) >= m_MaximumHSClusterSizeX or
257 m_clusterSize >= m_MaximumHSClusterSize or currentIndexX > m_nAngleSectors) return;
258
259 // The cell (currentIndexX, currentIndexY) is the current one has already been checked, so continue
260 if (lastIndexX == currentIndexX && lastIndexY == currentIndexY) continue;
261
262 // first check bounds to avoid out-of-bound array access
263 // as they are uints, they are always >= 0, and in case of an overflow they would be too large
264 if (currentIndexX < m_nAngleSectors and currentIndexY < m_nVerticalSectors) {
265
266 auto activeSector = m_activeSectors.find({currentIndexX, currentIndexY});
267 // Only continue searching if the current cluster is smaller than the maximum cluster size
268 if (activeSector != m_activeSectors.end() and activeSector->second.first < 0 /*and m_clusterSize < m_MaximumHSClusterSize*/) {
269 activeSector->second.first = m_clusterCount;
271
272 // No need to check whether currentIndex exists as a key in m_activeSectors as they were created at the same time
273 // so it's certain the key exists.
274 for (VXDHoughState* hit : activeSector->second.second) {
275 if (not TrackFindingCDC::is_in(hit, m_currentTrackCandidate)) {
276 m_currentTrackCandidate.emplace_back(hit);
277 }
278 }
279
280 // search in the next Hough Space cells...
281 DepthFirstSearch(currentIndexX, currentIndexY);
282 }
283 }
284 }
285 }
286}
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