From 0053203c2b5bc9d4a1403ce80df6282bca2a8c33 Mon Sep 17 00:00:00 2001
From: Dave Gaskell <gaskelld@jlab.org>
Date: Fri, 8 Jul 2022 15:35:45 -0400
Subject: [PATCH] Fixes for rho decay
---
rho_decay.f | 53 +++++++++++++++++++++++++++++++++++++++++++--------
rho_physics.f | 46 ++++++++++++++++++++++++++------------------
2 files changed, 72 insertions(+), 27 deletions(-)
diff --git a/rho_decay.f b/rho_decay.f
index 64db156..7ecf7e7 100644
--- a/rho_decay.f
+++ b/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
diff --git a/rho_physics.f b/rho_physics.f
index db6ff93..fd87c63 100644
--- a/rho_physics.f
+++ b/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
--
GitLab