pythiaeRHIC/pythia/pygbeh.f

0001
0002
0003 C*********************************************************************
0004
0005 C...PYGBEH
0006 C...Evaluates the Bethe-Heitler cross section for heavy flavour
0007 C...production.
0008 C...Adapted from SaSgam library, authors G.A. Schuler and T. Sjostrand.
0009
0010       SUBROUTINE PYGBEH(KF,X,Q2,P2,PM2,XPBH)
0011
0012 C...Double precision and integer declarations.
0013       IMPLICIT DOUBLE PRECISION(A-H, O-Z)
0014       IMPLICIT INTEGER(I-N)
0015       INTEGER PYK,PYCHGE,PYCOMP
0016
0017 C...Local data.
0018       DATA AEM2PI/0.0011614D0/
0019
0020 C...Reset output.
0021       XPBH=0D0
0022       SIGBH=0D0
0023
0024 C...Check kinematics limits.
0025       IF(X.GE.Q2/(4D0*PM2+Q2+P2)) RETURN
0026       W2=Q2*(1D0-X)/X-P2
0027       BETA2=1D0-4D0*PM2/W2
0028       IF(BETA2.LT.1D-10) RETURN
0029       BETA=SQRT(BETA2)
0030       RMQ=4D0*PM2/Q2
0031
0032 C...Simple case: P2 = 0.
0033       IF(P2.LT.1D-4) THEN
0034         IF(BETA.LT.0.99D0) THEN
0035           XBL=LOG((1D0+BETA)/(1D0-BETA))
0036         ELSE
0037           XBL=LOG((1D0+BETA)**2*W2/(4D0*PM2))
0038         ENDIF
0039         SIGBH=BETA*(8D0*X*(1D0-X)-1D0-RMQ*X*(1D0-X))+
0040      &  XBL*(X**2+(1D0-X)**2+RMQ*X*(1D0-3D0*X)-0.5D0*RMQ**2*X**2)
0041
0042 C...Complicated case: P2 > 0, based on approximation of
0043 C...C.T. Hill and G.G. Ross, Nucl. Phys. B148 (1979) 373
0044       ELSE
0045         RPQ=1D0-4D0*X**2*P2/Q2
0046         IF(RPQ.GT.1D-10) THEN
0047           RPBE=SQRT(RPQ*BETA2)
0048           IF(RPBE.LT.0.99D0) THEN
0049             XBL=LOG((1D0+RPBE)/(1D0-RPBE))
0050             XBI=2D0*RPBE/(1D0-RPBE**2)
0051           ELSE
0052             RPBESN=4D0*PM2/W2+(4D0*X**2*P2/Q2)*BETA2
0053             XBL=LOG((1D0+RPBE)**2/RPBESN)
0054             XBI=2D0*RPBE/RPBESN
0055           ENDIF
0056           SIGBH=BETA*(6D0*X*(1D0-X)-1D0)+
0057      &    XBL*(X**2+(1D0-X)**2+RMQ*X*(1D0-3D0*X)-0.5D0*RMQ**2*X**2)+
0058      &    XBI*(2D0*X/Q2)*(PM2*X*(2D0-RMQ)-P2*X)
0059         ENDIF
0060       ENDIF
0061
0062 C...Multiply by charge-squared etc. to get parton distribution.
0063       CHSQ=1D0/9D0
0064       IF(IABS(KF).EQ.2.OR.IABS(KF).EQ.4) CHSQ=4D0/9D0
0065       XPBH=3D0*CHSQ*AEM2PI*X*SIGBH
0066
0067       RETURN
0068       END