LorentzRotation.cxx

Go to the documentation of this file.

// @(#)root/mathcore:$Id$
// Authors: W. Brown, M. Fischler, L. Moneta    2005
 
/**********************************************************************
 *                                                                    *
 * Copyright (c) 2005 , LCG ROOT FNAL MathLib Team                    *
 *                                                                    *
 *                                                                    *
 **********************************************************************/
 
// Header file for class LorentzRotation, a 4x4 matrix representation of
// a general Lorentz transformation
//
// Created by: Mark Fischler Mon Aug 8  2005
//
 
#include "MathX/GenVectorX/GenVectorIO.h"
 
#include "MathX/GenVectorX/LorentzRotation.h"
#include "MathX/GenVectorX/LorentzVector.h"
#include "MathX/GenVectorX/PxPyPzE4D.h"
#include "MathX/GenVectorX/GenVector_exception.h"
 
#include <cmath>
#include <algorithm>
 
#include "MathX/GenVectorX/Rotation3D.h"
#include "MathX/GenVectorX/RotationX.h"
#include "MathX/GenVectorX/RotationY.h"
#include "MathX/GenVectorX/RotationZ.h"
 
#include "MathX/GenVectorX/AccHeaders.h"
 
#include "MathX/GenVectorX/MathHeaders.h"
 
