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