HIJING/hipyset1_35/pythia_hijing.F

   1 * $Id$
   2
   3 C*********************************************************************
   4
   5       SUBROUTINE PYTHIA_HIJING
   6
   7 C...Administers the generation of a high-pt event via calls to a number
   8 C...of subroutines; also computes cross-sections.
   9 #include "lujets_hijing.inc"
  10 #include "ludat1_hijing.inc"
  11 #include "ludat2_hijing.inc"
  12 #include "pysubs_hijing.inc"
  13 #include "pypars_hijing.inc"
  14 #include "pyint1_hijing.inc"
  15 #include "pyint2_hijing.inc"
  16 #include "pyint5_hijing.inc"
  17
  18 C...Loop over desired number of overlayed events (normally 1).
  19       MINT(7)=0
  20       MINT(8)=0
  21       NOVL=1
  22       IF(MSTP(131).NE.0) CALL PYOVLY_HIJING(2)
  23       IF(MSTP(131).NE.0) NOVL=MINT(81)
  24       MINT(83)=0
  25       MINT(84)=MSTP(126)
  26       MSTU(70)=0
  27       DO 190 IOVL=1,NOVL
  28       IF(MINT(84)+100.GE.MSTU(4)) THEN
  29         CALL LUERRM_HIJING(11,
  30      &        '(PYTHIA_HIJING:) no more space in LUJETS_HIJING '/
  31      $        /'for overlayed events')
  32         IF(MSTU(21).GE.1) GOTO 200
  33       ENDIF
  34       MINT(82)=IOVL
  35
  36 C...Generate variables of hard scattering.
  37   100 CONTINUE
  38       IF(IOVL.EQ.1) NGEN(0,2)=NGEN(0,2)+1
  39       MINT(31)=0
  40       MINT(51)=0
  41       CALL PYRAND_HIJING
  42       ISUB=MINT(1)
  43       IF(IOVL.EQ.1) THEN
  44         NGEN(ISUB,2)=NGEN(ISUB,2)+1
  45
  46 C...Store information on hard interaction.
  47         DO 110 J=1,200
  48         MSTI(J)=0
  49   110   PARI(J)=0.
  50         MSTI(1)=MINT(1)
  51         MSTI(2)=MINT(2)
  52         MSTI(11)=MINT(11)
  53         MSTI(12)=MINT(12)
  54         MSTI(15)=MINT(15)
  55         MSTI(16)=MINT(16)
  56         MSTI(17)=MINT(17)
  57         MSTI(18)=MINT(18)
  58         PARI(11)=VINT(1)
  59         PARI(12)=VINT(2)
  60         IF(ISUB.NE.95) THEN
  61           DO 120 J=13,22
  62   120     PARI(J)=VINT(30+J)
  63           PARI(33)=VINT(41)
  64           PARI(34)=VINT(42)
  65           PARI(35)=PARI(33)-PARI(34)
  66           PARI(36)=VINT(21)
  67           PARI(37)=VINT(22)
  68           PARI(38)=VINT(26)
  69           PARI(41)=VINT(23)
  70         ENDIF
  71       ENDIF
  72
  73       IF(MSTP(111).EQ.-1) GOTO 160
  74       IF(ISUB.LE.90.OR.ISUB.GE.95) THEN
  75 C...Hard scattering (including low-pT):
  76 C...reconstruct kinematics and colour flow of hard scattering.
  77         CALL PYSCAT_HIJING
  78         IF(MINT(51).EQ.1) GOTO 100
  79
  80 C...Showering of initial state partons (optional).
  81         IPU1=MINT(84)+1
  82         IPU2=MINT(84)+2
  83         IF(MSTP(61).GE.1.AND.MINT(43).NE.1.AND.ISUB.NE.95)
  84      &  CALL PYSSPA_HIJING(IPU1,IPU2)
  85         NSAV1=N
  86
  87 C...Multiple interactions.
  88         IF(MSTP(81).GE.1.AND.MINT(43).EQ.4.AND.ISUB.NE.95)
  89      &  CALL PYMULT_HIJING(6)
  90         MINT(1)=ISUB
  91         NSAV2=N
  92
  93 C...Hadron remnants and primordial kT.
  94         CALL PYREMN_HIJING(IPU1,IPU2)
  95         IF(MINT(51).EQ.1) GOTO 100
  96         NSAV3=N
  97
  98 C...Showering of final state partons (optional).
  99         IPU3=MINT(84)+3
 100         IPU4=MINT(84)+4
 101         IF(MSTP(71).GE.1.AND.ISUB.NE.95.AND.K(IPU3,1).GT.0.AND.
 102      &  K(IPU3,1).LE.10.AND.K(IPU4,1).GT.0.AND.K(IPU4,1).LE.10) THEN
 103           QMAX=SQRT(PARP(71)*VINT(52))
 104           IF(ISUB.EQ.5) QMAX=SQRT(PMAS(23,1)**2)
 105           IF(ISUB.EQ.8) QMAX=SQRT(PMAS(24,1)**2)
 106           CALL LUSHOW_HIJING(IPU3,IPU4,QMAX)
 107         ENDIF
 108
 109 C...Sum up transverse and longitudinal momenta.
 110         IF(IOVL.EQ.1) THEN
 111           PARI(65)=2.*PARI(17)
 112           DO 130 I=MSTP(126)+1,N
 113           IF(K(I,1).LE.0.OR.K(I,1).GT.10) GOTO 130
 114           PT=SQRT(P(I,1)**2+P(I,2)**2)
 115           PARI(69)=PARI(69)+PT
 116           IF(I.LE.NSAV1.OR.I.GT.NSAV3) PARI(66)=PARI(66)+PT
 117           IF(I.GT.NSAV1.AND.I.LE.NSAV2) PARI(68)=PARI(68)+PT
 118   130     CONTINUE
 119           PARI(67)=PARI(68)
 120           PARI(71)=VINT(151)
 121           PARI(72)=VINT(152)
 122           PARI(73)=VINT(151)
 123           PARI(74)=VINT(152)
 124         ENDIF
 125
 126 C...Decay of final state resonances.
 127         IF(MSTP(41).GE.1.AND.ISUB.NE.95) CALL PYRESD_HIJING
 128
 129       ELSE
 130 C...Diffractive and elastic scattering.
 131         CALL PYDIFF_HIJING
 132         IF(IOVL.EQ.1) THEN
 133           PARI(65)=2.*PARI(17)
 134           PARI(66)=PARI(65)
 135           PARI(69)=PARI(65)
 136         ENDIF
 137       ENDIF
 138
 139 C...Recalculate energies from momenta and masses (if desired).
 140       IF(MSTP(113).GE.1) THEN
 141         DO 140 I=MINT(83)+1,N
 142   140   IF(K(I,1).GT.0.AND.K(I,1).LE.10) P(I,4)=SQRT(P(I,1)**2+
 143      &  P(I,2)**2+P(I,3)**2+P(I,5)**2)
 144       ENDIF
 145
 146 C...Rearrange partons along strings, check invariant mass cuts.
 147       MSTU(28)=0
 148       CALL LUPREP_HIJING(MINT(84)+1)
 149       IF(MSTP(112).EQ.1.AND.MSTU(28).EQ.3) GOTO 100
 150       IF(MSTP(125).EQ.0.OR.MSTP(125).EQ.1) THEN
 151         DO 150 I=MINT(84)+1,N
 152         IF(K(I,2).NE.94) GOTO 150
 153         K(I+1,3)=MOD(K(I+1,4)/MSTU(5),MSTU(5))
 154         K(I+2,3)=MOD(K(I+2,4)/MSTU(5),MSTU(5))
 155   150   CONTINUE
 156         CALL LUEDIT_HIJING(12)
 157         CALL LUEDIT_HIJING(14)
 158         IF(MSTP(125).EQ.0) CALL LUEDIT_HIJING(15)
 159         IF(MSTP(125).EQ.0) MINT(4)=0
 160       ENDIF
 161
 162 C...Introduce separators between sections in LULIST_HIJING event listing.
 163       IF(IOVL.EQ.1.AND.MSTP(125).LE.0) THEN
 164         MSTU(70)=1
 165         MSTU(71)=N
 166       ELSEIF(IOVL.EQ.1) THEN
 167         MSTU(70)=3
 168         MSTU(71)=2
 169         MSTU(72)=MINT(4)
 170         MSTU(73)=N
 171       ENDIF
 172
 173 C...Perform hadronization (if desired).
 174       IF(MSTP(111).GE.1) CALL LUEXEC_HIJING
 175       IF(MSTP(125).EQ.0.OR.MSTP(125).EQ.1) CALL LUEDIT_HIJING(14)
 176
 177 C...Calculate Monte Carlo estimates of cross-sections.
 178   160 IF(IOVL.EQ.1) THEN
 179         IF(MSTP(111).NE.-1) NGEN(ISUB,3)=NGEN(ISUB,3)+1
 180         NGEN(0,3)=NGEN(0,3)+1
 181         XSEC(0,3)=0.
 182         DO 170 I=1,200
 183         IF(I.EQ.96) THEN
 184           XSEC(I,3)=0.
 185         ELSEIF(MSUB(95).EQ.1.AND.(I.EQ.11.OR.I.EQ.12.OR.I.EQ.13.OR.
 186      &  I.EQ.28.OR.I.EQ.53.OR.I.EQ.68)) THEN
 187           XSEC(I,3)=XSEC(96,2)*NGEN(I,3)/MAX(1.,FLOAT(NGEN(96,1))*
 188      &    FLOAT(NGEN(96,2)))
 189         ELSEIF(NGEN(I,1).EQ.0) THEN
 190           XSEC(I,3)=0.
 191         ELSEIF(NGEN(I,2).EQ.0) THEN
 192           XSEC(I,3)=XSEC(I,2)*NGEN(0,3)/(FLOAT(NGEN(I,1))*
 193      &    FLOAT(NGEN(0,2)))
 194         ELSE
 195           XSEC(I,3)=XSEC(I,2)*NGEN(I,3)/(FLOAT(NGEN(I,1))*
 196      &    FLOAT(NGEN(I,2)))
 197         ENDIF
 198   170   XSEC(0,3)=XSEC(0,3)+XSEC(I,3)
 199         IF(MSUB(95).EQ.1) THEN
 200           NGENS=NGEN(91,3)+NGEN(92,3)+NGEN(93,3)+NGEN(94,3)+NGEN(95,3)
 201           XSECS=XSEC(91,3)+XSEC(92,3)+XSEC(93,3)+XSEC(94,3)+XSEC(95,3)
 202           XMAXS=XSEC(95,1)
 203           IF(MSUB(91).EQ.1) XMAXS=XMAXS+XSEC(91,1)
 204           IF(MSUB(92).EQ.1) XMAXS=XMAXS+XSEC(92,1)
 205           IF(MSUB(93).EQ.1) XMAXS=XMAXS+XSEC(93,1)
 206           IF(MSUB(94).EQ.1) XMAXS=XMAXS+XSEC(94,1)
 207           FAC=1.
 208           IF(NGENS.LT.NGEN(0,3)) FAC=(XMAXS-XSECS)/(XSEC(0,3)-XSECS)
 209           XSEC(11,3)=FAC*XSEC(11,3)
 210           XSEC(12,3)=FAC*XSEC(12,3)
 211           XSEC(13,3)=FAC*XSEC(13,3)
 212           XSEC(28,3)=FAC*XSEC(28,3)
 213           XSEC(53,3)=FAC*XSEC(53,3)
 214           XSEC(68,3)=FAC*XSEC(68,3)
 215           XSEC(0,3)=XSEC(91,3)+XSEC(92,3)+XSEC(93,3)+XSEC(94,3)+
 216      &    XSEC(95,1)
 217         ENDIF
 218
 219 C...Store final information.
 220         MINT(5)=MINT(5)+1
 221         MSTI(3)=MINT(3)
 222         MSTI(4)=MINT(4)
 223         MSTI(5)=MINT(5)
 224         MSTI(6)=MINT(6)
 225         MSTI(7)=MINT(7)
 226         MSTI(8)=MINT(8)
 227         MSTI(13)=MINT(13)
 228         MSTI(14)=MINT(14)
 229         MSTI(21)=MINT(21)
 230         MSTI(22)=MINT(22)
 231         MSTI(23)=MINT(23)
 232         MSTI(24)=MINT(24)
 233         MSTI(25)=MINT(25)
 234         MSTI(26)=MINT(26)
 235         MSTI(31)=MINT(31)
 236         PARI(1)=XSEC(0,3)
 237         PARI(2)=XSEC(0,3)/MINT(5)
 238         PARI(31)=VINT(141)
 239         PARI(32)=VINT(142)
 240         IF(ISUB.NE.95.AND.MINT(7)*MINT(8).NE.0) THEN
 241           PARI(42)=2.*VINT(47)/VINT(1)
 242           DO 180 IS=7,8
 243           PARI(36+IS)=P(MINT(IS),3)/VINT(1)
 244           PARI(38+IS)=P(MINT(IS),4)/VINT(1)
 245           I=MINT(IS)
 246           PR=MAX(1E-20,P(I,5)**2+P(I,1)**2+P(I,2)**2)
 247           PARI(40+IS)=SIGN(LOG(MIN((SQRT(PR+P(I,3)**2)+ABS(P(I,3)))/
 248      &    SQRT(PR),1E20)),P(I,3))
 249           PR=MAX(1E-20,P(I,1)**2+P(I,2)**2)
 250           PARI(42+IS)=SIGN(LOG(MIN((SQRT(PR+P(I,3)**2)+ABS(P(I,3)))/
 251      &    SQRT(PR),1E20)),P(I,3))
 252           PARI(44+IS)=P(I,3)/SQRT(P(I,1)**2+P(I,2)**2+P(I,3)**2)
 253           PARI(46+IS)=ULANGL_HIJING(P(I,3),SQRT(P(I,1)**2+P(I,2)**2))
 254           PARI(48+IS)=ULANGL_HIJING(P(I,1),P(I,2))
 255   180     CONTINUE
 256         ENDIF
 257         PARI(61)=VINT(148)
 258         IF(ISET(ISUB).EQ.1.OR.ISET(ISUB).EQ.3) THEN
 259           MSTU(161)=MINT(21)
 260           MSTU(162)=0
 261         ELSE
 262           MSTU(161)=MINT(21)
 263           MSTU(162)=MINT(22)
 264         ENDIF
 265       ENDIF
 266
 267 C...Prepare to go to next overlayed event.
 268       MSTI(41)=IOVL
 269       IF(IOVL.GE.2.AND.IOVL.LE.10) MSTI(40+IOVL)=ISUB
 270       IF(MSTU(70).LT.10) THEN
 271         MSTU(70)=MSTU(70)+1
 272         MSTU(70+MSTU(70))=N
 273       ENDIF
 274       MINT(83)=N
 275       MINT(84)=N+MSTP(126)
 276   190 CONTINUE
 277
 278 C...Information on overlayed events.
 279       IF(MSTP(131).EQ.1.AND.MSTP(133).GE.1) THEN
 280         PARI(91)=VINT(132)
 281         PARI(92)=VINT(133)
 282         PARI(93)=VINT(134)
 283         IF(MSTP(133).EQ.2) PARI(93)=PARI(93)*XSEC(0,3)/VINT(131)
 284       ENDIF
 285
 286 C...Transform to the desired coordinate frame.
 287   200 CALL PYFRAM_HIJING(MSTP(124))
 288
 289       RETURN
 290       END