00001 SUBROUTINE PXCONE(MODE,NTRAK,ITKDM,PTRAK,CONER,EPSLON,OVLIM,
00002 + MXJET,NJET,PJET,IPASS,IJMUL,IERR)
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028
00029
00030
00031
00032
00033
00034
00035
00036
00037
00038
00039
00040
00041
00042
00043
00044
00045
00046
00047
00048
00049
00050
00051
00052
00053
00054
00055
00056
00057
00058
00059
00060
00061
00062
00063
00064
00065
00066
00067
00068
00069
00070
00071
00072
00073
00074
00075
00076
00077
00078
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
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
00089 DOUBLE PRECISION COSR,COS2R, VSEED(3), VEC1(3), VEC2(3),PTSQ,PPSQ,
00090 + COSVAL,PXMDPI
00091
00092 DOUBLE PRECISION RSEP
00093
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
00103
00104 RSEP = 2D0
00105
00106
00107 IERR=0
00108
00109
00110 IF(NCALL.LE.0) THEN
00111 ROLD = 0.
00112 EPSOLD = 0.
00113 OVOLD = 0.
00114 ENDIF
00115 NCALL = NCALL + 1
00116
00117
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
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
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
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
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
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
00213 COSR = 1-CONER**2
00214
00215 COS2R = 1-(RSEP*CONER)**2
00216
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
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
00235
00236
00237
00238
00239
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
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
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
00276
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
00288
00289 CALL PXORD(EPSLON,NJET,NTRAK,JETLIS,PJ)
00290
00291
00292 CALL PXOLAP(MODE,NJET,NTRAK,JETLIS,PJ,PP,OVLIM)
00293
00294
00295 CALL PXORD(EPSLON,NJET,NTRAK,JETLIS,PJ)
00296
00297
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
00329
00330
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
00345
00346
00347
00348
00349 SUBROUTINE PXOLAP(MODE,NJET,NTRAK,JETLIS,PJ,PP,OVLIM)
00350
00351
00352
00353
00354
00355
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
00371 DO 100 I = 2,NJET
00372
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
00386 IF (EOVER.GT.OVLIM*PJ(4,I)) THEN
00387
00388 DO 130 N = 1,NTRAK
00389 JETLIS(I,N) = .FALSE.
00390 130 CONTINUE
00391 ENDIF
00392 100 CONTINUE
00393
00394
00395
00396
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
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
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
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
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
00453
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
00462
00463
00464
00465
00466 SUBROUTINE PXORD(EPSLON,NJET,NTRAK,JETLIS,PJ)
00467
00468
00469
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
00479
00480
00481
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
00494 CALL PXSORV(NJET,ELIST,INDEX,'I')
00495
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
00505
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
00517
00518
00519
00520 SUBROUTINE PXSEAR(MODE,COSR,NTRAK,PU,PP,VSEED,NJET,
00521 + JETLIS,PJ,UNSTBL,IERR)
00522
00523
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
00533
00534
00535
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
00552 IF (.NOT.OK) THEN
00553 RETURN
00554 ENDIF
00555 IF(PXSAME(NEWLIS,OLDLIS,NTRAK)) THEN
00556
00557 IF (PXNEW(NEWLIS,JETLIS,NTRAK,NJET)) THEN
00558
00559
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
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
00587
00588
00589 SUBROUTINE PXSORV(N,A,K,OPT)
00590
00591
00592
00593
00594 INTEGER NMAX
00595 PARAMETER (NMAX=5000)
00596
00597
00598
00599 INTEGER N,I,J,K(NMAX),IL(NMAX),IR(NMAX)
00600 CHARACTER OPT
00601
00602
00603 DOUBLE PRECISION A(NMAX),B(NMAX)
00604
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
00645
00646
00647
00648
00649 SUBROUTINE PXTRY(MODE,COSR,NTRAK,PU,PP,OAXIS,NAXIS,
00650 + PNEW,NEWLIS,OK)
00651
00652
00653 INTEGER MXTRAK
00654 PARAMETER (MXTRAK=5000)
00655 INTEGER NTRAK,MODE
00656
00657
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
00663
00664
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
00701
00702
00703
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
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
00733
00734
00735
00736
00737 SUBROUTINE PXUVEC(NTRAK,PP,PU,IERR)
00738
00739
00740
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
00766
00767
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
00775
00776
00777
00778
00779
00780
00781
00782
00783
00784
00785
00786
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
00795
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
00804
00805
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
00816
00817 LOGICAL FUNCTION PXNEW(TSTLIS,JETLIS,NTRAK,NJET)
00818
00819 INTEGER MXTRAK,MXPROT
00820 PARAMETER (MXTRAK=5000,MXPROT=5000)
00821 INTEGER NTRAK,NJET
00822
00823
00824 LOGICAL TSTLIS(MXTRAK),JETLIS(MXPROT*MXTRAK)
00825
00826
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
00849
00850 LOGICAL FUNCTION PXSAME(LIST1,LIST2,N)
00851
00852 LOGICAL LIST1(*),LIST2(*)
00853 INTEGER N
00854
00855
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
00868
00869
00870 FUNCTION PXMDPI(PHI)
00871 IMPLICIT NONE
00872
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