pythiaeRHIC/pythia/pyedit.f

0001
0002 C*********************************************************************
0003
0004 C...PYEDIT
0005 C...Performs global manipulations on the event record, in particular
0006 C...to exclude unstable or undetectable partons/particles.
0007
0008       SUBROUTINE PYEDIT(MEDIT)
0009
0010 C...Double precision and integer declarations.
0011       IMPLICIT DOUBLE PRECISION(A-H, O-Z)
0012       IMPLICIT INTEGER(I-N)
0013       INTEGER PYK,PYCHGE,PYCOMP
0014 C...Parameter statement to help give large particle numbers.
0015       PARAMETER (KSUSY1=1000000,KSUSY2=2000000,KTECHN=3000000,
0016      &KEXCIT=4000000,KDIMEN=5000000)
0017 C...Commonblocks.
0018       COMMON/PYJETS/N,NPAD,K(4000,5),P(4000,5),V(4000,5)
0019       COMMON/PYDAT1/MSTU(200),PARU(200),MSTJ(200),PARJ(200)
0020       COMMON/PYDAT2/KCHG(500,4),PMAS(500,4),PARF(2000),VCKM(4,4)
0021       COMMON/PYCTAG/NCT,MCT(4000,2)
0022       SAVE /PYJETS/,/PYDAT1/,/PYDAT2/,/PYCTAG/
0023 C...Local arrays.
0024       DIMENSION NS(2),PTS(2),PLS(2)
0025
0026 C...Remove unwanted partons/particles.
0027       IF((MEDIT.GE.0.AND.MEDIT.LE.3).OR.MEDIT.EQ.5) THEN
0028         IMAX=N
0029         IF(MSTU(2).GT.0) IMAX=MSTU(2)
0030         I1=MAX(1,MSTU(1))-1
0031         DO 110 I=MAX(1,MSTU(1)),IMAX
0032           IF(K(I,1).EQ.0.OR.(K(I,1).GE.21.AND.K(I,1).LE.40)) GOTO 110
0033           IF(MEDIT.EQ.1) THEN
0034             IF(K(I,1).GT.10.AND.K(I,1).NE.41.AND.K(I,1).NE.42) GOTO 110
0035           ELSEIF(MEDIT.EQ.2) THEN
0036             IF(K(I,1).GT.10.AND.K(I,1).NE.41.AND.K(I,1).NE.42) GOTO 110
0037             KC=PYCOMP(K(I,2))
0038             IF(KC.EQ.0.OR.KC.EQ.12.OR.KC.EQ.14.OR.KC.EQ.16.OR.
0039      &      KC.EQ.18.OR.K(I,2).EQ.KSUSY1+22.OR.K(I,2).EQ.39.OR.
0040      &      K(I,2).EQ.KSUSY1+39) GOTO 110
0041           ELSEIF(MEDIT.EQ.3) THEN
0042             IF(K(I,1).GT.10.AND.K(I,1).NE.41.AND.K(I,1).NE.42) GOTO 110
0043             KC=PYCOMP(K(I,2))
0044             IF(KC.EQ.0) GOTO 110
0045             IF(KCHG(KC,2).EQ.0.AND.PYCHGE(K(I,2)).EQ.0) GOTO 110
0046           ELSEIF(MEDIT.EQ.5) THEN
0047             IF(K(I,1).EQ.13.OR.K(I,1).EQ.14.OR.K(I,1).EQ.52) GOTO 110
0048             KC=PYCOMP(K(I,2))
0049             IF(KC.EQ.0) GOTO 110
0050             IF(K(I,1).GT.10.AND.K(I,1).NE.41.AND.K(I,1).NE.42.AND.
0051      &      KCHG(KC,2).EQ.0) GOTO 110
0052           ENDIF
0053
0054 C...Pack remaining partons/particles. Origin no longer known.
0055           I1=I1+1
0056           DO 100 J=1,5
0057             K(I1,J)=K(I,J)
0058             P(I1,J)=P(I,J)
0059             V(I1,J)=V(I,J)
0060   100     CONTINUE
0061           K(I1,3)=0
0062   110   CONTINUE
0063         IF(I1.LT.N) MSTU(3)=0
0064         IF(I1.LT.N) MSTU(70)=0
0065         N=I1
0066
0067 C...Selective removal of class of entries. New position of retained.
0068       ELSEIF(MEDIT.GE.11.AND.MEDIT.LE.15) THEN
0069         I1=0
0070         DO 120 I=1,N
0071           K(I,3)=MOD(K(I,3),MSTU(5))
0072           IF(MEDIT.EQ.11.AND.K(I,1).LT.0) GOTO 120
0073           IF(MEDIT.EQ.12.AND.K(I,1).EQ.0) GOTO 120
0074           IF(MEDIT.EQ.13.AND.(K(I,1).EQ.11.OR.K(I,1).EQ.12.OR.
0075      &    K(I,1).EQ.15.OR.K(I,1).EQ.51).AND.K(I,2).NE.94) GOTO 120
0076           IF(MEDIT.EQ.14.AND.(K(I,1).EQ.13.OR.K(I,1).EQ.14.OR.
0077      &    K(I,1).EQ.52.OR.K(I,2).EQ.94)) GOTO 120
0078           IF(MEDIT.EQ.15.AND.K(I,1).GE.21.AND.K(I,1).LE.40) GOTO 120
0079           I1=I1+1
0080           K(I,3)=K(I,3)+MSTU(5)*I1
0081   120   CONTINUE
0082
0083 C...Find new event history information and replace old.
0084         DO 140 I=1,N
0085           IF(K(I,1).LE.0.OR.(K(I,1).GE.21.AND.K(I,1).LE.40).OR.
0086      &    K(I,3)/MSTU(5).EQ.0) GOTO 140
0087           ID=I
0088   130     IM=MOD(K(ID,3),MSTU(5))
0089           IF(MEDIT.EQ.13.AND.IM.GT.0.AND.IM.LE.N) THEN
0090             IF((K(IM,1).EQ.11.OR.K(IM,1).EQ.12.OR.K(IM,1).EQ.15.OR.
0091      &      K(IM,1).EQ.51).AND.K(IM,2).NE.94) THEN
0092               ID=IM
0093               GOTO 130
0094             ENDIF
0095           ELSEIF(MEDIT.EQ.14.AND.IM.GT.0.AND.IM.LE.N) THEN
0096             IF(K(IM,1).EQ.13.OR.K(IM,1).EQ.14.OR.K(IM,1).EQ.52.OR.
0097      &      K(IM,2).EQ.94) THEN
0098               ID=IM
0099               GOTO 130
0100             ENDIF
0101           ENDIF
0102           K(I,3)=MSTU(5)*(K(I,3)/MSTU(5))
0103           IF(IM.NE.0) K(I,3)=K(I,3)+K(IM,3)/MSTU(5)
0104           IF(K(I,1).NE.3.AND.K(I,1).NE.13.AND.K(I,1).NE.14.AND.
0105      &      K(I,1).NE.42.AND.K(I,1).NE.52) THEN
0106             IF(K(I,4).GT.0.AND.K(I,4).LE.MSTU(4)) K(I,4)=
0107      &      K(K(I,4),3)/MSTU(5)
0108             IF(K(I,5).GT.0.AND.K(I,5).LE.MSTU(4)) K(I,5)=
0109      &      K(K(I,5),3)/MSTU(5)
0110           ELSE
0111             KCM=MOD(K(I,4)/MSTU(5),MSTU(5))
0112             IF(KCM.GT.0.AND.KCM.LE.MSTU(4).AND.K(I,1).NE.42.AND.
0113      &      K(I,1).NE.52) KCM=K(KCM,3)/MSTU(5)
0114             KCD=MOD(K(I,4),MSTU(5))
0115             IF(KCD.GT.0.AND.KCD.LE.MSTU(4)) KCD=K(KCD,3)/MSTU(5)
0116             K(I,4)=MSTU(5)**2*(K(I,4)/MSTU(5)**2)+MSTU(5)*KCM+KCD
0117             KCM=MOD(K(I,5)/MSTU(5),MSTU(5))
0118             IF(KCM.GT.0.AND.KCM.LE.MSTU(4)) KCM=K(KCM,3)/MSTU(5)
0119             KCD=MOD(K(I,5),MSTU(5))
0120             IF(KCD.GT.0.AND.KCD.LE.MSTU(4)) KCD=K(KCD,3)/MSTU(5)
0121             K(I,5)=MSTU(5)**2*(K(I,5)/MSTU(5)**2)+MSTU(5)*KCM+KCD
0122           ENDIF
0123   140   CONTINUE
0124
0125 C...Pack remaining entries.
0126         I1=0
0127         MSTU90=MSTU(90)
0128         MSTU(90)=0
0129         DO 170 I=1,N
0130           IF(K(I,3)/MSTU(5).EQ.0) GOTO 170
0131           I1=I1+1
0132           DO 150 J=1,5
0133             K(I1,J)=K(I,J)
0134             P(I1,J)=P(I,J)
0135             V(I1,J)=V(I,J)
0136   150     CONTINUE
0137 C...Also update LHA1 colour tags
0138           MCT(I1,1)=MCT(I,1)
0139           MCT(I1,2)=MCT(I,2)
0140           K(I1,3)=MOD(K(I1,3),MSTU(5))
0141           DO 160 IZ=1,MSTU90
0142             IF(I.EQ.MSTU(90+IZ)) THEN
0143               MSTU(90)=MSTU(90)+1
0144               MSTU(90+MSTU(90))=I1
0145               PARU(90+MSTU(90))=PARU(90+IZ)
0146             ENDIF
0147   160     CONTINUE
0148   170   CONTINUE
0149         IF(I1.LT.N) MSTU(3)=0
0150         IF(I1.LT.N) MSTU(70)=0
0151         N=I1
0152
0153 C...Fill in some missing daughter pointers (lost in colour flow).
0154       ELSEIF(MEDIT.EQ.16) THEN
0155         DO 220 I=1,N
0156           IF(K(I,1).LE.10.OR.(K(I,1).GE.21.AND.K(I,1).LE.50)) GOTO 220
0157           IF(K(I,4).NE.0.OR.K(I,5).NE.0) GOTO 220
0158 C...Find daughters who point to mother.
0159           DO 180 I1=I+1,N
0160             IF(K(I1,3).NE.I) THEN
0161             ELSEIF(K(I,4).EQ.0) THEN
0162               K(I,4)=I1
0163             ELSE
0164               K(I,5)=I1
0165             ENDIF
0166   180     CONTINUE
0167           IF(K(I,5).EQ.0) K(I,5)=K(I,4)
0168           IF(K(I,4).NE.0) GOTO 220
0169 C...Find daughters who point to documentation version of mother.
0170           IM=K(I,3)
0171           IF(IM.LE.0.OR.IM.GE.I) GOTO 220
0172           IF(K(IM,1).LE.20.OR.K(IM,1).GT.30) GOTO 220
0173           IF(K(IM,2).NE.K(I,2).OR.ABS(P(IM,5)-P(I,5)).GT.1D-2) GOTO 220
0174           DO 190 I1=I+1,N
0175             IF(K(I1,3).NE.IM) THEN
0176             ELSEIF(K(I,4).EQ.0) THEN
0177               K(I,4)=I1
0178             ELSE
0179               K(I,5)=I1
0180             ENDIF
0181   190     CONTINUE
0182           IF(K(I,5).EQ.0) K(I,5)=K(I,4)
0183           IF(K(I,4).NE.0) GOTO 220
0184 C...Find daughters who point to documentation daughters who,
0185 C...in their turn, point to documentation mother.
0186           ID1=IM
0187           ID2=IM
0188           DO 200 I1=IM+1,I-1
0189             IF(K(I1,3).EQ.IM.AND.K(I1,1).GE.21.AND.K(I1,1).LE.30) THEN
0190               ID2=I1
0191               IF(ID1.EQ.IM) ID1=I1
0192             ENDIF
0193   200     CONTINUE
0194           DO 210 I1=I+1,N
0195             IF(K(I1,3).NE.ID1.AND.K(I1,3).NE.ID2) THEN
0196             ELSEIF(K(I,4).EQ.0) THEN
0197               K(I,4)=I1
0198             ELSE
0199               K(I,5)=I1
0200             ENDIF
0201   210     CONTINUE
0202           IF(K(I,5).EQ.0) K(I,5)=K(I,4)
0203   220   CONTINUE
0204
0205 C...Save top entries at bottom of PYJETS commonblock.
0206       ELSEIF(MEDIT.EQ.21) THEN
0207         IF(2*N.GE.MSTU(4)) THEN
0208           CALL PYERRM(11,'(PYEDIT:) no more memory left in PYJETS')
0209           RETURN
0210         ENDIF
0211         DO 240 I=1,N
0212           DO 230 J=1,5
0213             K(MSTU(4)-I,J)=K(I,J)
0214             P(MSTU(4)-I,J)=P(I,J)
0215             V(MSTU(4)-I,J)=V(I,J)
0216   230     CONTINUE
0217   240   CONTINUE
0218         MSTU(32)=N
0219
0220 C...Restore bottom entries of commonblock PYJETS to top.
0221       ELSEIF(MEDIT.EQ.22) THEN
0222         DO 260 I=1,MSTU(32)
0223           DO 250 J=1,5
0224             K(I,J)=K(MSTU(4)-I,J)
0225             P(I,J)=P(MSTU(4)-I,J)
0226             V(I,J)=V(MSTU(4)-I,J)
0227   250     CONTINUE
0228   260   CONTINUE
0229         N=MSTU(32)
0230
0231 C...Mark primary entries at top of commonblock PYJETS as untreated.
0232       ELSEIF(MEDIT.EQ.23) THEN
0233         I1=0
0234         DO 270 I=1,N
0235           KH=K(I,3)
0236           IF(KH.GE.1) THEN
0237             IF(K(KH,1).GE.21.AND.K(KH,1).LE.30) KH=0
0238           ENDIF
0239           IF(KH.NE.0) GOTO 280
0240           I1=I1+1
0241           IF(K(I,1).GE.11.AND.K(I,1).LE.20) K(I,1)=K(I,1)-10
0242           IF(K(I,1).GE.51.AND.K(I,1).LE.60) K(I,1)=K(I,1)-10
0243   270   CONTINUE
0244   280   N=I1
0245
0246 C...Place largest axis along z axis and second largest in xy plane.
0247       ELSEIF(MEDIT.EQ.31.OR.MEDIT.EQ.32) THEN
0248         CALL PYROBO(1,N+MSTU(3),0D0,-PYANGL(P(MSTU(61),1),
0249      &  P(MSTU(61),2)),0D0,0D0,0D0)
0250         CALL PYROBO(1,N+MSTU(3),-PYANGL(P(MSTU(61),3),
0251      &  P(MSTU(61),1)),0D0,0D0,0D0,0D0)
0252         CALL PYROBO(1,N+MSTU(3),0D0,-PYANGL(P(MSTU(61)+1,1),
0253      &  P(MSTU(61)+1,2)),0D0,0D0,0D0)
0254         IF(MEDIT.EQ.31) RETURN
0255
0256 C...Rotate to put slim jet along +z axis.
0257         DO 290 IS=1,2
0258           NS(IS)=0
0259           PTS(IS)=0D0
0260           PLS(IS)=0D0
0261   290   CONTINUE
0262         DO 300 I=1,N
0263           IF(K(I,1).LE.0.OR.K(I,1).GT.10) GOTO 300
0264           IF(MSTU(41).GE.2) THEN
0265             KC=PYCOMP(K(I,2))
0266             IF(KC.EQ.0.OR.KC.EQ.12.OR.KC.EQ.14.OR.KC.EQ.16.OR.
0267      &      KC.EQ.18.OR.K(I,2).EQ.KSUSY1+22.OR.K(I,2).EQ.39.OR.
0268      &      K(I,2).EQ.KSUSY1+39) GOTO 300
0269             IF(MSTU(41).GE.3.AND.KCHG(KC,2).EQ.0.AND.PYCHGE(K(I,2))
0270      &      .EQ.0) GOTO 300
0271           ENDIF
0272           IS=2D0-SIGN(0.5D0,P(I,3))
0273           NS(IS)=NS(IS)+1
0274           PTS(IS)=PTS(IS)+SQRT(P(I,1)**2+P(I,2)**2)
0275   300   CONTINUE
0276         IF(NS(1)*PTS(2)**2.LT.NS(2)*PTS(1)**2)
0277      &  CALL PYROBO(1,N+MSTU(3),PARU(1),0D0,0D0,0D0,0D0)
0278
0279 C...Rotate to put second largest jet into -z,+x quadrant.
0280         DO 310 I=1,N
0281           IF(P(I,3).GE.0D0) GOTO 310
0282           IF(K(I,1).LE.0.OR.K(I,1).GT.10) GOTO 310
0283           IF(MSTU(41).GE.2) THEN
0284             KC=PYCOMP(K(I,2))
0285             IF(KC.EQ.0.OR.KC.EQ.12.OR.KC.EQ.14.OR.KC.EQ.16.OR.
0286      &      KC.EQ.18.OR.K(I,2).EQ.KSUSY1+22.OR.K(I,2).EQ.39.OR.
0287      &      K(I,2).EQ.KSUSY1+39) GOTO 310
0288             IF(MSTU(41).GE.3.AND.KCHG(KC,2).EQ.0.AND.PYCHGE(K(I,2))
0289      &      .EQ.0) GOTO 310
0290           ENDIF
0291           IS=2D0-SIGN(0.5D0,P(I,1))
0292           PLS(IS)=PLS(IS)-P(I,3)
0293   310   CONTINUE
0294         IF(PLS(2).GT.PLS(1)) CALL PYROBO(1,N+MSTU(3),0D0,PARU(1),
0295      &  0D0,0D0,0D0)
0296       ENDIF
0297
0298       RETURN
0299       END