Belle II Software development
AxialHitQuadTreeProcessor.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/trackFindingCDC/legendre/quadtree/AxialHitQuadTreeProcessor.h>
9
10#include <tracking/trackFindingCDC/eventdata/hits/CDCWireHit.h>
11
12#include <array>
13#include <vector>
14
15#include <TF1.h>
16#include <TCanvas.h>
17#include <TGraph.h>
18#include <TAxis.h>
19
20using namespace Belle2;
21using namespace TrackFindingCDC;
22
23namespace {
24 bool sameSign(double n1, double n2, double n3, double n4)
25 {
26 return ((n1 > 0 && n2 > 0 && n3 > 0 && n4 > 0) || (n1 < 0 && n2 < 0 && n3 < 0 && n4 < 0));
27 }
28
29 using YSpan = AxialHitQuadTreeProcessor::YSpan;
30 YSpan splitCurvSpan(const YSpan& curvSpan, int nodeLevel, int lastLevel, int j)
31 {
32 const float meanCurv = curvSpan[0] + (curvSpan[1] - curvSpan[0]) / 2.0;
33 const std::array<float, 3> binBounds{curvSpan[0], meanCurv, curvSpan[1]};
34 const float binWidth = binBounds[j + 1] - binBounds[j];
35
36 const bool standardBinning = (nodeLevel <= lastLevel - 7) or (std::fabs(meanCurv) <= 0.005);
37 if (standardBinning) {
38 // B2INFO("Case 1* " << meanCurv << " " << (meanCurv <= 0.005));
39 float curv1 = binBounds[j];
40 float curv2 = binBounds[j + 1];
41
42 // Standard bin division
43 return {curv1, curv2};
44 }
45
46 // Non-standard binning
47 // For level 6 to 7 only expand 1 / 4, for higher levels expand 1 / 8.
48 // (assuming last level == 12)
49 if (nodeLevel < lastLevel - 5) {
50 // B2INFO("Case 2*");
51 float curv1 = binBounds[j] - binWidth / 4.;
52 float curv2 = binBounds[j + 1] + binWidth / 4.;
53 return {curv1, curv2};
54 } else {
55 // B2INFO("Case 3*");
56 float curv1 = binBounds[j] - binWidth / 8.;
57 float curv2 = binBounds[j + 1] + binWidth / 8.;
58 return {curv1, curv2};
59 }
60 }
61
62}
63
64std::vector<float> AxialHitQuadTreeProcessor::createCurvBound(YSpan curvSpan, int lastLevel)
65{
66 std::vector<YSpan> spans{{curvSpan}};
67
68 std::vector<YSpan> nextSpans;
69 for (int level = 1; level <= lastLevel; ++level) {
70 nextSpans.clear();
71 for (const YSpan& span : spans) {
72 nextSpans.push_back(splitCurvSpan(span, level, lastLevel, 0));
73 nextSpans.push_back(splitCurvSpan(span, level, lastLevel, 1));
74 }
75 spans.swap(nextSpans);
76 }
77
78 std::vector<float> bounds;
79 for (const YSpan& span : spans) {
80 bounds.push_back(span[0]);
81 bounds.push_back(span[1]);
82 }
83
84 assert(bounds.size() == std::pow(2, lastLevel));
85 return bounds;
86}
87
89{
90 static const int maxLevel = PrecisionUtil::getLookupGridLevel();
91 static const int nBins = std::pow(2, maxLevel);
92 static LookupTable<Vector2D> trigonometricLookUpTable(&Vector2D::Phi, nBins, -M_PI, M_PI);
93 return trigonometricLookUpTable;
94}
95
97 int seedLevel,
98 const XYSpans& ranges,
99 PrecisionUtil::PrecisionFunction precisionFunction)
100 : QuadTreeProcessor(lastLevel, seedLevel, ranges)
101 , m_precisionFunction(precisionFunction)
102 , m_localOrigin(0.0, 0.0)
104{
105 m_twoSidedPhaseSpace = m_quadTree->getYMin() * m_quadTree->getYMax() < 0;
106}
107
109 const YSpan& curvSpan,
110 const LookupTable<Vector2D>* cosSinLookupTable)
111 : QuadTreeProcessor(0, 0, { {0, cosSinLookupTable->getNPoints() - 1}, curvSpan})
112, m_localOrigin(localOrigin)
113, m_cosSinLookupTable(cosSinLookupTable)
114{
115 // Never use two sided mode in off origin extension
116 m_twoSidedPhaseSpace = false;
117}
118
119
121{
122 if (node->getLevel() <= 6) return false;
123 if (node->getLevel() >= getLastLevel()) return true;
124
125 const double nodeResolution = fabs(node->getYMin() - node->getYMax());
126 const double meanCurv = (node->getYMax() + node->getYMin()) / 2;
127
128 const double resolution = m_precisionFunction(meanCurv);
129 if (resolution >= nodeResolution) return true;
130
131 return false;
132}
133
136{
137 const int nodeLevel = node->getLevel();
138 const int lastLevel = getLastLevel();
139 const float meanCurv = std::fabs(node->getYMax() + node->getYMin()) / 2;
140
141 // Expand bins for all nodes 7 levels before the last level (for lastLevel = 12 starting at 6)
142 // but only in a curvature region higher than 0.005. Lower than that use always standard.
143 bool standardBinning = (nodeLevel <= lastLevel - 7) or (meanCurv <= 0.005);
144
145 if (standardBinning) {
146 float r1 = node->getYLowerBound(j);
147 float r2 = node->getYUpperBound(j);
148 long theta1 = node->getXLowerBound(i);
149 long theta2 = node->getXUpperBound(i);
150
151 // Standard bin division
152 return XYSpans({theta1, theta2}, {r1, r2});
153 }
154
155 // Non-standard binning
156 // For level 6 to 7 only expand 1 / 4, for higher levels expand 1 / 8.
157 // (assuming last level == 12)
158 if (nodeLevel < lastLevel - 5) {
159 float r1 = node->getYLowerBound(j) - node->getYBinWidth(j) / 4.;
160 float r2 = node->getYUpperBound(j) + node->getYBinWidth(j) / 4.;
161
162 // long extension = pow(2, lastLevel - nodeLevel) / 4; is same as:
163 long extension = pow(2, lastLevel - nodeLevel - 2);
164
165 long theta1 = node->getXLowerBound(i) - extension;
166 if (theta1 < 0) theta1 = 0;
167
168 long theta2 = node->getXUpperBound(i) + extension;
169 if (theta2 >= m_cosSinLookupTable->getNPoints()) {
170 theta2 = m_cosSinLookupTable->getNPoints() - 1;
171 }
172
173 return XYSpans({theta1, theta2}, {r1, r2});
174 } else {
175 float r1 = node->getYLowerBound(j) - node->getYBinWidth(j) / 8.;
176 float r2 = node->getYUpperBound(j) + node->getYBinWidth(j) / 8.;
177
178 // long extension = pow(2, lastLevel - nodeLevel) / 8; is same as
179 long extension = pow(2, lastLevel - nodeLevel - 3);
180
181 long theta1 = node->getXLowerBound(i) - extension;
182 if (theta1 < 0) theta1 = 0;
183
184 long theta2 = node->getXUpperBound(i) + extension;
185 if (theta2 >= m_cosSinLookupTable->getNPoints()) {
186 theta2 = m_cosSinLookupTable->getNPoints() - 1;
187 }
188
189 return XYSpans({theta1, theta2}, {r1, r2});
190 }
191}
192
194{
195 // Check whether the hit lies in the forward direction
196 if (node->getLevel() <= 4 and m_twoSidedPhaseSpace and node->getYMin() > -c_curlCurv and
197 node->getYMax() < c_curlCurv) {
198 if (not checkDerivative(node, wireHit)) return false;
199 }
200
201 const double& l = wireHit->getRefDriftLength();
202 const Vector2D& pos2D = wireHit->getRefPos2D() - m_localOrigin;
203 double r2 = pos2D.normSquared() - l * l;
204
205 using Quadlet = std::array<std::array<float, 2>, 2>;
206 Quadlet distRight{};
207 Quadlet distLeft{};
208
209 // get top and bottom borders of the node
210 float rMin = node->getYMin() * r2 / 2;
211 float rMax = node->getYMax() * r2 / 2;
212
213 // get left and right borders of the node
214 long thetaMin = node->getXMin();
215 long thetaMax = node->getXMax();
216
217 const Vector2D& thetaVecMin = m_cosSinLookupTable->at(thetaMin);
218 const Vector2D& thetaVecMax = m_cosSinLookupTable->at(thetaMax);
219
220 float rHitMin = thetaVecMin.dot(pos2D);
221 float rHitMax = thetaVecMax.dot(pos2D);
222
223 // compute sinograms at the left and right borders of the node
224 float rHitMinRight = rHitMin - l;
225 float rHitMaxRight = rHitMax - l;
226
227 float rHitMinLeft = rHitMin + l;
228 float rHitMaxLeft = rHitMax + l;
229
230 // Compute distance from the sinograms to bottom and top borders of the node
231 distRight[0][0] = rMin - rHitMinRight;
232 distRight[0][1] = rMin - rHitMaxRight;
233 distRight[1][0] = rMax - rHitMinRight;
234 distRight[1][1] = rMax - rHitMaxRight;
235
236 distLeft[0][0] = rMin - rHitMinLeft;
237 distLeft[0][1] = rMin - rHitMaxLeft;
238 distLeft[1][0] = rMax - rHitMinLeft;
239 distLeft[1][1] = rMax - rHitMaxLeft;
240
241 // Compare distance signs from sinograms to the node
242 // Check right
243 if (not sameSign(distRight[0][0], distRight[0][1], distRight[1][0], distRight[1][1])) {
244 return true;
245 }
246
247 // Check left
248 if (not sameSign(distLeft[0][0], distLeft[0][1], distLeft[1][0], distLeft[1][1])) {
249 return true;
250 }
251
252 // Check the extremum
253 float rHitMinExtr = thetaVecMin.cross(pos2D);
254 float rHitMaxExtr = thetaVecMax.cross(pos2D);
255 if (rHitMinExtr * rHitMaxExtr < 0.) return checkExtremum(node, wireHit);
256
257 // Not contained
258 return false;
259}
260
262{
263 const Vector2D& pos2D = wireHit->getRefPos2D() - m_localOrigin;
264
265 long thetaMin = node->getXMin();
266 long thetaMax = node->getXMax();
267
268 const Vector2D& thetaVecMin = m_cosSinLookupTable->at(thetaMin);
269 const Vector2D& thetaVecMax = m_cosSinLookupTable->at(thetaMax);
270
271 float rMinD = thetaVecMin.cross(pos2D);
272 float rMaxD = thetaVecMax.cross(pos2D);
273
274 // Does not really make sense...
275 if ((rMinD > 0) && (rMaxD * rMinD >= 0)) return true;
276 if ((rMaxD * rMinD < 0)) return true;
277 return false;
278}
279
281{
282 const double& l = wireHit->getRefDriftLength();
283 const Vector2D& pos2D = wireHit->getRefPos2D() - m_localOrigin;
284 double r2 = pos2D.normSquared() - l * l;
285
286 // get left and right borders of the node
287 long thetaMin = node->getXMin();
288 long thetaMax = node->getXMax();
289
290 const Vector2D& thetaVecMin = m_cosSinLookupTable->at(thetaMin);
291 const Vector2D& thetaVecMax = m_cosSinLookupTable->at(thetaMax);
292
293 if (not pos2D.isBetween(thetaVecMin, thetaVecMax)) return false;
294
295 // compute sinograms at the position
296 double r = pos2D.norm();
297 float rRight = r - l;
298 float rLeft = r + l;
299
300 // get top and bottom borders of the node
301 float rMin = node->getYMin() * r2 / 2;
302 float rMax = node->getYMax() * r2 / 2;
303
304 bool crossesRight = (rMin - rRight) * (rMax - rRight) < 0;
305 bool crossesLeft = (rMin - rLeft) * (rMax - rLeft) < 0;
306 return crossesRight or crossesLeft;
307}
308
309void AxialHitQuadTreeProcessor::drawHits(std::vector<const CDCWireHit*> hits, unsigned int color) const
310{
311 static int nevent(0);
312
313 TCanvas* canv = new TCanvas("canv", "legendre transform", 0, 0, 1200, 600);
314 canv->cd(1);
315 TGraph* dummyGraph = new TGraph();
316 dummyGraph->SetPoint(1, -M_PI, 0);
317 dummyGraph->SetPoint(2, M_PI, 0);
318 dummyGraph->Draw("AP");
319 dummyGraph->GetXaxis()->SetTitle("#theta");
320 dummyGraph->GetYaxis()->SetTitle("#rho");
321 dummyGraph->GetXaxis()->SetRangeUser(-M_PI, M_PI);
322 dummyGraph->GetYaxis()->SetRangeUser(-0.02, 0.15);
323
324 for (const CDCWireHit* wireHit : hits) {
325 const double& l = wireHit->getRefDriftLength();
326 const Vector2D& pos2D = wireHit->getRefPos2D() - m_localOrigin;
327 double x = pos2D.x();
328 double y = pos2D.y();
329 double r2 = pos2D.normSquared() - l * l;
330
331 TF1* concaveHitLegendre = new TF1("concaveHitLegendre", "2*([0]/[3])*cos(x) + 2*([1]/[3])*sin(x) + 2*([2]/[3])", -M_PI, M_PI);
332 TF1* convexHitLegendre = new TF1("convexHitLegendre", "2*([0]/[3])*cos(x) + 2*([1]/[3])*sin(x) - 2*([2]/[3])", -M_PI, M_PI);
333 concaveHitLegendre->SetLineWidth(1);
334 convexHitLegendre->SetLineWidth(1);
335 concaveHitLegendre->SetLineColor(color);
336 convexHitLegendre->SetLineColor(color);
337
338 concaveHitLegendre->SetParameters(x, y, l, r2);
339 convexHitLegendre->SetParameters(x, y, l, r2);
340 concaveHitLegendre->Draw("CSAME");
341 convexHitLegendre->Draw("CSAME");
342 }
343// canv->Print(Form("legendreHits_%i.root", nevent));
344// canv->Print(Form("legendreHits_%i.eps", nevent));
345 canv->Print(Form("legendreHits_%i.png", nevent));
346 delete canv;
347
348 nevent++;
349}
350
352{
353 std::vector<const CDCWireHit*> hits;
354 for (Item* item : node->getItems()) {
355 const CDCWireHit* wireHit = item->getPointer();
356 hits.push_back(wireHit);
357 }
358 drawHits(hits);
359}
XYSpans createChild(QuadTree *node, int i, int j) const final
Return the new ranges.
const double c_curlCurv
The curvature above which the trajectory is considered a curler.
bool checkDerivative(QuadTree *node, const CDCWireHit *wireHit) const
Check derivative of the sinogram.
void drawHits(std::vector< const CDCWireHit * > hits, unsigned int color=46) const
Draw QuadTree node.
bool m_twoSidedPhaseSpace
Indicator whether the two sided phases space insertion check should be used This option should automa...
Vector2D m_localOrigin
Local origin on which the phase space coordinates are centered.
const LookupTable< Vector2D > * m_cosSinLookupTable
Pinned lookup table for precomputed cosine and sine values.
bool isLeaf(QuadTree *node) const final
lastLevel depends on curvature of the track candidate
bool checkExtremum(QuadTree *node, const CDCWireHit *wireHit) const
Checks whether extreme point is located within QuadTree node's ranges.
AxialHitQuadTreeProcessor(int lastLevel, int seedLevel, const XYSpans &ranges, PrecisionUtil::PrecisionFunction precisionFunction)
Constructor.
PrecisionUtil::PrecisionFunction m_precisionFunction
Lambda which holds resolution function for the quadtree.
bool isInNode(QuadTree *node, const CDCWireHit *wireHit) const final
Check whether hit belongs to the quadtree node:
static const LookupTable< Vector2D > & getCosSinLookupTable()
Get the standard lookup table containing equally spaces unit vectors (cos, sin)
void drawNode(QuadTree *node) const
Draw QuadTree node.
static std::vector< float > createCurvBound(YSpan curvSpan, int lastLevel)
Constructs an array with the curvature bounds as generated by the default bin divisions.
Class representing a hit wire in the central drift chamber.
Definition CDCWireHit.h:55
Class which holds precomputed values of a function.
Definition LookupTable.h:50
int getNPoints() const
Return the number of finite sampling points in this lookup table.
static constexpr int getLookupGridLevel()
Returns desired deepness of the trigonometrical lookup table. Used as template parameter for the Trig...
std::function< double(double)> PrecisionFunction
Function type which is used for resolution calculations (resolution=f(curvature)) Takes a curvature v...
std::vector< AItem * > & getItems()
Get items from node.
AY getYBinWidth(int iBin)
Getter for the width of the iBin bin in "r" direction.
AX getXMax() const
Get maximal "Theta" value of the node.
AY getYMax() const
Get maximal "r" value of the node.
AY getYMin() const
Get minimal "r" value of the node.
int getLevel() const
Returns level of the node in tree (i.e., how much ancestors the node has)
AX getXLowerBound(int iBin) const
Get lower "Theta" value of given bin.
AY getYUpperBound(int iBin) const
Get upper "r" value of given bin.
AY getYLowerBound(int iBin) const
Get lower "r" value of given bin.
AX getXMin() const
Get minimal "Theta" value of the node.
AX getXUpperBound(int iBin) const
Get upper "Theta" value of given bin.
QuadTreeProcessor(int lastLevel, int seedLevel, const XYSpans &xySpans, bool debugOutput=false)
A two dimensional vector which is equipped with functions for correct handling of orientation relate...
Definition Vector2D.h:32
double dot(const Vector2D &rhs) const
Calculates the two dimensional dot product.
Definition Vector2D.h:158
double cross(const Vector2D &rhs) const
Calculated the two dimensional cross product.
Definition Vector2D.h:163
double x() const
Getter for the x coordinate.
Definition Vector2D.h:595
bool isBetween(const Vector2D &lower, const Vector2D &upper) const
Checks if this vector is between two other vectors Between means here that when rotating the lower ve...
Definition Vector2D.h:525
double normSquared() const
Calculates .
Definition Vector2D.h:169
double y() const
Getter for the y coordinate.
Definition Vector2D.h:605
double norm() const
Calculates the length of the vector.
Definition Vector2D.h:175
static Vector2D Phi(const double phi)
Constructs a unit vector with azimuth angle equal to phi.
Definition Vector2D.h:62
bool sameSign(double expected, double actual)
Predicate checking that two values have the same sign.
Abstract base class for different kinds of events.