C++ Interface to Tauola
photos-F/photos.f
1 /* copyright(c) 1991-2021 free software foundation, inc.
2  this file is part of the gnu c library.
3 
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.
8 
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.
13 
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/>. */
17 
18 
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. */
26 
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. */
33 
34 
35 
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:
39  - 56 emoji characters
40  - 285 hentaigana
41  - 3 additional Zanabazar Square characters */
42 
43 *///////////////////////////////////////////////////////////////////////
44 *//
45 *// !!!!!!! WARNING!!!!! This source may be agressive !!!!
46 *//
47 *// Due to short common block names it may owerwrite variables in other
48 *// parts of the code.
49 *//
50 *// One should add suffix c_Photos_ to names of all commons as soon as
51 *// possible!!
52 *///////////////////////////////////////////////////////////////////////
53 
54 C.----------------------------------------------------------------------
55 C.
56 C. PHOTOS: PHOtos CDE-s
57 C.
58 C. Purpose: Keep definitions for PHOTOS QED correction Monte Carlo.
59 C.
60 C. Input Parameters: None
61 C.
62 C. Output Parameters: None
63 C.
64 C. Author(s): Z. Was, B. van Eijk Created at: 29/11/89
65 C. Last Update: 10/08/93
66 C.
67 C. =========================================================
68 C. General Structure Information: =
69 C. =========================================================
70 C: ROUTINES:
71 C. 1) INITIALIZATION:
72 C. PHOCDE
73 C. PHOINI
74 C. PHOCIN
75 C. PHOINF
76 C. 2) GENERAL INTERFACE:
77 C. PHOTOS
78 C. PHOTOS_GET
79 C. PHOTOS_SET
80 C. PHOTOS_MAKE
81 C. PHOBOS
82 C. PHOIN
83 C. PHOTWO (specific interface
84 C. PHOOUT
85 C. PHOCHK
86 C. PHTYPE (specific interface
87 C. PHOMAK (specific interface
88 C. 3) QED PHOTON GENERATION:
89 C. PHINT
90 C. PHOBW
91 C. PHOPRE
92 C. PHOOMA
93 C. PHOENE
94 C. PHOCOR
95 C. PHOFAC
96 C. PHODO
97 C. 4) UTILITIES:
98 C. PHOTRI
99 C. PHOAN1
100 C. PHOAN2
101 C. PHOBO3
102 C. PHORO2
103 C. PHORO3
104 C. PHORIN
105 C. PHORAN
106 C. PHOCHA
107 C. PHOSPI
108 C. PHOERR
109 C. PHOREP
110 C. PHLUPA
111 C. PHCORK
112 C. IPHQRK
113 C. IPHEKL
114 C. COMMONS:
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.----------------------------------------------------------------------
134  SUBROUTINE PHOINI
135 C.----------------------------------------------------------------------
136 C.
137 C. PHOTOS: PHOton radiation in decays INItialisation
138 C.
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.
142 C.
143 C. Input Parameters: None
144 C.
145 C. Output Parameters: None
146 C.
147 C. Author(s): Z. Was, B. van Eijk Created at: 26/11/89
148 C. Last Update: 12/04/90
149 C.
150 C.----------------------------------------------------------------------
151  IMPLICIT NONE
152  INTEGER INIT,IDUM,IPHQRK,IPHEKL
153  SAVE INIT
154  DATA INIT/ 0/
155 C--
156 C-- Return if already initialized...
157 .NE. IF (INIT0) RETURN
158  INIT=1
159 C--
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
167  CALL PHCORK(1)
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
175 C--
176 C-- Preset parameters in PHOTOS commons
177  CALL PHOCIN
178 C--
179 C-- Print info
180  CALL PHOINF
181 
182 C--
183 C-- Initialize Marsaglia and Zaman random number generator
184  CALL PHORIN
185  RETURN
186  END
187  SUBROUTINE PHOCIN
188 C.----------------------------------------------------------------------
189 C.
190 C. PHOTOS: PHOton Common INitialisation
191 C.
192 C. Purpose: Initialisation of parameters in common blocks.
193 C.
194 C. Input Parameters: None
195 C.
196 C. Output Parameters: Commons /PHOLUN/, /PHOPHO/, /PHOCOP/, /PHPICO/
197 C. and /PHSEED/.
198 C.
199 C. Author(s): B. van Eijk Created at: 26/11/89
200 C. Z. Was Last Update: 29/01/05
201 C.
202 C.----------------------------------------------------------------------
203  IMPLICIT NONE
204  INTEGER d_h_NMXHEP
205  PARAMETER (d_h_NMXHEP=10000)
206  LOGICAL QEDRAD
207  COMMON/PHOQED/QEDRAD(d_h_NMXHEP)
208  INTEGER PHLUN
209  COMMON/PHOLUN/PHLUN
210  REAL*8 ALPHA,XPHCUT
211  COMMON/PHOCOP/ALPHA,XPHCUT
212  REAL*8 PI,TWOPI
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
217  INTEGER PHOMES
218  PARAMETER (PHOMES=10)
219  INTEGER STATUS
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
224  INTEGER INIT,I
225  SAVE INIT
226  DATA INIT/ 0/
227 C--
228 C-- Return if already initialized...
229 .NE. IF (INIT0) RETURN
230  INIT=1
231 C--
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...
234  DO 10 I=1,d_h_NMXHEP
235  10 QEDRAD(I)=.TRUE.
236 C--
237 C-- Logical output unit for printing of PHOTOS error messages
238  PHLUN=6
239 C--
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
243 C--
244 C-- Define some constants
245  ALPHA=0.00729735039D0
246  PI=3.14159265358979324D0
247  TWOPI=6.28318530717958648D0
248 C--
249 C-- Default seeds Marsaglia and Zaman random number generator
250  ISEED(1)=1802
251  ISEED(2)=9373
252 C--
253 C-- Iitialization for extra options
254 C-- (1)
255 C-- Interference weight now universal.
256  INTERF=.TRUE.
257 C-- (2)
258 C-- Second order - double photon switch
259  ISEC=.TRUE.
260 C-- Third/fourth order - triple (or quatric) photon switch,
261 C-- see dipswitch ifour
262  ITRE=.FALSE.
263 C-- Exponentiation on:
264  IEXP=.FALSE. !.TRUE.
265  IF (IEXP) THEN
266  ISEC=.FALSE.
267  ITRE=.FALSE.
268  CALL PHCORK(5) ! in case of exponentiation correction of ph space
269  ! is a default mandatory
270  XPHCUT=0.000 000 1
271  EXPEPS=1D-4
272  ENDIF
273 C-- (3)
274 C-- Emision in the hard process g g (q qbar) --> t tbar
275 C-- t --> W b
276  IFTOP=.TRUE.
277 C--
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
284  IF (INTERF) THEN
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
288  FINT=2.0D0
289  ELSE
290  FINT=1.0D0
291  ENDIF
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
296 C
297 C IF (INTERF) THEN
298 C FINT=1.80D0
299 C ELSE
300 C FINT=0.0D0
301 C ENDIF
302 C----------END VARIANTS A B -----------------------
303 
304 C-- Effects of initial state charge (in leptonic W decays)
305 C--
306  IFW=.TRUE.
307 C-- Initialise status counter for warning messages
308  DO 20 I=1,PHOMES
309  20 STATUS(I)=0
310  RETURN
311  END
312  SUBROUTINE PHOINF
313 C.----------------------------------------------------------------------
314 C.
315 C. PHOTOS: PHOton radiation in decays general INFo
316 C.
317 C. Purpose: Print PHOTOS info
318 C.
319 C. Input Parameters: PHOLUN
320 C.
321 C. Output Parameters: PHOVN1, PHOVN2
322 C.
323 C. Author(s): B. van Eijk Created at: 12/04/90
324 C. Last Update: 27/06/04
325 C.
326 C.----------------------------------------------------------------------
327  IMPLICIT NONE
328  INTEGER IV1,IV2,IV3
329  INTEGER PHOVN1,PHOVN2
330  COMMON/PHOVER/PHOVN1,PHOVN2
331  INTEGER PHLUN
332  COMMON/PHOLUN/PHLUN
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
336  REAL*8 ALPHA,XPHCUT
337  COMMON/PHOCOP/ALPHA,XPHCUT
338 C--
339 C-- PHOTOS version number and release date
340  PHOVN1=215
341  PHOVN2=111005
342 C--
343 C-- Print info
344  WRITE(PHLUN,9000)
345  WRITE(PHLUN,9020)
346  WRITE(PHLUN,9010)
347  WRITE(PHLUN,9030)
348  IV1=PHOVN1/100
349  IV2=PHOVN1-IV1*100
350  WRITE(PHLUN,9040) IV1,IV2
351  IV1=PHOVN2/10000
352  IV2=(PHOVN2-IV1*10000)/100
353  IV3=PHOVN2-IV1*10000-IV2*100
354  WRITE(PHLUN,9050) IV1,IV2,IV3
355  WRITE(PHLUN,9030)
356  WRITE(PHLUN,9010)
357  WRITE(PHLUN,9060)
358  WRITE(PHLUN,9010)
359  WRITE(PHLUN,9070)
360  WRITE(PHLUN,9010)
361  WRITE(PHLUN,9020)
362  WRITE(PHLUN,9010)
363  WRITE(PHLUN,9064) INTERF,ISEC,ITRE,IEXP,IFTOP,IFW,ALPHA,XPHCUT
364  WRITE(PHLUN,9010)
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)
371  WRITE(PHLUN,9080)
372  WRITE(PHLUN,9010)
373  WRITE(PHLUN,9020)
374  RETURN
375  9000 FORMAT(1H1)
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',
382  &T81,'*')
383  9061 FORMAT(1H ,'*',18X,'option with interference is active ',
384  &T81,'*')
385  9062 FORMAT(1H ,'*',18X,'option with double photons is active ',
386  &T81,'*')
387  9066 FORMAT(1H ,'*',18X,'option with triple/quatric photons is active',
388  &T81,'*')
389  9067 FORMAT(1H ,'*',18X,'option with exponentiation is active epsexp=',
390  &E10.4,T81,'*')
391  9063 FORMAT(1H ,'*',18X,'emision in t tbar production is active ',
392  &T81,'*')
393  9064 FORMAT(1H ,'*',18X,'internal input parameters:',T81,'*'
394  &,/, 1H ,'*',T81,'*'
395  &,/, 1H ,'*',18X,'interf=',L2,' isec=',L2,' itre=',L2,
396  & ' iexp=',L2,' iftop=',L2,
397  & ' ifw=',L2,T81,'*'
398  &,/, 1H ,'*',18X,'alpha_qed=',F8.5,' xphcut=',E8.3,T81,'*')
399  9065 FORMAT(1H ,'*',18X,'correction wt in decay of w is active ',
400  &T81,'*')
401  9070 FORMAT(1H ,'*',9X,
402  &'monte carlo Program - by e. barberio, b. van eijk and z. was',
403  & T81,'*',/,
404  & 1H ,'*',9X,'version 2.09 - by p. golonka and z.w.',T81,'*')
405  9080 FORMAT( 1H ,'*',9X,' ',T81,'*',/,
406  & 1H ,'*',9X,
407  & ' warning(1): /hepevt/ is not anymore the standard common block'
408  & ,T81,'*',/,
409  & 1H ,'*',9X,' ',T81,'*',/,
410  & 1H ,'*',9X,
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'
420  & ,T81,'*')
421  END
422  SUBROUTINE PHOTOS(ID)
423 C.----------------------------------------------------------------------
424 C.
425 C. PHOTOS: General search routine + _GET + _SET
426 C.
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.
430 C.
431 C.
432 C. Input Parameters: ID see routine PHOTOS_MAKE
433 C.
434 C. Output Parameters: None
435 C.
436 C. Author(s): Z. Was Created at: 21/07/98
437 C. Last Update: 21/07/98
438 C.
439 C.----------------------------------------------------------------------
440  COMMON /PHLUPY/ IPOIN,IPOINM
441  INTEGER IPOIN,IPOINM
442  COMMON /PHNUM/ IEV
443  INTEGER IEV
444  INTEGER PHLUN
445  COMMON/PHOLUN/PHLUN
446 
447 .GT..AND..LT. IF (1IPOINM1IPOIN ) THEN
448  WRITE(PHLUN,*) 'event nr=',IEV,
449  $ 'we are testing /hepevt/ at ipoint=1 (input)'
450  CALL PHODMP
451  ENDIF
452  CALL PHOTOS_GET
453  CALL PHOTOS_MAKE(ID)
454  CALL PHOTOS_SET
455 .GT..AND..LT. IF (1IPOINM1IPOIN ) THEN
456  WRITE(PHLUN,*) 'event nr=',IEV,
457  $ 'we are testing /hepevt/ at ipoint=1 (output)'
458  CALL PHODMP
459  ENDIF
460 
461  END
462 
463  SUBROUTINE PHOTOS_GET
464 C.----------------------------------------------------------------------
465 C.
466 C. Getter for PHOTOS:
467 C.
468 C. Purpose: Copies /HEPEVT/ into /PH_HEPEVT/
469 C.
470 C.
471 C. Input Parameters: None
472 C.
473 C. Output Parameters: None
474 C.
475 C. Author(s): Z. Was Created at: 21/07/98
476 C. Last Update: 21/07/98
477 C.
478 C.----------------------------------------------------------------------
479 
480  IMPLICIT NONE
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,
485  $ d_h_jdahep
486  COMMON /hepevt/
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 * ----------------------------------------------------------------------
496  LOGICAL d_h_qedrad
497  COMMON /phoqed/
498  $ d_h_qedrad(d_h_nmxhep) ! Photos flag
499  INTEGER NMXHEP
500  PARAMETER (NMXHEP=10000)
501  INTEGER IDHEP,ISTHEP,JDAHEP,JMOHEP,NEVHEP,NHEP
502  REAL*8 PHEP,VHEP
503  COMMON/PH_HEPEVT/NEVHEP,NHEP,ISTHEP(NMXHEP),IDHEP(NMXHEP),
504  &JMOHEP(2,NMXHEP),JDAHEP(2,NMXHEP),PHEP(5,NMXHEP),VHEP(4,NMXHEP)
505  LOGICAL QEDRAD
506  COMMON/PH_PHOQED/QEDRAD(NMXHEP)
507  integer k,l
508  nevhep= d_h_nevhep ! serial number
509  nhep = d_h_nhep ! number of particles
510  DO K=1,nhep
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
517  DO l=1,4
518  phep(l,k) =d_h_phep(l,k) ! four-momentum, mass [GeV]
519  vhep(l,k) =d_h_vhep(l,k) ! vertex [mm]
520  ENDDO
521  phep(5,k) =d_h_phep(5,k) ! four-momentum, mass [GeV]
522  qedrad(k) =d_h_qedrad(k) ! Photos special flag
523  ENDDO
524  END
525 
526 
527  SUBROUTINE PHOTOS_SET
528 C.----------------------------------------------------------------------
529 C.
530 C. Setter for PHOTOS:
531 C.
532 C. Purpose: Copies /PH_HEPEVT/ into /HEPEVT/
533 C.
534 C.
535 C. Input Parameters: None
536 C.
537 C. Output Parameters: None
538 C.
539 C. Author(s): Z. Was Created at: 21/07/98
540 C. Last Update: 21/07/98
541 C.
542 C.----------------------------------------------------------------------
543  IMPLICIT NONE
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,
548  $ d_h_jdahep
549  COMMON /hepevt/
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 * ----------------------------------------------------------------------
559  LOGICAL d_h_qedrad
560  COMMON /phoqed/
561  $ d_h_qedrad(d_h_nmxhep) ! Photos flag
562  INTEGER NMXHEP
563  PARAMETER (NMXHEP=10000)
564  INTEGER IDHEP,ISTHEP,JDAHEP,JMOHEP,NEVHEP,NHEP
565  REAL*8 PHEP,VHEP
566  COMMON/PH_HEPEVT/NEVHEP,NHEP,ISTHEP(NMXHEP),IDHEP(NMXHEP),
567  &JMOHEP(2,NMXHEP),JDAHEP(2,NMXHEP),PHEP(5,NMXHEP),VHEP(4,NMXHEP)
568  LOGICAL QEDRAD
569  COMMON/PH_PHOQED/QEDRAD(NMXHEP)
570  INTEGER K,L
571 
572  d_h_nevhep= nevhep ! serial number
573  d_h_nhep = nhep ! number of particles
574  DO K=1,nhep
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
581  DO l=1,4
582  d_h_phep(l,k) =phep(l,k) ! four-momentum, mass [GeV]
583  d_h_vhep(l,k) =vhep(l,k) ! vertex [mm]
584  ENDDO
585  d_h_phep(5,k) =phep(5,k) ! four-momentum, mass [GeV]
586  d_h_qedrad(k) =qedrad(k) ! Photos special flag
587  ENDDO
588  END
589  SUBROUTINE PHOTOS_MAKE(IPARR)
590 C.----------------------------------------------------------------------
591 C.
592 C. PHOTOS_MAKE: General search routine
593 C.
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.
599 C.
600 C. Input Parameter: IPPAR: Pointer to decaying particle in
601 C. /PH_HEPEVT/ and the common itself,
602 C.
603 C. Output Parameters: Common /PH_HEPEVT/, either with or without
604 C. new particles added.
605 C.
606 C. Author(s): Z. Was, B. van Eijk Created at: 26/11/89
607 C. Last Update: 30/08/93
608 C.
609 C.----------------------------------------------------------------------
610  IMPLICIT NONE
611  REAL*8 PHOTON(5)
612  INTEGER IP,IPARR,IPPAR,I,J,K,L,NLAST
613  DOUBLE PRECISION DATA
614  INTEGER MOTHER,POSPHO
615  LOGICAL CASCAD
616  INTEGER NMXHEP
617  PARAMETER (NMXHEP=10000)
618  INTEGER IDHEP,ISTHEP,JDAHEP,JMOHEP,NEVHEP,NHEP
619  REAL*8 PHEP,VHEP
620  COMMON/PH_HEPEVT/NEVHEP,NHEP,ISTHEP(NMXHEP),IDHEP(NMXHEP),
621  &JMOHEP(2,NMXHEP),JDAHEP(2,NMXHEP),PHEP(5,NMXHEP),VHEP(4,NMXHEP)
622  LOGICAL QEDRAD
623  COMMON/PH_PHOQED/QEDRAD(NMXHEP)
624  INTEGER NMXPHO
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)
629 C--
630  IPPAR=ABS(IPARR)
631 C-- Store pointers for cascade treatement...
632  IP=IPPAR
633  NLAST=NHEP
634  CASCAD=.FALSE.
635 C--
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
638 C--
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
642  ISTACK(0)=IPPAR
643 C-- NUMIT denotes number of secondary decay branches
644  NUMIT=0
645 C-- NTRY denotes number of secondary branches already checked for
646 C-- for existence of further branches
647  NTRY=0
648 C-- let-s search if IPARR does not prevent searching.
649 .GT. IF (IPARR0) THEN
650  30 CONTINUE
651  DO I=JDAHEP(1,IP),JDAHEP(2,IP)
652 .NE..AND..EQ. IF (JDAHEP(1,I)0JMOHEP(1,JDAHEP(1,I))I) THEN
653  NUMIT=NUMIT+1
654 .GT. IF (NUMITNMXPHO) THEN
655  DATA=NUMIT
656  CALL PHOERR(7,'photos',DATA)
657  ENDIF
658  ISTACK(NUMIT)=I
659  ENDIF
660  ENDDO
661 .GT. IF(NUMITNTRY) THEN
662  NTRY=NTRY+1
663  IP=ISTACK(NTRY)
664  GOTO 30
665  ENDIF
666  ENDIF
667 C-- let-s do generation
668  DO 25 KK=0,NUMIT
669  NA=NHEP
670  FIRST=JDAHEP(1,ISTACK(KK))
671  LAST=JDAHEP(2,ISTACK(KK))
672  DO II=1,LAST-FIRST+1
673  DO LL=1,5
674  PORIG(LL,II)=PHEP(LL,FIRST+II-1)
675  ENDDO
676  ENDDO
677 C--
678  CALL PHTYPE(ISTACK(KK))
679 C--
680 C-- Correct energy/momentum of cascade daughters
681 .GT. IF(NHEPNA) THEN
682  DO II=1,LAST-FIRST+1
683  IPP=FIRST+II-1
684  FIRSTA=JDAHEP(1,IPP)
685  LASTA=JDAHEP(2,IPP)
686 .EQ. IF(JMOHEP(1,IPP)ISTACK(KK))
687  $ CALL PHOBOS(IPP,PORIG(1,II),PHEP(1,IPP),FIRSTA,LASTA)
688  ENDDO
689  ENDIF
690  25 CONTINUE
691 C--
692 C-- rearrange /PH_HEPEVT/ to get correct order..
693 .GT. IF (NHEPNLAST) THEN
694  DO 160 I=NLAST+1,NHEP
695 C--
696 C-- Photon mother and position...
697  MOTHER=JMOHEP(1,I)
698  POSPHO=JDAHEP(2,MOTHER)+1
699 C-- Intermediate save of photon energy/momentum and pointers
700  DO 90 J=1,5
701  90 PHOTON(J)=PHEP(J,I)
702  ISPHO =ISTHEP(I)
703  IDPHO =IDHEP(I)
704  MOTHER2 =JMOHEP(2,I)
705  IDA1 =JDAHEP(1,I)
706  IDA2 =JDAHEP(2,I)
707 C--
708 C-- Exclude photon in sequence !
709 .NE. IF (POSPHONHEP) THEN
710 C--
711 C--
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)
716  IDHEP(K)=IDHEP(K-1)
717  DO 100 L=1,2
718  JMOHEP(L,K)=JMOHEP(L,K-1)
719  100 JDAHEP(L,K)=JDAHEP(L,K-1)
720  DO 110 L=1,5
721  110 PHEP(L,K)=PHEP(L,K-1)
722  DO 120 L=1,4
723  120 VHEP(L,K)=VHEP(L,K-1)
724 C--
725 C-- Correct pointers assuming most dirty /PH_HEPEVT/...
726  DO 130 K=1,NHEP
727  DO 130 L=1,2
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
732  130 CONTINUE
733 C--
734 C-- Store photon energy/momentum
735  DO 140 J=1,5
736  140 PHEP(J,POSPHO)=PHOTON(J)
737  ENDIF
738 C--
739 C-- Store pointers for the photon...
740  JDAHEP(2,MOTHER)=POSPHO
741  ISTHEP(POSPHO)=ISPHO
742  IDHEP(POSPHO)=IDPHO
743  JMOHEP(1,POSPHO)=MOTHER
744  JMOHEP(2,POSPHO)=MOTHER2
745  JDAHEP(1,POSPHO)=IDA1
746  JDAHEP(2,POSPHO)=IDA2
747 C--
748 C-- Get photon production vertex position
749  DO 150 J=1,4
750  150 VHEP(J,POSPHO)=VHEP(J,POSPHO-1)
751  160 CONTINUE
752  ENDIF
753  RETURN
754  END
755  SUBROUTINE PHOBOS(IP,PBOOS1,PBOOS2,FIRST,LAST)
756 C.----------------------------------------------------------------------
757 C.
758 C. PHOBOS: PHOton radiation in decays BOoSt routine
759 C.
760 C. Purpose: Boost particles in cascade decay to parent rest frame
761 C. and boost back with modified boost vector.
762 C.
763 C. Input Parameters: IP: pointer of particle starting chain
764 C. to be boosted
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/).
771 C.
772 C. Output Parameters: Common /PH_HEPEVT/.
773 C.
774 C. Author(s): B. van Eijk Created at: 13/02/90
775 C. Z. Was Last Update: 16/11/93
776 C.
777 C.----------------------------------------------------------------------
778  IMPLICIT NONE
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)
784  INTEGER NMXHEP
785  PARAMETER (NMXHEP=10000)
786  INTEGER IDHEP,ISTHEP,JDAHEP,JMOHEP,NEVHEP,NHEP
787  REAL*8 PHEP,VHEP
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
791  NSTACK=0
792  DO 10 J=1,3
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)
797 C--
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
802  DO 30 J=1,3
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
805 C--
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)
808  DO 40 J=1,3
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
812  NSTACK=NSTACK+1
813 C--
814 C-- Check on stack length...
815 .GT. IF (NSTACKMAXSTA) THEN
816  DATA=NSTACK
817  CALL PHOERR(7,'phobos',DATA)
818  ENDIF
819  STACK(NSTACK)=I
820  ENDIF
821  ENDIF
822  50 CONTINUE
823 .NE. IF (NSTACK0) THEN
824 C--
825 C-- Now go one step further in the decay tree...
826  FIRST=JDAHEP(1,STACK(NSTACK))
827  LAST=JDAHEP(2,STACK(NSTACK))
828  IP=STACK(NSTACK)
829  NSTACK=NSTACK-1
830  GOTO 20
831  ENDIF
832  RETURN
833  END
834  SUBROUTINE PHOIN(IP,BOOST,NHEP0)
835 C.----------------------------------------------------------------------
836 C.
837 C. PHOIN: PHOtos INput
838 C.
839 C. Purpose: copies IP branch of the common /PH_HEPEVT/ into /PHOEVT/
840 C. moves branch into its CMS system.
841 C.
842 C. Input Parameters: IP: pointer of particle starting branch
843 C. to be copied
844 C. BOOST: Flag whether boost to CMS was or was
845 C . not performed.
846 C.
847 C. Output Parameters: Commons: /PHOEVT/, /PHOCMS/
848 C.
849 C. Author(s): Z. Was Created at: 24/05/93
850 C. Last Update: 16/11/93
851 C.
852 C.----------------------------------------------------------------------
853  IMPLICIT NONE
854  INTEGER NMXHEP
855  PARAMETER (NMXHEP=10000)
856  INTEGER IDHEP,ISTHEP,JDAHEP,JMOHEP,NEVHEP,NHEP
857  REAL*8 PHEP,VHEP
858  COMMON/PH_HEPEVT/NEVHEP,NHEP,ISTHEP(NMXHEP),IDHEP(NMXHEP),
859  &JMOHEP(2,NMXHEP),JDAHEP(2,NMXHEP),PHEP(5,NMXHEP),VHEP(4,NMXHEP)
860  INTEGER NMXPHO
861  PARAMETER (NMXPHO=10000)
862  INTEGER IDPHO,ISTPHO,JDAPHO,JMOPHO,NEVPHO,NPHO
863  REAL*8 PPHO,VPHO
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
867  LOGICAL BOOST
868  INTEGER J,NHEP0
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
874 C--
875 C let-s calculate size of the little common entry
876  FIRST=JDAHEP(1,IP)
877  LAST =JDAHEP(2,IP)
878  NPHO=3+LAST-FIRST+NHEP-NHEP0
879  NEVPHO=NPHO
880 C let-s take in decaying particle
881  IDPHO(1)=IDHEP(IP)
882  JDAPHO(1,1)=3
883  JDAPHO(2,1)=3+LAST-FIRST
884  DO I=1,5
885  PPHO(I,1)=PHEP(I,IP)
886  ENDDO
887 C let-s take in eventual second mother
888  IP2=JMOHEP(2,JDAHEP(1,IP))
889 .NE..AND..NE. IF((IP20)(IP2IP)) THEN
890  IDPHO(2)=IDHEP(IP2)
891  JDAPHO(1,2)=3
892  JDAPHO(2,2)=3+LAST-FIRST
893  DO I=1,5
894  PPHO(I,2)=PHEP(I,IP2)
895  ENDDO
896  ELSE
897  IDPHO(2)=0
898  DO I=1,5
899  PPHO(I,2)=0.0D0
900  ENDDO
901  ENDIF
902 C let-s take in daughters
903  DO LL=0,LAST-FIRST
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
907  DO I=1,5
908  PPHO(I,3+LL)=PHEP(I,FIRST+LL)
909  ENDDO
910  ENDDO
911 .GT. IF (NHEPNHEP0) THEN
912 C let-s take in illegitimate daughters
913  NA=3+LAST-FIRST
914  DO LL=1,NHEP-NHEP0
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
918  DO I=1,5
919  PPHO(I,NA+LL)=PHEP(I,NHEP0+LL)
920  ENDDO
921  ENDDO
922 C-- there is NHEP-NHEP0 daugters more.
923  JDAPHO(2,1)=3+LAST-FIRST+NHEP-NHEP0
924  ENDIF
925 .EQ. IF(IDPHO(NPHO)22)CALL PHLUPA(100001)
926 .EQ.! IF(IDPHO(NPHO)22) stop
927  CALL PHCORK(0)
928 .EQ. IF(IDPHO(NPHO)22)CALL PHLUPA(100002)
929 C special case of t tbar production process
930  IF(IFTOP) CALL PHOTWO(0)
931  BOOST=.FALSE.
932 C-- Check whether parent is in its rest frame...
933 C ZBW ND 27.07.2009:
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.
940 
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
943 
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
946 
947  BOOST=.TRUE.
948 C--
949 C-- Boost daughter particles to rest frame of parent...
950 C-- Resultant neutral system already calculated in rest frame !
951  DO 10 J=1,3
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)
956  DO 20 J=1,3
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
960  I=1
961  PB=BET(1)*PPHO(1,I)+BET(2)*PPHO(2,I)+BET(3)*PPHO(3,I)
962  DO J=1,3
963  PPHO(J,I)=PPHO(J,I)+BET(J)*(PPHO(4,I)+PB/(GAM+1.D0))
964  ENDDO
965  PPHO(4,I)=GAM*PPHO(4,I)+PB
966  ENDIF
967 C special case of t tbar production process
968  IF(IFTOP) CALL PHOTWO(1)
969  CALL PHLUPA(2)
970 .EQ. IF(IDPHO(NPHO)22) CALL PHLUPA(10000)
971 .EQ.! IF(IDPHO(NPHO-1)22) stop
972  END
973  SUBROUTINE PHOTWO(MODE)
974 C.----------------------------------------------------------------------
975 C.
976 C. PHOTWO: PHOtos but TWO mothers allowed
977 C.
978 C. Purpose: Combines two mothers into one in /PHOEVT/
979 C. necessary eg in case of g g (q qbar) --> t tbar
980 C.
981 C. Input Parameters: Common /PHOEVT/ (/PHOCMS/)
982 C.
983 C. Output Parameters: Common /PHOEVT/, (stored mothers)
984 C.
985 C. Author(s): Z. Was Created at: 5/08/93
986 C. Last Update:10/08/93
987 C.
988 C.----------------------------------------------------------------------
989  IMPLICIT NONE
990  INTEGER NMXPHO
991  PARAMETER (NMXPHO=10000)
992  INTEGER IDPHO,ISTPHO,JDAPHO,JMOPHO,NEVPHO,NPHO
993  REAL*8 PPHO,VPHO
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
998  INTEGER I,MODE
999  REAL*8 MPASQR
1000  LOGICAL IFRAD
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
1005 C
1006 C t tbar case
1007 .EQ. IF(MODE0) THEN
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)
1010  IFRAD=IFRAD
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)
1015  IF(IFRAD) THEN
1016 c.....combining first and second mother
1017  DO I=1,4
1018  PPHO(I,1)=PPHO(I,1)+PPHO(I,2)
1019  ENDDO
1020  PPHO(5,1)=SQRT(MPASQR)
1021 c.....removing second mother,
1022  DO I=1,5
1023  PPHO(I,2)=0.0D0
1024  ENDDO
1025  ENDIF
1026  ELSE
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.
1030  ENDIF
1031  END
1032  SUBROUTINE PHOOUT(IP,BOOST,NHEP0)
1033 C.----------------------------------------------------------------------
1034 C.
1035 C. PHOOUT: PHOtos OUTput
1036 C.
1037 C. Purpose: copies back IP branch of the common /PH_HEPEVT/ from
1038 C. /PHOEVT/ moves branch back from its CMS system.
1039 C.
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
1043 C . not performed.
1044 C.
1045 C. Output Parameters: Common /PHOEVT/,
1046 C.
1047 C. Author(s): Z. Was Created at: 24/05/93
1048 C. Last Update:
1049 C.
1050 C.----------------------------------------------------------------------
1051  IMPLICIT NONE
1052  INTEGER NMXHEP
1053  PARAMETER (NMXHEP=10000)
1054  INTEGER IDHEP,ISTHEP,JDAHEP,JMOHEP,NEVHEP,NHEP
1055  REAL*8 PHEP,VHEP
1056  COMMON/PH_HEPEVT/NEVHEP,NHEP,ISTHEP(NMXHEP),IDHEP(NMXHEP),
1057  &JMOHEP(2,NMXHEP),JDAHEP(2,NMXHEP),PHEP(5,NMXHEP),VHEP(4,NMXHEP)
1058  INTEGER NMXPHO
1059  PARAMETER (NMXPHO=10000)
1060  INTEGER IDPHO,ISTPHO,JDAPHO,JMOPHO,NEVPHO,NPHO
1061  REAL*8 PPHO,VPHO
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
1065  LOGICAL BOOST
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...
1071  CALL PHLUPA(10)
1072  IF (BOOST) THEN
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)
1075  DO 100 K=1,3
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
1079  DO NN=NEVPHO+1,NPHO
1080  PB=-BET(1)*PPHO(1,NN)-BET(2)*PPHO(2,NN)-BET(3)*PPHO(3,NN)
1081  DO 120 K=1,3
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
1084  ENDDO
1085  ENDIF
1086  FIRST=JDAHEP(1,IP)
1087  LAST =JDAHEP(2,IP)
1088 C let-s take in original daughters
1089  DO LL=0,LAST-FIRST
1090  IDHEP(FIRST+LL) = IDPHO(3+LL)
1091  DO I=1,5
1092  PHEP(I,FIRST+LL) = PPHO(I,3+LL)
1093  ENDDO
1094  ENDDO
1095 C let-s take newcomers to the end of HEPEVT.
1096  NA=3+LAST-FIRST
1097  DO LL=1,NPHO-NA
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
1104  DO I=1,5
1105  PHEP(I,NHEP0+LL) = PPHO(I,NA+LL)
1106  ENDDO
1107  ENDDO
1108  NHEP=NHEP+NPHO-NEVPHO
1109  CALL PHLUPA(20)
1110  END
1111  SUBROUTINE PHOCHK(JFIRST)
1112 C.----------------------------------------------------------------------
1113 C.
1114 C. PHOCHK: checking branch.
1115 C.
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.
1120 C.
1121 C.
1122 C. Author(s): Z. Was Created at: 22/10/92
1123 C. Last Update: 11/12/00
1124 C.
1125 C.----------------------------------------------------------------------
1126 C ********************
1127  IMPLICIT NONE
1128  INTEGER NMXPHO
1129  PARAMETER (NMXPHO=10000)
1130  INTEGER IDPHO,ISTPHO,JDAPHO,JMOPHO,NEVPHO,NPHO
1131  REAL*8 PPHO,VPHO
1132  COMMON/PHOEVT/NEVPHO,NPHO,ISTPHO(NMXPHO),IDPHO(NMXPHO),
1133  &JMOPHO(2,NMXPHO),JDAPHO(2,NMXPHO),PPHO(5,NMXPHO),VPHO(4,NMXPHO)
1134  LOGICAL CHKIF
1135  COMMON/PHOIF/CHKIF(NMXPHO)
1136  INTEGER NMXHEP
1137  PARAMETER (NMXHEP=10000)
1138  LOGICAL QEDRAD
1139  COMMON/PH_PHOQED/QEDRAD(NMXHEP)
1140  INTEGER JFIRST
1141  LOGICAL F
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
1146  LOGICAL IFRAD
1147  INTEGER IDENT,K,IQRK,IPHQRK,IEKL,IPHEKL
1148 C these are OK .... if you do not like somebody else, add here.
1149  F(IDABS)=
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)
1156 C
1157  IQRK=IPHQRK(0) ! switch for emission from quark
1158  IEKL=IPHEKL(0)
1159  NLAST = NPHO
1160 C
1161  IPPAR=1
1162 C checking for good particles
1163  IFNPI0=.TRUE.
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)) )
1172 
1173 .AND..NOT. IFNPI0=(IFNPI0(IFKL))
1174  ENDIF
1175  DO 10 I=IPPAR,NLAST
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)
1181 .AND. & IFNPI0
1182  10 CONTINUE
1183 C--
1184 C now we go to special cases, where CHKIF(I) will be overwritten
1185 C--
1186  IF(IFTOP) THEN
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
1190  IDENT=K
1191  GOTO 15
1192  ENDIF
1193  ENDDO
1194  15 CONTINUE
1195 .EQ..AND..EQ. IFRAD=((IDPHO(1)21)(IDPHO(2)21))
1196 .OR..LE..AND..EQ. & ((ABS(IDPHO(1))6)((IDPHO(2))(-IDPHO(1))))
1197  IFRAD=IFRAD
1198 .AND..EQ..AND..EQ. & (ABS(IDPHO(3))6)((IDPHO(4))(-IDPHO(3)))
1199 .AND..EQ. & (IDENT4)
1200  IF(IFRAD) THEN
1201  DO 20 I=IPPAR,NLAST
1202  CHKIF(I)= .TRUE.
1203 .GT..AND. IF(I2) CHKIF(I)=CHKIF(I)QEDRAD(JFIRST+I-IPPAR-2)
1204  20 CONTINUE
1205  ENDIF
1206  ENDIF
1207 C--
1208 C--
1209  IF(IFTOP) THEN
1210 C special case of top decay
1211  DO K=JDAPHO(2,1),JDAPHO(1,1),-1
1212 .NE. IF(IDPHO(K)22) THEN
1213  IDENT=K
1214  GOTO 25
1215  ENDIF
1216  ENDDO
1217  25 CONTINUE
1218 .EQ..AND..EQ. IFRAD=((ABS(IDPHO(1))6)(IDPHO(2)0))
1219  IFRAD=IFRAD
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)
1223  IF(IFRAD) THEN
1224  DO 30 I=IPPAR,NLAST
1225  CHKIF(I)= .TRUE.
1226 .GT..AND. IF(I2) CHKIF(I)=CHKIF(I)QEDRAD(JFIRST+I-IPPAR-2)
1227  30 CONTINUE
1228  ENDIF
1229  ENDIF
1230 C--
1231 C--
1232  END
1233  SUBROUTINE PHTYPE(ID)
1234 C.----------------------------------------------------------------------
1235 C.
1236 C. PHTYPE: Central manadgement routine.
1237 C.
1238 C. Purpose: defines what kind of the
1239 C. actions will be performed at point ID.
1240 C.
1241 C. Input Parameters: ID: pointer of particle starting branch
1242 C. in /PH_HEPEVT/ to be treated.
1243 C.
1244 C. Output Parameters: Common /PH_HEPEVT/.
1245 C.
1246 C. Author(s): Z. Was Created at: 24/05/93
1247 C. P. Golonka Last Update: 27/06/04
1248 C.
1249 C.----------------------------------------------------------------------
1250  IMPLICIT NONE
1251  INTEGER NMXHEP
1252  PARAMETER (NMXHEP=10000)
1253  INTEGER IDHEP,ISTHEP,JDAHEP,JMOHEP,NEVHEP,NHEP
1254  REAL*8 PHEP,VHEP
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
1260  LOGICAL EXPINI
1261  INTEGER NX,K,NCHAN
1262  PARAMETER (NX=10)
1263  REAL*8 PRO,PRSUM,ESU
1264  COMMON /PHOEXP/ PRO(NX),NCHAN,EXPINI
1265 
1266  INTEGER ID,NHEP0
1267  LOGICAL IPAIR
1268  REAL*8 RN,PHORAN,SUM
1269  INTEGER WTDUM
1270  LOGICAL IFOUR
1271 C--
1272 .AND. IFOUR=(.TRUE.)(ITRE) ! we can make internal choice whether
1273  ! we want 3 or four photons at most.
1274  IPAIR=.TRUE.
1275 C-- Check decay multiplicity..
1276 .EQ. IF (JDAHEP(1,ID)0) RETURN
1277 .EQ.C IF (JDAHEP(1,ID)JDAHEP(2,ID)) RETURN
1278 C--
1279  NHEP0=NHEP
1280 C--
1281  IF (IEXP) THEN
1282  EXPINI=.TRUE. ! Initialization/cleaning
1283  DO NCHAN=1,NX
1284  PRO(NCHAN)=0.D0
1285  ENDDO
1286  NCHAN=0
1287 
1288  FSEC=1.0D0
1289  CALL PHOMAK(ID,NHEP0)! Initialization/crude formfactors into
1290  ! PRO(NCHAN)
1291  EXPINI=.FALSE.
1292  RN=PHORAN(WTDUM)
1293  PRSUM=0
1294  DO K=1,NX
1295  PRSUM=PRSUM+PRO(K)
1296  ENDDO
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
1301  ! at first for k=0
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.
1306  NCHAN=0
1307  CALL PHOMAK(ID,NHEP0) ! LOOPING
1308 .GT. IF(SUM1D0-EXPEPS) GOTO 100
1309  ENDDO
1310  100 CONTINUE
1311  ELSEIF(IFOUR) THEN
1312 C-- quatro photon emission
1313  FSEC=1.0D0
1314  RN=PHORAN(WTDUM)
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)
1325  ENDIF
1326  ELSEIF(ITRE) THEN
1327 C-- triple photon emission
1328  FSEC=1.0D0
1329  RN=PHORAN(WTDUM)
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)
1336  ENDIF
1337  ELSEIF(ISEC) THEN
1338 C-- double photon emission
1339  FSEC=1.0D0
1340  RN=PHORAN(WTDUM)
1341 .GE. IF (RN0.5D0) THEN
1342  CALL PHOMAK(ID,NHEP0)
1343  CALL PHOMAK(ID,NHEP0)
1344  ENDIF
1345  ELSE
1346 C-- single photon emission
1347  FSEC=1.0D0
1348  CALL PHOMAK(ID,NHEP0)
1349  ENDIF
1350 C--
1351 C-- electron positron pair (coomented out for a while
1352 C IF (IPAIR) CALL PHOPAR(ID,NHEP0)
1353  END
1354  SUBROUTINE PHOMAK(IPPAR,NHEP0)
1355 C.----------------------------------------------------------------------
1356 C.
1357 C. PHOMAK: PHOtos MAKe
1358 C.
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
1362 C. general tools.
1363 C.
1364 C. Input Parameter: IPPAR: Pointer to decaying particle in
1365 C. /PH_HEPEVT/ and the common itself
1366 C.
1367 C. Output Parameters: Common /PH_HEPEVT/, either with or without
1368 C. particles added.
1369 C.
1370 C. Author(s): Z. Was, Created at: 26/05/93
1371 C. Last Update: 29/01/05
1372 C.
1373 C.----------------------------------------------------------------------
1374  IMPLICIT NONE
1375  DOUBLE PRECISION DATA
1376  REAL*8 PHORAN
1377  INTEGER IP,IPPAR,NCHARG
1378  INTEGER WTDUM,IDUM,NHEP0
1379  INTEGER NCHARB,NEUDAU
1380  REAL*8 RN,WT,PHINT
1381  LOGICAL BOOST
1382  INTEGER NMXHEP
1383  PARAMETER (NMXHEP=10000)
1384  INTEGER IDHEP,ISTHEP,JDAHEP,JMOHEP,NEVHEP,NHEP
1385  REAL*8 PHEP,VHEP
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
1391 C--
1392  IP=IPPAR
1393  IDUM=1
1394  NCHARG=0
1395 C--
1396  CALL PHOIN(IP,BOOST,NHEP0)
1397  CALL PHOCHK(JDAHEP(1,IP))
1398  WT=0.0D0
1399  CALL PHOPRE(1,WT,NEUDAU,NCHARB)
1400 
1401 .EQ. IF (WT0.0D0) RETURN
1402  RN=PHORAN(WTDUM)
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
1408  DATA=WT
1409 .GT. IF (WT1.0D0) CALL PHOERR(3,'wt_int',DATA)
1410 C weighting
1411 .LE. IF (RNWT) THEN
1412  CALL PHOOUT(IP,BOOST,NHEP0)
1413  ENDIF
1414  RETURN
1415  END
1416  FUNCTION PHINT1(IDUM)
1417 C.----------------------------------------------------------------------
1418 C.
1419 C. PHINT: PHotos INTerference (Old version kept for tests only.
1420 C.
1421 C. Purpose: Calculates interference between emission of photons from
1422 C. different possible chaged daughters stored in
1423 C. the HEP common /PHOEVT/.
1424 C.
1425 C. Input Parameter: commons /PHOEVT/ /PHOMOM/ /PHOPHS/
1426 C.
1427 C.
1428 C. Output Parameters:
1429 C.
1430 C.
1431 C. Author(s): Z. Was, Created at: 10/08/93
1432 C. Last Update: 15/03/99
1433 C.
1434 C.----------------------------------------------------------------------
1435 
1436  IMPLICIT NONE
1437  REAL*8 PHINT,phint1
1438  REAL*8 PHOCHA
1439  INTEGER IDUM
1440  INTEGER NMXPHO
1441  PARAMETER (NMXPHO=10000)
1442  INTEGER IDPHO,ISTPHO,JDAPHO,JMOPHO,NEVPHO,NPHO
1443  REAL*8 PPHO,VPHO
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
1447  REAL*8 PNEUTR
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
1453  LOGICAL IFINT
1454  INTEGER K,IDENT
1455 C
1456  DO K=JDAPHO(2,1),JDAPHO(1,1),-1
1457 .NE. IF(IDPHO(K)22) THEN
1458  IDENT=K
1459  GOTO 20
1460  ENDIF
1461  ENDDO
1462  20 CONTINUE
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
1472  IF(IFINT) THEN
1473  MPASQR = PPHO(5,1)**2
1474  XX=4.D0*MCHSQR/MPASQR*(1.D0-XPHOTO)/(1.D0-XPHOTO+(MCHSQR-MNESQR)
1475  & /MPASQR)**2
1476  BETA=SQRT(1.D0-XX)
1477  PHINT = 2D0/(1D0+COSTHG**2*BETA**2)
1478  ELSE
1479  PHINT = 1D0
1480  ENDIF
1481  phint1=1
1482  END
1483 
1484  FUNCTION PHINT2(IDUM)
1485 C.----------------------------------------------------------------------
1486 C.
1487 C. PHINT: PHotos INTerference
1488 C.
1489 C. Purpose: Calculates interference between emission of photons from
1490 C. different possible chaged daughters stored in
1491 C. the HEP common /PHOEVT/.
1492 C.
1493 C. Input Parameter: commons /PHOEVT/ /PHOMOM/ /PHOPHS/
1494 C.
1495 C.
1496 C. Output Parameters:
1497 C.
1498 C.
1499 C. Author(s): Z. Was, Created at: 10/08/93
1500 C. Last Update:
1501 C.
1502 C.----------------------------------------------------------------------
1503 
1504  IMPLICIT NONE
1505  REAL*8 PHINT,PHINT1,PHINT2
1506  REAL*8 PHOCHA
1507  INTEGER IDUM
1508  INTEGER NMXPHO
1509  PARAMETER (NMXPHO=10000)
1510  INTEGER IDPHO,ISTPHO,JDAPHO,JMOPHO,NEVPHO,NPHO
1511  REAL*8 PPHO,VPHO
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
1515  REAL*8 PNEUTR
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
1522  LOGICAL IFINT
1523  INTEGER K,IDENT
1524 C
1525  DO K=JDAPHO(2,1),JDAPHO(1,1),-1
1526 .NE. IF(IDPHO(K)22) THEN
1527  IDENT=K
1528  GOTO 20
1529  ENDIF
1530  ENDDO
1531  20 CONTINUE
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
1541 .AND. IFINT= IFINT
1542 .gt. $ abs(PHOCHA(IDPHO(JDAPHO(1,1))))0.01D0
1543 C calculates interference weight contribution
1544  IF(IFINT) THEN
1545  MPASQR = PPHO(5,1)**2
1546  XX=4.D0*MCHSQR/MPASQR*(1.-XPHOTO)/(1.-XPHOTO+(MCHSQR-MNESQR)/
1547  & MPASQR)**2
1548  BETA=SQRT(1.D0-XX)
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
1552  PP =SQRT(PP2)
1553  E1 =SQRT(PP2+MCHSQR)
1554  E2 =SQRT(PP2+MNESQR)
1555  PHINT= (E1+E2)**2/((E2+COSTHG*PP)**2+(E1-COSTHG*PP)**2)
1556 C
1557  q1=PHOCHA(IDPHO(JDAPHO(1,1)))
1558  q2=PHOCHA(IDPHO(IDENT))
1559  do k=1,4
1560  pq1(k)=ppho(k,JDAPHO(1,1))
1561  pq2(k)=ppho(k,JDAPHO(1,1)+1)
1562  pphot(k)=ppho(k,npho)
1563  enddo
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)
1567 C
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
1573 
1574  PHINT= (q1*(E2+COSTHG*PP)-q2*(E1-COSTHG*PP))**2
1575  & /(q1**2*(E2+COSTHG*PP)**2+q2**2*(E1-COSTHG*PP)**2)
1576  ELSE
1577 
1578  PHINT= (q1*(E1-COSTHG*PP)-q2*(E2+COSTHG*PP))**2
1579  & /(q1**2*(E1-COSTHG*PP)**2+q2**2*(E2+COSTHG*PP)**2)
1580  ENDIF
1581  ELSE
1582  PHINT = 1D0
1583  ENDIF
1584  phint1=1
1585  phint2=1
1586  END
1587 
1588 
1589  SUBROUTINE PHOPRE(IPARR,WT,NEUDAU,NCHARB)
1590 C.----------------------------------------------------------------------
1591 C.
1592 C. PHOTOS: Photon radiation in decays
1593 C.
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.
1600 C.
1601 C. Input Parameter: IPPAR: Pointer to decaying particle in
1602 C. /PHOEVT/ and the common itself,
1603 C.
1604 C. Output Parameters: Common /PHOEVT/, either with or without a
1605 C. photon(s) added.
1606 C. WT weight of the configuration
1607 C.
1608 C. Author(s): Z. Was, B. van Eijk Created at: 26/11/89
1609 C. Last Update: 29/01/05
1610 C.
1611 C.----------------------------------------------------------------------
1612  IMPLICIT NONE
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
1617  INTEGER IDABS,IDUM
1618  INTEGER NCHARB,NEUDAU
1619  REAL*8 WT,WGT
1620  INTEGER NMXPHO
1621  PARAMETER (NMXPHO=10000)
1622  INTEGER IDPHO,ISTPHO,JDAPHO,JMOPHO,NEVPHO,NPHO
1623  REAL*8 PPHO,VPHO
1624  COMMON/PHOEVT/NEVPHO,NPHO,ISTPHO(NMXPHO),IDPHO(NMXPHO),
1625  &JMOPHO(2,NMXPHO),JDAPHO(2,NMXPHO),PPHO(5,NMXPHO),VPHO(4,NMXPHO)
1626  LOGICAL CHKIF
1627  COMMON/PHOIF/CHKIF(NMXPHO)
1628  INTEGER CHAPOI(NMXPHO)
1629  DOUBLE PRECISION MCHSQR,MNESQR
1630  REAL*8 PNEUTR
1631  COMMON/PHOMOM/MCHSQR,MNESQR,PNEUTR(5)
1632  DOUBLE PRECISION COSTHG,SINTHG
1633  REAL*8 XPHMAX,XPHOTO
1634  COMMON/PHOPHS/XPHMAX,XPHOTO,COSTHG,SINTHG
1635  REAL*8 ALPHA,XPHCUT
1636  COMMON/PHOCOP/ALPHA,XPHCUT
1637  INTEGER IREP
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
1644 
1645 C--
1646  IPPAR=IPARR
1647 C-- Store pointers for cascade treatement...
1648  IP=IPPAR
1649  NLAST=NPHO
1650  IDUM=1
1651 C--
1652 C-- Check decay multiplicity..
1653 .EQ. IF (JDAPHO(1,IP)0) RETURN
1654 C--
1655 C-- Loop over daughters, determine charge multiplicity
1656  10 NCHARG=0
1657  IREP=0
1658  MINMAS=0.D0
1659  MASSUM=0.D0
1660  DO 20 I=JDAPHO(1,IP),JDAPHO(2,IP)
1661 C--
1662 C--
1663 C-- Exclude marked particles, quarks and gluons etc...
1664  IDABS=ABS(IDPHO(I))
1665  IF (CHKIF(I-JDAPHO(1,IP)+3)) THEN
1666 .NE. IF (PHOCHA(IDPHO(I))0) THEN
1667  NCHARG=NCHARG+1
1668 .GT. IF (NCHARGNMXPHO) THEN
1669  DATA=NCHARG
1670  CALL PHOERR(1,'photos',DATA)
1671  ENDIF
1672  CHAPOI(NCHARG)=I
1673  ENDIF
1674  MINMAS=MINMAS+PPHO(5,I)**2
1675  ENDIF
1676  MASSUM=MASSUM+PPHO(5,I)
1677  20 CONTINUE
1678 .NE. IF (NCHARG0) THEN
1679 C--
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
1682 C--
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)
1686 C--
1687  30 CONTINUE
1688  DO 70 J=1,3
1689  70 PNEUTR(J)=-PPHO(J,CHAPOI(NCHARG))
1690  PNEUTR(4)=PPHO(5,IP)-PPHO(4,CHAPOI(NCHARG))
1691 C--
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
1696  NEUPOI=JDAPHO(1,IP)
1697 .EQ. IF (NEUPOICHAPOI(NCHARG)) NEUPOI=JDAPHO(2,IP)
1698  MNESQR=PPHO(5,NEUPOI)**2
1699  PNEUTR(5)=PPHO(5,NEUPOI)
1700  ELSE
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)
1704  ENDIF
1705 C--
1706 C-- Determine kinematical limit...
1707  XPHMAX=(MPASQR-(PNEUTR(5)+PPHO(5,CHAPOI(NCHARG)))**2)/MPASQR
1708 C--
1709 C-- Photon energy fraction...
1710  CALL PHOENE(MPASQR,MCHREN,BETA,BIGLOG,IDPHO(CHAPOI(NCHARG)))
1711 C--
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
1717 C--
1718 C-- No radiation was accepted, check for more daughters that may ra-
1719 C-- diate and correct radiation probability...
1720  NCHARG=NCHARG-1
1721 .GT. IF (NCHARG0) THEN
1722  IREP=IREP+1
1723  GOTO 30
1724  ENDIF
1725  ELSE
1726 C--
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)
1730 C--
1731 C-- Calculate sin(theta) and cos(theta) from interval variables
1732  DEL1=(2.D0-EPS)*(EPS/(2.D0-EPS))**PHORAN(THEDUM)
1733  DEL2=2.D0-DEL1
1734 
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
1741 C ELSE
1742 C COSTHG=-1D0+2*PHORAN(THEDUM)
1743 C SINTHG= SQRT(1D0-COSTHG**2)
1744 C ENDIF
1745 C
1746 .GT.C IF (FINT1.0D0) THEN
1747 C
1748 C WGT=1D0/(1D0-BETA*COSTHG)
1749 C WGT=WGT/(WGT+FINT)
1750 C ! WGT=1D0 ! ??
1751 C
1752 C ELSE
1753 C WGT=1D0
1754 C ENDIF
1755 C
1756 C ----------- END OF VARIANT B ------------------
1757 
1758 C ----------- VARIANT A ------------------
1759  COSTHG=(1.D0-DEL1)/BETA
1760  SINTHG=SQRT(DEL1*DEL2-MCHREN)/BETA
1761  WGT=1D0
1762 C ----------- END OF VARIANT A ------------------
1763 
1764 C--
1765 C-- Determine spin of particle and construct code for matrix element
1766  ME=2.D0*PHOSPI(IDPHO(CHAPOI(NCHARG)))+1.D0
1767 C--
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
1773  NEUDAU=I
1774  GOTO 51
1775  ENDIF
1776  ENDDO
1777 C--
1778 C-- Pointer not found...
1779  DATA=NCHARG
1780  CALL PHOERR(5,'phokin',DATA)
1781  51 CONTINUE
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
1786  ENDIF
1787  ELSE
1788  DATA=PPHO(5,IP)-MASSUM
1789  CALL PHOERR(10,'photos',DATA)
1790  ENDIF
1791  ENDIF
1792 C--
1793  RETURN
1794  END
1795  SUBROUTINE PHOOMA(IFIRST,ILAST,POINTR)
1796 C.----------------------------------------------------------------------
1797 C.
1798 C. PHOTOS: PHOton radiation in decays Order MAss vector
1799 C.
1800 C. Purpose: Order the contents of array 'pointr' according to the
1801 C. decreasing value in the array 'mass'.
1802 C.
1803 C. Input Parameters: IFIRST, ILAST: Pointers to the vector loca-
1804 C. tion be sorted,
1805 C. POINTR: Unsorted array with pointers to
1806 C. /PHOEVT/.
1807 C.
1808 C. Output Parameter: POINTR: Sorted arrays with respect to
1809 C. particle mass 'ppho(5,*)'.
1810 C.
1811 C. Author(s): B. van Eijk Created at: 28/11/89
1812 C. Last Update: 27/05/93
1813 C.
1814 C.----------------------------------------------------------------------
1815  IMPLICIT NONE
1816  INTEGER NMXPHO
1817  PARAMETER (NMXPHO=10000)
1818  INTEGER IDPHO,ISTPHO,JDAPHO,JMOPHO,NEVPHO,NPHO
1819  REAL*8 PPHO,VPHO
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
1825 C--
1826 C-- Copy particle masses
1827  DO 10 I=IFIRST,ILAST
1828  10 MASS(I)=PPHO(5,POINTR(I))
1829 C--
1830 C-- Order the masses in a decreasing series
1831  DO 30 I=IFIRST,ILAST-1
1832  DO 20 J=I+1,ILAST
1833 .LE. IF (MASS(J)MASS(I)) GOTO 20
1834  BUFPOI=POINTR(J)
1835  POINTR(J)=POINTR(I)
1836  POINTR(I)=BUFPOI
1837  BUFMAS=MASS(J)
1838  MASS(J)=MASS(I)
1839  MASS(I)=BUFMAS
1840  20 CONTINUE
1841  30 CONTINUE
1842  RETURN
1843  END
1844  SUBROUTINE PHOENE(MPASQR,MCHREN,BETA,BIGLOG,IDENT)
1845 C.----------------------------------------------------------------------
1846 C.
1847 C. PHOTOS: PHOton radiation in decays calculation of photon ENErgy
1848 C. fraction
1849 C.
1850 C. Purpose: Subroutine returns photon energy fraction (in (parent
1851 C. mass)/2 units) for the decay bremsstrahlung.
1852 C.
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.
1856 C.
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.
1861 C.
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
1864 C.
1865 C.----------------------------------------------------------------------
1866  IMPLICIT NONE
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
1871  REAL*8 PNEUTR
1872  INTEGER IDENT
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
1878  REAL*8 ALPHA,XPHCUT
1879  COMMON/PHOCOP/ALPHA,XPHCUT
1880  REAL*8 PI,TWOPI
1881  COMMON/PHPICO/PI,TWOPI
1882  INTEGER IREP
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
1888  INTEGER NX,NCHAN,K
1889  PARAMETER (NX=10)
1890  LOGICAL EXPINI
1891  REAL*8 PRO,PRSUM
1892  COMMON /PHOEXP/ PRO(NX),NCHAN,EXPINI
1893 C--
1894 .LE. IF (XPHMAXXPHCUT) THEN
1895  BETA=PHOFAC(-1) ! to zero counter, here beta is dummy
1896  XPHOTO=0.0D0
1897  RETURN
1898  ENDIF
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)
1902 
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 ------------------
1912 
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
1921  PRKILL=0d0
1922  IF (IEXP) THEN ! IEXP
1923  NCHAN=NCHAN+1
1924  IF (EXPINI) THEN ! EXPINI
1925  PRO(NCHAN)=PRHARD+0.05*(1.0+FINT) ! we store hard photon emission prob
1926  !for leg NCHAN
1927  PRHARD=0D0 ! to kill emission at initialization call
1928  PROBH=PRHARD
1929  ELSE ! EXPINI
1930  PRSUM=0
1931  DO K=NCHAN,NX
1932  PRSUM=PRSUM+PRO(K)
1933  ENDDO
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 !
1939 
1940  ENDIF ! EXPINI
1941  PRSOFT=1.D0-PRHARD
1942  ELSE ! IEXP
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.
1946  PROBH=PRHARD
1947  ENDIF ! IEXP
1948  PRSOFT=1.D0-PRHARD
1949 C--
1950 C-- Check on kinematical bounds
1951  IF (IEXP) THEN
1952 .LT. IF (PRSOFT-5.0D-8) THEN
1953  DATA=PRSOFT
1954  CALL PHOERR(2,'phoene',DATA)
1955  ENDIF
1956  ELSE
1957 .LT. IF (PRSOFT0.1D0) THEN
1958  DATA=PRSOFT
1959  CALL PHOERR(2,'phoene',DATA)
1960  ENDIF
1961  ENDIF
1962 
1963  RRR=PHORAN(IWT1)
1964 .LT. IF (RRRPRSOFT) THEN
1965 C--
1966 C-- No photon... (ie. photon too soft)
1967  XPHOTO=0.D0
1968 .LT. IF (RRRPRKILL) XPHOTO=-5d0 ! No photon...no further trials
1969  ELSE
1970 C--
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))
1976  & GOTO 10
1977  ENDIF
1978 C--
1979 C-- Calculate parameter for PHOFAC function
1980  XF=4.D0*MCHSQR*MPASQR/(MPASQR+MCHSQR-MNESQR)**2
1981  RETURN
1982  END
1983  FUNCTION PHOCOR(MPASQR,MCHREN,ME)
1984 C.----------------------------------------------------------------------
1985 C.
1986 C. PHOTOS: PHOton radiation in decays CORrection weight from
1987 C. matrix elements
1988 C.
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 !
1992 C.
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
1996 C.
1997 C. Output Parameter: Function value.
1998 C.
1999 C. Author(s): Z. Was, B. van Eijk Created at: 26/11/89
2000 C. Last Update: 21/03/93
2001 C.
2002 C.----------------------------------------------------------------------
2003  IMPLICIT NONE
2004  DOUBLE PRECISION MPASQR,MCHREN,BETA,XX,YY,DATA
2005  INTEGER ME
2006  REAL*8 PHOCOR,PHOFAC,WT1,WT2,WT3
2007  DOUBLE PRECISION MCHSQR,MNESQR
2008  REAL*8 PNEUTR
2009  COMMON/PHOMOM/MCHSQR,MNESQR,PNEUTR(5)
2010  DOUBLE PRECISION COSTHG,SINTHG
2011  REAL*8 XPHMAX,XPHOTO
2012  COMMON/PHOPHS/XPHMAX,XPHOTO,COSTHG,SINTHG
2013  INTEGER IREP
2014  REAL*8 PROBH,CORWT,XF
2015  COMMON/PHOPRO/PROBH,CORWT,XF,IREP
2016 C--
2017 C-- Shaping (modified by ZW)...
2018  XX=4.D0*MCHSQR/MPASQR*(1.D0-XPHOTO)/(1.D0-XPHOTO+(MCHSQR-MNESQR)/
2019  &MPASQR)**2
2020 .EQ. IF (ME1) THEN
2021  YY=1.D0
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))
2025  WT3=1.D0
2026 .EQ..OR..EQ..OR..EQ. ELSEIF ((ME3)(ME4)(ME5)) THEN
2027  YY=1.D0
2028  WT3=(1.D0+(1.D0-XPHOTO/XPHMAX)**2-(XPHOTO/XPHMAX)**3)/
2029  & (1.D0+(1.D0-XPHOTO/XPHMAX)** 2)
2030  ELSE
2031  DATA=(ME-1.D0)/2.D0
2032  CALL PHOERR(6,'phocor',DATA)
2033  YY=1.D0
2034  WT3=1.D0
2035  ENDIF
2036  BETA=SQRT(1.D0-XX)
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
2039  WT2=WT2*PHOFAC(1)
2040  PHOCOR=WT1*WT2*WT3
2041  CORWT=PHOCOR
2042 .GT. IF (PHOCOR1.D0) THEN
2043  DATA=PHOCOR
2044  CALL PHOERR(3,'phocor',DATA)
2045  ENDIF
2046  RETURN
2047  END
2048  FUNCTION PHOFAC(MODE)
2049 C.----------------------------------------------------------------------
2050 C.
2051 C. PHOTOS: PHOton radiation in decays control FACtor
2052 C.
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.
2069 C.
2070 C. Input Parameters: MODE, PROBH, XF
2071 C.
2072 C. Output Parameter: Function value
2073 C.
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
2076 C.
2077 C.----------------------------------------------------------------------
2078  IMPLICIT NONE
2079  REAL*8 PHOFAC,FF,PRX
2080  INTEGER MODE
2081  INTEGER IREP
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
2087  SAVE PRX,FF
2088  DATA PRX,FF/ 0.D0, 0.D0/
2089  IF (IEXP) THEN ! In case of exponentiation this routine is useles
2090  PHOFAC=1
2091  RETURN
2092  ENDIF
2093 .EQ. IF (MODE-1) THEN
2094  PRX=1.D0
2095  FF=1.D0
2096  PROBH=0.0
2097 .EQ. ELSEIF (MODE0) THEN
2098 .EQ. IF (IREP0) PRX=1.D0
2099  PRX=PRX/(1.D0-PROBH)
2100  FF=1.D0
2101 C--
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
2107  PHOFAC=FF*PRX
2108  ELSE
2109  PHOFAC=1.D0/FF
2110  ENDIF
2111  END
2112  SUBROUTINE PHOBW(WT)
2113 C.----------------------------------------------------------------------
2114 C.
2115 C. PHOTOS: PHOtos Boson W correction weight
2116 C.
2117 C. Purpose: calculates correction weight due to amplitudes of
2118 C. emission from W boson.
2119 C.
2120 C.
2121 C.
2122 C.
2123 C.
2124 C. Input Parameters: Common /PHOEVT/, with photon added.
2125 C. wt to be corrected
2126 C.
2127 C.
2128 C.
2129 C. Output Parameters: wt
2130 C.
2131 C. Author(s): G. Nanava, Z. Was Created at: 13/03/03
2132 C. Last Update: 13/03/03
2133 C.
2134 C.----------------------------------------------------------------------
2135  IMPLICIT NONE
2136  DOUBLE PRECISION WT
2137  INTEGER NMXPHO
2138  PARAMETER (NMXPHO=10000)
2139  INTEGER IDPHO,ISTPHO,JDAPHO,JMOPHO,NEVPHO,NPHO
2140  REAL*8 PPHO,VPHO
2141  COMMON/PHOEVT/NEVPHO,NPHO,ISTPHO(NMXPHO),IDPHO(NMXPHO),
2142  &JMOPHO(2,NMXPHO),JDAPHO(2,NMXPHO),PPHO(5,NMXPHO),VPHO(4,NMXPHO)
2143  INTEGER I
2144  DOUBLE PRECISION EMU,MCHREN,BETA,COSTHG,MPASQR,XPH
2145 C--
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
2151 
2152  IF(
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
2156  I=JDAPHO(1,1)
2157  ELSE
2158  I=JDAPHO(1,1)+1
2159  ENDIF
2160  EMU=PPHO(4,I)
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)
2168  MPASQR=PPHO(4,1)**2
2169  XPH=PPHO(4,NPHO)
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))
2173  ENDIF
2174 c write(*,*) IDPHO(1),IDPHO(JDAPHO(1,1)),IDPHO(JDAPHO(1,1)+1)
2175 c write(*,*) emu,xph,costhg,beta,mpasqr,mchren
2176 
2177  END
2178  SUBROUTINE PHODO(IP,NCHARB,NEUDAU)
2179 C.----------------------------------------------------------------------
2180 C.
2181 C. PHOTOS: PHOton radiation in decays DOing of KINematics
2182 C.
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.
2189 C.
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.
2196 C.
2197 C. Author(s): Z. Was, B. van Eijk Created at: 26/11/89
2198 C. Last Update: 27/05/93
2199 C.
2200 C.----------------------------------------------------------------------
2201  IMPLICIT NONE
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
2205  INTEGER NCHARB
2206  REAL*8 EPHOTO,PMAVIR,PHOTRI
2207  REAL*8 GNEUT,PHORAN,CCOSTH,SSINTH,PVEC(4)
2208  INTEGER NMXPHO
2209  PARAMETER (NMXPHO=10000)
2210  INTEGER IDPHO,ISTPHO,JDAPHO,JMOPHO,NEVPHO,NPHO
2211  REAL*8 PPHO,VPHO
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
2215  REAL*8 PNEUTR
2216  COMMON/PHOMOM/MCHSQR,MNESQR,PNEUTR(5)
2217  DOUBLE PRECISION COSTHG,SINTHG
2218  REAL*8 XPHMAX,XPHOTO
2219  COMMON/PHOPHS/XPHMAX,XPHOTO,COSTHG,SINTHG
2220  REAL*8 PI,TWOPI
2221  COMMON/PHPICO/PI,TWOPI
2222 C--
2223  EPHOTO=XPHOTO*PPHO(5,IP)/2.D0
2224  PMAVIR=SQRT(PPHO(5,IP)*(PPHO(5,IP)-2.D0*EPHOTO))
2225 C--
2226 C-- Reconstruct kinematics of charged particle and neutral system
2227  FI1=PHOAN1(PNEUTR(1),PNEUTR(2))
2228 C--
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))
2234 C--
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))
2242  QOLD=PNEUTR(3)
2243  GNEUT=(QNEW**2+QOLD**2+MNESQR)/(QNEW*QOLD+SQRT((QNEW**2+MNESQR)*
2244  &(QOLD**2+MNESQR)))
2245 .LT. IF (GNEUT1.D0) THEN
2246  DATA=0.D0
2247  CALL PHOERR(4,'phokin',DATA)
2248  ENDIF
2249  PARNE=GNEUT-SQRT(MAX(GNEUT**2-1.0D0,0.D0))
2250 C--
2251 C-- 2) ...reductive boost...
2252  CALL PHOBO3(PARNE,PNEUTR)
2253 C--
2254 C-- ...calculate photon energy in the reduced system...
2255  NPHO=NPHO+1
2256  ISTPHO(NPHO)=1
2257  IDPHO(NPHO) =22
2258 C-- Photon mother and daughter pointers !
2259  JMOPHO(1,NPHO)=IP
2260  JMOPHO(2,NPHO)=0
2261  JDAPHO(1,NPHO)=0
2262  JDAPHO(2,NPHO)=0
2263  PPHO(4,NPHO)=EPHOTO*PPHO(5,IP)/PMAVIR
2264 C--
2265 C-- ...and photon momenta
2266  CCOSTH=-COSTHG
2267  SSINTH=SINTHG
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)
2272 C--
2273 C-- Minus sign because axis opposite direction of charged particle !
2274  PPHO(3,NPHO)=-PPHO(4,NPHO)*COSTHG
2275  PPHO(5,NPHO)=0.D0
2276 C--
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)
2283 C--
2284 C-- Boost to the rest frame of decaying particle
2285  CALL PHOBO3(ANGLE,PNEUTR(1))
2286  CALL PHOBO3(ANGLE,PPHO(1,NPHO))
2287 C--
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))
2293 C--
2294  DO 60 I=2,4
2295  60 PVEC(I)=0.D0
2296  PVEC(1)=1.D0
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))
2305 C--
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
2316 C--
2317  FIRST=NEUDAU
2318  LAST=JDAPHO(2,IP)
2319  DO 70 I=FIRST,LAST
2320 .NE..AND..EQ. IF (INCHARB(JMOPHO(1,I)IP)) THEN
2321 C--
2322 C-- Reconstruct kinematics...
2323  CALL PHORO3(-FI1,PPHO(1,I))
2324  CALL PHORO2(-TH1,PPHO(1,I))
2325 C--
2326 C-- ...reductive boost
2327  CALL PHOBO3(PARNE,PPHO(1,I))
2328 C--
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))
2332 C--
2333 C-- Boost to the rest frame of decaying particle
2334  CALL PHOBO3(ANGLE,PPHO(1,I))
2335 C--
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))
2339 C--
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))
2344  ENDIF
2345  70 CONTINUE
2346  ELSE
2347 C--
2348 C-- ...only one 'neutral' particle in addition to photon!
2349  DO 80 J=1,4
2350  80 PPHO(J,NEUDAU)=PNEUTR(J)
2351  ENDIF
2352 C--
2353 C-- All 'neutrals' treated, fill /PHOEVT/ for charged particle...
2354  DO 90 J=1,3
2355  90 PPHO(J,NCHARB)=-(PPHO(J,NPHO)+PNEUTR(J))
2356  PPHO(4,NCHARB)=PPHO(5,IP)-(PPHO(4,NPHO)+PNEUTR(4))
2357 C--
2358  END
2359  FUNCTION PHOTRI(A,B,C)
2360 C.----------------------------------------------------------------------
2361 C.
2362 C. PHOTOS: PHOton radiation in decays calculation of TRIangle fie
2363 C.
2364 C. Purpose: Calculation of triangle function for phase space.
2365 C.
2366 C. Input Parameters: A, B, C (Virtual) particle masses.
2367 C.
2368 C. Output Parameter: Function value =
2369 C. SQRT(LAMBDA(A**2,B**2,C**2))/(2*A)
2370 C.
2371 C. Author(s): B. van Eijk Created at: 15/11/89
2372 C. Last Update: 02/01/90
2373 C.
2374 C.----------------------------------------------------------------------
2375  IMPLICIT NONE
2376  DOUBLE PRECISION DA,DB,DC,DAPB,DAMB,DTRIAN
2377  REAL*8 A,B,C,PHOTRI
2378  DA=A
2379  DB=B
2380  DC=C
2381  DAPB=DA+DB
2382  DAMB=DA-DB
2383  DTRIAN=SQRT((DAMB-DC)*(DAPB+DC)*(DAMB+DC)*(DAPB-DC))
2384  PHOTRI=DTRIAN/(DA+DA)
2385  RETURN
2386  END
2387  FUNCTION PHOAN1(X,Y)
2388 C.----------------------------------------------------------------------
2389 C.
2390 C. PHOTOS: PHOton radiation in decays calculation of ANgle '1'
2391 C.
2392 C. Purpose: Calculate angle from X and Y
2393 C.
2394 C. Input Parameters: X, Y
2395 C.
2396 C. Output Parameter: Function value
2397 C.
2398 C. Author(s): S. Jadach Created at: 01/01/89
2399 C. B. van Eijk Last Update: 02/01/90
2400 C.
2401 C.----------------------------------------------------------------------
2402  IMPLICIT NONE
2403  DOUBLE PRECISION PHOAN1
2404  REAL*8 X,Y
2405  REAL*8 PI,TWOPI
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
2410  ELSE
2411  PHOAN1=ACOS(X/SQRT(X**2+Y**2))
2412  ENDIF
2413 .LT. IF (Y0.D0) PHOAN1=TWOPI-PHOAN1
2414  RETURN
2415  END
2416  FUNCTION PHOAN2(X,Y)
2417 C.----------------------------------------------------------------------
2418 C.
2419 C. PHOTOS: PHOton radiation in decays calculation of ANgle '2'
2420 C.
2421 C. Purpose: Calculate angle from X and Y
2422 C.
2423 C. Input Parameters: X, Y
2424 C.
2425 C. Output Parameter: Function value
2426 C.
2427 C. Author(s): S. Jadach Created at: 01/01/89
2428 C. B. van Eijk Last Update: 02/01/90
2429 C.
2430 C.----------------------------------------------------------------------
2431  IMPLICIT NONE
2432  DOUBLE PRECISION PHOAN2
2433  REAL*8 X,Y
2434  REAL*8 PI,TWOPI
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
2439  ELSE
2440  PHOAN2=ACOS(X/SQRT(X**2+Y**2))
2441  ENDIF
2442  RETURN
2443  END
2444  SUBROUTINE PHOBO3(ANGLE,PVEC)
2445 C.----------------------------------------------------------------------
2446 C.
2447 C. PHOTOS: PHOton radiation in decays BOost routine '3'
2448 C.
2449 C. Purpose: Boost vector PVEC along z-axis where ANGLE = EXP(ETA),
2450 C. ETA is the hyperbolic velocity.
2451 C.
2452 C. Input Parameters: ANGLE, PVEC
2453 C.
2454 C. Output Parameter: PVEC
2455 C.
2456 C. Author(s): S. Jadach Created at: 01/01/89
2457 C. B. van Eijk Last Update: 02/01/90
2458 C.
2459 C.----------------------------------------------------------------------
2460  IMPLICIT NONE
2461  DOUBLE PRECISION QPL,QMI,ANGLE
2462  REAL*8 PVEC(4)
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
2467  RETURN
2468  END
2469  SUBROUTINE PHORO2(ANGLE,PVEC)
2470 C.----------------------------------------------------------------------
2471 C.
2472 C. PHOTOS: PHOton radiation in decays ROtation routine '2'
2473 C.
2474 C. Purpose: Rotate x and z components of vector PVEC around angle
2475 C. 'angle'.
2476 C.
2477 C. Input Parameters: ANGLE, PVEC
2478 C.
2479 C. Output Parameter: PVEC
2480 C.
2481 C. Author(s): S. Jadach Created at: 01/01/89
2482 C. B. van Eijk Last Update: 02/01/90
2483 C.
2484 C.----------------------------------------------------------------------
2485  IMPLICIT NONE
2486  DOUBLE PRECISION CS,SN,ANGLE
2487  REAL*8 PVEC(4)
2488  CS=COS(ANGLE)*PVEC(1)+SIN(ANGLE)*PVEC(3)
2489  SN=-SIN(ANGLE)*PVEC(1)+COS(ANGLE)*PVEC(3)
2490  PVEC(1)=CS
2491  PVEC(3)=SN
2492  RETURN
2493  END
2494  SUBROUTINE PHORO3(ANGLE,PVEC)
2495 C.----------------------------------------------------------------------
2496 C.
2497 C. PHOTOS: PHOton radiation in decays ROtation routine '3'
2498 C.
2499 C. Purpose: Rotate x and y components of vector PVEC around angle
2500 C. 'angle'.
2501 C.
2502 C. Input Parameters: ANGLE, PVEC
2503 C.
2504 C. Output Parameter: PVEC
2505 C.
2506 C. Author(s): S. Jadach Created at: 01/01/89
2507 C. B. van Eijk Last Update: 02/01/90
2508 C.
2509 C.----------------------------------------------------------------------
2510  IMPLICIT NONE
2511  DOUBLE PRECISION CS,SN,ANGLE
2512  REAL*8 PVEC(4)
2513  CS=COS(ANGLE)*PVEC(1)-SIN(ANGLE)*PVEC(2)
2514  SN=SIN(ANGLE)*PVEC(1)+COS(ANGLE)*PVEC(2)
2515  PVEC(1)=CS
2516  PVEC(2)=SN
2517  RETURN
2518  END
2519  SUBROUTINE PHORIN
2520 C.----------------------------------------------------------------------
2521 C.
2522 C. PHOTOS: PHOton radiation in decays RANdom number generator init
2523 C.
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.
2527 C.
2528 C. Input Parameters: ISEED(*)
2529 C.
2530 C. Output Parameters: URAN, CRAN, CDRAN, CMRAN, I97, J97
2531 C.
2532 C. Author(s): B. van Eijk and F. James Created at: 27/09/89
2533 C. Last Update: 22/02/90
2534 C.
2535 C.----------------------------------------------------------------------
2536  IMPLICIT NONE
2537  DOUBLE PRECISION DATA
2538  REAL*8 S,T
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
2543 C--
2544 C-- Check value range of seeds
2545 .LT..OR..GE. IF ((ISEED(1)0)(ISEED(1)31328)) THEN
2546  DATA=ISEED(1)
2547  CALL PHOERR(8,'phorin',DATA)
2548  ENDIF
2549 .LT..OR..GE. IF ((ISEED(2)0)(ISEED(2)30081)) THEN
2550  DATA=ISEED(2)
2551  CALL PHOERR(9,'phorin',DATA)
2552  ENDIF
2553 C--
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)
2559  DO 20 I=1,97
2560  S=0.D0
2561  T=0.5D0
2562  DO 10 J=1,24
2563  IS5=MOD (MOD(IS1*IS2,179)*IS3,179)
2564  IS1=IS2
2565  IS2=IS3
2566  IS3=IS5
2567  IS4=MOD(53*IS4+1,169)
2568 .GE. IF (MOD(IS4*IS5,64)32) S=S+T
2569  10 T=0.5D0*T
2570  20 URAN(I)=S
2571  CRAN=362436.D0/16777216.D0
2572  CDRAN=7654321.D0/16777216.D0
2573  CMRAN=16777213.D0/16777216.D0
2574  I97=97
2575  J97=33
2576  RETURN
2577  END
2578  FUNCTION PHORAN(IDUM)
2579 C.----------------------------------------------------------------------
2580 C.
2581 C. PHOTOS: PHOton radiation in decays RANdom number generator based
2582 C. on Marsaglia Algorithm
2583 C.
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.
2590 C.
2591 C. Input Parameters: IDUM (integer dummy)
2592 C.
2593 C. Output Parameters: Function value
2594 C.
2595 C. Author(s): B. van Eijk, G. Marsaglia and Created at: 27/09/89
2596 C. A. Zaman Last Update: 27/09/89
2597 C.
2598 C.----------------------------------------------------------------------
2599  IMPLICIT NONE
2600  REAL*8 PHORAN
2601  INTEGER IDUM
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
2607  URAN(I97)=PHORAN
2608  I97=I97-1
2609 .EQ. IF (I970) I97=97
2610  J97=J97-1
2611 .EQ. IF (J970) J97=97
2612  CRAN=CRAN-CDRAN
2613 .LT. IF (CRAN0.D0) CRAN=CRAN+CMRAN
2614  PHORAN=PHORAN-CRAN
2615 .LT. IF (PHORAN0.D0) PHORAN=PHORAN+1.D0
2616 .LE. IF (PHORAN0.D0) GOTO 10
2617  RETURN
2618  END
2619  FUNCTION PHOCHA(IDHEP)
2620 C.----------------------------------------------------------------------
2621 C.
2622 C. PHOTOS: PHOton radiation in decays CHArge determination
2623 C.
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.
2627 C.
2628 C. Input Parameter: IDHEP
2629 C.
2630 C. Output Parameter: Funtion value = charge of particle with code
2631 C. IDHEP
2632 C.
2633 C. Author(s): E. Barberio and B. van Eijk Created at: 29/11/89
2634 C. Last update: 02/01/90
2635 C.
2636 C.----------------------------------------------------------------------
2637  IMPLICIT NONE
2638  REAL*8 PHOCHA
2639  INTEGER IDHEP,IDABS,Q1,Q2,Q3
2640 C--
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)
2644  DATA CHARGE/ 0.D0,
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/
2649  IDABS=ABS(IDHEP)
2650 .LE. IF (IDABS100) THEN
2651 C--
2652 C-- Charge of quark, lepton, boson etc....
2653  PHOCHA = CHARGE(IDABS)
2654  ELSE
2655 C--
2656 C-- Check on particle build out of quarks, unpack its code...
2657  Q3=MOD(IDABS/1000,10)
2658  Q2=MOD(IDABS/100,10)
2659  Q1=MOD(IDABS/10,10)
2660 .EQ. IF (Q30) THEN
2661 C--
2662 C-- ...meson...
2663 .EQ. IF(MOD(Q2,2)0) THEN
2664  PHOCHA=CHARGE(Q2)-CHARGE(Q1)
2665  ELSE
2666  PHOCHA=CHARGE(Q1)-CHARGE(Q2)
2667  ENDIF
2668  ELSE
2669 C--
2670 C-- ...diquarks or baryon.
2671  PHOCHA=CHARGE(Q1)+CHARGE(Q2)+CHARGE(Q3)
2672  ENDIF
2673  ENDIF
2674 C--
2675 C-- Find the sign of the charge...
2676 .LT. IF (IDHEP0.D0) PHOCHA=-PHOCHA
2677 .lt. IF (PHOCHA**21d-6) PHOCHA=0.D0
2678  RETURN
2679  END
2680  FUNCTION PHOSPI(IDHEP)
2681 C.----------------------------------------------------------------------
2682 C.
2683 C. PHOTOS: PHOton radiation in decays function for SPIn determina-
2684 C. tion
2685 C.
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.
2689 C.
2690 C. Input Parameter: IDHEP
2691 C.
2692 C. Output Parameter: Funtion value = spin of particle with code
2693 C. IDHEP
2694 C.
2695 C. Author(s): E. Barberio and B. van Eijk Created at: 29/11/89
2696 C. Last update: 02/01/90
2697 C.
2698 C.----------------------------------------------------------------------
2699  IMPLICIT NONE
2700  REAL*8 PHOSPI
2701  INTEGER IDHEP,IDABS
2702 C--
2703 C-- Array 'spin' contains the spin of the first 100 particles accor-
2704 C-- ding to the PDG particle code...
2705  REAL*8 SPIN(100)
2706  DATA SPIN/ 8*.5D0, 1.D0, 0.D0, 8*.5D0, 2*0.D0, 4*1.D0, 76*0.D0/
2707  IDABS=ABS(IDHEP)
2708 C--
2709 C-- Spin of quark, lepton, boson etc....
2710 .LE. IF (IDABS100) THEN
2711  PHOSPI=SPIN(IDABS)
2712  ELSE
2713 C--
2714 C-- ...other particles, however...
2715  PHOSPI=(MOD(IDABS,10)-1.D0)/2.D0
2716 C--
2717 C-- ...K_short and K_long are special !!
2718  PHOSPI=MAX(PHOSPI,0.D0)
2719  ENDIF
2720  RETURN
2721  END
2722  SUBROUTINE PHOERR(IMES,TEXT,DATA)
2723 C.----------------------------------------------------------------------
2724 C.
2725 C. PHOTOS: PHOton radiation in decays ERRror handling
2726 C.
2727 C. Purpose: Inform user about (fatal) errors and warnings generated
2728 C. by either the user or the program.
2729 C.
2730 C. Input Parameters: IMES, TEXT, DATA
2731 C.
2732 C. Output Parameters: None
2733 C.
2734 C. Author(s): B. van Eijk Created at: 29/11/89
2735 C. Last Update: 10/01/92
2736 C.
2737 C.----------------------------------------------------------------------
2738  IMPLICIT NONE
2739  DOUBLE PRECISION DATA
2740  INTEGER IMES,IERROR
2741  REAL*8 SDATA
2742  INTEGER PHLUN
2743  COMMON/PHOLUN/PHLUN
2744  INTEGER PHOMES
2745  PARAMETER (PHOMES=10)
2746  INTEGER STATUS
2747  COMMON/PHOSTA/STATUS(PHOMES)
2748  CHARACTER TEXT*(*)
2749  SAVE IERROR
2750 C-- security STOP switch
2751  LOGICAL ISEC
2752  SAVE ISEC
2753  DATA ISEC /.TRUE./
2754  DATA IERROR/ 0/
2755 .LE. IF (IMESPHOMES) STATUS(IMES)=STATUS(IMES)+1
2756 C--
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
2760  SDATA=DATA
2761  WRITE(PHLUN,9000)
2762  WRITE(PHLUN,9120)
2763  GOTO (10,20,30,40,50,60,70,80,90,100),IMES
2764  WRITE(PHLUN,9130) IMES
2765  GOTO 120
2766  10 WRITE(PHLUN,9010) TEXT,INT(SDATA)
2767  GOTO 110
2768  20 WRITE(PHLUN,9020) TEXT,SDATA
2769  GOTO 110
2770  30 WRITE(PHLUN,9030) TEXT,SDATA
2771  GOTO 110
2772  40 WRITE(PHLUN,9040) TEXT
2773  GOTO 110
2774  50 WRITE(PHLUN,9050) TEXT,INT(SDATA)
2775  GOTO 110
2776  60 WRITE(PHLUN,9060) TEXT,SDATA
2777  GOTO 130
2778  70 WRITE(PHLUN,9070) TEXT,INT(SDATA)
2779  GOTO 110
2780  80 WRITE(PHLUN,9080) TEXT,INT(SDATA)
2781  GOTO 110
2782  90 WRITE(PHLUN,9090) TEXT,INT(SDATA)
2783  GOTO 110
2784  100 WRITE(PHLUN,9100) TEXT,SDATA
2785  GOTO 130
2786  110 CONTINUE
2787  WRITE(PHLUN,9140)
2788  WRITE(PHLUN,9120)
2789  WRITE(PHLUN,9000)
2790  IF (ISEC) THEN
2791  STOP
2792  ELSE
2793  GOTO 130
2794  ENDIF
2795  120 IERROR=IERROR+1
2796 .GE. IF (IERROR10) THEN
2797  WRITE(PHLUN,9150)
2798  WRITE(PHLUN,9120)
2799  WRITE(PHLUN,9000)
2800  IF (ISEC) THEN
2801  STOP
2802  ELSE
2803  GOTO 130
2804  ENDIF
2805  ENDIF
2806  130 WRITE(PHLUN,9120)
2807  WRITE(PHLUN,9000)
2808  RETURN
2809  9000 FORMAT(1H ,80('*'))
2810  9010 FORMAT(1H ,'* ',A,': too many charged particles, ncharg =',I6,T81,
2811  &'*')
2812  9020 FORMAT(1H ,'* ',A,': too much bremsstrahlung required, prsoft = ',
2813  &F15.6,T81,'*')
2814  9030 FORMAT(1H ,'* ',A,': combined weight is exceeding 1., weight = ',
2815  &F15.6,T81,'*')
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,
2822  &'*')
2823  9070 FORMAT(1H ,'* ',A, ': stack length exceeded, nstack = ',I5 ,T81,
2824  &'*')
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,
2831  &'*')
2832  9120 FORMAT(1H ,'*',T81,'*')
2833  9130 FORMAT(1H ,'* funny error message: ',I4,' ! What to do ?',T81,'*')
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,
2836  &'*')
2837  END
2838  SUBROUTINE phorep
2839 c.----------------------------------------------------------------------
2840 c.
2841 c. photos: photon radiation in decays run summary report
2842 c.
2843 c. purpose: inform user about success and/or restrictions of photos
2844 c. encountered during execution.
2845 c.
2846 c. input parameters: Common /phosta/
2847 c.
2848 c. output parameters: none
2849 c.
2850 c. author(s): b. van eijk created at: 10/01/92
2851 c. last update: 10/01/92
2852 c.
2853 c.----------------------------------------------------------------------
2854  IMPLICIT NONE
2855  INTEGER PHLUN
2856  common/pholun/phlun
2857  INTEGER PHOMES
2858  parameter(phomes=10)
2859  INTEGER STATUS
2860  common/phosta/status(phomes)
2861  INTEGER I
2862  LOGICAL ERROR
2863  error=.false.
2864  WRITE(phlun,9000)
2865  WRITE(phlun,9010)
2866  WRITE(phlun,9020)
2867  WRITE(phlun,9030)
2868  WRITE(phlun,9040)
2869  WRITE(phlun,9030)
2870  WRITE(phlun,9020)
2871  DO 10 i=1,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)
2875  ELSE
2876  error=.true.
2877  WRITE(phlun,9060) i,status(i)
2878  ENDIF
2879  10 CONTINUE
2880  IF (.NOT.error) WRITE(phlun,9070)
2881  WRITE(phlun,9020)
2882  WRITE(phlun,9010)
2883  RETURN
2884  9000 FORMAT(1h1)
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',
2892  &t81,'*')
2893  END
2894  SUBROUTINE phlupa(IPOINT)
2895  IMPLICIT NONE
2896 c.----------------------------------------------------------------------
2897 c.
2898 c. phlupa: debugging tool
2899 c.
2900 c. purpose: none, eventually may printout content of the
2901 c. /phoevt/ common
2902 c.
2903 c. input parameters: Common /phoevt/ and /phnum/
2904 c. latter may have number of the event.
2905 c.
2906 c. output parameters: none
2907 c.
2908 c. author(s): z. was created at: 30/05/93
2909 c. last update: 09/10/05
2910 c.
2911 c.----------------------------------------------------------------------
2912  INTEGER NMXPHO
2913  parameter(nmxpho=10000)
2914  INTEGER IDPHO,ISTPHO,JDAPHO,JMOPHO,NEVPHO,NPHO,I,J,IPOINT
2915  INTEGER IPOIN,IPOIN0,IPOINM,IEV
2916  INTEGER IOUT
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)
2920  COMMON /phnum/ iev
2921  INTEGER PHLUN
2922  common/pholun/phlun
2923  dimension sum(5)
2924  DATA ipoin0/ -5/
2925  COMMON /phlupy/ ipoin,ipoinm
2926  SAVE ipoin0
2927  IF (ipoin0.LT.0) THEN
2928  ipoin0=300 000 ! maximal no-print point
2929  ipoin =ipoin0
2930  ipoinm=300 001 ! minimal no-print point
2931  ENDIF
2932  IF (ipoint.LE.ipoinm.OR.ipoint.GE.ipoin ) RETURN
2933  iout=56
2934  IF (iev.LT.1000) THEN
2935  DO i=1,5
2936  sum(i)=0.0d0
2937  ENDDO
2938  WRITE(phlun,*) 'EVENT NR=',iev,
2939  $ 'WE ARE TESTING /PHOEVT/ at IPOINT=',ipoint
2940  WRITE(phlun,10)
2941  i=1
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)
2944  i=2
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)
2947  WRITE(phlun,*) ' '
2948  DO i=3,npho
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)
2951  DO j=1,4
2952  sum(j)=sum(j)+ppho(j,i)
2953  ENDDO
2954  ENDDO
2955  sum(5)=sqrt(abs(sum(4)**2-sum(1)**2-sum(2)**2-sum(3)**2))
2956  WRITE(phlun,30) sum
2957  10 FORMAT(1x,' ID ','p_x ','p_y ','p_z ',
2958  $ 'E ','m ',
2959  $ 'ID-MO_DA1','ID-MO DA2' )
2960  20 FORMAT(1x,i4,5(f14.9),2i9)
2961  30 FORMAT(1x,' SUM',5(f14.9))
2962  ENDIF
2963  END
2964 
2965 
2966 
2967  FUNCTION iphqrk(MODCOR)
2968  implicit none
2969 c.----------------------------------------------------------------------
2970 c.
2971 c. iphqrk: enables blocks emision from quarks
2972 c.
2973 c
2974 c. input parameters: modcor
2975 c. modcor >0 type of action
2976 c. =1 blocks
2977 c. =2 enables
2978 c. =0 execution mode(retrns stored value)
2979 c.
2980 c.
2981 c. author(s): z. was created at: 11/12/00
2982 c. modified :
2983 c.----------------------------------------------------------------------
2984  INTEGER IPHQRK,MODCOR,MODOP
2985  INTEGER PHLUN
2986  common/pholun/phlun
2987 
2988  SAVE modop
2989  DATA modop /0/
2990  IF (modcor.NE.0) THEN
2991 c initialization
2992  modop=modcor
2993 
2994  WRITE(phlun,*)
2995  $ 'Message from PHOTOS: IPHQRK(MODCOR):: (re)initialization'
2996  IF (modop.EQ.1) THEN
2997  WRITE(phlun,*)
2998  $ 'MODOP=1 -- blocks emission from light quarks: DEFAULT'
2999  ELSEIF (modop.EQ.2) THEN
3000  WRITE(phlun,*)
3001  $ 'MODOP=2 -- enables emission from light quarks: TEST '
3002  ELSE
3003  WRITE(phlun,*) 'IPHQRK wrong MODCOR=',modcor
3004  stop
3005  ENDIF
3006  RETURN
3007  ENDIF
3008 
3009  IF (modop.EQ.0.AND.modcor.EQ.0) THEN
3010  WRITE(phlun,*) 'IPHQRK lack of initialization'
3011  stop
3012  ENDIF
3013  iphqrk=modop
3014  END
3015 
3016 
3017  FUNCTION iphekl(MODCOR)
3018  implicit none
3019 c.----------------------------------------------------------------------
3020 c.
3021 c. iphekl: enables/blocks emision in: pi0 to gamma e+ e-
3022 c.
3023 c
3024 c. input parameters: modcor
3025 c. modcor >0 type of action
3026 c. =1 blocks
3027 c. =2 enables
3028 c. =0 execution mode(retrns stored value)
3029 c.
3030 c.
3031 c. author(s): z. was created at: 11/12/00
3032 c. modified :
3033 c.----------------------------------------------------------------------
3034  INTEGER IPHEKL,MODCOR,MODOP
3035  INTEGER PHLUN
3036  common/pholun/phlun
3037 
3038  SAVE modop
3039  DATA modop /0/
3040 
3041  IF (modcor.NE.0) THEN
3042 c initialization
3043  modop=modcor
3044 
3045  WRITE(phlun,*)
3046  $ 'Message from PHOTOS: IPHEKL(MODCOR):: (re)initialization'
3047  IF (modop.EQ.2) THEN
3048  WRITE(phlun,*)
3049  $ 'MODOP=2 -- blocks emission in pi0 to gamma e+e-: DEFAULT'
3050  WRITE(phlun,*)
3051  $ 'MODOP=2 -- blocks emission in Kl to gamma e+e-: DEFAULT'
3052  ELSEIF (modop.EQ.1) THEN
3053  WRITE(phlun,*)
3054  $ 'MODOP=1 -- enables emission in pi0 to gamma e+e- : TEST '
3055  WRITE(phlun,*)
3056  $ 'MODOP=1 -- enables emission in Kl to gamma e+e- : TEST '
3057  ELSE
3058  WRITE(phlun,*) 'IPHEKL wrong MODCOR=',modcor
3059  stop
3060  ENDIF
3061  RETURN
3062  ENDIF
3063 
3064  IF (modop.EQ.0.AND.modcor.EQ.0) THEN
3065  WRITE(phlun,*) 'IPHELK lack of initialization'
3066  stop
3067  ENDIF
3068  iphekl=modop
3069  END
3070 
3071  SUBROUTINE phcork(MODCOR)
3072  implicit none
3073 c.----------------------------------------------------------------------
3074 c.
3075 c. phcork: corrects kinmatics of subbranch needed if host program
3076 c. produces events with the shaky momentum conservation
3077 c
3078 c. input parameters: Common /phoevt/, modcor
3079 c. modcor >0 type of action
3080 c. =1 no 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
3089 c.
3090 c. output parameters: corrected /phoevt/
3091 c.
3092 c. author(s): p.golonka, z. was created at: 01/02/99
3093 c. modified : 08/02/99
3094 c.----------------------------------------------------------------------
3095  INTEGER NMXPHO
3096  parameter(nmxpho=10000)
3097 
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
3101  real*8 ppho,vpho
3102  common/phoevt/nevpho,npho,istpho(nmxpho),idpho(nmxpho),
3103  &jmopho(2,nmxpho),jdapho(2,nmxpho),ppho(5,nmxpho),vpho(4,nmxpho)
3104 
3105  INTEGER PHLUN
3106  common/pholun/phlun
3107 
3108  COMMON /phnum/ iev
3109  SAVE modop
3110  DATA modop /0/
3111  SAVE iprint
3112  DATA iprint /0/
3113  SAVE mcut
3114  IF (modcor.NE.0) THEN
3115 c initialization
3116  modop=modcor
3117 
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 '
3128  mcut=0.4
3129  WRITE(phlun,*) 'Mcut=',mcut,'GeV'
3130  ELSEIF (modop.EQ.5) THEN
3131  WRITE(phlun,*) 'MODOP=5 -- corrects Energy from mass+flow'
3132 
3133  ELSE
3134  WRITE(phlun,*) 'PHCORK wrong MODCOR=',modcor
3135  stop
3136  ENDIF
3137  RETURN
3138  ENDIF
3139 
3140  IF (modop.EQ.0.AND.modcor.EQ.0) THEN
3141  WRITE(phlun,*) 'PHCORK lack of initialization'
3142  stop
3143  ENDIF
3144 
3145 c execution mode
3146 c ==============
3147 c ==============
3148 
3149 
3150  px=0
3151  py=0
3152  pz=0
3153  e =0
3154 
3155  IF (modop.EQ.1) THEN
3156 c -----------------------
3157 c in this case we do nothing
3158  RETURN
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
3163 
3164  DO i=3,npho
3165 
3166  px=px+ppho(1,i)
3167  py=py+ppho(2,i)
3168  pz=pz+ppho(3,i)
3169 
3170  p2=ppho(1,i)**2+ppho(2,i)**2+ppho(3,i)**2
3171 
3172  en=sqrt( ppho(5,i)**2 + p2)
3173 
3174  IF (iprint.EQ.1)
3175  & WRITE(phlun,*) "CORRECTING ENERGY OF ",i,":",
3176  & ppho(4,i),"=>",en
3177 
3178  ppho(4,i)=en
3179  e = e+ppho(4,i)
3180 
3181  ENDDO
3182 
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
3187 
3188  DO i=3,npho
3189 
3190  px=px+ppho(1,i)
3191  py=py+ppho(2,i)
3192  pz=pz+ppho(3,i)
3193 
3194  p2=ppho(1,i)**2+ppho(2,i)**2+ppho(3,i)**2
3195 
3196  en=sqrt( ppho(5,i)**2 + p2)
3197 
3198  IF (iprint.EQ.1)
3199  & WRITE(phlun,*) "CORRECTING ENERGY OF ",i,":",
3200  & ppho(4,i),"=>",en
3201 
3202  ppho(4,i)=en
3203  e = e+ppho(4,i)
3204 
3205  ENDDO
3206  DO k=1,4
3207  ppho(k,1)=0d0
3208  DO i=3,npho
3209  ppho(k,1)=ppho(k,1)+ppho(k,i)
3210  ENDDO
3211  ENDDO
3212  xms=sqrt(ppho(4,1)**2-ppho(3,1)**2-ppho(2,1)**2-ppho(1,1)**2)
3213  ppho(5,1)=xms
3214  ELSEIF(modop.EQ.3) THEN
3215 c -----------------------
3216 
3217 cc lets loop thru all daughters and correct their masses
3218 cc according to e^2=p^2+m^2
3219 
3220  DO i=3,npho
3221 
3222  px=px+ppho(1,i)
3223  py=py+ppho(2,i)
3224  pz=pz+ppho(3,i)
3225  e = e+ppho(4,i)
3226 
3227  p2=ppho(1,i)**2+ppho(2,i)**2+ppho(3,i)**2
3228 
3229  m=sqrt(abs( ppho(4,i)**2 - p2))
3230 
3231  IF (iprint.EQ.1)
3232  & WRITE(phlun,*) "CORRECTING MASS OF ",i,":",
3233  & ppho(5,i),"=>",m
3234 
3235  ppho(5,i)=m
3236 
3237  ENDDO
3238 
3239 
3240  ELSEIF(modop.EQ.4) THEN
3241 c -----------------------
3242 
3243 cc lets loop thru all daughters and correct their masses
3244 cc or energies according to e^2=p^2+m^2
3245 
3246  DO i=3,npho
3247 
3248  px=px+ppho(1,i)
3249  py=py+ppho(2,i)
3250  pz=pz+ppho(3,i)
3251 
3252  p2=ppho(1,i)**2+ppho(2,i)**2+ppho(3,i)**2
3253 
3254  m=sqrt(abs( ppho(4,i)**2 - p2))
3255 
3256  IF (m.GT.mcut) THEN
3257  IF (iprint.EQ.1)
3258  & WRITE(phlun,*) "CORRECTING MASS OF ",i,":",
3259  & ppho(5,i),"=>",m
3260  ppho(5,i)=m
3261  e = e+ppho(4,i)
3262  ELSE
3263 
3264  en=sqrt( ppho(5,i)**2 + p2)
3265 
3266  IF (iprint.EQ.1)
3267  & WRITE(phlun,*) "CORRECTING ENERGY OF ",i,":",
3268  & ppho(4,i),"=>",en
3269 
3270  ppho(4,i)=en
3271  e = e+ppho(4,i)
3272  ENDIF
3273 
3274  ENDDO
3275  ENDIF
3276 c -----
3277 
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)
3284  ENDIF
3285 
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)
3290 
3291  p2=ppho(1,1)**2+ppho(2,1)**2+ppho(3,1)**2
3292 
3293  IF (ppho(4,1)**2.GT.p2) THEN
3294  m=sqrt( ppho(4,1)**2 - p2 )
3295  IF (iprint.EQ.1)
3296  & WRITE(phlun,*) " M:",ppho(5,1),"=>",m
3297  ppho(5,1)=m
3298  ENDIF
3299 
3300  CALL phlupa(25)
3301 
3302  END
3303 
3304 
3305 
3306  FUNCTION phint(IDUM)
3307 c --- can be used with variant a. for b use phint1 or 2 --------------
3308 c.----------------------------------------------------------------------
3309 c.
3310 c. phint: photos universal interference correction weight
3311 c.
3312 c. purpose: calculates correction weight as expressed by
3313 c formula(17) from cpc 79 (1994), 291.
3314 c.
3315 c. input parameters: Common /phoevt/, with photon added.
3316 c.
3317 c. output parameters: correction weight
3318 c.
3319 c. author(s): z. was, p.golonka created at: 19/01/05
3320 c. last update: 25/01/05
3321 c.
3322 c.----------------------------------------------------------------------
3323  IMPLICIT NONE
3324  real*8 phint,phint2
3325  INTEGER IDUM
3326  INTEGER NMXPHO
3327  parameter(nmxpho=10000)
3328  INTEGER IDPHO,ISTPHO,JDAPHO,JMOPHO,NEVPHO,NPHO
3329  real*8 ppho,vpho
3330  common/phoevt/nevpho,npho,istpho(nmxpho),idpho(nmxpho),
3331  &jmopho(2,nmxpho),jdapho(2,nmxpho),ppho(5,nmxpho),vpho(4,nmxpho)
3332  INTEGER I,K,L
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)
3336  real*8 phocha
3337 c--
3338 
3339 c calculate polarimetric vector: ph, eps1, eps2 are orthogonal
3340 
3341  DO k=1,4
3342  ph(k)=ppho(k,npho)
3343  eps2(k)=1d0
3344  ENDDO
3345 
3346  CALL phoeps(ph,eps2,eps1)
3347  CALL phoeps(ph,eps1,eps2)
3348 
3349 
3350  xnum1=0d0
3351  xnum2=0d0
3352  xdeno=0d0
3353 
3354  DO k=jdapho(1,1),npho-1 ! or JDAPHO(1,2)
3355 
3356 c momenta of charged particle in pl
3357  DO l=1,4
3358  pl(l)=ppho(l,k)
3359  ENDDO
3360 c scalar products: epsilon*p/k*p
3361 
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) )
3365 
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) )
3369 
3370 
3371 c accumulate the currents
3372  xnum1 = xnum1+xc1
3373  xnum2 = xnum2+xc2
3374 
3375  xdeno = xdeno + xc1**2 + xc2**2
3376 
3377  ENDDO
3378 
3379  phint=(xnum1**2 + xnum2**2) / xdeno
3380  phint2=phint
3381 
3382  END
3383 
3384 
3385  SUBROUTINE phoeps (VEC1, VEC2, EPS)
3386 c.----------------------------------------------------------------------
3387 c.
3388 c. phoeps: phoeps vector product(normalized to unity)
3389 c.
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.
3393 c.
3394 c. input parameters: vec1,vec2 - input 4-vectors
3395 c.
3396 c. output parameters: eps - normalized 4-vector, orthogonal to
3397 c vec1 and vec2
3398 c.
3399 c. author(s): z. was, p.golonka created at: 19/01/05
3400 c. last update: 25/01/05
3401 c.
3402 c.----------------------------------------------------------------------
3403 
3404  DOUBLE PRECISION VEC1(4), VEC2(4), EPS(4),XN
3405 
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)
3409  eps(4)=0d0
3410 
3411  xn=sqrt( eps(1)**2 +eps(2)**2 +eps(3)**2)
3412 
3413  eps(1)=eps(1)/xn
3414  eps(2)=eps(2)/xn
3415  eps(3)=eps(3)/xn
3416 
3417 
3418  END
3419  SUBROUTINE phodmp
3420 c.----------------------------------------------------------------------
3421 c.
3422 c. photos: photon radiation in decays event dump routine
3423 c.
3424 c. purpose: print event record.
3425 c.
3426 c. input parameters: Common /hepevt/
3427 c.
3428 c. output parameters: none
3429 c.
3430 c. author(s): b. van eijk created at: 05/06/90
3431 c. last update: 05/06/90
3432 c.
3433 c.----------------------------------------------------------------------
3434 c IMPLICIT NONE
3435  DOUBLE PRECISION SUMVEC(5)
3436  INTEGER I,J
3437 c this is the hepevt class in old style. No d_h_ class pre-name
3438  INTEGER NMXHEP
3439  parameter(nmxhep=10000)
3440  real*8 phep, vhep ! to be real*4/ *8 depending on host
3441  INTEGER nevhep,nhep,isthep,idhep,jmohep,
3442  $ jdahep
3443  COMMON /hepevt/
3444  $ nevhep, ! serial number
3445  $ nhep, ! number of particles
3446  $ isthep(nmxhep), ! status code
3447  $ idhep(nmxhep), ! particle ident KF
3448  $ jmohep(2,nmxhep), ! parent particles
3449  $ jdahep(2,nmxhep), ! childreen particles
3450  $ phep(5,nmxhep), ! four-momentum, mass [GeV]
3451  $ vhep(4,nmxhep) ! vertex [mm]
3452 * ----------------------------------------------------------------------
3453  LOGICAL qedrad
3454  COMMON /phoqed/
3455  $ qedrad(nmxhep) ! Photos flag
3456 * ----------------------------------------------------------------------
3457  SAVE hepevt,phoqed
3458  INTEGER PHLUN
3459  common/pholun/phlun
3460  DO 10 i=1,5
3461  10 sumvec(i)=0.
3462 c--
3463 c-- print event number...
3464  WRITE(phlun,9000)
3465  WRITE(phlun,9010) nevhep
3466  WRITE(phlun,9080)
3467  WRITE(phlun,9020)
3468  DO 30 i=1,nhep
3469 c--
3470 c-- for 'stable particle' calculate vector momentum sum
3471  IF (jdahep(1,i).EQ.0) THEN
3472  DO 20 j=1,4
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)
3476  ELSE
3477  WRITE(phlun,9040) i,idhep(i),jmohep(1,i),jmohep(2,i),(phep
3478  & (j,i),j=1,5)
3479  ENDIF
3480  ELSE
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)
3484  ELSE
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)
3487  ENDIF
3488  ENDIF
3489  30 CONTINUE
3490  sumvec(5)=sqrt(sumvec(4)**2-sumvec(1)**2-sumvec(2)**2-
3491  &sumvec(3)**2)
3492  WRITE(phlun,9070) (sumvec(j),j=1,5)
3493  RETURN
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')
3504  END