Fixes for rho decay

rho_decay.f

@@ -36,6 +36,10 @@ C Local declarations.
 	real*8 norm,dist
 	real*8 grnd
 
+	real*8 mmu,mt2,ml2,ml,ksi2
+	real*8 pol,rtmp,rtmp2,angle,costh
+	real*8 delta,fracl
+
 	logical success
 
 	parameter (ctaurho=1.31467e-13)
@@ -45,9 +49,29 @@ C Calulate R = sigmaL/sigmaT
 C Warning!  Note that if you change the parameterization of R in
 C rho-physics, do it here as well for self consistency!
 C This parameterization is taken from HERMES data...
+C DJG	R_rho = 0.33*(orig%Q2/mrho2)**(0.61)
+C DJG Update to use HepGen parameterization
+	mmu=105.66
+	mt2 = 0.62*mrho2
+	ml2 = 1.5*mt2
+	ml = sqrt(ml2)
+	ksi2 = 1.0
+
+
 
-	R_rho = 0.33*(orig%Q2/mrho2)**(0.61)
+	R_rho=(1.+orig%Q2/mt2)**2 * ksi2 *
+     >    (pi/2.*ml2/orig%Q2 - ml2*ml/(sqrt(orig%Q2)*(orig%Q2+ml2))
+     >    -ml2/orig%Q2*atan(ml/sqrt(orig%Q2)))**2
 
+	delta = 2.0*mmu**2*(1.0-epsilon)/orig%Q2
+	fracl = (epsilon+delta)*R_rho/(1.0+(epsilon+delta)*R_rho)
+
+	rtmp=grnd()
+	if(rtmp.lt.fracl) then
+	   pol=0.0 ! long.
+	else
+	   pol=1.0 ! transverse
+	endif
 
 	ph = orig%p%P
 	beta = ph/sqrt(ph**2+Mh2)
@@ -56,13 +80,26 @@ C This parameterization is taken from HERMES data...
 
 C Generate center/mass decay angles and momenta.
 	rph = grnd()*2.*pi
- 100	rth1 = grnd()*2.-1.
-	rth = acos(rth1)
 
+	rtmp2 = grnd()
+	if(pol.eq.0.0) then
+	   if(rtmp2.lt.0.5) then
+	      costh=-1.0*abs(2*rtmp2-1.)**(1./3.)
+	   elseif(rtmp2.gt.0.5) then
+	      costh=+1.0*abs(2*rtmp2-1.)**(1./3.)
+	   else
+	      costh=0.0
+	   endif
+	elseif(pol.eq.1.0) then
+	   angle=acos(abs(2*rtmp2-1.))
+	   if(rtmp2.lt.0.5) then
+	      costh=-2.0*cos((pi+angle)/3.)
+	   else
+	      costh=+2.0*cos((pi+angle)/3.)
+	   endif
+	endif
+	rth=acos(costh)
 
-	norm = (1.0+2.0*epsilon*R_rho)*grnd()
-	dist = sin(rth)**2+2.0*epsilon*R_rho*cos(rth)**2
-	if(dist.lt.norm) goto 100
 c       write(6,*) 'now decaying the rho',ctaurho,mh2
 	ntup%rhotheta=rth
 c	er = 384.65
@@ -117,7 +154,7 @@ C DJG If pion's heading for spectrometer, try to calculate xptar and yptar
 	   else
 	      th_inp = spec%p%theta+th_inp
 	   endif
-	elseif(hadron_arm.eq.2.or.hadron_arm.eq.4) then
+	elseif(hadron_arm.eq.2.or.hadron_arm.eq.4.or.hadron_arm.eq.5) then
 	   if(pyf.gt.0.0) then
 	      th_inp = th_inp - spec%p%theta
 	   else

rho_physics.f

@@ -238,7 +238,8 @@ c	  mtar_offshell = sqrt(efer**2-pfer**2)
 ! 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	sig0 = 41.263/(vertex%nu/1000.0)**0.4765   ! microbarns
+! DJG - replace w/HepGen (from Vardan)
 
 ! 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.
@@ -246,36 +247,43 @@ c	  mtar_offshell = sqrt(efer**2-pfer**2)
 ! DJG:  should change it in rho_decay also if you want to be self
 ! DJG:  consistent.
 
-	R = 0.33*(vertex%Q2/mrho2)**(0.61)
+c DJG	R = 0.33*(vertex%Q2/mrho2)**(0.61)
 ! PYB added this
-	if(R.lt.0.) write(222,'(''r'',3e12.3)')
-     >    r,vertex%Q2,mrho2
-	if(R.lt.0.) R=0.
+c DJG	if(R.lt.0.) write(222,'(''r'',3e12.3)')
+c DJG     >    r,vertex%Q2,mrho2
+c DJG	if(R.lt.0.) R=0.
 
 ! 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)
+C DJG	sigt = sig0*(1.0+epsi*R)*(mrho2/(vertex%Q2+mrho2))**(2.575)
 
 ! PYB added this
-	if(sigt.lt.0.) write(222,'(''sigt'',6e12.3)')
-     >    sigt, sig0, epsi,r,mrho2,vertex%q2
-	if(sigt<0.) sigt = 0.
-
+c	if(sigt.lt.0.) write(222,'(''sigt'',6e12.3)')
+c     >    sigt, sig0, epsi,r,mrho2,vertex%q2
+c	if(sigt<0.) sigt = 0.
 
 ! 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*log10(cdeltatau)
+! DJG: March 2022: HepGen uses a fixed value 
+c DJG	cdeltatau = hbarc/(sqrt(vertex%nu**2+vertex%Q2+mrho2)-vertex%nu) !in fm!
+c DJG	if(cdeltatau.lt.2.0) then
+c DJG	   brho = 4.4679 + 8.6106*log10(cdeltatau)
 ! PYB added this
-	   if(brho.lt.1.0) brho = 1.0
-	else
-	   brho = 7.0
-	endif
+c DJG	   if(brho.lt.1.0) brho = 1.0
+c DJG	else
+c DJG	   brho = 7.0
+c DJG	endif
 
-	sig219 = sigt*brho*exp(-brho*tprime)/2.0/pi !ub/GeV**2/rad
 
+c DJG HepGen model (from Vardan)
+	sig0 = 1.0D-3*27.4*(6.0/(vertex%Q2/1.0E6))**1.96   ! microbarns
+	sigt = sig0*2.33
+	if(sigt.lt.0.) write(222,'(''sigt'',6e12.3)')
+     >    sigt, sig0, epsi,r,mrho2,vertex%q2
+	if(sigt<0.) sigt = 0.
+	brho=5.0
+
+	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