ROOT   Reference Guide
Searching...
No Matches
Functions.h
Go to the documentation of this file.
1// @(#)root/smatrix:$Id$
2// Authors: T. Glebe, L. Moneta 2005
3
4#ifndef ROOT_Math_Functions
5#define ROOT_Math_Functions
6// ********************************************************************
7//
8// source:
9//
10// type: source code
11//
12// created: 16. Mar 2001
13//
14// author: Thorsten Glebe
15// HERA-B Collaboration
16// Max-Planck-Institut fuer Kernphysik
17// Saupfercheckweg 1
18// 69117 Heidelberg
19// Germany
20// E-mail: T.Glebe@mpi-hd.mpg.de
21//
22// Description: Functions which are not applied like a unary operator
23//
24// changes:
25// 16 Mar 2001 (TG) creation
27// 07 Apr 2001 (TG) Lmag2, Lmag added
28// 19 Apr 2001 (TG) added #include <cmath>
29// 24 Apr 2001 (TG) added sign()
30// 26 Jul 2001 (TG) round() added
31// 27 Sep 2001 (TG) added Expr declaration
32//
33// ********************************************************************
34#include <cmath>
35
36#include "Math/Expression.h"
37
38/**
39\defgroup TempFunction Generic Template Functions
40\ingroup SMatrixGroup
41
42 These functions apply for any type T, such as a scalar, a vector or a matrix.
43*/
44
45/**
46\defgroup VectFunction Vector Template Functions
47\ingroup SMatrixGroup
48
49These functions apply to SVector types (and also to Vector expressions) and can
50return a vector expression or
51a scalar, like in the Dot product, or a matrix, like in the Tensor product
52*/
53
54
55namespace ROOT {
56
57 namespace Math {
58
59
60
61template <class T, unsigned int D> class SVector;
62
63
64/** square
65 Template function to compute \f$x\cdot x \f$, for any type T returning a type T
66
67 @ingroup TempFunction
68 @author T. Glebe
69*/
70//==============================================================================
71// square: x*x
72//==============================================================================
73template <class T>
74inline const T Square(const T& x) { return x*x; }
75
76/** maximum.
77 Template to find max(a,b) where a,b are of type T
78
79 @ingroup TempFunction
80 @author T. Glebe
81*/
82//==============================================================================
83// maximum
84//==============================================================================
85template <class T>
86inline const T Maximum(const T& lhs, const T& rhs) {
87 return (lhs > rhs) ? lhs : rhs;
88}
89
90/** minimum.
91 Template to find min(a,b) where a,b are of type T
92
93 @ingroup TempFunction
94 @author T. Glebe
95*/
96//==============================================================================
97// minimum
98//==============================================================================
99template <class T>
100inline const T Minimum(const T& lhs, const T& rhs) {
101 return (lhs < rhs) ? lhs : rhs;
102}
103
104/** round.
105 Template to compute nearest integer value for any type T
106 @ingroup TempFunction
107 @author T. Glebe
108*/
109//==============================================================================
110// round
111//==============================================================================
112template <class T>
113inline int Round(const T& x) {
114 return (x-static_cast<int>(x) < 0.5) ? static_cast<int>(x) : static_cast<int>(x+1);
115}
116
117
118/** sign.
119 Template to compute the sign of a number
120
121 @ingroup TempFunction
122 @author T. Glebe
123*/
124//==============================================================================
125// sign
126//==============================================================================
127template <class T>
128inline int Sign(const T& x) { return (x==0)? 0 : (x<0)? -1 : 1; }
129
130//==============================================================================
131// meta_dot
132//==============================================================================
133template <unsigned int I>
134struct meta_dot {
135 template <class A, class B, class T>
136 static inline T f(const A& lhs, const B& rhs, const T& x) {
137 return lhs.apply(I) * rhs.apply(I) + meta_dot<I-1>::f(lhs,rhs,x);
138 }
139};
140
141
142//==============================================================================
143// meta_dot<0>
144//==============================================================================
145template <>
146struct meta_dot<0> {
147 template <class A, class B, class T>
148 static inline T f(const A& lhs, const B& rhs, const T& /*x */) {
149 return lhs.apply(0) * rhs.apply(0);
150 }
151};
152
153
154/**
155 Vector dot product.
156 Template to compute \f$\vec{a}\cdot\vec{b} = \sum_i a_i\cdot b_i \f$.
157
158 @ingroup VectFunction
159 @author T. Glebe
160*/
161//==============================================================================
162// dot
163//==============================================================================
164template <class T, unsigned int D>
165inline T Dot(const SVector<T,D>& lhs, const SVector<T,D>& rhs) {
166 return meta_dot<D-1>::f(lhs,rhs, T());
167}
168
169//==============================================================================
170// dot
171//==============================================================================
172template <class A, class T, unsigned int D>
173inline T Dot(const SVector<T,D>& lhs, const VecExpr<A,T,D>& rhs) {
174 return meta_dot<D-1>::f(lhs,rhs, T());
175}
176
177//==============================================================================
178// dot
179//==============================================================================
180template <class A, class T, unsigned int D>
181inline T Dot(const VecExpr<A,T,D>& lhs, const SVector<T,D>& rhs) {
182 return meta_dot<D-1>::f(lhs,rhs, T());
183}
184
185
186//==============================================================================
187// dot
188//==============================================================================
189template <class A, class B, class T, unsigned int D>
190inline T Dot(const VecExpr<A,T,D>& lhs, const VecExpr<B,T,D>& rhs) {
191 return meta_dot<D-1>::f(rhs,lhs, T());
192}
193
194
195//==============================================================================
196// meta_mag
197//==============================================================================
198template <unsigned int I>
199struct meta_mag {
200 template <class A, class T>
201 static inline T f(const A& rhs, const T& x) {
202 return Square(rhs.apply(I)) + meta_mag<I-1>::f(rhs, x);
203 }
204};
205
206
207//==============================================================================
208// meta_mag<0>
209//==============================================================================
210template <>
211struct meta_mag<0> {
212 template <class A, class T>
213 static inline T f(const A& rhs, const T& ) {
214 return Square(rhs.apply(0));
215 }
216};
217
218
219/**
220 Vector magnitude square
221 Template to compute \f$|\vec{v}|^2 = \sum_iv_i^2 \f$.
222
223 @ingroup VectFunction
224 @author T. Glebe
225*/
226//==============================================================================
227// mag2
228//==============================================================================
229template <class T, unsigned int D>
230inline T Mag2(const SVector<T,D>& rhs) {
231 return meta_mag<D-1>::f(rhs, T());
232}
233
234//==============================================================================
235// mag2
236//==============================================================================
237template <class A, class T, unsigned int D>
238inline T Mag2(const VecExpr<A,T,D>& rhs) {
239 return meta_mag<D-1>::f(rhs, T());
240}
241
242/**
243 Vector magnitude (Euclidian norm)
244 Compute : \f$|\vec{v}| = \sqrt{\sum_iv_i^2} \f$.
245
246 @ingroup VectFunction
247 @author T. Glebe
248*/
249//==============================================================================
250// mag
251//==============================================================================
252template <class T, unsigned int D>
253inline T Mag(const SVector<T,D>& rhs) {
254 return std::sqrt(Mag2(rhs));
255}
256
257//==============================================================================
258// mag
259//==============================================================================
260template <class A, class T, unsigned int D>
261inline T Mag(const VecExpr<A,T,D>& rhs) {
262 return std::sqrt(Mag2(rhs));
263}
264
265
266/** Lmag2: Square of Minkowski Lorentz-Vector norm (only for 4D Vectors)
267 Template to compute \f$|\vec{v}|^2 = v_0^2 - v_1^2 - v_2^2 -v_3^2 \f$.
268
269 @ingroup VectFunction
270 @author T. Glebe
271*/
272//==============================================================================
273// Lmag2
274//==============================================================================
275template <class T>
276inline T Lmag2(const SVector<T,4>& rhs) {
277 return Square(rhs[0]) - Square(rhs[1]) - Square(rhs[2]) - Square(rhs[3]);
278}
279
280//==============================================================================
281// Lmag2
282//==============================================================================
283template <class A, class T>
284inline T Lmag2(const VecExpr<A,T,4>& rhs) {
285 return Square(rhs.apply(0))
286 - Square(rhs.apply(1)) - Square(rhs.apply(2)) - Square(rhs.apply(3));
287}
288
289/** Lmag: Minkowski Lorentz-Vector norm (only for 4-dim vectors)
290 Length of a vector Lorentz-Vector:
291 \f$|\vec{v}| = \sqrt{v_0^2 - v_1^2 - v_2^2 -v_3^2} \f$.
292
293 @ingroup VectFunction
294 @author T. Glebe
295*/
296//==============================================================================
297// Lmag
298//==============================================================================
299template <class T>
300inline T Lmag(const SVector<T,4>& rhs) {
301 return std::sqrt(Lmag2(rhs));
302}
303
304//==============================================================================
305// Lmag
306//==============================================================================
307template <class A, class T>
308inline T Lmag(const VecExpr<A,T,4>& rhs) {
309 return std::sqrt(Lmag2(rhs));
310}
311
312
313/** Vector Cross Product (only for 3-dim vectors)
314 \f$\vec{c} = \vec{a}\times\vec{b} \f$.
315
316 @ingroup VectFunction
317 @author T. Glebe
318*/
319//==============================================================================
320// cross product
321//==============================================================================
322template <class T>
323inline SVector<T,3> Cross(const SVector<T,3>& lhs, const SVector<T,3>& rhs) {
324 return SVector<T,3>(lhs.apply(1)*rhs.apply(2) -
325 lhs.apply(2)*rhs.apply(1),
326 lhs.apply(2)*rhs.apply(0) -
327 lhs.apply(0)*rhs.apply(2),
328 lhs.apply(0)*rhs.apply(1) -
329 lhs.apply(1)*rhs.apply(0));
330}
331
332//==============================================================================
333// cross product
334//==============================================================================
335template <class A, class T>
336inline SVector<T,3> Cross(const VecExpr<A,T,3>& lhs, const SVector<T,3>& rhs) {
337 return SVector<T,3>(lhs.apply(1)*rhs.apply(2) -
338 lhs.apply(2)*rhs.apply(1),
339 lhs.apply(2)*rhs.apply(0) -
340 lhs.apply(0)*rhs.apply(2),
341 lhs.apply(0)*rhs.apply(1) -
342 lhs.apply(1)*rhs.apply(0));
343}
344
345//==============================================================================
346// cross product
347//==============================================================================
348template <class T, class A>
349inline SVector<T,3> Cross(const SVector<T,3>& lhs, const VecExpr<A,T,3>& rhs) {
350 return SVector<T,3>(lhs.apply(1)*rhs.apply(2) -
351 lhs.apply(2)*rhs.apply(1),
352 lhs.apply(2)*rhs.apply(0) -
353 lhs.apply(0)*rhs.apply(2),
354 lhs.apply(0)*rhs.apply(1) -
355 lhs.apply(1)*rhs.apply(0));
356}
357
358//==============================================================================
359// cross product
360//==============================================================================
361template <class A, class B, class T>
362inline SVector<T,3> Cross(const VecExpr<A,T,3>& lhs, const VecExpr<B,T,3>& rhs) {
363 return SVector<T,3>(lhs.apply(1)*rhs.apply(2) -
364 lhs.apply(2)*rhs.apply(1),
365 lhs.apply(2)*rhs.apply(0) -
366 lhs.apply(0)*rhs.apply(2),
367 lhs.apply(0)*rhs.apply(1) -
368 lhs.apply(1)*rhs.apply(0));
369}
370
371
372/** Unit.
373 Return a vector of unit length: \f$\vec{e}_v = \vec{v}/|\vec{v}| \f$.
374
375 @ingroup VectFunction
376 @author T. Glebe
377*/
378//==============================================================================
379// unit: returns a unit vector
380//==============================================================================
381template <class T, unsigned int D>
382inline SVector<T,D> Unit(const SVector<T,D>& rhs) {
383 return SVector<T,D>(rhs).Unit();
384}
385
386//==============================================================================
387// unit: returns a unit vector
388//==============================================================================
389template <class A, class T, unsigned int D>
390inline SVector<T,D> Unit(const VecExpr<A,T,D>& rhs) {
391 return SVector<T,D>(rhs).Unit();
392}
393
394#ifdef XXX
395//==============================================================================
396// unit: returns a unit vector (worse performance)
397//==============================================================================
398template <class T, unsigned int D>
399inline VecExpr<BinaryOp<DivOp<T>, SVector<T,D>, Constant<T>, T>, T, D>
400 unit(const SVector<T,D>& lhs) {
401 typedef BinaryOp<DivOp<T>, SVector<T,D>, Constant<T>, T> DivOpBinOp;
402 return VecExpr<DivOpBinOp,T,D>(DivOpBinOp(DivOp<T>(),lhs,Constant<T>(mag(lhs))));
403}
404
405//==============================================================================
406// unit: returns a unit vector (worse performance)
407//==============================================================================
408template <class A, class T, unsigned int D>
409inline VecExpr<BinaryOp<DivOp<T>, VecExpr<A,T,D>, Constant<T>, T>, T, D>
410 unit(const VecExpr<A,T,D>& lhs) {
411 typedef BinaryOp<DivOp<T>, VecExpr<A,T,D>, Constant<T>, T> DivOpBinOp;
412 return VecExpr<DivOpBinOp,T,D>(DivOpBinOp(DivOp<T>(),lhs,Constant<T>(mag(lhs))));
413}
414#endif
415
416
417 } // namespace Math
418
419} // namespace ROOT
420
421
422
423#endif /* ROOT_Math_Functions */
SVector: a generic fixed size Vector class.
Definition SVector.h:75
SVector< T, D > & Unit()
transform vector into a vector of length 1
Definition SVector.icc:477
T apply(unsigned int i) const
access the parse tree. Index starts from zero
Definition SVector.icc:537
Expression wrapper class for Vector objects.
Definition Expression.h:64
T apply(unsigned int i) const
Definition Expression.h:77
const T Minimum(const T &lhs, const T &rhs)
minimum.
Definition Functions.h:100
int Sign(const T &x)
sign.
Definition Functions.h:128
int Round(const T &x)
round.
Definition Functions.h:113
const T Maximum(const T &lhs, const T &rhs)
maximum.
Definition Functions.h:86
SVector< T, D > Unit(const SVector< T, D > &rhs)
Unit.
Definition Functions.h:382
SVector< T, 3 > Cross(const SVector< T, 3 > &lhs, const SVector< T, 3 > &rhs)
Vector Cross Product (only for 3-dim vectors) .
Definition Functions.h:323
T Lmag2(const SVector< T, 4 > &rhs)
Lmag2: Square of Minkowski Lorentz-Vector norm (only for 4D Vectors) Template to compute .
Definition Functions.h:276
T Mag2(const SVector< T, D > &rhs)
Vector magnitude square Template to compute .
Definition Functions.h:230
T Mag(const SVector< T, D > &rhs)
Vector magnitude (Euclidian norm) Compute : .
Definition Functions.h:253
T Lmag(const SVector< T, 4 > &rhs)
Lmag: Minkowski Lorentz-Vector norm (only for 4-dim vectors) Length of a vector Lorentz-Vector: .
Definition Functions.h:300
Double_t x[n]
Definition legend1.C:17
#define I(x, y, z)
Namespace for new Math classes and functions.
tbb::task_arena is an alias of tbb::interface7::task_arena, which doesn't allow to forward declare tb...
#define Dot(u, v)
Definition normal.c:49
static T f(const A &lhs, const B &rhs, const T &)
Definition Functions.h:148
static T f(const A &lhs, const B &rhs, const T &x)
Definition Functions.h:136
static T f(const A &rhs, const T &)
Definition Functions.h:213
static T f(const A &rhs, const T &x)
Definition Functions.h:201
#define Square(a, x, y)
Definition triangle.c:4812