pythiaeRHIC/pythia/pysghg.f

0001
0002 C*********************************************************************
0003
0004 C...PYSGHG
0005 C...Subprocess cross sections for Higgs processes,
0006 C...except Higgs pairs in PYSGSU, but including WW scattering.
0007 C...Auxiliary to PYSIGH.
0008
0009       SUBROUTINE PYSGHG(NCHN,SIGS)
0010
0011 C...Double precision and integer declarations
0012       IMPLICIT DOUBLE PRECISION(A-H, O-Z)
0013       IMPLICIT INTEGER(I-N)
0014       INTEGER PYK,PYCHGE,PYCOMP
0015 C...Parameter statement to help give large particle numbers.
0016       PARAMETER (KSUSY1=1000000,KSUSY2=2000000,KTECHN=3000000,
0017      &KEXCIT=4000000,KDIMEN=5000000)
0018 C...Commonblocks
0019       COMMON/PYDAT1/MSTU(200),PARU(200),MSTJ(200),PARJ(200)
0020       COMMON/PYDAT2/KCHG(500,4),PMAS(500,4),PARF(2000),VCKM(4,4)
0021       COMMON/PYDAT3/MDCY(500,3),MDME(8000,2),BRAT(8000),KFDP(8000,5)
0022       COMMON/PYPARS/MSTP(200),PARP(200),MSTI(200),PARI(200)
0023       COMMON/PYINT1/MINT(400),VINT(400)
0024       COMMON/PYINT2/ISET(500),KFPR(500,2),COEF(500,20),ICOL(40,4,2)
0025       COMMON/PYINT3/XSFX(2,-40:40),ISIG(1000,3),SIGH(1000)
0026       COMMON/PYINT4/MWID(500),WIDS(500,5)
0027       COMMON/PYSUBS/MSEL,MSELPD,MSUB(500),KFIN(2,-40:40),CKIN(200)
0028       COMMON/PYMSSM/IMSS(0:99),RMSS(0:99)
0029       COMMON/PYSGCM/ISUB,ISUBSV,MMIN1,MMAX1,MMIN2,MMAX2,MMINA,MMAXA,
0030      &KFAC(2,-40:40),COMFAC,FACK,FACA,SH,TH,UH,SH2,TH2,UH2,SQM3,SQM4,
0031      &SHR,SQPTH,TAUP,BE34,CTH,X(2),SQMZ,SQMW,GMMZ,GMMW,
0032      &AEM,AS,XW,XW1,XWC,XWV,POLL,POLR,POLLL,POLRR
0033       SAVE /PYDAT1/,/PYDAT2/,/PYDAT3/,/PYPARS/,/PYINT1/,/PYINT2/,
0034      &/PYINT3/,/PYINT4/,/PYSUBS/,/PYMSSM/,/PYSGCM/
0035 C...Local arrays and complex variables
0036       DIMENSION WDTP(0:400),WDTE(0:400,0:5)
0037       COMPLEX*16 A004,A204,A114,A00U,A20U,A11U
0038       COMPLEX*16 CIGTOT,CIZTOT,F0ALP,F1ALP,F2ALP,F0BET,F1BET,F2BET,FIF
0039
0040 C...Convert H or A process into equivalent h one
0041       IHIGG=1
0042       KFHIGG=25
0043       IF(ISUB.EQ.401.OR.ISUB.EQ.402) THEN
0044          KFHIGG=KFPR(ISUB,1)
0045       END IF
0046       IF((ISUB.GE.151.AND.ISUB.LE.160).OR.(ISUB.GE.171.AND.
0047      &ISUB.LE.190)) THEN
0048         IHIGG=2
0049         IF(MOD(ISUB-1,10).GE.5) IHIGG=3
0050         KFHIGG=33+IHIGG
0051         IF(ISUB.EQ.151.OR.ISUB.EQ.156) ISUB=3
0052         IF(ISUB.EQ.152.OR.ISUB.EQ.157) ISUB=102
0053         IF(ISUB.EQ.153.OR.ISUB.EQ.158) ISUB=103
0054         IF(ISUB.EQ.171.OR.ISUB.EQ.176) ISUB=24
0055         IF(ISUB.EQ.172.OR.ISUB.EQ.177) ISUB=26
0056         IF(ISUB.EQ.173.OR.ISUB.EQ.178) ISUB=123
0057         IF(ISUB.EQ.174.OR.ISUB.EQ.179) ISUB=124
0058         IF(ISUB.EQ.181.OR.ISUB.EQ.186) ISUB=121
0059         IF(ISUB.EQ.182.OR.ISUB.EQ.187) ISUB=122
0060         IF(ISUB.EQ.183.OR.ISUB.EQ.188) ISUB=111
0061         IF(ISUB.EQ.184.OR.ISUB.EQ.189) ISUB=112
0062         IF(ISUB.EQ.185.OR.ISUB.EQ.190) ISUB=113
0063       ENDIF
0064       SQMH=PMAS(KFHIGG,1)**2
0065       GMMH=PMAS(KFHIGG,1)*PMAS(KFHIGG,2)
0066
0067 C...Strongly interacting Z_L/W_L model of Dobado, Herrero, Terron
0068       IF((MSTP(46).GE.3.AND.MSTP(46).LE.6).AND.(ISUB.EQ.71.OR.ISUB.EQ.
0069      &72.OR.ISUB.EQ.73.OR.ISUB.EQ.76.OR.ISUB.EQ.77)) THEN
0070 C...Calculate M_R and N_R functions for Higgs-like and QCD-like models
0071         IF(MSTP(46).LE.4) THEN
0072           HDTLH=LOG(PMAS(25,1)/PARP(44))
0073           HDTMR=(4.5D0*PARU(1)/SQRT(3D0)-74D0/9D0)/8D0+HDTLH/12D0
0074           HDTNR=-1D0/18D0+HDTLH/6D0
0075         ELSE
0076           HDTNM=0.125D0*(1D0/(288D0*PARU(1)**2)+(PARP(47)/PARP(45))**2)
0077           HDTLQ=LOG(PARP(45)/PARP(44))
0078           HDTMR=-(4D0*PARU(1))**2*0.5D0*HDTNM+HDTLQ/12D0
0079           HDTNR=(4D0*PARU(1))**2*HDTNM+HDTLQ/6D0
0080         ENDIF
0081
0082 C...Calculate lowest and next-to-lowest order partial wave amplitudes
0083         HDTV=1D0/(16D0*PARU(1)*PARP(47)**2)
0084         A00L=DBLE(HDTV*SH)
0085         A20L=-0.5D0*A00L
0086         A11L=A00L/6D0
0087         HDTLS=LOG(SH/PARP(44)**2)
0088         A004=DBLE((HDTV*SH)**2/(4D0*PARU(1)))*
0089      &  CMPLX(DBLE((176D0*HDTMR+112D0*HDTNR)/3D0+11D0/27D0-
0090      &  (50D0/9D0)*HDTLS),DBLE(4D0*PARU(1)))
0091         A204=DBLE((HDTV*SH)**2/(4D0*PARU(1)))*
0092      &  CMPLX(DBLE(32D0*(HDTMR+2D0*HDTNR)/3D0+25D0/54D0-
0093      &  (20D0/9D0)*HDTLS),DBLE(PARU(1)))
0094         A114=DBLE((HDTV*SH)**2/(6D0*PARU(1)))*
0095      &  CMPLX(DBLE(4D0*(-2D0*HDTMR+HDTNR)-1D0/18D0),DBLE(PARU(1)/6D0))
0096
0097 C...Unitarize partial wave amplitudes with Pade or K-matrix method
0098         IF(MSTP(46).EQ.3.OR.MSTP(46).EQ.5) THEN
0099           A00U=A00L/(1D0-A004/A00L)
0100           A20U=A20L/(1D0-A204/A20L)
0101           A11U=A11L/(1D0-A114/A11L)
0102         ELSE
0103           A00U=(A00L+DBLE(A004))/(1D0-DCMPLX(0.D0,A00L+DBLE(A004)))
0104           A20U=(A20L+DBLE(A204))/(1D0-DCMPLX(0.D0,A20L+DBLE(A204)))
0105           A11U=(A11L+DBLE(A114))/(1D0-DCMPLX(0.D0,A11L+DBLE(A114)))
0106         ENDIF
0107       ENDIF
0108
0109 C...Differential cross section expressions.
0110
0111       IF(ISUB.LE.60) THEN
0112         IF(ISUB.EQ.3) THEN
0113 C...f + fbar -> h0 (or H0, or A0)
0114           CALL PYWIDT(KFHIGG,SH,WDTP,WDTE)
0115           HS=SHR*WDTP(0)
0116           FACBW=4D0*COMFAC/((SH-SQMH)**2+HS**2)
0117           IF(ABS(SHR-PMAS(KFHIGG,1)).GT.PARP(48)*PMAS(KFHIGG,2))
0118      &    FACBW=0D0
0119           HP=AEM/(8D0*XW)*SH/SQMW*SH
0120           HF=SHR*(WDTE(0,1)+WDTE(0,2)+WDTE(0,4))
0121           DO 100 I=MMINA,MMAXA
0122             IF(I.EQ.0.OR.KFAC(1,I)*KFAC(2,-I).EQ.0) GOTO 100
0123             IA=IABS(I)
0124             RMQ=PYMRUN(IA,SH)**2/SH
0125             HI=HP*RMQ
0126             IF(IA.LE.10) HI=HP*RMQ*FACA/3D0
0127             IF(MSTP(4).GE.1.OR.IHIGG.GE.2) THEN
0128               IKFI=1
0129               IF(IA.LE.10.AND.MOD(IA,2).EQ.0) IKFI=2
0130               IF(IA.GT.10) IKFI=3
0131               HI=HI*PARU(150+10*IHIGG+IKFI)**2
0132               IF(IMSS(1).NE.0.AND.IA.EQ.5) THEN
0133                 HI=HI/(1D0+RMSS(41))**2
0134                 IF(IHIGG.NE.3) THEN
0135                   HI=HI*(1D0+RMSS(41)*PARU(152+10*IHIGG)/
0136      &            PARU(151+10*IHIGG))**2
0137                 ENDIF
0138               ENDIF
0139             ENDIF
0140             NCHN=NCHN+1
0141             ISIG(NCHN,1)=I
0142             ISIG(NCHN,2)=-I
0143             ISIG(NCHN,3)=1
0144             SIGH(NCHN)=HI*FACBW*HF
0145   100     CONTINUE
0146
0147         ELSEIF(ISUB.EQ.5) THEN
0148 C...Z0 + Z0 -> h0
0149           CALL PYWIDT(25,SH,WDTP,WDTE)
0150           HS=SHR*WDTP(0)
0151           FACBW=4D0*COMFAC/((SH-SQMH)**2+HS**2)
0152           IF(ABS(SHR-PMAS(25,1)).GT.PARP(48)*PMAS(25,2)) FACBW=0D0
0153           HP=AEM/(8D0*XW)*SH/SQMW*SH
0154           HF=SHR*(WDTE(0,1)+WDTE(0,2)+WDTE(0,4))
0155           HI=HP/4D0
0156           FACI=8D0/(PARU(1)**2*XW1)*(AEM*XWC)**2
0157           DO 120 I=MMIN1,MMAX1
0158             IF(I.EQ.0.OR.KFAC(1,I).EQ.0) GOTO 120
0159             DO 110 J=MMIN2,MMAX2
0160               IF(J.EQ.0.OR.KFAC(2,J).EQ.0) GOTO 110
0161               EI=KCHG(IABS(I),1)/3D0
0162               AI=SIGN(1D0,EI)
0163               VI=AI-4D0*EI*XWV
0164               EJ=KCHG(IABS(J),1)/3D0
0165               AJ=SIGN(1D0,EJ)
0166               VJ=AJ-4D0*EJ*XWV
0167               NCHN=NCHN+1
0168               ISIG(NCHN,1)=I
0169               ISIG(NCHN,2)=J
0170               ISIG(NCHN,3)=1
0171               SIGH(NCHN)=FACI*(VI**2+AI**2)*(VJ**2+AJ**2)*HI*FACBW*HF
0172   110       CONTINUE
0173   120     CONTINUE
0174
0175         ELSEIF(ISUB.EQ.8) THEN
0176 C...W+ + W- -> h0
0177           CALL PYWIDT(25,SH,WDTP,WDTE)
0178           HS=SHR*WDTP(0)
0179           FACBW=4D0*COMFAC/((SH-SQMH)**2+HS**2)
0180           IF(ABS(SHR-PMAS(25,1)).GT.PARP(48)*PMAS(25,2)) FACBW=0D0
0181           HP=AEM/(8D0*XW)*SH/SQMW*SH
0182           HF=SHR*(WDTE(0,1)+WDTE(0,2)+WDTE(0,4))
0183           HI=HP/2D0
0184           FACI=1D0/(4D0*PARU(1)**2)*(AEM/XW)**2
0185           DO 140 I=MMIN1,MMAX1
0186             IF(I.EQ.0.OR.KFAC(1,I).EQ.0) GOTO 140
0187             EI=SIGN(1D0,DBLE(I))*KCHG(IABS(I),1)
0188             DO 130 J=MMIN2,MMAX2
0189               IF(J.EQ.0.OR.KFAC(2,J).EQ.0) GOTO 130
0190               EJ=SIGN(1D0,DBLE(J))*KCHG(IABS(J),1)
0191               IF(EI*EJ.GT.0D0) GOTO 130
0192               NCHN=NCHN+1
0193               ISIG(NCHN,1)=I
0194               ISIG(NCHN,2)=J
0195               ISIG(NCHN,3)=1
0196               SIGH(NCHN)=FACI*VINT(180+I)*VINT(180+J)*HI*FACBW*HF
0197   130       CONTINUE
0198   140     CONTINUE
0199
0200         ELSEIF(ISUB.EQ.24) THEN
0201 C...f + fbar -> Z0 + h0 (or H0, or A0)
0202 C...Propagators: Z0, h0 as simulated in PYOFSH and as desired
0203           HBW3=GMMZ/((SQM3-SQMZ)**2+GMMZ**2)
0204           CALL PYWIDT(23,SQM3,WDTP,WDTE)
0205           GMMZ3=SQRT(SQM3)*WDTP(0)
0206           HBW3C=GMMZ3/((SQM3-SQMZ)**2+GMMZ3**2)
0207           HBW4=GMMH/((SQM4-SQMH)**2+GMMH**2)
0208           CALL PYWIDT(KFHIGG,SQM4,WDTP,WDTE)
0209           GMMH4=SQRT(SQM4)*WDTP(0)
0210           HBW4C=GMMH4/((SQM4-SQMH)**2+GMMH4**2)
0211           THUH=MAX(TH*UH-SQM3*SQM4,SH*CKIN(3)**2)
0212           FACHZ=COMFAC*(HBW3C/HBW3)*(HBW4C/HBW4)*8D0*(AEM*XWC)**2*
0213      &    (THUH+2D0*SH*SQM3)/((SH-SQMZ)**2+GMMZ**2)
0214           FACHZ=FACHZ*WIDS(23,2)*WIDS(KFHIGG,2)
0215           IF(MSTP(4).GE.1.OR.IHIGG.GE.2) FACHZ=FACHZ*
0216      &    PARU(154+10*IHIGG)**2
0217           DO 150 I=MMINA,MMAXA
0218             IF(I.EQ.0.OR.KFAC(1,I)*KFAC(2,-I).EQ.0) GOTO 150
0219             EI=KCHG(IABS(I),1)/3D0
0220             AI=SIGN(1D0,EI)
0221             VI=AI-4D0*EI*XWV
0222             FCOI=1D0
0223             IF(IABS(I).LE.10) FCOI=FACA/3D0
0224             NCHN=NCHN+1
0225             ISIG(NCHN,1)=I
0226             ISIG(NCHN,2)=-I
0227             ISIG(NCHN,3)=1
0228             SIGH(NCHN)=FACHZ*FCOI*(VI**2+AI**2)
0229   150     CONTINUE
0230
0231         ELSEIF(ISUB.EQ.26) THEN
0232 C...f + fbar' -> W+/- + h0 (or H0, or A0)
0233 C...Propagators: W+-, h0 as simulated in PYOFSH and as desired
0234           HBW3=GMMW/((SQM3-SQMW)**2+GMMW**2)
0235           CALL PYWIDT(24,SQM3,WDTP,WDTE)
0236           GMMW3=SQRT(SQM3)*WDTP(0)
0237           HBW3C=GMMW3/((SQM3-SQMW)**2+GMMW3**2)
0238           HBW4=GMMH/((SQM4-SQMH)**2+GMMH**2)
0239           CALL PYWIDT(KFHIGG,SQM4,WDTP,WDTE)
0240           GMMH4=SQRT(SQM4)*WDTP(0)
0241           HBW4C=GMMH4/((SQM4-SQMH)**2+GMMH4**2)
0242           THUH=MAX(TH*UH-SQM3*SQM4,SH*CKIN(3)**2)
0243           FACHW=COMFAC*0.125D0*(AEM/XW)**2*(THUH+2D0*SH*SQM3)/
0244      &    ((SH-SQMW)**2+GMMW**2)*(HBW3C/HBW3)*(HBW4C/HBW4)
0245           FACHW=FACHW*WIDS(KFHIGG,2)
0246           IF(MSTP(4).GE.1.OR.IHIGG.GE.2) FACHW=FACHW*
0247      &    PARU(155+10*IHIGG)**2
0248           DO 170 I=MMIN1,MMAX1
0249             IA=IABS(I)
0250             IF(I.EQ.0.OR.IA.GT.20.OR.KFAC(1,I).EQ.0) GOTO 170
0251             DO 160 J=MMIN2,MMAX2
0252               JA=IABS(J)
0253               IF(J.EQ.0.OR.JA.GT.20.OR.KFAC(1,J).EQ.0) GOTO 160
0254               IF(I*J.GT.0.OR.MOD(IA+JA,2).EQ.0) GOTO 160
0255               IF((IA.LE.10.AND.JA.GT.10).OR.(IA.GT.10.AND.JA.LE.10))
0256      &        GOTO 160
0257               KCHW=(KCHG(IA,1)*ISIGN(1,I)+KCHG(JA,1)*ISIGN(1,J))/3
0258               FCKM=1D0
0259               IF(IA.LE.10) FCKM=VCKM((IA+1)/2,(JA+1)/2)
0260               FCOI=1D0
0261               IF(IA.LE.10) FCOI=FACA/3D0
0262               NCHN=NCHN+1
0263               ISIG(NCHN,1)=I
0264               ISIG(NCHN,2)=J
0265               ISIG(NCHN,3)=1
0266               SIGH(NCHN)=FACHW*FCOI*FCKM*WIDS(24,(5-KCHW)/2)
0267   160       CONTINUE
0268   170     CONTINUE
0269
0270         ELSEIF(ISUB.EQ.32) THEN
0271 C...f + g -> f + h0 (q + g -> q + h0 only)
0272           FHCQ=COMFAC*FACA*AS*AEM/XW*1D0/24D0
0273 C...H propagator: as simulated in PYOFSH and as desired
0274           SQMHC=PMAS(25,1)**2
0275           GMMHC=PMAS(25,1)*PMAS(25,2)
0276           HBW4=GMMHC/((SQM4-SQMHC)**2+GMMHC**2)
0277           CALL PYWIDT(25,SQM4,WDTP,WDTE)
0278           GMMHCC=SQRT(SQM4)*WDTP(0)
0279           HBW4C=GMMHCC/((SQM4-SQMHC)**2+GMMHCC**2)
0280           FHCQ=FHCQ*HBW4C/HBW4
0281           DO 190 I=MMINA,MMAXA
0282             IA=IABS(I)
0283             IF(IA.NE.5) GOTO 190
0284             SQML=PYMRUN(IA,SH)**2
0285             SQMQ=PMAS(IA,1)**2
0286             FACHCQ=FHCQ*SQML/SQMW*
0287      &      (SH/(SQMQ-UH)+2D0*SQMQ*(SQM4-UH)/(SQMQ-UH)**2+(SQMQ-UH)/SH-
0288      &      2D0*SQMQ/(SQMQ-UH)+2D0*(SQM4-UH)/(SQMQ-UH)*
0289      &      (SQM4-SQMQ-SH)/SH)
0290             DO 180 ISDE=1,2
0291               IF(ISDE.EQ.1.AND.KFAC(1,I)*KFAC(2,21).EQ.0) GOTO 180
0292               IF(ISDE.EQ.2.AND.KFAC(1,21)*KFAC(2,I).EQ.0) GOTO 180
0293               NCHN=NCHN+1
0294               ISIG(NCHN,ISDE)=I
0295               ISIG(NCHN,3-ISDE)=21
0296               ISIG(NCHN,3)=1
0297               SIGH(NCHN)=FACHCQ*WIDS(25,2)
0298   180       CONTINUE
0299   190     CONTINUE
0300         ENDIF
0301
0302       ELSEIF(ISUB.LE.80) THEN
0303         IF(ISUB.EQ.71) THEN
0304 C...Z0 + Z0 -> Z0 + Z0
0305           IF(SH.LE.4.01D0*SQMZ) GOTO 220
0306
0307           IF(MSTP(46).LE.2) THEN
0308 C...Exact scattering ME:s for on-mass-shell gauge bosons
0309             BE2=1D0-4D0*SQMZ/SH
0310             TH=-0.5D0*SH*BE2*(1D0-CTH)
0311             UH=-0.5D0*SH*BE2*(1D0+CTH)
0312             IF(MAX(TH,UH).GT.-1D0) GOTO 220
0313             SHANG=1D0/XW1*SQMW/SQMZ*(1D0+BE2)**2
0314             ASHRE=(SH-SQMH)/((SH-SQMH)**2+GMMH**2)*SHANG
0315             ASHIM=-GMMH/((SH-SQMH)**2+GMMH**2)*SHANG
0316             THANG=1D0/XW1*SQMW/SQMZ*(BE2-CTH)**2
0317             ATHRE=(TH-SQMH)/((TH-SQMH)**2+GMMH**2)*THANG
0318             ATHIM=-GMMH/((TH-SQMH)**2+GMMH**2)*THANG
0319             UHANG=1D0/XW1*SQMW/SQMZ*(BE2+CTH)**2
0320             AUHRE=(UH-SQMH)/((UH-SQMH)**2+GMMH**2)*UHANG
0321             AUHIM=-GMMH/((UH-SQMH)**2+GMMH**2)*UHANG
0322             FACZZ=COMFAC*1D0/(4096D0*PARU(1)**2*16D0*XW1**2)*
0323      &      (AEM/XW)**4*(SH/SQMW)**2*(SQMZ/SQMW)*SH2
0324             IF(MSTP(46).LE.0) FACZZ=FACZZ*(ASHRE**2+ASHIM**2)
0325             IF(MSTP(46).EQ.1) FACZZ=FACZZ*((ASHRE+ATHRE+AUHRE)**2+
0326      &      (ASHIM+ATHIM+AUHIM)**2)
0327             IF(MSTP(46).EQ.2) FACZZ=0D0
0328
0329           ELSE
0330 C...Strongly interacting Z_L/W_L model of Dobado, Herrero, Terron
0331             FACZZ=COMFAC*(AEM/(16D0*PARU(1)*XW*XW1))**2*(64D0/9D0)*
0332      &      ABS(A00U+2D0*A20U)**2
0333           ENDIF
0334           FACZZ=FACZZ*WIDS(23,1)
0335
0336           DO 210 I=MMIN1,MMAX1
0337             IF(I.EQ.0.OR.KFAC(1,I).EQ.0) GOTO 210
0338             EI=KCHG(IABS(I),1)/3D0
0339             AI=SIGN(1D0,EI)
0340             VI=AI-4D0*EI*XWV
0341             AVI=AI**2+VI**2
0342             DO 200 J=MMIN2,MMAX2
0343               IF(J.EQ.0.OR.KFAC(2,J).EQ.0) GOTO 200
0344               EJ=KCHG(IABS(J),1)/3D0
0345               AJ=SIGN(1D0,EJ)
0346               VJ=AJ-4D0*EJ*XWV
0347               AVJ=AJ**2+VJ**2
0348               NCHN=NCHN+1
0349               ISIG(NCHN,1)=I
0350               ISIG(NCHN,2)=J
0351               ISIG(NCHN,3)=1
0352               SIGH(NCHN)=0.5D0*FACZZ*AVI*AVJ
0353   200       CONTINUE
0354   210     CONTINUE
0355   220     CONTINUE
0356
0357         ELSEIF(ISUB.EQ.72) THEN
0358 C...Z0 + Z0 -> W+ + W-
0359           IF(SH.LE.4.01D0*SQMZ) GOTO 250
0360
0361           IF(MSTP(46).LE.2) THEN
0362 C...Exact scattering ME:s for on-mass-shell gauge bosons
0363             BE2=SQRT((1D0-4D0*SQMW/SH)*(1D0-4D0*SQMZ/SH))
0364             CTH2=CTH**2
0365             TH=-0.5D0*SH*(1D0-2D0*(SQMW+SQMZ)/SH-BE2*CTH)
0366             UH=-0.5D0*SH*(1D0-2D0*(SQMW+SQMZ)/SH+BE2*CTH)
0367             IF(MAX(TH,UH).GT.-1D0) GOTO 250
0368             SHANG=4D0*SQRT(SQMW/(SQMZ*XW1))*(1D0-2D0*SQMW/SH)*
0369      &      (1D0-2D0*SQMZ/SH)
0370             ASHRE=(SH-SQMH)/((SH-SQMH)**2+GMMH**2)*SHANG
0371             ASHIM=-GMMH/((SH-SQMH)**2+GMMH**2)*SHANG
0372             ATWRE=XW1/SQMZ*SH/(TH-SQMW)*((CTH-BE2)**2*(3D0/2D0+BE2/2D0*
0373      &      CTH-(SQMW+SQMZ)/SH+(SQMW-SQMZ)**2/(SH*SQMW))+4D0*
0374      &      ((SQMW+SQMZ)/SH*(1D0-3D0*CTH2)+8D0*SQMW*SQMZ/SH2*
0375      &      (2D0*CTH2-1D0)+4D0*(SQMW**2+SQMZ**2)/SH2*CTH2+
0376      &      2D0*(SQMW+SQMZ)/SH*BE2*CTH))
0377             ATWIM=0D0
0378             AUWRE=XW1/SQMZ*SH/(UH-SQMW)*((CTH+BE2)**2*(3D0/2D0-BE2/2D0*
0379      &      CTH-(SQMW+SQMZ)/SH+(SQMW-SQMZ)**2/(SH*SQMW))+4D0*
0380      &      ((SQMW+SQMZ)/SH*(1D0-3D0*CTH2)+8D0*SQMW*SQMZ/SH2*
0381      &      (2D0*CTH2-1D0)+4D0*(SQMW**2+SQMZ**2)/SH2*CTH2-
0382      &      2D0*(SQMW+SQMZ)/SH*BE2*CTH))
0383             AUWIM=0D0
0384             A4RE=2D0*XW1/SQMZ*(3D0-CTH2-4D0*(SQMW+SQMZ)/SH)
0385             A4IM=0D0
0386             FACWW=COMFAC*1D0/(4096D0*PARU(1)**2*16D0*XW1**2)*
0387      &      (AEM/XW)**4*(SH/SQMW)**2*(SQMZ/SQMW)*SH2
0388             IF(MSTP(46).LE.0) FACWW=FACWW*(ASHRE**2+ASHIM**2)
0389             IF(MSTP(46).EQ.1) FACWW=FACWW*((ASHRE+ATWRE+AUWRE+A4RE)**2+
0390      &      (ASHIM+ATWIM+AUWIM+A4IM)**2)
0391             IF(MSTP(46).EQ.2) FACWW=FACWW*((ATWRE+AUWRE+A4RE)**2+
0392      &      (ATWIM+AUWIM+A4IM)**2)
0393
0394           ELSE
0395 C...Strongly interacting Z_L/W_L model of Dobado, Herrero, Terron
0396             FACWW=COMFAC*(AEM/(16D0*PARU(1)*XW*XW1))**2*(64D0/9D0)*
0397      &      ABS(A00U-A20U)**2
0398           ENDIF
0399           FACWW=FACWW*WIDS(24,1)
0400
0401           DO 240 I=MMIN1,MMAX1
0402             IF(I.EQ.0.OR.KFAC(1,I).EQ.0) GOTO 240
0403             EI=KCHG(IABS(I),1)/3D0
0404             AI=SIGN(1D0,EI)
0405             VI=AI-4D0*EI*XWV
0406             AVI=AI**2+VI**2
0407             DO 230 J=MMIN2,MMAX2
0408               IF(J.EQ.0.OR.KFAC(2,J).EQ.0) GOTO 230
0409               EJ=KCHG(IABS(J),1)/3D0
0410               AJ=SIGN(1D0,EJ)
0411               VJ=AJ-4D0*EJ*XWV
0412               AVJ=AJ**2+VJ**2
0413               NCHN=NCHN+1
0414               ISIG(NCHN,1)=I
0415               ISIG(NCHN,2)=J
0416               ISIG(NCHN,3)=1
0417               SIGH(NCHN)=FACWW*AVI*AVJ
0418   230       CONTINUE
0419   240     CONTINUE
0420   250     CONTINUE
0421
0422         ELSEIF(ISUB.EQ.73) THEN
0423 C...Z0 + W+/- -> Z0 + W+/-
0424           IF(SH.LE.2D0*SQMZ+2D0*SQMW) GOTO 280
0425
0426           IF(MSTP(46).LE.2) THEN
0427 C...Exact scattering ME:s for on-mass-shell gauge bosons
0428             BE2=1D0-2D0*(SQMZ+SQMW)/SH+((SQMZ-SQMW)/SH)**2
0429             EP1=1D0-(SQMZ-SQMW)/SH
0430             EP2=1D0+(SQMZ-SQMW)/SH
0431             TH=-0.5D0*SH*BE2*(1D0-CTH)
0432             UH=(SQMZ-SQMW)**2/SH-0.5D0*SH*BE2*(1D0+CTH)
0433             IF(MAX(TH,UH).GT.-1D0) GOTO 280
0434             THANG=(BE2-EP1*CTH)*(BE2-EP2*CTH)
0435             ATHRE=(TH-SQMH)/((TH-SQMH)**2+GMMH**2)*THANG
0436             ATHIM=-GMMH/((TH-SQMH)**2+GMMH**2)*THANG
0437             ASWRE=-XW1/SQMZ*SH/(SH-SQMW)*(-BE2*(EP1+EP2)**4*CTH+
0438      &      1D0/4D0*(BE2+EP1*EP2)**2*((EP1-EP2)**2-4D0*BE2*CTH)+
0439      &      2D0*BE2*(BE2+EP1*EP2)*(EP1+EP2)**2*CTH-
0440      &      1D0/16D0*SH/SQMW*(EP1**2-EP2**2)**2*(BE2+EP1*EP2)**2)
0441             ASWIM=0D0
0442             AUWRE=XW1/SQMZ*SH/(UH-SQMW)*(-BE2*(EP2+EP1*CTH)*
0443      &      (EP1+EP2*CTH)*(BE2+EP1*EP2)+BE2*(EP2+EP1*CTH)*
0444      &      (BE2+EP1*EP2*CTH)*(2D0*EP2-EP2*CTH+EP1)-
0445      &      BE2*(EP2+EP1*CTH)**2*(BE2-EP2**2*CTH)-1D0/8D0*
0446      &      (BE2+EP1*EP2*CTH)**2*((EP1+EP2)**2+2D0*BE2*(1D0-CTH))+
0447      &      1D0/32D0*SH/SQMW*(BE2+EP1*EP2*CTH)**2*
0448      &      (EP1**2-EP2**2)**2-BE2*(EP1+EP2*CTH)*(EP2+EP1*CTH)*
0449      &      (BE2+EP1*EP2)+BE2*(EP1+EP2*CTH)*(BE2+EP1*EP2*CTH)*
0450      &      (2D0*EP1-EP1*CTH+EP2)-BE2*(EP1+EP2*CTH)**2*
0451      &      (BE2-EP1**2*CTH)-1D0/8D0*(BE2+EP1*EP2*CTH)**2*
0452      &      ((EP1+EP2)**2+2D0*BE2*(1D0-CTH))+1D0/32D0*SH/SQMW*
0453      &      (BE2+EP1*EP2*CTH)**2*(EP1**2-EP2**2)**2)
0454             AUWIM=0D0
0455             A4RE=XW1/SQMZ*(EP1**2*EP2**2*(CTH**2-1D0)-
0456      &      2D0*BE2*(EP1**2+EP2**2+EP1*EP2)*CTH-2D0*BE2*EP1*EP2)
0457             A4IM=0D0
0458             FACZW=COMFAC*1D0/(4096D0*PARU(1)**2*4D0*XW1)*(AEM/XW)**4*
0459      &      (SH/SQMW)**2*SQRT(SQMZ/SQMW)*SH2
0460             IF(MSTP(46).LE.0) FACZW=0D0
0461             IF(MSTP(46).EQ.1) FACZW=FACZW*((ATHRE+ASWRE+AUWRE+A4RE)**2+
0462      &      (ATHIM+ASWIM+AUWIM+A4IM)**2)
0463             IF(MSTP(46).EQ.2) FACZW=FACZW*((ASWRE+AUWRE+A4RE)**2+
0464      &      (ASWIM+AUWIM+A4IM)**2)
0465
0466           ELSE
0467 C...Strongly interacting Z_L/W_L model of Dobado, Herrero, Terron
0468             FACZW=COMFAC*AEM**2/(64D0*PARU(1)**2*XW**2*XW1)*16D0*
0469      &      ABS(A20U+3D0*A11U*DBLE(CTH))**2
0470           ENDIF
0471           FACZW=FACZW*WIDS(23,2)
0472
0473           DO 270 I=MMIN1,MMAX1
0474             IF(I.EQ.0.OR.KFAC(1,I).EQ.0) GOTO 270
0475             EI=KCHG(IABS(I),1)/3D0
0476             AI=SIGN(1D0,EI)
0477             VI=AI-4D0*EI*XWV
0478             AVI=AI**2+VI**2
0479             KCHWI=ISIGN(1,KCHG(IABS(I),1)*ISIGN(1,I))
0480             DO 260 J=MMIN2,MMAX2
0481               IF(J.EQ.0.OR.KFAC(2,J).EQ.0) GOTO 260
0482               EJ=KCHG(IABS(J),1)/3D0
0483               AJ=SIGN(1D0,EJ)
0484               VJ=AI-4D0*EJ*XWV
0485               AVJ=AJ**2+VJ**2
0486               KCHWJ=ISIGN(1,KCHG(IABS(J),1)*ISIGN(1,J))
0487               NCHN=NCHN+1
0488               ISIG(NCHN,1)=I
0489               ISIG(NCHN,2)=J
0490               ISIG(NCHN,3)=1
0491               SIGH(NCHN)=FACZW*AVI*VINT(180+J)*WIDS(24,(5-KCHWJ)/2)
0492               NCHN=NCHN+1
0493               ISIG(NCHN,1)=I
0494               ISIG(NCHN,2)=J
0495               ISIG(NCHN,3)=2
0496               SIGH(NCHN)=FACZW*VINT(180+I)*WIDS(24,(5-KCHWI)/2)*AVJ
0497   260       CONTINUE
0498   270     CONTINUE
0499   280     CONTINUE
0500
0501         ELSEIF(ISUB.EQ.75) THEN
0502 C...W+ + W- -> gamma + gamma
0503
0504         ELSEIF(ISUB.EQ.76) THEN
0505 C...W+ + W- -> Z0 + Z0
0506           IF(SH.LE.4.01D0*SQMZ) GOTO 310
0507
0508           IF(MSTP(46).LE.2) THEN
0509 C...Exact scattering ME:s for on-mass-shell gauge bosons
0510             BE2=SQRT((1D0-4D0*SQMW/SH)*(1D0-4D0*SQMZ/SH))
0511             CTH2=CTH**2
0512             TH=-0.5D0*SH*(1D0-2D0*(SQMW+SQMZ)/SH-BE2*CTH)
0513             UH=-0.5D0*SH*(1D0-2D0*(SQMW+SQMZ)/SH+BE2*CTH)
0514             IF(MAX(TH,UH).GT.-1D0) GOTO 310
0515             SHANG=4D0*SQRT(SQMW/(SQMZ*XW1))*(1D0-2D0*SQMW/SH)*
0516      &      (1D0-2D0*SQMZ/SH)
0517             ASHRE=(SH-SQMH)/((SH-SQMH)**2+GMMH**2)*SHANG
0518             ASHIM=-GMMH/((SH-SQMH)**2+GMMH**2)*SHANG
0519             ATWRE=XW1/SQMZ*SH/(TH-SQMW)*((CTH-BE2)**2*(3D0/2D0+BE2/2D0*
0520      &      CTH-(SQMW+SQMZ)/SH+(SQMW-SQMZ)**2/(SH*SQMW))+4D0*
0521      &      ((SQMW+SQMZ)/SH*(1D0-3D0*CTH2)+8D0*SQMW*SQMZ/SH2*
0522      &      (2D0*CTH2-1D0)+4D0*(SQMW**2+SQMZ**2)/SH2*CTH2+
0523      &      2D0*(SQMW+SQMZ)/SH*BE2*CTH))
0524             ATWIM=0D0
0525             AUWRE=XW1/SQMZ*SH/(UH-SQMW)*((CTH+BE2)**2*(3D0/2D0-BE2/2D0*
0526      &      CTH-(SQMW+SQMZ)/SH+(SQMW-SQMZ)**2/(SH*SQMW))+4D0*
0527      &      ((SQMW+SQMZ)/SH*(1D0-3D0*CTH2)+8D0*SQMW*SQMZ/SH2*
0528      &      (2D0*CTH2-1D0)+4D0*(SQMW**2+SQMZ**2)/SH2*CTH2-
0529      &      2D0*(SQMW+SQMZ)/SH*BE2*CTH))
0530             AUWIM=0D0
0531             A4RE=2D0*XW1/SQMZ*(3D0-CTH2-4D0*(SQMW+SQMZ)/SH)
0532             A4IM=0D0
0533             FACZZ=COMFAC*1D0/(4096D0*PARU(1)**2)*(AEM/XW)**4*
0534      &      (SH/SQMW)**2*SH2
0535             IF(MSTP(46).LE.0) FACZZ=FACZZ*(ASHRE**2+ASHIM**2)
0536             IF(MSTP(46).EQ.1) FACZZ=FACZZ*((ASHRE+ATWRE+AUWRE+A4RE)**2+
0537      &      (ASHIM+ATWIM+AUWIM+A4IM)**2)
0538             IF(MSTP(46).EQ.2) FACZZ=FACZZ*((ATWRE+AUWRE+A4RE)**2+
0539      &      (ATWIM+AUWIM+A4IM)**2)
0540
0541           ELSE
0542 C...Strongly interacting Z_L/W_L model of Dobado, Herrero, Terron
0543             FACZZ=COMFAC*(AEM/(4D0*PARU(1)*XW))**2*(64D0/9D0)*
0544      &      ABS(A00U-A20U)**2
0545           ENDIF
0546           FACZZ=FACZZ*WIDS(23,1)
0547
0548           DO 300 I=MMIN1,MMAX1
0549             IF(I.EQ.0.OR.KFAC(1,I).EQ.0) GOTO 300
0550             EI=SIGN(1D0,DBLE(I))*KCHG(IABS(I),1)
0551             DO 290 J=MMIN2,MMAX2
0552               IF(J.EQ.0.OR.KFAC(2,J).EQ.0) GOTO 290
0553               EJ=SIGN(1D0,DBLE(J))*KCHG(IABS(J),1)
0554               IF(EI*EJ.GT.0D0) GOTO 290
0555               NCHN=NCHN+1
0556               ISIG(NCHN,1)=I
0557               ISIG(NCHN,2)=J
0558               ISIG(NCHN,3)=1
0559               SIGH(NCHN)=0.5D0*FACZZ*VINT(180+I)*VINT(180+J)
0560   290       CONTINUE
0561   300     CONTINUE
0562   310     CONTINUE
0563
0564         ELSEIF(ISUB.EQ.77) THEN
0565 C...W+/- + W+/- -> W+/- + W+/-
0566           IF(SH.LE.4.01D0*SQMW) GOTO 340
0567
0568           IF(MSTP(46).LE.2) THEN
0569 C...Exact scattering ME:s for on-mass-shell gauge bosons
0570             BE2=1D0-4D0*SQMW/SH
0571             BE4=BE2**2
0572             CTH2=CTH**2
0573             CTH3=CTH**3
0574             TH=-0.5D0*SH*BE2*(1D0-CTH)
0575             UH=-0.5D0*SH*BE2*(1D0+CTH)
0576             IF(MAX(TH,UH).GT.-1D0) GOTO 340
0577             SHANG=(1D0+BE2)**2
0578             ASHRE=(SH-SQMH)/((SH-SQMH)**2+GMMH**2)*SHANG
0579             ASHIM=-GMMH/((SH-SQMH)**2+GMMH**2)*SHANG
0580             THANG=(BE2-CTH)**2
0581             ATHRE=(TH-SQMH)/((TH-SQMH)**2+GMMH**2)*THANG
0582             ATHIM=-GMMH/((TH-SQMH)**2+GMMH**2)*THANG
0583             UHANG=(BE2+CTH)**2
0584             AUHRE=(UH-SQMH)/((UH-SQMH)**2+GMMH**2)*UHANG
0585             AUHIM=-GMMH/((UH-SQMH)**2+GMMH**2)*UHANG
0586             SGZANG=1D0/SQMW*BE2*(3D0-BE2)**2*CTH
0587             ASGRE=XW*SGZANG
0588             ASGIM=0D0
0589             ASZRE=XW1*SH/(SH-SQMZ)*SGZANG
0590             ASZIM=0D0
0591             TGZANG=1D0/SQMW*(BE2*(4D0-2D0*BE2+BE4)+BE2*(4D0-10D0*BE2+
0592      &      BE4)*CTH+(2D0-11D0*BE2+10D0*BE4)*CTH2+BE2*CTH3)
0593             ATGRE=0.5D0*XW*SH/TH*TGZANG
0594             ATGIM=0D0
0595             ATZRE=0.5D0*XW1*SH/(TH-SQMZ)*TGZANG
0596             ATZIM=0D0
0597             UGZANG=1D0/SQMW*(BE2*(4D0-2D0*BE2+BE4)-BE2*(4D0-10D0*BE2+
0598      &      BE4)*CTH+(2D0-11D0*BE2+10D0*BE4)*CTH2-BE2*CTH3)
0599             AUGRE=0.5D0*XW*SH/UH*UGZANG
0600             AUGIM=0D0
0601             AUZRE=0.5D0*XW1*SH/(UH-SQMZ)*UGZANG
0602             AUZIM=0D0
0603             A4ARE=1D0/SQMW*(1D0+2D0*BE2-6D0*BE2*CTH-CTH2)
0604             A4AIM=0D0
0605             A4SRE=2D0/SQMW*(1D0+2D0*BE2-CTH2)
0606             A4SIM=0D0
0607             FWW=COMFAC*1D0/(4096D0*PARU(1)**2)*(AEM/XW)**4*
0608      &      (SH/SQMW)**2*SH2
0609             IF(MSTP(46).LE.0) THEN
0610               AWWARE=ASHRE
0611               AWWAIM=ASHIM
0612               AWWSRE=0D0
0613               AWWSIM=0D0
0614             ELSEIF(MSTP(46).EQ.1) THEN
0615               AWWARE=ASHRE+ATHRE+ASGRE+ASZRE+ATGRE+ATZRE+A4ARE
0616               AWWAIM=ASHIM+ATHIM+ASGIM+ASZIM+ATGIM+ATZIM+A4AIM
0617               AWWSRE=-ATHRE-AUHRE+ATGRE+ATZRE+AUGRE+AUZRE+A4SRE
0618               AWWSIM=-ATHIM-AUHIM+ATGIM+ATZIM+AUGIM+AUZIM+A4SIM
0619             ELSE
0620               AWWARE=ASGRE+ASZRE+ATGRE+ATZRE+A4ARE
0621               AWWAIM=ASGIM+ASZIM+ATGIM+ATZIM+A4AIM
0622               AWWSRE=ATGRE+ATZRE+AUGRE+AUZRE+A4SRE
0623               AWWSIM=ATGIM+ATZIM+AUGIM+AUZIM+A4SIM
0624             ENDIF
0625             AWWA2=AWWARE**2+AWWAIM**2
0626             AWWS2=AWWSRE**2+AWWSIM**2
0627
0628           ELSE
0629 C...Strongly interacting Z_L/W_L model of Dobado, Herrero, Terron
0630             FWWA=COMFAC*(AEM/(4D0*PARU(1)*XW))**2*(64D0/9D0)*
0631      &      ABS(A00U+0.5D0*A20U+4.5D0*A11U*DBLE(CTH))**2
0632             FWWS=COMFAC*(AEM/(4D0*PARU(1)*XW))**2*64D0*ABS(A20U)**2
0633           ENDIF
0634
0635           DO 330 I=MMIN1,MMAX1
0636             IF(I.EQ.0.OR.KFAC(1,I).EQ.0) GOTO 330
0637             EI=SIGN(1D0,DBLE(I))*KCHG(IABS(I),1)
0638             DO 320 J=MMIN2,MMAX2
0639               IF(J.EQ.0.OR.KFAC(2,J).EQ.0) GOTO 320
0640               EJ=SIGN(1D0,DBLE(J))*KCHG(IABS(J),1)
0641               IF(EI*EJ.LT.0D0) THEN
0642 C...W+W-
0643                 IF(MSTP(45).EQ.1) GOTO 320
0644                 IF(MSTP(46).LE.2) FACWW=FWW*AWWA2*WIDS(24,1)
0645                 IF(MSTP(46).GE.3) FACWW=FWWA*WIDS(24,1)
0646               ELSE
0647 C...W+W+/W-W-
0648                 IF(MSTP(45).EQ.2) GOTO 320
0649                 IF(MSTP(46).LE.2) FACWW=FWW*AWWS2
0650                 IF(MSTP(46).GE.3) FACWW=FWWS
0651                 IF(EI.GT.0D0) FACWW=FACWW*WIDS(24,4)
0652                 IF(EI.LT.0D0) FACWW=FACWW*WIDS(24,5)
0653               ENDIF
0654               NCHN=NCHN+1
0655               ISIG(NCHN,1)=I
0656               ISIG(NCHN,2)=J
0657               ISIG(NCHN,3)=1
0658               SIGH(NCHN)=FACWW*VINT(180+I)*VINT(180+J)
0659               IF(EI*EJ.GT.0D0) SIGH(NCHN)=0.5D0*SIGH(NCHN)
0660   320       CONTINUE
0661   330     CONTINUE
0662   340     CONTINUE
0663         ENDIF
0664
0665       ELSEIF(ISUB.LE.120) THEN
0666         IF(ISUB.EQ.102) THEN
0667 C...g + g -> h0 (or H0, or A0)
0668           CALL PYWIDT(KFHIGG,SH,WDTP,WDTE)
0669           WDTP13=0D0
0670           DO 345 IDC=MDCY(KFHIGG,2),MDCY(KFHIGG,2)+MDCY(KFHIGG,3)-1
0671             IF(KFDP(IDC,1).EQ.21.AND.KFDP(IDC,2).EQ.21.AND.
0672      &      KFDP(IDC,3).EQ.0) WDTP13=PMAS(KFHIGG,2)*BRAT(IDC)
0673   345     CONTINUE
0674           IF(WDTP13.EQ.0D0) CALL PYERRM(26,
0675      &    '(PYSGHG:) did not find Higgs -> g g channel')
0676           HS=SHR*WDTP(0)
0677           HF=SHR*(WDTE(0,1)+WDTE(0,2)+WDTE(0,4))
0678           FACBW=4D0*COMFAC/((SH-SQMH)**2+HS**2)
0679           IF(ABS(SHR-PMAS(KFHIGG,1)).GT.PARP(48)*PMAS(KFHIGG,2))
0680      &    FACBW=0D0
0681           HI=SHR*WDTP13/32D0
0682           IF(KFAC(1,21)*KFAC(2,21).EQ.0) GOTO 350
0683           NCHN=NCHN+1
0684           ISIG(NCHN,1)=21
0685           ISIG(NCHN,2)=21
0686           ISIG(NCHN,3)=1
0687           SIGH(NCHN)=HI*FACBW*HF
0688   350     CONTINUE
0689
0690         ELSEIF(ISUB.EQ.103) THEN
0691 C...gamma + gamma -> h0 (or H0, or A0)
0692           CALL PYWIDT(KFHIGG,SH,WDTP,WDTE)
0693           WDTP14=0D0
0694           DO 355 IDC=MDCY(KFHIGG,2),MDCY(KFHIGG,2)+MDCY(KFHIGG,3)-1
0695             IF(KFDP(IDC,1).EQ.22.AND.KFDP(IDC,2).EQ.22.AND.
0696      &      KFDP(IDC,3).EQ.0) WDTP14=PMAS(KFHIGG,2)*BRAT(IDC)
0697   355     CONTINUE
0698           IF(WDTP14.EQ.0D0) CALL PYERRM(26,
0699      &    '(PYSGHG:) did not find Higgs -> gamma gamma channel')
0700           HS=SHR*WDTP(0)
0701           HF=SHR*(WDTE(0,1)+WDTE(0,2)+WDTE(0,4))
0702           FACBW=4D0*COMFAC/((SH-SQMH)**2+HS**2)
0703           IF(ABS(SHR-PMAS(KFHIGG,1)).GT.PARP(48)*PMAS(KFHIGG,2))
0704      &    FACBW=0D0
0705           HI=SHR*WDTP14*2D0
0706           IF(KFAC(1,22)*KFAC(2,22).EQ.0) GOTO 360
0707           NCHN=NCHN+1
0708           ISIG(NCHN,1)=22
0709           ISIG(NCHN,2)=22
0710           ISIG(NCHN,3)=1
0711           SIGH(NCHN)=HI*FACBW*HF
0712   360     CONTINUE
0713
0714         ELSEIF(ISUB.EQ.110) THEN
0715 C...f + fbar -> gamma + h0
0716           THUH=MAX(TH*UH,SH*CKIN(3)**2)
0717           FACHG=COMFAC*(3D0*AEM**4)/(2D0*PARU(1)**2*XW*SQMW)*SH*THUH
0718           FACHG=FACHG*WIDS(KFHIGG,2)
0719 C...Calculate loop contributions for intermediate gamma* and Z0
0720           CIGTOT=DCMPLX(0D0,0D0)
0721           CIZTOT=DCMPLX(0D0,0D0)
0722           JMAX=3*MSTP(1)+1
0723           DO 370 J=1,JMAX
0724             IF(J.LE.2*MSTP(1)) THEN
0725               FNC=1D0
0726               EJ=KCHG(J,1)/3D0
0727               AJ=SIGN(1D0,EJ+0.1D0)
0728               VJ=AJ-4D0*EJ*XWV
0729               BALP=SQM4/(2D0*PMAS(J,1))**2
0730               BBET=SH/(2D0*PMAS(J,1))**2
0731             ELSEIF(J.LE.3*MSTP(1)) THEN
0732               FNC=3D0
0733               JL=2*(J-2*MSTP(1))-1
0734               EJ=KCHG(10+JL,1)/3D0
0735               AJ=SIGN(1D0,EJ+0.1D0)
0736               VJ=AJ-4D0*EJ*XWV
0737               BALP=SQM4/(2D0*PMAS(10+JL,1))**2
0738               BBET=SH/(2D0*PMAS(10+JL,1))**2
0739             ELSE
0740               BALP=SQM4/(2D0*PMAS(24,1))**2
0741               BBET=SH/(2D0*PMAS(24,1))**2
0742             ENDIF
0743             BABI=1D0/(BALP-BBET)
0744             IF(BALP.LT.1D0) THEN
0745               F0ALP=DCMPLX(DBLE(ASIN(SQRT(BALP))),0D0)
0746               F1ALP=F0ALP**2
0747             ELSE
0748               F0ALP=DCMPLX(DBLE(LOG(SQRT(BALP)+SQRT(BALP-1D0))),
0749      &        -DBLE(0.5D0*PARU(1)))
0750               F1ALP=-F0ALP**2
0751             ENDIF
0752             F2ALP=DBLE(SQRT(ABS(BALP-1D0)/BALP))*F0ALP
0753             IF(BBET.LT.1D0) THEN
0754               F0BET=DCMPLX(DBLE(ASIN(SQRT(BBET))),0D0)
0755               F1BET=F0BET**2
0756             ELSE
0757               F0BET=DCMPLX(DBLE(LOG(SQRT(BBET)+SQRT(BBET-1D0))),
0758      &        -DBLE(0.5D0*PARU(1)))
0759               F1BET=-F0BET**2
0760             ENDIF
0761             F2BET=DBLE(SQRT(ABS(BBET-1D0)/BBET))*F0BET
0762             IF(J.LE.3*MSTP(1)) THEN
0763               FIF=DBLE(0.5D0*BABI)+DBLE(BABI**2)*(DBLE(0.5D0*(1D0-BALP+
0764      &        BBET))*(F1BET-F1ALP)+DBLE(BBET)*(F2BET-F2ALP))
0765               CIGTOT=CIGTOT+DBLE(FNC*EJ**2)*FIF
0766               CIZTOT=CIZTOT+DBLE(FNC*EJ*VJ)*FIF
0767             ELSE
0768               TXW=XW/XW1
0769               CIGTOT=CIGTOT-0.5*(DBLE(BABI*(1.5D0+BALP))+DBLE(BABI**2)*
0770      &        (DBLE(1.5D0-3D0*BALP+4D0*BBET)*(F1BET-F1ALP)+
0771      &        DBLE(BBET*(2D0*BALP+3D0))*(F2BET-F2ALP)))
0772               CIZTOT=CIZTOT-DBLE(0.5D0*BABI*XW1)*(DBLE(5D0-TXW+2D0*BALP*
0773      &        (1D0-TXW))*(1D0+DBLE(2D0*BABI*BBET)*(F2BET-F2ALP))+
0774      &        DBLE(BABI*(4D0*BBET*(3D0-TXW)-(2D0*BALP-1D0)*(5D0-TXW)))*
0775      &        (F1BET-F1ALP))
0776             ENDIF
0777   370     CONTINUE
0778           CIGTOT=CIGTOT/DBLE(SH)
0779           CIZTOT=CIZTOT*DBLE(XWC)/DCMPLX(DBLE(SH-SQMZ),DBLE(GMMZ))
0780 C...Loop over initial flavours
0781           DO 380 I=MMINA,MMAXA
0782             IF(I.EQ.0.OR.KFAC(1,I)*KFAC(2,-I).EQ.0) GOTO 380
0783             EI=KCHG(IABS(I),1)/3D0
0784             AI=SIGN(1D0,EI)
0785             VI=AI-4D0*EI*XWV
0786             FCOI=1D0
0787             IF(IABS(I).LE.10) FCOI=FACA/3D0
0788             NCHN=NCHN+1
0789             ISIG(NCHN,1)=I
0790             ISIG(NCHN,2)=-I
0791             ISIG(NCHN,3)=1
0792             SIGH(NCHN)=FACHG*FCOI*(ABS(DBLE(EI)*CIGTOT+DBLE(VI)*
0793      &      CIZTOT)**2+AI**2*ABS(CIZTOT)**2)
0794   380     CONTINUE
0795
0796         ELSEIF(ISUB.EQ.111) THEN
0797 C...f + fbar -> g + h0 (q + qbar -> g + h0 only)
0798           IF(MSTP(38).NE.0) THEN
0799 C...Simple case: only do gg <-> h exactly.
0800           CALL PYWIDT(KFHIGG,SQM4,WDTP,WDTE)
0801           WDTP13=0D0
0802           DO 385 IDC=MDCY(KFHIGG,2),MDCY(KFHIGG,2)+MDCY(KFHIGG,3)-1
0803             IF(KFDP(IDC,1).EQ.21.AND.KFDP(IDC,2).EQ.21.AND.
0804      &      KFDP(IDC,3).EQ.0) WDTP13=PMAS(KFHIGG,2)*BRAT(IDC)
0805   385     CONTINUE
0806           IF(WDTP13.EQ.0D0) CALL PYERRM(26,
0807      &    '(PYSGHG:) did not find Higgs -> g g channel')
0808           FACGH=COMFAC*FACA*(2D0/9D0)*AS*(WDTP13/SQRT(SQM4))*
0809      &    (TH**2+UH**2)/(SH*SQM4)
0810 C...Propagators: as simulated in PYOFSH and as desired
0811           HBW4=GMMH/((SQM4-SQMH)**2+GMMH**2)
0812           GMMHC=SQRT(SQM4)*WDTP(0)
0813           HBW4C=SQRT(SQM4)*(WDTE(0,1)+WDTE(0,2)+WDTE(0,4))/
0814      &    ((SQM4-SQMH)**2+GMMHC**2)
0815           FACGH=FACGH*HBW4C/HBW4
0816           ELSE
0817 C...Messy case: do full loop integrals
0818           A5STUR=0D0
0819           A5STUI=0D0
0820           DO 390 I=1,2*MSTP(1)
0821             SQMQ=PMAS(I,1)**2
0822             EPSS=4D0*SQMQ/SH
0823             EPSH=4D0*SQMQ/SQMH
0824             CALL PYWAUX(1,EPSS,W1SR,W1SI)
0825             CALL PYWAUX(1,EPSH,W1HR,W1HI)
0826             CALL PYWAUX(2,EPSS,W2SR,W2SI)
0827             CALL PYWAUX(2,EPSH,W2HR,W2HI)
0828             A5STUR=A5STUR+EPSH*(1D0+SH/(TH+UH)*(W1SR-W1HR)+
0829      &      (0.25D0-SQMQ/(TH+UH))*(W2SR-W2HR))
0830             A5STUI=A5STUI+EPSH*(SH/(TH+UH)*(W1SI-W1HI)+
0831      &      (0.25D0-SQMQ/(TH+UH))*(W2SI-W2HI))
0832   390     CONTINUE
0833           FACGH=COMFAC*FACA/(144D0*PARU(1)**2)*AEM/XW*AS**3*SQMH/SQMW*
0834      &    SQMH/SH*(UH**2+TH**2)/(UH+TH)**2*(A5STUR**2+A5STUI**2)
0835           FACGH=FACGH*WIDS(25,2)
0836           ENDIF
0837           DO 400 I=MMINA,MMAXA
0838             IF(I.EQ.0.OR.IABS(I).GT.MSTP(58).OR.
0839      &      KFAC(1,I)*KFAC(2,-I).EQ.0) GOTO 400
0840             NCHN=NCHN+1
0841             ISIG(NCHN,1)=I
0842             ISIG(NCHN,2)=-I
0843             ISIG(NCHN,3)=1
0844             SIGH(NCHN)=FACGH
0845   400     CONTINUE
0846
0847         ELSEIF(ISUB.EQ.112) THEN
0848 C...f + g -> f + h0 (q + g -> q + h0 only)
0849           IF(MSTP(38).NE.0) THEN
0850 C...Simple case: only do gg <-> h exactly.
0851           CALL PYWIDT(KFHIGG,SQM4,WDTP,WDTE)
0852           WDTP13=0D0
0853           DO 405 IDC=MDCY(KFHIGG,2),MDCY(KFHIGG,2)+MDCY(KFHIGG,3)-1
0854             IF(KFDP(IDC,1).EQ.21.AND.KFDP(IDC,2).EQ.21.AND.
0855      &      KFDP(IDC,3).EQ.0) WDTP13=PMAS(KFHIGG,2)*BRAT(IDC)
0856   405     CONTINUE
0857           IF(WDTP13.EQ.0D0) CALL PYERRM(26,
0858      &    '(PYSGHG:) did not find Higgs -> g g channel')
0859           FACQH=COMFAC*FACA*(1D0/12D0)*AS*(WDTP13/SQRT(SQM4))*
0860      &    (SH**2+UH**2)/(-TH*SQM4)
0861 C...Propagators: as simulated in PYOFSH and as desired
0862           HBW4=GMMH/((SQM4-SQMH)**2+GMMH**2)
0863           GMMHC=SQRT(SQM4)*WDTP(0)
0864           HBW4C=SQRT(SQM4)*(WDTE(0,1)+WDTE(0,2)+WDTE(0,4))/
0865      &    ((SQM4-SQMH)**2+GMMHC**2)
0866           FACQH=FACQH*HBW4C/HBW4
0867           ELSE
0868 C...Messy case: do full loop integrals
0869           A5TSUR=0D0
0870           A5TSUI=0D0
0871           DO 410 I=1,2*MSTP(1)
0872             SQMQ=PMAS(I,1)**2
0873             EPST=4D0*SQMQ/TH
0874             EPSH=4D0*SQMQ/SQMH
0875             CALL PYWAUX(1,EPST,W1TR,W1TI)
0876             CALL PYWAUX(1,EPSH,W1HR,W1HI)
0877             CALL PYWAUX(2,EPST,W2TR,W2TI)
0878             CALL PYWAUX(2,EPSH,W2HR,W2HI)
0879             A5TSUR=A5TSUR+EPSH*(1D0+TH/(SH+UH)*(W1TR-W1HR)+
0880      &      (0.25D0-SQMQ/(SH+UH))*(W2TR-W2HR))
0881             A5TSUI=A5TSUI+EPSH*(TH/(SH+UH)*(W1TI-W1HI)+
0882      &      (0.25D0-SQMQ/(SH+UH))*(W2TI-W2HI))
0883   410     CONTINUE
0884           FACQH=COMFAC*FACA/(384D0*PARU(1)**2)*AEM/XW*AS**3*SQMH/SQMW*
0885      &    SQMH/(-TH)*(UH**2+SH**2)/(UH+SH)**2*(A5TSUR**2+A5TSUI**2)
0886           FACQH=FACQH*WIDS(25,2)
0887           ENDIF
0888           DO 430 I=MMINA,MMAXA
0889             IF(I.EQ.0.OR.IABS(I).GT.MSTP(58)) GOTO 430
0890             DO 420 ISDE=1,2
0891               IF(ISDE.EQ.1.AND.KFAC(1,I)*KFAC(2,21).EQ.0) GOTO 420
0892               IF(ISDE.EQ.2.AND.KFAC(1,21)*KFAC(2,I).EQ.0) GOTO 420
0893               NCHN=NCHN+1
0894               ISIG(NCHN,ISDE)=I
0895               ISIG(NCHN,3-ISDE)=21
0896               ISIG(NCHN,3)=1
0897               SIGH(NCHN)=FACQH
0898   420       CONTINUE
0899   430     CONTINUE
0900
0901         ELSEIF(ISUB.EQ.113) THEN
0902 C...g + g -> g + h0
0903           IF(MSTP(38).NE.0) THEN
0904 C...Simple case: only do gg <-> h exactly.
0905           CALL PYWIDT(KFHIGG,SQM4,WDTP,WDTE)
0906           WDTP13=0D0
0907           DO 435 IDC=MDCY(KFHIGG,2),MDCY(KFHIGG,2)+MDCY(KFHIGG,3)-1
0908             IF(KFDP(IDC,1).EQ.21.AND.KFDP(IDC,2).EQ.21.AND.
0909      &      KFDP(IDC,3).EQ.0) WDTP13=PMAS(KFHIGG,2)*BRAT(IDC)
0910   435     CONTINUE
0911           IF(WDTP13.EQ.0D0) CALL PYERRM(26,
0912      &    '(PYSGHG:) did not find Higgs -> g g channel')
0913           FACGH=COMFAC*FACA*(3D0/16D0)*AS*(WDTP13/SQRT(SQM4))*
0914      &    (SH**4+TH**4+UH**4+SQM4**4)/(SH*TH*UH*SQM4)
0915 C...Propagators: as simulated in PYOFSH and as desired
0916           HBW4=GMMH/((SQM4-SQMH)**2+GMMH**2)
0917           GMMHC=SQRT(SQM4)*WDTP(0)
0918           HBW4C=SQRT(SQM4)*(WDTE(0,1)+WDTE(0,2)+WDTE(0,4))/
0919      &    ((SQM4-SQMH)**2+GMMHC**2)
0920           FACGH=FACGH*HBW4C/HBW4
0921           ELSE
0922 C...Messy case: do full loop integrals
0923           A2STUR=0D0
0924           A2STUI=0D0
0925           A2USTR=0D0
0926           A2USTI=0D0
0927           A2TUSR=0D0
0928           A2TUSI=0D0
0929           A4STUR=0D0
0930           A4STUI=0D0
0931           DO 440 I=1,2*MSTP(1)
0932             SQMQ=PMAS(I,1)**2
0933             EPSS=4D0*SQMQ/SH
0934             EPST=4D0*SQMQ/TH
0935             EPSU=4D0*SQMQ/UH
0936             EPSH=4D0*SQMQ/SQMH
0937             IF(EPSH.LT.1D-6) GOTO 440
0938             CALL PYWAUX(1,EPSS,W1SR,W1SI)
0939             CALL PYWAUX(1,EPST,W1TR,W1TI)
0940             CALL PYWAUX(1,EPSU,W1UR,W1UI)
0941             CALL PYWAUX(1,EPSH,W1HR,W1HI)
0942             CALL PYWAUX(2,EPSS,W2SR,W2SI)
0943             CALL PYWAUX(2,EPST,W2TR,W2TI)
0944             CALL PYWAUX(2,EPSU,W2UR,W2UI)
0945             CALL PYWAUX(2,EPSH,W2HR,W2HI)
0946             CALL PYI3AU(EPSS,TH/UH,Y3STUR,Y3STUI)
0947             CALL PYI3AU(EPSS,UH/TH,Y3SUTR,Y3SUTI)
0948             CALL PYI3AU(EPST,SH/UH,Y3TSUR,Y3TSUI)
0949             CALL PYI3AU(EPST,UH/SH,Y3TUSR,Y3TUSI)
0950             CALL PYI3AU(EPSU,SH/TH,Y3USTR,Y3USTI)
0951             CALL PYI3AU(EPSU,TH/SH,Y3UTSR,Y3UTSI)
0952             CALL PYI3AU(EPSH,SQMH/SH*TH/UH,YHSTUR,YHSTUI)
0953             CALL PYI3AU(EPSH,SQMH/SH*UH/TH,YHSUTR,YHSUTI)
0954             CALL PYI3AU(EPSH,SQMH/TH*SH/UH,YHTSUR,YHTSUI)
0955             CALL PYI3AU(EPSH,SQMH/TH*UH/SH,YHTUSR,YHTUSI)
0956             CALL PYI3AU(EPSH,SQMH/UH*SH/TH,YHUSTR,YHUSTI)
0957             CALL PYI3AU(EPSH,SQMH/UH*TH/SH,YHUTSR,YHUTSI)
0958             W3STUR=YHSTUR-Y3STUR-Y3UTSR
0959             W3STUI=YHSTUI-Y3STUI-Y3UTSI
0960             W3SUTR=YHSUTR-Y3SUTR-Y3TUSR
0961             W3SUTI=YHSUTI-Y3SUTI-Y3TUSI
0962             W3TSUR=YHTSUR-Y3TSUR-Y3USTR
0963             W3TSUI=YHTSUI-Y3TSUI-Y3USTI
0964             W3TUSR=YHTUSR-Y3TUSR-Y3SUTR
0965             W3TUSI=YHTUSI-Y3TUSI-Y3SUTI
0966             W3USTR=YHUSTR-Y3USTR-Y3TSUR
0967             W3USTI=YHUSTI-Y3USTI-Y3TSUI
0968             W3UTSR=YHUTSR-Y3UTSR-Y3STUR
0969             W3UTSI=YHUTSI-Y3UTSI-Y3STUI
0970             B2STUR=SQMQ/SQMH**2*(SH*(UH-SH)/(SH+UH)+2D0*TH*UH*
0971      &      (UH+2D0*SH)/(SH+UH)**2*(W1TR-W1HR)+(SQMQ-SH/4D0)*
0972      &      (0.5D0*W2SR+0.5D0*W2HR-W2TR+W3STUR)+SH2*(2D0*SQMQ/
0973      &      (SH+UH)**2-0.5D0/(SH+UH))*(W2TR-W2HR)+0.5D0*TH*UH/SH*
0974      &      (W2HR-2D0*W2TR)+0.125D0*(SH-12D0*SQMQ-4D0*TH*UH/SH)*W3TSUR)
0975             B2STUI=SQMQ/SQMH**2*(2D0*TH*UH*(UH+2D0*SH)/(SH+UH)**2*
0976      &      (W1TI-W1HI)+(SQMQ-SH/4D0)*(0.5D0*W2SI+0.5D0*W2HI-W2TI+
0977      &      W3STUI)+SH2*(2D0*SQMQ/(SH+UH)**2-0.5D0/(SH+UH))*
0978      &      (W2TI-W2HI)+0.5D0*TH*UH/SH*(W2HI-2D0*W2TI)+0.125D0*
0979      &      (SH-12D0*SQMQ-4D0*TH*UH/SH)*W3TSUI)
0980             B2SUTR=SQMQ/SQMH**2*(SH*(TH-SH)/(SH+TH)+2D0*UH*TH*
0981      &      (TH+2D0*SH)/(SH+TH)**2*(W1UR-W1HR)+(SQMQ-SH/4D0)*
0982      &      (0.5D0*W2SR+0.5D0*W2HR-W2UR+W3SUTR)+SH2*(2D0*SQMQ/
0983      &      (SH+TH)**2-0.5D0/(SH+TH))*(W2UR-W2HR)+0.5D0*UH*TH/SH*
0984      &      (W2HR-2D0*W2UR)+0.125D0*(SH-12D0*SQMQ-4D0*UH*TH/SH)*W3USTR)
0985             B2SUTI=SQMQ/SQMH**2*(2D0*UH*TH*(TH+2D0*SH)/(SH+TH)**2*
0986      &      (W1UI-W1HI)+(SQMQ-SH/4D0)*(0.5D0*W2SI+0.5D0*W2HI-W2UI+
0987      &      W3SUTI)+SH2*(2D0*SQMQ/(SH+TH)**2-0.5D0/(SH+TH))*
0988      &      (W2UI-W2HI)+0.5D0*UH*TH/SH*(W2HI-2D0*W2UI)+0.125D0*
0989      &      (SH-12D0*SQMQ-4D0*UH*TH/SH)*W3USTI)
0990             B2TSUR=SQMQ/SQMH**2*(TH*(UH-TH)/(TH+UH)+2D0*SH*UH*
0991      &      (UH+2D0*TH)/(TH+UH)**2*(W1SR-W1HR)+(SQMQ-TH/4D0)*
0992      &      (0.5D0*W2TR+0.5D0*W2HR-W2SR+W3TSUR)+TH2*(2D0*SQMQ/
0993      &      (TH+UH)**2-0.5D0/(TH+UH))*(W2SR-W2HR)+0.5D0*SH*UH/TH*
0994      &      (W2HR-2D0*W2SR)+0.125D0*(TH-12D0*SQMQ-4D0*SH*UH/TH)*W3STUR)
0995             B2TSUI=SQMQ/SQMH**2*(2D0*SH*UH*(UH+2D0*TH)/(TH+UH)**2*
0996      &      (W1SI-W1HI)+(SQMQ-TH/4D0)*(0.5D0*W2TI+0.5D0*W2HI-W2SI+
0997      &      W3TSUI)+TH2*(2D0*SQMQ/(TH+UH)**2-0.5D0/(TH+UH))*
0998      &      (W2SI-W2HI)+0.5D0*SH*UH/TH*(W2HI-2D0*W2SI)+0.125D0*
0999      &      (TH-12D0*SQMQ-4D0*SH*UH/TH)*W3STUI)
1000             B2TUSR=SQMQ/SQMH**2*(TH*(SH-TH)/(TH+SH)+2D0*UH*SH*
1001      &      (SH+2D0*TH)/(TH+SH)**2*(W1UR-W1HR)+(SQMQ-TH/4D0)*
1002      &      (0.5D0*W2TR+0.5D0*W2HR-W2UR+W3TUSR)+TH2*(2D0*SQMQ/
1003      &      (TH+SH)**2-0.5D0/(TH+SH))*(W2UR-W2HR)+0.5D0*UH*SH/TH*
1004      &      (W2HR-2D0*W2UR)+0.125D0*(TH-12D0*SQMQ-4D0*UH*SH/TH)*W3UTSR)
1005             B2TUSI=SQMQ/SQMH**2*(2D0*UH*SH*(SH+2D0*TH)/(TH+SH)**2*
1006      &      (W1UI-W1HI)+(SQMQ-TH/4D0)*(0.5D0*W2TI+0.5D0*W2HI-W2UI+
1007      &      W3TUSI)+TH2*(2D0*SQMQ/(TH+SH)**2-0.5D0/(TH+SH))*
1008      &      (W2UI-W2HI)+0.5D0*UH*SH/TH*(W2HI-2D0*W2UI)+0.125D0*
1009      &      (TH-12D0*SQMQ-4D0*UH*SH/TH)*W3UTSI)
1010             B2USTR=SQMQ/SQMH**2*(UH*(TH-UH)/(UH+TH)+2D0*SH*TH*
1011      &      (TH+2D0*UH)/(UH+TH)**2*(W1SR-W1HR)+(SQMQ-UH/4D0)*
1012      &      (0.5D0*W2UR+0.5D0*W2HR-W2SR+W3USTR)+UH2*(2D0*SQMQ/
1013      &      (UH+TH)**2-0.5D0/(UH+TH))*(W2SR-W2HR)+0.5D0*SH*TH/UH*
1014      &      (W2HR-2D0*W2SR)+0.125D0*(UH-12D0*SQMQ-4D0*SH*TH/UH)*W3SUTR)
1015             B2USTI=SQMQ/SQMH**2*(2D0*SH*TH*(TH+2D0*UH)/(UH+TH)**2*
1016      &      (W1SI-W1HI)+(SQMQ-UH/4D0)*(0.5D0*W2UI+0.5D0*W2HI-W2SI+
1017      &      W3USTI)+UH2*(2D0*SQMQ/(UH+TH)**2-0.5D0/(UH+TH))*
1018      &      (W2SI-W2HI)+0.5D0*SH*TH/UH*(W2HI-2D0*W2SI)+0.125D0*
1019      &      (UH-12D0*SQMQ-4D0*SH*TH/UH)*W3SUTI)
1020             B2UTSR=SQMQ/SQMH**2*(UH*(SH-UH)/(UH+SH)+2D0*TH*SH*
1021      &      (SH+2D0*UH)/(UH+SH)**2*(W1TR-W1HR)+(SQMQ-UH/4D0)*
1022      &      (0.5D0*W2UR+0.5D0*W2HR-W2TR+W3UTSR)+UH2*(2D0*SQMQ/
1023      &      (UH+SH)**2-0.5D0/(UH+SH))*(W2TR-W2HR)+0.5D0*TH*SH/UH*
1024      &      (W2HR-2D0*W2TR)+0.125D0*(UH-12D0*SQMQ-4D0*TH*SH/UH)*W3TUSR)
1025             B2UTSI=SQMQ/SQMH**2*(2D0*TH*SH*(SH+2D0*UH)/(UH+SH)**2*
1026      &      (W1TI-W1HI)+(SQMQ-UH/4D0)*(0.5D0*W2UI+0.5D0*W2HI-W2TI+
1027      &      W3UTSI)+UH2*(2D0*SQMQ/(UH+SH)**2-0.5D0/(UH+SH))*
1028      &      (W2TI-W2HI)+0.5D0*TH*SH/UH*(W2HI-2D0*W2TI)+0.125D0*
1029      &      (UH-12D0*SQMQ-4D0*TH*SH/UH)*W3TUSI)
1030             B4STUR=0.25D0*EPSH*(-2D0/3D0+0.25D0*(EPSH-1D0)*
1031      &      (W2SR-W2HR+W3STUR))
1032             B4STUI=0.25D0*EPSH*0.25D0*(EPSH-1D0)*(W2SI-W2HI+W3STUI)
1033             B4TUSR=0.25D0*EPSH*(-2D0/3D0+0.25D0*(EPSH-1D0)*
1034      &      (W2TR-W2HR+W3TUSR))
1035             B4TUSI=0.25D0*EPSH*0.25D0*(EPSH-1D0)*(W2TI-W2HI+W3TUSI)
1036             B4USTR=0.25D0*EPSH*(-2D0/3D0+0.25D0*(EPSH-1D0)*
1037      &      (W2UR-W2HR+W3USTR))
1038             B4USTI=0.25D0*EPSH*0.25D0*(EPSH-1D0)*(W2UI-W2HI+W3USTI)
1039             A2STUR=A2STUR+B2STUR+B2SUTR
1040             A2STUI=A2STUI+B2STUI+B2SUTI
1041             A2USTR=A2USTR+B2USTR+B2UTSR
1042             A2USTI=A2USTI+B2USTI+B2UTSI
1043             A2TUSR=A2TUSR+B2TUSR+B2TSUR
1044             A2TUSI=A2TUSI+B2TUSI+B2TSUI
1045             A4STUR=A4STUR+B4STUR+B4USTR+B4TUSR
1046             A4STUI=A4STUI+B4STUI+B4USTI+B4TUSI
1047   440     CONTINUE
1048           FACGH=COMFAC*FACA*3D0/(128D0*PARU(1)**2)*AEM/XW*AS**3*
1049      &    SQMH/SQMW*SQMH**3/(SH*TH*UH)*(A2STUR**2+A2STUI**2+A2USTR**2+
1050      &    A2USTI**2+A2TUSR**2+A2TUSI**2+A4STUR**2+A4STUI**2)
1051           FACGH=FACGH*WIDS(25,2)
1052           ENDIF
1053           IF(KFAC(1,21)*KFAC(2,21).EQ.0) GOTO 450
1054           NCHN=NCHN+1
1055           ISIG(NCHN,1)=21
1056           ISIG(NCHN,2)=21
1057           ISIG(NCHN,3)=1
1058           SIGH(NCHN)=FACGH
1059   450     CONTINUE
1060         ENDIF
1061
1062       ELSEIF(ISUB.LE.170) THEN
1063         IF(ISUB.EQ.121) THEN
1064 C...g + g -> Q + Qbar + h0
1065           IF(KFAC(1,21)*KFAC(2,21).EQ.0) GOTO 460
1066           IA=KFPR(ISUBSV,2)
1067           PMF=PYMRUN(IA,SH)
1068           FACQQH=COMFAC*(4D0*PARU(1)*AEM/XW)*(4D0*PARU(1)*AS)**2*
1069      &    (0.5D0*PMF/PMAS(24,1))**2
1070           WID2=1D0
1071           IF(IA.EQ.6.OR.IA.EQ.7.OR.IA.EQ.8) WID2=WIDS(IA,1)
1072           FACQQH=FACQQH*WID2
1073           IF(MSTP(4).GE.1.OR.IHIGG.GE.2) THEN
1074             IKFI=1
1075             IF(IA.LE.10.AND.MOD(IA,2).EQ.0) IKFI=2
1076             IF(IA.GT.10) IKFI=3
1077             FACQQH=FACQQH*PARU(150+10*IHIGG+IKFI)**2
1078             IF(IMSS(1).NE.0.AND.IA.EQ.5) THEN
1079               FACQQH=FACQQH/(1D0+RMSS(41))**2
1080               IF(IHIGG.NE.3) THEN
1081                 FACQQH=FACQQH*(1D0+RMSS(41)*PARU(152+10*IHIGG)/
1082      &          PARU(151+10*IHIGG))**2
1083               ENDIF
1084             ENDIF
1085           ENDIF
1086           CALL PYQQBH(WTQQBH)
1087           CALL PYWIDT(KFHIGG,SH,WDTP,WDTE)
1088           HS=SHR*WDTP(0)
1089           HF=SHR*(WDTE(0,1)+WDTE(0,2)+WDTE(0,4))
1090           FACBW=(1D0/PARU(1))*VINT(2)*HF/((SH-SQMH)**2+HS**2)
1091           IF(ABS(SHR-PMAS(KFHIGG,1)).GT.PARP(48)*PMAS(KFHIGG,2))
1092      &    FACBW=0D0
1093           NCHN=NCHN+1
1094           ISIG(NCHN,1)=21
1095           ISIG(NCHN,2)=21
1096           ISIG(NCHN,3)=1
1097           SIGH(NCHN)=FACQQH*WTQQBH*FACBW
1098   460     CONTINUE
1099
1100         ELSEIF(ISUB.EQ.122) THEN
1101 C...q + qbar -> Q + Qbar + h0
1102           IA=KFPR(ISUBSV,2)
1103           PMF=PYMRUN(IA,SH)
1104           FACQQH=COMFAC*(4D0*PARU(1)*AEM/XW)*(4D0*PARU(1)*AS)**2*
1105      &    (0.5D0*PMF/PMAS(24,1))**2
1106           WID2=1D0
1107           IF(IA.EQ.6.OR.IA.EQ.7.OR.IA.EQ.8) WID2=WIDS(IA,1)
1108           FACQQH=FACQQH*WID2
1109           IF(MSTP(4).GE.1.OR.IHIGG.GE.2) THEN
1110             IKFI=1
1111             IF(IA.LE.10.AND.MOD(IA,2).EQ.0) IKFI=2
1112             IF(IA.GT.10) IKFI=3
1113             FACQQH=FACQQH*PARU(150+10*IHIGG+IKFI)**2
1114             IF(IMSS(1).NE.0.AND.IA.EQ.5) THEN
1115               FACQQH=FACQQH/(1D0+RMSS(41))**2
1116               IF(IHIGG.NE.3) THEN
1117                 FACQQH=FACQQH*(1D0+RMSS(41)*PARU(152+10*IHIGG)/
1118      &          PARU(151+10*IHIGG))**2
1119               ENDIF
1120             ENDIF
1121           ENDIF
1122           CALL PYQQBH(WTQQBH)
1123           CALL PYWIDT(KFHIGG,SH,WDTP,WDTE)
1124           HS=SHR*WDTP(0)
1125           HF=SHR*(WDTE(0,1)+WDTE(0,2)+WDTE(0,4))
1126           FACBW=(1D0/PARU(1))*VINT(2)*HF/((SH-SQMH)**2+HS**2)
1127           IF(ABS(SHR-PMAS(KFHIGG,1)).GT.PARP(48)*PMAS(KFHIGG,2))
1128      &    FACBW=0D0
1129           DO 470 I=MMINA,MMAXA
1130             IF(I.EQ.0.OR.IABS(I).GT.MSTP(58).OR.
1131      &      KFAC(1,I)*KFAC(2,-I).EQ.0) GOTO 470
1132             NCHN=NCHN+1
1133             ISIG(NCHN,1)=I
1134             ISIG(NCHN,2)=-I
1135             ISIG(NCHN,3)=1
1136             SIGH(NCHN)=FACQQH*WTQQBH*FACBW
1137   470     CONTINUE
1138
1139         ELSEIF(ISUB.EQ.123) THEN
1140 C...f + f' -> f + f' + h0 (or H0, or A0) (Z0 + Z0 -> h0 as
1141 C...inner process)
1142           FACNOR=COMFAC*(4D0*PARU(1)*AEM/(XW*XW1))**3*SQMZ/32D0
1143           IF(MSTP(4).GE.1.OR.IHIGG.GE.2) FACNOR=FACNOR*
1144      &    PARU(154+10*IHIGG)**2
1145           FACPRP=1D0/((VINT(215)-VINT(204)**2)*
1146      &    (VINT(216)-VINT(209)**2))**2
1147           FACZZ1=FACNOR*FACPRP*(0.5D0*TAUP*VINT(2))*VINT(219)
1148           FACZZ2=FACNOR*FACPRP*VINT(217)*VINT(218)
1149           CALL PYWIDT(KFHIGG,SH,WDTP,WDTE)
1150           HS=SHR*WDTP(0)
1151           HF=SHR*(WDTE(0,1)+WDTE(0,2)+WDTE(0,4))
1152           FACBW=(1D0/PARU(1))*VINT(2)*HF/((SH-SQMH)**2+HS**2)
1153           IF(ABS(SHR-PMAS(KFHIGG,1)).GT.PARP(48)*PMAS(KFHIGG,2))
1154      &    FACBW=0D0
1155           DO 490 I=MMIN1,MMAX1
1156             IF(I.EQ.0.OR.KFAC(1,I).EQ.0) GOTO 490
1157             IA=IABS(I)
1158             DO 480 J=MMIN2,MMAX2
1159               IF(J.EQ.0.OR.KFAC(2,J).EQ.0) GOTO 480
1160               JA=IABS(J)
1161               EI=KCHG(IA,1)*ISIGN(1,I)/3D0
1162               AI=SIGN(1D0,KCHG(IA,1)+0.5D0)*ISIGN(1,I)
1163               VI=AI-4D0*EI*XWV
1164               EJ=KCHG(JA,1)*ISIGN(1,J)/3D0
1165               AJ=SIGN(1D0,KCHG(JA,1)+0.5D0)*ISIGN(1,J)
1166               VJ=AJ-4D0*EJ*XWV
1167               FACLR1=(VI**2+AI**2)*(VJ**2+AJ**2)+4D0*VI*AI*VJ*AJ
1168               FACLR2=(VI**2+AI**2)*(VJ**2+AJ**2)-4D0*VI*AI*VJ*AJ
1169               NCHN=NCHN+1
1170               ISIG(NCHN,1)=I
1171               ISIG(NCHN,2)=J
1172               ISIG(NCHN,3)=1
1173               SIGH(NCHN)=(FACLR1*FACZZ1+FACLR2*FACZZ2)*FACBW
1174   480       CONTINUE
1175   490     CONTINUE
1176
1177         ELSEIF(ISUB.EQ.124) THEN
1178 C...f + f' -> f" + f"' + h0 (or H0, or A0) (W+ + W- -> h0 as
1179 C...inner process)
1180           FACNOR=COMFAC*(4D0*PARU(1)*AEM/XW)**3*SQMW
1181           IF(MSTP(4).GE.1.OR.IHIGG.GE.2) FACNOR=FACNOR*
1182      &    PARU(155+10*IHIGG)**2
1183           FACPRP=1D0/((VINT(215)-VINT(204)**2)*
1184      &    (VINT(216)-VINT(209)**2))**2
1185           FACWW=FACNOR*FACPRP*(0.5D0*TAUP*VINT(2))*VINT(219)
1186           CALL PYWIDT(KFHIGG,SH,WDTP,WDTE)
1187           HS=SHR*WDTP(0)
1188           HF=SHR*(WDTE(0,1)+WDTE(0,2)+WDTE(0,4))
1189           FACBW=(1D0/PARU(1))*VINT(2)*HF/((SH-SQMH)**2+HS**2)
1190           IF(ABS(SHR-PMAS(KFHIGG,1)).GT.PARP(48)*PMAS(KFHIGG,2))
1191      &    FACBW=0D0
1192           DO 510 I=MMIN1,MMAX1
1193             IF(I.EQ.0.OR.KFAC(1,I).EQ.0) GOTO 510
1194             EI=SIGN(1D0,DBLE(I))*KCHG(IABS(I),1)
1195             DO 500 J=MMIN2,MMAX2
1196               IF(J.EQ.0.OR.KFAC(2,J).EQ.0) GOTO 500
1197               EJ=SIGN(1D0,DBLE(J))*KCHG(IABS(J),1)
1198               IF(EI*EJ.GT.0D0) GOTO 500
1199               FACLR=VINT(180+I)*VINT(180+J)
1200               NCHN=NCHN+1
1201               ISIG(NCHN,1)=I
1202               ISIG(NCHN,2)=J
1203               ISIG(NCHN,3)=1
1204               SIGH(NCHN)=FACLR*FACWW*FACBW
1205   500       CONTINUE
1206   510     CONTINUE
1207
1208         ELSEIF(ISUB.EQ.143) THEN
1209 C...f + fbar' -> H+/-
1210           SQMHC=PMAS(37,1)**2
1211           CALL PYWIDT(37,SH,WDTP,WDTE)
1212           HS=SHR*WDTP(0)
1213           FACBW=4D0*COMFAC/((SH-SQMHC)**2+HS**2)
1214           HP=AEM/(8D0*XW)*SH/SQMW*SH
1215           DO 530 I=MMIN1,MMAX1
1216             IF(I.EQ.0.OR.KFAC(1,I).EQ.0) GOTO 530
1217             IA=IABS(I)
1218             IM=(MOD(IA,10)+1)/2
1219             DO 520 J=MMIN2,MMAX2
1220               IF(J.EQ.0.OR.KFAC(2,J).EQ.0) GOTO 520
1221               JA=IABS(J)
1222               JM=(MOD(JA,10)+1)/2
1223               IF(I*J.GT.0.OR.IA.EQ.JA.OR.IM.NE.JM) GOTO 520
1224               IF((IA.LE.10.AND.JA.GT.10).OR.(IA.GT.10.AND.JA.LE.10))
1225      &        GOTO 520
1226               IF(MOD(IA,2).EQ.0) THEN
1227                 IU=IA
1228                 IL=JA
1229               ELSE
1230                 IU=JA
1231                 IL=IA
1232               ENDIF
1233               RML=PYMRUN(IL,SH)**2/SH
1234               RMU=PYMRUN(IU,SH)**2/SH
1235               HI=HP*(RML*PARU(141)**2+RMU/PARU(141)**2)
1236               IF(IA.LE.10) HI=HI*FACA/3D0
1237               KCHHC=(KCHG(IA,1)*ISIGN(1,I)+KCHG(JA,1)*ISIGN(1,J))/3
1238               HF=SHR*(WDTE(0,1)+WDTE(0,(5-KCHHC)/2)+WDTE(0,4))
1239               NCHN=NCHN+1
1240               ISIG(NCHN,1)=I
1241               ISIG(NCHN,2)=J
1242               ISIG(NCHN,3)=1
1243               SIGH(NCHN)=HI*FACBW*HF
1244   520       CONTINUE
1245   530     CONTINUE
1246
1247         ELSEIF(ISUB.EQ.161) THEN
1248 C...f + g -> f' + H+/- (b + g -> t + H+/- only)
1249 C...(choice of only b and t to avoid kinematics problems)
1250           FHCQ=COMFAC*FACA*AS*AEM/XW*1D0/24
1251 C...H propagator: as simulated in PYOFSH and as desired
1252           SQMHC=PMAS(37,1)**2
1253           GMMHC=PMAS(37,1)*PMAS(37,2)
1254           HBW4=GMMHC/((SQM4-SQMHC)**2+GMMHC**2)
1255           CALL PYWIDT(37,SQM4,WDTP,WDTE)
1256           GMMHCC=SQRT(SQM4)*WDTP(0)
1257           HBW4C=GMMHCC/((SQM4-SQMHC)**2+GMMHCC**2)
1258           FHCQ=FHCQ*HBW4C/HBW4
1259           Q2RM=SH
1260           IF(MSTP(32).EQ.12) Q2RM=PARP(194)
1261           DO 550 I=MMINA,MMAXA
1262             IA=IABS(I)
1263             IF(IA.NE.5) GOTO 550
1264             SQML=PYMRUN(IA,Q2RM)**2
1265             IUA=IA+MOD(IA,2)
1266             SQMQ=PYMRUN(IUA,Q2RM)**2
1267             FACHCQ=FHCQ*(SQML*PARU(141)**2+SQMQ/PARU(141)**2)/SQMW*
1268      &      (SH/(SQMQ-UH)+2D0*SQMQ*(SQMHC-UH)/(SQMQ-UH)**2+(SQMQ-UH)/SH-
1269      &      2D0*SQMQ/(SQMQ-UH)+2D0*(SQMHC-UH)/(SQMQ-UH)*
1270      &      (SQMHC-SQMQ-SH)/SH)
1271             KCHHC=ISIGN(1,KCHG(IA,1)*ISIGN(1,I))
1272             DO 540 ISDE=1,2
1273               IF(ISDE.EQ.1.AND.KFAC(1,I)*KFAC(2,21).EQ.0) GOTO 540
1274               IF(ISDE.EQ.2.AND.KFAC(1,21)*KFAC(2,I).EQ.0) GOTO 540
1275               NCHN=NCHN+1
1276               ISIG(NCHN,ISDE)=I
1277               ISIG(NCHN,3-ISDE)=21
1278               ISIG(NCHN,3)=1
1279               SIGH(NCHN)=FACHCQ*WIDS(37,(5-KCHHC)/2)
1280               IF(IUA.EQ.6) SIGH(NCHN)=SIGH(NCHN)*WIDS(6,(5+KCHHC)/2)
1281   540       CONTINUE
1282   550     CONTINUE
1283         ENDIF
1284
1285       ELSEIF(ISUB.LE.402) THEN
1286         IF(ISUB.EQ.401) THEN
1287 C...  g + g -> t + bbar + H-
1288           IF(KFAC(1,21)*KFAC(2,21).EQ.0) GOTO 560
1289           IA=KFPR(ISUBSV,2)
1290           CALL PYSTBH(WTTBH)
1291           CALL PYWIDT(KFHIGG,SH,WDTP,WDTE)
1292           HS=SHR*WDTP(0)
1293           FACBW=(1D0/PARU(1))*VINT(2)*HS/((SH-SQMH)**2+HS**2)
1294           IF(ABS(SHR-PMAS(KFHIGG,1)).GT.PARP(48)*PMAS(KFHIGG,2))
1295      &       FACBW=0D0
1296           NCHN=NCHN+1
1297           ISIG(NCHN,1)=21
1298           ISIG(NCHN,2)=21
1299           ISIG(NCHN,3)=1
1300           SIGH(NCHN)=2d0*COMFAC*WTTBH*FACBW
1301 c     Since we don't know yet if H+ or H-, assume H+
1302 c     when calculating suppression due to closed channels.
1303           SIGH(NCHN)=SIGH(NCHN)*WIDS(37,2)*WIDS(6,3)
1304           IF(ABS(WIDS(37,2)-WIDS(37,3))
1305      &       .GE.1D-6*(WIDS(37,2)+WIDS(37,3)).OR.
1306      &       ABS(WIDS(6,2)-WIDS(6,3))
1307      &       .GE.1D-6*(WIDS(6,2)+WIDS(6,3))) THEN
1308             WRITE(*,*)'Error: Process 401 cannot handle different'
1309             WRITE(*,*)'decays for H+ and H- or t and tbar.'
1310             WRITE(*,*)'Execution stopped.'
1311             CALL PYSTOP(108)
1312           END IF
1313  560      CONTINUE
1314
1315         ELSEIF(ISUB.EQ.402) THEN
1316 C...  q + qbar -> t + bbar + H-
1317           IA=KFPR(ISUBSV,2)
1318           CALL PYSTBH(WTTBH)
1319           CALL PYWIDT(KFHIGG,SH,WDTP,WDTE)
1320           HS=SHR*WDTP(0)
1321           FACBW=(1D0/PARU(1))*VINT(2)*HS/((SH-SQMH)**2+HS**2)
1322           IF(ABS(SHR-PMAS(KFHIGG,1)).GT.PARP(48)*PMAS(KFHIGG,2))
1323      &       FACBW=0D0
1324           DO 570 I=MMINA,MMAXA
1325             IF(I.EQ.0.OR.IABS(I).GT.MSTP(58).OR.
1326      &         KFAC(1,I)*KFAC(2,-I).EQ.0) GOTO 570
1327             NCHN=NCHN+1
1328             ISIG(NCHN,1)=I
1329             ISIG(NCHN,2)=-I
1330             ISIG(NCHN,3)=1
1331             SIGH(NCHN)=2d0*COMFAC*WTTBH*FACBW
1332 c     Since we don't know yet if H+ or H-, assume H+
1333 c     when calculating suppression due to closed channels.
1334             SIGH(NCHN)=SIGH(NCHN)*WIDS(37,2)*WIDS(6,3)
1335             IF(ABS(WIDS(37,2)-WIDS(37,3))/(WIDS(37,2)+WIDS(37,3))
1336      &         .GE.1D-6.OR.
1337      &         ABS(WIDS(6,2)-WIDS(6,3))/(WIDS(6,2)+WIDS(6,3))
1338      &         .GE.1D-6) THEN
1339               WRITE(*,*)'Error: Process 402 cannot handle different'
1340               WRITE(*,*)'decays for H+ and H- or t and tbar.'
1341               WRITE(*,*)'Execution stopped.'
1342               CALL PYSTOP(108)
1343             END IF
1344  570      CONTINUE
1345         ENDIF
1346       ENDIF
1347
1348       RETURN
1349       END