Rivet is hosted by Hepforge, IPPP Durham
LorentzTrans.hh
Go to the documentation of this file.
00001 #ifndef RIVET_MATH_LORENTZTRANS
00002 #define RIVET_MATH_LORENTZTRANS
00003 
00004 #include "Rivet/Math/MathHeader.hh"
00005 #include "Rivet/Math/MathUtils.hh"
00006 #include "Rivet/Math/MatrixN.hh"
00007 #include "Rivet/Math/Matrix3.hh"
00008 #include "Rivet/Math/Vector4.hh"
00009 #include <iostream>
00010 
00011 namespace Rivet {
00012 
00013 
00014   /// @brief Object implementing Lorentz transform calculations and boosts.
00015   ///
00016   /// @note These boosts are defined actively, i.e. as modifications of vectors
00017   /// rather than frame transformations. So the boost vector is the opposite of
00018   /// what you might expect.
00019   ///
00020   /// @todo Review the active/passive convention choice. Seems counterintuitive now...
00021   class LorentzTransform {
00022   public:
00023 
00024     /// @name Simple Lorentz factor conversions
00025     //@{
00026 
00027     /// Calculate the \f$ \gamma \f$ factor from \f$ \beta \f$
00028     static double beta2gamma(double beta) {
00029       return 1.0 / sqrt(1 - sqr(beta));
00030     }
00031 
00032     /// Calculate \f$ \beta \f$ from the \f$ \gamma \f$ factor
00033     static double gamma2beta(double gamma) {
00034       return sqrt(1 - sqr(1/gamma));
00035     }
00036 
00037     // /// Calculate the \f$ \gamma \f$ factor from \f$ \bar{\beta}^2 = 1-\beta^2 \f$
00038     // static double betabarsq2gamma(double betabarsq) {
00039     //   return 1.0 / sqrt(betabarsq);
00040     // }
00041 
00042     // /// Calculate \f$ \bar{\beta}^2 = 1-\beta^2 \f$ from the \f$ \gamma \f$ factor
00043     // static double gamma2betabarsq(double gamma) {
00044     //   return 1.0 / sqr(gamma);
00045     // }
00046 
00047     //@}
00048 
00049 
00050     /// @name Construction
00051     //@{
00052 
00053     /// Default (identity) constructor
00054     LorentzTransform() {
00055       _boostMatrix = Matrix<4>::mkIdentity();
00056     }
00057 
00058 
00059     /// Make an LT for an active boost (i.e. object velocity += in boost direction)
00060     static LorentzTransform mkObjTransformFromBeta(const Vector3& vbeta) {
00061       LorentzTransform rtn;
00062       return rtn.setBetaVec(vbeta);
00063     }
00064 
00065     /// Make an LT for a passive boost (i.e. object velocity -= in boost direction)
00066     static LorentzTransform mkFrameTransformFromBeta(const Vector3& vbeta) {
00067       LorentzTransform rtn;
00068       return rtn.setBetaVec(-vbeta);
00069     }
00070 
00071     /// Make an LT for an active boost (i.e. object velocity += in boost direction)
00072     static LorentzTransform mkObjTransformFromGamma(const Vector3& vgamma) {
00073       LorentzTransform rtn;
00074       if (!vgamma.isZero()) rtn.setGammaVec(vgamma);
00075       return rtn;
00076     }
00077 
00078     /// Make an LT for a passive boost (i.e. object velocity -= in boost direction)
00079     static LorentzTransform mkFrameTransformFromGamma(const Vector3& vgamma) {
00080       LorentzTransform rtn;
00081       if (!vgamma.isZero()) rtn.setGammaVec(-vgamma);
00082       return rtn;
00083     }
00084 
00085     //@}
00086 
00087 
00088     /// @name Boost vector and beta/gamma factors
00089     //@{
00090 
00091     /// Set up an active Lorentz boost from the \f$ \vec\beta \f$ vector
00092     LorentzTransform& setBetaVec(const Vector3& vbeta) {
00093       assert(vbeta.mod2() < 1);
00094       const double beta = vbeta.mod();
00095       const double gamma = beta2gamma(beta);
00096       _boostMatrix = Matrix<4>::mkIdentity();
00097       _boostMatrix.set(0, 0, gamma);
00098       _boostMatrix.set(1, 1, gamma);
00099       _boostMatrix.set(0, 1, +beta*gamma); //< +ve coeff since active boost
00100       _boostMatrix.set(1, 0, +beta*gamma); //< +ve coeff since active boost
00101       if (beta > 0) _boostMatrix = rotate(Vector3::mkX(), vbeta)._boostMatrix;
00102       return *this;
00103     }
00104 
00105     /// Get the \f$ \vec\beta \f$ vector for an active Lorentz boost
00106     Vector3 betaVec() const {
00107       FourMomentum boost(_boostMatrix.getColumn(0)); //< @todo WRONG?!
00108       //cout << "!!!" << boost << endl;
00109       if (boost.isZero()) return Vector3();
00110       assert(boost.E() > 0);
00111       const double beta = boost.p3().mod() / boost.E();
00112       return boost.p3().unit() * beta;
00113     }
00114 
00115     /// Get the \f$ \beta \f$ factor
00116     double beta() const {
00117       return betaVec().mod();
00118     }
00119 
00120 
00121     /// Set up an active Lorentz boost from the \f$ \vec\gamma \f$ vector
00122     LorentzTransform& setGammaVec(const Vector3& vgamma) {
00123       const double gamma = vgamma.mod();
00124       const double beta = gamma2beta(gamma);
00125       _boostMatrix = Matrix<4>::mkIdentity();
00126       _boostMatrix.set(0, 0, gamma);
00127       _boostMatrix.set(1, 1, gamma);
00128       _boostMatrix.set(0, 1, +beta*gamma); //< +ve coeff since active boost
00129       _boostMatrix.set(1, 0, +beta*gamma); //< +ve coeff since active boost
00130       if (beta > 0) _boostMatrix = rotate(Vector3::mkX(), vgamma)._boostMatrix;
00131       return *this;
00132     }
00133 
00134     /// Get the \f$ \vec\gamma \f$ vector for an active Lorentz boost
00135     Vector3 gammaVec() const {
00136       FourMomentum boost(_boostMatrix.getColumn(0)); //< @todo WRONG?!
00137       if (boost.isZero()) return Vector3();
00138       assert(boost.E() > 0);
00139       const double beta = boost.p3().mod() / boost.E();
00140       return boost.p3().unit() * beta;
00141     }
00142 
00143     /// Get the \f$ \gamma \f$ factor
00144     double gamma() const {
00145       return beta2gamma(beta());
00146     }
00147 
00148     //@}
00149 
00150 
00151     /// Apply this transformation to the given 4-vector
00152     FourVector transform(const FourVector& v4) const {
00153       return multiply(_boostMatrix, v4);
00154     }
00155 
00156     /// Apply this transformation to the given 4-mometum
00157     FourMomentum transform(const FourMomentum& v4) const {
00158       return multiply(_boostMatrix, v4);
00159     }
00160 
00161 
00162     /// @name Transform modifications
00163     //@{
00164 
00165     /// Rotate the transformation cf. the difference between vectors @a from and @a to
00166     LorentzTransform rotate(const Vector3& from, const Vector3& to) const {
00167       return rotate(Matrix3(from, to));
00168     }
00169 
00170     /// Rotate the transformation by @a angle radians about @a axis
00171     LorentzTransform rotate(const Vector3& axis, double angle) const {
00172       return rotate(Matrix3(axis, angle));
00173     }
00174 
00175     /// Rotate the transformation by the 3D rotation matrix @a rot
00176     LorentzTransform rotate(const Matrix3& rot) const {
00177       LorentzTransform lt = *this;
00178       const Matrix4 rot4 = _mkMatrix4(rot);
00179       const Matrix4 newlt = rot4 * _boostMatrix * rot4.inverse();
00180       lt._boostMatrix = newlt;
00181       return lt;
00182     }
00183 
00184     /// Calculate the inverse transform
00185     LorentzTransform inverse() const {
00186       LorentzTransform rtn;
00187       rtn._boostMatrix = _boostMatrix.inverse();
00188       return rtn;
00189     }
00190 
00191     /// Combine LTs, treating @a this as the LH matrix.
00192     LorentzTransform combine(const LorentzTransform& lt) const {
00193       LorentzTransform rtn;
00194       rtn._boostMatrix = _boostMatrix * lt._boostMatrix;
00195       return rtn;
00196     }
00197 
00198     /// Operator combination of two LTs
00199     LorentzTransform operator*(const LorentzTransform& lt) const {
00200       return combine(lt);
00201     }
00202 
00203     /// Pre-multiply m3 by this LT
00204     LorentzTransform preMult(const Matrix3& m3) {
00205       _boostMatrix = multiply(_mkMatrix4(m3),_boostMatrix);
00206       return *this;
00207     }
00208 
00209     /// Post-multiply m3 by this LT
00210     LorentzTransform postMult(const Matrix3& m3) {
00211       _boostMatrix *= _mkMatrix4(m3);
00212       return *this;
00213     }
00214 
00215     //@}
00216 
00217 
00218     /// Return the matrix form
00219     Matrix4 toMatrix() const {
00220       return _boostMatrix;
00221     }
00222 
00223 
00224   private:
00225 
00226     Matrix4 _mkMatrix4(const Matrix3& m3) const {
00227       Matrix4 m4 = Matrix4::mkIdentity();
00228       for (size_t i = 0; i < 3; ++i) {
00229         for (size_t j = 0; j < 3; ++j) {
00230           m4.set(i+1, j+1, m3.get(i, j));
00231         }
00232       }
00233       return m4;
00234     }
00235 
00236     Matrix4 _boostMatrix;
00237 
00238   };
00239 
00240 
00241 
00242   inline LorentzTransform inverse(const LorentzTransform& lt) {
00243     return lt.inverse();
00244   }
00245 
00246   inline LorentzTransform combine(const LorentzTransform& a, const LorentzTransform& b) {
00247     return a.combine(b);
00248   }
00249 
00250   inline FourVector transform(const LorentzTransform& lt, const FourVector& v4) {
00251       return lt.transform(v4);
00252   }
00253 
00254 
00255   //////////////////////////
00256 
00257 
00258   inline string toString(const LorentzTransform& lt) {
00259     return toString(lt.toMatrix());
00260   }
00261 
00262   inline ostream& operator<<(std::ostream& out, const LorentzTransform& lt) {
00263     out << toString(lt);
00264     return out;
00265   }
00266 
00267 
00268 }
00269 
00270 #endif