110 INTEGER INIT,IDUM,IPHQRK,IPHEKL
115 IF (init.NE.0)
RETURN
164 parameter(d_h_nmxhep=2000)
165 #elif defined (CZTERYD)
166 parameter(d_h_nmxhep=4000)
168 parameter(d_h_nmxhep=2000)
170 parameter(d_h_nmxhep=2000)
171 #elif defined (DZIESD)
172 parameter(d_h_nmxhep=10000)
174 parameter(d_h_nmxhep=2000)
177 common/phoqed/qedrad(d_h_nmxhep)
181 common/phocop/alpha,xphcut
183 common/phpico/pi,twopi
184 INTEGER ISEED,I97,J97
185 real*8 uran,cran,cdran,cmran
186 common/phseed/iseed(2),i97,j97,uran(97),cran,cdran,cmran
190 common/phosta/status(phomes)
191 LOGICAL INTERF,ISEC,ITRE,IEXP,IFTOP,IFW
192 real*8 fint,fsec,expeps
193 COMMON /phokey/ fsec,fint,expeps,interf,isec,itre,iexp,iftop,ifw
199 IF (init.NE.0)
RETURN
215 alpha=0.00729735039d0
216 pi=3.14159265358979324d0
217 twopi=6.28318530717958648d0
299 INTEGER PHOVN1,PHOVN2
300 common/phover/phovn1,phovn2
303 LOGICAL INTERF,ISEC,ITRE,IEXP,IFTOP,IFW
304 real*8 fint,fsec,expeps
305 COMMON /phokey/ fsec,fint,expeps,interf,isec,itre,iexp,iftop,ifw
307 common/phocop/alpha,xphcut
320 WRITE(phlun,9040) iv1,iv2
322 iv2=(phovn2-iv1*10000)/100
323 iv3=phovn2-iv1*10000-iv2*100
324 WRITE(phlun,9050) iv1,iv2,iv3
333 WRITE(phlun,9064) interf,isec,itre,iexp,iftop,ifw,alpha,xphcut
335 IF (interf)
WRITE(phlun,9061)
336 IF (isec)
WRITE(phlun,9062)
337 IF (itre)
WRITE(phlun,9066)
338 IF (iexp)
WRITE(phlun,9067) expeps
339 IF (iftop)
WRITE(phlun,9063)
340 IF (ifw)
WRITE(phlun,9065)
346 9010
FORMAT(1h ,
'*',t81,
'*')
347 9020
FORMAT(1h ,80(
'*'))
348 9030
FORMAT(1h ,
'*',26x,26(
'='),t81,
'*')
349 9040
FORMAT(1h ,
'*',28x,
'PHOTOS, Version: ',i2,
'.',i2,t81,
'*')
350 9050
FORMAT(1h ,
'*',28x,
'Released at: ',i2,
'/',i2,
'/',i2,t81,
'*')
351 9060
FORMAT(1h ,
'*',18x,
'PHOTOS QED Corrections in Particle Decays',
353 9061
FORMAT(1h ,
'*',18x,
'option with interference is active ',
355 9062
FORMAT(1h ,
'*',18x,
'option with double photons is active ',
357 9066
FORMAT(1h ,
'*',18x,
'option with triple/quatric photons is active',
359 9067
FORMAT(1h ,
'*',18x,
'option with exponentiation is active EPSEXP=',
361 9063
FORMAT(1h ,
'*',18x,
'emision in t tbar production is active ',
363 9064
FORMAT(1h ,
'*',18x,
'Internal input parameters:',t81,
'*'
365 &,/, 1h ,
'*',18x,
'INTERF=',l2,
' ISEC=',l2,
' ITRE=',l2,
366 &
' IEXP=',l2,
' IFTOP=',l2,
368 &,/, 1h ,
'*',18x,
'ALPHA_QED=',f8.5,
' XPHCUT=',e8.3,t81,
'*')
369 9065
FORMAT(1h ,
'*',18x,
'correction wt in decay of W is active ',
371 9070
FORMAT(1h ,
'*',9x,
372 &
'Monte Carlo Program - by E. Barberio, B. van Eijk and Z. Was',
374 & 1h ,
'*',9x,
'Version 2.09 - by P. Golonka and Z.W.',t81,
'*')
375 9080
FORMAT( 1h ,
'*',9x,
' ',t81,
'*',/,
377 &
' WARNING (1): /HEPEVT/ is not anymore the standard common block'
379 & 1h ,
'*',9x,
' ',t81,
'*',/,
382 &
' PHOTOS expects /d_HEPEVT/ to have REAL*8 variabl. To change to'
383 & ,t81,
'*',/, 1h ,
'*',9x,
384 &
' REAL*4 modify its declaration in subr. PHOTOS_GET PHOTOS_SET:'
385 & ,t81,
'*',/, 1h ,
'*',9x,
386 &
' REAL*8 d_h_phep, d_h_vhep'
387 & ,t81,
'*',/, 1h ,
'*',9x,
388 &
' WARNING (2): check dims. of /d_hepevt/ /phoqed/ /ph_hepevt/.'
389 & ,t81,
'*',/, 1h ,
'*',9x,
390 &
' HERE: d_h_nmxhep=2000 and NMXHEP=10000'
391 #elif defined (CZTERYD)
392 &
' PHOTOS expects /HEPEVT/ to have REAL*8 variables. To change to'
393 & ,t81,
'*',/, 1h ,
'*',9x,
394 &
' REAL*4 modify its declaration in subr. PHOTOS_GET PHOTOS_SET:'
395 & ,t81,
'*',/, 1h ,
'*',9x,
396 &
' REAL*8 d_h_phep, d_h_vhep'
397 & ,t81,
'*',/, 1h ,
'*',9x,
398 &
' WARNING (2): check dims. of /hepevt/ /phoqed/ /ph_hepevt/.'
399 & ,t81,
'*',/, 1h ,
'*',9x,
400 &
' HERE: d_h_nmxhep=4000 and NMXHEP=10000'
402 &
' PHOTOS expects /HEPEVT/ to have REAL*4 variables. To change to'
403 & ,t81,
'*',/, 1h ,
'*',9x,
404 &
' REAL*8 modify its declaration in subr. PHOTOS_GET PHOTOS_SET:'
405 & ,t81,
'*',/, 1h ,
'*',9x,
406 &
' REAL*8 d_h_phep, d_h_vhep'
407 & ,t81,
'*',/, 1h ,
'*',9x,
408 &
' WARNING (2): check dims. of /hepevt/ /phoqed/ /ph_hepevt/.'
409 & ,t81,
'*',/, 1h ,
'*',9x,
410 &
' HERE: d_h_nmxhep=2000 and NMXHEP=10000'
412 &
' PHOTOS expects /HEPEVT/ to have REAL*8 variables. To change to'
413 & ,t81,
'*',/, 1h ,
'*',9x,
414 &
' REAL*4 modify its declaration in subr. PHOTOS_GET PHOTOS_SET:'
415 & ,t81,
'*',/, 1h ,
'*',9x,
416 &
' REAL*8 d_h_phep, d_h_vhep'
417 & ,t81,
'*',/, 1h ,
'*',9x,
418 &
' WARNING (2): check dims. of /hepevt/ /phoqed/ /ph_hepevt/.'
419 & ,t81,
'*',/, 1h ,
'*',9x,
420 &
' HERE: d_h_nmxhep=2000 and NMXHEP=10000'
421 #elif defined (DZIESD)
422 &
' PHOTOS expects /HEPEVT/ to have REAL*8 variables. To change to'
423 & ,t81,
'*',/, 1h ,
'*',9x,
424 &
' REAL*4 modify its declaration in subr. PHOTOS_GET PHOTOS_SET:'
425 & ,t81,
'*',/, 1h ,
'*',9x,
426 &
' REAL*8 d_h_phep, d_h_vhep'
427 & ,t81,
'*',/, 1h ,
'*',9x,
428 &
' WARNING (2): check dims. of /hepevt/ /phoqed/ /ph_hepevt/.'
429 & ,t81,
'*',/, 1h ,
'*',9x,
430 &
' HERE: d_h_nmxhep=10000 and NMXHEP=10000'
432 &
' PHOTOS expects /HEPEVT/ to have REAL*8 variables. To change to'
433 & ,t81,
'*',/, 1h ,
'*',9x,
434 &
' REAL*4 modify its declaration in subr. PHOTOS_GET PHOTOS_SET:'
435 & ,t81,
'*',/, 1h ,
'*',9x,
436 &
' REAL*8 d_h_phep, d_h_vhep'
437 & ,t81,
'*',/, 1h ,
'*',9x,
438 &
' WARNING (2): check dims. of /hepevt/ /phoqed/ /ph_hepevt/.'
439 & ,t81,
'*',/, 1h ,
'*',9x,
440 &
' HERE: d_h_nmxhep=2000 and NMXHEP=10000'
444 SUBROUTINE photos(ID)
462 COMMON /phlupy/ ipoin,ipoinm
469 IF (1.GT.ipoinm.AND.1.LT.ipoin )
THEN
470 WRITE(phlun,*)
'EVENT NR=',iev,
471 $
'WE ARE TESTING /HEPEVT/ at IPOINT=1 (input)'
477 IF (1.GT.ipoinm.AND.1.LT.ipoin )
THEN
478 WRITE(phlun,*)
'EVENT NR=',iev,
479 $
'WE ARE TESTING /HEPEVT/ at IPOINT=1 (output)'
485 SUBROUTINE photos_get
505 parameter(d_h_nmxhep=2000)
506 real*8 d_h_phep, d_h_vhep
507 INTEGER d_h_nevhep,d_h_nhep,d_h_isthep,d_h_idhep,d_h_jmohep,
510 #elif defined (CZTERYD)
511 parameter(d_h_nmxhep=4000)
512 real*8 d_h_phep, d_h_vhep
513 INTEGER d_h_nevhep,d_h_nhep,d_h_isthep,d_h_idhep,d_h_jmohep,
517 parameter(d_h_nmxhep=2000)
518 real*4 d_h_phep, d_h_vhep
519 INTEGER d_h_nevhep,d_h_nhep,d_h_isthep,d_h_idhep,d_h_jmohep,
523 parameter(d_h_nmxhep=2000)
524 real*8 d_h_phep, d_h_vhep
525 INTEGER d_h_nevhep,d_h_nhep,d_h_isthep,d_h_idhep,d_h_jmohep,
528 #elif defined (DZIESD)
529 parameter(d_h_nmxhep=10000)
530 real*8 d_h_phep, d_h_vhep
531 INTEGER d_h_nevhep,d_h_nhep,d_h_isthep,d_h_idhep,d_h_jmohep,
535 parameter(d_h_nmxhep=2000)
536 real*8 d_h_phep, d_h_vhep
537 INTEGER d_h_nevhep,d_h_nhep,d_h_isthep,d_h_idhep,d_h_jmohep,
543 $ d_h_isthep(d_h_nmxhep),
544 $ d_h_idhep(d_h_nmxhep),
545 $ d_h_jmohep(2,d_h_nmxhep),
546 $ d_h_jdahep(2,d_h_nmxhep),
547 $ d_h_phep(5,d_h_nmxhep),
548 $ d_h_vhep(4,d_h_nmxhep)
552 $ d_h_qedrad(d_h_nmxhep)
554 parameter(nmxhep=10000)
555 INTEGER IDHEP,ISTHEP,JDAHEP,JMOHEP,NEVHEP,NHEP
557 common/ph_hepevt/nevhep,nhep,isthep(nmxhep),idhep(nmxhep),
558 &jmohep(2,nmxhep),jdahep(2,nmxhep),phep(5,nmxhep),vhep(4,nmxhep)
560 common/ph_phoqed/qedrad(nmxhep)
565 isthep(k) =d_h_isthep(k)
566 idhep(k) =d_h_idhep(k)
567 jmohep(1,k) =d_h_jmohep(1,k)
568 jdahep(1,k) =d_h_jdahep(1,k)
569 jmohep(2,k) =d_h_jmohep(2,k)
570 jdahep(2,k) =d_h_jdahep(2,k)
572 phep(l,k) =d_h_phep(l,k)
573 vhep(l,k) =d_h_vhep(l,k)
575 phep(5,k) =d_h_phep(5,k)
576 qedrad(k) =d_h_qedrad(k)
581 SUBROUTINE photos_set
600 parameter(d_h_nmxhep=2000)
601 real*8 d_h_phep, d_h_vhep
602 INTEGER d_h_nevhep,d_h_nhep,d_h_isthep,d_h_idhep,d_h_jmohep,
605 #elif defined (CZTERYD)
606 parameter(d_h_nmxhep=4000)
607 real*8 d_h_phep, d_h_vhep
608 INTEGER d_h_nevhep,d_h_nhep,d_h_isthep,d_h_idhep,d_h_jmohep,
612 parameter(d_h_nmxhep=2000)
613 real*4 d_h_phep, d_h_vhep
614 INTEGER d_h_nevhep,d_h_nhep,d_h_isthep,d_h_idhep,d_h_jmohep,
618 parameter(d_h_nmxhep=2000)
619 real*8 d_h_phep, d_h_vhep
620 INTEGER d_h_nevhep,d_h_nhep,d_h_isthep,d_h_idhep,d_h_jmohep,
623 #elif defined (DZIESD)
624 parameter(d_h_nmxhep=10000)
625 real*8 d_h_phep, d_h_vhep
626 INTEGER d_h_nevhep,d_h_nhep,d_h_isthep,d_h_idhep,d_h_jmohep,
630 parameter(d_h_nmxhep=2000)
631 real*8 d_h_phep, d_h_vhep
632 INTEGER d_h_nevhep,d_h_nhep,d_h_isthep,d_h_idhep,d_h_jmohep,
638 $ d_h_isthep(d_h_nmxhep),
639 $ d_h_idhep(d_h_nmxhep),
640 $ d_h_jmohep(2,d_h_nmxhep),
641 $ d_h_jdahep(2,d_h_nmxhep),
642 $ d_h_phep(5,d_h_nmxhep),
643 $ d_h_vhep(4,d_h_nmxhep)
647 $ d_h_qedrad(d_h_nmxhep)
649 parameter(nmxhep=10000)
650 INTEGER IDHEP,ISTHEP,JDAHEP,JMOHEP,NEVHEP,NHEP
652 common/ph_hepevt/nevhep,nhep,isthep(nmxhep),idhep(nmxhep),
653 &jmohep(2,nmxhep),jdahep(2,nmxhep),phep(5,nmxhep),vhep(4,nmxhep)
655 common/ph_phoqed/qedrad(nmxhep)
661 d_h_isthep(k) =isthep(k)
662 d_h_idhep(k) =idhep(k)
663 d_h_jmohep(1,k) =jmohep(1,k)
664 d_h_jdahep(1,k) =jdahep(1,k)
665 d_h_jmohep(2,k) =jmohep(2,k)
666 d_h_jdahep(2,k) =jdahep(2,k)
668 d_h_phep(l,k) =phep(l,k)
669 d_h_vhep(l,k) =vhep(l,k)
671 d_h_phep(5,k) =phep(5,k)
672 d_h_qedrad(k) =qedrad(k)
675 SUBROUTINE photos_make(IPARR)
698 INTEGER IP,IPARR,IPPAR,I,J,K,L,NLAST
699 DOUBLE PRECISION DATA
700 INTEGER MOTHER,POSPHO
703 parameter(nmxhep=10000)
704 INTEGER IDHEP,ISTHEP,JDAHEP,JMOHEP,NEVHEP,NHEP
706 common/ph_hepevt/nevhep,nhep,isthep(nmxhep),idhep(nmxhep),
707 &jmohep(2,nmxhep),jdahep(2,nmxhep),phep(5,nmxhep),vhep(4,nmxhep)
709 common/ph_phoqed/qedrad(nmxhep)
711 parameter(nmxpho=10000)
712 INTEGER ISTACK(0:NMXPHO),NUMIT,NTRY,KK,LL,II,NA,FIRST,LAST
713 INTEGER FIRSTA,LASTA,IPP,IDA1,IDA2,MOTHER2,IDPHO,ISPHO
714 real*8 porig(5,nmxpho)
723 IF ((jdahep(1,ip).EQ.0).OR.(jmohep(1,jdahep(1,ip)).NE.ip))
RETURN
737 DO i=jdahep(1,ip),jdahep(2,ip)
738 IF (jdahep(1,i).NE.0.AND.jmohep(1,jdahep(1,i)).EQ.i)
THEN
740 IF (numit.GT.nmxpho)
THEN
742 CALL phoerr(7,
'PHOTOS',data)
747 IF(numit.GT.ntry)
THEN
756 first=jdahep(1,istack(kk))
757 last=jdahep(2,istack(kk))
760 porig(ll,ii)=phep(ll,first+ii-1)
764 CALL phtype(istack(kk))
772 IF(jmohep(1,ipp).EQ.istack(kk))
773 $
CALL phobos(ipp,porig(1,ii),phep(1,ipp),firsta,lasta)
779 IF (nhep.GT.nlast)
THEN
780 DO 160 i=nlast+1,nhep
784 pospho=jdahep(2,mother)+1
787 90 photon(j)=phep(j,i)
795 IF (pospho.NE.nhep)
THEN
799 DO 120 k=i,pospho+1,-1
800 isthep(k)=isthep(k-1)
801 qedrad(k)=qedrad(k-1)
804 jmohep(l,k)=jmohep(l,k-1)
805 100 jdahep(l,k)=jdahep(l,k-1)
807 110 phep(l,k)=phep(l,k-1)
809 120 vhep(l,k)=vhep(l,k-1)
814 IF ((jmohep(l,k).NE.0).AND.(jmohep(l,k).GE.
815 & pospho)) jmohep(l,k)=jmohep(l,k)+1
816 IF ((jdahep(l,k).NE.0).AND.(jdahep(l,k).GE.
817 & pospho)) jdahep(l,k)=jdahep(l,k)+1
822 140 phep(j,pospho)=photon(j)
826 jdahep(2,mother)=pospho
829 jmohep(1,pospho)=mother
830 jmohep(2,pospho)=mother2
831 jdahep(1,pospho)=ida1
832 jdahep(2,pospho)=ida2
836 150 vhep(j,pospho)=vhep(j,pospho-1)
841 SUBROUTINE phobos(IP,PBOOS1,PBOOS2,FIRST,LAST)
865 DOUBLE PRECISION BET1(3),BET2(3),GAM1,GAM2,PB,DATA
866 INTEGER I,J,FIRST,LAST,MAXSTA,NSTACK,IP
867 parameter(maxsta=10000)
868 INTEGER STACK(MAXSTA)
869 real*8 pboos1(5),pboos2(5)
871 parameter(nmxhep=10000)
872 INTEGER IDHEP,ISTHEP,JDAHEP,JMOHEP,NEVHEP,NHEP
874 common/ph_hepevt/nevhep,nhep,isthep(nmxhep),idhep(nmxhep),
875 &jmohep(2,nmxhep),jdahep(2,nmxhep),phep(5,nmxhep),vhep(4,nmxhep)
876 IF ((last.EQ.0).OR.(last.LT.first))
RETURN
879 bet1(j)=-pboos1(j)/pboos1(5)
880 10 bet2(j)=pboos2(j)/pboos2(5)
881 gam1=pboos1(4)/pboos1(5)
882 gam2=pboos2(4)/pboos2(5)
885 20
DO 50 i=first,last
886 pb=bet1(1)*phep(1,i)+bet1(2)*phep(2,i)+bet1(3)*phep(3,i)
887 IF (jmohep(1,i).EQ.ip)
THEN
889 30 phep(j,i)=phep(j,i)+bet1(j)*(phep(4,i)+pb/(gam1+1.d0))
890 phep(4,i)=gam1*phep(4,i)+pb
893 pb=bet2(1)*phep(1,i)+bet2(2)*phep(2,i)+bet2(3)*phep(3,i)
895 40 phep(j,i)=phep(j,i)+bet2(j)*(phep(4,i)+pb/(gam2+1.d0))
896 phep(4,i)=gam2*phep(4,i)+pb
897 IF (jdahep(1,i).NE.0)
THEN
901 IF (nstack.GT.maxsta)
THEN
903 CALL phoerr(7,
'PHOBOS',data)
909 IF (nstack.NE.0)
THEN
912 first=jdahep(1,stack(nstack))
913 last=jdahep(2,stack(nstack))
920 SUBROUTINE phoin(IP,BOOST,NHEP0)
941 parameter(nmxhep=10000)
942 INTEGER IDHEP,ISTHEP,JDAHEP,JMOHEP,NEVHEP,NHEP
944 common/ph_hepevt/nevhep,nhep,isthep(nmxhep),idhep(nmxhep),
945 &jmohep(2,nmxhep),jdahep(2,nmxhep),phep(5,nmxhep),vhep(4,nmxhep)
947 parameter(nmxpho=10000)
948 INTEGER IDPHO,ISTPHO,JDAPHO,JMOPHO,NEVPHO,NPHO
950 common/phoevt/nevpho,npho,istpho(nmxpho),idpho(nmxpho),
951 &jmopho(2,nmxpho),jdapho(2,nmxpho),ppho(5,nmxpho),vpho(4,nmxpho)
952 INTEGER IP,IP2,I,FIRST,LAST,LL,NA
955 DOUBLE PRECISION BET(3),GAM,PB
956 COMMON /phocms/ bet,gam
957 LOGICAL INTERF,ISEC,ITRE,IEXP,IFTOP,IFW
958 real*8 fint,fsec,expeps
959 COMMON /phokey/ fsec,fint,expeps,interf,isec,itre,iexp,iftop,ifw
964 npho=3+last-first+nhep-nhep0
969 jdapho(2,1)=3+last-first
974 ip2=jmohep(2,jdahep(1,ip))
975 IF((ip2.NE.0).AND.(ip2.NE.ip))
THEN
978 jdapho(2,2)=3+last-first
980 ppho(i,2)=phep(i,ip2)
990 idpho(3+ll)=idhep(first+ll)
991 jmopho(1,3+ll)=jmohep(1,first+ll)
992 IF (jmohep(1,first+ll).EQ.ip) jmopho(1,3+ll)=1
994 ppho(i,3+ll)=phep(i,first+ll)
997 IF (nhep.GT.nhep0)
THEN
1001 idpho(na+ll)=idhep(nhep0+ll)
1002 jmopho(1,na+ll)=jmohep(1,nhep0+ll)
1003 IF (jmohep(1,nhep0+ll).EQ.ip) jmopho(1,na+ll)=1
1005 ppho(i,na+ll)=phep(i,nhep0+ll)
1009 jdapho(2,1)=3+last-first+nhep-nhep0
1011 IF(idpho(npho).EQ.22)
CALL phlupa(100001)
1014 IF(idpho(npho).EQ.22)
CALL phlupa(100002)
1016 IF(iftop)
CALL photwo(0)
1030 IF ((abs(ppho(1,1)+abs(ppho(2,1))+abs(ppho(3,1))).GT.
1031 $ ppho(5,1)*1.d-8).AND.(ppho(5,1).NE.0))
THEN
1038 10 bet(j)=-ppho(j,1)/ppho(5,1)
1039 gam=ppho(4,1)/ppho(5,1)
1040 DO 30 i=jdapho(1,1),jdapho(2,1)
1041 pb=bet(1)*ppho(1,i)+bet(2)*ppho(2,i)+bet(3)*ppho(3,i)
1043 20 ppho(j,i)=ppho(j,i)+bet(j)*(ppho(4,i)+pb/(gam+1.d0))
1044 30 ppho(4,i)=gam*ppho(4,i)+pb
1047 pb=bet(1)*ppho(1,i)+bet(2)*ppho(2,i)+bet(3)*ppho(3,i)
1049 ppho(j,i)=ppho(j,i)+bet(j)*(ppho(4,i)+pb/(gam+1.d0))
1051 ppho(4,i)=gam*ppho(4,i)+pb
1054 IF(iftop)
CALL photwo(1)
1056 IF(idpho(npho).EQ.22)
CALL phlupa(10000)
1059 SUBROUTINE photwo(MODE)
1077 parameter(nmxpho=10000)
1078 INTEGER IDPHO,ISTPHO,JDAPHO,JMOPHO,NEVPHO,NPHO
1080 common/phoevt/nevpho,npho,istpho(nmxpho),idpho(nmxpho),
1081 &jmopho(2,nmxpho),jdapho(2,nmxpho),ppho(5,nmxpho),vpho(4,nmxpho)
1082 DOUBLE PRECISION BET(3),GAM
1083 COMMON /phocms/ bet,gam
1094 ifrad=(idpho(1).EQ.21).AND.(idpho(2).EQ.21)
1095 ifrad=ifrad.OR.(idpho(1).EQ.-idpho(2).AND.abs(idpho(1)).LE.6)
1097 & .AND.(abs(idpho(3)).EQ.6).AND.(abs(idpho(4)).EQ.6)
1098 mpasqr= (ppho(4,1)+ppho(4,2))**2-(ppho(3,1)+ppho(3,2))**2
1099 & -(ppho(2,1)+ppho(2,2))**2-(ppho(1,1)+ppho(1,2))**2
1100 ifrad=ifrad.AND.(mpasqr.GT.0.0d0)
1104 ppho(i,1)=ppho(i,1)+ppho(i,2)
1106 ppho(5,1)=sqrt(mpasqr)
1118 SUBROUTINE phoout(IP,BOOST,NHEP0)
1139 parameter(nmxhep=10000)
1140 INTEGER IDHEP,ISTHEP,JDAHEP,JMOHEP,NEVHEP,NHEP
1142 common/ph_hepevt/nevhep,nhep,isthep(nmxhep),idhep(nmxhep),
1143 &jmohep(2,nmxhep),jdahep(2,nmxhep),phep(5,nmxhep),vhep(4,nmxhep)
1145 parameter(nmxpho=10000)
1146 INTEGER IDPHO,ISTPHO,JDAPHO,JMOPHO,NEVPHO,NPHO
1148 common/phoevt/nevpho,npho,istpho(nmxpho),idpho(nmxpho),
1149 &jmopho(2,nmxpho),jdapho(2,nmxpho),ppho(5,nmxpho),vpho(4,nmxpho)
1150 INTEGER IP,LL,FIRST,LAST,I
1152 INTEGER NN,J,K,NHEP0,NA
1153 DOUBLE PRECISION BET(3),GAM,PB
1154 COMMON /phocms/ bet,gam
1155 IF(npho.EQ.nevpho)
RETURN
1159 DO 110 j=jdapho(1,1),jdapho(2,1)
1160 pb=-bet(1)*ppho(1,j)-bet(2)*ppho(2,j)-bet(3)*ppho(3,j)
1162 100 ppho(k,j)=ppho(k,j)-bet(k)*(ppho(4,j)+pb/(gam+1.d0))
1163 110 ppho(4,j)=gam*ppho(4,j)+pb
1166 pb=-bet(1)*ppho(1,nn)-bet(2)*ppho(2,nn)-bet(3)*ppho(3,nn)
1168 120 ppho(k,nn)=ppho(k,nn)-bet(k)*(ppho(4,nn)+pb/(gam+1.d0))
1169 ppho(4,nn)=gam*ppho(4,nn)+pb
1176 idhep(first+ll) = idpho(3+ll)
1178 phep(i,first+ll) = ppho(i,3+ll)
1184 idhep(nhep0+ll) = idpho(na+ll)
1185 isthep(nhep0+ll)=istpho(na+ll)
1186 jmohep(1,nhep0+ll)=ip
1187 jmohep(2,nhep0+ll)=jmohep(2,jdahep(1,ip))
1188 jdahep(1,nhep0+ll)=0
1189 jdahep(2,nhep0+ll)=0
1191 phep(i,nhep0+ll) = ppho(i,na+ll)
1194 nhep=nhep+npho-nevpho
1197 SUBROUTINE phochk(JFIRST)
1215 parameter(nmxpho=10000)
1216 INTEGER IDPHO,ISTPHO,JDAPHO,JMOPHO,NEVPHO,NPHO
1218 common/phoevt/nevpho,npho,istpho(nmxpho),idpho(nmxpho),
1219 &jmopho(2,nmxpho),jdapho(2,nmxpho),ppho(5,nmxpho),vpho(4,nmxpho)
1221 common/phoif/chkif(nmxpho)
1223 parameter(nmxhep=10000)
1225 common/ph_phoqed/qedrad(nmxhep)
1228 INTEGER IDABS,NLAST,I,IPPAR
1229 LOGICAL INTERF,ISEC,ITRE,IEXP,IFTOP,IFW,IFNPI0,IFKL
1230 real*8 fint,fsec,expeps
1231 COMMON /phokey/ fsec,fint,expeps,interf,isec,itre,iexp,iftop,ifw
1233 INTEGER IDENT,K,IQRK,IPHQRK,IEKL,IPHEKL
1236 & ( ((idabs.GT.9.OR.iqrk.NE.1).AND.(idabs.LE.40))
1237 & .OR.(idabs.GT.100) )
1238 & .AND.(idabs.NE.21)
1239 $ .AND.(idabs.NE.2101).AND.(idabs.NE.3101).AND.(idabs.NE.3201)
1240 & .AND.(idabs.NE.1103).AND.(idabs.NE.2103).AND.(idabs.NE.2203)
1241 & .AND.(idabs.NE.3103).AND.(idabs.NE.3203).AND.(idabs.NE.3303)
1252 ifnpi0= (idpho(1).NE.111)
1253 ifkl = ((idpho(1).EQ.130).AND.
1254 $ ((idpho(3).EQ.22).OR.(idpho(4).EQ.22).OR.
1255 $ (idpho(5).EQ.22)).AND.
1256 $ ((idpho(3).EQ.11).OR.(idpho(4).EQ.11).OR.
1257 $ (idpho(5).EQ.11)) )
1259 ifnpi0=(ifnpi0.AND.(.NOT.ifkl))
1262 idabs = abs(idpho(i))
1264 chkif(i)= f(idabs) .AND.f(abs(idpho(1)))
1265 & .AND. (idpho(2).EQ.0)
1266 IF(i.GT.2) chkif(i)=chkif(i).AND.qedrad(jfirst+i-ippar-2)
1274 DO k=jdapho(2,1),jdapho(1,1),-1
1275 IF(idpho(k).NE.22)
THEN
1281 ifrad=((idpho(1).EQ.21).AND.(idpho(2).EQ.21))
1282 & .OR. ((abs(idpho(1)).LE.6).AND.((idpho(2)).EQ.(-idpho(1))))
1284 & .AND.(abs(idpho(3)).EQ.6).AND.((idpho(4)).EQ.(-idpho(3)))
1289 IF(i.GT.2) chkif(i)=chkif(i).AND.qedrad(jfirst+i-ippar-2)
1297 DO k=jdapho(2,1),jdapho(1,1),-1
1298 IF(idpho(k).NE.22)
THEN
1304 ifrad=((abs(idpho(1)).EQ.6).AND.(idpho(2).EQ.0))
1306 & .AND.((abs(idpho(3)).EQ.24).AND.(abs(idpho(4)).EQ.5)
1307 & .OR.(abs(idpho(3)).EQ.5).AND.(abs(idpho(4)).EQ.24))
1312 IF(i.GT.2) chkif(i)=chkif(i).AND.qedrad(jfirst+i-ippar-2)
1319 SUBROUTINE phtype(ID)
1338 parameter(nmxhep=10000)
1339 INTEGER IDHEP,ISTHEP,JDAHEP,JMOHEP,NEVHEP,NHEP
1341 common/ph_hepevt/nevhep,nhep,isthep(nmxhep),idhep(nmxhep),
1342 &jmohep(2,nmxhep),jdahep(2,nmxhep),phep(5,nmxhep),vhep(4,nmxhep)
1343 LOGICAL INTERF,ISEC,ITRE,IEXP,IFTOP,IFW
1344 real*8 fint,fsec,expeps
1345 COMMON /phokey/ fsec,fint,expeps,interf,isec,itre,iexp,iftop,ifw
1349 real*8 pro,prsum,esu
1350 COMMON /phoexp/ pro(nx),nchan,expini
1354 real*8 rn,phoran,sum
1358 ifour=(.true.).AND.(itre)
1362 IF (jdahep(1,id).EQ.0)
RETURN
1375 CALL phomak(id,nhep0)
1389 IF(rn.LT.sum)
GOTO 100
1393 CALL phomak(id,nhep0)
1394 IF(sum.GT.1d0-expeps)
GOTO 100
1401 IF (rn.GE.23.d0/24d0)
THEN
1402 CALL phomak(id,nhep0)
1403 CALL phomak(id,nhep0)
1404 CALL phomak(id,nhep0)
1405 CALL phomak(id,nhep0)
1406 ELSEIF (rn.GE.17.d0/24d0)
THEN
1407 CALL phomak(id,nhep0)
1408 CALL phomak(id,nhep0)
1409 ELSEIF (rn.GE.9.d0/24d0)
THEN
1410 CALL phomak(id,nhep0)
1416 IF (rn.GE.5.d0/6d0)
THEN
1417 CALL phomak(id,nhep0)
1418 CALL phomak(id,nhep0)
1419 CALL phomak(id,nhep0)
1420 ELSEIF (rn.GE.2.d0/6d0)
THEN
1421 CALL phomak(id,nhep0)
1427 IF (rn.GE.0.5d0)
THEN
1428 CALL phomak(id,nhep0)
1429 CALL phomak(id,nhep0)
1434 CALL phomak(id,nhep0)
1440 SUBROUTINE phomak(IPPAR,NHEP0)
1461 DOUBLE PRECISION DATA
1463 INTEGER IP,IPPAR,NCHARG
1464 INTEGER WTDUM,IDUM,NHEP0
1465 INTEGER NCHARB,NEUDAU
1469 parameter(nmxhep=10000)
1470 INTEGER IDHEP,ISTHEP,JDAHEP,JMOHEP,NEVHEP,NHEP
1472 common/ph_hepevt/nevhep,nhep,isthep(nmxhep),idhep(nmxhep),
1473 &jmohep(2,nmxhep),jdahep(2,nmxhep),phep(5,nmxhep),vhep(4,nmxhep)
1474 LOGICAL INTERF,ISEC,ITRE,IEXP,IFTOP,IFW
1475 real*8 fint,fsec,expeps
1476 COMMON /phokey/ fsec,fint,expeps,interf,isec,itre,iexp,iftop,ifw
1482 CALL phoin(ip,boost,nhep0)
1483 CALL phochk(jdahep(1,ip))
1485 CALL phopre(1,wt,neudau,ncharb)
1487 IF (wt.EQ.0.0d0)
RETURN
1490 CALL phodo(1,ncharb,neudau)
1492 IF (interf) wt=wt*phint(idum) /fint
1493 IF (ifw)
CALL phobw(wt)
1495 IF (wt.GT.1.0d0)
CALL phoerr(3,
'WT_INT',data)
1498 CALL phoout(ip,boost,nhep0)
1502 FUNCTION phint1(IDUM)
1527 parameter(nmxpho=10000)
1528 INTEGER IDPHO,ISTPHO,JDAPHO,JMOPHO,NEVPHO,NPHO
1530 common/phoevt/nevpho,npho,istpho(nmxpho),idpho(nmxpho),
1531 &jmopho(2,nmxpho),jdapho(2,nmxpho),ppho(5,nmxpho),vpho(4,nmxpho)
1532 DOUBLE PRECISION MCHSQR,MNESQR
1534 common/phomom/mchsqr,mnesqr,pneutr(5)
1535 DOUBLE PRECISION COSTHG,SINTHG
1536 real*8 xphmax,xphoto
1537 common/phophs/xphmax,xphoto,costhg,sinthg
1538 real*8 mpasqr,xx,beta
1542 DO k=jdapho(2,1),jdapho(1,1),-1
1543 IF(idpho(k).NE.22)
THEN
1550 ifint= npho.GT.ident
1552 ifint= ifint.AND.(ident-jdapho(1,1)).EQ.1
1554 ifint= ifint.AND.idpho(jdapho(1,1)).EQ.-idpho(ident)
1556 ifint= ifint.AND.phocha(idpho(ident)).NE.0
1559 mpasqr = ppho(5,1)**2
1560 xx=4.d0*mchsqr/mpasqr*(1.d0-xphoto)/(1.d0-xphoto+(mchsqr-mnesqr)
1563 phint = 2d0/(1d0+costhg**2*beta**2)
1570 FUNCTION phint2(IDUM)
1591 real*8 phint,phint1,phint2
1595 parameter(nmxpho=10000)
1596 INTEGER IDPHO,ISTPHO,JDAPHO,JMOPHO,NEVPHO,NPHO
1598 common/phoevt/nevpho,npho,istpho(nmxpho),idpho(nmxpho),
1599 &jmopho(2,nmxpho),jdapho(2,nmxpho),ppho(5,nmxpho),vpho(4,nmxpho)
1600 DOUBLE PRECISION MCHSQR,MNESQR
1602 common/phomom/mchsqr,mnesqr,pneutr(5)
1603 DOUBLE PRECISION COSTHG,SINTHG
1604 real*8 xphmax,xphoto
1605 common/phophs/xphmax,xphoto,costhg,sinthg
1606 real*8 mpasqr,xx,beta,pq1(4),pq2(4),pphot(4)
1607 real*8 ss,pp2,pp,e1,e2,q1,q2,costhe
1611 DO k=jdapho(2,1),jdapho(1,1),-1
1612 IF(idpho(k).NE.22)
THEN
1619 ifint= npho.GT.ident
1621 ifint= ifint.AND.(ident-jdapho(1,1)).EQ.1
1625 ifint= ifint.AND.abs(phocha(idpho(ident))).GT.0.01d0
1628 $ abs(phocha(idpho(jdapho(1,1)))).gt.0.01d0
1631 mpasqr = ppho(5,1)**2
1632 xx=4.d0*mchsqr/mpasqr*(1.-xphoto)/(1.-xphoto+(mchsqr-mnesqr)/
1635 phint = 2d0/(1d0+costhg**2*beta**2)
1636 ss =mpasqr*(1.d0-xphoto)
1637 pp2=((ss-mchsqr-mnesqr)**2-4*mchsqr*mnesqr)/ss/4
1639 e1 =sqrt(pp2+mchsqr)
1640 e2 =sqrt(pp2+mnesqr)
1641 phint= (e1+e2)**2/((e2+costhg*pp)**2+(e1-costhg*pp)**2)
1643 q1=phocha(idpho(jdapho(1,1)))
1644 q2=phocha(idpho(ident))
1646 pq1(k)=ppho(k,jdapho(1,1))
1647 pq2(k)=ppho(k,jdapho(1,1)+1)
1648 pphot(k)=ppho(k,npho)
1650 costhe=(pphot(1)*pq1(1)+pphot(2)*pq1(2)+pphot(3)*pq1(3))
1651 costhe=costhe/sqrt(pq1(1)**2+pq1(2)**2+pq1(3)**2)
1652 costhe=costhe/sqrt(pphot(1)**2+pphot(2)**2+pphot(3)**2)
1658 IF (costhg*costhe.GT.0)
then
1660 phint= (q1*(e2+costhg*pp)-q2*(e1-costhg*pp))**2
1661 & /(q1**2*(e2+costhg*pp)**2+q2**2*(e1-costhg*pp)**2)
1664 phint= (q1*(e1-costhg*pp)-q2*(e2+costhg*pp))**2
1665 & /(q1**2*(e1-costhg*pp)**2+q2**2*(e2+costhg*pp)**2)
1675 SUBROUTINE phopre(IPARR,WT,NEUDAU,NCHARB)
1699 DOUBLE PRECISION MINMAS,MPASQR,MCHREN
1700 DOUBLE PRECISION BETA,EPS,DEL1,DEL2,DATA,BIGLOG
1701 real*8 phocha,phospi,phoran,phocor,massum
1702 INTEGER IP,IPARR,IPPAR,I,J,ME,NCHARG,NEUPOI,NLAST,THEDUM
1704 INTEGER NCHARB,NEUDAU
1707 parameter(nmxpho=10000)
1708 INTEGER IDPHO,ISTPHO,JDAPHO,JMOPHO,NEVPHO,NPHO
1710 common/phoevt/nevpho,npho,istpho(nmxpho),idpho(nmxpho),
1711 &jmopho(2,nmxpho),jdapho(2,nmxpho),ppho(5,nmxpho),vpho(4,nmxpho)
1713 common/phoif/chkif(nmxpho)
1714 INTEGER CHAPOI(NMXPHO)
1715 DOUBLE PRECISION MCHSQR,MNESQR
1717 common/phomom/mchsqr,mnesqr,pneutr(5)
1718 DOUBLE PRECISION COSTHG,SINTHG
1719 real*8 xphmax,xphoto
1720 common/phophs/xphmax,xphoto,costhg,sinthg
1722 common/phocop/alpha,xphcut
1724 real*8 probh,corwt,xf
1725 common/phopro/probh,corwt,xf,irep
1727 LOGICAL INTERF,ISEC,ITRE,IEXP,IFTOP,IFW
1728 real*8 fint,fsec,expeps
1729 COMMON /phokey/ fsec,fint,expeps,interf,isec,itre,iexp,iftop,ifw
1739 IF (jdapho(1,ip).EQ.0)
RETURN
1746 DO 20 i=jdapho(1,ip),jdapho(2,ip)
1751 IF (chkif(i-jdapho(1,ip)+3))
THEN
1752 IF (phocha(idpho(i)).NE.0)
THEN
1754 IF (ncharg.GT.nmxpho)
THEN
1756 CALL phoerr(1,
'PHOTOS',data)
1760 minmas=minmas+ppho(5,i)**2
1762 massum=massum+ppho(5,i)
1764 IF (ncharg.NE.0)
THEN
1767 IF ((ppho(5,ip)-massum)/ppho(5,ip).GT.2.d0*xphcut)
THEN
1771 IF (ncharg.GT.1)
CALL phooma(1,ncharg,chapoi)
1775 70 pneutr(j)=-ppho(j,chapoi(ncharg))
1776 pneutr(4)=ppho(5,ip)-ppho(4,chapoi(ncharg))
1779 mpasqr=ppho(5,ip)**2
1780 mchsqr=ppho(5,chapoi(ncharg))**2
1781 IF ((jdapho(2,ip)-jdapho(1,ip)).EQ.1)
THEN
1783 IF (neupoi.EQ.chapoi(ncharg)) neupoi=jdapho(2,ip)
1784 mnesqr=ppho(5,neupoi)**2
1785 pneutr(5)=ppho(5,neupoi)
1787 mnesqr=pneutr(4)**2-pneutr(1)**2-pneutr(2)**2-pneutr(3)**2
1788 mnesqr=max(mnesqr,minmas-mchsqr)
1789 pneutr(5)=sqrt(mnesqr)
1793 xphmax=(mpasqr-(pneutr(5)+ppho(5,chapoi(ncharg)))**2)/mpasqr
1796 CALL phoene(mpasqr,mchren,beta,biglog,idpho(chapoi(ncharg)))
1798 IF (xphoto.LT.-4d0)
THEN
1802 ELSEIF ((xphoto.LT.xphcut).OR.(xphoto.GT.xphmax))
THEN
1807 IF (ncharg.GT.0)
THEN
1815 eps=mchren/(1.d0+beta)
1818 del1=(2.d0-eps)*(eps/(2.d0-eps))**phoran(thedum)
1845 costhg=(1.d0-del1)/beta
1846 sinthg=sqrt(del1*del2-mchren)/beta
1852 me=2.d0*phospi(idpho(chapoi(ncharg)))+1.d0
1857 DO i=jdapho(1,ip),jdapho(2,ip)
1858 IF (i.NE.chapoi(ncharg))
THEN
1866 CALL phoerr(5,
'PHOKIN',data)
1868 ncharb=chapoi(ncharg)
1869 ncharb=ncharb-jdapho(1,ip)+3
1870 neudau=neudau-jdapho(1,ip)+3
1871 wt=phocor(mpasqr,mchren,me)*wgt
1874 data=ppho(5,ip)-massum
1875 CALL phoerr(10,
'PHOTOS',data)
1881 SUBROUTINE phooma(IFIRST,ILAST,POINTR)
1903 parameter(nmxpho=10000)
1904 INTEGER IDPHO,ISTPHO,JDAPHO,JMOPHO,NEVPHO,NPHO
1906 common/phoevt/nevpho,npho,istpho(nmxpho),idpho(nmxpho),
1907 &jmopho(2,nmxpho),jdapho(2,nmxpho),ppho(5,nmxpho),vpho(4,nmxpho)
1908 INTEGER IFIRST,ILAST,I,J,BUFPOI,POINTR(NMXPHO)
1909 real*8 bufmas,mass(nmxpho)
1910 IF (ifirst.EQ.ilast)
RETURN
1913 DO 10 i=ifirst,ilast
1914 10 mass(i)=ppho(5,pointr(i))
1917 DO 30 i=ifirst,ilast-1
1919 IF (mass(j).LE.mass(i))
GOTO 20
1930 SUBROUTINE phoene(MPASQR,MCHREN,BETA,BIGLOG,IDENT)
1953 DOUBLE PRECISION MPASQR,MCHREN,BIGLOG,BETA,DATA
1954 INTEGER IWT1,IRN,IWT2
1955 real*8 prsoft,prhard,phoran,phofac
1956 DOUBLE PRECISION MCHSQR,MNESQR
1959 real*8 phocha,prkill,rrr
1960 common/phomom/mchsqr,mnesqr,pneutr(5)
1961 DOUBLE PRECISION COSTHG,SINTHG
1962 real*8 xphmax,xphoto
1963 common/phophs/xphmax,xphoto,costhg,sinthg
1965 common/phocop/alpha,xphcut
1967 common/phpico/pi,twopi
1969 real*8 probh,corwt,xf
1970 common/phopro/probh,corwt,xf,irep
1971 LOGICAL INTERF,ISEC,ITRE,IEXP,IFTOP,IFW
1972 real*8 fint,fsec,expeps
1973 COMMON /phokey/ fsec,fint,expeps,interf,isec,itre,iexp,iftop,ifw
1978 COMMON /phoexp/ pro(nx),nchan,expini
1980 IF (xphmax.LE.xphcut)
THEN
1986 mchren=4.d0*mchsqr/mpasqr/(1.d0+mchsqr/mpasqr)**2
1987 beta=sqrt(1.d0-mchren)
2000 biglog=log(mpasqr/mchsqr*(1.d0+beta)**2/4.d0*
2001 & (1.d0+mchsqr/mpasqr)**2)
2002 prhard=alpha/pi*(1d0/beta*biglog)*
2003 &(log(xphmax/xphcut)-.75d0+xphcut/xphmax-.25d0*xphcut**2/xphmax**2)
2004 prhard=prhard*phocha(ident)**2*fsec*fint
2006 IF (irep.EQ.0) probh=0.d0
2011 pro(nchan)=prhard+0.05*(1.0+fint)
2024 prkill=pro(nchan)/prsum-prhard
2029 prhard=prhard*phofac(0)
2038 IF (prsoft.LT.-5.0d-8)
THEN
2040 CALL phoerr(2,
'PHOENE',data)
2043 IF (prsoft.LT.0.1d0)
THEN
2045 CALL phoerr(2,
'PHOENE',data)
2050 IF (rrr.LT.prsoft)
THEN
2054 IF (rrr.LT.prkill) xphoto=-5d0
2059 10 xphoto=exp(phoran(irn)*log(xphcut/xphmax))
2060 xphoto=xphoto*xphmax
2061 IF (phoran(iwt2).GT.((1.d0+(1.d0-xphoto/xphmax)**2)/2.d0))
2066 xf=4.d0*mchsqr*mpasqr/(mpasqr+mchsqr-mnesqr)**2
2069 FUNCTION phocor(MPASQR,MCHREN,ME)
2090 DOUBLE PRECISION MPASQR,MCHREN,BETA,XX,YY,DATA
2092 real*8 phocor,phofac,wt1,wt2,wt3
2093 DOUBLE PRECISION MCHSQR,MNESQR
2095 common/phomom/mchsqr,mnesqr,pneutr(5)
2096 DOUBLE PRECISION COSTHG,SINTHG
2097 real*8 xphmax,xphoto
2098 common/phophs/xphmax,xphoto,costhg,sinthg
2100 real*8 probh,corwt,xf
2101 common/phopro/probh,corwt,xf,irep
2104 xx=4.d0*mchsqr/mpasqr*(1.d0-xphoto)/(1.d0-xphoto+(mchsqr-mnesqr)/
2108 wt3=(1.d0-xphoto/xphmax)/((1.d0+(1.d0-xphoto/xphmax)**2)/2.d0)
2109 ELSEIF (me.EQ.2)
THEN
2110 yy=0.5d0*(1.d0-xphoto/xphmax+1.d0/(1.d0-xphoto/xphmax))
2112 ELSEIF ((me.EQ.3).OR.(me.EQ.4).OR.(me.EQ.5))
THEN
2114 wt3=(1.d0+(1.d0-xphoto/xphmax)**2-(xphoto/xphmax)**3)/
2115 & (1.d0+(1.d0-xphoto/xphmax)** 2)
2118 CALL phoerr(6,
'PHOCOR',data)
2123 wt1=(1.d0-costhg*sqrt(1.d0-mchren))/(1.d0-costhg*beta)
2124 wt2=(1.d0-xx/yy/(1.d0-beta**2*costhg**2))*(1.d0+costhg*beta)/2.d0
2128 IF (phocor.GT.1.d0)
THEN
2130 CALL phoerr(3,
'PHOCOR',data)
2134 FUNCTION phofac(MODE)
2165 real*8 phofac,ff,prx
2168 real*8 probh,corwt,xf
2169 common/phopro/probh,corwt,xf,irep
2170 LOGICAL INTERF,ISEC,ITRE,IEXP,IFTOP,IFW
2171 real*8 fint,fsec,expeps
2172 COMMON /phokey/ fsec,fint,expeps,interf,isec,itre,iexp,iftop,ifw
2174 DATA prx,ff/ 0.d0, 0.d0/
2179 IF (mode.EQ.-1)
THEN
2183 ELSEIF (mode.EQ.0)
THEN
2184 IF (irep.EQ.0) prx=1.d0
2185 prx=prx/(1.d0-probh)
2198 SUBROUTINE phobw(WT)
2224 parameter(nmxpho=10000)
2225 INTEGER IDPHO,ISTPHO,JDAPHO,JMOPHO,NEVPHO,NPHO
2227 common/phoevt/nevpho,npho,istpho(nmxpho),idpho(nmxpho),
2228 &jmopho(2,nmxpho),jdapho(2,nmxpho),ppho(5,nmxpho),vpho(4,nmxpho)
2230 DOUBLE PRECISION EMU,MCHREN,BETA,COSTHG,MPASQR,XPH
2232 IF (abs(idpho(1)).EQ.24.AND.
2233 $ abs(idpho(jdapho(1,1) )).GE.11.AND.
2234 $ abs(idpho(jdapho(1,1) )).LE.16.AND.
2235 $ abs(idpho(jdapho(1,1)+1)).GE.11.AND.
2236 $ abs(idpho(jdapho(1,1)+1)).LE.16 )
THEN
2239 $ abs(idpho(jdapho(1,1) )).EQ.11.OR.
2240 $ abs(idpho(jdapho(1,1) )).EQ.13.OR.
2241 $ abs(idpho(jdapho(1,1) )).EQ.15 )
THEN
2247 mchren=abs(ppho(4,i)**2-ppho(3,i)**2
2248 $ -ppho(2,i)**2-ppho(1,i)**2)
2249 beta=sqrt(1- mchren/ ppho(4,i)**2)
2250 costhg=(ppho(3,i)*ppho(3,npho)+ppho(2,i)*ppho(2,npho)
2251 $ +ppho(1,i)*ppho(1,npho))/
2252 $ sqrt(ppho(3,i)**2+ppho(2,i)**2+ppho(1,i)**2) /
2253 $ sqrt(ppho(3,npho)**2+ppho(2,npho)**2+ppho(1,npho)**2)
2256 wt=wt*(1-8*emu*xph*(1-costhg*beta)*
2257 $ (mchren+2*xph*sqrt(mpasqr))/
2258 $ mpasqr**2/(1-mchren/mpasqr)/(4-mchren/mpasqr))
2264 SUBROUTINE phodo(IP,NCHARB,NEUDAU)
2288 DOUBLE PRECISION PHOAN1,PHOAN2,ANGLE,FI1,FI3,FI4,FI5,TH1,TH3,TH4
2289 DOUBLE PRECISION PARNE,QNEW,QOLD,DATA
2290 INTEGER IP,FI3DUM,I,J,NEUDAU,FIRST,LAST
2292 real*8 ephoto,pmavir,photri
2293 real*8 gneut,phoran,ccosth,ssinth,pvec(4)
2295 parameter(nmxpho=10000)
2296 INTEGER IDPHO,ISTPHO,JDAPHO,JMOPHO,NEVPHO,NPHO
2298 common/phoevt/nevpho,npho,istpho(nmxpho),idpho(nmxpho),
2299 &jmopho(2,nmxpho),jdapho(2,nmxpho),ppho(5,nmxpho),vpho(4,nmxpho)
2300 DOUBLE PRECISION MCHSQR,MNESQR
2302 common/phomom/mchsqr,mnesqr,pneutr(5)
2303 DOUBLE PRECISION COSTHG,SINTHG
2304 real*8 xphmax,xphoto
2305 common/phophs/xphmax,xphoto,costhg,sinthg
2307 common/phpico/pi,twopi
2309 ephoto=xphoto*ppho(5,ip)/2.d0
2310 pmavir=sqrt(ppho(5,ip)*(ppho(5,ip)-2.d0*ephoto))
2313 fi1=phoan1(pneutr(1),pneutr(2))
2317 th1=phoan2(pneutr(3),sqrt(pneutr(1)**2+pneutr(2)**2))
2318 CALL phoro3(-fi1,pneutr(1))
2319 CALL phoro2(-th1,pneutr(1))
2327 qnew=photri(pmavir,pneutr(5),ppho(5,ncharb))
2329 gneut=(qnew**2+qold**2+mnesqr)/(qnew*qold+sqrt((qnew**2+mnesqr)*
2331 IF (gneut.LT.1.d0)
THEN
2333 CALL phoerr(4,
'PHOKIN',data)
2335 parne=gneut-sqrt(max(gneut**2-1.0d0,0.d0))
2338 CALL phobo3(parne,pneutr)
2349 ppho(4,npho)=ephoto*ppho(5,ip)/pmavir
2354 th3=phoan2(ccosth,ssinth)
2355 fi3=twopi*phoran(fi3dum)
2356 ppho(1,npho)=ppho(4,npho)*sinthg*cos(fi3)
2357 ppho(2,npho)=ppho(4,npho)*sinthg*sin(fi3)
2360 ppho(3,npho)=-ppho(4,npho)*costhg
2364 CALL phoro3(-fi3,pneutr(1))
2365 CALL phoro3(-fi3,ppho(1,npho))
2366 CALL phoro2(-th3,pneutr(1))
2367 CALL phoro2(-th3,ppho(1,npho))
2368 angle=ephoto/ppho(4,npho)
2371 CALL phobo3(angle,pneutr(1))
2372 CALL phobo3(angle,ppho(1,npho))
2375 fi4=phoan1(pneutr(1),pneutr(2))
2376 th4=phoan2(pneutr(3),sqrt(pneutr(1)**2+pneutr(2)**2))
2377 CALL phoro3(fi4,pneutr(1))
2378 CALL phoro3(fi4,ppho(1,npho))
2383 CALL phoro3(-fi3,pvec)
2384 CALL phoro2(-th3,pvec)
2385 CALL phobo3(angle,pvec)
2386 CALL phoro3(fi4,pvec)
2387 CALL phoro2(-th4,pneutr)
2388 CALL phoro2(-th4,ppho(1,npho))
2389 CALL phoro2(-th4,pvec)
2390 fi5=phoan1(pvec(1),pvec(2))
2393 CALL phoro3(-fi5,pneutr)
2394 CALL phoro3(-fi5,ppho(1,npho))
2395 CALL phoro2(th1,pneutr(1))
2396 CALL phoro2(th1,ppho(1,npho))
2397 CALL phoro3(fi1,pneutr)
2398 CALL phoro3(fi1,ppho(1,npho))
2400 IF ((jdapho(2,ip)-jdapho(1,ip)).GT.1)
THEN
2406 IF (i.NE.ncharb.AND.(jmopho(1,i).EQ.ip))
THEN
2409 CALL phoro3(-fi1,ppho(1,i))
2410 CALL phoro2(-th1,ppho(1,i))
2413 CALL phobo3(parne,ppho(1,i))
2416 CALL phoro3(-fi3,ppho(1,i))
2417 CALL phoro2(-th3,ppho(1,i))
2420 CALL phobo3(angle,ppho(1,i))
2423 CALL phoro3(fi4,ppho(1,i))
2424 CALL phoro2(-th4,ppho(1,i))
2427 CALL phoro3(-fi5,ppho(1,i))
2428 CALL phoro2(th1,ppho(1,i))
2429 CALL phoro3(fi1,ppho(1,i))
2436 80 ppho(j,neudau)=pneutr(j)
2441 90 ppho(j,ncharb)=-(ppho(j,npho)+pneutr(j))
2442 ppho(4,ncharb)=ppho(5,ip)-(ppho(4,npho)+pneutr(4))
2445 FUNCTION photri(A,B,C)
2462 DOUBLE PRECISION DA,DB,DC,DAPB,DAMB,DTRIAN
2469 dtrian=sqrt((damb-dc)*(dapb+dc)*(damb+dc)*(dapb-dc))
2470 photri=dtrian/(da+da)
2473 FUNCTION phoan1(X,Y)
2489 DOUBLE PRECISION PHOAN1
2492 common/phpico/pi,twopi
2493 IF (abs(y).LT.abs(x))
THEN
2494 phoan1=atan(abs(y/x))
2495 IF (x.LE.0.d0) phoan1=pi-phoan1
2497 phoan1=acos(x/sqrt(x**2+y**2))
2499 IF (y.LT.0.d0) phoan1=twopi-phoan1
2502 FUNCTION phoan2(X,Y)
2518 DOUBLE PRECISION PHOAN2
2521 common/phpico/pi,twopi
2522 IF (abs(y).LT.abs(x))
THEN
2523 phoan2=atan(abs(y/x))
2524 IF (x.LE.0.d0) phoan2=pi-phoan2
2526 phoan2=acos(x/sqrt(x**2+y**2))
2530 SUBROUTINE phobo3(ANGLE,PVEC)
2547 DOUBLE PRECISION QPL,QMI,ANGLE
2549 qpl=(pvec(4)+pvec(3))*angle
2550 qmi=(pvec(4)-pvec(3))/angle
2551 pvec(3)=(qpl-qmi)/2.d0
2552 pvec(4)=(qpl+qmi)/2.d0
2555 SUBROUTINE phoro2(ANGLE,PVEC)
2572 DOUBLE PRECISION CS,SN,ANGLE
2574 cs=cos(angle)*pvec(1)+sin(angle)*pvec(3)
2575 sn=-sin(angle)*pvec(1)+cos(angle)*pvec(3)
2580 SUBROUTINE phoro3(ANGLE,PVEC)
2597 DOUBLE PRECISION CS,SN,ANGLE
2599 cs=cos(angle)*pvec(1)-sin(angle)*pvec(2)
2600 sn=sin(angle)*pvec(1)+cos(angle)*pvec(2)
2605 #include "../randg/photos-random.h"
2606 FUNCTION phocha(IDHEP)
2626 INTEGER IDHEP,IDABS,Q1,Q2,Q3
2630 real*8 charge(0:100)
2632 &-0.3333333333d0, 0.6666666667d0, -0.3333333333d0, 0.6666666667d0,
2633 &-0.3333333333d0, 0.6666666667d0, -0.3333333333d0, 0.6666666667d0,
2634 & 2*0.d0, -1.d0, 0.d0, -1.d0, 0.d0, -1.d0, 0.d0, -1.d0, 6*0.d0,
2635 & 1.d0, 12*0.d0, 1.d0, 63*0.d0/
2637 IF (idabs.LE.100)
THEN
2640 phocha = charge(idabs)
2644 q3=mod(idabs/1000,10)
2645 q2=mod(idabs/100,10)
2650 IF(mod(q2,2).EQ.0)
THEN
2651 phocha=charge(q2)-charge(q1)
2653 phocha=charge(q1)-charge(q2)
2658 phocha=charge(q1)+charge(q2)+charge(q3)
2663 IF (idhep.LT.0.d0) phocha=-phocha
2664 IF (phocha**2.lt.1d-6) phocha=0.d0
2667 FUNCTION phospi(IDHEP)
2693 DATA spin/ 8*.5d0, 1.d0, 0.d0, 8*.5d0, 2*0.d0, 4*1.d0, 76*0.d0/
2697 IF (idabs.LE.100)
THEN
2702 phospi=(mod(idabs,10)-1.d0)/2.d0
2705 phospi=max(phospi,0.d0)
2709 SUBROUTINE phoerr(IMES,TEXT,DATA)
2726 DOUBLE PRECISION DATA
2732 parameter(phomes=10)
2734 common/phosta/status(phomes)
2742 IF (imes.LE.phomes) status(imes)=status(imes)+1
2745 IF ((imes.EQ. 6).AND.(status(imes).GE.2))
RETURN
2746 IF ((imes.EQ.10).AND.(status(imes).GE.2))
RETURN
2750 GOTO (10,20,30,40,50,60,70,80,90,100),imes
2751 WRITE(phlun,9130) imes
2753 10
WRITE(phlun,9010) text,int(sdata)
2755 20
WRITE(phlun,9020) text,sdata
2757 30
WRITE(phlun,9030) text,sdata
2759 40
WRITE(phlun,9040) text
2761 50
WRITE(phlun,9050) text,int(sdata)
2763 60
WRITE(phlun,9060) text,sdata
2765 70
WRITE(phlun,9070) text,int(sdata)
2767 80
WRITE(phlun,9080) text,int(sdata)
2769 90
WRITE(phlun,9090) text,int(sdata)
2771 100
WRITE(phlun,9100) text,sdata
2783 IF (ierror.GE.10)
THEN
2793 130
WRITE(phlun,9120)
2796 9000
FORMAT(1h ,80(
'*'))
2797 9010
FORMAT(1h ,
'* ',a,
': Too many charged Particles, NCHARG =',i6,t81,
2799 9020
FORMAT(1h ,
'* ',a,
': Too much Bremsstrahlung required, PRSOFT = ',
2801 9030
FORMAT(1h ,
'* ',a,
': Combined Weight is exceeding 1., Weight = ',
2803 9040
FORMAT(1h ,
'* ',a,
2804 &
': Error in Rescaling charged and neutral Vectors',t81,
'*')
2805 9050
FORMAT(1h ,
'* ',a,
2806 &
': Non matching charged Particle Pointer, NCHARG = ',i5,t81,
'*')
2807 9060
FORMAT(1h ,
'* ',a,
2808 &
': Do you really work with a Particle of Spin: ',f4.1,
' ?',t81,
2810 9070
FORMAT(1h ,
'* ',a,
': Stack Length exceeded, NSTACK = ',i5 ,t81,
2812 9080
FORMAT(1h ,
'* ',a,
2813 &
': Random Number Generator Seed(1) out of Range: ',i8,t81,
'*')
2814 9090
FORMAT(1h ,
'* ',a,
2815 &
': Random Number Generator Seed(2) out of Range: ',i8,t81,
'*')
2816 9100
FORMAT(1h ,
'* ',a,
2817 &
': Available Phase Space below Cut-off: ',f15.6,
' GeV/c^2',t81,
2819 9120
FORMAT(1h ,
'*',t81,
'*')
2820 9130
FORMAT(1h ,
'* Funny Error Message: ',i4,
' ! What to do ?',t81,
'*')
2821 9140
FORMAT(1h ,
'* Fatal Error Message, I stop this Run !',t81,
'*')
2822 9150
FORMAT(1h ,
'* 10 Error Messages generated, I stop this Run !',t81,
2845 parameter(phomes=10)
2847 common/phosta/status(phomes)
2859 IF (status(i).EQ.0)
GOTO 10
2860 IF ((i.EQ.6).OR.(i.EQ.10))
THEN
2861 WRITE(phlun,9050) i,status(i)
2864 WRITE(phlun,9060) i,status(i)
2867 IF (.NOT.error)
WRITE(phlun,9070)
2872 9010
FORMAT(1h ,80(
'*'))
2873 9020
FORMAT(1h ,
'*',t81,
'*')
2874 9030
FORMAT(1h ,
'*',26x,25(
'='),t81,
'*')
2875 9040
FORMAT(1h ,
'*',30x,
'PHOTOS Run Summary',t81,
'*')
2876 9050
FORMAT(1h ,
'*',22x,
'Warning #',i2,
' occured',i6,
' times',t81,
'*')
2877 9060
FORMAT(1h ,
'*',23x,
'Error #',i2,
' occured',i6,
' times',t81,
'*')
2878 9070
FORMAT(1h ,
'*',16x,
'PHOTOS Execution has successfully terminated',
2881 SUBROUTINE phlupa(IPOINT)
2900 parameter(nmxpho=10000)
2901 INTEGER IDPHO,ISTPHO,JDAPHO,JMOPHO,NEVPHO,NPHO,I,J,IPOINT
2902 INTEGER IPOIN,IPOIN0,IPOINM,IEV
2904 real*8 ppho,vpho,sum
2905 common/phoevt/nevpho,npho,istpho(nmxpho),idpho(nmxpho),
2906 &jmopho(2,nmxpho),jdapho(2,nmxpho),ppho(5,nmxpho),vpho(4,nmxpho)
2912 COMMON /phlupy/ ipoin,ipoinm
2914 IF (ipoin0.LT.0)
THEN
2919 IF (ipoint.LE.ipoinm.OR.ipoint.GE.ipoin )
RETURN
2921 IF (iev.LT.1000)
THEN
2925 WRITE(phlun,*)
'EVENT NR=',iev,
2926 $
'WE ARE TESTING /PHOEVT/ at IPOINT=',ipoint
2929 WRITE(phlun,20) idpho(i),ppho(1,i),ppho(2,i),ppho(3,i),
2930 $ ppho(4,i),ppho(5,i),jdapho(1,i),jdapho(2,i)
2932 WRITE(phlun,20) idpho(i),ppho(1,i),ppho(2,i),ppho(3,i),
2933 $ ppho(4,i),ppho(5,i),jdapho(1,i),jdapho(2,i)
2936 WRITE(phlun,20) idpho(i),ppho(1,i),ppho(2,i),ppho(3,i),
2937 $ ppho(4,i),ppho(5,i),jmopho(1,i),jmopho(2,i)
2939 sum(j)=sum(j)+ppho(j,i)
2942 sum(5)=sqrt(abs(sum(4)**2-sum(1)**2-sum(2)**2-sum(3)**2))
2944 10
FORMAT(1x,
' ID ',
'p_x ',
'p_y ',
'p_z ',
2946 $
'ID-MO_DA1',
'ID-MO DA2' )
2947 20
FORMAT(1x,i4,5(f14.9),2i9)
2948 30
FORMAT(1x,
' SUM',5(f14.9))
2954 FUNCTION iphqrk(MODCOR)
2971 INTEGER IPHQRK,MODCOR,MODOP
2977 IF (modcor.NE.0)
THEN
2982 $
'Message from PHOTOS: IPHQRK(MODCOR):: (re)initialization'
2983 IF (modop.EQ.1)
THEN
2985 $
'MODOP=1 -- blocks emission from light quarks: DEFAULT'
2986 ELSEIF (modop.EQ.2)
THEN
2988 $
'MODOP=2 -- enables emission from light quarks: TEST '
2990 WRITE(phlun,*)
'IPHQRK wrong MODCOR=',modcor
2996 IF (modop.EQ.0.AND.modcor.EQ.0)
THEN
2997 WRITE(phlun,*)
'IPHQRK lack of initialization'
3004 FUNCTION iphekl(MODCOR)
3021 INTEGER IPHEKL,MODCOR,MODOP
3028 IF (modcor.NE.0)
THEN
3033 $
'Message from PHOTOS: IPHEKL(MODCOR):: (re)initialization'
3034 IF (modop.EQ.2)
THEN
3036 $
'MODOP=2 -- blocks emission in pi0 to gamma e+e-: DEFAULT'
3038 $
'MODOP=2 -- blocks emission in Kl to gamma e+e-: DEFAULT'
3039 ELSEIF (modop.EQ.1)
THEN
3041 $
'MODOP=1 -- enables emission in pi0 to gamma e+e- : TEST '
3043 $
'MODOP=1 -- enables emission in Kl to gamma e+e- : TEST '
3045 WRITE(phlun,*)
'IPHEKL wrong MODCOR=',modcor
3051 IF (modop.EQ.0.AND.modcor.EQ.0)
THEN
3052 WRITE(phlun,*)
'IPHELK lack of initialization'
3058 SUBROUTINE phcork(MODCOR)
3083 parameter(nmxpho=10000)
3085 real*8 m,p2,px,py,pz,e,en,mcut,xms
3086 INTEGER MODCOR,MODOP,I,IEV,IPRINT,K
3087 INTEGER IDPHO,ISTPHO,JDAPHO,JMOPHO,NEVPHO,NPHO
3089 common/phoevt/nevpho,npho,istpho(nmxpho),idpho(nmxpho),
3090 &jmopho(2,nmxpho),jdapho(2,nmxpho),ppho(5,nmxpho),vpho(4,nmxpho)
3101 IF (modcor.NE.0)
THEN
3105 WRITE(phlun,*)
'Message from PHCORK(MODCOR):: initialization'
3106 IF (modop.EQ.1)
THEN
3107 WRITE(phlun,*)
'MODOP=1 -- no corrections on event: DEFAULT'
3108 ELSEIF (modop.EQ.2)
THEN
3109 WRITE(phlun,*)
'MODOP=2 -- corrects Energy from mass'
3110 ELSEIF (modop.EQ.3)
THEN
3111 WRITE(phlun,*)
'MODOP=3 -- corrects mass from Energy'
3112 ELSEIF (modop.EQ.4)
THEN
3113 WRITE(phlun,*)
'MODOP=4 -- corrects Energy from mass to Mcut'
3114 WRITE(phlun,*)
' and mass from energy above Mcut '
3116 WRITE(phlun,*)
'Mcut=',mcut,
'GeV'
3117 ELSEIF (modop.EQ.5)
THEN
3118 WRITE(phlun,*)
'MODOP=5 -- corrects Energy from mass+flow'
3121 WRITE(phlun,*)
'PHCORK wrong MODCOR=',modcor
3127 IF (modop.EQ.0.AND.modcor.EQ.0)
THEN
3128 WRITE(phlun,*)
'PHCORK lack of initialization'
3142 IF (modop.EQ.1)
THEN
3146 ELSEIF(modop.EQ.2)
THEN
3157 p2=ppho(1,i)**2+ppho(2,i)**2+ppho(3,i)**2
3159 en=sqrt( ppho(5,i)**2 + p2)
3162 &
WRITE(phlun,*)
"CORRECTING ENERGY OF ",i,
":",
3170 ELSEIF(modop.EQ.5)
THEN
3181 p2=ppho(1,i)**2+ppho(2,i)**2+ppho(3,i)**2
3183 en=sqrt( ppho(5,i)**2 + p2)
3186 &
WRITE(phlun,*)
"CORRECTING ENERGY OF ",i,
":",
3196 ppho(k,1)=ppho(k,1)+ppho(k,i)
3199 xms=sqrt(ppho(4,1)**2-ppho(3,1)**2-ppho(2,1)**2-ppho(1,1)**2)
3201 ELSEIF(modop.EQ.3)
THEN
3214 p2=ppho(1,i)**2+ppho(2,i)**2+ppho(3,i)**2
3216 m=sqrt(abs( ppho(4,i)**2 - p2))
3219 &
WRITE(phlun,*)
"CORRECTING MASS OF ",i,
":",
3227 ELSEIF(modop.EQ.4)
THEN
3239 p2=ppho(1,i)**2+ppho(2,i)**2+ppho(3,i)**2
3241 m=sqrt(abs( ppho(4,i)**2 - p2))
3245 &
WRITE(phlun,*)
"CORRECTING MASS OF ",i,
":",
3251 en=sqrt( ppho(5,i)**2 + p2)
3254 &
WRITE(phlun,*)
"CORRECTING ENERGY OF ",i,
":",
3265 IF (iprint.EQ.1)
THEN
3266 WRITE(phlun,*)
"CORRECTING MOTHER"
3267 WRITE(phlun,*)
"PX:",ppho(1,1),
"=>",px-ppho(1,2)
3268 WRITE(phlun,*)
"PY:",ppho(2,1),
"=>",py-ppho(2,2)
3269 WRITE(phlun,*)
"PZ:",ppho(3,1),
"=>",pz-ppho(3,2)
3270 WRITE(phlun,*)
" E:",ppho(4,1),
"=>",e-ppho(4,2)
3273 ppho(1,1)=px-ppho(1,2)
3274 ppho(2,1)=py-ppho(2,2)
3275 ppho(3,1)=pz-ppho(3,2)
3276 ppho(4,1)=e -ppho(4,2)
3278 p2=ppho(1,1)**2+ppho(2,1)**2+ppho(3,1)**2
3280 IF (ppho(4,1)**2.GT.p2)
THEN
3281 m=sqrt( ppho(4,1)**2 - p2 )
3283 &
WRITE(phlun,*)
" M:",ppho(5,1),
"=>",m
3293 FUNCTION phint(IDUM)
3314 parameter(nmxpho=10000)
3315 INTEGER IDPHO,ISTPHO,JDAPHO,JMOPHO,NEVPHO,NPHO
3317 common/phoevt/nevpho,npho,istpho(nmxpho),idpho(nmxpho),
3318 &jmopho(2,nmxpho),jdapho(2,nmxpho),ppho(5,nmxpho),vpho(4,nmxpho)
3320 DOUBLE PRECISION EMU,MCHREN,BETA,COSTHG,MPASQR,XPH, XC1, XC2,XDENO
3321 DOUBLE PRECISION XNUM1,XNUM2
3322 DOUBLE PRECISION EPS1(4),EPS2(4),PH(4),PL(4)
3333 CALL phoeps(ph,eps2,eps1)
3334 CALL phoeps(ph,eps1,eps2)
3341 DO k=jdapho(1,1),npho-1
3349 xc1 = - phocha(idpho(k)) *
3350 & ( pl(1)*eps1(1) + pl(2)*eps1(2) + pl(3)*eps1(3) ) /
3351 & ( ph(4)*pl(4) - ph(1)*pl(1) - ph(2)*pl(2) - ph(3)*pl(3) )
3353 xc2 = - phocha(idpho(k)) *
3354 & ( pl(1)*eps2(1) + pl(2)*eps2(2) + pl(3)*eps2(3) ) /
3355 & ( ph(4)*pl(4) - ph(1)*pl(1) - ph(2)*pl(2) - ph(3)*pl(3) )
3362 xdeno = xdeno + xc1**2 + xc2**2
3366 phint=(xnum1**2 + xnum2**2) / xdeno
3372 SUBROUTINE phoeps (VEC1, VEC2, EPS)
3391 DOUBLE PRECISION VEC1(4), VEC2(4), EPS(4),XN
3393 eps(1)=vec1(2)*vec2(3) - vec1(3)*vec2(2)
3394 eps(2)=vec1(3)*vec2(1) - vec1(1)*vec2(3)
3395 eps(3)=vec1(1)*vec2(2) - vec1(2)*vec2(1)
3398 xn=sqrt( eps(1)**2 +eps(2)**2 +eps(3)**2)
3422 DOUBLE PRECISION SUMVEC(5)
3424 #include
"../include/HEPEVT.h"
3432 WRITE(phlun,9010) nevhep
3438 IF (jdahep(1,i).EQ.0)
THEN
3440 20 sumvec(j)=sumvec(j)+phep(j,i)
3441 IF (jmohep(2,i).EQ.0)
THEN
3442 WRITE(phlun,9030) i,idhep(i),jmohep(1,i),(phep(j,i),j=1,5)
3444 WRITE(phlun,9040) i,idhep(i),jmohep(1,i),jmohep(2,i),(phep
3448 IF (jmohep(2,i).EQ.0)
THEN
3449 WRITE(phlun,9050) i,idhep(i),jmohep(1,i),jdahep(1,i),
3450 & jdahep(2,i),(phep(j,i),j=1,5)
3452 WRITE(phlun,9060) i,idhep(i),jmohep(1,i),jmohep(2,i),
3453 & jdahep(1,i),jdahep(2,i),(phep(j,i),j=1,5)
3457 sumvec(5)=sqrt(sumvec(4)**2-sumvec(1)**2-sumvec(2)**2-
3459 WRITE(phlun,9070) (sumvec(j),j=1,5)
3461 9000
FORMAT(1h0,80(
'='))
3462 9010
FORMAT(1h ,29x,
'Event No.:',i10)
3463 9020
FORMAT(1h0,1x,
'Nr',3x,
'Type',3x,
'Parent(s)',2x,
'Daughter(s)',6x,
3464 &
'Px',7x,
'Py',7x,
'Pz',7x,
'E',4x,
'Inv. M.')
3465 9030
FORMAT(1h ,i4,i7,3x,i4,9x,
'Stable',2x,5f9.2)
3466 9040
FORMAT(1h ,i4,i7,i4,
' - ',i4,5x,
'Stable',2x,5f9.2)
3467 9050
FORMAT(1h ,i4,i7,3x,i4,6x,i4,
' - ',i4,5f9.2)
3468 9060
FORMAT(1h ,i4,i7,i4,
' - ',i4,2x,i4,
' - ',i4,5f9.2)
3469 9070
FORMAT(1h0,23x,
'Vector Sum: ', 5f9.2)
3470 9080
FORMAT(1h0,6x,
'Particle Parameters')