/*************************************************************************** * * $Id$ * * Author: Brian Lasiuk, Thomas Ullrich, April 1998 *************************************************************************** * * Description: * * Remarks: Since not all compilers support member templates * we have to specialize the templated member on these * platforms. If member templates are not supported the * ST_NO_MEMBER_TEMPLATES flag has to be set. tu. * *************************************************************************** * * $Log$ * Revision 1.1.1.1 2007/03/07 10:14:49 mchojnacki * First version on CVS * * Revision 1.15 2005/09/22 20:09:20 fisyak * Make AliFmLorentzVector persistent * * Revision 1.14 2005/07/19 22:27:11 perev * Cleanup * * Revision 1.13 2005/07/06 18:49:57 fisyak * Replace AliFmHelixD, AliFmLorentzVectorD,AliFmLorentzVectorF,AliFmMatrixD,AliFmMatrixF,AliFmPhysicalHelixD,AliFmThreeVectorD,AliFmThreeVectorF by templated version * * Revision 1.12 2005/03/28 06:03:41 perev * Defence FPE added * * Revision 1.11 2004/12/02 20:07:32 fine * define the valid method for both flavor of AliFmThreeVector * * Revision 1.10 2003/10/30 20:06:46 perev * Check of quality added * * Revision 1.9 2003/09/02 17:59:35 perev * gcc 3.2 updates + WarnOff * * Revision 1.8 2002/06/21 17:47:37 genevb * Added pseudoProduct * * Revision 1.7 2000/01/04 19:56:05 ullrich * Added cpp macro for CINT. * * Revision 1.6 1999/12/21 15:14:31 ullrich * Modified to cope with new compiler version on Sun (CC5.0). * * Revision 1.5 1999/10/15 15:46:54 ullrich * Changed output format in operator<< * * Revision 1.4 1999/06/04 18:00:05 ullrich * Added new constructor which takes C-style array as argument. * New operators operator() and operator[] which can be used * as lvalues. * * Revision 1.3 1999/02/17 11:42:19 ullrich * Removed specialization for 'long double'. * * Revision 1.2 1999/02/14 23:11:48 fisyak * Fixes for Rootcint * * Revision 1.1 1999/01/30 03:59:05 fisyak * Root Version of AliFmarClassLibrary * * Revision 1.1 1999/01/23 00:28:04 ullrich * Initial Revision * **************************************************************************/ #ifndef ST_THREE_VECTOR_HH #define ST_THREE_VECTOR_HH #ifdef __ROOT__ #include "Rtypes.h" #endif #ifndef __CINT__ #include #include #include #ifdef GNU_GCC # include #endif #if defined (__SUNPRO_CC) && __SUNPRO_CC < 0x500 # include #endif #ifndef ST_NO_EXCEPTIONS # include # if !defined(ST_NO_NAMESPACES) using std::out_of_range; # endif #endif #endif // __CINT__ #ifdef WIN32 #include "gcc2vs.h" #endif class TRootIOCtor;//nic nie rozumiem using namespace std; template class AliFmThreeVector { public: AliFmThreeVector(T = 0, T = 0, T = 0); // ROOT_VERSION(5,03,01) #if ROOT_VERSION_CODE >= 328449 AliFmThreeVector(TRootIOCtor*) : mX1(0), mX2(0), mX3(0) {} #endif virtual ~AliFmThreeVector(); #ifndef ST_NO_MEMBER_TEMPLATES template AliFmThreeVector(const AliFmThreeVector&); template AliFmThreeVector(const X*); template AliFmThreeVector& operator=(const AliFmThreeVector&); // AliFmThreeVector(const AliFmThreeVector&); use default // AliFmThreeVector& operator=(const AliFmThreeVector&); use default #else AliFmThreeVector(const AliFmThreeVector&); AliFmThreeVector(const AliFmThreeVector&); AliFmThreeVector(const float*); AliFmThreeVector(const double*); AliFmThreeVector& operator=(const AliFmThreeVector&); AliFmThreeVector& operator=(const AliFmThreeVector&); #endif void setX(T); void setY(T); void setZ(T); void setPhi(T); void setTheta(T); void setMag(T); void setMagnitude(T); T x() const; T y() const; T z() const; T theta() const; T cosTheta() const; T phi() const; T perp() const; T perp2() const; T magnitude() const; T mag() const; T mag2() const; T pseudoRapidity() const; T operator() (size_t) const; T operator[] (size_t) const; T& operator() (size_t); T& operator[] (size_t); T massHypothesis(T mass) const; AliFmThreeVector unit() const; AliFmThreeVector orthogonal() const; void rotateX(T); void rotateY(T); void rotateZ(T); AliFmThreeVector operator- (); AliFmThreeVector operator+ (); AliFmThreeVector& operator*= (double); AliFmThreeVector& operator/= (double); AliFmThreeVector pseudoProduct(double,double,double) const; #ifndef ST_NO_MEMBER_TEMPLATES template T angle(const AliFmThreeVector&) const; template AliFmThreeVector cross(const AliFmThreeVector&) const; template T dot (const AliFmThreeVector&) const; template AliFmThreeVector pseudoProduct(const AliFmThreeVector&) const; template bool operator == (const AliFmThreeVector& v) const; template bool operator != (const AliFmThreeVector& v) const; template AliFmThreeVector& operator+= (const AliFmThreeVector&); template AliFmThreeVector& operator-= (const AliFmThreeVector&); #else T angle(const AliFmThreeVector&) const; AliFmThreeVector cross(const AliFmThreeVector&) const; T dot (const AliFmThreeVector&) const; AliFmThreeVector pseudoProduct(const AliFmThreeVector&) const; T angle(const AliFmThreeVector&) const; T dot (const AliFmThreeVector&) const; AliFmThreeVector cross(const AliFmThreeVector&) const; AliFmThreeVector pseudoProduct(const AliFmThreeVector&) const; bool operator == (const AliFmThreeVector& v) const; bool operator != (const AliFmThreeVector& v) const; AliFmThreeVector& operator+= (const AliFmThreeVector&); AliFmThreeVector& operator-= (const AliFmThreeVector&); bool operator == (const AliFmThreeVector& v) const; bool operator != (const AliFmThreeVector& v) const; AliFmThreeVector& operator+= (const AliFmThreeVector&); AliFmThreeVector& operator-= (const AliFmThreeVector&); #endif int valid(double world = 1.e+5) const; int bad(double world = 1.e+5) const; protected: T mX1, mX2, mX3; #ifdef __ROOT__ ClassDef(AliFmThreeVector,3) #endif /* __ROOT__ */ }; #ifndef __CINT__ // // Implementation of member functions // template inline AliFmThreeVector::AliFmThreeVector(T x, T y, T z) : mX1(x), mX2(y), mX3(z) {/* nop */} template inline AliFmThreeVector::~AliFmThreeVector() {/* nop */} template inline void AliFmThreeVector::setX(T x) {mX1 = x;} template inline void AliFmThreeVector::setY(T y) {mX2 = y;} template inline void AliFmThreeVector::setZ(T z) {mX3 = z;} template void AliFmThreeVector::setPhi(T angle) { double r = magnitude(); double th = theta(); mX1 = r*sin(th)*cos(angle); mX2 = r*sin(th)*sin(angle); } template void AliFmThreeVector::setTheta(T angle) { double r = magnitude(); double ph = phi(); mX1 = r*sin(angle)*cos(ph); mX2 = r*sin(angle)*sin(ph); mX3 = r*cos(angle); } template void AliFmThreeVector::setMagnitude(T r) { double th = theta(); double ph = phi(); mX1 = r*sin(th)*cos(ph); mX2 = r*sin(th)*sin(ph); mX3 = r*cos(th); } template void AliFmThreeVector::setMag(T mag) { setMagnitude(mag); } template inline T AliFmThreeVector::x() const {return mX1;} template inline T AliFmThreeVector::y() const {return mX2;} template inline T AliFmThreeVector::z() const {return mX3;} template inline T AliFmThreeVector::theta() const { return acos(cosTheta()); } template inline T AliFmThreeVector::cosTheta() const { return mX3/(mag()+1e-20); } template inline T AliFmThreeVector::phi() const { return atan2(mX2,mX1); } template inline T AliFmThreeVector::pseudoRapidity() const { // // change code to more optimal: // double m = mag(); // return 0.5*::log( (m+z())/(m-z()) ); double tmp = tan(theta()/2.); if (tmp <=0.) return 1e20; return -::log(tmp); } template inline AliFmThreeVector AliFmThreeVector::unit() const { double tmp = mag(); if (tmp<=0.) tmp = 1e-20; return *this/tmp; } template T AliFmThreeVector::massHypothesis(T mass) const { return ::sqrt((*this)*(*this) + mass*mass); } template AliFmThreeVector AliFmThreeVector::orthogonal() const { // Direct copy from CLHEP--it is probably better to // use your own dot/cross product code... double x = (mX1 < 0.0) ? -mX1 : mX1; double y = (mX2 < 0.0) ? -mX2 : mX2; double z = (mX3 < 0.0) ? -mX3 : mX3; if(x(0,mX3,-mX2) : AliFmThreeVector(mX2,-mX1,0); else return mX2 < mX3 ? AliFmThreeVector(-mX3,0,mX1) : AliFmThreeVector(mX2,-mX1,0); } template void AliFmThreeVector::rotateX(T angle) { // may in the future make use of the AliFmRotation class! double yPrime = cos(angle)*mX2 - sin(angle)*mX3; double zPrime = sin(angle)*mX2 + cos(angle)*mX3; mX2 = yPrime; mX3 = zPrime; } template void AliFmThreeVector::rotateY(T angle) { // may in the future make use of the AliFmRotation class! double zPrime = cos(angle)*mX3 - sin(angle)*mX1; double xPrime = sin(angle)*mX3 + cos(angle)*mX1; mX1 = xPrime; mX3 = zPrime; } template void AliFmThreeVector::rotateZ(T angle) { // may in the future make use of the AliFmRotation class! double xPrime = cos(angle)*mX1 - sin(angle)*mX2; double yPrime = sin(angle)*mX1 + cos(angle)*mX2; mX1 = xPrime; mX2 = yPrime; } template inline T AliFmThreeVector::perp() const { return ::sqrt(mX1*mX1+mX2*mX2); } template inline T AliFmThreeVector::perp2() const { return mX1*mX1+mX2*mX2; } template inline T AliFmThreeVector::magnitude() const { return mag(); } template inline T AliFmThreeVector::mag() const { return ::sqrt(mX1*mX1+mX2*mX2+mX3*mX3); } template inline T AliFmThreeVector::mag2() const { return mX1*mX1+mX2*mX2+mX3*mX3; } template inline T AliFmThreeVector::operator() (size_t i) const { if (0 <=i && i <= 2) return (&mX1)[i]; #ifndef ST_NO_EXCEPTIONS throw out_of_range("AliFmThreeVector::operator(): bad index"); #else cerr << "AliFmThreeVector::operator(): bad index" << endl; #endif return 0; } template inline T& AliFmThreeVector::operator() (size_t i) { if (0 <=i && i <= 2) return (&mX1)[i]; #ifndef ST_NO_EXCEPTIONS throw out_of_range("AliFmThreeVector::operator(): bad index"); #else cerr << "AliFmThreeVector::operator(): bad index" << endl; #endif return mX1; } template inline T AliFmThreeVector::operator[] (size_t i) const { if (0 <=i && i <= 2) return (&mX1)[i]; #ifndef ST_NO_EXCEPTIONS throw out_of_range("AliFmThreeVector::operator[]: bad index"); #else cerr << "AliFmThreeVector::operator[]: bad index" << endl; #endif return 0; } template inline T &AliFmThreeVector::operator[] (size_t i) { if (0 <=i && i <= 2) return (&mX1)[i]; #ifndef ST_NO_EXCEPTIONS throw out_of_range("AliFmThreeVector::operator[]: bad index"); #else cerr << "AliFmThreeVector::operator[]: bad index" << endl; #endif return mX1; } template inline AliFmThreeVector& AliFmThreeVector::operator*= (double c) { mX1 *= c; mX2 *= c; mX3 *= c; return *this; } template inline AliFmThreeVector& AliFmThreeVector::operator/= (double c) { mX1 /= c; mX2 /= c; mX3 /= c; return *this; } template inline AliFmThreeVector AliFmThreeVector::pseudoProduct(double x,double y,double z) const { return AliFmThreeVector(mX1*x,mX2*y,mX3*z); } template AliFmThreeVector AliFmThreeVector::operator- () { return AliFmThreeVector(-mX1, -mX2, -mX3); } template AliFmThreeVector AliFmThreeVector::operator+ () { return *this; } #ifndef ST_NO_MEMBER_TEMPLATES #ifndef WIN32 template template inline AliFmThreeVector::AliFmThreeVector(const AliFmThreeVector& v) : mX1(v.x()), mX2(v.y()), mX3(v.z()) {/* nop */} template template inline AliFmThreeVector::AliFmThreeVector(const X *a) { mX1 = a[0]; mX2 = a[1]; mX3 = a[2]; } template template inline AliFmThreeVector& AliFmThreeVector::operator=(const AliFmThreeVector& v) { mX1 = v.x(); mX2 = v.y(); mX3 = v.z(); return *this; } template template inline bool AliFmThreeVector::operator== (const AliFmThreeVector& v) const { return mX1 == v.x() && mX2 == v.y() && mX3 == v.z(); } template template inline bool AliFmThreeVector::operator!= (const AliFmThreeVector& v) const { return !(*this == v); } template template inline AliFmThreeVector& AliFmThreeVector::operator+= (const AliFmThreeVector& v) { mX1 += v.x(); mX2 += v.y(); mX3 += v.z(); return *this; } template template inline AliFmThreeVector& AliFmThreeVector::operator-= (const AliFmThreeVector& v) { mX1 -= v.x(); mX2 -= v.y(); mX3 -= v.z(); return *this; } template template inline T AliFmThreeVector::dot(const AliFmThreeVector& v) const { return mX1*v.x() + mX2*v.y() + mX3*v.z(); } template template inline AliFmThreeVector AliFmThreeVector::cross(const AliFmThreeVector& v) const { return AliFmThreeVector(mX2*v.z() - mX3*v.y(), mX3*v.x() - mX1*v.z(), mX1*v.y() - mX2*v.x()); } template template inline T AliFmThreeVector::angle(const AliFmThreeVector& vec) const { double norm = this->mag2()*vec.mag2(); return norm > 0 ? acos(this->dot(vec)/(::sqrt(norm))) : 0; } template template inline AliFmThreeVector AliFmThreeVector::pseudoProduct(const AliFmThreeVector& v) const { return this->pseudoProduct(v.x(),v.y(),v.z()); } #endif #else template inline AliFmThreeVector::AliFmThreeVector(const AliFmThreeVector& v) : mX1(v.x()), mX2(v.y()), mX3(v.z()) {/* nop */} template inline AliFmThreeVector::AliFmThreeVector(const AliFmThreeVector& v) : mX1(v.x()), mX2(v.y()), mX3(v.z()) {/* nop */} template inline AliFmThreeVector::AliFmThreeVector(const float *a) { mX1 = a[0]; mX2 = a[1]; mX3 = a[2]; } template inline AliFmThreeVector::AliFmThreeVector(const double *a) { mX1 = a[0]; mX2 = a[1]; mX3 = a[2]; } template inline AliFmThreeVector& AliFmThreeVector::operator=(const AliFmThreeVector& v) { mX1 = v.x(); mX2 = v.y(); mX3 = v.z(); return *this; } template inline AliFmThreeVector& AliFmThreeVector::operator=(const AliFmThreeVector& v) { mX1 = v.x(); mX2 = v.y(); mX3 = v.z(); return *this; } template inline bool AliFmThreeVector::operator== (const AliFmThreeVector& v) const { return mX1 == v.x() && mX2 == v.y() && mX3 == v.z(); } template inline bool AliFmThreeVector::operator== (const AliFmThreeVector& v) const { return mX1 == v.x() && mX2 == v.y() && mX3 == v.z(); } template inline bool AliFmThreeVector::operator!= (const AliFmThreeVector& v) const { return !(*this == v); } template inline bool AliFmThreeVector::operator!= (const AliFmThreeVector& v) const { return !(*this == v); } template inline AliFmThreeVector& AliFmThreeVector::operator+= (const AliFmThreeVector& v) { mX1 += v.x(); mX2 += v.y(); mX3 += v.z(); return *this; } template inline AliFmThreeVector& AliFmThreeVector::operator+= (const AliFmThreeVector& v) { mX1 += v.x(); mX2 += v.y(); mX3 += v.z(); return *this; } template inline AliFmThreeVector& AliFmThreeVector::operator-= (const AliFmThreeVector& v) { mX1 -= v.x(); mX2 -= v.y(); mX3 -= v.z(); return *this; } template inline AliFmThreeVector& AliFmThreeVector::operator-= (const AliFmThreeVector& v) { mX1 -= v.x(); mX2 -= v.y(); mX3 -= v.z(); return *this; } template inline T AliFmThreeVector::dot(const AliFmThreeVector& v) const { return mX1*v.x() + mX2*v.y() + mX3*v.z(); } template inline T AliFmThreeVector::dot(const AliFmThreeVector& v) const { return mX1*v.x() + mX2*v.y() + mX3*v.z(); } template inline AliFmThreeVector AliFmThreeVector::cross(const AliFmThreeVector& v) const { return AliFmThreeVector(mX2*v.z() - mX3*v.y(), mX3*v.x() - mX1*v.z(), mX1*v.y() - mX2*v.x()); } template inline AliFmThreeVector AliFmThreeVector::cross(const AliFmThreeVector& v) const { return AliFmThreeVector(mX2*v.z() - mX3*v.y(), mX3*v.x() - mX1*v.z(), mX1*v.y() - mX2*v.x()); } template inline T AliFmThreeVector::angle(const AliFmThreeVector& v) const { double tmp = mag()*v.mag(); if (tmp <=0) tmp = 1e-20; return acos(this->dot(v)/tmp); } template inline T AliFmThreeVector::angle(const AliFmThreeVector& v) const { double tmp = mag()*v.mag(); if (tmp <=0) tmp = 1e-20; return acos(this->dot(v)/tmp); } template inline AliFmThreeVector AliFmThreeVector::pseudoProduct(const AliFmThreeVector& v) const { return this->pseudoProduct(v.x(),v.y(),v.z()); } template inline AliFmThreeVector AliFmThreeVector::pseudoProduct(const AliFmThreeVector& v) const { return this->pseudoProduct(v.x(),v.y(),v.z()); } #endif // ST_NO_MEMBER_TEMPLATES template inline int AliFmThreeVector::valid(double world) const {return !bad(world);} template inline int AliFmThreeVector::bad(double world) const { for (int i=0;i<3;i++) { if (!finite((&mX1)[i]) ) return 10+i; if ( fabs ((&mX1)[i])>world) return 20+i; } return 0; } #endif /*! __CINT__ */ #ifdef __CINT__ template<> float abs(const AliFmThreeVector& v); template<> double abs(const AliFmThreeVector& v); template<> AliFmThreeVector cross_product(const AliFmThreeVector& v1, const AliFmThreeVector& v2); template<> AliFmThreeVector cross_product(const AliFmThreeVector& v1, const AliFmThreeVector& v2); template<> AliFmThreeVector cross_product(const AliFmThreeVector& v1, const AliFmThreeVector& v2); template<> AliFmThreeVector cross_product(const AliFmThreeVector& v1, const AliFmThreeVector& v2); template<> AliFmThreeVector operator+ (const AliFmThreeVector& v1, const AliFmThreeVector& v2); template<> AliFmThreeVector operator+ (const AliFmThreeVector& v1, const AliFmThreeVector& v2); template<> AliFmThreeVector operator+ (const AliFmThreeVector& v1, const AliFmThreeVector& v2); template<> AliFmThreeVector operator+ (const AliFmThreeVector& v1, const AliFmThreeVector& v2); template<> AliFmThreeVector operator- (const AliFmThreeVector& v1, const AliFmThreeVector& v2); template<> AliFmThreeVector operator- (const AliFmThreeVector& v1, const AliFmThreeVector& v2); template<> AliFmThreeVector operator- (const AliFmThreeVector& v1, const AliFmThreeVector& v2); template<> AliFmThreeVector operator- (const AliFmThreeVector& v1, const AliFmThreeVector& v2); template<> AliFmThreeVector operator* (const AliFmThreeVector& v1, const AliFmThreeVector& v2); template<> AliFmThreeVector operator* (const AliFmThreeVector& v1, const AliFmThreeVector& v2); template<> AliFmThreeVector operator* (const AliFmThreeVector& v1, const AliFmThreeVector& v2); template<> AliFmThreeVector operator* (const AliFmThreeVector& v1, const AliFmThreeVector& v2); template<> AliFmThreeVector operator* (const double v1, const AliFmThreeVector& v2); template<> AliFmThreeVector operator* (const AliFmThreeVector& v1, const double v2); template<> AliFmThreeVector operator* (const double v1, const AliFmThreeVector& v2); template<> AliFmThreeVector operator* (const AliFmThreeVector& v1, const double v2); template<> AliFmThreeVector operator/ (const AliFmThreeVector& v1, const AliFmThreeVector& v2); template<> AliFmThreeVector operator/ (const AliFmThreeVector& v1, const AliFmThreeVector& v2); template<> AliFmThreeVector operator/ (const AliFmThreeVector& v1, const AliFmThreeVector& v2); template<> AliFmThreeVector operator/ (const AliFmThreeVector& v1, const AliFmThreeVector& v2); template<> AliFmThreeVector operator/ (const double v1, const AliFmThreeVector& v2); template<> AliFmThreeVector operator/ (const double v1, const AliFmThreeVector& v2); template<> AliFmThreeVector operator/ (const AliFmThreeVector& v1, const double v2); template<> AliFmThreeVector operator/ (const AliFmThreeVector& v1, const double v2); template<> istream& operator>>(istream& is,const AliFmThreeVector& v); template<> istream& operator>>(istream& is,const AliFmThreeVector& v); template<> ostream& operator<<(ostream& os,const AliFmThreeVector& v); template<> ostream& operator<<(ostream& os,const AliFmThreeVector& v); #else // // Non-member functions // template inline T abs(const AliFmThreeVector& v) {return v.mag();} template inline AliFmThreeVector cross_product(const AliFmThreeVector& v1, const AliFmThreeVector& v2) { return v1.cross(v2); } // // Non-member operators // template inline AliFmThreeVector operator+ (const AliFmThreeVector& v1, const AliFmThreeVector& v2) { return AliFmThreeVector(v1) += v2; } template inline AliFmThreeVector operator- (const AliFmThreeVector& v1, const AliFmThreeVector& v2) { return AliFmThreeVector(v1) -= v2; } template inline T operator* (const AliFmThreeVector& v1, const AliFmThreeVector& v2) { return AliFmThreeVector(v1).dot(v2); } template inline AliFmThreeVector operator* (const AliFmThreeVector& v, double c) { return AliFmThreeVector