Return to rho_physics.f CVS log

Up to [HallC] / simc_semi

Changes for rho0 generation

	real*8 function peerho(vertex,main)

* Purpose:
* This routine calculates p(e,e'rho)p cross sections using my attempt
* to code in the form used in PYTHIA with modifications implmented
* in the HERMES Monte Carlo by Patricia Liebing (Thanks Patty!)
*
*
* variables:
*   input:
*	omega		!energy of virtual photon		(MeV)
*	theta_pq	!angle between pi and q			(rad)
*	Q2_g            !4-momentum of virtual photon, squared  (GeV/c)^2
*	epsilon		!epsilon				(dimensionless)
*
*   output:
*	sigma_eerho	!d3sigma/dEe'dOmegae'Omegapi	(microbarn/MeV/sr^2)

	implicit none
	include 'simulate.inc'

	record /event_main/ main
	record /event/	vertex

* NOTE: when we refer to the center of mass system, it always refers to the
* photon-NUCLEON center of mass, not the photon-NUCLEUS!  The model gives
* the cross section in the photon-nucleon center of mass frame.

	real*8 sigma_eerho		!final cross section (returned as peepi)
	real*8 sig219,sig,factor
	real*8 sigt		        !components of dsigma/dt
	real*8 s5lab			!dsigma/dE_e/dOmega_e/dOmega_rho (in some lab)
c	real*8 efer			!energy of target particle
	real*8 epsi			!epsilon of virtual photon
	real*8 gtpr			!gamma_t prime.
	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 bstar,bstarx,bstary,bstarz,gstar	!beta/gamma of C.M. system in lab
	real*8 ppix,ppiy,ppiz		!pion momentum in lab.
	real*8 thetacm,phicm		!C.M. scattering angles
	real*8 costcm,costo2cm
	real*8 sintcm,sinto2cm
	real*8 ppicm,ppicmx,ppicmy,ppicmz,epicm	 !pion E,p in C.M.
	real*8 qstar,qstarx,qstary,qstarz,nustar !c.m. photon momentum/energy
	real*8 zero

	real*8 s,t,Q2_g			!t,s, Q2 in (GeV/c)**2
	real*8 tmin, tprime, cdeltatau,brho
	real*8 nu			!equivilent photon energy (MeV)

	real*8 sig0,R

	real*8 new_x_x,new_x_y,new_x_z
	real*8 new_y_x,new_y_y,new_y_z
	real*8 new_z_x,new_z_y,new_z_z
	real*8 dummy,p_new_x,p_new_y,phiqn
	real*8 davesig,phipq,cospq,sinpq
	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
	real*8 jacobian

	real*8 pbeam,beam_newx,beam_newy,beam_newz
	real*8 pbeamcmx,pbeamcmy,pbeamcmz,ebeamcm,pbeamcm
	real*8 ppicm_newx,ppicm_newy,ppicm_newz

	real*8 dEcmdcos,dEcmdphi,dcoscmdcos,dcoscmdphi
	real*8 qx,qy,qz,px,py,pz


* Initialize some stuff.
	Q2_g = vertex.Q2/1000000.
	phipq = main.phi_pq
	cospq = cos(phipq)
	sinpq = sin(phipq)

* calculate energy of initial (struck) nucleon, using the same assumptions that
* go into calculating the pion angle/momentum (in event.f).  For A>1, the struck
* nucleon is off shell, the 2nd nucleon (always a neutron) is on shell, and has
* p = -p_fermi, and any additional nucleons are at rest.

	efer = sqrt(pfer**2+targ.Mtar_struck**2)
	if(doing_deutpi.or.doing_hepi) then
	  efer = targ.M-sqrt(Mn2+pfer**2)
	  if(doing_hepi)efer=efer-mp
c	  mtar_offshell = sqrt(efer**2-pfer**2)
	endif


* calculate some kinematical variables
* 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)

	epsi = 1./(1.+2*(1.+vertex.nu**2/vertex.Q2)*tan(vertex.e.theta/2.)**2)

	s = (vertex.nu+efer)**2-(vertex.q+pfer*tfcos)**2-(pfer*tfsin)**2
	nu = (s-(efer**2-pfer**2))/2./(efer-pfer*tfcos) !equiv pho energy(MeV)

	s = s/1.d6				!CONVERT TO (GeV)**2
	main.w = sqrt(s)

	t = vertex.Q2-Mrho2+2.*vertex.nu*vertex.p.E-2.*vertex.p.P*vertex.q*tcos
	t = t/1.d6				!CONVERT TO (GeV/c)**2
	main.t = t


