pyhisigh.f

Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file pyhisigh.f

    
C***********************************************************************    
    
      SUBROUTINE pyhisigh(NCHN,SIGS)  
    
C...Differential matrix elements for all included subprocesses. 
C...Note that what is coded is (disregarding the COMFAC factor) 
C...1) for 2 -> 1 processes: s-hat/pi*d(sigma-hat), where,  
C...when d(sigma-hat) is given in the zero-width limit, the delta   
C...function in tau is replaced by a Breit-Wigner:  
C...1/pi*(s*m_res*Gamma_res)/((s*tau-m_res^2)^2+(m_res*Gamma_res)^2);   
C...2) for 2 -> 2 processes: (s-hat)**2/pi*d(sigma-hat)/d(t-hat);   
C...i.e., dimensionless quantities. COMFAC contains the factor  
C...pi/s and the conversion factor from GeV^-2 to mb.   
      common/ludat1/mstu(200),paru(200),mstj(200),parj(200) 
      SAVE /ludat1/ 
      common/ludat2/kchg(500,3),pmas(500,4),parf(2000),vckm(4,4)    
      SAVE /ludat2/ 
      common/ludat3/mdcy(500,3),mdme(2000,2),brat(2000),kfdp(2000,5)    
      SAVE /ludat3/ 
      common/pyhisubs/msel,msub(200),kfin(2,-40:40),ckin(200) 
      SAVE /pyhisubs/ 
      common/pyhipars/mstp(200),parp(200),msti(200),pari(200) 
      SAVE /pyhipars/ 
      common/pyhiint1/mint(400),vint(400) 
      SAVE /pyhiint1/ 
      common/pyhiint2/iset(200),kfpr(200,2),coef(200,20),icol(40,4,2) 
      SAVE /pyhiint2/ 
      common/pyhiint3/xsfx(2,-40:40),isig(1000,3),sigh(1000)  
      SAVE /pyhiint3/ 
      common/pyhiint4/widp(21:40,0:40),wide(21:40,0:40),wids(21:40,3) 
      SAVE /pyhiint4/ 
      common/pyhiint5/ngen(0:200,3),xsec(0:200,3) 
      SAVE /pyhiint5/ 
      dimension x(2),xpq(-6:6),kfac(2,-40:40),wdtp(0:40),wdte(0:40,0:5) 
    
C...Reset number of channels and cross-section. 
      nchn=0    
      sigs=0.   
    
C...Read kinematical variables and limits.  
      isub=mint(1)  
      taumin=vint(11)   
      ystmin=vint(12)   
      ctnmin=vint(13)   
      ctpmin=vint(14)   
      xt2min=vint(15)   
      taupmn=vint(16)   
      tau=vint(21)  
      yst=vint(22)  
      cth=vint(23)  
      xt2=vint(25)  
      taup=vint(26) 
      taumax=vint(31)   
      ystmax=vint(32)   
      ctnmax=vint(33)   
      ctpmax=vint(34)   
      xt2max=vint(35)   
      taupmx=vint(36)   
    
C...Derive kinematical quantities.  
      IF(iset(isub).LE.2.OR.iset(isub).EQ.5) THEN   
        x(1)=sqrt(tau)*exp(yst) 
        x(2)=sqrt(tau)*exp(-yst)    
      ELSE  
        x(1)=sqrt(taup)*exp(yst)    
        x(2)=sqrt(taup)*exp(-yst)   
      ENDIF 
      IF(mint(43).EQ.4.AND.iset(isub).GE.1.AND. 
     &(x(1).GT.0.999.OR.x(2).GT.0.999)) RETURN  
      sh=tau*vint(2)    
      sqm3=vint(63) 
      sqm4=vint(64) 
      rm3=sqm3/sh   
      rm4=sqm4/sh   
      be34=sqrt((1.-rm3-rm4)**2-4.*rm3*rm4) 
      rpts=4.*vint(71)**2/sh    
      be34l=sqrt(max(0.,(1.-rm3-rm4)**2-4.*rm3*rm4-rpts))   
      rm34=2.*rm3*rm4   
      rsqm=1.+rm34  
      rthm=(4.*rm3*rm4+rpts)/(1.-rm3-rm4+be34l) 
      th=-0.5*sh*max(rthm,1.-rm3-rm4-be34*cth)  
      uh=-0.5*sh*max(rthm,1.-rm3-rm4+be34*cth)  
      sqpth=0.25*sh*be34**2*(1.-cth**2) 
      sh2=sh**2 
      th2=th**2 
      uh2=uh**2 
    
C...Choice of Q2 scale. 
      IF(iset(isub).EQ.1.OR.iset(isub).EQ.3) THEN   
        q2=sh   
      ELSEIF(mod(iset(isub),2).EQ.0.OR.iset(isub).EQ.5) THEN    
        IF(mstp(32).EQ.1) THEN  
          q2=2.*sh*th*uh/(sh**2+th**2+uh**2)    
        ELSEIF(mstp(32).EQ.2) THEN  
          q2=sqpth+0.5*(sqm3+sqm4)  
        ELSEIF(mstp(32).EQ.3) THEN  
          q2=min(-th,-uh)   
        ELSEIF(mstp(32).EQ.4) THEN  
          q2=sh 
        ENDIF   
        IF(iset(isub).EQ.5.AND.mstp(82).GE.2) q2=q2+parp(82)**2 
      ENDIF 
    
C...Store derived kinematical quantities.   
      vint(41)=x(1) 
      vint(42)=x(2) 
      vint(44)=sh   
      vint(43)=sqrt(sh) 
      vint(45)=th   
      vint(46)=uh   
      vint(48)=sqpth    
      vint(47)=sqrt(sqpth)  
      vint(50)=taup*vint(2) 
      vint(49)=sqrt(max(0.,vint(50)))   
      vint(52)=q2   
      vint(51)=sqrt(q2) 
    
C...Calculate parton structure functions.   
      IF(iset(isub).LE.0) goto 145  
      IF(mint(43).GE.2) THEN    
        q2sf=q2 
        IF(iset(isub).EQ.3.OR.iset(isub).EQ.4) THEN 
          q2sf=pmas(23,1)**2    
          IF(isub.EQ.8.OR.isub.EQ.76.OR.isub.EQ.77) q2sf=pmas(24,1)**2  
        ENDIF   
        DO 100 i=3-mint(41),mint(42)    
        xsf=x(i)    
        IF(iset(isub).EQ.5) xsf=x(i)/vint(142+i)    
        CALL pyhistfu(mint(10+i),xsf,q2sf,xpq,i)    
        DO 100 kfl=-6,6 
  100   xsfx(i,kfl)=xpq(kfl)
      ENDIF 
    
C...Calculate alpha_strong and K-factor.    
      IF(mstp(33).NE.3) as=ulalps(q2)   
      fack=1.   
      faca=1.   
      IF(mstp(33).EQ.1) THEN    
        fack=parp(31)   
      ELSEIF(mstp(33).EQ.2) THEN    
        fack=parp(31)   
        faca=parp(32)/parp(31)  
      ELSEIF(mstp(33).EQ.3) THEN    
        q2as=parp(33)*q2    
        IF(iset(isub).EQ.5.AND.mstp(82).GE.2) q2as=q2as+    
     &  paru(112)*parp(82)  
        as=ulalps(q2as) 
      ENDIF 
      radc=1.+as/paru(1)    
    
C...Set flags for allowed reacting partons/leptons. 
      DO 130 i=1,2  
      DO 110 j=-40,40   
  110 kfac(i,j)=0   
      IF(mint(40+i).EQ.1) THEN  
        kfac(i,mint(10+i))=1    
      ELSE  
        DO 120 j=-40,40 
        kfac(i,j)=kfin(i,j) 
        IF(abs(j).GT.mstp(54).AND.j.NE.21) kfac(i,j)=0  
        IF(abs(j).LE.6) THEN    
          IF(xsfx(i,j).LT.1.e-10) kfac(i,j)=0   
        ELSEIF(j.EQ.21) THEN    
          IF(xsfx(i,0).LT.1.e-10) kfac(i,21)=0  
        ENDIF   
  120   CONTINUE    
      ENDIF 
  130 CONTINUE  
    
C...Lower and upper limit for flavour loops.    
      min1=0    
      max1=0    
      min2=0    
      max2=0    
      DO 140 j=-20,20   
      IF(kfac(1,-j).EQ.1) min1=-j   
      IF(kfac(1,j).EQ.1) max1=j 
      IF(kfac(2,-j).EQ.1) min2=-j   
      IF(kfac(2,j).EQ.1) max2=j 
  140 CONTINUE  
      mina=min(min1,min2)   
      maxa=max(max1,max2)   
    
C...Common conversion factors (including Jacobian) for subprocesses.    
      sqmz=pmas(23,1)**2    
      gmmz=pmas(23,1)*pmas(23,2)    
      sqmw=pmas(24,1)**2    
      gmmw=pmas(24,1)*pmas(24,2)    
      sqmh=pmas(25,1)**2    
      gmmh=pmas(25,1)*pmas(25,2)    
      sqmzp=pmas(32,1)**2   
      gmmzp=pmas(32,1)*pmas(32,2)   
      sqmhc=pmas(37,1)**2   
      gmmhc=pmas(37,1)*pmas(37,2)   
      sqmr=pmas(40,1)**2    
      gmmr=pmas(40,1)*pmas(40,2)    
      aem=paru(101) 
      xw=paru(102)  
    
C...Phase space integral in tau and y*. 
      comfac=paru(1)*paru(5)/vint(2)    
      IF(mint(43).EQ.4) comfac=comfac*fack  
      IF((mint(43).GE.2.OR.iset(isub).EQ.3.OR.iset(isub).EQ.4).AND. 
     &iset(isub).NE.5) THEN 
        atau0=log(taumax/taumin)    
        atau1=(taumax-taumin)/(taumax*taumin)   
        h1=coef(isub,1)+(atau0/atau1)*coef(isub,2)/tau  
        IF(mint(72).GE.1) THEN  
          taur1=vint(73)    
          gamr1=vint(74)    
          atau2=log(taumax/taumin*(taumin+taur1)/(taumax+taur1))/taur1  
          atau3=(atan((taumax-taur1)/gamr1)-atan((taumin-taur1)/gamr1))/    
     &    gamr1 
          h1=h1+(atau0/atau2)*coef(isub,3)/(tau+taur1)+ 
     &    (atau0/atau3)*coef(isub,4)*tau/((tau-taur1)**2+gamr1**2)  
        ENDIF   
        IF(mint(72).EQ.2) THEN  
          taur2=vint(75)    
          gamr2=vint(76)    
          atau4=log(taumax/taumin*(taumin+taur2)/(taumax+taur2))/taur2  
          atau5=(atan((taumax-taur2)/gamr2)-atan((taumin-taur2)/gamr2))/    
     &    gamr2 
          h1=h1+(atau0/atau4)*coef(isub,5)/(tau+taur2)+ 
     &    (atau0/atau5)*coef(isub,6)*tau/((tau-taur2)**2+gamr2**2)  
        ENDIF   
        comfac=comfac*atau0/(tau*h1)    
      ENDIF 
      IF(mint(43).EQ.4.AND.iset(isub).NE.5) THEN    
        ayst0=ystmax-ystmin 
        ayst1=0.5*(ystmax-ystmin)**2    
        ayst2=ayst1 
        ayst3=2.*(atan(exp(ystmax))-atan(exp(ystmin)))  
        h2=(ayst0/ayst1)*coef(isub,7)*(yst-ystmin)+(ayst0/ayst2)*   
     &  coef(isub,8)*(ystmax-yst)+(ayst0/ayst3)*coef(isub,9)/cosh(yst)  
        comfac=comfac*ayst0/h2  
      ENDIF 
    
C...2 -> 1 processes: reduction in angular part of phase space integral 
C...for case of decaying resonance. 
      acth0=ctnmax-ctnmin+ctpmax-ctpmin 
      IF((iset(isub).EQ.1.OR.iset(isub).EQ.3).AND.  
     &mdcy(kfpr(isub,1),1).EQ.1) THEN   
        IF(kfpr(isub,1).EQ.25.OR.kfpr(isub,1).EQ.37) THEN   
          comfac=comfac*0.5*acth0   
        ELSE    
          comfac=comfac*0.125*(3.*acth0+ctnmax**3-ctnmin**3+    
     &    ctpmax**3-ctpmin**3)  
        ENDIF   
    
C...2 -> 2 processes: angular part of phase space integral. 
      ELSEIF(iset(isub).EQ.2.OR.iset(isub).EQ.4) THEN   
        acth1=log((max(rm34,rsqm-ctnmin)*max(rm34,rsqm-ctpmin))/    
     &  (max(rm34,rsqm-ctnmax)*max(rm34,rsqm-ctpmax)))  
        acth2=log((max(rm34,rsqm+ctnmax)*max(rm34,rsqm+ctpmax))/    
     &  (max(rm34,rsqm+ctnmin)*max(rm34,rsqm+ctpmin)))  
        acth3=1./max(rm34,rsqm-ctnmax)-1./max(rm34,rsqm-ctnmin)+    
     &  1./max(rm34,rsqm-ctpmax)-1./max(rm34,rsqm-ctpmin)   
        acth4=1./max(rm34,rsqm+ctnmin)-1./max(rm34,rsqm+ctnmax)+    
     &  1./max(rm34,rsqm+ctpmin)-1./max(rm34,rsqm+ctpmax)   
        h3=coef(isub,10)+   
     &  (acth0/acth1)*coef(isub,11)/max(rm34,rsqm-cth)+ 
     &  (acth0/acth2)*coef(isub,12)/max(rm34,rsqm+cth)+ 
     &  (acth0/acth3)*coef(isub,13)/max(rm34,rsqm-cth)**2+  
     &  (acth0/acth4)*coef(isub,14)/max(rm34,rsqm+cth)**2   
        comfac=comfac*acth0*0.5*be34/h3 
      ENDIF 
    
