1 SUBROUTINE PHODO(IP,NCHARB,NEUDAU)
2 C.----------------------------------------------------------------------
4 C. PHOTOS: PHOton radiation in decays DOing of KINematics
6 C. Purpose: Starting from the charged particle energy/momentum,
7 C. PNEUTR, photon energy fraction and photon angle with
8 C. respect to the axis formed by charged particle energy/
9 C. momentum vector and PNEUTR, scale the energy/momentum,
10 C. keeping the original direction of the neutral system in
11 C. the lab. frame untouched.
13 C. Input Parameters: IP: Pointer to decaying particle in
14 C. /PHOEVT/ and the common itself
15 C. NCHARB: pointer to the charged radiating
16 C. daughter in /PHOEVT/.
17 C. NEUDAU: pointer to the first neutral daughter
18 C. Output Parameters: Common /PHOEVT/, with photon added.
20 C. Author(s): Z. Was, B. van Eijk Created at: 26/11/89
21 C. Last Update: 27/05/93
23 C.----------------------------------------------------------------------
25 DOUBLE PRECISION PHOAN1,PHOAN2,ANGLE,FI1,FI3,FI4,FI5,TH1,TH3,TH4
26 DOUBLE PRECISION PARNE,QNEW,QOLD,DATA
27 INTEGER IP,FI3DUM,I,J,NEUDAU,FIRST,LAST
29 REAL*8 EPHOTO,PMAVIR,PHOTRI
30 REAL*8 GNEUT,PHORAN,CCOSTH,SSINTH,PVEC(4)
32 PARAMETER (NMXPHO=10000)
33 INTEGER IDPHO,ISTPHO,JDAPHO,JMOPHO,NEVPHO,NPHO
35 COMMON/PHOEVT/NEVPHO,NPHO,ISTPHO(NMXPHO),IDPHO(NMXPHO),
36 &JMOPHO(2,NMXPHO),JDAPHO(2,NMXPHO),PPHO(5,NMXPHO),VPHO(4,NMXPHO)
37 DOUBLE PRECISION MCHSQR,MNESQR
39 COMMON/PHOMOM/MCHSQR,MNESQR,PNEUTR(5)
40 DOUBLE PRECISION COSTHG,SINTHG
42 COMMON/PHOPHS/XPHMAX,XPHOTO,COSTHG,SINTHG
44 COMMON/PHPICO/PI,TWOPI
46 EPHOTO=XPHOTO*PPHO(5,IP)/2.D0
47 PMAVIR=SQRT(PPHO(5,IP)*(PPHO(5,IP)-2.D0*EPHOTO))
49 C-- Reconstruct kinematics of charged particle and neutral system
50 FI1=PHOAN1(PNEUTR(1),PNEUTR(2))
52 C-- Choose axis along z of PNEUTR, calculate angle between x and y
53 C-- components and z and x-y plane and perform Lorentz transform...
54 TH1=PHOAN2(PNEUTR(3),SQRT(PNEUTR(1)**2+PNEUTR(2)**2))
55 CALL PHORO3(-FI1,PNEUTR(1))
56 CALL PHORO2(-TH1,PNEUTR(1))
58 C-- Take away photon energy from charged particle and PNEUTR ! Thus
59 C-- the onshell charged particle decays into virtual charged particle
60 C-- and photon. The virtual charged particle mass becomes:
61 C-- SQRT(PPHO(5,IP)*(PPHO(5,IP)-2*EPHOTO)). Construct new PNEUTR mo-
62 C-- mentum in the rest frame of the parent:
63 C-- 1) Scaling parameters...
64 QNEW=PHOTRI(PMAVIR,PNEUTR(5),PPHO(5,NCHARB))
66 GNEUT=(QNEW**2+QOLD**2+MNESQR)/(QNEW*QOLD+SQRT((QNEW**2+MNESQR)*
68 IF (GNEUT.LT.1.D0) THEN
70 CALL PHOERR(4,'PHOKIN',DATA)
72 PARNE=GNEUT-SQRT(MAX(GNEUT**2-1.0D0,0.D0))
74 C-- 2) ...reductive boost...
75 CALL PHOBO3(PARNE,PNEUTR)
77 C-- ...calculate photon energy in the reduced system...
81 C-- Photon mother and daughter pointers !
86 PPHO(4,NPHO)=EPHOTO*PPHO(5,IP)/PMAVIR
88 C-- ...and photon momenta
91 TH3=PHOAN2(CCOSTH,SSINTH)
92 FI3=TWOPI*PHORAN(FI3DUM)
93 PPHO(1,NPHO)=PPHO(4,NPHO)*SINTHG*COS(FI3)
94 PPHO(2,NPHO)=PPHO(4,NPHO)*SINTHG*SIN(FI3)
96 C-- Minus sign because axis opposite direction of charged particle !
97 PPHO(3,NPHO)=-PPHO(4,NPHO)*COSTHG
100 C-- Rotate in order to get photon along z-axis
101 CALL PHORO3(-FI3,PNEUTR(1))
102 CALL PHORO3(-FI3,PPHO(1,NPHO))
103 CALL PHORO2(-TH3,PNEUTR(1))
104 CALL PHORO2(-TH3,PPHO(1,NPHO))
105 ANGLE=EPHOTO/PPHO(4,NPHO)
107 C-- Boost to the rest frame of decaying particle
108 CALL PHOBO3(ANGLE,PNEUTR(1))
109 CALL PHOBO3(ANGLE,PPHO(1,NPHO))
111 C-- Back in the parent rest frame but PNEUTR not yet oriented !
112 FI4=PHOAN1(PNEUTR(1),PNEUTR(2))
113 TH4=PHOAN2(PNEUTR(3),SQRT(PNEUTR(1)**2+PNEUTR(2)**2))
114 CALL PHORO3(FI4,PNEUTR(1))
115 CALL PHORO3(FI4,PPHO(1,NPHO))
120 CALL PHORO3(-FI3,PVEC)
121 CALL PHORO2(-TH3,PVEC)
122 CALL PHOBO3(ANGLE,PVEC)
123 CALL PHORO3(FI4,PVEC)
124 CALL PHORO2(-TH4,PNEUTR)
125 CALL PHORO2(-TH4,PPHO(1,NPHO))
126 CALL PHORO2(-TH4,PVEC)
127 FI5=PHOAN1(PVEC(1),PVEC(2))
129 C-- Charged particle restores original direction
130 CALL PHORO3(-FI5,PNEUTR)
131 CALL PHORO3(-FI5,PPHO(1,NPHO))
132 CALL PHORO2(TH1,PNEUTR(1))
133 CALL PHORO2(TH1,PPHO(1,NPHO))
134 CALL PHORO3(FI1,PNEUTR)
135 CALL PHORO3(FI1,PPHO(1,NPHO))
136 C-- See whether neutral system has multiplicity larger than 1...
137 IF ((JDAPHO(2,IP)-JDAPHO(1,IP)).GT.1) THEN
138 C-- Find pointers to components of 'neutral' system
143 IF (I.NE.NCHARB.AND.(JMOPHO(1,I).EQ.IP)) THEN
145 C-- Reconstruct kinematics...
146 CALL PHORO3(-FI1,PPHO(1,I))
147 CALL PHORO2(-TH1,PPHO(1,I))
149 C-- ...reductive boost
150 CALL PHOBO3(PARNE,PPHO(1,I))
152 C-- Rotate in order to get photon along z-axis
153 CALL PHORO3(-FI3,PPHO(1,I))
154 CALL PHORO2(-TH3,PPHO(1,I))
156 C-- Boost to the rest frame of decaying particle
157 CALL PHOBO3(ANGLE,PPHO(1,I))
159 C-- Back in the parent rest-frame but PNEUTR not yet oriented.
160 CALL PHORO3(FI4,PPHO(1,I))
161 CALL PHORO2(-TH4,PPHO(1,I))
163 C-- Charged particle restores original direction
164 CALL PHORO3(-FI5,PPHO(1,I))
165 CALL PHORO2(TH1,PPHO(1,I))
166 CALL PHORO3(FI1,PPHO(1,I))
171 C-- ...only one 'neutral' particle in addition to photon!
173 80 PPHO(J,NEUDAU)=PNEUTR(J)
176 C-- All 'neutrals' treated, fill /PHOEVT/ for charged particle...
178 90 PPHO(J,NCHARB)=-(PPHO(J,NPHO)+PNEUTR(J))
179 PPHO(4,NCHARB)=PPHO(5,IP)-(PPHO(4,NPHO)+PNEUTR(4))