VecGeom/base/LorentzVector.h

0001 /// \file LorentzVector.h
0002 /// \author Federico Carminati (based on existing Vector3D.h)
0003
0004 #ifndef VECGEOM_BASE_LORENTZVECTOR_H_
0005 #define VECGEOM_BASE_LORENTZVECTOR_H_
0006
0007 #include "VecGeom/base/Global.h"
0008 #include "VecGeom/base/Vector3D.h"
0009
0010 #include "VecGeom/base/AlignedBase.h"
0011 #include "VecGeom/base/Vector3D.h"
0012 #include "VecGeom/base/RNG.h"
0013
0014 #include <cstdlib>
0015 #include <ostream>
0016 #include <string>
0017 #include <iostream>
0018
0019 namespace vecgeom {
0020
0021 VECGEOM_DEVICE_FORWARD_DECLARE(template <typename T> class LorentzVector;);
0022
0023 inline namespace VECGEOM_IMPL_NAMESPACE {
0024
0025 template <typename T>
0026 class LorentzRotation;
0027
0028 /**
0029  * @brief Lorentz dimensional vector class supporting most arithmetic operations.
0030  * @details If vector acceleration is enabled, the scalar template instantiation
0031  *          will use vector instructions for operations when possible.
0032  */
0033 template <typename T>
0034 class LorentzVector : public AlignedBase {
0035
0036   typedef LorentzVector<T> VecType;
0037
0038 private:
0039   T fVec[4];
0040
0041 public:
0042   VECCORE_ATT_HOST_DEVICE
0043   VECGEOM_FORCE_INLINE
0044   LorentzVector(const T a, const T b, const T c, const T d) : fVec{a, b, c, d} {}
0045
0046   VECCORE_ATT_HOST_DEVICE
0047   VECGEOM_FORCE_INLINE
0048   LorentzVector() : fVec{0, 0, 0, 0} {}
0049
0050   VECCORE_ATT_HOST_DEVICE
0051   VECGEOM_FORCE_INLINE
0052   LorentzVector(const T a) : fVec{a, a, a, a} {}
0053
0054   template <typename U>
0055   VECGEOM_FORCE_INLINE
0056   VECCORE_ATT_HOST_DEVICE
0057   LorentzVector(LorentzVector<U> const &other) : fVec{other[0], other[1], other[2], other[3]}
0058   {
0059   }
0060
0061   template <typename U>
0062   VECGEOM_FORCE_INLINE
0063   VECCORE_ATT_HOST_DEVICE
0064   LorentzVector(Vector3D<U> const &other, T t) : fVec{other[0], other[1], other[2], t}
0065   {
0066   }
0067
0068   VECCORE_ATT_HOST_DEVICE
0069   VECGEOM_FORCE_INLINE
0070   LorentzVector &operator=(LorentzVector const &other)
0071   {
0072     if (this != &other) {
0073       fVec[0] = other[0];
0074       fVec[1] = other[1];
0075       fVec[2] = other[2];
0076       fVec[3] = other[3];
0077     }
0078     return *this;
0079   }
0080
0081   VECCORE_ATT_HOST_DEVICE
0082   VECGEOM_FORCE_INLINE
0083   operator Vector3D<T> &() { return reinterpret_cast<Vector3D<T> &>(*this); }
0084
0085   /**
0086    * Constructs a vector from an std::string of the same format as output by the
0087    * "<<"-operator for outstreams.
0088    * @param str String formatted as "(%d, %d, %d %d)".
0089    */
0090   VECCORE_ATT_HOST
0091   LorentzVector(std::string const &str)
0092   {
0093     int begin = str.find("(") + 1, end = str.find(",") - 1;
0094     fVec[0] = std::atof(str.substr(begin, end - begin + 1).c_str());
0095     begin   = end + 2;
0096     end     = str.find(",", begin) - 1;
0097     fVec[1] = std::atof(str.substr(begin, end - begin + 1).c_str());
0098     begin   = end + 2;
0099     end     = str.find(",", begin) - 1;
0100     fVec[2] = std::atof(str.substr(begin, end - begin + 1).c_str());
0101     begin   = end + 2;
0102     end     = str.find(")", begin) - 1;
0103     fVec[3] = std::atof(str.substr(begin, end - begin + 1).c_str());
0104   }
0105
0106   /**
0107    * Contains no check for correct indexing to avoid impairing performance.
0108    * @param index Index of content in the range [0-2].
0109    */
0110   VECCORE_ATT_HOST_DEVICE
0111   VECGEOM_FORCE_INLINE
0112   T &operator[](const int index) { return fVec[index]; }
0113
0114   /**
0115    * Contains no check for correct indexing to avoid impairing performance.
0116    * @param index Index of content in the range [0-2].
0117    */
0118   VECCORE_ATT_HOST_DEVICE
0119   VECGEOM_FORCE_INLINE
0120   T const &operator[](const int index) const { return fVec[index]; }
0121
0122   VECCORE_ATT_HOST_DEVICE
0123   VECGEOM_FORCE_INLINE
0124   T &x() { return fVec[0]; }
0125
0126   VECCORE_ATT_HOST_DEVICE
0127   VECGEOM_FORCE_INLINE
0128   T const &x() const { return fVec[0]; }
0129
0130   VECCORE_ATT_HOST_DEVICE
0131   VECGEOM_FORCE_INLINE
0132   T &y() { return fVec[1]; }
0133
0134   VECCORE_ATT_HOST_DEVICE
0135   VECGEOM_FORCE_INLINE
0136   T const &y() const { return fVec[1]; }
0137
0138   VECCORE_ATT_HOST_DEVICE
0139   VECGEOM_FORCE_INLINE
0140   T &z() { return fVec[2]; }
0141
0142   VECCORE_ATT_HOST_DEVICE
0143   VECGEOM_FORCE_INLINE
0144   T const &z() const { return fVec[2]; }
0145
0146   VECCORE_ATT_HOST_DEVICE
0147   VECGEOM_FORCE_INLINE
0148   T &t() { return fVec[3]; }
0149
0150   VECCORE_ATT_HOST_DEVICE
0151   VECGEOM_FORCE_INLINE
0152   T const &t() const { return fVec[3]; }
0153
0154   VECCORE_ATT_HOST_DEVICE
0155   VECGEOM_FORCE_INLINE
0156   T &px() { return fVec[0]; }
0157
0158   VECCORE_ATT_HOST_DEVICE
0159   VECGEOM_FORCE_INLINE
0160   T const &px() const { return fVec[0]; }
0161
0162   VECCORE_ATT_HOST_DEVICE
0163   VECGEOM_FORCE_INLINE
0164   T &py() { return fVec[1]; }
0165
0166   VECCORE_ATT_HOST_DEVICE
0167   VECGEOM_FORCE_INLINE
0168   T const &py() const { return fVec[1]; }
0169
0170   VECCORE_ATT_HOST_DEVICE
0171   VECGEOM_FORCE_INLINE
0172   T &pz() { return fVec[2]; }
0173
0174   VECCORE_ATT_HOST_DEVICE
0175   VECGEOM_FORCE_INLINE
0176   T const &pz() const { return fVec[2]; }
0177
0178   VECCORE_ATT_HOST_DEVICE
0179   VECGEOM_FORCE_INLINE
0180   T &e() { return fVec[3]; }
0181
0182   VECCORE_ATT_HOST_DEVICE
0183   VECGEOM_FORCE_INLINE
0184   T const e() const { return fVec[3]; }
0185
0186   VECCORE_ATT_HOST_DEVICE
0187   VECGEOM_FORCE_INLINE
0188   T const m2() const { return Mag2(); }
0189
0190   VECCORE_ATT_HOST_DEVICE
0191   VECGEOM_FORCE_INLINE
0192   T const m() const { return Mag(); }
0193
0194   VECCORE_ATT_HOST_DEVICE
0195   VECGEOM_FORCE_INLINE
0196   T const mt2() const { return e() * e() - pz() * pz(); }
0197
0198   VECCORE_ATT_HOST_DEVICE
0199   VECGEOM_FORCE_INLINE
0200   T const mt() const
0201   {
0202     T mmm = mt2();
0203     return mmm < 0 ? -Sqrt(-mmm) : Sqrt(mmm);
0204   }
0205
0206   VECCORE_ATT_HOST_DEVICE
0207   VECGEOM_FORCE_INLINE
0208   void Set(T const &a, T const &b, T const &c, T const &d)
0209   {
0210     fVec[0] = a;
0211     fVec[1] = b;
0212     fVec[2] = c;
0213     fVec[3] = d;
0214   }
0215
0216   template <typename U, typename V>
0217   VECGEOM_FORCE_INLINE
0218   VECCORE_ATT_HOST_DEVICE
0219   void Set(const Vector3D<U> &v, V t)
0220   {
0221     fVec[0] = v[0];
0222     fVec[1] = v[1];
0223     fVec[2] = v[2];
0224     fVec[3] = t;
0225   }
0226
0227   template <typename U>
0228   VECGEOM_FORCE_INLINE
0229   VECCORE_ATT_HOST_DEVICE
0230   void Set(const Vector3D<U> &v)
0231   {
0232     fVec[0] = v[0];
0233     fVec[1] = v[1];
0234     fVec[2] = v[2];
0235   }
0236
0237   VECCORE_ATT_HOST_DEVICE
0238   void Set(const T a) { Set(a, a, a, a); }
0239
0240   /// \return the length squared perpendicular to z direction
0241   VECCORE_ATT_HOST_DEVICE
0242   VECGEOM_FORCE_INLINE
0243   T Perp2() const { return fVec[0] * fVec[0] + fVec[1] * fVec[1]; }
0244
0245   /// \return the length perpendicular to z direction
0246   VECCORE_ATT_HOST_DEVICE
0247   VECGEOM_FORCE_INLINE
0248   T Perp() const { return Sqrt(Perp2()); }
0249
0250   /// \set the length perpendicular to z direction
0251   VECCORE_ATT_HOST_DEVICE
0252   VECGEOM_FORCE_INLINE
0253   T SetPerp(T perp) const
0254   {
0255     T factor = perp / NonZero(Perp());
0256     fVec[0] *= factor;
0257     fVec[1] *= factor;
0258   }
0259
0260   /// The dot product of two LorentzVector<T> objects
0261   /// \return T (where T is float, double, or various SIMD vector types)
0262   template <typename U>
0263   VECGEOM_FORCE_INLINE
0264   VECCORE_ATT_HOST_DEVICE
0265   static T Dot(LorentzVector<T> const &left, LorentzVector<U> const &right)
0266   {
0267     return -left[0] * right[0] - left[1] * right[1] - left[2] * right[2] + left[3] * right[3];
0268   }
0269
0270   /// The dot product of two LorentzVector<T> objects
0271   /// \return T (where T is float, double, or various SIMD vector types)
0272   template <typename U>
0273   VECGEOM_FORCE_INLINE
0274   VECCORE_ATT_HOST_DEVICE
0275   T Dot(LorentzVector<U> const &right) const
0276   {
0277     return Dot(*this, right);
0278   }
0279
0280   // For UVector3 compatibility. It is equal to normal multiplication.
0281   VECCORE_ATT_HOST_DEVICE
0282   VECGEOM_FORCE_INLINE
0283   VecType MultiplyByComponents(VecType const &other) const { return *this * other; }
0284
0285   /// \return Squared magnitude of the vector.
0286   VECCORE_ATT_HOST_DEVICE
0287   VECGEOM_FORCE_INLINE
0288   T Mag2() const { return Dot(*this, *this); }
0289
0290   /// \return Magnitude of the vector.
0291   VECCORE_ATT_HOST_DEVICE
0292   VECGEOM_FORCE_INLINE
0293   T Mag() const { return Sqrt(Mag2()); }
0294
0295   /// \return Azimuthal angle between -pi and pi.
0296   VECCORE_ATT_HOST_DEVICE
0297   VECGEOM_FORCE_INLINE
0298   T Phi() const { return ATan2(fVec[1], fVec[0]); }
0299
0300   /// \return Polar angle between 0 and pi.
0301   VECCORE_ATT_HOST_DEVICE
0302   VECGEOM_FORCE_INLINE
0303   T Theta() const { return SpaceVector<T>().Theta(); }
0304
0305   /// Maps each vector entry to a function that manipulates the entry type.
0306   /// \param f A function of type "T f(const T&)" to map over entries.
0307   VECCORE_ATT_HOST_DEVICE
0308   VECGEOM_FORCE_INLINE
0309   void Map(T (*f)(const T &))
0310   {
0311     fVec[0] = f(fVec[0]);
0312     fVec[1] = f(fVec[1]);
0313     fVec[2] = f(fVec[2]);
0314     fVec[3] = f(fVec[3]);
0315   }
0316
0317   VECCORE_ATT_HOST_DEVICE
0318   VECGEOM_FORCE_INLINE
0319   LorentzVector<T> Abs() const
0320   {
0321     using vecCore::math::Abs;
0322     return LorentzVector<T>(Abs(fVec[0]), Abs(fVec[1]), Abs(fVec[2]), Abs(fVec[3]));
0323   }
0324
0325   template <typename U>
0326   VECGEOM_FORCE_INLINE
0327   VECCORE_ATT_HOST_DEVICE
0328   void MaskedAssign(LorentzVector<U> const &condition, LorentzVector<T> const &value)
0329   {
0330     vecCore::MaskedAssign(fVec[0], condition[0], value[0]);
0331     vecCore::MaskedAssign(fVec[1], condition[1], value[1]);
0332     vecCore::MaskedAssign(fVec[2], condition[2], value[2]);
0333     vecCore::MaskedAssign(fVec[3], condition[3], value[3]);
0334   }
0335
0336   VECCORE_ATT_HOST_DEVICE
0337   VECGEOM_FORCE_INLINE
0338   static VecType FromCylindrical(T r, T phi, T z, T t) { return VecType(r * cos(phi), r * sin(phi), z, t); }
0339
0340   VECCORE_ATT_HOST_DEVICE
0341   VECGEOM_FORCE_INLINE
0342   VecType &FixZeroes()
0343   {
0344     using vecCore::math::Abs;
0345     for (int i = 0; i < 4; ++i) {
0346       vecCore::MaskedAssign(fVec[i], Abs(fVec[i]) < kTolerance, T(0.0));
0347     }
0348     return *this;
0349   }
0350
0351 // Inplace binary operators
0352
0353 #define LORENTZVECTOR_TEMPLATE_INPLACE_BINARY_OP(OPERATOR)                                               \
0354   VECCORE_ATT_HOST_DEVICE                                                                                \
0355   VECGEOM_FORCE_INLINE                                                                                   \
0356   VecType &operator OPERATOR(const VecType &other)                                                       \
0357   {                                                                                                      \
0358     fVec[0] OPERATOR other.fVec[0];                                                                      \
0359     fVec[1] OPERATOR other.fVec[1];                                                                      \
0360     fVec[2] OPERATOR other.fVec[2];                                                                      \
0361     fVec[3] OPERATOR other.fVec[3];                                                                      \
0362     return *this;                                                                                        \
0363   }                                                                                                      \
0364   template <typename V>                                                                                  \
0365   VECCORE_ATT_HOST_DEVICE VECGEOM_FORCE_INLINE VecType &operator OPERATOR(const LorentzVector<V> &other) \
0366   {                                                                                                      \
0367     fVec[0] OPERATOR other[0];                                                                           \
0368     fVec[1] OPERATOR other[1];                                                                           \
0369     fVec[2] OPERATOR other[2];                                                                           \
0370     fVec[3] OPERATOR other[3];                                                                           \
0371     return *this;                                                                                        \
0372   }                                                                                                      \
0373   VECCORE_ATT_HOST_DEVICE                                                                                \
0374   VECGEOM_FORCE_INLINE                                                                                   \
0375   VecType &operator OPERATOR(const T &scalar)                                                            \
0376   {                                                                                                      \
0377     fVec[0] OPERATOR scalar;                                                                             \
0378     fVec[1] OPERATOR scalar;                                                                             \
0379     fVec[2] OPERATOR scalar;                                                                             \
0380     fVec[3] OPERATOR scalar;                                                                             \
0381     return *this;                                                                                        \
0382   }
0383   LORENTZVECTOR_TEMPLATE_INPLACE_BINARY_OP(+=)
0384   LORENTZVECTOR_TEMPLATE_INPLACE_BINARY_OP(-=)
0385   LORENTZVECTOR_TEMPLATE_INPLACE_BINARY_OP(*=)
0386   LORENTZVECTOR_TEMPLATE_INPLACE_BINARY_OP(/=)
0387   LORENTZVECTOR_TEMPLATE_INPLACE_BINARY_OP(*)
0388
0389 #undef LORENTZVECTOR_TEMPLATE_INPLACE_BINARY_OP
0390
0391   VECCORE_ATT_HOST_DEVICE
0392   VECGEOM_FORCE_INLINE
0393   operator bool() const { return fVec[0] && fVec[1] && fVec[2] && fVec[3]; }
0394
0395   template <typename U>
0396   VECGEOM_FORCE_INLINE
0397   VECCORE_ATT_HOST_DEVICE
0398   Vector3D<U> SpaceVector() const
0399   {
0400     return Vector3D<U>(fVec[0], fVec[1], fVec[2]);
0401   }
0402
0403   VECGEOM_FORCE_INLINE
0404   VECCORE_ATT_HOST_DEVICE
0405   Vector3D<T> vect() const { return Vector3D<T>(fVec[0], fVec[1], fVec[2]); }
0406
0407   template <typename U>
0408   VECGEOM_FORCE_INLINE
0409   VECCORE_ATT_HOST_DEVICE
0410   VecType Boost(const Vector3D<U> &beta)
0411   {
0412     U beta2 = beta.Mag2();
0413     if (beta2 <= 1e-16) return LorentzVector(*this);
0414     U gamma  = Sqrt(1. / (1 - beta2));
0415     U gamma2 = (gamma - 1.0) / beta2;
0416     U bdotv  = beta.Dot(SpaceVector<U>());
0417     fVec[0] += (gamma2 * bdotv + gamma * fVec[3]) * beta[0];
0418     fVec[1] += (gamma2 * bdotv + gamma * fVec[3]) * beta[1];
0419     fVec[2] += (gamma2 * bdotv + gamma * fVec[3]) * beta[2];
0420     fVec[3] = gamma * (fVec[3] + bdotv);
0421     return *this;
0422   }
0423
0424   Vector3D<T> BoostVector() const
0425   {
0426     if (fVec[3] == 0) {
0427       if (vect().Mag2() < kTolerance) {
0428         return Vector3D<T>(0);
0429       } else {
0430         std::cerr << "LorentzVector::boostVector() - "
0431                   << "boostVector computed for LorentzVector with t=0 -- infinite result" << std::endl;
0432         return Vector3D<T>(InfinityLength<T>());
0433       }
0434     }
0435     if (m2() <= 0) {
0436       std::cerr << "LorentzVector::boostVector() - "
0437                 << "boostVector computed for a non-timelike LorentzVector " << std::endl;
0438       // result will make analytic sense but is physically meaningless
0439     }
0440     return (1. / fVec[3]) * vect();
0441   }
0442
0443   VECCORE_ATT_HOST_DEVICE
0444   VECGEOM_FORCE_INLINE
0445   T Beta2() const { return (SpaceVector<T>() / fVec[3]).Mag2(); }
0446
0447   VECCORE_ATT_HOST_DEVICE
0448   VECGEOM_FORCE_INLINE
0449   T Beta() const { return Sqrt(Beta2()); }
0450
0451   VECCORE_ATT_HOST_DEVICE
0452   VECGEOM_FORCE_INLINE
0453   T Gamma2() const { return 1 / (1 - Beta2()); }
0454
0455   VECCORE_ATT_HOST_DEVICE
0456   VECGEOM_FORCE_INLINE
0457   T Gamma() const { return Sqrt(Gamma2()); }
0458
0459   VECCORE_ATT_HOST_DEVICE
0460   VECGEOM_FORCE_INLINE
0461   T Rapidity() const { return 0.5 * Log((fVec[3] + fVec[2]) / (fVec[3] - fVec[2])); }
0462
0463   VECCORE_ATT_HOST_DEVICE
0464   VECGEOM_FORCE_INLINE
0465   T PseudoRapidity() const { return -Log(Tan(0.5 * Theta())); }
0466
0467   LorentzVector<T> &operator*=(const LorentzRotation<T> &m1);
0468
0469   LorentzVector<T> &transform(const LorentzRotation<T> &m1);
0470
0471   /**
0472    * This method generates a two body decay
0473    @param masses [input] a vector of lenght two with the masses of the products
0474    @param t1     [out]   LorentzVector of the first decay product
0475    @param t2     [out]   LorentzVector of the second decay product
0476    @param func   [in]    function defining the angular distribution of the particles if not uniform
0477   */
0478   template <typename U, typename V, typename W>
0479   VECCORE_ATT_HOST_DEVICE
0480   void TwoBodyDecay(const U masses[], LorentzVector<V> &t1, LorentzVector<W> &t2,
0481                     void (*func)(double &mthe, double &mphi) = 0) const
0482   {
0483     // Make two body decay
0484     double theta = 0;
0485     double phi   = 0;
0486     if (func)
0487       func(theta, phi);
0488     else {
0489       phi   = RNG::Instance().uniform() * kTwoPi;
0490       theta = RNG::Instance().uniform() * kPi;
0491     }
0492     double mass2                          = Mag2();
0493     if (mass2 > -1e-6 && mass2 < 0) mass2 = 0;
0494     double dmass                          = masses[0] + masses[1];
0495     if (mass2 * (1. + 1.e-7) < dmass * dmass || (mass2 == 0 && dmass == 0)) {
0496       if (dmass != 0)
0497         std::cout << "WARNING::LorentzVector::TwoBodyDecay: daughters' mass square " << dmass * dmass
0498                   << " larger than parent " << mass2 << ", no decay " << std::endl;
0499       t1.Set(0, 0, 0, masses[0]);
0500       t2.Set(0, 0, 0, masses[1]);
0501       return;
0502     }
0503     const double mass = sqrt(mass2);
0504
0505     double en1 = 0.5 * (mass * mass + masses[0] * masses[0] - masses[1] * masses[1]) / mass;
0506     double en2 = 0.5 * (mass * mass - masses[0] * masses[0] + masses[1] * masses[1]) / mass;
0507     if ((en1 - masses[0] < -1e-10) || (en2 - masses[1] < -1e-10)) {
0508       std::cout << "WARNING::LorentzVector::TwoBodyDecay en<0  en1-masses[0] " << en1 - masses[0] << " en2-masses[1] "
0509                 << en2 - masses[1] << std::endl;
0510     }
0511     double p1 = 0;
0512     if (en1 < masses[0])
0513       en1 = masses[0];
0514     else
0515       p1      = sqrt((en1 + masses[0]) * (en1 - masses[0]));
0516     double p2 = 0;
0517     if (en2 < masses[1])
0518       en2 = masses[1];
0519     else
0520       p2     = sqrt((en2 + masses[1]) * (en2 - masses[1]));
0521     double p = 0.5 * (p1 + p2);
0522     // debug
0523     if (fabs(p1 - p2) / mass > 1e-4)
0524       std::cout << "WARNING::LorentzVector::TwoBodyDecay: momentum imbalance err " << fabs(p1 - p2) / mass << " p1-p2 "
0525                 << p1 - p2 << " p1 " << p1 << " p2 " << p2 << " mass " << mass << std::endl;
0526     if (fabs(mass - en1 - en2) / (mass + en1 + en2) / 3 > 1e-7)
0527       std::cout << "WARNING::LorentzVector::TwoBodyDecay: mass imbalance mass - en1 - en2 " << mass - en1 - en2
0528                 << std::endl;
0529     // debug
0530
0531     // Set the two vectors
0532     const Vector3D<T> space(p1 * sin(theta) * cos(phi), p1 * sin(theta) * sin(phi), p1 * cos(theta));
0533     t1.Set(space, en1);
0534     t2.Set(-space, en2);
0535
0536     if (func) {
0537       // This part only if there is a matrix element, otherwise isotropic is isotropic
0538       const Transformation3D tr(SpaceVector<T>(), false);
0539       tr.Transform(space, static_cast<Vector3D<T> &>(t1));
0540       tr.Transform(-space, static_cast<Vector3D<T> &>(t2));
0541     }
0542
0543     Vector3D<T> beta(-fVec[0] / fVec[3], -fVec[1] / fVec[3], -fVec[2] / fVec[3]);
0544     t1 = t1.Boost(beta);
0545     t2 = t2.Boost(beta);
0546     // dedbug
0547     static const LorentzVector<T> metric(-1, -1, -1, 1);
0548     const LorentzVector<T> vdiff = (*this) - t1 - t2;
0549     const double mdiff =
0550         sqrt(vdiff.Dot(vdiff * metric) / (Dot((*this) * metric) + t1.Dot(t1 * metric) + t2.Dot(t2 * metric))) / 3;
0551     if (mdiff > 1e-7)
0552       std::cout << "WARNING::LorentzVector::TwoBodyDecay error " << mdiff
0553                 << " in momentum energy balance this=" << (*this) << " difference=" << (*this) - t1 - t2 << std::endl;
0554     // debug
0555   }
0556
0557   /**
0558    * This method generate recursively the phase space via two body decays
0559    * The method is robust, however it provides a biased distribution
0560    * this will be fixed in due time, for the moment it is good enough
0561    @param[in] masses the masses of the products
0562    @param[out] daughters a vector of LorentzVectors with the decay products
0563   */
0564   template <typename U>
0565   VECCORE_ATT_HOST_DEVICE
0566   void PhaseSpace(const Vector<U> &masses, Vector<LorentzVector<T>> &daughters)
0567   {
0568     // Relativistic phase space
0569     LorentzVector<T> pa(*this);
0570     LorentzVector<T> d1;
0571     LorentzVector<T> d2;
0572
0573     //
0574     int ndec = masses.size();
0575     daughters.clear();
0576     double mass2 = 0;
0577     for (int j = 0; j < ndec; ++j)
0578       mass2 += masses[j];
0579     if (mass2 > Mag()) {
0580       std::cout << "LorentzVector::PhaseSpace: cannot decay: parent mass " << Mag() << " sum of daughter mass " << mass2
0581                 << std::endl;
0582       return;
0583     }
0584     for (int j = 0; j < ndec - 1; ++j) {
0585       mass2 = masses[j + 1];
0586       if (j < ndec - 2) {
0587         for (int i = j + 2; i < ndec; ++i)
0588           mass2 += masses[i];
0589         double dm          = 0;
0590         double pam2        = pa.Mag2();
0591         if (pam2 < 0) pam2 = 0;
0592         double freen       = sqrt(pam2) - masses[j] - mass2;
0593         if (mass2 > 0)
0594           dm = freen *
0595                std::pow(RNG::Instance().uniform(), masses[j] / (mass2 - 0.5 * (masses[ndec - 1] + masses[ndec - 2])));
0596         while (dm < 1e-7 * freen)
0597           dm = freen * RNG::Instance().uniform();
0598         mass2 += dm;
0599       }
0600       double vmass[2] = {masses[j], mass2};
0601       pa.TwoBodyDecay(vmass, d1, d2);
0602       pa = d2;
0603       daughters.push_back(d1);
0604     }
0605     daughters.push_back(d2);
0606   }
0607 };
0608
0609 template <typename T>
0610 std::ostream &operator<<(std::ostream &os, LorentzVector<T> const &vec)
0611 {
0612   os << "(" << vec[0] << ", " << vec[1] << ", " << vec[2] << ", " << vec[3] << ")";
0613   return os;
0614 }
0615
0616 #define LORENTZVECTOR_BINARY_OP(OPERATOR, INPLACE)                                                             \
0617   template <typename T, typename V>                                                                            \
0618   VECGEOM_FORCE_INLINE VECCORE_ATT_HOST_DEVICE LorentzVector<T> operator OPERATOR(const LorentzVector<T> &lhs, \
0619                                                                                   const LorentzVector<V> &rhs) \
0620   {                                                                                                            \
0621     LorentzVector<T> result(lhs);                                                                              \
0622     result INPLACE rhs;                                                                                        \
0623     return result;                                                                                             \
0624   }                                                                                                            \
0625   template <typename T, typename V>                                                                            \
0626   VECGEOM_FORCE_INLINE VECCORE_ATT_HOST_DEVICE LorentzVector<T> operator OPERATOR(LorentzVector<T> const &lhs, \
0627                                                                                   const V rhs)                 \
0628   {                                                                                                            \
0629     LorentzVector<T> result(lhs);                                                                              \
0630     result INPLACE rhs;                                                                                        \
0631     return result;                                                                                             \
0632   }                                                                                                            \
0633   template <typename T, typename V>                                                                            \
0634   VECGEOM_FORCE_INLINE VECCORE_ATT_HOST_DEVICE LorentzVector<T> operator OPERATOR(const V lhs,                 \
0635                                                                                   LorentzVector<T> const &rhs) \
0636   {                                                                                                            \
0637     LorentzVector<T> result(lhs);                                                                              \
0638     result INPLACE rhs;                                                                                        \
0639     return result;                                                                                             \
0640   }
0641 LORENTZVECTOR_BINARY_OP(+, +=)
0642 LORENTZVECTOR_BINARY_OP(-, -=)
0643 LORENTZVECTOR_BINARY_OP(*, *=)
0644 LORENTZVECTOR_BINARY_OP(/, /=)
0645 #undef LORENTZVECTOR_BINARY_OP
0646
0647 VECGEOM_FORCE_INLINE
0648 VECCORE_ATT_HOST_DEVICE
0649 bool operator==(LorentzVector<Precision> const &lhs, LorentzVector<Precision> const &rhs)
0650 {
0651   return Abs(lhs[0] - rhs[0]) < kTolerance && Abs(lhs[1] - rhs[1]) < kTolerance && Abs(lhs[2] - rhs[2]) < kTolerance &&
0652          Abs(lhs[3] - rhs[3]) < kTolerance;
0653 }
0654
0655 VECGEOM_FORCE_INLINE
0656 VECCORE_ATT_HOST_DEVICE
0657 LorentzVector<bool> operator!=(LorentzVector<Precision> const &lhs, LorentzVector<Precision> const &rhs)
0658 {
0659   return !(lhs == rhs);
0660 }
0661
0662 template <typename T>
0663 VECGEOM_FORCE_INLINE
0664 VECCORE_ATT_HOST_DEVICE
0665 LorentzVector<T> operator-(LorentzVector<T> const &vec)
0666 {
0667   return LorentzVector<T>(-vec[0], -vec[1], -vec[2], -vec[3]);
0668 }
0669
0670 VECCORE_ATT_HOST_DEVICE
0671 VECGEOM_FORCE_INLINE
0672 LorentzVector<bool> operator!(LorentzVector<bool> const &vec)
0673 {
0674   return LorentzVector<bool>(!vec[0], !vec[1], !vec[2], !vec[3]);
0675 }
0676
0677 #pragma GCC diagnostic push
0678 #pragma GCC diagnostic ignored "-Weffc++"
0679 #define LORENTZVECTOR_SCALAR_BOOLEAN_LOGICAL_OP(OPERATOR)                                               \
0680   VECCORE_ATT_HOST_DEVICE                                                                               \
0681   VECGEOM_FORCE_INLINE                                                                                  \
0682   LorentzVector<bool> operator OPERATOR(LorentzVector<bool> const &lhs, LorentzVector<bool> const &rhs) \
0683   {                                                                                                     \
0684     return LorentzVector<bool>(lhs[0] OPERATOR rhs[0], lhs[1] OPERATOR rhs[1], lhs[2] OPERATOR rhs[2],  \
0685                                lhs[3] OPERATOR rhs[3]);                                                 \
0686   }
0687 LORENTZVECTOR_SCALAR_BOOLEAN_LOGICAL_OP(&&)
0688 LORENTZVECTOR_SCALAR_BOOLEAN_LOGICAL_OP(||)
0689 #undef LORENTZVECTOR_SCALAR_BOOLEAN_LOGICAL_OP
0690 #pragma GCC diagnostic pop
0691
0692 static const LorentzVector<double> LorentzMetric(-1, -1, -1, 1);
0693
0694 } // End inline namespace
0695 } // End global namespace
0696
0697 #endif // VECGEOM_BASE_LORENTZVECTOR_H_