C...2 -> 3, 4 processes: phace space integral in tau'.  
      IF(mint(43).GE.2.AND.(iset(isub).EQ.3.OR.iset(isub).EQ.4)) THEN   
        ataup0=log(taupmx/taupmn)   
        ataup1=((1.-tau/taupmx)**4-(1.-tau/taupmn)**4)/(4.*tau) 
        h4=coef(isub,15)+   
     &  ataup0/ataup1*coef(isub,16)/taup*(1.-tau/taup)**3   
        IF(1.-tau/taup.GT.1.e-4) THEN   
          fzw=(1.+tau/taup)*log(taup/tau)-2.*(1.-tau/taup)  
        ELSE    
          fzw=1./6.*(1.-tau/taup)**3*tau/taup   
        ENDIF   
        comfac=comfac*ataup0*fzw/h4 
      ENDIF 
    
C...Phase space integral for low-pT and multiple interactions.  
      IF(iset(isub).EQ.5) THEN  
        comfac=paru(1)*paru(5)*fack*0.5*vint(2)/sh2 
        atau0=log(2.*(1.+sqrt(1.-xt2))/xt2-1.)  
        atau1=2.*atan(1./xt2-1.)/sqrt(xt2)  
        h1=coef(isub,1)+(atau0/atau1)*coef(isub,2)/sqrt(tau)    
        comfac=comfac*atau0/h1  
        ayst0=ystmax-ystmin 
        ayst1=0.5*(ystmax-ystmin)**2    
        ayst3=2.*(atan(exp(ystmax))-atan(exp(ystmin)))  
        h2=(ayst0/ayst1)*coef(isub,7)*(yst-ystmin)+(ayst0/ayst1)*   
     &  coef(isub,8)*(ystmax-yst)+(ayst0/ayst3)*coef(isub,9)/cosh(yst)  
        comfac=comfac*ayst0/h2  
        IF(mstp(82).LE.1) comfac=comfac*xt2**2*(1./vint(149)-1.)    
C...For MSTP(82)>=2 an additional factor (xT2/(xT2+VINT(149))**2 is 
C...introduced to make cross-section finite for xT2 -> 0.   
        IF(mstp(82).GE.2) comfac=comfac*xt2**2/(vint(149)*  
     &  (1.+vint(149))) 
      ENDIF 
    
C...A: 2 -> 1, tree diagrams.   
    
  145 IF(isub.LE.10) THEN   
      IF(isub.EQ.1) THEN    
C...f + fb -> gamma*/Z0.    
        mint(61)=2  
        CALL pyhiwidt(23,sqrt(sh),wdtp,wdte)  
        facz=comfac*aem**2*4./3.    
        DO 150 i=mina,maxa  
        IF(i.EQ.0.OR.kfac(1,i)*kfac(2,-i).EQ.0) goto 150    
        ei=kchg(iabs(i),1)/3.   
        ai=sign(1.,ei)  
        vi=ai-4.*ei*xw  
        facf=1. 
        IF(iabs(i).LE.10) facf=faca/3.  
        nchn=nchn+1 
        isig(nchn,1)=i  
        isig(nchn,2)=-i 
        isig(nchn,3)=1  
        sigh(nchn)=facf*facz*(ei**2*vint(111)+ei*vi/(8.*xw*(1.-xw))*    
     &  sh*(sh-sqmz)/((sh-sqmz)**2+gmmz**2)*vint(112)+(vi**2+ai**2)/    
     &  (16.*xw*(1.-xw))**2*sh2/((sh-sqmz)**2+gmmz**2)*vint(114))   
  150   CONTINUE    
    
      ELSEIF(isub.EQ.2) THEN    
C...f + fb' -> W+/-.    
        CALL pyhiwidt(24,sqrt(sh),wdtp,wdte)  
        facw=comfac*(aem/xw)**2*1./24*sh2/((sh-sqmw)**2+gmmw**2)    
        DO 170 i=min1,max1  
        IF(i.EQ.0.OR.kfac(1,i).EQ.0) goto 170   
        ia=iabs(i)  
        DO 160 j=min2,max2  
        IF(j.EQ.0.OR.kfac(2,j).EQ.0) goto 160   
        ja=iabs(j)  
        IF(i*j.GT.0.OR.mod(ia+ja,2).EQ.0) goto 160  
        IF((ia.LE.10.AND.ja.GT.10).OR.(ia.GT.10.AND.ja.LE.10)) goto 160 
        kchw=(kchg(ia,1)*isign(1,i)+kchg(ja,1)*isign(1,j))/3    
        facf=1. 
        IF(ia.LE.10) facf=vckm((ia+1)/2,(ja+1)/2)*faca/3.   
        nchn=nchn+1 
        isig(nchn,1)=i  
        isig(nchn,2)=j  
        isig(nchn,3)=1  
        sigh(nchn)=facf*facw*(wdte(0,1)+wdte(0,(5-kchw)/2)+wdte(0,4))   
  160   CONTINUE    
  170   CONTINUE    
    
      ELSEIF(isub.EQ.3) THEN    
C...f + fb -> H0.   
        CALL pyhiwidt(25,sqrt(sh),wdtp,wdte)  
        fach=comfac*(aem/xw)**2*1./48.*(sh/sqmw)**2*    
     &  sh2/((sh-sqmh)**2+gmmh**2)*(wdte(0,1)+wdte(0,2)+wdte(0,4))  
        DO 180 i=mina,maxa  
        IF(i.EQ.0.OR.kfac(1,i)*kfac(2,-i).EQ.0) goto 180    
        rmq=pmas(iabs(i),1)**2/sh   
        nchn=nchn+1 
        isig(nchn,1)=i  
        isig(nchn,2)=-i 
        isig(nchn,3)=1  
        sigh(nchn)=fach*rmq*sqrt(max(0.,1.-4.*rmq)) 
  180   CONTINUE    
    
      ELSEIF(isub.EQ.4) THEN    
C...gamma + W+/- -> W+/-.   
    
      ELSEIF(isub.EQ.5) THEN    
C...Z0 + Z0 -> H0.  
        CALL pyhiwidt(25,sqrt(sh),wdtp,wdte)  
        fach=comfac*1./(128.*paru(1)**2*16.*(1.-xw)**3)*(aem/xw)**4*    
     &  (sh/sqmw)**2*sh2/((sh-sqmh)**2+gmmh**2)*    
     &  (wdte(0,1)+wdte(0,2)+wdte(0,4)) 
        DO 200 i=min1,max1  
        IF(i.EQ.0.OR.kfac(1,i).EQ.0) goto 200   
        DO 190 j=min2,max2  
        IF(j.EQ.0.OR.kfac(2,j).EQ.0) goto 190   
        ei=kchg(iabs(i),1)/3.   
        ai=sign(1.,ei)  
        vi=ai-4.*ei*xw  
        ej=kchg(iabs(j),1)/3.   
        aj=sign(1.,ej)  
        vj=aj-4.*ej*xw  
        nchn=nchn+1 
        isig(nchn,1)=i  
        isig(nchn,2)=j  
        isig(nchn,3)=1  
        sigh(nchn)=fach*(vi**2+ai**2)*(vj**2+aj**2) 
  190   CONTINUE    
  200   CONTINUE    
    
      ELSEIF(isub.EQ.6) THEN    
C...Z0 + W+/- -> W+/-.  
    
      ELSEIF(isub.EQ.7) THEN    
C...W+ + W- -> Z0.  
    
      ELSEIF(isub.EQ.8) THEN    
C...W+ + W- -> H0.  
        CALL pyhiwidt(25,sqrt(sh),wdtp,wdte)  
        fach=comfac*1./(128*paru(1)**2)*(aem/xw)**4*(sh/sqmw)**2*   
     &  sh2/((sh-sqmh)**2+gmmh**2)*(wdte(0,1)+wdte(0,2)+wdte(0,4))  
        DO 220 i=min1,max1  
        IF(i.EQ.0.OR.kfac(1,i).EQ.0) goto 220   
        ei=sign(1.,float(i))*kchg(iabs(i),1)    
        DO 210 j=min2,max2  
        IF(j.EQ.0.OR.kfac(2,j).EQ.0) goto 210   
        ej=sign(1.,float(j))*kchg(iabs(j),1)    
        IF(ei*ej.GT.0.) goto 210    
        nchn=nchn+1 
        isig(nchn,1)=i  
        isig(nchn,2)=j  
        isig(nchn,3)=1  
        sigh(nchn)=fach*vint(180+i)*vint(180+j) 
  210   CONTINUE    
  220   CONTINUE    
      ENDIF 
    
C...B: 2 -> 2, tree diagrams.   
    
      ELSEIF(isub.LE.20) THEN   
      IF(isub.EQ.11) THEN   
C...f + f' -> f + f'.   
        facqq1=comfac*as**2*4./9.*(sh2+uh2)/th2 
        facqqb=comfac*as**2*4./9.*((sh2+uh2)/th2*faca-  
     &  mstp(34)*2./3.*uh2/(sh*th)) 
        facqq2=comfac*as**2*4./9.*((sh2+th2)/uh2-   
     &  mstp(34)*2./3.*sh2/(th*uh)) 
        DO 240 i=min1,max1  
        IF(i.EQ.0.OR.kfac(1,i).EQ.0) goto 240   
        DO 230 j=min2,max2  
        IF(j.EQ.0.OR.kfac(2,j).EQ.0) goto 230   
        nchn=nchn+1 
        isig(nchn,1)=i  
        isig(nchn,2)=j  
        isig(nchn,3)=1  
        sigh(nchn)=facqq1   
        IF(i.EQ.-j) sigh(nchn)=facqqb   
        IF(i.EQ.j) THEN 
          sigh(nchn)=0.5*sigh(nchn) 
          nchn=nchn+1   
          isig(nchn,1)=i    
          isig(nchn,2)=j    
          isig(nchn,3)=2    
          sigh(nchn)=0.5*facqq2 
        ENDIF   
  230   CONTINUE    
  240   CONTINUE    
    
      ELSEIF(isub.EQ.12) THEN   
C...f + fb -> f' + fb' (q + qb -> q' + qb' only).   
        CALL pyhiwidt(21,sqrt(sh),wdtp,wdte)  
        facqqb=comfac*as**2*4./9.*(th2+uh2)/sh2*(wdte(0,1)+wdte(0,2)+   
     &  wdte(0,3)+wdte(0,4))    
        DO 250 i=mina,maxa  
        IF(i.EQ.0.OR.kfac(1,i)*kfac(2,-i).EQ.0) goto 250    
        nchn=nchn+1 
        isig(nchn,1)=i  
        isig(nchn,2)=-i 
        isig(nchn,3)=1  
        sigh(nchn)=facqqb   
  250   CONTINUE    
    
      ELSEIF(isub.EQ.13) THEN   
C...f + fb -> g + g (q + qb -> g + g only). 
        facgg1=comfac*as**2*32./27.*(uh/th-(2.+mstp(34)*1./4.)*uh2/sh2) 
        facgg2=comfac*as**2*32./27.*(th/uh-(2.+mstp(34)*1./4.)*th2/sh2) 
        DO 260 i=mina,maxa  
        IF(i.EQ.0.OR.kfac(1,i)*kfac(2,-i).EQ.0) goto 260    
        nchn=nchn+1 
        isig(nchn,1)=i  
        isig(nchn,2)=-i 
        isig(nchn,3)=1  
        sigh(nchn)=0.5*facgg1   
        nchn=nchn+1 
        isig(nchn,1)=i  
        isig(nchn,2)=-i 
        isig(nchn,3)=2  
        sigh(nchn)=0.5*facgg2   
  260   CONTINUE    
    
      ELSEIF(isub.EQ.14) THEN   
C...f + fb -> g + gamma (q + qb -> g + gamma only). 
        facgg=comfac*as*aem*8./9.*(th2+uh2)/(th*uh) 
        DO 270 i=mina,maxa  
        IF(i.EQ.0.OR.kfac(1,i)*kfac(2,-i).EQ.0) goto 270    
        ei=kchg(iabs(i),1)/3.   
        nchn=nchn+1 
        isig(nchn,1)=i  
        isig(nchn,2)=-i 
        isig(nchn,3)=1  
        sigh(nchn)=facgg*ei**2  
  270   CONTINUE    
    
      ELSEIF(isub.EQ.15) THEN   
C...f + fb -> g + Z0 (q + qb -> g + Z0 only).   
        faczg=comfac*as*aem/(xw*(1.-xw))*1./18.*    
     &  (th2+uh2+2.*sqm4*sh)/(th*uh)    
        faczg=faczg*wids(23,2)  
        DO 280 i=mina,maxa  
        IF(i.EQ.0.OR.kfac(1,i)*kfac(2,-i).EQ.0) goto 280    
        ei=kchg(iabs(i),1)/3.   
        ai=sign(1.,ei)  
        vi=ai-4.*ei*xw  
        nchn=nchn+1 
        isig(nchn,1)=i  
        isig(nchn,2)=-i 
        isig(nchn,3)=1  
        sigh(nchn)=faczg*(vi**2+ai**2)  
  280   CONTINUE    
    
      ELSEIF(isub.EQ.16) THEN   
C...f + fb' -> g + W+/- (q + qb' -> g + W+/- only). 
        facwg=comfac*as*aem/xw*2./9.*(th2+uh2+2.*sqm4*sh)/(th*uh)   
        DO 300 i=min1,max1  
        IF(i.EQ.0.OR.kfac(1,i).EQ.0) goto 300   
        ia=iabs(i)  
        DO 290 j=min2,max2  
        IF(j.EQ.0.OR.kfac(2,j).EQ.0) goto 290   
        ja=iabs(j)  
        IF(i*j.GT.0.OR.mod(ia+ja,2).EQ.0) goto 290  
        kchw=(kchg(ia,1)*isign(1,i)+kchg(ja,1)*isign(1,j))/3    
        fckm=1. 
        IF(mint(43).EQ.4) fckm=vckm((ia+1)/2,(ja+1)/2)  
        nchn=nchn+1 
        isig(nchn,1)=i  
        isig(nchn,2)=j  
        isig(nchn,3)=1  
        sigh(nchn)=facwg*fckm*wids(24,(5-kchw)/2)   
  290   CONTINUE    
  300   CONTINUE    
    
      ELSEIF(isub.EQ.17) THEN   
C...f + fb -> g + H0 (q + qb -> g + H0 only).   
    
      ELSEIF(isub.EQ.18) THEN   
