pythiaeRHIC/pythia/pyclus.f

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