|
fastjet 2.4.5
|
|
fastjet 2.4.5
|
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 ! GPS -- to get consistent behaviour (2-way midpnts) 00239 c DO 140 N1 = 1,NJTORG-1 ! GPS -- to get consistent behaviour (2-way midpnts) 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 ! GPS -- to get consistent behaviour 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 !VSEED(2)=VSEED(2)/2 00257 ! GPS 25/02/07 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!, Copy output into output arrays 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's energy is contained in 00353 *** higher energy jets, that jet is neglected. 00354 *** Particles assigned to the jet closest in angle (a la CDF, Snowmass). 00355 C+SEQ,DECLARE. 00356 INTEGER MXTRAK, MXPROT 00357 PARAMETER (MXTRAK=5000,MXPROT=5000) 00358 INTEGER NJET, NTRAK, MODE 00359 LOGICAL JETLIS(MXPROT,MXTRAK) 00360 DOUBLE PRECISION PJ(4,MXPROT),PP(4,MXTRAK),PXMDPI 00361 INTEGER I,J,N,MU 00362 LOGICAL OVELAP 00363 DOUBLE PRECISION EOVER 00364 DOUBLE PRECISION OVLIM 00365 INTEGER ITERR, IJMIN, IJET(MXPROT), NJ 00366 DOUBLE PRECISION VEC1(3), VEC2(3), COST, THET, THMIN 00367 DATA IJMIN/0/ 00368 * 00369 IF (NJET.LE.1) RETURN 00370 *** Look for jets with large overlaps with higher energy jets. 00371 DO 100 I = 2,NJET 00372 *** Find overlap energy between jets I and all higher energy jets. 00373 EOVER = 0.0 00374 DO 110 N = 1,NTRAK 00375 OVELAP = .FALSE. 00376 DO 120 J= 1,I-1 00377 IF (JETLIS(I,N).AND.JETLIS(J,N)) THEN 00378 OVELAP = .TRUE. 00379 ENDIF 00380 120 CONTINUE 00381 IF (OVELAP) THEN 00382 EOVER = EOVER + PP(4,N) 00383 ENDIF 00384 110 CONTINUE 00385 *** Is the fraction of energy shared larger than OVLIM? 00386 IF (EOVER.GT.OVLIM*PJ(4,I)) THEN 00387 *** De-assign all particles from Jet I 00388 DO 130 N = 1,NTRAK 00389 JETLIS(I,N) = .FALSE. 00390 130 CONTINUE 00391 ENDIF 00392 100 CONTINUE 00393 *** Now there are no big overlaps, assign every particle in 00394 *** more than 1 jet to the closet jet. 00395 *** Any particles now in more than 1 jet are assigned to the CLOSET 00396 *** jet (in angle). 00397 DO 140 I=1,NTRAK 00398 NJ=0 00399 DO 150 J=1, NJET 00400 IF(JETLIS(J,I)) THEN 00401 NJ=NJ+1 00402 IJET(NJ)=J 00403 ENDIF 00404 150 CONTINUE 00405 IF (NJ .GT. 1) THEN 00406 *** Particle in > 1 jet - calc angles... 00407 VEC1(1)=PP(1,I) 00408 VEC1(2)=PP(2,I) 00409 VEC1(3)=PP(3,I) 00410 THMIN=0. 00411 DO 160 J=1,NJ 00412 VEC2(1)=PJ(1,IJET(J)) 00413 VEC2(2)=PJ(2,IJET(J)) 00414 VEC2(3)=PJ(3,IJET(J)) 00415 IF (MODE.NE.2) THEN 00416 CALL PXANG3(VEC1,VEC2,COST,THET,ITERR) 00417 ELSE 00418 THET=(VEC1(1)-VEC2(1))**2+PXMDPI(VEC1(2)-VEC2(2))**2 00419 ENDIF 00420 IF (J .EQ. 1) THEN 00421 THMIN=THET 00422 IJMIN=IJET(J) 00423 ELSEIF (THET .LT. THMIN) THEN 00424 THMIN=THET 00425 IJMIN=IJET(J) 00426 ENDIF 00427 160 CONTINUE 00428 *** Assign track to IJMIN 00429 DO 170 J=1,NJET 00430 JETLIS(J,I) = .FALSE. 00431 170 CONTINUE 00432 JETLIS(IJMIN,I)=.TRUE. 00433 ENDIF 00434 140 CONTINUE 00435 *** Recompute PJ 00436 DO 200 I = 1,NJET 00437 DO 210 MU = 1,4 00438 PJ(MU,I) = 0.0 00439 210 CONTINUE 00440 DO 220 N = 1,NTRAK 00441 IF( JETLIS(I,N) ) THEN 00442 IF (MODE.NE.2) THEN 00443 DO 230 MU = 1,4 00444 PJ(MU,I) = PJ(MU,I) + PP(MU,N) 00445 230 CONTINUE 00446 ELSE 00447 PJ(1,I)=PJ(1,I) 00448 + + PP(4,N)/(PP(4,N)+PJ(4,I))*(PP(1,N)-PJ(1,I)) 00449 c GPS 25/02/07 00450 PJ(2,I)=PXMDPI(PJ(2,I) 00451 + + PP(4,N)/(PP(4,N)+PJ(4,I))*PXMDPI(PP(2,N)-PJ(2,I))) 00452 c PJ(2,I)=PJ(2,I) 00453 c + + PP(4,N)/(PP(4,N)+PJ(4,I))*PXMDPI(PP(2,N)-PJ(2,I)) 00454 PJ(4,I)=PJ(4,I)+PP(4,N) 00455 ENDIF 00456 ENDIF 00457 220 CONTINUE 00458 200 CONTINUE 00459 RETURN 00460 END 00461 *CMZ : 2.00/00 10/01/95 10.17.57 by P. Schleper 00462 *CMZ : 1.06/00 14/03/94 15.37.45 by P. Schleper 00463 *-- Author : 00464 * 00465 C+DECK,PXORD. 00466 SUBROUTINE PXORD(EPSLON,NJET,NTRAK,JETLIS,PJ) 00467 * 00468 *** Routine to put jets into order and eliminate tose less than EPSLON 00469 C+SEQ,DECLARE. 00470 INTEGER MXTRAK,MXPROT 00471 PARAMETER (MXTRAK=5000,MXPROT=5000) 00472 INTEGER I, J, INDEX(MXPROT) 00473 DOUBLE PRECISION PTEMP(4,MXPROT), ELIST(MXPROT) 00474 INTEGER NJET,NTRAK 00475 LOGICAL JETLIS(MXPROT,MXTRAK) 00476 LOGICAL LOGTMP(MXPROT,MXTRAK) 00477 DOUBLE PRECISION EPSLON,PJ(4,MXPROT) 00478 *** Puts jets in order of energy: 1 = highest energy etc. 00479 *** Then Eliminate jets with energy below EPSLON 00480 * 00481 *** Copy input arrays. 00482 DO 100 I=1,NJET 00483 DO 110 J=1,4 00484 PTEMP(J,I)=PJ(J,I) 00485 110 CONTINUE 00486 DO 120 J=1,NTRAK 00487 LOGTMP(I,J)=JETLIS(I,J) 00488 120 CONTINUE 00489 100 CONTINUE 00490 DO 150 I=1,NJET 00491 ELIST(I)=PJ(4,I) 00492 150 CONTINUE 00493 *** Sort the energies... 00494 CALL PXSORV(NJET,ELIST,INDEX,'I') 00495 *** Fill PJ and JETLIS according to sort ( sort is in ascending order!!) 00496 DO 200 I=1, NJET 00497 DO 210 J=1,4 00498 PJ(J,I)=PTEMP(J,INDEX(NJET+1-I)) 00499 210 CONTINUE 00500 DO 220 J=1,NTRAK 00501 JETLIS(I,J)=LOGTMP(INDEX(NJET+1-I),J) 00502 220 CONTINUE 00503 200 CONTINUE 00504 ** Jets are now in order 00505 *** Now eliminate jets with less than Epsilon energy 00506 DO 300, I=1, NJET 00507 IF (PJ(4,I) .LT. EPSLON) THEN 00508 NJET=NJET-1 00509 PJ(4,I)=0. 00510 ENDIF 00511 300 CONTINUE 00512 RETURN 00513 END 00514 00515 ******************************************************************** 00516 *CMZ : 2.00/00 10/01/95 10.17.57 by P. Schleper 00517 *CMZ : 1.06/00 14/03/94 15.37.44 by P. Schleper 00518 *-- Author : 00519 C+DECK,PXSEAR. 00520 SUBROUTINE PXSEAR(MODE,COSR,NTRAK,PU,PP,VSEED,NJET, 00521 + JETLIS,PJ,UNSTBL,IERR) 00522 * 00523 C+SEQ,DECLARE. 00524 INTEGER MXTRAK, MXPROT 00525 PARAMETER (MXTRAK=5000,MXPROT=5000) 00526 INTEGER NTRAK, IERR, MODE 00527 DOUBLE PRECISION COSR,PU(3,MXTRAK),PP(4,MXTRAK),VSEED(3) 00528 LOGICAL UNSTBL 00529 LOGICAL JETLIS(MXPROT,MXTRAK) 00530 INTEGER NJET 00531 DOUBLE PRECISION PJ(4,MXPROT) 00532 *** Using VSEED as a trial axis , look for a stable jet. 00533 *** Check stable jets against those already found and add to PJ. 00534 *** Will try up to MXITER iterations to get a stable set of particles 00535 *** in the cone. 00536 INTEGER MU,N,ITER 00537 LOGICAL PXSAME,PXNEW,OK 00538 LOGICAL NEWLIS(MXTRAK),OLDLIS(MXTRAK) 00539 DOUBLE PRECISION OAXIS(3),NAXIS(3),PNEW(4) 00540 INTEGER MXITER 00541 PARAMETER(MXITER = 30) 00542 * 00543 DO 100 MU=1,3 00544 OAXIS(MU) = VSEED(MU) 00545 100 CONTINUE 00546 DO 110 N = 1,NTRAK 00547 OLDLIS(N) = .FALSE. 00548 110 CONTINUE 00549 DO 120 ITER = 1,MXITER 00550 CALL PXTRY(MODE,COSR,NTRAK,PU,PP,OAXIS,NAXIS,PNEW,NEWLIS,OK) 00551 *** Return immediately if there were no particles in the cone. 00552 IF (.NOT.OK) THEN 00553 RETURN 00554 ENDIF 00555 IF(PXSAME(NEWLIS,OLDLIS,NTRAK)) THEN 00556 *** We have a stable jet. 00557 IF (PXNEW(NEWLIS,JETLIS,NTRAK,NJET)) THEN 00558 *** And the jet is a new one. So add it to our arrays. 00559 *** Check arrays are big anough... 00560 IF (NJET .EQ. MXPROT) THEN 00561 WRITE (6,*) ' PXCONE: Found more than MXPROT proto-jets' 00562 IERR = -1 00563 RETURN 00564 ENDIF 00565 NJET = NJET + 1 00566 DO 130 N = 1,NTRAK 00567 JETLIS(NJET,N) = NEWLIS(N) 00568 130 CONTINUE 00569 DO 140 MU=1,4 00570 PJ(MU,NJET)=PNEW(MU) 00571 140 CONTINUE 00572 ENDIF 00573 RETURN 00574 ENDIF 00575 *** The jet was not stable, so we iterate again 00576 DO 150 N=1,NTRAK 00577 OLDLIS(N)=NEWLIS(N) 00578 150 CONTINUE 00579 DO 160 MU=1,3 00580 OAXIS(MU)=NAXIS(MU) 00581 160 CONTINUE 00582 120 CONTINUE 00583 UNSTBL = .TRUE. 00584 RETURN 00585 END 00586 *CMZ : 1.06/00 28/02/94 15.44.44 by P. Schleper 00587 *-- Author : 00588 C----------------------------------------------------------------------- 00589 SUBROUTINE PXSORV(N,A,K,OPT) 00590 C Sort A(N) into ascending order 00591 C OPT = 'I' : return index array K only 00592 C OTHERWISE : return sorted A and index array K 00593 C----------------------------------------------------------------------- 00594 INTEGER NMAX 00595 PARAMETER (NMAX=5000) 00596 * 00597 * INTEGER N,I,J,K(N),IL(NMAX),IR(NMAX) 00598 *LUND 00599 INTEGER N,I,J,K(NMAX),IL(NMAX),IR(NMAX) 00600 CHARACTER OPT 00601 * 00602 * DOUBLE PRECISION A(N),B(NMAX) 00603 DOUBLE PRECISION A(NMAX),B(NMAX) 00604 *LUND 00605 IF (N.GT.NMAX) STOP 'Sorry, not enough room in Mike''s PXSORV' 00606 IL(1)=0 00607 IR(1)=0 00608 DO 10 I=2,N 00609 IL(I)=0 00610 IR(I)=0 00611 J=1 00612 2 IF(A(I).GT.A(J)) GO TO 5 00613 3 IF(IL(J).EQ.0) GO TO 4 00614 J=IL(J) 00615 GO TO 2 00616 4 IR(I)=-J 00617 IL(J)=I 00618 GO TO 10 00619 5 IF(IR(J).LE.0) GO TO 6 00620 J=IR(J) 00621 GO TO 2 00622 6 IR(I)=IR(J) 00623 IR(J)=I 00624 10 CONTINUE 00625 I=1 00626 J=1 00627 GO TO 8 00628 20 J=IL(J) 00629 8 IF(IL(J).GT.0) GO TO 20 00630 9 K(I)=J 00631 B(I)=A(J) 00632 I=I+1 00633 IF(IR(J)) 12,30,13 00634 13 J=IR(J) 00635 GO TO 8 00636 12 J=-IR(J) 00637 GO TO 9 00638 30 IF(OPT.EQ.'I') RETURN 00639 DO 31 I=1,N 00640 31 A(I)=B(I) 00641 999 END 00642 00643 ********************************************************************* 00644 *CMZ : 2.00/00 10/01/95 10.17.57 by P. Schleper 00645 *CMZ : 1.06/00 14/03/94 15.37.44 by P. Schleper 00646 *-- Author : 00647 * 00648 C+DECK,PXTRY. 00649 SUBROUTINE PXTRY(MODE,COSR,NTRAK,PU,PP,OAXIS,NAXIS, 00650 + PNEW,NEWLIS,OK) 00651 * 00652 C+SEQ,DECLARE. 00653 INTEGER MXTRAK 00654 PARAMETER (MXTRAK=5000) 00655 INTEGER NTRAK,MODE 00656 *** Note that although PU and PP are assumed to be 2d arrays, they 00657 *** are used as 1d in this routine for efficiency 00658 DOUBLE PRECISION COSR,PU(3*MXTRAK),PP(4*MXTRAK),OAXIS(3),PXMDPI 00659 LOGICAL OK 00660 LOGICAL NEWLIS(MXTRAK) 00661 DOUBLE PRECISION NAXIS(3),PNEW(4) 00662 *** Finds all particles in cone of size COSR about OAXIS direction. 00663 *** Calculates 4-momentum sum of all particles in cone (PNEW) , and 00664 *** returns this as new jet axis NAXIS (Both unit Vectors) 00665 INTEGER N,MU,NPU,NPP 00666 DOUBLE PRECISION COSVAL,NORMSQ,NORM 00667 * 00668 OK = .FALSE. 00669 DO 100 MU=1,4 00670 PNEW(MU)=0.0 00671 100 CONTINUE 00672 NPU=-3 00673 NPP=-4 00674 DO 110 N=1,NTRAK 00675 NPU=NPU+3 00676 NPP=NPP+4 00677 IF (MODE.NE.2) THEN 00678 COSVAL=0.0 00679 DO 120 MU=1,3 00680 COSVAL=COSVAL+OAXIS(MU)*PU(MU+NPU) 00681 120 CONTINUE 00682 ELSE 00683 IF (ABS(PU(1+NPU)).GE.20.OR.ABS(OAXIS(1)).GE.20) THEN 00684 COSVAL=-1000 00685 ELSE 00686 COSVAL=1- 00687 + ((OAXIS(1)-PU(1+NPU))**2+PXMDPI(OAXIS(2)-PU(2+NPU))**2) 00688 ENDIF 00689 ENDIF 00690 IF (COSVAL.GE.COSR)THEN 00691 NEWLIS(N) = .TRUE. 00692 OK = .TRUE. 00693 IF (MODE.NE.2) THEN 00694 DO 130 MU=1,4 00695 PNEW(MU) = PNEW(MU) + PP(MU+NPP) 00696 130 CONTINUE 00697 ELSE 00698 PNEW(1)=PNEW(1) 00699 + + PP(4+NPP)/(PP(4+NPP)+PNEW(4))*(PP(1+NPP)-PNEW(1)) 00700 c PNEW(2)=PNEW(2) 00701 c + + PP(4+NPP)/(PP(4+NPP)+PNEW(4)) 00702 c + *PXMDPI(PP(2+NPP)-PNEW(2)) 00703 ! GPS 25/02/07 00704 PNEW(2)=PXMDPI(PNEW(2) 00705 + + PP(4+NPP)/(PP(4+NPP)+PNEW(4)) 00706 + *PXMDPI(PP(2+NPP)-PNEW(2))) 00707 PNEW(4)=PNEW(4)+PP(4+NPP) 00708 ENDIF 00709 ELSE 00710 NEWLIS(N)=.FALSE. 00711 ENDIF 00712 110 CONTINUE 00713 *** If there are particles in the cone, calc new jet axis 00714 IF (OK) THEN 00715 IF (MODE.NE.2) THEN 00716 NORMSQ = 0.0 00717 DO 140 MU = 1,3 00718 NORMSQ = NORMSQ + PNEW(MU)**2 00719 140 CONTINUE 00720 NORM = SQRT(NORMSQ) 00721 ELSE 00722 NORM = 1 00723 ENDIF 00724 DO 150 MU=1,3 00725 NAXIS(MU) = PNEW(MU)/NORM 00726 150 CONTINUE 00727 ENDIF 00728 RETURN 00729 END 00730 00731 ********************************************************************* 00732 *CMZ : 2.00/00 10/01/95 10.17.57 by P. Schleper 00733 *CMZ : 1.06/00 28/02/94 15.44.44 by P. Schleper 00734 *-- Author : 00735 C+DECK,PXUVEC. 00736 * 00737 SUBROUTINE PXUVEC(NTRAK,PP,PU,IERR) 00738 * 00739 *** Routine to calculate unit vectors PU of all particles PP 00740 C+SEQ,DECLARE. 00741 INTEGER MXTRAK 00742 PARAMETER (MXTRAK=5000) 00743 INTEGER NTRAK, IERR 00744 DOUBLE PRECISION PP(4,MXTRAK) 00745 DOUBLE PRECISION PU(3,MXTRAK) 00746 INTEGER N,MU 00747 DOUBLE PRECISION MAG 00748 DO 100 N=1,NTRAK 00749 MAG=0.0 00750 DO 110 MU=1,3 00751 MAG=MAG+PP(MU,N)**2 00752 110 CONTINUE 00753 MAG=SQRT(MAG) 00754 IF (MAG.EQ.0.0) THEN 00755 WRITE(6,*)' PXCONE: An input particle has zero mod(p)' 00756 IERR=-1 00757 RETURN 00758 ENDIF 00759 DO 120 MU=1,3 00760 PU(MU,N)=PP(MU,N)/MAG 00761 120 CONTINUE 00762 100 CONTINUE 00763 RETURN 00764 END 00765 *CMZ : 1.06/00 28/02/94 15.44.44 by P. Schleper 00766 *-- Author : 00767 C----------------------------------------------------------------------- 00768 SUBROUTINE PXZERI(N,A) 00769 INTEGER I,N,A(N) 00770 DO 10 I=1,N 00771 A(I)=0 00772 10 CONTINUE 00773 END 00774 *CMZ : 1.06/00 28/02/94 15.44.44 by P. Schleper 00775 *-- Author : 00776 C----------------------------------------------------------------------- 00777 C This is a set of routines written by Mike Seymour to provide the 00778 C services presumably normally provided by standard OPAL routines 00779 C PXZERV zeroes a vector 00780 C PXZERI zeroes a vector of integers 00781 C PXNORV normalizes a vector 00782 C PXADDV adds two vectors 00783 C PXSORV sorts a vector (copied from HERWIG) 00784 C PXANG3 finds the angle (and its cosine) between two vectors 00785 C PXMDPI moves its argument onto the range [-pi,pi) 00786 C----------------------------------------------------------------------- 00787 SUBROUTINE PXZERV(N,A) 00788 INTEGER I,N 00789 DOUBLE PRECISION A(N) 00790 DO 10 I=1,N 00791 A(I)=0 00792 10 CONTINUE 00793 END 00794 *-- Author : 00795 C----------------------------------------------------------------------- 00796 SUBROUTINE PXADDV(N,A,B,C,ITERR) 00797 INTEGER I,N,ITERR 00798 DOUBLE PRECISION A(N),B(N),C(N) 00799 DO 10 I=1,N 00800 C(I)=A(I)+B(I) 00801 10 CONTINUE 00802 END 00803 *CMZ : 1.06/00 28/02/94 15.44.44 by P. Schleper 00804 *-- Author : 00805 C----------------------------------------------------------------------- 00806 SUBROUTINE PXANG3(A,B,COST,THET,ITERR) 00807 INTEGER ITERR 00808 DOUBLE PRECISION A(3),B(3),C,COST,THET 00809 C=(A(1)**2+A(2)**2+A(3)**2)*(B(1)**2+B(2)**2+B(3)**2) 00810 IF (C.LE.0) RETURN 00811 C=1/SQRT(C) 00812 COST=(A(1)*B(1)+A(2)*B(2)+A(3)*B(3))*C 00813 THET=ACOS(COST) 00814 END 00815 *CMZ : 1.06/00 14/03/94 15.41.57 by P. Schleper 00816 *-- Author : P. Schleper 28/02/94 00817 LOGICAL FUNCTION PXNEW(TSTLIS,JETLIS,NTRAK,NJET) 00818 * 00819 INTEGER MXTRAK,MXPROT 00820 PARAMETER (MXTRAK=5000,MXPROT=5000) 00821 INTEGER NTRAK,NJET 00822 *** Note that although JETLIS is assumed to be a 2d array, it 00823 *** it is used as 1d in this routine for efficiency 00824 LOGICAL TSTLIS(MXTRAK),JETLIS(MXPROT*MXTRAK) 00825 *** Checks to see if TSTLIS entries correspond to a jet already found 00826 *** and entered in JETLIS 00827 INTEGER N, I, IN 00828 LOGICAL MATCH 00829 * 00830 PXNEW = .TRUE. 00831 DO 100 I = 1,NJET 00832 MATCH = .TRUE. 00833 IN=I-MXPROT 00834 DO 110 N = 1,NTRAK 00835 IN=IN+MXPROT 00836 IF(TSTLIS(N).NEQV.JETLIS(IN)) THEN 00837 MATCH = .FALSE. 00838 GO TO 100 00839 ENDIF 00840 110 CONTINUE 00841 IF (MATCH) THEN 00842 PXNEW = .FALSE. 00843 RETURN 00844 ENDIF 00845 100 CONTINUE 00846 RETURN 00847 END 00848 *CMZ : 1.06/00 14/03/94 15.41.57 by P. Schleper 00849 *-- Author : P. Schleper 28/02/94 00850 LOGICAL FUNCTION PXSAME(LIST1,LIST2,N) 00851 * 00852 LOGICAL LIST1(*),LIST2(*) 00853 INTEGER N 00854 *** Returns T if the first N elements of LIST1 are the same as the 00855 *** first N elements of LIST2. 00856 INTEGER I 00857 * 00858 PXSAME = .TRUE. 00859 DO 100 I = 1,N 00860 IF ( LIST1(I).NEQV.LIST2(I) ) THEN 00861 PXSAME = .FALSE. 00862 RETURN 00863 ENDIF 00864 100 CONTINUE 00865 RETURN 00866 END 00867 *CMZ : 1.06/00 28/02/94 15.44.44 by P. Schleper 00868 *-- Author : 00869 C----------------------------------------------------------------------- 00870 FUNCTION PXMDPI(PHI) 00871 IMPLICIT NONE 00872 C---RETURNS PHI, MOVED ONTO THE RANGE [-PI,PI) 00873 DOUBLE PRECISION PXMDPI,PHI,PI,TWOPI,THRPI,EPS 00874 PARAMETER (PI=3.141592654,TWOPI=6.283185307,THRPI=9.424777961) 00875 PARAMETER (EPS=1E-15) 00876 PXMDPI=PHI 00877 IF (PXMDPI.LE.PI) THEN 00878 IF (PXMDPI.GT.-PI) THEN 00879 GOTO 100 00880 ELSEIF (PXMDPI.GT.-THRPI) THEN 00881 PXMDPI=PXMDPI+TWOPI 00882 ELSE 00883 PXMDPI=-MOD(PI-PXMDPI,TWOPI)+PI 00884 ENDIF 00885 ELSEIF (PXMDPI.LE.THRPI) THEN 00886 PXMDPI=PXMDPI-TWOPI 00887 ELSE 00888 PXMDPI=MOD(PI+PXMDPI,TWOPI)-PI 00889 ENDIF 00890 100 IF (ABS(PXMDPI).LT.EPS) PXMDPI=0 00891 END
1.7.4
