|
1
|
! tshift_kaon.for
|
|
2
|
|
|
3
|
c this program uses kinematics for (e,e'K) L/T separation experiment
|
|
4
|
c to calculate the implied shift in t if the MM is wrong by some amount
|
|
5
|
c - keeping Q2, W, and meson CM-angle constant, calculate the change in
|
|
6
|
c meson kinematics implied by the MMshift
|
|
7
|
c - meson momentum and meson lab angle changed slightly
|
|
8
|
|
|
9
|
c Conventions:
|
|
10
|
c MMshift = difference between correct and expt values
|
|
11
|
c MMexpt=MMcorrect+MMshift (i.e. opposite sign of the MMshift applied to data)
|
|
12
|
c tEXP=tCORRECT+tSHIFT (i.e. apply the opposite sign of tSHIFT to data)
|
|
13
|
|
|
14
|
implicit none
|
|
15
|
|
|
16
|
real*8 mkp,mkpg,mL,mLg,mp,mpg,me,pi
|
|
17
|
real*8 nu,nug,pgam,pgamg,pgamx,pgamy,pgamz,thgam
|
|
18
|
real*8 tinc,tincg,einc,eincg,escat,escatg,tscat,tscatg
|
|
19
|
real*8 pinc,pincx,pincz,pscat,thscat,pscatx,pscatz
|
|
20
|
real*8 thkp0,thkpcm,ekp0,ekpcm0,pkp0,pkpcm0,pkpx0,pkpy0,pkpz0
|
|
21
|
real*8 thkp1,ekp1,ekpcm1,pkp1,pkpcm1,pkpx1,pkpy1,pkpz1
|
|
22
|
real*8 mmshift,mmshgev,tshift,tshiftg
|
|
23
|
real*8 pkpg,phix,thkplab0,thkplab1
|
|
24
|
real*8 betacm,gammacm,pgamtemp,pkptemp
|
|
25
|
real*8 ep1cm,pp1cm,egamcm,pgamcm,p1plus,gamplus,piplus
|
|
26
|
real*8 t0,t1,w,wg,q2g,ang
|
|
27
|
|
|
28
|
c stuff for linux_suppl.inc
|
|
29
|
real*8 sind, cosd, tand
|
|
30
|
real*8 asind, acosd, atand, atan2d
|
|
31
|
external sind, cosd, tand
|
|
32
|
external asind, acosd, atand, atan2d
|
|
33
|
|
|
34
|
me = 0.511
|
|
35
|
mkp = 493.677
|
|
36
|
mkpg = mkp/1000.
|
|
37
|
mp = 938.27
|
|
38
|
mpg = mp/1000.
|
|
39
|
mL = 1115.683
|
|
40
|
mLg = mL/1000.
|
|
41
|
|
|
42
|
pi=3.14156
|
|
43
|
|
|
44
|
1000 write(6,100)
|
|
45
|
100 format(' Enter Q^2, W in GeV')
|
|
46
|
read(5,*)q2g,wg
|
|
47
|
if (q2g.eq.0.0) goto 200
|
|
48
|
w = wg*1000.
|
|
49
|
|
|
50
|
write(6,101)
|
|
51
|
101 format(' Enter theta_pq CM (deg)')
|
|
52
|
read(5,*)thkpcm
|
|
53
|
c thkpcm=0
|
|
54
|
|
|
55
|
if (thkpcm.eq.0.) then
|
|
56
|
write(6,102)
|
|
57
|
102 format(/' Kaon will be emitted in parallel kinematics')
|
|
58
|
phix=0.
|
|
59
|
elseif (thkpcm.ne.0.) then
|
|
60
|
write(6,105)
|
|
61
|
105 format(' Enter PHI angle in degrees (0->360)')
|
|
62
|
read(5,*)phix
|
|
63
|
if (phix.lt.90. .or. phix.ge. 270.) write(6,103)
|
|
64
|
103 format(/' Kaon will be emitted forward of q vector')
|
|
65
|
if (phix.ge.90. .and. phix.lt. 270.) write(6,104)
|
|
66
|
104 format(/' Kaon will be emitted rearward of q vector')
|
|
67
|
endif
|
|
68
|
|
|
69
|
write(6,110)
|
|
70
|
110 format(' Enter the MM shift (MeV) >0 means data is high')
|
|
71
|
read(5,*)mmshift
|
|
72
|
mmshgev = mmshift/1000.
|
|
73
|
|
|
74
|
nug = (wg**2 + q2g - mpg**2)/(2.*mpg)
|
|
75
|
nu = nug * 1000.
|
|
76
|
|
|
77
|
pgamg = sqrt( nug**2 + q2g)
|
|
78
|
pgam = pgamg * 1000.
|
|
79
|
|
|
80
|
! find kaon lab momentum
|
|
81
|
! first, find speed of virtual photon+proton c.m. frame
|
|
82
|
betacm = pgam/(nu+mp)
|
|
83
|
gammacm = (nu+mp)/w
|
|
84
|
|
|
85
|
ekpcm0 = (w**2 + mkp**2 - (mL**2)) / (2.*w)
|
|
86
|
ekpcm1 = (w**2 + mkp**2 - ((mL+mmshift)**2)) / (2.*w)
|
|
87
|
pkpcm0 = sqrt( ekpcm0**2 - mkp**2 )
|
|
88
|
pkpcm1 = sqrt( ekpcm1**2 - mkp**2 )
|
|
89
|
c pkpcmg = pkpcm/1.e3
|
|
90
|
|
|
91
|
! now transform to lab frame wrt q-vector
|
|
92
|
ekp0 = gammacm* (betacm*pkpcm0*cosd(thkpcm) + ekpcm0)
|
|
93
|
ekp1 = gammacm* (betacm*pkpcm1*cosd(thkpcm) + ekpcm1)
|
|
94
|
pkp0 = sqrt(ekp0**2 - mkp**2)
|
|
95
|
pkp1 = sqrt(ekp1**2 - mkp**2)
|
|
96
|
! pkpg=pkp/1000.
|
|
97
|
thkplab0 = asind(sind(thkpcm)*pkpcm0/pkp0)
|
|
98
|
thkplab1 = asind(sind(thkpcm)*pkpcm1/pkp1)
|
|
99
|
|
|
100
|
write(6,*)' Beam energy (MeV) '
|
|
101
|
read(5,*)tinc
|
|
102
|
tincg=tinc/1000.
|
|
103
|
einc=tinc+me
|
|
104
|
eincg=einc/1.e3
|
|
105
|
pinc=sqrt(einc**2-me**2)
|
|
106
|
pincz=pinc
|
|
107
|
pincx=0.
|
|
108
|
|
|
109
|
escat=einc-nu
|
|
110
|
escatg=escat/1.e3
|
|
111
|
if (escat.lt.me) goto 200
|
|
112
|
tscat=escat-me
|
|
113
|
tscatg=tscat/1000.
|
|
114
|
pscat=sqrt(escat**2-me**2)
|
|
115
|
|
|
116
|
ang=(pinc**2+pscat**2-pgam**2)/2./pinc/pscat
|
|
117
|
if (ang.lt.-1. .or. ang.gt.1.) goto 200
|
|
118
|
thscat=acosd(ang)
|
|
119
|
pscatz = pscat*cosd(thscat)
|
|
120
|
pscatx = pscat*sind(thscat)
|
|
121
|
|
|
122
|
pgamz = pincz-pscatz
|
|
123
|
pgamx = pincx-pscatx
|
|
124
|
pgamy = 0.
|
|
125
|
|
|
126
|
pgamtemp = sqrt(pgamx**2 + pgamy**2 + pgamz**2)
|
|
127
|
if (abs(pgamtemp-pgam).gt.0.2) write(6,150)pgam,
|
|
128
|
* pgamtemp
|
|
129
|
150 format(' PGam disagreement ',2f10.3)
|
|
130
|
thgam = atand(pgamx/pgamz)
|
|
131
|
|
|
132
|
! now rotate to lab frame wrt e- beam
|
|
133
|
pkpx0 = pkp0*( cosd(thkplab0)*sind(thgam) -
|
|
134
|
* sind(thkplab0)*cosd(thgam)*cosd(phix) )
|
|
135
|
pkpy0 = pkp0*sind(thkplab0)*sind(phix)
|
|
136
|
pkpz0 = pkp0*( cosd(thkplab0)*cosd(thgam) +
|
|
137
|
* sind(thkplab0)*sind(thgam)*cosd(phix) )
|
|
138
|
pkpx1 = pkp1*( cosd(thkplab1)*sind(thgam) -
|
|
139
|
* sind(thkplab1)*cosd(thgam)*cosd(phix) )
|
|
140
|
pkpy1 = pkp1*sind(thkplab1)*sind(phix)
|
|
141
|
pkpz1 = pkp1*( cosd(thkplab1)*cosd(thgam) +
|
|
142
|
* sind(thkplab1)*sind(thgam)*cosd(phix) )
|
|
143
|
|
|
144
|
pkptemp = sqrt(pkpx0**2 + pkpy0**2 + pkpz0**2)
|
|
145
|
if (abs(pkptemp-pkp0).gt.0.2) write(6,160)pkp0,pkptemp
|
|
146
|
160 format(' Pkp disagreement ',2f10.2)
|
|
147
|
|
|
148
|
! general equation for t for non-parallel kinematics
|
|
149
|
! equation is actually for -t
|
|
150
|
t0 = (pgamx-pkpx0)**2 +(pgamy-pkpy0)**2
|
|
151
|
* +(pgamz-pkpz0)**2 -(nu-ekp0)**2
|
|
152
|
t1 = (pgamx-pkpx1)**2 +(pgamy-pkpy1)**2
|
|
153
|
* +(pgamz-pkpz1)**2 -(nu-ekp0)**2
|
|
154
|
! tg = t/1.e6
|
|
155
|
tshift=t1-t0
|
|
156
|
tshiftg=tshift/1.e6
|
|
157
|
write(6,165)mmshift,tshiftg
|
|
158
|
165 format(' MM shift of ',f5.3,' MeV gives ',f8.5,' GeV^2 shift')
|
|
159
|
|
|
160
|
! equation for u: between initial proton and outgoing kaon
|
|
161
|
! u = mp**2 + mkp**2 -2.*mp*ekp
|
|
162
|
! ug= u/1.e6
|
|
163
|
|
|
164
|
200 continue
|
|
165
|
end
|
|
166
|
|
|
167
|
include 'linux_suppl.inc'
|