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