pythiaeRHIC/pythia/pyevol.f

0001
0002 C***********************************************************************
0003
0004 C...PYEVOL
0005 C...Handles intertwined pT-ordered spacelike initial-state parton
0006 C...and multiple interactions.
0007
0008       SUBROUTINE PYEVOL(MODE,PT2MAX,PT2MIN)
0009 C...Mode = -1 : Initialize first time. Determine MAX and MIN scales.
0010 C...MODE =  0 : (Re-)initialize ISR/MI evolution.
0011 C...Mode =  1 : Evolve event from PT2MAX to PT2MIN.
0012
0013 C...Double precision and integer declarations.
0014       IMPLICIT DOUBLE PRECISION(A-H, O-Z)
0015       IMPLICIT INTEGER(I-N)
0016       INTEGER PYK,PYCHGE,PYCOMP
0017 C...External
0018       EXTERNAL PYALPS
0019       DOUBLE PRECISION PYALPS
0020 C...Parameter statement for maximum size of showers.
0021       PARAMETER (MAXNUR=1000)
0022 C...Commonblocks.
0023       COMMON/PYPART/NPART,NPARTD,IPART(MAXNUR),PTPART(MAXNUR)
0024       COMMON/PYJETS/N,NPAD,K(4000,5),P(4000,5),V(4000,5)
0025       COMMON/PYDAT1/MSTU(200),PARU(200),MSTJ(200),PARJ(200)
0026       COMMON/PYDAT2/KCHG(500,4),PMAS(500,4),PARF(2000),VCKM(4,4)
0027       COMMON/PYPARS/MSTP(200),PARP(200),MSTI(200),PARI(200)
0028       COMMON/PYINT1/MINT(400),VINT(400)
0029       COMMON/PYINT2/ISET(500),KFPR(500,2),COEF(500,20),ICOL(40,4,2)
0030       COMMON/PYINT3/XSFX(2,-40:40),ISIG(1000,3),SIGH(1000)
0031       COMMON/PYINTM/KFIVAL(2,3),NMI(2),IMI(2,800,2),NVC(2,-6:6),
0032      &     XASSOC(2,-6:6,240),XPSVC(-6:6,-1:240),PVCTOT(2,-1:1),
0033      &     XMI(2,240),PT2MI(240),IMISEP(0:240)
0034       COMMON/PYCTAG/NCT,MCT(4000,2)
0035       COMMON/PYISMX/MIMX,JSMX,KFLAMX,KFLCMX,KFBEAM(2),NISGEN(2,240),
0036      &     PT2MX,PT2AMX,ZMX,RM2CMX,Q2BMX,PHIMX
0037       COMMON/PYISJN/MJN1MX,MJN2MX,MJOIND(2,240)
0038 C...Local arrays and saved variables.
0039       DIMENSION VINTSV(11:80),KSAV(4,5),PSAV(4,5),VSAV(4,5),SHAT(240)
0040       SAVE NSAV,NPARTS,M15SV,M16SV,M21SV,M22SV,VINTSV,SHAT,ISUBHD,ALAM3
0041      &     ,PSAV,KSAV,VSAV
0042
0043       SAVE /PYPART/,/PYJETS/,/PYDAT1/,/PYDAT2/,/PYPARS/,/PYINT1/,
0044      &     /PYINT2/,/PYINT3/,/PYINTM/,/PYCTAG/,/PYISMX/,/PYISJN/
0045
0046 C----------------------------------------------------------------------
0047 C...MODE=-1: Pre-initialization. Store info on hard scattering etc,
0048 C...done only once per event, while MODE=0 is repeated each time the
0049 C...evolution needs to be restarted.
0050       IF (MODE.EQ.-1) THEN
0051         ISUBHD=MINT(1)
0052         NSAV=N
0053         NPARTS=NPART
0054 C...Store hard scattering variables
0055         M15SV=MINT(15)
0056         M16SV=MINT(16)
0057         M21SV=MINT(21)
0058         M22SV=MINT(22)
0059         DO 100 J=11,80
0060           VINTSV(J)=VINT(J)
0061   100   CONTINUE
0062         DO 120 J=1,5
0063           DO 110 IS=1,4
0064             I=IS+MINT(84)
0065             PSAV(IS,J)=P(I,J)
0066             KSAV(IS,J)=K(I,J)
0067             VSAV(IS,J)=V(I,J)
0068   110     CONTINUE
0069   120   CONTINUE
0070
0071 C...Set shat for hardest scattering
0072         SHAT(1)=VINT(44)
0073         IF(ISET(ISUBHD).GE.3.AND.ISET(ISUBHD).LE.5) SHAT(1)=VINT(26)
0074      &       *VINT(2)
0075
0076 C...Compute 3-Flavour Lambda_QCD (sets absolute lowest PT scale below)
0077         RMC=PMAS(4,1)
0078         RMB=PMAS(5,1)
0079         ALAM4=PARP(61)
0080         IF(MSTU(112).LT.4) ALAM4=PARP(61)*(PARP(61)/RMC)**(2D0/25D0)
0081         IF(MSTU(112).GT.4) ALAM4=PARP(61)*(RMB/PARP(61))**(2D0/25D0)
0082         ALAM3=ALAM4*(RMC/ALAM4)**(2D0/27D0)
0083
0084 C----------------------------------------------------------------------
0085 C...MODE= 0: Initialize ISR/MI evolution, i.e. begin from hardest
0086 C...interaction initiators, with no previous evolution. Check the input
0087 C...PT2MAX and PT2MIN and impose extra constraints on minimum PT2 (e.g.
0088 C...must be larger than Lambda_QCD) and maximum PT2 (e.g. must be
0089 C...smaller than the CM energy / 2.)
0090       ELSEIF (MODE.EQ.0) THEN
0091 C...Reset counters and switches
0092         N=NSAV
0093         NPART=NPARTS
0094         MINT(30)=0
0095         MINT(31)=1
0096         MINT(36)=1
0097 C...Reset hard scattering variables
0098         MINT(1)=ISUBHD
0099         DO 130 J=11,80
0100           VINT(J)=VINTSV(J)
0101   130   CONTINUE
0102         DO 150 J=1,5
0103           DO 140 IS=1,4
0104             I=IS+MINT(84)
0105             P(I,J)=PSAV(IS,J)
0106             K(I,J)=KSAV(IS,J)
0107             V(I,J)=VSAV(IS,J)
0108             P(MINT(83)+4+IS,J)=PSAV(IS,J)
0109             V(MINT(83)+4+IS,J)=VSAV(IS,J)
0110   140     CONTINUE
0111   150   CONTINUE
0112 C...Reset statistics on activity in event.
0113         DO 160 J=351,359
0114           MINT(J)=0
0115           VINT(J)=0D0
0116   160   CONTINUE
0117 C...Reset extra companion reweighting factor
0118         VINT(140)=1D0
0119
0120 C...We do not generate MI for soft process (ISUB=95), but the
0121 C...initialization must be done regardless, for later purposes.
0122         MINT(36)=1
0123
0124 C...Initialize multiple interactions.
0125         CALL PYPTMI(-1,PTDUM1,PTDUM2,PTDUM3,IDUM)
0126         IF(MINT(51).NE.0) RETURN
0127
0128 C...Decide whether quarks in hard scattering were valence or sea
0129         PT2HD=VINT(54)
0130         DO 170 JS=1,2
0131           MINT(30)=JS
0132           CALL PYPTMI(2,PT2HD,PTDUM2,PTDUM3,IDUM)
0133           IF(MINT(51).NE.0) RETURN
0134   170   CONTINUE
0135
0136 C...Set lower cutoff for PT2 iteration and colour interference PT2 scale
0137         VINT(18)=0D0
0138         IF(MSTP(70).EQ.0) THEN
0139           PT20=PARP(62)**2
0140           PT2MIN=MAX(PT2MIN,PT20,(1.1D0*ALAM3)**2)
0141         ELSEIF(MSTP(70).EQ.1) THEN
0142           PT20=(PARP(81)*(VINT(1)/PARP(89))**PARP(90))**2
0143           PT2MIN=MAX(PT2MIN,PT20,(1.1D0*ALAM3)**2)
0144         ELSE
0145           VINT(18)=(PARP(82)*(VINT(1)/PARP(89))**PARP(90))**2
0146           PT2MIN=MAX(PT2MIN,(1.1D0*ALAM3)**2)
0147         ENDIF
0148 C...Also store PT2MIN in VINT(17).
0149   180   VINT(17)=PT2MIN
0150
0151 C...Set FS masses zero now.
0152         VINT(63)=0D0
0153         VINT(64)=0D0
0154
0155 C...Initialize IS showers with VINT(56) as max scale.
0156         PT2ISR=VINT(56)
0157         CALL PYPTIS(-1,PT2ISR,PT2MIN,PT2DUM,IFAIL)
0158         IF(MINT(51).NE.0) RETURN
0159
0160         RETURN
0161
0162 C----------------------------------------------------------------------
0163 C...MODE= 1: Evolve event from PTMAX to PTMIN.
0164       ELSEIF (MODE.EQ.1) THEN
0165
0166 C...Skip if no phase space.
0167   190   IF (PT2MAX.LE.PT2MIN) GOTO 330
0168
0169 C...Starting pT2 max scale (to be udpated successively).
0170         PT2CMX=PT2MAX
0171
0172 C...Evolve two sides of the event to find which branches at highest pT.
0173   200   JSMX=-1
0174         MIMX=0
0175         PT2MX=0D0
0176
0177 C...Loop over current shower initiators.
0178         IF (MSTP(61).GE.1) THEN
0179           DO 230 MI=1,MINT(31)
0180             IF (MI.GE.2.AND.MSTP(84).LE.0) GOTO 230
0181             ISUB=96
0182             IF (MI.EQ.1) ISUB=ISUBHD
0183             MINT(1)=ISUB
0184             MINT(36)=MI
0185 C...Set up shat, initiator x values, and x remaining in BR.
0186             VINT(44)=SHAT(MI)
0187             VINT(141)=XMI(1,MI)
0188             VINT(142)=XMI(2,MI)
0189             VINT(143)=1D0
0190             VINT(144)=1D0
0191             DO 210 JI=1,MINT(31)
0192               IF (JI.EQ.MINT(36)) GOTO 210
0193               VINT(143)=VINT(143)-XMI(1,JI)
0194               VINT(144)=VINT(144)-XMI(2,JI)
0195   210       CONTINUE
0196 C...Loop over sides.
0197 C...Generate trial branchings for this interaction. The hardest
0198 C...branching so far is automatically updated if necessary in /PYISMX/.
0199             DO 220 JS=1,2
0200               MINT(30)=JS
0201               CALL PYPTIS(0,PT2CMX,PT2MIN,PT2NEW,IFAIL)
0202               IF (MINT(51).NE.0) RETURN
0203   220       CONTINUE
0204   230     CONTINUE
0205         ENDIF
0206
0207 C...Generate trial additional interaction.
0208         MINT(36)=MINT(31)+1
0209   240   IF (MOD(MSTP(81),10).GE.1) THEN
0210           MINT(1)=96
0211 C...Set up X remaining in BR.
0212           VINT(143)=1D0
0213           VINT(144)=1D0
0214           DO 250 JI=1,MINT(31)
0215             VINT(143)=VINT(143)-XMI(1,JI)
0216             VINT(144)=VINT(144)-XMI(2,JI)
0217   250     CONTINUE
0218 C...Generate trial interaction
0219   260     CALL PYPTMI(0,PT2CMX,PT2MIN,PT2NEW,IFAIL)
0220           IF (MINT(51).EQ.1) RETURN
0221         ENDIF
0222
0223 C...And the winner is:
0224         IF (PT2MX.LT.PT2MIN) THEN
0225           GOTO 330
0226         ELSEIF (JSMX.EQ.0) THEN
0227 C...Accept additional interaction (may still fail).
0228           CALL PYPTMI(1,PT2NEW,PT2MIN,PT2DUM,IFAIL)
0229           IF(MINT(51).NE.0) RETURN
0230           IF (IFAIL.EQ.0) THEN
0231             SHAT(MINT(36))=VINT(44)
0232 C...Decide on flavours (valence/sea/companion).
0233             DO 270 JS=1,2
0234               MINT(30)=JS
0235               CALL PYPTMI(2,PT2NEW,PT2MIN,PT2DUM,IFAIL)
0236               IF(MINT(51).NE.0) RETURN
0237   270       CONTINUE
0238           ENDIF
0239         ELSEIF (JSMX.EQ.1.OR.JSMX.EQ.2) THEN
0240 C...Reconstruct kinematics of acceptable ISR branching.
0241 C...Set up shat, initiator x values, and x remaining in BR.
0242           MINT(30)=JSMX
0243           MINT(36)=MIMX
0244           VINT(44)=SHAT(MINT(36))
0245           VINT(141)=XMI(1,MINT(36))
0246           VINT(142)=XMI(2,MINT(36))
0247           VINT(143)=1D0
0248           VINT(144)=1D0
0249           DO 280 JI=1,MINT(31)
0250             IF (JI.EQ.MINT(36)) GOTO 280
0251             VINT(143)=VINT(143)-XMI(1,JI)
0252             VINT(144)=VINT(144)-XMI(2,JI)
0253   280     CONTINUE
0254           PT2NEW=PT2MX
0255           CALL PYPTIS(1,PT2NEW,PT2DM1,PT2DM2,IFAIL)
0256           IF (MINT(51).EQ.1) RETURN
0257         ELSEIF (JSMX.EQ.3.OR.JSMX.EQ.4) THEN
0258 C...Bookeep joining. Cannot (yet) be constructed kinematically.
0259           MINT(354)=MINT(354)+1
0260           VINT(354)=VINT(354)+SQRT(PT2MX)
0261           IF (MINT(354).EQ.1) VINT(359)=SQRT(PT2MX)
0262           MJOIND(JSMX-2,MJN1MX)=MJN2MX
0263           MJOIND(JSMX-2,MJN2MX)=MJN1MX
0264         ENDIF
0265
0266 C...Update PT2 iteration scale.
0267         PT2CMX=PT2MX
0268
0269 C...Loop back to continue evolution.
0270         IF(N.GT.MSTU(4)-MSTU(32)-10) THEN
0271           CALL PYERRM(11,'(PYEVOL:) no more memory left in PYJETS')
0272         ELSE
0273           IF (JSMX.GE.0.AND.PT2CMX.GE.PT2MIN) GOTO 200
0274         ENDIF
0275
0276 C----------------------------------------------------------------------
0277 C...MODE= 2: (Re-)store user information on hardest interaction etc.
0278       ELSEIF (MODE.EQ.2) THEN
0279
0280 C...Revert to "ordinary" meanings of some parameters.
0281   290   DO 310 JS=1,2
0282           MINT(12+JS)=K(IMI(JS,1,1),2)
0283           VINT(140+JS)=XMI(JS,1)
0284           IF(MINT(18+JS).EQ.1) VINT(140+JS)=VINT(154+JS)*XMI(JS,1)
0285           VINT(142+JS)=1D0
0286           DO 300 MI=1,MINT(31)
0287             VINT(142+JS)=VINT(142+JS)-XMI(JS,MI)
0288   300     CONTINUE
0289   310   CONTINUE
0290
0291 C...Restore saved quantities for hardest interaction.
0292         MINT(1)=ISUBHD
0293         MINT(15)=M15SV
0294         MINT(16)=M16SV
0295         MINT(21)=M21SV
0296         MINT(22)=M22SV
0297         DO 320 J=11,80
0298           VINT(J)=VINTSV(J)
0299   320   CONTINUE
0300
0301       ENDIF
0302
0303   330 RETURN
0304       END