C*********************************************************************
 
C...PYSTBH (and auxiliaries)
C.. Evaluates the matrix elements for t + b + H production.
 
      SUBROUTINE PYSTBH(WTTBH)
 
C...DOUBLE PRECISION AND INTEGER DECLARATIONS
      IMPLICIT DOUBLE PRECISION(A-H, O-Z)
      IMPLICIT INTEGER(I-N)
      INTEGER PYK,PYCHGE,PYCOMP
 
C...COMMONBLOCKS
      COMMON/PYDAT1/MSTU(200),PARU(200),MSTJ(200),PARJ(200)
      COMMON/PYDAT2/KCHG(500,4),PMAS(500,4),PARF(2000),VCKM(4,4)
      COMMON/PYPARS/MSTP(200),PARP(200),MSTI(200),PARI(200)
      COMMON/PYINT1/MINT(400),VINT(400)
      COMMON/PYINT2/ISET(500),KFPR(500,2),COEF(500,20),ICOL(40,4,2)
      COMMON/PYINT3/XSFX(2,-40:40),ISIG(1000,3),SIGH(1000)
      COMMON/PYINT4/MWID(500),WIDS(500,5)
      COMMON/PYSUBS/MSEL,MSELPD,MSUB(500),KFIN(2,-40:40),CKIN(200)
      COMMON/PYMSSM/IMSS(0:99),RMSS(0:99)
      COMMON/PYSGCM/ISUB,ISUBSV,MMIN1,MMAX1,MMIN2,MMAX2,MMINA,MMAXA,
     &KFAC(2,-40:40),COMFAC,FACK,FACA,SH,TH,UH,SH2,TH2,UH2,SQM3,SQM4,
     &SHR,SQPTH,TAUP,BE34,CTH,X(2),SQMZ,SQMW,GMMZ,GMMW,
     &AEM,AS,XW,XW1,XWC,XWV,POLL,POLR,POLLL,POLRR
      COMMON/PYCTBH/ ALPHA,ALPHAS,SW2,MW2,TANB,VTB,V,A
      DOUBLE PRECISION MW2
      SAVE /PYDAT1/,/PYDAT2/,/PYPARS/,/PYINT1/,/PYINT2/,/PYINT3/,
     &/PYINT4/,/PYSUBS/,/PYMSSM/,/PYSGCM/,/PYCTBH/
 
C...LOCAL ARRAYS AND COMPLEX VARIABLES
      DIMENSION QQ(4,2),PP(4,3)
      SAVE
      DATA QQ/8*0D0/
 
      WTTBH=0D0
 
C...KINEMATIC PARAMETERS.
      SHPR=SQRT(VINT(26))*VINT(1)
      PH=SQRT(VINT(21))*VINT(1)
      SPH=PH**2
 
C...SET UP OUTGOING KINEMATICS: 1=T, 2=TBAR, 3=H.
      DO 100 I=1,2
        PT=SQRT(MAX(0D0,VINT(197+5*I)))
        PP(1,I)=PT*COS(VINT(198+5*I))
        PP(2,I)=PT*SIN(VINT(198+5*I))
  100 CONTINUE
      PP(1,3)=-PP(1,1)-PP(1,2)
      PP(2,3)=-PP(2,1)-PP(2,2)
      PMS1=VINT(201)**2+PP(1,1)**2+PP(2,1)**2
      PMS2=VINT(206)**2+PP(1,2)**2+PP(2,2)**2
      PMS3=SPH+PP(1,3)**2+PP(2,3)**2
      PMT3=SQRT(PMS3)
      PP(3,3)=PMT3*SINH(VINT(211))
      PP(4,3)=PMT3*COSH(VINT(211))
      PMS12=(SHPR-PP(4,3))**2-PP(3,3)**2
      PP(3,1)=(-PP(3,3)*(PMS12+PMS1-PMS2)+
     &VINT(213)*(SHPR-PP(4,3))*VINT(220))/(2D0*PMS12)
      PP(3,2)=-PP(3,1)-PP(3,3)
      PP(4,1)=SQRT(PMS1+PP(3,1)**2)
      PP(4,2)=SQRT(PMS2+PP(3,2)**2)
 
C...CM SYSTEM, INGOING QUARKS/GLUONS
      QQ(3,1) = SHPR/2.D0
      QQ(4,1) = QQ(3,1)
      QQ(3,2) = -QQ(3,1)
      QQ(4,2) = QQ(4,1)
 
C...PARAMETERS FOR AMPLITUDE METHOD
      ALPHA = AEM
      ALPHAS = AS
      SW2 = PARU(102)
      MW2 = PMAS(24,1)**2
      TANB = PARU(141)
      VTB = VCKM(3,3)
      RMB=PYMRUN(5,VINT(52))
 
      ISUB=MINT(1)
 
      IF (ISUB.EQ.401) THEN
        CALL PYTBHG(QQ(1,1),QQ(1,2),PP(1,1),PP(1,2),PP(1,3),
     &  VINT(201),VINT(206),RMB,VINT(43),WTTBH)
      ELSE IF (ISUB.EQ.402) THEN
        CALL PYTBHQ(QQ(1,1),QQ(1,2),PP(1,1),PP(1,2),PP(1,3),
     &  VINT(201),VINT(206),RMB,VINT(43),WTTBH)
      END IF
 
      RETURN
      END