pythiaeRHIC/pythia/pysigh.f

0001
0002 C***********************************************************************
0003
0004 C...PYSIGH
0005 C...Differential matrix elements for all included subprocesses
0006 C...Note that what is coded is (disregarding the COMFAC factor)
0007 C...1) for 2 -> 1 processes: s-hat/pi*d(sigma-hat), where,
0008 C...when d(sigma-hat) is given in the zero-width limit, the delta
0009 C...function in tau is replaced by a (modified) Breit-Wigner:
0010 C...1/pi*s*H_res/((s*tau-m_res^2)^2+H_res^2),
0011 C...where H_res = s-hat/m_res*Gamma_res(s-hat);
0012 C...2) for 2 -> 2 processes: (s-hat)**2/pi*d(sigma-hat)/d(t-hat);
0013 C...i.e., dimensionless quantities
0014 C...3) for 2 -> 3 processes: abs(M)^2, where the total cross-section is
0015 C...Integral abs(M)^2/(2shat') * (prod_(i=1)^3 d^3p_i/((2pi)^3*2E_i)) *
0016 C...(2pi)^4 delta^4(P - sum p_i)
0017 C...COMFAC contains the factor pi/s (or equivalent) and
0018 C...the conversion factor from GeV^-2 to mb
0019
0020       SUBROUTINE PYSIGH(NCHN,SIGS)
0021
0022 C...Double precision and integer declarations
0023       IMPLICIT DOUBLE PRECISION(A-H, O-Z)
0024       IMPLICIT INTEGER(I-N)
0025       INTEGER PYK,PYCHGE,PYCOMP
0026 C...Parameter statement to help give large particle numbers.
0027       PARAMETER (KSUSY1=1000000,KSUSY2=2000000,KTECHN=3000000,
0028      &KEXCIT=4000000,KDIMEN=5000000)
0029 C...Commonblocks
0030       COMMON/PYJETS/N,NPAD,K(4000,5),P(4000,5),V(4000,5)
0031       COMMON/PYDAT1/MSTU(200),PARU(200),MSTJ(200),PARJ(200)
0032       COMMON/PYDAT2/KCHG(500,4),PMAS(500,4),PARF(2000),VCKM(4,4)
0033       COMMON/PYDAT3/MDCY(500,3),MDME(8000,2),BRAT(8000),KFDP(8000,5)
0034       COMMON/PYSUBS/MSEL,MSELPD,MSUB(500),KFIN(2,-40:40),CKIN(200)
0035       COMMON/PYPARS/MSTP(200),PARP(200),MSTI(200),PARI(200)
0036       COMMON/PYINT1/MINT(400),VINT(400)
0037       COMMON/PYINT2/ISET(500),KFPR(500,2),COEF(500,20),ICOL(40,4,2)
0038       COMMON/PYINT3/XSFX(2,-40:40),ISIG(1000,3),SIGH(1000)
0039       COMMON/PYINT4/MWID(500),WIDS(500,5)
0040       COMMON/PYINT5/NGENPD,NGEN(0:500,3),XSEC(0:500,3)
0041       COMMON/PYINT7/SIGT(0:6,0:6,0:5)
0042       COMMON/PYMSSM/IMSS(0:99),RMSS(0:99)
0043       COMMON/PYSSMT/ZMIX(4,4),UMIX(2,2),VMIX(2,2),SMZ(4),SMW(2),
0044      &SFMIX(16,4),ZMIXI(4,4),UMIXI(2,2),VMIXI(2,2)
0045       COMMON/PYTCSM/ITCM(0:99),RTCM(0:99)
0046       COMMON/PYSGCM/ISUB,ISUBSV,MMIN1,MMAX1,MMIN2,MMAX2,MMINA,MMAXA,
0047      &KFAC(2,-40:40),COMFAC,FACK,FACA,SH,TH,UH,SH2,TH2,UH2,SQM3,SQM4,
0048      &SHR,SQPTH,TAUP,BE34,CTH,X(2),SQMZ,SQMW,GMMZ,GMMW,
0049      &AEM,AS,XW,XW1,XWC,XWV,POLL,POLR,POLLL,POLRR
0050       SAVE /PYJETS/,/PYDAT1/,/PYDAT2/,/PYDAT3/,/PYSUBS/,/PYPARS/,
0051      &/PYINT1/,/PYINT2/,/PYINT3/,/PYINT4/,/PYINT5/,/PYINT7/,
0052      &/PYMSSM/,/PYSSMT/,/PYTCSM/,/PYSGCM/
0053       include "mc_set.inc"
0054 C...Local arrays and complex variables
0055       DIMENSION XPQ(-25:25)
0056
0057 C...Map of processes onto which routine to call
0058 C...in order to evaluate cross section:
0059 C...0 = not implemented;
0060 C...1 = standard QCD (including photons);
0061 C...2 = heavy flavours;
0062 C...3 = W/Z;
0063 C...4 = Higgs (2 doublets; including longitudinal W/Z scattering);
0064 C...5 = SUSY;
0065 C...6 = Technicolor;
0066 C...7 = exotics (Z'/W'/LQ/R/f*/H++/Z_R/W_R/G*).
0067       DIMENSION MAPPR(500)
0068       DATA (MAPPR(I),I=1,180)/
0069      &    3,  3,  4,  0,  4,  0,  0,  4,  0,  1,
0070      1    1,  1,  1,  1,  3,  3,  0,  1,  3,  3,
0071      2    0,  3,  3,  4,  3,  4,  0,  1,  1,  3,
0072      3    3,  4,  1,  1,  3,  3,  0,  0,  0,  0,
0073      4    0,  0,  0,  0,  0,  0,  0,  0,  0,  0,
0074      5    0,  0,  1,  1,  0,  0,  0,  1,  0,  0,
0075      6    0,  0,  0,  0,  0,  0,  0,  1,  3,  3,
0076      7    4,  4,  4,  0,  0,  4,  4,  0,  0,  1,
0077      8    2,  2,  2,  2,  2,  2,  2,  2,  2,  0,
0078      9    1,  1,  1,  1,  1,  1,  0,  0,  1,  0,
0079      &    0,  4,  4,  2,  2,  2,  2,  2,  0,  4,
0080      1    4,  4,  4,  1,  1,  0,  0,  0,  0,  0,
0081      2    4,  4,  4,  4,  0,  0,  0,  0,  0,  0,
0082      3    1,  1,  1,  1,  1,  1,  1,  1,  1,  1,
0083      4    7,  7,  4,  7,  7,  7,  7,  7,  6,  0,
0084      5    4,  4,  4,  0,  0,  4,  4,  4,  0,  0,
0085      6    4,  7,  7,  7,  6,  6,  7,  7,  7,  0,
0086      7    4,  4,  4,  4,  0,  4,  4,  4,  4,  0/
0087       DATA (MAPPR(I),I=181,500)/
0088      8    4,  4,  4,  4,  4,  4,  4,  4,  4,  4,
0089      9    6,  6,  6,  6,  6,  0,  0,  0,  0,  0,
0090      &    100*5,
0091      &    5,  0,  0,  0,  0,  0,  0,  0,  0,  0,
0092      1     30*0,
0093      4    7,  7,  7,  7,  7,  7,  7,  7,  7,  7,
0094      5    7,  7,  7,  7,  0,  0,  0,  0,  0,  0,
0095      6    6,  6,  6,  6,  6,  6,  6,  6,  0,  6,
0096      7    6,  6,  6,  6,  6,  6,  6,  0,  0,  0,
0097      8    6,  6,  6,  6,  6,  6,  6,  6,  0,  0,
0098      9    7,  7,  7,  7,  7,  0,  0,  0,  0,  0,
0099      &    4,  4,  18*0,
0100      2    2,  2,  2,  2,  2,  2,  2,  2,  2,  2,
0101      3    2,  2,  2,  2,  2,  2,  2,  2,  2,  0,
0102      4     20*0,
0103      6    2,  2,  2,  2,  2,  2,  2,  2,  2,  2,
0104      7    2,  2,  2,  2,  2,  2,  2,  2,  2,  0,
0105      8     20*0/
0106
0107 C...Reset number of channels and cross-section
0108       NCHN=0
0109       SIGS=0D0
0110
0111 C...Read process to consider.
0112       ISUB=MINT(1)
0113       ISUBSV=ISUB
0114       MAP=MAPPR(ISUB)
0115
0116 C...Read kinematical variables and limits
0117       ISTSB=ISET(ISUBSV)
0118       TAUMIN=VINT(11)
0119       YSTMIN=VINT(12)
0120       CTNMIN=VINT(13)
0121       CTPMIN=VINT(14)
0122       TAUPMN=VINT(16)
0123       TAU=VINT(21)
0124       YST=VINT(22)
0125       CTH=VINT(23)
0126       XT2=VINT(25)
0127       TAUP=VINT(26)
0128       TAUMAX=VINT(31)
0129       YSTMAX=VINT(32)
0130       CTNMAX=VINT(33)
0131       CTPMAX=VINT(34)
0132       TAUPMX=VINT(36)
0133
0134 C...Derive kinematical quantities
0135       TAUE=TAU
0136       IF(ISTSB.GE.3.AND.ISTSB.LE.5) TAUE=TAUP
0137       X(1)=SQRT(TAUE)*EXP(YST)
0138       X(2)=SQRT(TAUE)*EXP(-YST)
0139       IF(MINT(45).EQ.2.AND.ISTSB.GE.1) THEN
0140         IF(X(1).GT.1D0-1D-7) RETURN
0141       ELSEIF(MINT(45).EQ.3) THEN
0142         X(1)=MIN(1D0-1.1D-10,X(1))
0143       ENDIF
0144       IF(MINT(46).EQ.2.AND.ISTSB.GE.1) THEN
0145         IF(X(2).GT.1D0-1D-7) RETURN
0146       ELSEIF(MINT(46).EQ.3) THEN
0147         X(2)=MIN(1D0-1.1D-10,X(2))
0148       ENDIF
0149       SH=MAX(1D0,TAU*VINT(2))
0150       SQM3=VINT(63)
0151       SQM4=VINT(64)
0152       RM3=SQM3/SH
0153       RM4=SQM4/SH
0154       BE34=SQRT(MAX(0D0,(1D0-RM3-RM4)**2-4D0*RM3*RM4))
0155       RPTS=4D0*VINT(71)**2/SH
0156       BE34L=SQRT(MAX(0D0,(1D0-RM3-RM4)**2-4D0*RM3*RM4-RPTS))
0157       RM34=MAX(1D-20,2D0*RM3*RM4)
0158       RSQM=1D0+RM34
0159       IF(2D0*VINT(71)**2/MAX(1D0,VINT(21)*VINT(2)).LT.0.0001D0)
0160      &RM34=MAX(RM34,2D0*VINT(71)**2/MAX(1D0,VINT(21)*VINT(2)))
0161       RTHM=(4D0*RM3*RM4+RPTS)/(1D0-RM3-RM4+BE34L)
0162       IF(ISTSB.EQ.0) THEN
0163         TH=VINT(45)
0164         UH=-0.5D0*SH*MAX(RTHM,1D0-RM3-RM4+BE34*CTH)
0165         SQPTH=MAX(VINT(71)**2,0.25D0*SH*BE34**2*VINT(59)**2)
0166       ELSE
0167 C...Kinematics with incoming masses tricky: now depends on how
0168 C...subprocess has been set up w.r.t. order of incoming partons.
0169         RM1=0D0
0170         IF(MINT(15).EQ.22.AND.VINT(3).LT.0D0) RM1=-VINT(3)**2/SH
0171         RM2=0D0
0172         IF(MINT(16).EQ.22.AND.VINT(4).LT.0D0) RM2=-VINT(4)**2/SH
0173         IF(ISUB.EQ.35) THEN
0174           RM2=MIN(RM1,RM2)
0175           RM1=0D0
0176         ENDIF
0177         BE12=SQRT(MAX(0D0,(1D0-RM1-RM2)**2-4D0*RM1*RM2))
0178         TUCOM=(1D0-RM1-RM2)*(1D0-RM3-RM4)
0179         TH=-0.5D0*SH*MAX(RTHM,TUCOM-2D0*RM1*RM4-2D0*RM2*RM3-
0180      &  BE12*BE34*CTH)
0181         UH=-0.5D0*SH*MAX(RTHM,TUCOM-2D0*RM1*RM3-2D0*RM2*RM4+
0182      &  BE12*BE34*CTH)
0183         SQPTH=MAX(VINT(71)**2,0.25D0*SH*BE34**2*(1D0-CTH**2))
0184       ENDIF
0185       SHR=SQRT(SH)
0186       SH2=SH**2
0187       TH2=TH**2
0188       UH2=UH**2
0189
0190 C...Choice of Q2 scale for hard process (e.g. alpha_s).
0191       IF(ISTSB.EQ.1.OR.ISTSB.EQ.3.OR.ISTSB.EQ.5) THEN
0192         Q2=SH
0193       ELSEIF(ISTSB.EQ.8) THEN
0194         IF(MINT(107).EQ.4) Q2=VINT(307)
0195         IF(MINT(108).EQ.4) Q2=VINT(308)
0196       ELSEIF(MOD(ISTSB,2).EQ.0.OR.ISTSB.EQ.9) THEN
0197         Q2IN1=0D0
0198         IF(MINT(11).EQ.22.AND.VINT(3).LT.0D0) Q2IN1=VINT(3)**2
0199         Q2IN2=0D0
0200         IF(MINT(12).EQ.22.AND.VINT(4).LT.0D0) Q2IN2=VINT(4)**2
0201         IF(MSTP(32).EQ.1) THEN
0202           Q2=2D0*SH*TH*UH/(SH**2+TH**2+UH**2)
0203         ELSEIF(MSTP(32).EQ.2) THEN
0204           Q2=SQPTH+0.5D0*(SQM3+SQM4)
0205         ELSEIF(MSTP(32).EQ.3) THEN
0206           Q2=MIN(-TH,-UH)
0207         ELSEIF(MSTP(32).EQ.4) THEN
0208           Q2=SH
0209         ELSEIF(MSTP(32).EQ.5) THEN
0210           Q2=-TH
0211         ELSEIF(MSTP(32).EQ.6) THEN
0212           XSF1=X(1)
0213           IF(ISTSB.EQ.9) XSF1=X(1)/VINT(143)
0214           XSF2=X(2)
0215           IF(ISTSB.EQ.9) XSF2=X(2)/VINT(144)
0216           Q2=(1D0+XSF1*Q2IN1/SH+XSF2*Q2IN2/SH)*
0217      &    (SQPTH+0.5D0*(SQM3+SQM4))
0218         ELSEIF(MSTP(32).EQ.7) THEN
0219           Q2=(1D0+Q2IN1/SH+Q2IN2/SH)*(SQPTH+0.5D0*(SQM3+SQM4))
0220         ELSEIF(MSTP(32).EQ.8) THEN
0221           Q2=SQPTH+0.5D0*(Q2IN1+Q2IN2+SQM3+SQM4)
0222         ELSEIF(MSTP(32).EQ.9) THEN
0223           Q2=SQPTH+Q2IN1+Q2IN2+SQM3+SQM4
0224         ELSEIF(MSTP(32).EQ.10) THEN
0225           Q2=VINT(2)
0226 C..Begin JA 040914
0227         ELSEIF(MSTP(32).EQ.11) THEN
0228           Q2=0.25*(SQM3+SQM4+2*SQRT(SQM3*SQM4))
0229         ELSEIF(MSTP(32).EQ.12) THEN
0230           Q2=PARP(193)
0231 C..End JA
0232         ELSEIF(MSTP(32).EQ.13) THEN
0233           Q2=SQPTH
0234         ENDIF
0235         IF(MINT(35).LE.2.AND.ISTSB.EQ.9) Q2=SQPTH
0236         IF(ISTSB.EQ.9.AND.MSTP(82).GE.2) Q2=Q2+
0237      &  (PARP(82)*(VINT(1)/PARP(89))**PARP(90))**2
0238       ENDIF
0239
0240 C...Choice of Q2 scale for parton densities.
0241       Q2SF=Q2
0242 C..Begin JA 040914
0243       IF(MSTP(32).EQ.12.AND.(MOD(ISTSB,2).EQ.0.OR.ISTSB.EQ.9)
0244      &     .OR.MSTP(39).EQ.8.AND.(ISTSB.GE.3.AND.ISTSB.LE.5))
0245      &     Q2=PARP(194)
0246 C..End JA
0247       IF(ISTSB.GE.3.AND.ISTSB.LE.5) THEN
0248         Q2SF=PMAS(23,1)**2
0249         IF(ISUB.EQ.8.OR.ISUB.EQ.76.OR.ISUB.EQ.77.OR.ISUB.EQ.124.OR.
0250      &  ISUB.EQ.174.OR.ISUB.EQ.179.OR.ISUB.EQ.351) Q2SF=PMAS(24,1)**2
0251         IF(ISUB.EQ.352) Q2SF=PMAS(PYCOMP(9900024),1)**2
0252         IF(ISUB.EQ.121.OR.ISUB.EQ.122.OR.ISUB.EQ.181.OR.ISUB.EQ.182.OR.
0253      &  ISUB.EQ.186.OR.ISUB.EQ.187.OR.ISUB.EQ.401.OR.ISUB.EQ.402) THEN
0254           Q2SF=PMAS(PYCOMP(KFPR(ISUBSV,2)),1)**2
0255           IF(MSTP(39).EQ.2) Q2SF=
0256      &         MAX(VINT(201)**2+VINT(202),VINT(206)**2+VINT(207))
0257           IF(MSTP(39).EQ.3) Q2SF=SH
0258           IF(MSTP(39).EQ.4) Q2SF=VINT(26)*VINT(2)
0259           IF(MSTP(39).EQ.5) Q2SF=PMAS(PYCOMP(KFPR(ISUBSV,1)),1)**2
0260 C..Begin JA 040914
0261           IF(MSTP(39).EQ.6) Q2SF=0.25*(VINT(201)+SQRT(SH))**2
0262           IF(MSTP(39).EQ.7) Q2SF=
0263      &         (VINT(201)**2+VINT(202)+VINT(206)**2+VINT(207))/2d0
0264           IF(MSTP(39).EQ.8) Q2SF=PARP(193)
0265 C..End JA
0266         ENDIF
0267       ENDIF
0268       IF(MINT(35).GE.3.AND.ISTSB.EQ.9) Q2SF=SQPTH
0269
0270       Q2PS=Q2SF
0271       Q2SF=Q2SF*PARP(34)
0272       IF(MSTP(69).GE.1.AND.MINT(47).EQ.5) Q2SF=VINT(2)
0273       IF(MSTP(69).GE.2) Q2SF=VINT(2)
0274
0275 C...Identify to which class(es) subprocess belongs
0276       ISMECR=0
0277       ISQCD=0
0278       ISJETS=0
0279       IF (ISUBSV.EQ.1.OR.ISUBSV.EQ.2.OR.ISUBSV.EQ.3.OR.
0280      &     ISUBSV.EQ.102.OR.ISUBSV.EQ.141.OR.ISUBSV.EQ.142.OR.
0281      &     ISUBSV.EQ.144.OR.ISUBSV.EQ.151.OR.ISUBSV.EQ.152.OR.
0282      &     ISUBSV.EQ.156.OR.ISUBSV.EQ.157) ISMECR=1
0283       IF (ISUBSV.EQ.11.OR.ISUBSV.EQ.12.OR.ISUBSV.EQ.13.OR.
0284      &     ISUBSV.EQ.28.OR.ISUBSV.EQ.53.OR.ISUBSV.EQ.68) ISQCD=1
0285       IF ((ISUBSV.EQ.81.OR.ISUBSV.EQ.82).AND.MINT(55).LE.5) ISQCD=1
0286       IF (ISUBSV.GE.381.AND.ISUBSV.LE.386) ISQCD=1
0287       IF ((ISUBSV.EQ.387.OR.ISUBSV.EQ.388).AND.MINT(55).LE.5) ISQCD=1
0288       IF (ISTSB.EQ.9) ISQCD=1
0289       IF ((ISUBSV.GE.86.AND.ISUBSV.LE.89).OR.ISUBSV.EQ.107.OR.
0290      &     (ISUBSV.GE.14.AND.ISUBSV.LE.16).OR.(ISUBSV.GE.29.AND.
0291      &     ISUBSV.LE.32).OR.(ISUBSV.GE.111.AND.ISUBSV.LE.113).OR.
0292      &     ISUBSV.EQ.115.OR.(ISUBSV.GE.183.AND.ISUBSV.LE.185).OR.
0293      &     (ISUBSV.GE.188.AND.ISUBSV.LE.190).OR.ISUBSV.EQ.161.OR.
0294      &     ISUBSV.EQ.167.OR.ISUBSV.EQ.168.OR.(ISUBSV.GE.393.AND.
0295      &     ISUBSV.LE.395).OR.(ISUBSV.GE.421.AND.ISUBSV.LE.439).OR.
0296      &     (ISUBSV.GE.461.AND.ISUBSV.LE.479)) ISJETS=1
0297 C...WBF is special case of ISJETS
0298       IF (ISUBSV.EQ.5.OR.ISUBSV.EQ.8.OR.
0299      &    (ISUBSV.GE.71.AND.ISUBSV.LE.73).OR.
0300      &    ISUBSV.EQ.76.OR.ISUBSV.EQ.77.OR.
0301      &    (ISUBSV.GE.121.AND.ISUBSV.LE.124).OR.
0302      &    ISUBSV.EQ.173.OR.ISUBSV.EQ.174.OR.
0303      &    ISUBSV.EQ.178.OR.ISUBSV.EQ.179.OR.
0304      &    ISUBSV.EQ.181.OR.ISUBSV.EQ.182.OR.
0305      &    ISUBSV.EQ.186.OR.ISUBSV.EQ.187.OR.
0306      &    ISUBSV.EQ.351.OR.ISUBSV.EQ.352) ISJETS=2
0307 C...Some processes with photons also belong here.
0308       IF (ISUBSV.EQ.10.OR.(ISUBSV.GE.18.AND.ISUBSV.LE.20).OR.
0309      &     (ISUBSV.GE.33.AND.ISUBSV.LE.36).OR.ISUBSV.EQ.54.OR.
0310      &     ISUBSV.EQ.58.OR.ISUBSV.EQ.69.OR.ISUBSV.EQ.70.OR.
0311      &     ISUBSV.EQ.80.OR.(ISUBSV.GE.83.AND.ISUBSV.LE.85).OR.
0312      &     (ISUBSV.GE.106.AND.ISUBSV.LE.110).OR.ISUBSV.EQ.114.OR.
0313      &     (ISUBSV.GE.131.AND.ISUBSV.LE.140)) ISJETS=3
0314
0315 C...Choice of Q2 scale for parton-shower activity.
0316       IF(MSTP(22).GE.1.AND.(ISUB.EQ.10.OR.ISUB.EQ.83).AND.
0317      &(MINT(43).EQ.2.OR.MINT(43).EQ.3)) THEN
0318         XBJ=X(2)
0319         IF(MINT(43).EQ.3) XBJ=X(1)
0320         IF(MSTP(22).EQ.1) THEN
0321           Q2PS=-TH
0322         ELSEIF(MSTP(22).EQ.2) THEN
0323           Q2PS=((1D0-XBJ)/XBJ)*(-TH)
0324         ELSEIF(MSTP(22).EQ.3) THEN
0325           Q2PS=SQRT((1D0-XBJ)/XBJ)*(-TH)
0326         ELSE
0327           Q2PS=(1D0-XBJ)*MAX(1D0,-LOG(XBJ))*(-TH)
0328         ENDIF
0329       ENDIF
0330 C...For multiple interactions, start from scale defined above
0331 C...For all other QCD or "+jets"-type events, start shower from pThard.
0332       IF (ISJETS.EQ.1.OR.ISQCD.EQ.1.AND.ISTSB.NE.9) Q2PS=SQPTH
0333       IF((MSTP(68).EQ.1.OR.MSTP(68).EQ.3).AND.ISMECR.EQ.1) THEN
0334 C...Max shower scale = s for ME corrected processes.
0335 C...(pT-ordering: max pT2 is s/4)
0336         Q2PS=VINT(2)
0337         IF (MINT(35).GE.3) Q2PS=Q2PS*0.25D0
0338       ELSEIF(MSTP(68).GE.2.AND.ISQCD.EQ.0.AND.ISJETS.EQ.0) THEN
0339 C...Max shower scale = s for all non-QCD, non-"+ jet" type processes.
0340 C...(pT-ordering: max pT2 is s/4)
0341         Q2PS=VINT(2)
0342         IF (MINT(35).GE.3) Q2PS=Q2PS*0.25D0
0343       ENDIF
0344       IF(MINT(35).EQ.2.AND.ISTSB.EQ.9) Q2PS=SQPTH
0345
0346 C...Elastic and diffractive events not associated with scales so set 0.
0347       IF(ISUBSV.GE.91.AND.ISUBSV.LE.94) THEN
0348         Q2SF=0D0
0349         Q2PS=0D0
0350       ENDIF
0351
0352 C...Store derived kinematical quantities
0353       VINT(41)=X(1)
0354       VINT(42)=X(2)
0355       VINT(44)=SH
0356       VINT(43)=SQRT(SH)
0357       VINT(45)=TH
0358       VINT(46)=UH
0359       IF(ISTSB.NE.8) VINT(48)=SQPTH
0360       IF(ISTSB.NE.8) VINT(47)=SQRT(SQPTH)
0361       VINT(50)=TAUP*VINT(2)
0362       VINT(49)=SQRT(MAX(0D0,VINT(50)))
0363       VINT(52)=Q2
0364       VINT(51)=SQRT(Q2)
0365       VINT(54)=Q2SF
0366       VINT(53)=SQRT(Q2SF)
0367       VINT(56)=Q2PS
0368       VINT(55)=SQRT(Q2PS)
0369
0370 C...Set starting scale for multiple interactions
0371       IF (ISUBSV.EQ.95) THEN
0372         XT2GMX=0D0
0373       ELSEIF(MSTP(86).EQ.3.OR.(MSTP(86).EQ.2.AND.ISUBSV.NE.11.AND.
0374      &      ISUBSV.NE.12.AND.ISUBSV.NE.13.AND.ISUBSV.NE.28.AND.
0375      &      ISUBSV.NE.53.AND.ISUBSV.NE.68.AND.ISUBSV.NE.95.AND.
0376      &      ISUBSV.NE.96)) THEN
0377 C...All accessible phase space allowed.
0378         XT2GMX=(1D0-VINT(41))*(1D0-VINT(42))
0379       ELSE
0380 C...Scale of hard process sets limit.
0381 C...2 -> 1. Limit is tau = x1*x2.
0382 C...2 -> 2. Limit is XT2 for hard process + FS masses.
0383 C...2 -> n > 2. Limit is tau' = tau of outer process.
0384         XT2GMX=VINT(25)
0385         IF(ISTSB.EQ.1) XT2GMX=VINT(21)
0386         IF(ISTSB.EQ.2)
0387      &      XT2GMX=(4D0*VINT(48)+2D0*VINT(63)+2D0*VINT(64))/VINT(2)
0388         IF(ISTSB.GE.3.AND.ISTSB.LE.5) XT2GMX=VINT(26)
0389       ENDIF
0390       VINT(62)=0.25D0*XT2GMX*VINT(2)
0391       VINT(61)=SQRT(MAX(0D0,VINT(62)))
0392
0393 C...Calculate parton distributions
0394       IF(ISTSB.LE.0) GOTO 160
0395       IF(MINT(47).GE.2) THEN
0396         DO 110 I=3-MIN(2,MINT(45)),MIN(2,MINT(46))
0397           XSF=X(I)
0398           IF(ISTSB.EQ.9) XSF=X(I)/VINT(142+I)
0399           IF(ISUB.EQ.99) THEN
0400             IF(MINT(140+I).EQ.0) THEN
0401               XSF=VINT(309-I)/(VINT(2)+VINT(309-I)-VINT(I+2)**2)
0402             ELSE
0403               XSF=VINT(309-I)/(VINT(2)+VINT(307)+VINT(308))
0404             ENDIF
0405             VINT(40+I)=XSF
0406             Q2SF=VINT(309-I)
0407           ENDIF
0408           MINT(105)=MINT(102+I)
0409           MINT(109)=MINT(106+I)
0410           VINT(120)=VINT(2+I)
0411           IF(MSTP(57).LE.1) THEN
0412             CALL PYPDFU(MINT(10+I),XSF,Q2SF,XPQ)
0413           ELSE
0414             CALL PYPDFL(MINT(10+I),XSF,Q2SF,XPQ)
0415           ENDIF
0416 C...Safety margin against heavy flavour very close to threshold,
0417 C...e.g. caused by mismatch in c and b masses.
0418           IF(Q2SF.LT.1.1*PMAS(4,1)**2) THEN
0419             XPQ(4)=0D0
0420             XPQ(-4)=0D0
0421           ENDIF
0422           IF(Q2SF.LT.1.1*PMAS(5,1)**2) THEN
0423             XPQ(5)=0D0
0424             XPQ(-5)=0D0
0425           ENDIF
0426           DO 100 KFL=-25,25
0427             XSFX(I,KFL)=XPQ(KFL)
0428   100     CONTINUE
0429   110   CONTINUE
0430       ENDIF
0431
0432 C...Calculate alpha_em, alpha_strong and K-factor
0433       XW=PARU(102)
0434       XWV=XW
0435       IF(MSTP(8).GE.2.OR.(ISUB.GE.71.AND.ISUB.LE.77)) XW=
0436      &1D0-(PMAS(24,1)/PMAS(23,1))**2
0437       XW1=1D0-XW
0438       XWC=1D0/(16D0*XW*XW1)
0439       AEM=PYALEM(Q2)
0440       IF(MSTP(8).GE.1) AEM=SQRT(2D0)*PARU(105)*PMAS(24,1)**2*XW/PARU(1)
0441       IF(MSTP(33).NE.3) AS=PYALPS(PARP(34)*Q2)
0442       FACK=1D0
0443       FACA=1D0
0444       IF(MSTP(33).EQ.1) THEN
0445         FACK=PARP(31)
0446       ELSEIF(MSTP(33).EQ.2) THEN
0447         FACK=PARP(31)
0448         FACA=PARP(32)/PARP(31)
0449       ELSEIF(MSTP(33).EQ.3) THEN
0450         Q2AS=PARP(33)*Q2
0451         IF(ISTSB.EQ.9.AND.MSTP(82).GE.2) Q2AS=Q2AS+
0452      &  PARU(112)*PARP(82)*(VINT(1)/PARP(89))**PARP(90)
0453         AS=PYALPS(Q2AS)
0454       ENDIF
0455       VINT(138)=1D0
0456       VINT(57)=AEM
0457       VINT(58)=AS
0458
0459 C...Set flags for allowed reacting partons/leptons
0460       DO 140 I=1,2
0461         DO 120 J=-25,25
0462           KFAC(I,J)=0
0463   120   CONTINUE
0464         IF(MINT(44+I).EQ.1) THEN
0465           KFAC(I,MINT(10+I))=1
0466         ELSEIF(MINT(40+I).EQ.1.AND.MSTP(12).EQ.0) THEN
0467           KFAC(I,MINT(10+I))=1
0468           KFAC(I,22)=1
0469           KFAC(I,24)=1
0470           KFAC(I,-24)=1
0471         ELSE
0472           DO 130 J=-25,25
0473             KFAC(I,J)=KFIN(I,J)
0474             IF(IABS(J).GT.MSTP(58).AND.IABS(J).LE.10) KFAC(I,J)=0
0475             IF(XSFX(I,J).LT.1D-10) KFAC(I,J)=0
0476   130     CONTINUE
0477         ENDIF
0478   140 CONTINUE
0479
0480 C...Lower and upper limit for fermion flavour loops
0481       MMIN1=0
0482       MMAX1=0
0483       MMIN2=0
0484       MMAX2=0
0485       DO 150 J=-20,20
0486         IF(KFAC(1,-J).EQ.1) MMIN1=-J
0487         IF(KFAC(1,J).EQ.1) MMAX1=J
0488         IF(KFAC(2,-J).EQ.1) MMIN2=-J
0489         IF(KFAC(2,J).EQ.1) MMAX2=J
0490   150 CONTINUE
0491       MMINA=MIN(MMIN1,MMIN2)
0492       MMAXA=MAX(MMAX1,MMAX2)
0493
0494 C...Common resonance mass and width combinations
0495       SQMZ=PMAS(23,1)**2
0496       SQMW=PMAS(24,1)**2
0497       GMMZ=PMAS(23,1)*PMAS(23,2)
0498       GMMW=PMAS(24,1)*PMAS(24,2)
0499
0500 C...Polarization factors...implemented so far for W+W-(25)
0501       POLR=(1D0+PARJ(132))*(1D0-PARJ(131))
0502       POLL=(1D0-PARJ(132))*(1D0+PARJ(131))
0503       POLRR=(1D0+PARJ(132))*(1D0+PARJ(131))
0504       POLLL=(1D0-PARJ(132))*(1D0-PARJ(131))
0505
0506 C...Phase space integral in tau
0507       COMFAC=PARU(1)*PARU(5)/VINT(2)
0508       IF(MINT(41).EQ.2.AND.MINT(42).EQ.2) COMFAC=COMFAC*FACK
0509       IF((MINT(47).GE.2.OR.(ISTSB.GE.3.AND.ISTSB.LE.5)).AND.
0510      &ISTSB.NE.8.AND.ISTSB.NE.9) THEN
0511         ATAU1=LOG(TAUMAX/TAUMIN)
0512         ATAU2=(TAUMAX-TAUMIN)/(TAUMAX*TAUMIN)
0513         H1=COEF(ISUBSV,1)+(ATAU1/ATAU2)*COEF(ISUBSV,2)/TAU
0514         IF(MINT(72).GE.1) THEN
0515           TAUR1=VINT(73)
0516           GAMR1=VINT(74)
0517           ATAUD=LOG(TAUMAX/TAUMIN*(TAUMIN+TAUR1)/(TAUMAX+TAUR1))
0518           ATAU3=ATAUD/TAUR1
0519           IF(ATAUD.GT.1D-10) H1=H1+
0520      &    (ATAU1/ATAU3)*COEF(ISUBSV,3)/(TAU+TAUR1)
0521           ATAUD=ATAN((TAUMAX-TAUR1)/GAMR1)-ATAN((TAUMIN-TAUR1)/GAMR1)
0522           ATAU4=ATAUD/GAMR1
0523           IF(ATAUD.GT.1D-10) H1=H1+
0524      &    (ATAU1/ATAU4)*COEF(ISUBSV,4)*TAU/((TAU-TAUR1)**2+GAMR1**2)
0525         ENDIF
0526         IF(MINT(72).EQ.2) THEN
0527           TAUR2=VINT(75)
0528           GAMR2=VINT(76)
0529           ATAUD=LOG(TAUMAX/TAUMIN*(TAUMIN+TAUR2)/(TAUMAX+TAUR2))
0530           ATAU5=ATAUD/TAUR2
0531           IF(ATAUD.GT.1D-10) H1=H1+
0532      &    (ATAU1/ATAU5)*COEF(ISUBSV,5)/(TAU+TAUR2)
0533           ATAUD=ATAN((TAUMAX-TAUR2)/GAMR2)-ATAN((TAUMIN-TAUR2)/GAMR2)
0534           ATAU6=ATAUD/GAMR2
0535           IF(ATAUD.GT.1D-10) H1=H1+
0536      &    (ATAU1/ATAU6)*COEF(ISUBSV,6)*TAU/((TAU-TAUR2)**2+GAMR2**2)
0537         ENDIF
0538         IF(MINT(47).EQ.5.AND.(ISTSB.LE.2.OR.ISTSB.GE.5)) THEN
0539           ATAU7=LOG(MAX(2D-10,1D0-TAUMIN)/MAX(2D-10,1D0-TAUMAX))
0540           IF(ATAU7.GT.1D-10) H1=H1+(ATAU1/ATAU7)*COEF(ISUBSV,7)*TAU/
0541      &    MAX(2D-10,1D0-TAU)
0542         ELSEIF(MINT(47).GE.6.AND.(ISTSB.LE.2.OR.ISTSB.GE.5)) THEN
0543           ATAU7=LOG(MAX(1D-10,1D0-TAUMIN)/MAX(1D-10,1D0-TAUMAX))
0544           IF(ATAU7.GT.1D-10) H1=H1+(ATAU1/ATAU7)*COEF(ISUBSV,7)*TAU/
0545      &    MAX(1D-10,1D0-TAU)
0546         ENDIF
0547         COMFAC=COMFAC*ATAU1/(TAU*H1)
0548       ENDIF
0549
0550 C...Phase space integral in y*
0551       IF((MINT(47).EQ.4.OR.MINT(47).EQ.5).AND.ISTSB.NE.8.AND.ISTSB.NE.9)
0552      &THEN
0553         AYST0=YSTMAX-YSTMIN
0554         IF(AYST0.LT.1D-10) THEN
0555           COMFAC=0D0
0556         ELSE
0557           AYST1=0.5D0*(YSTMAX-YSTMIN)**2
0558           AYST2=AYST1
0559           AYST3=2D0*(ATAN(EXP(YSTMAX))-ATAN(EXP(YSTMIN)))
0560           H2=(AYST0/AYST1)*COEF(ISUBSV,8)*(YST-YSTMIN)+
0561      &    (AYST0/AYST2)*COEF(ISUBSV,9)*(YSTMAX-YST)+
0562      &    (AYST0/AYST3)*COEF(ISUBSV,10)/COSH(YST)
0563           IF(MINT(45).EQ.3) THEN
0564             YST0=-0.5D0*LOG(TAUE)
0565             AYST4=LOG(MAX(1D-10,EXP(YST0-YSTMIN)-1D0)/
0566      &      MAX(1D-10,EXP(YST0-YSTMAX)-1D0))
0567             IF(AYST4.GT.1D-10) H2=H2+(AYST0/AYST4)*COEF(ISUBSV,11)/
0568      &      MAX(1D-10,1D0-EXP(YST-YST0))
0569           ENDIF
0570           IF(MINT(46).EQ.3) THEN
0571             YST0=-0.5D0*LOG(TAUE)
0572             AYST5=LOG(MAX(1D-10,EXP(YST0+YSTMAX)-1D0)/
0573      &      MAX(1D-10,EXP(YST0+YSTMIN)-1D0))
0574             IF(AYST5.GT.1D-10) H2=H2+(AYST0/AYST5)*COEF(ISUBSV,12)/
0575      &      MAX(1D-10,1D0-EXP(-YST-YST0))
0576           ENDIF
0577           COMFAC=COMFAC*AYST0/H2
0578         ENDIF
0579       ENDIF
0580
0581 C...2 -> 1 processes: reduction in angular part of phase space integral
0582 C...for case of decaying resonance
0583       ACTH0=CTNMAX-CTNMIN+CTPMAX-CTPMIN
0584       IF((ISTSB.EQ.1.OR.ISTSB.EQ.3.OR.ISTSB.EQ.5)) THEN
0585         IF(MDCY(PYCOMP(KFPR(ISUBSV,1)),1).EQ.1) THEN
0586           IF(KFPR(ISUB,1).EQ.25.OR.KFPR(ISUB,1).EQ.37.OR.
0587      &    KFPR(ISUB,1).EQ.39) THEN
0588             COMFAC=COMFAC*0.5D0*ACTH0
0589           ELSE
0590             COMFAC=COMFAC*0.125D0*(3D0*ACTH0+CTNMAX**3-CTNMIN**3+
0591      &      CTPMAX**3-CTPMIN**3)
0592           ENDIF
0593         ENDIF
0594
0595 C...2 -> 2 processes: angular part of phase space integral
0596       ELSEIF(ISTSB.EQ.2.OR.ISTSB.EQ.4) THEN
0597         ACTH1=LOG((MAX(RM34,RSQM-CTNMIN)*MAX(RM34,RSQM-CTPMIN))/
0598      &  (MAX(RM34,RSQM-CTNMAX)*MAX(RM34,RSQM-CTPMAX)))
0599         ACTH2=LOG((MAX(RM34,RSQM+CTNMAX)*MAX(RM34,RSQM+CTPMAX))/
0600      &  (MAX(RM34,RSQM+CTNMIN)*MAX(RM34,RSQM+CTPMIN)))
0601         ACTH3=1D0/MAX(RM34,RSQM-CTNMAX)-1D0/MAX(RM34,RSQM-CTNMIN)+
0602      &  1D0/MAX(RM34,RSQM-CTPMAX)-1D0/MAX(RM34,RSQM-CTPMIN)
0603         ACTH4=1D0/MAX(RM34,RSQM+CTNMIN)-1D0/MAX(RM34,RSQM+CTNMAX)+
0604      &  1D0/MAX(RM34,RSQM+CTPMIN)-1D0/MAX(RM34,RSQM+CTPMAX)
0605         H3=COEF(ISUBSV,13)+
0606      &  (ACTH0/ACTH1)*COEF(ISUBSV,14)/MAX(RM34,RSQM-CTH)+
0607      &  (ACTH0/ACTH2)*COEF(ISUBSV,15)/MAX(RM34,RSQM+CTH)+
0608      &  (ACTH0/ACTH3)*COEF(ISUBSV,16)/MAX(RM34,RSQM-CTH)**2+
0609      &  (ACTH0/ACTH4)*COEF(ISUBSV,17)/MAX(RM34,RSQM+CTH)**2
0610         COMFAC=COMFAC*ACTH0*0.5D0*BE34/H3
0611
0612 C...2 -> 2 processes: take into account final state Breit-Wigners
0613         COMFAC=COMFAC*VINT(80)
0614       ENDIF
0615
0616 C...2 -> 3, 4 processes: phace space integral in tau'
0617       IF(MINT(47).GE.2.AND.ISTSB.GE.3.AND.ISTSB.LE.5) THEN
0618         ATAUP1=LOG(TAUPMX/TAUPMN)
0619         ATAUP2=((1D0-TAU/TAUPMX)**4-(1D0-TAU/TAUPMN)**4)/(4D0*TAU)
0620         H4=COEF(ISUBSV,18)+
0621      &  (ATAUP1/ATAUP2)*COEF(ISUBSV,19)*(1D0-TAU/TAUP)**3/TAUP
0622         IF(MINT(47).EQ.5) THEN
0623           ATAUP3=LOG(MAX(2D-10,1D0-TAUPMN)/MAX(2D-10,1D0-TAUPMX))
0624           H4=H4+(ATAUP1/ATAUP3)*COEF(ISUBSV,20)*TAUP/MAX(2D-10,1D0-TAUP)
0625         ELSEIF(MINT(47).GE.6) THEN
0626           ATAUP3=LOG(MAX(1D-10,1D0-TAUPMN)/MAX(1D-10,1D0-TAUPMX))
0627           H4=H4+(ATAUP1/ATAUP3)*COEF(ISUBSV,20)*TAUP/MAX(1D-10,1D0-TAUP)
0628         ENDIF
0629         COMFAC=COMFAC*ATAUP1/H4
0630       ENDIF
0631
0632 C...2 -> 3, 4 processes: effective W/Z parton distributions
0633       IF(ISTSB.EQ.3.OR.ISTSB.EQ.4) THEN
0634         IF(1D0-TAU/TAUP.GT.1D-4) THEN
0635           FZW=(1D0+TAU/TAUP)*LOG(TAUP/TAU)-2D0*(1D0-TAU/TAUP)
0636         ELSE
0637           FZW=1D0/6D0*(1D0-TAU/TAUP)**3*TAU/TAUP
0638         ENDIF
0639         COMFAC=COMFAC*FZW
0640       ENDIF
0641
0642 C...2 -> 3 processes: phase space integrals for pT1, pT2, y3, mirror
0643       IF(ISTSB.EQ.5) THEN
0644         COMFAC=COMFAC*VINT(205)*VINT(210)*VINT(212)*VINT(214)/
0645      &  (128D0*PARU(1)**4*VINT(220))*(TAU**2/TAUP)
0646       ENDIF
0647
0648 C...Phase space integral for low-pT and multiple interactions
0649       IF(ISTSB.EQ.9) THEN
0650         COMFAC=PARU(1)*PARU(5)*FACK*0.5D0*VINT(2)/SH2
0651         ATAU1=LOG(2D0*(1D0+SQRT(1D0-XT2))/XT2-1D0)
0652         ATAU2=2D0*ATAN(1D0/XT2-1D0)/SQRT(XT2)
0653         H1=COEF(ISUBSV,1)+(ATAU1/ATAU2)*COEF(ISUBSV,2)/SQRT(TAU)
0654         COMFAC=COMFAC*ATAU1/H1
0655         AYST0=YSTMAX-YSTMIN
0656         AYST1=0.5D0*(YSTMAX-YSTMIN)**2
0657         AYST3=2D0*(ATAN(EXP(YSTMAX))-ATAN(EXP(YSTMIN)))
0658         H2=(AYST0/AYST1)*COEF(ISUBSV,8)*(YST-YSTMIN)+
0659      &  (AYST0/AYST1)*COEF(ISUBSV,9)*(YSTMAX-YST)+
0660      &  (AYST0/AYST3)*COEF(ISUBSV,10)/COSH(YST)
0661         COMFAC=COMFAC*AYST0/H2
0662         IF(MSTP(82).LE.1) COMFAC=COMFAC*XT2**2*(1D0/VINT(149)-1D0)
0663 C...For MSTP(82)>=2 an additional factor (xT2/(xT2+VINT(149))**2 is
0664 C...introduced to make cross-section finite for xT2 -> 0
0665         IF(MSTP(82).GE.2) COMFAC=COMFAC*XT2**2/(VINT(149)*
0666      &  (1D0+VINT(149)))
0667       ENDIF
0668
0669 C...Real gamma + gamma: include factor 2 when different nature
0670   160 IF(MINT(11).EQ.22.AND.MINT(12).EQ.22.AND.MINT(123).GE.4.AND.
0671      &MSTP(14).LE.10) COMFAC=2D0*COMFAC
0672
0673 C...Extra factors to include the effects of
0674 C...longitudinal resolved photons (but not direct or DIS ones).
0675       DO 170 ISDE=1,2
0676         IF(MINT(10+ISDE).EQ.22.AND.MINT(106+ISDE).GE.1.AND.
0677      &  MINT(106+ISDE).LE.3) THEN
0678           VINT(314+ISDE)=1D0
0679           XY=PARP(166+ISDE)
0680           IF(MSTP(16).EQ.0) THEN
0681             IF(VINT(304+ISDE).GT.0D0.AND.VINT(304+ISDE).LT.1D0)
0682      &      XY=VINT(304+ISDE)
0683           ELSE
0684             IF(VINT(308+ISDE).GT.0D0.AND.VINT(308+ISDE).LT.1D0)
0685      &      XY=VINT(308+ISDE)
0686           ENDIF
0687           Q2GA=VINT(306+ISDE)
0688           IF(MSTP(17).GT.0.AND.XY.GT.0D0.AND.XY.LT.1D0.AND.
0689      &    Q2GA.GT.0D0) THEN
0690             REDUCE=0D0
0691             IF(MSTP(17).EQ.1) THEN
0692               REDUCE=4D0*Q2*Q2GA/(Q2+Q2GA)**2
0693             ELSEIF(MSTP(17).EQ.2) THEN
0694               REDUCE=4D0*Q2GA/(Q2+Q2GA)
0695             ELSEIF(MSTP(17).EQ.3) THEN
0696               PMVIRT=PMAS(PYCOMP(113),1)
0697               REDUCE=4D0*Q2GA/(PMVIRT**2+Q2GA)
0698             ELSEIF(MSTP(17).EQ.4.AND.MINT(106+ISDE).EQ.1) THEN
0699               PMVIRT=PMAS(PYCOMP(113),1)
0700               REDUCE=4D0*PMVIRT**2*Q2GA/(PMVIRT**2+Q2GA)**2
0701             ELSEIF(MSTP(17).EQ.4.AND.MINT(106+ISDE).EQ.2) THEN
0702               PMVIRT=PMAS(PYCOMP(113),1)
0703               REDUCE=4D0*PMVIRT**2*Q2GA/(PMVIRT**2+Q2GA)**2
0704             ELSEIF(MSTP(17).EQ.4.AND.MINT(106+ISDE).EQ.3) THEN
0705               PMVSMN=4D0*PARP(15)**2
0706               PMVSMX=4D0*VINT(154)**2
0707               REDTRA=1D0/(PMVSMN+Q2GA)-1D0/(PMVSMX+Q2GA)
0708               REDLON=(3D0*PMVSMN+Q2GA)/(PMVSMN+Q2GA)**3-
0709      &        (3D0*PMVSMX+Q2GA)/(PMVSMX+Q2GA)**3
0710               REDUCE=4D0*(Q2GA/6D0)*REDLON/REDTRA
0711             ELSEIF(MSTP(17).EQ.5.AND.MINT(106+ISDE).EQ.1) THEN
0712               PMVIRT=PMAS(PYCOMP(113),1)
0713               REDUCE=4D0*Q2GA/(PMVIRT**2+Q2GA)
0714             ELSEIF(MSTP(17).EQ.5.AND.MINT(106+ISDE).EQ.2) THEN
0715               PMVIRT=PMAS(PYCOMP(113),1)
0716               REDUCE=4D0*Q2GA/(PMVIRT**2+Q2GA)
0717             ELSEIF(MSTP(17).EQ.5.AND.MINT(106+ISDE).EQ.3) THEN
0718               PMVSMN=4D0*PARP(15)**2
0719               PMVSMX=4D0*VINT(154)**2
0720               REDTRA=1D0/(PMVSMN+Q2GA)-1D0/(PMVSMX+Q2GA)
0721               REDLON=1D0/(PMVSMN+Q2GA)**2-1D0/(PMVSMX+Q2GA)**2
0722               REDUCE=4D0*(Q2GA/2D0)*REDLON/REDTRA
0723 C ........Hermes version of R_VMD, this still needs work to
0724 C only apply it for process 91 and separate rho from the rest.
0725             ELSEIF(MSTP(17).EQ.6.AND.MINT(106+ISDE).EQ.2) THEN
0726              IF ((MINT(1).le.94).and.(MINT(1).gt.90)) then
0727               if (VINT(67).gt.0.0) then
0728                 if (VINT(67).eq.PMAS(PYCOMP(113),1)) then
0729                    REDUCE=(Q2GA/VINT(67)**2)**PARP(166)
0730                 else
0731                    REDUCE=Q2GA/VINT(67)**2
0732                 endif
0733               else
0734                 PMVIRT=PMAS(PYCOMP(113),1)
0735                 REDUCE=Q2GA/PMVIRT**2
0736               endif
0737              ENDIF
0738             ENDIF
0739             BEAMAS=PYMASS(11)
0740             IF(VINT(302+ISDE).GT.0D0) BEAMAS=VINT(302+ISDE)
0741             IF((MINT(11).EQ.22).and.
0742      &         (MINT(12).EQ.2212.or.MINT(12).EQ.2112)) THEN
0743               FRACLT=1D0/(1D0+(XY**2*(1D0-2D0*BEAMAS**2/Q2GA))/
0744      &         (2D0/(1D0+Q2GA/XY**2/ebeamEnucl**2)*(1D0-XY-
0745      &         (Q2GA/4D0/ebeamEnucl**2))))
0746             ELSE
0747               FRACLT=1D0/(1D0+XY**2/2D0/(1D0-XY)*
0748      &        (1D0-2D0*BEAMAS**2/Q2GA))
0749             ENDIF
0750             VINT(314+ISDE)=1D0+PARP(165)*REDUCE*FRACLT
0751           ENDIF
0752         ELSE
0753           VINT(314+ISDE)=1D0
0754         ENDIF
0755         COMFAC=COMFAC*VINT(314+ISDE)
0756   170 CONTINUE
0757
0758 C...Evaluate cross sections - done in separate routines by kind
0759 C...of physics, to keep PYSIGH of sensible size.
0760       IF(MAP.EQ.1) THEN
0761 C...Standard QCD (including photons).
0762         CALL PYSGQC(NCHN,SIGS)
0763       ELSEIF(MAP.EQ.2) THEN
0764 C...Heavy flavours.
0765         CALL PYSGHF(NCHN,SIGS)
0766       ELSEIF(MAP.EQ.3) THEN
0767 C...W/Z.
0768         CALL PYSGWZ(NCHN,SIGS)
0769       ELSEIF(MAP.EQ.4) THEN
0770 C...Higgs (2 doublets; including longitudinal W/Z scattering).
0771         CALL PYSGHG(NCHN,SIGS)
0772       ELSEIF(MAP.EQ.5) THEN
0773 C...SUSY.
0774         CALL PYSGSU(NCHN,SIGS)
0775       ELSEIF(MAP.EQ.6) THEN
0776 C...Technicolor.
0777         CALL PYSGTC(NCHN,SIGS)
0778       ELSEIF(MAP.EQ.7) THEN
0779 C...Exotics (Z'/W'/LQ/R/f*/H++/Z_R/W_R/G*).
0780         CALL PYSGEX(NCHN,SIGS)
0781       ENDIF
0782
0783 C...Multiply with parton distributions
0784       IF(ISUB.LE.90.OR.ISUB.GE.96) THEN
0785         DO 180 ICHN=1,NCHN
0786           IF(MINT(45).GE.2) THEN
0787             KFL1=ISIG(ICHN,1)
0788             SIGH(ICHN)=SIGH(ICHN)*XSFX(1,KFL1)
0789           ENDIF
0790           IF(MINT(46).GE.2) THEN
0791             KFL2=ISIG(ICHN,2)
0792             SIGH(ICHN)=SIGH(ICHN)*XSFX(2,KFL2)
0793           ENDIF
0794           SIGS=SIGS+SIGH(ICHN)
0795   180   CONTINUE
0796       ENDIF
0797
0798       RETURN
0799       END