simc_semi/rho_physics.f - annotate - 1.2

Return to rho_physics.f CVS log

Up to [HallC] / simc_semi

1 gaskelld 1.1 real8 function peerho(vertex,main) 2 3 Purpose: 4 * This routine calculates p(e,e'rho)p cross sections using my attempt 5 * to code in the form used in PYTHIA with modifications implmented 6 * in the HERMES Monte Carlo by Patricia Liebing (Thanks Patty!) 7 * 8 * 9 * variables: 10 * input: 11 * omega !energy of virtual photon (MeV) 12 * theta_pq !angle between pi and q (rad) 13 * Q2_g !4-momentum of virtual photon, squared (GeV/c)^2 14 * epsilon !epsilon (dimensionless) 15 * 16 * output: 17 * sigma_eerho !d3sigma/dEe'dOmegae'Omegapi (microbarn/MeV/sr^2) 18 19 implicit none 20 include 'simulate.inc' 21 22 gaskelld 1.1 record /event_main/ main 23 record /event/ vertex 24 25 * NOTE: when we refer to the center of mass system, it always refers to the 26 * photon-NUCLEON center of mass, not the photon-NUCLEUS! The model gives 27 * the cross section in the photon-nucleon center of mass frame. 28 29 real8 sigma_eerho !final cross section (returned as peepi) 30 real8 sig219,sig,factor 31 real8 sigt !components of dsigma/dt 32 real8 s5lab !dsigma/dE_e/dOmega_e/dOmega_rho (in some lab) 33 c real8 efer !energy of target particle 34 real8 epsi !epsilon of virtual photon 35 real8 gtpr !gamma_t prime. 36 real8 tcos,tsin !cos/sin of theta between ppi and q 37 real8 tfcos,tfsin !cos/sin of theta between pfermi and q 38 real8 bstar,bstarx,bstary,bstarz,gstar !beta/gamma of C.M. system in lab 39 real8 ppix,ppiy,ppiz !pion momentum in lab. 40 real8 thetacm,phicm !C.M. scattering angles 41 real8 costcm,costo2cm 42 real8 sintcm,sinto2cm 43 gaskelld 1.1 real8 ppicm,ppicmx,ppicmy,ppicmz,epicm !pion E,p in C.M. 44 real8 qstar,qstarx,qstary,qstarz,nustar !c.m. photon momentum/energy 45 real8 zero 46 47 real8 s,t,Q2_g !t,s, Q2 in (GeV/c)*2 48 real8 tmin, tprime, cdeltatau,brho 49 real8 nu !equivilent photon energy (MeV) 50 51 real8 sig0,R 52 53 real8 new_x_x,new_x_y,new_x_z 54 real8 new_y_x,new_y_y,new_y_z 55 real8 new_z_x,new_z_y,new_z_z 56 real8 dummy,p_new_x,p_new_y,phiqn 57 real8 davesig,phipq,cospq,sinpq 58 real8 square_root,dp_dcos_num,dp_dcos_den,dp_dcos 59 real8 dp_dphi_num,dp_dphi_den,dp_dphi 60 real8 dt_dcos_lab,dt_dphi_lab,psign 61 real8 dpxdphi,dpydphi,dpxdcos,dpydcos,dpzdcos,dpzdphi 62 real8 dpxnewdphi,dpynewdphi,dpxnewdcos,dpynewdcos 63 real8 dpznewdphi,dpznewdcos 64 gaskelld 1.1 real8 dphicmdphi,dphicmdcos 65 real8 jacobian 66 67 real8 pbeam,beam_newx,beam_newy,beam_newz 68 real8 pbeamcmx,pbeamcmy,pbeamcmz,ebeamcm,pbeamcm 69 real8 ppicm_newx,ppicm_newy,ppicm_newz 70 71 real8 dEcmdcos,dEcmdphi,dcoscmdcos,dcoscmdphi 72 real8 qx,qy,qz,px,py,pz 73 74 75 * Initialize some stuff. 76 Q2_g = vertex.Q2/1000000. 77 phipq = main.phi_pq 78 cospq = cos(phipq) 79 sinpq = sin(phipq) 80 81 * calculate energy of initial (struck) nucleon, using the same assumptions that 82 * go into calculating the pion angle/momentum (in event.f). For A>1, the struck 83 * nucleon is off shell, the 2nd nucleon (always a neutron) is on shell, and has 84 * p = -p_fermi, and any additional nucleons are at rest. 85 gaskelld 1.1 86 efer = sqrt(pfer2+targ.Mtar_struck2) 87 if(doing_deutpi.or.doing_hepi) then 88 efer = targ.M-sqrt(Mn2+pfer2) 89 if(doing_hepi)efer=efer-mp 90 c mtar_offshell = sqrt(efer2-pfer*2) 91 endif 92 93 94 calculate some kinematical variables 95 * f's and fer indicate fermi momenta, s, star or cm CM system 96 97 tcos = vertex.up.xvertex.uq.x+vertex.up.yvertex.uq.y+vertex.up.zvertex.uq.z 98 if(tcos-1..gt.0..and.tcos-1..lt.1.d-8)tcos=1.0 99 tsin=sqrt(1.-tcos2) 100 101 tfcos = pferxvertex.uq.x+pferyvertex.uq.y+pferzvertex.uq.z 102 if(tfcos-1..gt.0..and.tfcos-1..lt.1.d-8)tfcos=1.0 103 tfsin=sqrt(1.-tfcos*2) 104 105 epsi = 1./(1.+2(1.+vertex.nu*2/vertex.Q2)tan(vertex.e.theta/2.)2) 106 gaskelld 1.1 107 s = (vertex.nu+efer)2-(vertex.q+pfertfcos)2-(pfertfsin)2 108 nu = (s-(efer2-pfer*2))/2./(efer-pfertfcos) !equiv pho energy(MeV) 109 110 s = s/1.d6 !CONVERT TO (GeV)*2 111 main.w = sqrt(s) 112 113 t = vertex.Q2-Mrho2+2.vertex.nuvertex.p.E-2.vertex.p.Pvertex.qtcos 114 t = t/1.d6 !CONVERT TO (GeV/c)*2 115 main.t = t 116 117 118 Calculate velocity of PHOTON-NUCLEON C.M. system in the lab frame. Use beta 119 * and gamma of the cm system (bstar and gstar) to transform particles into 120 * c.m. frame. Define z along the direction of q, and x to be along the 121 * direction of the pion momentum perpendicular to q. 122 * DJG: Get pfer components in the lab "q" system 123 124 * JRA:The transformation is calculated starting in the coord. system used 125 * in the fpi/nucpi replay (see event.f), where x=right, y=down, z=along beam. 126 * We must convert from SIMC coords to these first. 127 gaskelld 1.1 qx = -vertex.uq.y !'right' 128 qy = vertex.uq.x !'down' 129 qz = vertex.uq.z 130 px = -pfery 131 py = pferx 132 pz = pferz 133 134 dummy=sqrt((qx2+qy2)(qx2+qy2+qz2)) 135 new_x_x = -qxqz/dummy 136 new_x_y = -qyqz/dummy 137 new_x_z = (qx2 + qy2)/dummy 138 139 dummy = sqrt(qx2 + qy2) 140 new_y_x = qy/dummy 141 new_y_y = -qx/dummy 142 new_y_z = 0.0 143 144 p_new_x = pfer(pxnew_x_x + pynew_x_y + pznew_x_z) 145 p_new_y = pfer(pxnew_y_x + pynew_y_y + pznew_y_z) 146 147 if(p_new_x.eq.0.)then 148 gaskelld 1.1 phiqn=0. 149 else 150 phiqn = atan2(p_new_y,p_new_x) 151 endif 152 if(phiqn.lt.0.)phiqn = phiqn+2.pi 153 154 * get beam in "q" system. 155 156 pbeam = sqrt(vertex.Ein2-me2) 157 beam_newx = pbeamnew_x_z 158 beam_newy = pbeamnew_y_z 159 beam_newz = pbeamvertex.uq.z 160 161 bstar=sqrt((vertex.q+pfertfcos)*2+(pfertfsin)2)/(efer+vertex.nu) 162 gstar=1./sqrt(1. - bstar2) 163 164 bstarz = (vertex.q+pfertfcos)/(efer+vertex.nu) 165 bstarx = p_new_x/(efer+vertex.nu) 166 bstary = p_new_y/(efer+vertex.nu) 167 168 DJG: Boost virtual photon to CM. 169 gaskelld 1.1 170 zero =0.d0 171 call loren(gstar,bstarx,bstary,bstarz,vertex.nu, 172 > zero,zero,vertex.q,nustar,qstarx,qstary,qstarz,qstar) 173 174 * DJG: Boost pion to CM. 175 176 ppiz = vertex.p.Ptcos 177 ppix = vertex.p.Ptsincospq 178 ppiy = vertex.p.Ptsinsinpq 179 call loren(gstar,bstarx,bstary,bstarz,vertex.p.E, 180 > ppix,ppiy,ppiz,epicm,ppicmx,ppicmy,ppicmz,ppicm) 181 thetacm = acos((ppicmxqstarx+ppicmyqstary+ppicmzqstarz)/ppicm/qstar) 182 main.pcm = ppicm 183 184 * DJG Boost the beam to CM. 185 186 call loren(gstar,bstarx,bstary,bstarz,vertex.Ein,beam_newx, 187 > beam_newy,beam_newz,ebeamcm,pbeamcmx,pbeamcmy,pbeamcmz,pbeamcm) 188 189 * Thetacm is defined as angle between ppicm and qstar. 190 gaskelld 1.1 * To get phicm, we need out of plane angle relative to scattering plane 191 * (plane defined by pbeamcm and qcm). For stationary target, this plane 192 * does not change. In general, the new coordinate system is defined such 193 * that the new y direction is given by (qcm x pbeamcm) and the new x 194 * is given by (qcm x pbeamcm) x qcm. 195 196 dummy = sqrt((qstarypbeamcmz-qstarzpbeamcmy)*2+ 197 > (qstarzpbeamcmx-qstarxpbeamcmz)2 198 > +(qstarxpbeamcmy-qstarypbeamcmx)2) 199 new_y_x = (qstarypbeamcmz-qstarzpbeamcmy)/dummy 200 new_y_y = (qstarzpbeamcmx-qstarxpbeamcmz)/dummy 201 new_y_z = (qstarxpbeamcmy-qstarypbeamcmx)/dummy 202 203 dummy = sqrt((new_y_yqstarz-new_y_zqstary)2 204 > +(new_y_zqstarx-new_y_xqstarz)2 205 > +(new_y_xqstary-new_y_yqstarx)2) 206 new_x_x = (new_y_yqstarz-new_y_zqstary)/dummy 207 new_x_y = (new_y_zqstarx-new_y_xqstarz)/dummy 208 new_x_z = (new_y_xqstary-new_y_yqstarx)/dummy 209 210 new_z_x = qstarx/qstar 211 gaskelld 1.1 new_z_y = qstary/qstar 212 new_z_z = qstarz/qstar 213 214 ppicm_newx = ppicmxnew_x_x + ppicmynew_x_y + ppicmznew_x_z 215 ppicm_newy = ppicmxnew_y_x + ppicmynew_y_y + ppicmznew_y_z 216 ppicm_newz = ppicmxnew_z_x + ppicmynew_z_y + ppicmznew_z_z 217 218 sintcm = sin(thetacm) 219 costcm = cos(thetacm) 220 sinto2cm = sin(thetacm/2.) 221 costo2cm = cos(thetacm/2.) 222 phicm = atan2(ppicm_newy,ppicm_newx) 223 if(phicm.lt.0.) phicm = 2.3.141592654+phicm 224 225 main.thetacm = thetacm 226 main.phicm = phicm 227 228 ! DJG need tprime and tmin here 229 ! DJG Need an overall minus sign since we calculate -t above. 230 231 tmin = -( ((-Q2_g-mrho2/1.d6-(targ.mtar_struck/1000.)2+ 232 gaskelld 1.1 1 (targ.mtar_struck/1000.)2)/(2.sqrt(s)))2- 233 2 ((qstar-ppicm)/1000.)2 ) 234 235 236 tprime = abs(t-tmin) 237 ! Put in some W dependence from photoproduction data
238 gaskelld 1.2 ! DJG: This is my rough fit to some old photoproduction data.
239 gaskelld 1.1 240 sig0 = 41.263/(vertex.nu/1000.0)*0.4765 ! microbarns 241 242 ! DJG: R is usually fit to the form c_0 (Q2/M2_rho)^c1 243 ! DJG: The c_0 and c_1 are taken from HERMES data. 244 ! DJG: WARNING!!!! If you change this parameterization, you really 245 ! DJG: should change it in rho_decay also if you want to be self 246 ! DJG: consistent. 247 248 R = 0.33(vertex.Q2/mrho2)*(0.61) 249 250 ! DJG: The Q2 dependence is usually given by (M2_rho/(Q2+M2_rho))^2 251 ! DJG: HERMES found that 2.575 works better than 2 in the exponent. 252 253 sigt = sig0(1.0+epsiR)(mrho2/(vertex.Q2+mrho2))(2.575) 254 255 ! DJG: Need to parameterize t-dependence with b parameter as a function of c-tau 256 257 cdeltatau = hbarc/(sqrt(vertex.nu2+vertex.Q2+mrho2)-vertex.nu) !in fm! 258 if(cdeltatau.lt.2.0) then 259 brho = 4.4679 + 8.6106dlog10(cdeltatau) 260 gaskelld 1.1 else 261 brho = 7.0 262 endif 263 264 sig219 = sigtbrhoexp(-brhotprime)/2.0/pi !ub/GeV2/rad 265 266 sig=sig219/1.d+06 !dsig/dtdphicm in microbarns/MeV2/rad 267
268 gaskelld 1.2 269 C DJG Convert to dsig/dOmega_cm using dt/d(costhetacm) = 2 qcm pcm 270 271 sig = sig2.qstar*ppicm 272
273 gaskelld 1.1 ******************************************************************************* 274 * GMH: Virtual photon to electron beam flux conversion factor 275 gtpr = alpha/2./(pi*2)vertex.e.E/vertex.Ein(s-(targ.Mtar_struck/1000.)2)/2./ 276 > ((efer-pfertfcos)/1000.)/Q2_g/(1.-epsi) 277 278 279 factor=targ.Mtar_struck/(targ.Mtar_struck+vertex.nu-vertex.qvertex.p.E/vertex.p.Ptcos) 280 s5lab = gtprsigvertex.qvertex.p.Pfactor/pi/1.d+06 281 282 * DJG The above assumes target at rest. We need a full blown Jacobian: 283 * DJG \| dt/dcos(theta) dphicm/dphi - dt/dphi dphicm/dcos(theta) \| 284 * DJG: Calculate dt/d(cos(theta)) and dt/dphi for the general case. 285 286 psign = cos(phiqn)cospq+sin(phiqn)sinpq 287 288 square_root = vertex.q + pfertfcos - vertex.p.Ptcos 289 dp_dcos_num = vertex.p.P + (vertex.p.P*2tcos - 290 > psignpfervertex.p.Ptfsintcos/tsin)/square_root 291 dp_dcos_den = ( (vertex.nu+efer-vertex.p.E)vertex.p.P/vertex.p.E + 292 > vertex.p.Ptsin*2-psignpfertfsintsin )/square_root - tcos 293 dp_dcos = dp_dcos_num/dp_dcos_den 294 gaskelld 1.1 295 dp_dphi_num = pfervertex.p.Ptsintfsin(cos(phiqn)sinpq- 296 > sin(phiqn)cospq)/square_root 297 dp_dphi_den = tcos + (pfertsintfsinpsign - vertex.p.Ptsin*2 298 > - (vertex.nu+efer-vertex.p.E)vertex.p.P/vertex.p.E)/square_root 299 dp_dphi = dp_dphi_num/dp_dphi_den 300 301 dt_dcos_lab = 2.(vertex.qvertex.p.P + 302 > (vertex.qtcos-vertex.nuvertex.p.P/vertex.p.E)dp_dcos) 303 304 dt_dphi_lab = 2.(vertex.qtcos-vertex.nuvertex.p.P/vertex.p.E)dp_dphi 305 306 DJG: Now calculate dphicm/dphi and dphicm/d(cos(theta)) in the most 307 * DJG: excruciating way possible. 308 309 dpxdphi = vertex.p.Ptsin(-sinpq+(gstar-1.)bstarx/bstar2 310 > (bstarycospq-bstarxsinpq) ) + 311 > ( (ppicmx+gstarbstarxvertex.p.E)/vertex.p.P - 312 > gstarbstarxvertex.p.P/vertex.p.E)dp_dphi 313 dpydphi = vertex.p.Ptsin(cospq+(gstar-1.)bstary/bstar*2 314 > (bstarycospq-bstarxsinpq) ) + 315 gaskelld 1.1 > ( (ppicmy+gstarbstaryvertex.p.E)/vertex.p.P - 316 > gstarbstaryvertex.p.P/vertex.p.E)dp_dphi 317 dpzdphi = vertex.p.P(gstar-1.)/bstar*2bstarztsin 318 > (bstarycospq-bstarxsinpq) + 319 > ((ppicmz+gstarbstarzvertex.p.E)/vertex.p.P- 320 > gstarbstarzvertex.p.P/vertex.p.E)dp_dphi 321 322 dpxdcos = -vertex.p.Ptcos/tsin(cospq+(gstar-1.)bstarx/bstar*2 323 > (bstarxcospq+bstarysinpq-bstarztsin/tcos)) + 324 > ( (ppicmx+gstarbstarxvertex.p.E)/vertex.p.P - 325 > gstarbstarxvertex.p.P/vertex.p.E)dp_dcos 326 dpydcos = -vertex.p.Ptcos/tsin(sinpq+(gstar-1.)bstary/bstar2 327 > (bstarxcospq+bstarysinpq-bstarztsin/tcos)) + 328 > ( (ppicmy+gstarbstaryvertex.p.E)/vertex.p.P - 329 > gstarbstaryvertex.p.P/vertex.p.E)dp_dcos 330 dpzdcos = vertex.p.P(1.-(gstar-1.)/bstar2bstarztcos/tsin 331 > (bstarxcospq+bstarysinpq-tsin/tcosbstarz)) 332 > +((ppicmz+gstarbstarzvertex.p.E)/vertex.p.P-gstarbstarz* 333 > vertex.p.P/vertex.p.E)dp_dcos 334 335 dpxnewdphi = dpxdphinew_x_x+dpydphinew_x_y+dpzdphinew_x_z 336 gaskelld 1.1 dpynewdphi = dpxdphinew_y_x+dpydphinew_y_y+dpzdphinew_y_z 337 dpznewdphi = dpxdphinew_z_x+dpydphinew_z_y+dpzdphinew_z_z 338 339 dphicmdphi = (dpynewdphippicm_newx - ppicm_newydpxnewdphi)/ 340 > (ppicm_newx2+ppicm_newy2) 341 342 dpxnewdcos = dpxdcosnew_x_x+dpydcosnew_x_y+dpzdcosnew_x_z 343 dpynewdcos = dpxdcosnew_y_x+dpydcosnew_y_y+dpzdcosnew_y_z 344 dpznewdcos = dpxdcosnew_z_x+dpydcosnew_z_y+dpzdcosnew_z_z 345 346 dphicmdcos = (dpynewdcosppicm_newx - ppicm_newydpxnewdcos) 347 > /(ppicm_newx2+ppicm_newy2) 348 349 dcoscmdcos = dpznewdcos/ppicm-ppicm_newzepicmdEcmdcos/ppicm1.5 350 dcoscmdphi = dpznewdphi/ppicm-ppicm_newzepicmdEcmdphi/ppicm1.5 351 352 jacobian = abs(dt_dcos_labdphicmdphi-dt_dphi_labdphicmdcos) 353 main.davejac = jacobian 354 355 main.johnjac = 2(efer-2pferzpfervertex.p.E/vertex.p.Ptcos)* 356 > (vertex.q+pferzpfer)vertex.p.P / 357 gaskelld 1.1 > ( efer+vertex.nu-(vertex.q+pferzpfer)vertex.p.E/vertex.p.Ptcos ) 358 > - 2vertex.p.P*pfer 359
360 gaskelld 1.2 c davesig = gtprsigjacobian 361 362 davesig = gtpr*sig
363 gaskelld 1.1 sigma_eerho = davesig/1.d3 !ub/GeV-sr --> ub/MeV-sr 364 peerho = sigma_eerho 365 ntup.sigcm = sig !sig_cm 366 367 202 format(/11X,f5.1/) 368 203 format(11X,f5.0) 369 204 format(6(/9X,7f8.3)) 370 371 return 372 end

Analyzer/Replay: Mark Jones, Documents: Stephen Wood

Powered by
ViewCVS 0.9.2-cvsgraph-1.4.0