pygbeh.f

Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file pygbeh.f

 
 
C*********************************************************************
 
C...PYGBEH
C...Evaluates the Bethe-Heitler cross section for heavy flavour
C...production.
C...Adapted from SaSgam library, authors G.A. Schuler and T. Sjostrand.
 
      SUBROUTINE pygbeh(KF,X,Q2,P2,PM2,XPBH)
 
C...Double precision and integer declarations.
      IMPLICIT DOUBLE PRECISION(a-h, o-z)
      IMPLICIT INTEGER(i-n)
      INTEGER pyk,pychge,pycomp
 
C...Local data.
      DATA aem2pi/0.0011614d0/
 
C...Reset output.
      xpbh=0d0
      sigbh=0d0
 
C...Check kinematics limits.
      IF(x.GE.q2/(4d0*pm2+q2+p2)) RETURN
      w2=q2*(1d0-x)/x-p2
      beta2=1d0-4d0*pm2/w2
      IF(beta2.LT.1d-10) RETURN
      beta=sqrt(beta2)
      rmq=4d0*pm2/q2
 
C...Simple case: P2 = 0.
      IF(p2.LT.1d-4) THEN
        IF(beta.LT.0.99d0) THEN
          xbl=log((1d0+beta)/(1d0-beta))
        ELSE
          xbl=log((1d0+beta)**2*w2/(4d0*pm2))
        ENDIF
        sigbh=beta*(8d0*x*(1d0-x)-1d0-rmq*x*(1d0-x))+
     &  xbl*(x**2+(1d0-x)**2+rmq*x*(1d0-3d0*x)-0.5d0*rmq**2*x**2)
 
C...Complicated case: P2 > 0, based on approximation of
C...C.T. Hill and G.G. Ross, Nucl. Phys. B148 (1979) 373
      ELSE
        rpq=1d0-4d0*x**2*p2/q2
        IF(rpq.GT.1d-10) THEN
          rpbe=sqrt(rpq*beta2)
          IF(rpbe.LT.0.99d0) THEN
            xbl=log((1d0+rpbe)/(1d0-rpbe))
            xbi=2d0*rpbe/(1d0-rpbe**2)
          ELSE
            rpbesn=4d0*pm2/w2+(4d0*x**2*p2/q2)*beta2
            xbl=log((1d0+rpbe)**2/rpbesn)
            xbi=2d0*rpbe/rpbesn
          ENDIF
          sigbh=beta*(6d0*x*(1d0-x)-1d0)+
     &    xbl*(x**2+(1d0-x)**2+rmq*x*(1d0-3d0*x)-0.5d0*rmq**2*x**2)+
     &    xbi*(2d0*x/q2)*(pm2*x*(2d0-rmq)-p2*x)
        ENDIF
      ENDIF
 
C...Multiply by charge-squared etc. to get parton distribution.
      chsq=1d0/9d0
      IF(iabs(kf).EQ.2.OR.iabs(kf).EQ.4) chsq=4d0/9d0
      xpbh=3d0*chsq*aem2pi*x*sigbh
 
      RETURN
      END