!Context: Tracking
!
!
!
!Author: Th. Ullrich
!Last update: 11/5/95
!
! Definition of 3-momentum and 4-momentum data types
! ==================================================
!
! Definition of 3-momentum (px, py, pz) and 4-momentum (p, E)
! and related operations: *, -, +, ==, /=, = , abs()
!
! type(mom3) :: q has components p%px, p%py, p%pz
! type(mom4) :: r has components r%E, r%p%px, r%p%py, r%p%pz
!
!
! In addition the following functions (kind=double) are defined:
! --------------------------------------------------------------
! thetaof(mom3 q) returns polar angle in rad (0 for p=0)
! phiof(mom3 q) returns aziumthal angle in rad
! ptof(mom3 q) returns pt
! etaof(mom3 q) returns pseudorapidity (0 for p=0)
!
! yof(mom4 q) returns rapidity
! mtof(mom4 q) returns mt
! massof(mom4 q) returns mass (same as abs(q))
!
!See also: track_m
!
module mom_m
use f90_kind
!
! Data type definitions
!
type mom3
real(kind=double) :: px, py, pz
end type mom3
type mom4
real(kind=double) :: E
type(mom3) :: p
end type mom4
!
! Interfaces
!
interface operator (+)
module procedure mom3_plus_mom3, mom4_plus_mom4
end interface
interface operator (-)
module procedure mom3_minus_mom3, minus_mom3, &
mom4_minus_mom4, minus_mom4
end interface
interface operator (*)
module procedure mom3_times_mom3, mom4_times_mom4, &
scalar_times_mom3, mom3_times_scalar
end interface
interface operator (==)
module procedure mom3_eq_mom3, mom4_eq_mom4
end interface
interface operator (/=)
module procedure mom3_ne_mom3, mom4_ne_mom4
end interface
interface assignment (=)
module procedure mom3_ass_mom3, mom4_ass_mom4
end interface
interface abs
module procedure abs_mom3, abs_mom4
end interface
contains
function mom3_times_mom3(p1, p2)
implicit none
real(kind=double) :: mom3_times_mom3
type(mom3), intent(in) :: p1, p2
mom3_times_mom3 = p1%px*p2%px + p1%py*p2%py + p1%pz*p2%pz
end function mom3_times_mom3
function mom3_times_scalar(scale, p1)
implicit none
type(mom3) :: mom3_times_scalar
type(mom3), intent(in) :: p1
real(kind=double), intent(in) :: scale
mom3_times_scalar%px = scale*p1%px
mom3_times_scalar%py = scale*p1%py
mom3_times_scalar%pz = scale*p1%pz
end function mom3_times_scalar
function scalar_times_mom3(p1, scale)
implicit none
type(mom3) :: scalar_times_mom3
type(mom3), intent(in) :: p1
real(kind=double), intent(in) :: scale
scalar_times_mom3%px = scale*p1%px
scalar_times_mom3%py = scale*p1%py
scalar_times_mom3%pz = scale*p1%pz
end function scalar_times_mom3
function mom4_times_mom4(p1, p2)
implicit none
real(kind=double) :: mom4_times_mom4
type(mom4), intent(in) :: p1, p2
mom4_times_mom4 = p1%E*p2%E - p1%p*p2%p
if (mom4_times_mom4 < 0_double ) mom4_times_mom4 = 0_double
end function mom4_times_mom4
function mom3_plus_mom3(p1, p2)
implicit none
type(mom3) :: mom3_plus_mom3
type(mom3), intent(in) :: p1, p2
mom3_plus_mom3%px = p1%px + p2%px
mom3_plus_mom3%py = p1%py + p2%py
mom3_plus_mom3%pz = p1%pz + p2%pz
end function mom3_plus_mom3
function mom4_plus_mom4(p1, p2)
implicit none
type(mom4) :: mom4_plus_mom4
type(mom4), intent(in) :: p1, p2
mom4_plus_mom4%E = p1%E + p2%E
mom4_plus_mom4%p = p1%p + p2%p
end function mom4_plus_mom4
function mom3_minus_mom3(p1, p2)
implicit none
type(mom3) :: mom3_minus_mom3
type(mom3), intent(in) :: p1, p2
mom3_minus_mom3%px = p1%px - p2%px
mom3_minus_mom3%py = p1%py - p2%py
mom3_minus_mom3%pz = p1%pz - p2%pz
end function mom3_minus_mom3
function minus_mom3(p1)
implicit none
type(mom3) :: minus_mom3
type(mom3), intent(in) :: p1
minus_mom3%px = -p1%px
minus_mom3%py = -p1%py
minus_mom3%pz = -p1%pz
end function minus_mom3
function minus_mom4(p1)
implicit none
type(mom4) :: minus_mom4
type(mom4), intent(in) :: p1
minus_mom4%E = -p1%E
minus_mom4%p = -p1%p
end function minus_mom4
function mom4_minus_mom4(p1, p2)
implicit none
type(mom4) :: mom4_minus_mom4
type(mom4), intent(in) :: p1, p2
mom4_minus_mom4%E = p1%E - p2%E
mom4_minus_mom4%p = p1%p - p2%p
end function mom4_minus_mom4
function abs_mom4(p)
implicit none
real(kind=double) :: abs_mom4
type(mom4), intent(in) :: p
abs_mom4 = sqrt(p*p)
end function abs_mom4
function abs_mom3(p)
implicit none
real(kind=double) :: abs_mom3
type(mom3), intent(in) :: p
abs_mom3 = sqrt(p*p)
end function abs_mom3
function mom3_eq_mom3(p, q)
implicit none
logical :: mom3_eq_mom3
type(mom3), intent(in) :: p, q
mom3_eq_mom3 = ( p%px == q%px .and. &
p%py == q%py .and. p%pz == q%pz )
end function mom3_eq_mom3
function mom4_eq_mom4(p, q)
implicit none
logical :: mom4_eq_mom4
type(mom4), intent(in) :: p, q
mom4_eq_mom4 = ( p%E == q%E .and. p%p%px == q%p%px .and. &
p%p%py == q%p%py .and. p%p%pz == q%p%pz )
end function mom4_eq_mom4
function mom3_ne_mom3(p, q)
implicit none
logical :: mom3_ne_mom3
type(mom3), intent(in) :: p, q
mom3_ne_mom3 = ( p%px /= q%px .or. &
p%py /= q%py .or. p%pz /= q%pz )
end function mom3_ne_mom3
function mom4_ne_mom4(p, q)
implicit none
logical :: mom4_ne_mom4
type(mom4), intent(in) :: p, q
mom4_ne_mom4 = ( p%E /= q%E .or. p%p%px /= q%p%px .or. &
p%p%py /= q%p%py .or. p%p%pz /= q%p%pz )
end function mom4_ne_mom4
subroutine mom3_ass_mom3(p, q)
implicit none
type(mom3), intent(out) :: p
type(mom3), intent(in) :: q
p%px = q%px
p%py = q%py
p%pz = q%pz
end subroutine mom3_ass_mom3
subroutine mom4_ass_mom4(p, q)
implicit none
type(mom4), intent(out) :: p
type(mom4), intent(in) :: q
p%E = q%E
p%p = q%p
end subroutine mom4_ass_mom4
function ptof(p)
implicit none
type(mom3), intent(in) :: p
real(kind=double) :: ptof
ptof = sqrt(p%px*p%px+p%py*p%py)
end function ptof
function thetaof(p)
implicit none
type(mom3), intent(in) :: p
real(kind=double) :: thetaof, abs_p
logical, save :: init = .true.
real(kind=double), save :: pi
if (init) then
pi = acos(-1.d0)
init = .false.
endif
abs_p = abs(p)
if (abs_p > 0) then
thetaof = acos(p%pz/abs_p)
else
thetaof = 0_double
endif
end function thetaof
function phiof(p)
implicit none
type(mom3), intent(in) :: p
real(kind=double) :: phiof
logical, save :: init = .true.
real(kind=double), save :: pi
if (init) then
pi = acos(-1.d0)
init = .false.
endif
phiof = atan2(p%py,p%px)+(pi-sign(pi,p%py))
end function phiof
function etaof(p)
implicit none
type(mom3), intent(in) :: p
real(kind=double) :: etaof, theta_p
theta_p = thetaof(p)
if (theta_p > 0) then
etaof = -log(tan(theta_p/2_double))
else
etaof = 0
endif
end function etaof
function yof(p)
implicit none
type(mom4), intent(in) :: p
real(kind=double) :: yof
yof = 0.5*log((p%E+p%p%pz)/(p%E-p%p%pz))
end function yof
function mtof(p)
implicit none
type(mom4), intent(in) :: p
real(kind=double) :: mtof
mtof = sqrt(p%E*p%E-p%p%pz*p%p%pz)
end function mtof
function massof(p)
implicit none
type(mom4), intent(in) :: p
real(kind=double) :: massof
massof = abs(p)
end function massof
end module mom_m