      SUBROUTINE CHECKINPUT_ST
      IMPLICIT REAL*8 (A-H,M-Z)
      IMPLICIT INTEGER (I,J)
      COMMON/CONST6/S,ENERGY,ALPHAS,MST1,MG,MT,MS,MST2,SIN2T
      COMMON/CONST2/SCALE,SCAFAC,ICOLL,ISCAPT
      COMMON/CONST5/ILO,INLO,IONLYLO
      COMMON/CUT1/PTMIN,PTMAX
      COMMON/CUT2/YMIN,YMAX
      COMMON/FLAVOR/IFLAVOR,ITOTAL
      COMMON/CHARCONJ/ICHARCONJ

C***  CHECK MST1, MG, MT, MS AND MST2 ( SHOULD BE ABOVE 1 GEV )
C***  CHECK SIN2T 
      IF (MST1.LT.1D0) THEN
         MST1 = 1D0
         PRINT *,'MST1 CHANGED TO 1.'
      END IF
      IF (MG.LT.1D0) THEN
         MG = 1D0
         PRINT *,'MG CHANGED TO 1.'
      END IF
      IF (MT.LT.1D0) THEN
         MT = 1D0
         PRINT *,'MT CHANGED TO 1.'
      END IF
      IF (MS.LT.1D0) THEN
         MS = 1D0
         PRINT *,'MS CHANGED TO 1.'
      END IF
      IF (MST2.LT.1D0) THEN
         MST2 = 1D0
         PRINT *,'MST2 CHANGED TO 1.'
      END IF
      IF (ABS(SIN2T).GT.1D0) THEN
         SIN2T = 0D0
         PRINT *,'SIN2T CHANGED TO 0.'
      END IF

C***  CHECK ICOLL
      IF ((ICOLL.LT.0).OR.(ICOLL.GT.1)) THEN
         PRINT *,'COLLIDER TYPE SPECIFIED INCORRECTLY; ICOLL = ',ICOLL
         PRINT *,'PROGRAM ABORTED'
         call HARD_STOP
      END IF

C***  CHECK ENERGY ( SHOULD BE ABOVE PRODUCTION THRESHOLD )
      IF (ENERGY.LT.2*MST1) THEN
         PRINT *,'ENERGY TOO SMALL TO PRODUCE STOP-ANTISTOP'
         PRINT *,'PROGRAM ABORTED'
         call HARD_STOP
      END IF

C***  CHECK IFLAVOR
      IF ((IFLAVOR.LT.0).OR.(IFLAVOR.GT.5)) THEN
         PRINT *,'IFLAVOR IN THE WRONG RANGE; IFLAVOR = ',IFLAVOR
         PRINT *,'IFLAVOR CHANGED TO 0 (ALL)'
         IFLAVOR = 0
      END IF

C***  CHECK ITOTAL
      IF ((ITOTAL.LT.0).OR.(ITOTAL.GT.1)) THEN
         PRINT *,'ITOTAL IN THE WRONG RANGE; ITOTAL = ',ITOTAL
         PRINT *,'ITOTAL CHANGED TO 1'
         ITOTAL = 1
      END IF

C***  CHECK ISCAPT AND SCAFAC
      IF ((ISCAPT.LT.0).OR.(ISCAPT.GT.1)) THEN
         PRINT *,'ISCAPT IN THE WRONG RANGE; ISCAPT = ',ISCAPT
         PRINT *,'ISCAPT CHANGED TO 0'
         ISCAPT = 0
      END IF
      IF (SCAFAC.LT.1D-2) THEN
         PRINT *, 'SCAFAC TOO SMALL'
         PRINT *, 'SCAFAC CHANGED TO 1.'
         SCAFAC = 1.D0
      END IF
      IF ((ITOTAL.EQ.1).AND.(ISCAPT.EQ.1)) THEN
         PRINT *,'WARNING: TRANSVERSE MASS IS NOT USED AS A SCALE'
      END IF

C***  CHECK PTMIN AND PTMAX
      IF (PTMIN.LT.0D0) THEN
         PTMIN = 0D0
         PRINT *,'PTMIN NEGATIVE'
         PRINT *,'PTMIN CHANGED TO 0.'
      END IF
      IF (PTMAX.LT.0D0) THEN
         PTMAX = 0D0
         PRINT *,'PTMAX NEGATIVE'
         PRINT *,'PTMAX CHANGED TO 0.'
      END IF
      IF (PTMAX.LT.PTMIN) THEN
         PTMAX = PTMIN
         PRINT *,'PTMAX SMALLER THAN PTMIN'
         PRINT *,'PTMAX CHANGED TO PTMIN'
      END IF
      IF ((PTMAX.EQ.PTMIN).AND.(PTMAX.EQ.0D0)) THEN
         PTMAX = 0.1D0
         PTMIN = 0.1D0
         PRINT *,'PTMAX = PTMIN CHANGED TO 0.1'
      END IF

C***  CHECK YMIN AND YMAX
      IF (YMIN.LT.0D0) THEN
         YMIN = 0D0
         PRINT *,'YMIN NEGATIVE'
         PRINT *,'YMIN CHANGED TO 0.'
      END IF
      IF (YMAX.LT.0D0) THEN
         YMAX = 0D0
         PRINT *,'YMAX NEGATIVE'
         PRINT *,'YMAX CHANGED TO 0.'
      END IF
      IF (YMAX.LT.YMIN) THEN
         YMAX = YMIN
         PRINT *,'YMAX SMALLER THAN YMIN'
         PRINT *,'YMAX CHANGED TO YMIN'
      END IF

C***  CHECK ICHACONJ
      IF ((ICHARCONJ.LT.-1).OR.(ICHARCONJ.GT.1)) THEN
         PRINT *,'ICHARCONJ IN THE WRONG RANGE; ICHARCONJ = ',ICHARCONJ
         PRINT *,'ICHARCONJ CHANGED TO 0'
         ICHARCONJ = 0
      END IF

C***  CHECK IONLYLO
      IF ((IONLYLO.LT.0).OR.(IONLYLO.GT.1)) THEN
         PRINT *,'IONLYLO IN THE WRONG RANGE'
         PRINT *,'IONLYLO CHANGED TO 0'
         IONLYLO = 0
      END IF

