1 /* copyright(c) 1991-2012 free software foundation, inc.
2 this file is part of the gnu c library.
4 the gnu c library is free software; you can redistribute it and/or
5 modify it under the terms of the gnu lesser general
Public
6 license as published by the free software foundation; either
7 version 2.1 of the license, or(at your option) any later version.
9 the gnu c library is distributed in the hope that it will be useful,
10 but without any warranty; without even the implied warranty of
11 merchantability or fitness for a particular purpose. see the gnu
12 lesser general
Public license for more details.
14 you should have received a copy of the gnu lesser general
Public
15 license along with the gnu c library;
if not, see
16 <http://www.gnu.org/licenses/>. */
19 /* this header is separate from features.h so that the compiler can
20 include it implicitly at the start of every compilation. it must
21 not itself include <features.h> or any other header that includes
22 <features.h> because the
implicit include comes before any feature
23 test macros that may be defined in a source file before it first
24 explicitly includes a system header. gcc knows the name of this
25 header in order to preinclude it. */
27 /* we
do support the iec 559 math functionality,
real and complex. */
29 /* wchar_t uses iso/iec 10646 (2nd ed., published 2011-03-15) /
32 /* we
do not support c11 <threads.h>. */
34 *///////////////////////////////////////////////////////////////////////
38 *// due to short
common block names it may owerwrite variables in other
39 *// parts of the code.
41 *// one should add suffix c_photos_ to names of all commons as soon as
43 *///////////////////////////////////////////////////////////////////////
45 c.----------------------------------------------------------------------
47 c. photos: photos cde-s
49 c. purpose: keep definitions for photos qed correction monte carlo.
51 c. input parameters: none
53 c. output parameters: none
55 c. author(s): z. was, b. van eijk created at: 29/11/89
56 c. last update: 10/08/93
58 c. =========================================================
59 c. general structure information: =
60 c. =========================================================
67 c. 2) general interface:
74 c. photwo(specific
interface
77 c. phtype(specific
interface
78 c. phomak(specific
interface
79 c. 3) qed photon generation:
106 c. name used in sect.
# OF OCC. Comment
107 c. phoqed 1) 2) 3 flags whether emisson to be gen.
108 c. pholun 1) 4) 6 output device number
109 c. phocop 1) 3) 4 photon coupling & min energy
110 c. phpico 1) 3) 4) 5 pi & 2*pi
111 c. phseed 1) 4) 3 rn seed
112 c. phosta 1) 4) 3 status information
113 c. phokey 1) 2) 3) 7 keys for nonstandard application
114 c. phover 1) 1 version info for outside
115 c. hepevt 2) 2 pdg
common
116 c. ph_hepevt2) 8 pdg
common internal
117 c. phoevt 2) 3) 10 pdg branch
118 c. phoif 2) 3) 2 emission flags for pdg branch
119 c. phomom 3) 5 param of char-neutr system
120 c. phophs 3) 5 photon momentum parameters
121 c. phopro 3) 4 var. for photon rep. (in branch)
122 c. phocms 2) 3 parameters of boost to branch cms
123 c. phnum 4) 1 event number from outside
124 c.----------------------------------------------------------------------
126 c.----------------------------------------------------------------------
128 c. photos: photon radiation in decays initialisation
130 c. purpose: initialisation routine for the photos qed radiation
131 c. package. should be called at least once before a call
132 c. to the steering
program 'PHOTOS' is made.
134 c. input parameters: none
136 c. output parameters: none
138 c. author(s): z. was, b. van eijk created at: 26/11/89
139 c. last update: 12/04/90
141 c.----------------------------------------------------------------------
143 INTEGER init,idum,iphqrk,iphekl
147 c--
Return if already initialized...
148 IF (init.NE.0)
RETURN
151 c-- all the following
parameter setters can be called after phoini.
152 c-- initialization of kinematic correction against rounding errors.
153 c-- the set values will be used later
if called wit zero.
154 c-- default
parameter is 1 (no correction) optionally 2, 3, 4
155 c-- in
case of exponentiation new version 5 is needed in most cases.
156 c-- definition given here will be thus overwritten in such a
case
157 c-- below in routine phocin
159 c-- blocks emission from quarks
if parameter is 1 (enables
if 2),
160 c-- physical treatment
161 c-- will be 3, option 2 is not realistic and for tests only,
163 c-- blocks emission in pi0 to gamma e+ e-
if parameter is gt.1
164 c-- (enables otherwise)
167 c-- preset parameters in photos commons
174 c-- initialize marsaglia and zaman random number generator
179 c.----------------------------------------------------------------------
181 c. photos: photon
Common initialisation
183 c. purpose: initialisation of parameters in
common blocks.
185 c. input parameters: none
187 c. output parameters: commons /pholun/, /phopho/, /phocop/, /phpico/
190 c. author(s): b. van eijk created at: 26/11/89
191 c. z. was last update: 29/01/05
193 c.----------------------------------------------------------------------
196 parameter(d_h_nmxhep=10000)
198 common/phoqed/qedrad(d_h_nmxhep)
202 common/phocop/alpha,xphcut
204 common/phpico/pi,twopi
205 INTEGER iseed,i97,j97
206 REAL*8 uran,cran,cdran,cmran
207 common/phseed/iseed(2),i97,j97,uran(97),cran,cdran,cmran
211 common/phosta/status(phomes)
212 LOGICAL interf,isec,itre,iexp,iftop,ifw
213 REAL*8 fint,fsec,expeps
214 COMMON /phokey/ fsec,fint,expeps,interf,isec,itre,iexp,iftop,ifw
219 c--
Return if already initialized...
220 IF (init.NE.0)
RETURN
223 c-- preset switch for photon emission to
'TRUE' for each particle in
224 c-- /ph_hepevt/, this
interface is needed for koralb and koralz...
228 c--
Logical output unit for printing of photos error messages
231 c-- set cut
parameter for photon radiation
235 c-- define some constants
236 alpha=0.00729735039d0
237 pi=3.14159265358979324d0
238 twopi=6.28318530717958648d0
240 c-- default seeds marsaglia and zaman random number generator
244 c-- iitialization for extra options
246 c-- interference weight now universal.
249 c-- second order - double photon switch
251 c-- third/fourth order - triple(or quatric) photon switch,
252 c-- see dipswitch ifour
254 c-- exponentiation on:
265 c-- emision in the hard process g g(q qbar) --> t tbar
269 c-- further initialization done automatically
270 c-- see places with - variant a - variant b - all over
271 c-- to switch between options.
272 c ----------- slower variant a, but stable ------------------
273 c --- it is limiting choice for small xphcut in fixed orer
274 c --- modes of operation
276 c-- best choice is
if fint=2**n
where n+1 is maximal number
277 c-- of charged daughters
278 c-- see report on overweihted events
283 c ----------- faster variant b ------------------
284 c -- it is good for tests of fixed order and small xphcut
285 c -- but is less promising for more
complex cases of interference
286 c -- sometimes fails because of that
293 c----------
END variants a b -----------------------
295 c-- effects of initial state charge(in leptonic w decays)
298 c-- initialise status counter for warning messages
304 c.----------------------------------------------------------------------
306 c. photos: photon radiation in decays general info
308 c. purpose: print photos info
310 c. input parameters: pholun
312 c. output parameters: phovn1, phovn2
314 c. author(s): b. van eijk created at: 12/04/90
315 c. last update: 27/06/04
317 c.----------------------------------------------------------------------
320 INTEGER phovn1,phovn2
321 common/phover/phovn1,phovn2
324 LOGICAL interf,isec,itre,iexp,iftop,ifw
325 REAL*8 fint,fsec,expeps
326 COMMON /phokey/ fsec,fint,expeps,interf,isec,itre,iexp,iftop,ifw
328 common/phocop/alpha,xphcut
330 c-- photos version number and release date
341 WRITE(phlun,9040) iv1,iv2
343 iv2=(phovn2-iv1*10000)/100
344 iv3=phovn2-iv1*10000-iv2*100
345 WRITE(phlun,9050) iv1,iv2,iv3
354 WRITE(phlun,9064) interf,isec,itre,iexp,iftop,ifw,alpha,xphcut
356 IF (interf)
WRITE(phlun,9061)
357 IF (isec)
WRITE(phlun,9062)
358 IF (itre)
WRITE(phlun,9066)
359 IF (iexp)
WRITE(phlun,9067) expeps
360 IF (iftop)
WRITE(phlun,9063)
361 IF (ifw)
WRITE(phlun,9065)
367 9010
FORMAT(1h ,
'*',t81,
'*')
368 9020
FORMAT(1h ,80(
'*'))
369 9030
FORMAT(1h ,
'*',26x,26(
'='),t81,
'*')
370 9040
FORMAT(1h ,
'*',28x,
'PHOTOS, Version: ',i2,
'.',i2,t81,
'*')
371 9050
FORMAT(1h ,
'*',28x,
'Released at: ',i2,
'/',i2,
'/',i2,t81,
'*')
372 9060
FORMAT(1h ,
'*',18x,
'PHOTOS QED Corrections in Particle Decays',
374 9061
FORMAT(1h ,
'*',18x,
'option with interference is active ',
376 9062
FORMAT(1h ,
'*',18x,
'option with double photons is active ',
378 9066
FORMAT(1h ,
'*',18x,
'option with triple/quatric photons is active',
380 9067
FORMAT(1h ,
'*',18x,
'option with exponentiation is active EPSEXP=',
382 9063
FORMAT(1h ,
'*',18x,
'emision in t tbar production is active ',
384 9064
FORMAT(1h ,
'*',18x,
'Internal input parameters:',t81,
'*'
386 &,/, 1h ,
'*',18x,
'INTERF=',l2,
' ISEC=',l2,
' ITRE=',l2,
387 &
' IEXP=',l2,
' IFTOP=',l2,
389 &,/, 1h ,
'*',18x,
'ALPHA_QED=',f8.5,
' XPHCUT=',e8.3,t81,
'*')
390 9065
FORMAT(1h ,
'*',18x,
'correction wt in decay of W is active ',
392 9070
FORMAT(1h ,
'*',9x,
393 &
'Monte Carlo Program - by E. Barberio, B. van Eijk and Z. Was',
395 & 1h ,
'*',9x,
'Version 2.09 - by P. Golonka and Z.W.',t81,
'*')
396 9080
FORMAT( 1h ,
'*',9x,
' ',t81,
'*',/,
398 &
' WARNING (1): /HEPEVT/ is not anymore the standard common block'
400 & 1h ,
'*',9x,
' ',t81,
'*',/,
402 &
' PHOTOS expects /HEPEVT/ to have REAL*8 variables. To change to'
403 & ,t81,
'*',/, 1h ,
'*',9x,
404 &
' REAL*4 modify its declaration in subr. PHOTOS_GET PHOTOS_SET:'
405 & ,t81,
'*',/, 1h ,
'*',9x,
406 &
' REAL*8 d_h_phep, d_h_vhep'
407 & ,t81,
'*',/, 1h ,
'*',9x,
408 &
' WARNING (2): check dims. of /hepevt/ /phoqed/ /ph_hepevt/.'
409 & ,t81,
'*',/, 1h ,
'*',9x,
410 &
' HERE: d_h_nmxhep=10000 and NMXHEP=10000'
413 SUBROUTINE photos(ID)
414 c.----------------------------------------------------------------------
416 c. photos: general search routine + _get + _set
418 c. purpose: /hepevt/ is not anymore a standard at least
419 c.
REAL*8 REAL*4 are in use. photos_get and photos_set
420 c. were to be introduced.
423 c. input parameters: id see routine photos_make
425 c. output parameters: none
427 c. author(s): z. was created at: 21/07/98
428 c. last update: 21/07/98
430 c.----------------------------------------------------------------------
431 COMMON /phlupy/ ipoin,ipoinm
438 IF (1.GT.ipoinm.AND.1.LT.ipoin )
THEN
439 WRITE(phlun,*)
'EVENT NR=',iev,
440 $
'WE ARE TESTING /HEPEVT/ at IPOINT=1 (input)'
446 IF (1.GT.ipoinm.AND.1.LT.ipoin )
THEN
447 WRITE(phlun,*)
'EVENT NR=',iev,
448 $
'WE ARE TESTING /HEPEVT/ at IPOINT=1 (output)'
454 SUBROUTINE photos_get
455 c.----------------------------------------------------------------------
457 c. getter for photos:
459 c. purpose: copies /hepevt/ into /ph_hepevt/
462 c. input parameters: none
464 c. output parameters: none
466 c. author(s): z. was created at: 21/07/98
467 c. last update: 21/07/98
469 c.----------------------------------------------------------------------
473 parameter(d_h_nmxhep=10000)
474 REAL*8 d_h_phep, d_h_vhep
475 INTEGER d_h_nevhep,d_h_nhep,d_h_isthep,d_h_idhep,d_h_jmohep,
480 $ d_h_isthep(d_h_nmxhep),
481 $ d_h_idhep(d_h_nmxhep),
482 $ d_h_jmohep(2,d_h_nmxhep),
483 $ d_h_jdahep(2,d_h_nmxhep),
484 $ d_h_phep(5,d_h_nmxhep),
485 $ d_h_vhep(4,d_h_nmxhep)
486 * ----------------------------------------------------------------------
489 $ d_h_qedrad(d_h_nmxhep)
491 parameter(nmxhep=10000)
492 INTEGER idhep,isthep,jdahep,jmohep,nevhep,nhep
494 common/ph_hepevt/nevhep,nhep,isthep(nmxhep),idhep(nmxhep),
495 &jmohep(2,nmxhep),jdahep(2,nmxhep),phep(5,nmxhep),vhep(4,nmxhep)
497 common/ph_phoqed/qedrad(nmxhep)
502 isthep(k) =d_h_isthep(k)
503 idhep(k) =d_h_idhep(k)
504 jmohep(1,k) =d_h_jmohep(1,k)
505 jdahep(1,k) =d_h_jdahep(1,k)
506 jmohep(2,k) =d_h_jmohep(2,k)
507 jdahep(2,k) =d_h_jdahep(2,k)
509 phep(l,k) =d_h_phep(l,k)
510 vhep(l,k) =d_h_vhep(l,k)
512 phep(5,k) =d_h_phep(5,k)
513 qedrad(k) =d_h_qedrad(k)
518 SUBROUTINE photos_set
519 c.----------------------------------------------------------------------
521 c. setter for photos:
523 c. purpose: copies /ph_hepevt/ into /hepevt/
526 c. input parameters: none
528 c. output parameters: none
530 c. author(s): z. was created at: 21/07/98
531 c. last update: 21/07/98
533 c.----------------------------------------------------------------------
536 parameter(d_h_nmxhep=10000)
537 REAL*8 d_h_phep, d_h_vhep
538 INTEGER d_h_nevhep,d_h_nhep,d_h_isthep,d_h_idhep,d_h_jmohep,
543 $ d_h_isthep(d_h_nmxhep),
544 $ d_h_idhep(d_h_nmxhep),
545 $ d_h_jmohep(2,d_h_nmxhep),
546 $ d_h_jdahep(2,d_h_nmxhep),
547 $ d_h_phep(5,d_h_nmxhep),
548 $ d_h_vhep(4,d_h_nmxhep)
549 * ----------------------------------------------------------------------
552 $ d_h_qedrad(d_h_nmxhep)
554 parameter(nmxhep=10000)
555 INTEGER idhep,isthep,jdahep,jmohep,nevhep,nhep
557 common/ph_hepevt/nevhep,nhep,isthep(nmxhep),idhep(nmxhep),
558 &jmohep(2,nmxhep),jdahep(2,nmxhep),phep(5,nmxhep),vhep(4,nmxhep)
560 common/ph_phoqed/qedrad(nmxhep)
566 d_h_isthep(k) =isthep(k)
567 d_h_idhep(k) =idhep(k)
568 d_h_jmohep(1,k) =jmohep(1,k)
569 d_h_jdahep(1,k) =jdahep(1,k)
570 d_h_jmohep(2,k) =jmohep(2,k)
571 d_h_jdahep(2,k) =jdahep(2,k)
573 d_h_phep(l,k) =phep(l,k)
574 d_h_vhep(l,k) =vhep(l,k)
576 d_h_phep(5,k) =phep(5,k)
577 d_h_qedrad(k) =qedrad(k)
580 SUBROUTINE photos_make(IPARR)
581 c.----------------------------------------------------------------------
583 c. photos_make: general search routine
585 c. purpose: search through the /ph_hepevt/ standard hep
common, sta-
586 c. rting from the ippar-th particle. whenevr branching
587 c. point is found routine phtype(ip) is called.
588 c. finally
if calls on phtype(ip) modified entries,
common
589 c /ph_hepevt/ is ordered.
591 c. input
Parameter: ippar:
Pointer to decaying particle in
592 c. /ph_hepevt/ and the
common itself,
594 c. output parameters:
Common /ph_hepevt/, either with or without
595 c. new particles added.
597 c. author(s): z. was, b. van eijk created at: 26/11/89
598 c. last update: 30/08/93
600 c.----------------------------------------------------------------------
603 INTEGER ip,iparr,ippar,i,j,k,l,nlast
604 DOUBLE PRECISION data
605 INTEGER mother,pospho
608 parameter(nmxhep=10000)
609 INTEGER idhep,isthep,jdahep,jmohep,nevhep,nhep
611 common/ph_hepevt/nevhep,nhep,isthep(nmxhep),idhep(nmxhep),
612 &jmohep(2,nmxhep),jdahep(2,nmxhep),phep(5,nmxhep),vhep(4,nmxhep)
614 common/ph_phoqed/qedrad(nmxhep)
616 parameter(nmxpho=10000)
617 INTEGER istack(0:nmxpho),numit,ntry,kk,ll,ii,na,first,last
618 INTEGER firsta,lasta,ipp,ida1,ida2,mother2,idpho,ispho
619 REAL*8 porig(5,nmxpho)
622 c-- store pointers for cascade treatement...
627 c-- check decay multiplicity and minimum of correctness..
628 IF ((jdahep(1,ip).EQ.0).OR.(jmohep(1,jdahep(1,ip)).NE.ip))
RETURN
630 c-- single branch mode
631 c-- we start looking for the decay points in the cascade
632 c-- ippar is original position
where the
program was called
634 c-- numit denotes number of secondary decay branches
636 c-- ntry denotes number of secondary branches already checked for
637 c-- for existence of further branches
639 c-- let-s search
if iparr does not prevent searching.
642 DO i=jdahep(1,ip),jdahep(2,ip)
643 IF (jdahep(1,i).NE.0.AND.jmohep(1,jdahep(1,i)).EQ.i)
THEN
645 IF (numit.GT.nmxpho)
THEN
647 CALL phoerr(7,
'PHOTOS',data)
652 IF(numit.GT.ntry)
THEN
658 c-- let-s
do generation
661 first=jdahep(1,istack(kk))
662 last=jdahep(2,istack(kk))
665 porig(ll,ii)=phep(ll,first+ii-1)
669 CALL phtype(istack(kk))
671 c-- correct energy/momentum of cascade daughters
677 IF(jmohep(1,ipp).EQ.istack(kk))
678 $ CALL phobos(ipp,porig(1,ii),phep(1,ipp),firsta,lasta)
683 c-- rearrange /ph_hepevt/ to get correct order..
684 IF (nhep.GT.nlast)
THEN
685 DO 160 i=nlast+1,nhep
687 c-- photon mother and position...
689 pospho=jdahep(2,mother)+1
690 c-- intermediate
save of photon energy/momentum and pointers
692 90 photon(j)=phep(j,i)
699 c-- exclude photon in sequence
700 IF (pospho.NE.nhep)
THEN
703 c-- order /ph_hepevt/
704 DO 120 k=i,pospho+1,-1
705 isthep(k)=isthep(k-1)
706 qedrad(k)=qedrad(k-1)
709 jmohep(l,k)=jmohep(l,k-1)
710 100 jdahep(l,k)=jdahep(l,k-1)
712 110 phep(l,k)=phep(l,k-1)
714 120 vhep(l,k)=vhep(l,k-1)
716 c-- correct pointers assuming most dirty /ph_hepevt/...
719 IF ((jmohep(l,k).NE.0).AND.(jmohep(l,k).GE.
720 & pospho)) jmohep(l,k)=jmohep(l,k)+1
721 IF ((jdahep(l,k).NE.0).AND.(jdahep(l,k).GE.
722 & pospho)) jdahep(l,k)=jdahep(l,k)+1
725 c-- store photon energy/momentum
727 140 phep(j,pospho)=photon(j)
730 c-- store pointers for the photon...
731 jdahep(2,mother)=pospho
734 jmohep(1,pospho)=mother
735 jmohep(2,pospho)=mother2
736 jdahep(1,pospho)=ida1
737 jdahep(2,pospho)=ida2
739 c-- get photon production vertex position
741 150 vhep(j,pospho)=vhep(j,pospho-1)
746 SUBROUTINE phobos(IP,PBOOS1,PBOOS2,FIRST,LAST)
747 c.----------------------------------------------------------------------
749 c. phobos: photon radiation in decays boost routine
751 c. purpose: boost particles in cascade decay to parent rest frame
752 c. and boost back with modified boost vector.
754 c. input parameters: ip:
pointer of particle starting chain
756 c. pboos1: boost vector to rest frame,
757 c. pboos2: boost vector to modified frame,
758 c. first:
Pointer to first particle to be boos-
760 c. last:
Pointer to last particle to be boos-
763 c. output parameters:
Common /ph_hepevt/.
765 c. author(s): b. van eijk created at: 13/02/90
766 c. z. was last update: 16/11/93
768 c.----------------------------------------------------------------------
770 DOUBLE PRECISION bet1(3),bet2(3),gam1,gam2,pb,data
771 INTEGER i,j,first,last,maxsta,nstack,ip
772 parameter(maxsta=10000)
773 INTEGER stack(maxsta)
774 REAL*8 pboos1(5),pboos2(5)
776 parameter(nmxhep=10000)
777 INTEGER idhep,isthep,jdahep,jmohep,nevhep,nhep
779 common/ph_hepevt/nevhep,nhep,isthep(nmxhep),idhep(nmxhep),
780 &jmohep(2,nmxhep),jdahep(2,nmxhep),phep(5,nmxhep),vhep(4,nmxhep)
781 IF ((last.EQ.0).OR.(last.LT.first))
RETURN
784 bet1(j)=-pboos1(j)/pboos1(5)
785 10 bet2(j)=pboos2(j)/pboos2(5)
786 gam1=pboos1(4)/pboos1(5)
787 gam2=pboos2(4)/pboos2(5)
789 c-- boost vector to parent rest frame...
790 20
DO 50 i=first,last
791 pb=bet1(1)*phep(1,i)+bet1(2)*phep(2,i)+bet1(3)*phep(3,i)
792 IF (jmohep(1,i).EQ.ip)
THEN
794 30 phep(j,i)=phep(j,i)+bet1(j)*(phep(4,i)+pb/(gam1+1.d0))
795 phep(4,i)=gam1*phep(4,i)+pb
797 c-- ...and boost back to modified parent frame.
798 pb=bet2(1)*phep(1,i)+bet2(2)*phep(2,i)+bet2(3)*phep(3,i)
800 40 phep(j,i)=phep(j,i)+bet2(j)*(phep(4,i)+pb/(gam2+1.d0))
801 phep(4,i)=gam2*phep(4,i)+pb
802 IF (jdahep(1,i).NE.0)
THEN
805 c-- check on stack length...
806 IF (nstack.GT.maxsta)
THEN
808 CALL phoerr(7,
'PHOBOS',data)
814 IF (nstack.NE.0)
THEN
816 c-- now go one step further in the decay tree...
817 first=jdahep(1,stack(nstack))
818 last=jdahep(2,stack(nstack))
825 SUBROUTINE phoin(IP,BOOST,NHEP0)
826 c.----------------------------------------------------------------------
828 c. phoin: photos input
830 c. purpose: copies ip branch of the
common /ph_hepevt/ into /phoevt/
831 c. moves branch into its cms system.
833 c. input parameters: ip:
pointer of particle starting branch
835 c. boost: flag whether boost to cms was or was
838 c. output parameters: commons: /phoevt/, /phocms/
840 c. author(s): z. was created at: 24/05/93
841 c. last update: 16/11/93
843 c.----------------------------------------------------------------------
846 parameter(nmxhep=10000)
847 INTEGER idhep,isthep,jdahep,jmohep,nevhep,nhep
849 common/ph_hepevt/nevhep,nhep,isthep(nmxhep),idhep(nmxhep),
850 &jmohep(2,nmxhep),jdahep(2,nmxhep),phep(5,nmxhep),vhep(4,nmxhep)
852 parameter(nmxpho=10000)
853 INTEGER idpho,istpho,jdapho,jmopho,nevpho,npho
855 common/phoevt/nevpho,npho,istpho(nmxpho),idpho(nmxpho),
856 &jmopho(2,nmxpho),jdapho(2,nmxpho),ppho(5,nmxpho),vpho(4,nmxpho)
857 INTEGER ip,ip2,i,first,last,ll,na
860 DOUBLE PRECISION bet(3),gam,pb
861 COMMON /phocms/ bet,gam
862 LOGICAL interf,isec,itre,iexp,iftop,ifw
863 REAL*8 fint,fsec,expeps
864 COMMON /phokey/ fsec,fint,expeps,interf,isec,itre,iexp,iftop,ifw
866 c let-s calculate
size of the little
common entry
869 npho=3+last-first+nhep-nhep0
871 c let-s take in decaying particle
874 jdapho(2,1)=3+last-first
878 c let-s take in eventual second mother
879 ip2=jmohep(2,jdahep(1,ip))
880 IF((ip2.NE.0).AND.(ip2.NE.ip))
THEN
883 jdapho(2,2)=3+last-first
885 ppho(i,2)=phep(i,ip2)
893 c let-s take in daughters
895 idpho(3+ll)=idhep(first+ll)
896 jmopho(1,3+ll)=jmohep(1,first+ll)
897 IF (jmohep(1,first+ll).EQ.ip) jmopho(1,3+ll)=1
899 ppho(i,3+ll)=phep(i,first+ll)
902 IF (nhep.GT.nhep0)
THEN
903 c let-s take in illegitimate daughters
906 idpho(na+ll)=idhep(nhep0+ll)
907 jmopho(1,na+ll)=jmohep(1,nhep0+ll)
908 IF (jmohep(1,nhep0+ll).EQ.ip) jmopho(1,na+ll)=1
910 ppho(i,na+ll)=phep(i,nhep0+ll)
913 c-- there is nhep-nhep0 daugters more.
914 jdapho(2,1)=3+last-first+nhep-nhep0
916 IF(idpho(npho).EQ.22)CALL phlupa(100001)
919 IF(idpho(npho).EQ.22)CALL phlupa(100002)
920 c special
case of t tbar production process
921 IF(iftop) CALL photwo(0)
923 c-- check whether parent is in its rest frame...
925 c bug reported by vladimir savinov localized and fixed.
926 c protection against rounding error was back-firing
if soft
927 c momentum of mother was physical. consequence was that phcork was
928 c messing up masses of final state particles in vertex of the decay.
929 c only configurations with previously generated photons of energy fraction
930 c smaller than 0.0001 were affected. effect was numerically insignificant.
932 c
IF ( (abs(ppho(4,1)-ppho(5,1)).GT.ppho(5,1)*1.d-8)
933 c $ .AND.(ppho(5,1).NE.0))
THEN
935 IF ((abs(ppho(1,1)+abs(ppho(2,1))+abs(ppho(3,1))).GT.
936 $ ppho(5,1)*1.d-8).AND.(ppho(5,1).NE.0))
THEN
940 c-- boost daughter particles to rest frame of parent...
941 c-- resultant neutral system already calculated in rest frame
943 10 bet(j)=-ppho(j,1)/ppho(5,1)
944 gam=ppho(4,1)/ppho(5,1)
945 DO 30 i=jdapho(1,1),jdapho(2,1)
946 pb=bet(1)*ppho(1,i)+bet(2)*ppho(2,i)+bet(3)*ppho(3,i)
948 20 ppho(j,i)=ppho(j,i)+bet(j)*(ppho(4,i)+pb/(gam+1.d0))
949 30 ppho(4,i)=gam*ppho(4,i)+pb
950 c-- finally boost mother as well
952 pb=bet(1)*ppho(1,i)+bet(2)*ppho(2,i)+bet(3)*ppho(3,i)
954 ppho(j,i)=ppho(j,i)+bet(j)*(ppho(4,i)+pb/(gam+1.d0))
956 ppho(4,i)=gam*ppho(4,i)+pb
958 c special
case of t tbar production process
959 IF(iftop) CALL photwo(1)
961 IF(idpho(npho).EQ.22) CALL phlupa(10000)
964 SUBROUTINE photwo(MODE)
965 c.----------------------------------------------------------------------
967 c. photwo: photos but two mothers allowed
969 c. purpose: combines two mothers into one in /phoevt/
970 c. necessary eg in
case of g g(q qbar) --> t tbar
972 c. input parameters:
Common /phoevt/ (/phocms/)
974 c. output parameters:
Common /phoevt/, (stored mothers)
976 c. author(s): z. was created at: 5/08/93
977 c. last update:10/08/93
979 c.----------------------------------------------------------------------
982 parameter(nmxpho=10000)
983 INTEGER idpho,istpho,jdapho,jmopho,nevpho,npho
985 common/phoevt/nevpho,npho,istpho(nmxpho),idpho(nmxpho),
986 &jmopho(2,nmxpho),jdapho(2,nmxpho),ppho(5,nmxpho),vpho(4,nmxpho)
987 DOUBLE PRECISION bet(3),gam
988 COMMON /phocms/ bet,gam
992 c
logical ifrad is used to tag cases when two mothers may be
993 c merged to the sole one.
994 c so far used in case:
995 c 1) of t tbar production
999 ifrad=(idpho(1).EQ.21).AND.(idpho(2).EQ.21)
1000 ifrad=ifrad.OR.(idpho(1).EQ.-idpho(2).AND.abs(idpho(1)).LE.6)
1002 & .AND.(abs(idpho(3)).EQ.6).AND.(abs(idpho(4)).EQ.6)
1003 mpasqr= (ppho(4,1)+ppho(4,2))**2-(ppho(3,1)+ppho(3,2))**2
1004 & -(ppho(2,1)+ppho(2,2))**2-(ppho(1,1)+ppho(1,2))**2
1005 ifrad=ifrad.AND.(mpasqr.GT.0.0d0)
1007 c.....combining first and second mother
1009 ppho(i,1)=ppho(i,1)+ppho(i,2)
1011 ppho(5,1)=sqrt(mpasqr)
1012 c.....removing second mother,
1018 c boosting of the mothers to the reaction frame not implemented yet.
1019 c to
do it in mode 0 original mothers have to be stored in new comon(?)
1020 c and in mode 1 boosted to cms.
1023 SUBROUTINE phoout(IP,BOOST,NHEP0)
1024 c.----------------------------------------------------------------------
1026 c. phoout: photos output
1028 c. purpose: copies back ip branch of the
common /ph_hepevt/ from
1029 c. /phoevt/ moves branch back from its cms system.
1031 c. input parameters: ip:
pointer of particle starting branch
1032 c. to be given back.
1033 c. boost: flag whether boost to cms was or was
1036 c. output parameters:
Common /phoevt/,
1038 c. author(s): z. was created at: 24/05/93
1041 c.----------------------------------------------------------------------
1044 parameter(nmxhep=10000)
1045 INTEGER idhep,isthep,jdahep,jmohep,nevhep,nhep
1047 common/ph_hepevt/nevhep,nhep,isthep(nmxhep),idhep(nmxhep),
1048 &jmohep(2,nmxhep),jdahep(2,nmxhep),phep(5,nmxhep),vhep(4,nmxhep)
1050 parameter(nmxpho=10000)
1051 INTEGER idpho,istpho,jdapho,jmopho,nevpho,npho
1053 common/phoevt/nevpho,npho,istpho(nmxpho),idpho(nmxpho),
1054 &jmopho(2,nmxpho),jdapho(2,nmxpho),ppho(5,nmxpho),vpho(4,nmxpho)
1055 INTEGER ip,ll,first,last,i
1057 INTEGER nn,j,k,nhep0,na
1058 DOUBLE PRECISION bet(3),gam,pb
1059 COMMON /phocms/ bet,gam
1060 IF(npho.EQ.nevpho)
RETURN
1061 c-- when parent was not in its rest-frame, boost back...
1064 DO 110 j=jdapho(1,1),jdapho(2,1)
1065 pb=-bet(1)*ppho(1,j)-bet(2)*ppho(2,j)-bet(3)*ppho(3,j)
1067 100 ppho(k,j)=ppho(k,j)-bet(k)*(ppho(4,j)+pb/(gam+1.d0))
1068 110 ppho(4,j)=gam*ppho(4,j)+pb
1069 c-- ...boost photon, or whatever
else has shown up
1071 pb=-bet(1)*ppho(1,nn)-bet(2)*ppho(2,nn)-bet(3)*ppho(3,nn)
1073 120 ppho(k,nn)=ppho(k,nn)-bet(k)*(ppho(4,nn)+pb/(gam+1.d0))
1074 ppho(4,nn)=gam*ppho(4,nn)+pb
1079 c let-s take in original daughters
1081 idhep(first+ll) = idpho(3+ll)
1083 phep(i,first+ll) = ppho(i,3+ll)
1086 c let-s take newcomers to the
end of hepevt.
1089 idhep(nhep0+ll) = idpho(na+ll)
1090 isthep(nhep0+ll)=istpho(na+ll)
1091 jmohep(1,nhep0+ll)=ip
1092 jmohep(2,nhep0+ll)=jmohep(2,jdahep(1,ip))
1093 jdahep(1,nhep0+ll)=0
1094 jdahep(2,nhep0+ll)=0
1096 phep(i,nhep0+ll) = ppho(i,na+ll)
1099 nhep=nhep+npho-nevpho
1102 SUBROUTINE phochk(JFIRST)
1103 c.----------------------------------------------------------------------
1105 c. phochk: checking branch.
1107 c. purpose: checks whether particles in the
common block /phoevt/
1108 c. can be served by phomak.
1109 c. jfirst is the position in /ph_hepevt/ (
1110 c. daughter of sub-branch under action.
1113 c. author(s): z. was created at: 22/10/92
1114 c. last update: 11/12/00
1116 c.----------------------------------------------------------------------
1117 c ********************
1120 parameter(nmxpho=10000)
1121 INTEGER idpho,istpho,jdapho,jmopho,nevpho,npho
1123 common/phoevt/nevpho,npho,istpho(nmxpho),idpho(nmxpho),
1124 &jmopho(2,nmxpho),jdapho(2,nmxpho),ppho(5,nmxpho),vpho(4,nmxpho)
1126 common/phoif/chkif(nmxpho)
1128 parameter(nmxhep=10000)
1130 common/ph_phoqed/qedrad(nmxhep)
1133 INTEGER idabs,nlast,i,ippar
1134 LOGICAL interf,isec,itre,iexp,iftop,ifw,ifnpi0,ifkl
1135 REAL*8 fint,fsec,expeps
1136 COMMON /phokey/ fsec,fint,expeps,interf,isec,itre,iexp,iftop,ifw
1138 INTEGER ident,k,iqrk,iphqrk,iekl,iphekl
1139 c these are ok ....
if you
do not like somebody
else, add here.
1141 & ( ((idabs.GT.9.OR.iqrk.NE.1).AND.(idabs.LE.40))
1142 & .OR.(idabs.GT.100) )
1143 & .AND.(idabs.NE.21)
1144 $ .AND.(idabs.NE.2101).AND.(idabs.NE.3101).AND.(idabs.NE.3201)
1145 & .AND.(idabs.NE.1103).AND.(idabs.NE.2103).AND.(idabs.NE.2203)
1146 & .AND.(idabs.NE.3103).AND.(idabs.NE.3203).AND.(idabs.NE.3303)
1153 c checking for good particles
1157 ifnpi0= (idpho(1).NE.111)
1158 ifkl = ((idpho(1).EQ.130).AND.
1159 $ ((idpho(3).EQ.22).OR.(idpho(4).EQ.22).OR.
1160 $ (idpho(5).EQ.22)).AND.
1161 $ ((idpho(3).EQ.11).OR.(idpho(4).EQ.11).OR.
1162 $ (idpho(5).EQ.11)) )
1164 ifnpi0=(ifnpi0.AND.(.NOT.ifkl))
1167 idabs = abs(idpho(i))
1168 c possibly call on phzode is a dead(to be omitted) code.
1169 chkif(i)= f(idabs) .AND.f(abs(idpho(1)))
1170 & .AND. (idpho(2).EQ.0)
1171 IF(i.GT.2) chkif(i)=chkif(i).AND.qedrad(jfirst+i-ippar-2)
1175 c now we go to special cases,
where chkif(i) will be overwritten
1178 c special
case of top pair production
1179 DO k=jdapho(2,1),jdapho(1,1),-1
1180 IF(idpho(k).NE.22)
THEN
1186 ifrad=((idpho(1).EQ.21).AND.(idpho(2).EQ.21))
1187 & .OR. ((abs(idpho(1)).LE.6).AND.((idpho(2)).EQ.(-idpho(1))))
1189 & .AND.(abs(idpho(3)).EQ.6).AND.((idpho(4)).EQ.(-idpho(3)))
1194 IF(i.GT.2) chkif(i)=chkif(i).AND.qedrad(jfirst+i-ippar-2)
1201 c special
case of top decay
1202 DO k=jdapho(2,1),jdapho(1,1),-1
1203 IF(idpho(k).NE.22)
THEN
1209 ifrad=((abs(idpho(1)).EQ.6).AND.(idpho(2).EQ.0))
1211 & .AND.((abs(idpho(3)).EQ.24).AND.(abs(idpho(4)).EQ.5)
1212 & .OR.(abs(idpho(3)).EQ.5).AND.(abs(idpho(4)).EQ.24))
1217 IF(i.GT.2) chkif(i)=chkif(i).AND.qedrad(jfirst+i-ippar-2)
1224 SUBROUTINE phtype(ID)
1225 c.----------------------------------------------------------------------
1227 c. phtype: central manadgement routine.
1229 c. purpose: defines what kind of the
1230 c. actions will be performed at point id.
1232 c. input parameters: id:
pointer of particle starting branch
1233 c. in /ph_hepevt/ to be treated.
1235 c. output parameters:
Common /ph_hepevt/.
1237 c. author(s): z. was created at: 24/05/93
1238 c. p. golonka last update: 27/06/04
1240 c.----------------------------------------------------------------------
1243 parameter(nmxhep=10000)
1244 INTEGER idhep,isthep,jdahep,jmohep,nevhep,nhep
1246 common/ph_hepevt/nevhep,nhep,isthep(nmxhep),idhep(nmxhep),
1247 &jmohep(2,nmxhep),jdahep(2,nmxhep),phep(5,nmxhep),vhep(4,nmxhep)
1248 LOGICAL interf,isec,itre,iexp,iftop,ifw
1249 REAL*8 fint,fsec,expeps
1250 COMMON /phokey/ fsec,fint,expeps,interf,isec,itre,iexp,iftop,ifw
1254 REAL*8 pro,prsum,esu
1255 COMMON /phoexp/ pro(nx),nchan,expini
1259 REAL*8 rn,phoran,sum
1263 ifour=(.true.).AND.(itre)
1266 c-- check decay multiplicity..
1267 IF (jdahep(1,id).EQ.0)
RETURN
1268 c
IF (jdahep(1,id).EQ.jdahep(2,id))
RETURN
1280 CALL phomak(id,nhep0)
1294 IF(rn.LT.sum) goto 100
1298 CALL phomak(id,nhep0)
1299 IF(sum.GT.1d0-expeps) goto 100
1303 c-- quatro photon emission
1306 IF (rn.GE.23.d0/24d0)
THEN
1307 CALL phomak(id,nhep0)
1308 CALL phomak(id,nhep0)
1309 CALL phomak(id,nhep0)
1310 CALL phomak(id,nhep0)
1311 ELSEIF (rn.GE.17.d0/24d0)
THEN
1312 CALL phomak(id,nhep0)
1313 CALL phomak(id,nhep0)
1314 ELSEIF (rn.GE.9.d0/24d0)
THEN
1315 CALL phomak(id,nhep0)
1318 c-- triple photon emission
1321 IF (rn.GE.5.d0/6d0)
THEN
1322 CALL phomak(id,nhep0)
1323 CALL phomak(id,nhep0)
1324 CALL phomak(id,nhep0)
1325 ELSEIF (rn.GE.2.d0/6d0)
THEN
1326 CALL phomak(id,nhep0)
1329 c-- double photon emission
1332 IF (rn.GE.0.5d0)
THEN
1333 CALL phomak(id,nhep0)
1334 CALL phomak(id,nhep0)
1337 c-- single photon emission
1339 CALL phomak(id,nhep0)
1342 c-- electron positron pair(coomented out for a
while
1343 c
IF (ipair) CALL phopar(id,nhep0)
1345 SUBROUTINE phomak(IPPAR,NHEP0)
1346 c.----------------------------------------------------------------------
1348 c. phomak: photos make
1350 c. purpose: single or double bremstrahlung radiative corrections
1351 c. are generated in the decay of the ippar-th particle in
1352 c. the hep
common /ph_hepevt/. example of the
use of
1355 c. input
Parameter: ippar:
Pointer to decaying particle in
1356 c. /ph_hepevt/ and the
common itself
1358 c. output parameters:
Common /ph_hepevt/, either with or without
1361 c. author(s): z. was, created at: 26/05/93
1362 c. last update: 29/01/05
1364 c.----------------------------------------------------------------------
1366 DOUBLE PRECISION data
1368 INTEGER ip,ippar,ncharg
1369 INTEGER wtdum,idum,nhep0
1370 INTEGER ncharb,neudau
1374 parameter(nmxhep=10000)
1375 INTEGER idhep,isthep,jdahep,jmohep,nevhep,nhep
1377 common/ph_hepevt/nevhep,nhep,isthep(nmxhep),idhep(nmxhep),
1378 &jmohep(2,nmxhep),jdahep(2,nmxhep),phep(5,nmxhep),vhep(4,nmxhep)
1379 LOGICAL interf,isec,itre,iexp,iftop,ifw
1380 REAL*8 fint,fsec,expeps
1381 COMMON /phokey/ fsec,fint,expeps,interf,isec,itre,iexp,iftop,ifw
1387 CALL phoin(ip,boost,nhep0)
1388 CALL phochk(jdahep(1,ip))
1390 CALL phopre(1,wt,neudau,ncharb)
1392 IF (wt.EQ.0.0d0)
RETURN
1394 c phodo is caling phoran, thus change of series
if it is moved before
if
1395 CALL phodo(1,ncharb,neudau)
1396 c we eliminate /fint in variant b.
1397 IF (interf) wt=wt*phint(idum) /fint
1398 IF (ifw) CALL phobw(wt)
1400 IF (wt.GT.1.0d0) CALL phoerr(3,
'WT_INT',data)
1403 CALL phoout(ip,boost,nhep0)
1407 FUNCTION phint1(IDUM)
1408 c.----------------------------------------------------------------------
1410 c. phint: photos interference(old version kept for tests only.
1412 c. purpose: calculates interference between emission of photons from
1413 c. different possible chaged daughters stored in
1414 c. the hep
common /phoevt/.
1416 c. input
Parameter: commons /phoevt/ /phomom/ /phophs/
1419 c. output parameters:
1422 c. author(s): z. was, created at: 10/08/93
1423 c. last update: 15/03/99
1425 c.----------------------------------------------------------------------
1432 parameter(nmxpho=10000)
1433 INTEGER idpho,istpho,jdapho,jmopho,nevpho,npho
1435 common/phoevt/nevpho,npho,istpho(nmxpho),idpho(nmxpho),
1436 &jmopho(2,nmxpho),jdapho(2,nmxpho),ppho(5,nmxpho),vpho(4,nmxpho)
1437 DOUBLE PRECISION mchsqr,mnesqr
1439 common/phomom/mchsqr,mnesqr,pneutr(5)
1440 DOUBLE PRECISION costhg,sinthg
1441 REAL*8 xphmax,xphoto
1442 common/phophs/xphmax,xphoto,costhg,sinthg
1443 REAL*8 mpasqr,xx,beta
1447 DO k=jdapho(2,1),jdapho(1,1),-1
1448 IF(idpho(k).NE.22)
THEN
1454 c check
if there is a photon
1455 ifint= npho.GT.ident
1456 c check
if it is two body + gammas reaction
1457 ifint= ifint.AND.(ident-jdapho(1,1)).EQ.1
1458 c check
if two body was particle antiparticle
1459 ifint= ifint.AND.idpho(jdapho(1,1)).EQ.-idpho(ident)
1460 c check
if particles were charged
1461 ifint= ifint.AND.phocha(idpho(ident)).NE.0
1462 c calculates interference weight contribution
1464 mpasqr = ppho(5,1)**2
1465 xx=4.d0*mchsqr/mpasqr*(1.d0-xphoto)/(1.d0-xphoto+(mchsqr-mnesqr)
1468 phint = 2d0/(1d0+costhg**2*beta**2)
1475 FUNCTION phint2(IDUM)
1476 c.----------------------------------------------------------------------
1478 c. phint: photos interference
1480 c. purpose: calculates interference between emission of photons from
1481 c. different possible chaged daughters stored in
1482 c. the hep
common /phoevt/.
1484 c. input
Parameter: commons /phoevt/ /phomom/ /phophs/
1487 c. output parameters:
1490 c. author(s): z. was, created at: 10/08/93
1493 c.----------------------------------------------------------------------
1496 REAL*8 phint,phint1,phint2
1500 parameter(nmxpho=10000)
1501 INTEGER idpho,istpho,jdapho,jmopho,nevpho,npho
1503 common/phoevt/nevpho,npho,istpho(nmxpho),idpho(nmxpho),
1504 &jmopho(2,nmxpho),jdapho(2,nmxpho),ppho(5,nmxpho),vpho(4,nmxpho)
1505 DOUBLE PRECISION mchsqr,mnesqr
1507 common/phomom/mchsqr,mnesqr,pneutr(5)
1508 DOUBLE PRECISION costhg,sinthg
1509 REAL*8 xphmax,xphoto
1510 common/phophs/xphmax,xphoto,costhg,sinthg
1511 REAL*8 mpasqr,xx,beta,pq1(4),pq2(4),pphot(4)
1512 REAL*8 ss,pp2,pp,e1,e2,q1,q2,costhe
1516 DO k=jdapho(2,1),jdapho(1,1),-1
1517 IF(idpho(k).NE.22)
THEN
1523 c check
if there is a photon
1524 ifint= npho.GT.ident
1525 c check
if it is two body + gammas reaction
1526 ifint= ifint.AND.(ident-jdapho(1,1)).EQ.1
1527 c check
if two body was particle antiparticle(we improve on it
1528 c ifint= ifint.AND.idpho(jdapho(1,1)).EQ.-idpho(ident)
1529 c check
if particles were charged
1530 ifint= ifint.AND.abs(phocha(idpho(ident))).GT.0.01d0
1531 c check
if they have both charge
1533 $ abs(phocha(idpho(jdapho(1,1)))).gt.0.01d0
1534 c calculates interference weight contribution
1536 mpasqr = ppho(5,1)**2
1537 xx=4.d0*mchsqr/mpasqr*(1.-xphoto)/(1.-xphoto+(mchsqr-mnesqr)/
1540 phint = 2d0/(1d0+costhg**2*beta**2)
1541 ss =mpasqr*(1.d0-xphoto)
1542 pp2=((ss-mchsqr-mnesqr)**2-4*mchsqr*mnesqr)/ss/4
1544 e1 =sqrt(pp2+mchsqr)
1545 e2 =sqrt(pp2+mnesqr)
1546 phint= (e1+e2)**2/((e2+costhg*pp)**2+(e1-costhg*pp)**2)
1548 q1=phocha(idpho(jdapho(1,1)))
1549 q2=phocha(idpho(ident))
1551 pq1(k)=ppho(k,jdapho(1,1))
1552 pq2(k)=ppho(k,jdapho(1,1)+1)
1553 pphot(k)=ppho(k,npho)
1555 costhe=(pphot(1)*pq1(1)+pphot(2)*pq1(2)+pphot(3)*pq1(3))
1556 costhe=costhe/sqrt(pq1(1)**2+pq1(2)**2+pq1(3)**2)
1557 costhe=costhe/sqrt(pphot(1)**2+pphot(2)**2+pphot(3)**2)
1563 IF (costhg*costhe.GT.0)
then
1565 phint= (q1*(e2+costhg*pp)-q2*(e1-costhg*pp))**2
1566 & /(q1**2*(e2+costhg*pp)**2+q2**2*(e1-costhg*pp)**2)
1569 phint= (q1*(e1-costhg*pp)-q2*(e2+costhg*pp))**2
1570 & /(q1**2*(e1-costhg*pp)**2+q2**2*(e2+costhg*pp)**2)
1580 SUBROUTINE phopre(IPARR,WT,NEUDAU,NCHARB)
1581 c.----------------------------------------------------------------------
1583 c. photos: photon radiation in decays
1585 c. purpose: order(alpha) radiative corrections are generated in
1586 c. the decay of the ippar-th particle in the hep-like
1587 c.
common /phoevt/. photon radiation takes place from one
1588 c. of the charged daughters of the decaying particle ippar
1589 c. wt is calculated, eventual rejection will be performed
1590 c. later after inclusion of interference weight.
1592 c. input
Parameter: ippar:
Pointer to decaying particle in
1593 c. /phoevt/ and the
common itself,
1595 c. output parameters:
Common /phoevt/, either with or without a
1597 c. wt weight of the configuration
1599 c. author(s): z. was, b. van eijk created at: 26/11/89
1600 c. last update: 29/01/05
1602 c.----------------------------------------------------------------------
1604 DOUBLE PRECISION minmas,mpasqr,mchren
1605 DOUBLE PRECISION beta,eps,del1,del2,data,biglog
1606 REAL*8 phocha,phospi,phoran,phocor,massum
1607 INTEGER ip,iparr,ippar,i,j,me,ncharg,neupoi,nlast,thedum
1609 INTEGER ncharb,neudau
1612 parameter(nmxpho=10000)
1613 INTEGER idpho,istpho,jdapho,jmopho,nevpho,npho
1615 common/phoevt/nevpho,npho,istpho(nmxpho),idpho(nmxpho),
1616 &jmopho(2,nmxpho),jdapho(2,nmxpho),ppho(5,nmxpho),vpho(4,nmxpho)
1618 common/phoif/chkif(nmxpho)
1619 INTEGER chapoi(nmxpho)
1620 DOUBLE PRECISION mchsqr,mnesqr
1622 common/phomom/mchsqr,mnesqr,pneutr(5)
1623 DOUBLE PRECISION costhg,sinthg
1624 REAL*8 xphmax,xphoto
1625 common/phophs/xphmax,xphoto,costhg,sinthg
1627 common/phocop/alpha,xphcut
1629 REAL*8 probh,corwt,xf
1630 common/phopro/probh,corwt,xf,irep
1631 c may be it is not the best place, but ...
1632 LOGICAL interf,isec,itre,iexp,iftop,ifw
1633 REAL*8 fint,fsec,expeps
1634 COMMON /phokey/ fsec,fint,expeps,interf,isec,itre,iexp,iftop,ifw
1638 c-- store pointers for cascade treatement...
1643 c-- check decay multiplicity..
1644 IF (jdapho(1,ip).EQ.0)
RETURN
1646 c-- loop over daughters, determine charge multiplicity
1651 DO 20 i=jdapho(1,ip),jdapho(2,ip)
1654 c-- exclude marked particles, quarks and gluons etc...
1656 IF (chkif(i-jdapho(1,ip)+3))
THEN
1657 IF (phocha(idpho(i)).NE.0)
THEN
1659 IF (ncharg.GT.nmxpho)
THEN
1661 CALL phoerr(1,
'PHOTOS',data)
1665 minmas=minmas+ppho(5,i)**2
1667 massum=massum+ppho(5,i)
1669 IF (ncharg.NE.0)
THEN
1671 c-- check that sum of daughter masses does not exceed parent mass
1672 IF ((ppho(5,ip)-massum)/ppho(5,ip).GT.2.d0*xphcut)
THEN
1674 c-- order charged particles according to decreasing mass, this to
1675 c-- increase efficiency(smallest mass is treated first).
1676 IF (ncharg.GT.1) CALL phooma(1,ncharg,chapoi)
1680 70 pneutr(j)=-ppho(j,chapoi(ncharg))
1681 pneutr(4)=ppho(5,ip)-ppho(4,chapoi(ncharg))
1683 c-- calculate invariant mass of
'neutral' etc. systems
1684 mpasqr=ppho(5,ip)**2
1685 mchsqr=ppho(5,chapoi(ncharg))**2
1686 IF ((jdapho(2,ip)-jdapho(1,ip)).EQ.1)
THEN
1688 IF (neupoi.EQ.chapoi(ncharg)) neupoi=jdapho(2,ip)
1689 mnesqr=ppho(5,neupoi)**2
1690 pneutr(5)=ppho(5,neupoi)
1692 mnesqr=pneutr(4)**2-pneutr(1)**2-pneutr(2)**2-pneutr(3)**2
1693 mnesqr=max(mnesqr,minmas-mchsqr)
1694 pneutr(5)=sqrt(mnesqr)
1697 c-- determine kinematical limit...
1698 xphmax=(mpasqr-(pneutr(5)+ppho(5,chapoi(ncharg)))**2)/mpasqr
1700 c-- photon energy fraction...
1701 CALL phoene(mpasqr,mchren,beta,biglog,idpho(chapoi(ncharg)))
1703 IF (xphoto.LT.-4d0)
THEN
1706 c-- energy fraction not too large(very seldom) ? define angle.
1707 ELSEIF ((xphoto.LT.xphcut).OR.(xphoto.GT.xphmax))
THEN
1709 c-- no radiation was accepted, check for more daughters that may ra-
1710 c-- diate and correct radiation probability...
1712 IF (ncharg.GT.0)
THEN
1718 c-- angle is generated in the frame defined by charged vector and
1719 c-- pneutr, distribution is taken in the infrared limit...
1720 eps=mchren/(1.d0+beta)
1722 c-- calculate sin(theta) and cos(theta) from interval variables
1723 del1=(2.d0-eps)*(eps/(2.d0-eps))**phoran(thedum)
1726 c ----------- variant b ------------------
1727 cc corrections for more efiicient interference correction,
1728 cc instead of doubling crude distribution, we add flat parallel channel
1729 c
IF (phoran(thedum).LT.biglog/beta/(biglog/beta+2*fint))
THEN
1730 c costhg=(1.d0-del1)/beta
1731 c sinthg=sqrt(del1*del2-mchren)/beta
1733 c costhg=-1d0+2*phoran(thedum)
1734 c sinthg= sqrt(1d0-costhg**2)
1737 c
IF (fint.GT.1.0d0)
THEN
1739 c wgt=1d0/(1d0-beta*costhg)
1740 c wgt=wgt/(wgt+fint)
1747 c -----------
END of variant b ------------------
1749 c ----------- variant a ------------------
1750 costhg=(1.d0-del1)/beta
1751 sinthg=sqrt(del1*del2-mchren)/beta
1753 c -----------
END of variant a ------------------
1756 c-- determine spin of particle and construct code for matrix element
1757 me=2.d0*phospi(idpho(chapoi(ncharg)))+1.d0
1759 c-- weighting procedure with
'exact' matrix element, reconstruct kine-
1760 c-- matics for photon, neutral and charged system and update /phoevt/.
1761 c-- find
pointer to the first component of
'neutral' system
1762 DO i=jdapho(1,ip),jdapho(2,ip)
1763 IF (i.NE.chapoi(ncharg))
THEN
1769 c--
Pointer not found...
1771 CALL phoerr(5,
'PHOKIN',data)
1773 ncharb=chapoi(ncharg)
1774 ncharb=ncharb-jdapho(1,ip)+3
1775 neudau=neudau-jdapho(1,ip)+3
1776 wt=phocor(mpasqr,mchren,me)*wgt
1779 data=ppho(5,ip)-massum
1780 CALL phoerr(10,
'PHOTOS',data)
1786 SUBROUTINE phooma(IFIRST,ILAST,POINTR)
1787 c.----------------------------------------------------------------------
1789 c. photos: photon radiation in decays order mass vector
1791 c. purpose: order the contents of array
'POINTR' according to the
1792 c. decreasing value in the array
'MASS'.
1794 c. input parameters: ifirst, ilast: pointers to the vector loca-
1796 c. pointr: unsorted array with pointers to
1799 c. output
Parameter: pointr: sorted arrays with respect to
1800 c. particle mass
'PPHO(5,*)'.
1802 c. author(s): b. van eijk created at: 28/11/89
1803 c. last update: 27/05/93
1805 c.----------------------------------------------------------------------
1808 parameter(nmxpho=10000)
1809 INTEGER idpho,istpho,jdapho,jmopho,nevpho,npho
1811 common/phoevt/nevpho,npho,istpho(nmxpho),idpho(nmxpho),
1812 &jmopho(2,nmxpho),jdapho(2,nmxpho),ppho(5,nmxpho),vpho(4,nmxpho)
1813 INTEGER ifirst,ilast,i,j,bufpoi,pointr(nmxpho)
1814 REAL*8 bufmas,mass(nmxpho)
1815 IF (ifirst.EQ.ilast)
RETURN
1817 c-- copy particle masses
1818 DO 10 i=ifirst,ilast
1819 10 mass(i)=ppho(5,pointr(i))
1821 c-- order the masses in a decreasing series
1822 DO 30 i=ifirst,ilast-1
1824 IF (mass(j).LE.mass(i)) goto 20
1835 SUBROUTINE phoene(MPASQR,MCHREN,BETA,BIGLOG,IDENT)
1836 c.----------------------------------------------------------------------
1838 c. photos: photon radiation in decays calculation of photon energy
1841 c. purpose:
Subroutine returns photon energy fraction (in (parent
1842 c. mass)/2 units) for the decay bremsstrahlung.
1844 c. input parameters: mpasqr: mass of decaying system squared,
1845 c. xphcut: minimum energy fraction of photon,
1846 c. xphmax: maximum energy fraction of photon.
1848 c. output
Parameter: mchren: renormalised mass squared,
1849 c. beta: beta factor due to renormalisation,
1850 c. xphoto: photon energy fraction,
1851 c. xf: correction factor for phofac.
1853 c. author(s): s. jadach, z. was created at: 01/01/89
1854 c. b. van eijk, p.golonka last update: 29/01/05
1856 c.----------------------------------------------------------------------
1858 DOUBLE PRECISION mpasqr,mchren,biglog,beta,data
1859 INTEGER iwt1,irn,iwt2
1860 REAL*8 prsoft,prhard,phoran,phofac
1861 DOUBLE PRECISION mchsqr,mnesqr
1864 REAL*8 phocha,prkill,rrr
1865 common/phomom/mchsqr,mnesqr,pneutr(5)
1866 DOUBLE PRECISION costhg,sinthg
1867 REAL*8 xphmax,xphoto
1868 common/phophs/xphmax,xphoto,costhg,sinthg
1870 common/phocop/alpha,xphcut
1872 common/phpico/pi,twopi
1874 REAL*8 probh,corwt,xf
1875 common/phopro/probh,corwt,xf,irep
1876 LOGICAL interf,isec,itre,iexp,iftop,ifw
1877 REAL*8 fint,fsec,expeps
1878 COMMON /phokey/ fsec,fint,expeps,interf,isec,itre,iexp,iftop,ifw
1883 COMMON /phoexp/ pro(nx),nchan,expini
1885 IF (xphmax.LE.xphcut)
THEN
1890 c-- probabilities for hard and soft bremstrahlung...
1891 mchren=4.d0*mchsqr/mpasqr/(1.d0+mchsqr/mpasqr)**2
1892 beta=sqrt(1.d0-mchren)
1894 c ----------- variant b ------------------
1895 cc we replace 1d0/beta*biglog with(1d0/beta*biglog+2*fint)
1896 cc for integral of new crude
1897 c biglog=log(mpasqr/mchsqr*(1.d0+beta)**2/4.d0*
1898 c & (1.d0+mchsqr/mpasqr)**2)
1899 c prhard=alpha/pi*(1d0/beta*biglog+2*fint)*(log(xphmax/xphcut)
1900 c &-.75d0+xphcut/xphmax-.25d0*xphcut**2/xphmax**2)
1901 c prhard=prhard*phocha(ident)**2*fsec
1902 c -----------
END of variant b ------------------
1904 c ----------- variant a ------------------
1905 biglog=log(mpasqr/mchsqr*(1.d0+beta)**2/4.d0*
1906 & (1.d0+mchsqr/mpasqr)**2)
1907 prhard=alpha/pi*(1d0/beta*biglog)*
1908 &(log(xphmax/xphcut)-.75d0+xphcut/xphmax-.25d0*xphcut**2/xphmax**2)
1909 prhard=prhard*phocha(ident)**2*fsec*fint
1910 c -----------
END of variant a ------------------
1911 IF (irep.EQ.0) probh=0.d0
1916 pro(nchan)=prhard+0.05*(1.0+fint)
1929 prkill=pro(nchan)/prsum-prhard
1934 prhard=prhard*phofac(0)
1941 c-- check on kinematical bounds
1943 IF (prsoft.LT.-5.0d-8)
THEN
1945 CALL phoerr(2,
'PHOENE',data)
1948 IF (prsoft.LT.0.1d0)
THEN
1950 CALL phoerr(2,
'PHOENE',data)
1955 IF (rrr.LT.prsoft)
THEN
1957 c-- no photon... (ie. photon too soft)
1959 IF (rrr.LT.prkill) xphoto=-5d0
1962 c-- hard photon... (ie. photon hard enough).
1963 c-- calculate altarelli-parisi kernel
1964 10 xphoto=exp(phoran(irn)*log(xphcut/xphmax))
1965 xphoto=xphoto*xphmax
1966 IF (phoran(iwt2).GT.((1.d0+(1.d0-xphoto/xphmax)**2)/2.d0))
1970 c-- calculate
parameter for phofac function
1971 xf=4.d0*mchsqr*mpasqr/(mpasqr+mchsqr-mnesqr)**2
1974 FUNCTION phocor(MPASQR,MCHREN,ME)
1975 c.----------------------------------------------------------------------
1977 c. photos: photon radiation in decays correction weight from
1980 c. purpose: calculate photon angle. the reshaping functions will
1981 c. have to depend on the spin s of the charged particle.
1982 c. we define: me = 2 * s + 1
1984 c. input parameters: mpasqr: parent mass squared,
1985 c. mchren: renormalised mass of charged system,
1986 c. me: 2 * spin + 1 determines matrix element
1988 c. output
Parameter:
Function value.
1990 c. author(s): z. was, b. van eijk created at: 26/11/89
1991 c. last update: 21/03/93
1993 c.----------------------------------------------------------------------
1995 DOUBLE PRECISION mpasqr,mchren,beta,xx,yy,data
1997 REAL*8 phocor,phofac,wt1,wt2,wt3
1998 DOUBLE PRECISION mchsqr,mnesqr
2000 common/phomom/mchsqr,mnesqr,pneutr(5)
2001 DOUBLE PRECISION costhg,sinthg
2002 REAL*8 xphmax,xphoto
2003 common/phophs/xphmax,xphoto,costhg,sinthg
2005 REAL*8 probh,corwt,xf
2006 common/phopro/probh,corwt,xf,irep
2008 c-- shaping(modified by zw)...
2009 xx=4.d0*mchsqr/mpasqr*(1.d0-xphoto)/(1.d0-xphoto+(mchsqr-mnesqr)/
2013 wt3=(1.d0-xphoto/xphmax)/((1.d0+(1.d0-xphoto/xphmax)**2)/2.d0)
2014 ELSEIF (me.EQ.2)
THEN
2015 yy=0.5d0*(1.d0-xphoto/xphmax+1.d0/(1.d0-xphoto/xphmax))
2017 ELSEIF ((me.EQ.3).OR.(me.EQ.4).OR.(me.EQ.5))
THEN
2019 wt3=(1.d0+(1.d0-xphoto/xphmax)**2-(xphoto/xphmax)**3)/
2020 & (1.d0+(1.d0-xphoto/xphmax)** 2)
2023 CALL phoerr(6,
'PHOCOR',data)
2028 wt1=(1.d0-costhg*sqrt(1.d0-mchren))/(1.d0-costhg*beta)
2029 wt2=(1.d0-xx/yy/(1.d0-beta**2*costhg**2))*(1.d0+costhg*beta)/2.d0
2033 IF (phocor.GT.1.d0)
THEN
2035 CALL phoerr(3,
'PHOCOR',data)
2039 FUNCTION phofac(MODE)
2040 c.----------------------------------------------------------------------
2042 c. photos: photon radiation in decays control factor
2044 c. purpose: this is the control
function for the photon spectrum and
2045 c. final weighting. it is called from phoene for genera-
2046 c. ting the raw photon energy spectrum(mode=0) and in pho-
2047 c. cor to scale the final weight(mode=1). the factor con-
2048 c. sists of 3 terms. addition of the factor ff which mul-
2049 c. tiplies phofac for mode=0 and divides phofac for mode=1,
2050 c. does not affect the results for the mc generation. an
2051 c. appropriate choice for ff can speed up the calculation.
2052 c. note that a too small value of ff may cause weight over-
2053 c. flow in phocor and will generate a warning, halting the
2054 c. execution. prx should be included for repeated calls
2055 c. for the same event, allowing more particles to radiate
2056 c. photons. at the first call irep=0, for more than 1
2057 c. charged decay products, irep >= 1. thus, prsoft(no
2058 c. photon radiation probability in the previous calls)
2059 c. appropriately scales the strength of the bremsstrahlung.
2061 c. input parameters: mode, probh, xf
2063 c. output
Parameter:
Function value
2065 c. author(s): s. jadach, z. was created at: 01/01/89
2066 c. b. van eijk, p.golonka last update: 26/06/04
2068 c.----------------------------------------------------------------------
2070 REAL*8 phofac,ff,prx
2073 REAL*8 probh,corwt,xf
2074 common/phopro/probh,corwt,xf,irep
2075 LOGICAL interf,isec,itre,iexp,iftop,ifw
2076 REAL*8 fint,fsec,expeps
2077 COMMON /phokey/ fsec,fint,expeps,interf,isec,itre,iexp,iftop,ifw
2079 DATA prx,ff/ 0.d0, 0.d0/
2084 IF (mode.EQ.-1)
THEN
2088 ELSEIF (mode.EQ.0)
THEN
2089 IF (irep.EQ.0) prx=1.d0
2090 prx=prx/(1.d0-probh)
2093 c-- following options are not considered for the time being...
2094 c-- (1) good choice, but does not
save very much time:
2095 c-- ff=(1.0d0-sqrt(xf)/2.0d0)/(1.0+sqrt(xf)/2.0d0)
2096 c-- (2) taken from the blue, but works without weight overflows...
2097 c-- ff=(1.d0-xf/(1-(1-sqrt(xf))**2))*(1+(1-sqrt(xf))/sqrt(1-xf))/2
2103 SUBROUTINE phobw(WT)
2104 c.----------------------------------------------------------------------
2106 c. photos: photos boson w correction weight
2108 c. purpose: calculates correction weight due to amplitudes of
2109 c. emission from w boson.
2115 c. input parameters:
Common /phoevt/, with photon added.
2116 c. wt to be corrected
2120 c. output parameters: wt
2122 c. author(s): g. nanava, z. was created at: 13/03/03
2123 c. last update: 13/03/03
2125 c.----------------------------------------------------------------------
2129 parameter(nmxpho=10000)
2130 INTEGER idpho,istpho,jdapho,jmopho,nevpho,npho
2132 common/phoevt/nevpho,npho,istpho(nmxpho),idpho(nmxpho),
2133 &jmopho(2,nmxpho),jdapho(2,nmxpho),ppho(5,nmxpho),vpho(4,nmxpho)
2135 DOUBLE PRECISION emu,mchren,beta,costhg,mpasqr,xph
2137 IF (abs(idpho(1)).EQ.24.AND.
2138 $ abs(idpho(jdapho(1,1) )).GE.11.AND.
2139 $ abs(idpho(jdapho(1,1) )).LE.16.AND.
2140 $ abs(idpho(jdapho(1,1)+1)).GE.11.AND.
2141 $ abs(idpho(jdapho(1,1)+1)).LE.16 )
THEN
2144 $ abs(idpho(jdapho(1,1) )).EQ.11.OR.
2145 $ abs(idpho(jdapho(1,1) )).EQ.13.OR.
2146 $ abs(idpho(jdapho(1,1) )).EQ.15 )
THEN
2152 mchren=abs(ppho(4,i)**2-ppho(3,i)**2
2153 $ -ppho(2,i)**2-ppho(1,i)**2)
2154 beta=sqrt(1- mchren/ ppho(4,i)**2)
2155 costhg=(ppho(3,i)*ppho(3,npho)+ppho(2,i)*ppho(2,npho)
2156 $ +ppho(1,i)*ppho(1,npho))/
2157 $ sqrt(ppho(3,i)**2+ppho(2,i)**2+ppho(1,i)**2) /
2158 $ sqrt(ppho(3,npho)**2+ppho(2,npho)**2+ppho(1,npho)**2)
2161 wt=wt*(1-8*emu*xph*(1-costhg*beta)*
2162 $ (mchren+2*xph*sqrt(mpasqr))/
2163 $ mpasqr**2/(1-mchren/mpasqr)/(4-mchren/mpasqr))
2165 c
write(*,*) idpho(1),idpho(jdapho(1,1)),idpho(jdapho(1,1)+1)
2166 c
write(*,*) emu,xph,costhg,beta,mpasqr,mchren
2169 SUBROUTINE phodo(IP,NCHARB,NEUDAU)
2170 c.----------------------------------------------------------------------
2172 c. photos: photon radiation in decays doing of kinematics
2174 c. purpose: starting from the charged particle energy/momentum,
2175 c. pneutr, photon energy fraction and photon angle with
2176 c. respect to the axis formed by charged particle energy/
2177 c. momentum vector and pneutr, scale the energy/momentum,
2178 c. keeping the original direction of the neutral system in
2179 c. the lab. frame untouched.
2181 c. input parameters: ip:
Pointer to decaying particle in
2182 c. /phoevt/ and the
common itself
2183 c. ncharb:
pointer to the charged radiating
2184 c. daughter in /phoevt/.
2185 c. neudau:
pointer to the first neutral daughter
2186 c. output parameters:
Common /phoevt/, with photon added.
2188 c. author(s): z. was, b. van eijk created at: 26/11/89
2189 c. last update: 27/05/93
2191 c.----------------------------------------------------------------------
2193 DOUBLE PRECISION phoan1,phoan2,angle,fi1,fi3,fi4,fi5,th1,th3,th4
2194 DOUBLE PRECISION parne,qnew,qold,data
2195 INTEGER ip,fi3dum,i,j,neudau,first,last
2197 REAL*8 ephoto,pmavir,photri
2198 REAL*8 gneut,phoran,ccosth,ssinth,pvec(4)
2200 parameter(nmxpho=10000)
2201 INTEGER idpho,istpho,jdapho,jmopho,nevpho,npho
2203 common/phoevt/nevpho,npho,istpho(nmxpho),idpho(nmxpho),
2204 &jmopho(2,nmxpho),jdapho(2,nmxpho),ppho(5,nmxpho),vpho(4,nmxpho)
2205 DOUBLE PRECISION mchsqr,mnesqr
2207 common/phomom/mchsqr,mnesqr,pneutr(5)
2208 DOUBLE PRECISION costhg,sinthg
2209 REAL*8 xphmax,xphoto
2210 common/phophs/xphmax,xphoto,costhg,sinthg
2212 common/phpico/pi,twopi
2214 ephoto=xphoto*ppho(5,ip)/2.d0
2215 pmavir=sqrt(ppho(5,ip)*(ppho(5,ip)-2.d0*ephoto))
2217 c-- reconstruct kinematics of charged particle and neutral system
2218 fi1=phoan1(pneutr(1),pneutr(2))
2220 c-- choose axis along z of pneutr, calculate angle between x and y
2221 c-- components and z and x-y plane and perform lorentz transform...
2222 th1=phoan2(pneutr(3),sqrt(pneutr(1)**2+pneutr(2)**2))
2223 CALL phoro3(-fi1,pneutr(1))
2224 CALL phoro2(-th1,pneutr(1))
2226 c-- take away photon energy from charged particle and pneutr
2227 c-- the onshell charged particle decays into virtual charged particle
2228 c-- and photon. the virtual charged particle mass becomes:
2229 c-- sqrt(ppho(5,ip)*(ppho(5,ip)-2*ephoto)). construct new pneutr mo-
2230 c-- mentum in the rest frame of the parent:
2231 c-- 1) scaling parameters...
2232 qnew=photri(pmavir,pneutr(5),ppho(5,ncharb))
2234 gneut=(qnew**2+qold**2+mnesqr)/(qnew*qold+sqrt((qnew**2+mnesqr)*
2236 IF (gneut.LT.1.d0)
THEN
2238 CALL phoerr(4,
'PHOKIN',data)
2240 parne=gneut-sqrt(max(gneut**2-1.0d0,0.d0))
2242 c-- 2) ...reductive boost...
2243 CALL phobo3(parne,pneutr)
2245 c-- ...calculate photon energy in the reduced system...
2249 c-- photon mother and daughter pointers
2254 ppho(4,npho)=ephoto*ppho(5,ip)/pmavir
2256 c-- ...and photon momenta
2259 th3=phoan2(ccosth,ssinth)
2260 fi3=twopi*phoran(fi3dum)
2261 ppho(1,npho)=ppho(4,npho)*sinthg*cos(fi3)
2262 ppho(2,npho)=ppho(4,npho)*sinthg*sin(fi3)
2264 c-- minus sign because axis opposite direction of charged particle
2265 ppho(3,npho)=-ppho(4,npho)*costhg
2268 c-- rotate in order to get photon along z-axis
2269 CALL phoro3(-fi3,pneutr(1))
2270 CALL phoro3(-fi3,ppho(1,npho))
2271 CALL phoro2(-th3,pneutr(1))
2272 CALL phoro2(-th3,ppho(1,npho))
2273 angle=ephoto/ppho(4,npho)
2275 c-- boost to the rest frame of decaying particle
2276 CALL phobo3(angle,pneutr(1))
2277 CALL phobo3(angle,ppho(1,npho))
2279 c-- back in the parent rest frame but pneutr not yet oriented
2280 fi4=phoan1(pneutr(1),pneutr(2))
2281 th4=phoan2(pneutr(3),sqrt(pneutr(1)**2+pneutr(2)**2))
2282 CALL phoro3(fi4,pneutr(1))
2283 CALL phoro3(fi4,ppho(1,npho))
2288 CALL phoro3(-fi3,pvec)
2289 CALL phoro2(-th3,pvec)
2290 CALL phobo3(angle,pvec)
2291 CALL phoro3(fi4,pvec)
2292 CALL phoro2(-th4,pneutr)
2293 CALL phoro2(-th4,ppho(1,npho))
2294 CALL phoro2(-th4,pvec)
2295 fi5=phoan1(pvec(1),pvec(2))
2297 c-- charged particle restores original direction
2298 CALL phoro3(-fi5,pneutr)
2299 CALL phoro3(-fi5,ppho(1,npho))
2300 CALL phoro2(th1,pneutr(1))
2301 CALL phoro2(th1,ppho(1,npho))
2302 CALL phoro3(fi1,pneutr)
2303 CALL phoro3(fi1,ppho(1,npho))
2304 c-- see whether neutral system has multiplicity larger than 1...
2305 IF ((jdapho(2,ip)-jdapho(1,ip)).GT.1)
THEN
2306 c-- find pointers to components of
'neutral' system
2311 IF (i.NE.ncharb.AND.(jmopho(1,i).EQ.ip))
THEN
2313 c-- reconstruct kinematics...
2314 CALL phoro3(-fi1,ppho(1,i))
2315 CALL phoro2(-th1,ppho(1,i))
2317 c-- ...reductive boost
2318 CALL phobo3(parne,ppho(1,i))
2320 c-- rotate in order to get photon along z-axis
2321 CALL phoro3(-fi3,ppho(1,i))
2322 CALL phoro2(-th3,ppho(1,i))
2324 c-- boost to the rest frame of decaying particle
2325 CALL phobo3(angle,ppho(1,i))
2327 c-- back in the parent rest-frame but pneutr not yet oriented.
2328 CALL phoro3(fi4,ppho(1,i))
2329 CALL phoro2(-th4,ppho(1,i))
2331 c-- charged particle restores original direction
2332 CALL phoro3(-fi5,ppho(1,i))
2333 CALL phoro2(th1,ppho(1,i))
2334 CALL phoro3(fi1,ppho(1,i))
2339 c-- ...only one
'neutral' particle in addition to photon
2341 80 ppho(j,neudau)=pneutr(j)
2344 c-- all
'neutrals' treated, fill /phoevt/ for charged particle...
2346 90 ppho(j,ncharb)=-(ppho(j,npho)+pneutr(j))
2347 ppho(4,ncharb)=ppho(5,ip)-(ppho(4,npho)+pneutr(4))
2350 FUNCTION photri(A,B,C)
2351 c.----------------------------------------------------------------------
2353 c. photos: photon radiation in decays calculation of triangle fie
2355 c. purpose: calculation of triangle
function for phase space.
2357 c. input parameters: a, b, c(virtual) particle masses.
2359 c. output
Parameter:
Function value =
2360 c. sqrt(lambda(a**2,b**2,c**2))/(2*a)
2362 c. author(s): b. van eijk created at: 15/11/89
2363 c. last update: 02/01/90
2365 c.----------------------------------------------------------------------
2367 DOUBLE PRECISION da,db,dc,dapb,damb,dtrian
2374 dtrian=sqrt((damb-dc)*(dapb+dc)*(damb+dc)*(dapb-dc))
2375 photri=dtrian/(da+da)
2378 FUNCTION phoan1(X,Y)
2379 c.----------------------------------------------------------------------
2381 c. photos: photon radiation in decays calculation of angle
'1'
2383 c. purpose: calculate angle from x and y
2385 c. input parameters: x, y
2387 c. output
Parameter:
Function value
2389 c. author(s): s. jadach created at: 01/01/89
2390 c. b. van eijk last update: 02/01/90
2392 c.----------------------------------------------------------------------
2394 DOUBLE PRECISION phoan1
2397 common/phpico/pi,twopi
2398 IF (abs(y).LT.abs(x))
THEN
2399 phoan1=atan(abs(y/x))
2400 IF (x.LE.0.d0) phoan1=pi-phoan1
2402 phoan1=acos(x/sqrt(x**2+y**2))
2404 IF (y.LT.0.d0) phoan1=twopi-phoan1
2407 FUNCTION phoan2(X,Y)
2408 c.----------------------------------------------------------------------
2410 c. photos: photon radiation in decays calculation of angle
'2'
2412 c. purpose: calculate angle from x and y
2414 c. input parameters: x, y
2416 c. output
Parameter:
Function value
2418 c. author(s): s. jadach created at: 01/01/89
2419 c. b. van eijk last update: 02/01/90
2421 c.----------------------------------------------------------------------
2423 DOUBLE PRECISION phoan2
2426 common/phpico/pi,twopi
2427 IF (abs(y).LT.abs(x))
THEN
2428 phoan2=atan(abs(y/x))
2429 IF (x.LE.0.d0) phoan2=pi-phoan2
2431 phoan2=acos(x/sqrt(x**2+y**2))
2435 SUBROUTINE phobo3(ANGLE,PVEC)
2436 c.----------------------------------------------------------------------
2438 c. photos: photon radiation in decays boost routine
'3'
2440 c. purpose: boost vector pvec along z-axis
where angle = exp(eta),
2441 c. eta is the hyperbolic velocity.
2443 c. input parameters: angle, pvec
2445 c. output
Parameter: pvec
2447 c. author(s): s. jadach created at: 01/01/89
2448 c. b. van eijk last update: 02/01/90
2450 c.----------------------------------------------------------------------
2452 DOUBLE PRECISION qpl,qmi,angle
2454 qpl=(pvec(4)+pvec(3))*angle
2455 qmi=(pvec(4)-pvec(3))/angle
2456 pvec(3)=(qpl-qmi)/2.d0
2457 pvec(4)=(qpl+qmi)/2.d0
2460 SUBROUTINE phoro2(ANGLE,PVEC)
2461 c.----------------------------------------------------------------------
2463 c. photos: photon radiation in decays rotation routine
'2'
2465 c. purpose: rotate x and z components of vector pvec around angle
2468 c. input parameters: angle, pvec
2470 c. output
Parameter: pvec
2472 c. author(s): s. jadach created at: 01/01/89
2473 c. b. van eijk last update: 02/01/90
2475 c.----------------------------------------------------------------------
2477 DOUBLE PRECISION cs,sn,angle
2479 cs=cos(angle)*pvec(1)+sin(angle)*pvec(3)
2480 sn=-sin(angle)*pvec(1)+cos(angle)*pvec(3)
2485 SUBROUTINE phoro3(ANGLE,PVEC)
2486 c.----------------------------------------------------------------------
2488 c. photos: photon radiation in decays rotation routine
'3'
2490 c. purpose: rotate x and y components of vector pvec around angle
2493 c. input parameters: angle, pvec
2495 c. output
Parameter: pvec
2497 c. author(s): s. jadach created at: 01/01/89
2498 c. b. van eijk last update: 02/01/90
2500 c.----------------------------------------------------------------------
2502 DOUBLE PRECISION cs,sn,angle
2504 cs=cos(angle)*pvec(1)-sin(angle)*pvec(2)
2505 sn=sin(angle)*pvec(1)+cos(angle)*pvec(2)
2511 c.----------------------------------------------------------------------
2513 c. photos: photon radiation in decays random number generator init
2515 c. purpose: initialse phoran with the user specified seeds in the
2516 c. array iseed. for details see also: f. james cern dd-
2517 c. report november 1988.
2519 c. input parameters: iseed(*)
2521 c. output parameters: uran, cran, cdran, cmran, i97, j97
2523 c. author(s): b. van eijk and f. james created at: 27/09/89
2524 c. last update: 22/02/90
2526 c.----------------------------------------------------------------------
2528 DOUBLE PRECISION data
2530 INTEGER i,is1,is2,is3,is4,is5,j
2531 INTEGER iseed,i97,j97
2532 REAL*8 uran,cran,cdran,cmran
2533 common/phseed/iseed(2),i97,j97,uran(97),cran,cdran,cmran
2535 c-- check value range of seeds
2536 IF ((iseed(1).LT.0).OR.(iseed(1).GE.31328))
THEN
2538 CALL phoerr(8,
'PHORIN',data)
2540 IF ((iseed(2).LT.0).OR.(iseed(2).GE.30081))
THEN
2542 CALL phoerr(9,
'PHORIN',data)
2545 c-- calculate marsaglia and zaman seeds(by f. james)
2546 is1=mod(iseed(1)/177,177)+2
2547 is2=mod(iseed(1),177)+2
2548 is3=mod(iseed(2)/169,178)+1
2549 is4=mod(iseed(2),169)
2554 is5=mod(mod(is1*is2,179)*is3,179)
2558 is4=mod(53*is4+1,169)
2559 IF (mod(is4*is5,64).GE.32) s=s+t
2562 cran=362436.d0/16777216.d0
2563 cdran=7654321.d0/16777216.d0
2564 cmran=16777213.d0/16777216.d0
2569 FUNCTION phoran(IDUM)
2570 c.----------------------------------------------------------------------
2572 c. photos: photon radiation in decays random number generator based
2573 c. on marsaglia algorithm
2575 c. purpose: generate uniformly distributed random numbers between
2576 c. 0 and 1. super long period: 2**144. see also:
2577 c. g. marsaglia and a. zaman, fsu-scr-87-50, for seed mo-
2578 c. difications to this version see: f. james dd-report,
2579 c. november 1988. the generator has to be initialized by
2580 c. a call to phorin.
2582 c. input parameters: idum(
integer dummy)
2584 c. output parameters:
Function value
2586 c. author(s): b. van eijk, g. marsaglia and created at: 27/09/89
2587 c. a. zaman last update: 27/09/89
2589 c.----------------------------------------------------------------------
2593 INTEGER iseed,i97,j97
2594 REAL*8 uran,cran,cdran,cmran
2595 common/phseed/iseed(2),i97,j97,uran(97),cran,cdran,cmran
2596 10 phoran=uran(i97)-uran(j97)
2597 IF (phoran.LT.0.d0) phoran=phoran+1.d0
2600 IF (i97.EQ.0) i97=97
2602 IF (j97.EQ.0) j97=97
2604 IF (cran.LT.0.d0) cran=cran+cmran
2606 IF (phoran.LT.0.d0) phoran=phoran+1.d0
2607 IF (phoran.LE.0.d0) goto 10
2610 FUNCTION phocha(IDHEP)
2611 c.----------------------------------------------------------------------
2613 c. photos: photon radiation in decays charge determination
2615 c. purpose: calculate the charge of particle with code idhep. the
2616 c. code of the particle is defined by the particle
Data
2617 c. group in phys. lett. b204(1988) 1.
2619 c. input
Parameter: idhep
2621 c. output
Parameter: funtion value = charge of particle with code
2624 c. author(s): e. barberio and b. van eijk created at: 29/11/89
2625 c. last update: 02/01/90
2627 c.----------------------------------------------------------------------
2630 INTEGER idhep,idabs,q1,q2,q3
2632 c-- array
'CHARGE' contains the charge of the first 101 particles ac-
2633 c-- cording to the pdg particle code... (0 is added for convenience)
2634 REAL*8 charge(0:100)
2636 &-0.3333333333d0, 0.6666666667d0, -0.3333333333d0, 0.6666666667d0,
2637 &-0.3333333333d0, 0.6666666667d0, -0.3333333333d0, 0.6666666667d0,
2638 & 2*0.d0, -1.d0, 0.d0, -1.d0, 0.d0, -1.d0, 0.d0, -1.d0, 6*0.d0,
2639 & 1.d0, 12*0.d0, 1.d0, 63*0.d0/
2641 IF (idabs.LE.100)
THEN
2643 c-- charge of quark, lepton, boson etc....
2644 phocha = charge(idabs)
2647 c-- check on particle build out of quarks, unpack its code...
2648 q3=mod(idabs/1000,10)
2649 q2=mod(idabs/100,10)
2654 IF(mod(q2,2).EQ.0)
THEN
2655 phocha=charge(q2)-charge(q1)
2657 phocha=charge(q1)-charge(q2)
2661 c-- ...diquarks or baryon.
2662 phocha=charge(q1)+charge(q2)+charge(q3)
2666 c-- find the sign of the charge...
2667 IF (idhep.LT.0.d0) phocha=-phocha
2668 IF (phocha**2.lt.1d-6) phocha=0.d0
2671 FUNCTION phospi(IDHEP)
2672 c.----------------------------------------------------------------------
2674 c. photos: photon radiation in decays
function for spin determina-
2677 c. purpose: calculate the spin of particle with code idhep. the
2678 c. code of the particle is defined by the particle
Data
2679 c. group in phys. lett. b204(1988) 1.
2681 c. input
Parameter: idhep
2683 c. output
Parameter: funtion value = spin of particle with code
2686 c. author(s): e. barberio and b. van eijk created at: 29/11/89
2687 c. last update: 02/01/90
2689 c.----------------------------------------------------------------------
2694 c-- array
'SPIN' contains the spin of the first 100 particles accor-
2695 c-- ding to the pdg particle code...
2697 DATA spin/ 8*.5d0, 1.d0, 0.d0, 8*.5d0, 2*0.d0, 4*1.d0, 76*0.d0/
2700 c-- spin of quark, lepton, boson etc....
2701 IF (idabs.LE.100)
THEN
2705 c-- ...other particles, however...
2706 phospi=(mod(idabs,10)-1.d0)/2.d0
2708 c-- ...k_short and k_long are special
2709 phospi=max(phospi,0.d0)
2713 SUBROUTINE phoerr(IMES,TEXT,DATA)
2714 c.----------------------------------------------------------------------
2716 c. photos: photon radiation in decays errror handling
2718 c. purpose: inform user about(fatal) errors and warnings generated
2719 c. by either the user or the program.
2721 c. input parameters: imes, text,
DATA
2723 c. output parameters: none
2725 c. author(s): b. van eijk created at: 29/11/89
2726 c. last update: 10/01/92
2728 c.----------------------------------------------------------------------
2730 DOUBLE PRECISION data
2736 parameter(phomes=10)
2738 common/phosta/status(phomes)
2741 c-- security stop switch
2746 IF (imes.LE.phomes) status(imes)=status(imes)+1
2748 c-- count number of non-fatal errors...
2749 IF ((imes.EQ. 6).AND.(status(imes).GE.2))
RETURN
2750 IF ((imes.EQ.10).AND.(status(imes).GE.2))
RETURN
2754 goto(10,20,30,40,50,60,70,80,90,100),imes
2755 WRITE(phlun,9130) imes
2757 10
WRITE(phlun,9010) text,int(sdata)
2759 20
WRITE(phlun,9020) text,sdata
2761 30
WRITE(phlun,9030) text,sdata
2763 40
WRITE(phlun,9040) text
2765 50
WRITE(phlun,9050) text,int(sdata)
2767 60
WRITE(phlun,9060) text,sdata
2769 70
WRITE(phlun,9070) text,int(sdata)
2771 80
WRITE(phlun,9080) text,int(sdata)
2773 90
WRITE(phlun,9090) text,int(sdata)
2775 100
WRITE(phlun,9100) text,sdata
2787 IF (ierror.GE.10)
THEN
2797 130
WRITE(phlun,9120)
2800 9000
FORMAT(1h ,80(
'*'))
2801 9010
FORMAT(1h ,
'* ',a,
': Too many charged Particles, NCHARG =',i6,t81,
2803 9020
FORMAT(1h ,
'* ',a,
': Too much Bremsstrahlung required, PRSOFT = ',
2805 9030
FORMAT(1h ,
'* ',a,
': Combined Weight is exceeding 1., Weight = ',
2807 9040
FORMAT(1h ,
'* ',a,
2808 &
': Error in Rescaling charged and neutral Vectors',t81,
'*')
2809 9050
FORMAT(1h ,
'* ',a,
2810 &
': Non matching charged Particle Pointer, NCHARG = ',i5,t81,
'*')
2811 9060
FORMAT(1h ,
'* ',a,
2812 &
': Do you really work with a Particle of Spin: ',f4.1,
' ?',t81,
2814 9070
FORMAT(1h ,
'* ',a,
': Stack Length exceeded, NSTACK = ',i5 ,t81,
2816 9080
FORMAT(1h ,
'* ',a,
2817 &
': Random Number Generator Seed(1) out of Range: ',i8,t81,
'*')
2818 9090
FORMAT(1h ,
'* ',a,
2819 &
': Random Number Generator Seed(2) out of Range: ',i8,t81,
'*')
2820 9100
FORMAT(1h ,
'* ',a,
2821 &
': Available Phase Space below Cut-off: ',f15.6,
' GeV/c^2',t81,
2823 9120
FORMAT(1h ,
'*',t81,
'*')
2824 9130
FORMAT(1h ,
'* Funny Error Message: ',i4,
' ! What to do ?',t81,
'*')
2825 9140
FORMAT(1h ,
'* Fatal Error Message, I stop this Run !',t81,
'*')
2826 9150
FORMAT(1h ,
'* 10 Error Messages generated, I stop this Run !',t81,
2830 c.----------------------------------------------------------------------
2832 c. photos: photon radiation in decays run summary report
2834 c. purpose: inform user about success and/or restrictions of photos
2835 c. encountered during execution.
2837 c. input parameters:
Common /phosta/
2839 c. output parameters: none
2841 c. author(s): b. van eijk created at: 10/01/92
2842 c. last update: 10/01/92
2844 c.----------------------------------------------------------------------
2849 parameter(phomes=10)
2851 common/phosta/status(phomes)
2863 IF (status(i).EQ.0) goto 10
2864 IF ((i.EQ.6).OR.(i.EQ.10))
THEN
2865 WRITE(phlun,9050) i,status(i)
2868 WRITE(phlun,9060) i,status(i)
2871 IF (.NOT.error)
WRITE(phlun,9070)
2876 9010
FORMAT(1h ,80(
'*'))
2877 9020
FORMAT(1h ,
'*',t81,
'*')
2878 9030
FORMAT(1h ,
'*',26x,25(
'='),t81,
'*')
2879 9040
FORMAT(1h ,
'*',30x,
'PHOTOS Run Summary',t81,
'*')
2880 9050
FORMAT(1h ,
'*',22x,
'Warning #',i2,
' occured',i6,
' times',t81,
'*')
2881 9060
FORMAT(1h ,
'*',23x,
'Error #',i2,
' occured',i6,
' times',t81,
'*')
2882 9070
FORMAT(1h ,
'*',16x,
'PHOTOS Execution has successfully terminated',
2885 SUBROUTINE phlupa(IPOINT)
2887 c.----------------------------------------------------------------------
2889 c. phlupa: debugging tool
2891 c. purpose: none, eventually may printout content of the
2894 c. input parameters:
Common /phoevt/ and /phnum/
2895 c. latter may have number of the event.
2897 c. output parameters: none
2899 c. author(s): z. was created at: 30/05/93
2900 c. last update: 09/10/05
2902 c.----------------------------------------------------------------------
2904 parameter(nmxpho=10000)
2905 INTEGER idpho,istpho,jdapho,jmopho,nevpho,npho,i,j,ipoint
2906 INTEGER ipoin,ipoin0,ipoinm,iev
2908 REAL*8 ppho,vpho,sum
2909 common/phoevt/nevpho,npho,istpho(nmxpho),idpho(nmxpho),
2910 &jmopho(2,nmxpho),jdapho(2,nmxpho),ppho(5,nmxpho),vpho(4,nmxpho)
2916 COMMON /phlupy/ ipoin,ipoinm
2918 IF (ipoin0.LT.0)
THEN
2923 IF (ipoint.LE.ipoinm.OR.ipoint.GE.ipoin )
RETURN
2925 IF (iev.LT.1000)
THEN
2929 WRITE(phlun,*)
'EVENT NR=',iev,
2930 $
'WE ARE TESTING /PHOEVT/ at IPOINT=',ipoint
2933 WRITE(phlun,20) idpho(i),ppho(1,i),ppho(2,i),ppho(3,i),
2934 $ ppho(4,i),ppho(5,i),jdapho(1,i),jdapho(2,i)
2936 WRITE(phlun,20) idpho(i),ppho(1,i),ppho(2,i),ppho(3,i),
2937 $ ppho(4,i),ppho(5,i),jdapho(1,i),jdapho(2,i)
2940 WRITE(phlun,20) idpho(i),ppho(1,i),ppho(2,i),ppho(3,i),
2941 $ ppho(4,i),ppho(5,i),jmopho(1,i),jmopho(2,i)
2943 sum(j)=sum(j)+ppho(j,i)
2946 sum(5)=sqrt(abs(sum(4)**2-sum(1)**2-sum(2)**2-sum(3)**2))
2948 10
FORMAT(1x,
' ID ',
'p_x ',
'p_y ',
'p_z ',
2950 $
'ID-MO_DA1',
'ID-MO DA2' )
2951 20
FORMAT(1x,i4,5(f14.9),2i9)
2952 30
FORMAT(1x,
' SUM',5(f14.9))
2958 FUNCTION iphqrk(MODCOR)
2960 c.----------------------------------------------------------------------
2962 c. iphqrk: enables blocks emision from quarks
2965 c. input parameters: modcor
2966 c. modcor >0
type of action
2969 c. =0 execution mode(retrns stored value)
2972 c. author(s): z. was created at: 11/12/00
2974 c.----------------------------------------------------------------------
2975 INTEGER iphqrk,modcor,modop
2981 IF (modcor.NE.0)
THEN
2986 $
'Message from PHOTOS: IPHQRK(MODCOR):: (re)initialization'
2987 IF (modop.EQ.1)
THEN
2989 $
'MODOP=1 -- blocks emission from light quarks: DEFAULT'
2990 ELSEIF (modop.EQ.2)
THEN
2992 $
'MODOP=2 -- enables emission from light quarks: TEST '
2994 WRITE(phlun,*)
'IPHQRK wrong MODCOR=',modcor
3000 IF (modop.EQ.0.AND.modcor.EQ.0)
THEN
3001 WRITE(phlun,*)
'IPHQRK lack of initialization'
3008 FUNCTION iphekl(MODCOR)
3010 c.----------------------------------------------------------------------
3012 c. iphekl: enables/blocks emision in: pi0 to gamma e+ e-
3015 c. input parameters: modcor
3016 c. modcor >0
type of action
3019 c. =0 execution mode(retrns stored value)
3022 c. author(s): z. was created at: 11/12/00
3024 c.----------------------------------------------------------------------
3025 INTEGER iphekl,modcor,modop
3032 IF (modcor.NE.0)
THEN
3037 $
'Message from PHOTOS: IPHEKL(MODCOR):: (re)initialization'
3038 IF (modop.EQ.2)
THEN
3040 $
'MODOP=2 -- blocks emission in pi0 to gamma e+e-: DEFAULT'
3042 $
'MODOP=2 -- blocks emission in Kl to gamma e+e-: DEFAULT'
3043 ELSEIF (modop.EQ.1)
THEN
3045 $
'MODOP=1 -- enables emission in pi0 to gamma e+e- : TEST '
3047 $
'MODOP=1 -- enables emission in Kl to gamma e+e- : TEST '
3049 WRITE(phlun,*)
'IPHEKL wrong MODCOR=',modcor
3055 IF (modop.EQ.0.AND.modcor.EQ.0)
THEN
3056 WRITE(phlun,*)
'IPHELK lack of initialization'
3062 SUBROUTINE phcork(MODCOR)
3064 c.----------------------------------------------------------------------
3066 c. phcork: corrects kinmatics of subbranch needed
if host
program
3067 c. produces events with the shaky momentum conservation
3069 c. input parameters:
Common /phoevt/, modcor
3070 c. modcor >0
type of action
3072 c. =2 corrects energy from mass
3073 c. =3 corrects mass from energy
3074 c. =4 corrects energy from mass for
3075 c. particles up to .4 gev mass,
3076 c. for heavier ones corrects mass,
3077 c. =5 most complete correct also of mother
3078 c. often necessary for exponentiation.
3079 c. =0 execution mode
3081 c. output parameters: corrected /phoevt/
3083 c. author(s): p.golonka, z. was created at: 01/02/99
3084 c. modified : 08/02/99
3085 c.----------------------------------------------------------------------
3087 parameter(nmxpho=10000)
3089 REAL*8 m,p2,px,py,pz,e,en,mcut,xms
3090 INTEGER modcor,modop,i,iev,iprint,k
3091 INTEGER idpho,istpho,jdapho,jmopho,nevpho,npho
3093 common/phoevt/nevpho,npho,istpho(nmxpho),idpho(nmxpho),
3094 &jmopho(2,nmxpho),jdapho(2,nmxpho),ppho(5,nmxpho),vpho(4,nmxpho)
3105 IF (modcor.NE.0)
THEN
3109 WRITE(phlun,*)
'Message from PHCORK(MODCOR):: initialization'
3110 IF (modop.EQ.1)
THEN
3111 WRITE(phlun,*)
'MODOP=1 -- no corrections on event: DEFAULT'
3112 ELSEIF (modop.EQ.2)
THEN
3113 WRITE(phlun,*)
'MODOP=2 -- corrects Energy from mass'
3114 ELSEIF (modop.EQ.3)
THEN
3115 WRITE(phlun,*)
'MODOP=3 -- corrects mass from Energy'
3116 ELSEIF (modop.EQ.4)
THEN
3117 WRITE(phlun,*)
'MODOP=4 -- corrects Energy from mass to Mcut'
3118 WRITE(phlun,*)
' and mass from energy above Mcut '
3120 WRITE(phlun,*)
'Mcut=',mcut,
'GeV'
3121 ELSEIF (modop.EQ.5)
THEN
3122 WRITE(phlun,*)
'MODOP=5 -- corrects Energy from mass+flow'
3125 WRITE(phlun,*)
'PHCORK wrong MODCOR=',modcor
3131 IF (modop.EQ.0.AND.modcor.EQ.0)
THEN
3132 WRITE(phlun,*)
'PHCORK lack of initialization'
3146 IF (modop.EQ.1)
THEN
3147 c -----------------------
3148 c in this
case we
do nothing
3150 ELSEIF(modop.EQ.2)
THEN
3151 c -----------------------
3152 cc lets loop thru all daughters and correct their energies
3153 cc according to e^2=p^2+m^2
3161 p2=ppho(1,i)**2+ppho(2,i)**2+ppho(3,i)**2
3163 en=sqrt( ppho(5,i)**2 + p2)
3166 &
WRITE(phlun,*)
"CORRECTING ENERGY OF ",i,
":",
3174 ELSEIF(modop.EQ.5)
THEN
3175 c -----------------------
3176 cc lets loop thru all daughters and correct their energies
3177 cc according to e^2=p^2+m^2
3185 p2=ppho(1,i)**2+ppho(2,i)**2+ppho(3,i)**2
3187 en=sqrt( ppho(5,i)**2 + p2)
3190 &
WRITE(phlun,*)
"CORRECTING ENERGY OF ",i,
":",
3200 ppho(k,1)=ppho(k,1)+ppho(k,i)
3203 xms=sqrt(ppho(4,1)**2-ppho(3,1)**2-ppho(2,1)**2-ppho(1,1)**2)
3205 ELSEIF(modop.EQ.3)
THEN
3206 c -----------------------
3208 cc lets loop thru all daughters and correct their masses
3209 cc according to e^2=p^2+m^2
3218 p2=ppho(1,i)**2+ppho(2,i)**2+ppho(3,i)**2
3220 m=sqrt(abs( ppho(4,i)**2 - p2))
3223 &
WRITE(phlun,*)
"CORRECTING MASS OF ",i,
":",
3231 ELSEIF(modop.EQ.4)
THEN
3232 c -----------------------
3234 cc lets loop thru all daughters and correct their masses
3235 cc or energies according to e^2=p^2+m^2
3243 p2=ppho(1,i)**2+ppho(2,i)**2+ppho(3,i)**2
3245 m=sqrt(abs( ppho(4,i)**2 - p2))
3249 &
WRITE(phlun,*)
"CORRECTING MASS OF ",i,
":",
3255 en=sqrt( ppho(5,i)**2 + p2)
3258 &
WRITE(phlun,*)
"CORRECTING ENERGY OF ",i,
":",
3269 IF (iprint.EQ.1)
THEN
3270 WRITE(phlun,*)
"CORRECTING MOTHER"
3271 WRITE(phlun,*)
"PX:",ppho(1,1),
"=>",px-ppho(1,2)
3272 WRITE(phlun,*)
"PY:",ppho(2,1),
"=>",py-ppho(2,2)
3273 WRITE(phlun,*)
"PZ:",ppho(3,1),
"=>",pz-ppho(3,2)
3274 WRITE(phlun,*)
" E:",ppho(4,1),
"=>",e-ppho(4,2)
3277 ppho(1,1)=px-ppho(1,2)
3278 ppho(2,1)=py-ppho(2,2)
3279 ppho(3,1)=pz-ppho(3,2)
3280 ppho(4,1)=e -ppho(4,2)
3282 p2=ppho(1,1)**2+ppho(2,1)**2+ppho(3,1)**2
3284 IF (ppho(4,1)**2.GT.p2)
THEN
3285 m=sqrt( ppho(4,1)**2 - p2 )
3287 &
WRITE(phlun,*)
" M:",ppho(5,1),
"=>",m
3297 FUNCTION phint(IDUM)
3298 c --- can be used with variant a. for b
use phint1 or 2 --------------
3299 c.----------------------------------------------------------------------
3301 c. phint: photos universal interference correction weight
3303 c. purpose: calculates correction weight as expressed by
3304 c formula(17) from cpc 79 (1994), 291.
3306 c. input parameters:
Common /phoevt/, with photon added.
3308 c. output parameters: correction weight
3310 c. author(s): z. was, p.golonka created at: 19/01/05
3311 c. last update: 25/01/05
3313 c.----------------------------------------------------------------------
3318 parameter(nmxpho=10000)
3319 INTEGER idpho,istpho,jdapho,jmopho,nevpho,npho
3321 common/phoevt/nevpho,npho,istpho(nmxpho),idpho(nmxpho),
3322 &jmopho(2,nmxpho),jdapho(2,nmxpho),ppho(5,nmxpho),vpho(4,nmxpho)
3324 DOUBLE PRECISION emu,mchren,beta,costhg,mpasqr,xph, xc1, xc2,xdeno
3325 DOUBLE PRECISION xnum1,xnum2
3326 DOUBLE PRECISION eps1(4),eps2(4),ph(4),pl(4)
3330 c calculate polarimetric vector: ph, eps1, eps2 are orthogonal
3337 CALL phoeps(ph,eps2,eps1)
3338 CALL phoeps(ph,eps1,eps2)
3345 DO k=jdapho(1,1),npho-1
3347 c momenta of charged particle in pl
3351 c scalar products: epsilon*p/k*p
3353 xc1 = - phocha(idpho(k)) *
3354 & ( pl(1)*eps1(1) + pl(2)*eps1(2) + pl(3)*eps1(3) ) /
3355 & ( ph(4)*pl(4) - ph(1)*pl(1) - ph(2)*pl(2) - ph(3)*pl(3) )
3357 xc2 = - phocha(idpho(k)) *
3358 & ( pl(1)*eps2(1) + pl(2)*eps2(2) + pl(3)*eps2(3) ) /
3359 & ( ph(4)*pl(4) - ph(1)*pl(1) - ph(2)*pl(2) - ph(3)*pl(3) )
3362 c accumulate the currents
3366 xdeno = xdeno + xc1**2 + xc2**2
3370 phint=(xnum1**2 + xnum2**2) / xdeno
3376 SUBROUTINE phoeps (VEC1, VEC2, EPS)
3377 c.----------------------------------------------------------------------
3379 c. phoeps: phoeps vector product(normalized to unity)
3381 c. purpose: calculates vector product,
then normalizes its length.
3382 c used to generate orthogonal vectors, i.e. to
3383 c generate polarimetric vectors for photons.
3385 c. input parameters: vec1,vec2 - input 4-vectors
3387 c. output parameters: eps - normalized 4-vector, orthogonal to
3390 c. author(s): z. was, p.golonka created at: 19/01/05
3391 c. last update: 25/01/05
3393 c.----------------------------------------------------------------------
3395 DOUBLE PRECISION vec1(4), vec2(4), eps(4),xn
3397 eps(1)=vec1(2)*vec2(3) - vec1(3)*vec2(2)
3398 eps(2)=vec1(3)*vec2(1) - vec1(1)*vec2(3)
3399 eps(3)=vec1(1)*vec2(2) - vec1(2)*vec2(1)
3402 xn=sqrt( eps(1)**2 +eps(2)**2 +eps(3)**2)
3411 c.----------------------------------------------------------------------
3413 c. photos: photon radiation in decays event dump routine
3415 c. purpose: print event record.
3417 c. input parameters:
Common /hepevt/
3419 c. output parameters: none
3421 c. author(s): b. van eijk created at: 05/06/90
3422 c. last update: 05/06/90
3424 c.----------------------------------------------------------------------
3426 DOUBLE PRECISION sumvec(5)
3428 c this is the hepevt class in old style. no d_h_ class pre-name
3430 parameter(nmxhep=10000)
3432 INTEGER nevhep,nhep,isthep,idhep,jmohep,
3443 * ----------------------------------------------------------------------
3447 * ----------------------------------------------------------------------
3454 c-- print event number...
3456 WRITE(phlun,9010) nevhep
3461 c-- for
'stable particle' calculate vector momentum sum
3462 IF (jdahep(1,i).EQ.0)
THEN
3464 20 sumvec(j)=sumvec(j)+phep(j,i)
3465 IF (jmohep(2,i).EQ.0)
THEN
3466 WRITE(phlun,9030) i,idhep(i),jmohep(1,i),(phep(j,i),j=1,5)
3468 WRITE(phlun,9040) i,idhep(i),jmohep(1,i),jmohep(2,i),(phep
3472 IF (jmohep(2,i).EQ.0)
THEN
3473 WRITE(phlun,9050) i,idhep(i),jmohep(1,i),jdahep(1,i),
3474 & jdahep(2,i),(phep(j,i),j=1,5)
3476 WRITE(phlun,9060) i,idhep(i),jmohep(1,i),jmohep(2,i),
3477 & jdahep(1,i),jdahep(2,i),(phep(j,i),j=1,5)
3481 sumvec(5)=sqrt(sumvec(4)**2-sumvec(1)**2-sumvec(2)**2-
3483 WRITE(phlun,9070) (sumvec(j),j=1,5)
3485 9000
FORMAT(1h0,80(
'='))
3486 9010
FORMAT(1h ,29x,
'Event No.:',i10)
3487 9020
FORMAT(1h0,1x,
'Nr',3x,
'Type',3x,
'Parent(s)',2x,
'Daughter(s)',6x,
3488 &
'Px',7x,
'Py',7x,
'Pz',7x,
'E',4x,
'Inv. M.')
3489 9030
FORMAT(1h ,i4,i7,3x,i4,9x,
'Stable',2x,5f9.2)
3490 9040
FORMAT(1h ,i4,i7,i4,
' - ',i4,5x,
'Stable',2x,5f9.2)
3491 9050
FORMAT(1h ,i4,i7,3x,i4,6x,i4,
' - ',i4,5f9.2)
3492 9060
FORMAT(1h ,i4,i7,i4,
' - ',i4,2x,i4,
' - ',i4,5f9.2)
3493 9070
FORMAT(1h0,23x,
'Vector Sum: ', 5f9.2)
3494 9080
FORMAT(1h0,6x,
'Particle Parameters')