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