C***  CHECK ILO, INLO
      IF (ILO.LT.10) THEN
         ILO = 10
         PRINT *,'ILO CHANGED TO 10'
      END IF
      IF (INLO.LT.10) THEN
         INLO = 10
         PRINT *,'INLO CHANGED TO 10'
      END IF

      RETURN
      END


C ======================================================================

      SUBROUTINE INTEGST(RESLO,ERRLO,RESNLO,ERRNLO)
      IMPLICIT REAL*8 (A-H,M-Z)
      IMPLICIT INTEGER (I,J)
      integer  ivegas(1:4)
      COMMON/CONST6/S,ENERGY,ALPHAS,MST1,MG,MT,MS,MST2,SIN2T
      COMMON/CONST5/ILO,INLO,IONLYLO

      EXTERNAL SIGLOST, SIGNLOST

      CALL CHECKINPUT_ST

      S = ENERGY**2

      CALL INILO(ivegas)
      idim = 3
      acc = 1.d-4
      call INTEG(SIGLOST,idim,ivegas,acc,reslo,errlo)

      IF (IONLYLO.EQ.0) THEN
         CALL ININLO(ivegas)
         idim = 5 
         call INTEG(SIGNLOST,idim,ivegas,acc,resnlo,errnlo)
      ELSE 
         RESNLO = 0D0
         ERRNLO = 0D0
      END IF

      RETURN
      END

C ======================================================================

C***  THE CROSS-SECTION IN LO

      REAL*8 FUNCTION SIGLOST(VAR)
      IMPLICIT REAL*8 (A-H,M-Z)
      IMPLICIT INTEGER (I,J)
      DIMENSION VAR(3)
      COMMON/CONST6/S,ENERGY,ALPHAS,MST1,MG,MT,MS,MST2,SIN2T
      COMMON/CONST2/SCALE,SCAFAC,ICOLL,ISCAPT
      COMMON/PTYST/PT,PREPT,MTRANS,YVAL,PREY
      COMMON/FLAVOR/IFLAVOR,ITOTAL
      COMMON/CHARCONJ/ICHARCONJ
      CALL DVEGLO_ST(VAR)
      IF (ITOTAL.EQ.0) THEN
         CALL DPTYST
         CALL DEFSCAST
         CALL DDILOST
         SIGLOST1 = ESTGGLO()
         SIGLOST2 = ESTQBLO()

C***     SUM OVER POSITIVE AND NEGATIVE RAPIDITY
         IF ((ICOLL.EQ.0).AND.(ICHARCONJ.NE.0)) THEN
            YVAL = -YVAL
            CALL DDILOST
            SIGLOST1 = 2 * SIGLOST1
            SIGLOST2 = SIGLOST2 + ESTQBLO()
         ELSE 
            SIGLOST1 = 2 * SIGLOST1
            SIGLOST2 = 2 * SIGLOST2
         END IF
      ELSE
         CALL DEFSCAST
         CALL DTOLOST
         SIGLOST1 = FSTGGLO()
         SIGLOST2 = FSTQBLO()
      END IF
      SIGLOST =  SUMLOST(SIGLOST1,SIGLOST2,IFLAVOR)
      RETURN
      END

C ======================================================================

C***  THE CROSS-SECTION IN NLO

      REAL*8 FUNCTION SIGNLOST(VAR)
      IMPLICIT REAL*8 (A-H,M-Z)
      IMPLICIT INTEGER (I,J)
      DIMENSION VAR(4)
      COMMON/CONST6/S,ENERGY,ALPHAS,MST1,MG,MT,MS,MST2,SIN2T
      COMMON/CONST2/SCALE,SCAFAC,ICOLL,ISCAPT
      COMMON/PTYST/PT,PREPT,MTRANS,YVAL,PREY
      COMMON/FLAVOR/IFLAVOR,ITOTAL
      COMMON/CHARCONJ/ICHARCONJ

      CALL DVEGNLO_ST(VAR)
      CALL DMINLST
      IF (ITOTAL.EQ.0) THEN
         CALL DPTYST
         CALL DEFSCAST
         CALL DDINLST
         IF ((IFLAVOR.EQ.0).OR.(IFLAVOR.EQ.1).OR.(IFLAVOR.EQ.4))      
     +        SIGNLOST1 = ESTGGNLO()
         IF ((IFLAVOR.EQ.0).OR.(IFLAVOR.EQ.2).OR.(IFLAVOR.EQ.4)) 
     +        SIGNLOST2 = ESTQBNLO()
         IF ((IFLAVOR.EQ.0).OR.(IFLAVOR.EQ.3).OR.(IFLAVOR.EQ.5))
     +        SIGNLOST3 = ESTGBNLO()

C***     SUM OVER POSITIVE AND NEGATIVE RAPIDITY
         IF ((ICOLL.EQ.0).AND.(ICHARCONJ.NE.0)) THEN
            SIGNLOST1 = 2 * SIGNLOST1
            YVAL = -YVAL
            CALL DDINLST
            IF ((IFLAVOR.EQ.0).OR.(IFLAVOR.EQ.2).OR.(IFLAVOR.EQ.4)) 
     +           SIGNLOST2 = SIGNLOST2 + ESTQBNLO()
            IF ((IFLAVOR.EQ.0).OR.(IFLAVOR.EQ.3).OR.(IFLAVOR.EQ.5))
     +           SIGNLOST3 = SIGNLOST3 + ESTGBNLO()
         ELSE 
            SIGNLOST1 = 2 * SIGNLOST1
            SIGNLOST2 = 2 * SIGNLOST2
            SIGNLOST3 = 2 * SIGNLOST3
         END IF
      ELSE
         CALL DEFSCAST
         CALL DTONLST
         IF ((IFLAVOR.EQ.0).OR.(IFLAVOR.EQ.1).OR.(IFLAVOR.EQ.4)) 
     +        SIGNLOST1 = FSTGGNLO()
         IF ((IFLAVOR.EQ.0).OR.(IFLAVOR.EQ.2).OR.(IFLAVOR.EQ.4)) 
     +        SIGNLOST2 = FSTQBNLO()
         IF ((IFLAVOR.EQ.0).OR.(IFLAVOR.EQ.3).OR.(IFLAVOR.EQ.5)) 
     +        SIGNLOST3 = FSTGBNLO()
      END IF
      SIGNLOST =  SUMNLOST(SIGNLOST1,SIGNLOST2,SIGNLOST3,IFLAVOR)
      RETURN
      END

