PYTHIA/jetset/lucell.F

   1
   2 C*********************************************************************
   3
   4       SUBROUTINE LUCELL(NJET)
   5
   6 C...Purpose: to provide a simple way of jet finding in an eta-phi-ET
   7 C...coordinate frame, as used for calorimeters at hadron colliders.
   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
  13 C...Loop over all particles. Find cell that was hit by given particle.
  14       PTLRAT=1./SINH(PARU(51))**2
  15       NP=0
  16       NC=N
  17       DO 110 I=1,N
  18       IF(K(I,1).LE.0.OR.K(I,1).GT.10) GOTO 110
  19       IF(P(I,1)**2+P(I,2)**2.LE.PTLRAT*P(I,3)**2) GOTO 110
  20       IF(MSTU(41).GE.2) THEN
  21         KC=LUCOMP(K(I,2))
  22         IF(KC.EQ.0.OR.KC.EQ.12.OR.KC.EQ.14.OR.KC.EQ.16.OR.
  23      &  KC.EQ.18) GOTO 110
  24         IF(MSTU(41).GE.3.AND.KCHG(KC,2).EQ.0.AND.LUCHGE(K(I,2)).EQ.0)
  25      &  GOTO 110
  26       ENDIF
  27       NP=NP+1
  28       PT=SQRT(P(I,1)**2+P(I,2)**2)
  29       ETA=SIGN(LOG((SQRT(PT**2+P(I,3)**2)+ABS(P(I,3)))/PT),P(I,3))
  30       IETA=MAX(1,MIN(MSTU(51),1+INT(MSTU(51)*0.5*(ETA/PARU(51)+1.))))
  31       PHI=ULANGL(P(I,1),P(I,2))
  32       IPHI=MAX(1,MIN(MSTU(52),1+INT(MSTU(52)*0.5*(PHI/PARU(1)+1.))))
  33       IETPH=MSTU(52)*IETA+IPHI
  34
  35 C...Add to cell already hit, or book new cell.
  36       DO 100 IC=N+1,NC
  37       IF(IETPH.EQ.K(IC,3)) THEN
  38         K(IC,4)=K(IC,4)+1
  39         P(IC,5)=P(IC,5)+PT
  40         GOTO 110
  41       ENDIF
  42   100 CONTINUE
  43       IF(NC.GE.MSTU(4)-MSTU(32)-5) THEN
  44         CALL LUERRM(11,'(LUCELL:) no more memory left in LUJETS')
  45         NJET=-2
  46         RETURN
  47       ENDIF
  48       NC=NC+1
  49       K(NC,3)=IETPH
  50       K(NC,4)=1
  51       K(NC,5)=2
  52       P(NC,1)=(PARU(51)/MSTU(51))*(2*IETA-1-MSTU(51))
  53       P(NC,2)=(PARU(1)/MSTU(52))*(2*IPHI-1-MSTU(52))
  54       P(NC,5)=PT
  55   110 CONTINUE
  56
  57 C...Smear true bin content by calorimeter resolution.
  58       IF(MSTU(53).GE.1) THEN
  59         DO 130 IC=N+1,NC
  60         PEI=P(IC,5)
  61         IF(MSTU(53).EQ.2) PEI=P(IC,5)*COSH(P(IC,1))
  62   120   PEF=PEI+PARU(55)*SQRT(-2.*LOG(MAX(1E-10,RLU(0)))*PEI)*
  63      &  COS(PARU(2)*RLU(0))
  64         IF(PEF.LT.0..OR.PEF.GT.PARU(56)*PEI) GOTO 120
  65         P(IC,5)=PEF
  66         IF(MSTU(53).EQ.2) P(IC,5)=PEF/COSH(P(IC,1))
  67   130   CONTINUE
  68       ENDIF
  69
  70 C...Remove cells below threshold.
  71       IF(PARU(58).GT.0.) THEN
  72         NCC=NC
  73         NC=N
  74         DO 140 IC=N+1,NCC
  75         IF(P(IC,5).GT.PARU(58)) THEN
  76           NC=NC+1
  77           K(NC,3)=K(IC,3)
  78           K(NC,4)=K(IC,4)
  79           K(NC,5)=K(IC,5)
  80           P(NC,1)=P(IC,1)
  81           P(NC,2)=P(IC,2)
  82           P(NC,5)=P(IC,5)
  83         ENDIF
  84   140   CONTINUE
  85       ENDIF
  86
  87 C...Find initiator cell: the one with highest pT of not yet used ones.
  88       NJ=NC
  89   150 ETMAX=0.
  90       DO 160 IC=N+1,NC
  91       IF(K(IC,5).NE.2) GOTO 160
  92       IF(P(IC,5).LE.ETMAX) GOTO 160
  93       ICMAX=IC
  94       ETA=P(IC,1)
  95       PHI=P(IC,2)
  96       ETMAX=P(IC,5)
  97   160 CONTINUE
  98       IF(ETMAX.LT.PARU(52)) GOTO 220
  99       IF(NJ.GE.MSTU(4)-MSTU(32)-5) THEN
 100         CALL LUERRM(11,'(LUCELL:) no more memory left in LUJETS')
 101         NJET=-2
 102         RETURN
 103       ENDIF
 104       K(ICMAX,5)=1
 105       NJ=NJ+1
 106       K(NJ,4)=0
 107       K(NJ,5)=1
 108       P(NJ,1)=ETA
 109       P(NJ,2)=PHI
 110       P(NJ,3)=0.
 111       P(NJ,4)=0.
 112       P(NJ,5)=0.
 113
 114 C...Sum up unused cells within required distance of initiator.
 115       DO 170 IC=N+1,NC
 116       IF(K(IC,5).EQ.0) GOTO 170
 117       IF(ABS(P(IC,1)-ETA).GT.PARU(54)) GOTO 170
 118       DPHIA=ABS(P(IC,2)-PHI)
 119       IF(DPHIA.GT.PARU(54).AND.DPHIA.LT.PARU(2)-PARU(54)) GOTO 170
 120       PHIC=P(IC,2)
 121       IF(DPHIA.GT.PARU(1)) PHIC=PHIC+SIGN(PARU(2),PHI)
 122       IF((P(IC,1)-ETA)**2+(PHIC-PHI)**2.GT.PARU(54)**2) GOTO 170
 123       K(IC,5)=-K(IC,5)
 124       K(NJ,4)=K(NJ,4)+K(IC,4)
 125       P(NJ,3)=P(NJ,3)+P(IC,5)*P(IC,1)
 126       P(NJ,4)=P(NJ,4)+P(IC,5)*PHIC
 127       P(NJ,5)=P(NJ,5)+P(IC,5)
 128   170 CONTINUE
 129
 130 C...Reject cluster below minimum ET, else accept.
 131       IF(P(NJ,5).LT.PARU(53)) THEN
 132         NJ=NJ-1
 133         DO 180 IC=N+1,NC
 134         IF(K(IC,5).LT.0) K(IC,5)=-K(IC,5)
 135   180   CONTINUE
 136       ELSEIF(MSTU(54).LE.2) THEN
 137         P(NJ,3)=P(NJ,3)/P(NJ,5)
 138         P(NJ,4)=P(NJ,4)/P(NJ,5)
 139         IF(ABS(P(NJ,4)).GT.PARU(1)) P(NJ,4)=P(NJ,4)-SIGN(PARU(2),
 140      &  P(NJ,4))
 141         DO 190 IC=N+1,NC
 142         IF(K(IC,5).LT.0) K(IC,5)=0
 143   190   CONTINUE
 144       ELSE
 145         DO 200 J=1,4
 146         P(NJ,J)=0.
 147   200   CONTINUE
 148         DO 210 IC=N+1,NC
 149         IF(K(IC,5).GE.0) GOTO 210
 150         P(NJ,1)=P(NJ,1)+P(IC,5)*COS(P(IC,2))
 151         P(NJ,2)=P(NJ,2)+P(IC,5)*SIN(P(IC,2))
 152         P(NJ,3)=P(NJ,3)+P(IC,5)*SINH(P(IC,1))
 153         P(NJ,4)=P(NJ,4)+P(IC,5)*COSH(P(IC,1))
 154         K(IC,5)=0
 155   210   CONTINUE
 156       ENDIF
 157       GOTO 150
 158
 159 C...Arrange clusters in falling ET sequence.
 160   220 DO 250 I=1,NJ-NC
 161       ETMAX=0.
 162       DO 230 IJ=NC+1,NJ
 163       IF(K(IJ,5).EQ.0) GOTO 230
 164       IF(P(IJ,5).LT.ETMAX) GOTO 230
 165       IJMAX=IJ
 166       ETMAX=P(IJ,5)
 167   230 CONTINUE
 168       K(IJMAX,5)=0
 169       K(N+I,1)=31
 170       K(N+I,2)=98
 171       K(N+I,3)=I
 172       K(N+I,4)=K(IJMAX,4)
 173       K(N+I,5)=0
 174       DO 240 J=1,5
 175       P(N+I,J)=P(IJMAX,J)
 176       V(N+I,J)=0.
 177   240 CONTINUE
 178   250 CONTINUE
 179       NJET=NJ-NC
 180
 181 C...Convert to massless or massive four-vectors.
 182       IF(MSTU(54).EQ.2) THEN
 183         DO 260 I=N+1,N+NJET
 184         ETA=P(I,3)
 185         P(I,1)=P(I,5)*COS(P(I,4))
 186         P(I,2)=P(I,5)*SIN(P(I,4))
 187         P(I,3)=P(I,5)*SINH(ETA)
 188         P(I,4)=P(I,5)*COSH(ETA)
 189         P(I,5)=0.
 190   260   CONTINUE
 191       ELSEIF(MSTU(54).GE.3) THEN
 192         DO 270 I=N+1,N+NJET
 193         P(I,5)=SQRT(MAX(0.,P(I,4)**2-P(I,1)**2-P(I,2)**2-P(I,3)**2))
 194   270   CONTINUE
 195       ENDIF
 196
 197 C...Information about storage.
 198       MSTU(61)=N+1
 199       MSTU(62)=NP
 200       MSTU(63)=NC-N
 201       IF(MSTU(43).LE.1) MSTU(3)=NJET
 202       IF(MSTU(43).GE.2) N=N+NJET
 203
 204       RETURN
 205       END