hijing/src/lucell.f

0001
0002 C*********************************************************************
0003
0004       SUBROUTINE LUCELL(NJET)
0005
0006 C...Purpose: to provide a simple way of jet finding in an eta-phi-ET
0007 C...coordinate frame, as used for calorimeters at hadron colliders.
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
0015 C...Loop over all particles. Find cell that was hit by given particle.
0016       NCE2=2*MSTU(51)*MSTU(52)
0017       PTLRAT=1./SINH(PARU(51))**2
0018       NP=0
0019       NC=N
0020       DO 110 I=1,N
0021       IF(K(I,1).LE.0.OR.K(I,1).GT.10) GOTO 110
0022       IF(P(I,1)**2+P(I,2)**2.LE.PTLRAT*P(I,3)**2) GOTO 110
0023       IF(MSTU(41).GE.2) THEN
0024         KC=LUCOMP(K(I,2))
0025         IF(KC.EQ.0.OR.KC.EQ.12.OR.KC.EQ.14.OR.KC.EQ.16.OR.
0026      &  KC.EQ.18) GOTO 110
0027         IF(MSTU(41).GE.3.AND.KCHG(KC,2).EQ.0.AND.LUCHGE(K(I,2)).EQ.0)
0028      &  GOTO 110
0029       ENDIF
0030       NP=NP+1
0031       PT=SQRT(P(I,1)**2+P(I,2)**2)
0032       ETA=SIGN(LOG((SQRT(PT**2+P(I,3)**2)+ABS(P(I,3)))/PT),P(I,3))
0033       IETA=MAX(1,MIN(MSTU(51),1+INT(MSTU(51)*0.5*(ETA/PARU(51)+1.))))
0034       PHI=ULANGL(P(I,1),P(I,2))
0035       IPHI=MAX(1,MIN(MSTU(52),1+INT(MSTU(52)*0.5*(PHI/PARU(1)+1.))))
0036       IETPH=MSTU(52)*IETA+IPHI
0037
0038 C...Add to cell already hit, or book new cell.
0039       DO 100 IC=N+1,NC
0040       IF(IETPH.EQ.K(IC,3)) THEN
0041         K(IC,4)=K(IC,4)+1
0042         P(IC,5)=P(IC,5)+PT
0043         GOTO 110
0044       ENDIF
0045   100 CONTINUE
0046       IF(NC.GE.MSTU(4)-MSTU(32)-5) THEN
0047         CALL LUERRM(11,'(LUCELL:) no more memory left in LUJETS')
0048         NJET=-2
0049         RETURN
0050       ENDIF
0051       NC=NC+1
0052       K(NC,3)=IETPH
0053       K(NC,4)=1
0054       K(NC,5)=2
0055       P(NC,1)=(PARU(51)/MSTU(51))*(2*IETA-1-MSTU(51))
0056       P(NC,2)=(PARU(1)/MSTU(52))*(2*IPHI-1-MSTU(52))
0057       P(NC,5)=PT
0058   110 CONTINUE
0059
0060 C...Smear true bin content by calorimeter resolution.
0061       IF(MSTU(53).GE.1) THEN
0062         DO 130 IC=N+1,NC
0063         PEI=P(IC,5)
0064         IF(MSTU(53).EQ.2) PEI=P(IC,5)/COSH(P(IC,1))
0065   120   PEF=PEI+PARU(55)*SQRT(-2.*LOG(MAX(1E-10,RLU(0)))*PEI)*
0066      &  COS(PARU(2)*RLU(0))
0067         IF(PEF.LT.0..OR.PEF.GT.PARU(56)*PEI) GOTO 120
0068         P(IC,5)=PEF
0069   130   IF(MSTU(53).EQ.2) P(IC,5)=PEF*COSH(P(IC,1))
0070       ENDIF
0071
0072 C...Find initiator cell: the one with highest pT of not yet used ones.
0073       NJ=NC
0074   140 ETMAX=0.
0075       DO 150 IC=N+1,NC
0076       IF(K(IC,5).NE.2) GOTO 150
0077       IF(P(IC,5).LE.ETMAX) GOTO 150
0078       ICMAX=IC
0079       ETA=P(IC,1)
0080       PHI=P(IC,2)
0081       ETMAX=P(IC,5)
0082   150 CONTINUE
0083       IF(ETMAX.LT.PARU(52)) GOTO 210
0084       IF(NJ.GE.MSTU(4)-MSTU(32)-5) THEN
0085         CALL LUERRM(11,'(LUCELL:) no more memory left in LUJETS')
0086         NJET=-2
0087         RETURN
0088       ENDIF
0089       K(ICMAX,5)=1
0090       NJ=NJ+1
0091       K(NJ,4)=0
0092       K(NJ,5)=1
0093       P(NJ,1)=ETA
0094       P(NJ,2)=PHI
0095       P(NJ,3)=0.
0096       P(NJ,4)=0.
0097       P(NJ,5)=0.
0098
0099 C...Sum up unused cells within required distance of initiator.
0100       DO 160 IC=N+1,NC
0101       IF(K(IC,5).EQ.0) GOTO 160
0102       IF(ABS(P(IC,1)-ETA).GT.PARU(54)) GOTO 160
0103       DPHIA=ABS(P(IC,2)-PHI)
0104       IF(DPHIA.GT.PARU(54).AND.DPHIA.LT.PARU(2)-PARU(54)) GOTO 160
0105       PHIC=P(IC,2)
0106       IF(DPHIA.GT.PARU(1)) PHIC=PHIC+SIGN(PARU(2),PHI)
0107       IF((P(IC,1)-ETA)**2+(PHIC-PHI)**2.GT.PARU(54)**2) GOTO 160
0108       K(IC,5)=-K(IC,5)
0109       K(NJ,4)=K(NJ,4)+K(IC,4)
0110       P(NJ,3)=P(NJ,3)+P(IC,5)*P(IC,1)
0111       P(NJ,4)=P(NJ,4)+P(IC,5)*PHIC
0112       P(NJ,5)=P(NJ,5)+P(IC,5)
0113   160 CONTINUE
0114
0115 C...Reject cluster below minimum ET, else accept.
0116       IF(P(NJ,5).LT.PARU(53)) THEN
0117         NJ=NJ-1
0118         DO 170 IC=N+1,NC
0119   170   IF(K(IC,5).LT.0) K(IC,5)=-K(IC,5)
0120       ELSEIF(MSTU(54).LE.2) THEN
0121         P(NJ,3)=P(NJ,3)/P(NJ,5)
0122         P(NJ,4)=P(NJ,4)/P(NJ,5)
0123         IF(ABS(P(NJ,4)).GT.PARU(1)) P(NJ,4)=P(NJ,4)-SIGN(PARU(2),
0124      &  P(NJ,4))
0125         DO 180 IC=N+1,NC
0126   180   IF(K(IC,1).LT.0) K(IC,1)=0
0127       ELSE
0128         DO 190 J=1,4
0129   190   P(NJ,J)=0.
0130         DO 200 IC=N+1,NC
0131         IF(K(IC,5).GE.0) GOTO 200
0132         P(NJ,1)=P(NJ,1)+P(IC,5)*COS(P(IC,2))
0133         P(NJ,2)=P(NJ,2)+P(IC,5)*SIN(P(IC,2))
0134         P(NJ,3)=P(NJ,3)+P(IC,5)*SINH(P(IC,1))
0135         P(NJ,4)=P(NJ,4)+P(IC,5)*COSH(P(IC,1))
0136         K(IC,5)=0
0137   200   CONTINUE
0138       ENDIF
0139       GOTO 140
0140
0141 C...Arrange clusters in falling ET sequence.
0142   210 DO 230 I=1,NJ-NC
0143       ETMAX=0.
0144       DO 220 IJ=NC+1,NJ
0145       IF(K(IJ,5).EQ.0) GOTO 220
0146       IF(P(IJ,5).LT.ETMAX) GOTO 220
0147       IJMAX=IJ
0148       ETMAX=P(IJ,5)
0149   220 CONTINUE
0150       K(IJMAX,5)=0
0151       K(N+I,1)=31
0152       K(N+I,2)=98
0153       K(N+I,3)=I
0154       K(N+I,4)=K(IJMAX,4)
0155       K(N+I,5)=0
0156       DO 230 J=1,5
0157       P(N+I,J)=P(IJMAX,J)
0158   230 V(N+I,J)=0.
0159       NJET=NJ-NC
0160
0161 C...Convert to massless or massive four-vectors.
0162       IF(MSTU(54).EQ.2) THEN
0163         DO 240 I=N+1,N+NJET
0164         ETA=P(I,3)
0165         P(I,1)=P(I,5)*COS(P(I,4))
0166         P(I,2)=P(I,5)*SIN(P(I,4))
0167         P(I,3)=P(I,5)*SINH(ETA)
0168         P(I,4)=P(I,5)*COSH(ETA)
0169   240   P(I,5)=0.
0170       ELSEIF(MSTU(54).GE.3) THEN
0171         DO 250 I=N+1,N+NJET
0172   250   P(I,5)=SQRT(MAX(0.,P(I,4)**2-P(I,1)**2-P(I,2)**2-P(I,3)**2))
0173       ENDIF
0174
0175 C...Information about storage.
0176       MSTU(61)=N+1
0177       MSTU(62)=NP
0178       MSTU(63)=NC-N
0179       IF(MSTU(43).LE.1) MSTU(3)=NJET
0180       IF(MSTU(43).GE.2) N=N+NJET
0181
0182       RETURN
0183       END