C...f + fb -> gamma + gamma.    
        facgg=comfac*faca*aem**2*1./3.*(th2+uh2)/(th*uh)    
        DO 310 i=mina,maxa  
        IF(i.EQ.0.OR.kfac(1,i)*kfac(2,-i).EQ.0) goto 310    
        ei=kchg(iabs(i),1)/3.   
        nchn=nchn+1 
        isig(nchn,1)=i  
        isig(nchn,2)=-i 
        isig(nchn,3)=1  
        sigh(nchn)=facgg*ei**4  
  310   CONTINUE    
    
      ELSEIF(isub.EQ.19) THEN   
C...f + fb -> gamma + Z0.   
        facgz=comfac*faca*aem**2/(xw*(1.-xw))*1./24.*   
     &  (th2+uh2+2.*sqm4*sh)/(th*uh)    
        facgz=facgz*wids(23,2)  
        DO 320 i=mina,maxa  
        IF(i.EQ.0.OR.kfac(1,i)*kfac(2,-i).EQ.0) goto 320    
        ei=kchg(iabs(i),1)/3.   
        ai=sign(1.,ei)  
        vi=ai-4.*ei*xw  
        nchn=nchn+1 
        isig(nchn,1)=i  
        isig(nchn,2)=-i 
        isig(nchn,3)=1  
        sigh(nchn)=facgz*ei**2*(vi**2+ai**2)    
  320   CONTINUE    
    
      ELSEIF(isub.EQ.20) THEN   
C...f + fb' -> gamma + W+/-.    
        facgw=comfac*faca*aem**2/xw*1./6.*  
     &  ((2.*uh-th)/(3.*(sh-sqm4)))**2*(th2+uh2+2.*sqm4*sh)/(th*uh) 
        DO 340 i=min1,max1  
        IF(i.EQ.0.OR.kfac(1,i).EQ.0) goto 340   
        ia=iabs(i)  
        DO 330 j=min2,max2  
        IF(j.EQ.0.OR.kfac(2,j).EQ.0) goto 330   
        ja=iabs(j)  
        IF(i*j.GT.0.OR.mod(ia+ja,2).EQ.0) goto 330  
        kchw=(kchg(ia,1)*isign(1,i)+kchg(ja,1)*isign(1,j))/3    
        fckm=1. 
        IF(mint(43).EQ.4) fckm=vckm((ia+1)/2,(ja+1)/2)  
        nchn=nchn+1 
        isig(nchn,1)=i  
        isig(nchn,2)=j  
        isig(nchn,3)=1  
        sigh(nchn)=facgw*fckm*wids(24,(5-kchw)/2)   
  330   CONTINUE    
  340   CONTINUE    
      ENDIF 
    
      ELSEIF(isub.LE.30) THEN   
      IF(isub.EQ.21) THEN   
C...f + fb -> gamma + H0.   
    
      ELSEIF(isub.EQ.22) THEN   
C...f + fb -> Z0 + Z0.  
        faczz=comfac*faca*(aem/(xw*(1.-xw)))**2*1./768.*    
     &  (uh/th+th/uh+2.*(sqm3+sqm4)*sh/(th*uh)- 
     &  sqm3*sqm4*(1./th2+1./uh2))  
        faczz=faczz*wids(23,1)  
        DO 350 i=mina,maxa  
        IF(i.EQ.0.OR.kfac(1,i)*kfac(2,-i).EQ.0) goto 350    
        ei=kchg(iabs(i),1)/3.   
        ai=sign(1.,ei)  
        vi=ai-4.*ei*xw  
        nchn=nchn+1 
        isig(nchn,1)=i  
        isig(nchn,2)=-i 
        isig(nchn,3)=1  
        sigh(nchn)=faczz*(vi**4+6.*vi**2*ai**2+ai**4)   
  350   CONTINUE    
    
      ELSEIF(isub.EQ.23) THEN   
C...f + fb' -> Z0 + W+/-.   
        faczw=comfac*faca*(aem/xw)**2*1./6. 
        faczw=faczw*wids(23,2)  
        thuh=max(th*uh-sqm3*sqm4,sh*ckin(3)**2) 
        DO 370 i=min1,max1  
        IF(i.EQ.0.OR.kfac(1,i).EQ.0) goto 370   
        ia=iabs(i)  
        DO 360 j=min2,max2  
        IF(j.EQ.0.OR.kfac(2,j).EQ.0) goto 360   
        ja=iabs(j)  
        IF(i*j.GT.0.OR.mod(ia+ja,2).EQ.0) goto 360  
        kchw=(kchg(ia,1)*isign(1,i)+kchg(ja,1)*isign(1,j))/3    
        ei=kchg(ia,1)/3.    
        ai=sign(1.,ei)  
        vi=ai-4.*ei*xw  
        ej=kchg(ja,1)/3.    
        aj=sign(1.,ej)  
        vj=aj-4.*ej*xw  
        IF(vi+ai.GT.0) THEN 
          visav=vi  
          aisav=ai  
          vi=vj 
          ai=aj 
          vj=visav  
          aj=aisav  
        ENDIF   
        fckm=1. 
        IF(mint(43).EQ.4) fckm=vckm((ia+1)/2,(ja+1)/2)  
        nchn=nchn+1 
        isig(nchn,1)=i  
        isig(nchn,2)=j  
        isig(nchn,3)=1  
        sigh(nchn)=faczw*fckm*(1./(sh-sqmw)**2* 
     &  ((9.-8.*xw)/4.*thuh+(8.*xw-6.)/4.*sh*(sqm3+sqm4))+  
     &  (thuh-sh*(sqm3+sqm4))/(2.*(sh-sqmw))*((vj+aj)/th-(vi+ai)/uh)+   
     &  thuh/(16.*(1.-xw))*((vj+aj)**2/th2+(vi+ai)**2/uh2)+ 
     &  sh*(sqm3+sqm4)/(8.*(1.-xw))*(vi+ai)*(vj+aj)/(th*uh))*   
     &  wids(24,(5-kchw)/2) 
  360   CONTINUE    
  370   CONTINUE    
    
      ELSEIF(isub.EQ.24) THEN   
C...f + fb -> Z0 + H0.  
        thuh=max(th*uh-sqm3*sqm4,sh*ckin(3)**2) 
        fachz=comfac*faca*(aem/(xw*(1.-xw)))**2*1./96.* 
     &  (thuh+2.*sh*sqmz)/(sh-sqmz)**2  
        fachz=fachz*wids(23,2)*wids(25,2)   
        DO 380 i=mina,maxa  
        IF(i.EQ.0.OR.kfac(1,i)*kfac(2,-i).EQ.0) goto 380    
        ei=kchg(iabs(i),1)/3.   
        ai=sign(1.,ei)  
        vi=ai-4.*ei*xw  
        nchn=nchn+1 
        isig(nchn,1)=i  
        isig(nchn,2)=-i 
        isig(nchn,3)=1  
        sigh(nchn)=fachz*(vi**2+ai**2)  
  380   CONTINUE    
    
      ELSEIF(isub.EQ.25) THEN   
C...f + fb -> W+ + W-.  
        facww=comfac*faca*(aem/xw)**2*1./12.    
        facww=facww*wids(24,1)  
        thuh=max(th*uh-sqm3*sqm4,sh*ckin(3)**2) 
        DO 390 i=mina,maxa  
        IF(i.EQ.0.OR.kfac(1,i)*kfac(2,-i).EQ.0) goto 390    
        ei=kchg(iabs(i),1)/3.   
        ai=sign(1.,ei)  
        vi=ai-4.*ei*xw  
        dsigww=thuh/sh2*(3.-(sh-3.*(sqm3+sqm4))/(sh-sqmz)*  
     &  (vi+ai)/(2.*ai*(1.-xw))+(sh/(sh-sqmz))**2*  
     &  (1.-2.*(sqm3+sqm4)/sh+12.*sqm3*sqm4/sh2)*(vi**2+ai**2)/ 
     &  (8.*(1.-xw)**2))-2.*sqmz/(sh-sqmz)*(vi+ai)/ai+  
     &  sqmz*sh/(sh-sqmz)**2*(1.-2.*(sqm3+sqm4)/sh)*(vi**2+ai**2)/  
     &  (2.*(1.-xw))    
        IF(kchg(iabs(i),1).LT.0) THEN   
          dsigww=dsigww+2.*(1.+sqmz/(sh-sqmz)*(vi+ai)/(2.*ai))* 
     &    (thuh/(sh*th)-(sqm3+sqm4)/th)+thuh/th2    
        ELSE    
          dsigww=dsigww+2.*(1.+sqmz/(sh-sqmz)*(vi+ai)/(2.*ai))* 
     &    (thuh/(sh*uh)-(sqm3+sqm4)/uh)+thuh/uh2    
        ENDIF   
        nchn=nchn+1 
        isig(nchn,1)=i  
        isig(nchn,2)=-i 
        isig(nchn,3)=1  
        sigh(nchn)=facww*dsigww 
  390   CONTINUE    
    
      ELSEIF(isub.EQ.26) THEN   
C...f + fb' -> W+/- + H0.   
        thuh=max(th*uh-sqm3*sqm4,sh*ckin(3)**2) 
        fachw=comfac*faca*(aem/xw)**2*1./24.*(thuh+2.*sh*sqmw)/ 
     &  (sh-sqmw)**2    
        fachw=fachw*wids(25,2)  
        DO 410 i=min1,max1  
        IF(i.EQ.0.OR.kfac(1,i).EQ.0) goto 410   
        ia=iabs(i)  
        DO 400 j=min2,max2  
        IF(j.EQ.0.OR.kfac(1,j).EQ.0) goto 400   
        ja=iabs(j)  
        IF(i*j.GT.0.OR.mod(ia+ja,2).EQ.0) goto 400  
        kchw=(kchg(ia,1)*isign(1,i)+kchg(ja,1)*isign(1,j))/3    
        fckm=1. 
        IF(mint(43).EQ.4) fckm=vckm((ia+1)/2,(ja+1)/2)  
        nchn=nchn+1 
        isig(nchn,1)=i  
        isig(nchn,2)=j  
        isig(nchn,3)=1  
        sigh(nchn)=fachw*fckm*wids(24,(5-kchw)/2)   
  400   CONTINUE    
  410   CONTINUE    
    
      ELSEIF(isub.EQ.27) THEN   
C...f + fb -> H0 + H0.  
    
      ELSEIF(isub.EQ.28) THEN   
C...f + g -> f + g (q + g -> q + g only).   
        facqg1=comfac*as**2*4./9.*((2.+mstp(34)*1./4.)*uh2/th2-uh/sh)*  
     &  faca    
        facqg2=comfac*as**2*4./9.*((2.+mstp(34)*1./4.)*sh2/th2-sh/uh)   
        DO 430 i=mina,maxa  
        IF(i.EQ.0) goto 430 
        DO 420 isde=1,2 
        IF(isde.EQ.1.AND.kfac(1,i)*kfac(2,21).EQ.0) goto 420    
        IF(isde.EQ.2.AND.kfac(1,21)*kfac(2,i).EQ.0) goto 420    
        nchn=nchn+1 
        isig(nchn,isde)=i   
        isig(nchn,3-isde)=21    
        isig(nchn,3)=1  
        sigh(nchn)=facqg1   
        nchn=nchn+1 
        isig(nchn,isde)=i   
        isig(nchn,3-isde)=21    
        isig(nchn,3)=2  
        sigh(nchn)=facqg2   
  420   CONTINUE    
  430   CONTINUE    
    
      ELSEIF(isub.EQ.29) THEN   
C...f + g -> f + gamma (q + g -> q + gamma only).   
        fgq=comfac*faca*as*aem*1./3.*(sh2+uh2)/(-sh*uh) 
        DO 450 i=mina,maxa  
        IF(i.EQ.0) goto 450 
        ei=kchg(iabs(i),1)/3.   
        facgq=fgq*ei**2 
        DO 440 isde=1,2 
        IF(isde.EQ.1.AND.kfac(1,i)*kfac(2,21).EQ.0) goto 440    
        IF(isde.EQ.2.AND.kfac(1,21)*kfac(2,i).EQ.0) goto 440    
        nchn=nchn+1 
        isig(nchn,isde)=i   
        isig(nchn,3-isde)=21    
        isig(nchn,3)=1  
        sigh(nchn)=facgq    
  440   CONTINUE    
  450   CONTINUE    
    
      ELSEIF(isub.EQ.30) THEN   
C...f + g -> f + Z0 (q + g -> q + Z0 only). 
        fzq=comfac*faca*as*aem/(xw*(1.-xw))*1./48.* 
     &  (sh2+uh2+2.*sqm4*th)/(-sh*uh)   
        fzq=fzq*wids(23,2)  
        DO 470 i=mina,maxa  
        IF(i.EQ.0) goto 470 
        ei=kchg(iabs(i),1)/3.   
        ai=sign(1.,ei)  
        vi=ai-4.*ei*xw  
        faczq=fzq*(vi**2+ai**2) 
        DO 460 isde=1,2 
        IF(isde.EQ.1.AND.kfac(1,i)*kfac(2,21).EQ.0) goto 460    
        IF(isde.EQ.2.AND.kfac(1,21)*kfac(2,i).EQ.0) goto 460    
        nchn=nchn+1 
        isig(nchn,isde)=i   
        isig(nchn,3-isde)=21    
        isig(nchn,3)=1  
        sigh(nchn)=faczq    
  460   CONTINUE    
  470   CONTINUE    
      ENDIF 
    
      ELSEIF(isub.LE.40) THEN   
      IF(isub.EQ.31) THEN   
C...f + g -> f' + W+/- (q + g -> q' + W+/- only).   
        facwq=comfac*faca*as*aem/xw*1./12.* 
     &  (sh2+uh2+2.*sqm4*th)/(-sh*uh)   
        DO 490 i=mina,maxa  
        IF(i.EQ.0) goto 490 
        ia=iabs(i)  
        kchw=isign(1,kchg(ia,1)*isign(1,i)) 
        DO 480 isde=1,2 
        IF(isde.EQ.1.AND.kfac(1,i)*kfac(2,21).EQ.0) goto 480    
        IF(isde.EQ.2.AND.kfac(1,21)*kfac(2,i).EQ.0) goto 480    
        nchn=nchn+1 
        isig(nchn,isde)=i   
        isig(nchn,3-isde)=21    
        isig(nchn,3)=1  
        sigh(nchn)=facwq*vint(180+i)*wids(24,(5-kchw)/2)    
  480   CONTINUE    
  490   CONTINUE    
    
      ELSEIF(isub.EQ.32) THEN   
