Revision 4bee12876976bec72820e6464bf8665e58e624d0 authored by Gerardo Ganis on 24 July 2009, 07:33:50 UTC, committed by Gerardo Ganis on 24 July 2009, 07:33:50 UTC
git-svn-id: http://root.cern.ch/svn/root/branches/v5-24-00-patches@29570 27541ba8-7e3a-0410-8455-c3a389f83636
1 parent 3517aba
Boost.cxx
// @(#)root/mathcore:$Id$
// Authors: M. Fischler 2005
/**********************************************************************
* *
* Copyright (c) 2005 , LCG ROOT FNAL MathLib Team *
* *
* *
**********************************************************************/
// Header file for class Boost, a 4x4 symmetric matrix representation of
// an axial Lorentz transformation
//
// Created by: Mark Fischler Mon Nov 1 2005
//
#include "Math/GenVector/Boost.h"
#include "Math/GenVector/LorentzVector.h"
#include "Math/GenVector/PxPyPzE4D.h"
#include "Math/GenVector/DisplacementVector3D.h"
#include "Math/GenVector/Cartesian3D.h"
#include "Math/GenVector/GenVector_exception.h"
#include <cmath>
#include <algorithm>
//#ifdef TEX
/**
A variable names bgamma appears in several places in this file. A few
words of elaboration are needed to make its meaning clear. On page 69
of Misner, Thorne and Wheeler, (Exercise 2.7) the elements of the matrix
for a general Lorentz boost are given as
\f[ \Lambda^{j'}_k = \Lambda^{k'}_j
= (\gamma - 1) n^j n^k + \delta^{jk} \f]
where the n^i are unit vectors in the direction of the three spatial
axes. Using the definitions, \f$ n^i = \beta_i/\beta \f$ , then, for example,
\f[ \Lambda_{xy} = (\gamma - 1) n_x n_y
= (\gamma - 1) \beta_x \beta_y/\beta^2 \f]
By definition, \f[ \gamma^2 = 1/(1 - \beta^2) \f]
so that \f[ \gamma^2 \beta^2 = \gamma^2 - 1 \f]
or \f[ \beta^2 = (\gamma^2 - 1)/\gamma^2 \f]
If we insert this into the expression for \f$ \Lambda_{xy} \f$, we get
\f[ \Lambda_{xy} = (\gamma - 1) \gamma^2/(\gamma^2 - 1) \beta_x \beta_y \f]
or, finally
\f[ \Lambda_{xy} = \gamma^2/(\gamma+1) \beta_x \beta_y \f]
The expression \f$ \gamma^2/(\gamma+1) \f$ is what we call <em>bgamma</em> in the code below.
\class ROOT::Math::Boost
*/
//#endif
namespace ROOT {
namespace Math {
void Boost::SetIdentity() {
// set identity boost
fM[kXX] = 1.0; fM[kXY] = 0.0; fM[kXZ] = 0.0; fM[kXT] = 0.0;
fM[kYY] = 1.0; fM[kYZ] = 0.0; fM[kYT] = 0.0;
fM[kZZ] = 1.0; fM[kZT] = 0.0;
fM[kTT] = 1.0;
}
void Boost::SetComponents (Scalar bx, Scalar by, Scalar bz) {
// set the boost beta as 3 components
Scalar bp2 = bx*bx + by*by + bz*bz;
if (bp2 >= 1) {
GenVector::Throw (
"Beta Vector supplied to set Boost represents speed >= c");
// SetIdentity();
return;
}
Scalar gamma = 1.0 / std::sqrt(1.0 - bp2);
Scalar bgamma = gamma * gamma / (1.0 + gamma);
fM[kXX] = 1.0 + bgamma * bx * bx;
fM[kYY] = 1.0 + bgamma * by * by;
fM[kZZ] = 1.0 + bgamma * bz * bz;
fM[kXY] = bgamma * bx * by;
fM[kXZ] = bgamma * bx * bz;
fM[kYZ] = bgamma * by * bz;
fM[kXT] = gamma * bx;
fM[kYT] = gamma * by;
fM[kZT] = gamma * bz;
fM[kTT] = gamma;
}
void Boost::GetComponents (Scalar& bx, Scalar& by, Scalar& bz) const {
// get beta of the boots as 3 components
Scalar gaminv = 1.0/fM[kTT];
bx = fM[kXT]*gaminv;
by = fM[kYT]*gaminv;
bz = fM[kZT]*gaminv;
}
DisplacementVector3D< Cartesian3D<Boost::Scalar> >
Boost::BetaVector() const {
// get boost beta vector
Scalar gaminv = 1.0/fM[kTT];
return DisplacementVector3D< Cartesian3D<Scalar> >
( fM[kXT]*gaminv, fM[kYT]*gaminv, fM[kZT]*gaminv );
}
void Boost::GetLorentzRotation (Scalar r[]) const {
// get Lorentz rotation corresponding to this boost as an array of 16 values
r[kLXX] = fM[kXX]; r[kLXY] = fM[kXY]; r[kLXZ] = fM[kXZ]; r[kLXT] = fM[kXT];
r[kLYX] = fM[kXY]; r[kLYY] = fM[kYY]; r[kLYZ] = fM[kYZ]; r[kLYT] = fM[kYT];
r[kLZX] = fM[kXZ]; r[kLZY] = fM[kYZ]; r[kLZZ] = fM[kZZ]; r[kLZT] = fM[kZT];
r[kLTX] = fM[kXT]; r[kLTY] = fM[kYT]; r[kLTZ] = fM[kZT]; r[kLTT] = fM[kTT];
}
void Boost::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.
if (fM[kTT] <= 0) {
GenVector::Throw (
"Attempt to rectify a boost with non-positive gamma");
return;
}
DisplacementVector3D< Cartesian3D<Scalar> > beta ( fM[kXT], fM[kYT], fM[kZT] );
beta /= fM[kTT];
if ( beta.mag2() >= 1 ) {
beta /= ( beta.R() * ( 1.0 + 1.0e-16 ) );
}
SetComponents ( beta );
}
LorentzVector< PxPyPzE4D<double> >
Boost::operator() (const LorentzVector< PxPyPzE4D<double> > & v) const {
// apply bosost to a PxPyPzE LorentzVector
Scalar x = v.Px();
Scalar y = v.Py();
Scalar z = v.Pz();
Scalar t = v.E();
return LorentzVector< PxPyPzE4D<double> >
( fM[kXX]*x + fM[kXY]*y + fM[kXZ]*z + fM[kXT]*t
, fM[kXY]*x + fM[kYY]*y + fM[kYZ]*z + fM[kYT]*t
, fM[kXZ]*x + fM[kYZ]*y + fM[kZZ]*z + fM[kZT]*t
, fM[kXT]*x + fM[kYT]*y + fM[kZT]*z + fM[kTT]*t );
}
void Boost::Invert() {
// invert in place boost (modifying the object)
fM[kXT] = -fM[kXT];
fM[kYT] = -fM[kYT];
fM[kZT] = -fM[kZT];
}
Boost Boost::Inverse() const {
// return inverse of boost
Boost tmp(*this);
tmp.Invert();
return tmp;
}
// ========== I/O =====================
std::ostream & operator<< (std::ostream & os, const Boost & b) {
// TODO - this will need changing for machine-readable issues
// and even the human readable form needs formatiing improvements
double m[16];
b.GetLorentzRotation(m);
os << "\n" << m[0] << " " << m[1] << " " << m[2] << " " << m[3];
os << "\n" << "\t" << " " << m[5] << " " << m[6] << " " << m[7];
os << "\n" << "\t" << " " << "\t" << " " << m[10] << " " << m[11];
os << "\n" << "\t" << " " << "\t" << " " << "\t" << " " << m[15] << "\n";
return os;
}
} //namespace Math
} //namespace ROOT
![swh spinner](/static/img/swh-spinner.gif)
Computing file changes ...