pythiaeRHIC/pythia/pyevnw.f

0001
0002 C*********************************************************************
0003
0004 C...PYEVNW
0005 C...Administers the generation of a high-pT event via calls to
0006 C...a number of subroutines for the new multiple interactions and
0007 C...showering framework.
0008
0009       SUBROUTINE PYEVNW
0010
0011 C...Double precision and integer declarations.
0012       IMPLICIT DOUBLE PRECISION(A-H, O-Z)
0013       IMPLICIT INTEGER(I-N)
0014       INTEGER PYK,PYCHGE,PYCOMP
0015 C...Commonblocks.
0016       COMMON/PYJETS/N,NPAD,K(4000,5),P(4000,5),V(4000,5)
0017       COMMON/PYCTAG/NCT,MCT(4000,2)
0018       COMMON/PYDAT1/MSTU(200),PARU(200),MSTJ(200),PARJ(200)
0019       COMMON/PYDAT2/KCHG(500,4),PMAS(500,4),PARF(2000),VCKM(4,4)
0020       COMMON/PYDAT3/MDCY(500,3),MDME(8000,2),BRAT(8000),KFDP(8000,5)
0021       COMMON/PYPARS/MSTP(200),PARP(200),MSTI(200),PARI(200)
0022       COMMON/PYINT1/MINT(400),VINT(400)
0023       COMMON/PYINT2/ISET(500),KFPR(500,2),COEF(500,20),ICOL(40,4,2)
0024       COMMON/PYINT4/MWID(500),WIDS(500,5)
0025       COMMON/PYINT5/NGENPD,NGEN(0:500,3),XSEC(0:500,3)
0026       COMMON/PYINTM/KFIVAL(2,3),NMI(2),IMI(2,800,2),NVC(2,-6:6),
0027      &     XASSOC(2,-6:6,240),XPSVC(-6:6,-1:240),PVCTOT(2,-1:1),
0028      &     XMI(2,240),PT2MI(240),IMISEP(0:240)
0029       SAVE /PYJETS/,/PYCTAG/,/PYDAT1/,/PYDAT2/,/PYDAT3/,
0030      &     /PYPARS/,/PYINT1/,/PYINT2/,/PYINT4/,/PYINT5/,/PYINTM/
0031 C...Local arrays.
0032       DIMENSION VTX(4)
0033
0034 C...Stop if no subprocesses on.
0035       IF(MINT(121).EQ.1.AND.MSTI(53).EQ.1) THEN
0036         WRITE(MSTU(11),5100)
0037         CALL PYSTOP(1)
0038       ENDIF
0039
0040 C...Initial values for some counters.
0041       MSTU(1)=0
0042       MSTU(2)=0
0043       N=0
0044       MINT(5)=MINT(5)+1
0045       MINT(7)=0
0046       MINT(8)=0
0047       MINT(30)=0
0048       MINT(83)=0
0049       MINT(84)=MSTP(126)
0050       MSTU(24)=0
0051       MSTU70=0
0052       MSTJ14=MSTJ(14)
0053 C...Normally, use K(I,4:5) colour info rather than /PYCT/.
0054       NCT=0
0055       MINT(33)=0
0056
0057 C...Let called routines know call is from PYEVNW (not PYEVNT).
0058       MINT(35)=3
0059
0060 C...If variable energies: redo incoming kinematics and cross-section.
0061       MSTI(61)=0
0062       IF(MSTP(171).EQ.1) THEN
0063         CALL PYINKI(1)
0064         IF(MSTI(61).EQ.1) THEN
0065           MINT(5)=MINT(5)-1
0066           RETURN
0067         ENDIF
0068         IF(MINT(121).GT.1) CALL PYSAVE(3,1)
0069         CALL PYXTOT
0070       ENDIF
0071
0072 C...Loop over number of pileup events; check space left.
0073       IF(MSTP(131).LE.0) THEN
0074         NPILE=1
0075       ELSE
0076         CALL PYPILE(2)
0077         NPILE=MINT(81)
0078       ENDIF
0079       DO 300 IPILE=1,NPILE
0080         IF(MINT(84)+100.GE.MSTU(4)) THEN
0081           CALL PYERRM(11,
0082      &    '(PYEVNW:) no more space in PYJETS for pileup events')
0083           IF(MSTU(21).GE.1) GOTO 310
0084         ENDIF
0085         MINT(82)=IPILE
0086
0087 C...Generate variables of hard scattering.
0088         MINT(51)=0
0089         MSTI(52)=0
0090   100   CONTINUE
0091         IF(MINT(51).NE.0.OR.MSTU(24).NE.0) MSTI(52)=MSTI(52)+1
0092         MINT(31)=0
0093         MINT(39)=0
0094         MINT(36)=0
0095         MINT(51)=0
0096         MINT(57)=0
0097         CALL PYRAND
0098         IF(MSTI(61).EQ.1) THEN
0099           MINT(5)=MINT(5)-1
0100           RETURN
0101         ENDIF
0102         IF(MINT(51).EQ.2) RETURN
0103         ISUB=MINT(1)
0104         IF(MSTP(111).EQ.-1) GOTO 290
0105
0106 C...Loopback point if PYPREP fails, especially for junction topologies.
0107         NPREP=0
0108         MNT31S=MINT(31)
0109   110   NPREP=NPREP+1
0110         MINT(31)=MNT31S
0111
0112         IF((ISUB.LE.90.OR.ISUB.GE.95).AND.ISUB.NE.99) THEN
0113 C...Hard scattering (including low-pT):
0114 C...reconstruct kinematics and colour flow of hard scattering.
0115           MINT31=MINT(31)
0116   120     MINT(31)=MINT31
0117           MINT(51)=0
0118           CALL PYSCAT
0119           IF(MINT(51).EQ.1) GOTO 100
0120           NPARTD=N
0121           NFIN=N
0122
0123 C...Intertwined initial state showers and multiple interactions.
0124 C...Force no IS showers if no pdfs defined: MSTP(61) -> 0 for PYEVOL.
0125 C...Force no MI if cross section not known: MSTP(81) -> 0 for PYEVOL.
0126           MSTP61=MSTP(61)
0127           IF (MINT(47).LT.2) MSTP(61)=0
0128           MSTP81=MSTP(81)
0129           IF (MINT(50).EQ.0) MSTP(81)=0
0130           IF ((MSTP(61).GE.1.OR.MOD(MSTP(81),10).GE.0).AND.
0131      &    MINT(111).NE.12) THEN
0132 C...Absolute max pT2 scale for evolution: phase space limit.
0133             PT2MXS=0.25D0*VINT(2)
0134 C...Check if more constrained by ISR and MI max scales:
0135             PT2MXS=MIN(PT2MXS,MAX(VINT(56),VINT(62)))
0136 C...Loopback point in case of failure in evolution.
0137             LOOP=0
0138   130       LOOP=LOOP+1
0139             MINT(51)=0
0140             IF(LOOP.GT.100) THEN
0141               CALL PYERRM(9,'(PYEVNW:) failed to evolve shower or '
0142      &             //'multiple interactions.')
0143               MINT(51)=1
0144               RETURN
0145             ENDIF
0146
0147 C...Pre-initialization of interleaved MI/ISR/JI evolution, only done
0148 C...once per event. (E.g. compute constants and save variables to be
0149 C...restored later in case of failure.)
0150             IF (LOOP.EQ.1) CALL PYEVOL(-1,DUMMY1,DUMMY2)
0151
0152 C...Initialize interleaved MI/ISR/JI evolution.
0153 C...PT2MAX: absolute upper limit for evolution - Initialization may
0154 C...        return a PT2MAX which is lower than this.
0155 C...PT2MIN: absolute lower limit for evolution - Initialization may
0156 C...        return a PT2MIN which is larger than this (e.g. Lambda_QCD).
0157             PT2MAX=PT2MXS
0158             PT2MIN=0D0
0159             CALL PYEVOL(0,PT2MAX,PT2MIN)
0160             IF (MINT(51).EQ.1) GOTO 130
0161
0162 C...Perform interleaved MI/ISR/JI evolution from PT2MAX to PT2MIN.
0163 C...In principle factorized, so can be stopped and restarted.
0164 C...Example: stop/start at pT=10 GeV. (Commented out for now.)
0165 C            PT2MED=MAX(10D0**2,PT2MIN)
0166 C            CALL PYEVOL(1,PT2MAX,PT2MED)
0167 C            IF (MINT(51).EQ.1) GOTO 160
0168 C            PT2MAX=PT2MED
0169             CALL PYEVOL(1,PT2MAX,PT2MIN)
0170             IF (MINT(51).EQ.1) GOTO 130
0171
0172 C...Finalize interleaved MI/ISR/JI evolution.
0173             CALL PYEVOL(2,PT2MAX,PT2MIN)
0174             IF (MINT(51).EQ.1) GOTO 130
0175
0176           ENDIF
0177           MSTP(61)=MSTP61
0178           MSTP(81)=MSTP81
0179           IF(MINT(51).EQ.1) GOTO 100
0180 C...(MINT(52) is actually obsolete in this routine. Set anyway
0181 C...to ensure PYDOCU stable.)
0182           MINT(52)=N
0183           MINT(53)=N
0184
0185 C...Beam remnants - new scheme.
0186   140     IF(MINT(50).EQ.1) THEN
0187             IF (ISUB.EQ.95) MINT(31)=1
0188
0189 C...Beam remnant flavour and colour assignments - new scheme.
0190             CALL PYMIHK
0191             IF(MINT(51).EQ.1.AND.MINT(57).GE.1.AND.MINT(57).LE.5)
0192      &           GOTO 120
0193             IF(MINT(51).EQ.1) GOTO 100
0194
0195 C...Primordial kT and beam remnant momentum sharing - new scheme.
0196             CALL PYMIRM
0197             IF(MINT(51).EQ.1.AND.MINT(57).GE.1.AND.MINT(57).LE.5)
0198      &      GOTO 120
0199             IF(MINT(51).EQ.1) GOTO 100
0200             IF (ISUB.EQ.95) MINT(31)=0
0201           ELSEIF(MINT(111).NE.12) THEN
0202 C...Hadron remnants and primordial kT - old model.
0203 C...Happens e.g. for direct photon on one side.
0204             IPU1=IMI(1,1,1)
0205             IPU2=IMI(2,1,1)
0206             CALL PYREMN(IPU1,IPU2)
0207             IF(MINT(51).EQ.1.AND.MINT(57).GE.1.AND.MINT(57).LE.5) GOTO
0208      &           110
0209             IF(MINT(51).EQ.1) GOTO 100
0210 C...PYREMN does not set colour tags for BRs, so needs to be done now.
0211             DO 160 I=MINT(53)+1,N
0212               DO 150 KCS=4,5
0213                 IDA=MOD(K(I,KCS),MSTU(5))
0214                 IF (IDA.NE.0) THEN
0215                   MCT(I,KCS-3)=MCT(IDA,6-KCS)
0216                 ELSE
0217                   MCT(I,KCS-3)=0
0218                 ENDIF
0219   150         CONTINUE
0220   160       CONTINUE
0221 C...Instruct PYPREP to use colour tags
0222             MINT(33)=1
0223
0224             DO 360 MQGST=1,2
0225               DO 350 I=MINT(84)+1,N
0226
0227 C...Look for coloured string endpoint, or (later) leftover gluon.
0228                 IF (K(I,1).NE.3) GOTO 350
0229                 KC=PYCOMP(K(I,2))
0230                 IF(KC.EQ.0) GOTO 350
0231                 KQ=KCHG(KC,2)
0232                 IF(KQ.EQ.0.OR.(MQGST.EQ.1.AND.KQ.EQ.2)) GOTO 350
0233
0234 C...  Pick up loose string end with no previous tag.
0235                 KCS=4
0236                 IF(KQ*ISIGN(1,K(I,2)).LT.0) KCS=5
0237                 IF(MCT(I,KCS-3).NE.0) GOTO 350
0238
0239                 CALL PYCTTR(I,KCS,I)
0240                 IF(MINT(51).NE.0) RETURN
0241
0242  350          CONTINUE
0243  360        CONTINUE
0244 C...Now delete any colour processing information if set (since partons
0245 C...otherwise not FS showered!)
0246             DO 170 I=MINT(84)+1,N
0247               IF (I.LE.N) THEN
0248                 K(I,4)=MOD(K(I,4),MSTU(5)**2)
0249                 K(I,5)=MOD(K(I,5),MSTU(5)**2)
0250               ENDIF
0251   170       CONTINUE
0252           ENDIF
0253
0254 C...Showering of final state partons (optional).
0255           ALAMSV=PARJ(81)
0256           PARJ(81)=PARP(72)
0257           IF(MSTP(71).GE.1.AND.ISET(ISUB).GE.1.AND.ISET(ISUB).LE.10)
0258      &    THEN
0259             QMAX=VINT(55)
0260             IF(ISET(ISUB).EQ.2) QMAX=SQRT(PARP(71))*VINT(55)
0261             CALL PYPTFS(1,QMAX,0D0,PTGEN)
0262 C...External processes: handle successive showers.
0263           ELSEIF(ISET(ISUB).EQ.11) THEN
0264             CALL PYADSH(NFIN)
0265           ENDIF
0266           PARJ(81)=ALAMSV
0267
0268 C...Allow possibility for user to abort event generation.
0269           IVETO=0
0270           IF(IPILE.EQ.1.AND.MSTP(143).EQ.1) CALL PYVETO(IVETO) ! sm
0271           IF(IVETO.EQ.1) GOTO 100
0272
0273
0274 C...Decay of final state resonances.
0275           MINT(32)=0
0276           IF(MSTP(41).GE.1.AND.ISET(ISUB).LE.10) THEN
0277             CALL PYRESD(0)
0278             IF(MINT(51).NE.0) GOTO 100
0279           ENDIF
0280
0281           IF(MINT(51).EQ.1) GOTO 100
0282
0283         ELSEIF(ISUB.NE.99) THEN
0284 C...Diffractive and elastic scattering.
0285           CALL PYDIFF
0286
0287         ELSE
0288 C...DIS scattering (photon flux external).
0289           CALL PYDISG
0290           IF(MINT(51).EQ.1) GOTO 100
0291         ENDIF
0292
0293 C...Check that no odd resonance left undecayed.
0294         MINT(54)=N
0295         IF(MSTP(111).GE.1) THEN
0296           NFIX=N
0297           DO 180 I=MINT(84)+1,NFIX
0298             IF(K(I,1).GE.1.AND.K(I,1).LE.10.AND.K(I,2).NE.21.AND.
0299      &      K(I,2).NE.22) THEN
0300               KCA=PYCOMP(K(I,2))
0301               IF(MWID(KCA).NE.0.AND.MDCY(KCA,1).GE.1) THEN
0302                 CALL PYRESD(I)
0303                 IF(MINT(51).EQ.1) GOTO 100
0304               ENDIF
0305             ENDIF
0306   180     CONTINUE
0307         ENDIF
0308
0309 C...Boost hadronic subsystem to overall rest frame.
0310 C..(Only relevant when photon inside lepton beam.)
0311         IF(MINT(141).NE.0.OR.MINT(142).NE.0) CALL PYGAGA(4,WTGAGA)
0312
0313 C...Recalculate energies from momenta and masses (if desired).
0314         IF(MSTP(113).GE.1) THEN
0315           DO 190 I=MINT(83)+1,N
0316             IF(K(I,1).GT.0.AND.K(I,1).LE.10) P(I,4)=SQRT(P(I,1)**2+
0317      &      P(I,2)**2+P(I,3)**2+P(I,5)**2)
0318   190     CONTINUE
0319           NRECAL=N
0320         ENDIF
0321
0322 C...Colour reconnection before string formation
0323         CALL PYFSCR(MINT(84)+1)
0324
0325 C...Rearrange partons along strings, check invariant mass cuts.
0326         MSTU(28)=0
0327         IF(MSTP(111).LE.0) MSTJ(14)=-1
0328         CALL PYPREP(MINT(84)+1)
0329         MSTJ(14)=MSTJ14
0330         IF(MINT(51).EQ.1.AND.MSTU(24).EQ.1) THEN
0331           MSTU(24)=0
0332           GOTO 100
0333         ENDIF
0334         IF(MINT(51).EQ.1) GOTO 110
0335         IF(MSTP(112).EQ.1.AND.MSTU(28).EQ.3) GOTO 100
0336         IF(MSTP(125).EQ.0.OR.MSTP(125).EQ.1) THEN
0337           DO 220 I=MINT(84)+1,N
0338             IF(K(I,2).EQ.94) THEN
0339               DO 210 I1=I+1,MIN(N,I+10)
0340                 IF(K(I1,3).EQ.I) THEN
0341                   K(I1,3)=MOD(K(I1,4)/MSTU(5),MSTU(5))
0342                   IF(K(I1,3).EQ.0) THEN
0343                     DO 200 II=MINT(84)+1,I-1
0344                         IF(K(II,2).EQ.K(I1,2)) THEN
0345                           IF(MOD(K(II,4),MSTU(5)).EQ.I1.OR.
0346      &                    MOD(K(II,5),MSTU(5)).EQ.I1) K(I1,3)=II
0347                         ENDIF
0348   200               CONTINUE
0349                     IF(K(I+1,3).EQ.0) K(I+1,3)=K(I,3)
0350                   ENDIF
0351                 ENDIF
0352   210         CONTINUE
0353             ENDIF
0354   220     CONTINUE
0355           CALL PYEDIT(12)
0356           CALL PYEDIT(14)
0357           IF(MSTP(125).EQ.0) CALL PYEDIT(15)
0358           IF(MSTP(125).EQ.0) MINT(4)=0
0359           DO 240 I=MINT(83)+1,N
0360             IF(K(I,1).EQ.11.AND.K(I,4).EQ.0.AND.K(I,5).EQ.0) THEN
0361               DO 230 I1=I+1,N
0362                 IF(K(I1,3).EQ.I.AND.K(I,4).EQ.0) K(I,4)=I1
0363                 IF(K(I1,3).EQ.I) K(I,5)=I1
0364   230         CONTINUE
0365             ENDIF
0366   240     CONTINUE
0367         ENDIF
0368
0369 C...Introduce separators between sections in PYLIST event listing.
0370         IF(IPILE.EQ.1.AND.MSTP(125).LE.0) THEN
0371           MSTU70=1
0372           MSTU(71)=N
0373         ELSEIF(IPILE.EQ.1) THEN
0374           MSTU70=3
0375           MSTU(71)=2
0376           MSTU(72)=MINT(4)
0377           MSTU(73)=N
0378         ENDIF
0379
0380 C...Go back to lab frame (needed for vertices, also in fragmentation).
0381         CALL PYFRAM(1)
0382
0383 C...Set nonvanishing production vertex (optional).
0384         IF(MSTP(151).EQ.1) THEN
0385           DO 250 J=1,4
0386             VTX(J)=PARP(150+J)*SQRT(-2D0*LOG(MAX(1D-10,PYR(0))))*
0387      &      SIN(PARU(2)*PYR(0))
0388   250     CONTINUE
0389           DO 270 I=MINT(83)+1,N
0390             DO 260 J=1,4
0391               V(I,J)=V(I,J)+VTX(J)
0392   260       CONTINUE
0393   270     CONTINUE
0394         ENDIF
0395
0396 C...Perform hadronization (if desired).
0397         IF(MSTP(111).GE.1) THEN
0398           CALL PYEXEC
0399           IF(MSTU(24).NE.0) GOTO 100
0400         ENDIF
0401         IF(MSTP(113).GE.1) THEN
0402           DO 280 I=NRECAL,N
0403             IF(P(I,5).GT.0D0) P(I,4)=SQRT(P(I,1)**2+
0404      &      P(I,2)**2+P(I,3)**2+P(I,5)**2)
0405   280     CONTINUE
0406         ENDIF
0407         IF(MSTP(125).EQ.0.OR.MSTP(125).EQ.1) CALL PYEDIT(14)
0408
0409 C...Store event information and calculate Monte Carlo estimates of
0410 C...subprocess cross-sections.
0411   290   IF(IPILE.EQ.1) CALL PYDOCU
0412
0413 C...Set counters for current pileup event and loop to next one.
0414         MSTI(41)=IPILE
0415         IF(IPILE.GE.2.AND.IPILE.LE.10) MSTI(40+IPILE)=ISUB
0416         IF(MSTU70.LT.10) THEN
0417           MSTU70=MSTU70+1
0418           MSTU(70+MSTU70)=N
0419         ENDIF
0420         MINT(83)=N
0421         MINT(84)=N+MSTP(126)
0422         IF(IPILE.LT.NPILE) CALL PYFRAM(2)
0423   300 CONTINUE
0424
0425 C...Generic information on pileup events. Reconstruct missing history.
0426       IF(MSTP(131).EQ.1.AND.MSTP(133).GE.1) THEN
0427         PARI(91)=VINT(132)
0428         PARI(92)=VINT(133)
0429         PARI(93)=VINT(134)
0430         IF(MSTP(133).GE.2) PARI(93)=PARI(93)*XSEC(0,3)/VINT(131)
0431       ENDIF
0432       CALL PYEDIT(16)
0433
0434 C...Transform to the desired coordinate frame.
0435   310 CALL PYFRAM(MSTP(124))
0436       MSTU(70)=MSTU70
0437       PARU(21)=VINT(1)
0438
0439 C...Error messages
0440  5100 FORMAT(1X,'Error: no subprocess switched on.'/
0441      &1X,'Execution stopped.')
0442
0443       RETURN
0444       END