1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
|
subroutine momntx(energy,mass,costh,phi , p)
c
c This subroutine sets up a four-momentum from the four inputs.
c
c input:
c real energy : energy
c real mass : mass
c real costh : cos(theta)
c real phi : azimuthal angle
c
c output:
c real p(0:3) : four-momentum
c
implicit none
double precision p(0:3),energy,mass,costh,phi,pp,sinth
double precision rZero, rOne
parameter( rZero = 0.0d0, rOne = 1.0d0 )
#ifdef HELAS_CHECK
double precision rPi, rTwo
parameter( rPi = 3.14159265358979323846d0, rTwo = 2.d0 )
integer stdo
parameter( stdo = 6 )
#endif
c
#ifdef HELAS_CHECK
if (energy.lt.mass) then
write(stdo,*)
& ' helas-error : energy in momntx is less than mass'
write(stdo,*)
& ' : energy = ',energy,' : mass = ',mass
endif
if (mass.lt.rZero) then
write(stdo,*) ' helas-error : mass in momntx is negative'
write(stdo,*) ' : mass = ',mass
endif
if (abs(costh).gt.rOne) then
write(stdo,*) ' helas-error : costh in momntx is improper'
write(stdo,*) ' : costh = ',costh
endif
if (phi.lt.rZero .or. phi.gt.rTwo*rPi) then
write(stdo,*)
& ' helas-warn : phi in momntx does not lie on 0.0 thru 2.0*pi'
write(stdo,*)
& ' : phi = ',phi
endif
#endif
p(0) = energy
if ( energy.eq.mass ) then
p(1) = rZero
p(2) = rZero
p(3) = rZero
else
pp = sqrt((energy-mass)*(energy+mass))
sinth = sqrt((rOne-costh)*(rOne+costh))
p(3) = pp*costh
if ( phi.eq.rZero ) then
p(1) = pp*sinth
p(2) = rZero
else
p(1) = pp*sinth*cos(phi)
p(2) = pp*sinth*sin(phi)
endif
endif
c
return
end
|