PWGCF/FEMTOSCOPY/AliFemto/AliFemtoThreeVector.h

   1 /***************************************************************************
   2  *
   3  * $Id$
   4  *
   5  * Author: Brian Lasiuk, Thomas Ullrich, April 1998
   6  ***************************************************************************
   7  *
   8  * Description:
   9  *
  10  * Remarks:   Since not all compilers support member templates
  11  *            we have to specialize the templated member on these
  12  *            platforms. If member templates are not supported the
  13  *            ST_NO_MEMBER_TEMPLATES flag has to be set. tu.
  14  *
  15  ***************************************************************************
  16  *
  17  * $Log$
  18  * Revision 1.1.1.1  2007/04/25 15:38:41  panos
  19  * Importing the HBT code dir
  20  *
  21  * Revision 1.1.1.1  2007/03/07 10:14:49  mchojnacki
  22  * First version on CVS
  23  *
  24  * Revision 1.15  2005/09/22 20:09:20  fisyak
  25  * Make AliFmLorentzVector persistent
  26  *
  27  * Revision 1.14  2005/07/19 22:27:11  perev
  28  * Cleanup
  29  *
  30  * Revision 1.13  2005/07/06 18:49:57  fisyak
  31  * Replace AliFmHelixD, AliFmLorentzVectorD,AliFmLorentzVectorF,AliFmMatrixD,AliFmMatrixF,AliFmPhysicalHelixD,AliFmThreeVectorD,AliFmThreeVectorF by templated version
  32  *
  33  * Revision 1.12  2005/03/28 06:03:41  perev
  34  * Defence FPE added
  35  *
  36  * Revision 1.11  2004/12/02 20:07:32  fine
  37  * define the valid method for both flavor of AliFmThreeVector
  38  *
  39  * Revision 1.10  2003/10/30 20:06:46  perev
  40  * Check of quality added
  41  *
  42  * Revision 1.9  2003/09/02 17:59:35  perev
  43  * gcc 3.2 updates + WarnOff
  44  *
  45  * Revision 1.8  2002/06/21 17:47:37  genevb
  46  * Added pseudoProduct
  47  *
  48  * Revision 1.7  2000/01/04 19:56:05  ullrich
  49  * Added cpp macro for CINT.
  50  *
  51  * Revision 1.6  1999/12/21 15:14:31  ullrich
  52  * Modified to cope with new compiler version on Sun (CC5.0).
  53  *
  54  * Revision 1.5  1999/10/15 15:46:54  ullrich
  55  * Changed output format in operator<<
  56  *
  57  * Revision 1.4  1999/06/04 18:00:05  ullrich
  58  * Added new constructor which takes C-style array as argument.
  59  * New operators operator() and operator[] which can be used
  60  * as lvalues.
  61  *
  62  * Revision 1.3  1999/02/17 11:42:19  ullrich
  63  * Removed specialization for 'long double'.
  64  *
  65  * Revision 1.2  1999/02/14 23:11:48  fisyak
  66  * Fixes for Rootcint
  67  *
  68  * Revision 1.1  1999/01/30 03:59:05  fisyak
  69  * Root Version of AliFmarClassLibrary
  70  *
  71  * Revision 1.1  1999/01/23 00:28:04  ullrich
  72  * Initial Revision
  73  *
  74  **************************************************************************/
  75 /*//
  76 //// General class for a three-vector
  77 ///*/
  78 #ifndef ST_THREE_VECTOR_HH
  79 #define ST_THREE_VECTOR_HH
  80 #ifdef __ROOT__
  81 #include "Rtypes.h"
  82 #endif
  83 #ifndef __CINT__
  84 #include <iostream>
  85 #include <fstream>
  86 #include <math.h>
  87 #ifdef GNU_GCC
  88 #    include <stddef.h>
  89 #endif
  90 #if defined (__SUNPRO_CC) && __SUNPRO_CC < 0x500
  91 #    include <stdcomp.h>
  92 #endif
  93 #ifndef ST_NO_EXCEPTIONS
  94 #    include <stdexcept>
  95 #    if !defined(ST_NO_NAMESPACES)
  96 using std::out_of_range;
  97 #    endif
  98 #endif
  99 #endif // __CINT__
 100
 101 #ifdef WIN32
 102 #include "gcc2vs.h"
 103 #endif
 104
 105 class TRootIOCtor;//nic nie rozumiem
 106 using namespace std;
 107
 108
 109 template<class T> class AliFmThreeVector {
 110 public:
 111     AliFmThreeVector(T = 0, T = 0, T = 0);
 112   //                     ROOT_VERSION(5,03,01)
 113 #if ROOT_VERSION_CODE >= 328449
 114    AliFmThreeVector(TRootIOCtor*) : mX1(0), mX2(0), mX3(0) {}
 115 #endif
 116     virtual ~AliFmThreeVector();
 117
 118 #ifndef ST_NO_MEMBER_TEMPLATES
 119     template<class X> AliFmThreeVector(const AliFmThreeVector<X>&);
 120     template<class X> AliFmThreeVector(const X*);
 121     template<class X> AliFmThreeVector<T>& operator=(const AliFmThreeVector<X>&);
 122     // AliFmThreeVector(const AliFmThreeVector<T>&);                use default
 123     // AliFmThreeVector<T>& operator=(const AliFmThreeVector<T>&);  use default
 124 #else
 125     AliFmThreeVector(const AliFmThreeVector<float>&);
 126     AliFmThreeVector(const AliFmThreeVector<double>&);
 127
 128     AliFmThreeVector(const float*);
 129     AliFmThreeVector(const double*);
 130
 131     AliFmThreeVector<T>& operator=(const AliFmThreeVector<float>&);
 132     AliFmThreeVector<T>& operator=(const AliFmThreeVector<double>&);
 133 #endif
 134
 135     void SetX(T);
 136     void SetY(T);
 137     void SetZ(T);
 138
 139     void SetPhi(T);
 140     void SetTheta(T);
 141     void SetMag(T);
 142     void SetMagnitude(T);
 143
 144     T   x()                        const;
 145     T   y()                        const;
 146     T   z()                        const;
 147     T   Theta()                    const;
 148     T   CosTheta()                 const;
 149     T   Phi()                      const;
 150     T   Perp()                     const;
 151     T   Perp2()                    const;
 152     T   Magnitude()                const;
 153     T   Mag()                      const;
 154     T   Mag2()                     const;
 155     T   PseudoRapidity()           const;
 156     T   operator() (size_t)        const;
 157     T   operator[] (size_t)        const;
 158
 159     T&  operator() (size_t);
 160     T&  operator[] (size_t);
 161
 162     T   MassHypothesis(T mass)     const;
 163
 164     AliFmThreeVector<T>  Unit()       const;
 165     AliFmThreeVector<T>  Orthogonal() const;
 166
 167     void  RotateX(T);
 168     void  RotateY(T);
 169     void  RotateZ(T);
 170
 171     AliFmThreeVector<T>  operator- ();
 172     AliFmThreeVector<T>  operator+ ();
 173     AliFmThreeVector<T>& operator*= (double);
 174     AliFmThreeVector<T>& operator/= (double);
 175     AliFmThreeVector<T>  PseudoProduct(double,double,double) const;
 176
 177 #ifndef ST_NO_MEMBER_TEMPLATES
 178     template<class X> T                Angle(const AliFmThreeVector<X>&) const;
 179     template<class X> AliFmThreeVector<T> Cross(const AliFmThreeVector<X>&) const;
 180     template<class X> T                Dot  (const AliFmThreeVector<X>&) const;
 181     template<class X> AliFmThreeVector<T> PseudoProduct(const AliFmThreeVector<X>&) const;
 182
 183     template<class X> bool operator == (const AliFmThreeVector<X>& v) const;
 184     template<class X> bool operator != (const AliFmThreeVector<X>& v) const;
 185
 186     template<class X> AliFmThreeVector<T>& operator+= (const AliFmThreeVector<X>&);
 187     template<class X> AliFmThreeVector<T>& operator-= (const AliFmThreeVector<X>&);
 188 #else
 189     T                Angle(const AliFmThreeVector<float>&) const;
 190     AliFmThreeVector<T> Cross(const AliFmThreeVector<float>&) const;
 191     T                Dot  (const AliFmThreeVector<float>&) const;
 192     AliFmThreeVector<T> PseudoProduct(const AliFmThreeVector<float>&) const;
 193
 194     T                Angle(const AliFmThreeVector<double>&) const;
 195     T                Dot  (const AliFmThreeVector<double>&) const;
 196     AliFmThreeVector<T> Cross(const AliFmThreeVector<double>&) const;
 197     AliFmThreeVector<T> PseudoProduct(const AliFmThreeVector<double>&) const;
 198
 199     bool operator == (const AliFmThreeVector<float>& v) const;
 200     bool operator != (const AliFmThreeVector<float>& v) const;
 201     AliFmThreeVector<T>& operator+= (const AliFmThreeVector<float>&);
 202     AliFmThreeVector<T>& operator-= (const AliFmThreeVector<float>&);
 203
 204     bool operator == (const AliFmThreeVector<double>& v) const;
 205     bool operator != (const AliFmThreeVector<double>& v) const;
 206     AliFmThreeVector<T>& operator+= (const AliFmThreeVector<double>&);
 207     AliFmThreeVector<T>& operator-= (const AliFmThreeVector<double>&);
 208 #endif
 209   int             valid(double world = 1.e+5) const;
 210     int               bad(double world = 1.e+5) const;
 211 protected:
 212     T    mX1, mX2, mX3;  // Three vector components
 213 #ifdef __ROOT__
 214   ClassDef(AliFmThreeVector,3)
 215 #endif /* __ROOT__ */
 216 };
 217
 218 #ifndef __CINT__
 219 //
 220 //        Implementation of member functions
 221 //
 222 template<class T>
 223 inline AliFmThreeVector<T>::AliFmThreeVector(T x, T y, T z)
 224     : mX1(x), mX2(y), mX3(z) {/* nop */}
 225 template<class T>
 226 inline AliFmThreeVector<T>::~AliFmThreeVector() {/* nop */}
 227
 228 template<class T>
 229 inline void AliFmThreeVector<T>::SetX(T x) {mX1 = x;}
 230
 231 template<class T>
 232 inline void AliFmThreeVector<T>::SetY(T y) {mX2 = y;}
 233
 234 template<class T>
 235 inline void AliFmThreeVector<T>::SetZ(T z) {mX3 = z;}
 236
 237 template<class T>
 238 void AliFmThreeVector<T>::SetPhi(T angle)
 239 {
 240     double  r = Magnitude();
 241     double th = Theta();
 242
 243     mX1 = r*sin(th)*cos(angle);
 244     mX2 = r*sin(th)*sin(angle);
 245 }
 246
 247 template <class T>
 248 void AliFmThreeVector<T>::SetTheta(T angle)
 249 {
 250     double r  = Magnitude();
 251     double ph = Phi();
 252
 253     mX1 = r*sin(angle)*cos(ph);
 254     mX2 = r*sin(angle)*sin(ph);
 255     mX3 = r*cos(angle);
 256 }
 257
 258 template <class T>
 259 void AliFmThreeVector<T>::SetMagnitude(T r)
 260 {
 261     double th = Theta();
 262     double ph = Phi();
 263
 264     mX1 = r*sin(th)*cos(ph);
 265     mX2 = r*sin(th)*sin(ph);
 266     mX3 = r*cos(th);
 267 }
 268
 269 template <class T>
 270 void AliFmThreeVector<T>::SetMag(T mag)
 271 {
 272     SetMagnitude(mag);
 273 }
 274
 275 template<class T>
 276 inline T AliFmThreeVector<T>::x() const {return mX1;}
 277
 278 template<class T>
 279 inline T AliFmThreeVector<T>::y() const {return mX2;}
 280
 281 template<class T>
 282 inline T AliFmThreeVector<T>::z() const {return mX3;}
 283
 284 template<class T>
 285 inline T AliFmThreeVector<T>::Theta() const
 286 {
 287   return acos(CosTheta());
 288 }
 289
 290 template<class T>
 291 inline T AliFmThreeVector<T>::CosTheta() const
 292 {
 293   return mX3/(Mag()+1e-20);
 294 }
 295
 296 template<class T>
 297 inline T AliFmThreeVector<T>::Phi() const
 298 {
 299     return atan2(mX2,mX1);
 300 }
 301
 302 template<class T>
 303 inline T AliFmThreeVector<T>::PseudoRapidity() const
 304 {
 305     //
 306     // change code to more optimal:
 307     // double m = Mag();
 308     // return 0.5*::log( (m+z())/(m-z()) );
 309     double tmp = tan(Theta()/2.); if (tmp <=0.) return 1e20;
 310     return -::log(tmp);
 311 }
 312
 313 template<class T>
 314 inline AliFmThreeVector<T> AliFmThreeVector<T>::Unit() const
 315 {
 316     double tmp = Mag(); if (tmp<=0.) tmp = 1e-20;
 317     return *this/tmp;
 318 }
 319
 320 template <class T>
 321 T AliFmThreeVector<T>::MassHypothesis(T mass) const
 322 {
 323     return ::sqrt((*this)*(*this) + mass*mass);
 324 }
 325
 326 template <class T>
 327 AliFmThreeVector<T> AliFmThreeVector<T>::Orthogonal() const
 328 {
 329     // Direct copy from CLHEP--it is probably better to
 330     // use your own dot/cross product code...
 331     double x = (mX1 < 0.0) ? -mX1 : mX1;
 332     double y = (mX2 < 0.0) ? -mX2 : mX2;
 333     double z = (mX3 < 0.0) ? -mX3 : mX3;
 334
 335     if(x<y)
 336         return x < z ? AliFmThreeVector<T>(0,mX3,-mX2) :  AliFmThreeVector<T>(mX2,-mX1,0);
 337     else
 338         return  mX2 < mX3 ? AliFmThreeVector<T>(-mX3,0,mX1) :  AliFmThreeVector<T>(mX2,-mX1,0);
 339 }
 340
 341 template <class T>
 342 void AliFmThreeVector<T>::RotateX(T angle)
 343 {
 344     // may in the future make use of the AliFmRotation class!
 345     double yPrime = cos(angle)*mX2 - sin(angle)*mX3;
 346     double zPrime = sin(angle)*mX2 + cos(angle)*mX3;
 347
 348     mX2 = yPrime;
 349     mX3 = zPrime;
 350 }
 351
 352 template <class T>
 353 void AliFmThreeVector<T>::RotateY(T angle)
 354 {
 355     // may in the future make use of the AliFmRotation class!
 356     double zPrime = cos(angle)*mX3 - sin(angle)*mX1;
 357     double xPrime = sin(angle)*mX3 + cos(angle)*mX1;
 358
 359     mX1 = xPrime;
 360     mX3 = zPrime;
 361 }
 362
 363 template <class T>
 364 void AliFmThreeVector<T>::RotateZ(T angle)
 365 {
 366     // may in the future make use of the AliFmRotation class!
 367     double xPrime = cos(angle)*mX1 - sin(angle)*mX2;
 368     double yPrime = sin(angle)*mX1 + cos(angle)*mX2;
 369
 370     mX1 = xPrime;
 371     mX2 = yPrime;
 372 }
 373
 374 template<class T>
 375 inline T AliFmThreeVector<T>::Perp() const
 376 {
 377     return ::sqrt(mX1*mX1+mX2*mX2);
 378 }
 379
 380 template<class T>
 381 inline T AliFmThreeVector<T>::Perp2() const
 382 {
 383     return mX1*mX1+mX2*mX2;
 384 }
 385
 386 template<class T>
 387 inline T AliFmThreeVector<T>::Magnitude() const
 388 {
 389     return Mag();
 390 }
 391
 392 template<class T>
 393 inline T AliFmThreeVector<T>::Mag() const
 394 {
 395     return ::sqrt(mX1*mX1+mX2*mX2+mX3*mX3);
 396 }
 397
 398 template<class T>
 399 inline T AliFmThreeVector<T>::Mag2() const
 400 {
 401     return mX1*mX1+mX2*mX2+mX3*mX3;
 402 }
 403
 404 template<class T>
 405 inline T AliFmThreeVector<T>::operator() (size_t i) const
 406 {
 407     if (0 <=i && i <= 2)  return (&mX1)[i];
 408 #ifndef ST_NO_EXCEPTIONS
 409     throw out_of_range("AliFmThreeVector<T>::operator(): bad index");
 410 #else
 411     cerr << "AliFmThreeVector<T>::operator(): bad index" << endl;
 412 #endif
 413     return 0;
 414 }
 415
 416 template<class T>
 417 inline T& AliFmThreeVector<T>::operator() (size_t i)
 418 {
 419     if (0 <=i && i <= 2)  return (&mX1)[i];
 420 #ifndef ST_NO_EXCEPTIONS
 421     throw out_of_range("AliFmThreeVector<T>::operator(): bad index");
 422 #else
 423     cerr << "AliFmThreeVector<T>::operator(): bad index" << endl;
 424 #endif
 425     return mX1;
 426 }
 427
 428 template<class T>
 429 inline T AliFmThreeVector<T>::operator[] (size_t i) const
 430 {
 431     if (0 <=i && i <= 2)  return (&mX1)[i];
 432 #ifndef ST_NO_EXCEPTIONS
 433       throw out_of_range("AliFmThreeVector<T>::operator[]: bad index");
 434 #else
 435       cerr << "AliFmThreeVector<T>::operator[]: bad index" << endl;
 436 #endif
 437       return 0;
 438 }
 439
 440 template<class T>
 441 inline T &AliFmThreeVector<T>::operator[] (size_t i)
 442 {
 443     if (0 <=i && i <= 2)  return (&mX1)[i];
 444 #ifndef ST_NO_EXCEPTIONS
 445       throw out_of_range("AliFmThreeVector<T>::operator[]: bad index");
 446 #else
 447       cerr << "AliFmThreeVector<T>::operator[]: bad index" << endl;
 448 #endif
 449       return mX1;
 450 }
 451
 452 template<class T>
 453 inline AliFmThreeVector<T>& AliFmThreeVector<T>::operator*= (double c)
 454 {
 455     mX1 *= c; mX2 *= c; mX3 *= c;
 456     return *this;
 457 }
 458
 459 template<class T>
 460 inline AliFmThreeVector<T>& AliFmThreeVector<T>::operator/= (double c)
 461 {
 462     mX1 /= c; mX2 /= c; mX3 /= c;
 463     return *this;
 464 }
 465
 466 template<class T>
 467 inline AliFmThreeVector<T>
 468 AliFmThreeVector<T>::PseudoProduct(double x,double y,double z) const
 469 {
 470     return AliFmThreeVector<T>(mX1*x,mX2*y,mX3*z);
 471 }
 472
 473 template<class T>
 474 AliFmThreeVector<T> AliFmThreeVector<T>::operator- ()
 475 {
 476     return AliFmThreeVector<T>(-mX1, -mX2, -mX3);
 477 }
 478
 479 template<class T>
 480 AliFmThreeVector<T> AliFmThreeVector<T>::operator+ ()
 481 {
 482     return *this;
 483 }
 484
 485 #ifndef ST_NO_MEMBER_TEMPLATES
 486 #ifndef WIN32
 487
 488 template<class T>
 489 template<class X>
 490 inline AliFmThreeVector<T>::AliFmThreeVector(const AliFmThreeVector<X>& v)
 491     : mX1(v.x()), mX2(v.y()), mX3(v.z()) {/* nop */}
 492
 493 template<class T>
 494 template<class X>
 495 inline AliFmThreeVector<T>::AliFmThreeVector(const X *a)
 496 {
 497     mX1 = a[0];
 498     mX2 = a[1];
 499     mX3 = a[2];
 500 }
 501
 502 template<class T>
 503 template<class X>
 504 inline AliFmThreeVector<T>&
 505 AliFmThreeVector<T>::operator=(const AliFmThreeVector<X>& v)
 506 {
 507   if(this != &v) {
 508     mX1 = v.x();  mX2 = v.y();  mX3 = v.z();
 509   }
 510   return *this;
 511 }
 512
 513 template<class T>
 514 template<class X>
 515 inline bool AliFmThreeVector<T>::operator== (const AliFmThreeVector<X>& v) const
 516 {
 517     return mX1 == v.x() && mX2 == v.y() && mX3 == v.z();
 518 }
 519
 520 template<class T>
 521 template<class X>
 522 inline bool AliFmThreeVector<T>::operator!= (const AliFmThreeVector<X>& v) const
 523 {
 524     return !(*this == v);
 525 }
 526
 527 template<class T>
 528 template<class X>
 529 inline AliFmThreeVector<T>&
 530 AliFmThreeVector<T>::operator+= (const AliFmThreeVector<X>& v)
 531 {
 532     mX1 += v.x(); mX2 += v.y(); mX3 += v.z();
 533     return *this;
 534 }
 535
 536 template<class T>
 537 template<class X>
 538 inline AliFmThreeVector<T>&
 539 AliFmThreeVector<T>::operator-= (const AliFmThreeVector<X>& v)
 540 {
 541     mX1 -= v.x(); mX2 -= v.y(); mX3 -= v.z();
 542     return *this;
 543 }
 544
 545 template<class T>
 546 template<class X>
 547 inline T AliFmThreeVector<T>::Dot(const AliFmThreeVector<X>& v) const
 548 {
 549     return mX1*v.x() + mX2*v.y() + mX3*v.z();
 550 }
 551
 552 template<class T>
 553 template<class X>
 554 inline AliFmThreeVector<T>
 555 AliFmThreeVector<T>::Cross(const AliFmThreeVector<X>& v) const
 556 {
 557     return AliFmThreeVector<T>(mX2*v.z() - mX3*v.y(),
 558                             mX3*v.x() - mX1*v.z(),
 559                             mX1*v.y() - mX2*v.x());
 560 }
 561
 562 template<class T>
 563 template<class X>
 564 inline T AliFmThreeVector<T>::Angle(const AliFmThreeVector<X>& vec) const
 565 {
 566     double norm = this->Mag2()*vec.Mag2();
 567
 568     return norm > 0 ? acos(this->Dot(vec)/(::sqrt(norm))) : 0;
 569 }
 570
 571 template<class T>
 572 template<class X>
 573 inline AliFmThreeVector<T>
 574 AliFmThreeVector<T>::PseudoProduct(const AliFmThreeVector<X>& v) const
 575 {
 576     return this->PseudoProduct(v.x(),v.y(),v.z());
 577 }
 578
 579 #endif
 580 #else
 581
 582 template<class T>
 583 inline AliFmThreeVector<T>::AliFmThreeVector(const AliFmThreeVector<float>& v)
 584     : mX1(v.x()), mX2(v.y()), mX3(v.z()) {/* nop */}
 585
 586 template<class T>
 587 inline AliFmThreeVector<T>::AliFmThreeVector(const AliFmThreeVector<double>& v)
 588     : mX1(v.x()), mX2(v.y()), mX3(v.z()) {/* nop */}
 589
 590 template<class T>
 591 inline AliFmThreeVector<T>::AliFmThreeVector(const float *a)
 592 {
 593     mX1 = a[0];
 594     mX2 = a[1];
 595     mX3 = a[2];
 596 }
 597
 598 template<class T>
 599 inline AliFmThreeVector<T>::AliFmThreeVector(const double *a)
 600 {
 601     mX1 = a[0];
 602     mX2 = a[1];
 603     mX3 = a[2];
 604 }
 605
 606 template<class T>
 607 inline AliFmThreeVector<T>&
 608 AliFmThreeVector<T>::operator=(const AliFmThreeVector<float>& v)
 609 {
 610     mX1 = v.x();  mX2 = v.y();  mX3 = v.z();
 611     return *this;
 612 }
 613
 614 template<class T>
 615 inline AliFmThreeVector<T>&
 616 AliFmThreeVector<T>::operator=(const AliFmThreeVector<double>& v)
 617 {
 618     mX1 = v.x();  mX2 = v.y();  mX3 = v.z();
 619     return *this;
 620 }
 621
 622 template<class T>
 623 inline bool
 624 AliFmThreeVector<T>::operator== (const AliFmThreeVector<float>& v) const
 625 {
 626     return mX1 == v.x() && mX2 == v.y() && mX3 == v.z();
 627 }
 628
 629 template<class T>
 630 inline bool
 631 AliFmThreeVector<T>::operator== (const AliFmThreeVector<double>& v) const
 632 {
 633     return mX1 == v.x() && mX2 == v.y() && mX3 == v.z();
 634 }
 635
 636 template<class T>
 637 inline bool
 638 AliFmThreeVector<T>::operator!= (const AliFmThreeVector<float>& v) const
 639 {
 640     return !(*this == v);
 641 }
 642
 643 template<class T>
 644 inline bool
 645 AliFmThreeVector<T>::operator!= (const AliFmThreeVector<double>& v) const
 646 {
 647     return !(*this == v);
 648 }
 649
 650 template<class T>
 651 inline AliFmThreeVector<T>&
 652 AliFmThreeVector<T>::operator+= (const AliFmThreeVector<float>& v)
 653 {
 654     mX1 += v.x(); mX2 += v.y(); mX3 += v.z();
 655     return *this;
 656 }
 657
 658 template<class T>
 659 inline AliFmThreeVector<T>&
 660 AliFmThreeVector<T>::operator+= (const AliFmThreeVector<double>& v)
 661 {
 662     mX1 += v.x(); mX2 += v.y(); mX3 += v.z();
 663     return *this;
 664 }
 665
 666 template<class T>
 667 inline AliFmThreeVector<T>&
 668 AliFmThreeVector<T>::operator-= (const AliFmThreeVector<float>& v)
 669 {
 670     mX1 -= v.x(); mX2 -= v.y(); mX3 -= v.z();
 671     return *this;
 672 }
 673
 674 template<class T>
 675 inline AliFmThreeVector<T>&
 676 AliFmThreeVector<T>::operator-= (const AliFmThreeVector<double>& v)
 677 {
 678     mX1 -= v.x(); mX2 -= v.y(); mX3 -= v.z();
 679     return *this;
 680 }
 681
 682 template<class T>
 683 inline T AliFmThreeVector<T>::Dot(const AliFmThreeVector<float>& v) const
 684 {
 685     return mX1*v.x() + mX2*v.y() + mX3*v.z();
 686 }
 687
 688 template<class T>
 689 inline T AliFmThreeVector<T>::Dot(const AliFmThreeVector<double>& v) const
 690 {
 691     return mX1*v.x() + mX2*v.y() + mX3*v.z();
 692 }
 693
 694 template<class T>
 695 inline AliFmThreeVector<T>
 696 AliFmThreeVector<T>::Cross(const AliFmThreeVector<float>& v) const
 697 {
 698     return AliFmThreeVector<T>(mX2*v.z() - mX3*v.y(),
 699                             mX3*v.x() - mX1*v.z(),
 700                             mX1*v.y() - mX2*v.x());
 701 }
 702
 703 template<class T>
 704 inline AliFmThreeVector<T>
 705 AliFmThreeVector<T>::Cross(const AliFmThreeVector<double>& v) const
 706 {
 707     return AliFmThreeVector<T>(mX2*v.z() - mX3*v.y(),
 708                             mX3*v.x() - mX1*v.z(),
 709                             mX1*v.y() - mX2*v.x());
 710 }
 711
 712 template<class T>
 713 inline T AliFmThreeVector<T>::Angle(const AliFmThreeVector<float>& v) const
 714 {
 715     double tmp = Mag()*v.Mag(); if (tmp <=0) tmp = 1e-20;
 716     return acos(this->Dot(v)/tmp);
 717 }
 718
 719 template<class T>
 720 inline T AliFmThreeVector<T>::Angle(const AliFmThreeVector<double>& v) const
 721 {
 722     double tmp = Mag()*v.Mag(); if (tmp <=0) tmp = 1e-20;
 723     return acos(this->Dot(v)/tmp);
 724 }
 725
 726 template<class T>
 727 inline AliFmThreeVector<T>
 728 AliFmThreeVector<T>::PseudoProduct(const AliFmThreeVector<float>& v) const
 729 {
 730     return this->PseudoProduct(v.x(),v.y(),v.z());
 731 }
 732
 733 template<class T>
 734 inline AliFmThreeVector<T>
 735 AliFmThreeVector<T>::PseudoProduct(const AliFmThreeVector<double>& v) const
 736 {
 737     return this->PseudoProduct(v.x(),v.y(),v.z());
 738 }
 739 #endif  // ST_NO_MEMBER_TEMPLATES
 740 template<class T>
 741 inline int
 742 AliFmThreeVector<T>::valid(double world) const  {return !bad(world);}
 743
 744 template<class T>
 745 inline int
 746 AliFmThreeVector<T>::bad(double world) const
 747 {
 748   for (int i=0;i<3;i++) {
 749           if (!finite((&mX1)[i])      ) return 10+i;
 750           if ( fabs  ((&mX1)[i])>world) return 20+i;
 751   }
 752   return 0;
 753 }
 754 #endif /*! __CINT__ */
 755 #ifdef __CINT__
 756 template<> float abs(const AliFmThreeVector<float>& v);
 757 template<> double abs(const AliFmThreeVector<double>& v);
 758 template<> AliFmThreeVector<double> cross_product(const AliFmThreeVector<double>& v1, const AliFmThreeVector<double>& v2);
 759 template<> AliFmThreeVector<float>  cross_product(const AliFmThreeVector<float>&  v1, const AliFmThreeVector<float>& v2);
 760 template<> AliFmThreeVector<double> cross_product(const AliFmThreeVector<float>&  v1, const AliFmThreeVector<double>& v2);
 761 template<> AliFmThreeVector<double> cross_product(const AliFmThreeVector<double>& v1, const AliFmThreeVector<float>& v2);
 762 template<> AliFmThreeVector<double> operator+ (const AliFmThreeVector<double>& v1, const AliFmThreeVector<double>& v2);
 763 template<> AliFmThreeVector<float>  operator+ (const AliFmThreeVector<float>&  v1, const AliFmThreeVector<float>& v2);
 764 template<> AliFmThreeVector<double> operator+ (const AliFmThreeVector<double>& v1, const AliFmThreeVector<float>& v2);
 765 template<> AliFmThreeVector<double> operator+ (const AliFmThreeVector<float>&  v1, const AliFmThreeVector<double>& v2);
 766 template<> AliFmThreeVector<double> operator- (const AliFmThreeVector<double>& v1, const AliFmThreeVector<double>& v2);
 767 template<> AliFmThreeVector<float>  operator- (const AliFmThreeVector<float>&  v1, const AliFmThreeVector<float>& v2);
 768 template<> AliFmThreeVector<double> operator- (const AliFmThreeVector<double>& v1, const AliFmThreeVector<float>& v2);
 769 template<> AliFmThreeVector<double> operator- (const AliFmThreeVector<float>&  v1, const AliFmThreeVector<double>& v2);
 770 template<> AliFmThreeVector<double> operator* (const AliFmThreeVector<double>& v1, const AliFmThreeVector<double>& v2);
 771 template<> AliFmThreeVector<float>  operator* (const AliFmThreeVector<float>&  v1, const AliFmThreeVector<float>& v2);
 772 template<> AliFmThreeVector<double> operator* (const AliFmThreeVector<double>& v1, const AliFmThreeVector<float>& v2);
 773 template<> AliFmThreeVector<double> operator* (const AliFmThreeVector<float>&  v1, const AliFmThreeVector<double>& v2);
 774 template<> AliFmThreeVector<double> operator* (const double                 v1, const AliFmThreeVector<float>& v2);
 775 template<> AliFmThreeVector<double> operator* (const AliFmThreeVector<float>&  v1, const double v2);
 776 template<> AliFmThreeVector<double> operator* (const double                 v1, const AliFmThreeVector<double>& v2);
 777 template<> AliFmThreeVector<double> operator* (const AliFmThreeVector<double>& v1, const double v2);
 778 template<> AliFmThreeVector<double> operator/ (const AliFmThreeVector<double>& v1, const AliFmThreeVector<double>& v2);
 779 template<> AliFmThreeVector<float>  operator/ (const AliFmThreeVector<float>&  v1, const AliFmThreeVector<float>& v2);
 780 template<> AliFmThreeVector<double> operator/ (const AliFmThreeVector<double>& v1, const AliFmThreeVector<float>& v2);
 781 template<> AliFmThreeVector<double> operator/ (const AliFmThreeVector<float>&  v1, const AliFmThreeVector<double>& v2);
 782 template<> AliFmThreeVector<double> operator/ (const                double  v1, const AliFmThreeVector<double>& v2);
 783 template<> AliFmThreeVector<float>  operator/ (const                double  v1, const AliFmThreeVector<float>& v2);
 784 template<> AliFmThreeVector<double> operator/ (const AliFmThreeVector<double>& v1, const double v2);
 785 template<> AliFmThreeVector<double> operator/ (const AliFmThreeVector<float>&  v1, const double v2);
 786 template<> istream&  operator>>(istream& is,const AliFmThreeVector<double>& v);
 787 template<> istream&  operator>>(istream& is,const AliFmThreeVector<float>& v);
 788 template<> ostream&  operator<<(ostream& os,const AliFmThreeVector<double>& v);
 789 template<> ostream&  operator<<(ostream& os,const AliFmThreeVector<float>& v);
 790 #else
 791 //
 792 //        Non-member functions
 793 //
 794 template<class T>
 795 inline T abs(const AliFmThreeVector<T>& v) {return v.Mag();}
 796
 797 template<class T, class X>
 798 inline AliFmThreeVector<T>
 799 cross_product(const AliFmThreeVector<T>& v1, const AliFmThreeVector<X>& v2)
 800 {
 801     return v1.Cross(v2);
 802 }
 803
 804
 805 //
 806 //        Non-member operators
 807 //
 808 template<class T, class X>
 809 inline AliFmThreeVector<T>
 810 operator+ (const AliFmThreeVector<T>& v1, const AliFmThreeVector<X>& v2)
 811 {
 812     return AliFmThreeVector<T>(v1) += v2;
 813 }
 814
 815 template<class T, class X>
 816 inline AliFmThreeVector<T>
 817 operator- (const AliFmThreeVector<T>& v1, const AliFmThreeVector<X>& v2)
 818 {
 819     return AliFmThreeVector<T>(v1) -= v2;
 820 }
 821
 822 template<class T, class X>
 823 inline T operator* (const AliFmThreeVector<T>& v1, const AliFmThreeVector<X>& v2)
 824 {
 825     return AliFmThreeVector<T>(v1).Dot(v2);
 826 }
 827
 828 template<class T>
 829 inline AliFmThreeVector<T> operator* (const AliFmThreeVector<T>& v, double c)
 830 {
 831     return AliFmThreeVector<T>(v) *= c;
 832 }
 833
 834 template<class T>
 835 inline AliFmThreeVector<T> operator* (double c, const AliFmThreeVector<T>& v)
 836 {
 837     return AliFmThreeVector<T>(v) *= c;
 838 }
 839
 840 template<class T, class X>
 841 inline AliFmThreeVector<T> operator/ (const AliFmThreeVector<T>& v, X c)
 842 {
 843     return AliFmThreeVector<T>(v) /= c;
 844 }
 845
 846 template<class T>
 847 ostream&  operator<<(ostream& os, const AliFmThreeVector<T>& v)
 848 {
 849     return os << v.x() << '\t' << v.y() << '\t' << v.z();
 850 }
 851
 852 template<class T>
 853 istream&  operator>>(istream& is, AliFmThreeVector<T>& v)
 854 {
 855     T  x, y, z;
 856     is >> x >> y >> z;
 857     v.SetX(x);
 858     v.SetY(y);
 859     v.SetZ(z);
 860     return is;
 861 }
 862 #endif /* ! __CINT__ */
 863 #endif