Project
Loading...
Searching...
No Matches
HelixHelper.h
Go to the documentation of this file.
1// Copyright 2019-2020 CERN and copyright holders of ALICE O2.
2// See https://alice-o2.web.cern.ch/copyright for details of the copyright holders.
3// All rights not expressly granted are reserved.
4//
5// This software is distributed under the terms of the GNU General Public
6// License v3 (GPL Version 3), copied verbatim in the file "COPYING".
7//
8// In applying this license CERN does not waive the privileges and immunities
9// granted to it by virtue of its status as an Intergovernmental Organization
10// or submit itself to any jurisdiction.
11
15
16#ifndef _ALICEO2_HELIX_HELPER_
17#define _ALICEO2_HELIX_HELPER_
18
20#include "MathUtils/Utils.h"
22
23namespace o2
24{
25namespace track
26{
27
29//< precalculated track radius, center, alpha sin,cos and their combinations
31 float c, s, cc, ss, cs; // cos ans sin of track alpha and their products
32
33 GPUdDefault() TrackAuxPar() = default;
34
35 template <typename T>
36 GPUd() TrackAuxPar(const T& trc, float bz)
37 {
38 set(trc, bz);
39 }
40 GPUd() float cosDif(const TrackAuxPar& t) const { return c * t.c + s * t.s; } // cos(alpha_this - alha_t)
41 GPUd() float sinDif(const TrackAuxPar& t) const { return s * t.c - c * t.s; } // sin(alpha_this - alha_t)
42
43 template <typename T>
44 GPUd() void set(const T& trc, float bz)
45 {
46 trc.getCircleParams(bz, *this, s, c);
47 cc = c * c;
48 ss = s * s;
49 cs = c * s;
50 }
52};
53
54//__________________________________________________________
55//< crossing coordinates of 2 circles
56struct CrossInfo {
57 static constexpr float MaxDistXYDef = 10.;
58 float xDCA[2] = {};
59 float yDCA[2] = {};
60 int nDCA = 0;
61
63 {
64 const auto& trcA = trax0.rC > trax1.rC ? trax0 : trax1; // designate the largest circle as A
65 const auto& trcB = trax0.rC > trax1.rC ? trax1 : trax0;
66 nDCA = 0;
67 float xDist = trcB.xC - trcA.xC, yDist = trcB.yC - trcA.yC;
68 float dist2 = xDist * xDist + yDist * yDist, dist = o2::gpu::GPUCommonMath::Sqrt(dist2), rsum = trcA.rC + trcB.rC;
69 if (dist < 1e-12) {
70 return nDCA; // circles are concentric?
71 }
72 if (dist > rsum) { // circles don't touch, chose a point in between
73 // the parametric equation of lines connecting the centers is
74 // x = x0 + t/dist * (x1-x0), y = y0 + t/dist * (y1-y0)
75 if (dist - rsum > maxDistXY) { // too large distance
76 return nDCA;
77 }
78 notTouchingXY(dist, xDist, yDist, trcA, trcB.rC, isCollinear);
79 } else if (auto dfr = dist + trcB.rC - trcA.rC; dfr < 0.) { // the small circle is nestled into large one w/o touching
80 if (dfr > -maxDistXY) {
81 // select the point of closest approach of 2 circles
82 notTouchingXY(dist, xDist, yDist, trcA, -trcB.rC, isCollinear);
83 } else {
84 return nDCA;
85 }
86 } else { // 2 intersection points
87 if (isCollinear) {
91 float r2r = trcA.rC + trcB.rC;
92 float r1_r = trcA.rC / r2r;
93 float r2_r = trcB.rC / r2r;
94 xDCA[0] = r2_r * trcA.xC + r1_r * trcB.xC;
95 yDCA[0] = r2_r * trcA.yC + r1_r * trcB.yC;
96 nDCA = 1;
97 } else if (o2::gpu::GPUCommonMath::Abs(xDist) < o2::gpu::GPUCommonMath::Abs(yDist)) {
98 // to simplify calculations, we move to new frame x->x+Xc0, y->y+Yc0, so that
99 // the 1st one is centered in origin
100 float a = (trcA.rC * trcA.rC - trcB.rC * trcB.rC + dist2) / (2. * yDist), b = -xDist / yDist, ab = a * b, bb = b * b;
101 float det = ab * ab - (1. + bb) * (a * a - trcA.rC * trcA.rC);
102 if (det > 0.) {
103 det = o2::gpu::GPUCommonMath::Sqrt(det);
104 xDCA[0] = (-ab + det) / (1. + b * b);
105 yDCA[0] = a + b * xDCA[0] + trcA.yC;
106 xDCA[0] += trcA.xC;
107 xDCA[1] = (-ab - det) / (1. + b * b);
108 yDCA[1] = a + b * xDCA[1] + trcA.yC;
109 xDCA[1] += trcA.xC;
110 nDCA = 2;
111 } else { // due to the finite precision the det<=0, i.e. the circles are barely touching, fall back to this special case
112 notTouchingXY(dist, xDist, yDist, trcA, trcB.rC);
113 }
114 } else {
115 float a = (trcA.rC * trcA.rC - trcB.rC * trcB.rC + dist2) / (2. * xDist), b = -yDist / xDist, ab = a * b, bb = b * b;
116 float det = ab * ab - (1. + bb) * (a * a - trcA.rC * trcA.rC);
117 if (det > 0.) {
118 det = o2::gpu::GPUCommonMath::Sqrt(det);
119 yDCA[0] = (-ab + det) / (1. + bb);
120 xDCA[0] = a + b * yDCA[0] + trcA.xC;
121 yDCA[0] += trcA.yC;
122 yDCA[1] = (-ab - det) / (1. + bb);
123 xDCA[1] = a + b * yDCA[1] + trcA.xC;
124 yDCA[1] += trcA.yC;
125 nDCA = 2;
126 } else { // due to the finite precision the det<=0, i.e. the circles are barely touching, fall back to this special case
127 notTouchingXY(dist, xDist, yDist, trcA, trcB.rC);
128 }
129 }
130 }
131 return nDCA;
132 }
133
134 GPUd() void notTouchingXY(float dist, float xDist, float yDist, const TrackAuxPar& trcA, float rBSign, bool isCollinear = false)
135 {
136 if (isCollinear) {
139 float r2r = trcA.rC + rBSign;
140 float r1_r = trcA.rC / r2r;
141 float r2_r = rBSign / r2r;
142 xDCA[0] = r2_r * trcA.xC + r1_r * (xDist + trcA.xC);
143 yDCA[0] = r2_r * trcA.yC + r1_r * (yDist + trcA.yC);
144 } else {
145 // fast method to calculate DCA between 2 circles, assuming that they don't touch each outer:
146 // the parametric equation of lines connecting the centers is x = xA + t/dist * xDist, y = yA + t/dist * yDist
147 // with xA,yY being the center of the circle A ( = trcA.xC, trcA.yC ), xDist = trcB.xC = trcA.xC ...
148 // There are 2 special cases:
149 // (a) small circle is inside the large one: provide rBSign as -trcB.rC
150 // (b) circle are side by side: provide rBSign as trcB.rC
151 auto t2d = (dist + trcA.rC - rBSign) / dist;
152 xDCA[0] = trcA.xC + 0.5 * (xDist * t2d);
153 yDCA[0] = trcA.yC + 0.5 * (yDist * t2d);
154 }
155 nDCA = 1;
156 }
157
158 template <typename T>
159 GPUd() int linesCrossInfo(const TrackAuxPar& trax0, const T& tr0,
160 const TrackAuxPar& trax1, const T& tr1, float maxDistXY = MaxDistXYDef)
161 {
175 nDCA = 0;
176 float dx = trax1.xC - trax0.xC; // for straight line TrackAuxPar stores lab coordinates at referene point!!!
177 float dy = trax1.yC - trax0.yC; //
178 float dz = tr1.getZ() - tr0.getZ();
179 auto csp0i2 = 1. / tr0.getCsp2(); // 1 / csp^2
180 auto csp0i = o2::gpu::GPUCommonMath::Sqrt(csp0i2);
181 auto tgp0 = tr0.getSnp() * csp0i;
182 float kx0 = trax0.c - trax0.s * tgp0;
183 float ky0 = trax0.s + trax0.c * tgp0;
184 float kz0 = tr0.getTgl() * csp0i;
185 auto csp1i2 = 1. / tr1.getCsp2(); // 1 / csp^2
186 auto csp1i = o2::gpu::GPUCommonMath::Sqrt(csp1i2);
187 auto tgp1 = tr1.getSnp() * o2::gpu::GPUCommonMath::Sqrt(csp1i2);
188 float kx1 = trax1.c - trax1.s * tgp1;
189 float ky1 = trax1.s + trax1.c * tgp1;
190 float kz1 = tr1.getTgl() * csp1i;
199 float a00 = (1.f + tr0.getTgl() * tr0.getTgl()) * csp0i2, a11 = (1.f + tr1.getTgl() * tr1.getTgl()) * csp1i2, a01 = -(kx0 * kx1 + ky0 * ky1 + kz0 * kz1);
200 float b0 = dx * kx0 + dy * ky0 + dz * kz0, b1 = -(dx * kx1 + dy * ky1 + dz * kz1);
201 float det = a00 * a11 - a01 * a01, det0 = b0 * a11 - b1 * a01, det1 = a00 * b1 - a01 * b0;
202 if (o2::gpu::GPUCommonMath::Sqrt(det) > o2::constants::math::Almost0) {
203 auto detI = 1. / det;
204 auto t0 = det0 * detI;
205 auto t1 = det1 * detI;
206 float addx0 = kx0 * t0, addy0 = ky0 * t0, addx1 = kx1 * t1, addy1 = ky1 * t1;
207 dx += addx1 - addx0; // recalculate XY distance at DCA
208 dy += addy1 - addy0;
209 if (dx * dx + dy * dy > maxDistXY * maxDistXY) {
210 return nDCA;
211 }
212 xDCA[0] = (trax0.xC + addx0 + trax1.xC + addx1) * 0.5;
213 yDCA[0] = (trax0.yC + addy0 + trax1.yC + addy1) * 0.5;
214 nDCA = 1;
215 }
216 return nDCA;
217 }
218
219 template <typename T>
220 GPUd() int circleLineCrossInfo(const TrackAuxPar& trax0, const T& tr0,
221 const TrackAuxPar& trax1, const T& tr1, float maxDistXY = MaxDistXYDef)
222 {
236
237 const auto& traxH = trax0.rC > trax1.rC ? trax0 : trax1; // circle (for the line rC is set to 0)
238 const auto& traxL = trax0.rC > trax1.rC ? trax1 : trax0; // line
239 const auto& trcL = trax0.rC > trax1.rC ? tr1 : tr0; // track of the line
240
241 // solve quadratic equation of line crossing the circle
242 float dx = traxL.xC - traxH.xC; // X distance between the line lab reference and circle center
243 float dy = traxL.yC - traxH.yC; // Y...
244 // t^2(kx^2+ky^2) + 2t(dx*kx+dy*ky) + dx^2 + dy^2 - r^2 = 0
245 auto cspi2 = 1. / trcL.getCsp2(); // 1 / csp^2 == kx^2 + ky^2
246 auto cspi = o2::gpu::GPUCommonMath::Sqrt(cspi2);
247 auto tgp = trcL.getSnp() * cspi;
248 float kx = traxL.c - traxL.s * tgp;
249 float ky = traxL.s + traxL.c * tgp;
250 double dk = dx * kx + dy * ky;
251 double det = dk * dk - cspi2 * (dx * dx + dy * dy - traxH.rC * traxH.rC);
252 if (det > 0) { // 2 crossings
253 det = o2::gpu::GPUCommonMath::Sqrt(det);
254 float t0 = (-dk + det) * cspi2;
255 float t1 = (-dk - det) * cspi2;
256 xDCA[0] = traxL.xC + kx * t0;
257 yDCA[0] = traxL.yC + ky * t0;
258 xDCA[1] = traxL.xC + kx * t1;
259 yDCA[1] = traxL.yC + ky * t1;
260 nDCA = 2;
261 } else {
262 // there is no crossing, find the point of the closest approach on the line which is closest to the circle center
263 float t = -dk * cspi2;
264 float xL = traxL.xC + kx * t, yL = traxL.yC + ky * t; // point on the line, need to average with point on the circle
265 float dxc = xL - traxH.xC, dyc = yL - traxH.yC, dist = o2::gpu::GPUCommonMath::Sqrt(dxc * dxc + dyc * dyc);
266 if (dist - traxH.rC > maxDistXY) { // too large distance
267 return nDCA;
268 }
269 float drcf = traxH.rC / dist; // radius / distance to circle center
270 float xH = traxH.xC + dxc * drcf, yH = traxH.yC + dyc * drcf;
271 xDCA[0] = (xL + xH) * 0.5;
272 yDCA[0] = (yL + yH) * 0.5;
273 nDCA = 1;
274 }
275 return nDCA;
276 }
277
278 template <typename T>
279 GPUd() int set(const TrackAuxPar& trax0, const T& tr0, const TrackAuxPar& trax1, const T& tr1, float maxDistXY = MaxDistXYDef, bool isCollinear = false)
280 {
281 // calculate up to 2 crossings between 2 circles
282 nDCA = 0;
283 if (trax0.rC > o2::constants::math::Almost0 && trax1.rC > o2::constants::math::Almost0) { // both are not straight lines
285 } else if (trax0.rC < o2::constants::math::Almost0 && trax1.rC < o2::constants::math::Almost0) { // both are straigt lines
287 } else {
289 }
290 //
291 return nDCA;
292 }
293
294 GPUdDefault() CrossInfo() = default;
295
296 template <typename T>
297 GPUd() CrossInfo(const TrackAuxPar& trax0, const T& tr0, const TrackAuxPar& trax1, const T& tr1, float maxDistXY = MaxDistXYDef, bool isCollinear = false)
298 {
299 set(trax0, tr0, trax1, tr1, maxDistXY, isCollinear);
300 }
302};
303
304} // namespace track
305} // namespace o2
306
307#endif
General auxilliary methods.
#define GPUdDefault()
const GPUTPCGMMerger::trackCluster & b1
const int16_t bb
useful math constants
Declarations of 2D primitives.
GLboolean GLboolean GLboolean b
Definition glcorearb.h:1233
typedef void(APIENTRYP PFNGLCULLFACEPROC)(GLenum mode)
GLuint GLfloat GLfloat GLfloat GLfloat GLfloat GLfloat GLfloat t0
Definition glcorearb.h:5034
GLboolean GLboolean GLboolean GLboolean a
Definition glcorearb.h:1233
GLuint GLfloat GLfloat GLfloat GLfloat GLfloat GLfloat GLfloat GLfloat GLfloat t1
Definition glcorearb.h:5034
constexpr float Almost0
GPUd() value_T BetheBlochSolid(value_T bg
Definition TrackUtils.h:150
a couple of static helper functions to create timestamp values for CCDB queries or override obsolete ...
int set(const TrackAuxPar &trax0, const T &tr0, const TrackAuxPar &trax1, const T &tr1, float maxDistXY=MaxDistXYDef)
int linesCrossInfo(const TrackAuxPar &trax0, const T &tr0, const TrackAuxPar &trax1, const T &tr1, float maxDistXY=MaxDistXYDef)
const TrackAuxPar & trax1
Definition HelixHelper.h:62
const TrackAuxPar float bool isCollinear
Definition HelixHelper.h:62
GPUd() int circlesCrossInfo(const TrackAuxPar &trax0
ClassDefNV(CrossInfo, 1)
const T const TrackAuxPar const T & tr1
int circlesCrossInfo(const TrackAuxPar &trax0, const TrackAuxPar &trax1, float maxDistXY=MaxDistXYDef)
Definition HelixHelper.h:62
const TrackAuxPar float maxDistXY
Definition HelixHelper.h:62
static constexpr float MaxDistXYDef
Definition HelixHelper.h:57
int circleLineCrossInfo(const TrackAuxPar &trax0, const T &tr0, const TrackAuxPar &trax1, const T &tr1, float maxDistXY=MaxDistXYDef)
void notTouchingXY(float dist, float xDist, float yDist, const TrackAuxPar &trcA, float rBSign)
float sinDif(const TrackAuxPar &t) const
Definition HelixHelper.h:41
GPUd() float cosDif(const TrackAuxPar &t) const
Definition HelixHelper.h:40
GPUd() void set(const T &trc
GPUdDefault() TrackAuxPar()=default
ClassDefNV(TrackAuxPar, 1)
void set(const T &trc, float bz)
Definition HelixHelper.h:44
GPUd() float sinDif(const TrackAuxPar &t) const
Definition HelixHelper.h:41
float cosDif(const TrackAuxPar &t) const
Definition HelixHelper.h:40