C...f + g -> f + H0 (q + g -> q + H0 only). 
    
      ELSEIF(isub.EQ.33) THEN   
C...f + gamma -> f + g (q + gamma -> q + g only).   
    
      ELSEIF(isub.EQ.34) THEN   
C...f + gamma -> f + gamma. 
    
      ELSEIF(isub.EQ.35) THEN   
C...f + gamma -> f + Z0.    
    
      ELSEIF(isub.EQ.36) THEN   
C...f + gamma -> f' + W+/-. 
    
      ELSEIF(isub.EQ.37) THEN   
C...f + gamma -> f + H0.    
    
      ELSEIF(isub.EQ.38) THEN   
C...f + Z0 -> f + g (q + Z0 -> q + g only). 
    
      ELSEIF(isub.EQ.39) THEN   
C...f + Z0 -> f + gamma.    
    
      ELSEIF(isub.EQ.40) THEN   
C...f + Z0 -> f + Z0.   
      ENDIF 
    
      ELSEIF(isub.LE.50) THEN   
      IF(isub.EQ.41) THEN   
C...f + Z0 -> f' + W+/-.    
    
      ELSEIF(isub.EQ.42) THEN   
C...f + Z0 -> f + H0.   
    
      ELSEIF(isub.EQ.43) THEN   
C...f + W+/- -> f' + g (q + W+/- -> q' + g only).   
    
      ELSEIF(isub.EQ.44) THEN   
C...f + W+/- -> f' + gamma. 
    
      ELSEIF(isub.EQ.45) THEN   
C...f + W+/- -> f' + Z0.    
    
      ELSEIF(isub.EQ.46) THEN   
C...f + W+/- -> f' + W+/-.  
    
      ELSEIF(isub.EQ.47) THEN   
C...f + W+/- -> f' + H0.    
    
      ELSEIF(isub.EQ.48) THEN   
C...f + H0 -> f + g (q + H0 -> q + g only). 
    
      ELSEIF(isub.EQ.49) THEN   
C...f + H0 -> f + gamma.    
    
      ELSEIF(isub.EQ.50) THEN   
C...f + H0 -> f + Z0.   
      ENDIF 
    
      ELSEIF(isub.LE.60) THEN   
      IF(isub.EQ.51) THEN   
C...f + H0 -> f' + W+/-.    
    
      ELSEIF(isub.EQ.52) THEN   
C...f + H0 -> f + H0.   
    
      ELSEIF(isub.EQ.53) THEN   
C...g + g -> f + fb (g + g -> q + qb only). 
        CALL pyhiwidt(21,sqrt(sh),wdtp,wdte)  
        facqq1=comfac*as**2*1./6.*(uh/th-(2.+mstp(34)*1./4.)*uh2/sh2)*  
     &  (wdte(0,1)+wdte(0,2)+wdte(0,3)+wdte(0,4))*faca  
        facqq2=comfac*as**2*1./6.*(th/uh-(2.+mstp(34)*1./4.)*th2/sh2)*  
     &  (wdte(0,1)+wdte(0,2)+wdte(0,3)+wdte(0,4))*faca  
        IF(kfac(1,21)*kfac(2,21).EQ.0) goto 500 
        nchn=nchn+1 
        isig(nchn,1)=21 
        isig(nchn,2)=21 
        isig(nchn,3)=1  
        sigh(nchn)=facqq1   
        nchn=nchn+1 
        isig(nchn,1)=21 
        isig(nchn,2)=21 
        isig(nchn,3)=2  
        sigh(nchn)=facqq2   
  500   CONTINUE    
    
      ELSEIF(isub.EQ.54) THEN   
C...g + gamma -> f + fb (g + gamma -> q + qb only). 
    
      ELSEIF(isub.EQ.55) THEN   
C...g + gamma -> f + fb (g + gamma -> q + qb only). 
    
      ELSEIF(isub.EQ.56) THEN   
C...g + gamma -> f + fb (g + gamma -> q + qb only). 
    
      ELSEIF(isub.EQ.57) THEN   
C...g + gamma -> f + fb (g + gamma -> q + qb only). 
    
      ELSEIF(isub.EQ.58) THEN   
C...gamma + gamma -> f + fb.    
    
      ELSEIF(isub.EQ.59) THEN   
C...gamma + Z0 -> f + fb.   
    
      ELSEIF(isub.EQ.60) THEN   
C...gamma + W+/- -> f + fb'.    
      ENDIF 
    
      ELSEIF(isub.LE.70) THEN   
      IF(isub.EQ.61) THEN   
C...gamma + H0 -> f + fb.   
    
      ELSEIF(isub.EQ.62) THEN   
C...Z0 + Z0 -> f + fb.  
    
      ELSEIF(isub.EQ.63) THEN   
C...Z0 + W+/- -> f + fb'.   
    
      ELSEIF(isub.EQ.64) THEN   
C...Z0 + H0 -> f + fb.  
    
      ELSEIF(isub.EQ.65) THEN   
C...W+ + W- -> f + fb.  
    
      ELSEIF(isub.EQ.66) THEN   
C...W+/- + H0 -> f + fb'.   
    
      ELSEIF(isub.EQ.67) THEN   
C...H0 + H0 -> f + fb.  
    
      ELSEIF(isub.EQ.68) THEN   
C...g + g -> g + g. 
        facgg1=comfac*as**2*9./4.*(sh2/th2+2.*sh/th+3.+2.*th/sh+    
     &  th2/sh2)*faca   
        facgg2=comfac*as**2*9./4.*(uh2/sh2+2.*uh/sh+3.+2.*sh/uh+    
     &  sh2/uh2)*faca   
        facgg3=comfac*as**2*9./4.*(th2/uh2+2.*th/uh+3+2.*uh/th+uh2/th2) 
        IF(kfac(1,21)*kfac(2,21).EQ.0) goto 510 
        nchn=nchn+1 
        isig(nchn,1)=21 
        isig(nchn,2)=21 
        isig(nchn,3)=1  
        sigh(nchn)=0.5*facgg1   
        nchn=nchn+1 
        isig(nchn,1)=21 
        isig(nchn,2)=21 
        isig(nchn,3)=2  
        sigh(nchn)=0.5*facgg2   
        nchn=nchn+1 
        isig(nchn,1)=21 
        isig(nchn,2)=21 
        isig(nchn,3)=3  
        sigh(nchn)=0.5*facgg3   
  510   CONTINUE    
    
      ELSEIF(isub.EQ.69) THEN   
C...gamma + gamma -> W+ + W-.   
    
      ELSEIF(isub.EQ.70) THEN   
C...gamma + W+/- -> gamma + W+/-.   
      ENDIF 
    
      ELSEIF(isub.LE.80) THEN   
      IF(isub.EQ.71) THEN   
C...Z0 + Z0 -> Z0 + Z0. 
        be2=1.-4.*sqmz/sh   
        th=-0.5*sh*be2*(1.-cth) 
        uh=-0.5*sh*be2*(1.+cth) 
        shang=1./(1.-xw)*sqmw/sqmz*(1.+be2)**2  
        ashre=(sh-sqmh)/((sh-sqmh)**2+gmmh**2)*shang    
        ashim=-gmmh/((sh-sqmh)**2+gmmh**2)*shang    
        thang=1./(1.-xw)*sqmw/sqmz*(be2-cth)**2 
        athre=(th-sqmh)/((th-sqmh)**2+gmmh**2)*thang    
        athim=-gmmh/((th-sqmh)**2+gmmh**2)*thang    
        uhang=1./(1.-xw)*sqmw/sqmz*(be2+cth)**2 
        auhre=(uh-sqmh)/((uh-sqmh)**2+gmmh**2)*uhang    
        auhim=-gmmh/((uh-sqmh)**2+gmmh**2)*uhang    
        fach=0.5*comfac*1./(4096.*paru(1)**2*16.*(1.-xw)**2)*   
     &  (aem/xw)**4*(sh/sqmw)**2*((ashre+athre+auhre)**2+   
     &  (ashim+athim+auhim)**2)*sqmz/sqmw   
        DO 530 i=min1,max1  
        IF(i.EQ.0.OR.kfac(1,i).EQ.0) goto 530   
        ei=kchg(iabs(i),1)/3.   
        ai=sign(1.,ei)  
        vi=ai-4.*ei*xw  
        avi=ai**2+vi**2 
        DO 520 j=min2,max2  
        IF(j.EQ.0.OR.kfac(2,j).EQ.0) goto 520   
        ej=kchg(iabs(j),1)/3.   
        aj=sign(1.,ej)  
        vj=aj-4.*ej*xw  
        avj=aj**2+vj**2 
        nchn=nchn+1 
        isig(nchn,1)=i  
        isig(nchn,2)=j  
        isig(nchn,3)=1  
        sigh(nchn)=fach*avi*avj 
  520   CONTINUE    
  530   CONTINUE    
    
      ELSEIF(isub.EQ.72) THEN   
C...Z0 + Z0 -> W+ + W-. 
        be2=sqrt((1.-4.*sqmw/sh)*(1.-4.*sqmz/sh))   
        cth2=cth**2 
        th=-0.5*sh*(1.-2.*(sqmw+sqmz)/sh-be2*cth)   
        uh=-0.5*sh*(1.-2.*(sqmw+sqmz)/sh+be2*cth)   
        shang=4.*sqrt(sqmw/(sqmz*(1.-xw)))*(1.-2.*sqmw/sh)* 
     &  (1.-2.*sqmz/sh) 
        ashre=(sh-sqmh)/((sh-sqmh)**2+gmmh**2)*shang    
        ashim=-gmmh/((sh-sqmh)**2+gmmh**2)*shang    
        atwre=(1.-xw)/sqmz*sh/(th-sqmw)*((cth-be2)**2*(3./2.+be2/2.*cth-    
     &  (sqmw+sqmz)/sh+(sqmw-sqmz)**2/(sh*sqmw))+4.*((sqmw+sqmz)/sh*    
     &  (1.-3.*cth2)+8.*sqmw*sqmz/sh2*(2.*cth2-1.)+ 
     &  4.*(sqmw**2+sqmz**2)/sh2*cth2+2.*(sqmw+sqmz)/sh*be2*cth))   
        atwim=0.    
        auwre=(1.-xw)/sqmz*sh/(uh-sqmw)*((cth+be2)**2*(3./2.-be2/2.*cth-    
     &  (sqmw+sqmz)/sh+(sqmw-sqmz)**2/(sh*sqmw))+4.*((sqmw+sqmz)/sh*    
     &  (1.-3.*cth2)+8.*sqmw*sqmz/sh2*(2.*cth2-1.)+ 
     &  4.*(sqmw**2+sqmz**2)/sh2*cth2-2.*(sqmw+sqmz)/sh*be2*cth))   
        auwim=0.    
        a4re=2.*(1.-xw)/sqmz*(3.-cth2-4.*(sqmw+sqmz)/sh)    
        a4im=0. 
        fach=comfac*1./(4096.*paru(1)**2*16.*(1.-xw)**2)*(aem/xw)**4*   
     &  (sh/sqmw)**2*((ashre+atwre+auwre+a4re)**2+  
     &  (ashim+atwim+auwim+a4im)**2)*sqmz/sqmw  
        DO 550 i=min1,max1  
        IF(i.EQ.0.OR.kfac(1,i).EQ.0) goto 550   
        ei=kchg(iabs(i),1)/3.   
        ai=sign(1.,ei)  
        vi=ai-4.*ei*xw  
        avi=ai**2+vi**2 
        DO 540 j=min2,max2  
        IF(j.EQ.0.OR.kfac(2,j).EQ.0) goto 540   
        ej=kchg(iabs(j),1)/3.   
        aj=sign(1.,ej)  
        vj=aj-4.*ej*xw  
        avj=aj**2+vj**2 
        nchn=nchn+1 
        isig(nchn,1)=i  
        isig(nchn,2)=j  
        isig(nchn,3)=1  
        sigh(nchn)=fach*avi*avj 
  540   CONTINUE    
  550   CONTINUE    
    
      ELSEIF(isub.EQ.73) THEN   