C ======================================================================

      SUBROUTINE DVEGLO_ST(VAR)
C***  FILLS UP THE COMMON BLOCK VARVEGAS WITH 3 VARIABLES
      IMPLICIT REAL*8 (A-H,M-Z)
      DIMENSION VAR(3)
      COMMON/VARVEGAS/Y1,Y2,Y3,Y4
      Y1 = VAR(1)
      Y2 = VAR(2)
      Y3 = VAR(3)
      RETURN
      END

C ======================================================================

      SUBROUTINE DVEGNLO_ST(VAR)
C***  FILLS UP THE COMMON BLOCK VARVEGAS WITH 5 VARIABLES
      IMPLICIT REAL*8 (A-H,M-Z)
      DIMENSION VAR(4)
      COMMON/VARVEGAS/Y1,Y2,Y3,Y4
      Y1 = VAR(1)
      Y2 = VAR(2)
      Y3 = VAR(3)
      Y4 = VAR(4)
      RETURN
      END

C ======================================================================

      SUBROUTINE SWAP_ST(X,Y)
C***  INTERCHANGES X AND Y
      IMPLICIT REAL*8 (A-H,M-Z)
      Z = X
      X = Y
      Y = Z
      RETURN
      END

C ======================================================================

      SUBROUTINE AVER_ST(X,Y)
C***  AVERAGES X AND Y
      IMPLICIT REAL*8 (A-H,M-Z)
      X = ( X + Y )/2D0
      Y = X
      RETURN
      END

C ======================================================================

      SUBROUTINE DEFSCAST()
C***  DEFINES THE RENORMALIZATION AND FACTORIZATION SCALE
C***  FILLS UP PART OF THE COMMON BLOCKS CONST6 AND CONST2
      IMPLICIT REAL*8 (A-H,M-Z)
      COMMON/CONST6/S,ENERGY,ALPHAS,MST1,MG,MT,MS,MST2,SIN2T
      COMMON/CONST2/SCALE,SCAFAC,ICOLL,ISCAPT
      COMMON/PTYST/PT,PREPT,MTRANS,YVAL,PREY
      COMMON/FLAVOR/IFLAVOR,ITOTAL

      IF ((ITOTAL.EQ.0).AND.(ISCAPT.EQ.1)) THEN
         SCALE = MTRANS * SCAFAC
      ELSE 
         SCALE = MST1 * SCAFAC
      END IF

C***  THE SCALE SHOULD BE BELOW 1 TEV, BECAUSE SOME PARTON DENSITIES
C***  ARE RESTRICTED TO THIS REGION

      IF (SCALE.GE.1000D0) SCALE = 1000D0

      IF (SCALE.LE.5D0) SCALE = 5D0

      ALPHAS = ALPS(SCALE)

      RETURN
      END

C ======================================================================

      SUBROUTINE DPTYST()
C***  DEFINES PT AND Y (INCLUDING CUTS) 
C***  FILLS UP THE COMMON BLOCK PTYST
      IMPLICIT REAL*8 (A-H,M-Z)
      COMMON/CONST6/S,ENERGY,ALPHAS,MST1,MG,MT,MS,MST2,SIN2T
      COMMON/CUT1/PTMINC,PTMAXC
      COMMON/CUT2/YMINC,YMAXC
      COMMON/VARVEGAS/Y1,Y2,Y3,Y4
      COMMON/PTYST/PT,PREPT,MTRANS,YVAL,PREY
      PT2MAX = S/4.D0 -MST1**2
      PTMAX = SQRT( DMIN1(PT2MAX,PTMAXC**2) )
      PTMIN = PTMINC
      IF (PTMAX.GT.PTMIN) THEN
         PT = (PTMAX-PTMIN) * Y1 +PTMIN
         PREPT = 2.D0 * PT * (PTMAX -PTMIN)
      ELSE IF (PTMAX.EQ.PTMIN) THEN
         PT = PTMIN
         PREPT = 2.D0 * PT
      ELSE IF (PTMAX.LT.PTMIN) THEN 
         PT = 0D0
         PREPT = 0D0
      END IF
      MTRANS = SQRT(PT**2 +MST1**2)
      YMAXK = +DACOSH(SQRT(S)/2.D0/MTRANS)
      YMINK = -DACOSH(SQRT(S)/2.D0/MTRANS)
      YMAX = DMIN1( YMAXK, YMAXC)
      YMIN = DMAX1( YMINK, YMINC, 0D0)

      IF (YMAX.GT.YMIN) THEN
         YVAL = (YMAX -YMIN) * Y2 + YMIN
         PREY = YMAX - YMIN
      ELSE IF (YMAX.EQ.YMIN) THEN         
         YVAL = YMIN
         PREY = 1.D0
      ELSE IF (YMAX.LT.YMIN) THEN
         YVAL = 0D0
         PREY = 0D0
      END IF
      RETURN
      END

C ======================================================================

      SUBROUTINE DDILOST()
C***  DEFINES THE KINEMATICS FOR DIFFERENTIAL CROSS-SECTIONS IN LO
C***  FILLS UP THE COMMON BLOCK DILOST AND PDLOST
      IMPLICIT REAL*8 (A-H,M-Z)
      COMMON/VARVEGAS/Y1,Y2,Y3,Y4
      COMMON/CONST6/S,ENERGY,ALPHAS,MST1,MG,MT,MS,MST2,SIN2T
      COMMON/PTYST/PT,PREPT,MTRANS,YVAL,PREY
      COMMON/DILOST/SHATB,T1HATB,U1HATB,PREFB,X1,X2B
      COMMON/PDLOST/PDGGB,PDQBB,PDBQB
      T1 = -SQRT(S)*MTRANS*DEXP(+YVAL)
      U1 = -SQRT(S)*MTRANS*DEXP(-YVAL)
      X1MIN = -T1/(S+U1)
      X1 = X1MIN ** Y3
      X2B= -X1*U1/(X1*S +T1)
      SHATB  = X1 *X2B*S
      T1HATB = X2B*T1
      U1HATB = X1 *U1
      PREFB = S * PREPT * PREY *LOG(X1MIN)*(-X1*X2B)/(-X1*U1)
      CALL PDFLOST(X1,X2B)
      RETURN
      END

