diff --git a/include/Rivet/Math/LorentzTrans.hh b/include/Rivet/Math/LorentzTrans.hh
--- a/include/Rivet/Math/LorentzTrans.hh
+++ b/include/Rivet/Math/LorentzTrans.hh
@@ -1,296 +1,303 @@
 #ifndef RIVET_MATH_LORENTZTRANS
 #define RIVET_MATH_LORENTZTRANS
 
 #include "Rivet/Math/MathHeader.hh"
 #include "Rivet/Math/MathUtils.hh"
 #include "Rivet/Math/MatrixN.hh"
 #include "Rivet/Math/Matrix3.hh"
 #include "Rivet/Math/Vector4.hh"
 #include <iostream>
 
 namespace Rivet {
 
 
   /// @brief Object implementing Lorentz transform calculations and boosts.
   ///
   /// @note These boosts are defined actively, i.e. as modifications of vectors
   /// rather than frame transformations. So the boost vector is the opposite of
   /// what you might expect.
   ///
   /// @todo Review the active/passive convention choice. Seems counterintuitive now...
   class LorentzTransform {
   public:
 
     /// @name Simple Lorentz factor conversions
     //@{
 
     /// Calculate the \f$ \gamma \f$ factor from \f$ \beta \f$
     static double beta2gamma(double beta) {
       return 1.0 / sqrt(1 - sqr(beta));
     }
 
     /// Calculate \f$ \beta \f$ from the \f$ \gamma \f$ factor
     static double gamma2beta(double gamma) {
       return sqrt(1 - sqr(1/gamma));
     }
 
     // /// Calculate the \f$ \gamma \f$ factor from \f$ \bar{\beta}^2 = 1-\beta^2 \f$
     // static double betabarsq2gamma(double betabarsq) {
     //   return 1.0 / sqrt(betabarsq);
     // }
 
     // /// Calculate \f$ \bar{\beta}^2 = 1-\beta^2 \f$ from the \f$ \gamma \f$ factor
     // static double gamma2betabarsq(double gamma) {
     //   return 1.0 / sqr(gamma);
     // }
 
     //@}
 
 
     /// @name Construction
     //@{
 
     /// Default (identity) constructor
     LorentzTransform() {
       _boostMatrix = Matrix<4>::mkIdentity();
     }
 
+    /// Constructor from a 4x4 matrix
+    LorentzTransform(const Matrix<4>& boostMatrix) {
+      _boostMatrix = boostMatrix;
+    }
 
     /// Make an LT for an active boost (i.e. object velocity += in boost direction)
     static LorentzTransform mkObjTransformFromBeta(const Vector3& vbeta) {
       LorentzTransform rtn;
       return rtn.setBetaVec(vbeta);
     }
 
     /// Make an LT for a passive boost (i.e. object velocity -= in boost direction)
     static LorentzTransform mkFrameTransformFromBeta(const Vector3& vbeta) {
       LorentzTransform rtn;
       return rtn.setBetaVec(-vbeta);
     }
 
     /// Make an LT for an active boost (i.e. object velocity += in boost direction)
     static LorentzTransform mkObjTransformFromGamma(const Vector3& vgamma) {
       LorentzTransform rtn;
       if (!vgamma.isZero()) rtn.setGammaVec(vgamma);
       return rtn;
     }
 
     /// Make an LT for a passive boost (i.e. object velocity -= in boost direction)
     static LorentzTransform mkFrameTransformFromGamma(const Vector3& vgamma) {
       LorentzTransform rtn;
       if (!vgamma.isZero()) rtn.setGammaVec(-vgamma);
       return rtn;
     }
 
     /// Make an LT for an active boost (i.e. object velocity += in boost direction)
     static LorentzTransform mkObjTransform(const FourMomentum& p4) {
       return mkObjTransformFromBeta(p4.betaVec());
     }
 
     /// Make an LT for a passive boost (i.e. object velocity -= in boost direction)
     static LorentzTransform mkFrameTransform(const FourMomentum& p4) {
       return mkObjTransformFromBeta(-p4.betaVec());
     }
 
     //@}
 
 
     /// @name Boost vector and beta/gamma factors
     //@{
 
     /// Set up an active Lorentz boost from the \f$ \vec\beta \f$ vector
     LorentzTransform& setBetaVec(const Vector3& vbeta) {
       // Set to identity for null boosts
       _boostMatrix = Matrix<4>::mkIdentity();
-      if (vbeta.mod2() == 0) return *this;
-      assert(vbeta.mod2() < 1);
+      if (isZero(vbeta.mod2())) return *this;
+      //assert(vbeta.mod2() < 1);
       //
       // It's along the x, y, or z axis if 2 Cartesian components are zero
-      const bool alongxyz = (int(isZero(vbeta.x())) + int(isZero(vbeta.y())) + int(isZero(vbeta.z())) == 2);
-      const int i = !isZero(vbeta.x()) ? 1 : !isZero(vbeta.y()) ? 2 : !isZero(vbeta.z()) ? 3 : 1;
+      const bool alongxyz = (int(vbeta.x() == 0) + int(vbeta.y() == 0) + int(vbeta.z() == 0) == 2);
+      const int i = (!alongxyz || vbeta.x() != 0) ? 1 : (vbeta.y() != 0) ? 2 : 3;
+      const int isign = !alongxyz ? 1 : sign(vbeta[i-1]);
+      // cout << boolalpha << alongxyz << ", i=" << i << endl;
+      // cout << vbeta << ": " << isZero(vbeta.x()) << ", " << isZero(vbeta.y()) << ", " << isZero(vbeta.z()) << endl;
+      //
       const double beta = vbeta.mod();
       const double gamma = beta2gamma(beta);
       _boostMatrix.set(0, 0, gamma);
       _boostMatrix.set(i, i, gamma);
-      _boostMatrix.set(0, i, +beta*gamma); //< +ve coeff since active boost
-      _boostMatrix.set(i, 0, +beta*gamma); //< +ve coeff since active boost
+      _boostMatrix.set(0, i, +isign*beta*gamma); //< +ve coeff since active boost
+      _boostMatrix.set(i, 0, +isign*beta*gamma); //< +ve coeff since active boost
       if (!alongxyz) _boostMatrix = rotate(Vector3::mkX(), vbeta)._boostMatrix;
       return *this;
     }
 
     /// Get the \f$ \vec\beta \f$ vector for an active Lorentz boost
     Vector3 betaVec() const {
       FourMomentum boost(_boostMatrix.getColumn(0)); //< @todo WRONG?!
       //cout << "!!!" << boost << endl;
       if (boost.isZero()) return Vector3();
       assert(boost.E() > 0);
       const double beta = boost.p3().mod() / boost.E();
       return boost.p3().unit() * beta;
     }
 
     /// Get the \f$ \beta \f$ factor
     double beta() const {
       return betaVec().mod();
     }
 
 
     /// Set up an active Lorentz boost from the \f$ \vec\gamma \f$ vector
     LorentzTransform& setGammaVec(const Vector3& vgamma) {
-      const double gamma = vgamma.mod();
-      const double beta = gamma2beta(gamma);
+      // Set to identity for null boosts
       _boostMatrix = Matrix<4>::mkIdentity();
-      _boostMatrix.set(0, 0, gamma);
-      _boostMatrix.set(1, 1, gamma);
-      _boostMatrix.set(0, 1, +beta*gamma); //< +ve coeff since active boost
-      _boostMatrix.set(1, 0, +beta*gamma); //< +ve coeff since active boost
-      if (beta > 0) _boostMatrix = rotate(Vector3::mkX(), vgamma)._boostMatrix;
-      return *this;
+      if (isZero(vgamma.mod2() - 1)) return *this;
+      //assert(vgamma.mod2() >= 1);
+      //
+      /// @todo Avoid loss of precision by using an internal beta *and* gamma function?
+      const double beta = gamma2beta(vgamma.mod());
+      return setBetaVec(beta * vgamma.unit());
     }
 
     /// Get the \f$ \vec\gamma \f$ vector for an active Lorentz boost
     Vector3 gammaVec() const {
       FourMomentum boost(_boostMatrix.getColumn(0)); //< @todo WRONG?!
       if (boost.isZero()) return Vector3();
       assert(boost.E() > 0);
       const double beta = boost.p3().mod() / boost.E();
       return boost.p3().unit() * beta;
     }
 
     /// Get the \f$ \gamma \f$ factor
     double gamma() const {
       return beta2gamma(beta());
     }
 
     //@}
 
 
     /// Apply this transformation to the given 4-vector
     FourVector transform(const FourVector& v4) const {
       return multiply(_boostMatrix, v4);
     }
 
     /// Apply this transformation to the given 4-mometum
     FourMomentum transform(const FourMomentum& v4) const {
       return multiply(_boostMatrix, v4);
     }
 
     /// Apply this transformation to the given 4-vector
     FourVector operator () (const FourVector& v4) const {
       return transform(v4);
     }
 
     /// Apply this transformation to the given 4-mometum
     FourMomentum operator () (const FourMomentum& v4) const {
       return transform(v4);
     }
 
 
     /// @name Transform modifications
     //@{
 
     /// Rotate the transformation cf. the difference between vectors @a from and @a to
     LorentzTransform rotate(const Vector3& from, const Vector3& to) const {
       return rotate(Matrix3(from, to));
     }
 
     /// Rotate the transformation by @a angle radians about @a axis
     LorentzTransform rotate(const Vector3& axis, double angle) const {
       return rotate(Matrix3(axis, angle));
     }
 
     /// Rotate the transformation by the 3D rotation matrix @a rot
     LorentzTransform rotate(const Matrix3& rot) const {
       LorentzTransform lt = *this;
       const Matrix4 rot4 = _mkMatrix4(rot);
       const Matrix4 newlt = rot4 * _boostMatrix * rot4.inverse();
       lt._boostMatrix = newlt;
       return lt;
     }
 
     /// Calculate the inverse transform
     LorentzTransform inverse() const {
       LorentzTransform rtn;
       rtn._boostMatrix = _boostMatrix.inverse();
       return rtn;
     }
 
     /// Combine LTs, treating @a this as the LH matrix.
     LorentzTransform combine(const LorentzTransform& lt) const {
       LorentzTransform rtn;
       rtn._boostMatrix = _boostMatrix * lt._boostMatrix;
       return rtn;
     }
 
     /// Operator combination of two LTs
     LorentzTransform operator*(const LorentzTransform& lt) const {
       return combine(lt);
     }
 
     /// Pre-multiply m3 by this LT
     LorentzTransform preMult(const Matrix3& m3) {
       _boostMatrix = multiply(_mkMatrix4(m3),_boostMatrix);
       return *this;
     }
 
     /// Post-multiply m3 by this LT
     LorentzTransform postMult(const Matrix3& m3) {
       _boostMatrix *= _mkMatrix4(m3);
       return *this;
     }
 
     //@}
 
 
     /// Return the matrix form
     Matrix4 toMatrix() const {
       return _boostMatrix;
     }
 
 
   private:
 
     Matrix4 _mkMatrix4(const Matrix3& m3) const {
       Matrix4 m4 = Matrix4::mkIdentity();
       for (size_t i = 0; i < 3; ++i) {
         for (size_t j = 0; j < 3; ++j) {
           m4.set(i+1, j+1, m3.get(i, j));
         }
       }
       return m4;
     }
 
     Matrix4 _boostMatrix;
 
   };
 
 
 
   inline LorentzTransform inverse(const LorentzTransform& lt) {
     return lt.inverse();
   }
 
   inline LorentzTransform combine(const LorentzTransform& a, const LorentzTransform& b) {
     return a.combine(b);
   }
 
   inline FourVector transform(const LorentzTransform& lt, const FourVector& v4) {
       return lt.transform(v4);
   }
 
 
   //////////////////////////
 
 
   inline string toString(const LorentzTransform& lt) {
     return toString(lt.toMatrix());
   }
 
   inline ostream& operator<<(std::ostream& out, const LorentzTransform& lt) {
     out << toString(lt);
     return out;
   }
 
 
 }
 
 #endif
