C++InterfacetoTauola
photos.F
1 *///////////////////////////////////////////////////////////////////////
2 *//
3 *// !!!!!!! WARNING!!!!! This source may be agressive !!!!
4 *//
5 *// Due to short common block names it may owerwrite variables in other
6 *// parts of the code.
7 *//
8 *// One should add suffix c_Photos_ to names of all commons as soon as
9 *// possible!!
10 *///////////////////////////////////////////////////////////////////////
11 
12 C.----------------------------------------------------------------------
13 C.
14 C. PHOTOS: PHOtos CDE-s
15 C.
16 C. Purpose: Keep definitions for PHOTOS QED correction Monte Carlo.
17 C.
18 C. Input Parameters: None
19 C.
20 C. Output Parameters: None
21 C.
22 C. Author(s): Z. Was, B. van Eijk Created at: 29/11/89
23 C. Last Update: 10/08/93
24 C.
25 C. =========================================================
26 C. General Structure Information: =
27 C. =========================================================
28 C: ROUTINES:
29 C. 1) INITIALIZATION:
30 C. PHOCDE
31 C. PHOINI
32 C. PHOCIN
33 C. PHOINF
34 C. 2) GENERAL INTERFACE:
35 C. PHOTOS
36 C. PHOTOS_GET
37 C. PHOTOS_SET
38 C. PHOTOS_MAKE
39 C. PHOBOS
40 C. PHOIN
41 C. PHOTWO (specific interface
42 C. PHOOUT
43 C. PHOCHK
44 C. PHTYPE (specific interface
45 C. PHOMAK (specific interface
46 C. 3) QED PHOTON GENERATION:
47 C. PHINT
48 C. PHOBW
49 C. PHOPRE
50 C. PHOOMA
51 C. PHOENE
52 C. PHOCOR
53 C. PHOFAC
54 C. PHODO
55 C. 4) UTILITIES:
56 C. PHOTRI
57 C. PHOAN1
58 C. PHOAN2
59 C. PHOBO3
60 C. PHORO2
61 C. PHORO3
62 C. PHORIN
63 C. PHORAN
64 C. PHOCHA
65 C. PHOSPI
66 C. PHOERR
67 C. PHOREP
68 C. PHLUPA
69 C. PHCORK
70 C. IPHQRK
71 C. IPHEKL
72 C. COMMONS:
73 C. NAME USED IN SECT. # OF OCC. Comment
74 C. PHOQED 1) 2) 3 Flags whether emisson to be gen.
75 C. PHOLUN 1) 4) 6 Output device number
76 C. PHOCOP 1) 3) 4 photon coupling & min energy
77 C. PHPICO 1) 3) 4) 5 PI & 2*PI
78 C. PHSEED 1) 4) 3 RN seed
79 C. PHOSTA 1) 4) 3 Status information
80 C. PHOKEY 1) 2) 3) 7 Keys for nonstandard application
81 C. PHOVER 1) 1 Version info for outside
82 C. HEPEVT 2) 2 PDG common
83 C. PH_HEPEVT2) 8 PDG common internal
84 C. PHOEVT 2) 3) 10 PDG branch
85 C. PHOIF 2) 3) 2 emission flags for PDG branch
86 C. PHOMOM 3) 5 param of char-neutr system
87 C. PHOPHS 3) 5 photon momentum parameters
88 C. PHOPRO 3) 4 var. for photon rep. (in branch)
89 C. PHOCMS 2) 3 parameters of boost to branch CMS
90 C. PHNUM 4) 1 event number from outside
91 C.----------------------------------------------------------------------
92  SUBROUTINE phoini
93 C.----------------------------------------------------------------------
94 C.
95 C. PHOTOS: PHOton radiation in decays INItialisation
96 C.
97 C. Purpose: Initialisation routine for the PHOTOS QED radiation
98 C. package. Should be called at least once before a call
99 C. to the steering program 'PHOTOS' is made.
100 C.
101 C. Input Parameters: None
102 C.
103 C. Output Parameters: None
104 C.
105 C. Author(s): Z. Was, B. van Eijk Created at: 26/11/89
106 C. Last Update: 12/04/90
107 C.
108 C.----------------------------------------------------------------------
109  IMPLICIT NONE
110  INTEGER init,idum,iphqrk,iphekl
111  SAVE init
112  DATA init/ 0/
113 C--
114 C-- Return if already initialized...
115  IF (init.NE.0) RETURN
116  init=1
117 C--
118 C-- all the following parameter setters can be called after PHOINI.
119 C-- Initialization of kinematic correction against rounding errors.
120 C-- The set values will be used later if called wit zero.
121 C-- Default parameter is 1 (no correction) optionally 2, 3, 4
122 C-- In case of exponentiation new version 5 is needed in most cases.
123 C-- Definition given here will be thus overwritten in such a case
124 C-- below in routine PHOCIN
125  CALL phcork(1)
126 C-- blocks emission from quarks if parameter is 1 (enables if 2),
127 C-- physical treatment
128 C-- will be 3, option 2 is not realistic and for tests only,
129  idum= iphqrk(1) ! default is 1
130 C-- blocks emission in pi0 to gamma e+ e- if parameter is gt.1
131 C-- (enables otherwise)
132  idum= iphekl(2) ! default is 1
133 C--
134 C-- Preset parameters in PHOTOS commons
135  CALL phocin
136 C--
137 C-- Print info
138  CALL phoinf
139 
140 C--
141 C-- Initialize Marsaglia and Zaman random number generator
142  CALL phorin
143  RETURN
144  END
145  SUBROUTINE phocin
146 C.----------------------------------------------------------------------
147 C.
148 C. PHOTOS: PHOton Common INitialisation
149 C.
150 C. Purpose: Initialisation of parameters in common blocks.
151 C.
152 C. Input Parameters: None
153 C.
154 C. Output Parameters: Commons /PHOLUN/, /PHOPHO/, /PHOCOP/, /PHPICO/
155 C. and /PHSEED/.
156 C.
157 C. Author(s): B. van Eijk Created at: 26/11/89
158 C. Z. Was Last Update: 29/01/05
159 C.
160 C.----------------------------------------------------------------------
161  IMPLICIT NONE
162  INTEGER d_h_nmxhep
163 #if defined (KORALK)
164  parameter(d_h_nmxhep=2000)
165 #elif defined (CZTERYD)
166  parameter(d_h_nmxhep=4000)
167 #elif defined (DWAR)
168  parameter(d_h_nmxhep=2000)
169 #elif defined (DWAD)
170  parameter(d_h_nmxhep=2000)
171 #elif defined (DZIESD)
172  parameter(d_h_nmxhep=10000)
173 #else
174  parameter(d_h_nmxhep=2000)
175 #endif
176  LOGICAL qedrad
177  common/phoqed/qedrad(d_h_nmxhep)
178  INTEGER phlun
179  common/pholun/phlun
180  REAL*8 alpha,xphcut
181  common/phocop/alpha,xphcut
182  REAL*8 pi,twopi
183  common/phpico/pi,twopi
184  INTEGER iseed,i97,j97
185  REAL*8 uran,cran,cdran,cmran
186  common/phseed/iseed(2),i97,j97,uran(97),cran,cdran,cmran
187  INTEGER phomes
188  parameter(phomes=10)
189  INTEGER status
190  common/phosta/status(phomes)
191  LOGICAL interf,isec,itre,iexp,iftop,ifw
192  REAL*8 fint,fsec,expeps
193  COMMON /phokey/ fsec,fint,expeps,interf,isec,itre,iexp,iftop,ifw
194  INTEGER init,i
195  SAVE init
196  DATA init/ 0/
197 C--
198 C-- Return if already initialized...
199  IF (init.NE.0) RETURN
200  init=1
201 C--
202 C-- Preset switch for photon emission to 'TRUE' for each particle in
203 C-- /PH_HEPEVT/, this interface is needed for KORALB and KORALZ...
204  DO 10 i=1,d_h_nmxhep
205  10 qedrad(i)=.true.
206 C--
207 C-- Logical output unit for printing of PHOTOS error messages
208  phlun=6
209 C--
210 C-- Set cut parameter for photon radiation
211  xphcut=0.01 d0 ! 0.0001D0! to go to low valuex (IEXP excepted)
212 C-- ! switch to - VARIANT B
213 C--
214 C-- Define some constants
215  alpha=0.00729735039d0
216  pi=3.14159265358979324d0
217  twopi=6.28318530717958648d0
218 C--
219 C-- Default seeds Marsaglia and Zaman random number generator
220  iseed(1)=1802
221  iseed(2)=9373
222 C--
223 C-- Iitialization for extra options
224 C-- (1)
225 C-- Interference weight now universal.
226  interf=.true.
227 C-- (2)
228 C-- Second order - double photon switch
229  isec=.true.
230 C-- Third/fourth order - triple (or quatric) photon switch,
231 C-- see dipswitch ifour
232  itre=.false.
233 C-- Exponentiation on:
234  iexp=.false. !.TRUE.
235  IF (iexp) THEN
236  isec=.false.
237  itre=.false.
238  CALL phcork(5) ! in case of exponentiation correction of ph space
239  ! is a default mandatory
240  xphcut=0.000 000 1
241  expeps=1d-4
242  ENDIF
243 C-- (3)
244 C-- Emision in the hard process g g (q qbar) --> t tbar
245 C-- t --> W b
246  iftop=.true.
247 C--
248 C-- further initialization done automatically
249 C-- see places with - VARIANT A - VARIANT B - all over
250 C-- to switch between options.
251 C ----------- SLOWER VARIANT A, but stable ------------------
252 C --- it is limiting choice for small XPHCUT in fixed orer
253 C --- modes of operation
254  IF (interf) THEN
255 C-- best choice is if FINT=2**N where N+1 is maximal number
256 C-- of charged daughters
257 C-- see report on overweihted events
258  fint=2.0d0
259  ELSE
260  fint=1.0d0
261  ENDIF
262 C ----------- FASTER VARIANT B ------------------
263 C -- it is good for tests of fixed order and small XPHCUT
264 C -- but is less promising for more complex cases of interference
265 C -- sometimes fails because of that
266 C
267 C IF (INTERF) THEN
268 C FINT=1.80D0
269 C ELSE
270 C FINT=0.0D0
271 C ENDIF
272 C----------END VARIANTS A B -----------------------
273 
274 C-- Effects of initial state charge (in leptonic W decays)
275 C--
276  ifw=.true.
277 C-- Initialise status counter for warning messages
278  DO 20 i=1,phomes
279  20 status(i)=0
280  RETURN
281  END
282  SUBROUTINE phoinf
283 C.----------------------------------------------------------------------
284 C.
285 C. PHOTOS: PHOton radiation in decays general INFo
286 C.
287 C. Purpose: Print PHOTOS info
288 C.
289 C. Input Parameters: PHOLUN
290 C.
291 C. Output Parameters: PHOVN1, PHOVN2
292 C.
293 C. Author(s): B. van Eijk Created at: 12/04/90
294 C. Last Update: 27/06/04
295 C.
296 C.----------------------------------------------------------------------
297  IMPLICIT NONE
298  INTEGER iv1,iv2,iv3
299  INTEGER phovn1,phovn2
300  common/phover/phovn1,phovn2
301  INTEGER phlun
302  common/pholun/phlun
303  LOGICAL interf,isec,itre,iexp,iftop,ifw
304  REAL*8 fint,fsec,expeps
305  COMMON /phokey/ fsec,fint,expeps,interf,isec,itre,iexp,iftop,ifw
306  REAL*8 alpha,xphcut
307  common/phocop/alpha,xphcut
308 C--
309 C-- PHOTOS version number and release date
310  phovn1=215
311  phovn2=111005
312 C--
313 C-- Print info
314  WRITE(phlun,9000)
315  WRITE(phlun,9020)
316  WRITE(phlun,9010)
317  WRITE(phlun,9030)
318  iv1=phovn1/100
319  iv2=phovn1-iv1*100
320  WRITE(phlun,9040) iv1,iv2
321  iv1=phovn2/10000
322  iv2=(phovn2-iv1*10000)/100
323  iv3=phovn2-iv1*10000-iv2*100
324  WRITE(phlun,9050) iv1,iv2,iv3
325  WRITE(phlun,9030)
326  WRITE(phlun,9010)
327  WRITE(phlun,9060)
328  WRITE(phlun,9010)
329  WRITE(phlun,9070)
330  WRITE(phlun,9010)
331  WRITE(phlun,9020)
332  WRITE(phlun,9010)
333  WRITE(phlun,9064) interf,isec,itre,iexp,iftop,ifw,alpha,xphcut
334  WRITE(phlun,9010)
335  IF (interf) WRITE(phlun,9061)
336  IF (isec) WRITE(phlun,9062)
337  IF (itre) WRITE(phlun,9066)
338  IF (iexp) WRITE(phlun,9067) expeps
339  IF (iftop) WRITE(phlun,9063)
340  IF (ifw) WRITE(phlun,9065)
341  WRITE(phlun,9080)
342  WRITE(phlun,9010)
343  WRITE(phlun,9020)
344  RETURN
345  9000 FORMAT(1h1)
346  9010 FORMAT(1h ,'*',t81,'*')
347  9020 FORMAT(1h ,80('*'))
348  9030 FORMAT(1h ,'*',26x,26('='),t81,'*')
349  9040 FORMAT(1h ,'*',28x,'PHOTOS, Version: ',i2,'.',i2,t81,'*')
350  9050 FORMAT(1h ,'*',28x,'Released at: ',i2,'/',i2,'/',i2,t81,'*')
351  9060 FORMAT(1h ,'*',18x,'PHOTOS QED Corrections in Particle Decays',
352  &t81,'*')
353  9061 FORMAT(1h ,'*',18x,'option with interference is active ',
354  &t81,'*')
355  9062 FORMAT(1h ,'*',18x,'option with double photons is active ',
356  &t81,'*')
357  9066 FORMAT(1h ,'*',18x,'option with triple/quatric photons is active',
358  &t81,'*')
359  9067 FORMAT(1h ,'*',18x,'option with exponentiation is active EPSEXP=',
360  &e10.4,t81,'*')
361  9063 FORMAT(1h ,'*',18x,'emision in t tbar production is active ',
362  &t81,'*')
363  9064 FORMAT(1h ,'*',18x,'Internal input parameters:',t81,'*'
364  &,/, 1h ,'*',t81,'*'
365  &,/, 1h ,'*',18x,'INTERF=',l2,' ISEC=',l2,' ITRE=',l2,
366  & ' IEXP=',l2,' IFTOP=',l2,
367  & ' IFW=',l2,t81,'*'
368  &,/, 1h ,'*',18x,'ALPHA_QED=',f8.5,' XPHCUT=',e8.3,t81,'*')
369  9065 FORMAT(1h ,'*',18x,'correction wt in decay of W is active ',
370  &t81,'*')
371  9070 FORMAT(1h ,'*',9x,
372  &'Monte Carlo Program - by E. Barberio, B. van Eijk and Z. Was',
373  & t81,'*',/,
374  & 1h ,'*',9x,'Version 2.09 - by P. Golonka and Z.W.',t81,'*')
375  9080 FORMAT( 1h ,'*',9x,' ',t81,'*',/,
376  & 1h ,'*',9x,
377  & ' WARNING (1): /HEPEVT/ is not anymore the standard common block'
378  & ,t81,'*',/,
379  & 1h ,'*',9x,' ',t81,'*',/,
380  & 1h ,'*',9x,
381 #if defined (KORALK)
382  & ' PHOTOS expects /d_HEPEVT/ to have REAL*8 variabl. To change to'
383  & ,t81,'*',/, 1h ,'*',9x,
384  & ' REAL*4 modify its declaration in subr. PHOTOS_GET PHOTOS_SET:'
385  & ,t81,'*',/, 1h ,'*',9x,
386  & ' REAL*8 d_h_phep, d_h_vhep'
387  & ,t81,'*',/, 1h ,'*',9x,
388  & ' WARNING (2): check dims. of /d_hepevt/ /phoqed/ /ph_hepevt/.'
389  & ,t81,'*',/, 1h ,'*',9x,
390  & ' HERE: d_h_nmxhep=2000 and NMXHEP=10000'
391 #elif defined (CZTERYD)
392  & ' PHOTOS expects /HEPEVT/ to have REAL*8 variables. To change to'
393  & ,t81,'*',/, 1h ,'*',9x,
394  & ' REAL*4 modify its declaration in subr. PHOTOS_GET PHOTOS_SET:'
395  & ,t81,'*',/, 1h ,'*',9x,
396  & ' REAL*8 d_h_phep, d_h_vhep'
397  & ,t81,'*',/, 1h ,'*',9x,
398  & ' WARNING (2): check dims. of /hepevt/ /phoqed/ /ph_hepevt/.'
399  & ,t81,'*',/, 1h ,'*',9x,
400  & ' HERE: d_h_nmxhep=4000 and NMXHEP=10000'
401 #elif defined (DWAR)
402  & ' PHOTOS expects /HEPEVT/ to have REAL*4 variables. To change to'
403  & ,t81,'*',/, 1h ,'*',9x,
404  & ' REAL*8 modify its declaration in subr. PHOTOS_GET PHOTOS_SET:'
405  & ,t81,'*',/, 1h ,'*',9x,
406  & ' REAL*8 d_h_phep, d_h_vhep'
407  & ,t81,'*',/, 1h ,'*',9x,
408  & ' WARNING (2): check dims. of /hepevt/ /phoqed/ /ph_hepevt/.'
409  & ,t81,'*',/, 1h ,'*',9x,
410  & ' HERE: d_h_nmxhep=2000 and NMXHEP=10000'
411 #elif defined (DWAD)
412  & ' PHOTOS expects /HEPEVT/ to have REAL*8 variables. To change to'
413  & ,t81,'*',/, 1h ,'*',9x,
414  & ' REAL*4 modify its declaration in subr. PHOTOS_GET PHOTOS_SET:'
415  & ,t81,'*',/, 1h ,'*',9x,
416  & ' REAL*8 d_h_phep, d_h_vhep'
417  & ,t81,'*',/, 1h ,'*',9x,
418  & ' WARNING (2): check dims. of /hepevt/ /phoqed/ /ph_hepevt/.'
419  & ,t81,'*',/, 1h ,'*',9x,
420  & ' HERE: d_h_nmxhep=2000 and NMXHEP=10000'
421 #elif defined (DZIESD)
422  & ' PHOTOS expects /HEPEVT/ to have REAL*8 variables. To change to'
423  & ,t81,'*',/, 1h ,'*',9x,
424  & ' REAL*4 modify its declaration in subr. PHOTOS_GET PHOTOS_SET:'
425  & ,t81,'*',/, 1h ,'*',9x,
426  & ' REAL*8 d_h_phep, d_h_vhep'
427  & ,t81,'*',/, 1h ,'*',9x,
428  & ' WARNING (2): check dims. of /hepevt/ /phoqed/ /ph_hepevt/.'
429  & ,t81,'*',/, 1h ,'*',9x,
430  & ' HERE: d_h_nmxhep=10000 and NMXHEP=10000'
431 #else
432  & ' PHOTOS expects /HEPEVT/ to have REAL*8 variables. To change to'
433  & ,t81,'*',/, 1h ,'*',9x,
434  & ' REAL*4 modify its declaration in subr. PHOTOS_GET PHOTOS_SET:'
435  & ,t81,'*',/, 1h ,'*',9x,
436  & ' REAL*8 d_h_phep, d_h_vhep'
437  & ,t81,'*',/, 1h ,'*',9x,
438  & ' WARNING (2): check dims. of /hepevt/ /phoqed/ /ph_hepevt/.'
439  & ,t81,'*',/, 1h ,'*',9x,
440  & ' HERE: d_h_nmxhep=2000 and NMXHEP=10000'
441 #endif
442  & ,t81,'*')
443  END
444  SUBROUTINE photos(ID)
445 C.----------------------------------------------------------------------
446 C.
447 C. PHOTOS: General search routine + _GET + _SET
448 C.
449 C. Purpose: /HEPEVT/ is not anymore a standard at least
450 C. REAL*8 REAL*4 are in use. PHOTOS_GET and PHOTOS_SET
451 C. were to be introduced.
452 C.
453 C.
454 C. Input Parameters: ID see routine PHOTOS_MAKE
455 C.
456 C. Output Parameters: None
457 C.
458 C. Author(s): Z. Was Created at: 21/07/98
459 C. Last Update: 21/07/98
460 C.
461 C.----------------------------------------------------------------------
462  COMMON /phlupy/ ipoin,ipoinm
463  INTEGER ipoin,ipoinm
464  COMMON /phnum/ iev
465  INTEGER iev
466  INTEGER phlun
467  common/pholun/phlun
468 
469  IF (1.GT.ipoinm.AND.1.LT.ipoin ) THEN
470  WRITE(phlun,*) 'EVENT NR=',iev,
471  $ 'WE ARE TESTING /HEPEVT/ at IPOINT=1 (input)'
472  CALL phodmp
473  ENDIF
474  CALL photos_get
475  CALL photos_make(id)
476  CALL photos_set
477  IF (1.GT.ipoinm.AND.1.LT.ipoin ) THEN
478  WRITE(phlun,*) 'EVENT NR=',iev,
479  $ 'WE ARE TESTING /HEPEVT/ at IPOINT=1 (output)'
480  CALL phodmp
481  ENDIF
482 
483  END
484 
485  SUBROUTINE photos_get
486 C.----------------------------------------------------------------------
487 C.
488 C. Getter for PHOTOS:
489 C.
490 C. Purpose: Copies /HEPEVT/ into /PH_HEPEVT/
491 C.
492 C.
493 C. Input Parameters: None
494 C.
495 C. Output Parameters: None
496 C.
497 C. Author(s): Z. Was Created at: 21/07/98
498 C. Last Update: 21/07/98
499 C.
500 C.----------------------------------------------------------------------
501 
502  IMPLICIT NONE
503  INTEGER d_h_nmxhep ! maximum number of particles
504 #if defined (KORALK)
505  parameter(d_h_nmxhep=2000)
506  REAL*8 d_h_phep, d_h_vhep ! to be real*4/ *8 depending on host
507  INTEGER d_h_nevhep,d_h_nhep,d_h_isthep,d_h_idhep,d_h_jmohep,
508  $ d_h_jdahep
509  COMMON /d_hepevt/
510 #elif defined (CZTERYD)
511  parameter(d_h_nmxhep=4000)
512  REAL*8 d_h_phep, d_h_vhep ! to be real*4/ *8 depending on host
513  INTEGER d_h_nevhep,d_h_nhep,d_h_isthep,d_h_idhep,d_h_jmohep,
514  $ d_h_jdahep
515  COMMON /hepevt/
516 #elif defined (DWAR)
517  parameter(d_h_nmxhep=2000)
518  REAL*4 d_h_phep, d_h_vhep ! to be real*4/ *8 depending on host
519  INTEGER d_h_nevhep,d_h_nhep,d_h_isthep,d_h_idhep,d_h_jmohep,
520  $ d_h_jdahep
521  COMMON /hepevt/
522 #elif defined (DWAD)
523  parameter(d_h_nmxhep=2000)
524  REAL*8 d_h_phep, d_h_vhep ! to be real*4/ *8 depending on host
525  INTEGER d_h_nevhep,d_h_nhep,d_h_isthep,d_h_idhep,d_h_jmohep,
526  $ d_h_jdahep
527  COMMON /hepevt/
528 #elif defined (DZIESD)
529  parameter(d_h_nmxhep=10000)
530  REAL*8 d_h_phep, d_h_vhep ! to be real*4/ *8 depending on host
531  INTEGER d_h_nevhep,d_h_nhep,d_h_isthep,d_h_idhep,d_h_jmohep,
532  $ d_h_jdahep
533  COMMON /hepevt/
534 #else
535  parameter(d_h_nmxhep=2000)
536  REAL*8 d_h_phep, d_h_vhep ! to be real*4/ *8 depending on host
537  INTEGER d_h_nevhep,d_h_nhep,d_h_isthep,d_h_idhep,d_h_jmohep,
538  $ d_h_jdahep
539  COMMON /hepevt/
540 #endif
541  $ d_h_nevhep, ! serial number
542  $ d_h_nhep, ! number of particles
543  $ d_h_isthep(d_h_nmxhep), ! status code
544  $ d_h_idhep(d_h_nmxhep), ! particle ident KF
545  $ d_h_jmohep(2,d_h_nmxhep), ! parent particles
546  $ d_h_jdahep(2,d_h_nmxhep), ! childreen particles
547  $ d_h_phep(5,d_h_nmxhep), ! four-momentum, mass [GeV]
548  $ d_h_vhep(4,d_h_nmxhep) ! vertex [mm]
549 * ----------------------------------------------------------------------
550  LOGICAL d_h_qedrad
551  COMMON /phoqed/
552  $ d_h_qedrad(d_h_nmxhep) ! Photos flag
553  INTEGER nmxhep
554  parameter(nmxhep=10000)
555  INTEGER idhep,isthep,jdahep,jmohep,nevhep,nhep
556  REAL*8 phep,vhep
557  common/ph_hepevt/nevhep,nhep,isthep(nmxhep),idhep(nmxhep),
558  &jmohep(2,nmxhep),jdahep(2,nmxhep),phep(5,nmxhep),vhep(4,nmxhep)
559  LOGICAL qedrad
560  common/ph_phoqed/qedrad(nmxhep)
561  integer k,l
562  nevhep= d_h_nevhep ! serial number
563  nhep = d_h_nhep ! number of particles
564  DO k=1,nhep
565  isthep(k) =d_h_isthep(k) ! status code
566  idhep(k) =d_h_idhep(k) ! particle ident KF
567  jmohep(1,k) =d_h_jmohep(1,k) ! parent particles
568  jdahep(1,k) =d_h_jdahep(1,k) ! childreen particles
569  jmohep(2,k) =d_h_jmohep(2,k) ! parent particles
570  jdahep(2,k) =d_h_jdahep(2,k) ! childreen particles
571  DO l=1,4
572  phep(l,k) =d_h_phep(l,k) ! four-momentum, mass [GeV]
573  vhep(l,k) =d_h_vhep(l,k) ! vertex [mm]
574  ENDDO
575  phep(5,k) =d_h_phep(5,k) ! four-momentum, mass [GeV]
576  qedrad(k) =d_h_qedrad(k) ! Photos special flag
577  ENDDO
578  END
579 
580 
581  SUBROUTINE photos_set
582 C.----------------------------------------------------------------------
583 C.
584 C. Setter for PHOTOS:
585 C.
586 C. Purpose: Copies /PH_HEPEVT/ into /HEPEVT/
587 C.
588 C.
589 C. Input Parameters: None
590 C.
591 C. Output Parameters: None
592 C.
593 C. Author(s): Z. Was Created at: 21/07/98
594 C. Last Update: 21/07/98
595 C.
596 C.----------------------------------------------------------------------
597  IMPLICIT NONE
598  INTEGER d_h_nmxhep ! maximum number of particles
599 #if defined (KORALK)
600  parameter(d_h_nmxhep=2000)
601  REAL*8 d_h_phep, d_h_vhep ! to be real*4/ *8 depending on host
602  INTEGER d_h_nevhep,d_h_nhep,d_h_isthep,d_h_idhep,d_h_jmohep,
603  $ d_h_jdahep
604  COMMON /d_hepevt/
605 #elif defined (CZTERYD)
606  parameter(d_h_nmxhep=4000)
607  REAL*8 d_h_phep, d_h_vhep ! to be real*4/ *8 depending on host
608  INTEGER d_h_nevhep,d_h_nhep,d_h_isthep,d_h_idhep,d_h_jmohep,
609  $ d_h_jdahep
610  COMMON /hepevt/
611 #elif defined (DWAR)
612  parameter(d_h_nmxhep=2000)
613  REAL*4 d_h_phep, d_h_vhep ! to be real*4/ *8 depending on host
614  INTEGER d_h_nevhep,d_h_nhep,d_h_isthep,d_h_idhep,d_h_jmohep,
615  $ d_h_jdahep
616  COMMON /hepevt/
617 #elif defined (DWAD)
618  parameter(d_h_nmxhep=2000)
619  REAL*8 d_h_phep, d_h_vhep ! to be real*4/ *8 depending on host
620  INTEGER d_h_nevhep,d_h_nhep,d_h_isthep,d_h_idhep,d_h_jmohep,
621  $ d_h_jdahep
622  COMMON /hepevt/
623 #elif defined (DZIESD)
624  parameter(d_h_nmxhep=10000)
625  REAL*8 d_h_phep, d_h_vhep ! to be real*4/ *8 depending on host
626  INTEGER d_h_nevhep,d_h_nhep,d_h_isthep,d_h_idhep,d_h_jmohep,
627  $ d_h_jdahep
628  COMMON /hepevt/
629 #else
630  parameter(d_h_nmxhep=2000)
631  REAL*8 d_h_phep, d_h_vhep ! to be real*4/ *8 depending on host
632  INTEGER d_h_nevhep,d_h_nhep,d_h_isthep,d_h_idhep,d_h_jmohep,
633  $ d_h_jdahep
634  COMMON /hepevt/
635 #endif
636  $ d_h_nevhep, ! serial number
637  $ d_h_nhep, ! number of particles
638  $ d_h_isthep(d_h_nmxhep), ! status code
639  $ d_h_idhep(d_h_nmxhep), ! particle ident KF
640  $ d_h_jmohep(2,d_h_nmxhep), ! parent particles
641  $ d_h_jdahep(2,d_h_nmxhep), ! childreen particles
642  $ d_h_phep(5,d_h_nmxhep), ! four-momentum, mass [GeV]
643  $ d_h_vhep(4,d_h_nmxhep) ! vertex [mm]
644 * ----------------------------------------------------------------------
645  LOGICAL d_h_qedrad
646  COMMON /phoqed/
647  $ d_h_qedrad(d_h_nmxhep) ! Photos flag
648  INTEGER nmxhep
649  parameter(nmxhep=10000)
650  INTEGER idhep,isthep,jdahep,jmohep,nevhep,nhep
651  REAL*8 phep,vhep
652  common/ph_hepevt/nevhep,nhep,isthep(nmxhep),idhep(nmxhep),
653  &jmohep(2,nmxhep),jdahep(2,nmxhep),phep(5,nmxhep),vhep(4,nmxhep)
654  LOGICAL qedrad
655  common/ph_phoqed/qedrad(nmxhep)
656  INTEGER k,l
657 
658  d_h_nevhep= nevhep ! serial number
659  d_h_nhep = nhep ! number of particles
660  DO k=1,nhep
661  d_h_isthep(k) =isthep(k) ! status code
662  d_h_idhep(k) =idhep(k) ! particle ident KF
663  d_h_jmohep(1,k) =jmohep(1,k) ! parent particles
664  d_h_jdahep(1,k) =jdahep(1,k) ! childreen particles
665  d_h_jmohep(2,k) =jmohep(2,k) ! parent particles
666  d_h_jdahep(2,k) =jdahep(2,k) ! childreen particles
667  DO l=1,4
668  d_h_phep(l,k) =phep(l,k) ! four-momentum, mass [GeV]
669  d_h_vhep(l,k) =vhep(l,k) ! vertex [mm]
670  ENDDO
671  d_h_phep(5,k) =phep(5,k) ! four-momentum, mass [GeV]
672  d_h_qedrad(k) =qedrad(k) ! Photos special flag
673  ENDDO
674  END
675  SUBROUTINE photos_make(IPARR)
676 C.----------------------------------------------------------------------
677 C.
678 C. PHOTOS_MAKE: General search routine
679 C.
680 C. Purpose: Search through the /PH_HEPEVT/ standard HEP common, sta-
681 C. rting from the IPPAR-th particle. Whenevr branching
682 C. point is found routine PHTYPE(IP) is called.
683 C. Finally if calls on PHTYPE(IP) modified entries, common
684 C /PH_HEPEVT/ is ordered.
685 C.
686 C. Input Parameter: IPPAR: Pointer to decaying particle in
687 C. /PH_HEPEVT/ and the common itself,
688 C.
689 C. Output Parameters: Common /PH_HEPEVT/, either with or without
690 C. new particles added.
691 C.
692 C. Author(s): Z. Was, B. van Eijk Created at: 26/11/89
693 C. Last Update: 30/08/93
694 C.
695 C.----------------------------------------------------------------------
696  IMPLICIT NONE
697  REAL*8 photon(5)
698  INTEGER ip,iparr,ippar,i,j,k,l,nlast
699  DOUBLE PRECISION data
700  INTEGER mother,pospho
701  LOGICAL cascad
702  INTEGER nmxhep
703  parameter(nmxhep=10000)
704  INTEGER idhep,isthep,jdahep,jmohep,nevhep,nhep
705  REAL*8 phep,vhep
706  common/ph_hepevt/nevhep,nhep,isthep(nmxhep),idhep(nmxhep),
707  &jmohep(2,nmxhep),jdahep(2,nmxhep),phep(5,nmxhep),vhep(4,nmxhep)
708  LOGICAL qedrad
709  common/ph_phoqed/qedrad(nmxhep)
710  INTEGER nmxpho
711  parameter(nmxpho=10000)
712  INTEGER istack(0:nmxpho),numit,ntry,kk,ll,ii,na,first,last
713  INTEGER firsta,lasta,ipp,ida1,ida2,mother2,idpho,ispho
714  REAL*8 porig(5,nmxpho)
715 C--
716  ippar=abs(iparr)
717 C-- Store pointers for cascade treatement...
718  ip=ippar
719  nlast=nhep
720  cascad=.false.
721 C--
722 C-- Check decay multiplicity and minimum of correctness..
723  IF ((jdahep(1,ip).EQ.0).OR.(jmohep(1,jdahep(1,ip)).NE.ip)) RETURN
724 C--
725 C-- single branch mode
726 C-- we start looking for the decay points in the cascade
727 C-- IPPAR is original position where the program was called
728  istack(0)=ippar
729 C-- NUMIT denotes number of secondary decay branches
730  numit=0
731 C-- NTRY denotes number of secondary branches already checked for
732 C-- for existence of further branches
733  ntry=0
734 C-- let-s search if IPARR does not prevent searching.
735  IF (iparr.GT.0) THEN
736  30 CONTINUE
737  DO i=jdahep(1,ip),jdahep(2,ip)
738  IF (jdahep(1,i).NE.0.AND.jmohep(1,jdahep(1,i)).EQ.i) THEN
739  numit=numit+1
740  IF (numit.GT.nmxpho) THEN
741  data=numit
742  CALL phoerr(7,'PHOTOS',data)
743  ENDIF
744  istack(numit)=i
745  ENDIF
746  ENDDO
747  IF(numit.GT.ntry) THEN
748  ntry=ntry+1
749  ip=istack(ntry)
750  goto 30
751  ENDIF
752  ENDIF
753 C-- let-s do generation
754  DO 25 kk=0,numit
755  na=nhep
756  first=jdahep(1,istack(kk))
757  last=jdahep(2,istack(kk))
758  DO ii=1,last-first+1
759  DO ll=1,5
760  porig(ll,ii)=phep(ll,first+ii-1)
761  ENDDO
762  ENDDO
763 C--
764  CALL phtype(istack(kk))
765 C--
766 C-- Correct energy/momentum of cascade daughters
767  IF(nhep.GT.na) THEN
768  DO ii=1,last-first+1
769  ipp=first+ii-1
770  firsta=jdahep(1,ipp)
771  lasta=jdahep(2,ipp)
772  IF(jmohep(1,ipp).EQ.istack(kk))
773  $ CALL phobos(ipp,porig(1,ii),phep(1,ipp),firsta,lasta)
774  ENDDO
775  ENDIF
776  25 CONTINUE
777 C--
778 C-- rearrange /PH_HEPEVT/ to get correct order..
779  IF (nhep.GT.nlast) THEN
780  DO 160 i=nlast+1,nhep
781 C--
782 C-- Photon mother and position...
783  mother=jmohep(1,i)
784  pospho=jdahep(2,mother)+1
785 C-- Intermediate save of photon energy/momentum and pointers
786  DO 90 j=1,5
787  90 photon(j)=phep(j,i)
788  ispho =isthep(i)
789  idpho =idhep(i)
790  mother2 =jmohep(2,i)
791  ida1 =jdahep(1,i)
792  ida2 =jdahep(2,i)
793 C--
794 C-- Exclude photon in sequence !
795  IF (pospho.NE.nhep) THEN
796 C--
797 C--
798 C-- Order /PH_HEPEVT/
799  DO 120 k=i,pospho+1,-1
800  isthep(k)=isthep(k-1)
801  qedrad(k)=qedrad(k-1)
802  idhep(k)=idhep(k-1)
803  DO 100 l=1,2
804  jmohep(l,k)=jmohep(l,k-1)
805  100 jdahep(l,k)=jdahep(l,k-1)
806  DO 110 l=1,5
807  110 phep(l,k)=phep(l,k-1)
808  DO 120 l=1,4
809  120 vhep(l,k)=vhep(l,k-1)
810 C--
811 C-- Correct pointers assuming most dirty /PH_HEPEVT/...
812  DO 130 k=1,nhep
813  DO 130 l=1,2
814  IF ((jmohep(l,k).NE.0).AND.(jmohep(l,k).GE.
815  & pospho)) jmohep(l,k)=jmohep(l,k)+1
816  IF ((jdahep(l,k).NE.0).AND.(jdahep(l,k).GE.
817  & pospho)) jdahep(l,k)=jdahep(l,k)+1
818  130 CONTINUE
819 C--
820 C-- Store photon energy/momentum
821  DO 140 j=1,5
822  140 phep(j,pospho)=photon(j)
823  ENDIF
824 C--
825 C-- Store pointers for the photon...
826  jdahep(2,mother)=pospho
827  isthep(pospho)=ispho
828  idhep(pospho)=idpho
829  jmohep(1,pospho)=mother
830  jmohep(2,pospho)=mother2
831  jdahep(1,pospho)=ida1
832  jdahep(2,pospho)=ida2
833 C--
834 C-- Get photon production vertex position
835  DO 150 j=1,4
836  150 vhep(j,pospho)=vhep(j,pospho-1)
837  160 CONTINUE
838  ENDIF
839  RETURN
840  END
841  SUBROUTINE phobos(IP,PBOOS1,PBOOS2,FIRST,LAST)
842 C.----------------------------------------------------------------------
843 C.
844 C. PHOBOS: PHOton radiation in decays BOoSt routine
845 C.
846 C. Purpose: Boost particles in cascade decay to parent rest frame
847 C. and boost back with modified boost vector.
848 C.
849 C. Input Parameters: IP: pointer of particle starting chain
850 C. to be boosted
851 C. PBOOS1: Boost vector to rest frame,
852 C. PBOOS2: Boost vector to modified frame,
853 C. FIRST: Pointer to first particle to be boos-
854 C. ted (/PH_HEPEVT/),
855 C. LAST: Pointer to last particle to be boos-
856 C. ted (/PH_HEPEVT/).
857 C.
858 C. Output Parameters: Common /PH_HEPEVT/.
859 C.
860 C. Author(s): B. van Eijk Created at: 13/02/90
861 C. Z. Was Last Update: 16/11/93
862 C.
863 C.----------------------------------------------------------------------
864  IMPLICIT NONE
865  DOUBLE PRECISION bet1(3),bet2(3),gam1,gam2,pb,data
866  INTEGER i,j,first,last,maxsta,nstack,ip
867  parameter(maxsta=10000)
868  INTEGER stack(maxsta)
869  REAL*8 pboos1(5),pboos2(5)
870  INTEGER nmxhep
871  parameter(nmxhep=10000)
872  INTEGER idhep,isthep,jdahep,jmohep,nevhep,nhep
873  REAL*8 phep,vhep
874  common/ph_hepevt/nevhep,nhep,isthep(nmxhep),idhep(nmxhep),
875  &jmohep(2,nmxhep),jdahep(2,nmxhep),phep(5,nmxhep),vhep(4,nmxhep)
876  IF ((last.EQ.0).OR.(last.LT.first)) RETURN
877  nstack=0
878  DO 10 j=1,3
879  bet1(j)=-pboos1(j)/pboos1(5)
880  10 bet2(j)=pboos2(j)/pboos2(5)
881  gam1=pboos1(4)/pboos1(5)
882  gam2=pboos2(4)/pboos2(5)
883 C--
884 C-- Boost vector to parent rest frame...
885  20 DO 50 i=first,last
886  pb=bet1(1)*phep(1,i)+bet1(2)*phep(2,i)+bet1(3)*phep(3,i)
887  IF (jmohep(1,i).EQ.ip) THEN
888  DO 30 j=1,3
889  30 phep(j,i)=phep(j,i)+bet1(j)*(phep(4,i)+pb/(gam1+1.d0))
890  phep(4,i)=gam1*phep(4,i)+pb
891 C--
892 C-- ...and boost back to modified parent frame.
893  pb=bet2(1)*phep(1,i)+bet2(2)*phep(2,i)+bet2(3)*phep(3,i)
894  DO 40 j=1,3
895  40 phep(j,i)=phep(j,i)+bet2(j)*(phep(4,i)+pb/(gam2+1.d0))
896  phep(4,i)=gam2*phep(4,i)+pb
897  IF (jdahep(1,i).NE.0) THEN
898  nstack=nstack+1
899 C--
900 C-- Check on stack length...
901  IF (nstack.GT.maxsta) THEN
902  data=nstack
903  CALL phoerr(7,'PHOBOS',data)
904  ENDIF
905  stack(nstack)=i
906  ENDIF
907  ENDIF
908  50 CONTINUE
909  IF (nstack.NE.0) THEN
910 C--
911 C-- Now go one step further in the decay tree...
912  first=jdahep(1,stack(nstack))
913  last=jdahep(2,stack(nstack))
914  ip=stack(nstack)
915  nstack=nstack-1
916  goto 20
917  ENDIF
918  RETURN
919  END
920  SUBROUTINE phoin(IP,BOOST,NHEP0)
921 C.----------------------------------------------------------------------
922 C.
923 C. PHOIN: PHOtos INput
924 C.
925 C. Purpose: copies IP branch of the common /PH_HEPEVT/ into /PHOEVT/
926 C. moves branch into its CMS system.
927 C.
928 C. Input Parameters: IP: pointer of particle starting branch
929 C. to be copied
930 C. BOOST: Flag whether boost to CMS was or was
931 C . not performed.
932 C.
933 C. Output Parameters: Commons: /PHOEVT/, /PHOCMS/
934 C.
935 C. Author(s): Z. Was Created at: 24/05/93
936 C. Last Update: 16/11/93
937 C.
938 C.----------------------------------------------------------------------
939  IMPLICIT NONE
940  INTEGER nmxhep
941  parameter(nmxhep=10000)
942  INTEGER idhep,isthep,jdahep,jmohep,nevhep,nhep
943  REAL*8 phep,vhep
944  common/ph_hepevt/nevhep,nhep,isthep(nmxhep),idhep(nmxhep),
945  &jmohep(2,nmxhep),jdahep(2,nmxhep),phep(5,nmxhep),vhep(4,nmxhep)
946  INTEGER nmxpho
947  parameter(nmxpho=10000)
948  INTEGER idpho,istpho,jdapho,jmopho,nevpho,npho
949  REAL*8 ppho,vpho
950  common/phoevt/nevpho,npho,istpho(nmxpho),idpho(nmxpho),
951  &jmopho(2,nmxpho),jdapho(2,nmxpho),ppho(5,nmxpho),vpho(4,nmxpho)
952  INTEGER ip,ip2,i,first,last,ll,na
953  LOGICAL boost
954  INTEGER j,nhep0
955  DOUBLE PRECISION bet(3),gam,pb
956  COMMON /phocms/ bet,gam
957  LOGICAL interf,isec,itre,iexp,iftop,ifw
958  REAL*8 fint,fsec,expeps
959  COMMON /phokey/ fsec,fint,expeps,interf,isec,itre,iexp,iftop,ifw
960 C--
961 C let-s calculate size of the little common entry
962  first=jdahep(1,ip)
963  last =jdahep(2,ip)
964  npho=3+last-first+nhep-nhep0
965  nevpho=npho
966 C let-s take in decaying particle
967  idpho(1)=idhep(ip)
968  jdapho(1,1)=3
969  jdapho(2,1)=3+last-first
970  DO i=1,5
971  ppho(i,1)=phep(i,ip)
972  ENDDO
973 C let-s take in eventual second mother
974  ip2=jmohep(2,jdahep(1,ip))
975  IF((ip2.NE.0).AND.(ip2.NE.ip)) THEN
976  idpho(2)=idhep(ip2)
977  jdapho(1,2)=3
978  jdapho(2,2)=3+last-first
979  DO i=1,5
980  ppho(i,2)=phep(i,ip2)
981  ENDDO
982  ELSE
983  idpho(2)=0
984  DO i=1,5
985  ppho(i,2)=0.0d0
986  ENDDO
987  ENDIF
988 C let-s take in daughters
989  DO ll=0,last-first
990  idpho(3+ll)=idhep(first+ll)
991  jmopho(1,3+ll)=jmohep(1,first+ll)
992  IF (jmohep(1,first+ll).EQ.ip) jmopho(1,3+ll)=1
993  DO i=1,5
994  ppho(i,3+ll)=phep(i,first+ll)
995  ENDDO
996  ENDDO
997  IF (nhep.GT.nhep0) THEN
998 C let-s take in illegitimate daughters
999  na=3+last-first
1000  DO ll=1,nhep-nhep0
1001  idpho(na+ll)=idhep(nhep0+ll)
1002  jmopho(1,na+ll)=jmohep(1,nhep0+ll)
1003  IF (jmohep(1,nhep0+ll).EQ.ip) jmopho(1,na+ll)=1
1004  DO i=1,5
1005  ppho(i,na+ll)=phep(i,nhep0+ll)
1006  ENDDO
1007  ENDDO
1008 C-- there is NHEP-NHEP0 daugters more.
1009  jdapho(2,1)=3+last-first+nhep-nhep0
1010  ENDIF
1011  IF(idpho(npho).EQ.22)CALL phlupa(100001)
1012 ! IF(IDPHO(NPHO).EQ.22) stop
1013  CALL phcork(0)
1014  IF(idpho(npho).EQ.22)CALL phlupa(100002)
1015 C special case of t tbar production process
1016  IF(iftop) CALL photwo(0)
1017  boost=.false.
1018 C-- Check whether parent is in its rest frame...
1019 C ZBW ND 27.07.2009:
1020 C bug reported by Vladimir Savinov localized and fixed.
1021 C protection against rounding error was back-firing if soft
1022 C momentum of mother was physical. Consequence was that PHCORK was
1023 C messing up masses of final state particles in vertex of the decay.
1024 C Only configurations with previously generated photons of energy fraction
1025 C smaller than 0.0001 were affected. Effect was numerically insignificant.
1026 
1027 C IF ( (ABS(PPHO(4,1)-PPHO(5,1)).GT.PPHO(5,1)*1.D-8)
1028 C $ .AND.(PPHO(5,1).NE.0)) THEN
1029 
1030  IF ((abs(ppho(1,1)+abs(ppho(2,1))+abs(ppho(3,1))).GT.
1031  $ ppho(5,1)*1.d-8).AND.(ppho(5,1).NE.0)) THEN
1032 
1033  boost=.true.
1034 C--
1035 C-- Boost daughter particles to rest frame of parent...
1036 C-- Resultant neutral system already calculated in rest frame !
1037  DO 10 j=1,3
1038  10 bet(j)=-ppho(j,1)/ppho(5,1)
1039  gam=ppho(4,1)/ppho(5,1)
1040  DO 30 i=jdapho(1,1),jdapho(2,1)
1041  pb=bet(1)*ppho(1,i)+bet(2)*ppho(2,i)+bet(3)*ppho(3,i)
1042  DO 20 j=1,3
1043  20 ppho(j,i)=ppho(j,i)+bet(j)*(ppho(4,i)+pb/(gam+1.d0))
1044  30 ppho(4,i)=gam*ppho(4,i)+pb
1045 C-- Finally boost mother as well
1046  i=1
1047  pb=bet(1)*ppho(1,i)+bet(2)*ppho(2,i)+bet(3)*ppho(3,i)
1048  DO j=1,3
1049  ppho(j,i)=ppho(j,i)+bet(j)*(ppho(4,i)+pb/(gam+1.d0))
1050  ENDDO
1051  ppho(4,i)=gam*ppho(4,i)+pb
1052  ENDIF
1053 C special case of t tbar production process
1054  IF(iftop) CALL photwo(1)
1055  CALL phlupa(2)
1056  IF(idpho(npho).EQ.22) CALL phlupa(10000)
1057 ! IF(IDPHO(NPHO-1).EQ.22) stop
1058  END
1059  SUBROUTINE photwo(MODE)
1060 C.----------------------------------------------------------------------
1061 C.
1062 C. PHOTWO: PHOtos but TWO mothers allowed
1063 C.
1064 C. Purpose: Combines two mothers into one in /PHOEVT/
1065 C. necessary eg in case of g g (q qbar) --> t tbar
1066 C.
1067 C. Input Parameters: Common /PHOEVT/ (/PHOCMS/)
1068 C.
1069 C. Output Parameters: Common /PHOEVT/, (stored mothers)
1070 C.
1071 C. Author(s): Z. Was Created at: 5/08/93
1072 C. Last Update:10/08/93
1073 C.
1074 C.----------------------------------------------------------------------
1075  IMPLICIT NONE
1076  INTEGER nmxpho
1077  parameter(nmxpho=10000)
1078  INTEGER idpho,istpho,jdapho,jmopho,nevpho,npho
1079  REAL*8 ppho,vpho
1080  common/phoevt/nevpho,npho,istpho(nmxpho),idpho(nmxpho),
1081  &jmopho(2,nmxpho),jdapho(2,nmxpho),ppho(5,nmxpho),vpho(4,nmxpho)
1082  DOUBLE PRECISION bet(3),gam
1083  COMMON /phocms/ bet,gam
1084  INTEGER i,mode
1085  REAL*8 mpasqr
1086  LOGICAL ifrad
1087 C logical IFRAD is used to tag cases when two mothers may be
1088 C merged to the sole one.
1089 C So far used in case:
1090 C 1) of t tbar production
1091 C
1092 C t tbar case
1093  IF(mode.EQ.0) THEN
1094  ifrad=(idpho(1).EQ.21).AND.(idpho(2).EQ.21)
1095  ifrad=ifrad.OR.(idpho(1).EQ.-idpho(2).AND.abs(idpho(1)).LE.6)
1096  ifrad=ifrad
1097  & .AND.(abs(idpho(3)).EQ.6).AND.(abs(idpho(4)).EQ.6)
1098  mpasqr= (ppho(4,1)+ppho(4,2))**2-(ppho(3,1)+ppho(3,2))**2
1099  & -(ppho(2,1)+ppho(2,2))**2-(ppho(1,1)+ppho(1,2))**2
1100  ifrad=ifrad.AND.(mpasqr.GT.0.0d0)
1101  IF(ifrad) THEN
1102 c.....combining first and second mother
1103  DO i=1,4
1104  ppho(i,1)=ppho(i,1)+ppho(i,2)
1105  ENDDO
1106  ppho(5,1)=sqrt(mpasqr)
1107 c.....removing second mother,
1108  DO i=1,5
1109  ppho(i,2)=0.0d0
1110  ENDDO
1111  ENDIF
1112  ELSE
1113 C boosting of the mothers to the reaction frame not implemented yet.
1114 C to do it in mode 0 original mothers have to be stored in new comon (?)
1115 C and in mode 1 boosted to cms.
1116  ENDIF
1117  END
1118  SUBROUTINE phoout(IP,BOOST,NHEP0)
1119 C.----------------------------------------------------------------------
1120 C.
1121 C. PHOOUT: PHOtos OUTput
1122 C.
1123 C. Purpose: copies back IP branch of the common /PH_HEPEVT/ from
1124 C. /PHOEVT/ moves branch back from its CMS system.
1125 C.
1126 C. Input Parameters: IP: pointer of particle starting branch
1127 C. to be given back.
1128 C. BOOST: Flag whether boost to CMS was or was
1129 C . not performed.
1130 C.
1131 C. Output Parameters: Common /PHOEVT/,
1132 C.
1133 C. Author(s): Z. Was Created at: 24/05/93
1134 C. Last Update:
1135 C.
1136 C.----------------------------------------------------------------------
1137  IMPLICIT NONE
1138  INTEGER nmxhep
1139  parameter(nmxhep=10000)
1140  INTEGER idhep,isthep,jdahep,jmohep,nevhep,nhep
1141  REAL*8 phep,vhep
1142  common/ph_hepevt/nevhep,nhep,isthep(nmxhep),idhep(nmxhep),
1143  &jmohep(2,nmxhep),jdahep(2,nmxhep),phep(5,nmxhep),vhep(4,nmxhep)
1144  INTEGER nmxpho
1145  parameter(nmxpho=10000)
1146  INTEGER idpho,istpho,jdapho,jmopho,nevpho,npho
1147  REAL*8 ppho,vpho
1148  common/phoevt/nevpho,npho,istpho(nmxpho),idpho(nmxpho),
1149  &jmopho(2,nmxpho),jdapho(2,nmxpho),ppho(5,nmxpho),vpho(4,nmxpho)
1150  INTEGER ip,ll,first,last,i
1151  LOGICAL boost
1152  INTEGER nn,j,k,nhep0,na
1153  DOUBLE PRECISION bet(3),gam,pb
1154  COMMON /phocms/ bet,gam
1155  IF(npho.EQ.nevpho) RETURN
1156 C-- When parent was not in its rest-frame, boost back...
1157  CALL phlupa(10)
1158  IF (boost) THEN
1159  DO 110 j=jdapho(1,1),jdapho(2,1)
1160  pb=-bet(1)*ppho(1,j)-bet(2)*ppho(2,j)-bet(3)*ppho(3,j)
1161  DO 100 k=1,3
1162  100 ppho(k,j)=ppho(k,j)-bet(k)*(ppho(4,j)+pb/(gam+1.d0))
1163  110 ppho(4,j)=gam*ppho(4,j)+pb
1164 C-- ...boost photon, or whatever else has shown up
1165  DO nn=nevpho+1,npho
1166  pb=-bet(1)*ppho(1,nn)-bet(2)*ppho(2,nn)-bet(3)*ppho(3,nn)
1167  DO 120 k=1,3
1168  120 ppho(k,nn)=ppho(k,nn)-bet(k)*(ppho(4,nn)+pb/(gam+1.d0))
1169  ppho(4,nn)=gam*ppho(4,nn)+pb
1170  ENDDO
1171  ENDIF
1172  first=jdahep(1,ip)
1173  last =jdahep(2,ip)
1174 C let-s take in original daughters
1175  DO ll=0,last-first
1176  idhep(first+ll) = idpho(3+ll)
1177  DO i=1,5
1178  phep(i,first+ll) = ppho(i,3+ll)
1179  ENDDO
1180  ENDDO
1181 C let-s take newcomers to the end of HEPEVT.
1182  na=3+last-first
1183  DO ll=1,npho-na
1184  idhep(nhep0+ll) = idpho(na+ll)
1185  isthep(nhep0+ll)=istpho(na+ll)
1186  jmohep(1,nhep0+ll)=ip
1187  jmohep(2,nhep0+ll)=jmohep(2,jdahep(1,ip))
1188  jdahep(1,nhep0+ll)=0
1189  jdahep(2,nhep0+ll)=0
1190  DO i=1,5
1191  phep(i,nhep0+ll) = ppho(i,na+ll)
1192  ENDDO
1193  ENDDO
1194  nhep=nhep+npho-nevpho
1195  CALL phlupa(20)
1196  END
1197  SUBROUTINE phochk(JFIRST)
1198 C.----------------------------------------------------------------------
1199 C.
1200 C. PHOCHK: checking branch.
1201 C.
1202 C. Purpose: checks whether particles in the common block /PHOEVT/
1203 C. can be served by PHOMAK.
1204 C. JFIRST is the position in /PH_HEPEVT/ (!) of the first
1205 C. daughter of sub-branch under action.
1206 C.
1207 C.
1208 C. Author(s): Z. Was Created at: 22/10/92
1209 C. Last Update: 11/12/00
1210 C.
1211 C.----------------------------------------------------------------------
1212 C ********************
1213  IMPLICIT NONE
1214  INTEGER nmxpho
1215  parameter(nmxpho=10000)
1216  INTEGER idpho,istpho,jdapho,jmopho,nevpho,npho
1217  REAL*8 ppho,vpho
1218  common/phoevt/nevpho,npho,istpho(nmxpho),idpho(nmxpho),
1219  &jmopho(2,nmxpho),jdapho(2,nmxpho),ppho(5,nmxpho),vpho(4,nmxpho)
1220  LOGICAL chkif
1221  common/phoif/chkif(nmxpho)
1222  INTEGER nmxhep
1223  parameter(nmxhep=10000)
1224  LOGICAL qedrad
1225  common/ph_phoqed/qedrad(nmxhep)
1226  INTEGER jfirst
1227  LOGICAL f
1228  INTEGER idabs,nlast,i,ippar
1229  LOGICAL interf,isec,itre,iexp,iftop,ifw,ifnpi0,ifkl
1230  REAL*8 fint,fsec,expeps
1231  COMMON /phokey/ fsec,fint,expeps,interf,isec,itre,iexp,iftop,ifw
1232  LOGICAL ifrad
1233  INTEGER ident,k,iqrk,iphqrk,iekl,iphekl
1234 C these are OK .... if you do not like somebody else, add here.
1235  f(idabs)=
1236  & ( ((idabs.GT.9.OR.iqrk.NE.1).AND.(idabs.LE.40))
1237  & .OR.(idabs.GT.100) )
1238  & .AND.(idabs.NE.21)
1239  $ .AND.(idabs.NE.2101).AND.(idabs.NE.3101).AND.(idabs.NE.3201)
1240  & .AND.(idabs.NE.1103).AND.(idabs.NE.2103).AND.(idabs.NE.2203)
1241  & .AND.(idabs.NE.3103).AND.(idabs.NE.3203).AND.(idabs.NE.3303)
1242 C
1243  iqrk=iphqrk(0) ! switch for emission from quark
1244  iekl=iphekl(0)
1245  nlast = npho
1246 C
1247  ippar=1
1248 C checking for good particles
1249  ifnpi0=.true.
1250  IF (iekl.GT.1) THEN ! exclude radiative corr in decay of pi0
1251 C ! and Kl --> ee gamma
1252  ifnpi0= (idpho(1).NE.111) ! pi0
1253  ifkl = ((idpho(1).EQ.130).AND. ! Kl --> ee gamma
1254  $ ((idpho(3).EQ.22).OR.(idpho(4).EQ.22).OR.
1255  $ (idpho(5).EQ.22)).AND.
1256  $ ((idpho(3).EQ.11).OR.(idpho(4).EQ.11).OR.
1257  $ (idpho(5).EQ.11)) )
1258 
1259  ifnpi0=(ifnpi0.AND.(.NOT.ifkl))
1260  ENDIF
1261  DO 10 i=ippar,nlast
1262  idabs = abs(idpho(i))
1263 C possibly call on PHZODE is a dead (to be omitted) code.
1264  chkif(i)= f(idabs) .AND.f(abs(idpho(1)))
1265  & .AND. (idpho(2).EQ.0)
1266  IF(i.GT.2) chkif(i)=chkif(i).AND.qedrad(jfirst+i-ippar-2)
1267  & .AND.ifnpi0
1268  10 CONTINUE
1269 C--
1270 C now we go to special cases, where CHKIF(I) will be overwritten
1271 C--
1272  IF(iftop) THEN
1273 C special case of top pair production
1274  DO k=jdapho(2,1),jdapho(1,1),-1
1275  IF(idpho(k).NE.22) THEN
1276  ident=k
1277  goto 15
1278  ENDIF
1279  ENDDO
1280  15 CONTINUE
1281  ifrad=((idpho(1).EQ.21).AND.(idpho(2).EQ.21))
1282  & .OR. ((abs(idpho(1)).LE.6).AND.((idpho(2)).EQ.(-idpho(1))))
1283  ifrad=ifrad
1284  & .AND.(abs(idpho(3)).EQ.6).AND.((idpho(4)).EQ.(-idpho(3)))
1285  & .AND.(ident.EQ.4)
1286  IF(ifrad) THEN
1287  DO 20 i=ippar,nlast
1288  chkif(i)= .true.
1289  IF(i.GT.2) chkif(i)=chkif(i).AND.qedrad(jfirst+i-ippar-2)
1290  20 CONTINUE
1291  ENDIF
1292  ENDIF
1293 C--
1294 C--
1295  IF(iftop) THEN
1296 C special case of top decay
1297  DO k=jdapho(2,1),jdapho(1,1),-1
1298  IF(idpho(k).NE.22) THEN
1299  ident=k
1300  goto 25
1301  ENDIF
1302  ENDDO
1303  25 CONTINUE
1304  ifrad=((abs(idpho(1)).EQ.6).AND.(idpho(2).EQ.0))
1305  ifrad=ifrad
1306  & .AND.((abs(idpho(3)).EQ.24).AND.(abs(idpho(4)).EQ.5)
1307  & .OR.(abs(idpho(3)).EQ.5).AND.(abs(idpho(4)).EQ.24))
1308  & .AND.(ident.EQ.4)
1309  IF(ifrad) THEN
1310  DO 30 i=ippar,nlast
1311  chkif(i)= .true.
1312  IF(i.GT.2) chkif(i)=chkif(i).AND.qedrad(jfirst+i-ippar-2)
1313  30 CONTINUE
1314  ENDIF
1315  ENDIF
1316 C--
1317 C--
1318  END
1319  SUBROUTINE phtype(ID)
1320 C.----------------------------------------------------------------------
1321 C.
1322 C. PHTYPE: Central manadgement routine.
1323 C.
1324 C. Purpose: defines what kind of the
1325 C. actions will be performed at point ID.
1326 C.
1327 C. Input Parameters: ID: pointer of particle starting branch
1328 C. in /PH_HEPEVT/ to be treated.
1329 C.
1330 C. Output Parameters: Common /PH_HEPEVT/.
1331 C.
1332 C. Author(s): Z. Was Created at: 24/05/93
1333 C. P. Golonka Last Update: 27/06/04
1334 C.
1335 C.----------------------------------------------------------------------
1336  IMPLICIT NONE
1337  INTEGER nmxhep
1338  parameter(nmxhep=10000)
1339  INTEGER idhep,isthep,jdahep,jmohep,nevhep,nhep
1340  REAL*8 phep,vhep
1341  common/ph_hepevt/nevhep,nhep,isthep(nmxhep),idhep(nmxhep),
1342  &jmohep(2,nmxhep),jdahep(2,nmxhep),phep(5,nmxhep),vhep(4,nmxhep)
1343  LOGICAL interf,isec,itre,iexp,iftop,ifw
1344  REAL*8 fint,fsec,expeps
1345  COMMON /phokey/ fsec,fint,expeps,interf,isec,itre,iexp,iftop,ifw
1346  LOGICAL expini
1347  INTEGER nx,k,nchan
1348  parameter(nx=10)
1349  REAL*8 pro,prsum,esu
1350  COMMON /phoexp/ pro(nx),nchan,expini
1351 
1352  INTEGER id,nhep0
1353  LOGICAL ipair
1354  REAL*8 rn,phoran,sum
1355  INTEGER wtdum
1356  LOGICAL ifour
1357 C--
1358  ifour=(.true.).AND.(itre) ! we can make internal choice whether
1359  ! we want 3 or four photons at most.
1360  ipair=.true.
1361 C-- Check decay multiplicity..
1362  IF (jdahep(1,id).EQ.0) RETURN
1363 C IF (JDAHEP(1,ID).EQ.JDAHEP(2,ID)) RETURN
1364 C--
1365  nhep0=nhep
1366 C--
1367  IF (iexp) THEN
1368  expini=.true. ! Initialization/cleaning
1369  DO nchan=1,nx
1370  pro(nchan)=0.d0
1371  ENDDO
1372  nchan=0
1373 
1374  fsec=1.0d0
1375  CALL phomak(id,nhep0)! Initialization/crude formfactors into
1376  ! PRO(NCHAN)
1377  expini=.false.
1378  rn=phoran(wtdum)
1379  prsum=0
1380  DO k=1,nx
1381  prsum=prsum+pro(k)
1382  ENDDO
1383  esu=exp(-prsum) ! exponent for crude Poissonian multiplicity
1384  ! distribution, will be later overwritten
1385  ! to give probability for k
1386  sum=esu ! distribuant for the crude Poissonian
1387  ! at first for k=0
1388  DO k=1,100 ! hard coded max (photon) multiplicity is 100
1389  IF(rn.LT.sum) goto 100
1390  esu=esu*prsum/k ! we get at K ESU=EXP(-PRSUM)*PRSUM**K/K!
1391  sum=sum+esu ! thus we get distribuant at K.
1392  nchan=0
1393  CALL phomak(id,nhep0) ! LOOPING
1394  IF(sum.GT.1d0-expeps) goto 100
1395  ENDDO
1396  100 CONTINUE
1397  ELSEIF(ifour) THEN
1398 C-- quatro photon emission
1399  fsec=1.0d0
1400  rn=phoran(wtdum)
1401  IF (rn.GE.23.d0/24d0) THEN
1402  CALL phomak(id,nhep0)
1403  CALL phomak(id,nhep0)
1404  CALL phomak(id,nhep0)
1405  CALL phomak(id,nhep0)
1406  ELSEIF (rn.GE.17.d0/24d0) THEN
1407  CALL phomak(id,nhep0)
1408  CALL phomak(id,nhep0)
1409  ELSEIF (rn.GE.9.d0/24d0) THEN
1410  CALL phomak(id,nhep0)
1411  ENDIF
1412  ELSEIF(itre) THEN
1413 C-- triple photon emission
1414  fsec=1.0d0
1415  rn=phoran(wtdum)
1416  IF (rn.GE.5.d0/6d0) THEN
1417  CALL phomak(id,nhep0)
1418  CALL phomak(id,nhep0)
1419  CALL phomak(id,nhep0)
1420  ELSEIF (rn.GE.2.d0/6d0) THEN
1421  CALL phomak(id,nhep0)
1422  ENDIF
1423  ELSEIF(isec) THEN
1424 C-- double photon emission
1425  fsec=1.0d0
1426  rn=phoran(wtdum)
1427  IF (rn.GE.0.5d0) THEN
1428  CALL phomak(id,nhep0)
1429  CALL phomak(id,nhep0)
1430  ENDIF
1431  ELSE
1432 C-- single photon emission
1433  fsec=1.0d0
1434  CALL phomak(id,nhep0)
1435  ENDIF
1436 C--
1437 C-- electron positron pair (coomented out for a while
1438 C IF (IPAIR) CALL PHOPAR(ID,NHEP0)
1439  END
1440  SUBROUTINE phomak(IPPAR,NHEP0)
1441 C.----------------------------------------------------------------------
1442 C.
1443 C. PHOMAK: PHOtos MAKe
1444 C.
1445 C. Purpose: Single or double bremstrahlung radiative corrections
1446 C. are generated in the decay of the IPPAR-th particle in
1447 C. the HEP common /PH_HEPEVT/. Example of the use of
1448 C. general tools.
1449 C.
1450 C. Input Parameter: IPPAR: Pointer to decaying particle in
1451 C. /PH_HEPEVT/ and the common itself
1452 C.
1453 C. Output Parameters: Common /PH_HEPEVT/, either with or without
1454 C. particles added.
1455 C.
1456 C. Author(s): Z. Was, Created at: 26/05/93
1457 C. Last Update: 29/01/05
1458 C.
1459 C.----------------------------------------------------------------------
1460  IMPLICIT NONE
1461  DOUBLE PRECISION data
1462  REAL*8 phoran
1463  INTEGER ip,ippar,ncharg
1464  INTEGER wtdum,idum,nhep0
1465  INTEGER ncharb,neudau
1466  REAL*8 rn,wt,phint
1467  LOGICAL boost
1468  INTEGER nmxhep
1469  parameter(nmxhep=10000)
1470  INTEGER idhep,isthep,jdahep,jmohep,nevhep,nhep
1471  REAL*8 phep,vhep
1472  common/ph_hepevt/nevhep,nhep,isthep(nmxhep),idhep(nmxhep),
1473  &jmohep(2,nmxhep),jdahep(2,nmxhep),phep(5,nmxhep),vhep(4,nmxhep)
1474  LOGICAL interf,isec,itre,iexp,iftop,ifw
1475  REAL*8 fint,fsec,expeps
1476  COMMON /phokey/ fsec,fint,expeps,interf,isec,itre,iexp,iftop,ifw
1477 C--
1478  ip=ippar
1479  idum=1
1480  ncharg=0
1481 C--
1482  CALL phoin(ip,boost,nhep0)
1483  CALL phochk(jdahep(1,ip))
1484  wt=0.0d0
1485  CALL phopre(1,wt,neudau,ncharb)
1486 
1487  IF (wt.EQ.0.0d0) RETURN
1488  rn=phoran(wtdum)
1489 C PHODO is caling PHORAN, thus change of series if it is moved before if
1490  CALL phodo(1,ncharb,neudau)
1491 C we eliminate /FINT in variant B.
1492  IF (interf) wt=wt*phint(idum) /fint ! FINT must be in variant A
1493  IF (ifw) CALL phobw(wt) ! extra weight for leptonic W decay
1494  data=wt
1495  IF (wt.GT.1.0d0) CALL phoerr(3,'WT_INT',data)
1496 C weighting
1497  IF (rn.LE.wt) THEN
1498  CALL phoout(ip,boost,nhep0)
1499  ENDIF
1500  RETURN
1501  END
1502  FUNCTION phint1(IDUM)
1503 C.----------------------------------------------------------------------
1504 C.
1505 C. PHINT: PHotos INTerference (Old version kept for tests only.
1506 C.
1507 C. Purpose: Calculates interference between emission of photons from
1508 C. different possible chaged daughters stored in
1509 C. the HEP common /PHOEVT/.
1510 C.
1511 C. Input Parameter: commons /PHOEVT/ /PHOMOM/ /PHOPHS/
1512 C.
1513 C.
1514 C. Output Parameters:
1515 C.
1516 C.
1517 C. Author(s): Z. Was, Created at: 10/08/93
1518 C. Last Update: 15/03/99
1519 C.
1520 C.----------------------------------------------------------------------
1521 
1522  IMPLICIT NONE
1523  REAL*8 phint,phint1
1524  REAL*8 phocha
1525  INTEGER idum
1526  INTEGER nmxpho
1527  parameter(nmxpho=10000)
1528  INTEGER idpho,istpho,jdapho,jmopho,nevpho,npho
1529  REAL*8 ppho,vpho
1530  common/phoevt/nevpho,npho,istpho(nmxpho),idpho(nmxpho),
1531  &jmopho(2,nmxpho),jdapho(2,nmxpho),ppho(5,nmxpho),vpho(4,nmxpho)
1532  DOUBLE PRECISION mchsqr,mnesqr
1533  REAL*8 pneutr
1534  common/phomom/mchsqr,mnesqr,pneutr(5)
1535  DOUBLE PRECISION costhg,sinthg
1536  REAL*8 xphmax,xphoto
1537  common/phophs/xphmax,xphoto,costhg,sinthg
1538  REAL*8 mpasqr,xx,beta
1539  LOGICAL ifint
1540  INTEGER k,ident
1541 C
1542  DO k=jdapho(2,1),jdapho(1,1),-1
1543  IF(idpho(k).NE.22) THEN
1544  ident=k
1545  goto 20
1546  ENDIF
1547  ENDDO
1548  20 CONTINUE
1549 C check if there is a photon
1550  ifint= npho.GT.ident
1551 C check if it is two body + gammas reaction
1552  ifint= ifint.AND.(ident-jdapho(1,1)).EQ.1
1553 C check if two body was particle antiparticle
1554  ifint= ifint.AND.idpho(jdapho(1,1)).EQ.-idpho(ident)
1555 C check if particles were charged
1556  ifint= ifint.AND.phocha(idpho(ident)).NE.0
1557 C calculates interference weight contribution
1558  IF(ifint) THEN
1559  mpasqr = ppho(5,1)**2
1560  xx=4.d0*mchsqr/mpasqr*(1.d0-xphoto)/(1.d0-xphoto+(mchsqr-mnesqr)
1561  & /mpasqr)**2
1562  beta=sqrt(1.d0-xx)
1563  phint = 2d0/(1d0+costhg**2*beta**2)
1564  ELSE
1565  phint = 1d0
1566  ENDIF
1567  phint1=1
1568  END
1569 
1570  FUNCTION phint2(IDUM)
1571 C.----------------------------------------------------------------------
1572 C.
1573 C. PHINT: PHotos INTerference
1574 C.
1575 C. Purpose: Calculates interference between emission of photons from
1576 C. different possible chaged daughters stored in
1577 C. the HEP common /PHOEVT/.
1578 C.
1579 C. Input Parameter: commons /PHOEVT/ /PHOMOM/ /PHOPHS/
1580 C.
1581 C.
1582 C. Output Parameters:
1583 C.
1584 C.
1585 C. Author(s): Z. Was, Created at: 10/08/93
1586 C. Last Update:
1587 C.
1588 C.----------------------------------------------------------------------
1589 
1590  IMPLICIT NONE
1591  REAL*8 phint,phint1,phint2
1592  REAL*8 phocha
1593  INTEGER idum
1594  INTEGER nmxpho
1595  parameter(nmxpho=10000)
1596  INTEGER idpho,istpho,jdapho,jmopho,nevpho,npho
1597  REAL*8 ppho,vpho
1598  common/phoevt/nevpho,npho,istpho(nmxpho),idpho(nmxpho),
1599  &jmopho(2,nmxpho),jdapho(2,nmxpho),ppho(5,nmxpho),vpho(4,nmxpho)
1600  DOUBLE PRECISION mchsqr,mnesqr
1601  REAL*8 pneutr
1602  common/phomom/mchsqr,mnesqr,pneutr(5)
1603  DOUBLE PRECISION costhg,sinthg
1604  REAL*8 xphmax,xphoto
1605  common/phophs/xphmax,xphoto,costhg,sinthg
1606  REAL*8 mpasqr,xx,beta,pq1(4),pq2(4),pphot(4)
1607  REAL*8 ss,pp2,pp,e1,e2,q1,q2,costhe
1608  LOGICAL ifint
1609  INTEGER k,ident
1610 C
1611  DO k=jdapho(2,1),jdapho(1,1),-1
1612  IF(idpho(k).NE.22) THEN
1613  ident=k
1614  goto 20
1615  ENDIF
1616  ENDDO
1617  20 CONTINUE
1618 C check if there is a photon
1619  ifint= npho.GT.ident
1620 C check if it is two body + gammas reaction
1621  ifint= ifint.AND.(ident-jdapho(1,1)).EQ.1
1622 C check if two body was particle antiparticle (we improve on it !
1623 C IFINT= IFINT.AND.IDPHO(JDAPHO(1,1)).EQ.-IDPHO(IDENT)
1624 C check if particles were charged
1625  ifint= ifint.AND.abs(phocha(idpho(ident))).GT.0.01d0
1626 C check if they have both charge
1627  ifint= ifint.AND.
1628  $ abs(phocha(idpho(jdapho(1,1)))).gt.0.01d0
1629 C calculates interference weight contribution
1630  IF(ifint) THEN
1631  mpasqr = ppho(5,1)**2
1632  xx=4.d0*mchsqr/mpasqr*(1.-xphoto)/(1.-xphoto+(mchsqr-mnesqr)/
1633  & mpasqr)**2
1634  beta=sqrt(1.d0-xx)
1635  phint = 2d0/(1d0+costhg**2*beta**2)
1636  ss =mpasqr*(1.d0-xphoto)
1637  pp2=((ss-mchsqr-mnesqr)**2-4*mchsqr*mnesqr)/ss/4
1638  pp =sqrt(pp2)
1639  e1 =sqrt(pp2+mchsqr)
1640  e2 =sqrt(pp2+mnesqr)
1641  phint= (e1+e2)**2/((e2+costhg*pp)**2+(e1-costhg*pp)**2)
1642 C
1643  q1=phocha(idpho(jdapho(1,1)))
1644  q2=phocha(idpho(ident))
1645  do k=1,4
1646  pq1(k)=ppho(k,jdapho(1,1))
1647  pq2(k)=ppho(k,jdapho(1,1)+1)
1648  pphot(k)=ppho(k,npho)
1649  enddo
1650  costhe=(pphot(1)*pq1(1)+pphot(2)*pq1(2)+pphot(3)*pq1(3))
1651  costhe=costhe/sqrt(pq1(1)**2+pq1(2)**2+pq1(3)**2)
1652  costhe=costhe/sqrt(pphot(1)**2+pphot(2)**2+pphot(3)**2)
1653 C
1654 ! --- this IF checks whether JDAPHO(1,1) was MCH or MNE.
1655 ! --- COSTHG angle (and in-generation variables) may be better choice
1656 ! --- than costhe. note that in the formulae below amplitudes were
1657 ! --- multiplied by (E2+COSTHG*PP)*(E1-COSTHG*PP).
1658  IF (costhg*costhe.GT.0) then
1659 
1660  phint= (q1*(e2+costhg*pp)-q2*(e1-costhg*pp))**2
1661  & /(q1**2*(e2+costhg*pp)**2+q2**2*(e1-costhg*pp)**2)
1662  ELSE
1663 
1664  phint= (q1*(e1-costhg*pp)-q2*(e2+costhg*pp))**2
1665  & /(q1**2*(e1-costhg*pp)**2+q2**2*(e2+costhg*pp)**2)
1666  ENDIF
1667  ELSE
1668  phint = 1d0
1669  ENDIF
1670  phint1=1
1671  phint2=1
1672  END
1673 
1674 
1675  SUBROUTINE phopre(IPARR,WT,NEUDAU,NCHARB)
1676 C.----------------------------------------------------------------------
1677 C.
1678 C. PHOTOS: Photon radiation in decays
1679 C.
1680 C. Purpose: Order (alpha) radiative corrections are generated in
1681 C. the decay of the IPPAR-th particle in the HEP-like
1682 C. common /PHOEVT/. Photon radiation takes place from one
1683 C. of the charged daughters of the decaying particle IPPAR
1684 C. WT is calculated, eventual rejection will be performed
1685 C. later after inclusion of interference weight.
1686 C.
1687 C. Input Parameter: IPPAR: Pointer to decaying particle in
1688 C. /PHOEVT/ and the common itself,
1689 C.
1690 C. Output Parameters: Common /PHOEVT/, either with or without a
1691 C. photon(s) added.
1692 C. WT weight of the configuration
1693 C.
1694 C. Author(s): Z. Was, B. van Eijk Created at: 26/11/89
1695 C. Last Update: 29/01/05
1696 C.
1697 C.----------------------------------------------------------------------
1698  IMPLICIT NONE
1699  DOUBLE PRECISION minmas,mpasqr,mchren
1700  DOUBLE PRECISION beta,eps,del1,del2,data,biglog
1701  REAL*8 phocha,phospi,phoran,phocor,massum
1702  INTEGER ip,iparr,ippar,i,j,me,ncharg,neupoi,nlast,thedum
1703  INTEGER idabs,idum
1704  INTEGER ncharb,neudau
1705  REAL*8 wt,wgt
1706  INTEGER nmxpho
1707  parameter(nmxpho=10000)
1708  INTEGER idpho,istpho,jdapho,jmopho,nevpho,npho
1709  REAL*8 ppho,vpho
1710  common/phoevt/nevpho,npho,istpho(nmxpho),idpho(nmxpho),
1711  &jmopho(2,nmxpho),jdapho(2,nmxpho),ppho(5,nmxpho),vpho(4,nmxpho)
1712  LOGICAL chkif
1713  common/phoif/chkif(nmxpho)
1714  INTEGER chapoi(nmxpho)
1715  DOUBLE PRECISION mchsqr,mnesqr
1716  REAL*8 pneutr
1717  common/phomom/mchsqr,mnesqr,pneutr(5)
1718  DOUBLE PRECISION costhg,sinthg
1719  REAL*8 xphmax,xphoto
1720  common/phophs/xphmax,xphoto,costhg,sinthg
1721  REAL*8 alpha,xphcut
1722  common/phocop/alpha,xphcut
1723  INTEGER irep
1724  REAL*8 probh,corwt,xf
1725  common/phopro/probh,corwt,xf,irep
1726 C may be it is not the best place, but ...
1727  LOGICAL interf,isec,itre,iexp,iftop,ifw
1728  REAL*8 fint,fsec,expeps
1729  COMMON /phokey/ fsec,fint,expeps,interf,isec,itre,iexp,iftop,ifw
1730 
1731 C--
1732  ippar=iparr
1733 C-- Store pointers for cascade treatement...
1734  ip=ippar
1735  nlast=npho
1736  idum=1
1737 C--
1738 C-- Check decay multiplicity..
1739  IF (jdapho(1,ip).EQ.0) RETURN
1740 C--
1741 C-- Loop over daughters, determine charge multiplicity
1742  10 ncharg=0
1743  irep=0
1744  minmas=0.d0
1745  massum=0.d0
1746  DO 20 i=jdapho(1,ip),jdapho(2,ip)
1747 C--
1748 C--
1749 C-- Exclude marked particles, quarks and gluons etc...
1750  idabs=abs(idpho(i))
1751  IF (chkif(i-jdapho(1,ip)+3)) THEN
1752  IF (phocha(idpho(i)).NE.0) THEN
1753  ncharg=ncharg+1
1754  IF (ncharg.GT.nmxpho) THEN
1755  data=ncharg
1756  CALL phoerr(1,'PHOTOS',data)
1757  ENDIF
1758  chapoi(ncharg)=i
1759  ENDIF
1760  minmas=minmas+ppho(5,i)**2
1761  ENDIF
1762  massum=massum+ppho(5,i)
1763  20 CONTINUE
1764  IF (ncharg.NE.0) THEN
1765 C--
1766 C-- Check that sum of daughter masses does not exceed parent mass
1767  IF ((ppho(5,ip)-massum)/ppho(5,ip).GT.2.d0*xphcut) THEN
1768 C--
1769 C-- Order charged particles according to decreasing mass, this to
1770 C-- increase efficiency (smallest mass is treated first).
1771  IF (ncharg.GT.1) CALL phooma(1,ncharg,chapoi)
1772 C--
1773  30 CONTINUE
1774  DO 70 j=1,3
1775  70 pneutr(j)=-ppho(j,chapoi(ncharg))
1776  pneutr(4)=ppho(5,ip)-ppho(4,chapoi(ncharg))
1777 C--
1778 C-- Calculate invariant mass of 'neutral' etc. systems
1779  mpasqr=ppho(5,ip)**2
1780  mchsqr=ppho(5,chapoi(ncharg))**2
1781  IF ((jdapho(2,ip)-jdapho(1,ip)).EQ.1) THEN
1782  neupoi=jdapho(1,ip)
1783  IF (neupoi.EQ.chapoi(ncharg)) neupoi=jdapho(2,ip)
1784  mnesqr=ppho(5,neupoi)**2
1785  pneutr(5)=ppho(5,neupoi)
1786  ELSE
1787  mnesqr=pneutr(4)**2-pneutr(1)**2-pneutr(2)**2-pneutr(3)**2
1788  mnesqr=max(mnesqr,minmas-mchsqr)
1789  pneutr(5)=sqrt(mnesqr)
1790  ENDIF
1791 C--
1792 C-- Determine kinematical limit...
1793  xphmax=(mpasqr-(pneutr(5)+ppho(5,chapoi(ncharg)))**2)/mpasqr
1794 C--
1795 C-- Photon energy fraction...
1796  CALL phoene(mpasqr,mchren,beta,biglog,idpho(chapoi(ncharg)))
1797 C--
1798  IF (xphoto.LT.-4d0) THEN
1799  ncharg=0 ! we really stop trials
1800  xphoto=0d0! in this case !!
1801 C-- Energy fraction not too large (very seldom) ? Define angle.
1802  ELSEIF ((xphoto.LT.xphcut).OR.(xphoto.GT.xphmax)) THEN
1803 C--
1804 C-- No radiation was accepted, check for more daughters that may ra-
1805 C-- diate and correct radiation probability...
1806  ncharg=ncharg-1
1807  IF (ncharg.GT.0) THEN
1808  irep=irep+1
1809  goto 30
1810  ENDIF
1811  ELSE
1812 C--
1813 C-- Angle is generated in the frame defined by charged vector and
1814 C-- PNEUTR, distribution is taken in the infrared limit...
1815  eps=mchren/(1.d0+beta)
1816 C--
1817 C-- Calculate sin(theta) and cos(theta) from interval variables
1818  del1=(2.d0-eps)*(eps/(2.d0-eps))**phoran(thedum)
1819  del2=2.d0-del1
1820 
1821 C ----------- VARIANT B ------------------
1822 CC corrections for more efiicient interference correction,
1823 CC instead of doubling crude distribution, we add flat parallel channel
1824 C IF (PHORAN(THEDUM).LT.BIGLOG/BETA/(BIGLOG/BETA+2*FINT)) THEN
1825 C COSTHG=(1.D0-DEL1)/BETA
1826 C SINTHG=SQRT(DEL1*DEL2-MCHREN)/BETA
1827 C ELSE
1828 C COSTHG=-1D0+2*PHORAN(THEDUM)
1829 C SINTHG= SQRT(1D0-COSTHG**2)
1830 C ENDIF
1831 C
1832 C IF (FINT.GT.1.0D0) THEN
1833 C
1834 C WGT=1D0/(1D0-BETA*COSTHG)
1835 C WGT=WGT/(WGT+FINT)
1836 C ! WGT=1D0 ! ??
1837 C
1838 C ELSE
1839 C WGT=1D0
1840 C ENDIF
1841 C
1842 C ----------- END OF VARIANT B ------------------
1843 
1844 C ----------- VARIANT A ------------------
1845  costhg=(1.d0-del1)/beta
1846  sinthg=sqrt(del1*del2-mchren)/beta
1847  wgt=1d0
1848 C ----------- END OF VARIANT A ------------------
1849 
1850 C--
1851 C-- Determine spin of particle and construct code for matrix element
1852  me=2.d0*phospi(idpho(chapoi(ncharg)))+1.d0
1853 C--
1854 C-- Weighting procedure with 'exact' matrix element, reconstruct kine-
1855 C-- matics for photon, neutral and charged system and update /PHOEVT/.
1856 C-- Find pointer to the first component of 'neutral' system
1857  DO i=jdapho(1,ip),jdapho(2,ip)
1858  IF (i.NE.chapoi(ncharg)) THEN
1859  neudau=i
1860  goto 51
1861  ENDIF
1862  ENDDO
1863 C--
1864 C-- Pointer not found...
1865  data=ncharg
1866  CALL phoerr(5,'PHOKIN',data)
1867  51 CONTINUE
1868  ncharb=chapoi(ncharg)
1869  ncharb=ncharb-jdapho(1,ip)+3
1870  neudau=neudau-jdapho(1,ip)+3
1871  wt=phocor(mpasqr,mchren,me)*wgt
1872  ENDIF
1873  ELSE
1874  data=ppho(5,ip)-massum
1875  CALL phoerr(10,'PHOTOS',data)
1876  ENDIF
1877  ENDIF
1878 C--
1879  RETURN
1880  END
1881  SUBROUTINE phooma(IFIRST,ILAST,POINTR)
1882 C.----------------------------------------------------------------------
1883 C.
1884 C. PHOTOS: PHOton radiation in decays Order MAss vector
1885 C.
1886 C. Purpose: Order the contents of array 'POINTR' according to the
1887 C. decreasing value in the array 'MASS'.
1888 C.
1889 C. Input Parameters: IFIRST, ILAST: Pointers to the vector loca-
1890 C. tion be sorted,
1891 C. POINTR: Unsorted array with pointers to
1892 C. /PHOEVT/.
1893 C.
1894 C. Output Parameter: POINTR: Sorted arrays with respect to
1895 C. particle mass 'PPHO(5,*)'.
1896 C.
1897 C. Author(s): B. van Eijk Created at: 28/11/89
1898 C. Last Update: 27/05/93
1899 C.
1900 C.----------------------------------------------------------------------
1901  IMPLICIT NONE
1902  INTEGER nmxpho
1903  parameter(nmxpho=10000)
1904  INTEGER idpho,istpho,jdapho,jmopho,nevpho,npho
1905  REAL*8 ppho,vpho
1906  common/phoevt/nevpho,npho,istpho(nmxpho),idpho(nmxpho),
1907  &jmopho(2,nmxpho),jdapho(2,nmxpho),ppho(5,nmxpho),vpho(4,nmxpho)
1908  INTEGER ifirst,ilast,i,j,bufpoi,pointr(nmxpho)
1909  REAL*8 bufmas,mass(nmxpho)
1910  IF (ifirst.EQ.ilast) RETURN
1911 C--
1912 C-- Copy particle masses
1913  DO 10 i=ifirst,ilast
1914  10 mass(i)=ppho(5,pointr(i))
1915 C--
1916 C-- Order the masses in a decreasing series
1917  DO 30 i=ifirst,ilast-1
1918  DO 20 j=i+1,ilast
1919  IF (mass(j).LE.mass(i)) goto 20
1920  bufpoi=pointr(j)
1921  pointr(j)=pointr(i)
1922  pointr(i)=bufpoi
1923  bufmas=mass(j)
1924  mass(j)=mass(i)
1925  mass(i)=bufmas
1926  20 CONTINUE
1927  30 CONTINUE
1928  RETURN
1929  END
1930  SUBROUTINE phoene(MPASQR,MCHREN,BETA,BIGLOG,IDENT)
1931 C.----------------------------------------------------------------------
1932 C.
1933 C. PHOTOS: PHOton radiation in decays calculation of photon ENErgy
1934 C. fraction
1935 C.
1936 C. Purpose: Subroutine returns photon energy fraction (in (parent
1937 C. mass)/2 units) for the decay bremsstrahlung.
1938 C.
1939 C. Input Parameters: MPASQR: Mass of decaying system squared,
1940 C. XPHCUT: Minimum energy fraction of photon,
1941 C. XPHMAX: Maximum energy fraction of photon.
1942 C.
1943 C. Output Parameter: MCHREN: Renormalised mass squared,
1944 C. BETA: Beta factor due to renormalisation,
1945 C. XPHOTO: Photon energy fraction,
1946 C. XF: Correction factor for PHOFAC.
1947 C.
1948 C. Author(s): S. Jadach, Z. Was Created at: 01/01/89
1949 C. B. van Eijk, P.Golonka Last Update: 29/01/05
1950 C.
1951 C.----------------------------------------------------------------------
1952  IMPLICIT NONE
1953  DOUBLE PRECISION mpasqr,mchren,biglog,beta,data
1954  INTEGER iwt1,irn,iwt2
1955  REAL*8 prsoft,prhard,phoran,phofac
1956  DOUBLE PRECISION mchsqr,mnesqr
1957  REAL*8 pneutr
1958  INTEGER ident
1959  REAL*8 phocha,prkill,rrr
1960  common/phomom/mchsqr,mnesqr,pneutr(5)
1961  DOUBLE PRECISION costhg,sinthg
1962  REAL*8 xphmax,xphoto
1963  common/phophs/xphmax,xphoto,costhg,sinthg
1964  REAL*8 alpha,xphcut
1965  common/phocop/alpha,xphcut
1966  REAL*8 pi,twopi
1967  common/phpico/pi,twopi
1968  INTEGER irep
1969  REAL*8 probh,corwt,xf
1970  common/phopro/probh,corwt,xf,irep
1971  LOGICAL interf,isec,itre,iexp,iftop,ifw
1972  REAL*8 fint,fsec,expeps
1973  COMMON /phokey/ fsec,fint,expeps,interf,isec,itre,iexp,iftop,ifw
1974  INTEGER nx,nchan,k
1975  parameter(nx=10)
1976  LOGICAL expini
1977  REAL*8 pro,prsum
1978  COMMON /phoexp/ pro(nx),nchan,expini
1979 C--
1980  IF (xphmax.LE.xphcut) THEN
1981  beta=phofac(-1) ! to zero counter, here beta is dummy
1982  xphoto=0.0d0
1983  RETURN
1984  ENDIF
1985 C-- Probabilities for hard and soft bremstrahlung...
1986  mchren=4.d0*mchsqr/mpasqr/(1.d0+mchsqr/mpasqr)**2
1987  beta=sqrt(1.d0-mchren)
1988 
1989 C ----------- VARIANT B ------------------
1990 CC we replace 1D0/BETA*BIGLOG with (1D0/BETA*BIGLOG+2*FINT)
1991 CC for integral of new crude
1992 C BIGLOG=LOG(MPASQR/MCHSQR*(1.D0+BETA)**2/4.D0*
1993 C & (1.D0+MCHSQR/MPASQR)**2)
1994 C PRHARD=ALPHA/PI*(1D0/BETA*BIGLOG+2*FINT)*(LOG(XPHMAX/XPHCUT)
1995 C &-.75D0+XPHCUT/XPHMAX-.25D0*XPHCUT**2/XPHMAX**2)
1996 C PRHARD=PRHARD*PHOCHA(IDENT)**2*FSEC
1997 C ----------- END OF VARIANT B ------------------
1998 
1999 C ----------- VARIANT A ------------------
2000  biglog=log(mpasqr/mchsqr*(1.d0+beta)**2/4.d0*
2001  & (1.d0+mchsqr/mpasqr)**2)
2002  prhard=alpha/pi*(1d0/beta*biglog)*
2003  &(log(xphmax/xphcut)-.75d0+xphcut/xphmax-.25d0*xphcut**2/xphmax**2)
2004  prhard=prhard*phocha(ident)**2*fsec*fint
2005 C ----------- END OF VARIANT A ------------------
2006  IF (irep.EQ.0) probh=0.d0
2007  prkill=0d0
2008  IF (iexp) THEN ! IEXP
2009  nchan=nchan+1
2010  IF (expini) THEN ! EXPINI
2011  pro(nchan)=prhard+0.05*(1.0+fint) ! we store hard photon emission prob
2012  !for leg NCHAN
2013  prhard=0d0 ! to kill emission at initialization call
2014  probh=prhard
2015  ELSE ! EXPINI
2016  prsum=0
2017  DO k=nchan,nx
2018  prsum=prsum+pro(k)
2019  ENDDO
2020  prhard=prhard/prsum ! note that PRHARD may be smaller than
2021  !PRO(NCHAN) because it is calculated
2022  ! for kinematical configuartion as is
2023  ! (with effects of previous photons)
2024  prkill=pro(nchan)/prsum-prhard !
2025 
2026  ENDIF ! EXPINI
2027  prsoft=1.d0-prhard
2028  ELSE ! IEXP
2029  prhard=prhard*phofac(0) ! PHOFAC is used to control eikonal
2030  ! formfactors for non exp version only
2031  ! here PHOFAC(0)=1 at least now.
2032  probh=prhard
2033  ENDIF ! IEXP
2034  prsoft=1.d0-prhard
2035 C--
2036 C-- Check on kinematical bounds
2037  IF (iexp) THEN
2038  IF (prsoft.LT.-5.0d-8) THEN
2039  data=prsoft
2040  CALL phoerr(2,'PHOENE',data)
2041  ENDIF
2042  ELSE
2043  IF (prsoft.LT.0.1d0) THEN
2044  data=prsoft
2045  CALL phoerr(2,'PHOENE',data)
2046  ENDIF
2047  ENDIF
2048 
2049  rrr=phoran(iwt1)
2050  IF (rrr.LT.prsoft) THEN
2051 C--
2052 C-- No photon... (ie. photon too soft)
2053  xphoto=0.d0
2054  IF (rrr.LT.prkill) xphoto=-5d0 ! No photon...no further trials
2055  ELSE
2056 C--
2057 C-- Hard photon... (ie. photon hard enough).
2058 C-- Calculate Altarelli-Parisi Kernel
2059  10 xphoto=exp(phoran(irn)*log(xphcut/xphmax))
2060  xphoto=xphoto*xphmax
2061  IF (phoran(iwt2).GT.((1.d0+(1.d0-xphoto/xphmax)**2)/2.d0))
2062  & goto 10
2063  ENDIF
2064 C--
2065 C-- Calculate parameter for PHOFAC function
2066  xf=4.d0*mchsqr*mpasqr/(mpasqr+mchsqr-mnesqr)**2
2067  RETURN
2068  END
2069  FUNCTION phocor(MPASQR,MCHREN,ME)
2070 C.----------------------------------------------------------------------
2071 C.
2072 C. PHOTOS: PHOton radiation in decays CORrection weight from
2073 C. matrix elements
2074 C.
2075 C. Purpose: Calculate photon angle. The reshaping functions will
2076 C. have to depend on the spin S of the charged particle.
2077 C. We define: ME = 2 * S + 1 !
2078 C.
2079 C. Input Parameters: MPASQR: Parent mass squared,
2080 C. MCHREN: Renormalised mass of charged system,
2081 C. ME: 2 * spin + 1 determines matrix element
2082 C.
2083 C. Output Parameter: Function value.
2084 C.
2085 C. Author(s): Z. Was, B. van Eijk Created at: 26/11/89
2086 C. Last Update: 21/03/93
2087 C.
2088 C.----------------------------------------------------------------------
2089  IMPLICIT NONE
2090  DOUBLE PRECISION mpasqr,mchren,beta,xx,yy,data
2091  INTEGER me
2092  REAL*8 phocor,phofac,wt1,wt2,wt3
2093  DOUBLE PRECISION mchsqr,mnesqr
2094  REAL*8 pneutr
2095  common/phomom/mchsqr,mnesqr,pneutr(5)
2096  DOUBLE PRECISION costhg,sinthg
2097  REAL*8 xphmax,xphoto
2098  common/phophs/xphmax,xphoto,costhg,sinthg
2099  INTEGER irep
2100  REAL*8 probh,corwt,xf
2101  common/phopro/probh,corwt,xf,irep
2102 C--
2103 C-- Shaping (modified by ZW)...
2104  xx=4.d0*mchsqr/mpasqr*(1.d0-xphoto)/(1.d0-xphoto+(mchsqr-mnesqr)/
2105  &mpasqr)**2
2106  IF (me.EQ.1) THEN
2107  yy=1.d0
2108  wt3=(1.d0-xphoto/xphmax)/((1.d0+(1.d0-xphoto/xphmax)**2)/2.d0)
2109  ELSEIF (me.EQ.2) THEN
2110  yy=0.5d0*(1.d0-xphoto/xphmax+1.d0/(1.d0-xphoto/xphmax))
2111  wt3=1.d0
2112  ELSEIF ((me.EQ.3).OR.(me.EQ.4).OR.(me.EQ.5)) THEN
2113  yy=1.d0
2114  wt3=(1.d0+(1.d0-xphoto/xphmax)**2-(xphoto/xphmax)**3)/
2115  & (1.d0+(1.d0-xphoto/xphmax)** 2)
2116  ELSE
2117  data=(me-1.d0)/2.d0
2118  CALL phoerr(6,'PHOCOR',data)
2119  yy=1.d0
2120  wt3=1.d0
2121  ENDIF
2122  beta=sqrt(1.d0-xx)
2123  wt1=(1.d0-costhg*sqrt(1.d0-mchren))/(1.d0-costhg*beta)
2124  wt2=(1.d0-xx/yy/(1.d0-beta**2*costhg**2))*(1.d0+costhg*beta)/2.d0
2125  wt2=wt2*phofac(1)
2126  phocor=wt1*wt2*wt3
2127  corwt=phocor
2128  IF (phocor.GT.1.d0) THEN
2129  data=phocor
2130  CALL phoerr(3,'PHOCOR',data)
2131  ENDIF
2132  RETURN
2133  END
2134  FUNCTION phofac(MODE)
2135 C.----------------------------------------------------------------------
2136 C.
2137 C. PHOTOS: PHOton radiation in decays control FACtor
2138 C.
2139 C. Purpose: This is the control function for the photon spectrum and
2140 C. final weighting. It is called from PHOENE for genera-
2141 C. ting the raw photon energy spectrum (MODE=0) and in PHO-
2142 C. COR to scale the final weight (MODE=1). The factor con-
2143 C. sists of 3 terms. Addition of the factor FF which mul-
2144 C. tiplies PHOFAC for MODE=0 and divides PHOFAC for MODE=1,
2145 C. does not affect the results for the MC generation. An
2146 C. appropriate choice for FF can speed up the calculation.
2147 C. Note that a too small value of FF may cause weight over-
2148 C. flow in PHOCOR and will generate a warning, halting the
2149 C. execution. PRX should be included for repeated calls
2150 C. for the same event, allowing more particles to radiate
2151 C. photons. At the first call IREP=0, for more than 1
2152 C. charged decay products, IREP >= 1. Thus, PRSOFT (no
2153 C. photon radiation probability in the previous calls)
2154 C. appropriately scales the strength of the bremsstrahlung.
2155 C.
2156 C. Input Parameters: MODE, PROBH, XF
2157 C.
2158 C. Output Parameter: Function value
2159 C.
2160 C. Author(s): S. Jadach, Z. Was Created at: 01/01/89
2161 C. B. van Eijk, P.Golonka Last Update: 26/06/04
2162 C.
2163 C.----------------------------------------------------------------------
2164  IMPLICIT NONE
2165  REAL*8 phofac,ff,prx
2166  INTEGER mode
2167  INTEGER irep
2168  REAL*8 probh,corwt,xf
2169  common/phopro/probh,corwt,xf,irep
2170  LOGICAL interf,isec,itre,iexp,iftop,ifw
2171  REAL*8 fint,fsec,expeps
2172  COMMON /phokey/ fsec,fint,expeps,interf,isec,itre,iexp,iftop,ifw
2173  SAVE prx,ff
2174  DATA prx,ff/ 0.d0, 0.d0/
2175  IF (iexp) THEN ! In case of exponentiation this routine is useles
2176  phofac=1
2177  RETURN
2178  ENDIF
2179  IF (mode.EQ.-1) THEN
2180  prx=1.d0
2181  ff=1.d0
2182  probh=0.0
2183  ELSEIF (mode.EQ.0) THEN
2184  IF (irep.EQ.0) prx=1.d0
2185  prx=prx/(1.d0-probh)
2186  ff=1.d0
2187 C--
2188 C-- Following options are not considered for the time being...
2189 C-- (1) Good choice, but does not save very much time:
2190 C-- FF=(1.0D0-SQRT(XF)/2.0D0)/(1.0+SQRT(XF)/2.0D0)
2191 C-- (2) Taken from the blue, but works without weight overflows...
2192 C-- FF=(1.D0-XF/(1-(1-SQRT(XF))**2))*(1+(1-SQRT(XF))/SQRT(1-XF))/2
2193  phofac=ff*prx
2194  ELSE
2195  phofac=1.d0/ff
2196  ENDIF
2197  END
2198  SUBROUTINE phobw(WT)
2199 C.----------------------------------------------------------------------
2200 C.
2201 C. PHOTOS: PHOtos Boson W correction weight
2202 C.
2203 C. Purpose: calculates correction weight due to amplitudes of
2204 C. emission from W boson.
2205 C.
2206 C.
2207 C.
2208 C.
2209 C.
2210 C. Input Parameters: Common /PHOEVT/, with photon added.
2211 C. wt to be corrected
2212 C.
2213 C.
2214 C.
2215 C. Output Parameters: wt
2216 C.
2217 C. Author(s): G. Nanava, Z. Was Created at: 13/03/03
2218 C. Last Update: 13/03/03
2219 C.
2220 C.----------------------------------------------------------------------
2221  IMPLICIT NONE
2222  DOUBLE PRECISION wt
2223  INTEGER nmxpho
2224  parameter(nmxpho=10000)
2225  INTEGER idpho,istpho,jdapho,jmopho,nevpho,npho
2226  REAL*8 ppho,vpho
2227  common/phoevt/nevpho,npho,istpho(nmxpho),idpho(nmxpho),
2228  &jmopho(2,nmxpho),jdapho(2,nmxpho),ppho(5,nmxpho),vpho(4,nmxpho)
2229  INTEGER i
2230  DOUBLE PRECISION emu,mchren,beta,costhg,mpasqr,xph
2231 C--
2232  IF (abs(idpho(1)).EQ.24.AND.
2233  $ abs(idpho(jdapho(1,1) )).GE.11.AND.
2234  $ abs(idpho(jdapho(1,1) )).LE.16.AND.
2235  $ abs(idpho(jdapho(1,1)+1)).GE.11.AND.
2236  $ abs(idpho(jdapho(1,1)+1)).LE.16 ) THEN
2237 
2238  IF(
2239  $ abs(idpho(jdapho(1,1) )).EQ.11.OR.
2240  $ abs(idpho(jdapho(1,1) )).EQ.13.OR.
2241  $ abs(idpho(jdapho(1,1) )).EQ.15 ) THEN
2242  i=jdapho(1,1)
2243  ELSE
2244  i=jdapho(1,1)+1
2245  ENDIF
2246  emu=ppho(4,i)
2247  mchren=abs(ppho(4,i)**2-ppho(3,i)**2
2248  $ -ppho(2,i)**2-ppho(1,i)**2)
2249  beta=sqrt(1- mchren/ ppho(4,i)**2)
2250  costhg=(ppho(3,i)*ppho(3,npho)+ppho(2,i)*ppho(2,npho)
2251  $ +ppho(1,i)*ppho(1,npho))/
2252  $ sqrt(ppho(3,i)**2+ppho(2,i)**2+ppho(1,i)**2) /
2253  $ sqrt(ppho(3,npho)**2+ppho(2,npho)**2+ppho(1,npho)**2)
2254  mpasqr=ppho(4,1)**2
2255  xph=ppho(4,npho)
2256  wt=wt*(1-8*emu*xph*(1-costhg*beta)*
2257  $ (mchren+2*xph*sqrt(mpasqr))/
2258  $ mpasqr**2/(1-mchren/mpasqr)/(4-mchren/mpasqr))
2259  ENDIF
2260 c write(*,*) IDPHO(1),IDPHO(JDAPHO(1,1)),IDPHO(JDAPHO(1,1)+1)
2261 c write(*,*) emu,xph,costhg,beta,mpasqr,mchren
2262 
2263  END
2264  SUBROUTINE phodo(IP,NCHARB,NEUDAU)
2265 C.----------------------------------------------------------------------
2266 C.
2267 C. PHOTOS: PHOton radiation in decays DOing of KINematics
2268 C.
2269 C. Purpose: Starting from the charged particle energy/momentum,
2270 C. PNEUTR, photon energy fraction and photon angle with
2271 C. respect to the axis formed by charged particle energy/
2272 C. momentum vector and PNEUTR, scale the energy/momentum,
2273 C. keeping the original direction of the neutral system in
2274 C. the lab. frame untouched.
2275 C.
2276 C. Input Parameters: IP: Pointer to decaying particle in
2277 C. /PHOEVT/ and the common itself
2278 C. NCHARB: pointer to the charged radiating
2279 C. daughter in /PHOEVT/.
2280 C. NEUDAU: pointer to the first neutral daughter
2281 C. Output Parameters: Common /PHOEVT/, with photon added.
2282 C.
2283 C. Author(s): Z. Was, B. van Eijk Created at: 26/11/89
2284 C. Last Update: 27/05/93
2285 C.
2286 C.----------------------------------------------------------------------
2287  IMPLICIT NONE
2288  DOUBLE PRECISION phoan1,phoan2,angle,fi1,fi3,fi4,fi5,th1,th3,th4
2289  DOUBLE PRECISION parne,qnew,qold,data
2290  INTEGER ip,fi3dum,i,j,neudau,first,last
2291  INTEGER ncharb
2292  REAL*8 ephoto,pmavir,photri
2293  REAL*8 gneut,phoran,ccosth,ssinth,pvec(4)
2294  INTEGER nmxpho
2295  parameter(nmxpho=10000)
2296  INTEGER idpho,istpho,jdapho,jmopho,nevpho,npho
2297  REAL*8 ppho,vpho
2298  common/phoevt/nevpho,npho,istpho(nmxpho),idpho(nmxpho),
2299  &jmopho(2,nmxpho),jdapho(2,nmxpho),ppho(5,nmxpho),vpho(4,nmxpho)
2300  DOUBLE PRECISION mchsqr,mnesqr
2301  REAL*8 pneutr
2302  common/phomom/mchsqr,mnesqr,pneutr(5)
2303  DOUBLE PRECISION costhg,sinthg
2304  REAL*8 xphmax,xphoto
2305  common/phophs/xphmax,xphoto,costhg,sinthg
2306  REAL*8 pi,twopi
2307  common/phpico/pi,twopi
2308 C--
2309  ephoto=xphoto*ppho(5,ip)/2.d0
2310  pmavir=sqrt(ppho(5,ip)*(ppho(5,ip)-2.d0*ephoto))
2311 C--
2312 C-- Reconstruct kinematics of charged particle and neutral system
2313  fi1=phoan1(pneutr(1),pneutr(2))
2314 C--
2315 C-- Choose axis along z of PNEUTR, calculate angle between x and y
2316 C-- components and z and x-y plane and perform Lorentz transform...
2317  th1=phoan2(pneutr(3),sqrt(pneutr(1)**2+pneutr(2)**2))
2318  CALL phoro3(-fi1,pneutr(1))
2319  CALL phoro2(-th1,pneutr(1))
2320 C--
2321 C-- Take away photon energy from charged particle and PNEUTR ! Thus
2322 C-- the onshell charged particle decays into virtual charged particle
2323 C-- and photon. The virtual charged particle mass becomes:
2324 C-- SQRT(PPHO(5,IP)*(PPHO(5,IP)-2*EPHOTO)). Construct new PNEUTR mo-
2325 C-- mentum in the rest frame of the parent:
2326 C-- 1) Scaling parameters...
2327  qnew=photri(pmavir,pneutr(5),ppho(5,ncharb))
2328  qold=pneutr(3)
2329  gneut=(qnew**2+qold**2+mnesqr)/(qnew*qold+sqrt((qnew**2+mnesqr)*
2330  &(qold**2+mnesqr)))
2331  IF (gneut.LT.1.d0) THEN
2332  data=0.d0
2333  CALL phoerr(4,'PHOKIN',data)
2334  ENDIF
2335  parne=gneut-sqrt(max(gneut**2-1.0d0,0.d0))
2336 C--
2337 C-- 2) ...reductive boost...
2338  CALL phobo3(parne,pneutr)
2339 C--
2340 C-- ...calculate photon energy in the reduced system...
2341  npho=npho+1
2342  istpho(npho)=1
2343  idpho(npho) =22
2344 C-- Photon mother and daughter pointers !
2345  jmopho(1,npho)=ip
2346  jmopho(2,npho)=0
2347  jdapho(1,npho)=0
2348  jdapho(2,npho)=0
2349  ppho(4,npho)=ephoto*ppho(5,ip)/pmavir
2350 C--
2351 C-- ...and photon momenta
2352  ccosth=-costhg
2353  ssinth=sinthg
2354  th3=phoan2(ccosth,ssinth)
2355  fi3=twopi*phoran(fi3dum)
2356  ppho(1,npho)=ppho(4,npho)*sinthg*cos(fi3)
2357  ppho(2,npho)=ppho(4,npho)*sinthg*sin(fi3)
2358 C--
2359 C-- Minus sign because axis opposite direction of charged particle !
2360  ppho(3,npho)=-ppho(4,npho)*costhg
2361  ppho(5,npho)=0.d0
2362 C--
2363 C-- Rotate in order to get photon along z-axis
2364  CALL phoro3(-fi3,pneutr(1))
2365  CALL phoro3(-fi3,ppho(1,npho))
2366  CALL phoro2(-th3,pneutr(1))
2367  CALL phoro2(-th3,ppho(1,npho))
2368  angle=ephoto/ppho(4,npho)
2369 C--
2370 C-- Boost to the rest frame of decaying particle
2371  CALL phobo3(angle,pneutr(1))
2372  CALL phobo3(angle,ppho(1,npho))
2373 C--
2374 C-- Back in the parent rest frame but PNEUTR not yet oriented !
2375  fi4=phoan1(pneutr(1),pneutr(2))
2376  th4=phoan2(pneutr(3),sqrt(pneutr(1)**2+pneutr(2)**2))
2377  CALL phoro3(fi4,pneutr(1))
2378  CALL phoro3(fi4,ppho(1,npho))
2379 C--
2380  DO 60 i=2,4
2381  60 pvec(i)=0.d0
2382  pvec(1)=1.d0
2383  CALL phoro3(-fi3,pvec)
2384  CALL phoro2(-th3,pvec)
2385  CALL phobo3(angle,pvec)
2386  CALL phoro3(fi4,pvec)
2387  CALL phoro2(-th4,pneutr)
2388  CALL phoro2(-th4,ppho(1,npho))
2389  CALL phoro2(-th4,pvec)
2390  fi5=phoan1(pvec(1),pvec(2))
2391 C--
2392 C-- Charged particle restores original direction
2393  CALL phoro3(-fi5,pneutr)
2394  CALL phoro3(-fi5,ppho(1,npho))
2395  CALL phoro2(th1,pneutr(1))
2396  CALL phoro2(th1,ppho(1,npho))
2397  CALL phoro3(fi1,pneutr)
2398  CALL phoro3(fi1,ppho(1,npho))
2399 C-- See whether neutral system has multiplicity larger than 1...
2400  IF ((jdapho(2,ip)-jdapho(1,ip)).GT.1) THEN
2401 C-- Find pointers to components of 'neutral' system
2402 C--
2403  first=neudau
2404  last=jdapho(2,ip)
2405  DO 70 i=first,last
2406  IF (i.NE.ncharb.AND.(jmopho(1,i).EQ.ip)) THEN
2407 C--
2408 C-- Reconstruct kinematics...
2409  CALL phoro3(-fi1,ppho(1,i))
2410  CALL phoro2(-th1,ppho(1,i))
2411 C--
2412 C-- ...reductive boost
2413  CALL phobo3(parne,ppho(1,i))
2414 C--
2415 C-- Rotate in order to get photon along z-axis
2416  CALL phoro3(-fi3,ppho(1,i))
2417  CALL phoro2(-th3,ppho(1,i))
2418 C--
2419 C-- Boost to the rest frame of decaying particle
2420  CALL phobo3(angle,ppho(1,i))
2421 C--
2422 C-- Back in the parent rest-frame but PNEUTR not yet oriented.
2423  CALL phoro3(fi4,ppho(1,i))
2424  CALL phoro2(-th4,ppho(1,i))
2425 C--
2426 C-- Charged particle restores original direction
2427  CALL phoro3(-fi5,ppho(1,i))
2428  CALL phoro2(th1,ppho(1,i))
2429  CALL phoro3(fi1,ppho(1,i))
2430  ENDIF
2431  70 CONTINUE
2432  ELSE
2433 C--
2434 C-- ...only one 'neutral' particle in addition to photon!
2435  DO 80 j=1,4
2436  80 ppho(j,neudau)=pneutr(j)
2437  ENDIF
2438 C--
2439 C-- All 'neutrals' treated, fill /PHOEVT/ for charged particle...
2440  DO 90 j=1,3
2441  90 ppho(j,ncharb)=-(ppho(j,npho)+pneutr(j))
2442  ppho(4,ncharb)=ppho(5,ip)-(ppho(4,npho)+pneutr(4))
2443 C--
2444  END
2445  FUNCTION photri(A,B,C)
2446 C.----------------------------------------------------------------------
2447 C.
2448 C. PHOTOS: PHOton radiation in decays calculation of TRIangle fie
2449 C.
2450 C. Purpose: Calculation of triangle function for phase space.
2451 C.
2452 C. Input Parameters: A, B, C (Virtual) particle masses.
2453 C.
2454 C. Output Parameter: Function value =
2455 C. SQRT(LAMBDA(A**2,B**2,C**2))/(2*A)
2456 C.
2457 C. Author(s): B. van Eijk Created at: 15/11/89
2458 C. Last Update: 02/01/90
2459 C.
2460 C.----------------------------------------------------------------------
2461  IMPLICIT NONE
2462  DOUBLE PRECISION da,db,dc,dapb,damb,dtrian
2463  REAL*8 a,b,c,photri
2464  da=a
2465  db=b
2466  dc=c
2467  dapb=da+db
2468  damb=da-db
2469  dtrian=sqrt((damb-dc)*(dapb+dc)*(damb+dc)*(dapb-dc))
2470  photri=dtrian/(da+da)
2471  RETURN
2472  END
2473  FUNCTION phoan1(X,Y)
2474 C.----------------------------------------------------------------------
2475 C.
2476 C. PHOTOS: PHOton radiation in decays calculation of ANgle '1'
2477 C.
2478 C. Purpose: Calculate angle from X and Y
2479 C.
2480 C. Input Parameters: X, Y
2481 C.
2482 C. Output Parameter: Function value
2483 C.
2484 C. Author(s): S. Jadach Created at: 01/01/89
2485 C. B. van Eijk Last Update: 02/01/90
2486 C.
2487 C.----------------------------------------------------------------------
2488  IMPLICIT NONE
2489  DOUBLE PRECISION phoan1
2490  REAL*8 x,y
2491  REAL*8 pi,twopi
2492  common/phpico/pi,twopi
2493  IF (abs(y).LT.abs(x)) THEN
2494  phoan1=atan(abs(y/x))
2495  IF (x.LE.0.d0) phoan1=pi-phoan1
2496  ELSE
2497  phoan1=acos(x/sqrt(x**2+y**2))
2498  ENDIF
2499  IF (y.LT.0.d0) phoan1=twopi-phoan1
2500  RETURN
2501  END
2502  FUNCTION phoan2(X,Y)
2503 C.----------------------------------------------------------------------
2504 C.
2505 C. PHOTOS: PHOton radiation in decays calculation of ANgle '2'
2506 C.
2507 C. Purpose: Calculate angle from X and Y
2508 C.
2509 C. Input Parameters: X, Y
2510 C.
2511 C. Output Parameter: Function value
2512 C.
2513 C. Author(s): S. Jadach Created at: 01/01/89
2514 C. B. van Eijk Last Update: 02/01/90
2515 C.
2516 C.----------------------------------------------------------------------
2517  IMPLICIT NONE
2518  DOUBLE PRECISION phoan2
2519  REAL*8 x,y
2520  REAL*8 pi,twopi
2521  common/phpico/pi,twopi
2522  IF (abs(y).LT.abs(x)) THEN
2523  phoan2=atan(abs(y/x))
2524  IF (x.LE.0.d0) phoan2=pi-phoan2
2525  ELSE
2526  phoan2=acos(x/sqrt(x**2+y**2))
2527  ENDIF
2528  RETURN
2529  END
2530  SUBROUTINE phobo3(ANGLE,PVEC)
2531 C.----------------------------------------------------------------------
2532 C.
2533 C. PHOTOS: PHOton radiation in decays BOost routine '3'
2534 C.
2535 C. Purpose: Boost vector PVEC along z-axis where ANGLE = EXP(ETA),
2536 C. ETA is the hyperbolic velocity.
2537 C.
2538 C. Input Parameters: ANGLE, PVEC
2539 C.
2540 C. Output Parameter: PVEC
2541 C.
2542 C. Author(s): S. Jadach Created at: 01/01/89
2543 C. B. van Eijk Last Update: 02/01/90
2544 C.
2545 C.----------------------------------------------------------------------
2546  IMPLICIT NONE
2547  DOUBLE PRECISION qpl,qmi,angle
2548  REAL*8 pvec(4)
2549  qpl=(pvec(4)+pvec(3))*angle
2550  qmi=(pvec(4)-pvec(3))/angle
2551  pvec(3)=(qpl-qmi)/2.d0
2552  pvec(4)=(qpl+qmi)/2.d0
2553  RETURN
2554  END
2555  SUBROUTINE phoro2(ANGLE,PVEC)
2556 C.----------------------------------------------------------------------
2557 C.
2558 C. PHOTOS: PHOton radiation in decays ROtation routine '2'
2559 C.
2560 C. Purpose: Rotate x and z components of vector PVEC around angle
2561 C. 'ANGLE'.
2562 C.
2563 C. Input Parameters: ANGLE, PVEC
2564 C.
2565 C. Output Parameter: PVEC
2566 C.
2567 C. Author(s): S. Jadach Created at: 01/01/89
2568 C. B. van Eijk Last Update: 02/01/90
2569 C.
2570 C.----------------------------------------------------------------------
2571  IMPLICIT NONE
2572  DOUBLE PRECISION cs,sn,angle
2573  REAL*8 pvec(4)
2574  cs=cos(angle)*pvec(1)+sin(angle)*pvec(3)
2575  sn=-sin(angle)*pvec(1)+cos(angle)*pvec(3)
2576  pvec(1)=cs
2577  pvec(3)=sn
2578  RETURN
2579  END
2580  SUBROUTINE phoro3(ANGLE,PVEC)
2581 C.----------------------------------------------------------------------
2582 C.
2583 C. PHOTOS: PHOton radiation in decays ROtation routine '3'
2584 C.
2585 C. Purpose: Rotate x and y components of vector PVEC around angle
2586 C. 'ANGLE'.
2587 C.
2588 C. Input Parameters: ANGLE, PVEC
2589 C.
2590 C. Output Parameter: PVEC
2591 C.
2592 C. Author(s): S. Jadach Created at: 01/01/89
2593 C. B. van Eijk Last Update: 02/01/90
2594 C.
2595 C.----------------------------------------------------------------------
2596  IMPLICIT NONE
2597  DOUBLE PRECISION cs,sn,angle
2598  REAL*8 pvec(4)
2599  cs=cos(angle)*pvec(1)-sin(angle)*pvec(2)
2600  sn=sin(angle)*pvec(1)+cos(angle)*pvec(2)
2601  pvec(1)=cs
2602  pvec(2)=sn
2603  RETURN
2604  END
2605 #include "../randg/photos-random.h"
2606  FUNCTION phocha(IDHEP)
2607 C.----------------------------------------------------------------------
2608 C.
2609 C. PHOTOS: PHOton radiation in decays CHArge determination
2610 C.
2611 C. Purpose: Calculate the charge of particle with code IDHEP. The
2612 C. code of the particle is defined by the Particle Data
2613 C. Group in Phys. Lett. B204 (1988) 1.
2614 C.
2615 C. Input Parameter: IDHEP
2616 C.
2617 C. Output Parameter: Funtion value = charge of particle with code
2618 C. IDHEP
2619 C.
2620 C. Author(s): E. Barberio and B. van Eijk Created at: 29/11/89
2621 C. Last update: 02/01/90
2622 C.
2623 C.----------------------------------------------------------------------
2624  IMPLICIT NONE
2625  REAL*8 phocha
2626  INTEGER idhep,idabs,q1,q2,q3
2627 C--
2628 C-- Array 'CHARGE' contains the charge of the first 101 particles ac-
2629 C-- cording to the PDG particle code... (0 is added for convenience)
2630  REAL*8 charge(0:100)
2631  DATA charge/ 0.d0,
2632  &-0.3333333333d0, 0.6666666667d0, -0.3333333333d0, 0.6666666667d0,
2633  &-0.3333333333d0, 0.6666666667d0, -0.3333333333d0, 0.6666666667d0,
2634  & 2*0.d0, -1.d0, 0.d0, -1.d0, 0.d0, -1.d0, 0.d0, -1.d0, 6*0.d0,
2635  & 1.d0, 12*0.d0, 1.d0, 63*0.d0/
2636  idabs=abs(idhep)
2637  IF (idabs.LE.100) THEN
2638 C--
2639 C-- Charge of quark, lepton, boson etc....
2640  phocha = charge(idabs)
2641  ELSE
2642 C--
2643 C-- Check on particle build out of quarks, unpack its code...
2644  q3=mod(idabs/1000,10)
2645  q2=mod(idabs/100,10)
2646  q1=mod(idabs/10,10)
2647  IF (q3.EQ.0) THEN
2648 C--
2649 C-- ...meson...
2650  IF(mod(q2,2).EQ.0) THEN
2651  phocha=charge(q2)-charge(q1)
2652  ELSE
2653  phocha=charge(q1)-charge(q2)
2654  ENDIF
2655  ELSE
2656 C--
2657 C-- ...diquarks or baryon.
2658  phocha=charge(q1)+charge(q2)+charge(q3)
2659  ENDIF
2660  ENDIF
2661 C--
2662 C-- Find the sign of the charge...
2663  IF (idhep.LT.0.d0) phocha=-phocha
2664  IF (phocha**2.lt.1d-6) phocha=0.d0
2665  RETURN
2666  END
2667  FUNCTION phospi(IDHEP)
2668 C.----------------------------------------------------------------------
2669 C.
2670 C. PHOTOS: PHOton radiation in decays function for SPIn determina-
2671 C. tion
2672 C.
2673 C. Purpose: Calculate the spin of particle with code IDHEP. The
2674 C. code of the particle is defined by the Particle Data
2675 C. Group in Phys. Lett. B204 (1988) 1.
2676 C.
2677 C. Input Parameter: IDHEP
2678 C.
2679 C. Output Parameter: Funtion value = spin of particle with code
2680 C. IDHEP
2681 C.
2682 C. Author(s): E. Barberio and B. van Eijk Created at: 29/11/89
2683 C. Last update: 02/01/90
2684 C.
2685 C.----------------------------------------------------------------------
2686  IMPLICIT NONE
2687  REAL*8 phospi
2688  INTEGER idhep,idabs
2689 C--
2690 C-- Array 'SPIN' contains the spin of the first 100 particles accor-
2691 C-- ding to the PDG particle code...
2692  REAL*8 spin(100)
2693  DATA spin/ 8*.5d0, 1.d0, 0.d0, 8*.5d0, 2*0.d0, 4*1.d0, 76*0.d0/
2694  idabs=abs(idhep)
2695 C--
2696 C-- Spin of quark, lepton, boson etc....
2697  IF (idabs.LE.100) THEN
2698  phospi=spin(idabs)
2699  ELSE
2700 C--
2701 C-- ...other particles, however...
2702  phospi=(mod(idabs,10)-1.d0)/2.d0
2703 C--
2704 C-- ...K_short and K_long are special !!
2705  phospi=max(phospi,0.d0)
2706  ENDIF
2707  RETURN
2708  END
2709  SUBROUTINE phoerr(IMES,TEXT,DATA)
2710 C.----------------------------------------------------------------------
2711 C.
2712 C. PHOTOS: PHOton radiation in decays ERRror handling
2713 C.
2714 C. Purpose: Inform user about (fatal) errors and warnings generated
2715 C. by either the user or the program.
2716 C.
2717 C. Input Parameters: IMES, TEXT, DATA
2718 C.
2719 C. Output Parameters: None
2720 C.
2721 C. Author(s): B. van Eijk Created at: 29/11/89
2722 C. Last Update: 10/01/92
2723 C.
2724 C.----------------------------------------------------------------------
2725  IMPLICIT NONE
2726  DOUBLE PRECISION data
2727  INTEGER imes,ierror
2728  REAL*8 sdata
2729  INTEGER phlun
2730  common/pholun/phlun
2731  INTEGER phomes
2732  parameter(phomes=10)
2733  INTEGER status
2734  common/phosta/status(phomes)
2735  CHARACTER text*(*)
2736  SAVE ierror
2737 C-- security STOP switch
2738  LOGICAL isec
2739  SAVE isec
2740  DATA isec /.true./
2741  DATA ierror/ 0/
2742  IF (imes.LE.phomes) status(imes)=status(imes)+1
2743 C--
2744 C-- Count number of non-fatal errors...
2745  IF ((imes.EQ. 6).AND.(status(imes).GE.2)) RETURN
2746  IF ((imes.EQ.10).AND.(status(imes).GE.2)) RETURN
2747  sdata=DATA
2748  WRITE(phlun,9000)
2749  WRITE(phlun,9120)
2750  goto(10,20,30,40,50,60,70,80,90,100),imes
2751  WRITE(phlun,9130) imes
2752  goto 120
2753  10 WRITE(phlun,9010) text,int(sdata)
2754  goto 110
2755  20 WRITE(phlun,9020) text,sdata
2756  goto 110
2757  30 WRITE(phlun,9030) text,sdata
2758  goto 110
2759  40 WRITE(phlun,9040) text
2760  goto 110
2761  50 WRITE(phlun,9050) text,int(sdata)
2762  goto 110
2763  60 WRITE(phlun,9060) text,sdata
2764  goto 130
2765  70 WRITE(phlun,9070) text,int(sdata)
2766  goto 110
2767  80 WRITE(phlun,9080) text,int(sdata)
2768  goto 110
2769  90 WRITE(phlun,9090) text,int(sdata)
2770  goto 110
2771  100 WRITE(phlun,9100) text,sdata
2772  goto 130
2773  110 CONTINUE
2774  WRITE(phlun,9140)
2775  WRITE(phlun,9120)
2776  WRITE(phlun,9000)
2777  IF (isec) THEN
2778  stop
2779  ELSE
2780  goto 130
2781  ENDIF
2782  120 ierror=ierror+1
2783  IF (ierror.GE.10) THEN
2784  WRITE(phlun,9150)
2785  WRITE(phlun,9120)
2786  WRITE(phlun,9000)
2787  IF (isec) THEN
2788  stop
2789  ELSE
2790  goto 130
2791  ENDIF
2792  ENDIF
2793  130 WRITE(phlun,9120)
2794  WRITE(phlun,9000)
2795  RETURN
2796  9000 FORMAT(1h ,80('*'))
2797  9010 FORMAT(1h ,'* ',a,': Too many charged Particles, NCHARG =',i6,t81,
2798  &'*')
2799  9020 FORMAT(1h ,'* ',a,': Too much Bremsstrahlung required, PRSOFT = ',
2800  &f15.6,t81,'*')
2801  9030 FORMAT(1h ,'* ',a,': Combined Weight is exceeding 1., Weight = ',
2802  &f15.6,t81,'*')
2803  9040 FORMAT(1h ,'* ',a,
2804  &': Error in Rescaling charged and neutral Vectors',t81,'*')
2805  9050 FORMAT(1h ,'* ',a,
2806  &': Non matching charged Particle Pointer, NCHARG = ',i5,t81,'*')
2807  9060 FORMAT(1h ,'* ',a,
2808  &': Do you really work with a Particle of Spin: ',f4.1,' ?',t81,
2809  &'*')
2810  9070 FORMAT(1h ,'* ',a, ': Stack Length exceeded, NSTACK = ',i5 ,t81,
2811  &'*')
2812  9080 FORMAT(1h ,'* ',a,
2813  &': Random Number Generator Seed(1) out of Range: ',i8,t81,'*')
2814  9090 FORMAT(1h ,'* ',a,
2815  &': Random Number Generator Seed(2) out of Range: ',i8,t81,'*')
2816  9100 FORMAT(1h ,'* ',a,
2817  &': Available Phase Space below Cut-off: ',f15.6,' GeV/c^2',t81,
2818  &'*')
2819  9120 FORMAT(1h ,'*',t81,'*')
2820  9130 FORMAT(1h ,'* Funny Error Message: ',i4,' ! What to do ?',t81,'*')
2821  9140 FORMAT(1h ,'* Fatal Error Message, I stop this Run !',t81,'*')
2822  9150 FORMAT(1h ,'* 10 Error Messages generated, I stop this Run !',t81,
2823  &'*')
2824  END
2825  SUBROUTINE phorep
2826 C.----------------------------------------------------------------------
2827 C.
2828 C. PHOTOS: PHOton radiation in decays run summary REPort
2829 C.
2830 C. Purpose: Inform user about success and/or restrictions of PHOTOS
2831 C. encountered during execution.
2832 C.
2833 C. Input Parameters: Common /PHOSTA/
2834 C.
2835 C. Output Parameters: None
2836 C.
2837 C. Author(s): B. van Eijk Created at: 10/01/92
2838 C. Last Update: 10/01/92
2839 C.
2840 C.----------------------------------------------------------------------
2841  IMPLICIT NONE
2842  INTEGER phlun
2843  common/pholun/phlun
2844  INTEGER phomes
2845  parameter(phomes=10)
2846  INTEGER status
2847  common/phosta/status(phomes)
2848  INTEGER i
2849  LOGICAL error
2850  error=.false.
2851  WRITE(phlun,9000)
2852  WRITE(phlun,9010)
2853  WRITE(phlun,9020)
2854  WRITE(phlun,9030)
2855  WRITE(phlun,9040)
2856  WRITE(phlun,9030)
2857  WRITE(phlun,9020)
2858  DO 10 i=1,phomes
2859  IF (status(i).EQ.0) goto 10
2860  IF ((i.EQ.6).OR.(i.EQ.10)) THEN
2861  WRITE(phlun,9050) i,status(i)
2862  ELSE
2863  error=.true.
2864  WRITE(phlun,9060) i,status(i)
2865  ENDIF
2866  10 CONTINUE
2867  IF (.NOT.error) WRITE(phlun,9070)
2868  WRITE(phlun,9020)
2869  WRITE(phlun,9010)
2870  RETURN
2871  9000 FORMAT(1h1)
2872  9010 FORMAT(1h ,80('*'))
2873  9020 FORMAT(1h ,'*',t81,'*')
2874  9030 FORMAT(1h ,'*',26x,25('='),t81,'*')
2875  9040 FORMAT(1h ,'*',30x,'PHOTOS Run Summary',t81,'*')
2876  9050 FORMAT(1h ,'*',22x,'Warning #',i2,' occured',i6,' times',t81,'*')
2877  9060 FORMAT(1h ,'*',23x,'Error #',i2,' occured',i6,' times',t81,'*')
2878  9070 FORMAT(1h ,'*',16x,'PHOTOS Execution has successfully terminated',
2879  &t81,'*')
2880  END
2881  SUBROUTINE phlupa(IPOINT)
2882  IMPLICIT NONE
2883 C.----------------------------------------------------------------------
2884 C.
2885 C. PHLUPA: debugging tool
2886 C.
2887 C. Purpose: NONE, eventually may printout content of the
2888 C. /PHOEVT/ common
2889 C.
2890 C. Input Parameters: Common /PHOEVT/ and /PHNUM/
2891 C. latter may have number of the event.
2892 C.
2893 C. Output Parameters: None
2894 C.
2895 C. Author(s): Z. Was Created at: 30/05/93
2896 C. Last Update: 09/10/05
2897 C.
2898 C.----------------------------------------------------------------------
2899  INTEGER nmxpho
2900  parameter(nmxpho=10000)
2901  INTEGER idpho,istpho,jdapho,jmopho,nevpho,npho,i,j,ipoint
2902  INTEGER ipoin,ipoin0,ipoinm,iev
2903  INTEGER iout
2904  REAL*8 ppho,vpho,sum
2905  common/phoevt/nevpho,npho,istpho(nmxpho),idpho(nmxpho),
2906  &jmopho(2,nmxpho),jdapho(2,nmxpho),ppho(5,nmxpho),vpho(4,nmxpho)
2907  COMMON /phnum/ iev
2908  INTEGER phlun
2909  common/pholun/phlun
2910  dimension sum(5)
2911  DATA ipoin0/ -5/
2912  COMMON /phlupy/ ipoin,ipoinm
2913  SAVE ipoin0
2914  IF (ipoin0.LT.0) THEN
2915  ipoin0=300 000 ! maximal no-print point
2916  ipoin =ipoin0
2917  ipoinm=300 001 ! minimal no-print point
2918  ENDIF
2919  IF (ipoint.LE.ipoinm.OR.ipoint.GE.ipoin ) RETURN
2920  iout=56
2921  IF (iev.LT.1000) THEN
2922  DO i=1,5
2923  sum(i)=0.0d0
2924  ENDDO
2925  WRITE(phlun,*) 'EVENT NR=',iev,
2926  $ 'WE ARE TESTING /PHOEVT/ at IPOINT=',ipoint
2927  WRITE(phlun,10)
2928  i=1
2929  WRITE(phlun,20) idpho(i),ppho(1,i),ppho(2,i),ppho(3,i),
2930  $ ppho(4,i),ppho(5,i),jdapho(1,i),jdapho(2,i)
2931  i=2
2932  WRITE(phlun,20) idpho(i),ppho(1,i),ppho(2,i),ppho(3,i),
2933  $ ppho(4,i),ppho(5,i),jdapho(1,i),jdapho(2,i)
2934  WRITE(phlun,*) ' '
2935  DO i=3,npho
2936  WRITE(phlun,20) idpho(i),ppho(1,i),ppho(2,i),ppho(3,i),
2937  $ ppho(4,i),ppho(5,i),jmopho(1,i),jmopho(2,i)
2938  DO j=1,4
2939  sum(j)=sum(j)+ppho(j,i)
2940  ENDDO
2941  ENDDO
2942  sum(5)=sqrt(abs(sum(4)**2-sum(1)**2-sum(2)**2-sum(3)**2))
2943  WRITE(phlun,30) sum
2944  10 FORMAT(1x,' ID ','p_x ','p_y ','p_z ',
2945  $ 'E ','m ',
2946  $ 'ID-MO_DA1','ID-MO DA2' )
2947  20 FORMAT(1x,i4,5(f14.9),2i9)
2948  30 FORMAT(1x,' SUM',5(f14.9))
2949  ENDIF
2950  END
2951 
2952 
2953 
2954  FUNCTION iphqrk(MODCOR)
2955  implicit none
2956 C.----------------------------------------------------------------------
2957 C.
2958 C. IPHQRK: enables blocks emision from quarks
2959 C.
2960 C
2961 C. Input Parameters: MODCOR
2962 C. MODCOR >0 type of action
2963 C. =1 blocks
2964 C. =2 enables
2965 C. =0 execution mode (retrns stored value)
2966 C.
2967 C.
2968 C. Author(s): Z. Was Created at: 11/12/00
2969 C. Modified :
2970 C.----------------------------------------------------------------------
2971  INTEGER iphqrk,modcor,modop
2972  INTEGER phlun
2973  common/pholun/phlun
2974 
2975  SAVE modop
2976  DATA modop /0/
2977  IF (modcor.NE.0) THEN
2978 C INITIALIZATION
2979  modop=modcor
2980 
2981  WRITE(phlun,*)
2982  $ 'Message from PHOTOS: IPHQRK(MODCOR):: (re)initialization'
2983  IF (modop.EQ.1) THEN
2984  WRITE(phlun,*)
2985  $ 'MODOP=1 -- blocks emission from light quarks: DEFAULT'
2986  ELSEIF (modop.EQ.2) THEN
2987  WRITE(phlun,*)
2988  $ 'MODOP=2 -- enables emission from light quarks: TEST '
2989  ELSE
2990  WRITE(phlun,*) 'IPHQRK wrong MODCOR=',modcor
2991  stop
2992  ENDIF
2993  RETURN
2994  ENDIF
2995 
2996  IF (modop.EQ.0.AND.modcor.EQ.0) THEN
2997  WRITE(phlun,*) 'IPHQRK lack of initialization'
2998  stop
2999  ENDIF
3000  iphqrk=modop
3001  END
3002 
3003 
3004  FUNCTION iphekl(MODCOR)
3005  implicit none
3006 C.----------------------------------------------------------------------
3007 C.
3008 C. IPHEKL: enables/blocks emision in: pi0 to gamma e+ e-
3009 C.
3010 C
3011 C. Input Parameters: MODCOR
3012 C. MODCOR >0 type of action
3013 C. =1 blocks
3014 C. =2 enables
3015 C. =0 execution mode (retrns stored value)
3016 C.
3017 C.
3018 C. Author(s): Z. Was Created at: 11/12/00
3019 C. Modified :
3020 C.----------------------------------------------------------------------
3021  INTEGER iphekl,modcor,modop
3022  INTEGER phlun
3023  common/pholun/phlun
3024 
3025  SAVE modop
3026  DATA modop /0/
3027 
3028  IF (modcor.NE.0) THEN
3029 C INITIALIZATION
3030  modop=modcor
3031 
3032  WRITE(phlun,*)
3033  $ 'Message from PHOTOS: IPHEKL(MODCOR):: (re)initialization'
3034  IF (modop.EQ.2) THEN
3035  WRITE(phlun,*)
3036  $ 'MODOP=2 -- blocks emission in pi0 to gamma e+e-: DEFAULT'
3037  WRITE(phlun,*)
3038  $ 'MODOP=2 -- blocks emission in Kl to gamma e+e-: DEFAULT'
3039  ELSEIF (modop.EQ.1) THEN
3040  WRITE(phlun,*)
3041  $ 'MODOP=1 -- enables emission in pi0 to gamma e+e- : TEST '
3042  WRITE(phlun,*)
3043  $ 'MODOP=1 -- enables emission in Kl to gamma e+e- : TEST '
3044  ELSE
3045  WRITE(phlun,*) 'IPHEKL wrong MODCOR=',modcor
3046  stop
3047  ENDIF
3048  RETURN
3049  ENDIF
3050 
3051  IF (modop.EQ.0.AND.modcor.EQ.0) THEN
3052  WRITE(phlun,*) 'IPHELK lack of initialization'
3053  stop
3054  ENDIF
3055  iphekl=modop
3056  END
3057 
3058  SUBROUTINE phcork(MODCOR)
3059  implicit none
3060 C.----------------------------------------------------------------------
3061 C.
3062 C. PHCORK: corrects kinmatics of subbranch needed if host program
3063 C. produces events with the shaky momentum conservation
3064 C
3065 C. Input Parameters: Common /PHOEVT/, MODCOR
3066 C. MODCOR >0 type of action
3067 C. =1 no action
3068 C. =2 corrects energy from mass
3069 C. =3 corrects mass from energy
3070 C. =4 corrects energy from mass for
3071 C. particles up to .4 GeV mass,
3072 C. for heavier ones corrects mass,
3073 C. =5 most complete correct also of mother
3074 C. often necessary for exponentiation.
3075 C. =0 execution mode
3076 C.
3077 C. Output Parameters: corrected /PHOEVT/
3078 C.
3079 C. Author(s): P.Golonka, Z. Was Created at: 01/02/99
3080 C. Modified : 08/02/99
3081 C.----------------------------------------------------------------------
3082  INTEGER nmxpho
3083  parameter(nmxpho=10000)
3084 
3085  REAL*8 m,p2,px,py,pz,e,en,mcut,xms
3086  INTEGER modcor,modop,i,iev,iprint,k
3087  INTEGER idpho,istpho,jdapho,jmopho,nevpho,npho
3088  REAL*8 ppho,vpho
3089  common/phoevt/nevpho,npho,istpho(nmxpho),idpho(nmxpho),
3090  &jmopho(2,nmxpho),jdapho(2,nmxpho),ppho(5,nmxpho),vpho(4,nmxpho)
3091 
3092  INTEGER phlun
3093  common/pholun/phlun
3094 
3095  COMMON /phnum/ iev
3096  SAVE modop
3097  DATA modop /0/
3098  SAVE iprint
3099  DATA iprint /0/
3100  SAVE mcut
3101  IF (modcor.NE.0) THEN
3102 C INITIALIZATION
3103  modop=modcor
3104 
3105  WRITE(phlun,*) 'Message from PHCORK(MODCOR):: initialization'
3106  IF (modop.EQ.1) THEN
3107  WRITE(phlun,*) 'MODOP=1 -- no corrections on event: DEFAULT'
3108  ELSEIF (modop.EQ.2) THEN
3109  WRITE(phlun,*) 'MODOP=2 -- corrects Energy from mass'
3110  ELSEIF (modop.EQ.3) THEN
3111  WRITE(phlun,*) 'MODOP=3 -- corrects mass from Energy'
3112  ELSEIF (modop.EQ.4) THEN
3113  WRITE(phlun,*) 'MODOP=4 -- corrects Energy from mass to Mcut'
3114  WRITE(phlun,*) ' and mass from energy above Mcut '
3115  mcut=0.4
3116  WRITE(phlun,*) 'Mcut=',mcut,'GeV'
3117  ELSEIF (modop.EQ.5) THEN
3118  WRITE(phlun,*) 'MODOP=5 -- corrects Energy from mass+flow'
3119 
3120  ELSE
3121  WRITE(phlun,*) 'PHCORK wrong MODCOR=',modcor
3122  stop
3123  ENDIF
3124  RETURN
3125  ENDIF
3126 
3127  IF (modop.EQ.0.AND.modcor.EQ.0) THEN
3128  WRITE(phlun,*) 'PHCORK lack of initialization'
3129  stop
3130  ENDIF
3131 
3132 C execution mode
3133 C ==============
3134 C ==============
3135 
3136 
3137  px=0
3138  py=0
3139  pz=0
3140  e =0
3141 
3142  IF (modop.EQ.1) THEN
3143 C -----------------------
3144 C In this case we do nothing
3145  RETURN
3146  ELSEIF(modop.EQ.2) THEN
3147 C -----------------------
3148 CC lets loop thru all daughters and correct their energies
3149 CC according to E^2=p^2+m^2
3150 
3151  DO i=3,npho
3152 
3153  px=px+ppho(1,i)
3154  py=py+ppho(2,i)
3155  pz=pz+ppho(3,i)
3156 
3157  p2=ppho(1,i)**2+ppho(2,i)**2+ppho(3,i)**2
3158 
3159  en=sqrt( ppho(5,i)**2 + p2)
3160 
3161  IF (iprint.EQ.1)
3162  & WRITE(phlun,*) "CORRECTING ENERGY OF ",i,":",
3163  & ppho(4,i),"=>",en
3164 
3165  ppho(4,i)=en
3166  e = e+ppho(4,i)
3167 
3168  ENDDO
3169 
3170  ELSEIF(modop.EQ.5) THEN
3171 C -----------------------
3172 CC lets loop thru all daughters and correct their energies
3173 CC according to E^2=p^2+m^2
3174 
3175  DO i=3,npho
3176 
3177  px=px+ppho(1,i)
3178  py=py+ppho(2,i)
3179  pz=pz+ppho(3,i)
3180 
3181  p2=ppho(1,i)**2+ppho(2,i)**2+ppho(3,i)**2
3182 
3183  en=sqrt( ppho(5,i)**2 + p2)
3184 
3185  IF (iprint.EQ.1)
3186  & WRITE(phlun,*) "CORRECTING ENERGY OF ",i,":",
3187  & ppho(4,i),"=>",en
3188 
3189  ppho(4,i)=en
3190  e = e+ppho(4,i)
3191 
3192  ENDDO
3193  DO k=1,4
3194  ppho(k,1)=0d0
3195  DO i=3,npho
3196  ppho(k,1)=ppho(k,1)+ppho(k,i)
3197  ENDDO
3198  ENDDO
3199  xms=sqrt(ppho(4,1)**2-ppho(3,1)**2-ppho(2,1)**2-ppho(1,1)**2)
3200  ppho(5,1)=xms
3201  ELSEIF(modop.EQ.3) THEN
3202 C -----------------------
3203 
3204 CC lets loop thru all daughters and correct their masses
3205 CC according to E^2=p^2+m^2
3206 
3207  DO i=3,npho
3208 
3209  px=px+ppho(1,i)
3210  py=py+ppho(2,i)
3211  pz=pz+ppho(3,i)
3212  e = e+ppho(4,i)
3213 
3214  p2=ppho(1,i)**2+ppho(2,i)**2+ppho(3,i)**2
3215 
3216  m=sqrt(abs( ppho(4,i)**2 - p2))
3217 
3218  IF (iprint.EQ.1)
3219  & WRITE(phlun,*) "CORRECTING MASS OF ",i,":",
3220  & ppho(5,i),"=>",m
3221 
3222  ppho(5,i)=m
3223 
3224  ENDDO
3225 
3226 
3227  ELSEIF(modop.EQ.4) THEN
3228 C -----------------------
3229 
3230 CC lets loop thru all daughters and correct their masses
3231 CC or energies according to E^2=p^2+m^2
3232 
3233  DO i=3,npho
3234 
3235  px=px+ppho(1,i)
3236  py=py+ppho(2,i)
3237  pz=pz+ppho(3,i)
3238 
3239  p2=ppho(1,i)**2+ppho(2,i)**2+ppho(3,i)**2
3240 
3241  m=sqrt(abs( ppho(4,i)**2 - p2))
3242 
3243  IF (m.GT.mcut) THEN
3244  IF (iprint.EQ.1)
3245  & WRITE(phlun,*) "CORRECTING MASS OF ",i,":",
3246  & ppho(5,i),"=>",m
3247  ppho(5,i)=m
3248  e = e+ppho(4,i)
3249  ELSE
3250 
3251  en=sqrt( ppho(5,i)**2 + p2)
3252 
3253  IF (iprint.EQ.1)
3254  & WRITE(phlun,*) "CORRECTING ENERGY OF ",i,":",
3255  & ppho(4,i),"=>",en
3256 
3257  ppho(4,i)=en
3258  e = e+ppho(4,i)
3259  ENDIF
3260 
3261  ENDDO
3262  ENDIF
3263 C -----
3264 
3265  IF (iprint.EQ.1) THEN
3266  WRITE(phlun,*) "CORRECTING MOTHER"
3267  WRITE(phlun,*) "PX:",ppho(1,1),"=>",px-ppho(1,2)
3268  WRITE(phlun,*) "PY:",ppho(2,1),"=>",py-ppho(2,2)
3269  WRITE(phlun,*) "PZ:",ppho(3,1),"=>",pz-ppho(3,2)
3270  WRITE(phlun,*) " E:",ppho(4,1),"=>",e-ppho(4,2)
3271  ENDIF
3272 
3273  ppho(1,1)=px-ppho(1,2)
3274  ppho(2,1)=py-ppho(2,2)
3275  ppho(3,1)=pz-ppho(3,2)
3276  ppho(4,1)=e -ppho(4,2)
3277 
3278  p2=ppho(1,1)**2+ppho(2,1)**2+ppho(3,1)**2
3279 
3280  IF (ppho(4,1)**2.GT.p2) THEN
3281  m=sqrt( ppho(4,1)**2 - p2 )
3282  IF (iprint.EQ.1)
3283  & WRITE(phlun,*) " M:",ppho(5,1),"=>",m
3284  ppho(5,1)=m
3285  ENDIF
3286 
3287  CALL phlupa(25)
3288 
3289  END
3290 
3291 
3292 
3293  FUNCTION phint(IDUM)
3294 C --- can be used with VARIANT A. For B use PHINT1 or 2 --------------
3295 C.----------------------------------------------------------------------
3296 C.
3297 C. PHINT: PHotos universal INTerference correction weight
3298 C.
3299 C. Purpose: calculates correction weight as expressed by
3300 C formula (17) from CPC 79 (1994), 291.
3301 C.
3302 C. Input Parameters: Common /PHOEVT/, with photon added.
3303 C.
3304 C. Output Parameters: correction weight
3305 C.
3306 C. Author(s): Z. Was, P.Golonka Created at: 19/01/05
3307 C. Last Update: 25/01/05
3308 C.
3309 C.----------------------------------------------------------------------
3310  IMPLICIT NONE
3311  REAL*8 phint,phint2
3312  INTEGER idum
3313  INTEGER nmxpho
3314  parameter(nmxpho=10000)
3315  INTEGER idpho,istpho,jdapho,jmopho,nevpho,npho
3316  REAL*8 ppho,vpho
3317  common/phoevt/nevpho,npho,istpho(nmxpho),idpho(nmxpho),
3318  &jmopho(2,nmxpho),jdapho(2,nmxpho),ppho(5,nmxpho),vpho(4,nmxpho)
3319  INTEGER i,k,l
3320  DOUBLE PRECISION emu,mchren,beta,costhg,mpasqr,xph, xc1, xc2,xdeno
3321  DOUBLE PRECISION xnum1,xnum2
3322  DOUBLE PRECISION eps1(4),eps2(4),ph(4),pl(4)
3323  REAL*8 phocha
3324 C--
3325 
3326 C Calculate polarimetric vector: ph, eps1, eps2 are orthogonal
3327 
3328  DO k=1,4
3329  ph(k)=ppho(k,npho)
3330  eps2(k)=1d0
3331  ENDDO
3332 
3333  CALL phoeps(ph,eps2,eps1)
3334  CALL phoeps(ph,eps1,eps2)
3335 
3336 
3337  xnum1=0d0
3338  xnum2=0d0
3339  xdeno=0d0
3340 
3341  DO k=jdapho(1,1),npho-1 ! or JDAPHO(1,2)
3342 
3343 C momenta of charged particle in PL
3344  DO l=1,4
3345  pl(l)=ppho(l,k)
3346  ENDDO
3347 C scalar products: epsilon*p/k*p
3348 
3349  xc1 = - phocha(idpho(k)) *
3350  & ( pl(1)*eps1(1) + pl(2)*eps1(2) + pl(3)*eps1(3) ) /
3351  & ( ph(4)*pl(4) - ph(1)*pl(1) - ph(2)*pl(2) - ph(3)*pl(3) )
3352 
3353  xc2 = - phocha(idpho(k)) *
3354  & ( pl(1)*eps2(1) + pl(2)*eps2(2) + pl(3)*eps2(3) ) /
3355  & ( ph(4)*pl(4) - ph(1)*pl(1) - ph(2)*pl(2) - ph(3)*pl(3) )
3356 
3357 
3358 C accumulate the currents
3359  xnum1 = xnum1+xc1
3360  xnum2 = xnum2+xc2
3361 
3362  xdeno = xdeno + xc1**2 + xc2**2
3363 
3364  ENDDO
3365 
3366  phint=(xnum1**2 + xnum2**2) / xdeno
3367  phint2=phint
3368 
3369  END
3370 
3371 
3372  SUBROUTINE phoeps (VEC1, VEC2, EPS)
3373 C.----------------------------------------------------------------------
3374 C.
3375 C. PHOEPS: PHOeps vector product (normalized to unity)
3376 C.
3377 C. Purpose: calculates vector product, then normalizes its length.
3378 C used to generate orthogonal vectors, i.e. to
3379 C generate polarimetric vectors for photons.
3380 C.
3381 C. Input Parameters: VEC1,VEC2 - input 4-vectors
3382 C.
3383 C. Output Parameters: EPS - normalized 4-vector, orthogonal to
3384 C VEC1 and VEC2
3385 C.
3386 C. Author(s): Z. Was, P.Golonka Created at: 19/01/05
3387 C. Last Update: 25/01/05
3388 C.
3389 C.----------------------------------------------------------------------
3390 
3391  DOUBLE PRECISION vec1(4), vec2(4), eps(4),xn
3392 
3393  eps(1)=vec1(2)*vec2(3) - vec1(3)*vec2(2)
3394  eps(2)=vec1(3)*vec2(1) - vec1(1)*vec2(3)
3395  eps(3)=vec1(1)*vec2(2) - vec1(2)*vec2(1)
3396  eps(4)=0d0
3397 
3398  xn=sqrt( eps(1)**2 +eps(2)**2 +eps(3)**2)
3399 
3400  eps(1)=eps(1)/xn
3401  eps(2)=eps(2)/xn
3402  eps(3)=eps(3)/xn
3403 
3404 
3405  END
3406  SUBROUTINE phodmp
3407 C.----------------------------------------------------------------------
3408 C.
3409 C. PHOTOS: PHOton radiation in decays event DuMP routine
3410 C.
3411 C. Purpose: Print event record.
3412 C.
3413 C. Input Parameters: Common /HEPEVT/
3414 C.
3415 C. Output Parameters: None
3416 C.
3417 C. Author(s): B. van Eijk Created at: 05/06/90
3418 C. Last Update: 05/06/90
3419 C.
3420 C.----------------------------------------------------------------------
3421 C IMPLICIT NONE
3422  DOUBLE PRECISION sumvec(5)
3423  INTEGER i,j
3424 #include "../include/HEPEVT.h"
3425  INTEGER phlun
3426  common/pholun/phlun
3427  DO 10 i=1,5
3428  10 sumvec(i)=0.
3429 C--
3430 C-- Print event number...
3431  WRITE(phlun,9000)
3432  WRITE(phlun,9010) nevhep
3433  WRITE(phlun,9080)
3434  WRITE(phlun,9020)
3435  DO 30 i=1,nhep
3436 C--
3437 C-- For 'stable particle' calculate vector momentum sum
3438  IF (jdahep(1,i).EQ.0) THEN
3439  DO 20 j=1,4
3440  20 sumvec(j)=sumvec(j)+phep(j,i)
3441  IF (jmohep(2,i).EQ.0) THEN
3442  WRITE(phlun,9030) i,idhep(i),jmohep(1,i),(phep(j,i),j=1,5)
3443  ELSE
3444  WRITE(phlun,9040) i,idhep(i),jmohep(1,i),jmohep(2,i),(phep
3445  & (j,i),j=1,5)
3446  ENDIF
3447  ELSE
3448  IF (jmohep(2,i).EQ.0) THEN
3449  WRITE(phlun,9050) i,idhep(i),jmohep(1,i),jdahep(1,i),
3450  & jdahep(2,i),(phep(j,i),j=1,5)
3451  ELSE
3452  WRITE(phlun,9060) i,idhep(i),jmohep(1,i),jmohep(2,i),
3453  & jdahep(1,i),jdahep(2,i),(phep(j,i),j=1,5)
3454  ENDIF
3455  ENDIF
3456  30 CONTINUE
3457  sumvec(5)=sqrt(sumvec(4)**2-sumvec(1)**2-sumvec(2)**2-
3458  &sumvec(3)**2)
3459  WRITE(phlun,9070) (sumvec(j),j=1,5)
3460  RETURN
3461  9000 FORMAT(1h0,80('='))
3462  9010 FORMAT(1h ,29x,'Event No.:',i10)
3463  9020 FORMAT(1h0,1x,'Nr',3x,'Type',3x,'Parent(s)',2x,'Daughter(s)',6x,
3464  &'Px',7x,'Py',7x,'Pz',7x,'E',4x,'Inv. M.')
3465  9030 FORMAT(1h ,i4,i7,3x,i4,9x,'Stable',2x,5f9.2)
3466  9040 FORMAT(1h ,i4,i7,i4,' - ',i4,5x,'Stable',2x,5f9.2)
3467  9050 FORMAT(1h ,i4,i7,3x,i4,6x,i4,' - ',i4,5f9.2)
3468  9060 FORMAT(1h ,i4,i7,i4,' - ',i4,2x,i4,' - ',i4,5f9.2)
3469  9070 FORMAT(1h0,23x,'Vector Sum: ', 5f9.2)
3470  9080 FORMAT(1h0,6x,'Particle Parameters')
3471  END