C ======================================================================

      SUBROUTINE DTOLOST()
C***  DEFINES THE KINEMATICS FOR TOTAL CROSS-SECTIONS IN LO
C***  FILLS UP THE COMMON BLOCK TOLOST AND PDLOST
      IMPLICIT REAL*8 (A-H,M-Z)
      COMMON/VARVEGAS/Y1,Y2,Y3,Y4
      COMMON/CONST6/S,ENERGY,ALPHAS,MST1,MG,MT,MS,MST2,SIN2T
      COMMON/TOLOST/SHAT,T1HAT,PREF,X1,X2
      COMMON/PDLOST/PDGGB,PDQBB,PDBQB
      X1 = (4*MST1**2/S) ** Y1
      X2 = (4*MST1**2/S/X1) ** Y2
      SHAT = X1*X2*S
      BETA = SQRT(1D0 -4*MST1**2/SHAT)
      T1HAT = -SHAT/2D0 +(Y3 -0.5D0)*SHAT*BETA
      PREF = SHAT*BETA *LOG(4*MST1**2/S)*LOG(4*MST1**2/S/X1)
      CALL PDFLOST(X1,X2)
      RETURN
      END

C ======================================================================

      REAL*8 FUNCTION ESTGGLO()
C***  CALCULATES THE DIFFERENTIAL CROSS-SECTION IN LO FOR G G
      IMPLICIT REAL*8 (A-H,M-Z)
      COMMON/PDLOST/PDGGB,PDQBB,PDBQB
      COMMON/CONST6/S,ENERGY,ALPHAS,MST1,MG,MT,MS,MST2,SIN2T
      COMMON/DILOST/SHATB,T1HATB,U1HATB,PREFB,X1,X2B
      ESTGGLO = 
     +     PDGGB *PREFB * DSTGGB(ALPHAS,SHATB,T1HATB,MST1,MG) 
      RETURN
      END

C ======================================================================

      REAL*8 FUNCTION ESTQBLO()
C***  CALCULATES THE DIFFERENTIAL CROSS-SECTION IN LO FOR Q QB
      IMPLICIT REAL*8 (A-H,M-Z)
      INTEGER ICHARCONJ
      COMMON/PDLOST/PDGGB,PDQBB,PDBQB
      COMMON/CONST6/S,ENERGY,ALPHAS,MST1,MG,MT,MS,MST2,SIN2T
      COMMON/DILOST/SHATB,T1HATB,U1HATB,PREFB,X1,X2B
      COMMON/CHARCONJ/ICHARCONJ
      
C***  DISTINGUISH BETWEEN STOPS AND ANTISTOPS IN THE FINAL STATE
C***  ICHARCONJ =  1  <-- STOPS
C***  ICHARCONJ = -1  <-- ANTISTOPS
C***  ICHARCONJ =  0  <-- AVERAGE OF STOPS AND ANTISTOPS

      IF (ICHARCONJ.EQ.1) THEN
         CALL SWAP_ST(PDQBB, PDBQB)
      ELSE IF (ICHARCONJ.EQ.0) THEN
         CALL AVER_ST(PDQBB, PDBQB)
      END IF

      ESTQBLO = 
     +     +  PDQBB * PREFB*DSTQBB(ALPHAS,SHATB,T1HATB,MST1,MG)
     +     +  PDBQB * PREFB*DSTQBB(ALPHAS,SHATB,U1HATB,MST1,MG) 

      RETURN
      END

C ======================================================================

      REAL*8 FUNCTION FSTGGLO()
C***  CALCULATES THE TOTAL CROSS-SECTION IN LO FOR G G
      IMPLICIT REAL*8 (A-H,M-Z)
      COMMON/PDLOST/PDGGB,PDQBB,PDBQB
      COMMON/CONST6/S,ENERGY,ALPHAS,MST1,MG,MT,MS,MST2,SIN2T
      COMMON/TOLOST/SHAT,T1HAT,PREF,X1,X2
      FSTGGLO = 
     +     PDGGB * PREF*DSTGGB(ALPHAS,SHAT,T1HAT,MST1,MG) 
      RETURN
      END

C ====================================================================== 
      REAL*8 FUNCTION FSTQBLO()
C***  CALCULATES THE TOTAL CROSS-SECTION IN LO FOR Q QB
      IMPLICIT REAL*8 (A-H,M-Z)
      COMMON/PDLOST/PDGGB,PDQBB,PDBQB
      COMMON/CONST6/S,ENERGY,ALPHAS,MST1,MG,MT,MS,MST2,SIN2T
      COMMON/TOLOST/SHAT,T1HAT,PREF,X1,X2
      FSTQBLO = 
     +     + (PDQBB +PDBQB ) * PREF*DSTQBB(ALPHAS,SHAT,T1HAT,MST1,MG) 
      RETURN
      END

C ======================================================================
 
      REAL*8 FUNCTION SUMLOST(SIGLOST1,SIGLOST2,IFLAVOR)
C***  SUMS UP DIFFERENT FLAVOR STATES FOR CROSS-SECTION IN LO
      IMPLICIT REAL*8 (A-H,M-Z)
      IMPLICIT INTEGER (I)
      IF ((IFLAVOR.EQ.0).OR.(IFLAVOR.EQ.4))
     +                  SUMLOST = SIGLOST1 + SIGLOST2
      IF (IFLAVOR.EQ.1) SUMLOST = SIGLOST1
      IF (IFLAVOR.EQ.2) SUMLOST = SIGLOST2
      IF ((IFLAVOR.EQ.3).OR.(IFLAVOR.EQ.5)) SUMLOST = 0D0
      RETURN
      END

C ======================================================================

      SUBROUTINE DMINLST()
