hijing/src/pyhirand.f

0001
0002 C*********************************************************************
0003
0004       SUBROUTINE PYHIRAND
0005
0006 C...Generates quantities characterizing the high-pT scattering at the
0007 C...parton level according to the matrix elements. Chooses incoming,
0008 C...reacting partons, their momentum fractions and one of the possible
0009 C...subprocesses.
0010       COMMON/LUDAT1/MSTU(200),PARU(200),MSTJ(200),PARJ(200)
0011       SAVE /LUDAT1/
0012       COMMON/LUDAT2/KCHG(500,3),PMAS(500,4),PARF(2000),VCKM(4,4)
0013       SAVE /LUDAT2/
0014       COMMON/PYHISUBS/MSEL,MSUB(200),KFIN(2,-40:40),CKIN(200)
0015       SAVE /PYHISUBS/
0016       COMMON/PYHIPARS/MSTP(200),PARP(200),MSTI(200),PARI(200)
0017       SAVE /PYHIPARS/
0018       COMMON/PYHIINT1/MINT(400),VINT(400)
0019       SAVE /PYHIINT1/
0020       COMMON/PYHIINT2/ISET(200),KFPR(200,2),COEF(200,20),ICOL(40,4,2)
0021       SAVE /PYHIINT2/
0022       COMMON/PYHIINT3/XSFX(2,-40:40),ISIG(1000,3),SIGH(1000)
0023       SAVE /PYHIINT3/
0024       COMMON/PYHIINT4/WIDP(21:40,0:40),WIDE(21:40,0:40),WIDS(21:40,3)
0025       SAVE /PYHIINT4/
0026       COMMON/PYHIINT5/NGEN(0:200,3),XSEC(0:200,3)
0027       SAVE /PYHIINT5/
0028
0029 C...Initial values, specifically for (first) semihard interaction.
0030       MINT(17)=0
0031       MINT(18)=0
0032       VINT(143)=1.
0033       VINT(144)=1.
0034       IF(MSUB(95).EQ.1.OR.MINT(82).GE.2) CALL PYHIMULT(2)
0035       ISUB=0
0036   100 MINT(51)=0
0037
0038 C...Choice of process type - first event of overlay.
0039       IF(MINT(82).EQ.1.AND.(ISUB.LE.90.OR.ISUB.GT.96)) THEN
0040         RSUB=XSEC(0,1)*RLU(0)
0041         DO 110 I=1,200
0042         IF(MSUB(I).NE.1) GOTO 110
0043         ISUB=I
0044         RSUB=RSUB-XSEC(I,1)
0045         IF(RSUB.LE.0.) GOTO 120
0046   110   CONTINUE
0047   120   IF(ISUB.EQ.95) ISUB=96
0048
0049 C...Choice of inclusive process type - overlayed events.
0050       ELSEIF(MINT(82).GE.2.AND.ISUB.EQ.0) THEN
0051         RSUB=VINT(131)*RLU(0)
0052         ISUB=96
0053         IF(RSUB.GT.VINT(106)) ISUB=93
0054         IF(RSUB.GT.VINT(106)+VINT(104)) ISUB=92
0055         IF(RSUB.GT.VINT(106)+VINT(104)+VINT(103)) ISUB=91
0056       ENDIF
0057       IF(MINT(82).EQ.1) NGEN(0,1)=NGEN(0,1)+1
0058       IF(MINT(82).EQ.1) NGEN(ISUB,1)=NGEN(ISUB,1)+1
0059       MINT(1)=ISUB
0060
0061 C...Find resonances (explicit or implicit in cross-section).
0062       MINT(72)=0
0063       KFR1=0
0064       IF(ISET(ISUB).EQ.1.OR.ISET(ISUB).EQ.3) THEN
0065         KFR1=KFPR(ISUB,1)
0066       ELSEIF(ISUB.GE.71.AND.ISUB.LE.77) THEN
0067         KFR1=25
0068       ENDIF
0069       IF(KFR1.NE.0) THEN
0070         TAUR1=PMAS(KFR1,1)**2/VINT(2)
0071         GAMR1=PMAS(KFR1,1)*PMAS(KFR1,2)/VINT(2)
0072         MINT(72)=1
0073         MINT(73)=KFR1
0074         VINT(73)=TAUR1
0075         VINT(74)=GAMR1
0076       ENDIF
0077       IF(ISUB.EQ.141) THEN
0078         KFR2=23
0079         TAUR2=PMAS(KFR2,1)**2/VINT(2)
0080         GAMR2=PMAS(KFR2,1)*PMAS(KFR2,2)/VINT(2)
0081         MINT(72)=2
0082         MINT(74)=KFR2
0083         VINT(75)=TAUR2
0084         VINT(76)=GAMR2
0085       ENDIF
0086
0087 C...Find product masses and minimum pT of process,
0088 C...optionally with broadening according to a truncated Breit-Wigner.
0089       VINT(63)=0.
0090       VINT(64)=0.
0091       MINT(71)=0
0092       VINT(71)=CKIN(3)
0093       IF(MINT(82).GE.2) VINT(71)=0.
0094       IF(ISET(ISUB).EQ.2.OR.ISET(ISUB).EQ.4) THEN
0095         DO 130 I=1,2
0096         IF(KFPR(ISUB,I).EQ.0) THEN
0097         ELSEIF(MSTP(42).LE.0) THEN
0098           VINT(62+I)=PMAS(KFPR(ISUB,I),1)**2
0099         ELSE
0100           VINT(62+I)=ULMASS(KFPR(ISUB,I))**2
0101         ENDIF
0102   130   CONTINUE
0103         IF(MIN(VINT(63),VINT(64)).LT.CKIN(6)**2) MINT(71)=1
0104         IF(MINT(71).EQ.1) VINT(71)=MAX(CKIN(3),CKIN(5))
0105       ENDIF
0106
0107       IF(ISET(ISUB).EQ.0) THEN
0108 C...Double or single diffractive, or elastic scattering:
0109 C...choose m^2 according to 1/m^2 (diffractive), constant (elastic)
0110         IS=INT(1.5+RLU(0))
0111         VINT(63)=VINT(3)**2
0112         VINT(64)=VINT(4)**2
0113         IF(ISUB.EQ.92.OR.ISUB.EQ.93) VINT(62+IS)=PARP(111)**2
0114         IF(ISUB.EQ.93) VINT(65-IS)=PARP(111)**2
0115         SH=VINT(2)
0116         SQM1=VINT(3)**2
0117         SQM2=VINT(4)**2
0118         SQM3=VINT(63)
0119         SQM4=VINT(64)
0120         SQLA12=(SH-SQM1-SQM2)**2-4.*SQM1*SQM2
0121         SQLA34=(SH-SQM3-SQM4)**2-4.*SQM3*SQM4
0122         THTER1=SQM1+SQM2+SQM3+SQM4-(SQM1-SQM2)*(SQM3-SQM4)/SH-SH
0123         THTER2=SQRT(MAX(0.,SQLA12))*SQRT(MAX(0.,SQLA34))/SH
0124         THL=0.5*(THTER1-THTER2)
0125         THU=0.5*(THTER1+THTER2)
0126         THM=MIN(MAX(THL,PARP(101)),THU)
0127         JTMAX=0
0128         IF(ISUB.EQ.92.OR.ISUB.EQ.93) JTMAX=ISUB-91
0129         DO 140 JT=1,JTMAX
0130         MINT(13+3*JT-IS*(2*JT-3))=1
0131         SQMMIN=VINT(59+3*JT-IS*(2*JT-3))
0132         SQMI=VINT(8-3*JT+IS*(2*JT-3))**2
0133         SQMJ=VINT(3*JT-1-IS*(2*JT-3))**2
0134         SQMF=VINT(68-3*JT+IS*(2*JT-3))
0135         SQUA=0.5*SH/SQMI*((1.+(SQMI-SQMJ)/SH)*THM+SQMI-SQMF-
0136      &  SQMJ**2/SH+(SQMI+SQMJ)*SQMF/SH+(SQMI-SQMJ)**2/SH**2*SQMF)
0137         QUAR=SH/SQMI*(THM*(THM+SH-SQMI-SQMJ-SQMF*(1.-(SQMI-SQMJ)/SH))+
0138      &  SQMI*SQMJ-SQMJ*SQMF*(1.+(SQMI-SQMJ-SQMF)/SH))
0139         SQMMAX=SQUA+SQRT(MAX(0.,SQUA**2-QUAR))
0140         IF(ABS(QUAR/SQUA**2).LT.1.E-06) SQMMAX=0.5*QUAR/SQUA
0141         SQMMAX=MIN(SQMMAX,(VINT(1)-SQRT(SQMF))**2)
0142         VINT(59+3*JT-IS*(2*JT-3))=SQMMIN*(SQMMAX/SQMMIN)**RLU(0)
0143   140   CONTINUE
0144 C...Choose t-hat according to exp(B*t-hat+C*t-hat^2).
0145         SQM3=VINT(63)
0146         SQM4=VINT(64)
0147         SQLA34=(SH-SQM3-SQM4)**2-4.*SQM3*SQM4
0148         THTER1=SQM1+SQM2+SQM3+SQM4-(SQM1-SQM2)*(SQM3-SQM4)/SH-SH
0149         THTER2=SQRT(MAX(0.,SQLA12))*SQRT(MAX(0.,SQLA34))/SH
0150         THL=0.5*(THTER1-THTER2)
0151         THU=0.5*(THTER1+THTER2)
0152         B=VINT(121)
0153         C=VINT(122)
0154         IF(ISUB.EQ.92.OR.ISUB.EQ.93) THEN
0155           B=0.5*B
0156           C=0.5*C
0157         ENDIF
0158         THM=MIN(MAX(THL,PARP(101)),THU)
0159         EXPTH=0.
0160         THARG=B*(THM-THU)
0161         IF(THARG.GT.-20.) EXPTH=EXP(THARG)
0162   150   TH=THU+LOG(EXPTH+(1.-EXPTH)*RLU(0))/B
0163         TH=MAX(THM,MIN(THU,TH))
0164         RATLOG=MIN((B+C*(TH+THM))*(TH-THM),(B+C*(TH+THU))*(TH-THU))
0165         IF(RATLOG.LT.LOG(RLU(0))) GOTO 150
0166         VINT(21)=1.
0167         VINT(22)=0.
0168         VINT(23)=MIN(1.,MAX(-1.,(2.*TH-THTER1)/THTER2))
0169
0170 C...Note: in the following, by In is meant the integral over the
0171 C...quantity multiplying coefficient cn.
0172 C...Choose tau according to h1(tau)/tau, where
0173 C...h1(tau) = c0 + I0/I1*c1*1/tau + I0/I2*c2*1/(tau+tau_R) +
0174 C...I0/I3*c3*tau/((s*tau-m^2)^2+(m*Gamma)^2) +
0175 C...I0/I4*c4*1/(tau+tau_R') +
0176 C...I0/I5*c5*tau/((s*tau-m'^2)^2+(m'*Gamma')^2), and
0177 C...c0 + c1 + c2 + c3 + c4 + c5 = 1
0178       ELSEIF(ISET(ISUB).GE.1.AND.ISET(ISUB).LE.4) THEN
0179         CALL PYHIKLIM(1)
0180         IF(MINT(51).NE.0) GOTO 100
0181         RTAU=RLU(0)
0182         MTAU=1
0183         IF(RTAU.GT.COEF(ISUB,1)) MTAU=2
0184         IF(RTAU.GT.COEF(ISUB,1)+COEF(ISUB,2)) MTAU=3
0185         IF(RTAU.GT.COEF(ISUB,1)+COEF(ISUB,2)+COEF(ISUB,3)) MTAU=4
0186         IF(RTAU.GT.COEF(ISUB,1)+COEF(ISUB,2)+COEF(ISUB,3)+COEF(ISUB,4))
0187      &  MTAU=5
0188         IF(RTAU.GT.COEF(ISUB,1)+COEF(ISUB,2)+COEF(ISUB,3)+COEF(ISUB,4)+
0189      &  COEF(ISUB,5)) MTAU=6
0190         CALL PYHIKMAP(1,MTAU,RLU(0))
0191
0192 C...2 -> 3, 4 processes:
0193 C...Choose tau' according to h4(tau,tau')/tau', where
0194 C...h4(tau,tau') = c0 + I0/I1*c1*(1 - tau/tau')^3/tau', and
0195 C...c0 + c1 = 1.
0196         IF(ISET(ISUB).EQ.3.OR.ISET(ISUB).EQ.4) THEN
0197           CALL PYHIKLIM(4)
0198           IF(MINT(51).NE.0) GOTO 100
0199           RTAUP=RLU(0)
0200           MTAUP=1
0201           IF(RTAUP.GT.COEF(ISUB,15)) MTAUP=2
0202           CALL PYHIKMAP(4,MTAUP,RLU(0))
0203         ENDIF
0204
0205 C...Choose y* according to h2(y*), where
0206 C...h2(y*) = I0/I1*c1*(y*-y*min) + I0/I2*c2*(y*max-y*) +
0207 C...I0/I3*c3*1/cosh(y*), I0 = y*max-y*min, and c1 + c2 + c3 = 1.
0208         CALL PYHIKLIM(2)
0209         IF(MINT(51).NE.0) GOTO 100
0210         RYST=RLU(0)
0211         MYST=1
0212         IF(RYST.GT.COEF(ISUB,7)) MYST=2
0213         IF(RYST.GT.COEF(ISUB,7)+COEF(ISUB,8)) MYST=3
0214         CALL PYHIKMAP(2,MYST,RLU(0))
0215
0216 C...2 -> 2 processes:
0217 C...Choose cos(theta-hat) (cth) according to h3(cth), where
0218 C...h3(cth) = c0 + I0/I1*c1*1/(A - cth) + I0/I2*c2*1/(A + cth) +
0219 C...I0/I3*c3*1/(A - cth)^2 + I0/I4*c4*1/(A + cth)^2,
0220 C...A = 1 + 2*(m3*m4/sh)^2 (= 1 for massless products),
0221 C...and c0 + c1 + c2 + c3 + c4 = 1.
0222         CALL PYHIKLIM(3)
0223         IF(MINT(51).NE.0) GOTO 100
0224         IF(ISET(ISUB).EQ.2.OR.ISET(ISUB).EQ.4) THEN
0225           RCTH=RLU(0)
0226           MCTH=1
0227           IF(RCTH.GT.COEF(ISUB,10)) MCTH=2
0228           IF(RCTH.GT.COEF(ISUB,10)+COEF(ISUB,11)) MCTH=3
0229           IF(RCTH.GT.COEF(ISUB,10)+COEF(ISUB,11)+COEF(ISUB,12)) MCTH=4
0230           IF(RCTH.GT.COEF(ISUB,10)+COEF(ISUB,11)+COEF(ISUB,12)+
0231      &    COEF(ISUB,13)) MCTH=5
0232           CALL PYHIKMAP(3,MCTH,RLU(0))
0233         ENDIF
0234
0235 C...Low-pT or multiple interactions (first semihard interaction).
0236       ELSEIF(ISET(ISUB).EQ.5) THEN
0237         CALL PYHIMULT(3)
0238         ISUB=MINT(1)
0239       ENDIF
0240
0241 C...Choose azimuthal angle.
0242       VINT(24)=PARU(2)*RLU(0)
0243
0244 C...Check against user cuts on kinematics at parton level.
0245       MINT(51)=0
0246       IF(ISUB.LE.90.OR.ISUB.GT.100) CALL PYHIKLIM(0)
0247       IF(MINT(51).NE.0) GOTO 100
0248       IF(MINT(82).EQ.1.AND.MSTP(141).GE.1) THEN
0249         MCUT=0
0250         IF(MSUB(91)+MSUB(92)+MSUB(93)+MSUB(94)+MSUB(95).EQ.0)
0251      &  CALL PYHIKCUT(MCUT)
0252         IF(MCUT.NE.0) GOTO 100
0253       ENDIF
0254
0255 C...Calculate differential cross-section for different subprocesses.
0256       CALL PYHISIGH(NCHN,SIGS)
0257
0258 C...Calculations for Monte Carlo estimate of all cross-sections.
0259       IF(MINT(82).EQ.1.AND.ISUB.LE.90.OR.ISUB.GE.96) THEN
0260         XSEC(ISUB,2)=XSEC(ISUB,2)+SIGS
0261       ELSEIF(MINT(82).EQ.1) THEN
0262         XSEC(ISUB,2)=XSEC(ISUB,2)+XSEC(ISUB,1)
0263       ENDIF
0264
0265 C...Multiple interactions: store results of cross-section calculation.
0266       IF(MINT(43).EQ.4.AND.MSTP(82).GE.3) THEN
0267         VINT(153)=SIGS
0268         CALL PYHIMULT(4)
0269       ENDIF
0270
0271 C...Weighting using estimate of maximum of differential cross-section.
0272       VIOL=SIGS/XSEC(ISUB,1)
0273       IF(VIOL.LT.RLU(0)) GOTO 100
0274
0275 C...Check for possible violation of estimated maximum of differential
0276 C...cross-section used in weighting.
0277       IF(MSTP(123).LE.0) THEN
0278         IF(VIOL.GT.1.) THEN
0279           WRITE(MSTU(11),1000) VIOL,NGEN(0,3)+1
0280           WRITE(MSTU(11),1100) ISUB,VINT(21),VINT(22),VINT(23),VINT(26)
0281           STOP
0282         ENDIF
0283       ELSEIF(MSTP(123).EQ.1) THEN
0284         IF(VIOL.GT.VINT(108)) THEN
0285           VINT(108)=VIOL
0286 C          IF(VIOL.GT.1.) THEN
0287 C            WRITE(MSTU(11),1200) VIOL,NGEN(0,3)+1
0288 C            WRITE(MSTU(11),1100) ISUB,VINT(21),VINT(22),VINT(23),
0289 C     &      VINT(26)
0290 C          ENDIF
0291         ENDIF
0292       ELSEIF(VIOL.GT.VINT(108)) THEN
0293         VINT(108)=VIOL
0294         IF(VIOL.GT.1.) THEN
0295           XDIF=XSEC(ISUB,1)*(VIOL-1.)
0296           XSEC(ISUB,1)=XSEC(ISUB,1)+XDIF
0297           IF(MSUB(ISUB).EQ.1.AND.(ISUB.LE.90.OR.ISUB.GT.96))
0298      &    XSEC(0,1)=XSEC(0,1)+XDIF
0299 C          WRITE(MSTU(11),1200) VIOL,NGEN(0,3)+1
0300 C          WRITE(MSTU(11),1100) ISUB,VINT(21),VINT(22),VINT(23),VINT(26)
0301 C          IF(ISUB.LE.9) THEN
0302 C            WRITE(MSTU(11),1300) ISUB,XSEC(ISUB,1)
0303 C          ELSEIF(ISUB.LE.99) THEN
0304 C            WRITE(MSTU(11),1400) ISUB,XSEC(ISUB,1)
0305 C          ELSE
0306 C            WRITE(MSTU(11),1500) ISUB,XSEC(ISUB,1)
0307 C          ENDIF
0308           VINT(108)=1.
0309         ENDIF
0310       ENDIF
0311
0312 C...Multiple interactions: choose impact parameter.
0313       VINT(148)=1.
0314       IF(MINT(43).EQ.4.AND.(ISUB.LE.90.OR.ISUB.GE.96).AND.MSTP(82).GE.3)
0315      &THEN
0316         CALL PYHIMULT(5)
0317         IF(VINT(150).LT.RLU(0)) GOTO 100
0318       ENDIF
0319       IF(MINT(82).EQ.1.AND.MSUB(95).EQ.1) THEN
0320         IF(ISUB.LE.90.OR.ISUB.GE.95) NGEN(95,1)=NGEN(95,1)+1
0321         IF(ISUB.LE.90.OR.ISUB.GE.96) NGEN(96,2)=NGEN(96,2)+1
0322       ENDIF
0323       IF(ISUB.LE.90.OR.ISUB.GE.96) MINT(31)=MINT(31)+1
0324
0325 C...Choose flavour of reacting partons (and subprocess).
0326       RSIGS=SIGS*RLU(0)
0327       QT2=VINT(48)
0328       RQQBAR=PARP(87)*(1.-(QT2/(QT2+(PARP(88)*PARP(82))**2))**2)
0329       IF(ISUB.NE.95.AND.(ISUB.NE.96.OR.MSTP(82).LE.1.OR.
0330      &RLU(0).GT.RQQBAR)) THEN
0331         DO 190 ICHN=1,NCHN
0332         KFL1=ISIG(ICHN,1)
0333         KFL2=ISIG(ICHN,2)
0334         MINT(2)=ISIG(ICHN,3)
0335         RSIGS=RSIGS-SIGH(ICHN)
0336         IF(RSIGS.LE.0.) GOTO 210
0337   190   CONTINUE
0338
0339 C...Multiple interactions: choose qqbar preferentially at small pT.
0340       ELSEIF(ISUB.EQ.96) THEN
0341         CALL PYHISPLI(MINT(11),21,KFL1,KFLDUM)
0342         CALL PYHISPLI(MINT(12),21,KFL2,KFLDUM)
0343         MINT(1)=11
0344         MINT(2)=1
0345         IF(KFL1.EQ.KFL2.AND.RLU(0).LT.0.5) MINT(2)=2
0346
0347 C...Low-pT: choose string drawing configuration.
0348       ELSE
0349         KFL1=21
0350         KFL2=21
0351         RSIGS=6.*RLU(0)
0352         MINT(2)=1
0353         IF(RSIGS.GT.1.) MINT(2)=2
0354         IF(RSIGS.GT.2.) MINT(2)=3
0355       ENDIF
0356
0357 C...Reassign QCD process. Partons before initial state radiation.
0358   210 IF(MINT(2).GT.10) THEN
0359         MINT(1)=MINT(2)/10
0360         MINT(2)=MOD(MINT(2),10)
0361       ENDIF
0362       MINT(15)=KFL1
0363       MINT(16)=KFL2
0364       MINT(13)=MINT(15)
0365       MINT(14)=MINT(16)
0366       VINT(141)=VINT(41)
0367       VINT(142)=VINT(42)
0368
0369 C...Format statements for differential cross-section maximum violations.
0370  1000 FORMAT(1X,'Error: maximum violated by',1P,E11.3,1X,
0371      &'in event',1X,I7,'.'/1X,'Execution stopped!')
0372  1100 FORMAT(1X,'ISUB = ',I3,'; Point of violation:'/1X,'tau =',1P,
0373      &E11.3,', y* =',E11.3,', cthe = ',0P,F11.7,', tau'' =',1P,E11.3)
0374  1200 FORMAT(1X,'Warning: maximum violated by',1P,E11.3,1X,
0375      &'in event',1X,I7)
0376  1300 FORMAT(1X,'XSEC(',I1,',1) increased to',1P,E11.3)
0377  1400 FORMAT(1X,'XSEC(',I2,',1) increased to',1P,E11.3)
0378  1500 FORMAT(1X,'XSEC(',I3,',1) increased to',1P,E11.3)
0379
0380       RETURN
0381       END