pythiaeRHIC/pythia/pyfscr.f

0001
0002 C*********************************************************************
0003
0004 C...PYFSCR
0005 C...Performs colour annealing.
0006 C...MSTP(95) : CR Type
0007 C...         = 1  : old cut-and-paste reconnections, handled in PYMIHK
0008 C...         = 2  : Type I(no gg loops); hadron-hadron only
0009 C...         = 3  : Type I(no gg loops); all beams
0010 C...         = 4  : Type II(gg loops)  ; hadron-hadron only
0011 C...         = 5  : Type II(gg loops)  ; all beams
0012 C...         = 6  : Type S             ; hadron-hadron only
0013 C...         = 7  : Type S             ; all beams
0014 C...Types I and II are described in Sandhoff+Skands, in hep-ph/0604120.
0015 C...Type S is driven by starting only from free triplets, not octets.
0016 C...A string piece remains unchanged with probability
0017 C...    PKEEP = (1-PARP(78))**N
0018 C...This scaling corresponds to each string piece having to go through
0019 C...N other ones, each with probability PARP(78) for reconnection, where
0020 C...N is here chosen simply as the number of multiple interactions,
0021 C...for a rough scaling with the general level of activity.
0022
0023       SUBROUTINE PYFSCR(IP)
0024 C...Double precision and integer declarations.
0025       IMPLICIT DOUBLE PRECISION(A-H, O-Z)
0026       INTEGER PYK,PYCHGE,PYCOMP
0027 C...Commonblocks.
0028       COMMON/PYJETS/N,NPAD,K(4000,5),P(4000,5),V(4000,5)
0029       COMMON/PYDAT1/MSTU(200),PARU(200),MSTJ(200),PARJ(200)
0030       COMMON/PYPARS/MSTP(200),PARP(200),MSTI(200),PARI(200)
0031       COMMON/PYDAT2/KCHG(500,4),PMAS(500,4),PARF(2000),VCKM(4,4)
0032       COMMON/PYINT1/MINT(400),VINT(400)
0033 C...The common block of colour tags.
0034       COMMON/PYCTAG/NCT,MCT(4000,2)
0035       SAVE /PYJETS/,/PYDAT1/,/PYDAT2/,/PYINT1/,/PYCTAG/,
0036      &/PYPARS/
0037 C...MCN: Temporary storage of new colour tags
0038       DOUBLE PRECISION MCN(4000,2)
0039
0040 C...Function to give four-product.
0041       FOUR(I,J)=P(I,4)*P(J,4)
0042      &          -P(I,1)*P(J,1)-P(I,2)*P(J,2)-P(I,3)*P(J,3)
0043
0044 C...Check valid range of MSTP(95), local copy
0045       IF (MSTP(95).LE.1.OR.MSTP(95).GE.8) RETURN
0046       MSTP95=MOD(MSTP(95),10)
0047 C...Set whether CR allowed inside resonance systems or not
0048 C...(not implemented yet)
0049 C      MRESCR=1
0050 C      IF (MSTP(95).GE.10) MRESCR=0
0051
0052 C...Check whether colour tags already defined
0053       IF (MINT(33).EQ.0) THEN
0054 C...Erase any existing colour tags for this event
0055         DO 100 I=1,N
0056           MCT(I,1)=0
0057           MCT(I,2)=0
0058   100   CONTINUE
0059 C...Create colour tags for this event
0060         DO 120 I=1,N
0061           IF (K(I,1).EQ.3) THEN
0062             DO 110 KCS=4,5
0063               KCSIN=KCS
0064               IF (MCT(I,KCSIN-3).EQ.0) THEN
0065                 CALL PYCTTR(I,KCSIN,I)
0066               ENDIF
0067   110       CONTINUE
0068           ENDIF
0069   120 CONTINUE
0070 C...Instruct PYPREP to use colour tags
0071         MINT(33)=1
0072       ENDIF
0073
0074 C...For MSTP(95) even, only apply to hadron-hadron
0075       IF (MOD(MSTP(95),2).EQ.0) THEN
0076          KA1=IABS(MINT(11))
0077          KA2=IABS(MINT(12))
0078          IF (KA1.LT.100.OR.KA2.LT.100) GOTO 9999
0079       ENDIF
0080
0081 C...Initialize new tag array (but do not delete old yet)
0082       LCT=NCT
0083       DO 130 I=MAX(1,IP),N
0084          MCN(I,1)=0
0085          MCN(I,2)=0
0086   130 CONTINUE
0087
0088 C...For each final-state dipole, check whether string should be
0089 C...preserved.
0090       DO 150 ICT=1,NCT
0091         IC=0
0092         IA=0
0093         DO 140 I=MAX(1,IP),N
0094           IF (K(I,1).EQ.3.AND.MCT(I,1).EQ.ICT) IC=I
0095           IF (K(I,1).EQ.3.AND.MCT(I,2).EQ.ICT) IA=I
0096   140   CONTINUE
0097         IF (IC.NE.0.AND.IA.NE.0) THEN
0098 C...Chiefly consider large strings.
0099           PKEEP=(1D0-PARP(78))**MINT(31)
0100           IF (PYR(0).LE.PKEEP) THEN
0101             LCT=LCT+1
0102             MCN(IC,1)=LCT
0103             MCN(IA,2)=LCT
0104           ENDIF
0105         ENDIF
0106   150 CONTINUE
0107
0108 C...Loop over event record, starting from IP
0109 C...(Ignore junctions for now.)
0110       NLOOP=0
0111   160 NLOOP=NLOOP+1
0112       MCIMAX=0
0113       MCJMAX=0
0114       RLMAX=0D0
0115       ILMAX=0
0116       JLMAX=0
0117       DO 230 I=MAX(1,IP),N
0118          IF (K(I,1).NE.3) GOTO 230
0119 C...Check colour charge
0120          MCI=KCHG(PYCOMP(K(I,2)),2)*ISIGN(1,K(I,2))
0121          IF (MCI.EQ.0) GOTO 230
0122 C...For Seattle algorithm, only start from partons with one dangling
0123 C...colour tag
0124          IF (MSTP(95).EQ.6.OR.MSTP(95).EQ.7) THEN
0125            IF (MCI.EQ.2.AND.MCN(I,1).EQ.0.AND.MCN(I,2).EQ.0) GOTO 230
0126          ENDIF
0127 C...  Find optimal partner
0128          JLOPT=0
0129          MCJOPT=0
0130          MBROPT=0
0131          MGGOPT=0
0132          RLOPT=1D19
0133 C...Loop over I colour/anticolour, check whether already connected
0134   170    DO 220 ICL=1,2
0135             IF (MCN(I,ICL).NE.0) GOTO 220
0136             IF (ICL.EQ.1.AND.MCI.EQ.-1) GOTO 220
0137             IF (ICL.EQ.2.AND.MCI.EQ.1) GOTO 220
0138 C...Check whether this is a dangling colour tag (ie to junction!)
0139             IFOUND=0
0140             DO 180 J=MAX(1,IP),N
0141                IF (K(J,1).EQ.3.AND.MCT(J,3-ICL).EQ.MCT(I,ICL)) IFOUND=1
0142   180       CONTINUE
0143             IF (IFOUND.EQ.0) GOTO 220
0144             DO 210 J=MAX(1,IP),N
0145                IF (K(J,1).NE.3.OR.I.EQ.J) GOTO 210
0146 C...Do not make direct connections between partons in same Beam Remnant
0147                MBRSTR=0
0148                IF (K(I,3).LE.2.AND.K(J,3).LE.2.AND.K(I,3).EQ.K(J,3))
0149      &              MBRSTR=1
0150 C...Check colour charge
0151                MCJ=KCHG(PYCOMP(K(J,2)),2)*ISIGN(1,K(J,2))
0152                IF (MCJ.EQ.0.OR.(MCJ.EQ.MCI.AND.MCI.NE.2)) GOTO 210
0153 C...Check for gluon loops
0154                MGGSTR=0
0155                IF (MCJ.EQ.2.AND.MCI.EQ.2) THEN
0156                  ICLA=3-ICL
0157                  IF (MCN(I,ICLA).EQ.MCN(J,ICL).AND.MSTP(95).LE.3.AND.
0158      &                MCN(I,ICLA).NE.0) MGGSTR=1
0159                ENDIF
0160 C...Loop over J colour/anticolour, check whether already connected
0161                DO 200 JCL=1,2
0162                   IF (MCN(J,JCL).NE.0) GOTO 200
0163                   IF (JCL.EQ.ICL) GOTO 200
0164                   IF (JCL.EQ.1.AND.MCJ.EQ.-1) GOTO 200
0165                   IF (JCL.EQ.2.AND.MCJ.EQ.1) GOTO 200
0166 C...Check whether this is a dangling colour tag (ie to junction!)
0167                   IFOUND=0
0168                   DO 190 J2=MAX(1,IP),N
0169                      IF (K(J2,1).EQ.3.AND.MCT(J2,3-JCL).EQ.MCT(J,JCL))
0170      &                    IFOUND=1
0171   190             CONTINUE
0172                   IF (IFOUND.EQ.0) GOTO 200
0173 C...Save connection with smallest lambda measure
0174 C...If best so far was a BR string and this is not, also save.
0175 C...If best so far was a gg string and this is not, also save.
0176                   RL=FOUR(I,J)
0177                   IF (RL.LT.RLOPT.OR.(RL.EQ.RLOPT.AND.PYR(0).LE.0.5D0)
0178      &                 .OR.(MBROPT.EQ.1.AND.MBRSTR.EQ.0)
0179      &                 .OR.(MGGOPT.EQ.1.AND.MGGSTR.EQ.0)) THEN
0180                      RLOPT=RL
0181                      JLOPT=J
0182                      ICOPT=ICL
0183                      JCOPT=JCL
0184                      MCJOPT=MCJ
0185                      MBROPT=MBRSTR
0186                      MGGOPT=MGGSTR
0187                   ENDIF
0188   200          CONTINUE
0189   210       CONTINUE
0190   220    CONTINUE
0191          IF (JLOPT.NE.0) THEN
0192 C...Save pair with largest RLOPT so far
0193             IF (RLOPT.GE.RLMAX) THEN
0194                RLMAX=RLOPT
0195                ILMAX=I
0196                JLMAX=JLOPT
0197                ICMAX=ICOPT
0198                JCMAX=JCOPT
0199                MCJMAX=MCJOPT
0200                MCIMAX=MCI
0201             ENDIF
0202          ENDIF
0203   230 CONTINUE
0204 C...Save and iterate
0205       IF (ILMAX.GT.0) THEN
0206          LCT=LCT+1
0207          MCN(ILMAX,ICMAX)=LCT
0208          MCN(JLMAX,JCMAX)=LCT
0209          IF (NLOOP.LE.2*(N-IP)) THEN
0210             GOTO 160
0211          ELSE
0212             CALL PYERRM(31,' PYFSCR: infinite loop in color annealing')
0213             CALL PYSTOP(11)
0214          ENDIF
0215       ELSE
0216 C...Save and exit. First check for leftover gluon(s)
0217          DO 260 I=MAX(1,IP),N
0218 C...Check colour charge
0219             MCI=KCHG(PYCOMP(K(I,2)),2)*ISIGN(1,K(I,2))
0220             IF (K(I,1).NE.3.OR.MCI.NE.2) GOTO 260
0221             IF(MCN(I,1).EQ.0.AND.MCN(I,2).EQ.0) THEN
0222 C...Decide where to put left-over gluon (minimal insertion)
0223                ILMAX=0
0224                RLMAX=1D19
0225                DO 250 KCT=NCT+1,LCT
0226                   DO 240 IT=MAX(1,IP),N
0227                      IF (IT.EQ.I.OR.K(IT,1).NE.3) GOTO 240
0228                      IF (MCN(IT,1).EQ.KCT) IC=IT
0229                      IF (MCN(IT,2).EQ.KCT) IA=IT
0230   240             CONTINUE
0231                   RL=FOUR(IC,I)*FOUR(IA,I)
0232                   IF (RL.LT.RLMAX) THEN
0233                      RLMAX=RL
0234                      ICMAX=IC
0235                      IAMAX=IA
0236                   ENDIF
0237   250          CONTINUE
0238                LCT=LCT+1
0239                MCN(I,1)=MCN(ICMAX,1)
0240                MCN(I,2)=LCT
0241                MCN(ICMAX,1)=LCT
0242             ENDIF
0243   260    CONTINUE
0244          DO 270 I=MAX(1,IP),N
0245 C...Do not erase parton shower colour history
0246             IF (K(I,1).NE.3) GOTO 270
0247 C...Check colour charge
0248             MCI=KCHG(PYCOMP(K(I,2)),2)*ISIGN(1,K(I,2))
0249             IF (MCI.EQ.0) GOTO 270
0250             IF (MCN(I,1).NE.0) MCT(I,1)=MCN(I,1)
0251             IF (MCN(I,2).NE.0) MCT(I,2)=MCN(I,2)
0252   270    CONTINUE
0253       ENDIF
0254
0255  9999 RETURN
0256       END