C***  DEFINES SOME MISCELLANEOUS VARIABLES FOR CROSS-SECTION IN NLO
C***  FILLS UP THE COMMON BLOCK MINLST
      IMPLICIT REAL*8 (A-H,M-Z)
      COMMON/CONST6/S,ENERGY,ALPHAS,MST1,MG,MT,MS,MST2,SIN2T
      COMMON/CONST2/SCALE,SCAFAC,ICOLL,ISCAPT
      COMMON/MINLST/DEL,SCA
      DEL = 1.D-4 * MST1**2 
      SCA = SCALE**2/MST1**2
      RETURN
      END

C ======================================================================

      SUBROUTINE DDINLST()
C***  DEFINES THE KINEMATICS FOR DIFFERENTIAL CROSS-SECTIONS IN NLO
C***  FILLS UP THE COMMON BLOCK DINLST1-2 AND PDLOST, PDNLST1
      IMPLICIT REAL*8 (A-H,M-Z)
      COMMON/VARVEGAS/Y1,Y2,Y3,Y4
      COMMON/CONST6/S,ENERGY,ALPHAS,MST1,MG,MT,MS,MST2,SIN2T
      COMMON/PTYST/PT,PREPT,MTRANS,YVAL,PREY
      COMMON/MINLST/DEL,SCA
      COMMON/DINLST1/SHATB,T1HATB,U1HATB,PREFB,X1,X2B
      COMMON/DINLST2/SHAT,T1HAT,U1HAT,S4,PREF,X2,S4MAX
      COMMON/PDLOST/PDGGB,PDQBB,PDBQB
      COMMON/PDNLST1/PDGG,PDQB,PDBQ,PDQG, PDGQ,PDBG, PDGB

      T1 = -SQRT(S)*MTRANS*DEXP(+YVAL)
      U1 = -SQRT(S)*MTRANS*DEXP(-YVAL)
      X1MIN = -T1/(S+U1)
      X1 = X1MIN ** Y3
      S4MAX = X1 * (S+U1) +T1
      S4 = (S4MAX -DEL)* Y4 +DEL
      X2 = (S4 -X1*U1)/(X1*S +T1)
      X2B= -X1*U1/(X1*S +T1)
      SHAT   = X1 *X2 *S
      SHATB  = X1 *X2B*S
      T1HAT  = X2 *T1
      T1HATB = X2B*T1
      U1HAT  = X1 *U1
      U1HATB = X1 *U1
      PREF  = S * PREPT * PREY *LOG(X1MIN)*(-X1*X2 )
     +     *(S4MAX -DEL)/(S4-X1*U1)
      PREFB = S * PREPT * PREY *LOG(X1MIN)*(-X1*X2B)/(-X1*U1)
      CALL PDFNLOST1(X1,X2,X2B)
      RETURN
      END

C ======================================================================

      SUBROUTINE DTONLST()
C***  DEFINES THE KINEMATICS FOR TOTAL CROSS-SECTIONS IN NLO
C***  FILLS UP THE COMMON BLOCKS DINLST1-2 AND PDNLSTT
      IMPLICIT REAL*8 (A-H,M-Z)
      COMMON/MINLST/DEL,SCA
      COMMON/CONST6/S,ENERGY,ALPHAS,MST1,MG,MT,MS,MST2,SIN2T
      COMMON/VARVEGAS/Y1,Y2,Y3,Y4
      COMMON/DINLST1/SHATB,T1HATB,U1HATB,PREFB,X1,X2B
      COMMON/DINLST2/SHAT,T1HAT,U1HAT,S4,PREF,X2,S4MAX
      COMMON/PDNLSTT/PDGGT, PDQBT, PDGQT
      X1 = (4*MST1**2/S) ** Y1
      X2 = (4*MST1**2/S/X1) ** Y2
      SHAT = X1*X2*S
      DEL = 1.D-4 * MST1**2 * (1D0 -4*MST1**2/SHAT)
      BETA = SQRT(1D0 -4*MST1**2/SHAT)
      T1HAT = 
     + -(SHAT -DEL)/2D0 +(Y3 -0.5D0)*SQRT((SHAT-DEL)**2-4*MST1**2*SHAT)
      S4MAX = SHAT +T1HAT +SHAT*MST1**2/T1HAT
      S4 = (S4MAX -DEL)*Y4 +DEL
      PREF  = SHAT*BETA *S4MAX*LOG(4*MST1**2/S)*LOG(4*MST1**2/S/X1)
      PREFB = SHAT*BETA *      LOG(4*MST1**2/S)*LOG(4*MST1**2/S/X1)
      CALL PDFNLOSTT(X1,X2)
      RETURN
      END

C ======================================================================

      REAL*8 FUNCTION ESTGGNLO()
C***  CALCULATES THE DIFFERENTIAL CROSS-SECTION IN NLO FOR G G
      IMPLICIT REAL*8 (A-H,M-Z)
      COMMON/MINLST/DEL,SCA
      COMMON/PDLOST/PDGGB,PDQBB,PDBQB
      COMMON/PDNLST1/PDGG,PDQB,PDBQ,PDQG, PDGQ,PDBG, PDGB
      COMMON/CONST6/S,ENERGY,ALPHAS,MST1,MG,MT,MS,MST2,SIN2T
      COMMON/DINLST1/SHATB,T1HATB,U1HATB,PREFB,X1,X2B
      COMMON/DINLST2/SHAT,T1HAT,U1HAT,S4,PREF,X2,S4MAX

      ESTGGNLO = 
     +     +PDGGB * PREFB * DSTGGB(ALPHAS,SHATB,T1HATB,MST1,MG) 
     +     +PDGGB * PREFB * 
     +     DSTGGV(ALPHAS,SHATB,T1HATB,MST1,MG,MT,MS,MST2,sin2t)
     +     +PDGGB * PREFB * DSTGG1(ALPHAS,SHATB,T1HATB,MST1,MG,SCA) 
     +     +PDGGB * PREFB * 
     +     DSTGGD(ALPHAS,SHATB,T1HATB,S4,MST1,MG,DEL,S4MAX)
     +     +PDGGB * PREFB * 
     +     DSTGG2(ALPHAS,SHATB,T1HATB,S4,MST1,MG,DEL,S4MAX,SCA)
     +     +PDGG  * PREF  * DSTGGH(ALPHAS,SHAT,T1HAT,S4,MST1,MG) 
     +     +PDGG  * PREF  * DSTGG3(ALPHAS,SHAT,T1HAT,S4,MST1,MG,SCA) 
      RETURN
      END