C...Z0 + W+/- -> Z0 + W+/-. 
        be2=1.-2.*(sqmz+sqmw)/sh+((sqmz-sqmw)/sh)**2    
        ep1=1.+(sqmz-sqmw)/sh   
        ep2=1.-(sqmz-sqmw)/sh   
        th=-0.5*sh*be2*(1.-cth) 
        uh=(sqmz-sqmw)**2/sh-0.5*sh*be2*(1.+cth)    
        thang=sqrt(sqmw/(sqmz*(1.-xw)))*(be2-ep1*cth)*(be2-ep2*cth) 
        athre=(th-sqmh)/((th-sqmh)**2+gmmh**2)*thang    
        athim=-gmmh/((th-sqmh)**2+gmmh**2)*thang    
        aswre=(1.-xw)/sqmz*sh/(sh-sqmw)*(-be2*(ep1+ep2)**4*cth+ 
     &  1./4.*(be2+ep1*ep2)**2*((ep1-ep2)**2-4.*be2*cth)+   
     &  2.*be2*(be2+ep1*ep2)*(ep1+ep2)**2*cth-  
     &  1./16.*sh/sqmw*(ep1**2-ep2**2)**2*(be2+ep1*ep2)**2) 
        aswim=0.    
        auwre=(1.-xw)/sqmz*sh/(uh-sqmw)*(-be2*(ep2+ep1*cth)*    
     &  (ep1+ep2*cth)*(be2+ep1*ep2)+be2*(ep2+ep1*cth)*  
     &  (be2+ep1*ep2*cth)*(2.*ep2-ep2*cth+ep1)-be2*(ep2+ep1*cth)**2*    
     &  (be2-ep2**2*cth)-1./8.*(be2+ep1*ep2*cth)**2*((ep1+ep2)**2+  
     &  2.*be2*(1.-cth))+1./32.*sh/sqmw*(be2+ep1*ep2*cth)**2*   
     &  (ep1**2-ep2**2)**2-be2*(ep1+ep2*cth)*(ep2+ep1*cth)* 
     &  (be2+ep1*ep2)+be2*(ep1+ep2*cth)*(be2+ep1*ep2*cth)*  
     &  (2.*ep1-ep1*cth+ep2)-be2*(ep1+ep2*cth)**2*(be2-ep1**2*cth)- 
     &  1./8.*(be2+ep1*ep2*cth)**2*((ep1+ep2)**2+2.*be2*(1.-cth))+  
     &  1./32.*sh/sqmw*(be2+ep1*ep2*cth)**2*(ep1**2-ep2**2)**2) 
        auwim=0.    
        a4re=(1.-xw)/sqmz*(ep1**2*ep2**2*(cth**2-1.)-   
     &  2.*be2*(ep1**2+ep2**2+ep1*ep2)*cth-2.*be2*ep1*ep2)  
        a4im=0. 
        fach=comfac*1./(4096.*paru(1)**2*4.*(1.-xw))*(aem/xw)**4*   
     &  (sh/sqmw)**2*((athre+aswre+auwre+a4re)**2+  
     &  (athim+aswim+auwim+a4im)**2)*sqrt(sqmz/sqmw)    
        DO 570 i=min1,max1  
        IF(i.EQ.0.OR.kfac(1,i).EQ.0) goto 570   
        ei=kchg(iabs(i),1)/3.   
        ai=sign(1.,ei)  
        vi=ai-4.*ei*xw  
        avi=ai**2+vi**2 
        DO 560 j=min2,max2  
        IF(j.EQ.0.OR.kfac(2,j).EQ.0) goto 560   
        ej=kchg(iabs(j),1)/3.   
        aj=sign(1.,ej)  
        vj=ai-4.*ej*xw  
        avj=aj**2+vj**2 
        nchn=nchn+1 
        isig(nchn,1)=i  
        isig(nchn,2)=j  
        isig(nchn,3)=1  
        sigh(nchn)=fach*(avi*vint(180+j)+vint(180+i)*avj)   
  560   CONTINUE    
  570   CONTINUE    
    
      ELSEIF(isub.EQ.75) THEN   
C...W+ + W- -> gamma + gamma.   
    
      ELSEIF(isub.EQ.76) THEN   
C...W+ + W- -> Z0 + Z0. 
        be2=sqrt((1.-4.*sqmw/sh)*(1.-4.*sqmz/sh))   
        cth2=cth**2 
        th=-0.5*sh*(1.-2.*(sqmw+sqmz)/sh-be2*cth)   
        uh=-0.5*sh*(1.-2.*(sqmw+sqmz)/sh+be2*cth)   
        shang=4.*sqrt(sqmw/(sqmz*(1.-xw)))*(1.-2.*sqmw/sh)* 
     &  (1.-2.*sqmz/sh) 
        ashre=(sh-sqmh)/((sh-sqmh)**2+gmmh**2)*shang    
        ashim=-gmmh/((sh-sqmh)**2+gmmh**2)*shang    
        atwre=(1.-xw)/sqmz*sh/(th-sqmw)*((cth-be2)**2*(3./2.+be2/2.*cth-    
     &  (sqmw+sqmz)/sh+(sqmw-sqmz)**2/(sh*sqmw))+4.*((sqmw+sqmz)/sh*    
     &  (1.-3.*cth2)+8.*sqmw*sqmz/sh2*(2.*cth2-1.)+ 
     &  4.*(sqmw**2+sqmz**2)/sh2*cth2+2.*(sqmw+sqmz)/sh*be2*cth))   
        atwim=0.    
        auwre=(1.-xw)/sqmz*sh/(uh-sqmw)*((cth+be2)**2*(3./2.-be2/2.*cth-    
     &  (sqmw+sqmz)/sh+(sqmw-sqmz)**2/(sh*sqmw))+4.*((sqmw+sqmz)/sh*    
     &  (1.-3.*cth2)+8.*sqmw*sqmz/sh2*(2.*cth2-1.)+ 
     &  4.*(sqmw**2+sqmz**2)/sh2*cth2-2.*(sqmw+sqmz)/sh*be2*cth))   
        auwim=0.    
        a4re=2.*(1.-xw)/sqmz*(3.-cth2-4.*(sqmw+sqmz)/sh)    
        a4im=0. 
        fach=0.5*comfac*1./(4096.*paru(1)**2)*(aem/xw)**4*(sh/sqmw)**2* 
     &  ((ashre+atwre+auwre+a4re)**2+(ashim+atwim+auwim+a4im)**2)   
        DO 590 i=min1,max1  
        IF(i.EQ.0.OR.kfac(1,i).EQ.0) goto 590   
        ei=sign(1.,float(i))*kchg(iabs(i),1)    
        DO 580 j=min2,max2  
        IF(j.EQ.0.OR.kfac(2,j).EQ.0) goto 580   
        ej=sign(1.,float(j))*kchg(iabs(j),1)    
        IF(ei*ej.GT.0.) goto 580    
        nchn=nchn+1 
        isig(nchn,1)=i  
        isig(nchn,2)=j  
        isig(nchn,3)=1  
        sigh(nchn)=fach*vint(180+i)*vint(180+j) 
  580   CONTINUE    
  590   CONTINUE    
    
      ELSEIF(isub.EQ.77) THEN   
C...W+/- + W+/- -> W+/- + W+/-. 
        be2=1.-4.*sqmw/sh   
        be4=be2**2  
        cth2=cth**2 
        cth3=cth**3 
        th=-0.5*sh*be2*(1.-cth) 
        uh=-0.5*sh*be2*(1.+cth) 
        shang=(1.+be2)**2   
        ashre=(sh-sqmh)/((sh-sqmh)**2+gmmh**2)*shang    
        ashim=-gmmh/((sh-sqmh)**2+gmmh**2)*shang    
        thang=(be2-cth)**2  
        athre=(th-sqmh)/((th-sqmh)**2+gmmh**2)*thang    
        athim=-gmmh/((th-sqmh)**2+gmmh**2)*thang    
        sgzang=1./sqmw*be2*(3.-be2)**2*cth  
        asgre=xw*sgzang 
        asgim=0.    
        aszre=(1.-xw)*sh/(sh-sqmz)*sgzang   
        aszim=0.    
        tgzang=1./sqmw*(be2*(4.-2.*be2+be4)+be2*(4.-10.*be2+be4)*cth+   
     &  (2.-11.*be2+10.*be4)*cth2+be2*cth3) 
        atgre=0.5*xw*sh/th*tgzang   
        atgim=0.    
        atzre=0.5*(1.-xw)*sh/(th-sqmz)*tgzang   
        atzim=0.    
        a4re=1./sqmw*(1.+2.*be2-6.*be2*cth-cth2)    
        a4im=0. 
        fach=comfac*1./(4096.*paru(1)**2)*(aem/xw)**4*(sh/sqmw)**2* 
     &  ((ashre+athre+asgre+aszre+atgre+atzre+a4re)**2+ 
     &  (ashim+athim+asgim+aszim+atgim+atzim+a4im)**2)  
        DO 610 i=min1,max1  
        IF(i.EQ.0.OR.kfac(1,i).EQ.0) goto 610   
        ei=sign(1.,float(i))*kchg(iabs(i),1)    
        DO 600 j=min2,max2  
        IF(j.EQ.0.OR.kfac(2,j).EQ.0) goto 600   
        ej=sign(1.,float(j))*kchg(iabs(j),1)    
        IF(ei*ej.GT.0.) goto 600    
        nchn=nchn+1 
        isig(nchn,1)=i  
        isig(nchn,2)=j  
        isig(nchn,3)=1  
        sigh(nchn)=fach*vint(180+i)*vint(180+j) 
  600   CONTINUE    
  610   CONTINUE    
    
      ELSEIF(isub.EQ.78) THEN   
C...W+/- + H0 -> W+/- + H0. 
    
      ELSEIF(isub.EQ.79) THEN   
C...H0 + H0 -> H0 + H0. 
    
      ENDIF 
    
C...C: 2 -> 2, tree diagrams with masses.   
    
      ELSEIF(isub.LE.90) THEN   
      IF(isub.EQ.81) THEN   
C...q + qb -> Q + QB.   
        facqqb=comfac*as**2*4./9.*(((th-sqm3)**2+   
     &  (uh-sqm3)**2)/sh2+2.*sqm3/sh)   
        IF(mstp(35).GE.1) THEN  
          IF(mstp(35).EQ.1) THEN    
            alssg=parp(35)  
          ELSE  
            mst115=mstu(115)    
            mstu(115)=mstp(36)  
            q2bn=sqrt(sqm3*((sqrt(sh)-2.*sqrt(sqm3))**2+parp(36)**2))   
            alssg=ulalps(q2bn)  
            mstu(115)=mst115    
          ENDIF 
          xrepu=paru(1)*alssg/(6.*sqrt(max(1e-20,1.-4.*sqm3/sh)))   
          frepu=xrepu/(exp(min(100.,xrepu))-1.) 
          pari(81)=frepu    
          facqqb=facqqb*frepu   
        ENDIF   
        DO 620 i=mina,maxa  
        IF(i.EQ.0.OR.kfac(1,i)*kfac(2,-i).EQ.0) goto 620    
        nchn=nchn+1 
        isig(nchn,1)=i  
        isig(nchn,2)=-i 
        isig(nchn,3)=1  
        sigh(nchn)=facqqb   
  620   CONTINUE    
    
      ELSEIF(isub.EQ.82) THEN   
C...g + g -> Q + QB.    
        facqq1=comfac*faca*as**2*1./6.*((uh-sqm3)/(th-sqm3)-    
     &  2.*(uh-sqm3)**2/sh2+4.*sqm3/sh*(th*uh-sqm3**2)/(th-sqm3)**2)    
        facqq2=comfac*faca*as**2*1./6.*((th-sqm3)/(uh-sqm3)-    
     &  2.*(th-sqm3)**2/sh2+4.*sqm3/sh*(th*uh-sqm3**2)/(uh-sqm3)**2)    
        IF(mstp(35).GE.1) THEN  
          IF(mstp(35).EQ.1) THEN    
            alssg=parp(35)  
          ELSE  
            mst115=mstu(115)    
            mstu(115)=mstp(36)  
            q2bn=sqrt(sqm3*((sqrt(sh)-2.*sqrt(sqm3))**2+parp(36)**2))   
            alssg=ulalps(q2bn)  
            mstu(115)=mst115    
          ENDIF 
          xattr=4.*paru(1)*alssg/(3.*sqrt(max(1e-20,1.-4.*sqm3/sh)))    
          fattr=xattr/(1.-exp(-min(100.,xattr)))    
          xrepu=paru(1)*alssg/(6.*sqrt(max(1e-20,1.-4.*sqm3/sh)))   
          frepu=xrepu/(exp(min(100.,xrepu))-1.) 
          fatre=(2.*fattr+5.*frepu)/7.  
          pari(81)=fatre    
          facqq1=facqq1*fatre   
          facqq2=facqq2*fatre   
        ENDIF   
        IF(kfac(1,21)*kfac(2,21).EQ.0) goto 630 
        nchn=nchn+1 
        isig(nchn,1)=21 
        isig(nchn,2)=21 
        isig(nchn,3)=1  
        sigh(nchn)=facqq1   
        nchn=nchn+1 
        isig(nchn,1)=21 
        isig(nchn,2)=21 
        isig(nchn,3)=2  
        sigh(nchn)=facqq2   
  630   CONTINUE    
    
      ENDIF 
    
C...D: Mimimum bias processes.  
    
      ELSEIF(isub.LE.100) THEN  
      IF(isub.EQ.91) THEN   
C...Elastic scattering. 
        sigs=xsec(isub,1)   
    
      ELSEIF(isub.EQ.92) THEN   
C...Single diffractive scattering.  
        sigs=xsec(isub,1)   
    
      ELSEIF(isub.EQ.93) THEN   
C...Double diffractive scattering.  
        sigs=xsec(isub,1)   
    
      ELSEIF(isub.EQ.94) THEN   
C...Central diffractive scattering. 
        sigs=xsec(isub,1)   
    
      ELSEIF(isub.EQ.95) THEN   
C...Low-pT scattering.  
        sigs=xsec(isub,1)   
    
      ELSEIF(isub.EQ.96) THEN   
C...Multiple interactions: sum of QCD processes.    
        CALL pyhiwidt(21,sqrt(sh),wdtp,wdte)  
    
C...q + q' -> q + q'.   
        facqq1=comfac*as**2*4./9.*(sh2+uh2)/th2 
        facqqb=comfac*as**2*4./9.*((sh2+uh2)/th2*faca-  
     &  mstp(34)*2./3.*uh2/(sh*th)) 
        facqq2=comfac*as**2*4./9.*((sh2+th2)/uh2-   
     &  mstp(34)*2./3.*sh2/(th*uh)) 
        DO 650 i=-3,3   
        IF(i.EQ.0) goto 650 
        DO 640 j=-3,3   
        IF(j.EQ.0) goto 640 
        nchn=nchn+1 
        isig(nchn,1)=i  
        isig(nchn,2)=j  
        isig(nchn,3)=111    
        sigh(nchn)=facqq1   
        IF(i.EQ.-j) sigh(nchn)=facqqb   
        IF(i.EQ.j) THEN 
          sigh(nchn)=0.5*sigh(nchn) 
          nchn=nchn+1   
          isig(nchn,1)=i    
          isig(nchn,2)=j    
          isig(nchn,3)=112  
          sigh(nchn)=0.5*facqq2 
        ENDIF   
  640   CONTINUE    
  650   CONTINUE    
    
