bin/BCVEGPY/phase_point.F

0001
0002 c*********************************************************************

0003
0004 c...generate allowed phase_space and right kinematical ranges(to ensure

0005 c...the subprocess has enough energy).

0006       subroutine phpoint(zup,wt)
0007 c...preamble: declarations.

0008       implicit double precision (a-h, o-z)
0009         implicit integer(i-n)
0010
0011 c...pythia common block.

0012         parameter (maxnup=500)
0013       common/hepeup/nup,idprup,xwgtup,scalup,aqedup,aqcdup,idup(maxnup),
0014      &istup(maxnup),mothup(2,maxnup),icolup(2,maxnup),pup(5,maxnup),
0015      &vtimup(maxnup),spinup(maxnup)
0016       save /hepeup/
0017
0018         common/ptpass/ptmin,ptmax,crossmax,etamin,etamax,
0019      &  smin,smax,ymin,ymax,psetamin,psetamax
0020
0021         double complex colmat,bundamp
0022         common/upcom/ecm,pmbc,pmb,pmc,fbcc,pmomup(5,8),
0023      &  colmat(10,64),bundamp(4),pmomzero(5,8)
0024
0025 c...to get the subprocess cross-section.

0026       common/subopen/subfactor,subenergy,isubonly
0027
0028       dimension zup(7),zzup(5)
0029
0030 c...initial value.

0031       wt=0.0d0
0032
0033 c--------------------------------------------

0034 c****constraints from the kinematical region.

0035
0036       taumin =((pmbc+pmb+pmc)/ecm)**2
0037         taumax =1.0d0
0038
0039 c------------------------------------------

0040       yymin = 0
0041       yymax = 0
0042       if(isubonly.eq.1) then
0043          x1=subfactor
0044            x2=subfactor
0045            if((x1*x2).lt.taumin) then
0046              write(*,'(a)') 'energy not high enough!'
0047              write(3,'(a)') 'energy not high enough!'
0048              stop
0049            end if
0050         else
0051            tau  = (taumax-taumin)*zup(6)+taumin
0052            yymin= dlog(dsqrt(tau))
0053            yymax=-dlog(dsqrt(tau))
0054            yy   = (yymax-yymin)*zup(7)+yymin
0055            x1=dsqrt(tau)*exp(yy)
0056            x2=dsqrt(tau)*exp(-yy)
0057         end if
0058
0059 c...in the c.m. system.

0060         pup(1,1) = 0.0d0
0061         pup(2,1) = 0.0d0
0062         pup(3,1) = dsqrt(x1*x2)*ecm/2.0d0
0063         pup(4,1) = dsqrt(x1*x2)*ecm/2.0d0
0064         pup(5,1) = 0.0d0
0065
0066         pup(1,2) = 0.0d0
0067         pup(2,2) = 0.0d0
0068         pup(3,2) =-dsqrt(x1*x2)*ecm/2.0d0
0069         pup(4,2) = dsqrt(x1*x2)*ecm/2.0d0
0070         pup(5,2) = 0.0d0
0071
0072 c----------------------------------------------

0073
0074       do i=1,5
0075            zzup(i)=zup(i)
0076         end do
0077
0078 c...calculate the jacobian of the final particles in c.m. frame.

0079         et =pup(4,1)+pup(4,2)
0080
0081       izero=0
0082
0083 c...increase the phase-space efficiency.

0084 321   call phase_gen(zzup,et,wt)
0085         if((wt .lt. 1.0d-16).and.(izero.lt.100000))then
0086           izero=izero+1
0087           goto 321
0088         end if
0089
0090 c--------------------------------------------

0091
0092         if(isubonly.eq.1) then
0093           return
0094         else
0095           wt=(taumax-taumin)*(yymax-yymin)*wt
0096           return
0097         end if
0098
0099 c...after calling phase_gen all the particle's momenta are

0100 c...constructed effectively, which is in the c.m.s.

0101
0102         end
0103