C ======================================================================

      REAL*8 FUNCTION ESTQBNLO()
C***  CALCULATES THE DIFFERENTIAL CROSS-SECTION IN NLO FOR Q QB
      IMPLICIT REAL*8 (A-H,M-Z)
      INTEGER ICHARCONJ
      COMMON/MINLST/DEL,SCA
      COMMON/PDLOST/PDGGB,PDQBB,PDBQB
      COMMON/PDNLST1/PDGG,PDQB,PDBQ,PDQG, PDGQ,PDBG, PDGB
      COMMON/CONST6/S,ENERGY,ALPHAS,MST1,MG,MT,MS,MST2,SIN2T
      COMMON/DINLST1/SHATB,T1HATB,U1HATB,PREFB,X1,X2B
      COMMON/DINLST2/SHAT,T1HAT,U1HAT,S4,PREF,X2,S4MAX
      COMMON/CHARCONJ/ICHARCONJ

C***  DISTINGUISH BETWEEN STOPS AND ANTISTOPS IN THE FINAL STATE
C***  ICHARCONJ =  1  <-- STOPS
C***  ICHARCONJ = -1  <-- ANTISTOPS
C***  ICHARCONJ =  0  <-- AVERAGE OF STOPS AND ANTISTOPS

      IF (ICHARCONJ.EQ.1) THEN
         CALL SWAP_ST(PDQBB, PDBQB)
         CALL SWAP_ST(PDQB , PDBQ )
      ELSE IF (ICHARCONJ.EQ.0) THEN
         CALL AVER_ST(PDQBB, PDBQB)
         CALL AVER_ST(PDQB , PDBQ )
      END IF

      ESTQBNLO = 
     +        +PDQBB * PREFB * DSTQBB(ALPHAS,SHATB,T1HATB,MST1,MG) 
     +        +PDQBB * PREFB * 
     +        DSTQBV(ALPHAS,SHATB,T1HATB,MST1,MG,MT,MS,MST2,SIN2T)
     +        +PDQBB * PREFB * DSTQB1(ALPHAS,SHATB,T1HATB,MST1,MG,SCA) 
     +        +PDQBB * PREFB * 
     +        DSTQBD(ALPHAS,SHATB,T1HATB,S4,MST1,MG,DEL,S4MAX)
     +        +PDQBB * PREFB * 
     +        DSTQB2(ALPHAS,SHATB,T1HATB,S4,MST1,MG,DEL,S4MAX,SCA)
     +        +PDQB  * PREF  * 
     +        DSTQBH(ALPHAS,SHAT,T1HAT,S4,MST1,MG) 
     +        +PDQB  * PREF  * 
     +        DSTQB3(ALPHAS,SHAT,T1HAT,S4,MST1,MG,SCA) 

      ESTQBNLO = ESTQBNLO
     +        +PDBQB * PREFB * DSTQBB(ALPHAS,SHATB,U1HATB,MST1,MG) 
     +        +PDBQB * PREFB * 
     +        DSTQBV(ALPHAS,SHATB,U1HATB,MST1,MG,MT,MS,MST2,SIN2T)
     +        +PDBQB * PREFB * DSTQB1(ALPHAS,SHATB,U1HATB,MST1,MG,SCA) 
     +        +PDBQB * PREFB * 
     +        DSTQBD(ALPHAS,SHATB,U1HATB,S4,MST1,MG,DEL,S4MAX)
     +        +PDBQB * PREFB * 
     +        DSTQB2(ALPHAS,SHATB,U1HATB,S4,MST1,MG,DEL,S4MAX,SCA)
     +        +PDBQ  * PREF  * 
     +        DSTQBH(ALPHAS,SHAT,U1HAT,S4,MST1,MG) 
     +        +PDBQ  * PREF  * 
     +        DSTQB3(ALPHAS,SHAT,U1HAT,S4,MST1,MG,SCA) 
      RETURN
      END

C ======================================================================

      REAL*8 FUNCTION ESTGBNLO()
C***  CALCULATES THE DIFFERENTIAL CROSS-SECTION IN NLO FOR G QB AND G Q
      IMPLICIT REAL*8 (A-H,M-Z)
      COMMON/MINLST/DEL,SCA
      COMMON/PDLOST/PDGGB,PDQBB,PDBQB
      COMMON/PDNLST1/PDGG,PDQB,PDBQ,PDQG, PDGQ,PDBG, PDGB
      COMMON/CONST6/S,ENERGY,ALPHAS,MST1,MG,MT,MS,MST2,SIN2T
      COMMON/DINLST1/SHATB,T1HATB,U1HATB,PREFB,X1,X2B
      COMMON/DINLST2/SHAT,T1HAT,U1HAT,S4,PREF,X2,S4MAX
      COMMON/CHARCONJ/ICHARCONJ

C***  DISTINGUISH BETWEEN STOPS AND ANTISTOPS IN THE FINAL STATE
C***  ICHARCONJ =  1  <-- STOPS
C***  ICHARCONJ = -1  <-- ANTISTOPS
C***  ICHARCONJ =  0  <-- AVERAGE OF STOPS AND ANTISTOPS

      IF (ICHARCONJ.EQ.1) THEN
         CALL SWAP_ST(PDGB, PDGQ)
         CALL SWAP_ST(PDBG, PDQG)
      ELSE IF (ICHARCONJ.EQ.0) THEN
         CALL AVER_ST(PDGB, PDGQ)
         CALL AVER_ST(PDBG, PDQG)
      END IF

      ESTGBNLO = 
     +        +PDGB  * PREF  * 
     +        DSTGBH(ALPHAS,SHAT,T1HAT,S4,MST1,MG) 
     +        +PDGB  * PREF  * 
     +        DSTGB3(ALPHAS,SHAT,T1HAT,S4,MST1,MG,SCA) 
     +        +PDBG  * PREF  * 
     +        DSTGBH(ALPHAS,SHAT,U1HAT,S4,MST1,MG) 
     +        +PDBG  * PREF  * 
     +        DSTGB3(ALPHAS,SHAT,U1HAT,S4,MST1,MG,SCA) 
     +        
     +        +PDQG  * PREF  * 
     +        DSTQGH(ALPHAS,SHAT,T1HAT,S4,MST1,MG) 
     +        +PDQG  * PREF  * 
     +        DSTQG3(ALPHAS,SHAT,T1HAT,S4,MST1,MG,SCA) 
     +        +PDGQ  * PREF  * 
     +        DSTQGH(ALPHAS,SHAT,U1HAT,S4,MST1,MG) 
     +        +PDGQ  * PREF  * 
     +        DSTQG3(ALPHAS,SHAT,U1HAT,S4,MST1,MG,SCA) 

      RETURN
      END

