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