Project
Loading...
Searching...
No Matches
trigonometric.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 MATHUTILS_INCLUDE_MATHUTILS_DETAIL_TRIGONOMETRIC_H_
17#define MATHUTILS_INCLUDE_MATHUTILS_DETAIL_TRIGONOMETRIC_H_
18
19#include "GPUCommonDef.h"
20#include "GPUCommonMath.h"
23
24#ifndef GPUCA_GPUCODE_DEVICE
25#include <cmath>
26#include <tuple>
27#include <array>
28#endif
29
30namespace o2
31{
32namespace math_utils
33{
34namespace detail
35{
36
37template <typename T>
38GPUhdi() T to02Pi(T phi)
39{
40 T res = phi;
41 // ensure angle in [0:2pi] for the input in [-pi:pi] or [0:pi]
42 if (res < 0.f) {
44 }
45 return res;
46}
47
48template <typename T>
50{
51 phi = to02Pi<T>(phi);
52}
53
54template <typename T>
56{
57 T res = phi;
58 // ensure angle in [0:2pi] for the any input angle
59 while (res < 0.f) {
61 }
64 }
65 return res;
66}
67
68template <typename T>
69inline void bringTo02PiGen(T& phi)
70{
71 phi = to02PiGen<T>(phi);
72}
73
74template <typename T>
76{
77 T res = phi;
78 // ensure angle in [-pi:pi] for the input in [-pi:pi] or [0:pi]
81 } else if (res < -o2::constants::math::PI) {
83 }
84 return res;
85}
86
87template <typename T>
89{
90 phi = toPMPi<T>(phi);
91}
92
93template <typename T>
95{
96 // ensure angle in [-pi:pi] for any input angle
97 T res = phi;
98 while (res < -o2::constants::math::PI) {
100 }
101 while (res > o2::constants::math::PI) {
103 }
104 return res;
105}
106
107template <typename T>
108inline void bringToPMPiGen(T& phi)
109{
110 phi = toPMPiGen<T>(phi);
111}
112
113#ifdef __OPENCL__ // TODO: get rid of that stupid workaround for OpenCL template address spaces
114template <typename T, typename S, typename U>
115GPUhdi() void sincos(T ang, S& s, U& c)
116{
117 return o2::gpu::GPUCommonMath::SinCos(ang, s, c);
118}
119#else
120template <typename T>
121GPUhdi() void sincos(T ang, T& s, T& c)
122{
123 return o2::gpu::GPUCommonMath::SinCos(ang, s, c);
124}
125template <>
126GPUhdi() void sincos(double ang, double& s, double& c)
127{
128 return o2::gpu::GPUCommonMath::SinCosd(ang, s, c);
129}
130#endif
131
132#ifndef GPUCA_GPUCODE_DEVICE
133
134template <typename T>
135GPUhdi() std::tuple<T, T> sincos(T ang)
136{
137 T sin = 0;
138 T cos = 0;
139 sincos<T>(ang, sin, cos);
140 return std::make_tuple(sin, cos);
141}
142
143template <typename T>
144inline std::tuple<T, T> rotateZ(T xL, T yL, T snAlp, T csAlp)
145{
146 // 2D rotation of the point by angle alpha (local to global)
147 return std::make_tuple(xL * csAlp - yL * snAlp, xL * snAlp + yL * csAlp);
148}
149
150template <typename T>
151GPUhdi() std::tuple<T, T> rotateZInv(T xG, T yG, T snAlp, T csAlp)
152{
153 // inverse 2D rotation of the point by angle alpha (global to local)
154 return rotateZ<T>(xG, yG, -snAlp, csAlp);
155}
156
157#endif
158
159template <typename T>
160GPUhdi() void rotateZ(std::array<T, 3>& xy, T alpha)
161{
162 // transforms vector in tracking frame alpha to global frame
163 T sin, cos;
164 sincos<T>(alpha, sin, cos);
165 const T x = xy[0];
166 xy[0] = x * cos - xy[1] * sin;
167 xy[1] = x * sin + xy[1] * cos;
168}
169
170template <typename T>
172{
173 xG = xL * csAlp - yL * snAlp;
174 yG = xL * snAlp + yL * csAlp;
175 ;
176 // std::tie(xG, yG) = rotateZ<T>(xL, yL, snAlp, csAlp); // must not use std:: - versions on GPU
177}
178
179template <typename T>
181{
182 // inverse 2D rotation of the point by angle alpha (global to local)
183 rotateZ<T>(xG, yG, xL, yL, -snAlp, csAlp);
184}
185
186template <typename T = float>
187GPUhdi() constexpr T sectorDAlpha()
188{
190}
191
192template <typename T>
194{
195 // convert angle to sector ID, phi can be either in 0:2pi or -pi:pi convention
197 if (phi < 0.f) {
199 }
200 return sect % o2::constants::math::NSectors;
201}
202
203template <typename T>
204GPUhdi() T sector2Angle(int sect)
205{
206 // convert sector to its angle center, in -pi:pi convention
207 T ang = o2::constants::math::SectorSpanRad * (0.5f + sect);
208 bringToPMPi<T>(ang);
209 return ang;
210}
211
212template <typename T>
214{
215 // convert angle to its sector alpha
216 return sector2Angle<T>(angle2Sector<T>(phi));
217}
218
219template <typename T>
220GPUhdi() constexpr bool okForPhiMin(T phiMin, T phi)
221{
222 // check if phi is above the phiMin, phi's must be in 0-2pi range
223 const T dphi = phi - phiMin;
224 return ((dphi > 0 && dphi < constants::math::PI) || dphi < -constants::math::PI) ? true : false;
225}
226
227template <typename T>
228GPUhdi() constexpr bool okForPhiMax(T phiMax, T phi)
229{
230 // check if phi is below the phiMax, phi's must be in 0-2pi range
231 const T dphi = phi - phiMax;
232 return ((dphi < 0 && dphi > -constants::math::PI) || dphi > constants::math::PI) ? true : false;
233}
234
235template <typename T>
236GPUhdi() constexpr T meanPhiSmall(T phi0, T phi1)
237{
238 // return mean phi, assume phis in 0:2pi
239 T phi;
240 if (!okForPhiMin(phi0, phi1)) {
241 phi = phi0;
242 phi0 = phi1;
243 phi1 = phi;
244 }
245 if (phi0 > phi1) {
246 phi = 0.5 * (phi1 - (constants::math::TwoPI - phi0)); // wrap
247 } else {
248 phi = 0.5 * (phi0 + phi1);
249 }
251 return phi;
252}
253
254template <typename T>
255GPUhdi() constexpr T deltaPhiSmall(T phi0, T phi1)
256{
257 // return delta phi, assume phi is in 0:2pi
258 T del;
259 if (!okForPhiMin(phi0, phi1)) {
260 del = phi0;
261 phi0 = phi1;
262 phi1 = del;
263 }
264 del = phi1 - phi0;
265 if (del < 0) {
267 }
268 return del;
269}
270
271template <typename T>
272GPUhdi() T fastATan2(T y, T x)
273{
274 // Fast atan2(y,x) for any angle [-Pi,Pi]
275 // Average inaccuracy: 0.00048
276 // Max inaccuracy: 0.00084
277 // Speed: 6.2 times faster than atan2f()
278 constexpr T Pi = 3.1415926535897932384626433832795;
279
280 auto atan = [](T a) -> T {
281 // returns the arctan for the angular range [-Pi/4, Pi/4]
282 // the polynomial coefficients are taken from:
283 // https://stackoverflow.com/questions/42537957/fast-accurate-atan-arctan-approximation-algorithm
284 constexpr T A = 0.0776509570923569;
285 constexpr T B = -0.287434475393028;
286 constexpr T C = (Pi / 4 - A - B);
287 const T a2 = a * a;
288 return ((A * a2 + B) * a2 + C) * a;
289 };
290
291 auto atan2P = [atan](T yy, T xx) -> T {
292 // fast atan2(yy,xx) for the angular range [0,+Pi]
293 constexpr T Pi025 = 1 * Pi / 4;
294 constexpr T Pi075 = 3 * Pi / 4;
295 const T x1 = xx + yy; // point p1 (x1,y1) = (xx,yy) - Pi/4
296 const T y1 = yy - xx;
297 T phi0 = 0;
298 T tan = 0;
299 if (xx < 0) { // p1 is in the range [Pi/4, 3*Pi/4]
300 phi0 = Pi075;
301 tan = -x1 / y1;
302 } else { // p1 is in the range [-Pi/4, Pi/4]
303 phi0 = Pi025;
304 tan = y1 / x1;
305 }
306 return phi0 + atan(tan);
307 };
308
309 // fast atan2(y,x) for any angle [-Pi,Pi]
310 return copysign<T>(atan2P(o2::gpu::GPUCommonMath::Abs(y), x), y);
311}
312
313template <class T>
314GPUdi() T asin(T x)
315{
316 return o2::gpu::GPUCommonMath::ASin(x);
317};
318
319template <class T>
320GPUdi() T acos(T x)
321{
322 return o2::gpu::GPUCommonMath::ACos(x);
323};
324
325template <class T>
326GPUdi() T atan(T x)
327{
328 return o2::gpu::GPUCommonMath::ATan(x);
329};
330
331template <class T>
332GPUdi() T atan2(T y, T x)
333{
334 return o2::gpu::GPUCommonMath::ATan2(y, x);
335};
336
337template <class T>
338GPUdi() T sin(T x)
339{
340 return o2::gpu::GPUCommonMath::Sin(x);
341};
342
343template <class T>
344GPUdi() T cos(T x)
345{
346 return o2::gpu::GPUCommonMath::Cos(x);
347};
348
349template <class T>
350GPUdi() T tan(T x)
351{
352 return o2::gpu::GPUCommonMath::Tan(x);
353};
354
355template <class T>
356GPUdi() T twoPi()
357{
358 return o2::gpu::GPUCommonMath::TwoPi();
359};
360
361template <>
362GPUdi() double twoPi()
363{
365};
366
367template <class T>
368GPUdi() T pi()
369{
370 return o2::gpu::GPUCommonMath::Pi();
371}
372
373template <>
374GPUdi() double pi()
375{
377}
378
379#ifndef GPUCA_GPUCODE_DEVICE
380template <>
381GPUdi() double tan(double x)
382{
383 return std::tan(x);
384};
385template <>
386GPUdi() double sin(double x)
387{
388 return std::sin(x);
389};
390template <>
391GPUdi() double atan2(double y, double x)
392{
393 return std::atan2(y, x);
394};
395template <>
396GPUdi() double atan(double x)
397{
398 return std::atan(x);
399};
400template <>
401GPUdi() double asin(double x)
402{
403 return std::asin(x);
404};
405template <>
406GPUdi() double acos(double x)
407{
408 return std::acos(x);
409};
410template <>
411GPUdi() double cos(double x)
412{
413 return std::cos(x);
414};
415#endif
416
417} // namespace detail
418} // namespace math_utils
419} // namespace o2
420
421#endif /* MATHUTILS_INCLUDE_MATHUTILS_DETAIL_TRIGONOMETRIC_H_ */
useful math constants
uint32_t res
Definition RawData.h:0
Definition A.h:16
Definition B.h:16
GLfloat GLfloat GLfloat alpha
Definition glcorearb.h:279
GLint GLenum GLint x
Definition glcorearb.h:403
GLuint GLfloat GLfloat GLfloat GLfloat y1
Definition glcorearb.h:5034
GLuint GLfloat GLfloat GLfloat x1
Definition glcorearb.h:5034
GLenum array
Definition glcorearb.h:4274
typedef void(APIENTRYP PFNGLCULLFACEPROC)(GLenum mode)
GLboolean GLboolean GLboolean GLboolean a
Definition glcorearb.h:1233
constexpr int NSectors
constexpr float SectorSpanDeg
constexpr float TwoPI
constexpr float SectorSpanRad
constexpr float PI
constexpr float Rad2Deg
void bringTo02PiGen(T &phi)
return((dphi > 0 &&dphi< constants::math::PI)||dphi< -constants::math::PI) ? true GPUhdi() const expr bool okForPhiMax(T phiMax
return((dphi< 0 &&dphi > -constants::math::PI)||dphi > constants::math::PI) ? true T phi1
void bringToPMPiGen(T &phi)
std::tuple< T, T > rotateZ(T xL, T yL, T snAlp, T csAlp)
return copysign< T >(atan2P(o2::gpu::GPUCommonMath::Abs(y), x), y)
GPUdi() int nint(T x)
Definition basicMath.h:67
float angle2Alpha(float phi)
Definition Utils.h:203
int angle2Sector(float phi)
Definition Utils.h:183
float toPMPi(float phi)
Definition Utils.h:90
std::tuple< float, float > rotateZInv(float xG, float yG, float snAlp, float csAlp)
Definition Utils.h:142
float sector2Angle(int sect)
Definition Utils.h:193
void bringToPMPi(float &phi)
Definition Utils.h:100
a couple of static helper functions to create timestamp values for CCDB queries or override obsolete ...
Defining DataPointCompositeObject explicitly as copiable.