File indexing completed on 2023-03-17 11:04:21
0001
0002
0003
0004
0005
0006
0007
0008
0009
0010
0011
0012
0013
0014
0015
0016 subroutine phase_gen(y,et,wt)
0017 implicit double precision (a-h,o-z)
0018 implicit integer (i-n)
0019
0020 double complex colmat,bundamp
0021 common/upcom/ecm,pmbc,pmb,pmc,fbcc,pmomup(5,8),
0022 & colmat(10,64),bundamp(4),pmomzero(5,8)
0023 common/rconst/pi
0024 dimension xm(3),p(4,3),q(4,3),r(4),z(5),
0025 & p2(3),xm2(3),e(3),v(3),y(5)
0026
0027
0028 acc=1.0d-14
0029 itmax=6
0030 in_it=0
0031 pi2log=1.837877d0
0032
0033 xm(1)=pmc
0034 xm(2)=pmb
0035 xm(3)=pmbc
0036
0037 n=3
0038
0039 if(in_it.eq.0) then
0040 in_it=1
0041 z(1)=0.0d0
0042 do k=2,5
0043 z(k)=z(k-1)+dlog(dfloat(k))
0044 end do
0045 end if
0046
0047
0048 xmt=0.0d0
0049 nm=0
0050 do i=1,n
0051 if(xm(i) .gt. 1.0d-6) nm=nm+1
0052 xmt=xmt+dble(abs(xm(i)))
0053 end do
0054
0055
0056 extra_weight=1.0d0
0057 do i=1,n-1
0058 if(i.eq.1)then
0059 c=2.0d0*y(i)-1.0d0
0060 s=dsqrt(1.0d0-c*c)
0061
0062
0063
0064
0065
0066
0067 f=2*pi*pyr(0)
0068 fac=y(i+1)
0069 else
0070 c=2.0d0*y(i+1)-1.0d0
0071 s=dsqrt(1.0d0-c*c)
0072 f=2.0d0*pi*y(i+2)
0073 fac=y(i+3)
0074 end if
0075 q(4,i)=-dlog(fac)
0076 extra_weight=extra_weight*q(4,i)
0077 q(3,i)=q(4,i)*c
0078 q(2,i)=q(4,i)*s*cos(f)
0079 q(1,i)=q(4,i)*s*sin(f)
0080 end do
0081
0082
0083 do i=1,4
0084 r(i)=0.0d0
0085 end do
0086 do i=1,n-1
0087 do k=1,4
0088 r(k)=r(k)+q(k,i)
0089 end do
0090 end do
0091
0092
0093 do k=1,3
0094 q(k,n)=-r(k)
0095 end do
0096 q(4,n)=dsqrt(q(1,n)**2+q(2,n)**2+q(3,n)**2)
0097 r(4)=r(4)+q(4,n)
0098 rmas=r(4)
0099 if(rmas .lt. 1.0d-6)then
0100 wt=0.d0
0101 return
0102 end if
0103 x=et/rmas
0104
0105
0106 do i=1,n
0107 do k=1,4
0108 p(k,i)=x*q(k,i)
0109 end do
0110 end do
0111
0112
0113 wt=pi2log*(n-1)-dlog(et-p(4,n))*2.0d0*(1-n)-dlog(et)-
0114 & dlog(2.0d0)-dlog(p(4,n))-z(2*n-3)
0115 if(nm.eq.0) then
0116 wt=extra_weight*dexp(wt)
0117 return
0118 end if
0119
0120
0121 xmax=dsqrt(1.0d0-(xmt/et)**2)
0122 do i=1,n
0123 xm2(i)=xm(i)**2
0124 p2(i) =p(4,i)**2
0125 end do
0126 iter=0
0127 x=xmax
0128 accu=et*acc
0129 do while (iter.le.itmax)
0130 f0=-et
0131 g0=0.d0
0132 x2=x*x
0133 do i=1,n
0134 e(i)=dsqrt(xm2(i)+x2*p2(i))
0135 f0=f0+e(i)
0136 if(e(i) .gt. 0.0d0)then
0137 g0=g0+p2(i)/e(i)
0138 end if
0139 end do
0140 if(dble(abs(f0)).gt.accu.and.iter.le.itmax) then
0141 iter=iter+1
0142 x=x-f0/(x*g0)
0143 else if(dble(abs(f0)).le.accu) then
0144 iter=itmax+1
0145 end if
0146 end do
0147
0148 do i=1,n
0149 v(i)=x*p(4,i)
0150 do k=1,3
0151 p(k,i)=x*p(k,i)
0152 end do
0153 p(4,i)=e(i)
0154 end do
0155
0156
0157 wt2=1.0d0
0158 wt3=0.0d0
0159 do i=1,n
0160 if(e(i) .gt. 0.0d0)then
0161 wt2=wt2*v(i)/e(i)
0162 wt3=wt3+v(i)**2/e(i)
0163 end if
0164 end do
0165 wtm=(2.0d0*n-3.0d0)*dlog(x)+dlog(wt2/wt3*et)
0166
0167
0168 wt=wt+wtm
0169 wt=dexp(wt)*extra_weight
0170 if(xmt .gt. et)then
0171 wt=0.0d0
0172 return
0173 end if
0174
0175
0176 pmomup(5,3)=xm(3)
0177 do j=1,4
0178 pmomup(j,3)=p(j,3)
0179 end do
0180
0181
0182 pmomup(5,4)=xm(2)
0183 do j=1,4
0184 pmomup(j,4)=p(j,2)
0185 end do
0186
0187
0188 pmomup(5,5)=xm(1)
0189 do j=1,4
0190 pmomup(j,5)=p(j,1)
0191 end do
0192
0193 return
0194 end
0195
0196
0197
0198
0199 subroutine lorentz(pb,pc,pl)
0200 implicit double precision (a-h,o-z)
0201 implicit integer(i-n)
0202
0203 dimension pb(4),pc(4),pl(4)
0204
0205 q =sqrt(pb(4)**2-pb(1)**2-pb(2)**2-pb(3)**2)
0206 pl(4)=(pb(4)*pc(4)+pb(3)*pc(3)+pb(2)*pc(2)+pb(1)*pc(1))/q
0207 f =(pl(4)+pc(4))/(q+pb(4))
0208 do j=1,3
0209 pl(j)=pc(j)+f*pb(j)
0210 end do
0211
0212 return
0213 end
0214