C...q + qb -> q' + qb' or g + g.    
        facqqb=comfac*as**2*4./9.*(th2+uh2)/sh2*(wdte(0,1)+wdte(0,2)+   
     &  wdte(0,3)+wdte(0,4))    
        facgg1=comfac*as**2*32./27.*(uh/th-(2.+mstp(34)*1./4.)*uh2/sh2) 
        facgg2=comfac*as**2*32./27.*(th/uh-(2.+mstp(34)*1./4.)*th2/sh2) 
        DO 660 i=-3,3   
        IF(i.EQ.0) goto 660 
        nchn=nchn+1 
        isig(nchn,1)=i  
        isig(nchn,2)=-i 
        isig(nchn,3)=121    
        sigh(nchn)=facqqb   
        nchn=nchn+1 
        isig(nchn,1)=i  
        isig(nchn,2)=-i 
        isig(nchn,3)=131    
        sigh(nchn)=0.5*facgg1   
        nchn=nchn+1 
        isig(nchn,1)=i  
        isig(nchn,2)=-i 
        isig(nchn,3)=132    
        sigh(nchn)=0.5*facgg2   
  660   CONTINUE    
    
C...q + g -> q + g. 
        facqg1=comfac*as**2*4./9.*((2.+mstp(34)*1./4.)*uh2/th2-uh/sh)*  
     &  faca    
        facqg2=comfac*as**2*4./9.*((2.+mstp(34)*1./4.)*sh2/th2-sh/uh)   
        DO 680 i=-3,3   
        IF(i.EQ.0) goto 680 
        DO 670 isde=1,2 
        nchn=nchn+1 
        isig(nchn,isde)=i   
        isig(nchn,3-isde)=21    
        isig(nchn,3)=281    
        sigh(nchn)=facqg1   
        nchn=nchn+1 
        isig(nchn,isde)=i   
        isig(nchn,3-isde)=21    
        isig(nchn,3)=282    
        sigh(nchn)=facqg2   
  670   CONTINUE    
  680   CONTINUE    
    
C...g + g -> q + qb or g + g.   
        facqq1=comfac*as**2*1./6.*(uh/th-(2.+mstp(34)*1./4.)*uh2/sh2)*  
     &  (wdte(0,1)+wdte(0,2)+wdte(0,3)+wdte(0,4))*faca  
        facqq2=comfac*as**2*1./6.*(th/uh-(2.+mstp(34)*1./4.)*th2/sh2)*  
     &  (wdte(0,1)+wdte(0,2)+wdte(0,3)+wdte(0,4))*faca  
        facgg1=comfac*as**2*9./4.*(sh2/th2+2.*sh/th+3.+2.*th/sh+    
     &  th2/sh2)*faca   
        facgg2=comfac*as**2*9./4.*(uh2/sh2+2.*uh/sh+3.+2.*sh/uh+    
     &  sh2/uh2)*faca   
        facgg3=comfac*as**2*9./4.*(th2/uh2+2.*th/uh+3+2.*uh/th+uh2/th2) 
        nchn=nchn+1 
        isig(nchn,1)=21 
        isig(nchn,2)=21 
        isig(nchn,3)=531    
        sigh(nchn)=facqq1   
        nchn=nchn+1 
        isig(nchn,1)=21 
        isig(nchn,2)=21 
        isig(nchn,3)=532    
        sigh(nchn)=facqq2   
        nchn=nchn+1 
        isig(nchn,1)=21 
        isig(nchn,2)=21 
        isig(nchn,3)=681    
        sigh(nchn)=0.5*facgg1   
        nchn=nchn+1 
        isig(nchn,1)=21 
        isig(nchn,2)=21 
        isig(nchn,3)=682    
        sigh(nchn)=0.5*facgg2   
        nchn=nchn+1 
        isig(nchn,1)=21 
        isig(nchn,2)=21 
        isig(nchn,3)=683    
        sigh(nchn)=0.5*facgg3   
      ENDIF 
    
C...E: 2 -> 1, loop diagrams.   
    
      ELSEIF(isub.LE.110) THEN  
      IF(isub.EQ.101) THEN  
C...g + g -> gamma*/Z0. 
    
      ELSEIF(isub.EQ.102) THEN  
C...g + g -> H0.    
        CALL pyhiwidt(25,sqrt(sh),wdtp,wdte)  
        etare=0.    
        etaim=0.    
        DO 690 i=1,2*mstp(1)    
        eps=4.*pmas(i,1)**2/sh  
        IF(eps.LE.1.) THEN  
          IF(eps.GT.1.e-4) THEN 
            root=sqrt(1.-eps)   
            rln=log((1.+root)/(1.-root))    
          ELSE  
            rln=log(4./eps-2.)  
          ENDIF 
          phire=0.25*(rln**2-paru(1)**2)    
          phiim=0.5*paru(1)*rln 
        ELSE    
          phire=-(asin(1./sqrt(eps)))**2    
          phiim=0.  
        ENDIF   
        etare=etare+0.5*eps*(1.+(eps-1.)*phire) 
        etaim=etaim+0.5*eps*(eps-1.)*phiim  
  690   CONTINUE    
        eta2=etare**2+etaim**2  
        fach=comfac*faca*(as/paru(1)*aem/xw)**2*1./512.*    
     &  (sh/sqmw)**2*eta2*sh2/((sh-sqmh)**2+gmmh**2)*   
     &  (wdte(0,1)+wdte(0,2)+wdte(0,4)) 
        IF(kfac(1,21)*kfac(2,21).EQ.0) goto 700 
        nchn=nchn+1 
        isig(nchn,1)=21 
        isig(nchn,2)=21 
        isig(nchn,3)=1  
        sigh(nchn)=fach 
  700   CONTINUE    
    
      ENDIF 
    
C...F: 2 -> 2, box diagrams.    
    
      ELSEIF(isub.LE.120) THEN  
      IF(isub.EQ.111) THEN  
C...f + fb -> g + H0 (q + qb -> g + H0 only).   
        a5stur=0.   
        a5stui=0.   
        DO 710 i=1,2*mstp(1)    
        sqmq=pmas(i,1)**2   
        epss=4.*sqmq/sh 
        epsh=4.*sqmq/sqmh   
        a5stur=a5stur+sqmq/sqmh*(4.+4.*sh/(th+uh)*(pyhiw1au(epss,1)-  
     &  pyhiw1au(epsh,1))+(1.-4.*sqmq/(th+uh))*(pyhiw2au(epss,1)-   
     &  pyhiw2au(epsh,1)))    
        a5stui=a5stui+sqmq/sqmh*(4.*sh/(th+uh)*(pyhiw1au(epss,2)- 
     &  pyhiw1au(epsh,2))+(1.-4.*sqmq/(th+uh))*(pyhiw2au(epss,2)-   
     &  pyhiw2au(epsh,2)))    
  710   CONTINUE    
        facgh=comfac*faca/(144.*paru(1)**2)*aem/xw*as**3*sqmh/sqmw* 
     &  sqmh/sh*(uh**2+th**2)/(uh+th)**2*(a5stur**2+a5stui**2)  
        facgh=facgh*wids(25,2)  
        DO 720 i=mina,maxa  
        IF(i.EQ.0.OR.kfac(1,i)*kfac(2,-i).EQ.0) goto 720    
        nchn=nchn+1 
        isig(nchn,1)=i  
        isig(nchn,2)=-i 
        isig(nchn,3)=1  
        sigh(nchn)=facgh    
  720   CONTINUE    
    
      ELSEIF(isub.EQ.112) THEN  
C...f + g -> f + H0 (q + g -> q + H0 only). 
        a5tsur=0.   
        a5tsui=0.   
        DO 730 i=1,2*mstp(1)    
        sqmq=pmas(i,1)**2   
        epst=4.*sqmq/th 
        epsh=4.*sqmq/sqmh   
        a5tsur=a5tsur+sqmq/sqmh*(4.+4.*th/(sh+uh)*(pyhiw1au(epst,1)-  
     &  pyhiw1au(epsh,1))+(1.-4.*sqmq/(sh+uh))*(pyhiw2au(epst,1)-   
     &  pyhiw2au(epsh,1)))    
        a5tsui=a5tsui+sqmq/sqmh*(4.*th/(sh+uh)*(pyhiw1au(epst,2)- 
     &  pyhiw1au(epsh,2))+(1.-4.*sqmq/(sh+uh))*(pyhiw2au(epst,2)-   
     &  pyhiw2au(epsh,2)))    
  730   CONTINUE    
        facqh=comfac*faca/(384.*paru(1)**2)*aem/xw*as**3*sqmh/sqmw* 
     &  sqmh/(-th)*(uh**2+sh**2)/(uh+sh)**2*(a5tsur**2+a5tsui**2)   
        facqh=facqh*wids(25,2)  
        DO 750 i=mina,maxa  
        IF(i.EQ.0) goto 750 
        DO 740 isde=1,2 
        IF(isde.EQ.1.AND.kfac(1,i)*kfac(2,21).EQ.0) goto 740    
        IF(isde.EQ.2.AND.kfac(1,21)*kfac(2,i).EQ.0) goto 740    
        nchn=nchn+1 
        isig(nchn,isde)=i   
        isig(nchn,3-isde)=21    
        isig(nchn,3)=1  
        sigh(nchn)=facqh    
  740   CONTINUE    
  750   CONTINUE    
    
      ELSEIF(isub.EQ.113) THEN  
C...g + g -> g + H0.    
        a2stur=0.   
        a2stui=0.   
        a2ustr=0.   
        a2usti=0.   
        a2tusr=0.   
        a2tusi=0.   
        a4stur=0.   
        a4stui=0.   
        DO 760 i=6,2*mstp(1)    
