Belle II Software  release-05-01-25
DATCONTrackingInterceptFinder.cc
1 /**************************************************************************
2  * BASF2 (Belle Analysis Framework 2) *
3  * Copyright(C) 2013 - Belle II Collaboration *
4  * *
5  * Author: The Belle II Collaboration *
6  * Contributors: Michael Schnell, Christian Wessel *
7  * *
8  * This software is provided "as is" without any warranty. *
9  **************************************************************************/
10 
11 #include <tracking/modules/DATCON/DATCONTrackingModule.h>
12 
13 using namespace std;
14 using namespace Belle2;
15 
16 /*
17 * Very simple sensor filter / layer filter / sensor layer filter.
18 */
19 bool
20 DATCONTrackingModule::layerFilter(bool* layer)
21 {
22  unsigned int layercount = 0;
23 
24  /* Count number of found layers */
25  for (int i = 0; i < 4; ++i) {
26  if (layer[i] == true) {
27  ++layercount;
28  }
29  }
30  // return true, if lines in sector are from at least minLines different layers
31  if (layercount >= m_minimumLines) {
32  return (true);
33  }
34 
35  return (false);
36 }
37 
38 /*
39 * Find the intercept in hough, or hess hough space. // does "hess hough space" here mean "conformal mapped hough space"?
40 * Return zero on success.
41 * Parameter for uSide:
42 * x = phi0
43 * y = r
44 * Parameter for v_side:
45 * x = m
46 * y = a
47 */
48 int
49 DATCONTrackingModule::fastInterceptFinder2d(houghMap& hits, bool uSide, TVector2 v1_s,
50  TVector2 v2_s,
51  TVector2 v4_s,
52  unsigned int iterations, unsigned int maxIterations)
53 {
54  int hitID;
55  unsigned int countLayer;
56  double unitX, unitY;
57  double y1, y2;
58  double m, a;
59  houghPair hp;
60  VxdID sensor;
61  vector<unsigned int> candidateIDList;
62  double xdiff, ydiff;
63  int sectorX, sectorY;
64  double baseX, baseY;
65  int maxSectorX, maxSectorY;
66  houghMap containedHits;
67 
68  TVector2 v1, v2, v3, v4;
69 
70  unitX = ((v2_s.X() - v1_s.X()) / 2.0);
71  unitY = ((v1_s.Y() - v4_s.Y()) / 2.0);
72 
73  if (uSide) {
74  baseX = -M_PI;
75  baseY = -m_rectSizeU;
76  maxSectorX = (int)pow(2, m_maxIterationsU + 1) - 1;
77  maxSectorY = (int)pow(2, m_maxIterationsU + 1) - 1;
78  } else {
79  baseX = -M_PI;
80  baseY = -m_rectSizeV;
81  maxSectorX = (int)pow(2, m_maxIterationsV + 1) - 1;
82  maxSectorY = (int)pow(2, m_maxIterationsV + 1) - 1;
83  }
84 
85  countLayer = 0;
86  for (int i = 0; i < 2 ; ++i) {
87  for (int j = 0; j < 2; ++j) {
88  v1.Set((v4_s.X() + (double) i * unitX), (v4_s.Y() + (double) j * unitY + unitY));
89  v2.Set((v4_s.X() + (double) i * unitX + unitX), (v4_s.Y() + (double) j * unitY + unitY));
90  v3.Set((v4_s.X() + (double) i * unitX + unitX), (v4_s.Y() + (double) j * unitY));
91  v4.Set((v4_s.X() + (double) i * unitX), (v4_s.Y() + (double) j * unitY));
92 
93  candidateIDList.clear();
94  bool layerHit[4] = {false}; /* For layer filter */
95  for (const auto& hit : hits) {
96  hitID = hit.first;
97  hp = hit.second;
98  sensor = hp.first;
99 
100  if (!uSide) {
101  m = hp.second.X();
102  a = hp.second.Y();
103  y1 = m * cos(v1.X()) + a * sin(v1.X());
104  y2 = m * cos(v2.X()) + a * sin(v2.X());
105  } else {
106  m = hp.second.X();
107  a = hp.second.Y();
108  y1 = 2. * (m * cos(v1.X()) + a * sin(v1.X()));
109  y2 = 2. * (m * cos(v2.X()) + a * sin(v2.X()));
110  if (m_usePhase2Simulation) {
111  m = hp.second.X();
112  a = hp.second.Y();
113  y1 = m * cos(v1.X()) + a * sin(v1.X());
114  y2 = m * cos(v2.X()) + a * sin(v2.X());
115  }
116  }
117 
118  /* Check if HS-parameter curve is inside (or outside) actual sub-HS */
119  if (y1 <= v1.Y() && y2 >= v3.Y() && y2 >= y1) {
120  if (iterations == maxIterations) {
121  candidateIDList.push_back(hitID);
122  }
123 
124  layerHit[sensor.getLayerNumber() - 3] = true; /* layer filter */
125  ++countLayer;
126  containedHits.insert(hit);
127  }
128  }
129 
130  if (countLayer >= m_minimumLines) {
131  if (layerFilter(layerHit)) {
132  // recursive / iterative call of fastInterceptFinder2d, until iterations = critIterations (critIterations-1),
133  // actual values for v1...v4 are new startingpoints
134  if (iterations != maxIterations) {
135  fastInterceptFinder2d(containedHits, uSide, v1, v2, v4,
136  iterations + 1, maxIterations);
137  } else {
138  if (!uSide) {
139  vHoughCand.push_back(DATCONHoughCand(candidateIDList, make_pair(v1, v3)));
140  if (m_useHoughSpaceClustering) {
141  xdiff = v4.X() - baseX;
142  ydiff = v4.Y() - baseY;
143  sectorX = round(xdiff / unitX);
144  sectorY = maxSectorY - round(ydiff / unitY);
145  if (sectorX > maxSectorX) {
146  sectorX = maxSectorX;
147  } else if (sectorX < 0) {
148  sectorX = 0;
149  }
150  if (sectorY > maxSectorY) {
151  sectorY = maxSectorY;
152  } else if (sectorY < 0) {
153  sectorY = 0;
154  }
155  ArrayOfActiveSectorsThetaHS[sectorY][sectorX] = -1;
156  vHoughSpaceClusterCand.push_back(DATCONHoughSpaceClusterCand(candidateIDList, TVector2(sectorX, sectorY)));
157 
158  activeSectorVectorTheta.push_back(true);
159  }
160 
161  } else {
162  uHoughCand.push_back(DATCONHoughCand(candidateIDList, make_pair(v1, v3)));
163  if (m_useHoughSpaceClustering) {
164  baseX = -M_PI;
165  baseY = -m_rectSizeU;
166  xdiff = v4.X() - baseX;
167  ydiff = v4.Y() - baseY;
168  maxSectorX = (int)pow(2, m_maxIterationsU + 1) - 1;
169  maxSectorY = (int)pow(2, m_maxIterationsU + 1) - 1;
170  sectorX = round(xdiff / unitX);
171  sectorY = maxSectorY - round(ydiff / unitY);
172  if (sectorX > maxSectorX) {
173  sectorX = maxSectorX;
174  } else if (sectorX < 0) {
175  sectorX = 0;
176  }
177  if (sectorY > maxSectorY) {
178  sectorY = maxSectorY;
179  } else if (sectorY < 0) {
180  sectorY = 0;
181  }
182  ArrayOfActiveSectorsPhiHS[sectorY][sectorX] = -1;
183  uHoughSpaceClusterCand.push_back(DATCONHoughSpaceClusterCand(candidateIDList, TVector2(sectorX, sectorY)));
184 
185  activeSectorVectorPhi.push_back(true);
186  }
187  }
188  }
189  }
190  } else {
191  if (m_useHoughSpaceClustering) {
192  if (uSide) {
193  activeSectorVectorPhi.push_back(false);
194  } else {
195  activeSectorVectorTheta.push_back(false);
196  }
197  }
198  }
199 
200  countLayer = 0;
201  }
202  }
203 
204  return 0;
205 }
206 
207 
208 /*
209 * Find the intercept in hough, or hess hough space. // does "hess hough space" here mean "conformal mapped hough space"?
210 * Return zero on success.
211 * Parameter for uSide:
212 * x = phi0
213 * y = r
214 * Parameter for v_side:
215 * x = m
216 * y = a
217 */
218 int
219 DATCONTrackingModule::slowInterceptFinder2d(houghMap& hits, bool uSide)
220 {
221  int hitID;
222  unsigned int countLayer;
223  double unitX, unitY;
224  double y1, y2;
225  double m, a;
226  houghPair hp;
227  VxdID sensor;
228  vector<unsigned int> candidateIDList;
229  unsigned short angleSectors, vertSectors;
230  double angleRange, vertRange;
231  double left, right, up, down;
232 
233  TVector2 v1, v2, v3, v4;
234  if (m_usePhase2Simulation) {
235  // ATTENTION TODO FIXME : This still has to be implemented!!!
236  // So far no phase 2 specific algorithms have been implemented and tested!
237  B2WARNING("This mode is not yet implemented, nothing will happen! Return...");
238  return 0;
239  } else {
240  if (uSide) {
241  angleSectors = m_nAngleSectorsU;
242  vertSectors = m_nVertSectorsU;
243  left = -M_PI;
244  right = M_PI;
245  up = m_rectSizeU;
246  down = -m_rectSizeU;
247  } else {
248  angleSectors = m_nAngleSectorsV;
249  vertSectors = m_nVertSectorsV;
250  left = -M_PI;
251  right = 0.;
252  up = m_rectSizeV;
253  down = -m_rectSizeV;
254  }
255  }
256 
257  angleRange = (right - left);
258  vertRange = (up - down);
259  unitX = angleRange / (double)(angleSectors);
260  unitY = vertRange / (double)(vertSectors);
261 
262  countLayer = 0;
263  for (int i = vertSectors - 1; i >= 0; i--) {
264  for (int j = 0; j < angleSectors; j++) {
265  v1.Set((left + (double) j * unitX), (up - (double) i * unitY));
266  v2.Set((left + (double) j * unitX + unitX), (up - (double) i * unitY));
267  v3.Set((left + (double) j * unitX + unitX), (up - (double) i * unitY - unitY));
268  v4.Set((left + (double) j * unitX), (up - (double) i * unitY - unitY));
269 
270 
271  candidateIDList.clear();
272  bool layerHit[4] = {false}; /* For layer filter */
273  for (const auto& hit : hits) {
274  hitID = hit.first;
275  hp = hit.second;
276  sensor = hp.first;
277 
278  if (!uSide) {
279  m = hp.second.X();
280  a = hp.second.Y();
281  y1 = m * cos(v1.X()) + a * sin(v1.X());
282  y2 = m * cos(v2.X()) + a * sin(v2.X());
283  } else {
284  m = hp.second.X();
285  a = hp.second.Y();
286  y1 = 2 * (m * cos(v1.X()) + a * sin(v1.X()));
287  y2 = 2 * (m * cos(v2.X()) + a * sin(v2.X()));
288  }
289 
290  if (y1 <= v1.Y() && y2 >= v3.Y() && y2 > y1) {
291  candidateIDList.push_back(hitID);
292  layerHit[sensor.getLayerNumber() - 3] = true; /* layer filter */
293  ++countLayer;
294  }
295  }
296  if (countLayer >= m_minimumLines) {
297  if (layerFilter(layerHit)) {
298  if (!uSide) {
299  vHoughCand.push_back(DATCONHoughCand(candidateIDList, make_pair(v1, v3)));
300  ArrayOfActiveSectorsThetaHS[i][j] = -1;
301  vHoughSpaceClusterCand.push_back(DATCONHoughSpaceClusterCand(candidateIDList, TVector2(j, i)));
302  activeSectorVectorTheta.push_back(true);
303  } else {
304  uHoughCand.push_back(DATCONHoughCand(candidateIDList, make_pair(v1, v3)));
305  ArrayOfActiveSectorsPhiHS[i][j] = -1;
306  uHoughSpaceClusterCand.push_back(DATCONHoughSpaceClusterCand(candidateIDList, TVector2(j, i)));
307  activeSectorVectorPhi.push_back(true);
308  }
309  }
310  } else {
311  if (uSide) {
312  activeSectorVectorPhi.push_back(false);
313  } else {
314  activeSectorVectorTheta.push_back(false);
315  }
316  }
317 
318  countLayer = 0;
319  } // for (int j ...
320  } // for (int i ...
321 
322  return 0;
323 }
Belle2::VxdID
Class to uniquely identify a any structure of the PXD and SVD.
Definition: VxdID.h:43
Belle2::houghMap
std::map< int, houghPair > houghMap
Map containing integer ID and corresponding houghPair for the HS TODO make this description better.
Definition: DATCONTrackingModule.h:73
Belle2::DATCONHoughCand
Hough Candidates class.
Definition: DATCONHoughCand.h:37
Belle2
Abstract base class for different kinds of events.
Definition: MillepedeAlgorithm.h:19
Belle2::DATCONHoughSpaceClusterCand
The Hough Space Cluster Candidate represents the candidates of clusters found in the Hough Space.
Definition: DATCONHoughSpaceClusterCand.h:35
Belle2::houghPair
std::pair< VxdID, TVector2 > houghPair
Hough Tuples.
Definition: DATCONTrackingModule.h:71