-
Notifications
You must be signed in to change notification settings - Fork 88
/
jacobians.f
282 lines (228 loc) · 10.7 KB
/
jacobians.f
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
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
subroutine transform_to_cm(vertex,main,
& gstar,bstar,bstarx,bstary,bstarz,
& nustar,qstar,qstarx,qstary,qstarz,
& ehadcm,phadcm,phadcmx,phadcmy,phadcmz,
& ebeamcm,pbeamcm,pbeamcmx,pbeamcmy,pbeamcmz,
& etarcm,ptarcm,ptarcmx,ptarcmy,ptarcmz,
& thetacm,phicm,phiqn,jacobian,jac_old)
* GENERATE_CM
* This routine determines the transformation from the lab variables to
* the PHOTON-NUCLEON CENTER-OF-MASS frame, and returns the center-of-mass
* four-vectors for the q-vector, beam, and produced hadron (pion or kaon).
USE structureModule
implicit none
include 'simulate.inc'
type(event_main):: main
type(event):: vertex
* Variables to be returned:
real*8 gstar,bstar,bstarx,bstary,bstarz !beta of boost toC.M.
real*8 nustar,qstar,qstarx,qstary,qstarz !q in C.M.
real*8 ehadcm,phadcm,phadcmx,phadcmy,phadcmz !p_hadron in C.M.
real*8 ebeamcm,pbeamcm,pbeamcmx,pbeamcmy,pbeamcmz !p_beam in C.M.
real*8 etarcm,ptarcm,ptarcmx,ptarcmy,ptarcmz !p_hadron in C.M.
real*8 thetacm,phicm,phiqn
real*8 jacobian,jac_old
* Temporary variables: Boost related
real*8 tcos,tsin !cos/sin of theta between ppi and q
real*8 tfcos,tfsin !cos/sin of theta between pfermi and q
real*8 cospq,sinpq
real*8 dummy,zero
real*8 qx,qy,qz,px,py,pz
real*8 tmp_x_x,tmp_x_y,tmp_x_z
real*8 tmp_y_x,tmp_y_y,tmp_y_z
real*8 p_tmp_x,p_tmp_y
real*8 pbeam,beam_tmpx,beam_tmpy,beam_tmpz !lab variables to be boosted.
real*8 phadx,phady,phadz ! " " " "
real*8 ptarx,ptary,ptarz ! " " " "
real*8 tmp2_x_x,tmp2_x_y,tmp2_x_z
real*8 tmp2_y_x,tmp2_y_y,tmp2_y_z
real*8 tmp2_z_x,tmp2_z_y,tmp2_z_z
real*8 phadcm_tmp2x,phadcm_tmp2y,phadcm_tmp2z
* Temporary variables: Jacobian related
real*8 square_root,dp_dcos_num,dp_dcos_den,dp_dcos
real*8 dp_dphi_num,dp_dphi_den,dp_dphi
real*8 dt_dcos_lab,dt_dphi_lab,psign
real*8 dpxdphi,dpydphi,dpxdcos,dpydcos,dpzdcos,dpzdphi
real*8 dpxnewdphi,dpynewdphi,dpxnewdcos,dpynewdcos
real*8 dpznewdphi,dpznewdcos
real*8 dphicmdphi,dphicmdcos
* f's and fer indicate fermi momenta, s, star or cm CM system
tcos = vertex%up%x*vertex%uq%x+vertex%up%y*vertex%uq%y+vertex%up%z*vertex%uq%z
if(tcos-1..gt.0..and.tcos-1..lt.1.e-8)tcos=1.0
tsin=sqrt(1.-tcos**2)
tfcos = pferx*vertex%uq%x+pfery*vertex%uq%y+pferz*vertex%uq%z
if(tfcos-1..gt.0..and.tfcos-1..lt.1.e-8)tfcos=1.0
tfsin=sqrt(1.-tfcos**2)
cospq = cos(main%phi_pq)
sinpq = sin(main%phi_pq)
* JRA:The transformation is calculated starting in the coord. system used
* in the fpi/nucpi replay (see event.f), where x=right, y=down, z=along beam.
* We must convert from SIMC coords to these first.
qx = -vertex%uq%y !'right'
qy = vertex%uq%x !'down'
qz = vertex%uq%z
px = -pfery
py = pferx
pz = pferz
dummy=sqrt((qx**2+qy**2)*(qx**2+qy**2+qz**2))
tmp_x_x = -qx*qz/dummy
tmp_x_y = -qy*qz/dummy
tmp_x_z = (qx**2 + qy**2)/dummy
dummy = sqrt(qx**2 + qy**2)
tmp_y_x = qy/dummy
tmp_y_y = -qx/dummy
tmp_y_z = 0.0
p_tmp_x = pfer*(px*tmp_x_x + py*tmp_x_y + pz*tmp_x_z)
p_tmp_y = pfer*(px*tmp_y_x + py*tmp_y_y + pz*tmp_y_z)
if(p_tmp_x.eq.0.)then
phiqn=0.
else
phiqn = atan2(p_tmp_y,p_tmp_x)
endif
if(phiqn.lt.0.)phiqn = phiqn+2.*pi
* get beam in "q" system.
pbeam = vertex%Ein
beam_tmpx = pbeam*tmp_x_z
beam_tmpy = pbeam*tmp_y_z
beam_tmpz = pbeam*vertex%uq%z
bstar=sqrt((vertex%q+pfer*tfcos)**2+(pfer*tfsin)**2)/(efer+vertex%nu)
if (abs(bstar).gt.1.) write(6,'(a,4f15.5)') 'q,nu,efer,beta=',vertex%q,vertex%nu,efer,bstar
! if (bstar.gt.1. .and. bstar.lt.1.001) bstar=0.9999999
! if (bstar.lt.-1. .and. bstar.gt.-1.001) bstar=-0.999999
gstar=1./sqrt(1. - bstar**2)
bstarz = (vertex%q+pfer*tfcos)/(efer+vertex%nu)
bstarx = p_tmp_x/(efer+vertex%nu)
bstary = p_tmp_y/(efer+vertex%nu)
* DJG Boost the beam to CM.
call loren(gstar,bstarx,bstary,bstarz,vertex%Ein,beam_tmpx,
> beam_tmpy,beam_tmpz,ebeamcm,pbeamcmx,pbeamcmy,pbeamcmz,pbeamcm)
* DJG: Boost virtual photon to CM.
zero =0.e0
call loren(gstar,bstarx,bstary,bstarz,vertex%nu,
> zero,zero,vertex%q,nustar,qstarx,qstary,qstarz,qstar)
* DJG: Boost pion to CM.
phadz = vertex%p%P*tcos
phadx = vertex%p%P*tsin*cospq
phady = vertex%p%P*tsin*sinpq
call loren(gstar,bstarx,bstary,bstarz,vertex%p%E,
> phadx,phady,phadz,ehadcm,phadcmx,phadcmy,phadcmz,phadcm)
thetacm = acos((phadcmx*qstarx+phadcmy*qstary+phadcmz*qstarz)/phadcm/qstar)
* DJG: Boost target nucleon to CM.
ptarz = pfer*tfcos
ptarx = p_tmp_x
ptary = p_tmp_y
call loren(gstar,bstarx,bstary,bstarz,efer,
> ptarx,ptary,ptarz,etarcm,ptarcmx,ptarcmy,ptarcmz,ptarcm)
* Thetacm is defined as angle between phadcm and qstar.
* To get phicm, we need out of plane angle relative to scattering plane
* (plane defined by pbeamcm and qcm). For stationary target, this plane
* does not change. In general, the new coordinate system is defined such
* that the new y direction is given by (qcm x pbeamcm) and the new x
* is given by (qcm x pbeamcm) x qcm.
dummy = sqrt((qstary*pbeamcmz-qstarz*pbeamcmy)**2+
> (qstarz*pbeamcmx-qstarx*pbeamcmz)**2
> +(qstarx*pbeamcmy-qstary*pbeamcmx)**2)
tmp2_y_x = (qstary*pbeamcmz-qstarz*pbeamcmy)/dummy
tmp2_y_y = (qstarz*pbeamcmx-qstarx*pbeamcmz)/dummy
tmp2_y_z = (qstarx*pbeamcmy-qstary*pbeamcmx)/dummy
dummy = sqrt((tmp2_y_y*qstarz-tmp2_y_z*qstary)**2
> +(tmp2_y_z*qstarx-tmp2_y_x*qstarz)**2
> +(tmp2_y_x*qstary-tmp2_y_y*qstarx)**2)
tmp2_x_x = (tmp2_y_y*qstarz-tmp2_y_z*qstary)/dummy
tmp2_x_y = (tmp2_y_z*qstarx-tmp2_y_x*qstarz)/dummy
tmp2_x_z = (tmp2_y_x*qstary-tmp2_y_y*qstarx)/dummy
tmp2_z_x = qstarx/qstar
tmp2_z_y = qstary/qstar
tmp2_z_z = qstarz/qstar
phadcm_tmp2x = phadcmx*tmp2_x_x + phadcmy*tmp2_x_y + phadcmz*tmp2_x_z
phadcm_tmp2y = phadcmx*tmp2_y_x + phadcmy*tmp2_y_y + phadcmz*tmp2_y_z
phadcm_tmp2z = phadcmx*tmp2_z_x + phadcmy*tmp2_z_y + phadcmz*tmp2_z_z
phicm = atan2(phadcm_tmp2y,phadcm_tmp2x)
if(phicm.lt.0.) phicm = 2.*pi+phicm
* While we're here, and have all of the above factors available, We'll
* calculate the Jacobian:
* Now, Henk's CM -> LAB transformation OLD VERSION OF THE JACOBIAN!!!!
* HPB: Brauel's cross section is expressed in an invariant way as dsigma/dtdphi,
* HPB: with a virtual photon flux factor based on a full cross section, written
* HPB: as: dsigma/dQ2dWdtdphi. One can convert this cross section to dsigma/dom+
* HPB: in the lab frame or the cm frame by multiplying with the factor q*p_pi/3+
* HPB: (with the variables in the lab frame, or the cm frame)
* HPB: this factor is built from dt/dcos(theta)=2pq, and a factor 1/2pi, which
* HPB: comes from the conversion of Brauel's virtual photon flux factor and the
* HPB: one we use, based on dsigma/dE_e'dOmega_e'dOmega_pi
* HPB: see Devenish and Lyth, Phys. Rev. C D5, 47(1972)
*
* JRA: Simpler (I hope) explanation: sig is d2sigma/dt/dphi_cm. To get
* d2sigma/domega_cm (=s2cm), we multiply by dt/domega_cm = 1/2/pi*dt/dcos(theta+
* = 1/pi*q_cm*p_cm. We can then get d5sigma/dE_e'/dOmega_e'/dOmega_pi_cm by
* multiplying by gammav. Either of these can be converted to the dOmega_lab
* by multiplying by dOmega_lab/dOmega_cm = dcos(theta_lab)/dcos(theta_cm)
* (dphi_cm/dphi_lab=1 since the frames are collinear). This is 'factor',
* which is equal to dt/dcos(theta_cm) / dt/dcos(theta_lab). Hence, the
* following cross sections can be calculated:
* s2cm = sig*qcm*ppicm /pi
* s5cm =gammav*sig*qcm*ppicm /pi
* s2lab= sig*q *ppi *factor/pi
* s5lab=gammav*sig*q *ppi *factor/pi
* DJG The above assumes target at rest. We need a full blown Jacobian:
* DJG: | dt/dcos(theta) dphicm/dphi - dt/dphi dphicm/dcos(theta) |
* DJG: Calculate dt/d(cos(theta)) and dt/dphi for the general case.
psign = cos(phiqn)*cospq+sin(phiqn)*sinpq
square_root = vertex%q + pfer*tfcos - vertex%p%P*tcos
dp_dcos_num = vertex%p%P + (vertex%p%P**2*tcos -
> psign*pfer*vertex%p%P*tfsin*tcos/tsin)/square_root
dp_dcos_den = ( (vertex%nu+efer-vertex%p%E)*vertex%p%P/vertex%p%E +
> vertex%p%P*tsin**2-psign*pfer*tfsin*tsin )/square_root - tcos
dp_dcos = dp_dcos_num/dp_dcos_den
dp_dphi_num = pfer*vertex%p%P*tsin*tfsin*(cos(phiqn)*sinpq-
> sin(phiqn)*cospq)/square_root
dp_dphi_den = tcos + (pfer*tsin*tfsin*psign - vertex%p%P*tsin**2
> - (vertex%nu+efer-vertex%p%E)*vertex%p%P/vertex%p%E)/square_root
dp_dphi = dp_dphi_num/dp_dphi_den
dt_dcos_lab = 2.*(vertex%q*vertex%p%P +
> (vertex%q*tcos-vertex%nu*vertex%p%P/vertex%p%E)*dp_dcos)
dt_dphi_lab = 2.*(vertex%q*tcos-vertex%nu*vertex%p%P/vertex%p%E)*dp_dphi
* DJG: Now calculate dphicm/dphi and dphicm/d(cos(theta)) in the most
* DJG: excruciating way possible.
dpxdphi = vertex%p%P*tsin*(-sinpq+(gstar-1.)*bstarx/bstar**2*
> (bstary*cospq-bstarx*sinpq) ) +
> ( (phadcmx+gstar*bstarx*vertex%p%E)/vertex%p%P -
> gstar*bstarx*vertex%p%P/vertex%p%E)*dp_dphi
dpydphi = vertex%p%P*tsin*(cospq+(gstar-1.)*bstary/bstar**2*
> (bstary*cospq-bstarx*sinpq) ) +
> ( (phadcmy+gstar*bstary*vertex%p%E)/vertex%p%P -
> gstar*bstary*vertex%p%P/vertex%p%E)*dp_dphi
dpzdphi = vertex%p%P*(gstar-1.)/bstar**2*bstarz*tsin*
> (bstary*cospq-bstarx*sinpq) +
> ((phadcmz+gstar*bstarz*vertex%p%E)/vertex%p%P-
> gstar*bstarz*vertex%p%P/vertex%p%E)*dp_dphi
dpxdcos = -vertex%p%P*tcos/tsin*(cospq+(gstar-1.)*bstarx/bstar**2*
> (bstarx*cospq+bstary*sinpq-bstarz*tsin/tcos)) +
> ( (phadcmx+gstar*bstarx*vertex%p%E)/vertex%p%P -
> gstar*bstarx*vertex%p%P/vertex%p%E)*dp_dcos
dpydcos = -vertex%p%P*tcos/tsin*(sinpq+(gstar-1.)*bstary/bstar**2*
> (bstarx*cospq+bstary*sinpq-bstarz*tsin/tcos)) +
> ( (phadcmy+gstar*bstary*vertex%p%E)/vertex%p%P -
> gstar*bstary*vertex%p%P/vertex%p%E)*dp_dcos
dpzdcos = vertex%p%P*(1.-(gstar-1.)/bstar**2*bstarz*tcos/tsin*
> (bstarx*cospq+bstary*sinpq-tsin/tcos*bstarz))
> +((phadcmz+gstar*bstarz*vertex%p%E)/vertex%p%P-gstar*bstarz*
> vertex%p%P/vertex%p%E)*dp_dcos
dpxnewdphi = dpxdphi*tmp2_x_x+dpydphi*tmp2_x_y+dpzdphi*tmp2_x_z
dpynewdphi = dpxdphi*tmp2_y_x+dpydphi*tmp2_y_y+dpzdphi*tmp2_y_z
dpznewdphi = dpxdphi*tmp2_z_x+dpydphi*tmp2_z_y+dpzdphi*tmp2_z_z
dphicmdphi = (dpynewdphi*phadcm_tmp2x - phadcm_tmp2y*dpxnewdphi)/
> (phadcm_tmp2x**2+phadcm_tmp2y**2)
dpxnewdcos = dpxdcos*tmp2_x_x+dpydcos*tmp2_x_y+dpzdcos*tmp2_x_z
dpynewdcos = dpxdcos*tmp2_y_x+dpydcos*tmp2_y_y+dpzdcos*tmp2_y_z
dpznewdcos = dpxdcos*tmp2_z_x+dpydcos*tmp2_z_y+dpzdcos*tmp2_z_z
dphicmdcos = (dpynewdcos*phadcm_tmp2x - phadcm_tmp2y*dpxnewdcos)
> /(phadcm_tmp2x**2+phadcm_tmp2y**2)
jacobian = abs(dt_dcos_lab*dphicmdphi-dt_dphi_lab*dphicmdcos)
* Old jacobian - assumes collinear boost (i.e. pfer along q vector).
jac_old = 2*(efer-2*pferz*pfer*vertex%p%E/vertex%p%P*tcos)*
> (vertex%q+pferz*pfer)*vertex%p%P /
> ( efer+vertex%nu-(vertex%q+pferz*pfer)*vertex%p%E/vertex%p%P*tcos )
> - 2*vertex%p%P*pfer
return
end