diff --git a/test/testMatVec.cc b/test/testMatVec.cc
--- a/test/testMatVec.cc
+++ b/test/testMatVec.cc
@@ -1,155 +1,197 @@
 #include "Rivet/Math/MathUtils.hh"
 #include "Rivet/Math/Vectors.hh"
 #include "Rivet/Math/Matrices.hh"
-//#include "Rivet/Math/MatrixDiag.hh"
+#include "Rivet/Math/MatrixDiag.hh"
 using namespace Rivet;
 
 #include <iostream>
 #include <limits>
 #include <cassert>
 using namespace std;
 
+
+/// Set up an active Lorentz boost in the z direction
+LorentzTransform mkLTz(const double beta) {
+  assert(abs(beta) < 1);
+  const double gamma = LorentzTransform::beta2gamma(beta);
+  Matrix<4> rtn = Matrix<4>::mkIdentity();
+  rtn.set(0, 0, gamma);
+  rtn.set(3, 3, gamma);
+  rtn.set(0, 3, +beta*gamma); //< +ve coeff since active boost
+  rtn.set(3, 0, +beta*gamma); //< +ve coeff since active boost
+  return LorentzTransform(rtn);
+}
+
+
 int main() {
 
   FourVector a(1,0,0,0);
   cout << a << ": interval = " << a.invariant() << endl;
   assert(fuzzyEquals(a.invariant(), 1));
   a.setZ(1);
   assert(isZero(a.invariant()));
   cout << a << ": interval = " << a.invariant() << endl;
   a.setY(2).setZ(3);
   cout << a << ": interval = " << a.invariant() << endl;
   assert(fuzzyEquals(a.invariant(), -12));
   cout << a << ": vector = " << a.vector3() << endl << endl;
 
   FourMomentum b(1,0,0,0);
   cout << b << ": mass = " << b.mass() << endl;
   assert(fuzzyEquals(b.mass2(), 1));
   b.setPz(1);
   cout << b << ": mass = " << b.mass() << endl;
   assert(isZero(b.mass2()));
   b.setPy(2).setPz(3).setE(6);
   cout << b << ": mass = " << b.mass2() << endl;
   assert(fuzzyEquals(b.mass2(), 23));
   cout << b << ": vector = " << b.vector3() << endl << endl;
 
   Matrix3 m;
   m.set(0, 0, 7/4.0);
   m.set(0, 1, 3 * sqrt(3)/4.0);
   m.set(1, 0, 3 * sqrt(3)/4.0);
   m.set(1, 1, 13/4.0);
   m.set(2, 2, 9);
   cout << m << endl << endl;
   // EigenSystem<3> es = diagonalize(m);
 
   cout << "Matrices:" << endl;
   cout << Matrix3() << endl;
   cout << Matrix3::mkIdentity() << endl;
   const Matrix3 I3 = Matrix3::mkIdentity();
   cout << Matrix3::mkIdentity() * m * I3 << endl;
   cout << "tr(0) & det(0): " << Matrix3().trace() << ", " << Matrix3().det() << endl;
   cout << "tr(I3) & det(I3): " << I3.trace() << ", " << I3.det() << endl;
   Matrix3 m1 = Matrix3::mkIdentity();
   Matrix3 m2 = m1;
   m1.setRow(1, Vector3(1,2,3));
   m2.setColumn(1, Vector3(3,2,1));
   Matrix3 m3 = Matrix3::mkZero();
   cout << m1 << " + " << m2 << " = " << m1 + m2 << endl;
   m3.setRow(0, Vector3(2,3,0)).setRow(1, Vector3(1,4,3)).setRow(2, Vector3(0,1,2));
   cout << m1+m2 << " == " << m3 << ": " << (m1+m2 == m3 ? "true" : "false") << endl;
   cout << endl;
 
 
   Vector3 v3(1,2,3);
   cout << "Vector: " << v3 << endl;
   cout << "Invert: " << v3 << " --> " << -v3 << endl;
   const Matrix3 rot90(Vector3(0,0,1), PI/2.0);
   const Matrix3 rot90m(Vector3(0,0,1), -PI/2.0);
   const Matrix3 rot180(Vector3(0,0,1), PI);
   const Matrix3 rot180m(Vector3(0,0,1), -PI);
   const Vector3 v3_90 = rot90*v3;
   cout << "Rot 90: " << v3 << " ---> " << v3_90 << endl;
   const Vector3 v3_90m = rot90m*v3;
   cout << "Rot -90: " << v3 << " ---> " << v3_90m << endl;
   const Vector3 v3_180 = rot180*v3;
   cout << "Rot 180: " << v3 << " ---> " << v3_180 << endl;
   const Vector3 v3_180m = rot180m*v3;
   cout << "Rot -180: " << v3 << " ---> " << v3_180m << endl;
   assert(fuzzyEquals(v3_180, v3_180m));
 
   const Vector3 v3_9090 = rot90*rot90*v3;
   cout << "Rot 2 x 90: " << v3 << " ---> " << v3_9090 << endl;
   assert(fuzzyEquals(v3_180, v3_9090));
 
   const Vector3 v3_90m90m = rot90m*rot90m*v3;
   cout << "Rot 2 x -90: " << v3 << " ---> " << v3_90m90m << endl;
   assert(fuzzyEquals(v3_180, v3_90m90m));
 
   const Vector3 v3_9090m = rot90*rot90m*v3;
   const Vector3 v3_90m90 = rot90m*rot90*v3;
   cout << "Rot 90*-90: "<< v3 << " ---> " << v3_9090m << endl;
   cout << "Rot -90*90: "<< v3 << " ---> " << v3_90m90 << endl;
   assert(fuzzyEquals(v3, v3_9090m));
   assert(fuzzyEquals(v3, v3_90m90));
 
   const Vector3 v3_90i = rot90.inverse()*v3;
   cout << "Rot (90)^-1: "<< v3 << " ---> " << v3_90i << endl;
   assert(fuzzyEquals(v3_90i, v3_90m));
 
   const Vector3 v3_9090i = rot90*rot90.inverse()*v3;
   const Vector3 v3_90i90 = rot90.inverse()*rot90*v3;
   cout << "Rot 90*(90)^-1: "<< v3 << " ---> " << v3_9090i << endl;
   cout << "Rot (90)^-1*90: "<< v3 << " ---> " << v3_90i90 << endl;
   assert(fuzzyEquals(v3, v3_9090i));
   assert(fuzzyEquals(v3, v3_90i90));
 
   const Matrix3 rot1(Vector3(0,1,0), PI/180.0);
   cout << "Rot 0 x 45 x 1: " << v3 << endl;
   for (size_t i = 0; i < 8; ++i) {
     for (size_t j = 0; j < 45; ++j) {
       v3 = rot1*v3;
     }
     cout << "Rot " << i+1 << " x 45 x 1: " << v3 << endl;
   }
   assert(fuzzyEquals(v3, Vector3(1,2,3)));
   cout << endl;
 
   cout << "Boosts:" << endl;
   LorentzTransform ltX = LorentzTransform::mkObjTransformFromBeta(Vector3(0.5, 0, 0));
   cout << "LTx: " << ltX << endl;
   cout << "I on LTx: " << ltX.rotate(Matrix3::mkIdentity()) << endl;
   cout << "Rot90 on LTx: " << ltX.rotate(rot90) << endl;
   cout << endl;
 
+  LorentzTransform ltY = LorentzTransform::mkObjTransformFromBeta(Vector3(0, 0.5, 0));
+  cout << "LTy: " << ltY << endl;
+  cout << "I on LTy: " << ltY.rotate(Matrix3::mkIdentity()) << endl;
+  cout << "Rot90 on LTy: " << ltY.rotate(rot90) << endl;
+  cout << endl;
+
+  LorentzTransform ltZ = LorentzTransform::mkObjTransformFromBeta(Vector3(0, 0, 0.5));
+  cout << "LTz: " << ltZ << endl;
+  cout << "I on LTz: " << ltZ.rotate(Matrix3::mkIdentity()) << endl;
+  cout << "Rot90 on LTz: " << ltZ.rotate(rot90) << endl;
+  cout << endl;
+
+  LorentzTransform ltZ2 = mkLTz(0.5);
+  cout << "LTz': " << ltZ2 << endl;
+  cout << endl;
+
+  LorentzTransform ltXYZ = LorentzTransform::mkObjTransformFromBeta(Vector3(0.1, -0.2, 0.3));
+  cout << "LTxyz: " << ltXYZ << endl;
+  cout << "I on LTxyz: " << ltXYZ.rotate(Matrix3::mkIdentity()) << endl;
+  cout << "Rot90 on LTxyz: " << ltXYZ.rotate(rot90) << endl;
+  cout << endl;
+
   cout << "X-boosts:" << endl;
   const FourMomentum p1 = FourMomentum(10, 0, 0, 1);
   const FourMomentum p2 = ltX.transform(p1);
   cout << p1 << " -> " << p2 << endl;
   cout << p2 << " -> " << ltX.inverse().transform(p2) << endl;
   //cout << p1.boostVector() << endl;
+  cout << "LT(p1) = " << LorentzTransform::mkFrameTransformFromBeta(p1.boostVector()) << endl;
   const FourMomentum p3 = LorentzTransform::mkFrameTransformFromBeta(p1.boostVector()).transform(p1);
   cout << p1 << " -> " << p3 << endl;
+  assert(isZero(p3.x()));
+  assert(isZero(p3.y()));
+  assert(isZero(p3.z()));
+  assert(p3.E() < p1.E());
+  assert(fuzzyEquals(p3.E(), p1.mass()));
   cout << endl;
 
-  LorentzTransform ltY = LorentzTransform::mkObjTransformFromGamma(Vector3(0, 1.2, 0));
+  // LorentzTransform ltY = LorentzTransform::mkObjTransformFromGamma(Vector3(0, 1.2, 0));
   cout << FourMomentum(1,0,0,1) << " -> " //<< "\n  "
        << (ltX * ltY).transform(FourMomentum(1,0,0,1)) << endl;
   cout << FourMomentum(1,0,0,1) << " -> " //<< "\n  "
        << (ltY * ltX).transform(FourMomentum(1,0,0,1)) << endl;
   cout << (ltX * ltY).betaVec() << endl;
   cout << (ltY * ltX).betaVec() << endl;
   cout << (ltX * ltX.inverse()).betaVec() << endl;
 
   // If we are already in the rest frame and there is no boost, then LT is trivial/identity
   LorentzTransform noBoost;
   cout << "Element  0,0 should be 1 and is " << noBoost.toMatrix().get(0,0) << endl;
   assert(noBoost.toMatrix().get(0,0)==1);
   cout << "Element  0,1 should be 0 and is " << noBoost.toMatrix().get(0,1) << endl;
   assert(noBoost.toMatrix().get(0,1)==0);
   cout << "Element  1,0 should be 0 and is " << noBoost.toMatrix().get(1,0) << endl;
   assert(noBoost.toMatrix().get(1,0)==0);
   cout << "Element  1,1 should be 1 and is " << noBoost.toMatrix().get(1,1) << endl;
   assert(noBoost.toMatrix().get(1,1)==1);
 
   return EXIT_SUCCESS;
 }