|
FastJet 3.0.4
|
|
FastJet 3.0.4
|
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
1.7.4
