PYTHIA/pythia/pygbeh.F

   1
   2 C*********************************************************************
   3
   4       SUBROUTINE PYGBEH(KF,X,Q2,P2,PM2,XPBH)
   5 C...Purpose: to evaluate the Bethe-Heitler cross section for
   6 C...heavy flavour production.
   7       DATA AEM2PI/0.0011614/
   8
   9 C...Reset output.
  10       XPBH=0.
  11       SIGBH=0.
  12
  13 C...Check kinematics limits.
  14       IF(X.GE.Q2/(4.*PM2+Q2+P2)) RETURN
  15       W2=Q2*(1.-X)/X-P2
  16       BETA2=1.-4.*PM2/W2
  17       IF(BETA2.LT.1E-10) RETURN
  18       RMQ=4.*PM2/Q2
  19
  20 C...Simple case: P2 = 0.
  21       IF(P2.LT.1E-4) THEN
  22         BETA=SQRT(BETA2)
  23         IF(BETA.LT.0.99) THEN
  24           XBL=LOG((1.+BETA)/(1.-BETA))
  25         ELSE
  26           XBL=LOG((1.+BETA)**2*W2/(4.*PM2))
  27         ENDIF
  28         SIGBH=BETA*(8.*X*(1.-X)-1.-RMQ*X*(1.-X))+
  29      &  XBL*(X**2+(1.-X)**2+RMQ*X*(1.-3.*X)-0.5*RMQ**2*X**2)
  30
  31 C...Complicated case: P2 > 0, based on approximation of
  32 C...C.T. Hill and G.G. Ross, Nucl. Phys. B148 (1979) 373
  33       ELSE
  34         RPQ=1.-4.*X**2*P2/Q2
  35         IF(RPQ.GT.1E-10) THEN
  36           RPBE=SQRT(RPQ*BETA2)
  37           IF(RPBE.LT.0.99) THEN
  38             XBL=LOG((1.+RPBE)/(1.-RPBE))
  39             XBI=2.*RPBE/(1.-RPBE**2)
  40           ELSE
  41             RPBESN=4.*PM2/W2+(4.*X**2*P2/Q2)*BETA2
  42             XBL=LOG((1.+RPBE)**2/RPBESN)
  43             XBI=2.*RPBE/RPBESN
  44           ENDIF
  45           SIGBH=BETA*(6.*X*(1.-X)-1.)+
  46      &    XBL*(X**2+(1.-X)**2+RMQ*X*(1.-3.*X)-0.5*RMQ**2*X**2)+
  47      &    XBI*(2.*X/Q2)*(PM2*X*(2.-RMQ)-P2*X)
  48         ENDIF
  49       ENDIF
  50
  51 C...Multiply by charge-squared etc. to get parton distribution.
  52       CHSQ=1./9.
  53       IF(IABS(KF).EQ.2.OR.IABS(KF).EQ.4) CHSQ=4./9.
  54       XPBH=3.*CHSQ*AEM2PI*X*SIGBH
  55
  56       RETURN
  57       END