pygbeh.F

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

      SUBROUTINE pygbeh(KF,X,Q2,P2,PM2,XPBH)
C...Purpose: to evaluate the Bethe-Heitler cross section for
C...heavy flavour production.
      DATA aem2pi/0.0011614/
 
C...Reset output.
      xpbh=0.
      sigbh=0.
 
C...Check kinematics limits.
      IF(x.GE.q2/(4.*pm2+q2+p2)) RETURN
      w2=q2*(1.-x)/x-p2
      beta2=1.-4.*pm2/w2
      IF(beta2.LT.1e-10) RETURN
      rmq=4.*pm2/q2
 
C...Simple case: P2 = 0.
      IF(p2.LT.1e-4) THEN
        beta=sqrt(beta2)
        IF(beta.LT.0.99) THEN
          xbl=log((1.+beta)/(1.-beta))
        ELSE
          xbl=log((1.+beta)**2*w2/(4.*pm2))
        ENDIF
        sigbh=beta*(8.*x*(1.-x)-1.-rmq*x*(1.-x))+
     &  xbl*(x**2+(1.-x)**2+rmq*x*(1.-3.*x)-0.5*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=1.-4.*x**2*p2/q2
        IF(rpq.GT.1e-10) THEN
          rpbe=sqrt(rpq*beta2)
          IF(rpbe.LT.0.99) THEN
            xbl=log((1.+rpbe)/(1.-rpbe))
            xbi=2.*rpbe/(1.-rpbe**2)
          ELSE
            rpbesn=4.*pm2/w2+(4.*x**2*p2/q2)*beta2
            xbl=log((1.+rpbe)**2/rpbesn)
            xbi=2.*rpbe/rpbesn
          ENDIF
          sigbh=beta*(6.*x*(1.-x)-1.)+
     &    xbl*(x**2+(1.-x)**2+rmq*x*(1.-3.*x)-0.5*rmq**2*x**2)+
     &    xbi*(2.*x/q2)*(pm2*x*(2.-rmq)-p2*x)
        ENDIF
      ENDIF
 
C...Multiply by charge-squared etc. to get parton distribution.
      chsq=1./9.
      IF(iabs(kf).EQ.2.OR.iabs(kf).EQ.4) chsq=4./9.
      xpbh=3.*chsq*aem2pi*x*sigbh
 
      RETURN
      END