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