C'''Only t-quarks yet included  
        sqmq=pmas(i,1)**2   
        epss=4.*sqmq/sh 
        epst=4.*sqmq/th 
        epsu=4.*sqmq/uh 
        epsh=4.*sqmq/sqmh   
        IF(epsh.LT.1.e-6) goto 760  
        bestu=0.5*(1.+sqrt(1.+epss*th/uh))  
        beust=0.5*(1.+sqrt(1.+epsu*sh/th))  
        betus=0.5*(1.+sqrt(1.+epst*uh/sh))  
        beuts=bestu 
        betsu=beust 
        besut=betus 
        w3stur=pyhii3au(bestu,epsh,1)-pyhii3au(bestu,epss,1)-   
     &  pyhii3au(bestu,epsu,1)    
        w3stui=pyhii3au(bestu,epsh,2)-pyhii3au(bestu,epss,2)-   
     &  pyhii3au(bestu,epsu,2)    
        w3sutr=pyhii3au(besut,epsh,1)-pyhii3au(besut,epss,1)-   
     &  pyhii3au(besut,epst,1)    
        w3suti=pyhii3au(besut,epsh,2)-pyhii3au(besut,epss,2)-   
     &  pyhii3au(besut,epst,2)    
        w3tsur=pyhii3au(betsu,epsh,1)-pyhii3au(betsu,epst,1)-   
     &  pyhii3au(betsu,epsu,1)    
        w3tsui=pyhii3au(betsu,epsh,2)-pyhii3au(betsu,epst,2)-   
     &  pyhii3au(betsu,epsu,2)    
        w3tusr=pyhii3au(betus,epsh,1)-pyhii3au(betus,epst,1)-   
     &  pyhii3au(betus,epss,1)    
        w3tusi=pyhii3au(betus,epsh,2)-pyhii3au(betus,epst,2)-   
     &  pyhii3au(betus,epss,2)    
        w3ustr=pyhii3au(beust,epsh,1)-pyhii3au(beust,epsu,1)-   
     &  pyhii3au(beust,epst,1)    
        w3usti=pyhii3au(beust,epsh,2)-pyhii3au(beust,epsu,2)-   
     &  pyhii3au(beust,epst,2)    
        w3utsr=pyhii3au(beuts,epsh,1)-pyhii3au(beuts,epsu,1)-   
     &  pyhii3au(beuts,epss,1)    
        w3utsi=pyhii3au(beuts,epsh,2)-pyhii3au(beuts,epsu,2)-   
     &  pyhii3au(beuts,epss,2)    
        b2stur=sqmq/sqmh**2*(sh*(uh-sh)/(sh+uh)+2.*th*uh*(uh+2.*sh)/    
     &  (sh+uh)**2*(pyhiw1au(epst,1)-pyhiw1au(epsh,1))+(sqmq-sh/4.)*    
     &  (0.5*pyhiw2au(epss,1)+0.5*pyhiw2au(epsh,1)-pyhiw2au(epst,1)+
     &  w3stur)+  
     &  sh**2*(2.*sqmq/(sh+uh)**2-0.5/(sh+uh))*(pyhiw2au(epst,1)- 
     &  pyhiw2au(epsh,1))+0.5*th*uh/sh*(pyhiw2au(epsh,1)-
     &  2.*pyhiw2au(epst,1))+    
     &  0.125*(sh-12.*sqmq-4.*th*uh/sh)*w3tsur) 
        b2stui=sqmq/sqmh**2*(2.*th*uh*(uh+2.*sh)/(sh+uh)**2*    
     &  (pyhiw1au(epst,2)-pyhiw1au(epsh,2))+(sqmq-sh/4.)*   
     &  (0.5*pyhiw2au(epss,2)+0.5*pyhiw2au(epsh,2)-pyhiw2au(epst,2)+
     &  w3stui)+  
     &  sh**2*(2.*sqmq/(sh+uh)**2-0.5/(sh+uh))*(pyhiw2au(epst,2)- 
     &  pyhiw2au(epsh,2))+0.5*th*uh/sh*(pyhiw2au(epsh,2)-
     &  2.*pyhiw2au(epst,2))+    
     &  0.125*(sh-12.*sqmq-4.*th*uh/sh)*w3tsui) 
        b2sutr=sqmq/sqmh**2*(sh*(th-sh)/(sh+th)+2.*uh*th*(th+2.*sh)/    
     &  (sh+th)**2*(pyhiw1au(epsu,1)-pyhiw1au(epsh,1))+(sqmq-sh/4.)*    
     &  (0.5*pyhiw2au(epss,1)+0.5*pyhiw2au(epsh,1)-pyhiw2au(epsu,1)+
     &  w3sutr)+  
     &  sh**2*(2.*sqmq/(sh+th)**2-0.5/(sh+th))*(pyhiw2au(epsu,1)- 
     &  pyhiw2au(epsh,1))+0.5*uh*th/sh*(pyhiw2au(epsh,1)-
     &  2.*pyhiw2au(epsu,1))+    
     &  0.125*(sh-12.*sqmq-4.*uh*th/sh)*w3ustr) 
        b2suti=sqmq/sqmh**2*(2.*uh*th*(th+2.*sh)/(sh+th)**2*    
     &  (pyhiw1au(epsu,2)-pyhiw1au(epsh,2))+(sqmq-sh/4.)*   
     &  (0.5*pyhiw2au(epss,2)+0.5*pyhiw2au(epsh,2)-pyhiw2au(epsu,2)+
     &  w3suti)+  
     &  sh**2*(2.*sqmq/(sh+th)**2-0.5/(sh+th))*(pyhiw2au(epsu,2)- 
     &  pyhiw2au(epsh,2))+0.5*uh*th/sh*(pyhiw2au(epsh,2)-
     &  2.*pyhiw2au(epsu,2))+    
     &  0.125*(sh-12.*sqmq-4.*uh*th/sh)*w3usti) 
        b2tsur=sqmq/sqmh**2*(th*(uh-th)/(th+uh)+2.*sh*uh*(uh+2.*th)/    
     &  (th+uh)**2*(pyhiw1au(epss,1)-pyhiw1au(epsh,1))+(sqmq-th/4.)*    
     &  (0.5*pyhiw2au(epst,1)+0.5*pyhiw2au(epsh,1)-pyhiw2au(epss,1)+
     &  w3tsur)+  
     &  th**2*(2.*sqmq/(th+uh)**2-0.5/(th+uh))*(pyhiw2au(epss,1)- 
     &  pyhiw2au(epsh,1))+0.5*sh*uh/th*(pyhiw2au(epsh,1)-
     &  2.*pyhiw2au(epss,1))+    
     &  0.125*(th-12.*sqmq-4.*sh*uh/th)*w3stur) 
        b2tsui=sqmq/sqmh**2*(2.*sh*uh*(uh+2.*th)/(th+uh)**2*    
     &  (pyhiw1au(epss,2)-pyhiw1au(epsh,2))+(sqmq-th/4.)*   
     &  (0.5*pyhiw2au(epst,2)+0.5*pyhiw2au(epsh,2)-pyhiw2au(epss,2)+
     &  w3tsui)+  
     &  th**2*(2.*sqmq/(th+uh)**2-0.5/(th+uh))*(pyhiw2au(epss,2)- 
     &  pyhiw2au(epsh,2))+0.5*sh*uh/th*(pyhiw2au(epsh,2)-
     &  2.*pyhiw2au(epss,2))+    
     &  0.125*(th-12.*sqmq-4.*sh*uh/th)*w3stui) 
        b2tusr=sqmq/sqmh**2*(th*(sh-th)/(th+sh)+2.*uh*sh*(sh+2.*th)/    
     &  (th+sh)**2*(pyhiw1au(epsu,1)-pyhiw1au(epsh,1))+(sqmq-th/4.)*    
     &  (0.5*pyhiw2au(epst,1)+0.5*pyhiw2au(epsh,1)-pyhiw2au(epsu,1)+
     &  w3tusr)+  
     &  th**2*(2.*sqmq/(th+sh)**2-0.5/(th+sh))*(pyhiw2au(epsu,1)- 
     &  pyhiw2au(epsh,1))+0.5*uh*sh/th*(pyhiw2au(epsh,1)-
     &  2.*pyhiw2au(epsu,1))+    
     &  0.125*(th-12.*sqmq-4.*uh*sh/th)*w3utsr) 
        b2tusi=sqmq/sqmh**2*(2.*uh*sh*(sh+2.*th)/(th+sh)**2*    
     &  (pyhiw1au(epsu,2)-pyhiw1au(epsh,2))+(sqmq-th/4.)*   
     &  (0.5*pyhiw2au(epst,2)+0.5*pyhiw2au(epsh,2)-pyhiw2au(epsu,2)+
     &  w3tusi)+  
     &  th**2*(2.*sqmq/(th+sh)**2-0.5/(th+sh))*(pyhiw2au(epsu,2)- 
     &  pyhiw2au(epsh,2))+0.5*uh*sh/th*(pyhiw2au(epsh,2)-
     &  2.*pyhiw2au(epsu,2))+    
     &  0.125*(th-12.*sqmq-4.*uh*sh/th)*w3utsi) 
        b2ustr=sqmq/sqmh**2*(uh*(th-uh)/(uh+th)+2.*sh*th*(th+2.*uh)/    
     &  (uh+th)**2*(pyhiw1au(epss,1)-pyhiw1au(epsh,1))+(sqmq-uh/4.)*    
     &  (0.5*pyhiw2au(epsu,1)+0.5*pyhiw2au(epsh,1)-pyhiw2au(epss,1)+
     &  w3ustr)+  
     &  uh**2*(2.*sqmq/(uh+th)**2-0.5/(uh+th))*(pyhiw2au(epss,1)- 
     &  pyhiw2au(epsh,1))+0.5*sh*th/uh*(pyhiw2au(epsh,1)-
     &  2.*pyhiw2au(epss,1))+    
     &  0.125*(uh-12.*sqmq-4.*sh*th/uh)*w3sutr) 
        b2usti=sqmq/sqmh**2*(2.*sh*th*(th+2.*uh)/(uh+th)**2*    
     &  (pyhiw1au(epss,2)-pyhiw1au(epsh,2))+(sqmq-uh/4.)*   
     &  (0.5*pyhiw2au(epsu,2)+0.5*pyhiw2au(epsh,2)-pyhiw2au(epss,2)+
     &  w3usti)+  
     &  uh**2*(2.*sqmq/(uh+th)**2-0.5/(uh+th))*(pyhiw2au(epss,2)- 
     &  pyhiw2au(epsh,2))+0.5*sh*th/uh*(pyhiw2au(epsh,2)-
     &  2.*pyhiw2au(epss,2))+    
     &  0.125*(uh-12.*sqmq-4.*sh*th/uh)*w3suti) 
        b2utsr=sqmq/sqmh**2*(uh*(sh-uh)/(uh+sh)+2.*th*sh*(sh+2.*uh)/    
     &  (uh+sh)**2*(pyhiw1au(epst,1)-pyhiw1au(epsh,1))+(sqmq-uh/4.)*    
     &  (0.5*pyhiw2au(epsu,1)+0.5*pyhiw2au(epsh,1)-pyhiw2au(epst,1)+
     &  w3utsr)+  
     &  uh**2*(2.*sqmq/(uh+sh)**2-0.5/(uh+sh))*(pyhiw2au(epst,1)- 
     &  pyhiw2au(epsh,1))+0.5*th*sh/uh*(pyhiw2au(epsh,1)-
     &  2.*pyhiw2au(epst,1))+    
     &  0.125*(uh-12.*sqmq-4.*th*sh/uh)*w3tusr) 
        b2utsi=sqmq/sqmh**2*(2.*th*sh*(sh+2.*uh)/(uh+sh)**2*    
     &  (pyhiw1au(epst,2)-pyhiw1au(epsh,2))+(sqmq-uh/4.)*   
     &  (0.5*pyhiw2au(epsu,2)+0.5*pyhiw2au(epsh,2)-pyhiw2au(epst,2)+
     &  w3utsi)+  
     &  uh**2*(2.*sqmq/(uh+sh)**2-0.5/(uh+sh))*(pyhiw2au(epst,2)- 
     &  pyhiw2au(epsh,2))+0.5*th*sh/uh*(pyhiw2au(epsh,2)-
     &  2.*pyhiw2au(epst,2))+    
     &  0.125*(uh-12.*sqmq-4.*th*sh/uh)*w3tusi) 
        b4stur=sqmq/sqmh*(-2./3.+(sqmq/sqmh-1./4.)*(pyhiw2au(epss,1)- 
     &  pyhiw2au(epsh,1)+w3stur)) 
        b4stui=sqmq/sqmh*(sqmq/sqmh-1./4.)*(pyhiw2au(epss,2)- 
     &  pyhiw2au(epsh,2)+w3stui)  
        b4tusr=sqmq/sqmh*(-2./3.+(sqmq/sqmh-1./4.)*(pyhiw2au(epst,1)- 
     &  pyhiw2au(epsh,1)+w3tusr)) 
        b4tusi=sqmq/sqmh*(sqmq/sqmh-1./4.)*(pyhiw2au(epst,2)- 
     &  pyhiw2au(epsh,2)+w3tusi)  
        b4ustr=sqmq/sqmh*(-2./3.+(sqmq/sqmh-1./4.)*(pyhiw2au(epsu,1)- 
     &  pyhiw2au(epsh,1)+w3ustr)) 
        b4usti=sqmq/sqmh*(sqmq/sqmh-1./4.)*(pyhiw2au(epsu,2)- 
     &  pyhiw2au(epsh,2)+w3usti)  
        a2stur=a2stur+b2stur+b2sutr 
        a2stui=a2stui+b2stui+b2suti 
        a2ustr=a2ustr+b2ustr+b2utsr 
        a2usti=a2usti+b2usti+b2utsi 
        a2tusr=a2tusr+b2tusr+b2tsur 
        a2tusi=a2tusi+b2tusi+b2tsui 
        a4stur=a4stur+b4stur+b4ustr+b4tusr  
        a4stui=a4stui+b4stui+b4usti+b4tusi  
  760   CONTINUE    
        facgh=comfac*faca*3./(128.*paru(1)**2)*aem/xw*as**3*    
     &  sqmh/sqmw*sqmh**3/(sh*th*uh)*(a2stur**2+a2stui**2+a2ustr**2+    
     &  a2usti**2+a2tusr**2+a2tusi**2+a4stur**2+a4stui**2)  
        facgh=facgh*wids(25,2)  
        IF(kfac(1,21)*kfac(2,21).EQ.0) goto 770 
        nchn=nchn+1 
        isig(nchn,1)=21 
        isig(nchn,2)=21 
        isig(nchn,3)=1  
        sigh(nchn)=facgh    
  770   CONTINUE    
    
      ELSEIF(isub.EQ.114) THEN  
C...g + g -> gamma + gamma. 
        asre=0. 
        asim=0. 
        DO 780 i=1,2*mstp(1)    
        ei=kchg(iabs(i),1)/3.   
        sqmq=pmas(i,1)**2   
        epss=4.*sqmq/sh 
        epst=4.*sqmq/th 
        epsu=4.*sqmq/uh 
        IF(epss+abs(epst)+abs(epsu).LT.3.e-6) THEN  
          a0stur=1.+(th-uh)/sh*log(th/uh)+0.5*(th2+uh2)/sh2*    
     &    (log(th/uh)**2+paru(1)**2)    
          a0stui=0. 
          a0tsur=1.+(sh-uh)/th*log(-sh/uh)+0.5*(sh2+uh2)/th2*   
     &    log(-sh/uh)**2    
          a0tsui=-paru(1)*((sh-uh)/th+(sh2+uh2)/th2*log(-sh/uh))    
          a0utsr=1.+(th-sh)/uh*log(-th/sh)+0.5*(th2+sh2)/uh2*   
     &    log(-th/sh)**2    
          a0utsi=paru(1)*((th-sh)/uh+(th2+sh2)/uh2*log(-th/sh)) 
          a1stur=-1.    
          a1stui=0. 
          a2stur=-1.    
          a2stui=0. 
        ELSE    
          bestu=0.5*(1.+sqrt(1.+epss*th/uh))    
          beust=0.5*(1.+sqrt(1.+epsu*sh/th))    
          betus=0.5*(1.+sqrt(1.+epst*uh/sh))    
          beuts=bestu   
          betsu=beust   
          besut=betus   
          a0stur=1.+(1.+2.*th/sh)*pyhiw1au(epst,1)+(1.+2.*uh/sh)* 
     &    pyhiw1au(epsu,1)+0.5*((th2+uh2)/sh2-epss)*(pyhiw2au(epst,1)+  
     &    pyhiw2au(epsu,1))-
     &    0.25*epst*(1.-0.5*epss)*(pyhii3au(besut,epss,1)+    
     &    pyhii3au(besut,epst,1))-0.25*epsu*(1.-0.5*epss)*    
     &    (pyhii3au(bestu,epss,1)+pyhii3au(bestu,epsu,1))+  
     &    0.25*(-2.*(th2+uh2)/sh2+4.*epss+epst+epsu+0.5*epst*epsu)* 
     &    (pyhii3au(betsu,epst,1)+pyhii3au(betsu,epsu,1))   
          a0stui=(1.+2.*th/sh)*pyhiw1au(epst,2)+(1.+2.*uh/sh)*    
     &    pyhiw1au(epsu,2)+0.5*((th2+uh2)/sh2-epss)*(pyhiw2au(epst,2)+  
     &    pyhiw2au(epsu,2))-
     &    0.25*epst*(1.-0.5*epss)*(pyhii3au(besut,epss,2)+    
     &    pyhii3au(besut,epst,2))-0.25*epsu*(1.-0.5*epss)*    
     &    (pyhii3au(bestu,epss,2)+pyhii3au(bestu,epsu,2))+  
     &    0.25*(-2.*(th2+uh2)/sh2+4.*epss+epst+epsu+0.5*epst*epsu)* 
     &    (pyhii3au(betsu,epst,2)+pyhii3au(betsu,epsu,2))   
          a0tsur=1.+(1.+2.*sh/th)*pyhiw1au(epss,1)+(1.+2.*uh/th)* 
     &    pyhiw1au(epsu,1)+0.5*((sh2+uh2)/th2-epst)*(pyhiw2au(epss,1)+  
     &    pyhiw2au(epsu,1))-
     &    0.25*epss*(1.-0.5*epst)*(pyhii3au(betus,epst,1)+    
     &    pyhii3au(betus,epss,1))-0.25*epsu*(1.-0.5*epst)*    
     &    (pyhii3au(betsu,epst,1)+pyhii3au(betsu,epsu,1))+  
     &    0.25*(-2.*(sh2+uh2)/th2+4.*epst+epss+epsu+0.5*epss*epsu)* 
     &    (pyhii3au(bestu,epss,1)+pyhii3au(bestu,epsu,1))   
          a0tsui=(1.+2.*sh/th)*pyhiw1au(epss,2)+(1.+2.*uh/th)*    
     &    pyhiw1au(epsu,2)+0.5*((sh2+uh2)/th2-epst)*(pyhiw2au(epss,2)+  
     &    pyhiw2au(epsu,2))-
     &    0.25*epss*(1.-0.5*epst)*(pyhii3au(betus,epst,2)+    
     &    pyhii3au(betus,epss,2))-0.25*epsu*(1.-0.5*epst)*    
     &    (pyhii3au(betsu,epst,2)+pyhii3au(betsu,epsu,2))+  
     &    0.25*(-2.*(sh2+uh2)/th2+4.*epst+epss+epsu+0.5*epss*epsu)* 
     &    (pyhii3au(bestu,epss,2)+pyhii3au(bestu,epsu,2))   
          a0utsr=1.+(1.+2.*th/uh)*pyhiw1au(epst,1)+(1.+2.*sh/uh)* 
     &    pyhiw1au(epss,1)+0.5*((th2+sh2)/uh2-epsu)*(pyhiw2au(epst,1)+  
     &    pyhiw2au(epss,1))-
     &    0.25*epst*(1.-0.5*epsu)*(pyhii3au(beust,epsu,1)+    
     &    pyhii3au(beust,epst,1))-0.25*epss*(1.-0.5*epsu)*    
     &    (pyhii3au(beuts,epsu,1)+pyhii3au(beuts,epss,1))+  
     &    0.25*(-2.*(th2+sh2)/uh2+4.*epsu+epst+epss+0.5*epst*epss)* 
     &    (pyhii3au(betus,epst,1)+pyhii3au(betus,epss,1))   
          a0utsi=(1.+2.*th/uh)*pyhiw1au(epst,2)+(1.+2.*sh/uh)*    
     &    pyhiw1au(epss,2)+0.5*((th2+sh2)/uh2-epsu)*(pyhiw2au(epst,2)+  
     &    pyhiw2au(epss,2))-
     &    0.25*epst*(1.-0.5*epsu)*(pyhii3au(beust,epsu,2)+    
     &    pyhii3au(beust,epst,2))-0.25*epss*(1.-0.5*epsu)*    
     &    (pyhii3au(beuts,epsu,2)+pyhii3au(beuts,epss,2))+  
     &    0.25*(-2.*(th2+sh2)/uh2+4.*epsu+epst+epss+0.5*epst*epss)* 
     &    (pyhii3au(betus,epst,2)+pyhii3au(betus,epss,2))   
          a1stur=-1.-0.25*(epss+epst+epsu)*(pyhiw2au(epss,1)+ 
     &    pyhiw2au(epst,1)+pyhiw2au(epsu,1))+0.25*(epsu+0.5*epss*epst)* 
     &    (pyhii3au(besut,epss,1)+pyhii3au(besut,epst,1))+  
     &    0.25*(epst+0.5*epss*epsu)*(pyhii3au(bestu,epss,1)+  
     &    pyhii3au(bestu,epsu,1))+0.25*(epss+0.5*epst*epsu)*  
     &    (pyhii3au(betsu,epst,1)+pyhii3au(betsu,epsu,1))   
          a1stui=-0.25*(epss+epst+epsu)*(pyhiw2au(epss,2)+
     &    pyhiw2au(epst,2)+ 
     &    pyhiw2au(epsu,2))+0.25*(epsu+0.5*epss*epst)*    
     &    (pyhii3au(besut,epss,2)+pyhii3au(besut,epst,2))+  
     &    0.25*(epst+0.5*epss*epsu)*(pyhii3au(bestu,epss,2)+  
     &    pyhii3au(bestu,epsu,2))+0.25*(epss+0.5*epst*epsu)*  
     &    (pyhii3au(betsu,epst,2)+pyhii3au(betsu,epsu,2))   
          a2stur=-1.+0.125*epss*epst*(pyhii3au(besut,epss,1)+ 
     &    pyhii3au(besut,epst,1))+
     &    0.125*epss*epsu*(pyhii3au(bestu,epss,1)+  
     &    pyhii3au(bestu,epsu,1))+
     &    0.125*epst*epsu*(pyhii3au(betsu,epst,1)+  
     &    pyhii3au(betsu,epsu,1)) 
          a2stui=0.125*epss*epst*(pyhii3au(besut,epss,2)+ 
     &    pyhii3au(besut,epst,2))+
     &    0.125*epss*epsu*(pyhii3au(bestu,epss,2)+  
     &    pyhii3au(bestu,epsu,2))+
     &    0.125*epst*epsu*(pyhii3au(betsu,epst,2)+  
     &    pyhii3au(betsu,epsu,2)) 
        ENDIF   
        asre=asre+ei**2*(a0stur+a0tsur+a0utsr+4.*a1stur+a2stur) 
        asim=asim+ei**2*(a0stui+a0tsui+a0utsi+4.*a1stui+a2stui) 
  780   CONTINUE    
        facgg=comfac*faca/(8.*paru(1)**2)*as**2*aem**2*(asre**2+asim**2)    
        IF(kfac(1,21)*kfac(2,21).EQ.0) goto 790 
        nchn=nchn+1 
        isig(nchn,1)=21 
        isig(nchn,2)=21 
        isig(nchn,3)=1  
        sigh(nchn)=facgg    
  790   CONTINUE    
    
      ELSEIF(isub.EQ.115) THEN  
C...g + g -> gamma + Z0.    
    
      ELSEIF(isub.EQ.116) THEN  
C...g + g -> Z0 + Z0.   
    
      ELSEIF(isub.EQ.117) THEN  
C...g + g -> W+ + W-.   
    
      ENDIF 
    
C...G: 2 -> 3, tree diagrams.   
    
      ELSEIF(isub.LE.140) THEN  
      IF(isub.EQ.121) THEN  
C...g + g -> f + fb + H0.   
    
      ENDIF 
    
C...H: 2 -> 1, tree diagrams, non-standard model processes. 
    
      ELSEIF(isub.LE.160) THEN  
      IF(isub.EQ.141) THEN  
C...f + fb -> gamma*/Z0/Z'0.    
        mint(61)=2  
        CALL pyhiwidt(32,sqrt(sh),wdtp,wdte)  
        faczp=comfac*aem**2*4./9.   
        DO 800 i=mina,maxa  
        IF(i.EQ.0.OR.kfac(1,i)*kfac(2,-i).EQ.0) goto 800    
        ei=kchg(iabs(i),1)/3.   
        ai=sign(1.,ei)  
        vi=ai-4.*ei*xw  
        api=sign(1.,ei) 
        vpi=api-4.*ei*xw    
        nchn=nchn+1 
        isig(nchn,1)=i  
        isig(nchn,2)=-i 
        isig(nchn,3)=1  
        sigh(nchn)=faczp*(ei**2*vint(111)+ei*vi/(8.*xw*(1.-xw))*    
     &  sh*(sh-sqmz)/((sh-sqmz)**2+gmmz**2)*vint(112)+ei*vpi/(8.*xw*    
     &  (1.-xw))*sh*(sh-sqmzp)/((sh-sqmzp)**2+gmmzp**2)*vint(113)+  
     &  (vi**2+ai**2)/(16.*xw*(1.-xw))**2*sh2/((sh-sqmz)**2+gmmz**2)*   
     &  vint(114)+2.*(vi*vpi+ai*api)/(16.*xw*(1.-xw))**2*sh2*   
     &  ((sh-sqmz)*(sh-sqmzp)+gmmz*gmmzp)/(((sh-sqmz)**2+gmmz**2)*  
     &  ((sh-sqmzp)**2+gmmzp**2))*vint(115)+(vpi**2+api**2)/    
     &  (16.*xw*(1.-xw))**2*sh2/((sh-sqmzp)**2+gmmzp**2)*vint(116)) 
  800   CONTINUE    
    
      ELSEIF(isub.EQ.142) THEN  
C...f + fb' -> H+/-.    
        CALL pyhiwidt(37,sqrt(sh),wdtp,wdte)  
        fhc=comfac*(aem/xw)**2*1./48.*(sh/sqmw)**2*sh2/ 
     &  ((sh-sqmhc)**2+gmmhc**2)    
C'''No construction yet for leptons 
        DO 840 i=1,mstp(54)/2   
        il=2*i-1    
        iu=2*i  
        rmql=pmas(il,1)**2/sh   
        rmqu=pmas(iu,1)**2/sh   
        fachc=fhc*((rmql*paru(121)+rmqu/paru(121))*(1.-rmql-rmqu)-  
     &  4.*rmql*rmqu)/sqrt(max(0.,(1.-rmql-rmqu)**2-4.*rmql*rmqu))  
        IF(kfac(1,il)*kfac(2,-iu).EQ.0) goto 810    
        kchhc=(kchg(il,1)-kchg(iu,1))/3 
        nchn=nchn+1 
        isig(nchn,1)=il 
        isig(nchn,2)=-iu    
        isig(nchn,3)=1  
        sigh(nchn)=fachc*(wdte(0,1)+wdte(0,(5-kchhc)/2)+wdte(0,4))  
  810   IF(kfac(1,-il)*kfac(2,iu).EQ.0) goto 820    
        kchhc=(-kchg(il,1)+kchg(iu,1))/3    
        nchn=nchn+1 
        isig(nchn,1)=-il    
        isig(nchn,2)=iu 
        isig(nchn,3)=1  
        sigh(nchn)=fachc*(wdte(0,1)+wdte(0,(5-kchhc)/2)+wdte(0,4))  
  820   IF(kfac(1,iu)*kfac(2,-il).EQ.0) goto 830    
        kchhc=(kchg(iu,1)-kchg(il,1))/3 
        nchn=nchn+1 
        isig(nchn,1)=iu 
        isig(nchn,2)=-il    
        isig(nchn,3)=1  
        sigh(nchn)=fachc*(wdte(0,1)+wdte(0,(5-kchhc)/2)+wdte(0,4))  
  830   IF(kfac(1,-iu)*kfac(2,il).EQ.0) goto 840    
        kchhc=(-kchg(iu,1)+kchg(il,1))/3    
        nchn=nchn+1 
        isig(nchn,1)=-iu    
        isig(nchn,2)=il 
        isig(nchn,3)=1  
        sigh(nchn)=fachc*(wdte(0,1)+wdte(0,(5-kchhc)/2)+wdte(0,4))  
  840   CONTINUE    
    
      ELSEIF(isub.EQ.143) THEN  
