HIJING/hipyset1_35/luprep_hijing.F

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