[Rivet] Rivet 2.1.2: LorentzTransform with 0 Boost should do nothing

[Michael Große]

Hi,

I'm using Rivet 2.1.2 and I noticed incorrect behavior of the class
LorentzTransform, if it is created with a (0,0,0) boost-vector. In
this case the resulting boost-matrix should be the unity-matrix, since
we are already in the restframe and no boosting should be happening.

My current workaround is to add an if-condition to the constructor
LorentzTransform& setBoost(const Vector3& boost) in LorentzTrans.hh,
but I'm not sure if this is the proper place to fix this, since I
suspect that the error is somewhere in the calculation of the angle
between a vector and a vector with length 0.
I have also added another unittest to testMatVec.cc to check for this.

I have attached both my files to this email.

Cheers,

Michael
-------------- next part --------------
#ifndef RIVET_MATH_LORENTZTRANS
#define RIVET_MATH_LORENTZTRANS

#include <iostream>

#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"

namespace Rivet {


  inline double lorentzGamma(const double beta) {
    return 1.0 / sqrt(1 - beta*beta);
  }


  /// @brief Object implementing Lorentz transform calculations and boosts.
  class LorentzTransform {
    friend string toString(const LorentzTransform& lt);

  public:
    LorentzTransform() {
      _boostMatrix = Matrix<4>::mkIdentity();
    }

    LorentzTransform(const Vector3& boost) {
      setBoost(boost);
    }

    LorentzTransform(const double betaX, const double betaY, const double betaZ) {
      setBoost(betaX, betaY, betaZ);
    }

    LorentzTransform& setBoost(const Vector3& boost) {
      assert(boost.mod2() < 1);
      const double beta = boost.mod();
      const double gamma = lorentzGamma(beta);
      _boostMatrix = Matrix<4>::mkIdentity();
      _boostMatrix.set(0, 0, gamma);
      _boostMatrix.set(1, 1, gamma);
      // Positive coeff since these are active boosts
      _boostMatrix.set(0, 1, +beta*gamma);
      _boostMatrix.set(1, 0, +beta*gamma);
      if (beta>0.0){
        _boostMatrix = rotate(Vector3::mkX(), boost)._boostMatrix;
      }
      return *this;
    }

  //   LorentzTransform& setBoost(const Vector3& boostdirection, const double beta) {
  //     const Vector3 boost = boostdirection.unit() * beta;
  //     return setBoost(boost);
  //   }

    LorentzTransform& setBoost(const double betaX, const double betaY, const double betaZ) {
      return setBoost(Vector3(betaX, betaY, betaZ));
    }

    Vector3 boost() const {
      FourMomentum boost(_boostMatrix.getColumn(0));
      //cout << "!!!" << boost << endl;
      if (boost.isZero()) return boost;
      assert(boost.E() > 0);
      const double beta = boost.p3().mod() / boost.E();
      return boost.p3().unit() * beta;
    }

    double beta() const {
      return boost().mod();
    }

    double gamma() const {
      return lorentzGamma(beta());
    }

    LorentzTransform rotate(const Vector3& from, const Vector3& to) const {
      return rotate(Matrix3(from, to));
    }

    LorentzTransform rotate(const Vector3& axis, const double angle) const {
      return rotate(Matrix3(axis, angle));
    }

    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;
    }

    FourVector transform(const FourVector& v4) const {
      return multiply(_boostMatrix, v4);
    }

    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;
    }

    Matrix4 toMatrix() const {
      return _boostMatrix;
    }


    LorentzTransform operator*(const LorentzTransform& lt) const {
      return combine(lt);
    }

    LorentzTransform preMult(const Matrix3 & m3) {
      _boostMatrix = multiply(mkMatrix4(m3),_boostMatrix);
      return *this;
    }

    LorentzTransform postMult(const Matrix3 & m3) {
      _boostMatrix *= mkMatrix4(m3);
      return *this;
    }

  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;
    }


  private:
    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._boostMatrix);
  }

  inline ostream& operator<<(std::ostream& out, const LorentzTransform& lt) {
    out << toString(lt);
    return out;
  }


}

#endif
-------------- next part --------------
#include <iostream>
#include <limits>
#include <cassert>

#include "Rivet/Math/MathUtils.hh"
#include "Rivet/Math/Vectors.hh"
#include "Rivet/Math/Matrices.hh"

using namespace std;
using namespace Rivet;

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(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;

  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;
  const FourMomentum p3 = LorentzTransform(-p1.boostVector()).transform(p1);
  cout << p1 << " -> " << p3 << endl;
  cout << endl;

  LorentzTransform ltY(0,0.4,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).boost() << endl;
  cout << (ltY * ltX).boost() << endl;
  cout << (ltX * ltX.inverse()).boost() << endl;
  
  
  // if we are already in the restframe and there is no boost, then a Lorentz Transformation should do nothing, hence _boostmatrix should be unity
  LorentzTransform noBoost(0.0,0.0,0.0);
  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;
}