PWG2/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     mX1 = v.x();  mX2 = v.y();  mX3 = v.z();
 508     return *this;
 509 }
 510
 511 template<class T>
 512 template<class X>
 513 inline bool AliFmThreeVector<T>::operator== (const AliFmThreeVector<X>& v) const
 514 {
 515     return mX1 == v.x() && mX2 == v.y() && mX3 == v.z();
 516 }
 517
 518 template<class T>
 519 template<class X>
 520 inline bool AliFmThreeVector<T>::operator!= (const AliFmThreeVector<X>& v) const
 521 {
 522     return !(*this == v);
 523 }
 524
 525 template<class T>
 526 template<class X>
 527 inline AliFmThreeVector<T>&
 528 AliFmThreeVector<T>::operator+= (const AliFmThreeVector<X>& v)
 529 {
 530     mX1 += v.x(); mX2 += v.y(); mX3 += v.z();
 531     return *this;
 532 }
 533
 534 template<class T>
 535 template<class X>
 536 inline AliFmThreeVector<T>&
 537 AliFmThreeVector<T>::operator-= (const AliFmThreeVector<X>& v)
 538 {
 539     mX1 -= v.x(); mX2 -= v.y(); mX3 -= v.z();
 540     return *this;
 541 }
 542
 543 template<class T>
 544 template<class X>
 545 inline T AliFmThreeVector<T>::Dot(const AliFmThreeVector<X>& v) const
 546 {
 547     return mX1*v.x() + mX2*v.y() + mX3*v.z();
 548 }
 549
 550 template<class T>
 551 template<class X>
 552 inline AliFmThreeVector<T>
 553 AliFmThreeVector<T>::Cross(const AliFmThreeVector<X>& v) const
 554 {
 555     return AliFmThreeVector<T>(mX2*v.z() - mX3*v.y(),
 556                             mX3*v.x() - mX1*v.z(),
 557                             mX1*v.y() - mX2*v.x());
 558 }
 559
 560 template<class T>
 561 template<class X>
 562 inline T AliFmThreeVector<T>::Angle(const AliFmThreeVector<X>& vec) const
 563 {
 564     double norm = this->Mag2()*vec.Mag2();
 565
 566     return norm > 0 ? acos(this->Dot(vec)/(::sqrt(norm))) : 0;
 567 }
 568
 569 template<class T>
 570 template<class X>
 571 inline AliFmThreeVector<T>
 572 AliFmThreeVector<T>::PseudoProduct(const AliFmThreeVector<X>& v) const
 573 {
 574     return this->PseudoProduct(v.x(),v.y(),v.z());
 575 }
 576
 577 #endif
 578 #else
 579
 580 template<class T>
 581 inline AliFmThreeVector<T>::AliFmThreeVector(const AliFmThreeVector<float>& v)
 582     : mX1(v.x()), mX2(v.y()), mX3(v.z()) {/* nop */}
 583
 584 template<class T>
 585 inline AliFmThreeVector<T>::AliFmThreeVector(const AliFmThreeVector<double>& v)
 586     : mX1(v.x()), mX2(v.y()), mX3(v.z()) {/* nop */}
 587
 588 template<class T>
 589 inline AliFmThreeVector<T>::AliFmThreeVector(const float *a)
 590 {
 591     mX1 = a[0];
 592     mX2 = a[1];
 593     mX3 = a[2];
 594 }
 595
 596 template<class T>
 597 inline AliFmThreeVector<T>::AliFmThreeVector(const double *a)
 598 {
 599     mX1 = a[0];
 600     mX2 = a[1];
 601     mX3 = a[2];
 602 }
 603
 604 template<class T>
 605 inline AliFmThreeVector<T>&
 606 AliFmThreeVector<T>::operator=(const AliFmThreeVector<float>& v)
 607 {
 608     mX1 = v.x();  mX2 = v.y();  mX3 = v.z();
 609     return *this;
 610 }
 611
 612 template<class T>
 613 inline AliFmThreeVector<T>&
 614 AliFmThreeVector<T>::operator=(const AliFmThreeVector<double>& v)
 615 {
 616     mX1 = v.x();  mX2 = v.y();  mX3 = v.z();
 617     return *this;
 618 }
 619
 620 template<class T>
 621 inline bool
 622 AliFmThreeVector<T>::operator== (const AliFmThreeVector<float>& v) const
 623 {
 624     return mX1 == v.x() && mX2 == v.y() && mX3 == v.z();
 625 }
 626
 627 template<class T>
 628 inline bool
 629 AliFmThreeVector<T>::operator== (const AliFmThreeVector<double>& v) const
 630 {
 631     return mX1 == v.x() && mX2 == v.y() && mX3 == v.z();
 632 }
 633
 634 template<class T>
 635 inline bool
 636 AliFmThreeVector<T>::operator!= (const AliFmThreeVector<float>& v) const
 637 {
 638     return !(*this == v);
 639 }
 640
 641 template<class T>
 642 inline bool
 643 AliFmThreeVector<T>::operator!= (const AliFmThreeVector<double>& v) const
 644 {
 645     return !(*this == v);
 646 }
 647
 648 template<class T>
 649 inline AliFmThreeVector<T>&
 650 AliFmThreeVector<T>::operator+= (const AliFmThreeVector<float>& v)
 651 {
 652     mX1 += v.x(); mX2 += v.y(); mX3 += v.z();
 653     return *this;
 654 }
 655
 656 template<class T>
 657 inline AliFmThreeVector<T>&
 658 AliFmThreeVector<T>::operator+= (const AliFmThreeVector<double>& v)
 659 {
 660     mX1 += v.x(); mX2 += v.y(); mX3 += v.z();
 661     return *this;
 662 }
 663
 664 template<class T>
 665 inline AliFmThreeVector<T>&
 666 AliFmThreeVector<T>::operator-= (const AliFmThreeVector<float>& v)
 667 {
 668     mX1 -= v.x(); mX2 -= v.y(); mX3 -= v.z();
 669     return *this;
 670 }
 671
 672 template<class T>
 673 inline AliFmThreeVector<T>&
 674 AliFmThreeVector<T>::operator-= (const AliFmThreeVector<double>& v)
 675 {
 676     mX1 -= v.x(); mX2 -= v.y(); mX3 -= v.z();
 677     return *this;
 678 }
 679
 680 template<class T>
 681 inline T AliFmThreeVector<T>::Dot(const AliFmThreeVector<float>& v) const
 682 {
 683     return mX1*v.x() + mX2*v.y() + mX3*v.z();
 684 }
 685
 686 template<class T>
 687 inline T AliFmThreeVector<T>::Dot(const AliFmThreeVector<double>& v) const
 688 {
 689     return mX1*v.x() + mX2*v.y() + mX3*v.z();
 690 }
 691
 692 template<class T>
 693 inline AliFmThreeVector<T>
 694 AliFmThreeVector<T>::Cross(const AliFmThreeVector<float>& v) const
 695 {
 696     return AliFmThreeVector<T>(mX2*v.z() - mX3*v.y(),
 697                             mX3*v.x() - mX1*v.z(),
 698                             mX1*v.y() - mX2*v.x());
 699 }
 700
 701 template<class T>
 702 inline AliFmThreeVector<T>
 703 AliFmThreeVector<T>::Cross(const AliFmThreeVector<double>& v) const
 704 {
 705     return AliFmThreeVector<T>(mX2*v.z() - mX3*v.y(),
 706                             mX3*v.x() - mX1*v.z(),
 707                             mX1*v.y() - mX2*v.x());
 708 }
 709
 710 template<class T>
 711 inline T AliFmThreeVector<T>::Angle(const AliFmThreeVector<float>& v) const
 712 {
 713     double tmp = Mag()*v.Mag(); if (tmp <=0) tmp = 1e-20;
 714     return acos(this->Dot(v)/tmp);
 715 }
 716
 717 template<class T>
 718 inline T AliFmThreeVector<T>::Angle(const AliFmThreeVector<double>& v) const
 719 {
 720     double tmp = Mag()*v.Mag(); if (tmp <=0) tmp = 1e-20;
 721     return acos(this->Dot(v)/tmp);
 722 }
 723
 724 template<class T>
 725 inline AliFmThreeVector<T>
 726 AliFmThreeVector<T>::PseudoProduct(const AliFmThreeVector<float>& v) const
 727 {
 728     return this->PseudoProduct(v.x(),v.y(),v.z());
 729 }
 730
 731 template<class T>
 732 inline AliFmThreeVector<T>
 733 AliFmThreeVector<T>::PseudoProduct(const AliFmThreeVector<double>& v) const
 734 {
 735     return this->PseudoProduct(v.x(),v.y(),v.z());
 736 }
 737 #endif  // ST_NO_MEMBER_TEMPLATES
 738 template<class T>
 739 inline int
 740 AliFmThreeVector<T>::valid(double world) const  {return !bad(world);}
 741
 742 template<class T>
 743 inline int
 744 AliFmThreeVector<T>::bad(double world) const
 745 {
 746   for (int i=0;i<3;i++) {
 747           if (!finite((&mX1)[i])      ) return 10+i;
 748           if ( fabs  ((&mX1)[i])>world) return 20+i;
 749   }
 750   return 0;
 751 }
 752 #endif /*! __CINT__ */
 753 #ifdef __CINT__
 754 template<> float abs(const AliFmThreeVector<float>& v);
 755 template<> double abs(const AliFmThreeVector<double>& v);
 756 template<> AliFmThreeVector<double> cross_product(const AliFmThreeVector<double>& v1, const AliFmThreeVector<double>& v2);
 757 template<> AliFmThreeVector<float>  cross_product(const AliFmThreeVector<float>&  v1, const AliFmThreeVector<float>& v2);
 758 template<> AliFmThreeVector<double> cross_product(const AliFmThreeVector<float>&  v1, const AliFmThreeVector<double>& v2);
 759 template<> AliFmThreeVector<double> cross_product(const AliFmThreeVector<double>& v1, const AliFmThreeVector<float>& v2);
 760 template<> AliFmThreeVector<double> operator+ (const AliFmThreeVector<double>& v1, const AliFmThreeVector<double>& v2);
 761 template<> AliFmThreeVector<float>  operator+ (const AliFmThreeVector<float>&  v1, const AliFmThreeVector<float>& v2);
 762 template<> AliFmThreeVector<double> operator+ (const AliFmThreeVector<double>& v1, const AliFmThreeVector<float>& v2);
 763 template<> AliFmThreeVector<double> operator+ (const AliFmThreeVector<float>&  v1, const AliFmThreeVector<double>& v2);
 764 template<> AliFmThreeVector<double> operator- (const AliFmThreeVector<double>& v1, const AliFmThreeVector<double>& v2);
 765 template<> AliFmThreeVector<float>  operator- (const AliFmThreeVector<float>&  v1, const AliFmThreeVector<float>& v2);
 766 template<> AliFmThreeVector<double> operator- (const AliFmThreeVector<double>& v1, const AliFmThreeVector<float>& v2);
 767 template<> AliFmThreeVector<double> operator- (const AliFmThreeVector<float>&  v1, const AliFmThreeVector<double>& v2);
 768 template<> AliFmThreeVector<double> operator* (const AliFmThreeVector<double>& v1, const AliFmThreeVector<double>& v2);
 769 template<> AliFmThreeVector<float>  operator* (const AliFmThreeVector<float>&  v1, const AliFmThreeVector<float>& v2);
 770 template<> AliFmThreeVector<double> operator* (const AliFmThreeVector<double>& v1, const AliFmThreeVector<float>& v2);
 771 template<> AliFmThreeVector<double> operator* (const AliFmThreeVector<float>&  v1, const AliFmThreeVector<double>& v2);
 772 template<> AliFmThreeVector<double> operator* (const double                 v1, const AliFmThreeVector<float>& v2);
 773 template<> AliFmThreeVector<double> operator* (const AliFmThreeVector<float>&  v1, const double v2);
 774 template<> AliFmThreeVector<double> operator* (const double                 v1, const AliFmThreeVector<double>& v2);
 775 template<> AliFmThreeVector<double> operator* (const AliFmThreeVector<double>& v1, const double v2);
 776 template<> AliFmThreeVector<double> operator/ (const AliFmThreeVector<double>& v1, const AliFmThreeVector<double>& v2);
 777 template<> AliFmThreeVector<float>  operator/ (const AliFmThreeVector<float>&  v1, const AliFmThreeVector<float>& v2);
 778 template<> AliFmThreeVector<double> operator/ (const AliFmThreeVector<double>& v1, const AliFmThreeVector<float>& v2);
 779 template<> AliFmThreeVector<double> operator/ (const AliFmThreeVector<float>&  v1, const AliFmThreeVector<double>& v2);
 780 template<> AliFmThreeVector<double> operator/ (const                double  v1, const AliFmThreeVector<double>& v2);
 781 template<> AliFmThreeVector<float>  operator/ (const                double  v1, const AliFmThreeVector<float>& v2);
 782 template<> AliFmThreeVector<double> operator/ (const AliFmThreeVector<double>& v1, const double v2);
 783 template<> AliFmThreeVector<double> operator/ (const AliFmThreeVector<float>&  v1, const double v2);
 784 template<> istream&  operator>>(istream& is,const AliFmThreeVector<double>& v);
 785 template<> istream&  operator>>(istream& is,const AliFmThreeVector<float>& v);
 786 template<> ostream&  operator<<(ostream& os,const AliFmThreeVector<double>& v);
 787 template<> ostream&  operator<<(ostream& os,const AliFmThreeVector<float>& v);
 788 #else
 789 //
 790 //        Non-member functions
 791 //
 792 template<class T>
 793 inline T abs(const AliFmThreeVector<T>& v) {return v.Mag();}
 794
 795 template<class T, class X>
 796 inline AliFmThreeVector<T>
 797 cross_product(const AliFmThreeVector<T>& v1, const AliFmThreeVector<X>& v2)
 798 {
 799     return v1.Cross(v2);
 800 }
 801
 802
 803 //
 804 //        Non-member operators
 805 //
 806 template<class T, class X>
 807 inline AliFmThreeVector<T>
 808 operator+ (const AliFmThreeVector<T>& v1, const AliFmThreeVector<X>& v2)
 809 {
 810     return AliFmThreeVector<T>(v1) += v2;
 811 }
 812
 813 template<class T, class X>
 814 inline AliFmThreeVector<T>
 815 operator- (const AliFmThreeVector<T>& v1, const AliFmThreeVector<X>& v2)
 816 {
 817     return AliFmThreeVector<T>(v1) -= v2;
 818 }
 819
 820 template<class T, class X>
 821 inline T operator* (const AliFmThreeVector<T>& v1, const AliFmThreeVector<X>& v2)
 822 {
 823     return AliFmThreeVector<T>(v1).Dot(v2);
 824 }
 825
 826 template<class T>
 827 inline AliFmThreeVector<T> operator* (const AliFmThreeVector<T>& v, double c)
 828 {
 829     return AliFmThreeVector<T>(v) *= c;
 830 }
 831
 832 template<class T>
 833 inline AliFmThreeVector<T> operator* (double c, const AliFmThreeVector<T>& v)
 834 {
 835     return AliFmThreeVector<T>(v) *= c;
 836 }
 837
 838 template<class T, class X>
 839 inline AliFmThreeVector<T> operator/ (const AliFmThreeVector<T>& v, X c)
 840 {
 841     return AliFmThreeVector<T>(v) /= c;
 842 }
 843
 844 template<class T>
 845 ostream&  operator<<(ostream& os, const AliFmThreeVector<T>& v)
 846 {
 847     return os << v.x() << '\t' << v.y() << '\t' << v.z();
 848 }
 849
 850 template<class T>
 851 istream&  operator>>(istream& is, AliFmThreeVector<T>& v)
 852 {
 853     T  x, y, z;
 854     is >> x >> y >> z;
 855     v.SetX(x);
 856     v.SetY(y);
 857     v.SetZ(z);
 858     return is;
 859 }
 860 #endif /* ! __CINT__ */
 861 #endif