PWG2/FEMTOSCOPY/AliFemto/AliFmThreeVector.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 #ifndef ST_THREE_VECTOR_HH
  76 #define ST_THREE_VECTOR_HH
  77 #ifdef __ROOT__
  78 #include "Rtypes.h"
  79 #endif
  80 #ifndef __CINT__
  81 #include <iostream>
  82 #include <fstream>
  83 #include <math.h>
  84 #ifdef GNU_GCC
  85 #    include <stddef.h>
  86 #endif
  87 #if defined (__SUNPRO_CC) && __SUNPRO_CC < 0x500
  88 #    include <stdcomp.h>
  89 #endif
  90 #ifndef ST_NO_EXCEPTIONS
  91 #    include <stdexcept>
  92 #    if !defined(ST_NO_NAMESPACES)
  93 using std::out_of_range;
  94 #    endif
  95 #endif
  96 #endif // __CINT__
  97
  98 #ifdef WIN32
  99 #include "gcc2vs.h"
 100 #endif
 101
 102 class TRootIOCtor;//nic nie rozumiem
 103 using namespace std;
 104
 105
 106 template<class T> class AliFmThreeVector {
 107 public:
 108     AliFmThreeVector(T = 0, T = 0, T = 0);
 109   //                     ROOT_VERSION(5,03,01)
 110 #if ROOT_VERSION_CODE >= 328449
 111    AliFmThreeVector(TRootIOCtor*) : mX1(0), mX2(0), mX3(0) {}
 112 #endif
 113     virtual ~AliFmThreeVector();
 114
 115 #ifndef ST_NO_MEMBER_TEMPLATES
 116     template<class X> AliFmThreeVector(const AliFmThreeVector<X>&);
 117     template<class X> AliFmThreeVector(const X*);
 118     template<class X> AliFmThreeVector<T>& operator=(const AliFmThreeVector<X>&);
 119     // AliFmThreeVector(const AliFmThreeVector<T>&);                use default
 120     // AliFmThreeVector<T>& operator=(const AliFmThreeVector<T>&);  use default
 121 #else
 122     AliFmThreeVector(const AliFmThreeVector<float>&);
 123     AliFmThreeVector(const AliFmThreeVector<double>&);
 124
 125     AliFmThreeVector(const float*);
 126     AliFmThreeVector(const double*);
 127
 128     AliFmThreeVector<T>& operator=(const AliFmThreeVector<float>&);
 129     AliFmThreeVector<T>& operator=(const AliFmThreeVector<double>&);
 130 #endif
 131
 132     void setX(T);
 133     void setY(T);
 134     void setZ(T);
 135
 136     void setPhi(T);
 137     void setTheta(T);
 138     void setMag(T);
 139     void setMagnitude(T);
 140
 141     T   x()                        const;
 142     T   y()                        const;
 143     T   z()                        const;
 144     T   theta()                    const;
 145     T   cosTheta()                 const;
 146     T   phi()                      const;
 147     T   perp()                     const;
 148     T   perp2()                    const;
 149     T   magnitude()                const;
 150     T   mag()                      const;
 151     T   mag2()                     const;
 152     T   pseudoRapidity()           const;
 153     T   operator() (size_t)        const;
 154     T   operator[] (size_t)        const;
 155
 156     T&  operator() (size_t);
 157     T&  operator[] (size_t);
 158
 159     T   massHypothesis(T mass)     const;
 160
 161     AliFmThreeVector<T>  unit()       const;
 162     AliFmThreeVector<T>  orthogonal() const;
 163
 164     void  rotateX(T);
 165     void  rotateY(T);
 166     void  rotateZ(T);
 167
 168     AliFmThreeVector<T>  operator- ();
 169     AliFmThreeVector<T>  operator+ ();
 170     AliFmThreeVector<T>& operator*= (double);
 171     AliFmThreeVector<T>& operator/= (double);
 172     AliFmThreeVector<T>  pseudoProduct(double,double,double) const;
 173
 174 #ifndef ST_NO_MEMBER_TEMPLATES
 175     template<class X> T                angle(const AliFmThreeVector<X>&) const;
 176     template<class X> AliFmThreeVector<T> cross(const AliFmThreeVector<X>&) const;
 177     template<class X> T                dot  (const AliFmThreeVector<X>&) const;
 178     template<class X> AliFmThreeVector<T> pseudoProduct(const AliFmThreeVector<X>&) const;
 179
 180     template<class X> bool operator == (const AliFmThreeVector<X>& v) const;
 181     template<class X> bool operator != (const AliFmThreeVector<X>& v) const;
 182
 183     template<class X> AliFmThreeVector<T>& operator+= (const AliFmThreeVector<X>&);
 184     template<class X> AliFmThreeVector<T>& operator-= (const AliFmThreeVector<X>&);
 185 #else
 186     T                angle(const AliFmThreeVector<float>&) const;
 187     AliFmThreeVector<T> cross(const AliFmThreeVector<float>&) const;
 188     T                dot  (const AliFmThreeVector<float>&) const;
 189     AliFmThreeVector<T> pseudoProduct(const AliFmThreeVector<float>&) const;
 190
 191     T                angle(const AliFmThreeVector<double>&) const;
 192     T                dot  (const AliFmThreeVector<double>&) const;
 193     AliFmThreeVector<T> cross(const AliFmThreeVector<double>&) const;
 194     AliFmThreeVector<T> pseudoProduct(const AliFmThreeVector<double>&) const;
 195
 196     bool operator == (const AliFmThreeVector<float>& v) const;
 197     bool operator != (const AliFmThreeVector<float>& v) const;
 198     AliFmThreeVector<T>& operator+= (const AliFmThreeVector<float>&);
 199     AliFmThreeVector<T>& operator-= (const AliFmThreeVector<float>&);
 200
 201     bool operator == (const AliFmThreeVector<double>& v) const;
 202     bool operator != (const AliFmThreeVector<double>& v) const;
 203     AliFmThreeVector<T>& operator+= (const AliFmThreeVector<double>&);
 204     AliFmThreeVector<T>& operator-= (const AliFmThreeVector<double>&);
 205 #endif
 206   int             valid(double world = 1.e+5) const;
 207     int               bad(double world = 1.e+5) const;
 208 protected:
 209     T    mX1, mX2, mX3;
 210 #ifdef __ROOT__
 211   ClassDef(AliFmThreeVector,3)
 212 #endif /* __ROOT__ */
 213 };
 214
 215 #ifndef __CINT__
 216 //
 217 //        Implementation of member functions
 218 //
 219 template<class T>
 220 inline AliFmThreeVector<T>::AliFmThreeVector(T x, T y, T z)
 221     : mX1(x), mX2(y), mX3(z) {/* nop */}
 222 template<class T>
 223 inline AliFmThreeVector<T>::~AliFmThreeVector() {/* nop */}
 224
 225 template<class T>
 226 inline void AliFmThreeVector<T>::setX(T x) {mX1 = x;}
 227
 228 template<class T>
 229 inline void AliFmThreeVector<T>::setY(T y) {mX2 = y;}
 230
 231 template<class T>
 232 inline void AliFmThreeVector<T>::setZ(T z) {mX3 = z;}
 233
 234 template<class T>
 235 void AliFmThreeVector<T>::setPhi(T angle)
 236 {
 237     double  r = magnitude();
 238     double th = theta();
 239
 240     mX1 = r*sin(th)*cos(angle);
 241     mX2 = r*sin(th)*sin(angle);
 242 }
 243
 244 template <class T>
 245 void AliFmThreeVector<T>::setTheta(T angle)
 246 {
 247     double r  = magnitude();
 248     double ph = phi();
 249
 250     mX1 = r*sin(angle)*cos(ph);
 251     mX2 = r*sin(angle)*sin(ph);
 252     mX3 = r*cos(angle);
 253 }
 254
 255 template <class T>
 256 void AliFmThreeVector<T>::setMagnitude(T r)
 257 {
 258     double th = theta();
 259     double ph = phi();
 260
 261     mX1 = r*sin(th)*cos(ph);
 262     mX2 = r*sin(th)*sin(ph);
 263     mX3 = r*cos(th);
 264 }
 265
 266 template <class T>
 267 void AliFmThreeVector<T>::setMag(T mag)
 268 {
 269     setMagnitude(mag);
 270 }
 271
 272 template<class T>
 273 inline T AliFmThreeVector<T>::x() const {return mX1;}
 274
 275 template<class T>
 276 inline T AliFmThreeVector<T>::y() const {return mX2;}
 277
 278 template<class T>
 279 inline T AliFmThreeVector<T>::z() const {return mX3;}
 280
 281 template<class T>
 282 inline T AliFmThreeVector<T>::theta() const
 283 {
 284   return acos(cosTheta());
 285 }
 286
 287 template<class T>
 288 inline T AliFmThreeVector<T>::cosTheta() const
 289 {
 290   return mX3/(mag()+1e-20);
 291 }
 292
 293 template<class T>
 294 inline T AliFmThreeVector<T>::phi() const
 295 {
 296     return atan2(mX2,mX1);
 297 }
 298
 299 template<class T>
 300 inline T AliFmThreeVector<T>::pseudoRapidity() const
 301 {
 302     //
 303     // change code to more optimal:
 304     // double m = mag();
 305     // return 0.5*::log( (m+z())/(m-z()) );
 306     double tmp = tan(theta()/2.); if (tmp <=0.) return 1e20;
 307     return -::log(tmp);
 308 }
 309
 310 template<class T>
 311 inline AliFmThreeVector<T> AliFmThreeVector<T>::unit() const
 312 {
 313     double tmp = mag(); if (tmp<=0.) tmp = 1e-20;
 314     return *this/tmp;
 315 }
 316
 317 template <class T>
 318 T AliFmThreeVector<T>::massHypothesis(T mass) const
 319 {
 320     return ::sqrt((*this)*(*this) + mass*mass);
 321 }
 322
 323 template <class T>
 324 AliFmThreeVector<T> AliFmThreeVector<T>::orthogonal() const
 325 {
 326     // Direct copy from CLHEP--it is probably better to
 327     // use your own dot/cross product code...
 328     double x = (mX1 < 0.0) ? -mX1 : mX1;
 329     double y = (mX2 < 0.0) ? -mX2 : mX2;
 330     double z = (mX3 < 0.0) ? -mX3 : mX3;
 331
 332     if(x<y)
 333         return x < z ? AliFmThreeVector<T>(0,mX3,-mX2) :  AliFmThreeVector<T>(mX2,-mX1,0);
 334     else
 335         return  mX2 < mX3 ? AliFmThreeVector<T>(-mX3,0,mX1) :  AliFmThreeVector<T>(mX2,-mX1,0);
 336 }
 337
 338 template <class T>
 339 void AliFmThreeVector<T>::rotateX(T angle)
 340 {
 341     // may in the future make use of the AliFmRotation class!
 342     double yPrime = cos(angle)*mX2 - sin(angle)*mX3;
 343     double zPrime = sin(angle)*mX2 + cos(angle)*mX3;
 344
 345     mX2 = yPrime;
 346     mX3 = zPrime;
 347 }
 348
 349 template <class T>
 350 void AliFmThreeVector<T>::rotateY(T angle)
 351 {
 352     // may in the future make use of the AliFmRotation class!
 353     double zPrime = cos(angle)*mX3 - sin(angle)*mX1;
 354     double xPrime = sin(angle)*mX3 + cos(angle)*mX1;
 355
 356     mX1 = xPrime;
 357     mX3 = zPrime;
 358 }
 359
 360 template <class T>
 361 void AliFmThreeVector<T>::rotateZ(T angle)
 362 {
 363     // may in the future make use of the AliFmRotation class!
 364     double xPrime = cos(angle)*mX1 - sin(angle)*mX2;
 365     double yPrime = sin(angle)*mX1 + cos(angle)*mX2;
 366
 367     mX1 = xPrime;
 368     mX2 = yPrime;
 369 }
 370
 371 template<class T>
 372 inline T AliFmThreeVector<T>::perp() const
 373 {
 374     return ::sqrt(mX1*mX1+mX2*mX2);
 375 }
 376
 377 template<class T>
 378 inline T AliFmThreeVector<T>::perp2() const
 379 {
 380     return mX1*mX1+mX2*mX2;
 381 }
 382
 383 template<class T>
 384 inline T AliFmThreeVector<T>::magnitude() const
 385 {
 386     return mag();
 387 }
 388
 389 template<class T>
 390 inline T AliFmThreeVector<T>::mag() const
 391 {
 392     return ::sqrt(mX1*mX1+mX2*mX2+mX3*mX3);
 393 }
 394
 395 template<class T>
 396 inline T AliFmThreeVector<T>::mag2() const
 397 {
 398     return mX1*mX1+mX2*mX2+mX3*mX3;
 399 }
 400
 401 template<class T>
 402 inline T AliFmThreeVector<T>::operator() (size_t i) const
 403 {
 404     if (0 <=i && i <= 2)  return (&mX1)[i];
 405 #ifndef ST_NO_EXCEPTIONS
 406     throw out_of_range("AliFmThreeVector<T>::operator(): bad index");
 407 #else
 408     cerr << "AliFmThreeVector<T>::operator(): bad index" << endl;
 409 #endif
 410     return 0;
 411 }
 412
 413 template<class T>
 414 inline T& AliFmThreeVector<T>::operator() (size_t i)
 415 {
 416     if (0 <=i && i <= 2)  return (&mX1)[i];
 417 #ifndef ST_NO_EXCEPTIONS
 418     throw out_of_range("AliFmThreeVector<T>::operator(): bad index");
 419 #else
 420     cerr << "AliFmThreeVector<T>::operator(): bad index" << endl;
 421 #endif
 422     return mX1;
 423 }
 424
 425 template<class T>
 426 inline T AliFmThreeVector<T>::operator[] (size_t i) const
 427 {
 428     if (0 <=i && i <= 2)  return (&mX1)[i];
 429 #ifndef ST_NO_EXCEPTIONS
 430       throw out_of_range("AliFmThreeVector<T>::operator[]: bad index");
 431 #else
 432       cerr << "AliFmThreeVector<T>::operator[]: bad index" << endl;
 433 #endif
 434       return 0;
 435 }
 436
 437 template<class T>
 438 inline T &AliFmThreeVector<T>::operator[] (size_t i)
 439 {
 440     if (0 <=i && i <= 2)  return (&mX1)[i];
 441 #ifndef ST_NO_EXCEPTIONS
 442       throw out_of_range("AliFmThreeVector<T>::operator[]: bad index");
 443 #else
 444       cerr << "AliFmThreeVector<T>::operator[]: bad index" << endl;
 445 #endif
 446       return mX1;
 447 }
 448
 449 template<class T>
 450 inline AliFmThreeVector<T>& AliFmThreeVector<T>::operator*= (double c)
 451 {
 452     mX1 *= c; mX2 *= c; mX3 *= c;
 453     return *this;
 454 }
 455
 456 template<class T>
 457 inline AliFmThreeVector<T>& AliFmThreeVector<T>::operator/= (double c)
 458 {
 459     mX1 /= c; mX2 /= c; mX3 /= c;
 460     return *this;
 461 }
 462
 463 template<class T>
 464 inline AliFmThreeVector<T>
 465 AliFmThreeVector<T>::pseudoProduct(double x,double y,double z) const
 466 {
 467     return AliFmThreeVector<T>(mX1*x,mX2*y,mX3*z);
 468 }
 469
 470 template<class T>
 471 AliFmThreeVector<T> AliFmThreeVector<T>::operator- ()
 472 {
 473     return AliFmThreeVector<T>(-mX1, -mX2, -mX3);
 474 }
 475
 476 template<class T>
 477 AliFmThreeVector<T> AliFmThreeVector<T>::operator+ ()
 478 {
 479     return *this;
 480 }
 481
 482 #ifndef ST_NO_MEMBER_TEMPLATES
 483 #ifndef WIN32
 484
 485 template<class T>
 486 template<class X>
 487 inline AliFmThreeVector<T>::AliFmThreeVector(const AliFmThreeVector<X>& v)
 488     : mX1(v.x()), mX2(v.y()), mX3(v.z()) {/* nop */}
 489
 490 template<class T>
 491 template<class X>
 492 inline AliFmThreeVector<T>::AliFmThreeVector(const X *a)
 493 {
 494     mX1 = a[0];
 495     mX2 = a[1];
 496     mX3 = a[2];
 497 }
 498
 499 template<class T>
 500 template<class X>
 501 inline AliFmThreeVector<T>&
 502 AliFmThreeVector<T>::operator=(const AliFmThreeVector<X>& v)
 503 {
 504     mX1 = v.x();  mX2 = v.y();  mX3 = v.z();
 505     return *this;
 506 }
 507
 508 template<class T>
 509 template<class X>
 510 inline bool AliFmThreeVector<T>::operator== (const AliFmThreeVector<X>& v) const
 511 {
 512     return mX1 == v.x() && mX2 == v.y() && mX3 == v.z();
 513 }
 514
 515 template<class T>
 516 template<class X>
 517 inline bool AliFmThreeVector<T>::operator!= (const AliFmThreeVector<X>& v) const
 518 {
 519     return !(*this == v);
 520 }
 521
 522 template<class T>
 523 template<class X>
 524 inline AliFmThreeVector<T>&
 525 AliFmThreeVector<T>::operator+= (const AliFmThreeVector<X>& v)
 526 {
 527     mX1 += v.x(); mX2 += v.y(); mX3 += v.z();
 528     return *this;
 529 }
 530
 531 template<class T>
 532 template<class X>
 533 inline AliFmThreeVector<T>&
 534 AliFmThreeVector<T>::operator-= (const AliFmThreeVector<X>& v)
 535 {
 536     mX1 -= v.x(); mX2 -= v.y(); mX3 -= v.z();
 537     return *this;
 538 }
 539
 540 template<class T>
 541 template<class X>
 542 inline T AliFmThreeVector<T>::dot(const AliFmThreeVector<X>& v) const
 543 {
 544     return mX1*v.x() + mX2*v.y() + mX3*v.z();
 545 }
 546
 547 template<class T>
 548 template<class X>
 549 inline AliFmThreeVector<T>
 550 AliFmThreeVector<T>::cross(const AliFmThreeVector<X>& v) const
 551 {
 552     return AliFmThreeVector<T>(mX2*v.z() - mX3*v.y(),
 553                             mX3*v.x() - mX1*v.z(),
 554                             mX1*v.y() - mX2*v.x());
 555 }
 556
 557 template<class T>
 558 template<class X>
 559 inline T AliFmThreeVector<T>::angle(const AliFmThreeVector<X>& vec) const
 560 {
 561     double norm = this->mag2()*vec.mag2();
 562
 563     return norm > 0 ? acos(this->dot(vec)/(::sqrt(norm))) : 0;
 564 }
 565
 566 template<class T>
 567 template<class X>
 568 inline AliFmThreeVector<T>
 569 AliFmThreeVector<T>::pseudoProduct(const AliFmThreeVector<X>& v) const
 570 {
 571     return this->pseudoProduct(v.x(),v.y(),v.z());
 572 }
 573
 574 #endif
 575 #else
 576
 577 template<class T>
 578 inline AliFmThreeVector<T>::AliFmThreeVector(const AliFmThreeVector<float>& v)
 579     : mX1(v.x()), mX2(v.y()), mX3(v.z()) {/* nop */}
 580
 581 template<class T>
 582 inline AliFmThreeVector<T>::AliFmThreeVector(const AliFmThreeVector<double>& v)
 583     : mX1(v.x()), mX2(v.y()), mX3(v.z()) {/* nop */}
 584
 585 template<class T>
 586 inline AliFmThreeVector<T>::AliFmThreeVector(const float *a)
 587 {
 588     mX1 = a[0];
 589     mX2 = a[1];
 590     mX3 = a[2];
 591 }
 592
 593 template<class T>
 594 inline AliFmThreeVector<T>::AliFmThreeVector(const double *a)
 595 {
 596     mX1 = a[0];
 597     mX2 = a[1];
 598     mX3 = a[2];
 599 }
 600
 601 template<class T>
 602 inline AliFmThreeVector<T>&
 603 AliFmThreeVector<T>::operator=(const AliFmThreeVector<float>& v)
 604 {
 605     mX1 = v.x();  mX2 = v.y();  mX3 = v.z();
 606     return *this;
 607 }
 608
 609 template<class T>
 610 inline AliFmThreeVector<T>&
 611 AliFmThreeVector<T>::operator=(const AliFmThreeVector<double>& v)
 612 {
 613     mX1 = v.x();  mX2 = v.y();  mX3 = v.z();
 614     return *this;
 615 }
 616
 617 template<class T>
 618 inline bool
 619 AliFmThreeVector<T>::operator== (const AliFmThreeVector<float>& v) const
 620 {
 621     return mX1 == v.x() && mX2 == v.y() && mX3 == v.z();
 622 }
 623
 624 template<class T>
 625 inline bool
 626 AliFmThreeVector<T>::operator== (const AliFmThreeVector<double>& v) const
 627 {
 628     return mX1 == v.x() && mX2 == v.y() && mX3 == v.z();
 629 }
 630
 631 template<class T>
 632 inline bool
 633 AliFmThreeVector<T>::operator!= (const AliFmThreeVector<float>& v) const
 634 {
 635     return !(*this == v);
 636 }
 637
 638 template<class T>
 639 inline bool
 640 AliFmThreeVector<T>::operator!= (const AliFmThreeVector<double>& v) const
 641 {
 642     return !(*this == v);
 643 }
 644
 645 template<class T>
 646 inline AliFmThreeVector<T>&
 647 AliFmThreeVector<T>::operator+= (const AliFmThreeVector<float>& v)
 648 {
 649     mX1 += v.x(); mX2 += v.y(); mX3 += v.z();
 650     return *this;
 651 }
 652
 653 template<class T>
 654 inline AliFmThreeVector<T>&
 655 AliFmThreeVector<T>::operator+= (const AliFmThreeVector<double>& v)
 656 {
 657     mX1 += v.x(); mX2 += v.y(); mX3 += v.z();
 658     return *this;
 659 }
 660
 661 template<class T>
 662 inline AliFmThreeVector<T>&
 663 AliFmThreeVector<T>::operator-= (const AliFmThreeVector<float>& v)
 664 {
 665     mX1 -= v.x(); mX2 -= v.y(); mX3 -= v.z();
 666     return *this;
 667 }
 668
 669 template<class T>
 670 inline AliFmThreeVector<T>&
 671 AliFmThreeVector<T>::operator-= (const AliFmThreeVector<double>& v)
 672 {
 673     mX1 -= v.x(); mX2 -= v.y(); mX3 -= v.z();
 674     return *this;
 675 }
 676
 677 template<class T>
 678 inline T AliFmThreeVector<T>::dot(const AliFmThreeVector<float>& v) const
 679 {
 680     return mX1*v.x() + mX2*v.y() + mX3*v.z();
 681 }
 682
 683 template<class T>
 684 inline T AliFmThreeVector<T>::dot(const AliFmThreeVector<double>& v) const
 685 {
 686     return mX1*v.x() + mX2*v.y() + mX3*v.z();
 687 }
 688
 689 template<class T>
 690 inline AliFmThreeVector<T>
 691 AliFmThreeVector<T>::cross(const AliFmThreeVector<float>& v) const
 692 {
 693     return AliFmThreeVector<T>(mX2*v.z() - mX3*v.y(),
 694                             mX3*v.x() - mX1*v.z(),
 695                             mX1*v.y() - mX2*v.x());
 696 }
 697
 698 template<class T>
 699 inline AliFmThreeVector<T>
 700 AliFmThreeVector<T>::cross(const AliFmThreeVector<double>& v) const
 701 {
 702     return AliFmThreeVector<T>(mX2*v.z() - mX3*v.y(),
 703                             mX3*v.x() - mX1*v.z(),
 704                             mX1*v.y() - mX2*v.x());
 705 }
 706
 707 template<class T>
 708 inline T AliFmThreeVector<T>::angle(const AliFmThreeVector<float>& v) const
 709 {
 710     double tmp = mag()*v.mag(); if (tmp <=0) tmp = 1e-20;
 711     return acos(this->dot(v)/tmp);
 712 }
 713
 714 template<class T>
 715 inline T AliFmThreeVector<T>::angle(const AliFmThreeVector<double>& v) const
 716 {
 717     double tmp = mag()*v.mag(); if (tmp <=0) tmp = 1e-20;
 718     return acos(this->dot(v)/tmp);
 719 }
 720
 721 template<class T>
 722 inline AliFmThreeVector<T>
 723 AliFmThreeVector<T>::pseudoProduct(const AliFmThreeVector<float>& v) const
 724 {
 725     return this->pseudoProduct(v.x(),v.y(),v.z());
 726 }
 727
 728 template<class T>
 729 inline AliFmThreeVector<T>
 730 AliFmThreeVector<T>::pseudoProduct(const AliFmThreeVector<double>& v) const
 731 {
 732     return this->pseudoProduct(v.x(),v.y(),v.z());
 733 }
 734 #endif  // ST_NO_MEMBER_TEMPLATES
 735 template<class T>
 736 inline int
 737 AliFmThreeVector<T>::valid(double world) const  {return !bad(world);}
 738
 739 template<class T>
 740 inline int
 741 AliFmThreeVector<T>::bad(double world) const
 742 {
 743   for (int i=0;i<3;i++) {
 744           if (!finite((&mX1)[i])      ) return 10+i;
 745           if ( fabs  ((&mX1)[i])>world) return 20+i;
 746   }
 747   return 0;
 748 }
 749 #endif /*! __CINT__ */
 750 #ifdef __CINT__
 751 template<> float abs(const AliFmThreeVector<float>& v);
 752 template<> double abs(const AliFmThreeVector<double>& v);
 753 template<> AliFmThreeVector<double> cross_product(const AliFmThreeVector<double>& v1, const AliFmThreeVector<double>& v2);
 754 template<> AliFmThreeVector<float>  cross_product(const AliFmThreeVector<float>&  v1, const AliFmThreeVector<float>& v2);
 755 template<> AliFmThreeVector<double> cross_product(const AliFmThreeVector<float>&  v1, const AliFmThreeVector<double>& v2);
 756 template<> AliFmThreeVector<double> cross_product(const AliFmThreeVector<double>& v1, const AliFmThreeVector<float>& v2);
 757 template<> AliFmThreeVector<double> operator+ (const AliFmThreeVector<double>& v1, const AliFmThreeVector<double>& v2);
 758 template<> AliFmThreeVector<float>  operator+ (const AliFmThreeVector<float>&  v1, const AliFmThreeVector<float>& v2);
 759 template<> AliFmThreeVector<double> operator+ (const AliFmThreeVector<double>& v1, const AliFmThreeVector<float>& v2);
 760 template<> AliFmThreeVector<double> operator+ (const AliFmThreeVector<float>&  v1, const AliFmThreeVector<double>& v2);
 761 template<> AliFmThreeVector<double> operator- (const AliFmThreeVector<double>& v1, const AliFmThreeVector<double>& v2);
 762 template<> AliFmThreeVector<float>  operator- (const AliFmThreeVector<float>&  v1, const AliFmThreeVector<float>& v2);
 763 template<> AliFmThreeVector<double> operator- (const AliFmThreeVector<double>& v1, const AliFmThreeVector<float>& v2);
 764 template<> AliFmThreeVector<double> operator- (const AliFmThreeVector<float>&  v1, const AliFmThreeVector<double>& v2);
 765 template<> AliFmThreeVector<double> operator* (const AliFmThreeVector<double>& v1, const AliFmThreeVector<double>& v2);
 766 template<> AliFmThreeVector<float>  operator* (const AliFmThreeVector<float>&  v1, const AliFmThreeVector<float>& v2);
 767 template<> AliFmThreeVector<double> operator* (const AliFmThreeVector<double>& v1, const AliFmThreeVector<float>& v2);
 768 template<> AliFmThreeVector<double> operator* (const AliFmThreeVector<float>&  v1, const AliFmThreeVector<double>& v2);
 769 template<> AliFmThreeVector<double> operator* (const double                 v1, const AliFmThreeVector<float>& v2);
 770 template<> AliFmThreeVector<double> operator* (const AliFmThreeVector<float>&  v1, const double v2);
 771 template<> AliFmThreeVector<double> operator* (const double                 v1, const AliFmThreeVector<double>& v2);
 772 template<> AliFmThreeVector<double> operator* (const AliFmThreeVector<double>& v1, const double v2);
 773 template<> AliFmThreeVector<double> operator/ (const AliFmThreeVector<double>& v1, const AliFmThreeVector<double>& v2);
 774 template<> AliFmThreeVector<float>  operator/ (const AliFmThreeVector<float>&  v1, const AliFmThreeVector<float>& v2);
 775 template<> AliFmThreeVector<double> operator/ (const AliFmThreeVector<double>& v1, const AliFmThreeVector<float>& v2);
 776 template<> AliFmThreeVector<double> operator/ (const AliFmThreeVector<float>&  v1, const AliFmThreeVector<double>& v2);
 777 template<> AliFmThreeVector<double> operator/ (const                double  v1, const AliFmThreeVector<double>& v2);
 778 template<> AliFmThreeVector<float>  operator/ (const                double  v1, const AliFmThreeVector<float>& v2);
 779 template<> AliFmThreeVector<double> operator/ (const AliFmThreeVector<double>& v1, const double v2);
 780 template<> AliFmThreeVector<double> operator/ (const AliFmThreeVector<float>&  v1, const double v2);
 781 template<> istream&  operator>>(istream& is,const AliFmThreeVector<double>& v);
 782 template<> istream&  operator>>(istream& is,const AliFmThreeVector<float>& v);
 783 template<> ostream&  operator<<(ostream& os,const AliFmThreeVector<double>& v);
 784 template<> ostream&  operator<<(ostream& os,const AliFmThreeVector<float>& v);
 785 #else
 786 //
 787 //        Non-member functions
 788 //
 789 template<class T>
 790 inline T abs(const AliFmThreeVector<T>& v) {return v.mag();}
 791
 792 template<class T, class X>
 793 inline AliFmThreeVector<T>
 794 cross_product(const AliFmThreeVector<T>& v1, const AliFmThreeVector<X>& v2)
 795 {
 796     return v1.cross(v2);
 797 }
 798
 799
 800 //
 801 //        Non-member operators
 802 //
 803 template<class T, class X>
 804 inline AliFmThreeVector<T>
 805 operator+ (const AliFmThreeVector<T>& v1, const AliFmThreeVector<X>& v2)
 806 {
 807     return AliFmThreeVector<T>(v1) += v2;
 808 }
 809
 810 template<class T, class X>
 811 inline AliFmThreeVector<T>
 812 operator- (const AliFmThreeVector<T>& v1, const AliFmThreeVector<X>& v2)
 813 {
 814     return AliFmThreeVector<T>(v1) -= v2;
 815 }
 816
 817 template<class T, class X>
 818 inline T operator* (const AliFmThreeVector<T>& v1, const AliFmThreeVector<X>& v2)
 819 {
 820     return AliFmThreeVector<T>(v1).dot(v2);
 821 }
 822
 823 template<class T>
 824 inline AliFmThreeVector<T> operator* (const AliFmThreeVector<T>& v, double c)
 825 {
 826     return AliFmThreeVector<T>(v) *= c;
 827 }
 828
 829 template<class T>
 830 inline AliFmThreeVector<T> operator* (double c, const AliFmThreeVector<T>& v)
 831 {
 832     return AliFmThreeVector<T>(v) *= c;
 833 }
 834
 835 template<class T, class X>
 836 inline AliFmThreeVector<T> operator/ (const AliFmThreeVector<T>& v, X c)
 837 {
 838     return AliFmThreeVector<T>(v) /= c;
 839 }
 840
 841 template<class T>
 842 ostream&  operator<<(ostream& os, const AliFmThreeVector<T>& v)
 843 {
 844     return os << v.x() << '\t' << v.y() << '\t' << v.z();
 845 }
 846
 847 template<class T>
 848 istream&  operator>>(istream& is, AliFmThreeVector<T>& v)
 849 {
 850     T  x, y, z;
 851     is >> x >> y >> z;
 852     v.setX(x);
 853     v.setY(y);
 854     v.setZ(z);
 855     return is;
 856 }
 857 #endif /* ! __CINT__ */
 858 #endif