FastJet 3.0.2
pxcone.f
00001       SUBROUTINE PXCONE(MODE,NTRAK,ITKDM,PTRAK,CONER,EPSLON,OVLIM,
00002      +                   MXJET,NJET,PJET,IPASS,IJMUL,IERR)
00003 *.*********************************************************
00004 *. ------
00005 *. PXCONE
00006 *. ------
00007 *.
00008 *. Code downloaded from the following web page 
00009 *.
00010 *.   http://aliceinfo.cern.ch/alicvs/viewvc/JETAN/pxcone.F?view=markup&pathrev=v4-05-04
00011 *. 
00012 *. on 17/10/2006 by G. Salam. Permission subsequently granted by Michael
00013 *. H. Seymour (on behalf of the PxCone authors) for this code to be
00014 *. distributed together with FastJet under the terms of the GNU Public
00015 *. License v2 (see the file COPYING in the main FastJet directory).
00016 *.
00017 *.********** Pre Release Version 26.2.93
00018 *.
00019 *. Driver for the Cone  Jet finding algorithm of L.A. del Pozo.
00020 *. Based on algorithm from D.E. Soper.
00021 *. Finds jets inside cone of half angle CONER with energy > EPSLON.
00022 *. Jets which receive more than a fraction OVLIM of their energy from
00023 *. overlaps with other jets are excluded.
00024 *. Output jets are ordered in energy.
00025 *. If MODE.EQ.2 momenta are stored as (eta,phi,<empty>,pt)
00026 *. Usage     :
00027 *.
00028 *.      INTEGER  ITKDM,MXTRK
00029 *.      PARAMETER  (ITKDM=4.or.more,MXTRK=1.or.more)
00030 *.      INTEGER  MXJET, MXTRAK, MXPROT
00031 *.      PARAMETER  (MXJET=10,MXTRAK=500,MXPROT=500)
00032 *.      INTEGER  IPASS (MXTRAK),IJMUL (MXJET)
00033 *.      INTEGER  NTRAK,NJET,IERR,MODE
00034 *.      DOUBLE PRECISION  PTRAK (ITKDM,MXTRK),PJET (5,MXJET)
00035 *.      DOUBLE PRECISION  CONER, EPSLON, OVLIM
00036 *.      NTRAK = 1.to.MXTRAK
00037 *.      CONER   = ...
00038 *.      EPSLON  = ...
00039 *.      OVLIM   = ...
00040 *.      CALL PXCONE (MODE,NTRAK,ITKDM,PTRAK,CONER,EPSLON,OVLIM,MXJET,
00041 *.     +             NJET,PJET,IPASS,IJMUL,IERR)
00042 *.
00043 *. INPUT     :  MODE      1=>e+e-, 2=>hadron-hadron
00044 *. INPUT     :  NTRAK     Number of particles
00045 *. INPUT     :  ITKDM     First dimension of PTRAK array
00046 *. INPUT     :  PTRAK     Array of particle 4-momenta (Px,Py,Pz,E)
00047 *. INPUT     :  CONER     Cone size (half angle) in radians
00048 *. INPUT     :  EPSLON    Minimum Jet energy (GeV)
00049 *. INPUT     :  OVLIM     Maximum fraction of overlap energy in a jet
00050 *. INPUT     :  MXJET     Maximum possible number of jets
00051 *. OUTPUT    :  NJET      Number of jets found
00052 *. OUTPUT    :  PJET      5-vectors of jets
00053 *. OUTPUT    :  IPASS(k)  Particle k belongs to jet number IPASS(k)
00054 *.                        IPASS = -1 if not assosciated to a jet
00055 *. OUTPUT    :  IJMUL(i)  Jet i contains IJMUL(i) particles
00056 *. OUTPUT    :  IERR      = 0 if all is OK ;   = -1 otherwise
00057 *.
00058 *. CALLS     : PXSEAR, PXSAME, PXNEW, PXTRY, PXORD, PXUVEC, PXOLAP
00059 *. CALLED    : User
00060 *.
00061 *. AUTHOR    :  L.A. del Pozo
00062 *. CREATED   :  26-Feb-93
00063 *. LAST MOD  :   2-Mar-93
00064 *.
00065 *. Modification Log.
00066 *. 25-Feb-07: G P Salam   - fix bugs concerning 2pi periodicity in eta phi mode
00067 *.                        - added commented code to get consistent behaviour
00068 *.                          regardless of particle order (replaces n-way 
00069 *.                          midpoints with 2-way midpoints however...)
00070 *. 2-Jan-97: M Wobisch    - fix bug concerning COS2R in eta phi mode
00071 *. 4-Apr-93: M H Seymour  - Change 2d arrays to 1d in PXTRY & PXNEW
00072 *. 2-Apr-93: M H Seymour  - Major changes to add boost-invariant mode
00073 *. 1-Apr-93: M H Seymour  - Increase all array sizes
00074 *. 30-Mar-93: M H Seymour - Change all REAL variables to DOUBLE PRECISION
00075 *. 30-Mar-93: M H Seymour - Change OVLIM into an input parameter
00076 *. 2-Mar-93: L A del Pozo - Fix Bugs in PXOLAP
00077 *. 1-Mar-93: L A del Pozo - Remove Cern library routine calls
00078 *. 1-Mar-93: L A del Pozo - Add Print out of welcome and R and Epsilon
00079 *.
00080 *.*********************************************************
00081 C+SEQ,DECLARE.
00082 *** External Arrays
00083       INTEGER  ITKDM,MXJET,NTRAK,NJET,IERR,MODE
00084       INTEGER  IPASS (NTRAK),IJMUL (MXJET)
00085       DOUBLE PRECISION  PTRAK (ITKDM,NTRAK),PJET (5,MXJET), 
00086      +     CONER, EPSLON, OVLIM
00087 *** Internal Arrays
00088       INTEGER MXPROT, MXTRAK
00089       PARAMETER (MXPROT=5000, MXTRAK=5000)
00090       DOUBLE PRECISION PP(4,MXTRAK), PU(3,MXTRAK), PJ(4,MXPROT)
00091       LOGICAL JETLIS(MXPROT,MXTRAK)
00092 *** Used in the routine.
00093       DOUBLE PRECISION COSR,COS2R, VSEED(3), VEC1(3), VEC2(3),PTSQ,PPSQ,
00094      +     COSVAL,PXMDPI
00095 cMWobisch
00096       DOUBLE PRECISION RSEP
00097 cMWobisch
00098       LOGICAL UNSTBL
00099       INTEGER I,J,N,MU,N1,N2, ITERR, NJTORG
00100       INTEGER NCALL, NPRINT
00101       DOUBLE PRECISION ROLD, EPSOLD, OVOLD
00102       SAVE NCALL,NPRINT,ROLD, EPSOLD, OVOLD
00103       DATA NCALL,NPRINT /0,0/
00104       DATA ROLD,EPSOLD,OVOLD/0.,0.,0./
00105 
00106 cMWobisch
00107 c***************************************
00108       RSEP  = 2D0
00109 c***************************************
00110 cMWobisch
00111       IERR=0
00112 *
00113 *** INITIALIZE
00114       IF(NCALL.LE.0)  THEN
00115          ROLD = 0.
00116          EPSOLD = 0.
00117          OVOLD = 0.
00118       ENDIF
00119       NCALL = NCALL + 1
00120 *
00121 *** Print welcome and Jetfinder parameters
00122       IF((CONER.NE.ROLD .OR. EPSLON.NE.EPSOLD .OR. OVLIM.NE.OVOLD)
00123      +     .AND. NPRINT.LE.10) THEN
00124          WRITE (6,*)
00125          WRITE (6,*) ' *********** PXCONE: Cone Jet-finder ***********'
00126          WRITE (6,*) '    Written by Luis Del Pozo of OPAL'
00127          WRITE (6,*) '    Modified for eta-phi by Mike Seymour'
00128          WRITE (6,*) '    Includes bug fixes by Wobisch, Salam'
00129          WRITE(6,1000)'   Cone Size R = ',CONER,' Radians'
00130          WRITE(6,1001)'   Min Jet energy Epsilon = ',EPSLON,' GeV'
00131          WRITE(6,1002)'   Overlap fraction parameter = ',OVLIM
00132          WRITE (6,*) '    PXCONE is not a supported product and is'
00133          WRITE (6,*) '    is provided for comparative purposes only'
00134          WRITE (6,*) ' ***********************************************'
00135 cMWobisch
00136          IF (RSEP .lt. 1.999) THEN
00137             WRITE(6,*) ' '
00138             WRITE (6,*) ' ******************************************'
00139             WRITE (6,*) ' ******************************************'
00140             WRITE(6,*) ' M Wobisch: private change !!!!!!!!!!!! '
00141             WRITE(6,*) '      Rsep is set to ',RSEP
00142             WRITE(6,*) ' this is ONLY meaningful in a NLO calculation'
00143             WRITE(6,*) '      ------------------------  '
00144             WRITE(6,*) '  please check what you''re doing!!'
00145             WRITE(6,*) '   or ask:  Markus.Wobisch@desy.de --'
00146             WRITE (6,*) ' ******************************************'
00147             WRITE (6,*) ' ******************************************'
00148             WRITE (6,*) ' ******************************************'
00149             WRITE(6,*) ' '
00150             WRITE(6,*) ' '
00151          ENDIF
00152 cMWobisch
00153 
00154          WRITE (6,*)
00155 1000     FORMAT(A18,F5.2,A10)
00156 1001     FORMAT(A29,F5.2,A5)
00157 1002     FORMAT(A33,F5.2)
00158          NPRINT = NPRINT + 1
00159          ROLD=CONER
00160          EPSOLD=EPSLON
00161          OVOLD=OVLIM
00162       ENDIF
00163 *
00164 *** Copy calling array PTRAK  to internal array PP(4,NTRAK)
00165 *
00166       IF (NTRAK .GT. MXTRAK) THEN
00167          WRITE (6,*) ' PXCONE: Ntrak too large: ',NTRAK
00168          IERR=-1
00169          RETURN
00170       ENDIF
00171       IF (MODE.NE.2) THEN
00172          DO  100 I=1, NTRAK
00173             DO  101 J=1,4
00174                PP(J,I)=PTRAK(J,I)
00175 101         CONTINUE
00176 100      CONTINUE
00177       ELSE
00178 *** Converting to eta,phi,pt if necessary
00179          DO  104 I=1,NTRAK
00180             PTSQ=PTRAK(1,I)**2+PTRAK(2,I)**2
00181             PPSQ=(SQRT(PTSQ+PTRAK(3,I)**2)+ABS(PTRAK(3,I)))**2
00182             IF (PTSQ.LE.4.25E-18*PPSQ) THEN
00183                PP(1,I)=20
00184             ELSE
00185                PP(1,I)=0.5*LOG(PPSQ/PTSQ)
00186             ENDIF
00187             PP(1,I)=SIGN(PP(1,I),PTRAK(3,I))
00188             IF (PTSQ.EQ.0) THEN
00189                PP(2,I)=0
00190             ELSE
00191                PP(2,I)=ATAN2(PTRAK(2,I),PTRAK(1,I))
00192             ENDIF
00193             PP(3,I)=0
00194             PP(4,I)=SQRT(PTSQ)
00195             PU(1,I)=PP(1,I)
00196             PU(2,I)=PP(2,I)
00197             PU(3,I)=PP(3,I)
00198 104      CONTINUE
00199       ENDIF
00200 *
00201 *** Zero output variables
00202 *
00203       NJET=0
00204       DO 102 I = 1, NTRAK
00205          DO 103 J = 1, MXPROT
00206            JETLIS(J,I) = .FALSE.
00207 103      CONTINUE
00208 102   CONTINUE
00209       CALL PXZERV(4*MXPROT,PJ)
00210       CALL PXZERI(MXJET,IJMUL)
00211 *
00212       IF (MODE.NE.2) THEN
00213          COSR = COS(CONER)
00214          COS2R = COS(CONER)
00215       ELSE
00216 *** Purely for convenience, work in terms of 1-R**2
00217          COSR = 1-CONER**2
00218 cMW -- select Rsep: 1-(Rsep*CONER)**2
00219          COS2R =  1-(RSEP*CONER)**2
00220 cORIGINAL         COS2R =  1-(2*CONER)**2
00221       ENDIF
00222       UNSTBL = .FALSE.
00223       IF (MODE.NE.2) THEN
00224          CALL PXUVEC(NTRAK,PP,PU,IERR)
00225          IF (IERR .NE. 0) RETURN
00226       ENDIF
00227 *** Look for jets using particle diretions as seed axes
00228 *
00229       DO 110 N = 1,NTRAK
00230         DO 120 MU = 1,3
00231           VSEED(MU) = PU(MU,N)
00232 120     CONTINUE
00233         CALL PXSEAR(MODE,COSR,NTRAK,PU,PP,VSEED,
00234      &                   NJET,JETLIS,PJ,UNSTBL,IERR)
00235          IF (IERR .NE. 0) RETURN
00236 110   CONTINUE
00237 
00238 cMW - for Rsep=1 goto 145
00239 c      GOTO 145
00240 
00241 *** Now look between all pairs of jets as seed axes.
00242 c      NJTORG = NJET           ! GPS -- to get consistent behaviour (2-way midpnts)
00243 c      DO 140 N1 = 1,NJTORG-1  ! GPS -- to get consistent behaviour (2-way midpnts)
00244       DO 140 N1 = 1,NJET-1  
00245          VEC1(1)=PJ(1,N1)
00246          VEC1(2)=PJ(2,N1)
00247          VEC1(3)=PJ(3,N1)
00248          IF (MODE.NE.2) CALL PXNORV(3,VEC1,VEC1,ITERR)
00249 C         DO 150 N2 = N1+1,NJTORG ! GPS -- to get consistent behaviour
00250          DO 150 N2 = N1+1,NJET
00251             VEC2(1)=PJ(1,N2)
00252             VEC2(2)=PJ(2,N2)
00253             VEC2(3)=PJ(3,N2)
00254             IF (MODE.NE.2) CALL PXNORV(3,VEC2,VEC2,ITERR)
00255             CALL PXADDV(3,VEC1,VEC2,VSEED,ITERR)
00256             IF (MODE.NE.2) THEN
00257                CALL PXNORV(3,VSEED,VSEED,ITERR)
00258             ELSE
00259                VSEED(1)=VSEED(1)/2
00260                !VSEED(2)=VSEED(2)/2
00261                ! GPS 25/02/07
00262                VSEED(2)=PXMDPI(VEC1(2)+0.5d0*PXMDPI(VEC2(2)-VEC1(2)))
00263             ENDIF
00264 C---ONLY BOTHER IF THEY ARE BETWEEN 1 AND 2 CONE RADII APART
00265             IF (MODE.NE.2) THEN
00266               COSVAL=VEC1(1)*VEC2(1)+VEC1(2)*VEC2(2)+VEC1(3)*VEC2(3)
00267             ELSE
00268                IF (ABS(VEC1(1)).GE.20.OR.ABS(VEC2(1)).GE.20) THEN
00269                   COSVAL=-1000
00270                ELSE
00271                   COSVAL=1-
00272      +              ((VEC1(1)-VEC2(1))**2+PXMDPI(VEC1(2)-VEC2(2))**2)
00273                ENDIF
00274             ENDIF
00275 
00276             IF (COSVAL.LE.COSR.AND.COSVAL.GE.COS2R)
00277      +           CALL PXSEAR(MODE,COSR,NTRAK,PU,PP,VSEED,NJET,
00278      +           JETLIS,PJ,UNSTBL,IERR)
00279 c            CALL PXSEAR(MODE,COSR,NTRAK,PU,PP,VSEED,NJET,
00280 c     +           JETLIS,PJ,UNSTBL,IERR)
00281             IF (IERR .NE. 0) RETURN
00282 150      CONTINUE
00283 140   CONTINUE
00284       IF (UNSTBL) THEN
00285         IERR=-1
00286         WRITE (6,*) ' PXCONE: Too many iterations to find a proto-jet'
00287         RETURN
00288       ENDIF
00289 
00290  145  CONTINUE
00291 *** Now put the jet list into order by jet energy, eliminating jets
00292 *** with energy less than EPSLON.
00293        CALL PXORD(EPSLON,NJET,NTRAK,JETLIS,PJ)
00294 *
00295 *** Take care of jet overlaps
00296        CALL PXOLAP(MODE,NJET,NTRAK,JETLIS,PJ,PP,OVLIM)
00297 *
00298 *** Order jets again as some have been eliminated, or lost energy.
00299        CALL PXORD(EPSLON,NJET,NTRAK,JETLIS,PJ)
00300 *
00301 *** All done!, Copy output into output arrays
00302       IF (NJET .GT. MXJET) THEN
00303          WRITE (6,*) ' PXCONE:  Found more than MXJET jets'
00304          IERR=-1
00305          GOTO 99
00306       ENDIF
00307       IF (MODE.NE.2) THEN
00308          DO 300 I=1, NJET
00309             DO 310 J=1,4
00310                PJET(J,I)=PJ(J,I)
00311 310         CONTINUE
00312 300      CONTINUE
00313       ELSE
00314          DO 315 I=1, NJET
00315             PJET(1,I)=PJ(4,I)*COS(PJ(2,I))
00316             PJET(2,I)=PJ(4,I)*SIN(PJ(2,I))
00317             PJET(3,I)=PJ(4,I)*SINH(PJ(1,I))
00318             PJET(4,I)=PJ(4,I)*COSH(PJ(1,I))
00319  315     CONTINUE
00320       ENDIF
00321       DO 320 I=1, NTRAK
00322          IPASS(I)=-1
00323          DO 330 J=1, NJET
00324             IF (JETLIS(J,I)) THEN
00325                IJMUL(J)=IJMUL(J)+1
00326                IPASS(I)=J
00327             ENDIF
00328 330      CONTINUE
00329 320   CONTINUE
00330 99    RETURN
00331       END
00332 *CMZ :  1.06/00 28/02/94  15.44.44  by  P. Schleper
00333 *-- Author :
00334 C-----------------------------------------------------------------------
00335       SUBROUTINE PXNORV(N,A,B,ITERR)
00336       INTEGER I,N,ITERR
00337       DOUBLE PRECISION A(N),B(N),C
00338       C=0
00339       DO 10 I=1,N
00340         C=C+A(I)**2
00341  10   CONTINUE
00342       IF (C.LE.0) RETURN
00343       C=1/SQRT(C)
00344       DO 20 I=1,N
00345         B(I)=A(I)*C
00346  20   CONTINUE
00347       END
00348 *CMZ :  2.00/00 10/01/95  10.17.57  by  P. Schleper
00349 *CMZ :  1.06/00 15/03/94  12.17.46  by  P. Schleper
00350 *-- Author :
00351 *
00352 C+DECK,PXOLAP.
00353       SUBROUTINE PXOLAP(MODE,NJET,NTRAK,JETLIS,PJ,PP,OVLIM)
00354 *
00355 *** Looks for particles assigned to more than 1 jet, and reassigns them
00356 *** If more than a fraction OVLIM of a jet's energy is contained in
00357 *** higher energy jets, that jet is neglected.
00358 *** Particles assigned to the jet closest in angle (a la CDF, Snowmass).
00359 C+SEQ,DECLARE.
00360       INTEGER MXTRAK, MXPROT
00361       PARAMETER (MXTRAK=5000,MXPROT=5000)
00362       INTEGER NJET, NTRAK, MODE
00363       LOGICAL JETLIS(MXPROT,MXTRAK)
00364       DOUBLE PRECISION PJ(4,MXPROT),PP(4,MXTRAK),PXMDPI
00365       INTEGER I,J,N,MU
00366       LOGICAL OVELAP
00367       DOUBLE PRECISION EOVER
00368       DOUBLE PRECISION OVLIM
00369       INTEGER ITERR, IJMIN, IJET(MXPROT), NJ
00370       DOUBLE PRECISION VEC1(3), VEC2(3), COST, THET, THMIN
00371       DATA IJMIN/0/
00372 *
00373       IF (NJET.LE.1) RETURN
00374 *** Look for jets with large overlaps with higher energy jets.
00375       DO 100 I = 2,NJET
00376 *** Find overlap energy between jets I and all higher energy jets.
00377        EOVER = 0.0
00378        DO 110 N = 1,NTRAK
00379          OVELAP = .FALSE.
00380          DO 120 J= 1,I-1
00381            IF (JETLIS(I,N).AND.JETLIS(J,N)) THEN
00382             OVELAP = .TRUE.
00383            ENDIF
00384 120      CONTINUE
00385          IF (OVELAP) THEN
00386            EOVER = EOVER + PP(4,N)
00387          ENDIF
00388 110     CONTINUE
00389 *** Is the fraction of energy shared larger than OVLIM?
00390         IF (EOVER.GT.OVLIM*PJ(4,I)) THEN
00391 *** De-assign all particles from Jet I
00392             DO 130 N = 1,NTRAK
00393               JETLIS(I,N) = .FALSE.
00394 130         CONTINUE
00395          ENDIF
00396 100   CONTINUE
00397 *** Now there are no big overlaps, assign every particle in
00398 *** more than 1 jet to the closet jet.
00399 *** Any particles now in more than 1 jet are assigned to the CLOSET
00400 *** jet (in angle).
00401       DO 140 I=1,NTRAK
00402          NJ=0
00403          DO 150 J=1, NJET
00404          IF(JETLIS(J,I)) THEN
00405             NJ=NJ+1
00406             IJET(NJ)=J
00407          ENDIF
00408 150      CONTINUE
00409          IF (NJ .GT. 1) THEN
00410 *** Particle in > 1 jet - calc angles...
00411             VEC1(1)=PP(1,I)
00412             VEC1(2)=PP(2,I)
00413             VEC1(3)=PP(3,I)
00414             THMIN=0.
00415             DO 160 J=1,NJ
00416                VEC2(1)=PJ(1,IJET(J))
00417                VEC2(2)=PJ(2,IJET(J))
00418                VEC2(3)=PJ(3,IJET(J))
00419                IF (MODE.NE.2) THEN
00420                   CALL PXANG3(VEC1,VEC2,COST,THET,ITERR)
00421                ELSE
00422                   THET=(VEC1(1)-VEC2(1))**2+PXMDPI(VEC1(2)-VEC2(2))**2
00423                ENDIF
00424                IF (J .EQ. 1) THEN
00425                   THMIN=THET
00426                   IJMIN=IJET(J)
00427                ELSEIF (THET .LT. THMIN) THEN
00428                   THMIN=THET
00429                   IJMIN=IJET(J)
00430                ENDIF
00431 160         CONTINUE
00432 *** Assign track to IJMIN
00433             DO 170 J=1,NJET
00434                JETLIS(J,I) = .FALSE.
00435 170         CONTINUE
00436             JETLIS(IJMIN,I)=.TRUE.
00437          ENDIF
00438 140   CONTINUE
00439 *** Recompute PJ
00440       DO 200 I = 1,NJET
00441         DO 210 MU = 1,4
00442           PJ(MU,I) = 0.0
00443 210     CONTINUE
00444         DO 220 N = 1,NTRAK
00445           IF( JETLIS(I,N) ) THEN
00446              IF (MODE.NE.2) THEN
00447                 DO 230 MU = 1,4
00448                    PJ(MU,I) = PJ(MU,I) + PP(MU,N)
00449 230             CONTINUE
00450              ELSE
00451                 PJ(1,I)=PJ(1,I)
00452      +               + PP(4,N)/(PP(4,N)+PJ(4,I))*(PP(1,N)-PJ(1,I))
00453 c GPS 25/02/07
00454                 PJ(2,I)=PXMDPI(PJ(2,I)
00455      +             + PP(4,N)/(PP(4,N)+PJ(4,I))*PXMDPI(PP(2,N)-PJ(2,I)))
00456 c                PJ(2,I)=PJ(2,I)
00457 c     +               + PP(4,N)/(PP(4,N)+PJ(4,I))*PXMDPI(PP(2,N)-PJ(2,I))
00458                 PJ(4,I)=PJ(4,I)+PP(4,N)
00459              ENDIF
00460           ENDIF
00461 220     CONTINUE
00462 200   CONTINUE
00463       RETURN
00464       END
00465 *CMZ :  2.00/00 10/01/95  10.17.57  by  P. Schleper
00466 *CMZ :  1.06/00 14/03/94  15.37.45  by  P. Schleper
00467 *-- Author :
00468 *
00469 C+DECK,PXORD.
00470        SUBROUTINE PXORD(EPSLON,NJET,NTRAK,JETLIS,PJ)
00471 *
00472 *** Routine to put jets into order and eliminate tose less than EPSLON
00473 C+SEQ,DECLARE.
00474       INTEGER MXTRAK,MXPROT
00475       PARAMETER (MXTRAK=5000,MXPROT=5000)
00476        INTEGER I, J, INDEX(MXPROT)
00477        DOUBLE PRECISION PTEMP(4,MXPROT), ELIST(MXPROT)
00478        INTEGER NJET,NTRAK
00479        LOGICAL JETLIS(MXPROT,MXTRAK)
00480        LOGICAL LOGTMP(MXPROT,MXTRAK)
00481        DOUBLE PRECISION EPSLON,PJ(4,MXPROT)
00482 *** Puts jets in order of energy: 1 = highest energy etc.
00483 *** Then Eliminate jets with energy below EPSLON
00484 *
00485 *** Copy input arrays.
00486       DO 100 I=1,NJET
00487          DO 110 J=1,4
00488             PTEMP(J,I)=PJ(J,I)
00489 110      CONTINUE
00490          DO 120 J=1,NTRAK
00491             LOGTMP(I,J)=JETLIS(I,J)
00492 120      CONTINUE
00493 100   CONTINUE
00494       DO 150 I=1,NJET
00495          ELIST(I)=PJ(4,I)
00496 150   CONTINUE
00497 *** Sort the energies...
00498       CALL PXSORV(NJET,ELIST,INDEX,'I')
00499 *** Fill PJ and JETLIS according to sort ( sort is in ascending order!!)
00500       DO 200 I=1, NJET
00501          DO 210 J=1,4
00502             PJ(J,I)=PTEMP(J,INDEX(NJET+1-I))
00503 210      CONTINUE
00504          DO 220 J=1,NTRAK
00505             JETLIS(I,J)=LOGTMP(INDEX(NJET+1-I),J)
00506 220      CONTINUE
00507 200   CONTINUE
00508 ** Jets are now in order
00509 *** Now eliminate jets with less than Epsilon energy
00510       DO 300, I=1, NJET
00511          IF (PJ(4,I) .LT. EPSLON) THEN
00512             NJET=NJET-1
00513             PJ(4,I)=0.
00514          ENDIF
00515 300   CONTINUE
00516       RETURN
00517       END
00518 
00519 ********************************************************************
00520 *CMZ :  2.00/00 10/01/95  10.17.57  by  P. Schleper
00521 *CMZ :  1.06/00 14/03/94  15.37.44  by  P. Schleper
00522 *-- Author :
00523 C+DECK,PXSEAR.
00524       SUBROUTINE PXSEAR(MODE,COSR,NTRAK,PU,PP,VSEED,NJET,
00525      +                JETLIS,PJ,UNSTBL,IERR)
00526 *
00527 C+SEQ,DECLARE.
00528       INTEGER MXTRAK, MXPROT
00529       PARAMETER (MXTRAK=5000,MXPROT=5000)
00530       INTEGER NTRAK, IERR, MODE
00531       DOUBLE PRECISION COSR,PU(3,MXTRAK),PP(4,MXTRAK),VSEED(3)
00532       LOGICAL UNSTBL
00533       LOGICAL JETLIS(MXPROT,MXTRAK)
00534       INTEGER NJET
00535       DOUBLE PRECISION  PJ(4,MXPROT)
00536 *** Using VSEED as a trial axis , look for a stable jet.
00537 *** Check stable jets against those already found and add to PJ.
00538 *** Will try up to MXITER iterations to get a stable set of particles
00539 *** in the cone.
00540       INTEGER MU,N,ITER
00541       LOGICAL PXSAME,PXNEW,OK
00542       LOGICAL NEWLIS(MXTRAK),OLDLIS(MXTRAK)
00543       DOUBLE PRECISION OAXIS(3),NAXIS(3),PNEW(4)
00544       INTEGER MXITER
00545       PARAMETER(MXITER = 30)
00546 *
00547       DO 100 MU=1,3
00548         OAXIS(MU) = VSEED(MU)
00549 100   CONTINUE
00550       DO 110 N = 1,NTRAK
00551         OLDLIS(N) = .FALSE.
00552 110   CONTINUE
00553       DO 120 ITER = 1,MXITER
00554         CALL PXTRY(MODE,COSR,NTRAK,PU,PP,OAXIS,NAXIS,PNEW,NEWLIS,OK)
00555 *** Return immediately if there were no particles in the cone.
00556        IF (.NOT.OK) THEN
00557          RETURN
00558        ENDIF
00559        IF(PXSAME(NEWLIS,OLDLIS,NTRAK)) THEN
00560 *** We have a stable jet.
00561              IF (PXNEW(NEWLIS,JETLIS,NTRAK,NJET)) THEN
00562 *** And the jet is a new one. So add it to our arrays.
00563 *** Check arrays are big anough...
00564              IF (NJET .EQ. MXPROT) THEN
00565              WRITE (6,*) ' PXCONE:  Found more than MXPROT proto-jets'
00566                 IERR = -1
00567                 RETURN
00568              ENDIF
00569                NJET = NJET + 1
00570                DO 130 N = 1,NTRAK
00571                  JETLIS(NJET,N) = NEWLIS(N)
00572 130            CONTINUE
00573                DO 140 MU=1,4
00574                  PJ(MU,NJET)=PNEW(MU)
00575 140          CONTINUE
00576              ENDIF
00577              RETURN
00578        ENDIF
00579 *** The jet was not stable, so we iterate again
00580        DO 150 N=1,NTRAK
00581          OLDLIS(N)=NEWLIS(N)
00582 150    CONTINUE
00583        DO 160 MU=1,3
00584          OAXIS(MU)=NAXIS(MU)
00585 160    CONTINUE
00586 120   CONTINUE
00587       UNSTBL = .TRUE.
00588       RETURN
00589       END
00590 *CMZ :  1.06/00 28/02/94  15.44.44  by  P. Schleper
00591 *-- Author :
00592 C-----------------------------------------------------------------------
00593       SUBROUTINE PXSORV(N,A,K,OPT)
00594 C     Sort A(N) into ascending order
00595 C     OPT = 'I' : return index array K only
00596 C     OTHERWISE : return sorted A and index array K
00597 C-----------------------------------------------------------------------
00598       INTEGER NMAX
00599       PARAMETER (NMAX=5000)
00600 *
00601 *      INTEGER N,I,J,K(N),IL(NMAX),IR(NMAX)
00602 *LUND
00603       INTEGER N,I,J,K(NMAX),IL(NMAX),IR(NMAX)
00604       CHARACTER OPT
00605 *
00606 *      DOUBLE PRECISION A(N),B(NMAX)
00607       DOUBLE PRECISION A(NMAX),B(NMAX)
00608 *LUND
00609       IF (N.GT.NMAX) STOP 'Sorry, not enough room in Mike''s PXSORV'
00610       IL(1)=0
00611       IR(1)=0
00612       DO 10 I=2,N
00613       IL(I)=0
00614       IR(I)=0
00615       J=1
00616    2  IF(A(I).GT.A(J)) GO TO 5
00617    3  IF(IL(J).EQ.0) GO TO 4
00618       J=IL(J)
00619       GO TO 2
00620    4  IR(I)=-J
00621       IL(J)=I
00622       GO TO 10
00623    5  IF(IR(J).LE.0) GO TO 6
00624       J=IR(J)
00625       GO TO 2
00626    6  IR(I)=IR(J)
00627       IR(J)=I
00628   10  CONTINUE
00629       I=1
00630       J=1
00631       GO TO 8
00632   20  J=IL(J)
00633    8  IF(IL(J).GT.0) GO TO 20
00634    9  K(I)=J
00635       B(I)=A(J)
00636       I=I+1
00637       IF(IR(J)) 12,30,13
00638   13  J=IR(J)
00639       GO TO 8
00640   12  J=-IR(J)
00641       GO TO 9
00642   30  IF(OPT.EQ.'I') RETURN
00643       DO 31 I=1,N
00644   31  A(I)=B(I)
00645  999  END
00646 
00647 *********************************************************************
00648 *CMZ :  2.00/00 10/01/95  10.17.57  by  P. Schleper
00649 *CMZ :  1.06/00 14/03/94  15.37.44  by  P. Schleper
00650 *-- Author :
00651 *
00652 C+DECK,PXTRY.
00653        SUBROUTINE PXTRY(MODE,COSR,NTRAK,PU,PP,OAXIS,NAXIS,
00654      +                  PNEW,NEWLIS,OK)
00655 *
00656 C+SEQ,DECLARE.
00657       INTEGER MXTRAK
00658       PARAMETER (MXTRAK=5000)
00659        INTEGER NTRAK,MODE
00660 *** Note that although PU and PP are assumed to be 2d arrays, they
00661 *** are used as 1d in this routine for efficiency
00662        DOUBLE PRECISION COSR,PU(3*MXTRAK),PP(4*MXTRAK),OAXIS(3),PXMDPI
00663        LOGICAL OK
00664        LOGICAL NEWLIS(MXTRAK)
00665        DOUBLE PRECISION NAXIS(3),PNEW(4)
00666 *** Finds all particles in cone of size COSR about OAXIS direction.
00667 *** Calculates 4-momentum sum of all particles in cone (PNEW) , and
00668 *** returns this as new jet axis NAXIS (Both unit Vectors)
00669        INTEGER N,MU,NPU,NPP
00670        DOUBLE PRECISION COSVAL,NORMSQ,NORM
00671 *
00672        OK = .FALSE.
00673        DO 100 MU=1,4
00674           PNEW(MU)=0.0
00675 100    CONTINUE
00676        NPU=-3
00677        NPP=-4
00678        DO 110 N=1,NTRAK
00679           NPU=NPU+3
00680           NPP=NPP+4
00681           IF (MODE.NE.2) THEN
00682              COSVAL=0.0
00683              DO 120 MU=1,3
00684                 COSVAL=COSVAL+OAXIS(MU)*PU(MU+NPU)
00685 120          CONTINUE
00686           ELSE
00687              IF (ABS(PU(1+NPU)).GE.20.OR.ABS(OAXIS(1)).GE.20) THEN
00688                 COSVAL=-1000
00689              ELSE
00690                 COSVAL=1-
00691      +           ((OAXIS(1)-PU(1+NPU))**2+PXMDPI(OAXIS(2)-PU(2+NPU))**2)
00692              ENDIF
00693           ENDIF
00694           IF (COSVAL.GE.COSR)THEN
00695              NEWLIS(N) = .TRUE.
00696              OK = .TRUE.
00697              IF (MODE.NE.2) THEN
00698                 DO 130 MU=1,4
00699                    PNEW(MU) = PNEW(MU) + PP(MU+NPP)
00700 130             CONTINUE
00701              ELSE
00702                 PNEW(1)=PNEW(1)
00703      +              + PP(4+NPP)/(PP(4+NPP)+PNEW(4))*(PP(1+NPP)-PNEW(1))
00704 c                PNEW(2)=PNEW(2)
00705 c     +              + PP(4+NPP)/(PP(4+NPP)+PNEW(4))
00706 c     +               *PXMDPI(PP(2+NPP)-PNEW(2))
00707 ! GPS 25/02/07
00708                 PNEW(2)=PXMDPI(PNEW(2)
00709      +              + PP(4+NPP)/(PP(4+NPP)+PNEW(4))
00710      +               *PXMDPI(PP(2+NPP)-PNEW(2)))
00711                 PNEW(4)=PNEW(4)+PP(4+NPP)
00712              ENDIF
00713           ELSE
00714              NEWLIS(N)=.FALSE.
00715           ENDIF
00716 110   CONTINUE
00717 *** If there are particles in the cone, calc new jet axis
00718        IF (OK) THEN
00719           IF (MODE.NE.2) THEN
00720              NORMSQ = 0.0
00721              DO 140 MU = 1,3
00722                 NORMSQ = NORMSQ + PNEW(MU)**2
00723 140          CONTINUE
00724              NORM = SQRT(NORMSQ)
00725           ELSE
00726              NORM = 1
00727           ENDIF
00728           DO 150 MU=1,3
00729              NAXIS(MU) = PNEW(MU)/NORM
00730 150       CONTINUE
00731        ENDIF
00732        RETURN
00733        END
00734 
00735 *********************************************************************
00736 *CMZ :  2.00/00 10/01/95  10.17.57  by  P. Schleper
00737 *CMZ :  1.06/00 28/02/94  15.44.44  by  P. Schleper
00738 *-- Author :
00739 C+DECK,PXUVEC.
00740 *
00741        SUBROUTINE PXUVEC(NTRAK,PP,PU,IERR)
00742 *
00743 *** Routine to calculate unit vectors PU of all particles PP
00744 C+SEQ,DECLARE.
00745       INTEGER MXTRAK
00746       PARAMETER (MXTRAK=5000)
00747       INTEGER NTRAK, IERR
00748       DOUBLE PRECISION PP(4,MXTRAK)
00749       DOUBLE PRECISION PU(3,MXTRAK)
00750       INTEGER N,MU
00751       DOUBLE PRECISION MAG
00752        DO 100 N=1,NTRAK
00753           MAG=0.0
00754           DO 110 MU=1,3
00755              MAG=MAG+PP(MU,N)**2
00756 110       CONTINUE
00757           MAG=SQRT(MAG)
00758           IF (MAG.EQ.0.0) THEN
00759              WRITE(6,*)' PXCONE: An input particle has zero mod(p)'
00760              IERR=-1
00761              RETURN
00762           ENDIF
00763           DO 120 MU=1,3
00764            PU(MU,N)=PP(MU,N)/MAG
00765 120       CONTINUE
00766 100    CONTINUE
00767        RETURN
00768        END
00769 *CMZ :  1.06/00 28/02/94  15.44.44  by  P. Schleper
00770 *-- Author :
00771 C-----------------------------------------------------------------------
00772       SUBROUTINE PXZERI(N,A)
00773       INTEGER I,N,A(N)
00774       DO 10 I=1,N
00775         A(I)=0
00776  10   CONTINUE
00777       END
00778 *CMZ :  1.06/00 28/02/94  15.44.44  by  P. Schleper
00779 *-- Author :
00780 C-----------------------------------------------------------------------
00781 C     This is a set of routines written by Mike Seymour to provide the
00782 C     services presumably normally provided by standard OPAL routines
00783 C     PXZERV zeroes a vector
00784 C     PXZERI zeroes a vector of integers
00785 C     PXNORV normalizes a vector
00786 C     PXADDV adds two vectors
00787 C     PXSORV sorts a vector (copied from HERWIG)
00788 C     PXANG3 finds the angle (and its cosine) between two vectors
00789 C     PXMDPI moves its argument onto the range [-pi,pi)
00790 C-----------------------------------------------------------------------
00791       SUBROUTINE PXZERV(N,A)
00792       INTEGER I,N
00793       DOUBLE PRECISION A(N)
00794       DO 10 I=1,N
00795         A(I)=0
00796  10   CONTINUE
00797       END
00798 *-- Author :
00799 C-----------------------------------------------------------------------
00800       SUBROUTINE PXADDV(N,A,B,C,ITERR)
00801       INTEGER I,N,ITERR
00802       DOUBLE PRECISION A(N),B(N),C(N)
00803       DO 10 I=1,N
00804         C(I)=A(I)+B(I)
00805  10   CONTINUE
00806       END
00807 *CMZ :  1.06/00 28/02/94  15.44.44  by  P. Schleper
00808 *-- Author :
00809 C-----------------------------------------------------------------------
00810       SUBROUTINE PXANG3(A,B,COST,THET,ITERR)
00811       INTEGER ITERR
00812       DOUBLE PRECISION A(3),B(3),C,COST,THET
00813       C=(A(1)**2+A(2)**2+A(3)**2)*(B(1)**2+B(2)**2+B(3)**2)
00814       IF (C.LE.0) RETURN
00815       C=1/SQRT(C)
00816       COST=(A(1)*B(1)+A(2)*B(2)+A(3)*B(3))*C
00817       THET=ACOS(COST)
00818       END
00819 *CMZ :  1.06/00 14/03/94  15.41.57  by  P. Schleper
00820 *-- Author :    P. Schleper   28/02/94
00821        LOGICAL FUNCTION PXNEW(TSTLIS,JETLIS,NTRAK,NJET)
00822 *
00823       INTEGER MXTRAK,MXPROT
00824       PARAMETER (MXTRAK=5000,MXPROT=5000)
00825        INTEGER NTRAK,NJET
00826 *** Note that although JETLIS is assumed to be a 2d array, it
00827 *** it is used as 1d in this routine for efficiency
00828        LOGICAL TSTLIS(MXTRAK),JETLIS(MXPROT*MXTRAK)
00829 *** Checks to see if TSTLIS entries correspond to a jet already found
00830 *** and entered in JETLIS
00831        INTEGER N, I, IN
00832        LOGICAL MATCH
00833 *
00834        PXNEW = .TRUE.
00835        DO 100 I = 1,NJET
00836           MATCH = .TRUE.
00837           IN=I-MXPROT
00838           DO 110 N = 1,NTRAK
00839             IN=IN+MXPROT
00840             IF(TSTLIS(N).NEQV.JETLIS(IN)) THEN
00841              MATCH = .FALSE.
00842              GO TO 100
00843             ENDIF
00844 110       CONTINUE
00845           IF (MATCH) THEN
00846            PXNEW = .FALSE.
00847            RETURN
00848           ENDIF
00849 100    CONTINUE
00850        RETURN
00851        END
00852 *CMZ :  1.06/00 14/03/94  15.41.57  by  P. Schleper
00853 *-- Author :    P. Schleper   28/02/94
00854        LOGICAL FUNCTION PXSAME(LIST1,LIST2,N)
00855 *
00856        LOGICAL LIST1(*),LIST2(*)
00857        INTEGER N
00858 *** Returns T if the first N elements of LIST1 are the same as the
00859 *** first N elements of LIST2.
00860        INTEGER I
00861 *
00862        PXSAME = .TRUE.
00863        DO 100 I = 1,N
00864         IF ( LIST1(I).NEQV.LIST2(I) ) THEN
00865           PXSAME = .FALSE.
00866           RETURN
00867         ENDIF
00868 100    CONTINUE
00869        RETURN
00870        END
00871 *CMZ :  1.06/00 28/02/94  15.44.44  by  P. Schleper
00872 *-- Author :
00873 C-----------------------------------------------------------------------
00874       FUNCTION PXMDPI(PHI)
00875       IMPLICIT NONE
00876 C---RETURNS PHI, MOVED ONTO THE RANGE [-PI,PI)
00877       DOUBLE PRECISION PXMDPI,PHI,PI,TWOPI,THRPI,EPS
00878       PARAMETER (PI=3.141592654,TWOPI=6.283185307,THRPI=9.424777961)
00879       PARAMETER (EPS=1E-15)
00880       PXMDPI=PHI
00881       IF (PXMDPI.LE.PI) THEN
00882         IF (PXMDPI.GT.-PI) THEN
00883           GOTO 100
00884         ELSEIF (PXMDPI.GT.-THRPI) THEN
00885           PXMDPI=PXMDPI+TWOPI
00886         ELSE
00887           PXMDPI=-MOD(PI-PXMDPI,TWOPI)+PI
00888         ENDIF
00889       ELSEIF (PXMDPI.LE.THRPI) THEN
00890         PXMDPI=PXMDPI-TWOPI
00891       ELSE
00892         PXMDPI=MOD(PI+PXMDPI,TWOPI)-PI
00893       ENDIF
00894  100  IF (ABS(PXMDPI).LT.EPS) PXMDPI=0
00895       END
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends