JETAN/pxcone.F

   1       SUBROUTINE PXCONE(MODE,NTRAK,ITKDM,PTRAK,CONER,EPSLON,OVLIM,
   2      +                   MXJET,NJET,PJET,IPASS,IJMUL,IERR)
   3 *.*********************************************************
   4 *. ------
   5 *. PXCONE
   6 *. ------
   7 *.
   8 *.********** Pre Release Version 26.2.93
   9 *.
  10 *. Driver for the Cone  Jet finding algorithm of L.A. del Pozo.
  11 *. Based on algorithm from D.E. Soper.
  12 *. Finds jets inside cone of half angle CONER with energy > EPSLON.
  13 *. Jets which receive more than a fraction OVLIM of their energy from
  14 *. overlaps with other jets are excluded.
  15 *. Output jets are ordered in energy.
  16 *. If MODE.EQ.2 momenta are stored as (eta,phi,<empty>,pt)
  17 *. Usage     :
  18 *.
  19 *.      INTEGER  ITKDM,MXTRK
  20 *.      PARAMETER  (ITKDM=4.or.more,MXTRK=1.or.more)
  21 *.      INTEGER  MXJET, MXTRAK, MXPROT
  22 *.      PARAMETER  (MXJET=10,MXTRAK=500,MXPROT=500)
  23 *.      INTEGER  IPASS (MXTRAK),IJMUL (MXJET)
  24 *.      INTEGER  NTRAK,NJET,IERR,MODE
  25 *.      DOUBLE PRECISION  PTRAK (ITKDM,MXTRK),PJET (5,MXJET)
  26 *.      DOUBLE PRECISION  CONER, EPSLON, OVLIM
  27 *.      NTRAK = 1.to.MXTRAK
  28 *.      CONER   = ...
  29 *.      EPSLON  = ...
  30 *.      OVLIM   = ...
  31 *.      CALL PXCONE (MODE,NTRAK,ITKDM,PTRAK,CONER,EPSLON,OVLIM,MXJET,
  32 *.     +             NJET,PJET,IPASS,IJMUL,IERR)
  33 *.
  34 *. INPUT     :  MODE      1=>e+e-, 2=>hadron-hadron
  35 *. INPUT     :  NTRAK     Number of particles
  36 *. INPUT     :  ITKDM     First dimension of PTRAK array
  37 *. INPUT     :  PTRAK     Array of particle 4-momenta (Px,Py,Pz,E)
  38 *. INPUT     :  CONER     Cone size (half angle) in radians
  39 *. INPUT     :  EPSLON    Minimum Jet energy (GeV)
  40 *. INPUT     :  OVLIM     Maximum fraction of overlap energy in a jet
  41 *. INPUT     :  MXJET     Maximum possible number of jets
  42 *. OUTPUT    :  NJET      Number of jets found
  43 *. OUTPUT    :  PJET      5-vectors of jets
  44 *. OUTPUT    :  IPASS(k)  Particle k belongs to jet number IPASS(k)
  45 *.                        IPASS = -1 if not assosciated to a jet
  46 *. OUTPUT    :  IJMUL(i)  Jet i contains IJMUL(i) particles
  47 *. OUTPUT    :  IERR      = 0 if all is OK ;   = -1 otherwise
  48 *.
  49 *. CALLS     : PXSEAR, PXSAME, PXNEW, PXTRY, PXORD, PXUVEC, PXOLAP
  50 *. CALLED    : User
  51 *.
  52 *. AUTHOR    :  L.A. del Pozo
  53 *. CREATED   :  26-Feb-93
  54 *. LAST MOD  :   2-Mar-93
  55 *.
  56 *. Modification Log.
  57 *. 2-Jan-97: M Wobisch    - fix bug concerning COS2R in eta phi mode
  58 *. 4-Apr-93: M H Seymour  - Change 2d arrays to 1d in PXTRY & PXNEW
  59 *. 2-Apr-93: M H Seymour  - Major changes to add boost-invariant mode
  60 *. 1-Apr-93: M H Seymour  - Increase all array sizes
  61 *. 30-Mar-93: M H Seymour - Change all REAL variables to DOUBLE PRECISION
  62 *. 30-Mar-93: M H Seymour - Change OVLIM into an input parameter
  63 *. 2-Mar-93: L A del Pozo - Fix Bugs in PXOLAP
  64 *. 1-Mar-93: L A del Pozo - Remove Cern library routine calls
  65 *. 1-Mar-93: L A del Pozo - Add Print out of welcome and R and Epsilon
  66 *.
  67 *.*********************************************************
  68 C+SEQ,DECLARE.
  69 *** External Arrays
  70       INTEGER  ITKDM,MXJET,NTRAK,NJET,IERR,MODE
  71       INTEGER  IPASS (NTRAK),IJMUL (MXJET)
  72       DOUBLE PRECISION  PTRAK (ITKDM,NTRAK),PJET (5,MXJET),
  73      +     CONER, EPSLON, OVLIM
  74 *** Internal Arrays
  75       INTEGER MXPROT, MXTRAK
  76       PARAMETER (MXPROT=5000, MXTRAK=5000)
  77       DOUBLE PRECISION PP(4,MXTRAK), PU(3,MXTRAK), PJ(4,MXPROT)
  78       LOGICAL JETLIS(MXPROT,MXTRAK)
  79 *** Used in the routine.
  80       DOUBLE PRECISION COSR,COS2R, VSEED(3), VEC1(3), VEC2(3),PTSQ,PPSQ,
  81      +     COSVAL,PXMDPI
  82 cMWobisch
  83       DOUBLE PRECISION RSEP
  84 cMWobisch
  85       LOGICAL UNSTBL
  86       INTEGER I,J,N,MU,N1,N2, ITERR
  87       INTEGER NCALL, NPRINT
  88       DOUBLE PRECISION ROLD, EPSOLD, OVOLD
  89       SAVE NCALL,NPRINT,ROLD, EPSOLD, OVOLD
  90       DATA NCALL,NPRINT /0,0/
  91       DATA ROLD,EPSOLD,OVOLD/0.,0.,0./
  92
  93 cMWobisch
  94 c***************************************
  95       RSEP  = 2D0
  96 c***************************************
  97 cMWobisch
  98       IERR=0
  99 *
 100 *** INITIALIZE
 101       IF(NCALL.LE.0)  THEN
 102          ROLD = 0.
 103          EPSOLD = 0.
 104          OVOLD = 0.
 105       ENDIF
 106       NCALL = NCALL + 1
 107 *
 108 *** Print welcome and Jetfinder parameters
 109       IF((CONER.NE.ROLD .OR. EPSLON.NE.EPSOLD .OR. OVLIM.NE.OVOLD)
 110      +     .AND. NPRINT.LE.10) THEN
 111          WRITE (6,*)
 112          WRITE (6,*) ' *********** PXCONE: Cone Jet-finder ***********'
 113          WRITE (6,*) '    Written by Luis Del Pozo of OPAL'
 114          WRITE (6,*) '    Modified for eta-phi by Mike Seymour'
 115          WRITE(6,1000)'   Cone Size R = ',CONER,' Radians'
 116          WRITE(6,1001)'   Min Jet energy Epsilon = ',EPSLON,' GeV'
 117          WRITE(6,1002)'   Overlap fraction parameter = ',OVLIM
 118          WRITE (6,*) ' ***********************************************'
 119 cMWobisch
 120          IF (RSEP .lt. 1.999) THEN
 121             WRITE(6,*) ' '
 122             WRITE (6,*) ' ******************************************'
 123             WRITE (6,*) ' ******************************************'
 124             WRITE(6,*) ' M Wobisch: private change !!!!!!!!!!!! '
 125             WRITE(6,*) '      Rsep is set to ',RSEP
 126             WRITE(6,*) ' this is ONLY meaningful in a NLO calculation'
 127             WRITE(6,*) '      ------------------------  '
 128             WRITE(6,*) '  please check what you''re doing!!'
 129             WRITE(6,*) '   or ask:  Markus.Wobisch@desy.de --'
 130             WRITE (6,*) ' ******************************************'
 131             WRITE (6,*) ' ******************************************'
 132             WRITE (6,*) ' ******************************************'
 133             WRITE(6,*) ' '
 134             WRITE(6,*) ' '
 135          ENDIF
 136 cMWobisch
 137
 138          WRITE (6,*)
 139 1000     FORMAT(A18,F5.2,A10)
 140 1001     FORMAT(A29,F5.2,A5)
 141 1002     FORMAT(A33,F5.2)
 142          NPRINT = NPRINT + 1
 143          ROLD=CONER
 144          EPSOLD=EPSLON
 145          OVOLD=OVLIM
 146       ENDIF
 147 *
 148 *** Copy calling array PTRAK  to internal array PP(4,NTRAK)
 149 *
 150       IF (NTRAK .GT. MXTRAK) THEN
 151          WRITE (6,*) ' PXCONE: Ntrak too large: ',NTRAK
 152          IERR=-1
 153          RETURN
 154       ENDIF
 155       IF (MODE.NE.2) THEN
 156          DO  100 I=1, NTRAK
 157             DO  101 J=1,4
 158                PP(J,I)=PTRAK(J,I)
 159 101         CONTINUE
 160 100      CONTINUE
 161       ELSE
 162 *** Converting to eta,phi,pt if necessary
 163          DO  104 I=1,NTRAK
 164             PTSQ=PTRAK(1,I)**2+PTRAK(2,I)**2
 165             PPSQ=(SQRT(PTSQ+PTRAK(3,I)**2)+ABS(PTRAK(3,I)))**2
 166             IF (PTSQ.LE.4.25E-18*PPSQ) THEN
 167                PP(1,I)=20
 168             ELSE
 169                PP(1,I)=0.5*LOG(PPSQ/PTSQ)
 170             ENDIF
 171             PP(1,I)=SIGN(PP(1,I),PTRAK(3,I))
 172             IF (PTSQ.EQ.0) THEN
 173                PP(2,I)=0
 174             ELSE
 175                PP(2,I)=ATAN2(PTRAK(2,I),PTRAK(1,I))
 176             ENDIF
 177             PP(3,I)=0
 178             PP(4,I)=SQRT(PTSQ)
 179             PU(1,I)=PP(1,I)
 180             PU(2,I)=PP(2,I)
 181             PU(3,I)=PP(3,I)
 182 104      CONTINUE
 183       ENDIF
 184 *
 185 *** Zero output variables
 186 *
 187       NJET=0
 188       DO 102 I = 1, NTRAK
 189          DO 103 J = 1, MXPROT
 190            JETLIS(J,I) = .FALSE.
 191 103      CONTINUE
 192 102   CONTINUE
 193       CALL PXZERV(4*MXPROT,PJ)
 194       CALL PXZERI(MXJET,IJMUL)
 195 *
 196       IF (MODE.NE.2) THEN
 197          COSR = COS(CONER)
 198          COS2R = COS(CONER)
 199       ELSE
 200 *** Purely for convenience, work in terms of 1-R**2
 201          COSR = 1-CONER**2
 202 cMW -- select Rsep: 1-(Rsep*CONER)**2
 203          COS2R =  1-(RSEP*CONER)**2
 204 cORIGINAL         COS2R =  1-(2*CONER)**2
 205       ENDIF
 206       UNSTBL = .FALSE.
 207       IF (MODE.NE.2) THEN
 208          CALL PXUVEC(NTRAK,PP,PU,IERR)
 209          IF (IERR .NE. 0) RETURN
 210       ENDIF
 211 *** Look for jets using particle diretions as seed axes
 212 *
 213       DO 110 N = 1,NTRAK
 214         DO 120 MU = 1,3
 215           VSEED(MU) = PU(MU,N)
 216 120     CONTINUE
 217         CALL PXSEAR(MODE,COSR,NTRAK,PU,PP,VSEED,
 218      &                   NJET,JETLIS,PJ,UNSTBL,IERR)
 219          IF (IERR .NE. 0) RETURN
 220 110   CONTINUE
 221
 222 cMW - for Rsep=1 goto 145
 223 c      GOTO 145
 224
 225 *** Now look between all pairs of jets as seed axes.
 226       DO 140 N1 = 1,NJET-1
 227          VEC1(1)=PJ(1,N1)
 228          VEC1(2)=PJ(2,N1)
 229          VEC1(3)=PJ(3,N1)
 230          IF (MODE.NE.2) CALL PXNORV(3,VEC1,VEC1,ITERR)
 231          DO 150 N2 = N1+1,NJET
 232             VEC2(1)=PJ(1,N2)
 233             VEC2(2)=PJ(2,N2)
 234             VEC2(3)=PJ(3,N2)
 235             IF (MODE.NE.2) CALL PXNORV(3,VEC2,VEC2,ITERR)
 236             CALL PXADDV(3,VEC1,VEC2,VSEED,ITERR)
 237             IF (MODE.NE.2) THEN
 238                CALL PXNORV(3,VSEED,VSEED,ITERR)
 239             ELSE
 240                VSEED(1)=VSEED(1)/2
 241                VSEED(2)=VSEED(2)/2
 242             ENDIF
 243 C---ONLY BOTHER IF THEY ARE BETWEEN 1 AND 2 CONE RADII APART
 244             IF (MODE.NE.2) THEN
 245               COSVAL=VEC1(1)*VEC2(1)+VEC1(2)*VEC2(2)+VEC1(3)*VEC2(3)
 246             ELSE
 247                IF (ABS(VEC1(1)).GE.20.OR.ABS(VEC2(1)).GE.20) THEN
 248                   COSVAL=-1000
 249                ELSE
 250                   COSVAL=1-
 251      +              ((VEC1(1)-VEC2(1))**2+PXMDPI(VEC1(2)-VEC2(2))**2)
 252                ENDIF
 253             ENDIF
 254
 255             IF (COSVAL.LE.COSR.AND.COSVAL.GE.COS2R)
 256      +           CALL PXSEAR(MODE,COSR,NTRAK,PU,PP,VSEED,NJET,
 257      +           JETLIS,PJ,UNSTBL,IERR)
 258 c            CALL PXSEAR(MODE,COSR,NTRAK,PU,PP,VSEED,NJET,
 259 c     +           JETLIS,PJ,UNSTBL,IERR)
 260             IF (IERR .NE. 0) RETURN
 261 150      CONTINUE
 262 140   CONTINUE
 263       IF (UNSTBL) THEN
 264         IERR=-1
 265         WRITE (6,*) ' PXCONE: Too many iterations to find a proto-jet'
 266         RETURN
 267       ENDIF
 268
 269  145  CONTINUE
 270 *** Now put the jet list into order by jet energy, eliminating jets
 271 *** with energy less than EPSLON.
 272        CALL PXORD(EPSLON,NJET,NTRAK,JETLIS,PJ)
 273 *
 274 *** Take care of jet overlaps
 275        CALL PXOLAP(MODE,NJET,NTRAK,JETLIS,PJ,PP,OVLIM)
 276 *
 277 *** Order jets again as some have been eliminated, or lost energy.
 278        CALL PXORD(EPSLON,NJET,NTRAK,JETLIS,PJ)
 279 *
 280 *** All done!, Copy output into output arrays
 281       IF (NJET .GT. MXJET) THEN
 282          WRITE (6,*) ' PXCONE:  Found more than MXJET jets'
 283          IERR=-1
 284          GOTO 99
 285       ENDIF
 286       IF (MODE.NE.2) THEN
 287          DO 300 I=1, NJET
 288             DO 310 J=1,4
 289                PJET(J,I)=PJ(J,I)
 290 310         CONTINUE
 291 300      CONTINUE
 292       ELSE
 293          DO 315 I=1, NJET
 294             PJET(1,I)=PJ(4,I)*COS(PJ(2,I))
 295             PJET(2,I)=PJ(4,I)*SIN(PJ(2,I))
 296             PJET(3,I)=PJ(4,I)*SINH(PJ(1,I))
 297             PJET(4,I)=PJ(4,I)*COSH(PJ(1,I))
 298  315     CONTINUE
 299       ENDIF
 300       DO 320 I=1, NTRAK
 301          IPASS(I)=-1
 302          DO 330 J=1, NJET
 303             IF (JETLIS(J,I)) THEN
 304                IJMUL(J)=IJMUL(J)+1
 305                IPASS(I)=J
 306             ENDIF
 307 330      CONTINUE
 308 320   CONTINUE
 309 99    RETURN
 310       END
 311 *CMZ :  1.06/00 28/02/94  15.44.44  by  P. Schleper
 312 *-- Author :
 313 C-----------------------------------------------------------------------
 314       SUBROUTINE PXNORV(N,A,B,ITERR)
 315       INTEGER I,N,ITERR
 316       DOUBLE PRECISION A(N),B(N),C
 317       C=0
 318       DO 10 I=1,N
 319         C=C+A(I)**2
 320  10   CONTINUE
 321       IF (C.LE.0) RETURN
 322       C=1/SQRT(C)
 323       DO 20 I=1,N
 324         B(I)=A(I)*C
 325  20   CONTINUE
 326       END
 327 *CMZ :  2.00/00 10/01/95  10.17.57  by  P. Schleper
 328 *CMZ :  1.06/00 15/03/94  12.17.46  by  P. Schleper
 329 *-- Author :
 330 *
 331 C+DECK,PXOLAP.
 332       SUBROUTINE PXOLAP(MODE,NJET,NTRAK,JETLIS,PJ,PP,OVLIM)
 333 *
 334 *** Looks for particles assigned to more than 1 jet, and reassigns them
 335 *** If more than a fraction OVLIM of a jet's energy is contained in
 336 *** higher energy jets, that jet is neglected.
 337 *** Particles assigned to the jet closest in angle (a la CDF, Snowmass).
 338 C+SEQ,DECLARE.
 339       INTEGER MXTRAK, MXPROT
 340       PARAMETER (MXTRAK=5000,MXPROT=5000)
 341       INTEGER NJET, NTRAK, MODE
 342       LOGICAL JETLIS(MXPROT,MXTRAK)
 343       DOUBLE PRECISION PJ(4,MXPROT),PP(4,MXTRAK),PXMDPI
 344       INTEGER I,J,N,MU
 345       LOGICAL OVELAP
 346       DOUBLE PRECISION EOVER
 347       DOUBLE PRECISION OVLIM
 348       INTEGER ITERR, IJMIN, IJET(MXPROT), NJ
 349       DOUBLE PRECISION VEC1(3), VEC2(3), COST, THET, THMIN
 350       DATA IJMIN/0/
 351 *
 352       IF (NJET.LE.1) RETURN
 353 *** Look for jets with large overlaps with higher energy jets.
 354       DO 100 I = 2,NJET
 355 *** Find overlap energy between jets I and all higher energy jets.
 356        EOVER = 0.0
 357        DO 110 N = 1,NTRAK
 358          OVELAP = .FALSE.
 359          DO 120 J= 1,I-1
 360            IF (JETLIS(I,N).AND.JETLIS(J,N)) THEN
 361             OVELAP = .TRUE.
 362            ENDIF
 363 120      CONTINUE
 364          IF (OVELAP) THEN
 365            EOVER = EOVER + PP(4,N)
 366          ENDIF
 367 110     CONTINUE
 368 *** Is the fraction of energy shared larger than OVLIM?
 369         IF (EOVER.GT.OVLIM*PJ(4,I)) THEN
 370 *** De-assign all particles from Jet I
 371             DO 130 N = 1,NTRAK
 372               JETLIS(I,N) = .FALSE.
 373 130         CONTINUE
 374          ENDIF
 375 100   CONTINUE
 376 *** Now there are no big overlaps, assign every particle in
 377 *** more than 1 jet to the closet jet.
 378 *** Any particles now in more than 1 jet are assigned to the CLOSET
 379 *** jet (in angle).
 380       DO 140 I=1,NTRAK
 381          NJ=0
 382          DO 150 J=1, NJET
 383          IF(JETLIS(J,I)) THEN
 384             NJ=NJ+1
 385             IJET(NJ)=J
 386          ENDIF
 387 150      CONTINUE
 388          IF (NJ .GT. 1) THEN
 389 *** Particle in > 1 jet - calc angles...
 390             VEC1(1)=PP(1,I)
 391             VEC1(2)=PP(2,I)
 392             VEC1(3)=PP(3,I)
 393             THMIN=0.
 394             DO 160 J=1,NJ
 395                VEC2(1)=PJ(1,IJET(J))
 396                VEC2(2)=PJ(2,IJET(J))
 397                VEC2(3)=PJ(3,IJET(J))
 398                IF (MODE.NE.2) THEN
 399                   CALL PXANG3(VEC1,VEC2,COST,THET,ITERR)
 400                ELSE
 401                   THET=(VEC1(1)-VEC2(1))**2+PXMDPI(VEC1(2)-VEC2(2))**2
 402                ENDIF
 403                IF (J .EQ. 1) THEN
 404                   THMIN=THET
 405                   IJMIN=IJET(J)
 406                ELSEIF (THET .LT. THMIN) THEN
 407                   THMIN=THET
 408                   IJMIN=IJET(J)
 409                ENDIF
 410 160         CONTINUE
 411 *** Assign track to IJMIN
 412             DO 170 J=1,NJET
 413                JETLIS(J,I) = .FALSE.
 414 170         CONTINUE
 415             JETLIS(IJMIN,I)=.TRUE.
 416          ENDIF
 417 140   CONTINUE
 418 *** Recompute PJ
 419       DO 200 I = 1,NJET
 420         DO 210 MU = 1,4
 421           PJ(MU,I) = 0.0
 422 210     CONTINUE
 423         DO 220 N = 1,NTRAK
 424           IF( JETLIS(I,N) ) THEN
 425              IF (MODE.NE.2) THEN
 426                 DO 230 MU = 1,4
 427                    PJ(MU,I) = PJ(MU,I) + PP(MU,N)
 428 230             CONTINUE
 429              ELSE
 430                 PJ(1,I)=PJ(1,I)
 431      +               + PP(4,N)/(PP(4,N)+PJ(4,I))*(PP(1,N)-PJ(1,I))
 432                 PJ(2,I)=PJ(2,I)
 433      +               + PP(4,N)/(PP(4,N)+PJ(4,I))*PXMDPI(PP(2,N)-PJ(2,I))
 434                 PJ(4,I)=PJ(4,I)+PP(4,N)
 435              ENDIF
 436           ENDIF
 437 220     CONTINUE
 438 200   CONTINUE
 439       RETURN
 440       END
 441 *CMZ :  2.00/00 10/01/95  10.17.57  by  P. Schleper
 442 *CMZ :  1.06/00 14/03/94  15.37.45  by  P. Schleper
 443 *-- Author :
 444 *
 445 C+DECK,PXORD.
 446        SUBROUTINE PXORD(EPSLON,NJET,NTRAK,JETLIS,PJ)
 447 *
 448 *** Routine to put jets into order and eliminate tose less than EPSLON
 449 C+SEQ,DECLARE.
 450       INTEGER MXTRAK,MXPROT
 451       PARAMETER (MXTRAK=5000,MXPROT=5000)
 452        INTEGER I, J, INDEX(MXPROT)
 453        DOUBLE PRECISION PTEMP(4,MXPROT), ELIST(MXPROT)
 454        INTEGER NJET,NTRAK
 455        LOGICAL JETLIS(MXPROT,MXTRAK)
 456        LOGICAL LOGTMP(MXPROT,MXTRAK)
 457        DOUBLE PRECISION EPSLON,PJ(4,MXPROT)
 458 *** Puts jets in order of energy: 1 = highest energy etc.
 459 *** Then Eliminate jets with energy below EPSLON
 460 *
 461 *** Copy input arrays.
 462       DO 100 I=1,NJET
 463          DO 110 J=1,4
 464             PTEMP(J,I)=PJ(J,I)
 465 110      CONTINUE
 466          DO 120 J=1,NTRAK
 467             LOGTMP(I,J)=JETLIS(I,J)
 468 120      CONTINUE
 469 100   CONTINUE
 470       DO 150 I=1,NJET
 471          ELIST(I)=PJ(4,I)
 472 150   CONTINUE
 473 *** Sort the energies...
 474       CALL PXSORV(NJET,ELIST,INDEX,'I')
 475 *** Fill PJ and JETLIS according to sort ( sort is in ascending order!!)
 476       DO 200 I=1, NJET
 477          DO 210 J=1,4
 478             PJ(J,I)=PTEMP(J,INDEX(NJET+1-I))
 479 210      CONTINUE
 480          DO 220 J=1,NTRAK
 481             JETLIS(I,J)=LOGTMP(INDEX(NJET+1-I),J)
 482 220      CONTINUE
 483 200   CONTINUE
 484 ** Jets are now in order
 485 *** Now eliminate jets with less than Epsilon energy
 486       DO 300, I=1, NJET
 487          IF (PJ(4,I) .LT. EPSLON) THEN
 488             NJET=NJET-1
 489             PJ(4,I)=0.
 490          ENDIF
 491 300   CONTINUE
 492       RETURN
 493       END
 494
 495 ********************************************************************
 496 *CMZ :  2.00/00 10/01/95  10.17.57  by  P. Schleper
 497 *CMZ :  1.06/00 14/03/94  15.37.44  by  P. Schleper
 498 *-- Author :
 499 C+DECK,PXSEAR.
 500       SUBROUTINE PXSEAR(MODE,COSR,NTRAK,PU,PP,VSEED,NJET,
 501      +                JETLIS,PJ,UNSTBL,IERR)
 502 *
 503 C+SEQ,DECLARE.
 504       INTEGER MXTRAK, MXPROT
 505       PARAMETER (MXTRAK=5000,MXPROT=5000)
 506       INTEGER NTRAK, IERR, MODE
 507       DOUBLE PRECISION COSR,PU(3,MXTRAK),PP(4,MXTRAK),VSEED(3)
 508       LOGICAL UNSTBL
 509       LOGICAL JETLIS(MXPROT,MXTRAK)
 510       INTEGER NJET
 511       DOUBLE PRECISION  PJ(4,MXPROT)
 512 *** Using VSEED as a trial axis , look for a stable jet.
 513 *** Check stable jets against those already found and add to PJ.
 514 *** Will try up to MXITER iterations to get a stable set of particles
 515 *** in the cone.
 516       INTEGER MU,N,ITER
 517       LOGICAL PXSAME,PXNEW,OK
 518       LOGICAL NEWLIS(MXTRAK),OLDLIS(MXTRAK)
 519       DOUBLE PRECISION OAXIS(3),NAXIS(3),PNEW(4)
 520       INTEGER MXITER
 521       PARAMETER(MXITER = 30)
 522 *
 523       DO 100 MU=1,3
 524         OAXIS(MU) = VSEED(MU)
 525 100   CONTINUE
 526       DO 110 N = 1,NTRAK
 527         OLDLIS(N) = .FALSE.
 528 110   CONTINUE
 529       DO 120 ITER = 1,MXITER
 530         CALL PXTRY(MODE,COSR,NTRAK,PU,PP,OAXIS,NAXIS,PNEW,NEWLIS,OK)
 531 *** Return immediately if there were no particles in the cone.
 532        IF (.NOT.OK) THEN
 533          RETURN
 534        ENDIF
 535        IF(PXSAME(NEWLIS,OLDLIS,NTRAK)) THEN
 536 *** We have a stable jet.
 537              IF (PXNEW(NEWLIS,JETLIS,NTRAK,NJET)) THEN
 538 *** And the jet is a new one. So add it to our arrays.
 539 *** Check arrays are big anough...
 540              IF (NJET .EQ. MXPROT) THEN
 541              WRITE (6,*) ' PXCONE:  Found more than MXPROT proto-jets'
 542                 IERR = -1
 543                 RETURN
 544              ENDIF
 545                NJET = NJET + 1
 546                DO 130 N = 1,NTRAK
 547                  JETLIS(NJET,N) = NEWLIS(N)
 548 130            CONTINUE
 549                DO 140 MU=1,4
 550                  PJ(MU,NJET)=PNEW(MU)
 551 140          CONTINUE
 552              ENDIF
 553              RETURN
 554        ENDIF
 555 *** The jet was not stable, so we iterate again
 556        DO 150 N=1,NTRAK
 557          OLDLIS(N)=NEWLIS(N)
 558 150    CONTINUE
 559        DO 160 MU=1,3
 560          OAXIS(MU)=NAXIS(MU)
 561 160    CONTINUE
 562 120   CONTINUE
 563       UNSTBL = .TRUE.
 564       RETURN
 565       END
 566 *CMZ :  1.06/00 28/02/94  15.44.44  by  P. Schleper
 567 *-- Author :
 568 C-----------------------------------------------------------------------
 569       SUBROUTINE PXSORV(N,A,K,OPT)
 570 C     Sort A(N) into ascending order
 571 C     OPT = 'I' : return index array K only
 572 C     OTHERWISE : return sorted A and index array K
 573 C-----------------------------------------------------------------------
 574       INTEGER NMAX
 575       PARAMETER (NMAX=5000)
 576 *
 577 *      INTEGER N,I,J,K(N),IL(NMAX),IR(NMAX)
 578 *LUND
 579       INTEGER N,I,J,K(NMAX),IL(NMAX),IR(NMAX)
 580       CHARACTER OPT
 581 *
 582 *      DOUBLE PRECISION A(N),B(NMAX)
 583       DOUBLE PRECISION A(NMAX),B(NMAX)
 584 *LUND
 585       IF (N.GT.NMAX) STOP 'Sorry, not enough room in Mike''s PXSORV'
 586       IL(1)=0
 587       IR(1)=0
 588       DO 10 I=2,N
 589       IL(I)=0
 590       IR(I)=0
 591       J=1
 592    2  IF(A(I).GT.A(J)) GO TO 5
 593    3  IF(IL(J).EQ.0) GO TO 4
 594       J=IL(J)
 595       GO TO 2
 596    4  IR(I)=-J
 597       IL(J)=I
 598       GO TO 10
 599    5  IF(IR(J).LE.0) GO TO 6
 600       J=IR(J)
 601       GO TO 2
 602    6  IR(I)=IR(J)
 603       IR(J)=I
 604   10  CONTINUE
 605       I=1
 606       J=1
 607       GO TO 8
 608   20  J=IL(J)
 609    8  IF(IL(J).GT.0) GO TO 20
 610    9  K(I)=J
 611       B(I)=A(J)
 612       I=I+1
 613       IF(IR(J)) 12,30,13
 614   13  J=IR(J)
 615       GO TO 8
 616   12  J=-IR(J)
 617       GO TO 9
 618   30  IF(OPT.EQ.'I') RETURN
 619       DO 31 I=1,N
 620   31  A(I)=B(I)
 621  999  END
 622
 623 *********************************************************************
 624 *CMZ :  2.00/00 10/01/95  10.17.57  by  P. Schleper
 625 *CMZ :  1.06/00 14/03/94  15.37.44  by  P. Schleper
 626 *-- Author :
 627 *
 628 C+DECK,PXTRY.
 629        SUBROUTINE PXTRY(MODE,COSR,NTRAK,PU,PP,OAXIS,NAXIS,
 630      +                  PNEW,NEWLIS,OK)
 631 *
 632 C+SEQ,DECLARE.
 633       INTEGER MXTRAK
 634       PARAMETER (MXTRAK=5000)
 635        INTEGER NTRAK,MODE
 636 *** Note that although PU and PP are assumed to be 2d arrays, they
 637 *** are used as 1d in this routine for efficiency
 638        DOUBLE PRECISION COSR,PU(3*MXTRAK),PP(4*MXTRAK),OAXIS(3),PXMDPI
 639        LOGICAL OK
 640        LOGICAL NEWLIS(MXTRAK)
 641        DOUBLE PRECISION NAXIS(3),PNEW(4)
 642 *** Finds all particles in cone of size COSR about OAXIS direction.
 643 *** Calculates 4-momentum sum of all particles in cone (PNEW) , and
 644 *** returns this as new jet axis NAXIS (Both unit Vectors)
 645        INTEGER N,MU,NPU,NPP
 646        DOUBLE PRECISION COSVAL,NORMSQ,NORM
 647 *
 648        OK = .FALSE.
 649        DO 100 MU=1,4
 650           PNEW(MU)=0.0
 651 100    CONTINUE
 652        NPU=-3
 653        NPP=-4
 654        DO 110 N=1,NTRAK
 655           NPU=NPU+3
 656           NPP=NPP+4
 657           IF (MODE.NE.2) THEN
 658              COSVAL=0.0
 659              DO 120 MU=1,3
 660                 COSVAL=COSVAL+OAXIS(MU)*PU(MU+NPU)
 661 120          CONTINUE
 662           ELSE
 663              IF (ABS(PU(1+NPU)).GE.20.OR.ABS(OAXIS(1)).GE.20) THEN
 664                 COSVAL=-1000
 665              ELSE
 666                 COSVAL=1-
 667      +           ((OAXIS(1)-PU(1+NPU))**2+PXMDPI(OAXIS(2)-PU(2+NPU))**2)
 668              ENDIF
 669           ENDIF
 670           IF (COSVAL.GE.COSR)THEN
 671              NEWLIS(N) = .TRUE.
 672              OK = .TRUE.
 673              IF (MODE.NE.2) THEN
 674                 DO 130 MU=1,4
 675                    PNEW(MU) = PNEW(MU) + PP(MU+NPP)
 676 130             CONTINUE
 677              ELSE
 678                 PNEW(1)=PNEW(1)
 679      +              + PP(4+NPP)/(PP(4+NPP)+PNEW(4))*(PP(1+NPP)-PNEW(1))
 680                 PNEW(2)=PNEW(2)
 681      +              + PP(4+NPP)/(PP(4+NPP)+PNEW(4))
 682      +               *PXMDPI(PP(2+NPP)-PNEW(2))
 683                 PNEW(4)=PNEW(4)+PP(4+NPP)
 684              ENDIF
 685           ELSE
 686              NEWLIS(N)=.FALSE.
 687           ENDIF
 688 110   CONTINUE
 689 *** If there are particles in the cone, calc new jet axis
 690        IF (OK) THEN
 691           IF (MODE.NE.2) THEN
 692              NORMSQ = 0.0
 693              DO 140 MU = 1,3
 694                 NORMSQ = NORMSQ + PNEW(MU)**2
 695 140          CONTINUE
 696              NORM = SQRT(NORMSQ)
 697           ELSE
 698              NORM = 1
 699           ENDIF
 700           DO 150 MU=1,3
 701              NAXIS(MU) = PNEW(MU)/NORM
 702 150       CONTINUE
 703        ENDIF
 704        RETURN
 705        END
 706
 707 *********************************************************************
 708 *CMZ :  2.00/00 10/01/95  10.17.57  by  P. Schleper
 709 *CMZ :  1.06/00 28/02/94  15.44.44  by  P. Schleper
 710 *-- Author :
 711 C+DECK,PXUVEC.
 712 *
 713        SUBROUTINE PXUVEC(NTRAK,PP,PU,IERR)
 714 *
 715 *** Routine to calculate unit vectors PU of all particles PP
 716 C+SEQ,DECLARE.
 717       INTEGER MXTRAK
 718       PARAMETER (MXTRAK=5000)
 719       INTEGER NTRAK, IERR
 720       DOUBLE PRECISION PP(4,MXTRAK)
 721       DOUBLE PRECISION PU(3,MXTRAK)
 722       INTEGER N,MU
 723       DOUBLE PRECISION MAG
 724        DO 100 N=1,NTRAK
 725           MAG=0.0
 726           DO 110 MU=1,3
 727              MAG=MAG+PP(MU,N)**2
 728 110       CONTINUE
 729           MAG=SQRT(MAG)
 730           IF (MAG.EQ.0.0) THEN
 731              WRITE(6,*)' PXCONE: An input particle has zero mod(p)'
 732              IERR=-1
 733              RETURN
 734           ENDIF
 735           DO 120 MU=1,3
 736            PU(MU,N)=PP(MU,N)/MAG
 737 120       CONTINUE
 738 100    CONTINUE
 739        RETURN
 740        END
 741 *CMZ :  1.06/00 28/02/94  15.44.44  by  P. Schleper
 742 *-- Author :
 743 C-----------------------------------------------------------------------
 744       SUBROUTINE PXZERI(N,A)
 745       INTEGER I,N,A(N)
 746       DO 10 I=1,N
 747         A(I)=0
 748  10   CONTINUE
 749       END
 750 *CMZ :  1.06/00 28/02/94  15.44.44  by  P. Schleper
 751 *-- Author :
 752 C-----------------------------------------------------------------------
 753 C     This is a set of routines written by Mike Seymour to provide the
 754 C     services presumably normally provided by standard OPAL routines
 755 C     PXZERV zeroes a vector
 756 C     PXZERI zeroes a vector of integers
 757 C     PXNORV normalizes a vector
 758 C     PXADDV adds two vectors
 759 C     PXSORV sorts a vector (copied from HERWIG)
 760 C     PXANG3 finds the angle (and its cosine) between two vectors
 761 C     PXMDPI moves its argument onto the range [-pi,pi)
 762 C-----------------------------------------------------------------------
 763       SUBROUTINE PXZERV(N,A)
 764       INTEGER I,N
 765       DOUBLE PRECISION A(N)
 766       DO 10 I=1,N
 767         A(I)=0
 768  10   CONTINUE
 769       END
 770 *-- Author :
 771 C-----------------------------------------------------------------------
 772       SUBROUTINE PXADDV(N,A,B,C,ITERR)
 773       INTEGER I,N,ITERR
 774       DOUBLE PRECISION A(N),B(N),C(N)
 775       DO 10 I=1,N
 776         C(I)=A(I)+B(I)
 777  10   CONTINUE
 778       END
 779 *CMZ :  1.06/00 28/02/94  15.44.44  by  P. Schleper
 780 *-- Author :
 781 C-----------------------------------------------------------------------
 782       SUBROUTINE PXANG3(A,B,COST,THET,ITERR)
 783       INTEGER ITERR
 784       DOUBLE PRECISION A(3),B(3),C,COST,THET
 785       C=(A(1)**2+A(2)**2+A(3)**2)*(B(1)**2+B(2)**2+B(3)**2)
 786       IF (C.LE.0) RETURN
 787       C=1/SQRT(C)
 788       COST=(A(1)*B(1)+A(2)*B(2)+A(3)*B(3))*C
 789       THET=ACOS(COST)
 790       END
 791 *CMZ :  1.06/00 14/03/94  15.41.57  by  P. Schleper
 792 *-- Author :    P. Schleper   28/02/94
 793        LOGICAL FUNCTION PXNEW(TSTLIS,JETLIS,NTRAK,NJET)
 794 *
 795       INTEGER MXTRAK,MXPROT
 796       PARAMETER (MXTRAK=5000,MXPROT=5000)
 797        INTEGER NTRAK,NJET
 798 *** Note that although JETLIS is assumed to be a 2d array, it
 799 *** it is used as 1d in this routine for efficiency
 800        LOGICAL TSTLIS(MXTRAK),JETLIS(MXPROT*MXTRAK)
 801 *** Checks to see if TSTLIS entries correspond to a jet already found
 802 *** and entered in JETLIS
 803        INTEGER N, I, IN
 804        LOGICAL MATCH
 805 *
 806        PXNEW = .TRUE.
 807        DO 100 I = 1,NJET
 808           MATCH = .TRUE.
 809           IN=I-MXPROT
 810           DO 110 N = 1,NTRAK
 811             IN=IN+MXPROT
 812             IF(TSTLIS(N).NEQV.JETLIS(IN)) THEN
 813              MATCH = .FALSE.
 814              GO TO 100
 815             ENDIF
 816 110       CONTINUE
 817           IF (MATCH) THEN
 818            PXNEW = .FALSE.
 819            RETURN
 820           ENDIF
 821 100    CONTINUE
 822        RETURN
 823        END
 824 *CMZ :  1.06/00 14/03/94  15.41.57  by  P. Schleper
 825 *-- Author :    P. Schleper   28/02/94
 826        LOGICAL FUNCTION PXSAME(LIST1,LIST2,N)
 827 *
 828        LOGICAL LIST1(*),LIST2(*)
 829        INTEGER N
 830 *** Returns T if the first N elements of LIST1 are the same as the
 831 *** first N elements of LIST2.
 832        INTEGER I
 833 *
 834        PXSAME = .TRUE.
 835        DO 100 I = 1,N
 836         IF ( LIST1(I).NEQV.LIST2(I) ) THEN
 837           PXSAME = .FALSE.
 838           RETURN
 839         ENDIF
 840 100    CONTINUE
 841        RETURN
 842        END
 843 *CMZ :  1.06/00 28/02/94  15.44.44  by  P. Schleper
 844 *-- Author :
 845 C-----------------------------------------------------------------------
 846       FUNCTION PXMDPI(PHI)
 847       IMPLICIT NONE
 848 C---RETURNS PHI, MOVED ONTO THE RANGE [-PI,PI)
 849       DOUBLE PRECISION PXMDPI,PHI,PI,TWOPI,THRPI,EPS
 850       PARAMETER (PI=3.141592654,TWOPI=6.283185307,THRPI=9.424777961)
 851       PARAMETER (EPS=1E-15)
 852       PXMDPI=PHI
 853       IF (PXMDPI.LE.PI) THEN
 854         IF (PXMDPI.GT.-PI) THEN
 855           GOTO 100
 856         ELSEIF (PXMDPI.GT.-THRPI) THEN
 857           PXMDPI=PXMDPI+TWOPI
 858         ELSE
 859           PXMDPI=-MOD(PI-PXMDPI,TWOPI)+PI
 860         ENDIF
 861       ELSEIF (PXMDPI.LE.THRPI) THEN
 862         PXMDPI=PXMDPI-TWOPI
 863       ELSE
 864         PXMDPI=MOD(PI+PXMDPI,TWOPI)-PI
 865       ENDIF
 866  100  IF (ABS(PXMDPI).LT.EPS) PXMDPI=0
 867       END