HIJING/hipyset1_35/luclus_hijing.F

   1 * $Id$
   2
   3 C*********************************************************************
   4
   5       SUBROUTINE LUCLUS_HIJING(NJET)
   6
   7 C...Purpose: to subdivide the particle content of an event into
   8 C...jets/clusters.
   9 #include "lujets_hijing.inc"
  10 #include "ludat1_hijing.inc"
  11 #include "ludat2_hijing.inc"
  12       DIMENSION PS(5)
  13       SAVE NSAV,NP,PS,PSS,RINIT,NPRE,NREM
  14
  15 C...Functions: distance measure in pT or (pseudo)mass.
  16       R2T(I1,I2)=(P(I1,5)*P(I2,5)-P(I1,1)*P(I2,1)-P(I1,2)*P(I2,2)-
  17      &P(I1,3)*P(I2,3))*2.*P(I1,5)*P(I2,5)/(0.0001+P(I1,5)+P(I2,5))**2
  18       R2M(I1,I2)=2.*P(I1,4)*P(I2,4)*(1.-(P(I1,1)*P(I2,1)+P(I1,2)*
  19      &P(I2,2)+P(I1,3)*P(I2,3))/(P(I1,5)*P(I2,5)))
  20
  21 C...If first time, reset. If reentering, skip preliminaries.
  22       IF(MSTU(48).LE.0) THEN
  23         NP=0
  24         DO 100 J=1,5
  25   100   PS(J)=0.
  26         PSS=0.
  27       ELSE
  28         NJET=NSAV
  29         IF(MSTU(43).GE.2) N=N-NJET
  30         DO 110 I=N+1,N+NJET
  31   110   P(I,5)=SQRT(P(I,1)**2+P(I,2)**2+P(I,3)**2)
  32         IF(MSTU(46).LE.3) R2ACC=PARU(44)**2
  33         IF(MSTU(46).GE.4) R2ACC=PARU(45)*PS(5)**2
  34         NLOOP=0
  35         GOTO 290
  36       ENDIF
  37
  38 C...Find which particles are to be considered in cluster search.
  39       DO 140 I=1,N
  40       IF(K(I,1).LE.0.OR.K(I,1).GT.10) GOTO 140
  41       IF(MSTU(41).GE.2) THEN
  42         KC=LUCOMP_HIJING(K(I,2))
  43         IF(KC.EQ.0.OR.KC.EQ.12.OR.KC.EQ.14.OR.KC.EQ.16.OR.
  44      &  KC.EQ.18) GOTO 140
  45         IF(MSTU(41).GE.3.AND.KCHG(KC,2).EQ.0.AND.LUCHGE_HIJING(K(I,2))
  46      $       .EQ.0)GOTO 140
  47       ENDIF
  48       IF(N+2*NP.GE.MSTU(4)-MSTU(32)-5) THEN
  49          CALL LUERRM_HIJING(11
  50      $        ,'(LUCLUS_HIJING:) no more memory left in LUJETS_HIJING')
  51         NJET=-1
  52         RETURN
  53       ENDIF
  54
  55 C...Take copy of these particles, with space left for jets later on.
  56       NP=NP+1
  57       K(N+NP,3)=I
  58       DO 120 J=1,5
  59   120 P(N+NP,J)=P(I,J)
  60       IF(MSTU(42).EQ.0) P(N+NP,5)=0.
  61       IF(MSTU(42).EQ.1.AND.K(I,2).NE.22) P(N+NP,5)=PMAS(101,1)
  62       P(N+NP,4)=SQRT(P(N+NP,5)**2+P(I,1)**2+P(I,2)**2+P(I,3)**2)
  63       P(N+NP,5)=SQRT(P(I,1)**2+P(I,2)**2+P(I,3)**2)
  64       DO 130 J=1,4
  65   130 PS(J)=PS(J)+P(N+NP,J)
  66       PSS=PSS+P(N+NP,5)
  67   140 CONTINUE
  68       DO 150 I=N+1,N+NP
  69       K(I+NP,3)=K(I,3)
  70       DO 150 J=1,5
  71   150 P(I+NP,J)=P(I,J)
  72       PS(5)=SQRT(MAX(0.,PS(4)**2-PS(1)**2-PS(2)**2-PS(3)**2))
  73
  74 C...Very low multiplicities not considered.
  75       IF(NP.LT.MSTU(47)) THEN
  76          CALL LUERRM_HIJING(8
  77      $        ,'(LUCLUS_HIJING:) too few particles for analysis')
  78         NJET=-1
  79         RETURN
  80       ENDIF
  81
  82 C...Find precluster configuration. If too few jets, make harder cuts.
  83       NLOOP=0
  84       IF(MSTU(46).LE.3) R2ACC=PARU(44)**2
  85       IF(MSTU(46).GE.4) R2ACC=PARU(45)*PS(5)**2
  86       RINIT=1.25*PARU(43)
  87       IF(NP.LE.MSTU(47)+2) RINIT=0.
  88   160 RINIT=0.8*RINIT
  89       NPRE=0
  90       NREM=NP
  91       DO 170 I=N+NP+1,N+2*NP
  92   170 K(I,4)=0
  93
  94 C...Sum up small momentum region. Jet if enough absolute momentum.
  95       IF(MSTU(46).LE.2) THEN
  96         DO 180 J=1,4
  97   180   P(N+1,J)=0.
  98         DO 200 I=N+NP+1,N+2*NP
  99         IF(P(I,5).GT.2.*RINIT) GOTO 200
 100         NREM=NREM-1
 101         K(I,4)=1
 102         DO 190 J=1,4
 103   190   P(N+1,J)=P(N+1,J)+P(I,J)
 104   200   CONTINUE
 105         P(N+1,5)=SQRT(P(N+1,1)**2+P(N+1,2)**2+P(N+1,3)**2)
 106         IF(P(N+1,5).GT.2.*RINIT) NPRE=1
 107         IF(RINIT.GE.0.2*PARU(43).AND.NPRE+NREM.LT.MSTU(47)) GOTO 160
 108       ENDIF
 109
 110 C...Find fastest remaining particle.
 111   210 NPRE=NPRE+1
 112       PMAX=0.
 113       DO 220 I=N+NP+1,N+2*NP
 114       IF(K(I,4).NE.0.OR.P(I,5).LE.PMAX) GOTO 220
 115       IMAX=I
 116       PMAX=P(I,5)
 117   220 CONTINUE
 118       DO 230 J=1,5
 119   230 P(N+NPRE,J)=P(IMAX,J)
 120       NREM=NREM-1
 121       K(IMAX,4)=NPRE
 122
 123 C...Sum up precluster around it according to pT separation.
 124       IF(MSTU(46).LE.2) THEN
 125         DO 250 I=N+NP+1,N+2*NP
 126         IF(K(I,4).NE.0) GOTO 250
 127         R2=R2T(I,IMAX)
 128         IF(R2.GT.RINIT**2) GOTO 250
 129         NREM=NREM-1
 130         K(I,4)=NPRE
 131         DO 240 J=1,4
 132   240   P(N+NPRE,J)=P(N+NPRE,J)+P(I,J)
 133   250   CONTINUE
 134         P(N+NPRE,5)=SQRT(P(N+NPRE,1)**2+P(N+NPRE,2)**2+P(N+NPRE,3)**2)
 135
 136 C...Sum up precluster around it according to mass separation.
 137       ELSE
 138   260   IMIN=0
 139         R2MIN=RINIT**2
 140         DO 270 I=N+NP+1,N+2*NP
 141         IF(K(I,4).NE.0) GOTO 270
 142         R2=R2M(I,N+NPRE)
 143         IF(R2.GE.R2MIN) GOTO 270
 144         IMIN=I
 145         R2MIN=R2
 146   270   CONTINUE
 147         IF(IMIN.NE.0) THEN
 148           DO 280 J=1,4
 149   280     P(N+NPRE,J)=P(N+NPRE,J)+P(IMIN,J)
 150           P(N+NPRE,5)=SQRT(P(N+NPRE,1)**2+P(N+NPRE,2)**2+P(N+NPRE,3)**2)
 151           NREM=NREM-1
 152           K(IMIN,4)=NPRE
 153           GOTO 260
 154         ENDIF
 155       ENDIF
 156
 157 C...Check if more preclusters to be found. Start over if too few.
 158       IF(RINIT.GE.0.2*PARU(43).AND.NPRE+NREM.LT.MSTU(47)) GOTO 160
 159       IF(NREM.GT.0) GOTO 210
 160       NJET=NPRE
 161
 162 C...Reassign all particles to nearest jet. Sum up new jet momenta.
 163   290 TSAV=0.
 164       PSJT=0.
 165   300 IF(MSTU(46).LE.1) THEN
 166         DO 310 I=N+1,N+NJET
 167         DO 310 J=1,4
 168   310   V(I,J)=0.
 169         DO 340 I=N+NP+1,N+2*NP
 170         R2MIN=PSS**2
 171         DO 320 IJET=N+1,N+NJET
 172         IF(P(IJET,5).LT.RINIT) GOTO 320
 173         R2=R2T(I,IJET)
 174         IF(R2.GE.R2MIN) GOTO 320
 175         IMIN=IJET
 176         R2MIN=R2
 177   320   CONTINUE
 178         K(I,4)=IMIN-N
 179         DO 330 J=1,4
 180   330   V(IMIN,J)=V(IMIN,J)+P(I,J)
 181   340   CONTINUE
 182         PSJT=0.
 183         DO 360 I=N+1,N+NJET
 184         DO 350 J=1,4
 185   350   P(I,J)=V(I,J)
 186         P(I,5)=SQRT(P(I,1)**2+P(I,2)**2+P(I,3)**2)
 187   360   PSJT=PSJT+P(I,5)
 188       ENDIF
 189
 190 C...Find two closest jets.
 191       R2MIN=2.*R2ACC
 192       DO 370 ITRY1=N+1,N+NJET-1
 193       DO 370 ITRY2=ITRY1+1,N+NJET
 194       IF(MSTU(46).LE.2) R2=R2T(ITRY1,ITRY2)
 195       IF(MSTU(46).GE.3) R2=R2M(ITRY1,ITRY2)
 196       IF(R2.GE.R2MIN) GOTO 370
 197       IMIN1=ITRY1
 198       IMIN2=ITRY2
 199       R2MIN=R2
 200   370 CONTINUE
 201
 202 C...If allowed, join two closest jets and start over.
 203       IF(NJET.GT.MSTU(47).AND.R2MIN.LT.R2ACC) THEN
 204         IREC=MIN(IMIN1,IMIN2)
 205         IDEL=MAX(IMIN1,IMIN2)
 206         DO 380 J=1,4
 207   380   P(IREC,J)=P(IMIN1,J)+P(IMIN2,J)
 208         P(IREC,5)=SQRT(P(IREC,1)**2+P(IREC,2)**2+P(IREC,3)**2)
 209         DO 390 I=IDEL+1,N+NJET
 210         DO 390 J=1,5
 211   390   P(I-1,J)=P(I,J)
 212         IF(MSTU(46).GE.2) THEN
 213           DO 400 I=N+NP+1,N+2*NP
 214           IORI=N+K(I,4)
 215           IF(IORI.EQ.IDEL) K(I,4)=IREC-N
 216   400     IF(IORI.GT.IDEL) K(I,4)=K(I,4)-1
 217         ENDIF
 218         NJET=NJET-1
 219         GOTO 290
 220
 221 C...Divide up broad jet if empty cluster in list of final ones.
 222       ELSEIF(NJET.EQ.MSTU(47).AND.MSTU(46).LE.1.AND.NLOOP.LE.2) THEN
 223         DO 410 I=N+1,N+NJET
 224   410   K(I,5)=0
 225         DO 420 I=N+NP+1,N+2*NP
 226   420   K(N+K(I,4),5)=K(N+K(I,4),5)+1
 227         IEMP=0
 228         DO 430 I=N+1,N+NJET
 229   430   IF(K(I,5).EQ.0) IEMP=I
 230         IF(IEMP.NE.0) THEN
 231           NLOOP=NLOOP+1
 232           ISPL=0
 233           R2MAX=0.
 234           DO 440 I=N+NP+1,N+2*NP
 235           IF(K(N+K(I,4),5).LE.1.OR.P(I,5).LT.RINIT) GOTO 440
 236           IJET=N+K(I,4)
 237           R2=R2T(I,IJET)
 238           IF(R2.LE.R2MAX) GOTO 440
 239           ISPL=I
 240           R2MAX=R2
 241   440     CONTINUE
 242           IF(ISPL.NE.0) THEN
 243             IJET=N+K(ISPL,4)
 244             DO 450 J=1,4
 245             P(IEMP,J)=P(ISPL,J)
 246   450       P(IJET,J)=P(IJET,J)-P(ISPL,J)
 247             P(IEMP,5)=P(ISPL,5)
 248             P(IJET,5)=SQRT(P(IJET,1)**2+P(IJET,2)**2+P(IJET,3)**2)
 249             IF(NLOOP.LE.2) GOTO 290
 250           ENDIF
 251         ENDIF
 252       ENDIF
 253
 254 C...If generalized thrust has not yet converged, continue iteration.
 255       IF(MSTU(46).LE.1.AND.NLOOP.LE.2.AND.PSJT/PSS.GT.TSAV+PARU(48))
 256      &THEN
 257         TSAV=PSJT/PSS
 258         GOTO 300
 259       ENDIF
 260
 261 C...Reorder jets according to energy.
 262       DO 460 I=N+1,N+NJET
 263       DO 460 J=1,5
 264   460 V(I,J)=P(I,J)
 265       DO 490 INEW=N+1,N+NJET
 266       PEMAX=0.
 267       DO 470 ITRY=N+1,N+NJET
 268       IF(V(ITRY,4).LE.PEMAX) GOTO 470
 269       IMAX=ITRY
 270       PEMAX=V(ITRY,4)
 271   470 CONTINUE
 272       K(INEW,1)=31
 273       K(INEW,2)=97
 274       K(INEW,3)=INEW-N
 275       K(INEW,4)=0
 276       DO 480 J=1,5
 277   480 P(INEW,J)=V(IMAX,J)
 278       V(IMAX,4)=-1.
 279   490 K(IMAX,5)=INEW
 280
 281 C...Clean up particle-jet assignments and jet information.
 282       DO 500 I=N+NP+1,N+2*NP
 283       IORI=K(N+K(I,4),5)
 284       K(I,4)=IORI-N
 285       IF(K(K(I,3),1).NE.3) K(K(I,3),4)=IORI-N
 286       K(IORI,4)=K(IORI,4)+1
 287   500 CONTINUE
 288       IEMP=0
 289       PSJT=0.
 290       DO 520 I=N+1,N+NJET
 291       K(I,5)=0
 292       PSJT=PSJT+P(I,5)
 293       P(I,5)=SQRT(MAX(P(I,4)**2-P(I,5)**2,0.))
 294       DO 510 J=1,5
 295   510 V(I,J)=0.
 296   520 IF(K(I,4).EQ.0) IEMP=I
 297
 298 C...Select storing option. Output variables. Check for failure.
 299       MSTU(61)=N+1
 300       MSTU(62)=NP
 301       MSTU(63)=NPRE
 302       PARU(61)=PS(5)
 303       PARU(62)=PSJT/PSS
 304       PARU(63)=SQRT(R2MIN)
 305       IF(NJET.LE.1) PARU(63)=0.
 306       IF(IEMP.NE.0) THEN
 307          CALL LUERRM_HIJING(8
 308      $        ,'(LUCLUS_HIJING:) failed to reconstruct as requested')
 309         NJET=-1
 310       ENDIF
 311       IF(MSTU(43).LE.1) MSTU(3)=NJET
 312       IF(MSTU(43).GE.2) N=N+NJET
 313       NSAV=NJET
 314
 315       RETURN
 316       END