namespace ROOT {
 
namespace ROOT_MATH_ARCH {
 
LorentzRotation::LorentzRotation()
{
   // constructor of an identity LR
   fM[kXX] = 1.0;
   fM[kXY] = 0.0;
   fM[kXZ] = 0.0;
   fM[kXT] = 0.0;
   fM[kYX] = 0.0;
   fM[kYY] = 1.0;
   fM[kYZ] = 0.0;
   fM[kYT] = 0.0;
   fM[kZX] = 0.0;
   fM[kZY] = 0.0;
   fM[kZZ] = 1.0;
   fM[kZT] = 0.0;
   fM[kTX] = 0.0;
   fM[kTY] = 0.0;
   fM[kTZ] = 0.0;
   fM[kTT] = 1.0;
}
 
LorentzRotation::LorentzRotation(Rotation3D const &r)
{
   // construct from  Rotation3D
   r.GetComponents(fM[kXX], fM[kXY], fM[kXZ], fM[kYX], fM[kYY], fM[kYZ], fM[kZX], fM[kZY], fM[kZZ]);
   fM[kXT] = 0.0;
   fM[kYT] = 0.0;
   fM[kZT] = 0.0;
   fM[kTX] = 0.0;
   fM[kTY] = 0.0;
   fM[kTZ] = 0.0;
   fM[kTT] = 1.0;
}
 
LorentzRotation::LorentzRotation(AxisAngle const &a)
{
   // construct from  AxisAngle
   const Rotation3D r(a);
   r.GetComponents(fM[kXX], fM[kXY], fM[kXZ], fM[kYX], fM[kYY], fM[kYZ], fM[kZX], fM[kZY], fM[kZZ]);
   fM[kXT] = 0.0;
   fM[kYT] = 0.0;
   fM[kZT] = 0.0;
   fM[kTX] = 0.0;
   fM[kTY] = 0.0;
   fM[kTZ] = 0.0;
   fM[kTT] = 1.0;
}
 
LorentzRotation::LorentzRotation(EulerAngles const &e)
{
   // construct from  EulerAngles
   const Rotation3D r(e);
   r.GetComponents(fM[kXX], fM[kXY], fM[kXZ], fM[kYX], fM[kYY], fM[kYZ], fM[kZX], fM[kZY], fM[kZZ]);
   fM[kXT] = 0.0;
   fM[kYT] = 0.0;
   fM[kZT] = 0.0;
   fM[kTX] = 0.0;
   fM[kTY] = 0.0;
   fM[kTZ] = 0.0;
   fM[kTT] = 1.0;
}
 
LorentzRotation::LorentzRotation(Quaternion const &q)
{
   // construct from Quaternion
   const Rotation3D r(q);
   r.GetComponents(fM[kXX], fM[kXY], fM[kXZ], fM[kYX], fM[kYY], fM[kYZ], fM[kZX], fM[kZY], fM[kZZ]);
   fM[kXT] = 0.0;
   fM[kYT] = 0.0;
   fM[kZT] = 0.0;
   fM[kTX] = 0.0;
   fM[kTY] = 0.0;
   fM[kTZ] = 0.0;
   fM[kTT] = 1.0;
}
 
LorentzRotation::LorentzRotation(RotationX const &r)
{
   // construct from  RotationX
   Scalar s = r.SinAngle();
   Scalar c = r.CosAngle();
   fM[kXX] = 1.0;
   fM[kXY] = 0.0;
   fM[kXZ] = 0.0;
   fM[kXT] = 0.0;
   fM[kYX] = 0.0;
   fM[kYY] = c;
   fM[kYZ] = -s;
   fM[kYT] = 0.0;
   fM[kZX] = 0.0;
   fM[kZY] = s;
   fM[kZZ] = c;
   fM[kZT] = 0.0;
   fM[kTX] = 0.0;
   fM[kTY] = 0.0;
   fM[kTZ] = 0.0;
   fM[kTT] = 1.0;
}
 
LorentzRotation::LorentzRotation(RotationY const &r)
{
   // construct from  RotationY
   Scalar s = r.SinAngle();
   Scalar c = r.CosAngle();
   fM[kXX] = c;
   fM[kXY] = 0.0;
   fM[kXZ] = s;
   fM[kXT] = 0.0;
   fM[kYX] = 0.0;
   fM[kYY] = 1.0;
   fM[kYZ] = 0.0;
   fM[kYT] = 0.0;
   fM[kZX] = -s;
   fM[kZY] = 0.0;
   fM[kZZ] = c;
   fM[kZT] = 0.0;
   fM[kTX] = 0.0;
   fM[kTY] = 0.0;
   fM[kTZ] = 0.0;
   fM[kTT] = 1.0;
}
 
LorentzRotation::LorentzRotation(RotationZ const &r)
{
   // construct from  RotationX
   Scalar s = r.SinAngle();
   Scalar c = r.CosAngle();
   fM[kXX] = c;
   fM[kXY] = -s;
   fM[kXZ] = 0.0;
   fM[kXT] = 0.0;
   fM[kYX] = s;
   fM[kYY] = c;
   fM[kYZ] = 0.0;
   fM[kYT] = 0.0;
   fM[kZX] = 0.0;
   fM[kZY] = 0.0;
   fM[kZZ] = 1.0;
   fM[kZT] = 0.0;
   fM[kTX] = 0.0;
   fM[kTY] = 0.0;
   fM[kTZ] = 0.0;
   fM[kTT] = 1.0;
}
 
void LorentzRotation::Rectify()
{
   // Assuming the representation of this is close to a true Lorentz Rotation,
   // but may have drifted due to round-off error from many operations,
   // this forms an "exact" orthosymplectic matrix for the Lorentz Rotation
   // again.
 
   typedef LorentzVector<PxPyPzE4D<Scalar>> FourVector;
   if (fM[kTT] <= 0) {
#if !defined(ROOT_MATH_SYCL) && !defined(ROOT_MATH_CUDA)
      GenVector_Throw("LorentzRotation:Rectify(): Non-positive TT component - cannot rectify");
#endif
      return;
   }
   FourVector t(fM[kTX], fM[kTY], fM[kTZ], fM[kTT]);
   Scalar m2 = t.M2();
   if (m2 <= 0) {
#if !defined(ROOT_MATH_SYCL) && !defined(ROOT_MATH_CUDA)
      GenVector_Throw("LorentzRotation:Rectify(): Non-timelike time row - cannot rectify");
#endif
      return;
   }
   t /= math_sqrt(m2);
   FourVector z(fM[kZX], fM[kZY], fM[kZZ], fM[kZT]);
   z = z - z.Dot(t) * t;
   m2 = z.M2();
   if (m2 >= 0) {
#if !defined(ROOT_MATH_SYCL) && !defined(ROOT_MATH_CUDA)
      GenVector_Throw("LorentzRotation:Rectify(): Non-spacelike Z row projection - "
                      "cannot rectify");
#endif
      return;
   }
   z /= math_sqrt(-m2);
   FourVector y(fM[kYX], fM[kYY], fM[kYZ], fM[kYT]);
   y = y - y.Dot(t) * t - y.Dot(z) * z;
   m2 = y.M2();
   if (m2 >= 0) {
#if !defined(ROOT_MATH_SYCL) && !defined(ROOT_MATH_CUDA)
      GenVector_Throw("LorentzRotation:Rectify(): Non-spacelike Y row projection - "
                      "cannot rectify");
#endif
      return;
   }
   y /= math_sqrt(-m2);
   FourVector x(fM[kXX], fM[kXY], fM[kXZ], fM[kXT]);
   x = x - x.Dot(t) * t - x.Dot(z) * z - x.Dot(y) * y;
   m2 = x.M2();
   if (m2 >= 0) {
#if !defined(ROOT_MATH_SYCL) && !defined(ROOT_MATH_CUDA)
      GenVector_Throw("LorentzRotation:Rectify(): Non-spacelike X row projection - "
                      "cannot rectify");
#endif
      return;
   }
   x /= math_sqrt(-m2);
}
 
void LorentzRotation::Invert()
{
   // invert modifying current content
   Scalar temp;
   temp = fM[kXY];
   fM[kXY] = fM[kYX];
   fM[kYX] = temp;
   temp = fM[kXZ];
   fM[kXZ] = fM[kZX];
   fM[kZX] = temp;
   temp = fM[kYZ];
   fM[kYZ] = fM[kZY];
   fM[kZY] = temp;
   temp = fM[kXT];
   fM[kXT] = -fM[kTX];
   fM[kTX] = -temp;
   temp = fM[kYT];
   fM[kYT] = -fM[kTY];
   fM[kTY] = -temp;
   temp = fM[kZT];
   fM[kZT] = -fM[kTZ];
   fM[kTZ] = -temp;
}
 
LorentzRotation LorentzRotation::Inverse() const
{
   // return an inverse LR
   return LorentzRotation(fM[kXX], fM[kYX], fM[kZX], -fM[kTX], fM[kXY], fM[kYY], fM[kZY], -fM[kTY], fM[kXZ], fM[kYZ],
                          fM[kZZ], -fM[kTZ], -fM[kXT], -fM[kYT], -fM[kZT], fM[kTT]);
}
 
LorentzRotation LorentzRotation::operator*(const LorentzRotation &r) const
{
   // combination with another LR
   return LorentzRotation(fM[kXX] * r.fM[kXX] + fM[kXY] * r.fM[kYX] + fM[kXZ] * r.fM[kZX] + fM[kXT] * r.fM[kTX],
                          fM[kXX] * r.fM[kXY] + fM[kXY] * r.fM[kYY] + fM[kXZ] * r.fM[kZY] + fM[kXT] * r.fM[kTY],
                          fM[kXX] * r.fM[kXZ] + fM[kXY] * r.fM[kYZ] + fM[kXZ] * r.fM[kZZ] + fM[kXT] * r.fM[kTZ],
                          fM[kXX] * r.fM[kXT] + fM[kXY] * r.fM[kYT] + fM[kXZ] * r.fM[kZT] + fM[kXT] * r.fM[kTT],
                          fM[kYX] * r.fM[kXX] + fM[kYY] * r.fM[kYX] + fM[kYZ] * r.fM[kZX] + fM[kYT] * r.fM[kTX],
                          fM[kYX] * r.fM[kXY] + fM[kYY] * r.fM[kYY] + fM[kYZ] * r.fM[kZY] + fM[kYT] * r.fM[kTY],
                          fM[kYX] * r.fM[kXZ] + fM[kYY] * r.fM[kYZ] + fM[kYZ] * r.fM[kZZ] + fM[kYT] * r.fM[kTZ],
                          fM[kYX] * r.fM[kXT] + fM[kYY] * r.fM[kYT] + fM[kYZ] * r.fM[kZT] + fM[kYT] * r.fM[kTT],
                          fM[kZX] * r.fM[kXX] + fM[kZY] * r.fM[kYX] + fM[kZZ] * r.fM[kZX] + fM[kZT] * r.fM[kTX],
                          fM[kZX] * r.fM[kXY] + fM[kZY] * r.fM[kYY] + fM[kZZ] * r.fM[kZY] + fM[kZT] * r.fM[kTY],
                          fM[kZX] * r.fM[kXZ] + fM[kZY] * r.fM[kYZ] + fM[kZZ] * r.fM[kZZ] + fM[kZT] * r.fM[kTZ],
                          fM[kZX] * r.fM[kXT] + fM[kZY] * r.fM[kYT] + fM[kZZ] * r.fM[kZT] + fM[kZT] * r.fM[kTT],
                          fM[kTX] * r.fM[kXX] + fM[kTY] * r.fM[kYX] + fM[kTZ] * r.fM[kZX] + fM[kTT] * r.fM[kTX],
                          fM[kTX] * r.fM[kXY] + fM[kTY] * r.fM[kYY] + fM[kTZ] * r.fM[kZY] + fM[kTT] * r.fM[kTY],
                          fM[kTX] * r.fM[kXZ] + fM[kTY] * r.fM[kYZ] + fM[kTZ] * r.fM[kZZ] + fM[kTT] * r.fM[kTZ],
                          fM[kTX] * r.fM[kXT] + fM[kTY] * r.fM[kYT] + fM[kTZ] * r.fM[kZT] + fM[kTT] * r.fM[kTT]);
}
 
#if !defined(ROOT_MATH_SYCL) && !defined(ROOT_MATH_CUDA)
 
std::ostream &operator<<(std::ostream &os, const LorentzRotation &r)
{
   // TODO - this will need changing for machine-readable issues
   //        and even the human readable form needs formatting improvements
   double m[16];
   r.GetComponents(m, m + 16);
   os << "\n" << m[0] << "  " << m[1] << "  " << m[2] << "  " << m[3];
   os << "\n" << m[4] << "  " << m[5] << "  " << m[6] << "  " << m[7];
   os << "\n" << m[8] << "  " << m[9] << "  " << m[10] << "  " << m[11];
   os << "\n" << m[12] << "  " << m[13] << "  " << m[14] << "  " << m[15] << "\n";
   return os;
}
#endif
 
} // namespace ROOT_MATH_ARCH
} // namespace ROOT