bin/BCVEGPY/phase_gen.F

0001 c********************************************************************

0002 c********************************************************************

0003 c*** for phase-space

0004 c*** phase_gen() ---- generate the phase-space

0005 c*** lorentz() ---- doing momentum transform.

0006 c********************************************************************

0007 c********************************************************************

0008
0009 c...from program rambos, which is a democratic multi-particle phase

0010 c...space generator.  authors:  s.d. ellis,  r. kleiss,  w.j. stirling

0011 c...this is version 1.0 - written by r. kleiss modified slightly by

0012 c...ian hinchliffe. here we make little changes and let it only adaptable

0013 c...to three final particals. the interesting reader can change back to

0014 c...it's original form easily.

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 c      data acc/1.0d-14/,itmax/6/,in_it/0/,pi2log/1.837877d0/

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 c count nonzero masses

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 c generate n massless momenta in infinite phase space

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 c...there is one redundant overall phi that can be chosing arbitrarily.

0062 c...by chosing the original f=0, we will always get

0063 c...q(1,1)=0, which corresponds to a specific receive direction of the

0064 c...final particles. this breaks the isotropy of the space, to avoid this,

0065 c...we may choose f=2*pi*pyr(0) .

0066 c          f=0.0d0

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 c calculate the parameters of the scale transformation

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 c calculate the n_th vector

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 c scale the q's to the p's

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 c return for weighted massless momenta

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 c massive particles: rescale the momenta by a factor x

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 c calculate the mass-effect weight factor

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 c return for weighted massive momenta

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 c...b_c

0176       pmomup(5,3)=xm(3)
0177       do j=1,4
0178            pmomup(j,3)=p(j,3)
0179       end do
0180
0181 c...b quark

0182       pmomup(5,4)=xm(2)
0183       do j=1,4
0184            pmomup(j,4)=p(j,2)
0185       end do
0186
0187 c...\bar{c} quark

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 c**********************************************************

0197
0198 c...this performs a lorentz transformation

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