root/Math/SVector.icc

0001 // @(#)root/smatrix:$Id$
0002 // Authors: T. Glebe, L. Moneta    2005
0003
0004 #ifndef ROOT_Math_SVector_icc
0005 #define ROOT_Math_SVector_icc
0006 // ********************************************************************
0007 //
0008 // source:
0009 //
0010 // type:      source code
0011 //
0012 // created:   21. 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: A fixed size Vector class
0023 //
0024 // changes:
0025 // 21 Mar 2001 (TG) creation
0026 // 26 Mar 2001 (TG) added place_at()
0027 // 06 Apr 2001 (TG) CTORS added
0028 // 07 Apr 2001 (TG) CTORS added
0029 // 22 Aug 2001 (TG) CTOR(T*,len) added
0030 // 14 Jan 2002 (TG) added operator==(), operator!=(), operator>(), operator<()
0031 //
0032 // ********************************************************************
0033
0034 #ifndef ROOT_Math_SVector
0035 #error "Do not use SVector.icc directly. #include \"Math/SVector.h\" instead."
0036 #endif // ROOT_Math_SVector
0037
0038 #include <iostream>
0039 #include <assert.h>
0040 #include <algorithm>
0041
0042 #include "Math/StaticCheck.h"
0043
0044 namespace ROOT {
0045
0046 namespace Math {
0047
0048
0049 //==============================================================================
0050 // Constructors
0051 //==============================================================================
0052 template <class T, unsigned int D>
0053 SVector<T,D>::SVector() {
0054    for(unsigned int i=0; i<D; ++i)
0055       fArray[i] = 0;
0056 }
0057
0058 template <class T, unsigned int D>
0059 template <class A>
0060 SVector<T,D>::SVector(const VecExpr<A,T,D>& rhs) {
0061    operator=(rhs);
0062 }
0063
0064 template <class T, unsigned int D>
0065 SVector<T,D>::SVector(const SVector<T,D>& rhs) {
0066    for(unsigned int i=0; i<D; ++i)
0067       fArray[i] = rhs.fArray[i];
0068 }
0069
0070
0071
0072 //==============================================================================
0073 // New Constructors from STL interfaces
0074 //==============================================================================
0075
0076 #ifdef LATER
0077 template <class T, unsigned int D>
0078 template <class InputIterator>
0079 SVector<T,D>::SVector(InputIterator begin, InputIterator end) {
0080    assert(begin + D == end);
0081    std::copy(begin, end, fArray);
0082 }
0083
0084 template <class T, unsigned int D>
0085 template <class InputIterator>
0086 SVector<T,D>::SVector(InputIterator begin, unsigned int size) {
0087    assert( size <= D);
0088    std::copy(begin, begin+size, fArray);
0089 }
0090
0091 #else
0092
0093 template <class T, unsigned int D>
0094 SVector<T,D>::SVector(const T* a, unsigned int len) {
0095    (void)len;
0096    assert(len == D);
0097    for(unsigned int i=0; i<D; ++i)
0098       fArray[i] = a[i];
0099 }
0100
0101 template <class T, unsigned int D>
0102 SVector<T,D>::SVector(const_iterator ibegin, const_iterator iend) {
0103    assert(ibegin + D == iend);
0104    std::copy(ibegin, iend, fArray);
0105 }
0106
0107
0108 #endif
0109
0110
0111 template <class T, unsigned int D>
0112 SVector<T,D>::SVector(const T& a1) {
0113    STATIC_CHECK( D == 1,SVector_dimension_not_right);
0114    fArray[0] = a1;
0115 }
0116
0117 template <class T, unsigned int D>
0118 SVector<T,D>::SVector(const T& a1, const T& a2) {
0119    STATIC_CHECK( D == 2,SVector_dimension_not_right);
0120    fArray[0] = a1; fArray[1] = a2;
0121 }
0122
0123 template <class T, unsigned int D>
0124 SVector<T,D>::SVector(const T& a1, const T& a2, const T& a3) {
0125    STATIC_CHECK( D == 3,SVector_dimension_not_right);
0126    fArray[0] = a1; fArray[1] = a2; fArray[2] = a3;
0127 }
0128
0129 template <class T, unsigned int D>
0130 SVector<T,D>::SVector(const T& a1, const T& a2, const T& a3, const T& a4) {
0131    STATIC_CHECK( D == 4,SVector_dimension_not_right);
0132    fArray[0] = a1; fArray[1] = a2; fArray[2] = a3; fArray[3] = a4;
0133 }
0134
0135 template <class T, unsigned int D>
0136 SVector<T,D>::SVector(const T& a1, const T& a2, const T& a3, const T& a4,
0137                       const T& a5) {
0138    STATIC_CHECK( D == 5,SVector_dimension_not_right);
0139    fArray[0] = a1; fArray[1] = a2; fArray[2] = a3; fArray[3] = a4;
0140    fArray[4] = a5;
0141 }
0142
0143 template <class T, unsigned int D>
0144 SVector<T,D>::SVector(const T& a1, const T& a2, const T& a3, const T& a4,
0145                       const T& a5, const T& a6) {
0146    STATIC_CHECK( D == 6,SVector_dimension_not_right);
0147    fArray[0] = a1; fArray[1] = a2; fArray[2] = a3; fArray[3] = a4;
0148    fArray[4] = a5; fArray[5] = a6;
0149 }
0150
0151 template <class T, unsigned int D>
0152 SVector<T,D>::SVector(const T& a1, const T& a2, const T& a3, const T& a4,
0153                       const T& a5, const T& a6, const T& a7) {
0154    STATIC_CHECK( D == 7,SVector_dimension_not_right);
0155    fArray[0] = a1; fArray[1] = a2; fArray[2] = a3; fArray[3] = a4;
0156    fArray[4] = a5; fArray[5] = a6; fArray[6] = a7;
0157 }
0158
0159 template <class T, unsigned int D>
0160 SVector<T,D>::SVector(const T& a1, const T& a2, const T& a3, const T& a4,
0161                       const T& a5, const T& a6, const T& a7, const T& a8) {
0162    STATIC_CHECK( D == 8,SVector_dimension_not_right);
0163    fArray[0] = a1; fArray[1] = a2; fArray[2] = a3; fArray[3] = a4;
0164    fArray[4] = a5; fArray[5] = a6; fArray[6] = a7; fArray[7] = a8;
0165 }
0166
0167 template <class T, unsigned int D>
0168 SVector<T,D>::SVector(const T& a1, const T& a2, const T& a3, const T& a4,
0169                       const T& a5, const T& a6, const T& a7, const T& a8,
0170                       const T& a9) {
0171    STATIC_CHECK( D == 9,SVector_dimension_not_right);
0172    fArray[0] = a1; fArray[1] = a2; fArray[2] = a3; fArray[3] = a4;
0173    fArray[4] = a5; fArray[5] = a6; fArray[6] = a7; fArray[7] = a8;
0174    fArray[8] = a9;
0175 }
0176
0177 template <class T, unsigned int D>
0178 SVector<T,D>::SVector(const T& a1, const T& a2, const T& a3, const T& a4,
0179                       const T& a5, const T& a6, const T& a7, const T& a8,
0180                       const T& a9, const T& a10) {
0181    STATIC_CHECK( D == 10,SVector_dimension_not_right);
0182    fArray[0] = a1; fArray[1] = a2; fArray[2] = a3; fArray[3] = a4;
0183    fArray[4] = a5; fArray[5] = a6; fArray[6] = a7; fArray[7] = a8;
0184    fArray[8] = a9; fArray[9] = a10;
0185 }
0186
0187 //==============================================================================
0188 // operator=
0189 //==============================================================================
0190 template <class T, unsigned int D>
0191 SVector<T,D>& SVector<T,D>::operator=(const T& a1) {
0192    // operator = for size 1 vectors
0193    STATIC_CHECK( D == 1,SVector_dimension_not_right);
0194    fArray[0] = a1;
0195    return *this;
0196 }
0197
0198 template <class T, unsigned int D>
0199 SVector<T,D>& SVector<T,D>::operator=(const SVector<T,D>& rhs) {
0200    for(unsigned int i=0; i<D; ++i)
0201       fArray[i] = rhs.fArray[i];
0202    return *this;
0203 }
0204
0205 template <class T, unsigned int D>
0206 template <class A>
0207 SVector<T,D>& SVector<T,D>::operator=(const VecExpr<A,T,D>& rhs) {
0208    if (! rhs.IsInUse(fArray) ) {
0209       for(unsigned int i=0; i<D; ++i) {
0210          fArray[i] = rhs.apply(i);
0211       }
0212    }
0213    else {
0214       // otherwise we need to create a temporary object
0215       T tmp[D];
0216       for(unsigned int i=0; i<D; ++i) {
0217          tmp[i] = rhs.apply(i);
0218       }
0219       for(unsigned int i=0; i<D; ++i) {
0220          fArray[i] = tmp[i];
0221       }
0222    }
0223    return *this;
0224 }
0225
0226 //==============================================================================
0227 // operator==
0228 //==============================================================================
0229 template <class T, unsigned int D>
0230 bool SVector<T,D>::operator==(const T& rhs) const {
0231    bool rc = true;
0232    for(unsigned int i=0; i<D; ++i) {
0233       rc = rc && (fArray[i] == rhs);
0234    }
0235    return rc;
0236 }
0237
0238 template <class T, unsigned int D>
0239 bool SVector<T,D>::operator==(const SVector<T,D>& rhs) const {
0240    bool rc = true;
0241    for(unsigned int i=0; i<D; ++i) {
0242       rc = rc && (fArray[i] == rhs.apply(i));
0243    }
0244    return rc;
0245 }
0246
0247 template <class T, unsigned int D>
0248 template <class A>
0249 bool SVector<T,D>::operator==(const VecExpr<A,T,D>& rhs) const {
0250    bool rc = true;
0251    for(unsigned int i=0; i<D; ++i) {
0252       rc = rc && (fArray[i] == rhs.apply(i));
0253    }
0254    return rc;
0255 }
0256
0257 //==============================================================================
0258 // operator!=
0259 //==============================================================================
0260 template <class T, unsigned int D>
0261 inline bool SVector<T,D>::operator!=(const T& rhs) const {
0262    return !operator==(rhs);
0263 }
0264
0265 template <class T, unsigned int D>
0266 inline bool SVector<T,D>::operator!=(const SVector<T,D>& rhs) const {
0267    return !operator==(rhs);
0268 }
0269
0270 template <class T, unsigned int D>
0271 template <class A>
0272 inline bool SVector<T,D>::operator!=(const VecExpr<A,T,D>& rhs) const {
0273    return !operator==(rhs);
0274 }
0275
0276 //==============================================================================
0277 // operator>
0278 //==============================================================================
0279 template <class T, unsigned int D>
0280 bool SVector<T,D>::operator>(const T& rhs) const {
0281    bool rc = true;
0282    for(unsigned int i=0; i<D; ++i) {
0283       rc = rc && (fArray[i] > rhs);
0284    }
0285    return rc;
0286 }
0287
0288 template <class T, unsigned int D>
0289 bool SVector<T,D>::operator>(const SVector<T,D>& rhs) const {
0290    bool rc = true;
0291    for(unsigned int i=0; i<D; ++i) {
0292       rc = rc && (fArray[i] > rhs.apply(i));
0293    }
0294    return rc;
0295 }
0296
0297 template <class T, unsigned int D>
0298 template <class A>
0299 bool SVector<T,D>::operator>(const VecExpr<A,T,D>& rhs) const {
0300    bool rc = true;
0301    for(unsigned int i=0; i<D; ++i) {
0302       rc = rc && (fArray[i] > rhs.apply(i));
0303    }
0304    return rc;
0305 }
0306
0307 //==============================================================================
0308 // operator<
0309 //==============================================================================
0310 template <class T, unsigned int D>
0311 bool SVector<T,D>::operator<(const T& rhs) const {
0312    bool rc = true;
0313    for(unsigned int i=0; i<D; ++i) {
0314       rc = rc && (fArray[i] < rhs);
0315    }
0316    return rc;
0317 }
0318
0319 template <class T, unsigned int D>
0320 bool SVector<T,D>::operator<(const SVector<T,D>& rhs) const {
0321    bool rc = true;
0322    for(unsigned int i=0; i<D; ++i) {
0323       rc = rc && (fArray[i] < rhs.apply(i));
0324    }
0325    return rc;
0326 }
0327
0328 template <class T, unsigned int D>
0329 template <class A>
0330 bool SVector<T,D>::operator<(const VecExpr<A,T,D>& rhs) const {
0331    bool rc = true;
0332    for(unsigned int i=0; i<D; ++i) {
0333       rc = rc && (fArray[i] < rhs.apply(i));
0334    }
0335    return rc;
0336 }
0337
0338 //==============================================================================
0339 // operator+=
0340 //==============================================================================
0341 #ifdef NEW_IMPL
0342 template <class T, unsigned int D>
0343 template<class A>
0344 SVector<T,D>& SVector<T,D>::operator+=(const  A& rhs) {
0345    return operator=(*this + rhs);
0346 }
0347
0348 template <class T, unsigned int D>
0349 template<class A>
0350 SVector<T,D>& SVector<T,D>::operator-=(const  A& rhs) {
0351    // self subtraction
0352    return operator=(*this - rhs);
0353 }
0354
0355 template <class T, unsigned int D>
0356 template<class A>
0357 SVector<T,D>& SVector<T,D>::operator*=(const  A& rhs) {
0358    // self multiplication
0359    return operator=(*this * rhs);
0360 }
0361
0362 template <class T, unsigned int D>
0363 template<class A>
0364 SVector<T,D>& SVector<T,D>::operator/=(const  A& rhs) {
0365    // self division
0366    return operator=(*this / rhs);
0367 }
0368 #endif
0369
0370 template <class T, unsigned int D>
0371 SVector<T,D>& SVector<T,D>::operator+=(const  T& rhs) {
0372    for(unsigned int i=0; i<D; ++i) {
0373       fArray[i] += rhs;
0374    }
0375    return *this;
0376 }
0377
0378 template <class T, unsigned int D>
0379 SVector<T,D>& SVector<T,D>::operator+=(const  SVector<T,D>& rhs) {
0380    for(unsigned int i=0; i<D; ++i) {
0381       fArray[i] += rhs.apply(i);
0382    }
0383    return *this;
0384 }
0385
0386
0387 template <class T, unsigned int D>
0388 template <class A>
0389 SVector<T,D>& SVector<T,D>::operator+=(const VecExpr<A,T,D>& rhs) {
0390    for(unsigned int i=0; i<D; ++i) {
0391       fArray[i] += rhs.apply(i);
0392    }
0393    return *this;
0394 }
0395
0396 //==============================================================================
0397 // operator-=
0398 //==============================================================================
0399 template <class T, unsigned int D>
0400 SVector<T,D>& SVector<T,D>::operator-=(const  T& rhs) {
0401    for(unsigned int i=0; i<D; ++i) {
0402       fArray[i] -= rhs;
0403    }
0404    return *this;
0405 }
0406
0407 template <class T, unsigned int D>
0408 SVector<T,D>& SVector<T,D>::operator-=(const SVector<T,D>& rhs) {
0409    for(unsigned int i=0; i<D; ++i) {
0410       fArray[i] -= rhs.apply(i);
0411    }
0412    return *this;
0413 }
0414
0415 template <class T, unsigned int D>
0416 template <class A>
0417 SVector<T,D>& SVector<T,D>::operator-=(const VecExpr<A,T,D>& rhs) {
0418    for(unsigned int i=0; i<D; ++i) {
0419       fArray[i] -= rhs.apply(i);
0420    }
0421    return *this;
0422 }
0423
0424 //==============================================================================
0425 // operator*= (only scalar values)
0426 //==============================================================================
0427 template <class T, unsigned int D>
0428 SVector<T,D>& SVector<T,D>::operator*=(const T& rhs) {
0429    for(unsigned int i=0; i<D; ++i) {
0430       fArray[i] *= rhs;
0431    }
0432    return *this;
0433 }
0434 #ifdef OLD_IMPL
0435 template <class T, unsigned int D>
0436 template <class A>
0437 SVector<T,D>& SVector<T,D>::operator*=(const VecExpr<A,T,D>& rhs) {
0438    for(unsigned int i=0; i<D; ++i) {
0439       fArray[i] *= rhs.apply(i);
0440    }
0441    return *this;
0442 }
0443
0444 //==============================================================================
0445 // operator/=
0446 //==============================================================================
0447 template <class T, unsigned int D>
0448 SVector<T,D>& SVector<T,D>::operator/=(const SVector<T,D>& rhs) {
0449    for(unsigned int i=0; i<D; ++i) {
0450       fArray[i] /= rhs.apply(i);
0451    }
0452    return *this;
0453 }
0454
0455 template <class T, unsigned int D>
0456 template <class A>
0457 SVector<T,D>& SVector<T,D>::operator/=(const VecExpr<A,T,D>& rhs) {
0458    for(unsigned int i=0; i<D; ++i) {
0459       fArray[i] /= rhs.apply(i);
0460    }
0461    return *this;
0462 }
0463 #endif
0464 template <class T, unsigned int D>
0465 SVector<T,D>& SVector<T,D>::operator/=(const T& rhs) {
0466    for(unsigned int i=0; i<D; ++i) {
0467       fArray[i] /= rhs;
0468    }
0469    return *this;
0470 }
0471
0472
0473 //==============================================================================
0474 // unit
0475 //==============================================================================
0476 template <class T, unsigned int D>
0477 SVector<T,D>& SVector<T,D>::Unit() {
0478    const T len = Mag(*this);
0479    for(unsigned int i=0; i<D; ++i) {
0480       fArray[i] /= len;
0481    }
0482    return *this;
0483 }
0484
0485 //==============================================================================
0486 // place_at
0487 //==============================================================================
0488 template <class T, unsigned int D>
0489 template <unsigned int D2>
0490 SVector<T,D>& SVector<T,D>::Place_at(const SVector<T,D2>& rhs, unsigned int row) {
0491
0492    assert(row+D2 <= D);
0493    //  Sassert(end <= D);
0494
0495    for(unsigned int i=row, j=0; j<D2; ++i,++j)
0496       fArray[i] = rhs.apply(j);
0497
0498    return *this;
0499 }
0500
0501
0502 //==============================================================================
0503 // place_at
0504 //==============================================================================
0505 template <class T, unsigned int D>
0506 template <class A, unsigned int D2>
0507 SVector<T,D>& SVector<T,D>::Place_at(const VecExpr<A,T,D2>& rhs, unsigned int row) {
0508
0509    assert(row+D2 <= D);
0510
0511    for(unsigned int i=row, j=0; j<D2; ++i,++j)
0512       fArray[i] = rhs.apply(j);
0513
0514    return *this;
0515 }
0516
0517 //==============================================================================
0518 // print
0519 //==============================================================================
0520 template <class T, unsigned int D>
0521 std::ostream& SVector<T,D>::Print(std::ostream& os) const {
0522    const std::ios_base::fmtflags prevFmt = os.setf(std::ios::right,std::ios::adjustfield);
0523    //  os.setf(ios::fixed);
0524
0525    for (unsigned int i = 0; i < D; ++i ) {
0526       os << fArray[i];
0527       if (i != D-1) os << ", ";
0528    }
0529    if (prevFmt != os.flags() ) os.setf(prevFmt, std::ios::adjustfield);
0530    return os;
0531 }
0532
0533 //==============================================================================
0534 // Access functions
0535 //==============================================================================
0536 template <class T, unsigned int D>
0537 inline T SVector<T,D>::apply(unsigned int i) const { return fArray[i]; }
0538
0539 template <class T, unsigned int D>
0540 inline const T* SVector<T,D>::Array() const { return fArray; }
0541
0542 template <class T, unsigned int D>
0543 inline T* SVector<T,D>::Array() { return fArray; }
0544
0545
0546 //==============================================================================
0547 // STL interface
0548 //==============================================================================
0549 template <class T, unsigned int D>
0550 inline T* SVector<T,D>::begin() { return fArray; }
0551
0552 template <class T, unsigned int D>
0553 inline const T* SVector<T,D>::begin() const { return fArray; }
0554
0555 template <class T, unsigned int D>
0556 inline T* SVector<T,D>::end() { return fArray + Dim(); }
0557
0558 template <class T, unsigned int D>
0559 inline const T* SVector<T,D>::end() const { return fArray + Dim(); }
0560
0561 template <class T, unsigned int D>
0562 template <class InputIterator>
0563 void SVector<T,D>::SetElements(InputIterator ibegin, InputIterator iend) {
0564    // iterator size must match vector size
0565    assert( ibegin + D == iend);
0566    std::copy(ibegin, iend, fArray);
0567 }
0568
0569 template <class T, unsigned int D>
0570 template <class InputIterator>
0571 void SVector<T,D>::SetElements(InputIterator ibegin, unsigned int size) {
0572    // size <= vector size
0573    assert( size <= D);
0574    std::copy(ibegin, ibegin+size, fArray);
0575 }
0576
0577
0578 //==============================================================================
0579 // Operators
0580 //==============================================================================
0581 template <class T, unsigned int D>
0582 inline const T& SVector<T,D>::operator[](unsigned int i) const { return fArray[i]; }
0583
0584 template <class T, unsigned int D>
0585 inline const T& SVector<T,D>::operator()(unsigned int i) const { return fArray[i]; }
0586
0587 template <class T, unsigned int D>
0588 inline T& SVector<T,D>::operator[](unsigned int i) { return fArray[i]; }
0589
0590 template <class T, unsigned int D>
0591 inline T& SVector<T,D>::operator()(unsigned int i) { return fArray[i]; }
0592 //==============================================================================
0593 // Element access with At()
0594 //==============================================================================
0595 template <class T, unsigned int D>
0596 inline const T& SVector<T,D>::At(unsigned int i) const {
0597    assert(i < D);
0598    return fArray[i];
0599 }
0600
0601 template <class T, unsigned int D>
0602 inline T& SVector<T,D>::At(unsigned int i) {
0603    assert(i < D);
0604    return fArray[i];
0605 }
0606
0607 //==============================================================================
0608 // SubVector
0609 //==============================================================================
0610 template <class T, unsigned int D>
0611 template <class SubVector>
0612 SubVector SVector<T,D>::Sub(unsigned int row) const {
0613
0614    STATIC_CHECK( SubVector::kSize <= D,SVector_dimension_too_small);
0615
0616    assert(row + SubVector::kSize <= D);
0617
0618    SubVector tmp;
0619    // need to use std::copy ??
0620    for(unsigned int i=0; i < SubVector::kSize; ++i) {
0621       tmp[i] = fArray[i+row];
0622    }
0623    return tmp;
0624 }
0625
0626 // check if the given passed pointer is the same contained in the vector
0627 template <class T, unsigned int D>
0628 bool SVector<T,D>::IsInUse( const T * p) const {
0629    return p == fArray;
0630 }
0631
0632
0633 //==============================================================================
0634 // operator<<
0635 //==============================================================================
0636 template <class T, unsigned int D>
0637 inline std::ostream& operator<<(std::ostream& os, const SVector<T,D>& rhs) {
0638    return rhs.Print(os);
0639 }
0640
0641
0642
0643
0644 }  // namespace Math
0645
0646 }  // namespace ROOT
0647
0648
0649 #endif