C...f + fb -> R.    
        CALL pyhiwidt(40,sqrt(sh),wdtp,wdte)  
        facr=comfac*(aem/xw)**2*1./9.*sh2/((sh-sqmr)**2+gmmr**2)    
        DO 860 i=min1,max1  
        IF(i.EQ.0.OR.kfac(1,i).EQ.0) goto 860   
        ia=iabs(i)  
        DO 850 j=min2,max2  
        IF(j.EQ.0.OR.kfac(2,j).EQ.0) goto 850   
        ja=iabs(j)  
        IF(i*j.GT.0.OR.iabs(ia-ja).NE.2) goto 850   
        nchn=nchn+1 
        isig(nchn,1)=i  
        isig(nchn,2)=j  
        isig(nchn,3)=1  
        sigh(nchn)=facr*(wdte(0,1)+wdte(0,(10-(i+j))/4)+wdte(0,4))  
  850   CONTINUE    
  860   CONTINUE    
    
      ENDIF 
    
C...I: 2 -> 2, tree diagrams, non-standard model processes. 
    
      ELSE  
      IF(isub.EQ.161) THEN  
C...f + g -> f' + H+/- (q + g -> q' + H+/- only).   
        fhcq=comfac*faca*as*aem/xw*1./24    
        DO 900 i=1,mstp(54) 
        iu=i+mod(i,2)   
        sqmq=pmas(iu,1)**2  
        fachcq=fhcq/paru(121)*sqmq/sqmw*(sh/(sqmq-uh)+  
     &  2.*sqmq*(sqmhc-uh)/(sqmq-uh)**2+(sqmq-uh)/sh+   
     &  2.*sqmq/(sqmq-uh)+2.*(sqmhc-uh)/(sqmq-uh)*(sqmhc-sqmq-sh)/sh)   
        IF(kfac(1,-i)*kfac(2,21).EQ.0) goto 870 
        kchhc=isign(1,-kchg(i,1))   
        nchn=nchn+1 
        isig(nchn,1)=-i 
        isig(nchn,2)=21 
        isig(nchn,3)=1  
        sigh(nchn)=fachcq*(wdte(0,1)+wdte(0,(5-kchhc)/2)+wdte(0,4)) 
  870   IF(kfac(1,i)*kfac(2,21).EQ.0) goto 880  
        kchhc=isign(1,kchg(i,1))    
        nchn=nchn+1 
        isig(nchn,1)=i  
        isig(nchn,2)=21 
        isig(nchn,3)=1  
        sigh(nchn)=fachcq*(wdte(0,1)+wdte(0,(5-kchhc)/2)+wdte(0,4)) 
  880   IF(kfac(1,21)*kfac(2,-i).EQ.0) goto 890 
        kchhc=isign(1,-kchg(i,1))   
        nchn=nchn+1 
        isig(nchn,1)=21 
        isig(nchn,2)=-i 
        isig(nchn,3)=1  
        sigh(nchn)=fachcq*(wdte(0,1)+wdte(0,(5-kchhc)/2)+wdte(0,4)) 
  890   IF(kfac(1,21)*kfac(2,i).EQ.0) goto 900  
        kchhc=isign(1,kchg(i,1))    
        nchn=nchn+1 
        isig(nchn,1)=21 
        isig(nchn,2)=i  
        isig(nchn,3)=1  
        sigh(nchn)=fachcq*(wdte(0,1)+wdte(0,(5-kchhc)/2)+wdte(0,4)) 
  900   CONTINUE    
    
      ENDIF 
      ENDIF 
    
C...Multiply with structure functions.  
      IF(isub.LE.90.OR.isub.GE.96) THEN 
        DO 910 ichn=1,nchn  
        IF(mint(41).EQ.2) THEN  
          kfl1=isig(ichn,1) 
          IF(kfl1.EQ.21) kfl1=0 
          sigh(ichn)=sigh(ichn)*xsfx(1,kfl1)    
        ENDIF   
        IF(mint(42).EQ.2) THEN  
          kfl2=isig(ichn,2) 
          IF(kfl2.EQ.21) kfl2=0 
          sigh(ichn)=sigh(ichn)*xsfx(2,kfl2)    
        ENDIF   
  910   sigs=sigs+sigh(ichn)    
      ENDIF 
    
      RETURN    
      END