Logo ROOT  
Reference Guide
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
26 // 03 Apr 2001 (TG) minimum added, doc++ comments added
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 
49 These functions apply to SVector types (and also to Vector expressions) and can
50 return a vector expression or
51 a scalar, like in the Dot product, or a matrix, like in the Tensor product
52 */
53 
54 
55 namespace ROOT {
56 
57  namespace Math {
58 
59 
60 
61 template <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 //==============================================================================
73 template <class T>
74 inline 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 //==============================================================================
85 template <class T>
86 inline 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 //==============================================================================
99 template <class T>
100 inline 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 //==============================================================================
112 template <class T>
113 inline 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 //==============================================================================
127 template <class T>
128 inline int Sign(const T& x) { return (x==0)? 0 : (x<0)? -1 : 1; }
129 
130 //==============================================================================
131 // meta_dot
132 //==============================================================================
133 template <unsigned int I>
134 struct 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 //==============================================================================
145 template <>
146 struct 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 //==============================================================================
164 template <class T, unsigned int D>
165 inline 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 //==============================================================================
172 template <class A, class T, unsigned int D>
173 inline 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 //==============================================================================
180 template <class A, class T, unsigned int D>
181 inline 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 //==============================================================================
189 template <class A, class B, class T, unsigned int D>
190 inline 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 //==============================================================================
198 template <unsigned int I>
199 struct 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 //==============================================================================
210 template <>
211 struct 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 //==============================================================================
229 template <class T, unsigned int D>
230 inline T Mag2(const SVector<T,D>& rhs) {
231  return meta_mag<D-1>::f(rhs, T());
232 }
233 
234 //==============================================================================
235 // mag2
236 //==============================================================================
237 template <class A, class T, unsigned int D>
238 inline 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 //==============================================================================
252 template <class T, unsigned int D>
253 inline T Mag(const SVector<T,D>& rhs) {
254  return std::sqrt(Mag2(rhs));
255 }
256 
257 //==============================================================================
258 // mag
259 //==============================================================================
260 template <class A, class T, unsigned int D>
261 inline 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 //==============================================================================
275 template <class T>
276 inline 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 //==============================================================================
283 template <class A, class T>
284 inline 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 //==============================================================================
299 template <class T>
300 inline T Lmag(const SVector<T,4>& rhs) {
301  return std::sqrt(Lmag2(rhs));
302 }
303 
304 //==============================================================================
305 // Lmag
306 //==============================================================================
307 template <class A, class T>
308 inline 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 //==============================================================================
322 template <class T>
323 inline 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 //==============================================================================
335 template <class A, class T>
336 inline 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 //==============================================================================
348 template <class T, class A>
349 inline 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 //==============================================================================
361 template <class A, class B, class T>
362 inline 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 //==============================================================================
381 template <class T, unsigned int D>
382 inline 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 //==============================================================================
389 template <class A, class T, unsigned int D>
390 inline 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 //==============================================================================
398 template <class T, unsigned int D>
399 inline 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 //==============================================================================
408 template <class A, class T, unsigned int D>
409 inline 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 */
ROOT::Math::SVector
SVector: a generic fixed size Vector class.
Definition: SVector.h:75
ROOT::Math::Sign
int Sign(const T &x)
sign.
Definition: Functions.h:128
ROOT::Math::meta_mag::f
static T f(const A &rhs, const T &x)
Definition: Functions.h:201
ROOT::Math::VecExpr::apply
T apply(unsigned int i) const
Definition: Expression.h:77
x
Double_t x[n]
Definition: legend1.C:17
ROOT::Math::Cephes::A
static double A[]
Definition: SpecFuncCephes.cxx:170
ROOT::Math::Maximum
const T Maximum(const T &lhs, const T &rhs)
maximum.
Definition: Functions.h:86
ROOT::Math::Round
int Round(const T &x)
round.
Definition: Functions.h:113
ROOT::Math::Square
const T Square(const T &x)
square Template function to compute , for any type T returning a type T
Definition: Functions.h:74
ROOT::Math::Mag
T Mag(const SVector< T, D > &rhs)
Vector magnitude (Euclidian norm) Compute : .
Definition: Functions.h:253
ROOT::Math::VecExpr
Expression wrapper class for Vector objects.
Definition: Expression.h:64
ROOT::Math::meta_mag< 0 >::f
static T f(const A &rhs, const T &)
Definition: Functions.h:213
ROOT::Math::Unit
SVector< T, D > Unit(const SVector< T, D > &rhs)
Unit.
Definition: Functions.h:382
ROOT::Math::meta_dot
Definition: Functions.h:134
sqrt
double sqrt(double)
ROOT::Math::SVector::apply
T apply(unsigned int i) const
access the parse tree. Index starts from zero
Definition: SVector.icc:537
ROOT::Math::SVector::Unit
SVector< T, D > & Unit()
transform vector into a vector of length 1
Definition: SVector.icc:477
ROOT::Math::Mag2
T Mag2(const SVector< T, D > &rhs)
Vector magnitude square Template to compute .
Definition: Functions.h:230
ROOT::Math::Cross
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
ROOT::Math::meta_mag
Definition: Functions.h:199
ROOT::Math::Lmag2
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
ROOT::Math::Chebyshev::T
double T(double x)
Definition: ChebyshevPol.h:34
I
#define I(x, y, z)
ROOT::Math::meta_dot< 0 >::f
static T f(const A &lhs, const B &rhs, const T &)
Definition: Functions.h:148
ROOT::Math::Minimum
const T Minimum(const T &lhs, const T &rhs)
minimum.
Definition: Functions.h:100
ROOT::Math::meta_dot::f
static T f(const A &lhs, const B &rhs, const T &x)
Definition: Functions.h:136
ROOT::Math::Cephes::B
static double B[]
Definition: SpecFuncCephes.cxx:178
Expression.h
ROOT::Math::Dot
T Dot(const SVector< T, D > &lhs, const SVector< T, D > &rhs)
Vector dot product.
Definition: Functions.h:165
ROOT
tbb::task_arena is an alias of tbb::interface7::task_arena, which doesn't allow to forward declare tb...
Definition: EExecutionPolicy.hxx:4
Math
Namespace for new Math classes and functions.
ROOT::Math::Lmag
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