ISAJET/code/ispjet.F

   1 #include "isajet/pilot.h"
   2       SUBROUTINE ISPJET(DRCUT,ETCUT,NPJ,PJPT,PJPHI,PJETA)
   3 C----------------------------------------------------------------------
   4 C-
   5 C-   Purpose and Methods : COMBINES PARTONS INTO PARTON JETS
   6 C-                         based on PJCONE
   7 C-   Inputs
   8 C-        DRCUT  - dR=sqrt(dETA**2+dPHI**2) cut around Leading Partons.
   9 C-        ETCUT  - Transverse Energy cut (minimum for defining a JET ).
  10 C-
  11 C-   Outputs :
  12 C-        NPJ       = No. of Parton Jets found.
  13 C-        PJPT(NPJ) = pt of partons
  14 C-        PJPHI(NPJ)= phi "
  15 C-        PJETA(NPJ)= eta "
  16 C-
  17 C-   created  16-APR-1996   Serban D. Protopopescu
  18 C-   Updated  16-JUN-1998   F. Paige
  19 C-    Copy of ISAZEB routine ISA_PJETS to be used by IPARTNS
  20 C-
  21 C----------------------------------------------------------------------
  22 #if defined(CERNLIB_IMPNONE)
  23       IMPLICIT NONE
  24 #endif
  25 #include "isajet/pjets.inc"
  26 #include "isajet/jetset.inc"
  27 C
  28       INTEGER NPJ
  29       REAL    DRCUT
  30       REAL    PJPHI(*), PJETA(*), PJPT(*)
  31       INTEGER NP, JP, J, JO, JOP1, JOP2, JP1, JP2, ISKP, IP
  32       REAL    X1, Y1, PHI1, PHI2,TH
  33       REAL    DETA, DPHI, DR,  ETCUT
  34       INTEGER NPMAX
  35       PARAMETER (NPMAX=50)
  36       INTEGER JIORD(NPMAX), JDORD(NPMAX), JCNN(NPMAX,NPMAX)
  37       INTEGER JSKP(NPMAX)
  38       INTEGER I, JJ, K
  39       REAL    PJIN(4,NPMAX), PINPHI(NPMAX), PINETA(NPMAX)
  40       REAL    PINPT(NPMAX),PDMPT(NPMAX)
  41       REAL    PJ(4,NPMAX)
  42       REAL    EPS
  43       DOUBLE PRECISION PI, TWOPI, HALFPI, RADIAN
  44 C
  45 C last significant (correctly rounded) decimal place on VAX:
  46 C                                                |
  47 C                                                V
  48         PARAMETER (PI=        3.1415 92653 58979 32384 6 D0)
  49         PARAMETER (TWOPI=     6.2831 85307 17958 64769 3 D0)
  50         PARAMETER (HALFPI=    1.5707 96326 79489 66192 3 D0)
  51         PARAMETER (RADIAN= 0.0174532 92519 94329 5769237 D0)
  52 C
  53       PARAMETER( EPS = 1.0E-5 )
  54 C----------------------------------------------------------------------
  55 C
  56       NP=0
  57       DO 10 I=1,NJSET
  58         IF(JDCAY(I).EQ.0.AND.IABS(JTYPE(I)).LT.10) THEN
  59           NP = NP + 1
  60           DO 11 K=1,4
  61             PJIN(K,NP)=PJSET(K,I)
  62    11     CONTINUE
  63           PINPT(NP)  = SQRT( PJIN(1,NP)**2+PJIN(2,NP)**2 )
  64           PINPHI(NP) = ATAN2 (PJIN(2,NP),PJIN(1,NP)+EPS)
  65           IF(PINPHI(NP).LT.0.)PINPHI(NP)=PINPHI(NP)+TWOPI
  66           TH = ATAN2 (PINPT(NP),PJIN(3,NP)+EPS)
  67           PINETA(NP) = -ALOG ( ABS(TAN(TH/2.)) + EPS )
  68           IF(NP.GE.NPMAX) GOTO 35
  69         ENDIF
  70    10 CONTINUE
  71    35 CONTINUE    ! jump here if more than NPMAX partons
  72 C
  73 C ... Order partons in pt
  74 C
  75       DO 100 JP = 1 , NP
  76         JIORD(JP) = JP
  77   100 PDMPT(JP)=PINPT(JP)
  78       CALL ISASRT(PDMPT(1),NP,JIORD)
  79       DO 105 JP = 1 , NP
  80   105 JDORD(JP) = JIORD(NP-JP+1)
  81 C
  82 C ... Combine partons close in r space
  83 C
  84       DO 110 J = 1 , NP
  85         JO=JDORD(J)
  86   110 JCNN(JO,1)=0
  87       ISKP=0
  88       DO 120 JP1 = 1 , NP-1
  89         JOP1=JDORD(JP1)
  90 C ... Check if parton already connected to other one
  91         IF ( JCNN(JOP1,1).EQ.-1 ) GOTO 120
  92         DO 130 JP2 = JP1+1 , NP
  93           JOP2=JDORD(JP2)
  94 C ... Check if parton already connected to other one
  95           IF ( JCNN(JOP2,1).EQ.-1 ) GOTO 130
  96           DETA = PINETA(JOP1) - PINETA(JOP2)
  97           PHI1 = PINPHI(JOP1)
  98           PHI2 = PINPHI(JOP2)
  99           X1 = COS(PHI2-PHI1)
 100           Y1 = SIN(PHI2-PHI1)
 101           IF(X1.EQ.0.0) THEN
 102             DPHI = HALFPI
 103           ELSE
 104             DPHI = ATAN2(Y1,X1)
 105           END IF
 106           DR = SQRT(DETA**2+DPHI**2)
 107 C --- Criterion for combining partons
 108           IF ( DR.LT.DRCUT ) THEN
 109             JCNN(JOP1,1)=JCNN(JOP1,1)+1
 110             JCNN(JOP2,1)=-1
 111             JCNN(JOP1,JCNN(JOP1,1)+1)=JOP2
 112             ISKP=ISKP+JCNN(JOP1,1)
 113             JSKP(ISKP)=JOP2
 114           ELSE
 115             GOTO 130
 116           ENDIF
 117   130   CONTINUE
 118   120 CONTINUE
 119 C
 120 C ... Bookkeeping for parton jets
 121 C
 122       DO 150 IP = 1 , NPJ
 123         PJPHI(IP)=0.
 124         PJETA(IP)=0.
 125         PJPT(IP) =0.
 126   150 CONTINUE
 127       NPJ=0
 128       DO 200 JP1 = 1 , NP
 129         JOP1=JDORD(JP1)
 130 C ... Already connected, single parton, or has others to connect to
 131         IF ( JCNN(JOP1,1).GE.0 ) THEN
 132           NPJ=NPJ+1
 133           DO 151 K=1,4
 134             PJ(K,NPJ)=PJIN(K,JOP1)
 135   151     CONTINUE
 136           PJPHI(NPJ) = PINPHI(JOP1)
 137           PJETA(NPJ) = PINETA(JOP1)
 138           PJPT(NPJ) = PINPT(JOP1)
 139           IF ( JCNN(JOP1,1).EQ.0 ) GOTO 205
 140           DO 210 JJ = 1 , JCNN(JOP1,1)
 141             PJ(1,NPJ) = PJ(1,NPJ) + PJIN(1,JCNN(JOP1,JJ+1))
 142             PJ(2,NPJ) = PJ(2,NPJ) + PJIN(2,JCNN(JOP1,JJ+1))
 143             PJ(3,NPJ) = PJ(3,NPJ) + PJIN(3,JCNN(JOP1,JJ+1))
 144             PJ(4,NPJ) = PJ(4,NPJ) + PJIN(4,JCNN(JOP1,JJ+1))
 145   210     CONTINUE
 146           PJPT(NPJ) = SQRT( PJ(1,NPJ)**2 + PJ(2,NPJ)**2 )
 147           PJPHI(NPJ) = ATAN2 (PJ(2,NPJ),PJ(1,NPJ)+EPS)
 148           IF(PJPHI(NPJ).LT.0.)PJPHI(NPJ)=PJPHI(NPJ)+TWOPI
 149           TH = ATAN2 (PJPT(NPJ),PJ(3,NPJ)+EPS)
 150           PJETA(NPJ) = -ALOG ( ABS(TAN(TH/2.)) + EPS )
 151 C ... Criterion for dropping a parton jet ( et < etcut )
 152   205     IF ( PJPT(NPJ).GT.ETCUT ) GOTO 200
 153           NPJ=NPJ-1
 154         ENDIF
 155   200 CONTINUE
 156 C
 157 C ... Order pjets in pt
 158 C
 159       DO 300 JP = 1 , NPJ
 160         JIORD(JP) = JP
 161   300 PDMPT(JP)=PJPT(JP)
 162       CALL ISASRT(PDMPT(1),NPJ,JIORD)
 163       DO 305 JP = 1 , NPJ
 164         PINPT(JP)=PJPT(JP)
 165         PINETA(JP)=PJETA(JP)
 166         PINPHI(JP)=PJPHI(JP)
 167   305 JDORD(JP) = JIORD(NPJ-JP+1)
 168       DO 306 JP = 1 , NPJ
 169         PJPT(JP)=PINPT(JDORD(JP))
 170         PJETA(JP)=PINETA(JDORD(JP))
 171         PJPHI(JP)=PINPHI(JDORD(JP))
 172   306 CONTINUE
 173 C-
 174   999 RETURN
 175       END