1 /* copyright(c) 1991-2021 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
48 *// parts of the code.
50 *// One should add suffix c_Photos_ to names of all commons as soon as
52 *///////////////////////////////////////////////////////////////////////
54 C.----------------------------------------------------------------------
56 C. PHOTOS: PHOtos CDE-s
58 C. Purpose: Keep definitions for PHOTOS QED correction Monte Carlo.
60 C. Input Parameters: None
62 C. Output Parameters: None
64 C. Author(s): Z. Was, B. van Eijk Created at: 29/11/89
65 C. Last Update: 10/08/93
67 C. =========================================================
68 C. General Structure Information: =
69 C. =========================================================
76 C. 2) GENERAL INTERFACE:
83 C. PHOTWO (specific interface
86 C. PHTYPE (specific interface
87 C. PHOMAK (specific interface
88 C. 3) QED PHOTON GENERATION:
115 C. NAME USED IN SECT. # OF OCC. Comment
116 C. PHOQED 1) 2) 3 Flags whether emisson to be gen.
117 C. PHOLUN 1) 4) 6 Output device number
118 C. PHOCOP 1) 3) 4 photon coupling & min energy
119 C. PHPICO 1) 3) 4) 5 PI & 2*PI
120 C. PHSEED 1) 4) 3 RN seed
121 C. PHOSTA 1) 4) 3 Status information
122 C. PHOKEY 1) 2) 3) 7 Keys for nonstandard application
123 C. PHOVER 1) 1 Version info for outside
124 C. HEPEVT 2) 2 PDG common
125 C. PH_HEPEVT2) 8 PDG common internal
126 C. PHOEVT 2) 3) 10 PDG branch
127 C. PHOIF 2) 3) 2 emission flags for PDG branch
128 C. PHOMOM 3) 5 param of char-neutr system
129 C. PHOPHS 3) 5 photon momentum parameters
130 C. PHOPRO 3) 4 var. for photon rep. (in branch)
131 C. PHOCMS 2) 3 parameters of boost to branch CMS
132 C. PHNUM 4) 1 event number from outside
133 C.----------------------------------------------------------------------
135 C.----------------------------------------------------------------------
137 C. PHOTOS: PHOton radiation in decays INItialisation
139 C. Purpose: Initialisation routine for the PHOTOS QED radiation
140 C. package. Should be called at least once before a call
141 C. to the steering program 'photos
' is made.
143 C. Input Parameters: None
145 C. Output Parameters: None
147 C. Author(s): Z. Was, B. van Eijk Created at: 26/11/89
148 C. Last Update: 12/04/90
150 C.----------------------------------------------------------------------
152 INTEGER INIT,IDUM,IPHQRK,IPHEKL
156 C-- Return if already initialized...
157 .NE.
IF (INIT0) RETURN
160 C-- all the following parameter setters can be called after PHOINI.
161 C-- Initialization of kinematic correction against rounding errors.
162 C-- The set values will be used later if called wit zero.
163 C-- Default parameter is 1 (no correction) optionally 2, 3, 4
164 C-- In case of exponentiation new version 5 is needed in most cases.
165 C-- Definition given here will be thus overwritten in such a case
166 C-- below in routine PHOCIN
168 C-- blocks emission from quarks if parameter is 1 (enables if 2),
169 C-- physical treatment
170 C-- will be 3, option 2 is not realistic and for tests only,
171 IDUM= IPHQRK(1) ! default is 1
172 C-- blocks emission in pi0 to gamma e+ e- if parameter is gt.1
173 C-- (enables otherwise)
174 IDUM= IPHEKL(2) ! default is 1
176 C-- Preset parameters in PHOTOS commons
183 C-- Initialize Marsaglia and Zaman random number generator
188 C.----------------------------------------------------------------------
190 C. PHOTOS: PHOton Common INitialisation
192 C. Purpose: Initialisation of parameters in common blocks.
194 C. Input Parameters: None
196 C. Output Parameters: Commons /PHOLUN/, /PHOPHO/, /PHOCOP/, /PHPICO/
199 C. Author(s): B. van Eijk Created at: 26/11/89
200 C. Z. Was Last Update: 29/01/05
202 C.----------------------------------------------------------------------
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
228 C-- Return if already initialized...
229 .NE.
IF (INIT0) RETURN
232 C-- Preset switch for photon emission to 'true
' for each particle in
233 C-- /PH_HEPEVT/, this interface is needed for KORALB and KORALZ...
237 C-- Logical output unit for printing of PHOTOS error messages
240 C-- Set cut parameter for photon radiation
241 XPHCUT=0.01 D0 ! 0.0001D0! to go to low valuex (IEXP excepted)
242 C-- ! switch to - VARIANT B
244 C-- Define some constants
245 ALPHA=0.00729735039D0
246 PI=3.14159265358979324D0
247 TWOPI=6.28318530717958648D0
249 C-- Default seeds Marsaglia and Zaman random number generator
253 C-- Iitialization for extra options
255 C-- Interference weight now universal.
258 C-- Second order - double photon switch
260 C-- Third/fourth order - triple (or quatric) photon switch,
261 C-- see dipswitch ifour
263 C-- Exponentiation on:
268 CALL PHCORK(5) ! in case of exponentiation correction of ph space
269 ! is a default mandatory
274 C-- Emision in the hard process g g (q qbar) --> t tbar
278 C-- further initialization done automatically
279 C-- see places with - VARIANT A - VARIANT B - all over
280 C-- to switch between options.
281 C ----------- SLOWER VARIANT A, but stable ------------------
282 C --- it is limiting choice for small XPHCUT in fixed orer
283 C --- modes of operation
285 C-- best choice is if FINT=2**N where N+1 is maximal number
286 C-- of charged daughters
287 C-- see report on overweihted events
292 C ----------- FASTER VARIANT B ------------------
293 C -- it is good for tests of fixed order and small XPHCUT
294 C -- but is less promising for more complex cases of interference
295 C -- sometimes fails because of that
302 C----------END VARIANTS A B -----------------------
304 C-- Effects of initial state charge (in leptonic W decays)
307 C-- Initialise status counter for warning messages
313 C.----------------------------------------------------------------------
315 C. PHOTOS: PHOton radiation in decays general INFo
317 C. Purpose: Print PHOTOS info
319 C. Input Parameters: PHOLUN
321 C. Output Parameters: PHOVN1, PHOVN2
323 C. Author(s): B. van Eijk Created at: 12/04/90
324 C. Last Update: 27/06/04
326 C.----------------------------------------------------------------------
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
339 C-- 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)
423 C.----------------------------------------------------------------------
425 C. PHOTOS: General search routine + _GET + _SET
427 C. Purpose: /HEPEVT/ is not anymore a standard at least
428 C. REAL*8 REAL*4 are in use. PHOTOS_GET and PHOTOS_SET
429 C. were to be introduced.
432 C. Input Parameters: ID see routine PHOTOS_MAKE
434 C. Output Parameters: None
436 C. Author(s): Z. Was Created at: 21/07/98
437 C. Last Update: 21/07/98
439 C.----------------------------------------------------------------------
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
464 C.----------------------------------------------------------------------
466 C. Getter for PHOTOS:
468 C. Purpose: Copies /HEPEVT/ into /PH_HEPEVT/
471 C. Input Parameters: None
473 C. Output Parameters: None
475 C. Author(s): Z. Was Created at: 21/07/98
476 C. Last Update: 21/07/98
478 C.----------------------------------------------------------------------
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
528 C.----------------------------------------------------------------------
530 C. Setter for PHOTOS:
532 C. Purpose: Copies /PH_HEPEVT/ into /HEPEVT/
535 C. Input Parameters: None
537 C. Output Parameters: None
539 C. Author(s): Z. Was Created at: 21/07/98
540 C. Last Update: 21/07/98
542 C.----------------------------------------------------------------------
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)
590 C.----------------------------------------------------------------------
592 C. PHOTOS_MAKE: General search routine
594 C. Purpose: Search through the /PH_HEPEVT/ standard HEP common, sta-
595 C. rting from the IPPAR-th particle. Whenevr branching
596 C. point is found routine PHTYPE(IP) is called.
597 C. Finally if calls on PHTYPE(IP) modified entries, common
598 C /PH_HEPEVT/ is ordered.
600 C. Input Parameter: IPPAR: Pointer to decaying particle in
601 C. /PH_HEPEVT/ and the common itself,
603 C. Output Parameters: Common /PH_HEPEVT/, either with or without
604 C. new particles added.
606 C. Author(s): Z. Was, B. van Eijk Created at: 26/11/89
607 C. Last Update: 30/08/93
609 C.----------------------------------------------------------------------
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)
631 C-- Store pointers for cascade treatement...
636 C-- Check decay multiplicity and minimum of correctness..
637 .EQ..OR..NE.
IF ((JDAHEP(1,IP)0)(JMOHEP(1,JDAHEP(1,IP))IP)) RETURN
639 C-- single branch mode
640 C-- we start looking for the decay points in the cascade
641 C-- IPPAR is original position where the program was called
643 C-- NUMIT denotes number of secondary decay branches
645 C-- NTRY denotes number of secondary branches already checked for
646 C-- for existence of further branches
648 C-- let-s search if IPARR does not prevent searching.
649 .GT.
IF (IPARR0) THEN
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
667 C-- 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))
680 C-- 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)
692 C-- rearrange /PH_HEPEVT/ to get correct order..
693 .GT.
IF (NHEPNLAST) THEN
694 DO 160 I=NLAST+1,NHEP
696 C-- Photon mother and position...
698 POSPHO=JDAHEP(2,MOTHER)+1
699 C-- Intermediate save of photon energy/momentum and pointers
701 90 PHOTON(J)=PHEP(J,I)
708 C-- Exclude photon in sequence !
709 .NE.
IF (POSPHONHEP) THEN
712 C-- Order /PH_HEPEVT/
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)
725 C-- 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
734 C-- Store photon energy/momentum
736 140 PHEP(J,POSPHO)=PHOTON(J)
739 C-- 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
748 C-- Get photon production vertex position
750 150 VHEP(J,POSPHO)=VHEP(J,POSPHO-1)
755 SUBROUTINE PHOBOS(IP,PBOOS1,PBOOS2,FIRST,LAST)
756 C.----------------------------------------------------------------------
758 C. PHOBOS: PHOton radiation in decays BOoSt routine
760 C. Purpose: Boost particles in cascade decay to parent rest frame
761 C. and boost back with modified boost vector.
763 C. Input Parameters: IP: pointer of particle starting chain
765 C. PBOOS1: Boost vector to rest frame,
766 C. PBOOS2: Boost vector to modified frame,
767 C. FIRST: Pointer to first particle to be boos-
768 C. ted (/PH_HEPEVT/),
769 C. LAST: Pointer to last particle to be boos-
770 C. ted (/PH_HEPEVT/).
772 C. Output Parameters: Common /PH_HEPEVT/.
774 C. Author(s): B. van Eijk Created at: 13/02/90
775 C. Z. Was Last Update: 16/11/93
777 C.----------------------------------------------------------------------
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)
798 C-- 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
806 C-- ...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
814 C-- Check on stack length...
815 .GT.
IF (NSTACKMAXSTA) THEN
817 CALL PHOERR(7,'phobos
',DATA)
823 .NE.
IF (NSTACK0) THEN
825 C-- 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)
835 C.----------------------------------------------------------------------
837 C. PHOIN: PHOtos INput
839 C. Purpose: copies IP branch of the common /PH_HEPEVT/ into /PHOEVT/
840 C. moves branch into its CMS system.
842 C. Input Parameters: IP: pointer of particle starting branch
844 C. BOOST: Flag whether boost to CMS was or was
847 C. Output Parameters: Commons: /PHOEVT/, /PHOCMS/
849 C. Author(s): Z. Was Created at: 24/05/93
850 C. Last Update: 16/11/93
852 C.----------------------------------------------------------------------
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
875 C let-s calculate size of the little common entry
878 NPHO=3+LAST-FIRST+NHEP-NHEP0
880 C let-s take in decaying particle
883 JDAPHO(2,1)=3+LAST-FIRST
887 C 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)
902 C 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
912 C 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)
922 C-- 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)
929 C special case of t tbar production process
930 IF(IFTOP) CALL PHOTWO(0)
932 C-- Check whether parent is in its rest frame...
934 C bug reported by Vladimir Savinov localized and fixed.
935 C protection against rounding error was back-firing if soft
936 C momentum of mother was physical. Consequence was that PHCORK was
937 C messing up masses of final state particles in vertex of the decay.
938 C Only configurations with previously generated photons of energy fraction
939 C 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
949 C-- Boost daughter particles to rest frame of parent...
950 C-- 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
959 C-- 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
967 C 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)
974 C.----------------------------------------------------------------------
976 C. PHOTWO: PHOtos but TWO mothers allowed
978 C. Purpose: Combines two mothers into one in /PHOEVT/
979 C. necessary eg in case of g g (q qbar) --> t tbar
981 C. Input Parameters: Common /PHOEVT/ (/PHOCMS/)
983 C. Output Parameters: Common /PHOEVT/, (stored mothers)
985 C. Author(s): Z. Was Created at: 5/08/93
986 C. Last Update:10/08/93
988 C.----------------------------------------------------------------------
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
1001 C logical IFRAD is used to tag cases when two mothers may be
1002 C merged to the sole one.
1003 C So far used in case:
1004 C 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)
1016 c.....combining first and second mother
1018 PPHO(I,1)=PPHO(I,1)+PPHO(I,2)
1020 PPHO(5,1)=SQRT(MPASQR)
1021 c.....removing second mother,
1027 C boosting of the mothers to the reaction frame not implemented yet.
1028 C to do it in mode 0 original mothers have to be stored in new comon (?)
1029 C and in mode 1 boosted to cms.
1032 SUBROUTINE PHOOUT(IP,BOOST,NHEP0)
1033 C.----------------------------------------------------------------------
1035 C. PHOOUT: PHOtos OUTput
1037 C. Purpose: copies back IP branch of the common /PH_HEPEVT/ from
1038 C. /PHOEVT/ moves branch back from its CMS system.
1040 C. Input Parameters: IP: pointer of particle starting branch
1041 C. to be given back.
1042 C. BOOST: Flag whether boost to CMS was or was
1045 C. Output Parameters: Common /PHOEVT/,
1047 C. Author(s): Z. Was Created at: 24/05/93
1050 C.----------------------------------------------------------------------
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
1070 C-- 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
1078 C-- ...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
1088 C let-s take in original daughters
1090 IDHEP(FIRST+LL) = IDPHO(3+LL)
1092 PHEP(I,FIRST+LL) = PPHO(I,3+LL)
1095 C 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)
1112 C.----------------------------------------------------------------------
1114 C. PHOCHK: checking branch.
1116 C. Purpose: checks whether particles in the common block /PHOEVT/
1117 C. can be served by PHOMAK.
1118 C. JFIRST is the position in /PH_HEPEVT/ (!) of the first
1119 C. daughter of sub-branch under action.
1122 C. Author(s): Z. Was Created at: 22/10/92
1123 C. Last Update: 11/12/00
1125 C.----------------------------------------------------------------------
1126 C ********************
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
1148 C 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
1162 C checking for good particles
1164 .GT.
IF (IEKL1) THEN ! exclude radiative corr in decay of pi0
1165 C ! 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))
1177 C 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)
1184 C now we go to special cases, where CHKIF(I) will be overwritten
1187 C 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)))
1199 .AND..EQ.
& (IDENT4)
1203 .GT..AND.
IF(I2) CHKIF(I)=CHKIF(I)QEDRAD(JFIRST+I-IPPAR-2)
1210 C 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))
1222 .AND..EQ.
& (IDENT4)
1226 .GT..AND.
IF(I2) CHKIF(I)=CHKIF(I)QEDRAD(JFIRST+I-IPPAR-2)
1233 SUBROUTINE PHTYPE(ID)
1234 C.----------------------------------------------------------------------
1236 C. PHTYPE: Central manadgement routine.
1238 C. Purpose: defines what kind of the
1239 C. actions will be performed at point ID.
1241 C. Input Parameters: ID: pointer of particle starting branch
1242 C. in /PH_HEPEVT/ to be treated.
1244 C. Output Parameters: Common /PH_HEPEVT/.
1246 C. Author(s): Z. Was Created at: 24/05/93
1247 C. P. Golonka Last Update: 27/06/04
1249 C.----------------------------------------------------------------------
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.
1275 C-- 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
1312 C-- 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)
1327 C-- 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)
1338 C-- double photon emission
1341 .GE.
IF (RN0.5D0) THEN
1342 CALL PHOMAK(ID,NHEP0)
1343 CALL PHOMAK(ID,NHEP0)
1346 C-- single photon emission
1348 CALL PHOMAK(ID,NHEP0)
1351 C-- electron positron pair (coomented out for a while
1352 C IF (IPAIR) CALL PHOPAR(ID,NHEP0)
1354 SUBROUTINE PHOMAK(IPPAR,NHEP0)
1355 C.----------------------------------------------------------------------
1357 C. PHOMAK: PHOtos MAKe
1359 C. Purpose: Single or double bremstrahlung radiative corrections
1360 C. are generated in the decay of the IPPAR-th particle in
1361 C. the HEP common /PH_HEPEVT/. Example of the use of
1364 C. Input Parameter: IPPAR: Pointer to decaying particle in
1365 C. /PH_HEPEVT/ and the common itself
1367 C. Output Parameters: Common /PH_HEPEVT/, either with or without
1370 C. Author(s): Z. Was, Created at: 26/05/93
1371 C. Last Update: 29/01/05
1373 C.----------------------------------------------------------------------
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
1403 C PHODO is caling PHORAN, thus change of series if it is moved before if
1404 CALL PHODO(1,NCHARB,NEUDAU)
1405 C 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)
1417 C.----------------------------------------------------------------------
1419 C. PHINT: PHotos INTerference (Old version kept for tests only.
1421 C. Purpose: Calculates interference between emission of photons from
1422 C. different possible chaged daughters stored in
1423 C. the HEP common /PHOEVT/.
1425 C. Input Parameter: commons /PHOEVT/ /PHOMOM/ /PHOPHS/
1428 C. Output Parameters:
1431 C. Author(s): Z. Was, Created at: 10/08/93
1432 C. Last Update: 15/03/99
1434 C.----------------------------------------------------------------------
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
1463 C check if there is a photon
1464 .GT.
IFINT= NPHOIDENT
1465 C check if it is two body + gammas reaction
1466 .AND..EQ.
IFINT= IFINT(IDENT-JDAPHO(1,1))1
1467 C check if two body was particle antiparticle
1468 .AND..EQ.
IFINT= IFINTIDPHO(JDAPHO(1,1))-IDPHO(IDENT)
1469 C check if particles were charged
1470 .AND..NE.
IFINT= IFINTPHOCHA(IDPHO(IDENT))0
1471 C 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)
1485 C.----------------------------------------------------------------------
1487 C. PHINT: PHotos INTerference
1489 C. Purpose: Calculates interference between emission of photons from
1490 C. different possible chaged daughters stored in
1491 C. the HEP common /PHOEVT/.
1493 C. Input Parameter: commons /PHOEVT/ /PHOMOM/ /PHOPHS/
1496 C. Output Parameters:
1499 C. Author(s): Z. Was, Created at: 10/08/93
1502 C.----------------------------------------------------------------------
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
1532 C check if there is a photon
1533 .GT.
IFINT= NPHOIDENT
1534 C check if it is two body + gammas reaction
1535 .AND..EQ.
IFINT= IFINT(IDENT-JDAPHO(1,1))1
1536 C check if two body was particle antiparticle (we improve on it !
1537 .AND..EQ.
C IFINT= IFINTIDPHO(JDAPHO(1,1))-IDPHO(IDENT)
1538 C check if particles were charged
1539 .AND..GT.
IFINT= IFINTabs(PHOCHA(IDPHO(IDENT)))0.01D0
1540 C check if they have both charge
1542 .gt.
$ abs(PHOCHA(IDPHO(JDAPHO(1,1))))0.01D0
1543 C 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)
1590 C.----------------------------------------------------------------------
1592 C. PHOTOS: Photon radiation in decays
1594 C. Purpose: Order (alpha) radiative corrections are generated in
1595 C. the decay of the IPPAR-th particle in the HEP-like
1596 C. common /PHOEVT/. Photon radiation takes place from one
1597 C. of the charged daughters of the decaying particle IPPAR
1598 C. WT is calculated, eventual rejection will be performed
1599 C. later after inclusion of interference weight.
1601 C. Input Parameter: IPPAR: Pointer to decaying particle in
1602 C. /PHOEVT/ and the common itself,
1604 C. Output Parameters: Common /PHOEVT/, either with or without a
1606 C. WT weight of the configuration
1608 C. Author(s): Z. Was, B. van Eijk Created at: 26/11/89
1609 C. Last Update: 29/01/05
1611 C.----------------------------------------------------------------------
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
1640 C 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
1647 C-- Store pointers for cascade treatement...
1652 C-- Check decay multiplicity..
1653 .EQ.
IF (JDAPHO(1,IP)0) RETURN
1655 C-- Loop over daughters, determine charge multiplicity
1660 DO 20 I=JDAPHO(1,IP),JDAPHO(2,IP)
1663 C-- 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
1680 C-- 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
1683 C-- Order charged particles according to decreasing mass, this to
1684 C-- 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))
1692 C-- 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)
1706 C-- Determine kinematical limit...
1707 XPHMAX=(MPASQR-(PNEUTR(5)+PPHO(5,CHAPOI(NCHARG)))**2)/MPASQR
1709 C-- 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 !!
1715 C-- Energy fraction not too large (very seldom) ? Define angle.
1716 .LT..OR..GT.
ELSEIF ((XPHOTOXPHCUT)(XPHOTOXPHMAX)) THEN
1718 C-- No radiation was accepted, check for more daughters that may ra-
1719 C-- diate and correct radiation probability...
1721 .GT.
IF (NCHARG0) THEN
1727 C-- Angle is generated in the frame defined by charged vector and
1728 C-- PNEUTR, distribution is taken in the infrared limit...
1729 EPS=MCHREN/(1.D0+BETA)
1731 C-- Calculate sin(theta) and cos(theta) from interval variables
1732 DEL1=(2.D0-EPS)*(EPS/(2.D0-EPS))**PHORAN(THEDUM)
1735 C ----------- VARIANT B ------------------
1736 CC corrections for more efiicient interference correction,
1737 CC instead of doubling crude distribution, we add flat parallel channel
1738 .LT.
C IF (PHORAN(THEDUM)BIGLOG/BETA/(BIGLOG/BETA+2*FINT)) THEN
1739 C COSTHG=(1.D0-DEL1)/BETA
1740 C SINTHG=SQRT(DEL1*DEL2-MCHREN)/BETA
1742 C COSTHG=-1D0+2*PHORAN(THEDUM)
1743 C SINTHG= SQRT(1D0-COSTHG**2)
1746 .GT.
C IF (FINT1.0D0) THEN
1748 C WGT=1D0/(1D0-BETA*COSTHG)
1749 C WGT=WGT/(WGT+FINT)
1756 C ----------- END OF VARIANT B ------------------
1758 C ----------- VARIANT A ------------------
1759 COSTHG=(1.D0-DEL1)/BETA
1760 SINTHG=SQRT(DEL1*DEL2-MCHREN)/BETA
1762 C ----------- END OF VARIANT A ------------------
1765 C-- Determine spin of particle and construct code for matrix element
1766 ME=2.D0*PHOSPI(IDPHO(CHAPOI(NCHARG)))+1.D0
1768 C-- Weighting procedure with 'exact
' matrix element, reconstruct kine-
1769 C-- matics for photon, neutral and charged system and update /PHOEVT/.
1770 C-- Find pointer to the first component of 'neutral
' system
1771 DO I=JDAPHO(1,IP),JDAPHO(2,IP)
1772 .NE.
IF (ICHAPOI(NCHARG)) THEN
1778 C-- 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)
1796 C.----------------------------------------------------------------------
1798 C. PHOTOS: PHOton radiation in decays Order MAss vector
1800 C. Purpose: Order the contents of array 'pointr
' according to the
1801 C. decreasing value in the array 'mass
'.
1803 C. Input Parameters: IFIRST, ILAST: Pointers to the vector loca-
1805 C. POINTR: Unsorted array with pointers to
1808 C. Output Parameter: POINTR: Sorted arrays with respect to
1809 C. particle mass 'ppho(5,*)
'.
1811 C. Author(s): B. van Eijk Created at: 28/11/89
1812 C. Last Update: 27/05/93
1814 C.----------------------------------------------------------------------
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
1826 C-- Copy particle masses
1827 DO 10 I=IFIRST,ILAST
1828 10 MASS(I)=PPHO(5,POINTR(I))
1830 C-- 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)
1845 C.----------------------------------------------------------------------
1847 C. PHOTOS: PHOton radiation in decays calculation of photon ENErgy
1850 C. Purpose: Subroutine returns photon energy fraction (in (parent
1851 C. mass)/2 units) for the decay bremsstrahlung.
1853 C. Input Parameters: MPASQR: Mass of decaying system squared,
1854 C. XPHCUT: Minimum energy fraction of photon,
1855 C. XPHMAX: Maximum energy fraction of photon.
1857 C. Output Parameter: MCHREN: Renormalised mass squared,
1858 C. BETA: Beta factor due to renormalisation,
1859 C. XPHOTO: Photon energy fraction,
1860 C. XF: Correction factor for PHOFAC.
1862 C. Author(s): S. Jadach, Z. Was Created at: 01/01/89
1863 C. B. van Eijk, P.Golonka Last Update: 29/01/05
1865 C.----------------------------------------------------------------------
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
1899 C-- Probabilities for hard and soft bremstrahlung...
1900 MCHREN=4.D0*MCHSQR/MPASQR/(1.D0+MCHSQR/MPASQR)**2
1901 BETA=SQRT(1.D0-MCHREN)
1903 C ----------- VARIANT B ------------------
1904 CC we replace 1D0/BETA*BIGLOG with (1D0/BETA*BIGLOG+2*FINT)
1905 CC for integral of new crude
1906 C BIGLOG=LOG(MPASQR/MCHSQR*(1.D0+BETA)**2/4.D0*
1907 C & (1.D0+MCHSQR/MPASQR)**2)
1908 C PRHARD=ALPHA/PI*(1D0/BETA*BIGLOG+2*FINT)*(LOG(XPHMAX/XPHCUT)
1909 C &-.75D0+XPHCUT/XPHMAX-.25D0*XPHCUT**2/XPHMAX**2)
1910 C PRHARD=PRHARD*PHOCHA(IDENT)**2*FSEC
1911 C ----------- END OF VARIANT B ------------------
1913 C ----------- 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
1919 C ----------- 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.
1950 C-- 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
1966 C-- No photon... (ie. photon too soft)
1968 .LT.
IF (RRRPRKILL) XPHOTO=-5d0 ! No photon...no further trials
1971 C-- Hard photon... (ie. photon hard enough).
1972 C-- 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))
1979 C-- Calculate parameter for PHOFAC function
1980 XF=4.D0*MCHSQR*MPASQR/(MPASQR+MCHSQR-MNESQR)**2
1983 FUNCTION PHOCOR(MPASQR,MCHREN,ME)
1984 C.----------------------------------------------------------------------
1986 C. PHOTOS: PHOton radiation in decays CORrection weight from
1989 C. Purpose: Calculate photon angle. The reshaping functions will
1990 C. have to depend on the spin S of the charged particle.
1991 C. We define: ME = 2 * S + 1 !
1993 C. Input Parameters: MPASQR: Parent mass squared,
1994 C. MCHREN: Renormalised mass of charged system,
1995 C. ME: 2 * spin + 1 determines matrix element
1997 C. Output Parameter: Function value.
1999 C. Author(s): Z. Was, B. van Eijk Created at: 26/11/89
2000 C. Last Update: 21/03/93
2002 C.----------------------------------------------------------------------
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
2017 C-- 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)
2049 C.----------------------------------------------------------------------
2051 C. PHOTOS: PHOton radiation in decays control FACtor
2053 C. Purpose: This is the control function for the photon spectrum and
2054 C. final weighting. It is called from PHOENE for genera-
2055 C. ting the raw photon energy spectrum (MODE=0) and in PHO-
2056 C. COR to scale the final weight (MODE=1). The factor con-
2057 C. sists of 3 terms. Addition of the factor FF which mul-
2058 C. tiplies PHOFAC for MODE=0 and divides PHOFAC for MODE=1,
2059 C. does not affect the results for the MC generation. An
2060 C. appropriate choice for FF can speed up the calculation.
2061 C. Note that a too small value of FF may cause weight over-
2062 C. flow in PHOCOR and will generate a warning, halting the
2063 C. execution. PRX should be included for repeated calls
2064 C. for the same event, allowing more particles to radiate
2065 C. photons. At the first call IREP=0, for more than 1
2066 C. charged decay products, IREP >= 1. Thus, PRSOFT (no
2067 C. photon radiation probability in the previous calls)
2068 C. appropriately scales the strength of the bremsstrahlung.
2070 C. Input Parameters: MODE, PROBH, XF
2072 C. Output Parameter: Function value
2074 C. Author(s): S. Jadach, Z. Was Created at: 01/01/89
2075 C. B. van Eijk, P.Golonka Last Update: 26/06/04
2077 C.----------------------------------------------------------------------
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)
2102 C-- Following options are not considered for the time being...
2103 C-- (1) Good choice, but does not save very much time:
2104 C-- FF=(1.0D0-SQRT(XF)/2.0D0)/(1.0+SQRT(XF)/2.0D0)
2105 C-- (2) Taken from the blue, but works without weight overflows...
2106 C-- FF=(1.D0-XF/(1-(1-SQRT(XF))**2))*(1+(1-SQRT(XF))/SQRT(1-XF))/2
2112 SUBROUTINE PHOBW(WT)
2113 C.----------------------------------------------------------------------
2115 C. PHOTOS: PHOtos Boson W correction weight
2117 C. Purpose: calculates correction weight due to amplitudes of
2118 C. emission from W boson.
2124 C. Input Parameters: Common /PHOEVT/, with photon added.
2125 C. wt to be corrected
2129 C. Output Parameters: wt
2131 C. Author(s): G. Nanava, Z. Was Created at: 13/03/03
2132 C. Last Update: 13/03/03
2134 C.----------------------------------------------------------------------
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))
2174 c write(*,*) IDPHO(1),IDPHO(JDAPHO(1,1)),IDPHO(JDAPHO(1,1)+1)
2175 c write(*,*) emu,xph,costhg,beta,mpasqr,mchren
2178 SUBROUTINE PHODO(IP,NCHARB,NEUDAU)
2179 C.----------------------------------------------------------------------
2181 C. PHOTOS: PHOton radiation in decays DOing of KINematics
2183 C. Purpose: Starting from the charged particle energy/momentum,
2184 C. PNEUTR, photon energy fraction and photon angle with
2185 C. respect to the axis formed by charged particle energy/
2186 C. momentum vector and PNEUTR, scale the energy/momentum,
2187 C. keeping the original direction of the neutral system in
2188 C. the lab. frame untouched.
2190 C. Input Parameters: IP: Pointer to decaying particle in
2191 C. /PHOEVT/ and the common itself
2192 C. NCHARB: pointer to the charged radiating
2193 C. daughter in /PHOEVT/.
2194 C. NEUDAU: pointer to the first neutral daughter
2195 C. Output Parameters: Common /PHOEVT/, with photon added.
2197 C. Author(s): Z. Was, B. van Eijk Created at: 26/11/89
2198 C. Last Update: 27/05/93
2200 C.----------------------------------------------------------------------
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))
2226 C-- Reconstruct kinematics of charged particle and neutral system
2227 FI1=PHOAN1(PNEUTR(1),PNEUTR(2))
2229 C-- Choose axis along z of PNEUTR, calculate angle between x and y
2230 C-- 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))
2235 C-- Take away photon energy from charged particle and PNEUTR ! Thus
2236 C-- the onshell charged particle decays into virtual charged particle
2237 C-- and photon. The virtual charged particle mass becomes:
2238 C-- SQRT(PPHO(5,IP)*(PPHO(5,IP)-2*EPHOTO)). Construct new PNEUTR mo-
2239 C-- mentum in the rest frame of the parent:
2240 C-- 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))
2251 C-- 2) ...reductive boost...
2252 CALL PHOBO3(PARNE,PNEUTR)
2254 C-- ...calculate photon energy in the reduced system...
2258 C-- Photon mother and daughter pointers !
2263 PPHO(4,NPHO)=EPHOTO*PPHO(5,IP)/PMAVIR
2265 C-- ...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)
2273 C-- Minus sign because axis opposite direction of charged particle !
2274 PPHO(3,NPHO)=-PPHO(4,NPHO)*COSTHG
2277 C-- 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)
2284 C-- Boost to the rest frame of decaying particle
2285 CALL PHOBO3(ANGLE,PNEUTR(1))
2286 CALL PHOBO3(ANGLE,PPHO(1,NPHO))
2288 C-- 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))
2306 C-- 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))
2313 C-- See whether neutral system has multiplicity larger than 1...
2314 .GT.
IF ((JDAPHO(2,IP)-JDAPHO(1,IP))1) THEN
2315 C-- Find pointers to components of 'neutral
' system
2320 .NE..AND..EQ.
IF (INCHARB(JMOPHO(1,I)IP)) THEN
2322 C-- Reconstruct kinematics...
2323 CALL PHORO3(-FI1,PPHO(1,I))
2324 CALL PHORO2(-TH1,PPHO(1,I))
2326 C-- ...reductive boost
2327 CALL PHOBO3(PARNE,PPHO(1,I))
2329 C-- Rotate in order to get photon along z-axis
2330 CALL PHORO3(-FI3,PPHO(1,I))
2331 CALL PHORO2(-TH3,PPHO(1,I))
2333 C-- Boost to the rest frame of decaying particle
2334 CALL PHOBO3(ANGLE,PPHO(1,I))
2336 C-- 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))
2340 C-- 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))
2348 C-- ...only one 'neutral
' particle in addition to photon!
2350 80 PPHO(J,NEUDAU)=PNEUTR(J)
2353 C-- 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)
2360 C.----------------------------------------------------------------------
2362 C. PHOTOS: PHOton radiation in decays calculation of TRIangle fie
2364 C. Purpose: Calculation of triangle function for phase space.
2366 C. Input Parameters: A, B, C (Virtual) particle masses.
2368 C. Output Parameter: Function value =
2369 C. SQRT(LAMBDA(A**2,B**2,C**2))/(2*A)
2371 C. Author(s): B. van Eijk Created at: 15/11/89
2372 C. Last Update: 02/01/90
2374 C.----------------------------------------------------------------------
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)
2388 C.----------------------------------------------------------------------
2390 C. PHOTOS: PHOton radiation in decays calculation of ANgle '1
'
2392 C. Purpose: Calculate angle from X and Y
2394 C. Input Parameters: X, Y
2396 C. Output Parameter: Function value
2398 C. Author(s): S. Jadach Created at: 01/01/89
2399 C. B. van Eijk Last Update: 02/01/90
2401 C.----------------------------------------------------------------------
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)
2417 C.----------------------------------------------------------------------
2419 C. PHOTOS: PHOton radiation in decays calculation of ANgle '2
'
2421 C. Purpose: Calculate angle from X and Y
2423 C. Input Parameters: X, Y
2425 C. Output Parameter: Function value
2427 C. Author(s): S. Jadach Created at: 01/01/89
2428 C. B. van Eijk Last Update: 02/01/90
2430 C.----------------------------------------------------------------------
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)
2445 C.----------------------------------------------------------------------
2447 C. PHOTOS: PHOton radiation in decays BOost routine '3
'
2449 C. Purpose: Boost vector PVEC along z-axis where ANGLE = EXP(ETA),
2450 C. ETA is the hyperbolic velocity.
2452 C. Input Parameters: ANGLE, PVEC
2454 C. Output Parameter: PVEC
2456 C. Author(s): S. Jadach Created at: 01/01/89
2457 C. B. van Eijk Last Update: 02/01/90
2459 C.----------------------------------------------------------------------
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)
2470 C.----------------------------------------------------------------------
2472 C. PHOTOS: PHOton radiation in decays ROtation routine '2
'
2474 C. Purpose: Rotate x and z components of vector PVEC around angle
2477 C. Input Parameters: ANGLE, PVEC
2479 C. Output Parameter: PVEC
2481 C. Author(s): S. Jadach Created at: 01/01/89
2482 C. B. van Eijk Last Update: 02/01/90
2484 C.----------------------------------------------------------------------
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)
2495 C.----------------------------------------------------------------------
2497 C. PHOTOS: PHOton radiation in decays ROtation routine '3
'
2499 C. Purpose: Rotate x and y components of vector PVEC around angle
2502 C. Input Parameters: ANGLE, PVEC
2504 C. Output Parameter: PVEC
2506 C. Author(s): S. Jadach Created at: 01/01/89
2507 C. B. van Eijk Last Update: 02/01/90
2509 C.----------------------------------------------------------------------
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)
2520 C.----------------------------------------------------------------------
2522 C. PHOTOS: PHOton radiation in decays RANdom number generator init
2524 C. Purpose: Initialse PHORAN with the user specified seeds in the
2525 C. array ISEED. For details see also: F. James CERN DD-
2526 C. Report November 1988.
2528 C. Input Parameters: ISEED(*)
2530 C. Output Parameters: URAN, CRAN, CDRAN, CMRAN, I97, J97
2532 C. Author(s): B. van Eijk and F. James Created at: 27/09/89
2533 C. Last Update: 22/02/90
2535 C.----------------------------------------------------------------------
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
2544 C-- 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)
2554 C-- 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)
2579 C.----------------------------------------------------------------------
2581 C. PHOTOS: PHOton radiation in decays RANdom number generator based
2582 C. on Marsaglia Algorithm
2584 C. Purpose: Generate uniformly distributed random numbers between
2585 C. 0 and 1. Super long period: 2**144. See also:
2586 C. G. Marsaglia and A. Zaman, FSU-SCR-87-50, for seed mo-
2587 C. difications to this version see: F. James DD-Report,
2588 C. November 1988. The generator has to be initialized by
2589 C. a call to PHORIN.
2591 C. Input Parameters: IDUM (integer dummy)
2593 C. Output Parameters: Function value
2595 C. Author(s): B. van Eijk, G. Marsaglia and Created at: 27/09/89
2596 C. A. Zaman Last Update: 27/09/89
2598 C.----------------------------------------------------------------------
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)
2620 C.----------------------------------------------------------------------
2622 C. PHOTOS: PHOton radiation in decays CHArge determination
2624 C. Purpose: Calculate the charge of particle with code IDHEP. The
2625 C. code of the particle is defined by the Particle Data
2626 C. Group in Phys. Lett. B204 (1988) 1.
2628 C. Input Parameter: IDHEP
2630 C. Output Parameter: Funtion value = charge of particle with code
2633 C. Author(s): E. Barberio and B. van Eijk Created at: 29/11/89
2634 C. Last update: 02/01/90
2636 C.----------------------------------------------------------------------
2639 INTEGER IDHEP,IDABS,Q1,Q2,Q3
2641 C-- Array 'charge
' contains the charge of the first 101 particles ac-
2642 C-- 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
2652 C-- Charge of quark, lepton, boson etc....
2653 PHOCHA = CHARGE(IDABS)
2656 C-- 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)
2670 C-- ...diquarks or baryon.
2671 PHOCHA=CHARGE(Q1)+CHARGE(Q2)+CHARGE(Q3)
2675 C-- 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)
2681 C.----------------------------------------------------------------------
2683 C. PHOTOS: PHOton radiation in decays function for SPIn determina-
2686 C. Purpose: Calculate the spin of particle with code IDHEP. The
2687 C. code of the particle is defined by the Particle Data
2688 C. Group in Phys. Lett. B204 (1988) 1.
2690 C. Input Parameter: IDHEP
2692 C. Output Parameter: Funtion value = spin of particle with code
2695 C. Author(s): E. Barberio and B. van Eijk Created at: 29/11/89
2696 C. Last update: 02/01/90
2698 C.----------------------------------------------------------------------
2703 C-- Array 'spin
' contains the spin of the first 100 particles accor-
2704 C-- 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/
2709 C-- Spin of quark, lepton, boson etc....
2710 .LE.
IF (IDABS100) THEN
2714 C-- ...other particles, however...
2715 PHOSPI=(MOD(IDABS,10)-1.D0)/2.D0
2717 C-- ...K_short and K_long are special !!
2718 PHOSPI=MAX(PHOSPI,0.D0)
2722 SUBROUTINE PHOERR(IMES,TEXT,DATA)
2723 C.----------------------------------------------------------------------
2725 C. PHOTOS: PHOton radiation in decays ERRror handling
2727 C. Purpose: Inform user about (fatal) errors and warnings generated
2728 C. by either the user or the program.
2730 C. Input Parameters: IMES, TEXT, DATA
2732 C. Output Parameters: None
2734 C. Author(s): B. van Eijk Created at: 29/11/89
2735 C. Last Update: 10/01/92
2737 C.----------------------------------------------------------------------
2739 DOUBLE PRECISION DATA
2745 PARAMETER (PHOMES=10)
2747 COMMON/PHOSTA/STATUS(PHOMES)
2750 C-- security STOP switch
2755 .LE.
IF (IMESPHOMES) STATUS(IMES)=STATUS(IMES)+1
2757 C-- 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,
2839 c.----------------------------------------------------------------------
2841 c. photos: photon radiation in decays run summary report
2843 c. purpose: inform user about success and/or restrictions of photos
2844 c. encountered during execution.
2846 c. input parameters:
Common /phosta/
2848 c. output parameters: none
2850 c. author(s): b. van eijk created at: 10/01/92
2851 c. last update: 10/01/92
2853 c.----------------------------------------------------------------------
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)
2896 c.----------------------------------------------------------------------
2898 c. phlupa: debugging tool
2900 c. purpose: none, eventually may printout content of the
2903 c. input parameters:
Common /phoevt/ and /phnum/
2904 c. latter may have number of the event.
2906 c. output parameters: none
2908 c. author(s): z. was created at: 30/05/93
2909 c. last update: 09/10/05
2911 c.----------------------------------------------------------------------
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)
2969 c.----------------------------------------------------------------------
2971 c. iphqrk: enables blocks emision from quarks
2974 c. input parameters: modcor
2975 c. modcor >0
type of action
2978 c. =0 execution mode(retrns stored
value)
2981 c. author(s): z. was created at: 11/12/00
2983 c.----------------------------------------------------------------------
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)
3019 c.----------------------------------------------------------------------
3021 c. iphekl: enables/blocks emision in: pi0 to gamma e+ e-
3024 c. input parameters: modcor
3025 c. modcor >0
type of action
3028 c. =0 execution mode(retrns stored
value)
3031 c. author(s): z. was created at: 11/12/00
3033 c.----------------------------------------------------------------------
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)
3073 c.----------------------------------------------------------------------
3075 c. phcork: corrects kinmatics of subbranch needed
if host
program
3076 c. produces events with the shaky momentum conservation
3078 c. input parameters:
Common /phoevt/, modcor
3079 c. modcor >0
type of action
3081 c. =2 corrects energy from mass
3082 c. =3 corrects mass from energy
3083 c. =4 corrects energy from mass for
3084 c. particles up to .4 gev mass,
3085 c. for heavier ones corrects mass,
3086 c. =5 most complete correct also of mother
3087 c. often necessary for exponentiation.
3088 c. =0 execution mode
3090 c. output parameters: corrected /phoevt/
3092 c. author(s): p.golonka, z. was created at: 01/02/99
3093 c. modified : 08/02/99
3094 c.----------------------------------------------------------------------
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
3156 c -----------------------
3157 c in this
case we
do nothing
3159 ELSEIF(modop.EQ.2)
THEN
3160 c -----------------------
3161 cc lets loop thru all daughters and correct their energies
3162 cc 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
3184 c -----------------------
3185 cc lets loop thru all daughters and correct their energies
3186 cc 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
3215 c -----------------------
3217 cc lets loop thru all daughters and correct their masses
3218 cc 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
3241 c -----------------------
3243 cc lets loop thru all daughters and correct their masses
3244 cc 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)
3307 c --- can be used with variant a. for b
use phint1 or 2 --------------
3308 c.----------------------------------------------------------------------
3310 c. phint: photos universal interference correction weight
3312 c. purpose: calculates correction weight as expressed by
3313 c formula(17) from cpc 79 (1994), 291.
3315 c. input parameters:
Common /phoevt/, with photon added.
3317 c. output parameters: correction weight
3319 c. author(s): z. was, p.golonka created at: 19/01/05
3320 c. last update: 25/01/05
3322 c.----------------------------------------------------------------------
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)
3339 c 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
3356 c momenta of charged particle in pl
3360 c 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) )
3371 c accumulate the currents
3375 xdeno = xdeno + xc1**2 + xc2**2
3379 phint=(xnum1**2 + xnum2**2) / xdeno
3385 SUBROUTINE phoeps (VEC1, VEC2, EPS)
3386 c.----------------------------------------------------------------------
3388 c. phoeps: phoeps vector product(normalized to unity)
3390 c. purpose: calculates vector product,
then normalizes its length.
3391 c used to generate orthogonal vectors, i.e. to
3392 c generate polarimetric vectors for photons.
3394 c. input parameters: vec1,vec2 - input 4-vectors
3396 c. output parameters: eps - normalized 4-vector, orthogonal to
3399 c. author(s): z. was, p.golonka created at: 19/01/05
3400 c. last update: 25/01/05
3402 c.----------------------------------------------------------------------
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)
3420 c.----------------------------------------------------------------------
3422 c. photos: photon radiation in decays event dump routine
3424 c. purpose: print event record.
3426 c. input parameters:
Common /hepevt/
3428 c. output parameters: none
3430 c. author(s): b. van eijk created at: 05/06/90
3431 c. last update: 05/06/90
3433 c.----------------------------------------------------------------------
3435 DOUBLE PRECISION SUMVEC(5)
3437 c 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 * ----------------------------------------------------------------------
3463 c-- print event number...
3465 WRITE(phlun,9010) nevhep
3470 c-- 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')