Hi,
I was looking at the documentation of TLorentzVector on the web site
of ROOT as the version of root I had here on disk (2.21) did not have
these classes. I will check them and make sure that nothing has
changed, but I suspect the problem I mention is still there as the
class internals would need to be completely redesigned to change this
problem.
Rene Brun writes:
> Hi Nick, Hi Jeff,
> As both of you know, we have introduced in 2.22 new versions of the
> classes
> TVector3 and TLorentzVector. We had several iterations with the FNAL
> team
> supporting ZOOM & CLHEP. We finally agreed on an implementation that
> you can find in 2.22/09 released yesterday.
I am not surprised that in HEP these problems are not present,
probably at HEP energies for most particles you can treat the masses
as irrelevant (E approx eq. p).
> Since both of you do not specify which version of Root you are referring
> to,
> I am not sure that you are talking of the same thing.
> During the discussions with FNAL, we investigated the point if
> TLorentzVector
> should derive from TVector3 or not. There are pros/cons with the two
> approaches.
I don't have too much an opinion on this, I could understand arguments
for and against it.
> The idea was to converge with the CLHEP versions of LorentzVector.
> We are convinced that for frequently used classes such as
> TLorentzVector,
> it is not a good idea that everybody implements his own version.
> We had many comments from users that a convergence should be possible
> and it is a very frustating situation to see different implementations.
Yes, this is after all why OO programming is supposed to be so much
better. If you ask me, the jury is still out ;-)
> I would like to understand what are precisely the criticisms with the
> current implementation.
> Jeff, please give more details on the problems.
The problem is when Lorentz vector is used to represent the momentum
of a real particle e.g. an electron or proton. This particle has an
invariant characteristic which is its rest mass, which is also the
norm of the four-vector. I have found that if one takes the total
relativistic energy and the three momentum components of a particle as
the "member data" of a four-vector class, the particle can lose its
identity through boosts and transformations, i.e. you will not get the
same value for the invariant square (or rest mass) before and after
transformations. This is particularly true for high-energy
electrons. Unfortunately this can cause problems in kinematics
calculations where different pieces of a relation are computed in
different frames for simplicity. I don't have any examples on hand,
sorry (I fixed them all apparently).
So to reiterate, I believe that four a Four-vector object representing
the momentum of a particle, the rest mass should be chosen as one of
the four member data instead of the energy, since the particle should
always preserve its identity even under boosts and such.
I include my Python implementation below for comparison, as well as an
example program which uses the implementation. The Python
implementation of Vec4 does not derive from Vec3 but it uses a Vec3
object to represent the spatial component of the four-vector.
I hope the attachments below work as advertised a few days ago by Matt
Langston.
Jeff
## -*- mode: python -*-
## Calculate d(e,e'p) kinematics
## for fixed value of neutron momentum (= p_m) and fixed Q^2,
## generate kinematics for angular distribution of neutron w.r.t. vector q
## formalism can be found in MRN I:
from vec4 import *
from math import *
from masstable import *
## "less fixed" parameters
BEAMEN =4745.0 # MeV
pn = 500.0 # MeV/c
Q2 = 2.00 # GeV/c
Q2 = Q2 * 1.0E6 # convert to MeV/c
## given above, below doesn't change:
e0 = Vec4(pz=BEAMEN, m0=mass['e'])
targ = Vec4(m0=mass['d'])
## specify neutron momentum and calculate stuff that doesn't change w/ angle
En = sqrt(pn**2 + mass['n']**2)
sigm = mass['d'] - En
lamd = (Q2/2 + En*mass['d']) - (mass['d']**2 + mass['n']**2 - mass['p']**2)/2
## print labels
print "Beam Energy %5.3f GeV, Q2 = %5.3f (GeV/c)**2, p_n = %5.1f MeV/c" % \
(BEAMEN/1000, Q2/1.0E6, pn)
print "th n/q theta_e mom_e theta_p mom_p T_p th_q theta_pq"
print "(deg) (deg) (MeV/c) (deg) (MeV/c) (MeV) (deg) (deg)"
## loop over neutron angle
for nangl in range(0,181,10) :
thet = nangl * pi/180.0
denom = (sigm**2 - (pn*cos(thet))**2)
b = -lamd*sigm / denom
c = (lamd**2 - Q2*(pn*cos(thet))**2) / denom
if nangl == 90 :
omega = -b # roundoff problems with sqrt avoided
elif nangl > 90 :
omega = -b + sqrt(b**2 - c)
else :
omega = -b - sqrt(b**2 - c)
q3 = sqrt(Q2 + omega**2) # vector q
arg = 1 - (Q2 + 2*mass['e']**2)/(2*BEAMEN*(BEAMEN-omega))
ethet = acos(arg)
scatang = 180*ethet/pi
## construct scattered-electron fourvec and fourvec q, check calc.
ef = Vec4(Vec3(mag=BEAMEN-omega,theta=scatang,phi=-90),m0=mass['e'])
q = e0 - ef
## construct neutron four-vector in lab
dstar = q + targ
angn = q.thet() - nangl
if angn > 0 :
phin = q.phi()
else :
angn = -angn
phin = 180 + q.phi()
neut = Vec4(Vec3(mag=pn,theta=angn,phi=phin),m0=mass['n'])
prot = dstar - neut
cth = prot.p.hat() * q.p.hat()
## code below handles problems with range over/underflow due to
## roundoff
if abs( abs(cth) - 1.0) < 1.0E-06 :
if cth < 0.0 :
cth = -1.0
else :
cth = 1.0
thetapq = acos( cth )
thetapq = thetapq * 180/pi
print \
"%3d %5.2f %7.2f %5.2f %7.2f %7.2f %6.2f %6.2f" % \
(nangl, scatang, ef.P(), prot.thet(), prot.P(), prot.T(),
q.thet(), thetapq)
# Vec4 class: implementation of four-vectors
# Original version Fourvec: J. A. Templon, NIKHEF-K, Amsterdam
# New version: same person but now at University of Georgia
#
# You can use any units you want for the masses, momenta etc. as long as
# they are consistent, and don't require multiplication by c to go
# back and forth. An example which I use all the time is for all masses
# in MeV/c^2, momenta in MeV/c, and energies in MeV.
#
# Angles must be entered in degrees
#
# Usage: initialize with
#
# from vec4 import Vec4
#
# a) your_fourvec = Vec4(3vec,m0=restmass) - give rest mass
# b) your_fourvec = Vec4(3vec,W2=4vsquare) - give invariant square
# c) your_fourvec = Vec4(3vec,Etot=total_energy) - give total energy
# d) any of the three above can also be used when you substitute
# px= , py=, pz= for the 3vec part. Omitting the 3vec means
# a particle at rest.
#
# 3vec is a three-momentum object (see class vec3)
#
# example:
# e0 = Vec4(Vec3(x=0.0, y=0.0, z=855.0), m0=0.510999)
# or
# e0 = Vec4(pz=855.00, m0=0.510999)
#
# creates a four-momentum for an electron travelling at 855.0 MeV/c
# along the z-axis.
#
# attributes:
#
# vec4.W2 - the invariant square of the four-momentum (MeV/c)**2
# vec4.p - the three momentum as a Vec3 object
#
# methods returning output:
#
# vec4.m0() - returns the particle rest mass in MeV/c^2
# vec4.E() - returns the particle total energy (relativistic) in MeV/c^2
# vec4.T() - returns the particle kinetic energy in MeV
# vec4.P() - returns the absolute magnitude of the 3-momentum (MeV/c)
# vec4.thet() - returns the polar angle of the momentum 3-vector (degrees)
# vec4.phi() - returns the azim. angle of the momentum 3-vector (degrees)
# vec4.gamma() - returns the Lorentz factor
# vec4.beta() - returns the reduced velocity
#
# operators:
#
# "-" (__sub__) - subtraction of two four-vectors
# "+" (__add__) - addition of two four-vectors
# "*" (__mul__) - scalar product of two four-vectors
# __repr__ - printing Vec4 objects
import math, string
from vec3 import *
DEGRAD = 180.0/math.pi # conversion radians -> degrees
class Vec4:
def __init__(self, vec3=0.0, px=0.0, py=0.0, pz=0.0,
m0=0.0, W2=0.0, Etot=0.0) :
## first get the vec3 part
if vec3 <> 0.0 :
self.p = vec3
elif px <> 0.0 or py <> 0.0 or pz <> 0.0 :
self.p = Vec3(cart=(px, py, pz))
else :
self.p = Vec3( )
##-> set the invariant square
if m0 <> 0.0 :
self.W2 = m0**2
elif W2 <> 0.0 :
self.W2 = W2
elif Etot <> 0.0 :
self.W2 = Etot**2 - self.p*self.p
else :
self.W2 = 0.0 # assume it's a photon
# end of __init__
def P(self) : # returns the magnitude of the 3-momentum
return self.p.norm()
def m0(self) : # returns the rest mass
if self.W2 < 0.0 : # then the rest mass is not a valid concept
return -1.0E+60 # hopefully this will cause an error somewhere
else :
return math.sqrt(self.W2)
def E(self) : # returns the total energy
return math.sqrt(self.W2 + self.p*self.p)
def T(self) :
ekin = self.E() - self.m0()
return ekin
def gamma(self) : # return the Lorentz factor
return self.E() / self.m0()
def beta(self) : # return the reduced velocity
invgam2 = 1.0 / ( pow(self.gamma(),2) )
return math.sqrt( 1.0 - invgam2 )
def thet(self) :
return self.p.thet()
def phi(self) :
return self.p.phi()
#####
def __add__(self, x) :
rp = self.p + x.p
rW2 = self.W2 + 2*(self*x) + x.W2
return Vec4(rp,W2=rW2)
def __sub__(self, x) :
rp = self.p - x.p
rW2 = self.W2 - 2*(self*x) + x.W2
return Vec4(rp,W2=rW2)
def __mul__(self, x) :
return self.E()*x.E() - self.p * x.p
def __repr__(self) :
return '[Vec4: E %-11.5g, px %-11.5g, py %-11.5g, pz %-11.5g]' % (
self.E(), self.p.x, self.p.y, self.p.z)
# END of Vec4 class
# Vec3 class: implementation of three vectors
#
# J. A. Templon, Department of Physics and Astronomy, University of Georgia
#
# Usage: There are a number of methods to create a three-vector
#
# from vec3 import *
#
# your_threevec = Vec3(x=<val>,y=<val>,z=<val>)
# your_threevec = Vec3(mag=<val>,theta=<val>,phi=<val>)
# your_threevec = Vec3(cart=(<val>,<val>,<val>)
# your_threevec = Vec3(polar=(<val>,<val>,<val>)
#
# cart= specifies a cartesian-representation for the three vector
# (x, y, z)
# polar= specifies a spherical-polar representation for the three vector
# (mag, theta, phi)
#
# mag is the norm of the vector
# theta and phi should be specified in degrees
#
# if no arguments are specified, you get back a vector of length zero
# (x=y=z=0)
#
# attributes:
#
# vec3.x - the x component
# vec3.y - the y component
# vec3.z - the z component
#
# methods returning output:
#
# vec3.norm() - returns the norm (length)
# vec3.thet() - returns the polar angle of the momentum 3-vector (degrees)
# vec3.phi() - returns the azim. angle of the momentum 3-vector (degrees)
# vec3.hat() - returns unit 3-vector in direction of vec3
#
# operators:
#
# "-" (__sub__) - subtraction of two three-vectors
# "+" (__add__) - addition of two three-vectors
# "*" (__mul__) - multiplication of two three-vectors (dot product)
import math, string
DEGRAD = 180.0/math.pi # conversion radians -> degrees
class Vec3:
def __init__(self, # uses optional arguments to allow
x=0.0, y=0.0, z=0.0, # multiple definition styles
mag=0.0, theta=0.0, phi=0.0,
cart=(0.0, 0.0, 0.0),
polar=(0.0, 0.0, 0.0)) :
## figure out what the user specified. Make sure only one (or none)
## set was specified. If users specifies for example x, y, and theta,
## raise an error. If nothing is specified, return a zero vector.
cart1 = x <> 0.0 or y <> 0.0 or z <> 0.0 # logical variables
pol1 = mag <> 0.0 or theta <> 0.0 or phi <> 0.0
cart2 = cart <> (0.0, 0.0, 0.0)
pol2 = polar <> (0.0, 0.0, 0.0)
sum = cart1 + pol1 + cart2 + pol2 # number of sets (partially) passed
if sum < 1 : cart1 = 1 # nothing (or 0,0,0) was specified
## trap and raise when more than one spec method used.
if sum > 1 :
raise Vec3Error, 'Ambiguous three-vector specification'
##-> initialize the momentum vector and magnitude
if cart1 :
self.x = x ; self.y = y; self.z = z
elif pol1 :
thetarad = theta / DEGRAD
phirad = phi / DEGRAD
self.x = mag * math.sin(thetarad) * math.cos(phirad)
self.y = mag * math.sin(thetarad) * math.sin(phirad)
self.z = mag * math.cos(thetarad)
elif cart2 :
self.x, self.y, self.z = cart
elif pol2 :
thetarad = polar[1] / DEGRAD
phirad = polar[2] / DEGRAD
mag = polar[0]
self.x = mag * math.sin(thetarad) * math.cos(phirad)
self.y = mag * math.sin(thetarad) * math.sin(phirad)
self.z = mag * math.cos(thetarad)
else :
raise Vec3Error, 'Three-vector init: this can\'t happen!'
# end of __init__
def norm(self) : # returns the magnitude of the 3-momentum
p_mag = self.x*self.x + self.y*self.y + self.z*self.z
return math.sqrt(p_mag)
def thet(self) :
return DEGRAD * math.acos( self.z / self.norm() )
def phi(self) :
if self.x == 0 and self.y == 0 :
return 0.0
else :
return DEGRAD * math.atan2( self.y, self.x )
def hat(self) :
N = self.norm()
return Vec3(x=self.x/N, y=self.y/N, z=self.z/N)
def __add__(self, other) :
rx = self.x + other.x
ry = self.y + other.y
rz = self.z + other.z
return Vec3(cart=(rx, ry, rz))
def __sub__(self, other) :
rx = self.x - other.x
ry = self.y - other.y
rz = self.z - other.z
return Vec3(cart=(rx, ry, rz))
def __mul__(self, other) :
return (self.x * other.x) + (self.y * other.y) + (self.z * other.z)
# END of Vec3 class
# Module masstable
# masses for nuclear, subnuclear particles
# from Tuli's book (NNDC) 5th edition (1995)
# masses in MeV/c^2
mass = { }
mass['e'] = 0.51099906
mass['n'] = 939.56563
mass['p'] = 938.27231
mass['d'] = 1875.614
mass['t'] = 2808.922
mass['3He'] = 2808.392
mass['4He'] = 3727.380
This archive was generated by hypermail 2b29 : Tue Jan 04 2000 - 00:43:35 MET