Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2023-10-25 09:48:25 

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

This page was automatically generated by the 2.2.1 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!