PYTHIA/jetset/luclus.F

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