VecGeom/base/LorentzRotation.h

0001 /// \file LorentzVector.h
0002 /// \author Andrei Gheata (based on CLHEP/Vector/LorentzRotation.h)
0003
0004 #ifndef VECGEOM_BASE_LORENTZROTATION_H_
0005 #define VECGEOM_BASE_LORENTZROTATION_H_
0006
0007 #include "VecGeom/base/AlignedBase.h"
0008 #include "VecGeom/base/LorentzVector.h"
0009
0010 namespace vecgeom {
0011
0012 VECGEOM_DEVICE_FORWARD_DECLARE(template <typename T> class LorentzRotation;);
0013
0014 inline namespace VECGEOM_IMPL_NAMESPACE {
0015
0016 /**
0017  * @brief Lorentz rotation class for performing rotations and boosts on Lorentz vectors
0018  * @details If vector acceleration is enabled, the scalar template instantiation
0019  *          will use vector instructions for operations when possible.
0020  */
0021 template <typename T>
0022 class LorentzRotation : public AlignedBase {
0023
0024   typedef LorentzVector<T> VecType;
0025
0026 private:
0027   T fxx, fxy, fxz, fxt, fyx, fyy, fyz, fyt, fzx, fzy, fzz, fzt, ftx, fty, ftz, ftt; // The matrix elements.
0028
0029 public:
0030   VECCORE_ATT_HOST_DEVICE
0031   VECGEOM_FORCE_INLINE
0032   LorentzRotation()
0033       : fxx(1.0), fxy(0.0), fxz(0.0), fxt(0.0), fyx(0.0), fyy(1.0), fyz(0.0), fyt(0.0), fzx(0.0), fzy(0.0), fzz(1.0),
0034         fzt(0.0), ftx(0.0), fty(0.0), ftz(0.0), ftt(1.0)
0035   {
0036   }
0037
0038   template <typename U>
0039   VECGEOM_FORCE_INLINE
0040   VECCORE_ATT_HOST_DEVICE
0041   LorentzRotation(LorentzRotation<U> const &r)
0042       : fxx(r.fxx), fxy(r.fxy), fxz(r.fxz), fxt(r.fxt), fyx(r.fyx), fyy(r.fyy), fyz(r.fyz), fyt(r.fyt), fzx(r.fzx),
0043         fzy(r.fzy), fzz(r.fzz), fzt(r.fzt), ftx(r.ftx), fty(r.fty), ftz(r.ftz), ftt(r.ftt)
0044   {
0045   }
0046
0047   VECCORE_ATT_HOST_DEVICE
0048   VECGEOM_FORCE_INLINE
0049   LorentzRotation &operator=(LorentzRotation const &r)
0050   {
0051     if (this != &r) {
0052       fxx = r.fxx;
0053       fxy = r.fxy;
0054       fxz = r.fxz;
0055       fxt = r.mxt;
0056       fyx = r.fyx;
0057       fyy = r.fyy;
0058       fyz = r.fyz;
0059       fyt = r.myt;
0060       fzx = r.fzx;
0061       fzy = r.fzy;
0062       fzz = r.fzz;
0063       fzt = r.mzt;
0064       ftx = r.ftx;
0065       fty = r.fty;
0066       ftz = r.ftz;
0067       ftt = r.mtt;
0068     }
0069     return *this;
0070   }
0071
0072   LorentzVector<T> vectorMultiplication(const LorentzVector<T> &p) const
0073   {
0074     T x(p.x());
0075     T y(p.y());
0076     T z(p.z());
0077     T t(p.t());
0078     return LorentzVector<T>(fxx * x + fxy * y + fxz * z + fxt * t, fyx * x + fyy * y + fyz * z + fyt * t,
0079                             fzx * x + fzy * y + fzz * z + fzt * t, ftx * x + fty * y + ftz * z + ftt * t);
0080   }
0081
0082   VECCORE_ATT_HOST_DEVICE
0083   VECGEOM_FORCE_INLINE
0084   LorentzVector<T> col1() const { return LorentzVector<T>(fxx, fyx, fzx, ftx); }
0085
0086   VECCORE_ATT_HOST_DEVICE
0087   VECGEOM_FORCE_INLINE
0088   LorentzVector<T> col2() const { return LorentzVector<T>(fxy, fyy, fzy, fty); }
0089
0090   VECCORE_ATT_HOST_DEVICE
0091   VECGEOM_FORCE_INLINE
0092   LorentzVector<T> col3() const { return LorentzVector<T>(fxz, fyz, fzz, ftz); }
0093
0094   VECCORE_ATT_HOST_DEVICE
0095   VECGEOM_FORCE_INLINE
0096   LorentzVector<T> col4() const { return LorentzVector<T>(fxt, fyt, fzt, ftt); }
0097
0098   VECCORE_ATT_HOST_DEVICE
0099   VECGEOM_FORCE_INLINE
0100   LorentzVector<T> row1() const { return LorentzVector<T>(fxx, fxy, fxz, fxt); }
0101
0102   VECCORE_ATT_HOST_DEVICE
0103   VECGEOM_FORCE_INLINE
0104   LorentzVector<T> row2() const { return LorentzVector<T>(fyx, fyy, fyz, fyt); }
0105
0106   VECCORE_ATT_HOST_DEVICE
0107   VECGEOM_FORCE_INLINE
0108   LorentzVector<T> row3() const { return LorentzVector<T>(fzx, fzy, fzz, fzt); }
0109
0110   VECCORE_ATT_HOST_DEVICE
0111   VECGEOM_FORCE_INLINE
0112   LorentzVector<T> row4() const { return LorentzVector<T>(ftx, fty, ftz, ftt); }
0113
0114   VECCORE_ATT_HOST_DEVICE
0115   LorentzRotation<T> &rotateX(T delta)
0116   {
0117     T c1                  = Cos(delta);
0118     T s1                  = Sin(delta);
0119     LorentzVector<T> rowy = row2();
0120     LorentzVector<T> rowz = row3();
0121     LorentzVector<T> r2   = c1 * rowy - s1 * rowz;
0122     LorentzVector<T> r3   = s1 * rowy + c1 * rowz;
0123     fyx                   = r2.x();
0124     fyy                   = r2.y();
0125     fyz                   = r2.z();
0126     fyt                   = r2.t();
0127     fzx                   = r3.x();
0128     fzy                   = r3.y();
0129     fzz                   = r3.z();
0130     fzt                   = r3.t();
0131     return *this;
0132   }
0133
0134   LorentzRotation<T> &rotateY(T delta)
0135   {
0136     T c1                  = std::cos(delta);
0137     T s1                  = std::sin(delta);
0138     LorentzVector<T> rowx = row1();
0139     LorentzVector<T> rowz = row3();
0140     LorentzVector<T> r1   = c1 * rowx + s1 * rowz;
0141     LorentzVector<T> r3   = -s1 * rowx + c1 * rowz;
0142     fxx                   = r1.x();
0143     fxy                   = r1.y();
0144     fxz                   = r1.z();
0145     fxt                   = r1.t();
0146     fzx                   = r3.x();
0147     fzy                   = r3.y();
0148     fzz                   = r3.z();
0149     fzt                   = r3.t();
0150     return *this;
0151   }
0152
0153   LorentzRotation<T> &rotateZ(T delta)
0154   {
0155     T c1                  = std::cos(delta);
0156     T s1                  = std::sin(delta);
0157     LorentzVector<T> rowx = row1();
0158     LorentzVector<T> rowy = row2();
0159     LorentzVector<T> r1   = c1 * rowx - s1 * rowy;
0160     LorentzVector<T> r2   = s1 * rowx + c1 * rowy;
0161     fxx                   = r1.x();
0162     fxy                   = r1.y();
0163     fxz                   = r1.z();
0164     fxt                   = r1.t();
0165     fyx                   = r2.x();
0166     fyy                   = r2.y();
0167     fyz                   = r2.z();
0168     fyt                   = r2.t();
0169     return *this;
0170   }
0171
0172   LorentzRotation<T> &boostX(T beta)
0173   {
0174     T b2                  = beta * beta;
0175     T g1                  = 1.0 / Sqrt(1.0 - b2);
0176     T bg                  = beta * g1;
0177     LorentzVector<T> rowx = row1();
0178     LorentzVector<T> rowt = row4();
0179     LorentzVector<T> r1   = g1 * rowx + bg * rowt;
0180     LorentzVector<T> r4   = bg * rowx + g1 * rowt;
0181     fxx                   = r1.x();
0182     fxy                   = r1.y();
0183     fxz                   = r1.z();
0184     fxt                   = r1.t();
0185     ftx                   = r4.x();
0186     fty                   = r4.y();
0187     ftz                   = r4.z();
0188     ftt                   = r4.t();
0189     return *this;
0190   }
0191
0192   LorentzRotation<T> &boostY(T beta)
0193   {
0194     T b2                  = beta * beta;
0195     T g1                  = 1.0 / std::sqrt(1.0 - b2);
0196     T bg                  = beta * g1;
0197     LorentzVector<T> rowy = row2();
0198     LorentzVector<T> rowt = row4();
0199     LorentzVector<T> r2   = g1 * rowy + bg * rowt;
0200     LorentzVector<T> r4   = bg * rowy + g1 * rowt;
0201     fyx                   = r2.x();
0202     fyy                   = r2.y();
0203     fyz                   = r2.z();
0204     fyt                   = r2.t();
0205     ftx                   = r4.x();
0206     fty                   = r4.y();
0207     ftz                   = r4.z();
0208     ftt                   = r4.t();
0209     return *this;
0210   }
0211
0212   LorentzRotation<T> &boostZ(T beta)
0213   {
0214     T b2                  = beta * beta;
0215     T g1                  = 1.0 / std::sqrt(1.0 - b2);
0216     T bg                  = beta * g1;
0217     LorentzVector<T> rowz = row3();
0218     LorentzVector<T> rowt = row4();
0219     LorentzVector<T> r3   = g1 * rowz + bg * rowt;
0220     LorentzVector<T> r4   = bg * rowz + g1 * rowt;
0221     ftx                   = r4.x();
0222     fty                   = r4.y();
0223     ftz                   = r4.z();
0224     ftt                   = r4.t();
0225     fzx                   = r3.x();
0226     fzy                   = r3.y();
0227     fzz                   = r3.z();
0228     fzt                   = r3.t();
0229     return *this;
0230   }
0231
0232   /** @brief Matrix multiplication. Note a *= b; <=> a = a * b; while a.transform(b); <=> a = b * a; */
0233   // LorentzRotation & operator *= (const LorentzRotation & r);
0234   // LorentzRotation & transform   (const LorentzRotation & r);
0235
0236 protected:
0237   VECCORE_ATT_HOST_DEVICE
0238   VECGEOM_FORCE_INLINE
0239   LorentzRotation(T rxx, T rxy, T rxz, T rxt, T ryx, T ryy, T ryz, T ryt, T rzx, T rzy, T rzz, T rzt, T rtx, T rty,
0240                   T rtz, T rtt)
0241       : fxx(rxx), fxy(rxy), fxz(rxz), fxt(rxt), fyx(ryx), fyy(ryy), fyz(ryz), fyt(ryt), fzx(rzx), fzy(rzy), fzz(rzz),
0242         fzt(rzt), ftx(rtx), fty(rty), ftz(rtz), ftt(rtt)
0243   {
0244   }
0245 };
0246
0247 } // End inline namespace
0248 } // End global namespace
0249
0250 #endif // VECGEOM_BASE_LORENTZVECTOR_H_