* Calculate velocity of PHOTON-NUCLEON C.M. system in the lab frame. Use beta
* and gamma of the cm system (bstar and gstar) to transform particles into
* c.m. frame.  Define z along the direction of q, and x to be along the
* direction of the pion momentum perpendicular to q.
* DJG: Get pfer components in the lab "q" system

* 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))
	new_x_x = -qx*qz/dummy
	new_x_y = -qy*qz/dummy
	new_x_z = (qx**2 + qy**2)/dummy

	dummy   = sqrt(qx**2 + qy**2)
	new_y_x =  qy/dummy
	new_y_y = -qx/dummy
	new_y_z =  0.0

	p_new_x = pfer*(px*new_x_x + py*new_x_y + pz*new_x_z)
	p_new_y = pfer*(px*new_y_x + py*new_y_y + pz*new_y_z)

	if(p_new_x.eq.0.)then
	  phiqn=0.
	else
	  phiqn = atan2(p_new_y,p_new_x)
	endif
	if(phiqn.lt.0.)phiqn = phiqn+2.*pi

* get beam in "q" system.

	pbeam = sqrt(vertex.Ein**2-me**2)
	beam_newx = pbeam*new_x_z
	beam_newy = pbeam*new_y_z
	beam_newz = pbeam*vertex.uq.z

	bstar=sqrt((vertex.q+pfer*tfcos)**2+(pfer*tfsin)**2)/(efer+vertex.nu)
	gstar=1./sqrt(1. - bstar**2)

	bstarz = (vertex.q+pfer*tfcos)/(efer+vertex.nu)
	bstarx = p_new_x/(efer+vertex.nu)
	bstary = p_new_y/(efer+vertex.nu)

* 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.
	
	ppiz = vertex.p.P*tcos
	ppix = vertex.p.P*tsin*cospq
	ppiy = vertex.p.P*tsin*sinpq
	call loren(gstar,bstarx,bstary,bstarz,vertex.p.E,
     >		ppix,ppiy,ppiz,epicm,ppicmx,ppicmy,ppicmz,ppicm)
	thetacm = acos((ppicmx*qstarx+ppicmy*qstary+ppicmz*qstarz)/ppicm/qstar)
	main.pcm = ppicm

* DJG Boost the beam to CM.

	call loren(gstar,bstarx,bstary,bstarz,vertex.Ein,beam_newx,
     >		beam_newy,beam_newz,ebeamcm,pbeamcmx,pbeamcmy,pbeamcmz,pbeamcm)

* Thetacm is defined as angle between ppicm 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)
	new_y_x = (qstary*pbeamcmz-qstarz*pbeamcmy)/dummy
	new_y_y = (qstarz*pbeamcmx-qstarx*pbeamcmz)/dummy
	new_y_z = (qstarx*pbeamcmy-qstary*pbeamcmx)/dummy

	dummy = sqrt((new_y_y*qstarz-new_y_z*qstary)**2
     >		+(new_y_z*qstarx-new_y_x*qstarz)**2
     >		+(new_y_x*qstary-new_y_y*qstarx)**2)
	new_x_x = (new_y_y*qstarz-new_y_z*qstary)/dummy
	new_x_y = (new_y_z*qstarx-new_y_x*qstarz)/dummy
	new_x_z = (new_y_x*qstary-new_y_y*qstarx)/dummy

	new_z_x = qstarx/qstar
	new_z_y = qstary/qstar
	new_z_z = qstarz/qstar

	ppicm_newx = ppicmx*new_x_x + ppicmy*new_x_y + ppicmz*new_x_z
	ppicm_newy = ppicmx*new_y_x + ppicmy*new_y_y + ppicmz*new_y_z
	ppicm_newz = ppicmx*new_z_x + ppicmy*new_z_y + ppicmz*new_z_z

	sintcm = sin(thetacm)
	costcm = cos(thetacm)
	sinto2cm = sin(thetacm/2.)
	costo2cm = cos(thetacm/2.)
	phicm = atan2(ppicm_newy,ppicm_newx)
	if(phicm.lt.0.) phicm = 2.*3.141592654+phicm

	main.thetacm = thetacm
	main.phicm = phicm

! DJG need tprime and tmin here
! DJG Need an overall minus sign since we calculate -t above.

	tmin = -( ((-Q2_g-mrho2/1.d6-(targ.mtar_struck/1000.)**2+
	1    (targ.mtar_struck/1000.)**2)/(2.*sqrt(s)))**2-
	2    ((qstar-ppicm)/1000.)**2 )


	tprime = abs(t-tmin)