C ======================================================================

      REAL*8 FUNCTION FSTGGNLO()
C***  CALCULATES THE TOTAL CROSS-SECTION IN NLO FOR G G
      IMPLICIT REAL*8 (A-H,M-Z)
      COMMON/CONST6/S,ENERGY,ALPHAS,MST1,MG,MT,MS,MST2,SIN2T
      COMMON/PDNLSTT/PDGGT, PDQBT, PDGQT
      COMMON/MINLST/DEL,SCA
      COMMON/DINLST1/SHATB,T1HATB,U1HATB,PREFB,X1,X2B
      COMMON/DINLST2/SHAT,T1HAT,U1HAT,S4,PREF,X2,S4MAX

      FSTGGNLO = 
     +        +PDGGT * PREFB * DSTGGB(ALPHAS,SHAT,T1HAT,MST1,MG) 
     +        +PDGGT * PREFB * 
     +        DSTGGV(ALPHAS,SHAT,T1HAT,MST1,MG,MT,MS,MST2,SIN2T)
     +        +PDGGT * PREFB * DSTGG1(ALPHAS,SHAT,T1HAT,MST1,MG,SCA) 
     +        +PDGGT * PREFB * 
     +        DSTGGD(ALPHAS,SHAT,T1HAT,S4,MST1,MG,DEL,S4MAX)
     +        +PDGGT * PREFB * 
     +        DSTGG2(ALPHAS,SHAT,T1HAT,S4,MST1,MG,DEL,S4MAX,SCA)
     +        +PDGGT * PREF  * DSTGGH(ALPHAS,SHAT,T1HAT,S4,MST1,MG) 
     +        +PDGGT * PREF  * DSTGG3(ALPHAS,SHAT,T1HAT,S4,MST1,MG,SCA) 
      RETURN
      END

C ======================================================================

      REAL*8 FUNCTION FSTQBNLO()
C***  CALCULATES THE TOTAL CROSS-SECTION IN NLO FOR Q QB
      IMPLICIT REAL*8 (A-H,M-Z)
      COMMON/PDNLSTT/PDGGT, PDQBT, PDGQT
      COMMON/CONST6/S,ENERGY,ALPHAS,MST1,MG,MT,MS,MST2,SIN2T
      COMMON/MINLST/DEL,SCA
      COMMON/DINLST1/SHATB,T1HATB,U1HATB,PREFB,X1,X2B
      COMMON/DINLST2/SHAT,T1HAT,U1HAT,S4,PREF,X2,S4MAX

      FSTQBNLO = 
     +        +PDQBT * PREFB * DSTQBB(ALPHAS,SHAT,T1HAT,MST1,MG) 
     +        +PDQBT * PREFB * 
     +        DSTQBV(ALPHAS,SHAT,T1HAT,MST1,MG,MT,MS,MST2,SIN2T)
     +        +PDQBT * PREFB * DSTQB1(ALPHAS,SHAT,T1HAT,MST1,MG,SCA) 
     +        +PDQBT * PREFB * 
     +        DSTQBD(ALPHAS,SHAT,T1HAT,S4,MST1,MG,DEL,S4MAX)
     +        +PDQBT * PREFB * 
     +        DSTQB2(ALPHAS,SHAT,T1HAT,S4,MST1,MG,DEL,S4MAX,SCA)
     +        +PDQBT * PREF  * 
     +        DSTQBH(ALPHAS,SHAT,T1HAT,S4,MST1,MG) 
     +        +PDQBT * PREF  * 
     +        DSTQB3(ALPHAS,SHAT,T1HAT,S4,MST1,MG,SCA) 
      RETURN
      END

C ======================================================================

      REAL*8 FUNCTION FSTGBNLO()
C***  CALCULATES THE TOTAL CROSS-SECTION IN NLO FOR G QB AND G Q 
      IMPLICIT REAL*8 (A-H,M-Z)
      COMMON/CONST6/S,ENERGY,ALPHAS,MST1,MG,MT,MS,MST2,SIN2T
      COMMON/PDNLSTT/PDGGT, PDQBT, PDGQT
      COMMON/MINLST/DEL,SCA
      COMMON/DINLST1/SHATB,T1HATB,U1HATB,PREFB,X1,X2B
      COMMON/DINLST2/SHAT,T1HAT,U1HAT,S4,PREF,X2,S4MAX

      FSTGBNLO = 
     +        +PDGQT * PREF  * 
     +        DSTGBH(ALPHAS,SHAT,T1HAT,S4,MST1,MG) 
     +        +PDGQT * PREF  * 
     +        DSTGB3(ALPHAS,SHAT,T1HAT,S4,MST1,MG,SCA) 
      RETURN
      END

C ======================================================================

      REAL*8 FUNCTION SUMNLOST(SIGNLOST1,SIGNLOST2,SIGNLOST3,IFLAVOR)
C***  SUMS UP DIFFERENT FLAVOR STATES FOR CROSS-SECTION IN NLO
      IMPLICIT REAL*8 (A-H,M-Z)
      IMPLICIT INTEGER (I)
      IF (IFLAVOR.EQ.0) SUMNLOST = SIGNLOST1 +SIGNLOST2 +SIGNLOST3
      IF (IFLAVOR.EQ.1) SUMNLOST = SIGNLOST1
      IF (IFLAVOR.EQ.2) SUMNLOST = SIGNLOST2
      IF (IFLAVOR.EQ.3) SUMNLOST = SIGNLOST3
      IF (IFLAVOR.EQ.4) SUMNLOST = SIGNLOST1 +SIGNLOST2
      IF (IFLAVOR.EQ.5) SUMNLOST = SIGNLOST3
      IF (.NOT.(ABS(SUMNLOST).LT.1.D40)) THEN
         SUMNLOST = 0.D0
         PRINT *,'NAN'
      END IF
      RETURN
      END
      
