|
1
|
! tshift_piff.for
|
|
2
|
|
|
3
|
c this program uses kinematics for (e,e'pi) 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 mpi,mpig,mn,mng,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 thpi0,thpicm,epi0,epicm0,ppi0,ppicm0,ppix0,ppiy0,ppiz0
|
|
21
|
real*8 thpi1,epi1,epicm1,ppi1,ppicm1,ppix1,ppiy1,ppiz1
|
|
22
|
real*8 mmshift,mmshgev,tshift,tshiftg
|
|
23
|
real*8 ppig,phix,thpilab0,thpilab1
|
|
24
|
real*8 betacm,gammacm,pgamtemp,ppitemp
|
|
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
|
mpi = 139.57
|
|
36
|
mpig = mpi/1000.
|
|
37
|
mp = 938.27
|
|
38
|
mpg = mp/1000.
|
|
39
|
mn = 939.57
|
|
40
|
mng = mn/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,*)thpicm
|
|
53
|
c thpicm=0
|
|
54
|
|
|
55
|
if (thpicm.eq.0.) then
|
|
56
|
write(6,102)
|
|
57
|
102 format(/' Pion will be emitted in parallel kinematics')
|
|
58
|
phix=0.
|
|
59
|
elseif (thpicm.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(/' Pion will be emitted forward of q vector')
|
|
65
|
if (phix.ge.90. .and. phix.lt. 270.) write(6,104)
|
|
66
|
104 format(/' Pion 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 pion 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
|
epicm0 = (w**2 + mpi**2 - (mn**2)) / (2.*w)
|
|
86
|
epicm1 = (w**2 + mpi**2 - ((mn+mmshift)**2)) / (2.*w)
|
|
87
|
ppicm0 = sqrt( epicm0**2 - mpi**2 )
|
|
88
|
ppicm1 = sqrt( epicm1**2 - mpi**2 )
|
|
89
|
c ppicmg = ppicm/1.e3
|
|
90
|
|
|
91
|
! now transform to lab frame wrt q-vector
|
|
92
|
epi0 = gammacm* (betacm*ppicm0*cosd(thpicm) + epicm0)
|
|
93
|
epi1 = gammacm* (betacm*ppicm1*cosd(thpicm) + epicm1)
|
|
94
|
ppi0 = sqrt(epi0**2 - mpi**2)
|
|
95
|
ppi1 = sqrt(epi1**2 - mpi**2)
|
|
96
|
! ppig=ppi/1000.
|
|
97
|
thpilab0 = asind(sind(thpicm)*ppicm0/ppi0)
|
|
98
|
thpilab1 = asind(sind(thpicm)*ppicm1/ppi1)
|
|
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
|
ppix0 = ppi0*( cosd(thpilab0)*sind(thgam) -
|
|
134
|
* sind(thpilab0)*cosd(thgam)*cosd(phix) )
|
|
135
|
ppiy0 = ppi0*sind(thpilab0)*sind(phix)
|
|
136
|
ppiz0 = ppi0*( cosd(thpilab0)*cosd(thgam) +
|
|
137
|
* sind(thpilab0)*sind(thgam)*cosd(phix) )
|
|
138
|
ppix1 = ppi1*( cosd(thpilab1)*sind(thgam) -
|
|
139
|
* sind(thpilab1)*cosd(thgam)*cosd(phix) )
|
|
140
|
ppiy1 = ppi1*sind(thpilab1)*sind(phix)
|
|
141
|
ppiz1 = ppi1*( cosd(thpilab1)*cosd(thgam) +
|
|
142
|
* sind(thpilab1)*sind(thgam)*cosd(phix) )
|
|
143
|
|
|
144
|
ppitemp = sqrt(ppix0**2 + ppiy0**2 + ppiz0**2)
|
|
145
|
if (abs(ppitemp-ppi0).gt.0.2) write(6,160)ppi0,ppitemp
|
|
146
|
160 format(' Ppi disagreement ',2f10.2)
|
|
147
|
|
|
148
|
! general equation for t for non-parallel kinematics
|
|
149
|
! equation is actually for -t
|
|
150
|
t0 = (pgamx-ppix0)**2 +(pgamy-ppiy0)**2
|
|
151
|
* +(pgamz-ppiz0)**2 -(nu-epi0)**2
|
|
152
|
t1 = (pgamx-ppix1)**2 +(pgamy-ppiy1)**2
|
|
153
|
* +(pgamz-ppiz1)**2 -(nu-epi0)**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 pion
|
|
161
|
! u = mp**2 + mpi**2 -2.*mp*epi
|
|
162
|
! ug= u/1.e6
|
|
163
|
|
|
164
|
200 continue
|
|
165
|
end
|
|
166
|
|
|
167
|
include 'linux_suppl.inc'
|