! Put in some W dependence from photoproduction data
! DJG:  This is my rough fit to some old photoproduction data.	

	sig0 = 41.263/(vertex.nu/1000.0)**0.4765   ! microbarns

! DJG:  R is usually fit to the form c_0 (Q2/M2_rho)^c1
! DJG:  The c_0 and c_1 are taken from HERMES data.
! DJG:  WARNING!!!! If you change this parameterization, you really
! DJG:  should change it in rho_decay also if you want to be self
! DJG:  consistent.

	R = 0.33*(vertex.Q2/mrho2)**(0.61)

! DJG:  The Q2 dependence is usually given by (M2_rho/(Q2+M2_rho))^2
! DJG:  HERMES found that 2.575 works better than 2 in the exponent.

	sigt = sig0*(1.0+epsi*R)*(mrho2/(vertex.Q2+mrho2))**(2.575)

! DJG:  Need to parameterize t-dependence with b parameter as a function of c-tau

	cdeltatau = hbarc/(sqrt(vertex.nu**2+vertex.Q2+mrho2)-vertex.nu) !in fm!
	if(cdeltatau.lt.2.0) then
	   brho = 4.4679 + 8.6106*dlog10(cdeltatau)
	else
	   brho = 7.0
	endif

	sig219 = sigt*brho*exp(-brho*tprime)/2.0/pi !ub/GeV**2/rad

	sig=sig219/1.d+06	!dsig/dtdphicm in microbarns/MeV**2/rad


C DJG Convert to dsig/dOmega_cm using dt/d(costhetacm) = 2 qcm pcm

	sig = sig*2.*qstar*ppicm

*******************************************************************************
* GMH: Virtual photon to electron beam flux conversion factor
	gtpr = alpha/2./(pi**2)*vertex.e.E/vertex.Ein*(s-(targ.Mtar_struck/1000.)**2)/2./
     >		((efer-pfer*tfcos)/1000.)/Q2_g/(1.-epsi)


	factor=targ.Mtar_struck/(targ.Mtar_struck+vertex.nu-vertex.q*vertex.p.E/vertex.p.P*tcos)
	s5lab = gtpr*sig*vertex.q*vertex.p.P*factor/pi/1.d+06

* 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) ) +
     >		( (ppicmx+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) ) +
     >		( (ppicmy+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) +
     > 		((ppicmz+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)) +
     >		( (ppicmx+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)) +
     >		( (ppicmy+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))
     >		+((ppicmz+gstar*bstarz*vertex.p.E)/vertex.p.P-gstar*bstarz*
     >		vertex.p.P/vertex.p.E)*dp_dcos

	dpxnewdphi = dpxdphi*new_x_x+dpydphi*new_x_y+dpzdphi*new_x_z
	dpynewdphi = dpxdphi*new_y_x+dpydphi*new_y_y+dpzdphi*new_y_z
	dpznewdphi = dpxdphi*new_z_x+dpydphi*new_z_y+dpzdphi*new_z_z

	dphicmdphi = (dpynewdphi*ppicm_newx - ppicm_newy*dpxnewdphi)/
     >			(ppicm_newx**2+ppicm_newy**2)

	dpxnewdcos = dpxdcos*new_x_x+dpydcos*new_x_y+dpzdcos*new_x_z
	dpynewdcos = dpxdcos*new_y_x+dpydcos*new_y_y+dpzdcos*new_y_z
	dpznewdcos = dpxdcos*new_z_x+dpydcos*new_z_y+dpzdcos*new_z_z

	dphicmdcos = (dpynewdcos*ppicm_newx - ppicm_newy*dpxnewdcos)
     >			/(ppicm_newx**2+ppicm_newy**2)

	dcoscmdcos = dpznewdcos/ppicm-ppicm_newz*epicm*dEcmdcos/ppicm**1.5
	dcoscmdphi = dpznewdphi/ppicm-ppicm_newz*epicm*dEcmdphi/ppicm**1.5

	jacobian = abs(dt_dcos_lab*dphicmdphi-dt_dphi_lab*dphicmdcos)
	main.davejac = jacobian

	main.johnjac = 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

c	davesig = gtpr*sig*jacobian

	davesig = gtpr*sig
	sigma_eerho = davesig/1.d3	!ub/GeV-sr --> ub/MeV-sr
	peerho = sigma_eerho
	ntup.sigcm = sig	!sig_cm

202	format(/11X,f5.1/)
203	format(11X,f5.0)
204	format(6(/9X,7f8.3))

	return
	end