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