1/* copyright(c) 1991-2025 free software foundation, inc.
2 this file is part of the gnu c library.
4 the gnu c library is free software; you can redistribute it and/or
5 modify it under the terms of the gnu lesser general
Public
6 license as published by the free software foundation; either
7 version 2.1 of the license, or(at your option) any later version.
9 the gnu c library is distributed in the hope that it will be useful,
10 but without any warranty; without even the implied warranty of
11 merchantability or fitness for a particular purpose. see the gnu
12 lesser general
Public license for more details.
14 you should have received a copy of the gnu lesser general
Public
15 license along with the gnu c library;
if not, see
16 <https://www.gnu.org/licenses/>. */
19/* this header is separate from features.h so that the compiler can
20 include it implicitly at the start of every compilation. it must
21 not itself include <features.h> or any other header that includes
22 <features.h> because the
implicit include comes before any feature
23 test macros that may be defined in a source file before it first
24 explicitly includes a system header. gcc knows the name of this
25 header in order to preinclude it. */
27/* glibc
's intent is to support the IEC 559 math functionality, real
28 and complex. If the GCC (4.9 and later) predefined macros
29 specifying compiler intent are available, use them to determine
30 whether the overall intent is to support these features; otherwise,
31 presume an older compiler has intent to support these features and
32 define these macros by default. */
36/* wchar_t uses Unicode 10.0.0. Version 10.0 of the Unicode Standard is
37 synchronized with ISO/IEC 10646:2017, fifth edition, plus
38 the following additions from Amendment 1 to the fifth edition:
41 - 3 additional Zanabazar Square characters */
43*///////////////////////////////////////////////////////////////////////
45*// !!!!!!! WARNING!!!!! This source may be agressive !!!!
47*// Due to short common block names it may owerwrite variables in other
50*// One should add suffix c_Photos_ to names of all commons as soon as
52*///////////////////////////////////////////////////////////////////////
54C.----------------------------------------------------------------------
56C. PHOTOS: PHOtos CDE-s
58C. Purpose: Keep definitions for PHOTOS QED correction Monte Carlo.
60C. Input Parameters: None
62C. Output Parameters: None
64C. Author(s): Z. Was, B. van Eijk Created at: 29/11/89
65C. Last Update: 10/08/93
67C. =========================================================
68C. General Structure Information: =
69C. =========================================================
76C. 2) GENERAL INTERFACE:
83C. PHOTWO (specific interface
86C. PHTYPE (specific interface
87C. PHOMAK (specific interface
88C. 3) QED PHOTON GENERATION:
115C. NAME USED IN SECT. # OF OCC. Comment
116C. PHOQED 1) 2) 3 Flags whether emisson to be gen.
117C. PHOLUN 1) 4) 6 Output device number
118C. PHOCOP 1) 3) 4 photon coupling & min energy
119C. PHPICO 1) 3) 4) 5 PI & 2*PI
120C. PHSEED 1) 4) 3 RN seed
121C. PHOSTA 1) 4) 3 Status information
122C. PHOKEY 1) 2) 3) 7 Keys for nonstandard application
123C. PHOVER 1) 1 Version info for outside
124C. HEPEVT 2) 2 PDG common
125C. PH_HEPEVT2) 8 PDG common internal
126C. PHOEVT 2) 3) 10 PDG branch
127C. PHOIF 2) 3) 2 emission flags for PDG branch
128C. PHOMOM 3) 5 param of char-neutr system
129C. PHOPHS 3) 5 photon momentum parameters
130C. PHOPRO 3) 4 var. for photon rep. (in branch)
131C. PHOCMS 2) 3 parameters of boost to branch CMS
132C. PHNUM 4) 1 event number from outside
133C.----------------------------------------------------------------------
135C.----------------------------------------------------------------------
137C. PHOTOS: PHOton radiation in decays INItialisation
139C. Purpose: Initialisation routine for the PHOTOS QED radiation
140C. package. Should be called at least once before a call
141C. to the steering program 'photos
' is made.
143C. Input Parameters: None
145C. Output Parameters: None
147C. Author(s): Z. Was, B. van Eijk Created at: 26/11/89
148C. Last Update: 12/04/90
150C.----------------------------------------------------------------------
152 INTEGER INIT,IDUM,IPHQRK,IPHEKL
156C-- Return if already initialized...
157.NE.
IF (INIT0) RETURN
160C-- all the following parameter setters can be called after PHOINI.
161C-- Initialization of kinematic correction against rounding errors.
162C-- The set values will be used later if called wit zero.
163C-- Default parameter is 1 (no correction) optionally 2, 3, 4
164C-- In case of exponentiation new version 5 is needed in most cases.
165C-- Definition given here will be thus overwritten in such a case
166C-- below in routine PHOCIN
168C-- blocks emission from quarks if parameter is 1 (enables if 2),
169C-- physical treatment
170C-- will be 3, option 2 is not realistic and for tests only,
171 IDUM= IPHQRK(1) ! default is 1
172C-- blocks emission in pi0 to gamma e+ e- if parameter is gt.1
173C-- (enables otherwise)
174 IDUM= IPHEKL(2) ! default is 1
176C-- Preset parameters in PHOTOS commons
183C-- Initialize Marsaglia and Zaman random number generator
188C.----------------------------------------------------------------------
190C. PHOTOS: PHOton Common INitialisation
192C. Purpose: Initialisation of parameters in common blocks.
194C. Input Parameters: None
196C. Output Parameters: Commons /PHOLUN/, /PHOPHO/, /PHOCOP/, /PHPICO/
199C. Author(s): B. van Eijk Created at: 26/11/89
200C. Z. Was Last Update: 29/01/05
202C.----------------------------------------------------------------------
205 PARAMETER (d_h_NMXHEP=10000)
207 COMMON/PHOQED/QEDRAD(d_h_NMXHEP)
211 COMMON/PHOCOP/ALPHA,XPHCUT
213 COMMON/PHPICO/PI,TWOPI
214 INTEGER ISEED,I97,J97
215 REAL*8 URAN,CRAN,CDRAN,CMRAN
216 COMMON/PHSEED/ISEED(2),I97,J97,URAN(97),CRAN,CDRAN,CMRAN
218 PARAMETER (PHOMES=10)
220 COMMON/PHOSTA/STATUS(PHOMES)
221 LOGICAL INTERF,ISEC,ITRE,IEXP,IFTOP,IFW
222 REAL*8 FINT,FSEC,EXPEPS
223 COMMON /PHOKEY/ FSEC,FINT,EXPEPS,INTERF,ISEC,ITRE,IEXP,IFTOP,IFW
228C-- Return if already initialized...
229.NE.
IF (INIT0) RETURN
232C-- Preset switch for photon emission to 'true
' for each particle in
233C-- /PH_HEPEVT/, this interface is needed for KORALB and KORALZ...
237C-- Logical output unit for printing of PHOTOS error messages
240C-- Set cut parameter for photon radiation
241 XPHCUT=0.01 D0 ! 0.0001D0! to go to low valuex (IEXP excepted)
242C-- ! switch to - VARIANT B
244C-- Define some constants
245 ALPHA=0.00729735039D0
246 PI=3.14159265358979324D0
247 TWOPI=6.28318530717958648D0
249C-- Default seeds Marsaglia and Zaman random number generator
253C-- Iitialization for extra options
255C-- Interference weight now universal.
258C-- Second order - double photon switch
260C-- Third/fourth order - triple (or quatric) photon switch,
261C-- see dipswitch ifour
263C-- Exponentiation on:
268 CALL PHCORK(5) ! in case of exponentiation correction of ph space
269 ! is a default mandatory
274C-- Emision in the hard process g g (q qbar) --> t tbar
278C-- further initialization done automatically
279C-- see places with - VARIANT A - VARIANT B - all over
280C-- to switch between options.
281C ----------- SLOWER VARIANT A, but stable ------------------
282C --- it is limiting choice for small XPHCUT in fixed orer
283C --- modes of operation
285C-- best choice is if FINT=2**N where N+1 is maximal number
286C-- of charged daughters
287C-- see report on overweihted events
292C ----------- FASTER VARIANT B ------------------
293C -- it is good for tests of fixed order and small XPHCUT
294C -- but is less promising for more complex cases of interference
295C -- sometimes fails because of that
302C----------END VARIANTS A B -----------------------
304C-- Effects of initial state charge (in leptonic W decays)
307C-- Initialise status counter for warning messages
313C.----------------------------------------------------------------------
315C. PHOTOS: PHOton radiation in decays general INFo
317C. Purpose: Print PHOTOS info
319C. Input Parameters: PHOLUN
321C. Output Parameters: PHOVN1, PHOVN2
323C. Author(s): B. van Eijk Created at: 12/04/90
324C. Last Update: 27/06/04
326C.----------------------------------------------------------------------
329 INTEGER PHOVN1,PHOVN2
330 COMMON/PHOVER/PHOVN1,PHOVN2
333 LOGICAL INTERF,ISEC,ITRE,IEXP,IFTOP,IFW
334 REAL*8 FINT,FSEC,EXPEPS
335 COMMON /PHOKEY/ FSEC,FINT,EXPEPS,INTERF,ISEC,ITRE,IEXP,IFTOP,IFW
337 COMMON/PHOCOP/ALPHA,XPHCUT
339C-- PHOTOS version number and release date
350 WRITE(PHLUN,9040) IV1,IV2
352 IV2=(PHOVN2-IV1*10000)/100
353 IV3=PHOVN2-IV1*10000-IV2*100
354 WRITE(PHLUN,9050) IV1,IV2,IV3
363 WRITE(PHLUN,9064) INTERF,ISEC,ITRE,IEXP,IFTOP,IFW,ALPHA,XPHCUT
365 IF (INTERF) WRITE(PHLUN,9061)
366 IF (ISEC) WRITE(PHLUN,9062)
367 IF (ITRE) WRITE(PHLUN,9066)
368 IF (IEXP) WRITE(PHLUN,9067) EXPEPS
369 IF (IFTOP) WRITE(PHLUN,9063)
370 IF (IFW) WRITE(PHLUN,9065)
376 9010 FORMAT(1H ,'*
',T81,'*
')
377 9020 FORMAT(1H ,80('*
'))
378 9030 FORMAT(1H ,'*
',26X,26('=
'),T81,'*
')
379 9040 FORMAT(1H ,'*
',28X,'photos, version:
',I2,'.
',I2,T81,'*
')
380 9050 FORMAT(1H ,'*
',28X,'released at:
',I2,'/
',I2,'/
',I2,T81,'*
')
381 9060 FORMAT(1H ,'*
',18X,'photos qed corrections in particle decays
',
383 9061 FORMAT(1H ,'*
',18X,'option with interference is active
',
385 9062 FORMAT(1H ,'*
',18X,'option with double photons is active
',
387 9066 FORMAT(1H ,'*
',18X,'option with triple/quatric photons is active
',
389 9067 FORMAT(1H ,'*
',18X,'option with exponentiation is active epsexp=
',
391 9063 FORMAT(1H ,'*
',18X,'emision in t tbar production is active
',
393 9064 FORMAT(1H ,'*
',18X,'internal input parameters:
',T81,'*
'
395 &,/, 1H ,'*
',18X,'interf=
',L2,' isec=
',L2,' itre=
',L2,
396 & ' iexp=
',L2,' iftop=
',L2,
398 &,/, 1H ,'*
',18X,'alpha_qed=
',F8.5,' xphcut=
',E8.3,T81,'*
')
399 9065 FORMAT(1H ,'*
',18X,'correction wt in decay of w is active
',
401 9070 FORMAT(1H ,'*
',9X,
402 &'monte carlo
Program - by e. barberio, b. van eijk and z. was
',
404 & 1H ,'*
',9X,'version 2.09 - by p. golonka and z.w.
',T81,'*
')
405 9080 FORMAT( 1H ,'*
',9X,' ',T81,'*
',/,
407 & ' warning(1): /hepevt/ is not anymore the standard
common block
'
409 & 1H ,'*
',9X,' ',T81,'*
',/,
411 & ' photos expects /hepevt/ to have real*8 variables. to change to
'
412 & ,T81,'*
',/, 1H ,'*
',9X,
413 & ' real*4 modify its declaration in subr. photos_get photos_set:
'
414 & ,T81,'*
',/, 1H ,'*
',9X,
415 & ' real*8 d_h_phep, d_h_vhep
'
416 & ,T81,'*
',/, 1H ,'*
',9X,
417 & ' warning(2): check dims. of /hepevt/ /phoqed/ /ph_hepevt/.
'
418 & ,T81,'*
',/, 1H ,'*
',9X,
419 & ' here: d_h_nmxhep=10000 and nmxhep=10000
'
422 SUBROUTINE PHOTOS(ID)
423C.----------------------------------------------------------------------
425C. PHOTOS: General search routine + _GET + _SET
427C. Purpose: /HEPEVT/ is not anymore a standard at least
428C. REAL*8 REAL*4 are in use. PHOTOS_GET and PHOTOS_SET
429C. were to be introduced.
432C. Input Parameters: ID see routine PHOTOS_MAKE
434C. Output Parameters: None
436C. Author(s): Z. Was Created at: 21/07/98
437C. Last Update: 21/07/98
439C.----------------------------------------------------------------------
440 COMMON /PHLUPY/ IPOIN,IPOINM
447.GT..AND..LT.
IF (1IPOINM1IPOIN ) THEN
448 WRITE(PHLUN,*) 'event nr=
',IEV,
449 $ 'we are testing /hepevt/ at ipoint=1 (input)
'
455.GT..AND..LT.
IF (1IPOINM1IPOIN ) THEN
456 WRITE(PHLUN,*) 'event nr=
',IEV,
457 $ 'we are testing /hepevt/ at ipoint=1 (output)
'
463 SUBROUTINE PHOTOS_GET
464C.----------------------------------------------------------------------
468C. Purpose: Copies /HEPEVT/ into /PH_HEPEVT/
471C. Input Parameters: None
473C. Output Parameters: None
475C. Author(s): Z. Was Created at: 21/07/98
476C. Last Update: 21/07/98
478C.----------------------------------------------------------------------
481 INTEGER d_h_nmxhep ! maximum number of particles
482 PARAMETER (d_h_NMXHEP=10000)
483 REAL*8 d_h_phep, d_h_vhep ! to be real*4/ *8 depending on host
484 INTEGER d_h_nevhep,d_h_nhep,d_h_isthep,d_h_idhep,d_h_jmohep,
487 $ d_h_nevhep, ! serial number
488 $ d_h_nhep, ! number of particles
489 $ d_h_isthep(d_h_nmxhep), ! status code
490 $ d_h_idhep(d_h_nmxhep), ! particle ident KF
491 $ d_h_jmohep(2,d_h_nmxhep), ! parent particles
492 $ d_h_jdahep(2,d_h_nmxhep), ! childreen particles
493 $ d_h_phep(5,d_h_nmxhep), ! four-momentum, mass [GeV]
494 $ d_h_vhep(4,d_h_nmxhep) ! vertex [mm]
495* ----------------------------------------------------------------------
498 $ d_h_qedrad(d_h_nmxhep) ! Photos flag
500 PARAMETER (NMXHEP=10000)
501 INTEGER IDHEP,ISTHEP,JDAHEP,JMOHEP,NEVHEP,NHEP
503 COMMON/PH_HEPEVT/NEVHEP,NHEP,ISTHEP(NMXHEP),IDHEP(NMXHEP),
504 &JMOHEP(2,NMXHEP),JDAHEP(2,NMXHEP),PHEP(5,NMXHEP),VHEP(4,NMXHEP)
506 COMMON/PH_PHOQED/QEDRAD(NMXHEP)
508 nevhep= d_h_nevhep ! serial number
509 nhep = d_h_nhep ! number of particles
511 isthep(k) =d_h_isthep(k) ! status code
512 idhep(k) =d_h_idhep(k) ! particle ident KF
513 jmohep(1,k) =d_h_jmohep(1,k) ! parent particles
514 jdahep(1,k) =d_h_jdahep(1,k) ! childreen particles
515 jmohep(2,k) =d_h_jmohep(2,k) ! parent particles
516 jdahep(2,k) =d_h_jdahep(2,k) ! childreen particles
518 phep(l,k) =d_h_phep(l,k) ! four-momentum, mass [GeV]
519 vhep(l,k) =d_h_vhep(l,k) ! vertex [mm]
521 phep(5,k) =d_h_phep(5,k) ! four-momentum, mass [GeV]
522 qedrad(k) =d_h_qedrad(k) ! Photos special flag
527 SUBROUTINE PHOTOS_SET
528C.----------------------------------------------------------------------
532C. Purpose: Copies /PH_HEPEVT/ into /HEPEVT/
535C. Input Parameters: None
537C. Output Parameters: None
539C. Author(s): Z. Was Created at: 21/07/98
540C. Last Update: 21/07/98
542C.----------------------------------------------------------------------
544 INTEGER d_h_nmxhep ! maximum number of particles
545 PARAMETER (d_h_NMXHEP=10000)
546 REAL*8 d_h_phep, d_h_vhep ! to be real*4/ *8 depending on host
547 INTEGER d_h_nevhep,d_h_nhep,d_h_isthep,d_h_idhep,d_h_jmohep,
550 $ d_h_nevhep, ! serial number
551 $ d_h_nhep, ! number of particles
552 $ d_h_isthep(d_h_nmxhep), ! status code
553 $ d_h_idhep(d_h_nmxhep), ! particle ident KF
554 $ d_h_jmohep(2,d_h_nmxhep), ! parent particles
555 $ d_h_jdahep(2,d_h_nmxhep), ! childreen particles
556 $ d_h_phep(5,d_h_nmxhep), ! four-momentum, mass [GeV]
557 $ d_h_vhep(4,d_h_nmxhep) ! vertex [mm]
558* ----------------------------------------------------------------------
561 $ d_h_qedrad(d_h_nmxhep) ! Photos flag
563 PARAMETER (NMXHEP=10000)
564 INTEGER IDHEP,ISTHEP,JDAHEP,JMOHEP,NEVHEP,NHEP
566 COMMON/PH_HEPEVT/NEVHEP,NHEP,ISTHEP(NMXHEP),IDHEP(NMXHEP),
567 &JMOHEP(2,NMXHEP),JDAHEP(2,NMXHEP),PHEP(5,NMXHEP),VHEP(4,NMXHEP)
569 COMMON/PH_PHOQED/QEDRAD(NMXHEP)
572 d_h_nevhep= nevhep ! serial number
573 d_h_nhep = nhep ! number of particles
575 d_h_isthep(k) =isthep(k) ! status code
576 d_h_idhep(k) =idhep(k) ! particle ident KF
577 d_h_jmohep(1,k) =jmohep(1,k) ! parent particles
578 d_h_jdahep(1,k) =jdahep(1,k) ! childreen particles
579 d_h_jmohep(2,k) =jmohep(2,k) ! parent particles
580 d_h_jdahep(2,k) =jdahep(2,k) ! childreen particles
582 d_h_phep(l,k) =phep(l,k) ! four-momentum, mass [GeV]
583 d_h_vhep(l,k) =vhep(l,k) ! vertex [mm]
585 d_h_phep(5,k) =phep(5,k) ! four-momentum, mass [GeV]
586 d_h_qedrad(k) =qedrad(k) ! Photos special flag
589 SUBROUTINE PHOTOS_MAKE(IPARR)
590C.----------------------------------------------------------------------
592C. PHOTOS_MAKE: General search routine
594C. Purpose: Search through the /PH_HEPEVT/ standard HEP common, sta-
595C. rting from the IPPAR-th particle. Whenevr branching
596C. point is found routine PHTYPE(IP) is called.
597C. Finally if calls on PHTYPE(IP) modified entries, common
598C /PH_HEPEVT/ is ordered.
600C. Input Parameter: IPPAR: Pointer to decaying particle in
601C. /PH_HEPEVT/ and the common itself,
603C. Output Parameters: Common /PH_HEPEVT/, either with or without
604C. new particles added.
606C. Author(s): Z. Was, B. van Eijk Created at: 26/11/89
607C. Last Update: 30/08/93
609C.----------------------------------------------------------------------
612 INTEGER IP,IPARR,IPPAR,I,J,K,L,NLAST
613 DOUBLE PRECISION DATA
614 INTEGER MOTHER,POSPHO
617 PARAMETER (NMXHEP=10000)
618 INTEGER IDHEP,ISTHEP,JDAHEP,JMOHEP,NEVHEP,NHEP
620 COMMON/PH_HEPEVT/NEVHEP,NHEP,ISTHEP(NMXHEP),IDHEP(NMXHEP),
621 &JMOHEP(2,NMXHEP),JDAHEP(2,NMXHEP),PHEP(5,NMXHEP),VHEP(4,NMXHEP)
623 COMMON/PH_PHOQED/QEDRAD(NMXHEP)
625 PARAMETER (NMXPHO=10000)
626 INTEGER ISTACK(0:NMXPHO),NUMIT,NTRY,KK,LL,II,NA,FIRST,LAST
627 INTEGER FIRSTA,LASTA,IPP,IDA1,IDA2,MOTHER2,IDPHO,ISPHO
628 REAL*8 PORIG(5,NMXPHO)
631C-- Store pointers for cascade treatement...
636C-- Check decay multiplicity and minimum of correctness..
637.EQ..OR..NE.
IF ((JDAHEP(1,IP)0)(JMOHEP(1,JDAHEP(1,IP))IP)) RETURN
639C-- single branch mode
640C-- we start looking for the decay points in the cascade
641C-- IPPAR is original position where the program was called
643C-- NUMIT denotes number of secondary decay branches
645C-- NTRY denotes number of secondary branches already checked for
646C-- for existence of further branches
648C-- let-s search if IPARR does not prevent searching.
651 DO I=JDAHEP(1,IP),JDAHEP(2,IP)
652.NE..AND..EQ.
IF (JDAHEP(1,I)0JMOHEP(1,JDAHEP(1,I))I) THEN
654.GT.
IF (NUMITNMXPHO) THEN
656 CALL PHOERR(7,'photos
',DATA)
661.GT.
IF(NUMITNTRY) THEN
667C-- let-s do generation
670 FIRST=JDAHEP(1,ISTACK(KK))
671 LAST=JDAHEP(2,ISTACK(KK))
674 PORIG(LL,II)=PHEP(LL,FIRST+II-1)
678 CALL PHTYPE(ISTACK(KK))
680C-- Correct energy/momentum of cascade daughters
686.EQ.
IF(JMOHEP(1,IPP)ISTACK(KK))
687 $ CALL PHOBOS(IPP,PORIG(1,II),PHEP(1,IPP),FIRSTA,LASTA)
692C-- rearrange /PH_HEPEVT/ to get correct order..
693.GT.
IF (NHEPNLAST) THEN
694 DO 160 I=NLAST+1,NHEP
696C-- Photon mother and position...
698 POSPHO=JDAHEP(2,MOTHER)+1
699C-- Intermediate save of photon energy/momentum and pointers
701 90 PHOTON(J)=PHEP(J,I)
708C-- Exclude photon in sequence !
709.NE.
IF (POSPHONHEP) THEN
713 DO 120 K=I,POSPHO+1,-1
714 ISTHEP(K)=ISTHEP(K-1)
715 QEDRAD(K)=QEDRAD(K-1)
718 JMOHEP(L,K)=JMOHEP(L,K-1)
719 100 JDAHEP(L,K)=JDAHEP(L,K-1)
721 110 PHEP(L,K)=PHEP(L,K-1)
723 120 VHEP(L,K)=VHEP(L,K-1)
725C-- Correct pointers assuming most dirty /PH_HEPEVT/...
728.NE..AND..GE.
IF ((JMOHEP(L,K)0)(JMOHEP(L,K)
729 & POSPHO)) JMOHEP(L,K)=JMOHEP(L,K)+1
730.NE..AND..GE.
IF ((JDAHEP(L,K)0)(JDAHEP(L,K)
731 & POSPHO)) JDAHEP(L,K)=JDAHEP(L,K)+1
734C-- Store photon energy/momentum
736 140 PHEP(J,POSPHO)=PHOTON(J)
739C-- Store pointers for the photon...
740 JDAHEP(2,MOTHER)=POSPHO
743 JMOHEP(1,POSPHO)=MOTHER
744 JMOHEP(2,POSPHO)=MOTHER2
745 JDAHEP(1,POSPHO)=IDA1
746 JDAHEP(2,POSPHO)=IDA2
748C-- Get photon production vertex position
750 150 VHEP(J,POSPHO)=VHEP(J,POSPHO-1)
755 SUBROUTINE PHOBOS(IP,PBOOS1,PBOOS2,FIRST,LAST)
756C.----------------------------------------------------------------------
758C. PHOBOS: PHOton radiation in decays BOoSt routine
760C. Purpose: Boost particles in cascade decay to parent rest frame
761C. and boost back with modified boost vector.
763C. Input Parameters: IP: pointer of particle starting chain
765C. PBOOS1: Boost vector to rest frame,
766C. PBOOS2: Boost vector to modified frame,
767C. FIRST: Pointer to first particle to be boos-
769C. LAST: Pointer to last particle to be boos-
772C. Output Parameters: Common /PH_HEPEVT/.
774C. Author(s): B. van Eijk Created at: 13/02/90
775C. Z. Was Last Update: 16/11/93
777C.----------------------------------------------------------------------
779 DOUBLE PRECISION BET1(3),BET2(3),GAM1,GAM2,PB,DATA
780 INTEGER I,J,FIRST,LAST,MAXSTA,NSTACK,IP
781 PARAMETER (MAXSTA=10000)
782 INTEGER STACK(MAXSTA)
783 REAL*8 PBOOS1(5),PBOOS2(5)
785 PARAMETER (NMXHEP=10000)
786 INTEGER IDHEP,ISTHEP,JDAHEP,JMOHEP,NEVHEP,NHEP
788 COMMON/PH_HEPEVT/NEVHEP,NHEP,ISTHEP(NMXHEP),IDHEP(NMXHEP),
789 &JMOHEP(2,NMXHEP),JDAHEP(2,NMXHEP),PHEP(5,NMXHEP),VHEP(4,NMXHEP)
790.EQ..OR..LT.
IF ((LAST0)(LASTFIRST)) RETURN
793 BET1(J)=-PBOOS1(J)/PBOOS1(5)
794 10 BET2(J)=PBOOS2(J)/PBOOS2(5)
795 GAM1=PBOOS1(4)/PBOOS1(5)
796 GAM2=PBOOS2(4)/PBOOS2(5)
798C-- Boost vector to parent rest frame...
799 20 DO 50 I=FIRST,LAST
800 PB=BET1(1)*PHEP(1,I)+BET1(2)*PHEP(2,I)+BET1(3)*PHEP(3,I)
801.EQ.
IF (JMOHEP(1,I)IP) THEN
803 30 PHEP(J,I)=PHEP(J,I)+BET1(J)*(PHEP(4,I)+PB/(GAM1+1.D0))
804 PHEP(4,I)=GAM1*PHEP(4,I)+PB
806C-- ...and boost back to modified parent frame.
807 PB=BET2(1)*PHEP(1,I)+BET2(2)*PHEP(2,I)+BET2(3)*PHEP(3,I)
809 40 PHEP(J,I)=PHEP(J,I)+BET2(J)*(PHEP(4,I)+PB/(GAM2+1.D0))
810 PHEP(4,I)=GAM2*PHEP(4,I)+PB
811.NE.
IF (JDAHEP(1,I)0) THEN
814C-- Check on stack length...
815.GT.
IF (NSTACKMAXSTA) THEN
817 CALL PHOERR(7,'phobos
',DATA)
823.NE.
IF (NSTACK0) THEN
825C-- Now go one step further in the decay tree...
826 FIRST=JDAHEP(1,STACK(NSTACK))
827 LAST=JDAHEP(2,STACK(NSTACK))
834 SUBROUTINE PHOIN(IP,BOOST,NHEP0)
835C.----------------------------------------------------------------------
837C. PHOIN: PHOtos INput
839C. Purpose: copies IP branch of the common /PH_HEPEVT/ into /PHOEVT/
840C. moves branch into its CMS system.
842C. Input Parameters: IP: pointer of particle starting branch
844C. BOOST: Flag whether boost to CMS was or was
847C. Output Parameters: Commons: /PHOEVT/, /PHOCMS/
849C. Author(s): Z. Was Created at: 24/05/93
850C. Last Update: 16/11/93
852C.----------------------------------------------------------------------
855 PARAMETER (NMXHEP=10000)
856 INTEGER IDHEP,ISTHEP,JDAHEP,JMOHEP,NEVHEP,NHEP
858 COMMON/PH_HEPEVT/NEVHEP,NHEP,ISTHEP(NMXHEP),IDHEP(NMXHEP),
859 &JMOHEP(2,NMXHEP),JDAHEP(2,NMXHEP),PHEP(5,NMXHEP),VHEP(4,NMXHEP)
861 PARAMETER (NMXPHO=10000)
862 INTEGER IDPHO,ISTPHO,JDAPHO,JMOPHO,NEVPHO,NPHO
864 COMMON/PHOEVT/NEVPHO,NPHO,ISTPHO(NMXPHO),IDPHO(NMXPHO),
865 &JMOPHO(2,NMXPHO),JDAPHO(2,NMXPHO),PPHO(5,NMXPHO),VPHO(4,NMXPHO)
866 INTEGER IP,IP2,I,FIRST,LAST,LL,NA
869 DOUBLE PRECISION BET(3),GAM,PB
870 COMMON /PHOCMS/ BET,GAM
871 LOGICAL INTERF,ISEC,ITRE,IEXP,IFTOP,IFW
872 REAL*8 FINT,FSEC,EXPEPS
873 COMMON /PHOKEY/ FSEC,FINT,EXPEPS,INTERF,ISEC,ITRE,IEXP,IFTOP,IFW
875C let-s calculate size of the little common entry
878 NPHO=3+LAST-FIRST+NHEP-NHEP0
880C let-s take in decaying particle
883 JDAPHO(2,1)=3+LAST-FIRST
887C let-s take in eventual second mother
888 IP2=JMOHEP(2,JDAHEP(1,IP))
889.NE..AND..NE.
IF((IP20)(IP2IP)) THEN
892 JDAPHO(2,2)=3+LAST-FIRST
894 PPHO(I,2)=PHEP(I,IP2)
902C let-s take in daughters
904 IDPHO(3+LL)=IDHEP(FIRST+LL)
905 JMOPHO(1,3+LL)=JMOHEP(1,FIRST+LL)
906.EQ.
IF (JMOHEP(1,FIRST+LL)IP) JMOPHO(1,3+LL)=1
908 PPHO(I,3+LL)=PHEP(I,FIRST+LL)
911.GT.
IF (NHEPNHEP0) THEN
912C let-s take in illegitimate daughters
915 IDPHO(NA+LL)=IDHEP(NHEP0+LL)
916 JMOPHO(1,NA+LL)=JMOHEP(1,NHEP0+LL)
917.EQ.
IF (JMOHEP(1,NHEP0+LL)IP) JMOPHO(1,NA+LL)=1
919 PPHO(I,NA+LL)=PHEP(I,NHEP0+LL)
922C-- there is NHEP-NHEP0 daugters more.
923 JDAPHO(2,1)=3+LAST-FIRST+NHEP-NHEP0
925.EQ.
IF(IDPHO(NPHO)22)CALL PHLUPA(100001)
926.EQ.
! IF(IDPHO(NPHO)22) stop
928.EQ.
IF(IDPHO(NPHO)22)CALL PHLUPA(100002)
929C special case of t tbar production process
930 IF(IFTOP) CALL PHOTWO(0)
932C-- Check whether parent is in its rest frame...
934C bug reported by Vladimir Savinov localized and fixed.
935C protection against rounding error was back-firing if soft
936C momentum of mother was physical. Consequence was that PHCORK was
937C messing up masses of final state particles in vertex of the decay.
938C Only configurations with previously generated photons of energy fraction
939C smaller than 0.0001 were affected. Effect was numerically insignificant.
941.GT.
C IF ( (ABS(PPHO(4,1)-PPHO(5,1))PPHO(5,1)*1.D-8)
942.AND..NE.
C $ (PPHO(5,1)0)) THEN
944.GT.
IF ((ABS(PPHO(1,1)+ABS(PPHO(2,1))+ABS(PPHO(3,1)))
945.AND..NE.
$ PPHO(5,1)*1.D-8)(PPHO(5,1)0)) THEN
949C-- Boost daughter particles to rest frame of parent...
950C-- Resultant neutral system already calculated in rest frame !
952 10 BET(J)=-PPHO(J,1)/PPHO(5,1)
953 GAM=PPHO(4,1)/PPHO(5,1)
954 DO 30 I=JDAPHO(1,1),JDAPHO(2,1)
955 PB=BET(1)*PPHO(1,I)+BET(2)*PPHO(2,I)+BET(3)*PPHO(3,I)
957 20 PPHO(J,I)=PPHO(J,I)+BET(J)*(PPHO(4,I)+PB/(GAM+1.D0))
958 30 PPHO(4,I)=GAM*PPHO(4,I)+PB
959C-- Finally boost mother as well
961 PB=BET(1)*PPHO(1,I)+BET(2)*PPHO(2,I)+BET(3)*PPHO(3,I)
963 PPHO(J,I)=PPHO(J,I)+BET(J)*(PPHO(4,I)+PB/(GAM+1.D0))
965 PPHO(4,I)=GAM*PPHO(4,I)+PB
967C special case of t tbar production process
968 IF(IFTOP) CALL PHOTWO(1)
970.EQ.
IF(IDPHO(NPHO)22) CALL PHLUPA(10000)
971.EQ.
! IF(IDPHO(NPHO-1)22) stop
973 SUBROUTINE PHOTWO(MODE)
974C.----------------------------------------------------------------------
976C. PHOTWO: PHOtos but TWO mothers allowed
978C. Purpose: Combines two mothers into one in /PHOEVT/
979C. necessary eg in case of g g (q qbar) --> t tbar
981C. Input Parameters: Common /PHOEVT/ (/PHOCMS/)
983C. Output Parameters: Common /PHOEVT/, (stored mothers)
985C. Author(s): Z. Was Created at: 5/08/93
986C. Last Update:10/08/93
988C.----------------------------------------------------------------------
991 PARAMETER (NMXPHO=10000)
992 INTEGER IDPHO,ISTPHO,JDAPHO,JMOPHO,NEVPHO,NPHO
994 COMMON/PHOEVT/NEVPHO,NPHO,ISTPHO(NMXPHO),IDPHO(NMXPHO),
995 &JMOPHO(2,NMXPHO),JDAPHO(2,NMXPHO),PPHO(5,NMXPHO),VPHO(4,NMXPHO)
996 DOUBLE PRECISION BET(3),GAM
997 COMMON /PHOCMS/ BET,GAM
1001C logical IFRAD is used to tag cases when two mothers may be
1002C merged to the sole one.
1003C So far used in case:
1004C 1) of t tbar production
1008.EQ..AND..EQ.
IFRAD=(IDPHO(1)21)(IDPHO(2)21)
1009.OR..EQ..AND..LE.
IFRAD=IFRAD(IDPHO(1)-IDPHO(2)ABS(IDPHO(1))6)
1011.AND..EQ..AND..EQ.
& (ABS(IDPHO(3))6)(ABS(IDPHO(4))6)
1012 MPASQR= (PPHO(4,1)+PPHO(4,2))**2-(PPHO(3,1)+PPHO(3,2))**2
1013 & -(PPHO(2,1)+PPHO(2,2))**2-(PPHO(1,1)+PPHO(1,2))**2
1014.AND..GT.
IFRAD=IFRAD(MPASQR0.0D0)
1016c.....combining first and second mother
1018 PPHO(I,1)=PPHO(I,1)+PPHO(I,2)
1020 PPHO(5,1)=SQRT(MPASQR)
1021c.....removing second mother,
1027C boosting of the mothers to the reaction frame not implemented yet.
1028C to do it in mode 0 original mothers have to be stored in new comon (?)
1029C and in mode 1 boosted to cms.
1032 SUBROUTINE PHOOUT(IP,BOOST,NHEP0)
1033C.----------------------------------------------------------------------
1035C. PHOOUT: PHOtos OUTput
1037C. Purpose: copies back IP branch of the common /PH_HEPEVT/ from
1038C. /PHOEVT/ moves branch back from its CMS system.
1040C. Input Parameters: IP: pointer of particle starting branch
1042C. BOOST: Flag whether boost to CMS was or was
1045C. Output Parameters: Common /PHOEVT/,
1047C. Author(s): Z. Was Created at: 24/05/93
1050C.----------------------------------------------------------------------
1053 PARAMETER (NMXHEP=10000)
1054 INTEGER IDHEP,ISTHEP,JDAHEP,JMOHEP,NEVHEP,NHEP
1056 COMMON/PH_HEPEVT/NEVHEP,NHEP,ISTHEP(NMXHEP),IDHEP(NMXHEP),
1057 &JMOHEP(2,NMXHEP),JDAHEP(2,NMXHEP),PHEP(5,NMXHEP),VHEP(4,NMXHEP)
1059 PARAMETER (NMXPHO=10000)
1060 INTEGER IDPHO,ISTPHO,JDAPHO,JMOPHO,NEVPHO,NPHO
1062 COMMON/PHOEVT/NEVPHO,NPHO,ISTPHO(NMXPHO),IDPHO(NMXPHO),
1063 &JMOPHO(2,NMXPHO),JDAPHO(2,NMXPHO),PPHO(5,NMXPHO),VPHO(4,NMXPHO)
1064 INTEGER IP,LL,FIRST,LAST,I
1066 INTEGER NN,J,K,NHEP0,NA
1067 DOUBLE PRECISION BET(3),GAM,PB
1068 COMMON /PHOCMS/ BET,GAM
1069.EQ.
IF(NPHONEVPHO) RETURN
1070C-- When parent was not in its rest-frame, boost back...
1073 DO 110 J=JDAPHO(1,1),JDAPHO(2,1)
1074 PB=-BET(1)*PPHO(1,J)-BET(2)*PPHO(2,J)-BET(3)*PPHO(3,J)
1076 100 PPHO(K,J)=PPHO(K,J)-BET(K)*(PPHO(4,J)+PB/(GAM+1.D0))
1077 110 PPHO(4,J)=GAM*PPHO(4,J)+PB
1078C-- ...boost photon, or whatever else has shown up
1080 PB=-BET(1)*PPHO(1,NN)-BET(2)*PPHO(2,NN)-BET(3)*PPHO(3,NN)
1082 120 PPHO(K,NN)=PPHO(K,NN)-BET(K)*(PPHO(4,NN)+PB/(GAM+1.D0))
1083 PPHO(4,NN)=GAM*PPHO(4,NN)+PB
1088C let-s take in original daughters
1090 IDHEP(FIRST+LL) = IDPHO(3+LL)
1092 PHEP(I,FIRST+LL) = PPHO(I,3+LL)
1095C let-s take newcomers to the end of HEPEVT.
1098 IDHEP(NHEP0+LL) = IDPHO(NA+LL)
1099 ISTHEP(NHEP0+LL)=ISTPHO(NA+LL)
1100 JMOHEP(1,NHEP0+LL)=IP
1101 JMOHEP(2,NHEP0+LL)=JMOHEP(2,JDAHEP(1,IP))
1102 JDAHEP(1,NHEP0+LL)=0
1103 JDAHEP(2,NHEP0+LL)=0
1105 PHEP(I,NHEP0+LL) = PPHO(I,NA+LL)
1108 NHEP=NHEP+NPHO-NEVPHO
1111 SUBROUTINE PHOCHK(JFIRST)
1112C.----------------------------------------------------------------------
1114C. PHOCHK: checking branch.
1116C. Purpose: checks whether particles in the common block /PHOEVT/
1117C. can be served by PHOMAK.
1118C. JFIRST is the position in /PH_HEPEVT/ (!) of the first
1119C. daughter of sub-branch under action.
1122C. Author(s): Z. Was Created at: 22/10/92
1123C. Last Update: 11/12/00
1125C.----------------------------------------------------------------------
1126C ********************
1129 PARAMETER (NMXPHO=10000)
1130 INTEGER IDPHO,ISTPHO,JDAPHO,JMOPHO,NEVPHO,NPHO
1132 COMMON/PHOEVT/NEVPHO,NPHO,ISTPHO(NMXPHO),IDPHO(NMXPHO),
1133 &JMOPHO(2,NMXPHO),JDAPHO(2,NMXPHO),PPHO(5,NMXPHO),VPHO(4,NMXPHO)
1135 COMMON/PHOIF/CHKIF(NMXPHO)
1137 PARAMETER (NMXHEP=10000)
1139 COMMON/PH_PHOQED/QEDRAD(NMXHEP)
1142 INTEGER IDABS,NLAST,I,IPPAR
1143 LOGICAL INTERF,ISEC,ITRE,IEXP,IFTOP,IFW,IFNPI0,IFKL
1144 REAL*8 FINT,FSEC,EXPEPS
1145 COMMON /PHOKEY/ FSEC,FINT,EXPEPS,INTERF,ISEC,ITRE,IEXP,IFTOP,IFW
1147 INTEGER IDENT,K,IQRK,IPHQRK,IEKL,IPHEKL
1148C these are OK .... if you do not like somebody else, add here.
1150.GT..OR..NE..AND..LE.
& ( ((IDABS9IQRK1)(IDABS40))
1151.OR..GT.
& (IDABS100) )
1152.AND..NE.
& (IDABS21)
1153.AND..NE..AND..NE..AND..NE.
$ (IDABS2101)(IDABS3101)(IDABS3201)
1154.AND..NE..AND..NE..AND..NE.
& (IDABS1103)(IDABS2103)(IDABS2203)
1155.AND..NE..AND..NE..AND..NE.
& (IDABS3103)(IDABS3203)(IDABS3303)
1157 IQRK=IPHQRK(0) ! switch for emission from quark
1162C checking for good particles
1164.GT.
IF (IEKL1) THEN ! exclude radiative corr in decay of pi0
1165C ! and Kl --> ee gamma
1166.NE.
IFNPI0= (IDPHO(1)111) ! pi0
1167.EQ..AND.
IFKL = ((IDPHO(1)130) ! Kl --> ee gamma
1168.EQ..OR..EQ..OR.
$ ((IDPHO(3)22)(IDPHO(4)22)
1169.EQ..AND.
$ (IDPHO(5)22))
1170.EQ..OR..EQ..OR.
$ ((IDPHO(3)11)(IDPHO(4)11)
1171.EQ.
$ (IDPHO(5)11)) )
1173.AND..NOT.
IFNPI0=(IFNPI0(IFKL))
1176 IDABS = ABS(IDPHO(I))
1177C possibly call on PHZODE is a dead (to be omitted) code.
1178.AND.
CHKIF(I)= F(IDABS) F(ABS(IDPHO(1)))
1179.AND..EQ.
& (IDPHO(2)0)
1180.GT..AND.
IF(I2) CHKIF(I)=CHKIF(I)QEDRAD(JFIRST+I-IPPAR-2)
1184C now we go to special cases, where CHKIF(I) will be overwritten
1187C special case of top pair production
1188 DO K=JDAPHO(2,1),JDAPHO(1,1),-1
1189.NE.
IF(IDPHO(K)22) THEN
1195.EQ..AND..EQ.
IFRAD=((IDPHO(1)21)(IDPHO(2)21))
1196.OR..LE..AND..EQ.
& ((ABS(IDPHO(1))6)((IDPHO(2))(-IDPHO(1))))
1198.AND..EQ..AND..EQ.
& (ABS(IDPHO(3))6)((IDPHO(4))(-IDPHO(3)))
1203.GT..AND.
IF(I2) CHKIF(I)=CHKIF(I)QEDRAD(JFIRST+I-IPPAR-2)
1210C special case of top decay
1211 DO K=JDAPHO(2,1),JDAPHO(1,1),-1
1212.NE.
IF(IDPHO(K)22) THEN
1218.EQ..AND..EQ.
IFRAD=((ABS(IDPHO(1))6)(IDPHO(2)0))
1220.AND..EQ..AND..EQ.
& ((ABS(IDPHO(3))24)(ABS(IDPHO(4))5)
1221.OR..EQ..AND..EQ.
& (ABS(IDPHO(3))5)(ABS(IDPHO(4))24))
1226.GT..AND.
IF(I2) CHKIF(I)=CHKIF(I)QEDRAD(JFIRST+I-IPPAR-2)
1233 SUBROUTINE PHTYPE(ID)
1234C.----------------------------------------------------------------------
1236C. PHTYPE: Central manadgement routine.
1238C. Purpose: defines what kind of the
1239C. actions will be performed at point ID.
1241C. Input Parameters: ID: pointer of particle starting branch
1242C. in /PH_HEPEVT/ to be treated.
1244C. Output Parameters: Common /PH_HEPEVT/.
1246C. Author(s): Z. Was Created at: 24/05/93
1247C. P. Golonka Last Update: 27/06/04
1249C.----------------------------------------------------------------------
1252 PARAMETER (NMXHEP=10000)
1253 INTEGER IDHEP,ISTHEP,JDAHEP,JMOHEP,NEVHEP,NHEP
1255 COMMON/PH_HEPEVT/NEVHEP,NHEP,ISTHEP(NMXHEP),IDHEP(NMXHEP),
1256 &JMOHEP(2,NMXHEP),JDAHEP(2,NMXHEP),PHEP(5,NMXHEP),VHEP(4,NMXHEP)
1257 LOGICAL INTERF,ISEC,ITRE,IEXP,IFTOP,IFW
1258 REAL*8 FINT,FSEC,EXPEPS
1259 COMMON /PHOKEY/ FSEC,FINT,EXPEPS,INTERF,ISEC,ITRE,IEXP,IFTOP,IFW
1263 REAL*8 PRO,PRSUM,ESU
1264 COMMON /PHOEXP/ PRO(NX),NCHAN,EXPINI
1268 REAL*8 RN,PHORAN,SUM
1272.AND.
IFOUR=(.TRUE.)(ITRE) ! we can make internal choice whether
1273 ! we want 3 or four photons at most.
1275C-- Check decay multiplicity..
1276.EQ.
IF (JDAHEP(1,ID)0) RETURN
1277.EQ.
C IF (JDAHEP(1,ID)JDAHEP(2,ID)) RETURN
1282 EXPINI=.TRUE. ! Initialization/cleaning
1289 CALL PHOMAK(ID,NHEP0)! Initialization/crude formfactors into
1297 ESU=EXP(-PRSUM) ! exponent for crude Poissonian multiplicity
1298 ! distribution, will be later overwritten
1299 ! to give probability for k
1300 SUM=ESU ! distribuant for the crude Poissonian
1302 DO K=1,100 ! hard coded max (photon) multiplicity is 100
1303.LT.
IF(RNSUM) GOTO 100
1304 ESU=ESU*PRSUM/K ! we get at K ESU=EXP(-PRSUM)*PRSUM**K/K!
1305 SUM=SUM+ESU ! thus we get distribuant at K.
1307 CALL PHOMAK(ID,NHEP0) ! LOOPING
1308.GT.
IF(SUM1D0-EXPEPS) GOTO 100
1312C-- quatro photon emission
1315.GE.
IF (RN23.D0/24D0) THEN
1316 CALL PHOMAK(ID,NHEP0)
1317 CALL PHOMAK(ID,NHEP0)
1318 CALL PHOMAK(ID,NHEP0)
1319 CALL PHOMAK(ID,NHEP0)
1320.GE.
ELSEIF (RN17.D0/24D0) THEN
1321 CALL PHOMAK(ID,NHEP0)
1322 CALL PHOMAK(ID,NHEP0)
1323.GE.
ELSEIF (RN9.D0/24D0) THEN
1324 CALL PHOMAK(ID,NHEP0)
1327C-- triple photon emission
1330.GE.
IF (RN5.D0/6D0) THEN
1331 CALL PHOMAK(ID,NHEP0)
1332 CALL PHOMAK(ID,NHEP0)
1333 CALL PHOMAK(ID,NHEP0)
1334.GE.
ELSEIF (RN2.D0/6D0) THEN
1335 CALL PHOMAK(ID,NHEP0)
1338C-- double photon emission
1341.GE.
IF (RN0.5D0) THEN
1342 CALL PHOMAK(ID,NHEP0)
1343 CALL PHOMAK(ID,NHEP0)
1346C-- single photon emission
1348 CALL PHOMAK(ID,NHEP0)
1351C-- electron positron pair (coomented out for a while
1352C IF (IPAIR) CALL PHOPAR(ID,NHEP0)
1354 SUBROUTINE PHOMAK(IPPAR,NHEP0)
1355C.----------------------------------------------------------------------
1357C. PHOMAK: PHOtos MAKe
1359C. Purpose: Single or double bremstrahlung radiative corrections
1360C. are generated in the decay of the IPPAR-th particle in
1361C. the HEP common /PH_HEPEVT/. Example of the use of
1364C. Input Parameter: IPPAR: Pointer to decaying particle in
1365C. /PH_HEPEVT/ and the common itself
1367C. Output Parameters: Common /PH_HEPEVT/, either with or without
1370C. Author(s): Z. Was, Created at: 26/05/93
1371C. Last Update: 29/01/05
1373C.----------------------------------------------------------------------
1375 DOUBLE PRECISION DATA
1377 INTEGER IP,IPPAR,NCHARG
1378 INTEGER WTDUM,IDUM,NHEP0
1379 INTEGER NCHARB,NEUDAU
1383 PARAMETER (NMXHEP=10000)
1384 INTEGER IDHEP,ISTHEP,JDAHEP,JMOHEP,NEVHEP,NHEP
1386 COMMON/PH_HEPEVT/NEVHEP,NHEP,ISTHEP(NMXHEP),IDHEP(NMXHEP),
1387 &JMOHEP(2,NMXHEP),JDAHEP(2,NMXHEP),PHEP(5,NMXHEP),VHEP(4,NMXHEP)
1388 LOGICAL INTERF,ISEC,ITRE,IEXP,IFTOP,IFW
1389 REAL*8 FINT,FSEC,EXPEPS
1390 COMMON /PHOKEY/ FSEC,FINT,EXPEPS,INTERF,ISEC,ITRE,IEXP,IFTOP,IFW
1396 CALL PHOIN(IP,BOOST,NHEP0)
1397 CALL PHOCHK(JDAHEP(1,IP))
1399 CALL PHOPRE(1,WT,NEUDAU,NCHARB)
1401.EQ.
IF (WT0.0D0) RETURN
1403C PHODO is caling PHORAN, thus change of series if it is moved before if
1404 CALL PHODO(1,NCHARB,NEUDAU)
1405C we eliminate /FINT in variant B.
1406 IF (INTERF) WT=WT*PHINT(IDUM) /FINT ! FINT must be in variant A
1407 IF (IFW) CALL PHOBW(WT) ! extra weight for leptonic W decay
1409.GT.
IF (WT1.0D0) CALL PHOERR(3,'wt_int
',DATA)
1412 CALL PHOOUT(IP,BOOST,NHEP0)
1416 FUNCTION PHINT1(IDUM)
1417C.----------------------------------------------------------------------
1419C. PHINT: PHotos INTerference (Old version kept for tests only.
1421C. Purpose: Calculates interference between emission of photons from
1422C. different possible chaged daughters stored in
1423C. the HEP common /PHOEVT/.
1425C. Input Parameter: commons /PHOEVT/ /PHOMOM/ /PHOPHS/
1428C. Output Parameters:
1431C. Author(s): Z. Was, Created at: 10/08/93
1432C. Last Update: 15/03/99
1434C.----------------------------------------------------------------------
1441 PARAMETER (NMXPHO=10000)
1442 INTEGER IDPHO,ISTPHO,JDAPHO,JMOPHO,NEVPHO,NPHO
1444 COMMON/PHOEVT/NEVPHO,NPHO,ISTPHO(NMXPHO),IDPHO(NMXPHO),
1445 &JMOPHO(2,NMXPHO),JDAPHO(2,NMXPHO),PPHO(5,NMXPHO),VPHO(4,NMXPHO)
1446 DOUBLE PRECISION MCHSQR,MNESQR
1448 COMMON/PHOMOM/MCHSQR,MNESQR,PNEUTR(5)
1449 DOUBLE PRECISION COSTHG,SINTHG
1450 REAL*8 XPHMAX,XPHOTO
1451 COMMON/PHOPHS/XPHMAX,XPHOTO,COSTHG,SINTHG
1452 REAL*8 MPASQR,XX,BETA
1456 DO K=JDAPHO(2,1),JDAPHO(1,1),-1
1457.NE.
IF(IDPHO(K)22) THEN
1463C check if there is a photon
1464.GT.
IFINT= NPHOIDENT
1465C check if it is two body + gammas reaction
1466.AND..EQ.
IFINT= IFINT(IDENT-JDAPHO(1,1))1
1467C check if two body was particle antiparticle
1468.AND..EQ.
IFINT= IFINTIDPHO(JDAPHO(1,1))-IDPHO(IDENT)
1469C check if particles were charged
1470.AND..NE.
IFINT= IFINTPHOCHA(IDPHO(IDENT))0
1471C calculates interference weight contribution
1473 MPASQR = PPHO(5,1)**2
1474 XX=4.D0*MCHSQR/MPASQR*(1.D0-XPHOTO)/(1.D0-XPHOTO+(MCHSQR-MNESQR)
1477 PHINT = 2D0/(1D0+COSTHG**2*BETA**2)
1484 FUNCTION PHINT2(IDUM)
1485C.----------------------------------------------------------------------
1487C. PHINT: PHotos INTerference
1489C. Purpose: Calculates interference between emission of photons from
1490C. different possible chaged daughters stored in
1491C. the HEP common /PHOEVT/.
1493C. Input Parameter: commons /PHOEVT/ /PHOMOM/ /PHOPHS/
1496C. Output Parameters:
1499C. Author(s): Z. Was, Created at: 10/08/93
1502C.----------------------------------------------------------------------
1505 REAL*8 PHINT,PHINT1,PHINT2
1509 PARAMETER (NMXPHO=10000)
1510 INTEGER IDPHO,ISTPHO,JDAPHO,JMOPHO,NEVPHO,NPHO
1512 COMMON/PHOEVT/NEVPHO,NPHO,ISTPHO(NMXPHO),IDPHO(NMXPHO),
1513 &JMOPHO(2,NMXPHO),JDAPHO(2,NMXPHO),PPHO(5,NMXPHO),VPHO(4,NMXPHO)
1514 DOUBLE PRECISION MCHSQR,MNESQR
1516 COMMON/PHOMOM/MCHSQR,MNESQR,PNEUTR(5)
1517 DOUBLE PRECISION COSTHG,SINTHG
1518 REAL*8 XPHMAX,XPHOTO
1519 COMMON/PHOPHS/XPHMAX,XPHOTO,COSTHG,SINTHG
1520 REAL*8 MPASQR,XX,BETA,PQ1(4),PQ2(4),PPHOT(4)
1521 REAL*8 SS,PP2,PP,E1,E2,Q1,Q2,COSTHE
1525 DO K=JDAPHO(2,1),JDAPHO(1,1),-1
1526.NE.
IF(IDPHO(K)22) THEN
1532C check if there is a photon
1533.GT.
IFINT= NPHOIDENT
1534C check if it is two body + gammas reaction
1535.AND..EQ.
IFINT= IFINT(IDENT-JDAPHO(1,1))1
1536C check if two body was particle antiparticle (we improve on it !
1537.AND..EQ.
C IFINT= IFINTIDPHO(JDAPHO(1,1))-IDPHO(IDENT)
1538C check if particles were charged
1539.AND..GT.
IFINT= IFINTabs(PHOCHA(IDPHO(IDENT)))0.01D0
1540C check if they have both charge
1542.gt.
$ abs(PHOCHA(IDPHO(JDAPHO(1,1))))0.01D0
1543C calculates interference weight contribution
1545 MPASQR = PPHO(5,1)**2
1546 XX=4.D0*MCHSQR/MPASQR*(1.-XPHOTO)/(1.-XPHOTO+(MCHSQR-MNESQR)/
1549 PHINT = 2D0/(1D0+COSTHG**2*BETA**2)
1550 SS =MPASQR*(1.D0-XPHOTO)
1551 PP2=((SS-MCHSQR-MNESQR)**2-4*MCHSQR*MNESQR)/SS/4
1553 E1 =SQRT(PP2+MCHSQR)
1554 E2 =SQRT(PP2+MNESQR)
1555 PHINT= (E1+E2)**2/((E2+COSTHG*PP)**2+(E1-COSTHG*PP)**2)
1557 q1=PHOCHA(IDPHO(JDAPHO(1,1)))
1558 q2=PHOCHA(IDPHO(IDENT))
1560 pq1(k)=ppho(k,JDAPHO(1,1))
1561 pq2(k)=ppho(k,JDAPHO(1,1)+1)
1562 pphot(k)=ppho(k,npho)
1564 costhe=(pphot(1)*pq1(1)+pphot(2)*pq1(2)+pphot(3)*pq1(3))
1565 costhe=costhe/sqrt(pq1(1)**2+pq1(2)**2+pq1(3)**2)
1566 costhe=costhe/sqrt(pphot(1)**2+pphot(2)**2+pphot(3)**2)
1568! --- this IF checks whether JDAPHO(1,1) was MCH or MNE.
1569! --- COSTHG angle (and in-generation variables) may be better choice
1570! --- than costhe. note that in the formulae below amplitudes were
1571! --- multiplied by (E2+COSTHG*PP)*(E1-COSTHG*PP).
1572.GT.
IF (costhg*costhe0) then
1574 PHINT= (q1*(E2+COSTHG*PP)-q2*(E1-COSTHG*PP))**2
1575 & /(q1**2*(E2+COSTHG*PP)**2+q2**2*(E1-COSTHG*PP)**2)
1578 PHINT= (q1*(E1-COSTHG*PP)-q2*(E2+COSTHG*PP))**2
1579 & /(q1**2*(E1-COSTHG*PP)**2+q2**2*(E2+COSTHG*PP)**2)
1589 SUBROUTINE PHOPRE(IPARR,WT,NEUDAU,NCHARB)
1590C.----------------------------------------------------------------------
1592C. PHOTOS: Photon radiation in decays
1594C. Purpose: Order (alpha) radiative corrections are generated in
1595C. the decay of the IPPAR-th particle in the HEP-like
1596C. common /PHOEVT/. Photon radiation takes place from one
1597C. of the charged daughters of the decaying particle IPPAR
1598C. WT is calculated, eventual rejection will be performed
1599C. later after inclusion of interference weight.
1601C. Input Parameter: IPPAR: Pointer to decaying particle in
1602C. /PHOEVT/ and the common itself,
1604C. Output Parameters: Common /PHOEVT/, either with or without a
1606C. WT weight of the configuration
1608C. Author(s): Z. Was, B. van Eijk Created at: 26/11/89
1609C. Last Update: 29/01/05
1611C.----------------------------------------------------------------------
1613 DOUBLE PRECISION MINMAS,MPASQR,MCHREN
1614 DOUBLE PRECISION BETA,EPS,DEL1,DEL2,DATA,BIGLOG
1615 REAL*8 PHOCHA,PHOSPI,PHORAN,PHOCOR,MASSUM
1616 INTEGER IP,IPARR,IPPAR,I,J,ME,NCHARG,NEUPOI,NLAST,THEDUM
1618 INTEGER NCHARB,NEUDAU
1621 PARAMETER (NMXPHO=10000)
1622 INTEGER IDPHO,ISTPHO,JDAPHO,JMOPHO,NEVPHO,NPHO
1624 COMMON/PHOEVT/NEVPHO,NPHO,ISTPHO(NMXPHO),IDPHO(NMXPHO),
1625 &JMOPHO(2,NMXPHO),JDAPHO(2,NMXPHO),PPHO(5,NMXPHO),VPHO(4,NMXPHO)
1627 COMMON/PHOIF/CHKIF(NMXPHO)
1628 INTEGER CHAPOI(NMXPHO)
1629 DOUBLE PRECISION MCHSQR,MNESQR
1631 COMMON/PHOMOM/MCHSQR,MNESQR,PNEUTR(5)
1632 DOUBLE PRECISION COSTHG,SINTHG
1633 REAL*8 XPHMAX,XPHOTO
1634 COMMON/PHOPHS/XPHMAX,XPHOTO,COSTHG,SINTHG
1636 COMMON/PHOCOP/ALPHA,XPHCUT
1638 REAL*8 PROBH,CORWT,XF
1639 COMMON/PHOPRO/PROBH,CORWT,XF,IREP
1640C may be it is not the best place, but ...
1641 LOGICAL INTERF,ISEC,ITRE,IEXP,IFTOP,IFW
1642 REAL*8 FINT,FSEC,EXPEPS
1643 COMMON /PHOKEY/ FSEC,FINT,EXPEPS,INTERF,ISEC,ITRE,IEXP,IFTOP,IFW
1647C-- Store pointers for cascade treatement...
1652C-- Check decay multiplicity..
1653.EQ.
IF (JDAPHO(1,IP)0) RETURN
1655C-- Loop over daughters, determine charge multiplicity
1660 DO 20 I=JDAPHO(1,IP),JDAPHO(2,IP)
1663C-- Exclude marked particles, quarks and gluons etc...
1665 IF (CHKIF(I-JDAPHO(1,IP)+3)) THEN
1666.NE.
IF (PHOCHA(IDPHO(I))0) THEN
1668.GT.
IF (NCHARGNMXPHO) THEN
1670 CALL PHOERR(1,'photos
',DATA)
1674 MINMAS=MINMAS+PPHO(5,I)**2
1676 MASSUM=MASSUM+PPHO(5,I)
1678.NE.
IF (NCHARG0) THEN
1680C-- Check that sum of daughter masses does not exceed parent mass
1681.GT.
IF ((PPHO(5,IP)-MASSUM)/PPHO(5,IP)2.D0*XPHCUT) THEN
1683C-- Order charged particles according to decreasing mass, this to
1684C-- increase efficiency (smallest mass is treated first).
1685.GT.
IF (NCHARG1) CALL PHOOMA(1,NCHARG,CHAPOI)
1689 70 PNEUTR(J)=-PPHO(J,CHAPOI(NCHARG))
1690 PNEUTR(4)=PPHO(5,IP)-PPHO(4,CHAPOI(NCHARG))
1692C-- Calculate invariant mass of 'neutral
' etc. systems
1693 MPASQR=PPHO(5,IP)**2
1694 MCHSQR=PPHO(5,CHAPOI(NCHARG))**2
1695.EQ.
IF ((JDAPHO(2,IP)-JDAPHO(1,IP))1) THEN
1697.EQ.
IF (NEUPOICHAPOI(NCHARG)) NEUPOI=JDAPHO(2,IP)
1698 MNESQR=PPHO(5,NEUPOI)**2
1699 PNEUTR(5)=PPHO(5,NEUPOI)
1701 MNESQR=PNEUTR(4)**2-PNEUTR(1)**2-PNEUTR(2)**2-PNEUTR(3)**2
1702 MNESQR=MAX(MNESQR,MINMAS-MCHSQR)
1703 PNEUTR(5)=SQRT(MNESQR)
1706C-- Determine kinematical limit...
1707 XPHMAX=(MPASQR-(PNEUTR(5)+PPHO(5,CHAPOI(NCHARG)))**2)/MPASQR
1709C-- Photon energy fraction...
1710 CALL PHOENE(MPASQR,MCHREN,BETA,BIGLOG,IDPHO(CHAPOI(NCHARG)))
1712.LT.
IF (XPHOTO-4D0) THEN
1713 NCHARG=0 ! we really stop trials
1714 XPHOTO=0d0! in this case !!
1715C-- Energy fraction not too large (very seldom) ? Define angle.
1716.LT..OR..GT.
ELSEIF ((XPHOTOXPHCUT)(XPHOTOXPHMAX)) THEN
1718C-- No radiation was accepted, check for more daughters that may ra-
1719C-- diate and correct radiation probability...
1721.GT.
IF (NCHARG0) THEN
1727C-- Angle is generated in the frame defined by charged vector and
1728C-- PNEUTR, distribution is taken in the infrared limit...
1729 EPS=MCHREN/(1.D0+BETA)
1731C-- Calculate sin(theta) and cos(theta) from interval variables
1732 DEL1=(2.D0-EPS)*(EPS/(2.D0-EPS))**PHORAN(THEDUM)
1735C ----------- VARIANT B ------------------
1736CC corrections for more efiicient interference correction,
1737CC instead of doubling crude distribution, we add flat parallel channel
1738.LT.
C IF (PHORAN(THEDUM)BIGLOG/BETA/(BIGLOG/BETA+2*FINT)) THEN
1739C COSTHG=(1.D0-DEL1)/BETA
1740C SINTHG=SQRT(DEL1*DEL2-MCHREN)/BETA
1742C COSTHG=-1D0+2*PHORAN(THEDUM)
1743C SINTHG= SQRT(1D0-COSTHG**2)
1746.GT.
C IF (FINT1.0D0) THEN
1748C WGT=1D0/(1D0-BETA*COSTHG)
1756C ----------- END OF VARIANT B ------------------
1758C ----------- VARIANT A ------------------
1759 COSTHG=(1.D0-DEL1)/BETA
1760 SINTHG=SQRT(DEL1*DEL2-MCHREN)/BETA
1762C ----------- END OF VARIANT A ------------------
1765C-- Determine spin of particle and construct code for matrix element
1766 ME=2.D0*PHOSPI(IDPHO(CHAPOI(NCHARG)))+1.D0
1768C-- Weighting procedure with 'exact
' matrix element, reconstruct kine-
1769C-- matics for photon, neutral and charged system and update /PHOEVT/.
1770C-- Find pointer to the first component of 'neutral
' system
1771 DO I=JDAPHO(1,IP),JDAPHO(2,IP)
1772.NE.
IF (ICHAPOI(NCHARG)) THEN
1778C-- Pointer not found...
1780 CALL PHOERR(5,'phokin
',DATA)
1782 NCHARB=CHAPOI(NCHARG)
1783 NCHARB=NCHARB-JDAPHO(1,IP)+3
1784 NEUDAU=NEUDAU-JDAPHO(1,IP)+3
1785 WT=PHOCOR(MPASQR,MCHREN,ME)*WGT
1788 DATA=PPHO(5,IP)-MASSUM
1789 CALL PHOERR(10,'photos
',DATA)
1795 SUBROUTINE PHOOMA(IFIRST,ILAST,POINTR)
1796C.----------------------------------------------------------------------
1798C. PHOTOS: PHOton radiation in decays Order MAss vector
1800C. Purpose: Order the contents of array 'pointr
' according to the
1801C. decreasing value in the array 'mass
'.
1803C. Input Parameters: IFIRST, ILAST: Pointers to the vector loca-
1805C. POINTR: Unsorted array with pointers to
1808C. Output Parameter: POINTR: Sorted arrays with respect to
1809C. particle mass 'ppho(5,*)
'.
1811C. Author(s): B. van Eijk Created at: 28/11/89
1812C. Last Update: 27/05/93
1814C.----------------------------------------------------------------------
1817 PARAMETER (NMXPHO=10000)
1818 INTEGER IDPHO,ISTPHO,JDAPHO,JMOPHO,NEVPHO,NPHO
1820 COMMON/PHOEVT/NEVPHO,NPHO,ISTPHO(NMXPHO),IDPHO(NMXPHO),
1821 &JMOPHO(2,NMXPHO),JDAPHO(2,NMXPHO),PPHO(5,NMXPHO),VPHO(4,NMXPHO)
1822 INTEGER IFIRST,ILAST,I,J,BUFPOI,POINTR(NMXPHO)
1823 REAL*8 BUFMAS,MASS(NMXPHO)
1824.EQ.
IF (IFIRSTILAST) RETURN
1826C-- Copy particle masses
1827 DO 10 I=IFIRST,ILAST
1828 10 MASS(I)=PPHO(5,POINTR(I))
1830C-- Order the masses in a decreasing series
1831 DO 30 I=IFIRST,ILAST-1
1833.LE.
IF (MASS(J)MASS(I)) GOTO 20
1844 SUBROUTINE PHOENE(MPASQR,MCHREN,BETA,BIGLOG,IDENT)
1845C.----------------------------------------------------------------------
1847C. PHOTOS: PHOton radiation in decays calculation of photon ENErgy
1850C. Purpose: Subroutine returns photon energy fraction (in (parent
1851C. mass)/2 units) for the decay bremsstrahlung.
1853C. Input Parameters: MPASQR: Mass of decaying system squared,
1854C. XPHCUT: Minimum energy fraction of photon,
1855C. XPHMAX: Maximum energy fraction of photon.
1857C. Output Parameter: MCHREN: Renormalised mass squared,
1858C. BETA: Beta factor due to renormalisation,
1859C. XPHOTO: Photon energy fraction,
1860C. XF: Correction factor for PHOFAC.
1862C. Author(s): S. Jadach, Z. Was Created at: 01/01/89
1863C. B. van Eijk, P.Golonka Last Update: 29/01/05
1865C.----------------------------------------------------------------------
1867 DOUBLE PRECISION MPASQR,MCHREN,BIGLOG,BETA,DATA
1868 INTEGER IWT1,IRN,IWT2
1869 REAL*8 PRSOFT,PRHARD,PHORAN,PHOFAC
1870 DOUBLE PRECISION MCHSQR,MNESQR
1873 REAL*8 PHOCHA,PRKILL,RRR
1874 COMMON/PHOMOM/MCHSQR,MNESQR,PNEUTR(5)
1875 DOUBLE PRECISION COSTHG,SINTHG
1876 REAL*8 XPHMAX,XPHOTO
1877 COMMON/PHOPHS/XPHMAX,XPHOTO,COSTHG,SINTHG
1879 COMMON/PHOCOP/ALPHA,XPHCUT
1881 COMMON/PHPICO/PI,TWOPI
1883 REAL*8 PROBH,CORWT,XF
1884 COMMON/PHOPRO/PROBH,CORWT,XF,IREP
1885 LOGICAL INTERF,ISEC,ITRE,IEXP,IFTOP,IFW
1886 REAL*8 FINT,FSEC,EXPEPS
1887 COMMON /PHOKEY/ FSEC,FINT,EXPEPS,INTERF,ISEC,ITRE,IEXP,IFTOP,IFW
1892 COMMON /PHOEXP/ PRO(NX),NCHAN,EXPINI
1894.LE.
IF (XPHMAXXPHCUT) THEN
1895 BETA=PHOFAC(-1) ! to zero counter, here beta is dummy
1899C-- Probabilities for hard and soft bremstrahlung...
1900 MCHREN=4.D0*MCHSQR/MPASQR/(1.D0+MCHSQR/MPASQR)**2
1901 BETA=SQRT(1.D0-MCHREN)
1903C ----------- VARIANT B ------------------
1904CC we replace 1D0/BETA*BIGLOG with (1D0/BETA*BIGLOG+2*FINT)
1905CC for integral of new crude
1906C BIGLOG=LOG(MPASQR/MCHSQR*(1.D0+BETA)**2/4.D0*
1907C & (1.D0+MCHSQR/MPASQR)**2)
1908C PRHARD=ALPHA/PI*(1D0/BETA*BIGLOG+2*FINT)*(LOG(XPHMAX/XPHCUT)
1909C &-.75D0+XPHCUT/XPHMAX-.25D0*XPHCUT**2/XPHMAX**2)
1910C PRHARD=PRHARD*PHOCHA(IDENT)**2*FSEC
1911C ----------- END OF VARIANT B ------------------
1913C ----------- VARIANT A ------------------
1914 BIGLOG=LOG(MPASQR/MCHSQR*(1.D0+BETA)**2/4.D0*
1915 & (1.D0+MCHSQR/MPASQR)**2)
1916 PRHARD=ALPHA/PI*(1D0/BETA*BIGLOG)*
1917 &(LOG(XPHMAX/XPHCUT)-.75D0+XPHCUT/XPHMAX-.25D0*XPHCUT**2/XPHMAX**2)
1918 PRHARD=PRHARD*PHOCHA(IDENT)**2*FSEC*FINT
1919C ----------- END OF VARIANT A ------------------
1920.EQ.
IF (IREP0) PROBH=0.D0
1922 IF (IEXP) THEN ! IEXP
1924 IF (EXPINI) THEN ! EXPINI
1925 PRO(NCHAN)=PRHARD+0.05*(1.0+FINT) ! we store hard photon emission prob
1927 PRHARD=0D0 ! to kill emission at initialization call
1934 PRHARD=PRHARD/PRSUM ! note that PRHARD may be smaller than
1935 !PRO(NCHAN) because it is calculated
1936 ! for kinematical configuartion as is
1937 ! (with effects of previous photons)
1938 PRKILL=PRO(NCHAN)/PRSUM-PRHARD !
1943 PRHARD=PRHARD*PHOFAC(0) ! PHOFAC is used to control eikonal
1944 ! formfactors for non exp version only
1945 ! here PHOFAC(0)=1 at least now.
1950C-- Check on kinematical bounds
1952.LT.
IF (PRSOFT-5.0D-8) THEN
1954 CALL PHOERR(2,'phoene
',DATA)
1957.LT.
IF (PRSOFT0.1D0) THEN
1959 CALL PHOERR(2,'phoene
',DATA)
1964.LT.
IF (RRRPRSOFT) THEN
1966C-- No photon... (ie. photon too soft)
1968.LT.
IF (RRRPRKILL) XPHOTO=-5d0 ! No photon...no further trials
1971C-- Hard photon... (ie. photon hard enough).
1972C-- Calculate Altarelli-Parisi Kernel
1973 10 XPHOTO=EXP(PHORAN(IRN)*LOG(XPHCUT/XPHMAX))
1974 XPHOTO=XPHOTO*XPHMAX
1975.GT.
IF (PHORAN(IWT2)((1.D0+(1.D0-XPHOTO/XPHMAX)**2)/2.D0))
1979C-- Calculate parameter for PHOFAC function
1980 XF=4.D0*MCHSQR*MPASQR/(MPASQR+MCHSQR-MNESQR)**2
1983 FUNCTION PHOCOR(MPASQR,MCHREN,ME)
1984C.----------------------------------------------------------------------
1986C. PHOTOS: PHOton radiation in decays CORrection weight from
1989C. Purpose: Calculate photon angle. The reshaping functions will
1990C. have to depend on the spin S of the charged particle.
1991C. We define: ME = 2 * S + 1 !
1993C. Input Parameters: MPASQR: Parent mass squared,
1994C. MCHREN: Renormalised mass of charged system,
1995C. ME: 2 * spin + 1 determines matrix element
1997C. Output Parameter: Function value.
1999C. Author(s): Z. Was, B. van Eijk Created at: 26/11/89
2000C. Last Update: 21/03/93
2002C.----------------------------------------------------------------------
2004 DOUBLE PRECISION MPASQR,MCHREN,BETA,XX,YY,DATA
2006 REAL*8 PHOCOR,PHOFAC,WT1,WT2,WT3
2007 DOUBLE PRECISION MCHSQR,MNESQR
2009 COMMON/PHOMOM/MCHSQR,MNESQR,PNEUTR(5)
2010 DOUBLE PRECISION COSTHG,SINTHG
2011 REAL*8 XPHMAX,XPHOTO
2012 COMMON/PHOPHS/XPHMAX,XPHOTO,COSTHG,SINTHG
2014 REAL*8 PROBH,CORWT,XF
2015 COMMON/PHOPRO/PROBH,CORWT,XF,IREP
2017C-- Shaping (modified by ZW)...
2018 XX=4.D0*MCHSQR/MPASQR*(1.D0-XPHOTO)/(1.D0-XPHOTO+(MCHSQR-MNESQR)/
2022 WT3=(1.D0-XPHOTO/XPHMAX)/((1.D0+(1.D0-XPHOTO/XPHMAX)**2)/2.D0)
2023.EQ.
ELSEIF (ME2) THEN
2024 YY=0.5D0*(1.D0-XPHOTO/XPHMAX+1.D0/(1.D0-XPHOTO/XPHMAX))
2026.EQ..OR..EQ..OR..EQ.
ELSEIF ((ME3)(ME4)(ME5)) THEN
2028 WT3=(1.D0+(1.D0-XPHOTO/XPHMAX)**2-(XPHOTO/XPHMAX)**3)/
2029 & (1.D0+(1.D0-XPHOTO/XPHMAX)** 2)
2032 CALL PHOERR(6,'phocor
',DATA)
2037 WT1=(1.D0-COSTHG*SQRT(1.D0-MCHREN))/(1.D0-COSTHG*BETA)
2038 WT2=(1.D0-XX/YY/(1.D0-BETA**2*COSTHG**2))*(1.D0+COSTHG*BETA)/2.D0
2042.GT.
IF (PHOCOR1.D0) THEN
2044 CALL PHOERR(3,'phocor
',DATA)
2048 FUNCTION PHOFAC(MODE)
2049C.----------------------------------------------------------------------
2051C. PHOTOS: PHOton radiation in decays control FACtor
2053C. Purpose: This is the control function for the photon spectrum and
2054C. final weighting. It is called from PHOENE for genera-
2055C. ting the raw photon energy spectrum (MODE=0) and in PHO-
2056C. COR to scale the final weight (MODE=1). The factor con-
2057C. sists of 3 terms. Addition of the factor FF which mul-
2058C. tiplies PHOFAC for MODE=0 and divides PHOFAC for MODE=1,
2059C. does not affect the results for the MC generation. An
2060C. appropriate choice for FF can speed up the calculation.
2061C. Note that a too small value of FF may cause weight over-
2062C. flow in PHOCOR and will generate a warning, halting the
2063C. execution. PRX should be included for repeated calls
2064C. for the same event, allowing more particles to radiate
2065C. photons. At the first call IREP=0, for more than 1
2066C. charged decay products, IREP >= 1. Thus, PRSOFT (no
2067C. photon radiation probability in the previous calls)
2068C. appropriately scales the strength of the bremsstrahlung.
2070C. Input Parameters: MODE, PROBH, XF
2072C. Output Parameter: Function value
2074C. Author(s): S. Jadach, Z. Was Created at: 01/01/89
2075C. B. van Eijk, P.Golonka Last Update: 26/06/04
2077C.----------------------------------------------------------------------
2079 REAL*8 PHOFAC,FF,PRX
2082 REAL*8 PROBH,CORWT,XF
2083 COMMON/PHOPRO/PROBH,CORWT,XF,IREP
2084 LOGICAL INTERF,ISEC,ITRE,IEXP,IFTOP,IFW
2085 REAL*8 FINT,FSEC,EXPEPS
2086 COMMON /PHOKEY/ FSEC,FINT,EXPEPS,INTERF,ISEC,ITRE,IEXP,IFTOP,IFW
2088 DATA PRX,FF/ 0.D0, 0.D0/
2089 IF (IEXP) THEN ! In case of exponentiation this routine is useles
2093.EQ.
IF (MODE-1) THEN
2097.EQ.
ELSEIF (MODE0) THEN
2098.EQ.
IF (IREP0) PRX=1.D0
2099 PRX=PRX/(1.D0-PROBH)
2102C-- Following options are not considered for the time being...
2103C-- (1) Good choice, but does not save very much time:
2104C-- FF=(1.0D0-SQRT(XF)/2.0D0)/(1.0+SQRT(XF)/2.0D0)
2105C-- (2) Taken from the blue, but works without weight overflows...
2106C-- FF=(1.D0-XF/(1-(1-SQRT(XF))**2))*(1+(1-SQRT(XF))/SQRT(1-XF))/2
2112 SUBROUTINE PHOBW(WT)
2113C.----------------------------------------------------------------------
2115C. PHOTOS: PHOtos Boson W correction weight
2117C. Purpose: calculates correction weight due to amplitudes of
2118C. emission from W boson.
2124C. Input Parameters: Common /PHOEVT/, with photon added.
2125C. wt to be corrected
2129C. Output Parameters: wt
2131C. Author(s): G. Nanava, Z. Was Created at: 13/03/03
2132C. Last Update: 13/03/03
2134C.----------------------------------------------------------------------
2138 PARAMETER (NMXPHO=10000)
2139 INTEGER IDPHO,ISTPHO,JDAPHO,JMOPHO,NEVPHO,NPHO
2141 COMMON/PHOEVT/NEVPHO,NPHO,ISTPHO(NMXPHO),IDPHO(NMXPHO),
2142 &JMOPHO(2,NMXPHO),JDAPHO(2,NMXPHO),PPHO(5,NMXPHO),VPHO(4,NMXPHO)
2144 DOUBLE PRECISION EMU,MCHREN,BETA,COSTHG,MPASQR,XPH
2146.EQ..AND.
IF (ABS(IDPHO(1))24
2147.GE..AND.
$ ABS(IDPHO(JDAPHO(1,1) ))11
2148.LE..AND.
$ ABS(IDPHO(JDAPHO(1,1) ))16
2149.GE..AND.
$ ABS(IDPHO(JDAPHO(1,1)+1))11
2150.LE.
$ ABS(IDPHO(JDAPHO(1,1)+1))16 ) THEN
2153.EQ..OR.
$ ABS(IDPHO(JDAPHO(1,1) ))11
2154.EQ..OR.
$ ABS(IDPHO(JDAPHO(1,1) ))13
2155.EQ.
$ ABS(IDPHO(JDAPHO(1,1) ))15 ) THEN
2161 MCHREN=ABS(PPHO(4,I)**2-PPHO(3,I)**2
2162 $ -PPHO(2,I)**2-PPHO(1,I)**2)
2163 BETA=SQRT(1- MCHREN/ PPHO(4,I)**2)
2164 COSTHG=(PPHO(3,I)*PPHO(3,NPHO)+PPHO(2,I)*PPHO(2,NPHO)
2165 $ +PPHO(1,I)*PPHO(1,NPHO))/
2166 $ SQRT(PPHO(3,I)**2+PPHO(2,I)**2+PPHO(1,I)**2) /
2167 $ SQRT(PPHO(3,NPHO)**2+PPHO(2,NPHO)**2+PPHO(1,NPHO)**2)
2170 WT=WT*(1-8*EMU*XPH*(1-COSTHG*BETA)*
2171 $ (MCHREN+2*XPH*SQRT(MPASQR))/
2172 $ MPASQR**2/(1-MCHREN/MPASQR)/(4-MCHREN/MPASQR))
2174c write(*,*) IDPHO(1),IDPHO(JDAPHO(1,1)),IDPHO(JDAPHO(1,1)+1)
2175c write(*,*) emu,xph,costhg,beta,mpasqr,mchren
2178 SUBROUTINE PHODO(IP,NCHARB,NEUDAU)
2179C.----------------------------------------------------------------------
2181C. PHOTOS: PHOton radiation in decays DOing of KINematics
2183C. Purpose: Starting from the charged particle energy/momentum,
2184C. PNEUTR, photon energy fraction and photon angle with
2185C. respect to the axis formed by charged particle energy/
2186C. momentum vector and PNEUTR, scale the energy/momentum,
2187C. keeping the original direction of the neutral system in
2188C. the lab. frame untouched.
2190C. Input Parameters: IP: Pointer to decaying particle in
2191C. /PHOEVT/ and the common itself
2192C. NCHARB: pointer to the charged radiating
2193C. daughter in /PHOEVT/.
2194C. NEUDAU: pointer to the first neutral daughter
2195C. Output Parameters: Common /PHOEVT/, with photon added.
2197C. Author(s): Z. Was, B. van Eijk Created at: 26/11/89
2198C. Last Update: 27/05/93
2200C.----------------------------------------------------------------------
2202 DOUBLE PRECISION PHOAN1,PHOAN2,ANGLE,FI1,FI3,FI4,FI5,TH1,TH3,TH4
2203 DOUBLE PRECISION PARNE,QNEW,QOLD,DATA
2204 INTEGER IP,FI3DUM,I,J,NEUDAU,FIRST,LAST
2206 REAL*8 EPHOTO,PMAVIR,PHOTRI
2207 REAL*8 GNEUT,PHORAN,CCOSTH,SSINTH,PVEC(4)
2209 PARAMETER (NMXPHO=10000)
2210 INTEGER IDPHO,ISTPHO,JDAPHO,JMOPHO,NEVPHO,NPHO
2212 COMMON/PHOEVT/NEVPHO,NPHO,ISTPHO(NMXPHO),IDPHO(NMXPHO),
2213 &JMOPHO(2,NMXPHO),JDAPHO(2,NMXPHO),PPHO(5,NMXPHO),VPHO(4,NMXPHO)
2214 DOUBLE PRECISION MCHSQR,MNESQR
2216 COMMON/PHOMOM/MCHSQR,MNESQR,PNEUTR(5)
2217 DOUBLE PRECISION COSTHG,SINTHG
2218 REAL*8 XPHMAX,XPHOTO
2219 COMMON/PHOPHS/XPHMAX,XPHOTO,COSTHG,SINTHG
2221 COMMON/PHPICO/PI,TWOPI
2223 EPHOTO=XPHOTO*PPHO(5,IP)/2.D0
2224 PMAVIR=SQRT(PPHO(5,IP)*(PPHO(5,IP)-2.D0*EPHOTO))
2226C-- Reconstruct kinematics of charged particle and neutral system
2227 FI1=PHOAN1(PNEUTR(1),PNEUTR(2))
2229C-- Choose axis along z of PNEUTR, calculate angle between x and y
2230C-- components and z and x-y plane and perform Lorentz transform...
2231 TH1=PHOAN2(PNEUTR(3),SQRT(PNEUTR(1)**2+PNEUTR(2)**2))
2232 CALL PHORO3(-FI1,PNEUTR(1))
2233 CALL PHORO2(-TH1,PNEUTR(1))
2235C-- Take away photon energy from charged particle and PNEUTR ! Thus
2236C-- the onshell charged particle decays into virtual charged particle
2237C-- and photon. The virtual charged particle mass becomes:
2238C-- SQRT(PPHO(5,IP)*(PPHO(5,IP)-2*EPHOTO)). Construct new PNEUTR mo-
2239C-- mentum in the rest frame of the parent:
2240C-- 1) Scaling parameters...
2241 QNEW=PHOTRI(PMAVIR,PNEUTR(5),PPHO(5,NCHARB))
2243 GNEUT=(QNEW**2+QOLD**2+MNESQR)/(QNEW*QOLD+SQRT((QNEW**2+MNESQR)*
2245.LT.
IF (GNEUT1.D0) THEN
2247 CALL PHOERR(4,'phokin
',DATA)
2249 PARNE=GNEUT-SQRT(MAX(GNEUT**2-1.0D0,0.D0))
2251C-- 2) ...reductive boost...
2252 CALL PHOBO3(PARNE,PNEUTR)
2254C-- ...calculate photon energy in the reduced system...
2258C-- Photon mother and daughter pointers !
2263 PPHO(4,NPHO)=EPHOTO*PPHO(5,IP)/PMAVIR
2265C-- ...and photon momenta
2268 TH3=PHOAN2(CCOSTH,SSINTH)
2269 FI3=TWOPI*PHORAN(FI3DUM)
2270 PPHO(1,NPHO)=PPHO(4,NPHO)*SINTHG*COS(FI3)
2271 PPHO(2,NPHO)=PPHO(4,NPHO)*SINTHG*SIN(FI3)
2273C-- Minus sign because axis opposite direction of charged particle !
2274 PPHO(3,NPHO)=-PPHO(4,NPHO)*COSTHG
2277C-- Rotate in order to get photon along z-axis
2278 CALL PHORO3(-FI3,PNEUTR(1))
2279 CALL PHORO3(-FI3,PPHO(1,NPHO))
2280 CALL PHORO2(-TH3,PNEUTR(1))
2281 CALL PHORO2(-TH3,PPHO(1,NPHO))
2282 ANGLE=EPHOTO/PPHO(4,NPHO)
2284C-- Boost to the rest frame of decaying particle
2285 CALL PHOBO3(ANGLE,PNEUTR(1))
2286 CALL PHOBO3(ANGLE,PPHO(1,NPHO))
2288C-- Back in the parent rest frame but PNEUTR not yet oriented !
2289 FI4=PHOAN1(PNEUTR(1),PNEUTR(2))
2290 TH4=PHOAN2(PNEUTR(3),SQRT(PNEUTR(1)**2+PNEUTR(2)**2))
2291 CALL PHORO3(FI4,PNEUTR(1))
2292 CALL PHORO3(FI4,PPHO(1,NPHO))
2297 CALL PHORO3(-FI3,PVEC)
2298 CALL PHORO2(-TH3,PVEC)
2299 CALL PHOBO3(ANGLE,PVEC)
2300 CALL PHORO3(FI4,PVEC)
2301 CALL PHORO2(-TH4,PNEUTR)
2302 CALL PHORO2(-TH4,PPHO(1,NPHO))
2303 CALL PHORO2(-TH4,PVEC)
2304 FI5=PHOAN1(PVEC(1),PVEC(2))
2306C-- Charged particle restores original direction
2307 CALL PHORO3(-FI5,PNEUTR)
2308 CALL PHORO3(-FI5,PPHO(1,NPHO))
2309 CALL PHORO2(TH1,PNEUTR(1))
2310 CALL PHORO2(TH1,PPHO(1,NPHO))
2311 CALL PHORO3(FI1,PNEUTR)
2312 CALL PHORO3(FI1,PPHO(1,NPHO))
2313C-- See whether neutral system has multiplicity larger than 1...
2314.GT.
IF ((JDAPHO(2,IP)-JDAPHO(1,IP))1) THEN
2315C-- Find pointers to components of 'neutral
' system
2320.NE..AND..EQ.
IF (INCHARB(JMOPHO(1,I)IP)) THEN
2322C-- Reconstruct kinematics...
2323 CALL PHORO3(-FI1,PPHO(1,I))
2324 CALL PHORO2(-TH1,PPHO(1,I))
2326C-- ...reductive boost
2327 CALL PHOBO3(PARNE,PPHO(1,I))
2329C-- Rotate in order to get photon along z-axis
2330 CALL PHORO3(-FI3,PPHO(1,I))
2331 CALL PHORO2(-TH3,PPHO(1,I))
2333C-- Boost to the rest frame of decaying particle
2334 CALL PHOBO3(ANGLE,PPHO(1,I))
2336C-- Back in the parent rest-frame but PNEUTR not yet oriented.
2337 CALL PHORO3(FI4,PPHO(1,I))
2338 CALL PHORO2(-TH4,PPHO(1,I))
2340C-- Charged particle restores original direction
2341 CALL PHORO3(-FI5,PPHO(1,I))
2342 CALL PHORO2(TH1,PPHO(1,I))
2343 CALL PHORO3(FI1,PPHO(1,I))
2348C-- ...only one 'neutral
' particle in addition to photon!
2350 80 PPHO(J,NEUDAU)=PNEUTR(J)
2353C-- All 'neutrals
' treated, fill /PHOEVT/ for charged particle...
2355 90 PPHO(J,NCHARB)=-(PPHO(J,NPHO)+PNEUTR(J))
2356 PPHO(4,NCHARB)=PPHO(5,IP)-(PPHO(4,NPHO)+PNEUTR(4))
2359 FUNCTION PHOTRI(A,B,C)
2360C.----------------------------------------------------------------------
2362C. PHOTOS: PHOton radiation in decays calculation of TRIangle fie
2364C. Purpose: Calculation of triangle function for phase space.
2366C. Input Parameters: A, B, C (Virtual) particle masses.
2368C. Output Parameter: Function value =
2369C. SQRT(LAMBDA(A**2,B**2,C**2))/(2*A)
2371C. Author(s): B. van Eijk Created at: 15/11/89
2372C. Last Update: 02/01/90
2374C.----------------------------------------------------------------------
2376 DOUBLE PRECISION DA,DB,DC,DAPB,DAMB,DTRIAN
2383 DTRIAN=SQRT((DAMB-DC)*(DAPB+DC)*(DAMB+DC)*(DAPB-DC))
2384 PHOTRI=DTRIAN/(DA+DA)
2387 FUNCTION PHOAN1(X,Y)
2388C.----------------------------------------------------------------------
2390C. PHOTOS: PHOton radiation in decays calculation of ANgle '1
'
2392C. Purpose: Calculate angle from X and Y
2394C. Input Parameters: X, Y
2396C. Output Parameter: Function value
2398C. Author(s): S. Jadach Created at: 01/01/89
2399C. B. van Eijk Last Update: 02/01/90
2401C.----------------------------------------------------------------------
2403 DOUBLE PRECISION PHOAN1
2406 COMMON/PHPICO/PI,TWOPI
2407.LT.
IF (ABS(Y)ABS(X)) THEN
2408 PHOAN1=ATAN(ABS(Y/X))
2409.LE.
IF (X0.D0) PHOAN1=PI-PHOAN1
2411 PHOAN1=ACOS(X/SQRT(X**2+Y**2))
2413.LT.
IF (Y0.D0) PHOAN1=TWOPI-PHOAN1
2416 FUNCTION PHOAN2(X,Y)
2417C.----------------------------------------------------------------------
2419C. PHOTOS: PHOton radiation in decays calculation of ANgle '2
'
2421C. Purpose: Calculate angle from X and Y
2423C. Input Parameters: X, Y
2425C. Output Parameter: Function value
2427C. Author(s): S. Jadach Created at: 01/01/89
2428C. B. van Eijk Last Update: 02/01/90
2430C.----------------------------------------------------------------------
2432 DOUBLE PRECISION PHOAN2
2435 COMMON/PHPICO/PI,TWOPI
2436.LT.
IF (ABS(Y)ABS(X)) THEN
2437 PHOAN2=ATAN(ABS(Y/X))
2438.LE.
IF (X0.D0) PHOAN2=PI-PHOAN2
2440 PHOAN2=ACOS(X/SQRT(X**2+Y**2))
2444 SUBROUTINE PHOBO3(ANGLE,PVEC)
2445C.----------------------------------------------------------------------
2447C. PHOTOS: PHOton radiation in decays BOost routine '3
'
2449C. Purpose: Boost vector PVEC along z-axis where ANGLE = EXP(ETA),
2450C. ETA is the hyperbolic velocity.
2452C. Input Parameters: ANGLE, PVEC
2454C. Output Parameter: PVEC
2456C. Author(s): S. Jadach Created at: 01/01/89
2457C. B. van Eijk Last Update: 02/01/90
2459C.----------------------------------------------------------------------
2461 DOUBLE PRECISION QPL,QMI,ANGLE
2463 QPL=(PVEC(4)+PVEC(3))*ANGLE
2464 QMI=(PVEC(4)-PVEC(3))/ANGLE
2465 PVEC(3)=(QPL-QMI)/2.D0
2466 PVEC(4)=(QPL+QMI)/2.D0
2469 SUBROUTINE PHORO2(ANGLE,PVEC)
2470C.----------------------------------------------------------------------
2472C. PHOTOS: PHOton radiation in decays ROtation routine '2
'
2474C. Purpose: Rotate x and z components of vector PVEC around angle
2477C. Input Parameters: ANGLE, PVEC
2479C. Output Parameter: PVEC
2481C. Author(s): S. Jadach Created at: 01/01/89
2482C. B. van Eijk Last Update: 02/01/90
2484C.----------------------------------------------------------------------
2486 DOUBLE PRECISION CS,SN,ANGLE
2488 CS=COS(ANGLE)*PVEC(1)+SIN(ANGLE)*PVEC(3)
2489 SN=-SIN(ANGLE)*PVEC(1)+COS(ANGLE)*PVEC(3)
2494 SUBROUTINE PHORO3(ANGLE,PVEC)
2495C.----------------------------------------------------------------------
2497C. PHOTOS: PHOton radiation in decays ROtation routine '3
'
2499C. Purpose: Rotate x and y components of vector PVEC around angle
2502C. Input Parameters: ANGLE, PVEC
2504C. Output Parameter: PVEC
2506C. Author(s): S. Jadach Created at: 01/01/89
2507C. B. van Eijk Last Update: 02/01/90
2509C.----------------------------------------------------------------------
2511 DOUBLE PRECISION CS,SN,ANGLE
2513 CS=COS(ANGLE)*PVEC(1)-SIN(ANGLE)*PVEC(2)
2514 SN=SIN(ANGLE)*PVEC(1)+COS(ANGLE)*PVEC(2)
2520C.----------------------------------------------------------------------
2522C. PHOTOS: PHOton radiation in decays RANdom number generator init
2524C. Purpose: Initialse PHORAN with the user specified seeds in the
2525C. array ISEED. For details see also: F. James CERN DD-
2526C. Report November 1988.
2528C. Input Parameters: ISEED(*)
2530C. Output Parameters: URAN, CRAN, CDRAN, CMRAN, I97, J97
2532C. Author(s): B. van Eijk and F. James Created at: 27/09/89
2533C. Last Update: 22/02/90
2535C.----------------------------------------------------------------------
2537 DOUBLE PRECISION DATA
2539 INTEGER I,IS1,IS2,IS3,IS4,IS5,J
2540 INTEGER ISEED,I97,J97
2541 REAL*8 URAN,CRAN,CDRAN,CMRAN
2542 COMMON/PHSEED/ISEED(2),I97,J97,URAN(97),CRAN,CDRAN,CMRAN
2544C-- Check value range of seeds
2545.LT..OR..GE.
IF ((ISEED(1)0)(ISEED(1)31328)) THEN
2547 CALL PHOERR(8,'phorin
',DATA)
2549.LT..OR..GE.
IF ((ISEED(2)0)(ISEED(2)30081)) THEN
2551 CALL PHOERR(9,'phorin
',DATA)
2554C-- Calculate Marsaglia and Zaman seeds (by F. James)
2555 IS1=MOD(ISEED(1)/177,177)+2
2556 IS2=MOD(ISEED(1),177)+2
2557 IS3=MOD(ISEED(2)/169,178)+1
2558 IS4=MOD(ISEED(2),169)
2563 IS5=MOD (MOD(IS1*IS2,179)*IS3,179)
2567 IS4=MOD(53*IS4+1,169)
2568.GE.
IF (MOD(IS4*IS5,64)32) S=S+T
2571 CRAN=362436.D0/16777216.D0
2572 CDRAN=7654321.D0/16777216.D0
2573 CMRAN=16777213.D0/16777216.D0
2578 FUNCTION PHORAN(IDUM)
2579C.----------------------------------------------------------------------
2581C. PHOTOS: PHOton radiation in decays RANdom number generator based
2582C. on Marsaglia Algorithm
2584C. Purpose: Generate uniformly distributed random numbers between
2585C. 0 and 1. Super long period: 2**144. See also:
2586C. G. Marsaglia and A. Zaman, FSU-SCR-87-50, for seed mo-
2587C. difications to this version see: F. James DD-Report,
2588C. November 1988. The generator has to be initialized by
2591C. Input Parameters: IDUM (integer dummy)
2593C. Output Parameters: Function value
2595C. Author(s): B. van Eijk, G. Marsaglia and Created at: 27/09/89
2596C. A. Zaman Last Update: 27/09/89
2598C.----------------------------------------------------------------------
2602 INTEGER ISEED,I97,J97
2603 REAL*8 URAN,CRAN,CDRAN,CMRAN
2604 COMMON/PHSEED/ISEED(2),I97,J97,URAN(97),CRAN,CDRAN,CMRAN
2605 10 PHORAN=URAN(I97)-URAN(J97)
2606.LT.
IF (PHORAN0.D0) PHORAN=PHORAN+1.D0
2609.EQ.
IF (I970) I97=97
2611.EQ.
IF (J970) J97=97
2613.LT.
IF (CRAN0.D0) CRAN=CRAN+CMRAN
2615.LT.
IF (PHORAN0.D0) PHORAN=PHORAN+1.D0
2616.LE.
IF (PHORAN0.D0) GOTO 10
2619 FUNCTION PHOCHA(IDHEP)
2620C.----------------------------------------------------------------------
2622C. PHOTOS: PHOton radiation in decays CHArge determination
2624C. Purpose: Calculate the charge of particle with code IDHEP. The
2625C. code of the particle is defined by the Particle Data
2626C. Group in Phys. Lett. B204 (1988) 1.
2628C. Input Parameter: IDHEP
2630C. Output Parameter: Funtion value = charge of particle with code
2633C. Author(s): E. Barberio and B. van Eijk Created at: 29/11/89
2634C. Last update: 02/01/90
2636C.----------------------------------------------------------------------
2639 INTEGER IDHEP,IDABS,Q1,Q2,Q3
2641C-- Array 'charge
' contains the charge of the first 101 particles ac-
2642C-- cording to the PDG particle code... (0 is added for convenience)
2643 REAL*8 CHARGE(0:100)
2645 &-0.3333333333D0, 0.6666666667D0, -0.3333333333D0, 0.6666666667D0,
2646 &-0.3333333333D0, 0.6666666667D0, -0.3333333333D0, 0.6666666667D0,
2647 & 2*0.D0, -1.D0, 0.D0, -1.D0, 0.D0, -1.D0, 0.D0, -1.D0, 6*0.D0,
2648 & 1.D0, 12*0.D0, 1.D0, 63*0.D0/
2650.LE.
IF (IDABS100) THEN
2652C-- Charge of quark, lepton, boson etc....
2653 PHOCHA = CHARGE(IDABS)
2656C-- Check on particle build out of quarks, unpack its code...
2657 Q3=MOD(IDABS/1000,10)
2658 Q2=MOD(IDABS/100,10)
2663.EQ.
IF(MOD(Q2,2)0) THEN
2664 PHOCHA=CHARGE(Q2)-CHARGE(Q1)
2666 PHOCHA=CHARGE(Q1)-CHARGE(Q2)
2670C-- ...diquarks or baryon.
2671 PHOCHA=CHARGE(Q1)+CHARGE(Q2)+CHARGE(Q3)
2675C-- Find the sign of the charge...
2676.LT.
IF (IDHEP0.D0) PHOCHA=-PHOCHA
2677.lt.
IF (PHOCHA**21d-6) PHOCHA=0.D0
2680 FUNCTION PHOSPI(IDHEP)
2681C.----------------------------------------------------------------------
2683C. PHOTOS: PHOton radiation in decays function for SPIn determina-
2686C. Purpose: Calculate the spin of particle with code IDHEP. The
2687C. code of the particle is defined by the Particle Data
2688C. Group in Phys. Lett. B204 (1988) 1.
2690C. Input Parameter: IDHEP
2692C. Output Parameter: Funtion value = spin of particle with code
2695C. Author(s): E. Barberio and B. van Eijk Created at: 29/11/89
2696C. Last update: 02/01/90
2698C.----------------------------------------------------------------------
2703C-- Array 'spin
' contains the spin of the first 100 particles accor-
2704C-- ding to the PDG particle code...
2706 DATA SPIN/ 8*.5D0, 1.D0, 0.D0, 8*.5D0, 2*0.D0, 4*1.D0, 76*0.D0/
2709C-- Spin of quark, lepton, boson etc....
2710.LE.
IF (IDABS100) THEN
2714C-- ...other particles, however...
2715 PHOSPI=(MOD(IDABS,10)-1.D0)/2.D0
2717C-- ...K_short and K_long are special !!
2718 PHOSPI=MAX(PHOSPI,0.D0)
2722 SUBROUTINE PHOERR(IMES,TEXT,DATA)
2723C.----------------------------------------------------------------------
2725C. PHOTOS: PHOton radiation in decays ERRror handling
2727C. Purpose: Inform user about (fatal) errors and warnings generated
2728C. by either the user or the program.
2730C. Input Parameters: IMES, TEXT, DATA
2732C. Output Parameters: None
2734C. Author(s): B. van Eijk Created at: 29/11/89
2735C. Last Update: 10/01/92
2737C.----------------------------------------------------------------------
2739 DOUBLE PRECISION DATA
2745 PARAMETER (PHOMES=10)
2747 COMMON/PHOSTA/STATUS(PHOMES)
2750C-- security STOP switch
2755.LE.
IF (IMESPHOMES) STATUS(IMES)=STATUS(IMES)+1
2757C-- Count number of non-fatal errors...
2758.EQ..AND..GE.
IF ((IMES 6)(STATUS(IMES)2)) RETURN
2759.EQ..AND..GE.
IF ((IMES10)(STATUS(IMES)2)) RETURN
2763 GOTO (10,20,30,40,50,60,70,80,90,100),IMES
2764 WRITE(PHLUN,9130) IMES
2766 10 WRITE(PHLUN,9010) TEXT,INT(SDATA)
2768 20 WRITE(PHLUN,9020) TEXT,SDATA
2770 30 WRITE(PHLUN,9030) TEXT,SDATA
2772 40 WRITE(PHLUN,9040) TEXT
2774 50 WRITE(PHLUN,9050) TEXT,INT(SDATA)
2776 60 WRITE(PHLUN,9060) TEXT,SDATA
2778 70 WRITE(PHLUN,9070) TEXT,INT(SDATA)
2780 80 WRITE(PHLUN,9080) TEXT,INT(SDATA)
2782 90 WRITE(PHLUN,9090) TEXT,INT(SDATA)
2784 100 WRITE(PHLUN,9100) TEXT,SDATA
2796.GE.
IF (IERROR10) THEN
2806 130 WRITE(PHLUN,9120)
2809 9000 FORMAT(1H ,80('*
'))
2810 9010 FORMAT(1H ,'*
',A,': too many charged particles, ncharg =
',I6,T81,
2812 9020 FORMAT(1H ,'*
',A,': too much bremsstrahlung required, prsoft =
',
2814 9030 FORMAT(1H ,'*
',A,': combined weight is exceeding 1., weight =
',
2816 9040 FORMAT(1H ,'*
',A,
2817 &': error in rescaling charged and neutral vectors
',T81,'*
')
2818 9050 FORMAT(1H ,'*
',A,
2819 &': non matching charged particle
Pointer, ncharg =
',I5,T81,'*
')
2820 9060 FORMAT(1H ,'*
',A,
2821 &':
Do you really work with a particle of spin:
',F4.1,' ?
',T81,
2823 9070 FORMAT(1H ,'*
',A, ': stack length exceeded, nstack =
',I5 ,T81,
2825 9080 FORMAT(1H ,'*
',A,
2826 &': random number generator seed(1) out of range:
',I8,T81,'*
')
2827 9090 FORMAT(1H ,'*
',A,
2828 &': random number generator seed(2) out of range:
',I8,T81,'*
')
2829 9100 FORMAT(1H ,'*
',A,
2830 &': available phase space below cut-off:
',F15.6,' gev/c^2
',T81,
2832 9120 FORMAT(1H ,'*
',T81,'*
')
2833 9130 FORMAT(1H ,'* funny error message:
',I4,'
2834 9140
FORMAT(1h ,
'* Fatal Error Message, I stop this Run !',t81,
'*')
2835 9150
FORMAT(1h ,
'* 10 Error Messages generated, I stop this Run !',t81,
2839c.----------------------------------------------------------------------
2841c. photos: photon radiation in decays run summary report
2843c. purpose: inform user about success and/or restrictions of photos
2844c. encountered during execution.
2846c. input parameters:
Common /phosta/
2848c. output parameters: none
2850c. author(s): b. van eijk created at: 10/01/92
2851c. last update: 10/01/92
2853c.----------------------------------------------------------------------
2858 parameter(phomes=10)
2860 common/phosta/status(phomes)
2872 IF (status(i).EQ.0)
GOTO 10
2873 IF ((i.EQ.6).OR.(i.EQ.10))
THEN
2874 WRITE(phlun,9050) i,status(i)
2877 WRITE(phlun,9060) i,status(i)
2880 IF (.NOT.error)
WRITE(phlun,9070)
2885 9010
FORMAT(1h ,80(
'*'))
2886 9020
FORMAT(1h ,
'*',t81,
'*')
2887 9030
FORMAT(1h ,
'*',26x,25(
'='),t81,
'*')
2888 9040
FORMAT(1h ,
'*',30x,
'PHOTOS Run Summary',t81,
'*')
2889 9050
FORMAT(1h ,
'*',22x,
'Warning #',i2,
' occured',i6,
' times',t81,
'*')
2890 9060
FORMAT(1h ,
'*',23x,
'Error #',i2,
' occured',i6,
' times',t81,
'*')
2891 9070
FORMAT(1h ,
'*',16x,
'PHOTOS Execution has successfully terminated',
2894 SUBROUTINE phlupa(IPOINT)
2896c.----------------------------------------------------------------------
2898c. phlupa: debugging tool
2900c. purpose: none, eventually may printout content of the
2903c. input parameters:
Common /phoevt/ and /phnum/
2904c. latter may have number of the event.
2906c. output parameters: none
2908c. author(s): z. was created at: 30/05/93
2909c. last update: 09/10/05
2911c.----------------------------------------------------------------------
2913 parameter(nmxpho=10000)
2914 INTEGER IDPHO,ISTPHO,JDAPHO,JMOPHO,NEVPHO,NPHO,I,J,IPOINT
2915 INTEGER IPOIN,IPOIN0,IPOINM,IEV
2917 real*8 ppho,vpho,sum
2918 common/phoevt/nevpho,npho,istpho(nmxpho),idpho(nmxpho),
2919 &jmopho(2,nmxpho),jdapho(2,nmxpho),ppho(5,nmxpho),vpho(4,nmxpho)
2925 COMMON /phlupy/ ipoin,ipoinm
2927 IF (ipoin0.LT.0)
THEN
2932 IF (ipoint.LE.ipoinm.OR.ipoint.GE.ipoin )
RETURN
2934 IF (iev.LT.1000)
THEN
2938 WRITE(phlun,*)
'EVENT NR=',iev,
2939 $
'WE ARE TESTING /PHOEVT/ at IPOINT=',ipoint
2942 WRITE(phlun,20) idpho(i),ppho(1,i),ppho(2,i),ppho(3,i),
2943 $ ppho(4,i),ppho(5,i),jdapho(1,i),jdapho(2,i)
2945 WRITE(phlun,20) idpho(i),ppho(1,i),ppho(2,i),ppho(3,i),
2946 $ ppho(4,i),ppho(5,i),jdapho(1,i),jdapho(2,i)
2949 WRITE(phlun,20) idpho(i),ppho(1,i),ppho(2,i),ppho(3,i),
2950 $ ppho(4,i),ppho(5,i),jmopho(1,i),jmopho(2,i)
2952 sum(j)=sum(j)+ppho(j,i)
2955 sum(5)=sqrt(abs(sum(4)**2-sum(1)**2-sum(2)**2-sum(3)**2))
2957 10
FORMAT(1x,
' ID ',
'p_x ',
'p_y ',
'p_z ',
2959 $
'ID-MO_DA1',
'ID-MO DA2' )
2960 20
FORMAT(1x,i4,5(f14.9),2i9)
2961 30
FORMAT(1x,
' SUM',5(f14.9))
2967 FUNCTION iphqrk(MODCOR)
2969c.----------------------------------------------------------------------
2971c. iphqrk: enables blocks emision from quarks
2974c. input parameters: modcor
2975c. modcor >0
type of action
2978c. =0 execution mode(retrns stored
value)
2981c. author(s): z. was created at: 11/12/00
2983c.----------------------------------------------------------------------
2984 INTEGER IPHQRK,MODCOR,MODOP
2990 IF (modcor.NE.0)
THEN
2995 $
'Message from PHOTOS: IPHQRK(MODCOR):: (re)initialization'
2996 IF (modop.EQ.1)
THEN
2998 $
'MODOP=1 -- blocks emission from light quarks: DEFAULT'
2999 ELSEIF (modop.EQ.2)
THEN
3001 $
'MODOP=2 -- enables emission from light quarks: TEST '
3003 WRITE(phlun,*)
'IPHQRK wrong MODCOR=',modcor
3009 IF (modop.EQ.0.AND.modcor.EQ.0)
THEN
3010 WRITE(phlun,*)
'IPHQRK lack of initialization'
3017 FUNCTION iphekl(MODCOR)
3019c.----------------------------------------------------------------------
3021c. iphekl: enables/blocks emision in: pi0 to gamma e+ e-
3024c. input parameters: modcor
3025c. modcor >0
type of action
3028c. =0 execution mode(retrns stored
value)
3031c. author(s): z. was created at: 11/12/00
3033c.----------------------------------------------------------------------
3034 INTEGER IPHEKL,MODCOR,MODOP
3041 IF (modcor.NE.0)
THEN
3046 $
'Message from PHOTOS: IPHEKL(MODCOR):: (re)initialization'
3047 IF (modop.EQ.2)
THEN
3049 $
'MODOP=2 -- blocks emission in pi0 to gamma e+e-: DEFAULT'
3051 $
'MODOP=2 -- blocks emission in Kl to gamma e+e-: DEFAULT'
3052 ELSEIF (modop.EQ.1)
THEN
3054 $
'MODOP=1 -- enables emission in pi0 to gamma e+e- : TEST '
3056 $
'MODOP=1 -- enables emission in Kl to gamma e+e- : TEST '
3058 WRITE(phlun,*)
'IPHEKL wrong MODCOR=',modcor
3064 IF (modop.EQ.0.AND.modcor.EQ.0)
THEN
3065 WRITE(phlun,*)
'IPHELK lack of initialization'
3071 SUBROUTINE phcork(MODCOR)
3073c.----------------------------------------------------------------------
3075c. phcork: corrects kinmatics of subbranch needed
if host
program
3076c. produces events with the shaky momentum conservation
3078c. input parameters:
Common /phoevt/, modcor
3079c. modcor >0
type of action
3081c. =2 corrects energy from mass
3082c. =3 corrects mass from energy
3083c. =4 corrects energy from mass for
3084c. particles up to .4 gev mass,
3085c. for heavier ones corrects mass,
3086c. =5 most complete correct also of mother
3087c. often necessary for exponentiation.
3090c. output parameters: corrected /phoevt/
3092c. author(s): p.golonka, z. was created at: 01/02/99
3093c. modified : 08/02/99
3094c.----------------------------------------------------------------------
3096 parameter(nmxpho=10000)
3098 real*8 m,p2,px,py,pz,e,en,mcut,xms
3099 INTEGER MODCOR,MODOP,I,IEV,IPRINT,K
3100 INTEGER IDPHO,ISTPHO,JDAPHO,JMOPHO,NEVPHO,NPHO
3102 common/phoevt/nevpho,npho,istpho(nmxpho),idpho(nmxpho),
3103 &jmopho(2,nmxpho),jdapho(2,nmxpho),ppho(5,nmxpho),vpho(4,nmxpho)
3114 IF (modcor.NE.0)
THEN
3118 WRITE(phlun,*)
'Message from PHCORK(MODCOR):: initialization'
3119 IF (modop.EQ.1)
THEN
3120 WRITE(phlun,*)
'MODOP=1 -- no corrections on event: DEFAULT'
3121 ELSEIF (modop.EQ.2)
THEN
3122 WRITE(phlun,*)
'MODOP=2 -- corrects Energy from mass'
3123 ELSEIF (modop.EQ.3)
THEN
3124 WRITE(phlun,*)
'MODOP=3 -- corrects mass from Energy'
3125 ELSEIF (modop.EQ.4)
THEN
3126 WRITE(phlun,*)
'MODOP=4 -- corrects Energy from mass to Mcut'
3127 WRITE(phlun,*)
' and mass from energy above Mcut '
3129 WRITE(phlun,*)
'Mcut=',mcut,
'GeV'
3130 ELSEIF (modop.EQ.5)
THEN
3131 WRITE(phlun,*)
'MODOP=5 -- corrects Energy from mass+flow'
3134 WRITE(phlun,*)
'PHCORK wrong MODCOR=',modcor
3140 IF (modop.EQ.0.AND.modcor.EQ.0)
THEN
3141 WRITE(phlun,*)
'PHCORK lack of initialization'
3155 IF (modop.EQ.1)
THEN
3156c -----------------------
3157c in this
case we
do nothing
3159 ELSEIF(modop.EQ.2)
THEN
3160c -----------------------
3161cc lets loop thru all daughters and correct their energies
3162cc according to e^2=p^2+m^2
3170 p2=ppho(1,i)**2+ppho(2,i)**2+ppho(3,i)**2
3172 en=sqrt( ppho(5,i)**2 + p2)
3175 &
WRITE(phlun,*)
"CORRECTING ENERGY OF ",i,
":",
3183 ELSEIF(modop.EQ.5)
THEN
3184c -----------------------
3185cc lets loop thru all daughters and correct their energies
3186cc according to e^2=p^2+m^2
3194 p2=ppho(1,i)**2+ppho(2,i)**2+ppho(3,i)**2
3196 en=sqrt( ppho(5,i)**2 + p2)
3199 &
WRITE(phlun,*)
"CORRECTING ENERGY OF ",i,
":",
3209 ppho(k,1)=ppho(k,1)+ppho(k,i)
3212 xms=sqrt(ppho(4,1)**2-ppho(3,1)**2-ppho(2,1)**2-ppho(1,1)**2)
3214 ELSEIF(modop.EQ.3)
THEN
3215c -----------------------
3217cc lets loop thru all daughters and correct their masses
3218cc according to e^2=p^2+m^2
3227 p2=ppho(1,i)**2+ppho(2,i)**2+ppho(3,i)**2
3229 m=sqrt(abs( ppho(4,i)**2 - p2))
3232 &
WRITE(phlun,*)
"CORRECTING MASS OF ",i,
":",
3240 ELSEIF(modop.EQ.4)
THEN
3241c -----------------------
3243cc lets loop thru all daughters and correct their masses
3244cc or energies according to e^2=p^2+m^2
3252 p2=ppho(1,i)**2+ppho(2,i)**2+ppho(3,i)**2
3254 m=sqrt(abs( ppho(4,i)**2 - p2))
3258 &
WRITE(phlun,*)
"CORRECTING MASS OF ",i,
":",
3264 en=sqrt( ppho(5,i)**2 + p2)
3267 &
WRITE(phlun,*)
"CORRECTING ENERGY OF ",i,
":",
3278 IF (iprint.EQ.1)
THEN
3279 WRITE(phlun,*)
"CORRECTING MOTHER"
3280 WRITE(phlun,*)
"PX:",ppho(1,1),
"=>",px-ppho(1,2)
3281 WRITE(phlun,*)
"PY:",ppho(2,1),
"=>",py-ppho(2,2)
3282 WRITE(phlun,*)
"PZ:",ppho(3,1),
"=>",pz-ppho(3,2)
3283 WRITE(phlun,*)
" E:",ppho(4,1),
"=>",e-ppho(4,2)
3286 ppho(1,1)=px-ppho(1,2)
3287 ppho(2,1)=py-ppho(2,2)
3288 ppho(3,1)=pz-ppho(3,2)
3289 ppho(4,1)=e -ppho(4,2)
3291 p2=ppho(1,1)**2+ppho(2,1)**2+ppho(3,1)**2
3293 IF (ppho(4,1)**2.GT.p2)
THEN
3294 m=sqrt( ppho(4,1)**2 - p2 )
3296 &
WRITE(phlun,*)
" M:",ppho(5,1),
"=>",m
3306 FUNCTION phint(IDUM)
3307c --- can be used with variant a. for b
use phint1 or 2 --------------
3308c.----------------------------------------------------------------------
3310c. phint: photos universal interference correction weight
3312c. purpose: calculates correction weight as expressed by
3313c formula(17) from cpc 79 (1994), 291.
3315c. input parameters:
Common /phoevt/, with photon added.
3317c. output parameters: correction weight
3319c. author(s): z. was, p.golonka created at: 19/01/05
3320c. last update: 25/01/05
3322c.----------------------------------------------------------------------
3327 parameter(nmxpho=10000)
3328 INTEGER IDPHO,ISTPHO,JDAPHO,JMOPHO,NEVPHO,NPHO
3330 common/phoevt/nevpho,npho,istpho(nmxpho),idpho(nmxpho),
3331 &jmopho(2,nmxpho),jdapho(2,nmxpho),ppho(5,nmxpho),vpho(4,nmxpho)
3333 DOUBLE PRECISION EMU,MCHREN,BETA,COSTHG,MPASQR,XPH, XC1, XC2,XDENO
3334 DOUBLE PRECISION XNUM1,XNUM2
3335 DOUBLE PRECISION EPS1(4),EPS2(4),PH(4),PL(4)
3339c calculate polarimetric
vector: ph, eps1, eps2 are orthogonal
3346 CALL phoeps(ph,eps2,eps1)
3347 CALL phoeps(ph,eps1,eps2)
3354 DO k=jdapho(1,1),npho-1
3356c momenta of charged particle in pl
3360c scalar products: epsilon*p/k*p
3362 xc1 = - phocha(idpho(k)) *
3363 & ( pl(1)*eps1(1) + pl(2)*eps1(2) + pl(3)*eps1(3) ) /
3364 & ( ph(4)*pl(4) - ph(1)*pl(1) - ph(2)*pl(2) - ph(3)*pl(3) )
3366 xc2 = - phocha(idpho(k)) *
3367 & ( pl(1)*eps2(1) + pl(2)*eps2(2) + pl(3)*eps2(3) ) /
3368 & ( ph(4)*pl(4) - ph(1)*pl(1) - ph(2)*pl(2) - ph(3)*pl(3) )
3371c accumulate the currents
3375 xdeno = xdeno + xc1**2 + xc2**2
3379 phint=(xnum1**2 + xnum2**2) / xdeno
3385 SUBROUTINE phoeps (VEC1, VEC2, EPS)
3386c.----------------------------------------------------------------------
3388c. phoeps: phoeps
vector product(normalized to unity)
3390c. purpose: calculates
vector product,
then normalizes its length.
3391c used to generate orthogonal vectors, i.e. to
3392c generate polarimetric vectors for photons.
3394c. input parameters: vec1,vec2 - input 4-vectors
3396c. output parameters: eps - normalized 4-
vector, orthogonal to
3399c. author(s): z. was, p.golonka created at: 19/01/05
3400c. last update: 25/01/05
3402c.----------------------------------------------------------------------
3404 DOUBLE PRECISION VEC1(4), VEC2(4), EPS(4),XN
3406 eps(1)=vec1(2)*vec2(3) - vec1(3)*vec2(2)
3407 eps(2)=vec1(3)*vec2(1) - vec1(1)*vec2(3)
3408 eps(3)=vec1(1)*vec2(2) - vec1(2)*vec2(1)
3411 xn=sqrt( eps(1)**2 +eps(2)**2 +eps(3)**2)
3420c.----------------------------------------------------------------------
3422c. photos: photon radiation in decays event dump routine
3424c. purpose: print event record.
3426c. input parameters:
Common /hepevt/
3428c. output parameters: none
3430c. author(s): b. van eijk created at: 05/06/90
3431c. last update: 05/06/90
3433c.----------------------------------------------------------------------
3435 DOUBLE PRECISION SUMVEC(5)
3437c this is the hepevt
class in old style. No d_h_
class pre-name
3439 parameter(nmxhep=10000)
3441 INTEGER nevhep,nhep,isthep,idhep,jmohep,
3452* ----------------------------------------------------------------------
3456* ----------------------------------------------------------------------
3463c-- print event number...
3465 WRITE(phlun,9010) nevhep
3470c-- for
'stable particle' calculate
vector momentum sum
3471 IF (jdahep(1,i).EQ.0)
THEN
3473 20 sumvec(j)=sumvec(j)+phep(j,i)
3474 IF (jmohep(2,i).EQ.0)
THEN
3475 WRITE(phlun,9030) i,idhep(i),jmohep(1,i),(phep(j,i),j=1,5)
3477 WRITE(phlun,9040) i,idhep(i),jmohep(1,i),jmohep(2,i),(phep
3481 IF (jmohep(2,i).EQ.0)
THEN
3482 WRITE(phlun,9050) i,idhep(i),jmohep(1,i),jdahep(1,i),
3483 & jdahep(2,i),(phep(j,i),j=1,5)
3485 WRITE(phlun,9060) i,idhep(i),jmohep(1,i),jmohep(2,i),
3486 & jdahep(1,i),jdahep(2,i),(phep(j,i),j=1,5)
3490 sumvec(5)=sqrt(sumvec(4)**2-sumvec(1)**2-sumvec(2)**2-
3492 WRITE(phlun,9070) (sumvec(j),j=1,5)
3494 9000
FORMAT(1h0,80(
'='))
3495 9010
FORMAT(1h ,29x,
'Event No.:',i10)
3496 9020
FORMAT(1h0,1x,
'Nr',3x,
'Type',3x,
'Parent(s)',2x,
'Daughter(s)',6x,
3497 &
'Px',7x,
'Py',7x,
'Pz',7x,
'E',4x,
'Inv. M.')
3498 9030
FORMAT(1h ,i4,i7,3x,i4,9x,
'Stable',2x,5f9.2)
3499 9040
FORMAT(1h ,i4,i7,i4,
' - ',i4,5x,
'Stable',2x,5f9.2)
3500 9050
FORMAT(1h ,i4,i7,3x,i4,6x,i4,
' - ',i4,5f9.2)
3501 9060
FORMAT(1h ,i4,i7,i4,
' - ',i4,2x,i4,
' - ',i4,5f9.2)
3502 9070
FORMAT(1h0,23x,
'Vector Sum: ', 5f9.2)
3503 9080
FORMAT(1h0,6x,
'Particle Parameters')