PHOS/shaker/lucell.f

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