Back to home page

EIC code displayed by LXR

 
 

    


File indexing completed on 2025-01-18 10:10:16

0001 // @(#)root/smatrix:$Id$
0002 // Authors: T. Glebe, L. Moneta    2005
0003 
0004 #ifndef ROOT_Math_Functions
0005 #define ROOT_Math_Functions
0006 // ********************************************************************
0007 //
0008 // source:
0009 //
0010 // type:      source code
0011 //
0012 // created:   16. Mar 2001
0013 //
0014 // author:    Thorsten Glebe
0015 //            HERA-B Collaboration
0016 //            Max-Planck-Institut fuer Kernphysik
0017 //            Saupfercheckweg 1
0018 //            69117 Heidelberg
0019 //            Germany
0020 //            E-mail: T.Glebe@mpi-hd.mpg.de
0021 //
0022 // Description: Functions which are not applied like a unary operator
0023 //
0024 // changes:
0025 // 16 Mar 2001 (TG) creation
0026 // 03 Apr 2001 (TG) minimum added, doc++ comments added
0027 // 07 Apr 2001 (TG) Lmag2, Lmag added
0028 // 19 Apr 2001 (TG) added #include <cmath>
0029 // 24 Apr 2001 (TG) added sign()
0030 // 26 Jul 2001 (TG) round() added
0031 // 27 Sep 2001 (TG) added Expr declaration
0032 //
0033 // ********************************************************************
0034 #include <cmath>
0035 
0036 #include "Math/Expression.h"
0037 
0038 /**
0039 \defgroup TempFunction Generic Template Functions
0040 \ingroup SMatrixGroup
0041 
0042    These functions apply for any type T, such as a scalar, a vector or a matrix.
0043 */
0044 
0045 /**
0046 \defgroup VectFunction Vector Template Functions
0047 \ingroup SMatrixGroup
0048 
0049 These functions apply to SVector types (and also to Vector expressions) and can
0050 return a vector expression or
0051 a scalar, like in the Dot product, or a matrix, like in the Tensor product
0052 */
0053 
0054 
0055 namespace ROOT {
0056 
0057   namespace Math {
0058 
0059 
0060 
0061 template <class T, unsigned int D> class SVector;
0062 
0063 
0064 /** square
0065     Template function to compute \f$x\cdot x \f$, for any type T returning a type T
0066 
0067     @ingroup TempFunction
0068     @author T. Glebe
0069 */
0070 //==============================================================================
0071 // square: x*x
0072 //==============================================================================
0073 template <class T>
0074 inline const T Square(const T& x) { return x*x; }
0075 
0076 /** maximum.
0077     Template to find max(a,b) where a,b are of type T
0078 
0079     @ingroup TempFunction
0080     @author T. Glebe
0081 */
0082 //==============================================================================
0083 // maximum
0084 //==============================================================================
0085 template <class T>
0086 inline const T Maximum(const T& lhs, const T& rhs) {
0087   return (lhs > rhs) ? lhs : rhs;
0088 }
0089 
0090 /** minimum.
0091     Template to find min(a,b) where a,b are of type T
0092 
0093     @ingroup TempFunction
0094     @author T. Glebe
0095 */
0096 //==============================================================================
0097 // minimum
0098 //==============================================================================
0099 template <class T>
0100 inline const T Minimum(const T& lhs, const T& rhs) {
0101   return (lhs < rhs) ? lhs : rhs;
0102 }
0103 
0104 /** round.
0105     Template to compute nearest integer value for any type T
0106     @ingroup TempFunction
0107     @author T. Glebe
0108 */
0109 //==============================================================================
0110 // round
0111 //==============================================================================
0112 template <class T>
0113 inline int Round(const T& x) {
0114   return (x-static_cast<int>(x) < 0.5) ? static_cast<int>(x) : static_cast<int>(x+1);
0115 }
0116 
0117 
0118 /** sign.
0119     Template to compute the sign of a number
0120 
0121     @ingroup TempFunction
0122     @author T. Glebe
0123 */
0124 //==============================================================================
0125 // sign
0126 //==============================================================================
0127 template <class T>
0128 inline int Sign(const T& x) { return (x==0)? 0 : (x<0)? -1 : 1; }
0129 
0130 //==============================================================================
0131 // meta_dot
0132 //==============================================================================
0133 template <unsigned int I>
0134 struct meta_dot {
0135   template <class A, class B, class T>
0136   static inline T f(const A& lhs, const B& rhs, const T& x) {
0137     return lhs.apply(I) * rhs.apply(I) + meta_dot<I-1>::f(lhs,rhs,x);
0138   }
0139 };
0140 
0141 
0142 //==============================================================================
0143 // meta_dot<0>
0144 //==============================================================================
0145 template <>
0146 struct meta_dot<0> {
0147   template <class A, class B, class T>
0148   static inline T f(const A& lhs, const B& rhs, const T& /*x */) {
0149     return lhs.apply(0) * rhs.apply(0);
0150   }
0151 };
0152 
0153 
0154 /**
0155     Vector dot product.
0156     Template to compute \f$\vec{a}\cdot\vec{b} = \sum_i a_i\cdot b_i \f$.
0157 
0158     @ingroup VectFunction
0159     @author T. Glebe
0160 */
0161 //==============================================================================
0162 // dot
0163 //==============================================================================
0164 template <class T, unsigned int D>
0165 inline T Dot(const SVector<T,D>& lhs, const SVector<T,D>& rhs) {
0166   return meta_dot<D-1>::f(lhs,rhs, T());
0167 }
0168 
0169 //==============================================================================
0170 // dot
0171 //==============================================================================
0172 template <class A, class T, unsigned int D>
0173 inline T Dot(const SVector<T,D>& lhs, const VecExpr<A,T,D>& rhs) {
0174   return meta_dot<D-1>::f(lhs,rhs, T());
0175 }
0176 
0177 //==============================================================================
0178 // dot
0179 //==============================================================================
0180 template <class A, class T, unsigned int D>
0181 inline T Dot(const VecExpr<A,T,D>& lhs, const SVector<T,D>& rhs) {
0182   return meta_dot<D-1>::f(lhs,rhs, T());
0183 }
0184 
0185 
0186 //==============================================================================
0187 // dot
0188 //==============================================================================
0189 template <class A, class B, class T, unsigned int D>
0190 inline T Dot(const VecExpr<A,T,D>& lhs, const VecExpr<B,T,D>& rhs) {
0191   return meta_dot<D-1>::f(rhs,lhs, T());
0192 }
0193 
0194 
0195 //==============================================================================
0196 // meta_mag
0197 //==============================================================================
0198 template <unsigned int I>
0199 struct meta_mag {
0200   template <class A, class T>
0201   static inline T f(const A& rhs, const T& x) {
0202     return Square(rhs.apply(I)) + meta_mag<I-1>::f(rhs, x);
0203   }
0204 };
0205 
0206 
0207 //==============================================================================
0208 // meta_mag<0>
0209 //==============================================================================
0210 template <>
0211 struct meta_mag<0> {
0212   template <class A, class T>
0213   static inline T f(const A& rhs, const T& ) {
0214     return Square(rhs.apply(0));
0215   }
0216 };
0217 
0218 
0219 /**
0220     Vector magnitude square
0221     Template to compute \f$|\vec{v}|^2 = \sum_iv_i^2 \f$.
0222 
0223     @ingroup VectFunction
0224     @author T. Glebe
0225 */
0226 //==============================================================================
0227 // mag2
0228 //==============================================================================
0229 template <class T, unsigned int D>
0230 inline T Mag2(const SVector<T,D>& rhs) {
0231   return meta_mag<D-1>::f(rhs, T());
0232 }
0233 
0234 //==============================================================================
0235 // mag2
0236 //==============================================================================
0237 template <class A, class T, unsigned int D>
0238 inline T Mag2(const VecExpr<A,T,D>& rhs) {
0239   return meta_mag<D-1>::f(rhs, T());
0240 }
0241 
0242 /**
0243     Vector magnitude (Euclidean norm)
0244     Compute : \f$ |\vec{v}| = \sqrt{\sum_iv_i^2} \f$.
0245 
0246     @ingroup VectFunction
0247     @author T. Glebe
0248 */
0249 //==============================================================================
0250 // mag
0251 //==============================================================================
0252 template <class T, unsigned int D>
0253 inline T Mag(const SVector<T,D>& rhs) {
0254   return std::sqrt(Mag2(rhs));
0255 }
0256 
0257 //==============================================================================
0258 // mag
0259 //==============================================================================
0260 template <class A, class T, unsigned int D>
0261 inline T Mag(const VecExpr<A,T,D>& rhs) {
0262   return std::sqrt(Mag2(rhs));
0263 }
0264 
0265 
0266 /** Lmag2: Square of Minkowski Lorentz-Vector norm (only for 4D Vectors)
0267     Template to compute \f$ |\vec{v}|^2 = v_0^2 - v_1^2 - v_2^2 -v_3^2 \f$.
0268 
0269     @ingroup VectFunction
0270     @author T. Glebe
0271 */
0272 //==============================================================================
0273 // Lmag2
0274 //==============================================================================
0275 template <class T>
0276 inline T Lmag2(const SVector<T,4>& rhs) {
0277   return Square(rhs[0]) - Square(rhs[1]) - Square(rhs[2]) - Square(rhs[3]);
0278 }
0279 
0280 //==============================================================================
0281 // Lmag2
0282 //==============================================================================
0283 template <class A, class T>
0284 inline T Lmag2(const VecExpr<A,T,4>& rhs) {
0285   return Square(rhs.apply(0))
0286     - Square(rhs.apply(1)) - Square(rhs.apply(2)) - Square(rhs.apply(3));
0287 }
0288 
0289 /** Lmag: Minkowski Lorentz-Vector norm (only for 4-dim vectors)
0290     Length of a vector Lorentz-Vector:
0291     \f$ |\vec{v}| = \sqrt{v_0^2 - v_1^2 - v_2^2 -v_3^2} \f$.
0292 
0293     @ingroup VectFunction
0294     @author T. Glebe
0295 */
0296 //==============================================================================
0297 // Lmag
0298 //==============================================================================
0299 template <class T>
0300 inline T Lmag(const SVector<T,4>& rhs) {
0301   return std::sqrt(Lmag2(rhs));
0302 }
0303 
0304 //==============================================================================
0305 // Lmag
0306 //==============================================================================
0307 template <class A, class T>
0308 inline T Lmag(const VecExpr<A,T,4>& rhs) {
0309   return std::sqrt(Lmag2(rhs));
0310 }
0311 
0312 
0313 /** Vector Cross Product (only for 3-dim vectors)
0314     \f$ \vec{c} = \vec{a}\times\vec{b} \f$.
0315 
0316     @ingroup VectFunction
0317     @author T. Glebe
0318 */
0319 //==============================================================================
0320 // cross product
0321 //==============================================================================
0322 template <class T>
0323 inline SVector<T,3> Cross(const SVector<T,3>& lhs, const SVector<T,3>& rhs) {
0324   return SVector<T,3>(lhs.apply(1)*rhs.apply(2) -
0325                       lhs.apply(2)*rhs.apply(1),
0326                       lhs.apply(2)*rhs.apply(0) -
0327                       lhs.apply(0)*rhs.apply(2),
0328                       lhs.apply(0)*rhs.apply(1) -
0329                       lhs.apply(1)*rhs.apply(0));
0330 }
0331 
0332 //==============================================================================
0333 // cross product
0334 //==============================================================================
0335 template <class A, class T>
0336 inline SVector<T,3> Cross(const VecExpr<A,T,3>& lhs, const SVector<T,3>& rhs) {
0337   return SVector<T,3>(lhs.apply(1)*rhs.apply(2) -
0338                       lhs.apply(2)*rhs.apply(1),
0339                       lhs.apply(2)*rhs.apply(0) -
0340                       lhs.apply(0)*rhs.apply(2),
0341                       lhs.apply(0)*rhs.apply(1) -
0342                       lhs.apply(1)*rhs.apply(0));
0343 }
0344 
0345 //==============================================================================
0346 // cross product
0347 //==============================================================================
0348 template <class T, class A>
0349 inline SVector<T,3> Cross(const SVector<T,3>& lhs, const VecExpr<A,T,3>& rhs) {
0350   return SVector<T,3>(lhs.apply(1)*rhs.apply(2) -
0351                       lhs.apply(2)*rhs.apply(1),
0352                       lhs.apply(2)*rhs.apply(0) -
0353                       lhs.apply(0)*rhs.apply(2),
0354                       lhs.apply(0)*rhs.apply(1) -
0355                       lhs.apply(1)*rhs.apply(0));
0356 }
0357 
0358 //==============================================================================
0359 // cross product
0360 //==============================================================================
0361 template <class A, class B, class T>
0362 inline SVector<T,3> Cross(const VecExpr<A,T,3>& lhs, const VecExpr<B,T,3>& rhs) {
0363   return SVector<T,3>(lhs.apply(1)*rhs.apply(2) -
0364                       lhs.apply(2)*rhs.apply(1),
0365                       lhs.apply(2)*rhs.apply(0) -
0366                       lhs.apply(0)*rhs.apply(2),
0367                       lhs.apply(0)*rhs.apply(1) -
0368                       lhs.apply(1)*rhs.apply(0));
0369 }
0370 
0371 
0372 /** Unit.
0373     Return a vector of unit length: \f$ \vec{e}_v = \vec{v}/|\vec{v}| \f$.
0374 
0375     @ingroup VectFunction
0376     @author T. Glebe
0377 */
0378 //==============================================================================
0379 // unit: returns a unit vector
0380 //==============================================================================
0381 template <class T, unsigned int D>
0382 inline SVector<T,D> Unit(const SVector<T,D>& rhs) {
0383   return SVector<T,D>(rhs).Unit();
0384 }
0385 
0386 //==============================================================================
0387 // unit: returns a unit vector
0388 //==============================================================================
0389 template <class A, class T, unsigned int D>
0390 inline SVector<T,D> Unit(const VecExpr<A,T,D>& rhs) {
0391   return SVector<T,D>(rhs).Unit();
0392 }
0393 
0394 #ifdef XXX
0395 //==============================================================================
0396 // unit: returns a unit vector (worse performance)
0397 //==============================================================================
0398 template <class T, unsigned int D>
0399 inline VecExpr<BinaryOp<DivOp<T>, SVector<T,D>, Constant<T>, T>, T, D>
0400  unit(const SVector<T,D>& lhs) {
0401   typedef BinaryOp<DivOp<T>, SVector<T,D>, Constant<T>, T> DivOpBinOp;
0402   return VecExpr<DivOpBinOp,T,D>(DivOpBinOp(DivOp<T>(),lhs,Constant<T>(mag(lhs))));
0403 }
0404 
0405 //==============================================================================
0406 // unit: returns a unit vector (worse performance)
0407 //==============================================================================
0408 template <class A, class T, unsigned int D>
0409 inline VecExpr<BinaryOp<DivOp<T>, VecExpr<A,T,D>, Constant<T>, T>, T, D>
0410  unit(const VecExpr<A,T,D>& lhs) {
0411   typedef BinaryOp<DivOp<T>, VecExpr<A,T,D>, Constant<T>, T> DivOpBinOp;
0412   return VecExpr<DivOpBinOp,T,D>(DivOpBinOp(DivOp<T>(),lhs,Constant<T>(mag(lhs))));
0413 }
0414 #endif
0415 
0416 
0417   }  // namespace Math
0418 
0419 }  // namespace ROOT
0420 
0421 
0422 
0423 #endif   /* ROOT_Math_Functions */