hijing/src/luclus.f

0001
0002 C*********************************************************************
0003
0004       SUBROUTINE LUCLUS(NJET)
0005
0006 C...Purpose: to subdivide the particle content of an event into
0007 C...jets/clusters.
0008       COMMON/LUJETS/N,K(9000,5),P(9000,5),V(9000,5)
0009       SAVE /LUJETS/
0010       COMMON/LUDAT1/MSTU(200),PARU(200),MSTJ(200),PARJ(200)
0011       SAVE /LUDAT1/
0012       COMMON/LUDAT2/KCHG(500,3),PMAS(500,4),PARF(2000),VCKM(4,4)
0013       SAVE /LUDAT2/
0014       DIMENSION PS(5)
0015       SAVE NSAV,NP,PS,PSS,RINIT,NPRE,NREM
0016
0017 C...Functions: distance measure in pT or (pseudo)mass.
0018       R2T(I1,I2)=(P(I1,5)*P(I2,5)-P(I1,1)*P(I2,1)-P(I1,2)*P(I2,2)-
0019      &P(I1,3)*P(I2,3))*2.*P(I1,5)*P(I2,5)/(0.0001+P(I1,5)+P(I2,5))**2
0020       R2M(I1,I2)=2.*P(I1,4)*P(I2,4)*(1.-(P(I1,1)*P(I2,1)+P(I1,2)*
0021      &P(I2,2)+P(I1,3)*P(I2,3))/(P(I1,5)*P(I2,5)))
0022
0023 C...If first time, reset. If reentering, skip preliminaries.
0024       IF(MSTU(48).LE.0) THEN
0025         NP=0
0026         DO 100 J=1,5
0027   100   PS(J)=0.
0028         PSS=0.
0029       ELSE
0030         NJET=NSAV
0031         IF(MSTU(43).GE.2) N=N-NJET
0032         DO 110 I=N+1,N+NJET
0033   110   P(I,5)=SQRT(P(I,1)**2+P(I,2)**2+P(I,3)**2)
0034         IF(MSTU(46).LE.3) R2ACC=PARU(44)**2
0035         IF(MSTU(46).GE.4) R2ACC=PARU(45)*PS(5)**2
0036         NLOOP=0
0037         GOTO 290
0038       ENDIF
0039
0040 C...Find which particles are to be considered in cluster search.
0041       DO 140 I=1,N
0042       IF(K(I,1).LE.0.OR.K(I,1).GT.10) GOTO 140
0043       IF(MSTU(41).GE.2) THEN
0044         KC=LUCOMP(K(I,2))
0045         IF(KC.EQ.0.OR.KC.EQ.12.OR.KC.EQ.14.OR.KC.EQ.16.OR.
0046      &  KC.EQ.18) GOTO 140
0047         IF(MSTU(41).GE.3.AND.KCHG(KC,2).EQ.0.AND.LUCHGE(K(I,2)).EQ.0)
0048      &  GOTO 140
0049       ENDIF
0050       IF(N+2*NP.GE.MSTU(4)-MSTU(32)-5) THEN
0051         CALL LUERRM(11,'(LUCLUS:) no more memory left in LUJETS')
0052         NJET=-1
0053         RETURN
0054       ENDIF
0055
0056 C...Take copy of these particles, with space left for jets later on.
0057       NP=NP+1
0058       K(N+NP,3)=I
0059       DO 120 J=1,5
0060   120 P(N+NP,J)=P(I,J)
0061       IF(MSTU(42).EQ.0) P(N+NP,5)=0.
0062       IF(MSTU(42).EQ.1.AND.K(I,2).NE.22) P(N+NP,5)=PMAS(101,1)
0063       P(N+NP,4)=SQRT(P(N+NP,5)**2+P(I,1)**2+P(I,2)**2+P(I,3)**2)
0064       P(N+NP,5)=SQRT(P(I,1)**2+P(I,2)**2+P(I,3)**2)
0065       DO 130 J=1,4
0066   130 PS(J)=PS(J)+P(N+NP,J)
0067       PSS=PSS+P(N+NP,5)
0068   140 CONTINUE
0069       DO 150 I=N+1,N+NP
0070       K(I+NP,3)=K(I,3)
0071       DO 150 J=1,5
0072   150 P(I+NP,J)=P(I,J)
0073       PS(5)=SQRT(MAX(0.,PS(4)**2-PS(1)**2-PS(2)**2-PS(3)**2))
0074
0075 C...Very low multiplicities not considered.
0076       IF(NP.LT.MSTU(47)) THEN
0077         CALL LUERRM(8,'(LUCLUS:) too few particles for analysis')
0078         NJET=-1
0079         RETURN
0080       ENDIF
0081
0082 C...Find precluster configuration. If too few jets, make harder cuts.
0083       NLOOP=0
0084       IF(MSTU(46).LE.3) R2ACC=PARU(44)**2
0085       IF(MSTU(46).GE.4) R2ACC=PARU(45)*PS(5)**2
0086       RINIT=1.25*PARU(43)
0087       IF(NP.LE.MSTU(47)+2) RINIT=0.
0088   160 RINIT=0.8*RINIT
0089       NPRE=0
0090       NREM=NP
0091       DO 170 I=N+NP+1,N+2*NP
0092   170 K(I,4)=0
0093
0094 C...Sum up small momentum region. Jet if enough absolute momentum.
0095       IF(MSTU(46).LE.2) THEN
0096         DO 180 J=1,4
0097   180   P(N+1,J)=0.
0098         DO 200 I=N+NP+1,N+2*NP
0099         IF(P(I,5).GT.2.*RINIT) GOTO 200
0100         NREM=NREM-1
0101         K(I,4)=1
0102         DO 190 J=1,4
0103   190   P(N+1,J)=P(N+1,J)+P(I,J)
0104   200   CONTINUE
0105         P(N+1,5)=SQRT(P(N+1,1)**2+P(N+1,2)**2+P(N+1,3)**2)
0106         IF(P(N+1,5).GT.2.*RINIT) NPRE=1
0107         IF(RINIT.GE.0.2*PARU(43).AND.NPRE+NREM.LT.MSTU(47)) GOTO 160
0108       ENDIF
0109
0110 C...Find fastest remaining particle.
0111   210 NPRE=NPRE+1
0112       PMAX=0.
0113       DO 220 I=N+NP+1,N+2*NP
0114       IF(K(I,4).NE.0.OR.P(I,5).LE.PMAX) GOTO 220
0115       IMAX=I
0116       PMAX=P(I,5)
0117   220 CONTINUE
0118       DO 230 J=1,5
0119   230 P(N+NPRE,J)=P(IMAX,J)
0120       NREM=NREM-1
0121       K(IMAX,4)=NPRE
0122
0123 C...Sum up precluster around it according to pT separation.
0124       IF(MSTU(46).LE.2) THEN
0125         DO 250 I=N+NP+1,N+2*NP
0126         IF(K(I,4).NE.0) GOTO 250
0127         R2=R2T(I,IMAX)
0128         IF(R2.GT.RINIT**2) GOTO 250
0129         NREM=NREM-1
0130         K(I,4)=NPRE
0131         DO 240 J=1,4
0132   240   P(N+NPRE,J)=P(N+NPRE,J)+P(I,J)
0133   250   CONTINUE
0134         P(N+NPRE,5)=SQRT(P(N+NPRE,1)**2+P(N+NPRE,2)**2+P(N+NPRE,3)**2)
0135
0136 C...Sum up precluster around it according to mass separation.
0137       ELSE
0138   260   IMIN=0
0139         R2MIN=RINIT**2
0140         DO 270 I=N+NP+1,N+2*NP
0141         IF(K(I,4).NE.0) GOTO 270
0142         R2=R2M(I,N+NPRE)
0143         IF(R2.GE.R2MIN) GOTO 270
0144         IMIN=I
0145         R2MIN=R2
0146   270   CONTINUE
0147         IF(IMIN.NE.0) THEN
0148           DO 280 J=1,4
0149   280     P(N+NPRE,J)=P(N+NPRE,J)+P(IMIN,J)
0150           P(N+NPRE,5)=SQRT(P(N+NPRE,1)**2+P(N+NPRE,2)**2+P(N+NPRE,3)**2)
0151           NREM=NREM-1
0152           K(IMIN,4)=NPRE
0153           GOTO 260
0154         ENDIF
0155       ENDIF
0156
0157 C...Check if more preclusters to be found. Start over if too few.
0158       IF(RINIT.GE.0.2*PARU(43).AND.NPRE+NREM.LT.MSTU(47)) GOTO 160
0159       IF(NREM.GT.0) GOTO 210
0160       NJET=NPRE
0161
0162 C...Reassign all particles to nearest jet. Sum up new jet momenta.
0163   290 TSAV=0.
0164       PSJT=0.
0165   300 IF(MSTU(46).LE.1) THEN
0166         DO 310 I=N+1,N+NJET
0167         DO 310 J=1,4
0168   310   V(I,J)=0.
0169         DO 340 I=N+NP+1,N+2*NP
0170         R2MIN=PSS**2
0171         DO 320 IJET=N+1,N+NJET
0172         IF(P(IJET,5).LT.RINIT) GOTO 320
0173         R2=R2T(I,IJET)
0174         IF(R2.GE.R2MIN) GOTO 320
0175         IMIN=IJET
0176         R2MIN=R2
0177   320   CONTINUE
0178         K(I,4)=IMIN-N
0179         DO 330 J=1,4
0180   330   V(IMIN,J)=V(IMIN,J)+P(I,J)
0181   340   CONTINUE
0182         PSJT=0.
0183         DO 360 I=N+1,N+NJET
0184         DO 350 J=1,4
0185   350   P(I,J)=V(I,J)
0186         P(I,5)=SQRT(P(I,1)**2+P(I,2)**2+P(I,3)**2)
0187   360   PSJT=PSJT+P(I,5)
0188       ENDIF
0189
0190 C...Find two closest jets.
0191       R2MIN=2.*R2ACC
0192       DO 370 ITRY1=N+1,N+NJET-1
0193       DO 370 ITRY2=ITRY1+1,N+NJET
0194       IF(MSTU(46).LE.2) R2=R2T(ITRY1,ITRY2)
0195       IF(MSTU(46).GE.3) R2=R2M(ITRY1,ITRY2)
0196       IF(R2.GE.R2MIN) GOTO 370
0197       IMIN1=ITRY1
0198       IMIN2=ITRY2
0199       R2MIN=R2
0200   370 CONTINUE
0201
0202 C...If allowed, join two closest jets and start over.
0203       IF(NJET.GT.MSTU(47).AND.R2MIN.LT.R2ACC) THEN
0204         IREC=MIN(IMIN1,IMIN2)
0205         IDEL=MAX(IMIN1,IMIN2)
0206         DO 380 J=1,4
0207   380   P(IREC,J)=P(IMIN1,J)+P(IMIN2,J)
0208         P(IREC,5)=SQRT(P(IREC,1)**2+P(IREC,2)**2+P(IREC,3)**2)
0209         DO 390 I=IDEL+1,N+NJET
0210         DO 390 J=1,5
0211   390   P(I-1,J)=P(I,J)
0212         IF(MSTU(46).GE.2) THEN
0213           DO 400 I=N+NP+1,N+2*NP
0214           IORI=N+K(I,4)
0215           IF(IORI.EQ.IDEL) K(I,4)=IREC-N
0216   400     IF(IORI.GT.IDEL) K(I,4)=K(I,4)-1
0217         ENDIF
0218         NJET=NJET-1
0219         GOTO 290
0220
0221 C...Divide up broad jet if empty cluster in list of final ones.
0222       ELSEIF(NJET.EQ.MSTU(47).AND.MSTU(46).LE.1.AND.NLOOP.LE.2) THEN
0223         DO 410 I=N+1,N+NJET
0224   410   K(I,5)=0
0225         DO 420 I=N+NP+1,N+2*NP
0226   420   K(N+K(I,4),5)=K(N+K(I,4),5)+1
0227         IEMP=0
0228         DO 430 I=N+1,N+NJET
0229   430   IF(K(I,5).EQ.0) IEMP=I
0230         IF(IEMP.NE.0) THEN
0231           NLOOP=NLOOP+1
0232           ISPL=0
0233           R2MAX=0.
0234           DO 440 I=N+NP+1,N+2*NP
0235           IF(K(N+K(I,4),5).LE.1.OR.P(I,5).LT.RINIT) GOTO 440
0236           IJET=N+K(I,4)
0237           R2=R2T(I,IJET)
0238           IF(R2.LE.R2MAX) GOTO 440
0239           ISPL=I
0240           R2MAX=R2
0241   440     CONTINUE
0242           IF(ISPL.NE.0) THEN
0243             IJET=N+K(ISPL,4)
0244             DO 450 J=1,4
0245             P(IEMP,J)=P(ISPL,J)
0246   450       P(IJET,J)=P(IJET,J)-P(ISPL,J)
0247             P(IEMP,5)=P(ISPL,5)
0248             P(IJET,5)=SQRT(P(IJET,1)**2+P(IJET,2)**2+P(IJET,3)**2)
0249             IF(NLOOP.LE.2) GOTO 290
0250           ENDIF
0251         ENDIF
0252       ENDIF
0253
0254 C...If generalized thrust has not yet converged, continue iteration.
0255       IF(MSTU(46).LE.1.AND.NLOOP.LE.2.AND.PSJT/PSS.GT.TSAV+PARU(48))
0256      &THEN
0257         TSAV=PSJT/PSS
0258         GOTO 300
0259       ENDIF
0260
0261 C...Reorder jets according to energy.
0262       DO 460 I=N+1,N+NJET
0263       DO 460 J=1,5
0264   460 V(I,J)=P(I,J)
0265       DO 490 INEW=N+1,N+NJET
0266       PEMAX=0.
0267       DO 470 ITRY=N+1,N+NJET
0268       IF(V(ITRY,4).LE.PEMAX) GOTO 470
0269       IMAX=ITRY
0270       PEMAX=V(ITRY,4)
0271   470 CONTINUE
0272       K(INEW,1)=31
0273       K(INEW,2)=97
0274       K(INEW,3)=INEW-N
0275       K(INEW,4)=0
0276       DO 480 J=1,5
0277   480 P(INEW,J)=V(IMAX,J)
0278       V(IMAX,4)=-1.
0279   490 K(IMAX,5)=INEW
0280
0281 C...Clean up particle-jet assignments and jet information.
0282       DO 500 I=N+NP+1,N+2*NP
0283       IORI=K(N+K(I,4),5)
0284       K(I,4)=IORI-N
0285       IF(K(K(I,3),1).NE.3) K(K(I,3),4)=IORI-N
0286       K(IORI,4)=K(IORI,4)+1
0287   500 CONTINUE
0288       IEMP=0
0289       PSJT=0.
0290       DO 520 I=N+1,N+NJET
0291       K(I,5)=0
0292       PSJT=PSJT+P(I,5)
0293       P(I,5)=SQRT(MAX(P(I,4)**2-P(I,5)**2,0.))
0294       DO 510 J=1,5
0295   510 V(I,J)=0.
0296   520 IF(K(I,4).EQ.0) IEMP=I
0297
0298 C...Select storing option. Output variables. Check for failure.
0299       MSTU(61)=N+1
0300       MSTU(62)=NP
0301       MSTU(63)=NPRE
0302       PARU(61)=PS(5)
0303       PARU(62)=PSJT/PSS
0304       PARU(63)=SQRT(R2MIN)
0305       IF(NJET.LE.1) PARU(63)=0.
0306       IF(IEMP.NE.0) THEN
0307         CALL LUERRM(8,'(LUCLUS:) failed to reconstruct as requested')
0308         NJET=-1
0309       ENDIF
0310       IF(MSTU(43).LE.1) MSTU(3)=NJET
0311       IF(MSTU(43).GE.2) N=N+NJET
0312       NSAV=NJET
0313
0314       RETURN
0315       END