pythiaeRHIC/pythia/pycell.f

0001
0002 C*********************************************************************
0003
0004 C...PYCELL
0005 C...Provides a simple way of jet finding in eta-phi-ET coordinates,
0006 C...as used for calorimeters at hadron colliders.
0007
0008       SUBROUTINE PYCELL(NJET)
0009
0010 C...Double precision and integer declarations.
0011       IMPLICIT DOUBLE PRECISION(A-H, O-Z)
0012       IMPLICIT INTEGER(I-N)
0013       INTEGER PYK,PYCHGE,PYCOMP
0014 C...Parameter statement to help give large particle numbers.
0015       PARAMETER (KSUSY1=1000000,KSUSY2=2000000,KTECHN=3000000,
0016      &KEXCIT=4000000,KDIMEN=5000000)
0017 C...Commonblocks.
0018       COMMON/PYJETS/N,NPAD,K(4000,5),P(4000,5),V(4000,5)
0019       COMMON/PYDAT1/MSTU(200),PARU(200),MSTJ(200),PARJ(200)
0020       COMMON/PYDAT2/KCHG(500,4),PMAS(500,4),PARF(2000),VCKM(4,4)
0021       SAVE /PYJETS/,/PYDAT1/,/PYDAT2/
0022
0023 C...Loop over all particles. Find cell that was hit by given particle.
0024       PTLRAT=1D0/SINH(PARU(51))**2
0025       NP=0
0026       NC=N
0027       DO 110 I=1,N
0028         IF(K(I,1).LE.0.OR.K(I,1).GT.10) GOTO 110
0029         IF(P(I,1)**2+P(I,2)**2.LE.PTLRAT*P(I,3)**2) GOTO 110
0030         IF(MSTU(41).GE.2) THEN
0031           KC=PYCOMP(K(I,2))
0032           IF(KC.EQ.0.OR.KC.EQ.12.OR.KC.EQ.14.OR.KC.EQ.16.OR.
0033      &    KC.EQ.18.OR.K(I,2).EQ.KSUSY1+22.OR.K(I,2).EQ.39.OR.
0034      &    K(I,2).EQ.KSUSY1+39) GOTO 110
0035           IF(MSTU(41).GE.3.AND.KCHG(KC,2).EQ.0.AND.PYCHGE(K(I,2)).EQ.0)
0036      &    GOTO 110
0037         ENDIF
0038         NP=NP+1
0039         PT=SQRT(P(I,1)**2+P(I,2)**2)
0040         ETA=SIGN(LOG((SQRT(PT**2+P(I,3)**2)+ABS(P(I,3)))/PT),P(I,3))
0041         IETA=MAX(1,MIN(MSTU(51),1+INT(MSTU(51)*0.5D0*
0042      &  (ETA/PARU(51)+1D0))))
0043         PHI=PYANGL(P(I,1),P(I,2))
0044         IPHI=MAX(1,MIN(MSTU(52),1+INT(MSTU(52)*0.5D0*
0045      &  (PHI/PARU(1)+1D0))))
0046         IETPH=MSTU(52)*IETA+IPHI
0047
0048 C...Add to cell already hit, or book new cell.
0049         DO 100 IC=N+1,NC
0050           IF(IETPH.EQ.K(IC,3)) THEN
0051             K(IC,4)=K(IC,4)+1
0052             P(IC,5)=P(IC,5)+PT
0053             GOTO 110
0054           ENDIF
0055   100   CONTINUE
0056         IF(NC.GE.MSTU(4)-MSTU(32)-5) THEN
0057           CALL PYERRM(11,'(PYCELL:) no more memory left in PYJETS')
0058           NJET=-2
0059           RETURN
0060         ENDIF
0061         NC=NC+1
0062         K(NC,3)=IETPH
0063         K(NC,4)=1
0064         K(NC,5)=2
0065         P(NC,1)=(PARU(51)/MSTU(51))*(2*IETA-1-MSTU(51))
0066         P(NC,2)=(PARU(1)/MSTU(52))*(2*IPHI-1-MSTU(52))
0067         P(NC,5)=PT
0068   110 CONTINUE
0069
0070 C...Smear true bin content by calorimeter resolution.
0071       IF(MSTU(53).GE.1) THEN
0072         DO 130 IC=N+1,NC
0073           PEI=P(IC,5)
0074           IF(MSTU(53).EQ.2) PEI=P(IC,5)*COSH(P(IC,1))
0075   120     PEF=PEI+PARU(55)*SQRT(-2D0*LOG(MAX(1D-10,PYR(0)))*PEI)*
0076      &    COS(PARU(2)*PYR(0))
0077           IF(PEF.LT.0D0.OR.PEF.GT.PARU(56)*PEI) GOTO 120
0078           P(IC,5)=PEF
0079           IF(MSTU(53).EQ.2) P(IC,5)=PEF/COSH(P(IC,1))
0080   130   CONTINUE
0081       ENDIF
0082
0083 C...Remove cells below threshold.
0084       IF(PARU(58).GT.0D0) THEN
0085         NCC=NC
0086         NC=N
0087         DO 140 IC=N+1,NCC
0088           IF(P(IC,5).GT.PARU(58)) THEN
0089             NC=NC+1
0090             K(NC,3)=K(IC,3)
0091             K(NC,4)=K(IC,4)
0092             K(NC,5)=K(IC,5)
0093             P(NC,1)=P(IC,1)
0094             P(NC,2)=P(IC,2)
0095             P(NC,5)=P(IC,5)
0096           ENDIF
0097   140   CONTINUE
0098       ENDIF
0099
0100 C...Find initiator cell: the one with highest pT of not yet used ones.
0101       NJ=NC
0102   150 ETMAX=0D0
0103       DO 160 IC=N+1,NC
0104         IF(K(IC,5).NE.2) GOTO 160
0105         IF(P(IC,5).LE.ETMAX) GOTO 160
0106         ICMAX=IC
0107         ETA=P(IC,1)
0108         PHI=P(IC,2)
0109         ETMAX=P(IC,5)
0110   160 CONTINUE
0111       IF(ETMAX.LT.PARU(52)) GOTO 220
0112       IF(NJ.GE.MSTU(4)-MSTU(32)-5) THEN
0113         CALL PYERRM(11,'(PYCELL:) no more memory left in PYJETS')
0114         NJET=-2
0115         RETURN
0116       ENDIF
0117       K(ICMAX,5)=1
0118       NJ=NJ+1
0119       K(NJ,4)=0
0120       K(NJ,5)=1
0121       P(NJ,1)=ETA
0122       P(NJ,2)=PHI
0123       P(NJ,3)=0D0
0124       P(NJ,4)=0D0
0125       P(NJ,5)=0D0
0126
0127 C...Sum up unused cells within required distance of initiator.
0128       DO 170 IC=N+1,NC
0129         IF(K(IC,5).EQ.0) GOTO 170
0130         IF(ABS(P(IC,1)-ETA).GT.PARU(54)) GOTO 170
0131         DPHIA=ABS(P(IC,2)-PHI)
0132         IF(DPHIA.GT.PARU(54).AND.DPHIA.LT.PARU(2)-PARU(54)) GOTO 170
0133         PHIC=P(IC,2)
0134         IF(DPHIA.GT.PARU(1)) PHIC=PHIC+SIGN(PARU(2),PHI)
0135         IF((P(IC,1)-ETA)**2+(PHIC-PHI)**2.GT.PARU(54)**2) GOTO 170
0136         K(IC,5)=-K(IC,5)
0137         K(NJ,4)=K(NJ,4)+K(IC,4)
0138         P(NJ,3)=P(NJ,3)+P(IC,5)*P(IC,1)
0139         P(NJ,4)=P(NJ,4)+P(IC,5)*PHIC
0140         P(NJ,5)=P(NJ,5)+P(IC,5)
0141   170 CONTINUE
0142
0143 C...Reject cluster below minimum ET, else accept.
0144       IF(P(NJ,5).LT.PARU(53)) THEN
0145         NJ=NJ-1
0146         DO 180 IC=N+1,NC
0147           IF(K(IC,5).LT.0) K(IC,5)=-K(IC,5)
0148   180   CONTINUE
0149       ELSEIF(MSTU(54).LE.2) THEN
0150         P(NJ,3)=P(NJ,3)/P(NJ,5)
0151         P(NJ,4)=P(NJ,4)/P(NJ,5)
0152         IF(ABS(P(NJ,4)).GT.PARU(1)) P(NJ,4)=P(NJ,4)-SIGN(PARU(2),
0153      &  P(NJ,4))
0154         DO 190 IC=N+1,NC
0155           IF(K(IC,5).LT.0) K(IC,5)=0
0156   190   CONTINUE
0157       ELSE
0158         DO 200 J=1,4
0159           P(NJ,J)=0D0
0160   200   CONTINUE
0161         DO 210 IC=N+1,NC
0162           IF(K(IC,5).GE.0) GOTO 210
0163           P(NJ,1)=P(NJ,1)+P(IC,5)*COS(P(IC,2))
0164           P(NJ,2)=P(NJ,2)+P(IC,5)*SIN(P(IC,2))
0165           P(NJ,3)=P(NJ,3)+P(IC,5)*SINH(P(IC,1))
0166           P(NJ,4)=P(NJ,4)+P(IC,5)*COSH(P(IC,1))
0167           K(IC,5)=0
0168   210   CONTINUE
0169       ENDIF
0170       GOTO 150
0171
0172 C...Arrange clusters in falling ET sequence.
0173   220 DO 250 I=1,NJ-NC
0174         ETMAX=0D0
0175         DO 230 IJ=NC+1,NJ
0176           IF(K(IJ,5).EQ.0) GOTO 230
0177           IF(P(IJ,5).LT.ETMAX) GOTO 230
0178           IJMAX=IJ
0179           ETMAX=P(IJ,5)
0180   230   CONTINUE
0181         K(IJMAX,5)=0
0182         K(N+I,1)=31
0183         K(N+I,2)=98
0184         K(N+I,3)=I
0185         K(N+I,4)=K(IJMAX,4)
0186         K(N+I,5)=0
0187         DO 240 J=1,5
0188           P(N+I,J)=P(IJMAX,J)
0189           V(N+I,J)=0D0
0190   240   CONTINUE
0191   250 CONTINUE
0192       NJET=NJ-NC
0193
0194 C...Convert to massless or massive four-vectors.
0195       IF(MSTU(54).EQ.2) THEN
0196         DO 260 I=N+1,N+NJET
0197           ETA=P(I,3)
0198           P(I,1)=P(I,5)*COS(P(I,4))
0199           P(I,2)=P(I,5)*SIN(P(I,4))
0200           P(I,3)=P(I,5)*SINH(ETA)
0201           P(I,4)=P(I,5)*COSH(ETA)
0202           P(I,5)=0D0
0203   260   CONTINUE
0204       ELSEIF(MSTU(54).GE.3) THEN
0205         DO 270 I=N+1,N+NJET
0206           P(I,5)=SQRT(MAX(0D0,P(I,4)**2-P(I,1)**2-P(I,2)**2-P(I,3)**2))
0207   270   CONTINUE
0208       ENDIF
0209
0210 C...Information about storage.
0211       MSTU(61)=N+1
0212       MSTU(62)=NP
0213       MSTU(63)=NC-N
0214       IF(MSTU(43).LE.1) MSTU(3)=MAX(0,NJET)
0215       IF(MSTU(43).GE.2) N=N+MAX(0,NJET)
0216
0217       RETURN
0218       END