00001
00002
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
00080
00081
00082
00083
00084
00085
00086
00087
00088
00089
00090
00091
00092 SUBROUTINE PHOINI
00093
00094
00095
00096
00097
00098
00099
00100
00101
00102
00103
00104
00105
00106
00107
00108
00109 IMPLICIT NONE
00110 INTEGER INIT,IDUM,IPHQRK,IPHEKL
00111 SAVE INIT
00112 DATA INIT/ 0/
00113
00114
00115 IF (INIT.NE.0) RETURN
00116 INIT=1
00117
00118
00119
00120
00121
00122
00123
00124
00125 CALL PHCORK(1)
00126
00127
00128
00129 IDUM= IPHQRK(1)
00130
00131
00132 IDUM= IPHEKL(2)
00133
00134
00135 CALL PHOCIN
00136
00137
00138 CALL PHOINF
00139
00140
00141
00142 CALL PHORIN
00143 RETURN
00144 END
00145 SUBROUTINE PHOCIN
00146
00147
00148
00149
00150
00151
00152
00153
00154
00155
00156
00157
00158
00159
00160
00161 IMPLICIT NONE
00162 INTEGER d_h_NMXHEP
00163 #if defined (KORALK)
00164 PARAMETER (d_h_NMXHEP=2000)
00165 #elif defined (CZTERYD)
00166 PARAMETER (d_h_NMXHEP=4000)
00167 #elif defined (DWAR)
00168 PARAMETER (d_h_NMXHEP=2000)
00169 #elif defined (DWAD)
00170 PARAMETER (d_h_NMXHEP=2000)
00171 #elif defined (DZIESD)
00172 PARAMETER (d_h_NMXHEP=10000)
00173 #else
00174 PARAMETER (d_h_NMXHEP=2000)
00175 #endif
00176 LOGICAL QEDRAD
00177 COMMON/PHOQED/QEDRAD(d_h_NMXHEP)
00178 INTEGER PHLUN
00179 COMMON/PHOLUN/PHLUN
00180 REAL*8 ALPHA,XPHCUT
00181 COMMON/PHOCOP/ALPHA,XPHCUT
00182 REAL*8 PI,TWOPI
00183 COMMON/PHPICO/PI,TWOPI
00184 INTEGER ISEED,I97,J97
00185 REAL*8 URAN,CRAN,CDRAN,CMRAN
00186 COMMON/PHSEED/ISEED(2),I97,J97,URAN(97),CRAN,CDRAN,CMRAN
00187 INTEGER PHOMES
00188 PARAMETER (PHOMES=10)
00189 INTEGER STATUS
00190 COMMON/PHOSTA/STATUS(PHOMES)
00191 LOGICAL INTERF,ISEC,ITRE,IEXP,IFTOP,IFW
00192 REAL*8 FINT,FSEC,EXPEPS
00193 COMMON /PHOKEY/ FSEC,FINT,EXPEPS,INTERF,ISEC,ITRE,IEXP,IFTOP,IFW
00194 INTEGER INIT,I
00195 SAVE INIT
00196 DATA INIT/ 0/
00197
00198
00199 IF (INIT.NE.0) RETURN
00200 INIT=1
00201
00202
00203
00204 DO 10 I=1,d_h_NMXHEP
00205 10 QEDRAD(I)=.TRUE.
00206
00207
00208 PHLUN=6
00209
00210
00211 XPHCUT=0.01 D0
00212
00213
00214
00215 ALPHA=0.00729735039D0
00216 PI=3.14159265358979324D0
00217 TWOPI=6.28318530717958648D0
00218
00219
00220 ISEED(1)=1802
00221 ISEED(2)=9373
00222
00223
00224
00225
00226 INTERF=.TRUE.
00227
00228
00229 ISEC=.TRUE.
00230
00231
00232 ITRE=.FALSE.
00233
00234 IEXP=.FALSE.
00235 IF (IEXP) THEN
00236 ISEC=.FALSE.
00237 ITRE=.FALSE.
00238 CALL PHCORK(5)
00239
00240 XPHCUT=0.000 000 1
00241 EXPEPS=1D-4
00242 ENDIF
00243
00244
00245
00246 IFTOP=.TRUE.
00247
00248
00249
00250
00251
00252
00253
00254 IF (INTERF) THEN
00255
00256
00257
00258 FINT=2.0D0
00259 ELSE
00260 FINT=1.0D0
00261 ENDIF
00262
00263
00264
00265
00266
00267
00268
00269
00270
00271
00272
00273
00274
00275
00276 IFW=.TRUE.
00277
00278 DO 20 I=1,PHOMES
00279 20 STATUS(I)=0
00280 RETURN
00281 END
00282 SUBROUTINE PHOINF
00283
00284
00285
00286
00287
00288
00289
00290
00291
00292
00293
00294
00295
00296
00297 IMPLICIT NONE
00298 INTEGER IV1,IV2,IV3
00299 INTEGER PHOVN1,PHOVN2
00300 COMMON/PHOVER/PHOVN1,PHOVN2
00301 INTEGER PHLUN
00302 COMMON/PHOLUN/PHLUN
00303 LOGICAL INTERF,ISEC,ITRE,IEXP,IFTOP,IFW
00304 REAL*8 FINT,FSEC,EXPEPS
00305 COMMON /PHOKEY/ FSEC,FINT,EXPEPS,INTERF,ISEC,ITRE,IEXP,IFTOP,IFW
00306 REAL*8 ALPHA,XPHCUT
00307 COMMON/PHOCOP/ALPHA,XPHCUT
00308
00309
00310 PHOVN1=215
00311 PHOVN2=111005
00312
00313
00314 WRITE(PHLUN,9000)
00315 WRITE(PHLUN,9020)
00316 WRITE(PHLUN,9010)
00317 WRITE(PHLUN,9030)
00318 IV1=PHOVN1/100
00319 IV2=PHOVN1-IV1*100
00320 WRITE(PHLUN,9040) IV1,IV2
00321 IV1=PHOVN2/10000
00322 IV2=(PHOVN2-IV1*10000)/100
00323 IV3=PHOVN2-IV1*10000-IV2*100
00324 WRITE(PHLUN,9050) IV1,IV2,IV3
00325 WRITE(PHLUN,9030)
00326 WRITE(PHLUN,9010)
00327 WRITE(PHLUN,9060)
00328 WRITE(PHLUN,9010)
00329 WRITE(PHLUN,9070)
00330 WRITE(PHLUN,9010)
00331 WRITE(PHLUN,9020)
00332 WRITE(PHLUN,9010)
00333 WRITE(PHLUN,9064) INTERF,ISEC,ITRE,IEXP,IFTOP,IFW,ALPHA,XPHCUT
00334 WRITE(PHLUN,9010)
00335 IF (INTERF) WRITE(PHLUN,9061)
00336 IF (ISEC) WRITE(PHLUN,9062)
00337 IF (ITRE) WRITE(PHLUN,9066)
00338 IF (IEXP) WRITE(PHLUN,9067) EXPEPS
00339 IF (IFTOP) WRITE(PHLUN,9063)
00340 IF (IFW) WRITE(PHLUN,9065)
00341 WRITE(PHLUN,9080)
00342 WRITE(PHLUN,9010)
00343 WRITE(PHLUN,9020)
00344 RETURN
00345 9000 FORMAT(1H1)
00346 9010 FORMAT(1H ,'*',T81,'*')
00347 9020 FORMAT(1H ,80('*'))
00348 9030 FORMAT(1H ,'*',26X,26('='),T81,'*')
00349 9040 FORMAT(1H ,'*',28X,'PHOTOS, Version: ',I2,'.',I2,T81,'*')
00350 9050 FORMAT(1H ,'*',28X,'Released at: ',I2,'/',I2,'/',I2,T81,'*')
00351 9060 FORMAT(1H ,'*',18X,'PHOTOS QED Corrections in Particle Decays',
00352 &T81,'*')
00353 9061 FORMAT(1H ,'*',18X,'option with interference is active ',
00354 &T81,'*')
00355 9062 FORMAT(1H ,'*',18X,'option with double photons is active ',
00356 &T81,'*')
00357 9066 FORMAT(1H ,'*',18X,'option with triple/quatric photons is active',
00358 &T81,'*')
00359 9067 FORMAT(1H ,'*',18X,'option with exponentiation is active EPSEXP=',
00360 &E10.4,T81,'*')
00361 9063 FORMAT(1H ,'*',18X,'emision in t tbar production is active ',
00362 &T81,'*')
00363 9064 FORMAT(1H ,'*',18X,'Internal input parameters:',T81,'*'
00364 &,/, 1H ,'*',T81,'*'
00365 &,/, 1H ,'*',18X,'INTERF=',L2,' ISEC=',L2,' ITRE=',L2,
00366 & ' IEXP=',L2,' IFTOP=',L2,
00367 & ' IFW=',L2,T81,'*'
00368 &,/, 1H ,'*',18X,'ALPHA_QED=',F8.5,' XPHCUT=',E8.3,T81,'*')
00369 9065 FORMAT(1H ,'*',18X,'correction wt in decay of W is active ',
00370 &T81,'*')
00371 9070 FORMAT(1H ,'*',9X,
00372 &'Monte Carlo Program - by E. Barberio, B. van Eijk and Z. Was',
00373 & T81,'*',/,
00374 & 1H ,'*',9X,'Version 2.09 - by P. Golonka and Z.W.',T81,'*')
00375 9080 FORMAT( 1H ,'*',9X,' ',T81,'*',/,
00376 & 1H ,'*',9X,
00377 & ' WARNING (1): /HEPEVT/ is not anymore the standard common block'
00378 & ,T81,'*',/,
00379 & 1H ,'*',9X,' ',T81,'*',/,
00380 & 1H ,'*',9X,
00381 #if defined (KORALK)
00382 & ' PHOTOS expects /d_HEPEVT/ to have REAL*8 variabl. To change to'
00383 & ,T81,'*',/, 1H ,'*',9X,
00384 & ' REAL*4 modify its declaration in subr. PHOTOS_GET PHOTOS_SET:'
00385 & ,T81,'*',/, 1H ,'*',9X,
00386 & ' REAL*8 d_h_phep, d_h_vhep'
00387 & ,T81,'*',/, 1H ,'*',9X,
00388 & ' WARNING (2): check dims. of /d_hepevt/ /phoqed/ /ph_hepevt/.'
00389 & ,T81,'*',/, 1H ,'*',9X,
00390 & ' HERE: d_h_nmxhep=2000 and NMXHEP=10000'
00391 #elif defined (CZTERYD)
00392 & ' PHOTOS expects /HEPEVT/ to have REAL*8 variables. To change to'
00393 & ,T81,'*',/, 1H ,'*',9X,
00394 & ' REAL*4 modify its declaration in subr. PHOTOS_GET PHOTOS_SET:'
00395 & ,T81,'*',/, 1H ,'*',9X,
00396 & ' REAL*8 d_h_phep, d_h_vhep'
00397 & ,T81,'*',/, 1H ,'*',9X,
00398 & ' WARNING (2): check dims. of /hepevt/ /phoqed/ /ph_hepevt/.'
00399 & ,T81,'*',/, 1H ,'*',9X,
00400 & ' HERE: d_h_nmxhep=4000 and NMXHEP=10000'
00401 #elif defined (DWAR)
00402 & ' PHOTOS expects /HEPEVT/ to have REAL*4 variables. To change to'
00403 & ,T81,'*',/, 1H ,'*',9X,
00404 & ' REAL*8 modify its declaration in subr. PHOTOS_GET PHOTOS_SET:'
00405 & ,T81,'*',/, 1H ,'*',9X,
00406 & ' REAL*8 d_h_phep, d_h_vhep'
00407 & ,T81,'*',/, 1H ,'*',9X,
00408 & ' WARNING (2): check dims. of /hepevt/ /phoqed/ /ph_hepevt/.'
00409 & ,T81,'*',/, 1H ,'*',9X,
00410 & ' HERE: d_h_nmxhep=2000 and NMXHEP=10000'
00411 #elif defined (DWAD)
00412 & ' PHOTOS expects /HEPEVT/ to have REAL*8 variables. To change to'
00413 & ,T81,'*',/, 1H ,'*',9X,
00414 & ' REAL*4 modify its declaration in subr. PHOTOS_GET PHOTOS_SET:'
00415 & ,T81,'*',/, 1H ,'*',9X,
00416 & ' REAL*8 d_h_phep, d_h_vhep'
00417 & ,T81,'*',/, 1H ,'*',9X,
00418 & ' WARNING (2): check dims. of /hepevt/ /phoqed/ /ph_hepevt/.'
00419 & ,T81,'*',/, 1H ,'*',9X,
00420 & ' HERE: d_h_nmxhep=2000 and NMXHEP=10000'
00421 #elif defined (DZIESD)
00422 & ' PHOTOS expects /HEPEVT/ to have REAL*8 variables. To change to'
00423 & ,T81,'*',/, 1H ,'*',9X,
00424 & ' REAL*4 modify its declaration in subr. PHOTOS_GET PHOTOS_SET:'
00425 & ,T81,'*',/, 1H ,'*',9X,
00426 & ' REAL*8 d_h_phep, d_h_vhep'
00427 & ,T81,'*',/, 1H ,'*',9X,
00428 & ' WARNING (2): check dims. of /hepevt/ /phoqed/ /ph_hepevt/.'
00429 & ,T81,'*',/, 1H ,'*',9X,
00430 & ' HERE: d_h_nmxhep=10000 and NMXHEP=10000'
00431 #else
00432 & ' PHOTOS expects /HEPEVT/ to have REAL*8 variables. To change to'
00433 & ,T81,'*',/, 1H ,'*',9X,
00434 & ' REAL*4 modify its declaration in subr. PHOTOS_GET PHOTOS_SET:'
00435 & ,T81,'*',/, 1H ,'*',9X,
00436 & ' REAL*8 d_h_phep, d_h_vhep'
00437 & ,T81,'*',/, 1H ,'*',9X,
00438 & ' WARNING (2): check dims. of /hepevt/ /phoqed/ /ph_hepevt/.'
00439 & ,T81,'*',/, 1H ,'*',9X,
00440 & ' HERE: d_h_nmxhep=2000 and NMXHEP=10000'
00441 #endif
00442 & ,T81,'*')
00443 END
00444 SUBROUTINE PHOTOS(ID)
00445
00446
00447
00448
00449
00450
00451
00452
00453
00454
00455
00456
00457
00458
00459
00460
00461
00462 COMMON /PHLUPY/ IPOIN,IPOINM
00463 INTEGER IPOIN,IPOINM
00464 COMMON /PHNUM/ IEV
00465 INTEGER IEV
00466 INTEGER PHLUN
00467 COMMON/PHOLUN/PHLUN
00468
00469 IF (1.GT.IPOINM.AND.1.LT.IPOIN ) THEN
00470 WRITE(PHLUN,*) 'EVENT NR=',IEV,
00471 $ 'WE ARE TESTING /HEPEVT/ at IPOINT=1 (input)'
00472 CALL PHODMP
00473 ENDIF
00474 CALL PHOTOS_GET
00475 CALL PHOTOS_MAKE(ID)
00476 CALL PHOTOS_SET
00477 IF (1.GT.IPOINM.AND.1.LT.IPOIN ) THEN
00478 WRITE(PHLUN,*) 'EVENT NR=',IEV,
00479 $ 'WE ARE TESTING /HEPEVT/ at IPOINT=1 (output)'
00480 CALL PHODMP
00481 ENDIF
00482
00483 END
00484
00485 SUBROUTINE PHOTOS_GET
00486
00487
00488
00489
00490
00491
00492
00493
00494
00495
00496
00497
00498
00499
00500
00501
00502 IMPLICIT NONE
00503 INTEGER d_h_nmxhep
00504 #if defined (KORALK)
00505 PARAMETER (d_h_NMXHEP=2000)
00506 REAL*8 d_h_phep, d_h_vhep
00507 INTEGER d_h_nevhep,d_h_nhep,d_h_isthep,d_h_idhep,d_h_jmohep,
00508 $ d_h_jdahep
00509 COMMON /d_hepevt/
00510 #elif defined (CZTERYD)
00511 PARAMETER (d_h_NMXHEP=4000)
00512 REAL*8 d_h_phep, d_h_vhep
00513 INTEGER d_h_nevhep,d_h_nhep,d_h_isthep,d_h_idhep,d_h_jmohep,
00514 $ d_h_jdahep
00515 COMMON /hepevt/
00516 #elif defined (DWAR)
00517 PARAMETER (d_h_NMXHEP=2000)
00518 REAL*4 d_h_phep, d_h_vhep
00519 INTEGER d_h_nevhep,d_h_nhep,d_h_isthep,d_h_idhep,d_h_jmohep,
00520 $ d_h_jdahep
00521 COMMON /hepevt/
00522 #elif defined (DWAD)
00523 PARAMETER (d_h_NMXHEP=2000)
00524 REAL*8 d_h_phep, d_h_vhep
00525 INTEGER d_h_nevhep,d_h_nhep,d_h_isthep,d_h_idhep,d_h_jmohep,
00526 $ d_h_jdahep
00527 COMMON /hepevt/
00528 #elif defined (DZIESD)
00529 PARAMETER (d_h_NMXHEP=10000)
00530 REAL*8 d_h_phep, d_h_vhep
00531 INTEGER d_h_nevhep,d_h_nhep,d_h_isthep,d_h_idhep,d_h_jmohep,
00532 $ d_h_jdahep
00533 COMMON /hepevt/
00534 #else
00535 PARAMETER (d_h_NMXHEP=2000)
00536 REAL*8 d_h_phep, d_h_vhep
00537 INTEGER d_h_nevhep,d_h_nhep,d_h_isthep,d_h_idhep,d_h_jmohep,
00538 $ d_h_jdahep
00539 COMMON /hepevt/
00540 #endif
00541 $ d_h_nevhep,
00542 $ d_h_nhep,
00543 $ d_h_isthep(d_h_nmxhep),
00544 $ d_h_idhep(d_h_nmxhep),
00545 $ d_h_jmohep(2,d_h_nmxhep),
00546 $ d_h_jdahep(2,d_h_nmxhep),
00547 $ d_h_phep(5,d_h_nmxhep),
00548 $ d_h_vhep(4,d_h_nmxhep)
00549
00550 LOGICAL d_h_qedrad
00551 COMMON /phoqed/
00552 $ d_h_qedrad(d_h_nmxhep)
00553 INTEGER NMXHEP
00554 PARAMETER (NMXHEP=10000)
00555 INTEGER IDHEP,ISTHEP,JDAHEP,JMOHEP,NEVHEP,NHEP
00556 REAL*8 PHEP,VHEP
00557 COMMON/PH_HEPEVT/NEVHEP,NHEP,ISTHEP(NMXHEP),IDHEP(NMXHEP),
00558 &JMOHEP(2,NMXHEP),JDAHEP(2,NMXHEP),PHEP(5,NMXHEP),VHEP(4,NMXHEP)
00559 LOGICAL QEDRAD
00560 COMMON/PH_PHOQED/QEDRAD(NMXHEP)
00561 integer k,l
00562 nevhep= d_h_nevhep
00563 nhep = d_h_nhep
00564 DO K=1,nhep
00565 isthep(k) =d_h_isthep(k)
00566 idhep(k) =d_h_idhep(k)
00567 jmohep(1,k) =d_h_jmohep(1,k)
00568 jdahep(1,k) =d_h_jdahep(1,k)
00569 jmohep(2,k) =d_h_jmohep(2,k)
00570 jdahep(2,k) =d_h_jdahep(2,k)
00571 DO l=1,4
00572 phep(l,k) =d_h_phep(l,k)
00573 vhep(l,k) =d_h_vhep(l,k)
00574 ENDDO
00575 phep(5,k) =d_h_phep(5,k)
00576 qedrad(k) =d_h_qedrad(k)
00577 ENDDO
00578 END
00579
00580
00581 SUBROUTINE PHOTOS_SET
00582
00583
00584
00585
00586
00587
00588
00589
00590
00591
00592
00593
00594
00595
00596
00597 IMPLICIT NONE
00598 INTEGER d_h_nmxhep
00599 #if defined (KORALK)
00600 PARAMETER (d_h_NMXHEP=2000)
00601 REAL*8 d_h_phep, d_h_vhep
00602 INTEGER d_h_nevhep,d_h_nhep,d_h_isthep,d_h_idhep,d_h_jmohep,
00603 $ d_h_jdahep
00604 COMMON /d_hepevt/
00605 #elif defined (CZTERYD)
00606 PARAMETER (d_h_NMXHEP=4000)
00607 REAL*8 d_h_phep, d_h_vhep
00608 INTEGER d_h_nevhep,d_h_nhep,d_h_isthep,d_h_idhep,d_h_jmohep,
00609 $ d_h_jdahep
00610 COMMON /hepevt/
00611 #elif defined (DWAR)
00612 PARAMETER (d_h_NMXHEP=2000)
00613 REAL*4 d_h_phep, d_h_vhep
00614 INTEGER d_h_nevhep,d_h_nhep,d_h_isthep,d_h_idhep,d_h_jmohep,
00615 $ d_h_jdahep
00616 COMMON /hepevt/
00617 #elif defined (DWAD)
00618 PARAMETER (d_h_NMXHEP=2000)
00619 REAL*8 d_h_phep, d_h_vhep
00620 INTEGER d_h_nevhep,d_h_nhep,d_h_isthep,d_h_idhep,d_h_jmohep,
00621 $ d_h_jdahep
00622 COMMON /hepevt/
00623 #elif defined (DZIESD)
00624 PARAMETER (d_h_NMXHEP=10000)
00625 REAL*8 d_h_phep, d_h_vhep
00626 INTEGER d_h_nevhep,d_h_nhep,d_h_isthep,d_h_idhep,d_h_jmohep,
00627 $ d_h_jdahep
00628 COMMON /hepevt/
00629 #else
00630 PARAMETER (d_h_NMXHEP=2000)
00631 REAL*8 d_h_phep, d_h_vhep
00632 INTEGER d_h_nevhep,d_h_nhep,d_h_isthep,d_h_idhep,d_h_jmohep,
00633 $ d_h_jdahep
00634 COMMON /hepevt/
00635 #endif
00636 $ d_h_nevhep,
00637 $ d_h_nhep,
00638 $ d_h_isthep(d_h_nmxhep),
00639 $ d_h_idhep(d_h_nmxhep),
00640 $ d_h_jmohep(2,d_h_nmxhep),
00641 $ d_h_jdahep(2,d_h_nmxhep),
00642 $ d_h_phep(5,d_h_nmxhep),
00643 $ d_h_vhep(4,d_h_nmxhep)
00644
00645 LOGICAL d_h_qedrad
00646 COMMON /phoqed/
00647 $ d_h_qedrad(d_h_nmxhep)
00648 INTEGER NMXHEP
00649 PARAMETER (NMXHEP=10000)
00650 INTEGER IDHEP,ISTHEP,JDAHEP,JMOHEP,NEVHEP,NHEP
00651 REAL*8 PHEP,VHEP
00652 COMMON/PH_HEPEVT/NEVHEP,NHEP,ISTHEP(NMXHEP),IDHEP(NMXHEP),
00653 &JMOHEP(2,NMXHEP),JDAHEP(2,NMXHEP),PHEP(5,NMXHEP),VHEP(4,NMXHEP)
00654 LOGICAL QEDRAD
00655 COMMON/PH_PHOQED/QEDRAD(NMXHEP)
00656 INTEGER K,L
00657
00658 d_h_nevhep= nevhep
00659 d_h_nhep = nhep
00660 DO K=1,nhep
00661 d_h_isthep(k) =isthep(k)
00662 d_h_idhep(k) =idhep(k)
00663 d_h_jmohep(1,k) =jmohep(1,k)
00664 d_h_jdahep(1,k) =jdahep(1,k)
00665 d_h_jmohep(2,k) =jmohep(2,k)
00666 d_h_jdahep(2,k) =jdahep(2,k)
00667 DO l=1,4
00668 d_h_phep(l,k) =phep(l,k)
00669 d_h_vhep(l,k) =vhep(l,k)
00670 ENDDO
00671 d_h_phep(5,k) =phep(5,k)
00672 d_h_qedrad(k) =qedrad(k)
00673 ENDDO
00674 END
00675 SUBROUTINE PHOTOS_MAKE(IPARR)
00676
00677
00678
00679
00680
00681
00682
00683
00684
00685
00686
00687
00688
00689
00690
00691
00692
00693
00694
00695
00696 IMPLICIT NONE
00697 REAL*8 PHOTON(5)
00698 INTEGER IP,IPARR,IPPAR,I,J,K,L,NLAST
00699 DOUBLE PRECISION DATA
00700 INTEGER MOTHER,POSPHO
00701 LOGICAL CASCAD
00702 INTEGER NMXHEP
00703 PARAMETER (NMXHEP=10000)
00704 INTEGER IDHEP,ISTHEP,JDAHEP,JMOHEP,NEVHEP,NHEP
00705 REAL*8 PHEP,VHEP
00706 COMMON/PH_HEPEVT/NEVHEP,NHEP,ISTHEP(NMXHEP),IDHEP(NMXHEP),
00707 &JMOHEP(2,NMXHEP),JDAHEP(2,NMXHEP),PHEP(5,NMXHEP),VHEP(4,NMXHEP)
00708 LOGICAL QEDRAD
00709 COMMON/PH_PHOQED/QEDRAD(NMXHEP)
00710 INTEGER NMXPHO
00711 PARAMETER (NMXPHO=10000)
00712 INTEGER ISTACK(0:NMXPHO),NUMIT,NTRY,KK,LL,II,NA,FIRST,LAST
00713 INTEGER FIRSTA,LASTA,IPP,IDA1,IDA2,MOTHER2,IDPHO,ISPHO
00714 REAL*8 PORIG(5,NMXPHO)
00715
00716 IPPAR=ABS(IPARR)
00717
00718 IP=IPPAR
00719 NLAST=NHEP
00720 CASCAD=.FALSE.
00721
00722
00723 IF ((JDAHEP(1,IP).EQ.0).OR.(JMOHEP(1,JDAHEP(1,IP)).NE.IP)) RETURN
00724
00725
00726
00727
00728 ISTACK(0)=IPPAR
00729
00730 NUMIT=0
00731
00732
00733 NTRY=0
00734
00735 IF (IPARR.GT.0) THEN
00736 30 CONTINUE
00737 DO I=JDAHEP(1,IP),JDAHEP(2,IP)
00738 IF (JDAHEP(1,I).NE.0.AND.JMOHEP(1,JDAHEP(1,I)).EQ.I) THEN
00739 NUMIT=NUMIT+1
00740 IF (NUMIT.GT.NMXPHO) THEN
00741 DATA=NUMIT
00742 CALL PHOERR(7,'PHOTOS',DATA)
00743 ENDIF
00744 ISTACK(NUMIT)=I
00745 ENDIF
00746 ENDDO
00747 IF(NUMIT.GT.NTRY) THEN
00748 NTRY=NTRY+1
00749 IP=ISTACK(NTRY)
00750 GOTO 30
00751 ENDIF
00752 ENDIF
00753
00754 DO 25 KK=0,NUMIT
00755 NA=NHEP
00756 FIRST=JDAHEP(1,ISTACK(KK))
00757 LAST=JDAHEP(2,ISTACK(KK))
00758 DO II=1,LAST-FIRST+1
00759 DO LL=1,5
00760 PORIG(LL,II)=PHEP(LL,FIRST+II-1)
00761 ENDDO
00762 ENDDO
00763
00764 CALL PHTYPE(ISTACK(KK))
00765
00766
00767 IF(NHEP.GT.NA) THEN
00768 DO II=1,LAST-FIRST+1
00769 IPP=FIRST+II-1
00770 FIRSTA=JDAHEP(1,IPP)
00771 LASTA=JDAHEP(2,IPP)
00772 IF(JMOHEP(1,IPP).EQ.ISTACK(KK))
00773 $ CALL PHOBOS(IPP,PORIG(1,II),PHEP(1,IPP),FIRSTA,LASTA)
00774 ENDDO
00775 ENDIF
00776 25 CONTINUE
00777
00778
00779 IF (NHEP.GT.NLAST) THEN
00780 DO 160 I=NLAST+1,NHEP
00781
00782
00783 MOTHER=JMOHEP(1,I)
00784 POSPHO=JDAHEP(2,MOTHER)+1
00785
00786 DO 90 J=1,5
00787 90 PHOTON(J)=PHEP(J,I)
00788 ISPHO =ISTHEP(I)
00789 IDPHO =IDHEP(I)
00790 MOTHER2 =JMOHEP(2,I)
00791 IDA1 =JDAHEP(1,I)
00792 IDA2 =JDAHEP(2,I)
00793
00794
00795 IF (POSPHO.NE.NHEP) THEN
00796
00797
00798
00799 DO 120 K=I,POSPHO+1,-1
00800 ISTHEP(K)=ISTHEP(K-1)
00801 QEDRAD(K)=QEDRAD(K-1)
00802 IDHEP(K)=IDHEP(K-1)
00803 DO 100 L=1,2
00804 JMOHEP(L,K)=JMOHEP(L,K-1)
00805 100 JDAHEP(L,K)=JDAHEP(L,K-1)
00806 DO 110 L=1,5
00807 110 PHEP(L,K)=PHEP(L,K-1)
00808 DO 120 L=1,4
00809 120 VHEP(L,K)=VHEP(L,K-1)
00810
00811
00812 DO 130 K=1,NHEP
00813 DO 130 L=1,2
00814 IF ((JMOHEP(L,K).NE.0).AND.(JMOHEP(L,K).GE.
00815 & POSPHO)) JMOHEP(L,K)=JMOHEP(L,K)+1
00816 IF ((JDAHEP(L,K).NE.0).AND.(JDAHEP(L,K).GE.
00817 & POSPHO)) JDAHEP(L,K)=JDAHEP(L,K)+1
00818 130 CONTINUE
00819
00820
00821 DO 140 J=1,5
00822 140 PHEP(J,POSPHO)=PHOTON(J)
00823 ENDIF
00824
00825
00826 JDAHEP(2,MOTHER)=POSPHO
00827 ISTHEP(POSPHO)=ISPHO
00828 IDHEP(POSPHO)=IDPHO
00829 JMOHEP(1,POSPHO)=MOTHER
00830 JMOHEP(2,POSPHO)=MOTHER2
00831 JDAHEP(1,POSPHO)=IDA1
00832 JDAHEP(2,POSPHO)=IDA2
00833
00834
00835 DO 150 J=1,4
00836 150 VHEP(J,POSPHO)=VHEP(J,POSPHO-1)
00837 160 CONTINUE
00838 ENDIF
00839 RETURN
00840 END
00841 SUBROUTINE PHOBOS(IP,PBOOS1,PBOOS2,FIRST,LAST)
00842
00843
00844
00845
00846
00847
00848
00849
00850
00851
00852
00853
00854
00855
00856
00857
00858
00859
00860
00861
00862
00863
00864 IMPLICIT NONE
00865 DOUBLE PRECISION BET1(3),BET2(3),GAM1,GAM2,PB,DATA
00866 INTEGER I,J,FIRST,LAST,MAXSTA,NSTACK,IP
00867 PARAMETER (MAXSTA=10000)
00868 INTEGER STACK(MAXSTA)
00869 REAL*8 PBOOS1(5),PBOOS2(5)
00870 INTEGER NMXHEP
00871 PARAMETER (NMXHEP=10000)
00872 INTEGER IDHEP,ISTHEP,JDAHEP,JMOHEP,NEVHEP,NHEP
00873 REAL*8 PHEP,VHEP
00874 COMMON/PH_HEPEVT/NEVHEP,NHEP,ISTHEP(NMXHEP),IDHEP(NMXHEP),
00875 &JMOHEP(2,NMXHEP),JDAHEP(2,NMXHEP),PHEP(5,NMXHEP),VHEP(4,NMXHEP)
00876 IF ((LAST.EQ.0).OR.(LAST.LT.FIRST)) RETURN
00877 NSTACK=0
00878 DO 10 J=1,3
00879 BET1(J)=-PBOOS1(J)/PBOOS1(5)
00880 10 BET2(J)=PBOOS2(J)/PBOOS2(5)
00881 GAM1=PBOOS1(4)/PBOOS1(5)
00882 GAM2=PBOOS2(4)/PBOOS2(5)
00883
00884
00885 20 DO 50 I=FIRST,LAST
00886 PB=BET1(1)*PHEP(1,I)+BET1(2)*PHEP(2,I)+BET1(3)*PHEP(3,I)
00887 IF (JMOHEP(1,I).EQ.IP) THEN
00888 DO 30 J=1,3
00889 30 PHEP(J,I)=PHEP(J,I)+BET1(J)*(PHEP(4,I)+PB/(GAM1+1.D0))
00890 PHEP(4,I)=GAM1*PHEP(4,I)+PB
00891
00892
00893 PB=BET2(1)*PHEP(1,I)+BET2(2)*PHEP(2,I)+BET2(3)*PHEP(3,I)
00894 DO 40 J=1,3
00895 40 PHEP(J,I)=PHEP(J,I)+BET2(J)*(PHEP(4,I)+PB/(GAM2+1.D0))
00896 PHEP(4,I)=GAM2*PHEP(4,I)+PB
00897 IF (JDAHEP(1,I).NE.0) THEN
00898 NSTACK=NSTACK+1
00899
00900
00901 IF (NSTACK.GT.MAXSTA) THEN
00902 DATA=NSTACK
00903 CALL PHOERR(7,'PHOBOS',DATA)
00904 ENDIF
00905 STACK(NSTACK)=I
00906 ENDIF
00907 ENDIF
00908 50 CONTINUE
00909 IF (NSTACK.NE.0) THEN
00910
00911
00912 FIRST=JDAHEP(1,STACK(NSTACK))
00913 LAST=JDAHEP(2,STACK(NSTACK))
00914 IP=STACK(NSTACK)
00915 NSTACK=NSTACK-1
00916 GOTO 20
00917 ENDIF
00918 RETURN
00919 END
00920 SUBROUTINE PHOIN(IP,BOOST,NHEP0)
00921
00922
00923
00924
00925
00926
00927
00928
00929
00930
00931
00932
00933
00934
00935
00936
00937
00938
00939 IMPLICIT NONE
00940 INTEGER NMXHEP
00941 PARAMETER (NMXHEP=10000)
00942 INTEGER IDHEP,ISTHEP,JDAHEP,JMOHEP,NEVHEP,NHEP
00943 REAL*8 PHEP,VHEP
00944 COMMON/PH_HEPEVT/NEVHEP,NHEP,ISTHEP(NMXHEP),IDHEP(NMXHEP),
00945 &JMOHEP(2,NMXHEP),JDAHEP(2,NMXHEP),PHEP(5,NMXHEP),VHEP(4,NMXHEP)
00946 INTEGER NMXPHO
00947 PARAMETER (NMXPHO=10000)
00948 INTEGER IDPHO,ISTPHO,JDAPHO,JMOPHO,NEVPHO,NPHO
00949 REAL*8 PPHO,VPHO
00950 COMMON/PHOEVT/NEVPHO,NPHO,ISTPHO(NMXPHO),IDPHO(NMXPHO),
00951 &JMOPHO(2,NMXPHO),JDAPHO(2,NMXPHO),PPHO(5,NMXPHO),VPHO(4,NMXPHO)
00952 INTEGER IP,IP2,I,FIRST,LAST,LL,NA
00953 LOGICAL BOOST
00954 INTEGER J,NHEP0
00955 DOUBLE PRECISION BET(3),GAM,PB
00956 COMMON /PHOCMS/ BET,GAM
00957 LOGICAL INTERF,ISEC,ITRE,IEXP,IFTOP,IFW
00958 REAL*8 FINT,FSEC,EXPEPS
00959 COMMON /PHOKEY/ FSEC,FINT,EXPEPS,INTERF,ISEC,ITRE,IEXP,IFTOP,IFW
00960
00961
00962 FIRST=JDAHEP(1,IP)
00963 LAST =JDAHEP(2,IP)
00964 NPHO=3+LAST-FIRST+NHEP-NHEP0
00965 NEVPHO=NPHO
00966
00967 IDPHO(1)=IDHEP(IP)
00968 JDAPHO(1,1)=3
00969 JDAPHO(2,1)=3+LAST-FIRST
00970 DO I=1,5
00971 PPHO(I,1)=PHEP(I,IP)
00972 ENDDO
00973
00974 IP2=JMOHEP(2,JDAHEP(1,IP))
00975 IF((IP2.NE.0).AND.(IP2.NE.IP)) THEN
00976 IDPHO(2)=IDHEP(IP2)
00977 JDAPHO(1,2)=3
00978 JDAPHO(2,2)=3+LAST-FIRST
00979 DO I=1,5
00980 PPHO(I,2)=PHEP(I,IP2)
00981 ENDDO
00982 ELSE
00983 IDPHO(2)=0
00984 DO I=1,5
00985 PPHO(I,2)=0.0D0
00986 ENDDO
00987 ENDIF
00988
00989 DO LL=0,LAST-FIRST
00990 IDPHO(3+LL)=IDHEP(FIRST+LL)
00991 JMOPHO(1,3+LL)=JMOHEP(1,FIRST+LL)
00992 IF (JMOHEP(1,FIRST+LL).EQ.IP) JMOPHO(1,3+LL)=1
00993 DO I=1,5
00994 PPHO(I,3+LL)=PHEP(I,FIRST+LL)
00995 ENDDO
00996 ENDDO
00997 IF (NHEP.GT.NHEP0) THEN
00998
00999 NA=3+LAST-FIRST
01000 DO LL=1,NHEP-NHEP0
01001 IDPHO(NA+LL)=IDHEP(NHEP0+LL)
01002 JMOPHO(1,NA+LL)=JMOHEP(1,NHEP0+LL)
01003 IF (JMOHEP(1,NHEP0+LL).EQ.IP) JMOPHO(1,NA+LL)=1
01004 DO I=1,5
01005 PPHO(I,NA+LL)=PHEP(I,NHEP0+LL)
01006 ENDDO
01007 ENDDO
01008
01009 JDAPHO(2,1)=3+LAST-FIRST+NHEP-NHEP0
01010 ENDIF
01011 IF(IDPHO(NPHO).EQ.22)CALL PHLUPA(100001)
01012
01013 CALL PHCORK(0)
01014 IF(IDPHO(NPHO).EQ.22)CALL PHLUPA(100002)
01015
01016 IF(IFTOP) CALL PHOTWO(0)
01017 BOOST=.FALSE.
01018
01019
01020
01021
01022
01023
01024
01025
01026
01027
01028
01029
01030 IF ((ABS(PPHO(1,1)+ABS(PPHO(2,1))+ABS(PPHO(3,1))).GT.
01031 $ PPHO(5,1)*1.D-8).AND.(PPHO(5,1).NE.0)) THEN
01032
01033 BOOST=.TRUE.
01034
01035
01036
01037 DO 10 J=1,3
01038 10 BET(J)=-PPHO(J,1)/PPHO(5,1)
01039 GAM=PPHO(4,1)/PPHO(5,1)
01040 DO 30 I=JDAPHO(1,1),JDAPHO(2,1)
01041 PB=BET(1)*PPHO(1,I)+BET(2)*PPHO(2,I)+BET(3)*PPHO(3,I)
01042 DO 20 J=1,3
01043 20 PPHO(J,I)=PPHO(J,I)+BET(J)*(PPHO(4,I)+PB/(GAM+1.D0))
01044 30 PPHO(4,I)=GAM*PPHO(4,I)+PB
01045
01046 I=1
01047 PB=BET(1)*PPHO(1,I)+BET(2)*PPHO(2,I)+BET(3)*PPHO(3,I)
01048 DO J=1,3
01049 PPHO(J,I)=PPHO(J,I)+BET(J)*(PPHO(4,I)+PB/(GAM+1.D0))
01050 ENDDO
01051 PPHO(4,I)=GAM*PPHO(4,I)+PB
01052 ENDIF
01053
01054 IF(IFTOP) CALL PHOTWO(1)
01055 CALL PHLUPA(2)
01056 IF(IDPHO(NPHO).EQ.22) CALL PHLUPA(10000)
01057
01058 END
01059 SUBROUTINE PHOTWO(MODE)
01060
01061
01062
01063
01064
01065
01066
01067
01068
01069
01070
01071
01072
01073
01074
01075 IMPLICIT NONE
01076 INTEGER NMXPHO
01077 PARAMETER (NMXPHO=10000)
01078 INTEGER IDPHO,ISTPHO,JDAPHO,JMOPHO,NEVPHO,NPHO
01079 REAL*8 PPHO,VPHO
01080 COMMON/PHOEVT/NEVPHO,NPHO,ISTPHO(NMXPHO),IDPHO(NMXPHO),
01081 &JMOPHO(2,NMXPHO),JDAPHO(2,NMXPHO),PPHO(5,NMXPHO),VPHO(4,NMXPHO)
01082 DOUBLE PRECISION BET(3),GAM
01083 COMMON /PHOCMS/ BET,GAM
01084 INTEGER I,MODE
01085 REAL*8 MPASQR
01086 LOGICAL IFRAD
01087
01088
01089
01090
01091
01092
01093 IF(MODE.EQ.0) THEN
01094 IFRAD=(IDPHO(1).EQ.21).AND.(IDPHO(2).EQ.21)
01095 IFRAD=IFRAD.OR.(IDPHO(1).EQ.-IDPHO(2).AND.ABS(IDPHO(1)).LE.6)
01096 IFRAD=IFRAD
01097 & .AND.(ABS(IDPHO(3)).EQ.6).AND.(ABS(IDPHO(4)).EQ.6)
01098 MPASQR= (PPHO(4,1)+PPHO(4,2))**2-(PPHO(3,1)+PPHO(3,2))**2
01099 & -(PPHO(2,1)+PPHO(2,2))**2-(PPHO(1,1)+PPHO(1,2))**2
01100 IFRAD=IFRAD.AND.(MPASQR.GT.0.0D0)
01101 IF(IFRAD) THEN
01102
01103 DO I=1,4
01104 PPHO(I,1)=PPHO(I,1)+PPHO(I,2)
01105 ENDDO
01106 PPHO(5,1)=SQRT(MPASQR)
01107
01108 DO I=1,5
01109 PPHO(I,2)=0.0D0
01110 ENDDO
01111 ENDIF
01112 ELSE
01113
01114
01115
01116 ENDIF
01117 END
01118 SUBROUTINE PHOOUT(IP,BOOST,NHEP0)
01119
01120
01121
01122
01123
01124
01125
01126
01127
01128
01129
01130
01131
01132
01133
01134
01135
01136
01137 IMPLICIT NONE
01138 INTEGER NMXHEP
01139 PARAMETER (NMXHEP=10000)
01140 INTEGER IDHEP,ISTHEP,JDAHEP,JMOHEP,NEVHEP,NHEP
01141 REAL*8 PHEP,VHEP
01142 COMMON/PH_HEPEVT/NEVHEP,NHEP,ISTHEP(NMXHEP),IDHEP(NMXHEP),
01143 &JMOHEP(2,NMXHEP),JDAHEP(2,NMXHEP),PHEP(5,NMXHEP),VHEP(4,NMXHEP)
01144 INTEGER NMXPHO
01145 PARAMETER (NMXPHO=10000)
01146 INTEGER IDPHO,ISTPHO,JDAPHO,JMOPHO,NEVPHO,NPHO
01147 REAL*8 PPHO,VPHO
01148 COMMON/PHOEVT/NEVPHO,NPHO,ISTPHO(NMXPHO),IDPHO(NMXPHO),
01149 &JMOPHO(2,NMXPHO),JDAPHO(2,NMXPHO),PPHO(5,NMXPHO),VPHO(4,NMXPHO)
01150 INTEGER IP,LL,FIRST,LAST,I
01151 LOGICAL BOOST
01152 INTEGER NN,J,K,NHEP0,NA
01153 DOUBLE PRECISION BET(3),GAM,PB
01154 COMMON /PHOCMS/ BET,GAM
01155 IF(NPHO.EQ.NEVPHO) RETURN
01156
01157 CALL PHLUPA(10)
01158 IF (BOOST) THEN
01159 DO 110 J=JDAPHO(1,1),JDAPHO(2,1)
01160 PB=-BET(1)*PPHO(1,J)-BET(2)*PPHO(2,J)-BET(3)*PPHO(3,J)
01161 DO 100 K=1,3
01162 100 PPHO(K,J)=PPHO(K,J)-BET(K)*(PPHO(4,J)+PB/(GAM+1.D0))
01163 110 PPHO(4,J)=GAM*PPHO(4,J)+PB
01164
01165 DO NN=NEVPHO+1,NPHO
01166 PB=-BET(1)*PPHO(1,NN)-BET(2)*PPHO(2,NN)-BET(3)*PPHO(3,NN)
01167 DO 120 K=1,3
01168 120 PPHO(K,NN)=PPHO(K,NN)-BET(K)*(PPHO(4,NN)+PB/(GAM+1.D0))
01169 PPHO(4,NN)=GAM*PPHO(4,NN)+PB
01170 ENDDO
01171 ENDIF
01172 FIRST=JDAHEP(1,IP)
01173 LAST =JDAHEP(2,IP)
01174
01175 DO LL=0,LAST-FIRST
01176 IDHEP(FIRST+LL) = IDPHO(3+LL)
01177 DO I=1,5
01178 PHEP(I,FIRST+LL) = PPHO(I,3+LL)
01179 ENDDO
01180 ENDDO
01181
01182 NA=3+LAST-FIRST
01183 DO LL=1,NPHO-NA
01184 IDHEP(NHEP0+LL) = IDPHO(NA+LL)
01185 ISTHEP(NHEP0+LL)=ISTPHO(NA+LL)
01186 JMOHEP(1,NHEP0+LL)=IP
01187 JMOHEP(2,NHEP0+LL)=JMOHEP(2,JDAHEP(1,IP))
01188 JDAHEP(1,NHEP0+LL)=0
01189 JDAHEP(2,NHEP0+LL)=0
01190 DO I=1,5
01191 PHEP(I,NHEP0+LL) = PPHO(I,NA+LL)
01192 ENDDO
01193 ENDDO
01194 NHEP=NHEP+NPHO-NEVPHO
01195 CALL PHLUPA(20)
01196 END
01197 SUBROUTINE PHOCHK(JFIRST)
01198
01199
01200
01201
01202
01203
01204
01205
01206
01207
01208
01209
01210
01211
01212
01213 IMPLICIT NONE
01214 INTEGER NMXPHO
01215 PARAMETER (NMXPHO=10000)
01216 INTEGER IDPHO,ISTPHO,JDAPHO,JMOPHO,NEVPHO,NPHO
01217 REAL*8 PPHO,VPHO
01218 COMMON/PHOEVT/NEVPHO,NPHO,ISTPHO(NMXPHO),IDPHO(NMXPHO),
01219 &JMOPHO(2,NMXPHO),JDAPHO(2,NMXPHO),PPHO(5,NMXPHO),VPHO(4,NMXPHO)
01220 LOGICAL CHKIF
01221 COMMON/PHOIF/CHKIF(NMXPHO)
01222 INTEGER NMXHEP
01223 PARAMETER (NMXHEP=10000)
01224 LOGICAL QEDRAD
01225 COMMON/PH_PHOQED/QEDRAD(NMXHEP)
01226 INTEGER JFIRST
01227 LOGICAL F
01228 INTEGER IDABS,NLAST,I,IPPAR
01229 LOGICAL INTERF,ISEC,ITRE,IEXP,IFTOP,IFW,IFNPI0,IFKL
01230 REAL*8 FINT,FSEC,EXPEPS
01231 COMMON /PHOKEY/ FSEC,FINT,EXPEPS,INTERF,ISEC,ITRE,IEXP,IFTOP,IFW
01232 LOGICAL IFRAD
01233 INTEGER IDENT,K,IQRK,IPHQRK,IEKL,IPHEKL
01234
01235 F(IDABS)=
01236 & ( ((IDABS.GT.9.OR.IQRK.NE.1).AND.(IDABS.LE.40))
01237 & .OR.(IDABS.GT.100) )
01238 & .AND.(IDABS.NE.21)
01239 $ .AND.(IDABS.NE.2101).AND.(IDABS.NE.3101).AND.(IDABS.NE.3201)
01240 & .AND.(IDABS.NE.1103).AND.(IDABS.NE.2103).AND.(IDABS.NE.2203)
01241 & .AND.(IDABS.NE.3103).AND.(IDABS.NE.3203).AND.(IDABS.NE.3303)
01242
01243 IQRK=IPHQRK(0)
01244 IEKL=IPHEKL(0)
01245 NLAST = NPHO
01246
01247 IPPAR=1
01248
01249 IFNPI0=.TRUE.
01250 IF (IEKL.GT.1) THEN
01251
01252 IFNPI0= (IDPHO(1).NE.111)
01253 IFKL = ((IDPHO(1).EQ.130).AND.
01254 $ ((IDPHO(3).EQ.22).OR.(IDPHO(4).EQ.22).OR.
01255 $ (IDPHO(5).EQ.22)).AND.
01256 $ ((IDPHO(3).EQ.11).OR.(IDPHO(4).EQ.11).OR.
01257 $ (IDPHO(5).EQ.11)) )
01258
01259 IFNPI0=(IFNPI0.AND.(.NOT.IFKL))
01260 ENDIF
01261 DO 10 I=IPPAR,NLAST
01262 IDABS = ABS(IDPHO(I))
01263
01264 CHKIF(I)= F(IDABS) .AND.F(ABS(IDPHO(1)))
01265 & .AND. (IDPHO(2).EQ.0)
01266 IF(I.GT.2) CHKIF(I)=CHKIF(I).AND.QEDRAD(JFIRST+I-IPPAR-2)
01267 & .AND.IFNPI0
01268 10 CONTINUE
01269
01270
01271
01272 IF(IFTOP) THEN
01273
01274 DO K=JDAPHO(2,1),JDAPHO(1,1),-1
01275 IF(IDPHO(K).NE.22) THEN
01276 IDENT=K
01277 GOTO 15
01278 ENDIF
01279 ENDDO
01280 15 CONTINUE
01281 IFRAD=((IDPHO(1).EQ.21).AND.(IDPHO(2).EQ.21))
01282 & .OR. ((ABS(IDPHO(1)).LE.6).AND.((IDPHO(2)).EQ.(-IDPHO(1))))
01283 IFRAD=IFRAD
01284 & .AND.(ABS(IDPHO(3)).EQ.6).AND.((IDPHO(4)).EQ.(-IDPHO(3)))
01285 & .AND.(IDENT.EQ.4)
01286 IF(IFRAD) THEN
01287 DO 20 I=IPPAR,NLAST
01288 CHKIF(I)= .TRUE.
01289 IF(I.GT.2) CHKIF(I)=CHKIF(I).AND.QEDRAD(JFIRST+I-IPPAR-2)
01290 20 CONTINUE
01291 ENDIF
01292 ENDIF
01293
01294
01295 IF(IFTOP) THEN
01296
01297 DO K=JDAPHO(2,1),JDAPHO(1,1),-1
01298 IF(IDPHO(K).NE.22) THEN
01299 IDENT=K
01300 GOTO 25
01301 ENDIF
01302 ENDDO
01303 25 CONTINUE
01304 IFRAD=((ABS(IDPHO(1)).EQ.6).AND.(IDPHO(2).EQ.0))
01305 IFRAD=IFRAD
01306 & .AND.((ABS(IDPHO(3)).EQ.24).AND.(ABS(IDPHO(4)).EQ.5)
01307 & .OR.(ABS(IDPHO(3)).EQ.5).AND.(ABS(IDPHO(4)).EQ.24))
01308 & .AND.(IDENT.EQ.4)
01309 IF(IFRAD) THEN
01310 DO 30 I=IPPAR,NLAST
01311 CHKIF(I)= .TRUE.
01312 IF(I.GT.2) CHKIF(I)=CHKIF(I).AND.QEDRAD(JFIRST+I-IPPAR-2)
01313 30 CONTINUE
01314 ENDIF
01315 ENDIF
01316
01317
01318 END
01319 SUBROUTINE PHTYPE(ID)
01320
01321
01322
01323
01324
01325
01326
01327
01328
01329
01330
01331
01332
01333
01334
01335
01336 IMPLICIT NONE
01337 INTEGER NMXHEP
01338 PARAMETER (NMXHEP=10000)
01339 INTEGER IDHEP,ISTHEP,JDAHEP,JMOHEP,NEVHEP,NHEP
01340 REAL*8 PHEP,VHEP
01341 COMMON/PH_HEPEVT/NEVHEP,NHEP,ISTHEP(NMXHEP),IDHEP(NMXHEP),
01342 &JMOHEP(2,NMXHEP),JDAHEP(2,NMXHEP),PHEP(5,NMXHEP),VHEP(4,NMXHEP)
01343 LOGICAL INTERF,ISEC,ITRE,IEXP,IFTOP,IFW
01344 REAL*8 FINT,FSEC,EXPEPS
01345 COMMON /PHOKEY/ FSEC,FINT,EXPEPS,INTERF,ISEC,ITRE,IEXP,IFTOP,IFW
01346 LOGICAL EXPINI
01347 INTEGER NX,K,NCHAN
01348 PARAMETER (NX=10)
01349 REAL*8 PRO,PRSUM,ESU
01350 COMMON /PHOEXP/ PRO(NX),NCHAN,EXPINI
01351
01352 INTEGER ID,NHEP0
01353 LOGICAL IPAIR
01354 REAL*8 RN,PHORAN,SUM
01355 INTEGER WTDUM
01356 LOGICAL IFOUR
01357
01358 IFOUR=(.TRUE.).AND.(ITRE)
01359
01360 IPAIR=.TRUE.
01361
01362 IF (JDAHEP(1,ID).EQ.0) RETURN
01363
01364
01365 NHEP0=NHEP
01366
01367 IF (IEXP) THEN
01368 EXPINI=.TRUE.
01369 DO NCHAN=1,NX
01370 PRO(NCHAN)=0.D0
01371 ENDDO
01372 NCHAN=0
01373
01374 FSEC=1.0D0
01375 CALL PHOMAK(ID,NHEP0)
01376
01377 EXPINI=.FALSE.
01378 RN=PHORAN(WTDUM)
01379 PRSUM=0
01380 DO K=1,NX
01381 PRSUM=PRSUM+PRO(K)
01382 ENDDO
01383 ESU=EXP(-PRSUM)
01384
01385
01386 SUM=ESU
01387
01388 DO K=1,100
01389 IF(RN.LT.SUM) GOTO 100
01390 ESU=ESU*PRSUM/K
01391 SUM=SUM+ESU
01392 NCHAN=0
01393 CALL PHOMAK(ID,NHEP0)
01394 IF(SUM.GT.1D0-EXPEPS) GOTO 100
01395 ENDDO
01396 100 CONTINUE
01397 ELSEIF(IFOUR) THEN
01398
01399 FSEC=1.0D0
01400 RN=PHORAN(WTDUM)
01401 IF (RN.GE.23.D0/24D0) THEN
01402 CALL PHOMAK(ID,NHEP0)
01403 CALL PHOMAK(ID,NHEP0)
01404 CALL PHOMAK(ID,NHEP0)
01405 CALL PHOMAK(ID,NHEP0)
01406 ELSEIF (RN.GE.17.D0/24D0) THEN
01407 CALL PHOMAK(ID,NHEP0)
01408 CALL PHOMAK(ID,NHEP0)
01409 ELSEIF (RN.GE.9.D0/24D0) THEN
01410 CALL PHOMAK(ID,NHEP0)
01411 ENDIF
01412 ELSEIF(ITRE) THEN
01413
01414 FSEC=1.0D0
01415 RN=PHORAN(WTDUM)
01416 IF (RN.GE.5.D0/6D0) THEN
01417 CALL PHOMAK(ID,NHEP0)
01418 CALL PHOMAK(ID,NHEP0)
01419 CALL PHOMAK(ID,NHEP0)
01420 ELSEIF (RN.GE.2.D0/6D0) THEN
01421 CALL PHOMAK(ID,NHEP0)
01422 ENDIF
01423 ELSEIF(ISEC) THEN
01424
01425 FSEC=1.0D0
01426 RN=PHORAN(WTDUM)
01427 IF (RN.GE.0.5D0) THEN
01428 CALL PHOMAK(ID,NHEP0)
01429 CALL PHOMAK(ID,NHEP0)
01430 ENDIF
01431 ELSE
01432
01433 FSEC=1.0D0
01434 CALL PHOMAK(ID,NHEP0)
01435 ENDIF
01436
01437
01438
01439 END
01440 SUBROUTINE PHOMAK(IPPAR,NHEP0)
01441
01442
01443
01444
01445
01446
01447
01448
01449
01450
01451
01452
01453
01454
01455
01456
01457
01458
01459
01460 IMPLICIT NONE
01461 DOUBLE PRECISION DATA
01462 REAL*8 PHORAN
01463 INTEGER IP,IPPAR,NCHARG
01464 INTEGER WTDUM,IDUM,NHEP0
01465 INTEGER NCHARB,NEUDAU
01466 REAL*8 RN,WT,PHINT
01467 LOGICAL BOOST
01468 INTEGER NMXHEP
01469 PARAMETER (NMXHEP=10000)
01470 INTEGER IDHEP,ISTHEP,JDAHEP,JMOHEP,NEVHEP,NHEP
01471 REAL*8 PHEP,VHEP
01472 COMMON/PH_HEPEVT/NEVHEP,NHEP,ISTHEP(NMXHEP),IDHEP(NMXHEP),
01473 &JMOHEP(2,NMXHEP),JDAHEP(2,NMXHEP),PHEP(5,NMXHEP),VHEP(4,NMXHEP)
01474 LOGICAL INTERF,ISEC,ITRE,IEXP,IFTOP,IFW
01475 REAL*8 FINT,FSEC,EXPEPS
01476 COMMON /PHOKEY/ FSEC,FINT,EXPEPS,INTERF,ISEC,ITRE,IEXP,IFTOP,IFW
01477
01478 IP=IPPAR
01479 IDUM=1
01480 NCHARG=0
01481
01482 CALL PHOIN(IP,BOOST,NHEP0)
01483 CALL PHOCHK(JDAHEP(1,IP))
01484 WT=0.0D0
01485 CALL PHOPRE(1,WT,NEUDAU,NCHARB)
01486
01487 IF (WT.EQ.0.0D0) RETURN
01488 RN=PHORAN(WTDUM)
01489
01490 CALL PHODO(1,NCHARB,NEUDAU)
01491
01492 IF (INTERF) WT=WT*PHINT(IDUM) /FINT
01493 IF (IFW) CALL PHOBW(WT)
01494 DATA=WT
01495 IF (WT.GT.1.0D0) CALL PHOERR(3,'WT_INT',DATA)
01496
01497 IF (RN.LE.WT) THEN
01498 CALL PHOOUT(IP,BOOST,NHEP0)
01499 ENDIF
01500 RETURN
01501 END
01502 FUNCTION PHINT1(IDUM)
01503
01504
01505
01506
01507
01508
01509
01510
01511
01512
01513
01514
01515
01516
01517
01518
01519
01520
01521
01522 IMPLICIT NONE
01523 REAL*8 PHINT,phint1
01524 REAL*8 PHOCHA
01525 INTEGER IDUM
01526 INTEGER NMXPHO
01527 PARAMETER (NMXPHO=10000)
01528 INTEGER IDPHO,ISTPHO,JDAPHO,JMOPHO,NEVPHO,NPHO
01529 REAL*8 PPHO,VPHO
01530 COMMON/PHOEVT/NEVPHO,NPHO,ISTPHO(NMXPHO),IDPHO(NMXPHO),
01531 &JMOPHO(2,NMXPHO),JDAPHO(2,NMXPHO),PPHO(5,NMXPHO),VPHO(4,NMXPHO)
01532 DOUBLE PRECISION MCHSQR,MNESQR
01533 REAL*8 PNEUTR
01534 COMMON/PHOMOM/MCHSQR,MNESQR,PNEUTR(5)
01535 DOUBLE PRECISION COSTHG,SINTHG
01536 REAL*8 XPHMAX,XPHOTO
01537 COMMON/PHOPHS/XPHMAX,XPHOTO,COSTHG,SINTHG
01538 REAL*8 MPASQR,XX,BETA
01539 LOGICAL IFINT
01540 INTEGER K,IDENT
01541
01542 DO K=JDAPHO(2,1),JDAPHO(1,1),-1
01543 IF(IDPHO(K).NE.22) THEN
01544 IDENT=K
01545 GOTO 20
01546 ENDIF
01547 ENDDO
01548 20 CONTINUE
01549
01550 IFINT= NPHO.GT.IDENT
01551
01552 IFINT= IFINT.AND.(IDENT-JDAPHO(1,1)).EQ.1
01553
01554 IFINT= IFINT.AND.IDPHO(JDAPHO(1,1)).EQ.-IDPHO(IDENT)
01555
01556 IFINT= IFINT.AND.PHOCHA(IDPHO(IDENT)).NE.0
01557
01558 IF(IFINT) THEN
01559 MPASQR = PPHO(5,1)**2
01560 XX=4.D0*MCHSQR/MPASQR*(1.D0-XPHOTO)/(1.D0-XPHOTO+(MCHSQR-MNESQR)
01561 & /MPASQR)**2
01562 BETA=SQRT(1.D0-XX)
01563 PHINT = 2D0/(1D0+COSTHG**2*BETA**2)
01564 ELSE
01565 PHINT = 1D0
01566 ENDIF
01567 phint1=1
01568 END
01569
01570 FUNCTION PHINT2(IDUM)
01571
01572
01573
01574
01575
01576
01577
01578
01579
01580
01581
01582
01583
01584
01585
01586
01587
01588
01589
01590 IMPLICIT NONE
01591 REAL*8 PHINT,PHINT1,PHINT2
01592 REAL*8 PHOCHA
01593 INTEGER IDUM
01594 INTEGER NMXPHO
01595 PARAMETER (NMXPHO=10000)
01596 INTEGER IDPHO,ISTPHO,JDAPHO,JMOPHO,NEVPHO,NPHO
01597 REAL*8 PPHO,VPHO
01598 COMMON/PHOEVT/NEVPHO,NPHO,ISTPHO(NMXPHO),IDPHO(NMXPHO),
01599 &JMOPHO(2,NMXPHO),JDAPHO(2,NMXPHO),PPHO(5,NMXPHO),VPHO(4,NMXPHO)
01600 DOUBLE PRECISION MCHSQR,MNESQR
01601 REAL*8 PNEUTR
01602 COMMON/PHOMOM/MCHSQR,MNESQR,PNEUTR(5)
01603 DOUBLE PRECISION COSTHG,SINTHG
01604 REAL*8 XPHMAX,XPHOTO
01605 COMMON/PHOPHS/XPHMAX,XPHOTO,COSTHG,SINTHG
01606 REAL*8 MPASQR,XX,BETA,PQ1(4),PQ2(4),PPHOT(4)
01607 REAL*8 SS,PP2,PP,E1,E2,Q1,Q2,COSTHE
01608 LOGICAL IFINT
01609 INTEGER K,IDENT
01610
01611 DO K=JDAPHO(2,1),JDAPHO(1,1),-1
01612 IF(IDPHO(K).NE.22) THEN
01613 IDENT=K
01614 GOTO 20
01615 ENDIF
01616 ENDDO
01617 20 CONTINUE
01618
01619 IFINT= NPHO.GT.IDENT
01620
01621 IFINT= IFINT.AND.(IDENT-JDAPHO(1,1)).EQ.1
01622
01623
01624
01625 IFINT= IFINT.AND.abs(PHOCHA(IDPHO(IDENT))).GT.0.01D0
01626
01627 IFINT= IFINT.AND.
01628 $ abs(PHOCHA(IDPHO(JDAPHO(1,1)))).gt.0.01D0
01629
01630 IF(IFINT) THEN
01631 MPASQR = PPHO(5,1)**2
01632 XX=4.D0*MCHSQR/MPASQR*(1.-XPHOTO)/(1.-XPHOTO+(MCHSQR-MNESQR)/
01633 & MPASQR)**2
01634 BETA=SQRT(1.D0-XX)
01635 PHINT = 2D0/(1D0+COSTHG**2*BETA**2)
01636 SS =MPASQR*(1.D0-XPHOTO)
01637 PP2=((SS-MCHSQR-MNESQR)**2-4*MCHSQR*MNESQR)/SS/4
01638 PP =SQRT(PP2)
01639 E1 =SQRT(PP2+MCHSQR)
01640 E2 =SQRT(PP2+MNESQR)
01641 PHINT= (E1+E2)**2/((E2+COSTHG*PP)**2+(E1-COSTHG*PP)**2)
01642
01643 q1=PHOCHA(IDPHO(JDAPHO(1,1)))
01644 q2=PHOCHA(IDPHO(IDENT))
01645 do k=1,4
01646 pq1(k)=ppho(k,JDAPHO(1,1))
01647 pq2(k)=ppho(k,JDAPHO(1,1)+1)
01648 pphot(k)=ppho(k,npho)
01649 enddo
01650 costhe=(pphot(1)*pq1(1)+pphot(2)*pq1(2)+pphot(3)*pq1(3))
01651 costhe=costhe/sqrt(pq1(1)**2+pq1(2)**2+pq1(3)**2)
01652 costhe=costhe/sqrt(pphot(1)**2+pphot(2)**2+pphot(3)**2)
01653
01654
01655
01656
01657
01658 IF (costhg*costhe.GT.0) then
01659
01660 PHINT= (q1*(E2+COSTHG*PP)-q2*(E1-COSTHG*PP))**2
01661 & /(q1**2*(E2+COSTHG*PP)**2+q2**2*(E1-COSTHG*PP)**2)
01662 ELSE
01663
01664 PHINT= (q1*(E1-COSTHG*PP)-q2*(E2+COSTHG*PP))**2
01665 & /(q1**2*(E1-COSTHG*PP)**2+q2**2*(E2+COSTHG*PP)**2)
01666 ENDIF
01667 ELSE
01668 PHINT = 1D0
01669 ENDIF
01670 phint1=1
01671 phint2=1
01672 END
01673
01674
01675 SUBROUTINE PHOPRE(IPARR,WT,NEUDAU,NCHARB)
01676
01677
01678
01679
01680
01681
01682
01683
01684
01685
01686
01687
01688
01689
01690
01691
01692
01693
01694
01695
01696
01697
01698 IMPLICIT NONE
01699 DOUBLE PRECISION MINMAS,MPASQR,MCHREN
01700 DOUBLE PRECISION BETA,EPS,DEL1,DEL2,DATA,BIGLOG
01701 REAL*8 PHOCHA,PHOSPI,PHORAN,PHOCOR,MASSUM
01702 INTEGER IP,IPARR,IPPAR,I,J,ME,NCHARG,NEUPOI,NLAST,THEDUM
01703 INTEGER IDABS,IDUM
01704 INTEGER NCHARB,NEUDAU
01705 REAL*8 WT,WGT
01706 INTEGER NMXPHO
01707 PARAMETER (NMXPHO=10000)
01708 INTEGER IDPHO,ISTPHO,JDAPHO,JMOPHO,NEVPHO,NPHO
01709 REAL*8 PPHO,VPHO
01710 COMMON/PHOEVT/NEVPHO,NPHO,ISTPHO(NMXPHO),IDPHO(NMXPHO),
01711 &JMOPHO(2,NMXPHO),JDAPHO(2,NMXPHO),PPHO(5,NMXPHO),VPHO(4,NMXPHO)
01712 LOGICAL CHKIF
01713 COMMON/PHOIF/CHKIF(NMXPHO)
01714 INTEGER CHAPOI(NMXPHO)
01715 DOUBLE PRECISION MCHSQR,MNESQR
01716 REAL*8 PNEUTR
01717 COMMON/PHOMOM/MCHSQR,MNESQR,PNEUTR(5)
01718 DOUBLE PRECISION COSTHG,SINTHG
01719 REAL*8 XPHMAX,XPHOTO
01720 COMMON/PHOPHS/XPHMAX,XPHOTO,COSTHG,SINTHG
01721 REAL*8 ALPHA,XPHCUT
01722 COMMON/PHOCOP/ALPHA,XPHCUT
01723 INTEGER IREP
01724 REAL*8 PROBH,CORWT,XF
01725 COMMON/PHOPRO/PROBH,CORWT,XF,IREP
01726
01727 LOGICAL INTERF,ISEC,ITRE,IEXP,IFTOP,IFW
01728 REAL*8 FINT,FSEC,EXPEPS
01729 COMMON /PHOKEY/ FSEC,FINT,EXPEPS,INTERF,ISEC,ITRE,IEXP,IFTOP,IFW
01730
01731
01732 IPPAR=IPARR
01733
01734 IP=IPPAR
01735 NLAST=NPHO
01736 IDUM=1
01737
01738
01739 IF (JDAPHO(1,IP).EQ.0) RETURN
01740
01741
01742 10 NCHARG=0
01743 IREP=0
01744 MINMAS=0.D0
01745 MASSUM=0.D0
01746 DO 20 I=JDAPHO(1,IP),JDAPHO(2,IP)
01747
01748
01749
01750 IDABS=ABS(IDPHO(I))
01751 IF (CHKIF(I-JDAPHO(1,IP)+3)) THEN
01752 IF (PHOCHA(IDPHO(I)).NE.0) THEN
01753 NCHARG=NCHARG+1
01754 IF (NCHARG.GT.NMXPHO) THEN
01755 DATA=NCHARG
01756 CALL PHOERR(1,'PHOTOS',DATA)
01757 ENDIF
01758 CHAPOI(NCHARG)=I
01759 ENDIF
01760 MINMAS=MINMAS+PPHO(5,I)**2
01761 ENDIF
01762 MASSUM=MASSUM+PPHO(5,I)
01763 20 CONTINUE
01764 IF (NCHARG.NE.0) THEN
01765
01766
01767 IF ((PPHO(5,IP)-MASSUM)/PPHO(5,IP).GT.2.D0*XPHCUT) THEN
01768
01769
01770
01771 IF (NCHARG.GT.1) CALL PHOOMA(1,NCHARG,CHAPOI)
01772
01773 30 CONTINUE
01774 DO 70 J=1,3
01775 70 PNEUTR(J)=-PPHO(J,CHAPOI(NCHARG))
01776 PNEUTR(4)=PPHO(5,IP)-PPHO(4,CHAPOI(NCHARG))
01777
01778
01779 MPASQR=PPHO(5,IP)**2
01780 MCHSQR=PPHO(5,CHAPOI(NCHARG))**2
01781 IF ((JDAPHO(2,IP)-JDAPHO(1,IP)).EQ.1) THEN
01782 NEUPOI=JDAPHO(1,IP)
01783 IF (NEUPOI.EQ.CHAPOI(NCHARG)) NEUPOI=JDAPHO(2,IP)
01784 MNESQR=PPHO(5,NEUPOI)**2
01785 PNEUTR(5)=PPHO(5,NEUPOI)
01786 ELSE
01787 MNESQR=PNEUTR(4)**2-PNEUTR(1)**2-PNEUTR(2)**2-PNEUTR(3)**2
01788 MNESQR=MAX(MNESQR,MINMAS-MCHSQR)
01789 PNEUTR(5)=SQRT(MNESQR)
01790 ENDIF
01791
01792
01793 XPHMAX=(MPASQR-(PNEUTR(5)+PPHO(5,CHAPOI(NCHARG)))**2)/MPASQR
01794
01795
01796 CALL PHOENE(MPASQR,MCHREN,BETA,BIGLOG,IDPHO(CHAPOI(NCHARG)))
01797
01798 IF (XPHOTO.LT.-4D0) THEN
01799 NCHARG=0
01800 XPHOTO=0d0
01801
01802 ELSEIF ((XPHOTO.LT.XPHCUT).OR.(XPHOTO.GT.XPHMAX)) THEN
01803
01804
01805
01806 NCHARG=NCHARG-1
01807 IF (NCHARG.GT.0) THEN
01808 IREP=IREP+1
01809 GOTO 30
01810 ENDIF
01811 ELSE
01812
01813
01814
01815 EPS=MCHREN/(1.D0+BETA)
01816
01817
01818 DEL1=(2.D0-EPS)*(EPS/(2.D0-EPS))**PHORAN(THEDUM)
01819 DEL2=2.D0-DEL1
01820
01821
01822
01823
01824
01825
01826
01827
01828
01829
01830
01831
01832
01833
01834
01835
01836
01837
01838
01839
01840
01841
01842
01843
01844
01845 COSTHG=(1.D0-DEL1)/BETA
01846 SINTHG=SQRT(DEL1*DEL2-MCHREN)/BETA
01847 WGT=1D0
01848
01849
01850
01851
01852 ME=2.D0*PHOSPI(IDPHO(CHAPOI(NCHARG)))+1.D0
01853
01854
01855
01856
01857 DO I=JDAPHO(1,IP),JDAPHO(2,IP)
01858 IF (I.NE.CHAPOI(NCHARG)) THEN
01859 NEUDAU=I
01860 GOTO 51
01861 ENDIF
01862 ENDDO
01863
01864
01865 DATA=NCHARG
01866 CALL PHOERR(5,'PHOKIN',DATA)
01867 51 CONTINUE
01868 NCHARB=CHAPOI(NCHARG)
01869 NCHARB=NCHARB-JDAPHO(1,IP)+3
01870 NEUDAU=NEUDAU-JDAPHO(1,IP)+3
01871 WT=PHOCOR(MPASQR,MCHREN,ME)*WGT
01872 ENDIF
01873 ELSE
01874 DATA=PPHO(5,IP)-MASSUM
01875 CALL PHOERR(10,'PHOTOS',DATA)
01876 ENDIF
01877 ENDIF
01878
01879 RETURN
01880 END
01881 SUBROUTINE PHOOMA(IFIRST,ILAST,POINTR)
01882
01883
01884
01885
01886
01887
01888
01889
01890
01891
01892
01893
01894
01895
01896
01897
01898
01899
01900
01901 IMPLICIT NONE
01902 INTEGER NMXPHO
01903 PARAMETER (NMXPHO=10000)
01904 INTEGER IDPHO,ISTPHO,JDAPHO,JMOPHO,NEVPHO,NPHO
01905 REAL*8 PPHO,VPHO
01906 COMMON/PHOEVT/NEVPHO,NPHO,ISTPHO(NMXPHO),IDPHO(NMXPHO),
01907 &JMOPHO(2,NMXPHO),JDAPHO(2,NMXPHO),PPHO(5,NMXPHO),VPHO(4,NMXPHO)
01908 INTEGER IFIRST,ILAST,I,J,BUFPOI,POINTR(NMXPHO)
01909 REAL*8 BUFMAS,MASS(NMXPHO)
01910 IF (IFIRST.EQ.ILAST) RETURN
01911
01912
01913 DO 10 I=IFIRST,ILAST
01914 10 MASS(I)=PPHO(5,POINTR(I))
01915
01916
01917 DO 30 I=IFIRST,ILAST-1
01918 DO 20 J=I+1,ILAST
01919 IF (MASS(J).LE.MASS(I)) GOTO 20
01920 BUFPOI=POINTR(J)
01921 POINTR(J)=POINTR(I)
01922 POINTR(I)=BUFPOI
01923 BUFMAS=MASS(J)
01924 MASS(J)=MASS(I)
01925 MASS(I)=BUFMAS
01926 20 CONTINUE
01927 30 CONTINUE
01928 RETURN
01929 END
01930 SUBROUTINE PHOENE(MPASQR,MCHREN,BETA,BIGLOG,IDENT)
01931
01932
01933
01934
01935
01936
01937
01938
01939
01940
01941
01942
01943
01944
01945
01946
01947
01948
01949
01950
01951
01952 IMPLICIT NONE
01953 DOUBLE PRECISION MPASQR,MCHREN,BIGLOG,BETA,DATA
01954 INTEGER IWT1,IRN,IWT2
01955 REAL*8 PRSOFT,PRHARD,PHORAN,PHOFAC
01956 DOUBLE PRECISION MCHSQR,MNESQR
01957 REAL*8 PNEUTR
01958 INTEGER IDENT
01959 REAL*8 PHOCHA,PRKILL,RRR
01960 COMMON/PHOMOM/MCHSQR,MNESQR,PNEUTR(5)
01961 DOUBLE PRECISION COSTHG,SINTHG
01962 REAL*8 XPHMAX,XPHOTO
01963 COMMON/PHOPHS/XPHMAX,XPHOTO,COSTHG,SINTHG
01964 REAL*8 ALPHA,XPHCUT
01965 COMMON/PHOCOP/ALPHA,XPHCUT
01966 REAL*8 PI,TWOPI
01967 COMMON/PHPICO/PI,TWOPI
01968 INTEGER IREP
01969 REAL*8 PROBH,CORWT,XF
01970 COMMON/PHOPRO/PROBH,CORWT,XF,IREP
01971 LOGICAL INTERF,ISEC,ITRE,IEXP,IFTOP,IFW
01972 REAL*8 FINT,FSEC,EXPEPS
01973 COMMON /PHOKEY/ FSEC,FINT,EXPEPS,INTERF,ISEC,ITRE,IEXP,IFTOP,IFW
01974 INTEGER NX,NCHAN,K
01975 PARAMETER (NX=10)
01976 LOGICAL EXPINI
01977 REAL*8 PRO,PRSUM
01978 COMMON /PHOEXP/ PRO(NX),NCHAN,EXPINI
01979
01980 IF (XPHMAX.LE.XPHCUT) THEN
01981 BETA=PHOFAC(-1)
01982 XPHOTO=0.0D0
01983 RETURN
01984 ENDIF
01985
01986 MCHREN=4.D0*MCHSQR/MPASQR/(1.D0+MCHSQR/MPASQR)**2
01987 BETA=SQRT(1.D0-MCHREN)
01988
01989
01990
01991
01992
01993
01994
01995
01996
01997
01998
01999
02000 BIGLOG=LOG(MPASQR/MCHSQR*(1.D0+BETA)**2/4.D0*
02001 & (1.D0+MCHSQR/MPASQR)**2)
02002 PRHARD=ALPHA/PI*(1D0/BETA*BIGLOG)*
02003 &(LOG(XPHMAX/XPHCUT)-.75D0+XPHCUT/XPHMAX-.25D0*XPHCUT**2/XPHMAX**2)
02004 PRHARD=PRHARD*PHOCHA(IDENT)**2*FSEC*FINT
02005
02006 IF (IREP.EQ.0) PROBH=0.D0
02007 PRKILL=0d0
02008 IF (IEXP) THEN
02009 NCHAN=NCHAN+1
02010 IF (EXPINI) THEN
02011 PRO(NCHAN)=PRHARD+0.05*(1.0+FINT)
02012
02013 PRHARD=0D0
02014 PROBH=PRHARD
02015 ELSE
02016 PRSUM=0
02017 DO K=NCHAN,NX
02018 PRSUM=PRSUM+PRO(K)
02019 ENDDO
02020 PRHARD=PRHARD/PRSUM
02021
02022
02023
02024 PRKILL=PRO(NCHAN)/PRSUM-PRHARD
02025
02026 ENDIF
02027 PRSOFT=1.D0-PRHARD
02028 ELSE
02029 PRHARD=PRHARD*PHOFAC(0)
02030
02031
02032 PROBH=PRHARD
02033 ENDIF
02034 PRSOFT=1.D0-PRHARD
02035
02036
02037 IF (IEXP) THEN
02038 IF (PRSOFT.LT.-5.0D-8) THEN
02039 DATA=PRSOFT
02040 CALL PHOERR(2,'PHOENE',DATA)
02041 ENDIF
02042 ELSE
02043 IF (PRSOFT.LT.0.1D0) THEN
02044 DATA=PRSOFT
02045 CALL PHOERR(2,'PHOENE',DATA)
02046 ENDIF
02047 ENDIF
02048
02049 RRR=PHORAN(IWT1)
02050 IF (RRR.LT.PRSOFT) THEN
02051
02052
02053 XPHOTO=0.D0
02054 IF (RRR.LT.PRKILL) XPHOTO=-5d0
02055 ELSE
02056
02057
02058
02059 10 XPHOTO=EXP(PHORAN(IRN)*LOG(XPHCUT/XPHMAX))
02060 XPHOTO=XPHOTO*XPHMAX
02061 IF (PHORAN(IWT2).GT.((1.D0+(1.D0-XPHOTO/XPHMAX)**2)/2.D0))
02062 & GOTO 10
02063 ENDIF
02064
02065
02066 XF=4.D0*MCHSQR*MPASQR/(MPASQR+MCHSQR-MNESQR)**2
02067 RETURN
02068 END
02069 FUNCTION PHOCOR(MPASQR,MCHREN,ME)
02070
02071
02072
02073
02074
02075
02076
02077
02078
02079
02080
02081
02082
02083
02084
02085
02086
02087
02088
02089 IMPLICIT NONE
02090 DOUBLE PRECISION MPASQR,MCHREN,BETA,XX,YY,DATA
02091 INTEGER ME
02092 REAL*8 PHOCOR,PHOFAC,WT1,WT2,WT3
02093 DOUBLE PRECISION MCHSQR,MNESQR
02094 REAL*8 PNEUTR
02095 COMMON/PHOMOM/MCHSQR,MNESQR,PNEUTR(5)
02096 DOUBLE PRECISION COSTHG,SINTHG
02097 REAL*8 XPHMAX,XPHOTO
02098 COMMON/PHOPHS/XPHMAX,XPHOTO,COSTHG,SINTHG
02099 INTEGER IREP
02100 REAL*8 PROBH,CORWT,XF
02101 COMMON/PHOPRO/PROBH,CORWT,XF,IREP
02102
02103
02104 XX=4.D0*MCHSQR/MPASQR*(1.D0-XPHOTO)/(1.D0-XPHOTO+(MCHSQR-MNESQR)/
02105 &MPASQR)**2
02106 IF (ME.EQ.1) THEN
02107 YY=1.D0
02108 WT3=(1.D0-XPHOTO/XPHMAX)/((1.D0+(1.D0-XPHOTO/XPHMAX)**2)/2.D0)
02109 ELSEIF (ME.EQ.2) THEN
02110 YY=0.5D0*(1.D0-XPHOTO/XPHMAX+1.D0/(1.D0-XPHOTO/XPHMAX))
02111 WT3=1.D0
02112 ELSEIF ((ME.EQ.3).OR.(ME.EQ.4).OR.(ME.EQ.5)) THEN
02113 YY=1.D0
02114 WT3=(1.D0+(1.D0-XPHOTO/XPHMAX)**2-(XPHOTO/XPHMAX)**3)/
02115 & (1.D0+(1.D0-XPHOTO/XPHMAX)** 2)
02116 ELSE
02117 DATA=(ME-1.D0)/2.D0
02118 CALL PHOERR(6,'PHOCOR',DATA)
02119 YY=1.D0
02120 WT3=1.D0
02121 ENDIF
02122 BETA=SQRT(1.D0-XX)
02123 WT1=(1.D0-COSTHG*SQRT(1.D0-MCHREN))/(1.D0-COSTHG*BETA)
02124 WT2=(1.D0-XX/YY/(1.D0-BETA**2*COSTHG**2))*(1.D0+COSTHG*BETA)/2.D0
02125 WT2=WT2*PHOFAC(1)
02126 PHOCOR=WT1*WT2*WT3
02127 CORWT=PHOCOR
02128 IF (PHOCOR.GT.1.D0) THEN
02129 DATA=PHOCOR
02130 CALL PHOERR(3,'PHOCOR',DATA)
02131 ENDIF
02132 RETURN
02133 END
02134 FUNCTION PHOFAC(MODE)
02135
02136
02137
02138
02139
02140
02141
02142
02143
02144
02145
02146
02147
02148
02149
02150
02151
02152
02153
02154
02155
02156
02157
02158
02159
02160
02161
02162
02163
02164 IMPLICIT NONE
02165 REAL*8 PHOFAC,FF,PRX
02166 INTEGER MODE
02167 INTEGER IREP
02168 REAL*8 PROBH,CORWT,XF
02169 COMMON/PHOPRO/PROBH,CORWT,XF,IREP
02170 LOGICAL INTERF,ISEC,ITRE,IEXP,IFTOP,IFW
02171 REAL*8 FINT,FSEC,EXPEPS
02172 COMMON /PHOKEY/ FSEC,FINT,EXPEPS,INTERF,ISEC,ITRE,IEXP,IFTOP,IFW
02173 SAVE PRX,FF
02174 DATA PRX,FF/ 0.D0, 0.D0/
02175 IF (IEXP) THEN
02176 PHOFAC=1
02177 RETURN
02178 ENDIF
02179 IF (MODE.EQ.-1) THEN
02180 PRX=1.D0
02181 FF=1.D0
02182 PROBH=0.0
02183 ELSEIF (MODE.EQ.0) THEN
02184 IF (IREP.EQ.0) PRX=1.D0
02185 PRX=PRX/(1.D0-PROBH)
02186 FF=1.D0
02187
02188
02189
02190
02191
02192
02193 PHOFAC=FF*PRX
02194 ELSE
02195 PHOFAC=1.D0/FF
02196 ENDIF
02197 END
02198 SUBROUTINE PHOBW(WT)
02199
02200
02201
02202
02203
02204
02205
02206
02207
02208
02209
02210
02211
02212
02213
02214
02215
02216
02217
02218
02219
02220
02221 IMPLICIT NONE
02222 DOUBLE PRECISION WT
02223 INTEGER NMXPHO
02224 PARAMETER (NMXPHO=10000)
02225 INTEGER IDPHO,ISTPHO,JDAPHO,JMOPHO,NEVPHO,NPHO
02226 REAL*8 PPHO,VPHO
02227 COMMON/PHOEVT/NEVPHO,NPHO,ISTPHO(NMXPHO),IDPHO(NMXPHO),
02228 &JMOPHO(2,NMXPHO),JDAPHO(2,NMXPHO),PPHO(5,NMXPHO),VPHO(4,NMXPHO)
02229 INTEGER I
02230 DOUBLE PRECISION EMU,MCHREN,BETA,COSTHG,MPASQR,XPH
02231
02232 IF (ABS(IDPHO(1)).EQ.24.AND.
02233 $ ABS(IDPHO(JDAPHO(1,1) )).GE.11.AND.
02234 $ ABS(IDPHO(JDAPHO(1,1) )).LE.16.AND.
02235 $ ABS(IDPHO(JDAPHO(1,1)+1)).GE.11.AND.
02236 $ ABS(IDPHO(JDAPHO(1,1)+1)).LE.16 ) THEN
02237
02238 IF(
02239 $ ABS(IDPHO(JDAPHO(1,1) )).EQ.11.OR.
02240 $ ABS(IDPHO(JDAPHO(1,1) )).EQ.13.OR.
02241 $ ABS(IDPHO(JDAPHO(1,1) )).EQ.15 ) THEN
02242 I=JDAPHO(1,1)
02243 ELSE
02244 I=JDAPHO(1,1)+1
02245 ENDIF
02246 EMU=PPHO(4,I)
02247 MCHREN=ABS(PPHO(4,I)**2-PPHO(3,I)**2
02248 $ -PPHO(2,I)**2-PPHO(1,I)**2)
02249 BETA=SQRT(1- MCHREN/ PPHO(4,I)**2)
02250 COSTHG=(PPHO(3,I)*PPHO(3,NPHO)+PPHO(2,I)*PPHO(2,NPHO)
02251 $ +PPHO(1,I)*PPHO(1,NPHO))/
02252 $ SQRT(PPHO(3,I)**2+PPHO(2,I)**2+PPHO(1,I)**2) /
02253 $ SQRT(PPHO(3,NPHO)**2+PPHO(2,NPHO)**2+PPHO(1,NPHO)**2)
02254 MPASQR=PPHO(4,1)**2
02255 XPH=PPHO(4,NPHO)
02256 WT=WT*(1-8*EMU*XPH*(1-COSTHG*BETA)*
02257 $ (MCHREN+2*XPH*SQRT(MPASQR))/
02258 $ MPASQR**2/(1-MCHREN/MPASQR)/(4-MCHREN/MPASQR))
02259 ENDIF
02260
02261
02262
02263 END
02264 SUBROUTINE PHODO(IP,NCHARB,NEUDAU)
02265
02266
02267
02268
02269
02270
02271
02272
02273
02274
02275
02276
02277
02278
02279
02280
02281
02282
02283
02284
02285
02286
02287 IMPLICIT NONE
02288 DOUBLE PRECISION PHOAN1,PHOAN2,ANGLE,FI1,FI3,FI4,FI5,TH1,TH3,TH4
02289 DOUBLE PRECISION PARNE,QNEW,QOLD,DATA
02290 INTEGER IP,FI3DUM,I,J,NEUDAU,FIRST,LAST
02291 INTEGER NCHARB
02292 REAL*8 EPHOTO,PMAVIR,PHOTRI
02293 REAL*8 GNEUT,PHORAN,CCOSTH,SSINTH,PVEC(4)
02294 INTEGER NMXPHO
02295 PARAMETER (NMXPHO=10000)
02296 INTEGER IDPHO,ISTPHO,JDAPHO,JMOPHO,NEVPHO,NPHO
02297 REAL*8 PPHO,VPHO
02298 COMMON/PHOEVT/NEVPHO,NPHO,ISTPHO(NMXPHO),IDPHO(NMXPHO),
02299 &JMOPHO(2,NMXPHO),JDAPHO(2,NMXPHO),PPHO(5,NMXPHO),VPHO(4,NMXPHO)
02300 DOUBLE PRECISION MCHSQR,MNESQR
02301 REAL*8 PNEUTR
02302 COMMON/PHOMOM/MCHSQR,MNESQR,PNEUTR(5)
02303 DOUBLE PRECISION COSTHG,SINTHG
02304 REAL*8 XPHMAX,XPHOTO
02305 COMMON/PHOPHS/XPHMAX,XPHOTO,COSTHG,SINTHG
02306 REAL*8 PI,TWOPI
02307 COMMON/PHPICO/PI,TWOPI
02308
02309 EPHOTO=XPHOTO*PPHO(5,IP)/2.D0
02310 PMAVIR=SQRT(PPHO(5,IP)*(PPHO(5,IP)-2.D0*EPHOTO))
02311
02312
02313 FI1=PHOAN1(PNEUTR(1),PNEUTR(2))
02314
02315
02316
02317 TH1=PHOAN2(PNEUTR(3),SQRT(PNEUTR(1)**2+PNEUTR(2)**2))
02318 CALL PHORO3(-FI1,PNEUTR(1))
02319 CALL PHORO2(-TH1,PNEUTR(1))
02320
02321
02322
02323
02324
02325
02326
02327 QNEW=PHOTRI(PMAVIR,PNEUTR(5),PPHO(5,NCHARB))
02328 QOLD=PNEUTR(3)
02329 GNEUT=(QNEW**2+QOLD**2+MNESQR)/(QNEW*QOLD+SQRT((QNEW**2+MNESQR)*
02330 &(QOLD**2+MNESQR)))
02331 IF (GNEUT.LT.1.D0) THEN
02332 DATA=0.D0
02333 CALL PHOERR(4,'PHOKIN',DATA)
02334 ENDIF
02335 PARNE=GNEUT-SQRT(MAX(GNEUT**2-1.0D0,0.D0))
02336
02337
02338 CALL PHOBO3(PARNE,PNEUTR)
02339
02340
02341 NPHO=NPHO+1
02342 ISTPHO(NPHO)=1
02343 IDPHO(NPHO) =22
02344
02345 JMOPHO(1,NPHO)=IP
02346 JMOPHO(2,NPHO)=0
02347 JDAPHO(1,NPHO)=0
02348 JDAPHO(2,NPHO)=0
02349 PPHO(4,NPHO)=EPHOTO*PPHO(5,IP)/PMAVIR
02350
02351
02352 CCOSTH=-COSTHG
02353 SSINTH=SINTHG
02354 TH3=PHOAN2(CCOSTH,SSINTH)
02355 FI3=TWOPI*PHORAN(FI3DUM)
02356 PPHO(1,NPHO)=PPHO(4,NPHO)*SINTHG*COS(FI3)
02357 PPHO(2,NPHO)=PPHO(4,NPHO)*SINTHG*SIN(FI3)
02358
02359
02360 PPHO(3,NPHO)=-PPHO(4,NPHO)*COSTHG
02361 PPHO(5,NPHO)=0.D0
02362
02363
02364 CALL PHORO3(-FI3,PNEUTR(1))
02365 CALL PHORO3(-FI3,PPHO(1,NPHO))
02366 CALL PHORO2(-TH3,PNEUTR(1))
02367 CALL PHORO2(-TH3,PPHO(1,NPHO))
02368 ANGLE=EPHOTO/PPHO(4,NPHO)
02369
02370
02371 CALL PHOBO3(ANGLE,PNEUTR(1))
02372 CALL PHOBO3(ANGLE,PPHO(1,NPHO))
02373
02374
02375 FI4=PHOAN1(PNEUTR(1),PNEUTR(2))
02376 TH4=PHOAN2(PNEUTR(3),SQRT(PNEUTR(1)**2+PNEUTR(2)**2))
02377 CALL PHORO3(FI4,PNEUTR(1))
02378 CALL PHORO3(FI4,PPHO(1,NPHO))
02379
02380 DO 60 I=2,4
02381 60 PVEC(I)=0.D0
02382 PVEC(1)=1.D0
02383 CALL PHORO3(-FI3,PVEC)
02384 CALL PHORO2(-TH3,PVEC)
02385 CALL PHOBO3(ANGLE,PVEC)
02386 CALL PHORO3(FI4,PVEC)
02387 CALL PHORO2(-TH4,PNEUTR)
02388 CALL PHORO2(-TH4,PPHO(1,NPHO))
02389 CALL PHORO2(-TH4,PVEC)
02390 FI5=PHOAN1(PVEC(1),PVEC(2))
02391
02392
02393 CALL PHORO3(-FI5,PNEUTR)
02394 CALL PHORO3(-FI5,PPHO(1,NPHO))
02395 CALL PHORO2(TH1,PNEUTR(1))
02396 CALL PHORO2(TH1,PPHO(1,NPHO))
02397 CALL PHORO3(FI1,PNEUTR)
02398 CALL PHORO3(FI1,PPHO(1,NPHO))
02399
02400 IF ((JDAPHO(2,IP)-JDAPHO(1,IP)).GT.1) THEN
02401
02402
02403 FIRST=NEUDAU
02404 LAST=JDAPHO(2,IP)
02405 DO 70 I=FIRST,LAST
02406 IF (I.NE.NCHARB.AND.(JMOPHO(1,I).EQ.IP)) THEN
02407
02408
02409 CALL PHORO3(-FI1,PPHO(1,I))
02410 CALL PHORO2(-TH1,PPHO(1,I))
02411
02412
02413 CALL PHOBO3(PARNE,PPHO(1,I))
02414
02415
02416 CALL PHORO3(-FI3,PPHO(1,I))
02417 CALL PHORO2(-TH3,PPHO(1,I))
02418
02419
02420 CALL PHOBO3(ANGLE,PPHO(1,I))
02421
02422
02423 CALL PHORO3(FI4,PPHO(1,I))
02424 CALL PHORO2(-TH4,PPHO(1,I))
02425
02426
02427 CALL PHORO3(-FI5,PPHO(1,I))
02428 CALL PHORO2(TH1,PPHO(1,I))
02429 CALL PHORO3(FI1,PPHO(1,I))
02430 ENDIF
02431 70 CONTINUE
02432 ELSE
02433
02434
02435 DO 80 J=1,4
02436 80 PPHO(J,NEUDAU)=PNEUTR(J)
02437 ENDIF
02438
02439
02440 DO 90 J=1,3
02441 90 PPHO(J,NCHARB)=-(PPHO(J,NPHO)+PNEUTR(J))
02442 PPHO(4,NCHARB)=PPHO(5,IP)-(PPHO(4,NPHO)+PNEUTR(4))
02443
02444 END
02445 FUNCTION PHOTRI(A,B,C)
02446
02447
02448
02449
02450
02451
02452
02453
02454
02455
02456
02457
02458
02459
02460
02461 IMPLICIT NONE
02462 DOUBLE PRECISION DA,DB,DC,DAPB,DAMB,DTRIAN
02463 REAL*8 A,B,C,PHOTRI
02464 DA=A
02465 DB=B
02466 DC=C
02467 DAPB=DA+DB
02468 DAMB=DA-DB
02469 DTRIAN=SQRT((DAMB-DC)*(DAPB+DC)*(DAMB+DC)*(DAPB-DC))
02470 PHOTRI=DTRIAN/(DA+DA)
02471 RETURN
02472 END
02473 FUNCTION PHOAN1(X,Y)
02474
02475
02476
02477
02478
02479
02480
02481
02482
02483
02484
02485
02486
02487
02488 IMPLICIT NONE
02489 DOUBLE PRECISION PHOAN1
02490 REAL*8 X,Y
02491 REAL*8 PI,TWOPI
02492 COMMON/PHPICO/PI,TWOPI
02493 IF (ABS(Y).LT.ABS(X)) THEN
02494 PHOAN1=ATAN(ABS(Y/X))
02495 IF (X.LE.0.D0) PHOAN1=PI-PHOAN1
02496 ELSE
02497 PHOAN1=ACOS(X/SQRT(X**2+Y**2))
02498 ENDIF
02499 IF (Y.LT.0.D0) PHOAN1=TWOPI-PHOAN1
02500 RETURN
02501 END
02502 FUNCTION PHOAN2(X,Y)
02503
02504
02505
02506
02507
02508
02509
02510
02511
02512
02513
02514
02515
02516
02517 IMPLICIT NONE
02518 DOUBLE PRECISION PHOAN2
02519 REAL*8 X,Y
02520 REAL*8 PI,TWOPI
02521 COMMON/PHPICO/PI,TWOPI
02522 IF (ABS(Y).LT.ABS(X)) THEN
02523 PHOAN2=ATAN(ABS(Y/X))
02524 IF (X.LE.0.D0) PHOAN2=PI-PHOAN2
02525 ELSE
02526 PHOAN2=ACOS(X/SQRT(X**2+Y**2))
02527 ENDIF
02528 RETURN
02529 END
02530 SUBROUTINE PHOBO3(ANGLE,PVEC)
02531
02532
02533
02534
02535
02536
02537
02538
02539
02540
02541
02542
02543
02544
02545
02546 IMPLICIT NONE
02547 DOUBLE PRECISION QPL,QMI,ANGLE
02548 REAL*8 PVEC(4)
02549 QPL=(PVEC(4)+PVEC(3))*ANGLE
02550 QMI=(PVEC(4)-PVEC(3))/ANGLE
02551 PVEC(3)=(QPL-QMI)/2.D0
02552 PVEC(4)=(QPL+QMI)/2.D0
02553 RETURN
02554 END
02555 SUBROUTINE PHORO2(ANGLE,PVEC)
02556
02557
02558
02559
02560
02561
02562
02563
02564
02565
02566
02567
02568
02569
02570
02571 IMPLICIT NONE
02572 DOUBLE PRECISION CS,SN,ANGLE
02573 REAL*8 PVEC(4)
02574 CS=COS(ANGLE)*PVEC(1)+SIN(ANGLE)*PVEC(3)
02575 SN=-SIN(ANGLE)*PVEC(1)+COS(ANGLE)*PVEC(3)
02576 PVEC(1)=CS
02577 PVEC(3)=SN
02578 RETURN
02579 END
02580 SUBROUTINE PHORO3(ANGLE,PVEC)
02581
02582
02583
02584
02585
02586
02587
02588
02589
02590
02591
02592
02593
02594
02595
02596 IMPLICIT NONE
02597 DOUBLE PRECISION CS,SN,ANGLE
02598 REAL*8 PVEC(4)
02599 CS=COS(ANGLE)*PVEC(1)-SIN(ANGLE)*PVEC(2)
02600 SN=SIN(ANGLE)*PVEC(1)+COS(ANGLE)*PVEC(2)
02601 PVEC(1)=CS
02602 PVEC(2)=SN
02603 RETURN
02604 END
02605 #include "../randg/photos-random.h"
02606 FUNCTION PHOCHA(IDHEP)
02607
02608
02609
02610
02611
02612
02613
02614
02615
02616
02617
02618
02619
02620
02621
02622
02623
02624 IMPLICIT NONE
02625 REAL*8 PHOCHA
02626 INTEGER IDHEP,IDABS,Q1,Q2,Q3
02627
02628
02629
02630 REAL*8 CHARGE(0:100)
02631 DATA CHARGE/ 0.D0,
02632 &-0.3333333333D0, 0.6666666667D0, -0.3333333333D0, 0.6666666667D0,
02633 &-0.3333333333D0, 0.6666666667D0, -0.3333333333D0, 0.6666666667D0,
02634 & 2*0.D0, -1.D0, 0.D0, -1.D0, 0.D0, -1.D0, 0.D0, -1.D0, 6*0.D0,
02635 & 1.D0, 12*0.D0, 1.D0, 63*0.D0/
02636 IDABS=ABS(IDHEP)
02637 IF (IDABS.LE.100) THEN
02638
02639
02640 PHOCHA = CHARGE(IDABS)
02641 ELSE
02642
02643
02644 Q3=MOD(IDABS/1000,10)
02645 Q2=MOD(IDABS/100,10)
02646 Q1=MOD(IDABS/10,10)
02647 IF (Q3.EQ.0) THEN
02648
02649
02650 IF(MOD(Q2,2).EQ.0) THEN
02651 PHOCHA=CHARGE(Q2)-CHARGE(Q1)
02652 ELSE
02653 PHOCHA=CHARGE(Q1)-CHARGE(Q2)
02654 ENDIF
02655 ELSE
02656
02657
02658 PHOCHA=CHARGE(Q1)+CHARGE(Q2)+CHARGE(Q3)
02659 ENDIF
02660 ENDIF
02661
02662
02663 IF (IDHEP.LT.0.D0) PHOCHA=-PHOCHA
02664 IF (PHOCHA**2.lt.1d-6) PHOCHA=0.D0
02665 RETURN
02666 END
02667 FUNCTION PHOSPI(IDHEP)
02668
02669
02670
02671
02672
02673
02674
02675
02676
02677
02678
02679
02680
02681
02682
02683
02684
02685
02686 IMPLICIT NONE
02687 REAL*8 PHOSPI
02688 INTEGER IDHEP,IDABS
02689
02690
02691
02692 REAL*8 SPIN(100)
02693 DATA SPIN/ 8*.5D0, 1.D0, 0.D0, 8*.5D0, 2*0.D0, 4*1.D0, 76*0.D0/
02694 IDABS=ABS(IDHEP)
02695
02696
02697 IF (IDABS.LE.100) THEN
02698 PHOSPI=SPIN(IDABS)
02699 ELSE
02700
02701
02702 PHOSPI=(MOD(IDABS,10)-1.D0)/2.D0
02703
02704
02705 PHOSPI=MAX(PHOSPI,0.D0)
02706 ENDIF
02707 RETURN
02708 END
02709 SUBROUTINE PHOERR(IMES,TEXT,DATA)
02710
02711
02712
02713
02714
02715
02716
02717
02718
02719
02720
02721
02722
02723
02724
02725 IMPLICIT NONE
02726 DOUBLE PRECISION DATA
02727 INTEGER IMES,IERROR
02728 REAL*8 SDATA
02729 INTEGER PHLUN
02730 COMMON/PHOLUN/PHLUN
02731 INTEGER PHOMES
02732 PARAMETER (PHOMES=10)
02733 INTEGER STATUS
02734 COMMON/PHOSTA/STATUS(PHOMES)
02735 CHARACTER TEXT*(*)
02736 SAVE IERROR
02737
02738 LOGICAL ISEC
02739 SAVE ISEC
02740 DATA ISEC /.TRUE./
02741 DATA IERROR/ 0/
02742 IF (IMES.LE.PHOMES) STATUS(IMES)=STATUS(IMES)+1
02743
02744
02745 IF ((IMES.EQ. 6).AND.(STATUS(IMES).GE.2)) RETURN
02746 IF ((IMES.EQ.10).AND.(STATUS(IMES).GE.2)) RETURN
02747 SDATA=DATA
02748 WRITE(PHLUN,9000)
02749 WRITE(PHLUN,9120)
02750 GOTO (10,20,30,40,50,60,70,80,90,100),IMES
02751 WRITE(PHLUN,9130) IMES
02752 GOTO 120
02753 10 WRITE(PHLUN,9010) TEXT,INT(SDATA)
02754 GOTO 110
02755 20 WRITE(PHLUN,9020) TEXT,SDATA
02756 GOTO 110
02757 30 WRITE(PHLUN,9030) TEXT,SDATA
02758 GOTO 110
02759 40 WRITE(PHLUN,9040) TEXT
02760 GOTO 110
02761 50 WRITE(PHLUN,9050) TEXT,INT(SDATA)
02762 GOTO 110
02763 60 WRITE(PHLUN,9060) TEXT,SDATA
02764 GOTO 130
02765 70 WRITE(PHLUN,9070) TEXT,INT(SDATA)
02766 GOTO 110
02767 80 WRITE(PHLUN,9080) TEXT,INT(SDATA)
02768 GOTO 110
02769 90 WRITE(PHLUN,9090) TEXT,INT(SDATA)
02770 GOTO 110
02771 100 WRITE(PHLUN,9100) TEXT,SDATA
02772 GOTO 130
02773 110 CONTINUE
02774 WRITE(PHLUN,9140)
02775 WRITE(PHLUN,9120)
02776 WRITE(PHLUN,9000)
02777 IF (ISEC) THEN
02778 STOP
02779 ELSE
02780 GOTO 130
02781 ENDIF
02782 120 IERROR=IERROR+1
02783 IF (IERROR.GE.10) THEN
02784 WRITE(PHLUN,9150)
02785 WRITE(PHLUN,9120)
02786 WRITE(PHLUN,9000)
02787 IF (ISEC) THEN
02788 STOP
02789 ELSE
02790 GOTO 130
02791 ENDIF
02792 ENDIF
02793 130 WRITE(PHLUN,9120)
02794 WRITE(PHLUN,9000)
02795 RETURN
02796 9000 FORMAT(1H ,80('*'))
02797 9010 FORMAT(1H ,'* ',A,': Too many charged Particles, NCHARG =',I6,T81,
02798 &'*')
02799 9020 FORMAT(1H ,'* ',A,': Too much Bremsstrahlung required, PRSOFT = ',
02800 &F15.6,T81,'*')
02801 9030 FORMAT(1H ,'* ',A,': Combined Weight is exceeding 1., Weight = ',
02802 &F15.6,T81,'*')
02803 9040 FORMAT(1H ,'* ',A,
02804 &': Error in Rescaling charged and neutral Vectors',T81,'*')
02805 9050 FORMAT(1H ,'* ',A,
02806 &': Non matching charged Particle Pointer, NCHARG = ',I5,T81,'*')
02807 9060 FORMAT(1H ,'* ',A,
02808 &': Do you really work with a Particle of Spin: ',F4.1,' ?',T81,
02809 &'*')
02810 9070 FORMAT(1H ,'* ',A, ': Stack Length exceeded, NSTACK = ',I5 ,T81,
02811 &'*')
02812 9080 FORMAT(1H ,'* ',A,
02813 &': Random Number Generator Seed(1) out of Range: ',I8,T81,'*')
02814 9090 FORMAT(1H ,'* ',A,
02815 &': Random Number Generator Seed(2) out of Range: ',I8,T81,'*')
02816 9100 FORMAT(1H ,'* ',A,
02817 &': Available Phase Space below Cut-off: ',F15.6,' GeV/c^2',T81,
02818 &'*')
02819 9120 FORMAT(1H ,'*',T81,'*')
02820 9130 FORMAT(1H ,'* Funny Error Message: ',I4,' ! What to do ?',T81,'*')
02821 9140 FORMAT(1H ,'* Fatal Error Message, I stop this Run !',T81,'*')
02822 9150 FORMAT(1H ,'* 10 Error Messages generated, I stop this Run !',T81,
02823 &'*')
02824 END
02825 SUBROUTINE PHOREP
02826
02827
02828
02829
02830
02831
02832
02833
02834
02835
02836
02837
02838
02839
02840
02841 IMPLICIT NONE
02842 INTEGER PHLUN
02843 COMMON/PHOLUN/PHLUN
02844 INTEGER PHOMES
02845 PARAMETER (PHOMES=10)
02846 INTEGER STATUS
02847 COMMON/PHOSTA/STATUS(PHOMES)
02848 INTEGER I
02849 LOGICAL ERROR
02850 ERROR=.FALSE.
02851 WRITE(PHLUN,9000)
02852 WRITE(PHLUN,9010)
02853 WRITE(PHLUN,9020)
02854 WRITE(PHLUN,9030)
02855 WRITE(PHLUN,9040)
02856 WRITE(PHLUN,9030)
02857 WRITE(PHLUN,9020)
02858 DO 10 I=1,PHOMES
02859 IF (STATUS(I).EQ.0) GOTO 10
02860 IF ((I.EQ.6).OR.(I.EQ.10)) THEN
02861 WRITE(PHLUN,9050) I,STATUS(I)
02862 ELSE
02863 ERROR=.TRUE.
02864 WRITE(PHLUN,9060) I,STATUS(I)
02865 ENDIF
02866 10 CONTINUE
02867 IF (.NOT.ERROR) WRITE(PHLUN,9070)
02868 WRITE(PHLUN,9020)
02869 WRITE(PHLUN,9010)
02870 RETURN
02871 9000 FORMAT(1H1)
02872 9010 FORMAT(1H ,80('*'))
02873 9020 FORMAT(1H ,'*',T81,'*')
02874 9030 FORMAT(1H ,'*',26X,25('='),T81,'*')
02875 9040 FORMAT(1H ,'*',30X,'PHOTOS Run Summary',T81,'*')
02876 9050 FORMAT(1H ,'*',22X,'Warning #',I2,' occured',I6,' times',T81,'*')
02877 9060 FORMAT(1H ,'*',23X,'Error #',I2,' occured',I6,' times',T81,'*')
02878 9070 FORMAT(1H ,'*',16X,'PHOTOS Execution has successfully terminated',
02879 &T81,'*')
02880 END
02881 SUBROUTINE PHLUPA(IPOINT)
02882 IMPLICIT NONE
02883
02884
02885
02886
02887
02888
02889
02890
02891
02892
02893
02894
02895
02896
02897
02898
02899 INTEGER NMXPHO
02900 PARAMETER (NMXPHO=10000)
02901 INTEGER IDPHO,ISTPHO,JDAPHO,JMOPHO,NEVPHO,NPHO,I,J,IPOINT
02902 INTEGER IPOIN,IPOIN0,IPOINM,IEV
02903 INTEGER IOUT
02904 REAL*8 PPHO,VPHO,SUM
02905 COMMON/PHOEVT/NEVPHO,NPHO,ISTPHO(NMXPHO),IDPHO(NMXPHO),
02906 &JMOPHO(2,NMXPHO),JDAPHO(2,NMXPHO),PPHO(5,NMXPHO),VPHO(4,NMXPHO)
02907 COMMON /PHNUM/ IEV
02908 INTEGER PHLUN
02909 COMMON/PHOLUN/PHLUN
02910 DIMENSION SUM(5)
02911 DATA IPOIN0/ -5/
02912 COMMON /PHLUPY/ IPOIN,IPOINM
02913 SAVE IPOIN0
02914 IF (IPOIN0.LT.0) THEN
02915 IPOIN0=300 000
02916 IPOIN =IPOIN0
02917 IPOINM=300 001
02918 ENDIF
02919 IF (IPOINT.LE.IPOINM.OR.IPOINT.GE.IPOIN ) RETURN
02920 IOUT=56
02921 IF (IEV.LT.1000) THEN
02922 DO I=1,5
02923 SUM(I)=0.0D0
02924 ENDDO
02925 WRITE(PHLUN,*) 'EVENT NR=',IEV,
02926 $ 'WE ARE TESTING /PHOEVT/ at IPOINT=',IPOINT
02927 WRITE(PHLUN,10)
02928 I=1
02929 WRITE(PHLUN,20) IDPHO(I),PPHO(1,I),PPHO(2,I),PPHO(3,I),
02930 $ PPHO(4,I),PPHO(5,I),JDAPHO(1,I),JDAPHO(2,I)
02931 I=2
02932 WRITE(PHLUN,20) IDPHO(I),PPHO(1,I),PPHO(2,I),PPHO(3,I),
02933 $ PPHO(4,I),PPHO(5,I),JDAPHO(1,I),JDAPHO(2,I)
02934 WRITE(PHLUN,*) ' '
02935 DO I=3,NPHO
02936 WRITE(PHLUN,20) IDPHO(I),PPHO(1,I),PPHO(2,I),PPHO(3,I),
02937 $ PPHO(4,I),PPHO(5,I),JMOPHO(1,I),JMOPHO(2,I)
02938 DO J=1,4
02939 SUM(J)=SUM(J)+PPHO(J,I)
02940 ENDDO
02941 ENDDO
02942 SUM(5)=SQRT(ABS(SUM(4)**2-SUM(1)**2-SUM(2)**2-SUM(3)**2))
02943 WRITE(PHLUN,30) SUM
02944 10 FORMAT(1X,' ID ','p_x ','p_y ','p_z ',
02945 $ 'E ','m ',
02946 $ 'ID-MO_DA1','ID-MO DA2' )
02947 20 FORMAT(1X,I4,5(F14.9),2I9)
02948 30 FORMAT(1X,' SUM',5(F14.9))
02949 ENDIF
02950 END
02951
02952
02953
02954 FUNCTION IPHQRK(MODCOR)
02955 implicit none
02956
02957
02958
02959
02960
02961
02962
02963
02964
02965
02966
02967
02968
02969
02970
02971 INTEGER IPHQRK,MODCOR,MODOP
02972 INTEGER PHLUN
02973 COMMON/PHOLUN/PHLUN
02974
02975 SAVE MODOP
02976 DATA MODOP /0/
02977 IF (MODCOR.NE.0) THEN
02978
02979 MODOP=MODCOR
02980
02981 WRITE(PHLUN,*)
02982 $ 'Message from PHOTOS: IPHQRK(MODCOR):: (re)initialization'
02983 IF (MODOP.EQ.1) THEN
02984 WRITE(PHLUN,*)
02985 $ 'MODOP=1 -- blocks emission from light quarks: DEFAULT'
02986 ELSEIF (MODOP.EQ.2) THEN
02987 WRITE(PHLUN,*)
02988 $ 'MODOP=2 -- enables emission from light quarks: TEST '
02989 ELSE
02990 WRITE(PHLUN,*) 'IPHQRK wrong MODCOR=',MODCOR
02991 STOP
02992 ENDIF
02993 RETURN
02994 ENDIF
02995
02996 IF (MODOP.EQ.0.AND.MODCOR.EQ.0) THEN
02997 WRITE(PHLUN,*) 'IPHQRK lack of initialization'
02998 STOP
02999 ENDIF
03000 IPHQRK=MODOP
03001 END
03002
03003
03004 FUNCTION IPHEKL(MODCOR)
03005 implicit none
03006
03007
03008
03009
03010
03011
03012
03013
03014
03015
03016
03017
03018
03019
03020
03021 INTEGER IPHEKL,MODCOR,MODOP
03022 INTEGER PHLUN
03023 COMMON/PHOLUN/PHLUN
03024
03025 SAVE MODOP
03026 DATA MODOP /0/
03027
03028 IF (MODCOR.NE.0) THEN
03029
03030 MODOP=MODCOR
03031
03032 WRITE(PHLUN,*)
03033 $ 'Message from PHOTOS: IPHEKL(MODCOR):: (re)initialization'
03034 IF (MODOP.EQ.2) THEN
03035 WRITE(PHLUN,*)
03036 $ 'MODOP=2 -- blocks emission in pi0 to gamma e+e-: DEFAULT'
03037 WRITE(PHLUN,*)
03038 $ 'MODOP=2 -- blocks emission in Kl to gamma e+e-: DEFAULT'
03039 ELSEIF (MODOP.EQ.1) THEN
03040 WRITE(PHLUN,*)
03041 $ 'MODOP=1 -- enables emission in pi0 to gamma e+e- : TEST '
03042 WRITE(PHLUN,*)
03043 $ 'MODOP=1 -- enables emission in Kl to gamma e+e- : TEST '
03044 ELSE
03045 WRITE(PHLUN,*) 'IPHEKL wrong MODCOR=',MODCOR
03046 STOP
03047 ENDIF
03048 RETURN
03049 ENDIF
03050
03051 IF (MODOP.EQ.0.AND.MODCOR.EQ.0) THEN
03052 WRITE(PHLUN,*) 'IPHELK lack of initialization'
03053 STOP
03054 ENDIF
03055 IPHEKL=MODOP
03056 END
03057
03058 SUBROUTINE PHCORK(MODCOR)
03059 implicit none
03060
03061
03062
03063
03064
03065
03066
03067
03068
03069
03070
03071
03072
03073
03074
03075
03076
03077
03078
03079
03080
03081
03082 INTEGER NMXPHO
03083 PARAMETER (NMXPHO=10000)
03084
03085 REAL*8 M,P2,PX,PY,PZ,E,EN,MCUT,XMS
03086 INTEGER MODCOR,MODOP,I,IEV,IPRINT,K
03087 INTEGER IDPHO,ISTPHO,JDAPHO,JMOPHO,NEVPHO,NPHO
03088 REAL*8 PPHO,VPHO
03089 COMMON/PHOEVT/NEVPHO,NPHO,ISTPHO(NMXPHO),IDPHO(NMXPHO),
03090 &JMOPHO(2,NMXPHO),JDAPHO(2,NMXPHO),PPHO(5,NMXPHO),VPHO(4,NMXPHO)
03091
03092 INTEGER PHLUN
03093 COMMON/PHOLUN/PHLUN
03094
03095 COMMON /PHNUM/ IEV
03096 SAVE MODOP
03097 DATA MODOP /0/
03098 SAVE IPRINT
03099 DATA IPRINT /0/
03100 SAVE MCUT
03101 IF (MODCOR.NE.0) THEN
03102
03103 MODOP=MODCOR
03104
03105 WRITE(PHLUN,*) 'Message from PHCORK(MODCOR):: initialization'
03106 IF (MODOP.EQ.1) THEN
03107 WRITE(PHLUN,*) 'MODOP=1 -- no corrections on event: DEFAULT'
03108 ELSEIF (MODOP.EQ.2) THEN
03109 WRITE(PHLUN,*) 'MODOP=2 -- corrects Energy from mass'
03110 ELSEIF (MODOP.EQ.3) THEN
03111 WRITE(PHLUN,*) 'MODOP=3 -- corrects mass from Energy'
03112 ELSEIF (MODOP.EQ.4) THEN
03113 WRITE(PHLUN,*) 'MODOP=4 -- corrects Energy from mass to Mcut'
03114 WRITE(PHLUN,*) ' and mass from energy above Mcut '
03115 MCUT=0.4
03116 WRITE(PHLUN,*) 'Mcut=',MCUT,'GeV'
03117 ELSEIF (MODOP.EQ.5) THEN
03118 WRITE(PHLUN,*) 'MODOP=5 -- corrects Energy from mass+flow'
03119
03120 ELSE
03121 WRITE(PHLUN,*) 'PHCORK wrong MODCOR=',MODCOR
03122 STOP
03123 ENDIF
03124 RETURN
03125 ENDIF
03126
03127 IF (MODOP.EQ.0.AND.MODCOR.EQ.0) THEN
03128 WRITE(PHLUN,*) 'PHCORK lack of initialization'
03129 STOP
03130 ENDIF
03131
03132
03133
03134
03135
03136
03137 PX=0
03138 PY=0
03139 PZ=0
03140 E =0
03141
03142 IF (MODOP.EQ.1) THEN
03143
03144
03145 RETURN
03146 ELSEIF(MODOP.EQ.2) THEN
03147
03148
03149
03150
03151 DO I=3,NPHO
03152
03153 PX=PX+PPHO(1,I)
03154 PY=PY+PPHO(2,I)
03155 PZ=PZ+PPHO(3,I)
03156
03157 P2=PPHO(1,I)**2+PPHO(2,I)**2+PPHO(3,I)**2
03158
03159 EN=SQRT( PPHO(5,I)**2 + P2)
03160
03161 IF (IPRINT.EQ.1)
03162 & WRITE(PHLUN,*) "CORRECTING ENERGY OF ",I,":",
03163 & PPHO(4,I),"=>",EN
03164
03165 PPHO(4,I)=EN
03166 E = E+PPHO(4,I)
03167
03168 ENDDO
03169
03170 ELSEIF(MODOP.EQ.5) THEN
03171
03172
03173
03174
03175 DO I=3,NPHO
03176
03177 PX=PX+PPHO(1,I)
03178 PY=PY+PPHO(2,I)
03179 PZ=PZ+PPHO(3,I)
03180
03181 P2=PPHO(1,I)**2+PPHO(2,I)**2+PPHO(3,I)**2
03182
03183 EN=SQRT( PPHO(5,I)**2 + P2)
03184
03185 IF (IPRINT.EQ.1)
03186 & WRITE(PHLUN,*) "CORRECTING ENERGY OF ",I,":",
03187 & PPHO(4,I),"=>",EN
03188
03189 PPHO(4,I)=EN
03190 E = E+PPHO(4,I)
03191
03192 ENDDO
03193 DO K=1,4
03194 PPHO(K,1)=0d0
03195 DO I=3,NPHO
03196 PPHO(K,1)=PPHO(K,1)+PPHO(K,I)
03197 ENDDO
03198 ENDDO
03199 XMS=SQRT(PPHO(4,1)**2-PPHO(3,1)**2-PPHO(2,1)**2-PPHO(1,1)**2)
03200 PPHO(5,1)=XMS
03201 ELSEIF(MODOP.EQ.3) THEN
03202
03203
03204
03205
03206
03207 DO I=3,NPHO
03208
03209 PX=PX+PPHO(1,I)
03210 PY=PY+PPHO(2,I)
03211 PZ=PZ+PPHO(3,I)
03212 E = E+PPHO(4,I)
03213
03214 P2=PPHO(1,I)**2+PPHO(2,I)**2+PPHO(3,I)**2
03215
03216 M=SQRT(ABS( PPHO(4,I)**2 - P2))
03217
03218 IF (IPRINT.EQ.1)
03219 & WRITE(PHLUN,*) "CORRECTING MASS OF ",I,":",
03220 & PPHO(5,I),"=>",M
03221
03222 PPHO(5,I)=M
03223
03224 ENDDO
03225
03226
03227 ELSEIF(MODOP.EQ.4) THEN
03228
03229
03230
03231
03232
03233 DO I=3,NPHO
03234
03235 PX=PX+PPHO(1,I)
03236 PY=PY+PPHO(2,I)
03237 PZ=PZ+PPHO(3,I)
03238
03239 P2=PPHO(1,I)**2+PPHO(2,I)**2+PPHO(3,I)**2
03240
03241 M=SQRT(ABS( PPHO(4,I)**2 - P2))
03242
03243 IF (M.GT.MCUT) THEN
03244 IF (IPRINT.EQ.1)
03245 & WRITE(PHLUN,*) "CORRECTING MASS OF ",I,":",
03246 & PPHO(5,I),"=>",M
03247 PPHO(5,I)=M
03248 E = E+PPHO(4,I)
03249 ELSE
03250
03251 EN=SQRT( PPHO(5,I)**2 + P2)
03252
03253 IF (IPRINT.EQ.1)
03254 & WRITE(PHLUN,*) "CORRECTING ENERGY OF ",I,":",
03255 & PPHO(4,I),"=>",EN
03256
03257 PPHO(4,I)=EN
03258 E = E+PPHO(4,I)
03259 ENDIF
03260
03261 ENDDO
03262 ENDIF
03263
03264
03265 IF (IPRINT.EQ.1) THEN
03266 WRITE(PHLUN,*) "CORRECTING MOTHER"
03267 WRITE(PHLUN,*) "PX:",PPHO(1,1),"=>",PX-PPHO(1,2)
03268 WRITE(PHLUN,*) "PY:",PPHO(2,1),"=>",PY-PPHO(2,2)
03269 WRITE(PHLUN,*) "PZ:",PPHO(3,1),"=>",PZ-PPHO(3,2)
03270 WRITE(PHLUN,*) " E:",PPHO(4,1),"=>",E-PPHO(4,2)
03271 ENDIF
03272
03273 PPHO(1,1)=PX-PPHO(1,2)
03274 PPHO(2,1)=PY-PPHO(2,2)
03275 PPHO(3,1)=PZ-PPHO(3,2)
03276 PPHO(4,1)=E -PPHO(4,2)
03277
03278 P2=PPHO(1,1)**2+PPHO(2,1)**2+PPHO(3,1)**2
03279
03280 IF (PPHO(4,1)**2.GT.P2) THEN
03281 M=SQRT( PPHO(4,1)**2 - P2 )
03282 IF (IPRINT.EQ.1)
03283 & WRITE(PHLUN,*) " M:",PPHO(5,1),"=>",M
03284 PPHO(5,1)=M
03285 ENDIF
03286
03287 CALL PHLUPA(25)
03288
03289 END
03290
03291
03292
03293 FUNCTION PHINT(IDUM)
03294
03295
03296
03297
03298
03299
03300
03301
03302
03303
03304
03305
03306
03307
03308
03309
03310 IMPLICIT NONE
03311 REAL*8 PHINT,PHINT2
03312 INTEGER IDUM
03313 INTEGER NMXPHO
03314 PARAMETER (NMXPHO=10000)
03315 INTEGER IDPHO,ISTPHO,JDAPHO,JMOPHO,NEVPHO,NPHO
03316 REAL*8 PPHO,VPHO
03317 COMMON/PHOEVT/NEVPHO,NPHO,ISTPHO(NMXPHO),IDPHO(NMXPHO),
03318 &JMOPHO(2,NMXPHO),JDAPHO(2,NMXPHO),PPHO(5,NMXPHO),VPHO(4,NMXPHO)
03319 INTEGER I,K,L
03320 DOUBLE PRECISION EMU,MCHREN,BETA,COSTHG,MPASQR,XPH, XC1, XC2,XDENO
03321 DOUBLE PRECISION XNUM1,XNUM2
03322 DOUBLE PRECISION EPS1(4),EPS2(4),PH(4),PL(4)
03323 REAL*8 PHOCHA
03324
03325
03326
03327
03328 DO K=1,4
03329 PH(K)=PPHO(K,NPHO)
03330 EPS2(K)=1D0
03331 ENDDO
03332
03333 CALL PHOEPS(PH,EPS2,EPS1)
03334 CALL PHOEPS(PH,EPS1,EPS2)
03335
03336
03337 XNUM1=0D0
03338 XNUM2=0D0
03339 XDENO=0D0
03340
03341 DO K=JDAPHO(1,1),NPHO-1
03342
03343
03344 DO L=1,4
03345 PL(L)=PPHO(L,K)
03346 ENDDO
03347
03348
03349 XC1 = - PHOCHA(IDPHO(K)) *
03350 & ( PL(1)*EPS1(1) + PL(2)*EPS1(2) + PL(3)*EPS1(3) ) /
03351 & ( PH(4)*PL(4) - PH(1)*PL(1) - PH(2)*PL(2) - PH(3)*PL(3) )
03352
03353 XC2 = - PHOCHA(IDPHO(K)) *
03354 & ( PL(1)*EPS2(1) + PL(2)*EPS2(2) + PL(3)*EPS2(3) ) /
03355 & ( PH(4)*PL(4) - PH(1)*PL(1) - PH(2)*PL(2) - PH(3)*PL(3) )
03356
03357
03358
03359 XNUM1 = XNUM1+XC1
03360 XNUM2 = XNUM2+XC2
03361
03362 XDENO = XDENO + XC1**2 + XC2**2
03363
03364 ENDDO
03365
03366 PHINT=(XNUM1**2 + XNUM2**2) / XDENO
03367 PHINT2=PHINT
03368
03369 END
03370
03371
03372 SUBROUTINE PHOEPS (VEC1, VEC2, EPS)
03373
03374
03375
03376
03377
03378
03379
03380
03381
03382
03383
03384
03385
03386
03387
03388
03389
03390
03391 DOUBLE PRECISION VEC1(4), VEC2(4), EPS(4),XN
03392
03393 EPS(1)=VEC1(2)*VEC2(3) - VEC1(3)*VEC2(2)
03394 EPS(2)=VEC1(3)*VEC2(1) - VEC1(1)*VEC2(3)
03395 EPS(3)=VEC1(1)*VEC2(2) - VEC1(2)*VEC2(1)
03396 EPS(4)=0D0
03397
03398 XN=SQRT( EPS(1)**2 +EPS(2)**2 +EPS(3)**2)
03399
03400 EPS(1)=EPS(1)/XN
03401 EPS(2)=EPS(2)/XN
03402 EPS(3)=EPS(3)/XN
03403
03404
03405 END
03406 SUBROUTINE PHODMP
03407
03408
03409
03410
03411
03412
03413
03414
03415
03416
03417
03418
03419
03420
03421
03422 DOUBLE PRECISION SUMVEC(5)
03423 INTEGER I,J
03424 #include "../include/HEPEVT.h"
03425 INTEGER PHLUN
03426 COMMON/PHOLUN/PHLUN
03427 DO 10 I=1,5
03428 10 SUMVEC(I)=0.
03429
03430
03431 WRITE(PHLUN,9000)
03432 WRITE(PHLUN,9010) NEVHEP
03433 WRITE(PHLUN,9080)
03434 WRITE(PHLUN,9020)
03435 DO 30 I=1,NHEP
03436
03437
03438 IF (JDAHEP(1,I).EQ.0) THEN
03439 DO 20 J=1,4
03440 20 SUMVEC(J)=SUMVEC(J)+PHEP(J,I)
03441 IF (JMOHEP(2,I).EQ.0) THEN
03442 WRITE(PHLUN,9030) I,IDHEP(I),JMOHEP(1,I),(PHEP(J,I),J=1,5)
03443 ELSE
03444 WRITE(PHLUN,9040) I,IDHEP(I),JMOHEP(1,I),JMOHEP(2,I),(PHEP
03445 & (J,I),J=1,5)
03446 ENDIF
03447 ELSE
03448 IF (JMOHEP(2,I).EQ.0) THEN
03449 WRITE(PHLUN,9050) I,IDHEP(I),JMOHEP(1,I),JDAHEP(1,I),
03450 & JDAHEP(2,I),(PHEP(J,I),J=1,5)
03451 ELSE
03452 WRITE(PHLUN,9060) I,IDHEP(I),JMOHEP(1,I),JMOHEP(2,I),
03453 & JDAHEP(1,I),JDAHEP(2,I),(PHEP(J,I),J=1,5)
03454 ENDIF
03455 ENDIF
03456 30 CONTINUE
03457 SUMVEC(5)=SQRT(SUMVEC(4)**2-SUMVEC(1)**2-SUMVEC(2)**2-
03458 &SUMVEC(3)**2)
03459 WRITE(PHLUN,9070) (SUMVEC(J),J=1,5)
03460 RETURN
03461 9000 FORMAT(1H0,80('='))
03462 9010 FORMAT(1H ,29X,'Event No.:',I10)
03463 9020 FORMAT(1H0,1X,'Nr',3X,'Type',3X,'Parent(s)',2X,'Daughter(s)',6X,
03464 &'Px',7X,'Py',7X,'Pz',7X,'E',4X,'Inv. M.')
03465 9030 FORMAT(1H ,I4,I7,3X,I4,9X,'Stable',2X,5F9.2)
03466 9040 FORMAT(1H ,I4,I7,I4,' - ',I4,5X,'Stable',2X,5F9.2)
03467 9050 FORMAT(1H ,I4,I7,3X,I4,6X,I4,' - ',I4,5F9.2)
03468 9060 FORMAT(1H ,I4,I7,I4,' - ',I4,2X,I4,' - ',I4,5F9.2)
03469 9070 FORMAT(1H0,23X,'Vector Sum: ', 5F9.2)
03470 9080 FORMAT(1H0,6X,'Particle Parameters')
03471 END