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). implicit none include 'simulate.inc' record /event_main/ main record /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.d-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.d-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.d0 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