hijing/src/pyhithia.f

0001
0002 C*********************************************************************
0003
0004       SUBROUTINE PYHITHIA
0005
0006 C...Administers the generation of a high-pt event via calls to a number
0007 C...of subroutines; also computes cross-sections.
0008       COMMON/LUJETS/N,K(9000,5),P(9000,5),V(9000,5)
0009       SAVE /LUJETS/
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/PYHIINT5/NGEN(0:200,3),XSEC(0:200,3)
0023       SAVE /PYHIINT5/
0024
0025 C...Loop over desired number of overlayed events (normally 1).
0026       MINT(7)=0
0027       MINT(8)=0
0028       NOVL=1
0029       IF(MSTP(131).NE.0) CALL PYHIOVLY(2)
0030       IF(MSTP(131).NE.0) NOVL=MINT(81)
0031       MINT(83)=0
0032       MINT(84)=MSTP(126)
0033       MSTU(70)=0
0034       DO 190 IOVL=1,NOVL
0035       IF(MINT(84)+100.GE.MSTU(4)) THEN
0036         CALL LUERRM(11,
0037      &  '(PYHITHIA:) no more space in LUJETS for overlayed events')
0038         IF(MSTU(21).GE.1) GOTO 200
0039       ENDIF
0040       MINT(82)=IOVL
0041
0042 C...Generate variables of hard scattering.
0043   100 CONTINUE
0044       IF(IOVL.EQ.1) NGEN(0,2)=NGEN(0,2)+1
0045       MINT(31)=0
0046       MINT(51)=0
0047       CALL PYHIRAND
0048       ISUB=MINT(1)
0049       IF(IOVL.EQ.1) THEN
0050         NGEN(ISUB,2)=NGEN(ISUB,2)+1
0051
0052 C...Store information on hard interaction.
0053         DO 110 J=1,200
0054         MSTI(J)=0
0055   110   PARI(J)=0.
0056         MSTI(1)=MINT(1)
0057         MSTI(2)=MINT(2)
0058         MSTI(11)=MINT(11)
0059         MSTI(12)=MINT(12)
0060         MSTI(15)=MINT(15)
0061         MSTI(16)=MINT(16)
0062         MSTI(17)=MINT(17)
0063         MSTI(18)=MINT(18)
0064         PARI(11)=VINT(1)
0065         PARI(12)=VINT(2)
0066         IF(ISUB.NE.95) THEN
0067           DO 120 J=13,22
0068   120     PARI(J)=VINT(30+J)
0069           PARI(33)=VINT(41)
0070           PARI(34)=VINT(42)
0071           PARI(35)=PARI(33)-PARI(34)
0072           PARI(36)=VINT(21)
0073           PARI(37)=VINT(22)
0074           PARI(38)=VINT(26)
0075           PARI(41)=VINT(23)
0076         ENDIF
0077       ENDIF
0078
0079       IF(MSTP(111).EQ.-1) GOTO 160
0080       IF(ISUB.LE.90.OR.ISUB.GE.95) THEN
0081 C...Hard scattering (including low-pT):
0082 C...reconstruct kinematics and colour flow of hard scattering.
0083         CALL PYHISCAT
0084         IF(MINT(51).EQ.1) GOTO 100
0085
0086 C...Showering of initial state partons (optional).
0087         IPU1=MINT(84)+1
0088         IPU2=MINT(84)+2
0089         IF(MSTP(61).GE.1.AND.MINT(43).NE.1.AND.ISUB.NE.95)
0090      &  CALL PYHISSPA(IPU1,IPU2)
0091         NSAV1=N
0092
0093 C...Multiple interactions.
0094         IF(MSTP(81).GE.1.AND.MINT(43).EQ.4.AND.ISUB.NE.95)
0095      &  CALL PYHIMULT(6)
0096         MINT(1)=ISUB
0097         NSAV2=N
0098
0099 C...Hadron remnants and primordial kT.
0100         CALL PYHIREMN(IPU1,IPU2)
0101         IF(MINT(51).EQ.1) GOTO 100
0102         NSAV3=N
0103
0104 C...Showering of final state partons (optional).
0105         IPU3=MINT(84)+3
0106         IPU4=MINT(84)+4
0107         IF(MSTP(71).GE.1.AND.ISUB.NE.95.AND.K(IPU3,1).GT.0.AND.
0108      &  K(IPU3,1).LE.10.AND.K(IPU4,1).GT.0.AND.K(IPU4,1).LE.10) THEN
0109           QMAX=SQRT(PARP(71)*VINT(52))
0110           IF(ISUB.EQ.5) QMAX=SQRT(PMAS(23,1)**2)
0111           IF(ISUB.EQ.8) QMAX=SQRT(PMAS(24,1)**2)
0112           CALL LUSHOW(IPU3,IPU4,QMAX)
0113         ENDIF
0114
0115 C...Sum up transverse and longitudinal momenta.
0116         IF(IOVL.EQ.1) THEN
0117           PARI(65)=2.*PARI(17)
0118           DO 130 I=MSTP(126)+1,N
0119           IF(K(I,1).LE.0.OR.K(I,1).GT.10) GOTO 130
0120           PT=SQRT(P(I,1)**2+P(I,2)**2)
0121           PARI(69)=PARI(69)+PT
0122           IF(I.LE.NSAV1.OR.I.GT.NSAV3) PARI(66)=PARI(66)+PT
0123           IF(I.GT.NSAV1.AND.I.LE.NSAV2) PARI(68)=PARI(68)+PT
0124   130     CONTINUE
0125           PARI(67)=PARI(68)
0126           PARI(71)=VINT(151)
0127           PARI(72)=VINT(152)
0128           PARI(73)=VINT(151)
0129           PARI(74)=VINT(152)
0130         ENDIF
0131
0132 C...Decay of final state resonances.
0133         IF(MSTP(41).GE.1.AND.ISUB.NE.95) CALL PYHIRESD
0134
0135       ELSE
0136 C...Diffractive and elastic scattering.
0137         CALL PYHIDIFF
0138         IF(IOVL.EQ.1) THEN
0139           PARI(65)=2.*PARI(17)
0140           PARI(66)=PARI(65)
0141           PARI(69)=PARI(65)
0142         ENDIF
0143       ENDIF
0144
0145 C...Recalculate energies from momenta and masses (if desired).
0146       IF(MSTP(113).GE.1) THEN
0147         DO 140 I=MINT(83)+1,N
0148   140   IF(K(I,1).GT.0.AND.K(I,1).LE.10) P(I,4)=SQRT(P(I,1)**2+
0149      &  P(I,2)**2+P(I,3)**2+P(I,5)**2)
0150       ENDIF
0151
0152 C...Rearrange partons along strings, check invariant mass cuts.
0153       MSTU(28)=0
0154       CALL LUPREP(MINT(84)+1)
0155       IF(MSTP(112).EQ.1.AND.MSTU(28).EQ.3) GOTO 100
0156       IF(MSTP(125).EQ.0.OR.MSTP(125).EQ.1) THEN
0157         DO 150 I=MINT(84)+1,N
0158         IF(K(I,2).NE.94) GOTO 150
0159         K(I+1,3)=MOD(K(I+1,4)/MSTU(5),MSTU(5))
0160         K(I+2,3)=MOD(K(I+2,4)/MSTU(5),MSTU(5))
0161   150   CONTINUE
0162         CALL LUEDIT(12)
0163         CALL LUEDIT(14)
0164         IF(MSTP(125).EQ.0) CALL LUEDIT(15)
0165         IF(MSTP(125).EQ.0) MINT(4)=0
0166       ENDIF
0167
0168 C...Introduce separators between sections in LULIST event listing.
0169       IF(IOVL.EQ.1.AND.MSTP(125).LE.0) THEN
0170         MSTU(70)=1
0171         MSTU(71)=N
0172       ELSEIF(IOVL.EQ.1) THEN
0173         MSTU(70)=3
0174         MSTU(71)=2
0175         MSTU(72)=MINT(4)
0176         MSTU(73)=N
0177       ENDIF
0178
0179 C...Perform hadronization (if desired).
0180       IF(MSTP(111).GE.1) CALL LUEXEC
0181       IF(MSTP(125).EQ.0.OR.MSTP(125).EQ.1) CALL LUEDIT(14)
0182
0183 C...Calculate Monte Carlo estimates of cross-sections.
0184   160 IF(IOVL.EQ.1) THEN
0185         IF(MSTP(111).NE.-1) NGEN(ISUB,3)=NGEN(ISUB,3)+1
0186         NGEN(0,3)=NGEN(0,3)+1
0187         XSEC(0,3)=0.
0188         DO 170 I=1,200
0189         IF(I.EQ.96) THEN
0190           XSEC(I,3)=0.
0191         ELSEIF(MSUB(95).EQ.1.AND.(I.EQ.11.OR.I.EQ.12.OR.I.EQ.13.OR.
0192      &  I.EQ.28.OR.I.EQ.53.OR.I.EQ.68)) THEN
0193           XSEC(I,3)=XSEC(96,2)*NGEN(I,3)/MAX(1.,FLOAT(NGEN(96,1))*
0194      &    FLOAT(NGEN(96,2)))
0195         ELSEIF(NGEN(I,1).EQ.0) THEN
0196           XSEC(I,3)=0.
0197         ELSEIF(NGEN(I,2).EQ.0) THEN
0198           XSEC(I,3)=XSEC(I,2)*NGEN(0,3)/(FLOAT(NGEN(I,1))*
0199      &    FLOAT(NGEN(0,2)))
0200         ELSE
0201           XSEC(I,3)=XSEC(I,2)*NGEN(I,3)/(FLOAT(NGEN(I,1))*
0202      &    FLOAT(NGEN(I,2)))
0203         ENDIF
0204   170   XSEC(0,3)=XSEC(0,3)+XSEC(I,3)
0205         IF(MSUB(95).EQ.1) THEN
0206           NGENS=NGEN(91,3)+NGEN(92,3)+NGEN(93,3)+NGEN(94,3)+NGEN(95,3)
0207           XSECS=XSEC(91,3)+XSEC(92,3)+XSEC(93,3)+XSEC(94,3)+XSEC(95,3)
0208           XMAXS=XSEC(95,1)
0209           IF(MSUB(91).EQ.1) XMAXS=XMAXS+XSEC(91,1)
0210           IF(MSUB(92).EQ.1) XMAXS=XMAXS+XSEC(92,1)
0211           IF(MSUB(93).EQ.1) XMAXS=XMAXS+XSEC(93,1)
0212           IF(MSUB(94).EQ.1) XMAXS=XMAXS+XSEC(94,1)
0213           FAC=1.
0214           IF(NGENS.LT.NGEN(0,3)) FAC=(XMAXS-XSECS)/(XSEC(0,3)-XSECS)
0215           XSEC(11,3)=FAC*XSEC(11,3)
0216           XSEC(12,3)=FAC*XSEC(12,3)
0217           XSEC(13,3)=FAC*XSEC(13,3)
0218           XSEC(28,3)=FAC*XSEC(28,3)
0219           XSEC(53,3)=FAC*XSEC(53,3)
0220           XSEC(68,3)=FAC*XSEC(68,3)
0221           XSEC(0,3)=XSEC(91,3)+XSEC(92,3)+XSEC(93,3)+XSEC(94,3)+
0222      &    XSEC(95,1)
0223         ENDIF
0224
0225 C...Store final information.
0226         MINT(5)=MINT(5)+1
0227         MSTI(3)=MINT(3)
0228         MSTI(4)=MINT(4)
0229         MSTI(5)=MINT(5)
0230         MSTI(6)=MINT(6)
0231         MSTI(7)=MINT(7)
0232         MSTI(8)=MINT(8)
0233         MSTI(13)=MINT(13)
0234         MSTI(14)=MINT(14)
0235         MSTI(21)=MINT(21)
0236         MSTI(22)=MINT(22)
0237         MSTI(23)=MINT(23)
0238         MSTI(24)=MINT(24)
0239         MSTI(25)=MINT(25)
0240         MSTI(26)=MINT(26)
0241         MSTI(31)=MINT(31)
0242         PARI(1)=XSEC(0,3)
0243         PARI(2)=XSEC(0,3)/MINT(5)
0244         PARI(31)=VINT(141)
0245         PARI(32)=VINT(142)
0246         IF(ISUB.NE.95.AND.MINT(7)*MINT(8).NE.0) THEN
0247           PARI(42)=2.*VINT(47)/VINT(1)
0248           DO 180 IS=7,8
0249           PARI(36+IS)=P(MINT(IS),3)/VINT(1)
0250           PARI(38+IS)=P(MINT(IS),4)/VINT(1)
0251           I=MINT(IS)
0252           PR=MAX(1E-20,P(I,5)**2+P(I,1)**2+P(I,2)**2)
0253           PARI(40+IS)=SIGN(LOG(MIN((SQRT(PR+P(I,3)**2)+ABS(P(I,3)))/
0254      &    SQRT(PR),1E20)),P(I,3))
0255           PR=MAX(1E-20,P(I,1)**2+P(I,2)**2)
0256           PARI(42+IS)=SIGN(LOG(MIN((SQRT(PR+P(I,3)**2)+ABS(P(I,3)))/
0257      &    SQRT(PR),1E20)),P(I,3))
0258           PARI(44+IS)=P(I,3)/SQRT(P(I,1)**2+P(I,2)**2+P(I,3)**2)
0259           PARI(46+IS)=ULANGL(P(I,3),SQRT(P(I,1)**2+P(I,2)**2))
0260           PARI(48+IS)=ULANGL(P(I,1),P(I,2))
0261   180     CONTINUE
0262         ENDIF
0263         PARI(61)=VINT(148)
0264         IF(ISET(ISUB).EQ.1.OR.ISET(ISUB).EQ.3) THEN
0265           MSTU(161)=MINT(21)
0266           MSTU(162)=0
0267         ELSE
0268           MSTU(161)=MINT(21)
0269           MSTU(162)=MINT(22)
0270         ENDIF
0271       ENDIF
0272
0273 C...Prepare to go to next overlayed event.
0274       MSTI(41)=IOVL
0275       IF(IOVL.GE.2.AND.IOVL.LE.10) MSTI(40+IOVL)=ISUB
0276       IF(MSTU(70).LT.10) THEN
0277         MSTU(70)=MSTU(70)+1
0278         MSTU(70+MSTU(70))=N
0279       ENDIF
0280       MINT(83)=N
0281       MINT(84)=N+MSTP(126)
0282   190 CONTINUE
0283
0284 C...Information on overlayed events.
0285       IF(MSTP(131).EQ.1.AND.MSTP(133).GE.1) THEN
0286         PARI(91)=VINT(132)
0287         PARI(92)=VINT(133)
0288         PARI(93)=VINT(134)
0289         IF(MSTP(133).EQ.2) PARI(93)=PARI(93)*XSEC(0,3)/VINT(131)
0290       ENDIF
0291
0292 C...Transform to the desired coordinate frame.
0293   200 CALL PYHIFRAM(MSTP(124))
0294
0295       RETURN
0296       END