Rivet/Math/LorentzTrans.hh

0001 #ifndef RIVET_MATH_LORENTZTRANS
0002 #define RIVET_MATH_LORENTZTRANS
0003
0004 #include "Rivet/Math/MathConstants.hh"
0005 #include "Rivet/Math/MathUtils.hh"
0006 #include "Rivet/Math/MatrixN.hh"
0007 #include "Rivet/Math/Matrix3.hh"
0008 #include "Rivet/Math/Vector4.hh"
0009 #include <ostream>
0010
0011 namespace Rivet {
0012
0013
0014   /// @brief Object implementing Lorentz transform calculations and boosts.
0015   ///
0016   /// @note These boosts are defined actively, i.e. as modifications of vectors
0017   /// rather than frame transformations. So the boost vector is the opposite of
0018   /// what you might expect.
0019   ///
0020   /// @todo Review the active/passive convention choice. Seems counterintuitive now...
0021   class LorentzTransform {
0022   public:
0023
0024     /// @name Simple Lorentz factor conversions
0025     /// @{
0026
0027     /// Calculate the \f$ \gamma \f$ factor from \f$ \beta \f$
0028     static double beta2gamma(double beta) {
0029       return 1.0 / sqrt(1 - sqr(beta));
0030     }
0031
0032     /// Calculate \f$ \beta \f$ from the \f$ \gamma \f$ factor
0033     static double gamma2beta(double gamma) {
0034       return sqrt(1 - sqr(1/gamma));
0035     }
0036
0037     /// @}
0038
0039
0040     /// @name Construction
0041     /// @{
0042
0043     /// Default (identity) constructor
0044     LorentzTransform() {
0045       _boostMatrix = Matrix<4>::mkIdentity();
0046     }
0047
0048     /// Constructor from a 4x4 matrix
0049     LorentzTransform(const Matrix<4>& boostMatrix) {
0050       _boostMatrix = boostMatrix;
0051     }
0052
0053     /// Make an LT for an active boost (i.e. object velocity += in boost direction)
0054     static LorentzTransform mkObjTransformFromBeta(const Vector3& vbeta) {
0055       LorentzTransform rtn;
0056       return rtn.setBetaVec(vbeta);
0057     }
0058
0059     /// Make an LT for a passive boost (i.e. object velocity -= in boost direction)
0060     static LorentzTransform mkFrameTransformFromBeta(const Vector3& vbeta) {
0061       LorentzTransform rtn;
0062       return rtn.setBetaVec(-vbeta);
0063     }
0064
0065     /// Make an LT for an active boost (i.e. object velocity += in boost direction)
0066     static LorentzTransform mkObjTransformFromGamma(const Vector3& vgamma) {
0067       LorentzTransform rtn;
0068       if (!vgamma.isZero()) rtn.setGammaVec(vgamma);
0069       return rtn;
0070     }
0071
0072     /// Make an LT for a passive boost (i.e. object velocity -= in boost direction)
0073     static LorentzTransform mkFrameTransformFromGamma(const Vector3& vgamma) {
0074       LorentzTransform rtn;
0075       if (!vgamma.isZero()) rtn.setGammaVec(-vgamma);
0076       return rtn;
0077     }
0078
0079     /// Make an LT for an active boost (i.e. object velocity += in boost direction)
0080     static LorentzTransform mkObjTransform(const FourMomentum& p4) {
0081       return mkObjTransformFromBeta(p4.betaVec());
0082     }
0083
0084     /// Make an LT for a passive boost (i.e. object velocity -= in boost direction)
0085     static LorentzTransform mkFrameTransform(const FourMomentum& p4) {
0086       return mkObjTransformFromBeta(-p4.betaVec());
0087     }
0088
0089     /// @}
0090
0091
0092     /// @name Boost vector and beta/gamma factors
0093     /// @{
0094
0095     /// Set up an active Lorentz boost from a (unit) direction vector, and \f$ \beta \f$ & \f$ \gamma \f$ factors
0096     LorentzTransform& _setBoost(const Vector3& vec, double beta, double gamma) {
0097       // Set to identity for null boosts
0098       _boostMatrix = Matrix<4>::mkIdentity();
0099       if (isZero(beta)) return *this;
0100       //
0101       // It's along the x, y, or z axis if 2 Cartesian components are zero
0102       const bool alongxyz = (int(vec.x() == 0) + int(vec.y() == 0) + int(vec.z() == 0) == 2);
0103       const int i = (!alongxyz || vec.x() != 0) ? 1 : (vec.y() != 0) ? 2 : 3;
0104       const int isign = !alongxyz ? 1 : sign(vec[i-1]);
0105       //
0106       _boostMatrix.set(0, 0, gamma);
0107       _boostMatrix.set(i, i, gamma);
0108       _boostMatrix.set(0, i, +isign*beta*gamma); //< +ve coeff since active boost
0109       _boostMatrix.set(i, 0, +isign*beta*gamma); //< +ve coeff since active boost
0110       //
0111       if (!alongxyz) _boostMatrix = rotate(Vector3::mkX(), vec)._boostMatrix;
0112       return *this;
0113     }
0114
0115     /// Set up an active Lorentz boost from the \f$ \vec\beta \f$ vector
0116     LorentzTransform& setBetaVec(const Vector3& vbeta) {
0117       // Set to identity for null boosts
0118       _boostMatrix = Matrix<4>::mkIdentity();
0119       if (isZero(vbeta.mod2())) return *this;
0120       const double beta = vbeta.mod();
0121       const double gamma = beta2gamma(beta);
0122       return _setBoost(vbeta.unit(), beta, gamma);
0123     }
0124
0125     /// Get the \f$ \vec\beta \f$ vector for an active Lorentz boost
0126     Vector3 betaVec() const {
0127       FourMomentum boost(_boostMatrix.getColumn(0)); //< @todo WRONG?!
0128       //cout << "!!!" << boost << '\n';
0129       if (boost.isZero()) return Vector3();
0130       assert(boost.E() > 0);
0131       const double beta = boost.p3().mod() / boost.E();
0132       return boost.p3().unit() * beta;
0133     }
0134
0135     /// Get the \f$ \beta \f$ factor
0136     double beta() const {
0137       return betaVec().mod();
0138     }
0139
0140
0141     /// Set up an active Lorentz boost from the \f$ \vec\gamma \f$ vector
0142     LorentzTransform& setGammaVec(const Vector3& vgamma) {
0143       // Set to identity for null boosts
0144       _boostMatrix = Matrix<4>::mkIdentity();
0145       if (isZero(vgamma.mod2() - 1)) return *this;
0146       const double gamma = vgamma.mod();
0147       const double beta = gamma2beta(gamma);
0148       return _setBoost(vgamma.unit(), beta, gamma);
0149     }
0150
0151     /// Get the \f$ \vec\gamma \f$ vector for an active Lorentz boost
0152     Vector3 gammaVec() const {
0153       FourMomentum boost(_boostMatrix.getColumn(0)); //< @todo WRONG?!
0154       if (boost.isZero()) return Vector3();
0155       assert(boost.E() > 0);
0156       const double beta = boost.p3().mod() / boost.E();
0157       return boost.p3().unit() * beta;
0158     }
0159
0160     /// Get the \f$ \gamma \f$ factor
0161     double gamma() const {
0162       return beta2gamma(beta());
0163     }
0164
0165     /// @}
0166
0167
0168     /// Apply this transformation to the given 4-vector
0169     FourVector transform(const FourVector& v4) const {
0170       return multiply(_boostMatrix, v4);
0171     }
0172
0173     /// Apply this transformation to the given 4-mometum
0174     FourMomentum transform(const FourMomentum& v4) const {
0175       return multiply(_boostMatrix, v4);
0176     }
0177
0178     /// Apply this transformation to the given 4-vector
0179     FourVector operator () (const FourVector& v4) const {
0180       return transform(v4);
0181     }
0182
0183     /// Apply this transformation to the given 4-mometum
0184     FourMomentum operator () (const FourMomentum& v4) const {
0185       return transform(v4);
0186     }
0187
0188
0189     /// @name Transform modifications
0190     /// @{
0191
0192     /// Rotate the transformation cf. the difference between vectors @a from and @a to
0193     LorentzTransform rotate(const Vector3& from, const Vector3& to) const {
0194       return rotate(Matrix3(from, to));
0195     }
0196
0197     /// Rotate the transformation by @a angle radians about @a axis
0198     LorentzTransform rotate(const Vector3& axis, double angle) const {
0199       return rotate(Matrix3(axis, angle));
0200     }
0201
0202     /// Rotate the transformation by the 3D rotation matrix @a rot
0203     LorentzTransform rotate(const Matrix3& rot) const {
0204       LorentzTransform lt = *this;
0205       const Matrix4 rot4 = _mkMatrix4(rot);
0206       const Matrix4 newlt = rot4 * _boostMatrix * rot4.inverse();
0207       lt._boostMatrix = newlt;
0208       return lt;
0209     }
0210
0211     /// Calculate the inverse transform
0212     LorentzTransform inverse() const {
0213       LorentzTransform rtn;
0214       rtn._boostMatrix = _boostMatrix.inverse();
0215       return rtn;
0216     }
0217
0218     /// Combine LTs, treating @a this as the LH matrix.
0219     LorentzTransform combine(const LorentzTransform& lt) const {
0220       LorentzTransform rtn;
0221       rtn._boostMatrix = _boostMatrix * lt._boostMatrix;
0222       return rtn;
0223     }
0224
0225     /// Operator combination of two LTs
0226     LorentzTransform operator*(const LorentzTransform& lt) const {
0227       return combine(lt);
0228     }
0229
0230     /// Pre-multiply m3 by this LT
0231     LorentzTransform preMult(const Matrix3& m3) {
0232       _boostMatrix = multiply(_mkMatrix4(m3),_boostMatrix);
0233       return *this;
0234     }
0235
0236     /// Post-multiply m3 by this LT
0237     LorentzTransform postMult(const Matrix3& m3) {
0238       _boostMatrix *= _mkMatrix4(m3);
0239       return *this;
0240     }
0241
0242     /// @}
0243
0244
0245     /// Return the matrix form
0246     Matrix4 toMatrix() const {
0247       return _boostMatrix;
0248     }
0249
0250
0251   private:
0252
0253     Matrix4 _mkMatrix4(const Matrix3& m3) const {
0254       Matrix4 m4 = Matrix4::mkIdentity();
0255       for (size_t i = 0; i < 3; ++i) {
0256         for (size_t j = 0; j < 3; ++j) {
0257           m4.set(i+1, j+1, m3.get(i, j));
0258         }
0259       }
0260       return m4;
0261     }
0262
0263     Matrix4 _boostMatrix;
0264
0265   };
0266
0267
0268
0269   inline LorentzTransform inverse(const LorentzTransform& lt) {
0270     return lt.inverse();
0271   }
0272
0273   inline LorentzTransform combine(const LorentzTransform& a, const LorentzTransform& b) {
0274     return a.combine(b);
0275   }
0276
0277   inline FourVector transform(const LorentzTransform& lt, const FourVector& v4) {
0278       return lt.transform(v4);
0279   }
0280
0281
0282   //////////////////////////
0283
0284
0285   inline string toString(const LorentzTransform& lt) {
0286     return toString(lt.toMatrix());
0287   }
0288
0289   inline std::ostream& operator<<(std::ostream& out, const LorentzTransform& lt) {
0290     out << toString(lt);
0291     return out;
0292   }
0293
0294
0295 }
0296
0297 #endif