|
||||
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 */
[ Source navigation ] | [ Diff markup ] | [ Identifier search ] | [ general search ] |
This page was automatically generated by the 2.3.7 LXR engine. The LXR team |