hijing/src/luprep.f

0001
0002 C*********************************************************************
0003
0004       SUBROUTINE LUPREP(IP)
0005
0006 C...Purpose: to rearrange partons along strings, to allow small systems
0007 C...to collapse into one or two particles and to check flavours.
0008       IMPLICIT DOUBLE PRECISION(D)
0009       COMMON/LUJETS/N,K(9000,5),P(9000,5),V(9000,5)
0010       SAVE /LUJETS/
0011       COMMON/LUDAT1/MSTU(200),PARU(200),MSTJ(200),PARJ(200)
0012       SAVE /LUDAT1/
0013       COMMON/LUDAT2/KCHG(500,3),PMAS(500,4),PARF(2000),VCKM(4,4)
0014       SAVE /LUDAT2/
0015       COMMON/LUDAT3/MDCY(500,3),MDME(2000,2),BRAT(2000),KFDP(2000,5)
0016       SAVE /LUDAT3/
0017       DIMENSION DPS(5),DPC(5),UE(3)
0018
0019 C...Rearrange parton shower product listing along strings: begin loop.
0020       I1=N
0021       DO 130 MQGST=1,2
0022       DO 120 I=MAX(1,IP),N
0023       IF(K(I,1).NE.3) GOTO 120
0024       KC=LUCOMP(K(I,2))
0025       IF(KC.EQ.0) GOTO 120
0026       KQ=KCHG(KC,2)
0027       IF(KQ.EQ.0.OR.(MQGST.EQ.1.AND.KQ.EQ.2)) GOTO 120
0028
0029 C...Pick up loose string end.
0030       KCS=4
0031       IF(KQ*ISIGN(1,K(I,2)).LT.0) KCS=5
0032       IA=I
0033       NSTP=0
0034   100 NSTP=NSTP+1
0035       IF(NSTP.GT.4*N) THEN
0036         CALL LUERRM(14,'(LUPREP:) caught in infinite loop')
0037         RETURN
0038       ENDIF
0039
0040 C...Copy undecayed parton.
0041       IF(K(IA,1).EQ.3) THEN
0042         IF(I1.GE.MSTU(4)-MSTU(32)-5) THEN
0043           CALL LUERRM(11,'(LUPREP:) no more memory left in LUJETS')
0044           RETURN
0045         ENDIF
0046         I1=I1+1
0047         K(I1,1)=2
0048         IF(NSTP.GE.2.AND.IABS(K(IA,2)).NE.21) K(I1,1)=1
0049         K(I1,2)=K(IA,2)
0050         K(I1,3)=IA
0051         K(I1,4)=0
0052         K(I1,5)=0
0053         DO 110 J=1,5
0054         P(I1,J)=P(IA,J)
0055   110   V(I1,J)=V(IA,J)
0056         K(IA,1)=K(IA,1)+10
0057         IF(K(I1,1).EQ.1) GOTO 120
0058       ENDIF
0059
0060 C...Go to next parton in colour space.
0061       IB=IA
0062       IF(MOD(K(IB,KCS)/MSTU(5)**2,2).EQ.0.AND.MOD(K(IB,KCS),MSTU(5)).
0063      &NE.0) THEN
0064         IA=MOD(K(IB,KCS),MSTU(5))
0065         K(IB,KCS)=K(IB,KCS)+MSTU(5)**2
0066         MREV=0
0067       ELSE
0068         IF(K(IB,KCS).GE.2*MSTU(5)**2.OR.MOD(K(IB,KCS)/MSTU(5),MSTU(5)).
0069      &  EQ.0) KCS=9-KCS
0070         IA=MOD(K(IB,KCS)/MSTU(5),MSTU(5))
0071         K(IB,KCS)=K(IB,KCS)+2*MSTU(5)**2
0072         MREV=1
0073       ENDIF
0074       IF(IA.LE.0.OR.IA.GT.N) THEN
0075         CALL LUERRM(12,'(LUPREP:) colour rearrangement failed')
0076         RETURN
0077       ENDIF
0078       IF(MOD(K(IA,4)/MSTU(5),MSTU(5)).EQ.IB.OR.MOD(K(IA,5)/MSTU(5),
0079      &MSTU(5)).EQ.IB) THEN
0080         IF(MREV.EQ.1) KCS=9-KCS
0081         IF(MOD(K(IA,KCS)/MSTU(5),MSTU(5)).NE.IB) KCS=9-KCS
0082         K(IA,KCS)=K(IA,KCS)+2*MSTU(5)**2
0083       ELSE
0084         IF(MREV.EQ.0) KCS=9-KCS
0085         IF(MOD(K(IA,KCS),MSTU(5)).NE.IB) KCS=9-KCS
0086         K(IA,KCS)=K(IA,KCS)+MSTU(5)**2
0087       ENDIF
0088       IF(IA.NE.I) GOTO 100
0089       K(I1,1)=1
0090   120 CONTINUE
0091   130 CONTINUE
0092       N=I1
0093
0094 C...Find lowest-mass colour singlet jet system, OK if above threshold.
0095       IF(MSTJ(14).LE.0) GOTO 320
0096       NS=N
0097   140 NSIN=N-NS
0098       PDM=1.+PARJ(32)
0099       IC=0
0100       DO 190 I=MAX(1,IP),NS
0101       IF(K(I,1).NE.1.AND.K(I,1).NE.2) THEN
0102       ELSEIF(K(I,1).EQ.2.AND.IC.EQ.0) THEN
0103         NSIN=NSIN+1
0104         IC=I
0105         DO 150 J=1,4
0106   150   DPS(J)=P(I,J)
0107         MSTJ(93)=1
0108         DPS(5)=ULMASS(K(I,2))
0109       ELSEIF(K(I,1).EQ.2) THEN
0110         DO 160 J=1,4
0111   160   DPS(J)=DPS(J)+P(I,J)
0112       ELSEIF(IC.NE.0.AND.KCHG(LUCOMP(K(I,2)),2).NE.0) THEN
0113         DO 170 J=1,4
0114   170   DPS(J)=DPS(J)+P(I,J)
0115         MSTJ(93)=1
0116         DPS(5)=DPS(5)+ULMASS(K(I,2))
0117         PD=SQRT(MAX(0D0,DPS(4)**2-DPS(1)**2-DPS(2)**2-DPS(3)**2))-DPS(5)
0118         IF(PD.LT.PDM) THEN
0119           PDM=PD
0120           DO 180 J=1,5
0121   180     DPC(J)=DPS(J)
0122           IC1=IC
0123           IC2=I
0124         ENDIF
0125         IC=0
0126       ELSE
0127         NSIN=NSIN+1
0128       ENDIF
0129   190 CONTINUE
0130       IF(PDM.GE.PARJ(32)) GOTO 320
0131
0132 C...Fill small-mass system as cluster.
0133       NSAV=N
0134       PECM=SQRT(MAX(0D0,DPC(4)**2-DPC(1)**2-DPC(2)**2-DPC(3)**2))
0135       K(N+1,1)=11
0136       K(N+1,2)=91
0137       K(N+1,3)=IC1
0138       K(N+1,4)=N+2
0139       K(N+1,5)=N+3
0140       P(N+1,1)=DPC(1)
0141       P(N+1,2)=DPC(2)
0142       P(N+1,3)=DPC(3)
0143       P(N+1,4)=DPC(4)
0144       P(N+1,5)=PECM
0145
0146 C...Form two particles from flavours of lowest-mass system, if feasible.
0147       K(N+2,1)=1
0148       K(N+3,1)=1
0149       IF(MSTU(16).NE.2) THEN
0150         K(N+2,3)=N+1
0151         K(N+3,3)=N+1
0152       ELSE
0153         K(N+2,3)=IC1
0154         K(N+3,3)=IC2
0155       ENDIF
0156       K(N+2,4)=0
0157       K(N+3,4)=0
0158       K(N+2,5)=0
0159       K(N+3,5)=0
0160       IF(IABS(K(IC1,2)).NE.21) THEN
0161         KC1=LUCOMP(K(IC1,2))
0162         KC2=LUCOMP(K(IC2,2))
0163         IF(KC1.EQ.0.OR.KC2.EQ.0) GOTO 320
0164         KQ1=KCHG(KC1,2)*ISIGN(1,K(IC1,2))
0165         KQ2=KCHG(KC2,2)*ISIGN(1,K(IC2,2))
0166         IF(KQ1+KQ2.NE.0) GOTO 320
0167   200   CALL LUKFDI(K(IC1,2),0,KFLN,K(N+2,2))
0168         CALL LUKFDI(K(IC2,2),-KFLN,KFLDMP,K(N+3,2))
0169         IF(K(N+2,2).EQ.0.OR.K(N+3,2).EQ.0) GOTO 200
0170       ELSE
0171         IF(IABS(K(IC2,2)).NE.21) GOTO 320
0172   210   CALL LUKFDI(1+INT((2.+PARJ(2))*RLU(0)),0,KFLN,KFDMP)
0173         CALL LUKFDI(KFLN,0,KFLM,K(N+2,2))
0174         CALL LUKFDI(-KFLN,-KFLM,KFLDMP,K(N+3,2))
0175         IF(K(N+2,2).EQ.0.OR.K(N+3,2).EQ.0) GOTO 210
0176       ENDIF
0177       P(N+2,5)=ULMASS(K(N+2,2))
0178       P(N+3,5)=ULMASS(K(N+3,2))
0179       IF(P(N+2,5)+P(N+3,5)+PARJ(64).GE.PECM.AND.NSIN.EQ.1) GOTO 320
0180       IF(P(N+2,5)+P(N+3,5)+PARJ(64).GE.PECM) GOTO 260
0181
0182 C...Perform two-particle decay of jet system, if possible.
0183       IF(PECM.GE.0.02*DPC(4)) THEN
0184         PA=SQRT((PECM**2-(P(N+2,5)+P(N+3,5))**2)*(PECM**2-
0185      &  (P(N+2,5)-P(N+3,5))**2))/(2.*PECM)
0186         UE(3)=2.*RLU(0)-1.
0187         PHI=PARU(2)*RLU(0)
0188         UE(1)=SQRT(1.-UE(3)**2)*COS(PHI)
0189         UE(2)=SQRT(1.-UE(3)**2)*SIN(PHI)
0190         DO 220 J=1,3
0191         P(N+2,J)=PA*UE(J)
0192   220   P(N+3,J)=-PA*UE(J)
0193         P(N+2,4)=SQRT(PA**2+P(N+2,5)**2)
0194         P(N+3,4)=SQRT(PA**2+P(N+3,5)**2)
0195         CALL LUDBRB(N+2,N+3,0.,0.,DPC(1)/DPC(4),DPC(2)/DPC(4),
0196      &  DPC(3)/DPC(4))
0197       ELSE
0198         NP=0
0199         DO 230 I=IC1,IC2
0200   230   IF(K(I,1).EQ.1.OR.K(I,1).EQ.2) NP=NP+1
0201         HA=P(IC1,4)*P(IC2,4)-P(IC1,1)*P(IC2,1)-P(IC1,2)*P(IC2,2)-
0202      &  P(IC1,3)*P(IC2,3)
0203         IF(NP.GE.3.OR.HA.LE.1.25*P(IC1,5)*P(IC2,5)) GOTO 260
0204         HD1=0.5*(P(N+2,5)**2-P(IC1,5)**2)
0205         HD2=0.5*(P(N+3,5)**2-P(IC2,5)**2)
0206         HR=SQRT(MAX(0.,((HA-HD1-HD2)**2-(P(N+2,5)*P(N+3,5))**2)/
0207      &  (HA**2-(P(IC1,5)*P(IC2,5))**2)))-1.
0208         HC=P(IC1,5)**2+2.*HA+P(IC2,5)**2
0209         HK1=((P(IC2,5)**2+HA)*HR+HD1-HD2)/HC
0210         HK2=((P(IC1,5)**2+HA)*HR+HD2-HD1)/HC
0211         DO 240 J=1,4
0212         P(N+2,J)=(1.+HK1)*P(IC1,J)-HK2*P(IC2,J)
0213   240   P(N+3,J)=(1.+HK2)*P(IC2,J)-HK1*P(IC1,J)
0214       ENDIF
0215       DO 250 J=1,4
0216       V(N+1,J)=V(IC1,J)
0217       V(N+2,J)=V(IC1,J)
0218   250 V(N+3,J)=V(IC2,J)
0219       V(N+1,5)=0.
0220       V(N+2,5)=0.
0221       V(N+3,5)=0.
0222       N=N+3
0223       GOTO 300
0224
0225 C...Else form one particle from the flavours available, if possible.
0226   260 K(N+1,5)=N+2
0227       IF(IABS(K(IC1,2)).GT.100.AND.IABS(K(IC2,2)).GT.100) THEN
0228         GOTO 320
0229       ELSEIF(IABS(K(IC1,2)).NE.21) THEN
0230         CALL LUKFDI(K(IC1,2),K(IC2,2),KFLDMP,K(N+2,2))
0231       ELSE
0232         KFLN=1+INT((2.+PARJ(2))*RLU(0))
0233         CALL LUKFDI(KFLN,-KFLN,KFLDMP,K(N+2,2))
0234       ENDIF
0235       IF(K(N+2,2).EQ.0) GOTO 260
0236       P(N+2,5)=ULMASS(K(N+2,2))
0237
0238 C...Find parton/particle which combines to largest extra mass.
0239       IR=0
0240       HA=0.
0241       DO 280 MCOMB=1,3
0242       IF(IR.NE.0) GOTO 280
0243       DO 270 I=MAX(1,IP),N
0244       IF(K(I,1).LE.0.OR.K(I,1).GT.10.OR.(I.GE.IC1.AND.I.LE.IC2.
0245      &AND.K(I,1).GE.1.AND.K(I,1).LE.2)) GOTO 270
0246       IF(MCOMB.EQ.1) KCI=LUCOMP(K(I,2))
0247       IF(MCOMB.EQ.1.AND.KCI.EQ.0) GOTO 270
0248       IF(MCOMB.EQ.1.AND.KCHG(KCI,2).EQ.0.AND.I.LE.NS) GOTO 270
0249       IF(MCOMB.EQ.2.AND.IABS(K(I,2)).GT.10.AND.IABS(K(I,2)).LE.100)
0250      &GOTO 270
0251       HCR=DPC(4)*P(I,4)-DPC(1)*P(I,1)-DPC(2)*P(I,2)-DPC(3)*P(I,3)
0252       IF(HCR.GT.HA) THEN
0253         IR=I
0254         HA=HCR
0255       ENDIF
0256   270 CONTINUE
0257   280 CONTINUE
0258
0259 C...Shuffle energy and momentum to put new particle on mass shell.
0260       HB=PECM**2+HA
0261       HC=P(N+2,5)**2+HA
0262       HD=P(IR,5)**2+HA
0263 C******************CHANGES BY HIJING************
0264       HK2=0.0
0265       IF(HA**2-(PECM*P(IR,5))**2.EQ.0.0.OR.HB+HD.EQ.0.0) GO TO 285
0266 C******************
0267       HK2=0.5*(HB*SQRT(((HB+HC)**2-4.*(HB+HD)*P(N+2,5)**2)/
0268      &(HA**2-(PECM*P(IR,5))**2))-(HB+HC))/(HB+HD)
0269   285 HK1=(0.5*(P(N+2,5)**2-PECM**2)+HD*HK2)/HB
0270       DO 290 J=1,4
0271       P(N+2,J)=(1.+HK1)*DPC(J)-HK2*P(IR,J)
0272       P(IR,J)=(1.+HK2)*P(IR,J)-HK1*DPC(J)
0273       V(N+1,J)=V(IC1,J)
0274   290 V(N+2,J)=V(IC1,J)
0275       V(N+1,5)=0.
0276       V(N+2,5)=0.
0277       N=N+2
0278
0279 C...Mark collapsed system and store daughter pointers. Iterate.
0280   300 DO 310 I=IC1,IC2
0281       IF((K(I,1).EQ.1.OR.K(I,1).EQ.2).AND.KCHG(LUCOMP(K(I,2)),2).NE.0)
0282      &THEN
0283         K(I,1)=K(I,1)+10
0284         IF(MSTU(16).NE.2) THEN
0285           K(I,4)=NSAV+1
0286           K(I,5)=NSAV+1
0287         ELSE
0288           K(I,4)=NSAV+2
0289           K(I,5)=N
0290         ENDIF
0291       ENDIF
0292   310 CONTINUE
0293       IF(N.LT.MSTU(4)-MSTU(32)-5) GOTO 140
0294
0295 C...Check flavours and invariant masses in parton systems.
0296   320 NP=0
0297       KFN=0
0298       KQS=0
0299       DO 330 J=1,5
0300   330 DPS(J)=0.
0301       DO 360 I=MAX(1,IP),N
0302       IF(K(I,1).LE.0.OR.K(I,1).GT.10) GOTO 360
0303       KC=LUCOMP(K(I,2))
0304       IF(KC.EQ.0) GOTO 360
0305       KQ=KCHG(KC,2)*ISIGN(1,K(I,2))
0306       IF(KQ.EQ.0) GOTO 360
0307       NP=NP+1
0308       IF(KQ.NE.2) THEN
0309         KFN=KFN+1
0310         KQS=KQS+KQ
0311         MSTJ(93)=1
0312         DPS(5)=DPS(5)+ULMASS(K(I,2))
0313       ENDIF
0314       DO 340 J=1,4
0315   340 DPS(J)=DPS(J)+P(I,J)
0316       IF(K(I,1).EQ.1) THEN
0317         IF(NP.NE.1.AND.(KFN.EQ.1.OR.KFN.GE.3.OR.KQS.NE.0)) CALL
0318      &  LUERRM(2,'(LUPREP:) unphysical flavour combination')
0319         IF(NP.NE.1.AND.DPS(4)**2-DPS(1)**2-DPS(2)**2-DPS(3)**2.LT.
0320      &  (0.9*PARJ(32)+DPS(5))**2) CALL LUERRM(3,
0321      &  '(LUPREP:) too small mass in jet system')
0322         NP=0
0323         KFN=0
0324         KQS=0
0325         DO 350 J=1,5
0326   350   DPS(J)=0.
0327       ENDIF
0328   360 CONTINUE
0329
0330       RETURN
0331       END