PYTHIA/jetset/luprep.F

   1
   2 C*********************************************************************
   3
   4       SUBROUTINE LUPREP(IP)
   5
   6 C...Purpose: to rearrange partons along strings, to allow small systems
   7 C...to collapse into one or two particles and to check flavours.
   8       IMPLICIT DOUBLE PRECISION(D)
   9       COMMON/LUJETS/N,K(4000,5),P(4000,5),V(4000,5)
  10       COMMON/LUDAT1/MSTU(200),PARU(200),MSTJ(200),PARJ(200)
  11       COMMON/LUDAT2/KCHG(500,3),PMAS(500,4),PARF(2000),VCKM(4,4)
  12       COMMON/LUDAT3/MDCY(500,3),MDME(2000,2),BRAT(2000),KFDP(2000,5)
  13       SAVE /LUJETS/,/LUDAT1/,/LUDAT2/,/LUDAT3/
  14       DIMENSION DPS(5),DPC(5),UE(3)
  15
  16 C...Rearrange parton shower product listing along strings: begin loop.
  17       I1=N
  18       DO 130 MQGST=1,2
  19       DO 120 I=MAX(1,IP),N
  20       IF(K(I,1).NE.3) GOTO 120
  21       KC=LUCOMP(K(I,2))
  22       IF(KC.EQ.0) GOTO 120
  23       KQ=KCHG(KC,2)
  24       IF(KQ.EQ.0.OR.(MQGST.EQ.1.AND.KQ.EQ.2)) GOTO 120
  25
  26 C...Pick up loose string end.
  27       KCS=4
  28       IF(KQ*ISIGN(1,K(I,2)).LT.0) KCS=5
  29       IA=I
  30       NSTP=0
  31   100 NSTP=NSTP+1
  32       IF(NSTP.GT.4*N) THEN
  33         CALL LUERRM(14,'(LUPREP:) caught in infinite loop')
  34         RETURN
  35       ENDIF
  36
  37 C...Copy undecayed parton.
  38       IF(K(IA,1).EQ.3) THEN
  39         IF(I1.GE.MSTU(4)-MSTU(32)-5) THEN
  40           CALL LUERRM(11,'(LUPREP:) no more memory left in LUJETS')
  41           RETURN
  42         ENDIF
  43         I1=I1+1
  44         K(I1,1)=2
  45         IF(NSTP.GE.2.AND.IABS(K(IA,2)).NE.21) K(I1,1)=1
  46         K(I1,2)=K(IA,2)
  47         K(I1,3)=IA
  48         K(I1,4)=0
  49         K(I1,5)=0
  50         DO 110 J=1,5
  51         P(I1,J)=P(IA,J)
  52         V(I1,J)=V(IA,J)
  53   110   CONTINUE
  54         K(IA,1)=K(IA,1)+10
  55         IF(K(I1,1).EQ.1) GOTO 120
  56       ENDIF
  57
  58 C...Go to next parton in colour space.
  59       IB=IA
  60       IF(MOD(K(IB,KCS)/MSTU(5)**2,2).EQ.0.AND.MOD(K(IB,KCS),MSTU(5))
  61      &.NE.0) THEN
  62         IA=MOD(K(IB,KCS),MSTU(5))
  63         K(IB,KCS)=K(IB,KCS)+MSTU(5)**2
  64         MREV=0
  65       ELSE
  66         IF(K(IB,KCS).GE.2*MSTU(5)**2.OR.MOD(K(IB,KCS)/MSTU(5),MSTU(5))
  67      &  .EQ.0) KCS=9-KCS
  68         IA=MOD(K(IB,KCS)/MSTU(5),MSTU(5))
  69         K(IB,KCS)=K(IB,KCS)+2*MSTU(5)**2
  70         MREV=1
  71       ENDIF
  72       IF(IA.LE.0.OR.IA.GT.N) THEN
  73         CALL LUERRM(12,'(LUPREP:) colour rearrangement failed')
  74         RETURN
  75       ENDIF
  76       IF(MOD(K(IA,4)/MSTU(5),MSTU(5)).EQ.IB.OR.MOD(K(IA,5)/MSTU(5),
  77      &MSTU(5)).EQ.IB) THEN
  78         IF(MREV.EQ.1) KCS=9-KCS
  79         IF(MOD(K(IA,KCS)/MSTU(5),MSTU(5)).NE.IB) KCS=9-KCS
  80         K(IA,KCS)=K(IA,KCS)+2*MSTU(5)**2
  81       ELSE
  82         IF(MREV.EQ.0) KCS=9-KCS
  83         IF(MOD(K(IA,KCS),MSTU(5)).NE.IB) KCS=9-KCS
  84         K(IA,KCS)=K(IA,KCS)+MSTU(5)**2
  85       ENDIF
  86       IF(IA.NE.I) GOTO 100
  87       K(I1,1)=1
  88   120 CONTINUE
  89   130 CONTINUE
  90       N=I1
  91       IF(MSTJ(14).LT.0) RETURN
  92
  93 C...Find lowest-mass colour singlet jet system, OK if above threshold.
  94       IF(MSTJ(14).EQ.0) GOTO 320
  95       NS=N
  96   140 NSIN=N-NS
  97       PDM=1.+PARJ(32)
  98       IC=0
  99       DO 190 I=MAX(1,IP),NS
 100       IF(K(I,1).NE.1.AND.K(I,1).NE.2) THEN
 101       ELSEIF(K(I,1).EQ.2.AND.IC.EQ.0) THEN
 102         NSIN=NSIN+1
 103         IC=I
 104         DO 150 J=1,4
 105         DPS(J)=P(I,J)
 106   150   CONTINUE
 107         MSTJ(93)=1
 108         DPS(5)=ULMASS(K(I,2))
 109       ELSEIF(K(I,1).EQ.2) THEN
 110         DO 160 J=1,4
 111         DPS(J)=DPS(J)+P(I,J)
 112   160   CONTINUE
 113       ELSEIF(IC.NE.0.AND.KCHG(LUCOMP(K(I,2)),2).NE.0) THEN
 114         DO 170 J=1,4
 115         DPS(J)=DPS(J)+P(I,J)
 116   170   CONTINUE
 117         MSTJ(93)=1
 118         DPS(5)=DPS(5)+ULMASS(K(I,2))
 119         PD=SQRT(MAX(0D0,DPS(4)**2-DPS(1)**2-DPS(2)**2-DPS(3)**2))-DPS(5)
 120         IF(PD.LT.PDM) THEN
 121           PDM=PD
 122           DO 180 J=1,5
 123           DPC(J)=DPS(J)
 124   180     CONTINUE
 125           IC1=IC
 126           IC2=I
 127         ENDIF
 128         IC=0
 129       ELSE
 130         NSIN=NSIN+1
 131       ENDIF
 132   190 CONTINUE
 133       IF(PDM.GE.PARJ(32)) GOTO 320
 134
 135 C...Fill small-mass system as cluster.
 136       NSAV=N
 137       PECM=SQRT(MAX(0D0,DPC(4)**2-DPC(1)**2-DPC(2)**2-DPC(3)**2))
 138       K(N+1,1)=11
 139       K(N+1,2)=91
 140       K(N+1,3)=IC1
 141       K(N+1,4)=N+2
 142       K(N+1,5)=N+3
 143       P(N+1,1)=DPC(1)
 144       P(N+1,2)=DPC(2)
 145       P(N+1,3)=DPC(3)
 146       P(N+1,4)=DPC(4)
 147       P(N+1,5)=PECM
 148
 149 C...Form two particles from flavours of lowest-mass system, if feasible.
 150       K(N+2,1)=1
 151       K(N+3,1)=1
 152       IF(MSTU(16).NE.2) THEN
 153         K(N+2,3)=N+1
 154         K(N+3,3)=N+1
 155       ELSE
 156         K(N+2,3)=IC1
 157         K(N+3,3)=IC2
 158       ENDIF
 159       K(N+2,4)=0
 160       K(N+3,4)=0
 161       K(N+2,5)=0
 162       K(N+3,5)=0
 163       IF(IABS(K(IC1,2)).NE.21) THEN
 164         KC1=LUCOMP(K(IC1,2))
 165         KC2=LUCOMP(K(IC2,2))
 166         IF(KC1.EQ.0.OR.KC2.EQ.0) GOTO 320
 167         KQ1=KCHG(KC1,2)*ISIGN(1,K(IC1,2))
 168         KQ2=KCHG(KC2,2)*ISIGN(1,K(IC2,2))
 169         IF(KQ1+KQ2.NE.0) GOTO 320
 170   200   CALL LUKFDI(K(IC1,2),0,KFLN,K(N+2,2))
 171         CALL LUKFDI(K(IC2,2),-KFLN,KFLDMP,K(N+3,2))
 172         IF(K(N+2,2).EQ.0.OR.K(N+3,2).EQ.0) GOTO 200
 173       ELSE
 174         IF(IABS(K(IC2,2)).NE.21) GOTO 320
 175   210   CALL LUKFDI(1+INT((2.+PARJ(2))*RLU(0)),0,KFLN,KFDMP)
 176         CALL LUKFDI(KFLN,0,KFLM,K(N+2,2))
 177         CALL LUKFDI(-KFLN,-KFLM,KFLDMP,K(N+3,2))
 178         IF(K(N+2,2).EQ.0.OR.K(N+3,2).EQ.0) GOTO 210
 179       ENDIF
 180       P(N+2,5)=ULMASS(K(N+2,2))
 181       P(N+3,5)=ULMASS(K(N+3,2))
 182       IF(P(N+2,5)+P(N+3,5)+PARJ(64).GE.PECM.AND.NSIN.EQ.1) GOTO 320
 183       IF(P(N+2,5)+P(N+3,5)+PARJ(64).GE.PECM) GOTO 260
 184
 185 C...Perform two-particle decay of jet system, if possible.
 186       IF(PECM.GE.0.02*DPC(4)) THEN
 187         PA=SQRT((PECM**2-(P(N+2,5)+P(N+3,5))**2)*(PECM**2-
 188      &  (P(N+2,5)-P(N+3,5))**2))/(2.*PECM)
 189         UE(3)=2.*RLU(0)-1.
 190         PHI=PARU(2)*RLU(0)
 191         UE(1)=SQRT(1.-UE(3)**2)*COS(PHI)
 192         UE(2)=SQRT(1.-UE(3)**2)*SIN(PHI)
 193         DO 220 J=1,3
 194         P(N+2,J)=PA*UE(J)
 195         P(N+3,J)=-PA*UE(J)
 196   220   CONTINUE
 197         P(N+2,4)=SQRT(PA**2+P(N+2,5)**2)
 198         P(N+3,4)=SQRT(PA**2+P(N+3,5)**2)
 199         MSTU(33)=1
 200         CALL LUDBRB(N+2,N+3,0.,0.,DPC(1)/DPC(4),DPC(2)/DPC(4),
 201      &  DPC(3)/DPC(4))
 202       ELSE
 203         NP=0
 204         DO 230 I=IC1,IC2
 205         IF(K(I,1).EQ.1.OR.K(I,1).EQ.2) NP=NP+1
 206   230   CONTINUE
 207         HA=P(IC1,4)*P(IC2,4)-P(IC1,1)*P(IC2,1)-P(IC1,2)*P(IC2,2)-
 208      &  P(IC1,3)*P(IC2,3)
 209         IF(NP.GE.3.OR.HA.LE.1.25*P(IC1,5)*P(IC2,5)) GOTO 260
 210         HD1=0.5*(P(N+2,5)**2-P(IC1,5)**2)
 211         HD2=0.5*(P(N+3,5)**2-P(IC2,5)**2)
 212         HR=SQRT(MAX(0.,((HA-HD1-HD2)**2-(P(N+2,5)*P(N+3,5))**2)/
 213      &  (HA**2-(P(IC1,5)*P(IC2,5))**2)))-1.
 214         HC=P(IC1,5)**2+2.*HA+P(IC2,5)**2
 215         HK1=((P(IC2,5)**2+HA)*HR+HD1-HD2)/HC
 216         HK2=((P(IC1,5)**2+HA)*HR+HD2-HD1)/HC
 217         DO 240 J=1,4
 218         P(N+2,J)=(1.+HK1)*P(IC1,J)-HK2*P(IC2,J)
 219         P(N+3,J)=(1.+HK2)*P(IC2,J)-HK1*P(IC1,J)
 220   240   CONTINUE
 221       ENDIF
 222       DO 250 J=1,4
 223       V(N+1,J)=V(IC1,J)
 224       V(N+2,J)=V(IC1,J)
 225       V(N+3,J)=V(IC2,J)
 226   250 CONTINUE
 227       V(N+1,5)=0.
 228       V(N+2,5)=0.
 229       V(N+3,5)=0.
 230       N=N+3
 231       GOTO 300
 232
 233 C...Else form one particle from the flavours available, if possible.
 234   260 K(N+1,5)=N+2
 235       IF(IABS(K(IC1,2)).GT.100.AND.IABS(K(IC2,2)).GT.100) THEN
 236         GOTO 320
 237       ELSEIF(IABS(K(IC1,2)).NE.21) THEN
 238         CALL LUKFDI(K(IC1,2),K(IC2,2),KFLDMP,K(N+2,2))
 239       ELSE
 240         KFLN=1+INT((2.+PARJ(2))*RLU(0))
 241         CALL LUKFDI(KFLN,-KFLN,KFLDMP,K(N+2,2))
 242       ENDIF
 243       IF(K(N+2,2).EQ.0) GOTO 260
 244       P(N+2,5)=ULMASS(K(N+2,2))
 245
 246 C...Find parton/particle which combines to largest extra mass.
 247       IR=0
 248       HA=0.
 249       HSM=0.
 250       DO 280 MCOMB=1,3
 251       IF(IR.NE.0) GOTO 280
 252       DO 270 I=MAX(1,IP),N
 253       IF(K(I,1).LE.0.OR.K(I,1).GT.10.OR.(I.GE.IC1.AND.I.LE.IC2
 254      &.AND.K(I,1).GE.1.AND.K(I,1).LE.2)) GOTO 270
 255       IF(MCOMB.EQ.1) KCI=LUCOMP(K(I,2))
 256       IF(MCOMB.EQ.1.AND.KCI.EQ.0) GOTO 270
 257       IF(MCOMB.EQ.1.AND.KCHG(KCI,2).EQ.0.AND.I.LE.NS) GOTO 270
 258       IF(MCOMB.EQ.2.AND.IABS(K(I,2)).GT.10.AND.IABS(K(I,2)).LE.100)
 259      &GOTO 270
 260       HCR=DPC(4)*P(I,4)-DPC(1)*P(I,1)-DPC(2)*P(I,2)-DPC(3)*P(I,3)
 261       HSR=2.*HCR+PECM**2-P(N+2,5)**2-2.*P(N+2,5)*P(I,5)
 262       IF(HSR.GT.HSM) THEN
 263         IR=I
 264         HA=HCR
 265         HSM=HSR
 266       ENDIF
 267   270 CONTINUE
 268   280 CONTINUE
 269
 270 C...Shuffle energy and momentum to put new particle on mass shell.
 271       IF(IR.NE.0) THEN
 272         HB=PECM**2+HA
 273         HC=P(N+2,5)**2+HA
 274         HD=P(IR,5)**2+HA
 275         HK2=0.5*(HB*SQRT(MAX(0.,((HB+HC)**2-4.*(HB+HD)*P(N+2,5)**2)/
 276      &  (HA**2-(PECM*P(IR,5))**2)))-(HB+HC))/(HB+HD)
 277         HK1=(0.5*(P(N+2,5)**2-PECM**2)+HD*HK2)/HB
 278         DO 290 J=1,4
 279         P(N+2,J)=(1.+HK1)*DPC(J)-HK2*P(IR,J)
 280         P(IR,J)=(1.+HK2)*P(IR,J)-HK1*DPC(J)
 281         V(N+1,J)=V(IC1,J)
 282         V(N+2,J)=V(IC1,J)
 283   290   CONTINUE
 284         V(N+1,5)=0.
 285         V(N+2,5)=0.
 286         N=N+2
 287       ELSE
 288         CALL LUERRM(3,'(LUPREP:) no match for collapsing cluster')
 289         RETURN
 290       ENDIF
 291
 292 C...Mark collapsed system and store daughter pointers. Iterate.
 293   300 DO 310 I=IC1,IC2
 294       IF((K(I,1).EQ.1.OR.K(I,1).EQ.2).AND.KCHG(LUCOMP(K(I,2)),2).NE.0)
 295      &THEN
 296         K(I,1)=K(I,1)+10
 297         IF(MSTU(16).NE.2) THEN
 298           K(I,4)=NSAV+1
 299           K(I,5)=NSAV+1
 300         ELSE
 301           K(I,4)=NSAV+2
 302           K(I,5)=N
 303         ENDIF
 304       ENDIF
 305   310 CONTINUE
 306       IF(N.LT.MSTU(4)-MSTU(32)-5) GOTO 140
 307
 308 C...Check flavours and invariant masses in parton systems.
 309   320 NP=0
 310       KFN=0
 311       KQS=0
 312       DO 330 J=1,5
 313       DPS(J)=0.
 314   330 CONTINUE
 315       DO 360 I=MAX(1,IP),N
 316       IF(K(I,1).LE.0.OR.K(I,1).GT.10) GOTO 360
 317       KC=LUCOMP(K(I,2))
 318       IF(KC.EQ.0) GOTO 360
 319       KQ=KCHG(KC,2)*ISIGN(1,K(I,2))
 320       IF(KQ.EQ.0) GOTO 360
 321       NP=NP+1
 322       IF(KQ.NE.2) THEN
 323         KFN=KFN+1
 324         KQS=KQS+KQ
 325         MSTJ(93)=1
 326         DPS(5)=DPS(5)+ULMASS(K(I,2))
 327       ENDIF
 328       DO 340 J=1,4
 329       DPS(J)=DPS(J)+P(I,J)
 330   340 CONTINUE
 331       IF(K(I,1).EQ.1) THEN
 332         IF(NP.NE.1.AND.(KFN.EQ.1.OR.KFN.GE.3.OR.KQS.NE.0)) CALL
 333      &  LUERRM(2,'(LUPREP:) unphysical flavour combination')
 334         IF(NP.NE.1.AND.DPS(4)**2-DPS(1)**2-DPS(2)**2-DPS(3)**2.LT.
 335      &  (0.9*PARJ(32)+DPS(5))**2) CALL LUERRM(3,
 336      &  '(LUPREP:) too small mass in jet system')
 337         NP=0
 338         KFN=0
 339         KQS=0
 340         DO 350 J=1,5
 341         DPS(J)=0.
 342   350   CONTINUE
 343       ENDIF
 344   360 CONTINUE
 345
 346       RETURN
 347       END