C ======================================================================

      SUBROUTINE PDFLOST(X1,X2)
      IMPLICIT REAL*8 (A-H,M-Z)
      IMPLICIT INTEGER (I,J)
      DIMENSION XPDF(-6:6)
      DIMENSION YPDF(-6:6)
      COMMON/CONST2/SCALE,SCAFAC,ICOLL,ISCAPT
      COMMON/PDLOST/PDGGB,PDQBB,PDBQB

      CALL PARTONDF(X1,SCALE,XPDF)
      CALL PARTONDF(X2,SCALE,YPDF)

      PDGGB = XPDF(0) * YPDF(0)
      PDQBB = 0D0
      PDBQB = 0D0
      IF (ICOLL.EQ.0) THEN
         DO 1000 I= 1,5
            PDQBB = PDQBB + XPDF( I) * YPDF( I) 
            PDBQB = PDBQB + XPDF(-I) * YPDF(-I) 
 1000    CONTINUE
      ELSE IF (ICOLL.EQ.1) THEN
         DO 1010 I= 1,5
            PDQBB = PDQBB + XPDF( I) * YPDF(-I) 
            PDBQB = PDBQB + XPDF(-I) * YPDF( I) 
 1010    CONTINUE
      END IF
      RETURN
      END

C ======================================================================

      SUBROUTINE PDFNLOSTT(X1,X2)
      IMPLICIT REAL*8 (A-H,M-Z)
      IMPLICIT INTEGER (I,J)
      DIMENSION XPDF(-6:6)
      DIMENSION YPDF(-6:6)
      COMMON/CONST2/SCALE,SCAFAC,ICOLL,ISCAPT
      COMMON/PDNLSTT/PDGGT, PDQBT, PDGQT

      CALL PARTONDF(X1,SCALE,XPDF)
      CALL PARTONDF(X2,SCALE,YPDF)

      PDGGT = XPDF(0) * YPDF(0)
      PDQBT = 0D0
      PDGQT = 0D0

      IF (ICOLL.EQ.0) THEN
         DO 1000 I= 1,5
            PDQBT = PDQBT + XPDF( I) * YPDF( I) + XPDF(-I) * YPDF(-I) 
            PDGQT = PDGQT + XPDF(0)  * YPDF(-I) + XPDF( I) * YPDF( 0) 
     +                    + XPDF(0)  * YPDF(I)  + XPDF(-I) * YPDF( 0) 
 1000    CONTINUE
      ELSE IF (ICOLL.EQ.1) THEN
         DO 1010 I= 1,5
            PDQBT = PDQBT + XPDF( I) * YPDF(-I) + XPDF(-I) * YPDF(I) 
            PDGQT = PDGQT + XPDF(0)  * YPDF(-I) + XPDF( I) * YPDF(0) 
     +                    + XPDF(0)  * YPDF(I)  + XPDF(-I) * YPDF(0) 
 1010    CONTINUE
      END IF
      RETURN
      END

C ======================================================================

      SUBROUTINE PDFNLOST1(X1,X2,X2B)
      IMPLICIT REAL*8 (A-H,M-Z)
      IMPLICIT INTEGER (I,J)
      DIMENSION XPDF(-6:6)
      DIMENSION YPDF(-6:6)
      DIMENSION ZPDF(-6:6)
      COMMON/CONST2/SCALE,SCAFAC,ICOLL,ISCAPT
      COMMON/PDLOST/PDGGB,PDQBB,PDBQB
      COMMON/PDNLST1/PDGG,PDQB,PDBQ,PDQG, PDGQ,PDBG, PDGB

      CALL PARTONDF(X1 ,SCALE,XPDF)
      CALL PARTONDF(X2B,SCALE,YPDF)
      CALL PARTONDF(X2 ,SCALE,ZPDF)

      PDGGB = XPDF(0) * YPDF(0)
      PDGG  = XPDF(0) * ZPDF(0)
      PDQBB = 0D0
      PDQB  = 0D0
      PDBQB = 0D0
      PDBQ  = 0D0
      PDQG  = 0D0
      PDGQ  = 0D0
      PDBG  = 0D0
      PDGB  = 0D0

      IF (ICOLL.EQ.0) THEN
         DO 1000 I= 1,5
            PDQBB = PDQBB + XPDF( I) * YPDF( I) 
            PDQB  = PDQB  + XPDF( I) * ZPDF( I) 
            PDBQB = PDBQB + XPDF(-I) * YPDF(-I) 
            PDBQ  = PDBQ  + XPDF(-I) * ZPDF(-I) 
            PDQG  = PDQG  + XPDF(+I) * ZPDF( 0) 
            PDGQ  = PDGQ  + XPDF( 0) * ZPDF(-I) 
            PDBG  = PDBG  + XPDF(-I) * ZPDF( 0) 
            PDGB  = PDGB  + XPDF( 0) * ZPDF( I) 
            
 1000    CONTINUE
      ELSE IF (ICOLL.EQ.1) THEN
         DO 1010 I= 1,5
            PDQBB = PDQBB + XPDF( I) * YPDF(-I) 
            PDQB  = PDQB  + XPDF( I) * ZPDF(-I) 
            PDBQB = PDBQB + XPDF(-I) * YPDF( I) 
            PDBQ  = PDBQ  + XPDF(-I) * ZPDF( I) 
            PDQG  = PDQG  + XPDF( I) * ZPDF( 0) 
            PDGQ  = PDGQ  + XPDF( 0) * ZPDF( I) 
            PDBG  = PDBG  + XPDF(-I) * ZPDF( 0) 
            PDGB  = PDGB  + XPDF( 0) * ZPDF(-I) 
 1010    CONTINUE
      END IF
      RETURN
      END

C ======================================================================

