c
c         main program of e+ + e- -> ee+g*+g --> ee + X

         IMPLICIT REAL*8(A-H,O-Z)
      COMMON/qcevnt/partn(10,1000),listp(10,1000),NTOP
      common/prcss/ecm,xbmin,nprc
      common/qcwght/wght0,fail,wght,accptx
      common/qcpar/roots,q2max,q02,q2min,flamd,qmass(10),nflav
      common/photon/alpha,rtx,rtph,rtvm,nfl,nfvm
      common/cross/sigma0,sigmax,sigma
      common/prcgg/q2xmax,q2xmin,zmin,ame2,ecmx0
      common/prcxx/xfail
         dimension leaves(1000)

c
c
c             
c      Steps for generation of
c      virtual photon and real photon events in e+ e- process
 
c      step 1: pick up Q^2
c      step 2: pick up energy of real photon
c      step 3: make shower of electron and photon system
c      step 4: make boost of Lab. system
c

c
c       subroutine for input of parameters  
                call qcinpt(neven)
c       subroutine for preparation for shower
                call qcprep
c

c       subroutine for analysis
                call anal(0)
c
              wght0 =0.0d0
              wght0x =0.0d0
              fail = 0.0d0
             xfail = 0.0d0
c
                do 1 i=1,neven
  5     continue
            if( nprc. eq. 2)then
c        subroutine to generate a real photon 
                call ggq2(q2,ecmx,nfailx)
                   if( nfailx .eq. 1)then
                      fail = fail + 1.0d0
                      goto 5
                                     endif
                q2max = q2
                ecm = ecmx
cx        write(6,*)' q2=',q2max,'  ecm=',ecm
c
                           endif
c      subroutine for generation of shower event
                call qcgen(nfail)
                  if( nfail .eq. 1)then
                      fail = fail + 1.0d0
                       goto 5
                                   endif
c
              if(nprc .eq.2)then
c       subroutine to go to CM system
                call labsys
c
           do 47 np=1,100
 47        leaves(np)=0
             call mklist(leaves, nprtn)
                if( i .le. 20)then
      write(6,*)'color flow, address  in order'
      write(6,*)'number of final partons',nprtn
c
c      In this case leaves(1) is address of scatterd lepton.
c      Please neglect leaves(1)
              do 46 np =0,nprtn/5
 46          write(6,461)(leaves(npi+np*5+1),npi=0,4)
 461  format('  ',4(I5,2X))
                              endif
                            endif
                if( i .le. 50)call qclist
c
              wght0 =wght0 + wght
c       subroutine for analysis
                call anal(1)
c
  1           continue

      write(6,100)
      write(6,101)
 100  format('  ')
 101  format(' end of generation')
          sigma = sigmax*wght0/(float(neven)+fail)
          write(6,110)sigma
 110  format(' the total cross section =',G15.5,' pb')
cx           write(6,*)' xfail =',xfail
c       subroutine for analysis
                call anal(2)
               stop
               end
c
c       input 
c
        subroutine qcinpt(neven)
             implicit real*8(a-h,o-z)
c
      common/color/cf,ca,tr,pi,beta0
      common/qcpar/roots,q2max,q02,q2min,flamd,qmass(10),nflav
      common/prcss/ecm,xbmin,nprc
      common/qcflg/kzap,nao
      COMMON /QCRND/ NSEED
      common/prcgg/q2xmax,q2xmin,zmin,ame2,ecmx0
c
c++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
c        users have to control following four parameters
c        for any process.  
c      event number of generation.
             neven = 50

c        center mass energy
             roots =58.0d0
c++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
c        in this case users should fix the minimum value of
c        virtual mass squared of the virtual photon.
               q2xmin = 20.0d0
c++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
c
c       choice of process
c      0 : calculate parton distributions in photon.
c      1 : process of real photon and electron scattering
c      2 : process of virtual photon and real photon 
c          of electron positron scattering
              nprc = 2
c++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
c
c
cx             q2max = roots*roots
             q02   = 0.1d0
             q2min = 0.25d0
              flamd = 0.23d0
c
              qmass(1) = 0.0d0
              qmass(2) = 0.0d0
              qmass(3) = 0.5d0
              qmass(4) = 1.5d0
              qmass(5) = 5.0d0
              qmass(6) = 130.0d0
c
               nflav = 3
c      color
             cf =4.0d0/3.0d0
             ca = 3.0d0
             tr = 0.5d0 
             pi = 4.0d0*atan(1.0d0)
             beta0 = (ca*11.0d0 -4.0d0*tr*nflav)/3.0d0
c
c
             nseed  = 1234567
c
             kzap = 0
               ecm = roots
               s = ecm*ecm
               xbmin = 1.0d-5
               zmin = 1.0d-5
               ame2 =(0.5d-3)**2
              if( nprc .eq. 1)then
c
c              real prcess of electron and photon 
c           ecm : total energy of c.m. system of electron and photon
c
               ecm = roots
               s = ecm*ecm
               q2max  = s/4.0d0
               q2xmax = q2max 
               xbmin = q2max / s
                              else
               if( nprc .eq. 2)then
c
c            virtual and real photon  process of electron and positron
c            scattering
c
               ecmx0 = roots
               ecm = ecmx0
                s = roots*roots
               q2max  = s/4.0d0
               q2xmax = q2max 
               zmin = q2xmin/s
               ame2 =(0.5d-3)**2
c
                              else
c
c           generate parton distributions of photon
c
                 xbmin = 0.0d0
               q2max = q2xmin
               q2xmax = q2max 
                              endif
                              endif
c
c
       WRITE(6,101)
 101   FORMAT('   Parton shower in  electron photon scattering'/
     & ,'(virtual photon and real photon process)')
       WRITE(6,110)nprc
 110   FORMAT(' real pro(ee):2, real pro.(e-gamma).:1',
     & 'parton distri.:0 ',I3)
       WRITE(6,115)ecm
 115   FORMAT(' total energy',G15.5)
       write(6,120)-q2max
 120   format(' maximum virtual mass squared of the photon  =',G15.5)
       write(6,125)-q2xmin
 125   format(' minimum virtual mass squared of the photon  =',G15.5)
       write(6,130)flamd,nflav
 130   format('  QCD param.  Lambda=',G15.5,'  no. of flavor =',I5)
       write(6,140)q02,q2min
 140   format(' cutoff mass squared =',G15.5,' initial Q**2 =',G15.5)
       write(6,150)neven
 150   format('  event number =',I8)
       write(6,160)kzap  
 160   format('  argument of alpha-s is K^2(0) or Pt^2(1) :',I5)
    
         if( nprc .eq. 2)then
        write(6,201)q2xmin,q2xmax
 201    format('  minimum and maximum of virtual photon',2G15.5)
        write(6,211)zmin
 211    format('mini. ener. frac. of real photon in positron ', G15.5)
                          endif
             return
             end
c
                subroutine ggq2(q2,ecmx,nfailx)
         IMPLICIT REAL*8(A-H,O-Z)
        COMMON/qcevnt/partn(10,1000),listp(10,1000),NTOP
      common/prcss/ecm,xbmin,nprc
      common/qcwght/wght0,fail,wght,accptx
      common/qcpar/roots,q2max,q02,q2min,flamd,qmass(10),nflav
      common/photon/alpha,rtx,rtph,rtvm,nfl,nfvm
      common/prcgg/q2xmax,q2xmin,zmin,ame2,ecmx0
      
               nfailx = 0
c
cx 5      continue
              theta  = 1.0d-1
c
c
c       pick up Q^2
c
           rm = 1.0d0/q2xmin - 1.0d0/q2xmax
c
       rq = rm*random(1) + 1.0d0/q2xmax
             q2 = 1.0d0/rq
c
           zmax = 1.0d0
          zlmin = log(zmin)
          zlmax = log(zmax)
          ame = sqrt(ame2)
          
          fxm=2.0d0*log(ecmx0*theta/ame/zmin)
           do 201 i=1,100
          zl = (zlmax-zlmin)* random(1)+zlmin
          z = exp(zl)
          fx = (1.0d0+(1.0d0-z)**2)
          fx = fx * log((1.0d0-z)/z*ecmx0*theta/ame)
c
           if( fx .ge. fxm*random(1))goto 211
 201      continue
 211      continue
cxxxx
cxxx              z =0.03
c
cx           write(6,*)'zmin,zmax,z =',zmin,zmax,z
          ecmx2 = z* ecmx0**2
          ecmx = sqrt(ecmx2)
c
         xbmin = q2/ecmx**2
cx          write(6,*)'q2,ecmx,xbmin,ecmx0= ',q2,ecmx,xbmin,ecmx0
cx           write(6,*)'zmin,zmax,z =',zmin,zmax,z
c
         if( xbmin .ge. 1.0d0)then
              nfailx = 1
ccx               goto 5
                              endif
          return
          end
c
                subroutine labsys
         IMPLICIT REAL*8(A-H,O-Z)
      common/color/cf,ca,tr,pi,beta0
      COMMON/qcevnt/partn(10,1000),listp(10,1000),NTOP
      common/prcss/ecm,xbmin,nprc
      common/qcwght/wght0,fail,wght,accptx
      common/qcpar/roots,q2max,q02,q2min,flamd,qmass(10),nflav
      common/photon/alpha,rtx,rtph,rtvm,nfl,nfvm
      common/prcgg/q2xmax,q2xmin,zmin,ame2,ecmx0
      dimension pppp(4,100),pref(4)
c
              z= ecm**2/ecmx0**2
c
c       write(6,*)' ecmx0,ecm,z =',ecmx0,ecm,z
c
           l1 = ntop + 1
         partn(2,L1) = 0.0d0
         partn(3,L1) =  0.0d0
         partn(4,L1) =  0.0d0
         partn(5,L1) =  0.0d0
         partn(6,l1) =  ecmx0/2.0d0
         partn(7,l1) =  ecmx0/2.0d0

         listp(1,L1) = -11
         listp(5,L1) = 2
         listp(6,L1) = -2

          
           l2 = ntop + 2
         partn(2,L2) = 0.0d0
         partn(3,L2) =  0.0d0
         partn(4,L2) =  0.0d0
         partn(5,L2) =  0.0d0
         partn(6,l2) = (1.0d0-z)*ecmx0/2.0d0
         partn(7,l2) = (1.0d0-z)*ecmx0/2.0d0

         listp(1,L2) = -11
         listp(5,L2) = 2
         listp(6,L2) = 2
c
c          call qclist
c         boost to center mass system of electron and positron
          pref(1)=0.0d0
          pref(2)=0.0d0
          pref(3)=-(1.0d0 - z)*ecmx0/2.0d0
          pref(4)= (1.0d0+z)*ecmx0/2.0d0
         do 201 i=1,ntop
         do 201 k=1,4
  201   pppp(k,i)=partn(3+k,i)
         nnn1 = 1
         nnn2 = ntop
      CALL BOOST( PPPP , PREF , NNN1 , NNN2 )
c
         do 211 i=1,ntop
         do 211 k=1,4
  211   partn(3+k,i) = pppp(k,i)
c
          ntop = ntop + 2
c
c         symmetric treatment
c
          xx = random(1)
          if( xx .ge. 0.5d0)then
          do 221 i=1,ntop
         partn(6,i) = -partn(6,i)
         id=listp(1,i)
         if(abs(id).eq.11)listp(1,i)=-id
 221     continue
                           endif
c         symmetric treatment on the azimuthal angle of lepton
c
          phi = 2.0d0*pi*random(1)
          cs = cos(phi)
          sn = sin(phi)
          do 231 i=1,ntop
         px =partn(4,i) 
         py =partn(5,i) 
         partn(4,i) = px*cs -py*sn
         partn(5,i) = px*sn +py*cs
 231     continue

          return
          end
c
              subroutine qcgen(nfailx)
         IMPLICIT REAL*8(A-H,O-Z)
        COMMON/qcevnt/partn(10,1000),listp(10,1000),NTOP
      common/prcss/ecm,xbmin,nprc
      common/qcwght/wght0,fail,wght,fail2
      common/prcgg/q2xmax,q2xmin,zmin,ame2,ecmx0
      common/prcxx/xfail
     common/qcpar/roots,q2max,q02,q2min,flamd,qmass(10),nflav
c
             nfailx = 0 
  5   continue
c
             do 1 i=1,100
             do 1 j=1,10
             partn(j,i)=0.0d0
             listp(j,i)=0
 1       continue
c
             ntop = 0
             call rslv(nfailr)
c
              if( nfailr .eq. 1)then
                 if( nprc.ne.2)then
                  write(6,*)'nfailr,nprc =',nfailr,nprc
                               endif
                   xfail = xfail + 1.0d0
                   nfailx = 1
                   return
                              endif
             call psint
c
             if( nprc .ge. 1)then
c
             call scatt(nfail)
c
              if( nfail .eq. 1)then
              if( nprc .eq. 2)then
                   nfailx = 1
                   return
                              endif
                  fail = fail + 1.0d0
                            goto 5
                                endif
c
c
             call psfnl
c
             call cms(nfail)
c
              if( nfail .eq. 1)then
        if( nprc .eq. 2)then
                   nfailx = 1
                   return
                        endif
                  fail = fail + 1.0d0
                  fail2 = fail2 + 1.0d0
                            goto 5
                                endif
c
                             endif
c
             return
             end
c
            subroutine rslv(nfail)
c       resolved photon or quark pair production
c
C   Q-EVOLUTION USING ALTARELLI-PARISI EQUATION
         IMPLICIT REAL*8(A-H,O-Z)
      COMMON /QCRND/ NSEED
      COMMON/qcevnt/partn(10,1000),listp(10,1000),NTOP
      common/color/cf,ca,tr,pi,beta0
      common/qcpar/roots,q2max,q02,q2min,flamd,qmass(10),nflav
      common/photon/alpha,rtx,rtph,rtvm,nfl,nfvm
      common/phovmu/vmuaa,vmua,vmub,vmuc,vxmu,vmunr
      common/phovmg/vmgaa,vmga,vmgb,vmgc,vxmg,vmgnr
      common/prcss/ecm,xbmin,nprc
      common/prcgg/q2xmax,q2xmin,zmin,ame2,ecmx0
      common/rslx/fph0,fph0x
c
cx           write(6,*)' q2max,q2min in resl =',q2max,q2min
          nfail = 0
c
         if( nprc .eq. 0)then
        L1= 1
          partn(1,L1) = 1.0d0
          listp(1,L1) =  22
          listp(2,L1) = 0
          listp(5,L1) = 0
          listp(6,L1) = -1

          ntop = 1
                        else
        L1= 1
          partn(1,L1) = 1.0d0
          partn(1,l1) = 0.0d0
          listp(1,L1) =  11
          listp(2,L1) = 0
          listp(5,L1) = 1
          listp(6,L1) = -2

        L1= 2
          partn(1,L1) = 1.0d0
          listp(1,L1) =  22
          listp(2,L1) = 0
          listp(5,L1) = 0
          listp(6,L1) = -1

          ntop = 2
                        endif
c

           alam2 = flamd**2
           qq0 = q2min
cx1           write(6,*)'q2max,q2min,nprc =',q2max,q2min,nprc  
         if( nprc .lt. 2)then
           qq  = q2max
                        else
           qq = q2xmax
                        endif
c
         rn = random(1)
          if( rn .le. rtph)then
c       resolved to quark-antiquark
c
c    fph0x : u-,d-,s-quark
c    fph0  : u-,d-,s-,c-quark
           if( random(1) .ge. fph0x/fph0)then
                id0 =4
         qqx = qmass(4)**2 * exp( rn/rtph *log(qq/qmass(4)**2))
         y0 = log(qqx/alam2)
                                         else
         xrn = random(1)*rtph
cx         qqx = qq0 * exp( rn/rtph *log(qq/qq0))
         qqx = qq0 * exp( xrn/rtph *log(qq/qq0))
           if( qqx .ge. qq)write(6,*)' qqx > qq, qqx,qq =',qqx,qq
         y0 = log(qqx/alam2)
c       nfl = 3 (u,d,s) case
c         select flavor according to chrage squared
c        charm quark included
              reu2 =4.0/(4.0+1.0+1.0)
              red2 = reu2+ 1.0/(4.0+1.0+1.0)
        xid =  random(1)
         if( xid .le. reu2)then
                          id0 = 1
                          elseif( xid.le.red2)then
                                 id0 =2
                          else
                                 id0= 3
                           endif
                                         endif
cx         write(6,*)' id0,qqx,mass**2(c)=',id0,qqx,qmass(4)**2
ccc    if charm is selected, check the kinematics
         if( id0.eq.4  .and. qqx.lt. qmass(4)**2)then
       write(6,*)' id0,qqx,mass**2(c)=',id0,qqx,qmass(4)**2
            stop
                                                endif

        if( id0 .gt. nfl)write(6,*)'id0,nfl= ',id0,nfl
        if( random(1) .ge. 0.5) id0 = -id0
          fxm = 1.0d0
          xmx =qqx/(qqx+qmass(abs(id0))**2)

         do 101 i=1,100
         x = xmx*random(1)
         fx = x*(x**2 + (1.0d0-x)**2)
         if( fx .ge. fxm)write(6,*)' fx > fxm in rslv.',fx,fxm
         if( fx .ge. fxm*random(1))goto 102
 101     continue
 102     xb  = x
c
c    check the kinematics
c
          xxx =qqx/xb - qmass(abs(id0))**2/(1.0d0-xb)
          if(xxx .lt. 0.0d0)then
c       cancel resolve
                   nfail =1
cx         write(6,*)'fail at rslv, qqx,xb,xxx,id0',qqx,xb,xxx,id0
cx          if(abs(id0).eq.4)then
cx          if(abs(id0).ge.3)then
cx          xb=0.5
cx            nfail= 0
cx           goto 881
cx                           endif
                   return
                           endif



c
          nad = ntop
          L1= NTOP+1
          partn(1,L1) =xb 

          listp(1,L1) = ID0
          listp(2,L1) = L1- NAD

          L2 = NTOP+2
          partn(1,L2) =(1.0D0-Xb)

          listp(1,L2) = -id0
          listp(2,L2) = L2- NAD
          listp(3,L2) = 0

          listp(3,NAD) = 2
          listp(4,NAD) = L1 - NAD
          listp(5,NAD) = 1

          partn(2,L1) = - qqx
          partn(2,L2) = qmass(abs(id0))**2
          listp(3,L1) = 0
          listp(3,L2) = 0
          listp(4,L1) = 0
          listp(4,L2) = 0
          listp(5,L1) = 0
          listp(5,L2) = 1
          listp(6,L1) = -1
          listp(6,L2) = 1

          ntop =ntop + 2

                           else
c     change to vector meson
c
         qqx = qq0 
         y0 = log(qqx/alam2)
          if( random(1) .le. vmunr/(vmunr+vmgnr))then
c     pick up quark
        id0 = nfvm * random(1) + 1
        if( id0 .gt. nfvm)write(6,*)'id0,nfvm= ',id0,nfvm
        if( random(1) .ge. 0.5d0) id0 = -id0
           do 111 i=1,100
           x = random(1)
c   quark distribution
        vx = x**vmua*(vmuaa+vmub*sqrt(x))*(1.0d0-x)**vmuc 
         if(  vx .ge. vxmu)write(6,*)' vx > vm ',vx,vxmu
         if( vx .ge. vxmu*random(1))goto 112
 111     continue
 112     xb  = x
                                             else
c   pick up gluon
         id0 = 0
           do 113 i=1,100
           x = random(1)
c   gluon distribution
        vx = x**vmga*(vmgaa+vmgb*sqrt(x))*(1.0d0-x)**vmgc 
         if(  vx .ge. vxmg)write(6,*)' vx > vm ',vx,vxmg
         if( vx .ge. vxmg*random(1))goto 114
 113     continue
 114     xb  = x
                                             endif
c
          nad = ntop
          L1= NTOP+1
          partn(1,L1) =xb 

          listp(1,L1) = ID0
          listp(2,L1) = L1- NAD

          L2 = NTOP+2
          partn(1,L2) =(1.0D0-Xb)

          listp(1,L2) = -id0
          listp(2,L2) = L2- NAD
          listp(3,L2) = 0

          listp(3,NAD) = 2
          listp(4,NAD) = L1 - NAD
          listp(5,NAD) = 1

          partn(2,L1) = - qqx
          partn(2,L2) = 0.0d0
          listp(3,L1) = 0
          listp(3,L2) = 0
          listp(4,L1) = 0
          listp(4,L2) = 0
          listp(5,L1) = 0
          listp(5,L2) = 1
          listp(6,L1) = -1
          listp(6,L2) = 1

          ntop =ntop + 2
c
               endif
c
c
c
              q2 = qqx
              xb0 = xb
              ifl = id0
              return
              end
c
c
         subroutine scatt(nfail)
         IMPLICIT REAL*8(A-H,O-Z)
        COMMON/qcevnt/partn(10,1000),listp(10,1000),NTOP
      common/prcss/ecm,xbmin,nprc
      common/qcwght/wght0,fail,wght,fail2
c
           nfail = 0
           weight = 1.0d0
c
           do 101 i=1,ntop
         if(listp(3,i).eq.0 .and. listp(6,i).eq.-1)then
                  lp =i
                  goto 102
                                                   endif
 101     continue
         write(6,*)' error in scatt, not fixed lp'
 102     continue
           id = listp(1,lp)
c
c      check kinematical contraint  xb > xbmin
c
           if( partn(1,lp) .le. xbmin)then
                 nfail = 1
                 return
                                       endif
         if( id .eq. 0)then
c      spacelike parton is gluon ... not couple with photon
              nfail = 1
              return
                      else
        ida = abs(id)
        if( ida.eq.1 .or. ida.eq.4  .or. ida.eq.6)then
          weight = 1.0d0
                                                  else
         weight = 0.25d0
          if( weight .le. random(1))then
                   nfail = 1
                   return
                                    endif
                                                  endif
                     endif
c
c      the scattered electron
c
           l1 = ntop + 1
c
           listp(1,l1) = 11
           listp(2,l1) = l1 - 1
           listp(3,l1) =  0
           listp(4,l1) =  0
           listp(5,l1) =  1
           listp(6,l1) =  2
c
           listp(3,1) = 1
           listp(4,1) = l1-1
c
           partn(1,l1) = 1.0d0
c
           ntop = ntop + 1
c
c      the scattered parton 
c
           l1 = ntop + 1
c
           listp(1,l1) =  id
           listp(2,l1) = l1 - lp
           listp(3,l1) =  0
           listp(4,l1) =  0
           listp(5,l1) =  0
           listp(6,l1) =  1
c
           listp(3,lp) = 1
           listp(4,lp) = l1-lp
c
           partn(1,l1) = 1.0d0
c
           ntop = ntop + 1
c
cx        write(6,*)' end of scatt. '
           return
           end
C------------------------------------------------------------------
C         RANDOM
C------------------------------------------------------------------
C         RANDOM NUMBER  ( 0- 1.0 , UNIFORM )
C------------------------------------------------------------------
C
      FUNCTION RANDOM( IDUM  )
      IMPLICIT REAL*8(A-H,O-Z)
C
      COMMON /QCRND/ NSEED
      NSEED = NSEED *48828125
      IF(NSEED.LT.0) NSEED=(NSEED+2147483647) + 1
      RANDOM=NSEED*0.4656613D-9
      RETURN
      END
         SUBROUTINE qclist 

         IMPLICIT REAL*8(A-H,O-Z)
        DIMENSION PSUM(3)
        COMMON/qcevnt/partn(10,1000),listp(10,1000),NTOP
             common/color/cf,ca,tr,pi,beta0
             common/qcpar/roots,q2max,q02,q2min,flamd,qmass(10),nflav

C
        WRITE(6,101)NTOP
 101    FORMAT(' LIST OF OUTPUT, # OF LISTS=',I5)
        WRITE(6,102)
 102    FORMAT('     ID',6X,'PX',8X,'PY',8X,
     &  'PZ',8X,'E',5X,'# of ch. ',' sp.(tm.)')
cx        WRITE(6,102)
cx 102    FORMAT('     ID',6X,' X ',6X,'Q2',8X,'PT2',8X,'PX',8X,'PY',8X,
cx     &  'PZ',8X,'E',7X,'listp(I,*)=',
cx     &  ' I=2',' I=3',' I=4',' I=5',' I=6')


        DO 110 I =1,NTOP
        WRITE(6,103)I,listp(1,I),(partn(K,I),K=4,7),
     &              listp(3,I),listp(6,i)
cx        WRITE(6,103)I,listp(1,I),(partn(K,I),K=1,7),(listp(K,I),K=2,7)
C
 110    CONTINUE

 103    FORMAT(' ',I3,I3,2X,4(G9.2,1X),2X,I3,2X,I3)
cx 103    FORMAT(' ',I3,I3,2X,7(G9.2,1X),3X,6I3)
C
        DO 203 I=1,3
 203    PSUM(I)=0.0d0
C
        DO 201  J=1,NTOP
      IF( listp(3,J).EQ.0 .AND. listp(6,J) .GT. 0)THEN
c        write(6,*)' J in qclist',j
        DO 202 I=1,3
 202    PSUM(I)=PSUM(I) + partn(3+I,J)
                                                  ENDIF
 201   CONTINUE
C
c      subtract  initial momentum
cx       psum(1)=psum(1)-partn(4,1)-partn(4,2)
cx       psum(2)=psum(2)-partn(5,1)-partn(5,2)
cx       psum(3)=psum(3)-partn(6,1)-partn(6,2)
       PSUM2 = 0.0
        DO 204 I=1,3
 204    PSUM2  =PSUM(I)**2 + PSUM2
C
        IF( PSUM2 .GT. 1.0d-5)THEN
        WRITE(6,*)' MOMENTUM CONSERVATION ???'
        WRITE(6,*)'  MOMENTUM SUM =',(PSUM(I),I=1,3)
                              ENDIF
C

        RETURN
        END
          SUBROUTINE   psint 
         IMPLICIT REAL*8(A-H,O-Z)
        COMMON/qcevnt/partn(10,1000),listp(10,1000),NTOP
             common/color/cf,ca,tr,pi,beta0
             common/qcpar/roots,q2max,q02,q2min,flamd,qmass(10),nflav
        common/prcss/ecm,xbmin,nprc
        common/qcflg/kzap,nao

C
C      branching

cxx         WRITE(6,*)'  ntop in sub. shower before branch =',NTOP
cx         CALL qclist
          NAD = 1 

   2      CONTINUE
c          write(6,*)' listp(5:6,nad)=',listp(5,nad),listp(6,nad)
          IF(listp(5,NAD).EQ.0 .AND. listp(6,NAD).LT.0)THEN
          CALL PSGEN1(NAD)
                                   ENDIF
           NAD = NAD + 1
          IF( NAD .LE. NTOP) GOTO 2


C      WRITE(6,111)
 111   FORMAT('   END OF PSINT ')
         RETURN
         END


          SUBROUTINE PSGEN1(NAD)
         IMPLICIT REAL*8(A-H,O-Z)

        COMMON/qcevnt/partn(10,1000),listp(10,1000),NTOP
             common/color/cf,ca,tr,pi,beta0
             common/qcpar/roots,q2max,q02,q2min,flamd,qmass(10),nflav
            common/tpnes/pnes(0:100,0:10),tmax,tmin,ndim
            common/prcss/ecm,xbmin,nprc
        common/qcflg/kzap,nao
 
         REAL*8 K2


C      start of the parton shower

         Q2MX = Q2MAX
         q2mn =  abs(partn(2,nad))
cxxxx
cx         q2mx = q2mn
cxxxx
         ID   = listp(1,NAD)
c
         flamd2 = flamd**2

cx       WRITE(6,*)'  Q2MX,Q2MN = ',Q2MX,Q2MN
         if( abs(id) .gt. 4)then
        write(6,*)' flavor id =',id
               call qclist
                stop 
                           endif
c       write(6,*)' size of array =',ndim
c        WRITE(6,113)
c       WRITE(6,112)(pnes(I,abs(id)),I=1,ndim)
c 112   FORMAT('  ',4(G12.3,3X))
c 113   FORMAT('  ARRAY OF THE SUDAKOV FORM FACTOR')
  5      CONTINUE
         T1 = LOG(Q2MX / flamd2 )
         T2 = LOG(Q2MN / flamd2 )
c       WRITE(6,*)'  T1,T2     = ',T1,T2
c       WRITE(6,*)'TMIN,TMAX   = ',TMIN,TMAX

         LT1 = (T1 - TMIN )/( TMAX - TMIN )*ndim
         LT2 = (T2 - TMIN )/( TMAX - TMIN )*ndim + 1
c       WRITE(6,*)'  LT1,LT2     = ',LT1,LT2
         psdk1 = pnes(lt1,abs(id))
         psdk2 = pnes(lt2,abs(id))
c
         IF( LT1 .LE. LT2 ) GOTO 101

          P0 =  ( PSDK2 - PSDK1 )
c       WRITE(6,*)'PSDK1,PSDK2 = ',PSDK1,PSDK2
c       ****  pr < 0 *****
          Pr = LOG( RANDOM(1) )
c      WRITE(6,*)' P0,   Pr=  ',P0,Pr
C
          IF ( Pr .LE. P0 ) GOTO 101
           Px =   PSDK2 - pr

          DO 201 LT =LT2, LT1
          IF( pnes( LT ,abs(id)) .ge. Px )GOTO 202
  201     CONTINUE

  202     CONTINUE


c
          CLT = LT-1 +(Px -pnes( LT-1,abs(id)) )
     &       /(pnes(LT,abs(id))-pnes(LT-1,abs(id)))
          T = CLT*((TMAX-TMIN)/ndim ) + TMIN
          K2 = flamd2*EXP(T)
c            WRITE(6,*)'LT1,LT2         = ',LT1,LT2
c            WRITE(6,*)'LT,CLT,T,k2 = ',lt,clt,t,k2
         IF( K2 .LE. Q2min)THEN
            WRITE(6,*)'LT1,LT2         = ',LT1,LT2
            WRITE(6,*)'LT,CLT,T,Q02,P2 = ',LT,CLT,T,Q02,P2
                             ENDIF
         IF( ABS(Q2min/K2). LE. 1.E-20)THEN
            WRITE(6,*)'Q02,K2          = ',Q02,K2
            WRITE(6,*)'LT1,LT2         = ',LT1,LT2
            WRITE(6,*)'LT,CLT,T,Q02,P2 = ',LT,CLT,T,Q02,P2
                             ENDIF

C
C
cxxx       FUNCTION    IDSEL(Y,YC, ID ,NF)
c      determine flavor of spacelike child
c
        nf = nflav
          if( k2 .le. 0.0d0)then
                 write(6,*)' k2 =',k2
                              stop 
                             endif
          y = log(k2/flamd2)
          yc = log(q02/flamd2)

          IF(  ID.EQ. 0)THEN
          if(kzap .eq. 0)then
          FY1=Y -  YC  -11.0d0/12.0d0
                         else
          FY1=log(y) - log(yc) -11.0d0/12.0d0
                         endif
          FY2=     1.0d0/3.0d0
          FY0= 6.0d0*FY1 + FY2*NF*2.0d0*TR
CX     FACTOR 2 COMES FROM Q-SPACE AND Q-BAR-SPACE
          FY1=FY1*6.0d0/FY0
         IF( FY1  .GE.  RANDOM(1) )THEN
                  IDNEW = 0
                                   ELSE
         IF( 0.5d0  .GE.  RANDOM(1) )THEN
                  IDNEW = NF*RANDOM( 1 ) +1
                                   ELSE
                  IDNEW = -(  NF*RANDOM( 1 ) +1 )
                                   ENDIF
C         WRITE(6,*)' NEW ID =',IDNEW
                                   ENDIF
                       ELSE
          if(kzap .eq. 0)then
          FY1= 2.0d0*Y -2.0d0*YC-17.d0/6.0d0
                         else
          FY1= 2.0d0*log(y) -2.0d0*log(yc)-17.0d0/6.0d0
                         endif
          FY2=4.0d0/3.0d0
          FY0=  (FY1 + FY2 )
C         WRITE(6,*)'  Y , YC ,LOG(Y) =',Y,YC,LOG(Y)
C         WRITE(6,*)' FY1, FY2,  FY0 =',FY1,FY2,FY0
          FY1=FY1  /FY0
         IF( FY1  .GE.  RANDOM(1) )THEN
                  IDNEW=ID
                                   ELSE
                  IDNEW =0
                                   ENDIF
                   ENDIF

c     idnew : flavor of spacelike child
c
c      determine x
c        FUNCTION  SPLIT(Y,YC,ID,IDNEW)
C
          IF( ID.EQ.0  )THEN
c     gluon branching
                       IF(   IDNEW.EQ.0   )THEN
c     gluon branches to gluon(s) and gluon(t)
           XLXMX =Y -YC
           DO 603 I = 1 , 100
             XLX = XLXMX*RANDOM(1)
                  X = EXP(-XLX)
       if(kzap .eq. 0)then
       AP =1.0d0-2.0d0*X*(1.0d0-X)+X*X*(1.0d0-X)**2
              IF( AP .GE. 1.2d0*RANDOM(1) ) GOTO 630
                      else
        AP =1.0d0 -2.0d0*X*(1.0d0-X)+X*X*(1.0d0-X)**2
        AP =1.0d0/log(x*k2/flamd2)*ap
        apx =1.2d0/log(q02/flamd2)
              IF( AP .GE. apx*RANDOM(1) ) GOTO 630
                      endif
  603    CONTINUE
  630   CONTINUE
                  SPLIT   = 1.0d0 - X
                                           ELSE
c      gluon branches to quark-antquark
           DO 604 I = 1 , 100
                  X = RANDOM(1)
        AP = X**3 +(1.0d0-X)**2*X
              IF( AP .GE.   RANDOM(1) ) GOTO 640
  604    CONTINUE
  640    CONTINUE
                  SPLIT   = X
                                           ENDIF
                       ELSE
c
c       quark branching
                       IF(   IDNEW.NE.0   )THEN
c       quark branches into quark(s) and gluon(t)
           XLXMX =Y -YC
           DO 605 I = 1 , 100
             XLX = XLXMX*RANDOM(1)
                  X = EXP(-XLX)
        if(kzap .eq. 0)then
        AP =(1.0d0 + (1.0d0-X)**2 )*(1.d0-X)
              IF( AP .GE.   2.0d0*RANDOM(1) ) GOTO 650
                     else
        AP =(1.0d0 + (1.0d0-X)**2 )*(1.d0-X)
        AP =1.0d0/log(x*k2/flamd2)*ap
        apx =2.0d0/log(q02/flamd2)
              IF( AP .GE. apx*RANDOM(1) ) GOTO 650
                     endif
  605    CONTINUE
  650    CONTINUE

                  SPLIT   = 1.0d0 - X
                                           ELSE
c       quark branches into quark(t) and gluon(s)
           DO 606 I = 1 , 100
             X   =  RANDOM(1)
        AP = 2.d0 - 2.d0*X  + X**2
              IF( AP .GE.   2*RANDOM(1) ) GOTO 660
  606    CONTINUE
  660    CONTINUE
                  SPLIT   = X
                                           ENDIF
                   ENDIF
              x = split
C------------------------------------------------------------------

C         check the constraint is satisfied or not

             if( nprc .ge. 1)then
          IF(  (1.0D0-X) .LT. K2/Q2MAX ) THEN
cx               write(6,*)'x<t, x,k2=',x,k2
                                         Q2MN = K2
                                         GOTO 5
                                        ENDIF
                             endif                                
           idx = id - idnew
          if( idx .ne. 0)then
             xxx=-q2mn+k2/x-qmass(abs(idx))**2/(1.0d0-x)
          IF(  xxx .LT. 0.0d0 ) THEN
                                 Q2MN = K2
                                 GOTO 5
                                ENDIF
                          endif


C         write partn and listp of the children


          L1= NTOP+1
          partn(1,L1) = X*partn(1,NAD)

          listp(1,L1) = IDnew
          listp(2,L1) = L1- NAD

          L2 = NTOP+2
          partn(1,L2) =(1.0D0-X)*partn(1,NAD)
cxx          if( idnew .eq. id)then
cxx          listp(1,L2) = 0
cxx                            else
cxx          if(id .eq. 0)then
cxx           listp(1,l2) = -idnew
cxx                       else
cxx            listp(1,l2) = id
cxx                        endif
            listp(1,l2) = id - idnew
          listp(2,L2) = L2- NAD
          listp(3,L2) = 0

          listp(3,NAD) = 2
          listp(4,NAD) = L1 - NAD
          listp(5,NAD) = 1

          partn(2,L1) = - K2
            idx = id - idnew
              if( idx .eq.0)then
          partn(2,L2) = 0.0d0
                             else
          partn(2,L2) = qmass(abs(idx))**2
                             endif
          listp(3,L1) = 0
          listp(3,L2) = 0
          listp(4,L1) = 0
          listp(4,L2) = 0
          listp(5,L1) = 0
          listp(5,L2) = 0
          listp(6,L1) = -1
          listp(6,L2) = 1

           NTOP = NTOP + 2
c       CALL qclist
c       WRITE(6,111)
          RETURN

C         not branch
 101      CONTINUE
c        WRITE(6,*)' NOT BRANCH'

          listp(3,NAD) = 0
          listp(4,NAD) = 0
          listp(5,NAD) = 1
c       WRITE(6,111)
 111   FORMAT('   END OF PSGEN1')

          RETURN
          END


          SUBROUTINE   psfnl
         IMPLICIT REAL*8(A-H,O-Z)
        COMMON/qcevnt/partn(10,1000),listp(10,1000),NTOP
        common/color/cf,ca,tr,pi,beta0
        common/qcpar/roots,q2max,q02,q2min,flamd,qmass(10),nflav
        common/prcss/ecm,xbmin,nprc
        common/qcflg/kzap,nao

C
C      branching


          NAD = 1
   4      CONTINUE

          IF(listp(5,NAD).EQ.0 .AND. listp(6,NAD).eq.1)THEN
          CALL LIMQ2(Q2MX,NAD)
cxx         WRITE(6,*)'  Q2MX IN PSFNL ',Q2MX,'  NAD =',NAD
          CALL PSGEN2(Q2MX,NAD)
                                   ENDIF
           NAD = NAD + 1
          IF( NAD .LE. NTOP) GOTO 4


cx      WRITE(6,111)
 111   FORMAT('   END OF PSfnl ')
         RETURN
         END


          SUBROUTINE PSGEN2(q2mx0,NAD)
         IMPLICIT REAL*8(A-H,O-Z)
        COMMON/qcevnt/partn(10,1000),listp(10,1000),NTOP
        common/color/cf,ca,tr,pi,beta0
        common/qcpar/roots,q2max,q02,q2min,flamd,qmass(10),nflav
        common/tpnet/pnet(0:100,0:10),tmax,tmin,ndim
        common/prcss/ecm,xbmin,nprc
        common/qcflg/kzap,nao
 
         REAL*8 K2


C      start of the parton shower

         Q2MX = q2mx0
         q2mn =  q2min
         ID   = listp(1,NAD)
c
         flamd2 = flamd**2
c
          nbrn = 0
           if( nbrn .eq. 1)goto 101
        if(nbrn.eq.2.and.listp(6,nad-listp(2,nad)).eq.1)goto 101
       if( q2mx .le. q2mn)goto 101
cx       WRITE(6,*)'  Q2MX,Q2MN = ',Q2MX,Q2MN
         if( abs(id) .gt. 4)then
        write(6,*)' flavor id =',id
               call qclist
                stop 
                           endif
c       write(6,*)' size of array =',ndim
c        WRITE(6,113)
c       WRITE(6,112)(pnet(I,abs(id)),I=1,ndim)
c 112   FORMAT('  ',4(G12.3,3X))
c 113   FORMAT('  ARRAY OF THE SUDAKOV FORM FATOR')
  5      CONTINUE
         T1 = LOG(Q2MX / flamd2 )
         T2 = LOG(Q2MN / flamd2 )
cx       WRITE(6,*)'  T1,T2     = ',T1,T2
cx       WRITE(6,*)'TMIN,TMAX   = ',TMIN,TMAX,'  ndim =',ndim

         LT1 = (T1 - TMIN )/( TMAX - TMIN )*ndim
         LT2 = (T2 - TMIN )/( TMAX - TMIN )*ndim + 1
cx       WRITE(6,*)'  LT1,LT2     = ',LT1,LT2
         psdk1 = pnet(lt1,abs(id))
         psdk2 = pnet(lt2,abs(id))
c
         IF( LT1 .LE. LT2 ) GOTO 101

          P0 =  ( PSDK2 - PSDK1 )
cx       WRITE(6,*)'PSDK1,PSDK2 = ',PSDK1,PSDK2
c       ****  pr < 0 *****
          Pr = LOG( RANDOM(1) )
cx      WRITE(6,*)' P0,   Pr=  ',P0,Pr
C
          IF ( Pr .LE. P0 ) GOTO 101
           Px =   PSDK1 + pr

          DO 201 LT =LT2, LT1
          IF( pnet( LT ,abs(id)) .ge. Px )GOTO 202
  201     CONTINUE

  202     CONTINUE


c
          CLT = LT-1 +(Px -pnet( LT-1,abs(id)) )
     &       /(pnet(LT,abs(id))-pnet(LT-1,abs(id)))
          T = CLT*((TMAX-TMIN)/ndim ) + TMIN
          K2 = flamd2*EXP(T)
cx            WRITE(6,*)'LT1,LT2         = ',LT1,LT2
cx            WRITE(6,*)'LT,CLT,T,k2 = ',lt,clt,t,k2
         IF( K2 .LE. Q2min)THEN
            WRITE(6,*)'LT1,LT2         = ',LT1,LT2
            WRITE(6,*)'LT,CLT,T,Q02,P2 = ',LT,CLT,T,Q02,P2
                             ENDIF
         IF( ABS(Q2min/K2). LE. 1.E-20)THEN
            WRITE(6,*)'Q02,K2          = ',Q02,K2
            WRITE(6,*)'LT1,LT2         = ',LT1,LT2
            WRITE(6,*)'LT,CLT,T,Q02,P2 = ',LT,CLT,T,Q02,P2
                             ENDIF

C
C
cxxx       FUNCTION    IDSEL(Y,YC, ID ,NF)
c      determine flavor of  child
c
        nf = nflav
          if( k2 .le. 0.0d0)then
                 write(6,*)' k2 =',k2
                              stop 
                             endif
          y = log(k2/flamd2)
          yc = log(q02/flamd2)

          IF(  ID.EQ. 0)THEN
          FY1=Y -  YC  -11.0d0/12.0d0
          FY2=     2.0d0/3.0d0
          FY0= 6.0d0*FY1 + FY2*NF*2.0d0*TR
          FY1=FY1*6.0d0/FY0
         IF( FY1  .GE.  RANDOM(1) )THEN
                  ID1 = 0
                  ID2 = 0
                                   ELSE
                  ID1 = NF*RANDOM( 1 ) +1
                  ID2 =  -id1
C         WRITE(6,*)' NEW ID =',IDNEW
                                   ENDIF
                       ELSE
                  ID1=ID
                  ID2 =0
                       ENDIF

c     idnew : flavor of  child
c
c      determine x
c        FUNCTION  SPLIT(Y,YC,ID,IDNEW)
C
          IF( ID.EQ.0  )THEN
c     gluon branching
           if( id1 .eq. 0)then
c      gluon deacy into two gluons
           XLXMX =Y -YC
           DO 603 I = 1 , 100
             XLX = XLXMX*RANDOM(1)
                  X = EXP(-XLX)
          if( x .ge. 0.5d0)goto 630
        if(kzap .eq. 0)then
        AP =1.0d0+x/(1.0d0-x) -2.0d0*x+x*x*(1.d0-X)
              IF( AP .GE.   2.0d0*RANDOM(1) ) GOTO 630
                     else
        AP =1.0d0+x/(1.0d0-x) -2.0d0*x+x*x*(1.d0-X)
        AP =1.0d0/log(x*k2/flamd2)*ap
        apx =1.2d0/log(q02/flamd2)
              IF( AP .GE. apx*RANDOM(1) ) GOTO 630
                     endif
cxx              IF( AP .GE. 1.2d0*RANDOM(1) ) GOTO 630
  603    CONTINUE
  630   CONTINUE
                                           ELSE
c      gluon branches to quark-antiquark
           DO 604 I = 1 , 100
                  X = RANDOM(1)
        AP = X**2 +(1.0d0-X)**2
              IF( AP .GE.   RANDOM(1) ) GOTO 640
  604    CONTINUE
  640    CONTINUE
                                           ENDIF
                       ELSE
c
c       quark branching
           XLXMX =Y -YC
           DO 605 I = 1 , 100
             XLX = XLXMX*RANDOM(1)
                  X = EXP(-XLX)
cxx        AP = 1.0d0 + (1.0d0-x)**2
cxx              IF( AP .GE.   2*RANDOM(1) ) GOTO 650
        if(kzap .eq. 0)then
        AP =(1.0d0 + (1.0d0-X)**2 )*(1.d0-X)
              IF( AP .GE.   2.0d0*RANDOM(1) ) GOTO 650
                     else
        AP =(1.0d0 + (1.0d0-X)**2 )*(1.d0-X)
        AP =1.0d0/log(x*k2/flamd2)*ap
        apx =2.0d0/log(q02/flamd2)
              IF( AP .GE. apx*RANDOM(1) ) GOTO 650
                     endif
  605    CONTINUE
  650    CONTINUE
            x = 1.0d0 - x
                   ENDIF
C------------------------------------------------------------------

C         check the constraint is satisfied or not

cx        write(6,*)'x,k2,q2max in psgen2=',x,k2,q2max
          IF(  (1.0D0-X) .LT. K2/Q2MAX ) THEN
cx        write(6,*)'x<t,x,k2 in psgen2=',x,k2
                                         Q2Mx = K2
                                         GOTO 5
                                        ENDIF
c
c         if gluon branch, impose the angular ordering
            npr = nad - listp(2,nad)
            xpr = partn(1,nad) / partn(1,npr)
            qxx = partn(2,npr)
              th0 = qxx/xpr/(1.0d0-xpr) 
              th1 = k2/x/(1.0d0-x)/xpr/xpr
       if(     id.eq.0 
     & .and.   id1.eq.0
     & .and.   listp(6,npr) .ge. 0
     & .and.  (1.0d0-xpr) .ge. xpr*(1.0d0-x)
     & .and.   th0 .le. th1  )then
c        Write(6,*)'break the A.O.,x,k2,xpr,qxx,th0,th1 in psgen2=',
c     &   x,k2,xpr,qxx ,th0,th1
                                         Q2Mx = K2
                                         GOTO 5
           endif
c
c   if gluon branch to q-q pair, check kinematics

       if(   id.eq.0 
     & .and. id1.ne.0
     & .and. k2.le.qmass(abs(id1))**2/x/(1.0d0-x))then
                                         Q2Mx = K2
                                         GOTO 5
           endif
c
c   if quark branch , check kinematics

       if(   id.ne.0 
     & .and. k2.le.qmass(abs(id))**2/x)then
                                         Q2Mx = K2
                                         GOTO 5
           endif
C         write partn and listp of the children
C         write partn and listp of the children

          L1= NTOP+1
          partn(1,L1) = X*partn(1,NAD)

          listp(1,L1) = ID1
          listp(2,L1) = L1- NAD

          L2 = NTOP+2
          partn(1,L2) =(1.0D0-X)*partn(1,NAD)
          listp(1,L2) = id2
          listp(2,L2) = L2- NAD

          partn(2,nad) = k2
          listp(3,NAD) = 2
          listp(4,NAD) = L1 - NAD
          listp(5,NAD) = 1

          listp(5,L1) = 0
          listp(5,L2) = 0
          listp(6,L1) = 1
          listp(6,L2) = 1

           NTOP = NTOP + 2
c       CALL qclist
c       WRITE(6,111)
          RETURN

C         not branch
 101      CONTINUE
c        WRITE(6,*)' NOT BRANCH'

          if( id .eq.0)then
          partn(2,nad) = 0.0d0
                                else
          partn(2,nad) = qmass(abs(id))**2
                                endif
          listp(3,NAD) = 0
          listp(4,NAD) = 0
          listp(5,NAD) = 1
c       WRITE(6,111)
 111   FORMAT('   END OF PSGEN2')

          RETURN
          END


         SUBROUTINE LIMQ2(Q2MX,NAD)
         IMPLICIT REAL*8(A-H,O-Z)
        COMMON/qcevnt/partn(10,1000),listp(10,1000),NTOP
             common/color/cf,ca,tr,pi,beta0
             common/qcpar/roots,q2max,q02,q2min,flamd,qmass(10),nflav
          NPR = NAD - listp(2,NAD)
        IF( listp(6,NPR).eq.-1 )THEN
           IF( listp(3,NPR).EQ.1 )THEN
             Q2MX=Q2MAX
                                  RETURN
                                  ELSE
        L1 = listp(4,NPR) + NPR
              IF( L1.EQ.NAD)L1=L1+1
          X = partn(1,NAD)/partn(1,NPR)
          Q2MX = X*(partn(2,NPR) - partn(2,L1)/(1.0D0-X) )
                  RETURN
                                ENDIF
                              ENDIF
            Q2X = partn(2,NPR)
        L1 = listp(4,NPR) + NPR
              IF( L1.EQ.NAD)L1=L1+1
            X   = partn(1,NAD)/partn(1,NPR)
          qms = qmass(abs(listp(1,l1)))**2
          qms = max(qms, partn(2,l1))
c          Q2MX = X*(Q2X - partn(2,l1)/(1.0d0-x))
          Q2MX = X*(Q2X - qms/(1.0d0-x))
          return
          end
c          calculation of non-branching probability
c
             subroutine qcprep
c
             implicit real*8(a-h,o-z)
c
             common/color/cf,ca,tr,pi,beta0
             common/qcpar/roots,q2max,q02,q2min,flamd,qmass(10),nflav
             common/qcflg/kzap,nao
             common/tpnes/pnes(0:100,0:10),tmaxs,tmins,ndims
             common/tpnet/pnet(0:100,0:10),tmaxt,tmint,ndimt
      common/photon/alpha,rtx,rtph,rtvm,nfl,nfvm
      common/phovmu/vmuaa,vmua,vmub,vmuc,vxmu,vmunr
      common/phovmg/vmgaa,vmga,vmgb,vmgc,vxmg,vmgnr
      common/qcwght/wght0,fail,wght,accptx
      common/cross/sigma0,sigmax,sigma
      common/prcgg/q2xmax,q2xmin,zmin,ame2,ecmx0
      common/rslx/fph0,fph0x
      common/prcss/ecm,xbmin,nprc

   

c
             ndim = 100
             ndims = ndim
             ndimt = ndim
c
             flamd2 = flamd**2
c
                fac = 2.0d0/beta0
                ymax = log(q2max/flamd2)
                ymin = log(q2min/flamd2)
                ycut = log(q02/flamd2)
                yycut =  log(ycut)
                dy = ( ymax - ymin )/ndim
c
               tmaxs = ymax
               tmins = ymin
               tmaxt = ymax
               tmint = ymin

c
              do 101 k=0,ndim
              y = ymin + k*dy
c
              if( kzap .eq. 0)then
              fys1= y - (ycut+11.0d0/12.0d0)*log(y)
              fyt1= y - (ycut+11.0d0/12.0d0)*log(y)
                              else
              fys1= y*(log(y)-1.0d0) - (yycut+11.0d0/12.0d0)*log(y)
              fyt1= y*(log(y)-1.0d0) - (yycut+11.0d0/12.0d0)*log(y)
                              endif
              fys2 = 1.0d0/3.0d0*log(y)
              fyt2 =2.0d0/3.0d0*log(y)
              pnes(k,0) = fac*(ca*2.0d0*fys1+nflav*2.0d0*tr*fys2)
              pnet(k,0) = fac*(ca*2.0d0*fyt1+nflav*2.0d0*tr*fyt2)
 101         continue
              pnsx = pnes(ndim,0)
              pntx = pnet(ndim,0)
              do 102 k=0,ndim
              pnes(k,0)=pnes(k,0) - pnsx
              pnet(k,0)=pnet(k,0) - pntx
 102          continue
c
c
              do 110 nf=1,nflav
              do 111 k=0,ndim
              y = ymin + k*dy
c

              y = ymin + k*dy
c
              if(kzap .eq. 0)then
              fys1= 2.0d0*y - (2.0d0*ycut+17.0d0/6.0d0)*log(y)
              fyt1= 2.0d0*y - (2.0d0*ycut+4.0d0/3.0d0)*log(y)
                             else
      fys1= 2.0d0*y*(log(y)-1.0d0)-(2.0d0*yycut+17.0d0/6.0d0)*log(y)
      fyt1= 2.0d0*y*(log(y)-1.0d0)-(2.0d0*yycut+4.0d0/3.0d0)*log(y)
                             endif
              fys2 = 4.0d0/3.0d0*log(y)
              pnes(k,nf) = fac*cf*(fys1+fys2)
              pnet(k,nf) = fac*cf*(fyt1)
 111  continue
              pnsx = pnes(ndim,nf)
              pntx = pnet(ndim,nf)
              do 112 k=0,ndim
              pnes(k,nf)=pnes(k,nf) - pnsx
              pnet(k,nf)=pnet(k,nf) - pntx
 112          continue
 110          continue
c
          ntbl =0
            if(ntbl.eq.1)then
        write(6,*)' table of PNE, gluon,quark(space), gluon,quark(time)'
            do 115 i=0,ndim
            write(6,116)pnes(i,0),pnes(i,1),pnet(i,0),pnet(i,1) 
 115        continue
 116        format(' ',4(G12.4,2X))
                           endif
c
c     parameters of resolved photon and quark pair created photon
c 
        alpha = 1.0d0/137.0d0
        nfl = 4
        nfvm = 2
c       eeav : charge average
c      1/3   :integration of x(x**2+(1-x)**2)
c      3  : color factor
c      2  : kinematical one
c
cx          eeav=(4.0/9.0+1.0/9.0+1.0/9.0)
          eeav=(4.0d0+2.0d0)/9.0d0
cx        spph = alpha/2.0/pi*(1.0/3.0)*2.0*3.0*eeav
        spph = alpha/pi*eeav
cx        spphc = alpha/2.0/pi*(1.0/3.0)*2.0*3.0*(4.0/9.0)
        spphc = alpha/pi*(4.0d0/9.0d0)
c
c     paramaters of hadronic components
c     ( Gluck .., P.R.D46(1992)1973   )
        vmuaa =0.235d0
        vmua = 0.5d0
        vmub = 0.082d0
        vmuc = 0.354d0
        vxmu = 0.5d0
       vmunr =  0.131d0 * nfvm*2
c
        vmgaa =0.535d0
        vmga = 0.462d0
        vmgb = 0.354d0
        vmgc = 0.233d0
        vxmg =  1.0d0
       vmgnr =  0.390
c
           alam2 = flamd**2
           qq0 = q2min
           qq  = q2max
           YC=    LOG(Q02 /ALAM2)
           Y0=    LOG(qq0/ALAM2)
           YS=    LOG(QQ/ALAM2)
c        WRITE(6,*)'Q-EVOLUTION USING ALTARELLI-PARISI EQUATION '
c        WRITE(6,*)'    CUTOFF Q02  =',Q02,' qq0 ',q2min
c        WRITE(6,*)'  Q**2=',QQ
c        WRITE(6,*)'  photon structure function '
c
         if(qq.le.qmass(4)**2)write(6,*)' qq < charm mass**2 '
        fph0x = spph*log(qq/qq0)
        fph0 = spph*log(qq/qq0)+ spphc*log(qq/qmass(4)**2)
         write(6,*)' fph0x(uds),fph0(udsc)=',fph0x,fph0

c
         alpi = alpha/2.0d0/pi
c     hadronic component
        fvm = alpha* ( vmunr + vmgnr)

cx         fvm = 0.0
c         fph0 = 0.0
        fint = fph0 + fvm
c
           walph = alpha
             rtph = fph0/walph
             rtvm = fvm/walph
       rtx =  rtph + rtvm
         wght0 = rtx
c           write(6,*)'  rtx,rtph,rtvm = ',rtx,rtph,rtvm
             rtph =rtph/rtx
             rtvm =rtvm/rtx
cx           write(6,*)'  rtx,rtph,rtvm = ',rtx,rtph,rtvm
cx        write(6,*)' fph0/w,fph0(c)/w = ',fph0/walph,(fph0-fph0x)/walph
cc           write(6,*)' sm,fint,fint/walph = ',sm,fint,fint/walph
c           write(6,*)' fph0,fvm =',fph0,fvm
cx             write(6,*)' q2max,q2min in prep.f =',q2max,q2min
c
c       tentative total cross section
c       reference cross section (QED muon pair)
c       sigma0=4pi*alpha**2/s = 86.856/s(GeV^2) (nb)
c      4pi*alpha**2 int dQ^2/Q^4 Fgamma/e(z,E) int dx xf(x)
c
       if( nprc .eq. 2)then
        esum = alpha*rtx*4.0d0/9.0d0 
        sigma0 = 86.856d3
         eng = ecmx0/2.0d0
         ame = sqrt(ame2)
         theta = 0.1d0
           rm = 1.0d0/q2xmin - 1.0d0/q2xmax
cx          write(6,*)' zmin,eng,ame= ',zmin,eng,ame
        fsum=ft(zmin,eng,ame,theta)*alpha/pi
        sigmax = sigma0*rm*fsum*esum
c
ccccc       factor 2 : e+ large angle scattering  as well as e-
              sigmax = sigmax*2.0d0 
        write(6,*)'sigma0,rm,fsum,esum=',sigma0,rm,fsum,esum 
        write(6,*)' sigmax =',sigmax
                     endif

              return
              end

c
         real*8 function ft(z,eng,ame,theta)
         IMPLICIT REAL*8(A-H,O-Z)
c
          zl=log(z)
          z1 = 1.0d0 - z
c
      f1= 2.0d0*(-zl-z1+0.25d0*(1.0d0-z*z))*log(eng*theta/ame)
c
         f2=zl**2+2.0d0*(-1.0d0-z*(zl-1.0d0))
     &     -0.25d0*(-1.0d0-2.0d0*z**2*zl-1.0d0)
c
         ft = f1+f2
cx         write(6,*)' z, f1,f2,ft=',z,f1,f2,ft
         return
         end
         SUBROUTINE ANAL(MODE)       
        IMPLICIT REAL*8(A-H,O-Z)    
      COMMON/qcevnt/partn(10,1000),listp(10,1000),NTOP
      common/color/cf,ca,tr,pi,beta0
      common/qcpar/roots,q2max,q02,q2min,flamd,qmass(10),nflav
      common/qcwght/wght0,fail,wght,accptx
      common/prcss/ecm,xbmin,nprc
      common/photon/alpha,rtx,rtph,rtvm,nfl,nfvm
      common/cross/sigma0,sigmax,sigma
c
        COMMON/cxxxx/accept,ncount,naccpt,nbin,nawgt
        COMMON/axxxx/xbin(10,100),ybin(-4:4,10,100),dx(10)
        common/bxxxx/xsum(-4:4),thmin
        IF( MODE .EQ. 0 )THEN  
      NCOUNT  = 0             
      NACCPT  = 0            
      ACCPT  = 0.0
      nawgt = 0
C                
C               
C              
        NBIN = 50   
C                  
       DX(1) = 1.0/NBIN 
       DO 101 J = 1,NBIN       
        XBIN(1,J) =DX(1)*(J-1)
      do 102 nf=-3,3
 102     YBIN(nf,1,J) =0.0   
 101  CONTINUE             
         return
                        endif
c
          if( mode .eq. 1)then
c
      naccpt = naccpt + 1
                                
c

       do 201 i=1,ntop
       if(listp(6,i).ge.1 .and. listp(3,i).eq.0)then
        xb = 2.0d0*partn(7,i)/ecm
        nx = xb/dx(1) + 1
        if( nx .le.0)nx=1
        if( nx .ge. nbin)nx=nbin
        ybin(0,1,nx)=ybin(0,1,nx)+wght
                                              endif
 201    continue
cccccccccccccccccccccccccccccc
                         RETURN    
                         ENDIF    
C                                
                               
         esum =rtx
                              
         DO 300 J=1,NBIN   
         YBIN(0,1,J)=YBIN(0,1,J)/DX(1)/wght0
 300     CONTINUE            
c
                                                          
          return
          END                                                                   
C   16/11/87 801271332  MEMBER NAME  MKLIST@  (FORT)     M  FORTRAN
C
C ------------------------------------------------------------
C  MKLIST
C ------------------------------------------------------------
C  MAKE LIST OF FINAL PARTONS ORDERED AS IN PLANAR TREE
C ------------------------------------------------------------
C
      SUBROUTINE MKLIST( LEAVES , NPRTN )
C
      IMPLICIT REAL*8(A-H,O-Z)
      INTEGER*4 STACK,PNTR
C
      COMMON/QCEVNT/PARTN(10,1000),LISTP(10,1000),NATOP
      DIMENSION STACK(1000),LEAVES(1000)
      DIMENSION MC(3),LC(3)
C  STACK : WORK AREA, USED AS A STACK
C  PNTR  : POINTER OF STACK
C  LEAVES : ORDERED LIST, KEEP ABSOLUTE ADDRESS IN PARTN AND LISTP
C  NPRTN  : NUMBER OF PARTONS IN FINAL STATE,  I.E., SIZE OF LEAVES
C
      NPRTN=0
c
      do 2000 ixx=1,2
      PNTR=1
      STACK(PNTR)=ixx
      IX=1
       if( ixx.eq.2)ix=0
C
 1000 N = STACK(PNTR)
      NC = LISTP(3,N)
      IF(NC.EQ.0) GO TO 100
C
      NA=LISTP(4,N)
      DO 1 I=1,NC
        MC(I)=LISTP(1, N+NA+I-1)
    1 CONTINUE
      CALL MKLSUB(MC,NC,LC,IX)
      IX=0
       if( ixx.eq.2)ix=1
C
      PNTR=PNTR-1
      DO 10 I=1,NC
        PNTR=PNTR+1
        STACK(PNTR)=N+NA+(LC(I)-1)
   10 CONTINUE
      NN =STACK(PNTR)
      STACK(PNTR)=STACK(PNTR-NC+1)
      STACK(PNTR-NC+1)=NN
      GO TO 1000
C
  100 NPRTN=NPRTN+1
      LEAVES(NPRTN)=STACK(PNTR)
      PNTR=PNTR-1
      IF(PNTR.LT.1) GO TO 2000
      GO TO 1000
C
 2000    continue
 9000 continue
cx         write(6,*)' # of final particles =',nprtn
cx         do 30 i=1,nprtn
cx 30          write(6,*)' i=',i,'  address =',leaves(i)
        nstop = 0
        if(nstop .eq. 1)stop
       RETURN
      END
C
C
C ------------------------------------------------------------
C  MKLSUB
C ------------------------------------------------------------
C  DETERMINE BARNCH TYPE AND ORDERING OF CHILD-PARTONS
C ------------------------------------------------------------
      SUBROUTINE MKLSUB( MC , NC , LC , IX )
C
      IMPLICIT REAL*8(A-H,O-Z)
      DIMENSION MC(3),LC(3)
C
C ALLOWED ORDER :
C TWO PARTONS  ; Q+G   , G+QBAR   , G+G      , QBAR+Q
C THREE PARTONS; Q+G+G , G+G+QBAR , Q+QBAR+Q , QBAR+Q+QBAR ,
C                G+G+G , G+QBAR+Q
C PRIMARY VERTEX (WHEN IX=1); Q+QBAR , Q+G+QBAR
      DIMENSION MAP2(3,3),MAP3(3,3,3),MAPL(3,6)
      DATA MAP2/ 0 , 1 , 2 , 2 , 0 , 1 , 1 , 2 , 1 /
      DATA MAP3/ 0 , 2 , 0 , 1 , 3 , 1 , 0 , 2 , 4 ,
     *3 , 1 , 3 , 2 , 0 , 0 , 6 , 0 , 1 ,
     *0 , 5 , 2 , 4 , 0 , 3 , 1 , 4 , 1 /
      DATA MAPL/ 1 , 2 , 3 ,
     *2 , 1 , 3 ,
     *1 , 3 , 2 ,
     *3 , 2 , 1 ,
     *3 , 1 , 2 ,
     *2 , 3 , 1 /
C
      if( nc .eq. 1)then
         lc(1)=1
             return
                     endif
      IF(NC.ge.3 ) GO TO 900
C  REDEFINE MC; Q=1, QBAR=2, G=3
      DO 1 I=1,NC
        M= MC(I)
        IF(M.GT.0) MC(I)=1
        IF(M.LT.0) MC(I)=2
        IF(M.EQ.0) MC(I)=3
    1 CONTINUE
C
      IF(NC.EQ.2) THEN
        LF=MAP2(MC(1),MC(2))
      ELSE
        LF=MAP3(MC(1),MC(2),MC(3))
      END IF
      IF(LF.EQ.0) GO TO 900
C
      DO 2 I=1,NC
        LC(I)=MAPL(I,LF)
    2 CONTINUE
C
      IF(IX.EQ.1) THEN
        IF(NC.EQ.2) THEN
          LL=LC(1)
          LC(1)=LC(2)
          LC(2)=LL
        ELSE
          LL=LC(1)
          LC(1)=LC(3)
          LC(3)=LC(2)
          LC(2)=LL
        END IF
      END IF
C
      RETURN
C
  900 WRITE(6,*) ' TROUBLE IN MKLSUB, NC=',NC
      WRITE(6,*) '                    MC=',MC
      STOP
      END
       subroutine cms(nfail)
         IMPLICIT REAL*8(A-H,O-Z)
      common/qcevnt/partn(10,1000),listp(10,1000),NTOP
      common/color/cf,ca,tr,pi,beta0
      common/qcpar/roots,q2max,q02,q2min,flamd,qmass(10),nflav
      common/prcss/ecm,xbmin,nprc
      common/cms1/cos2,phi2x
      common/qcwght/wght0,fail,wght,fail2
      dimension itm(2,100),pppp(4,100),pref(4)
c
             nfail = 0
c
c          assign momentum of initial photon and electron
c
          l1 = 2
          if(listp(1,l1) .ne. 22)then
                 write(6,*)' error in photon  in cms'
                  call qclist
                                stop
                                 endif
         partn(2,L1) = 0.0d0
         partn(3,L1) =  0.0d0
         partn(4,L1) =  0.0d0
         partn(5,L1) =  0.0d0
         partn(6,l1) = ecm/2.0d0
         partn(7,l1) =  ecm/2.0d0
         partn(8,L1)=( partn(7,L1)+partn(6,L1) )/SQRT(2.0d0)
         partn(9,L1)=( partn(7,L1)-partn(6,L1) )/SQRT(2.0d0)

         listp(5,L1) = 2

c    cms of spacelike partons
c
          do 100 i=3,ntop
          if( listp(5,i).ne.2 .and. listp(6,i).eq.-1)then
          call cmsint(i)
                                                     endif
 100   continue
cx        call qclist
c
c
c        cms of scattered electron and parton 
          call cmssct(nfailx)
          if(nfailx .eq.1 )then
                  nfail = 1
                  return
                          endif
c
cx        call qclist
c      cms of timelike partons
c
          do 200 i=3,ntop
        if( listp(5,i).ne.2 .and. listp(6,i).eq.1)then
          call cmsfnl(i)
                                                     endif
 200   continue
c
c         call qclist
       nstop = 0
      if(nstop.eq. 1)return
         do 210 i=3,ntop
          npr = i-listp(2,i)
        if(listp(6,i).eq.1 .and.listp(3,npr).eq.1)then
             itop = i
             goto 211
                                                    endif
 210     continue
 211     continue
cx            write(6,*)'scatt. par.=',itop
         itm(1,1) = itop
         itm(2,1) = 1
          nit = 1
          ni = 1
 220     continue  
         if( ni .gt. nit)goto 222
         if(itm(2,ni).eq.1)then
           nd =itm(1,ni)
         if(listp(3,nd) .eq. 2)then
          l1 = nd+listp(4,nd)
          l2 = l1+1
         itm(1,nit+1) = l1
         itm(1,nit+2) = l2
         itm(2,nit+1) = 1
         itm(2,nit+2) = 1
            nit = nit + 2
                              endif
          itm(2,ni) = 2
                           endif
            ni=ni+1
              goto 220
c
 222        continue
cx          write(6,*)' address ',(itm(1,i),i=1,nit)
c          write(6,*)' state ',(itm(2,i),i=1,nit)
        phi1 = 0.0d0
        cost = cos2
        phi2 = phi2x
         nnn1=1
         nnn2 = nit
         ndim = 4
        do 231 i=1,nit
         i0 =itm(1,i)
         do 231 k=1,3
  231   pppp(k,i)=partn(3+k,i0)
c
cx       write(6,*)' nit =',niu 
cx       write(6,*)'cost ,phi1,phi2=',cost,phi1,phi2
cx       write(6,*)'parton bf rt',(pppp(k,1),k=1,3)
        call ROTAT( PPPP , PHI1, COST , PHI2 , nnn1,nnn2 , NDIM )
cx       write(6,*)'parton af rt',(pppp(k,1),k=1,3)
c
        do 232 i=1,nit
         i0 =itm(1,i)
         do 232 k=1,3
  232   partn(3+k,i0) = pppp(k,i)
            
c
c         boost to center mass system of electron and photon
          pref(1)=-(partn(4,1)+partn(4,2))
          pref(2)=-(partn(5,1)+partn(5,2))
          pref(3)=-(partn(6,1)+partn(6,2))
          pref(4)=partn(7,1)+partn(7,2)
         do 241 i=1,ntop
         do 241 k=1,4
  241   pppp(k,i)=partn(3+k,i)
         nnn1 = 1
         nnn2 = ntop
      CALL BOOST( PPPP , PREF , NNN1 , NNN2 )
c
       px = pppp(1,2)
       py = pppp(2,2)
       pz = pppp(3,2)
       pt2 =px**2+py**2
        cos2 = pz/sqrt(pz**2 + pt2 )
        if( px .eq. 0)then
           phi1 = pi/2.0d0
           if(py .le. 0.0d0)phi1=-phi1
                      else
           phi1 = atan(py/px)
              if(px.lt. 0.0d0)then
             phi1 = phi1 + pi
                              endif
                      endif 
          cost =  cos2
          phi1 = pi - phi1
          phi2 =  0.0d0
        call ROTAT( PPPP , PHI1, COST , PHI2 , nnn1,nnn2 , NDIM )

        do 242 i=1,ntop
         do 242 k=1,4
  242   partn(3+k,i) = pppp(k,i)
             return
             end

       SUBROUTINE  cmssct(nfail)
         IMPLICIT REAL*8(A-H,O-Z)
      common/qcevnt/partn(10,1000),listp(10,1000),NTOP
      common/color/cf,ca,tr,pi,beta0
      common/qcpar/roots,q2max,q02,q2min,flamd,qmass(10),nflav
      common/prcss/ecm,xbmin,nprc
      common/cms1/cos2,phi2x
      common/qcwght/wght0,fail,wght,fail2
      common/prcgg/q2xmax,q2xmin,zmin,ame2,ecmx0


c
         nfail = 0
c
cx        write(6,*)' enter to cmssct'
           do 100 i=3,ntop
         if( listp(3,i).eq.1  .and. listp(6,i).eq.-1)goto 101
 100      continue
           write(6,*)' error ; not find scatted parton in cms'
c
 101     continue
          l0 = 2
          lp1 = 1
          lp2 = i
          l1 = lp1 + listp(4,lp1)
          l2 = lp2 + listp(4,lp2)
c
         smp2 = partn(2,lp2)
         sm2  = partn(2,l2)
         q2= -q2max
         ptx= partn(4,lp2)
         pty= partn(5,lp2)
         pt2 = ptx**2+pty**2
         if(pt2 .eq. 0.0d0)then
            cs= 1.0d0
            sn = 0.0d0
                           else
            cs =ptx/sqrt(pt2)
            sn =pty/sqrt(pt2)
                           endif
         x0 = partn(1,lp2)
    
         dd =(1.0d0+(smp2-sm2)/q2)**2-4.0d0*(smp2+pt2)/q2
         if(dd .lt. 0.0d0)then
       write(6,*)' dd<0 at cmssct,sm2,smp2,pt2=',sm2,smp2,pt2
           call qclist
                 stop
                         endif
         y =0.5d0*(1.0d0+(smp2-sm2)/q2+sqrt(dd))
        if(abs(y-1.0d0).ge.1.5d0)then
               if( nprc .eq.2)then
cx              write(6,*)'ecm ,q2 =',ecm,q2max
                              endif
cx             write(6,*)'y >< 1,y=',y,' x0=',x0
cx       write(6,*)' sm2,smp2,pt2=',sm2,smp2,pt2
cx         call qclist
                                 endif
             xb = x0/y  
          if( nprc .eq. 2)then
        wght = 1.0d0/xb*(1.0d0-q2max/ecm**2/xb)
c        write(6,*)' wght in cmssct =',wght
                         else
           wght = 1.0d0
                         endif
c
cxx        s = roots**2
        s = ecm**2
        pkt2=xb*s*(xb*s+q2)/(-q2)
cx          write(6,*)' pkt2=',pkt2
         if(pkt2 .lt. 0.0d0)then
cxc           write(6,*)'pkt2<0,pkt2,s,xb=',pkt2,s,xb
c            call qclist
              nfail = 1
              return
                            endif
c
c     P+ of the initial photon 
c
       ppls = partn(8,l0)
c
c       momentu of virtual(transfer) photon
c
       qpls = - xb*ppls
       qmns = q2/2.0d0/qpls
c
c        write(6,*)' qpls,qmns at cmssct',qpls,qmns
c      momentum assign of initial and scattered electrons
c
        partn(9,lp1) = s/2.0d0/ppls 
        partn(8,lp1) = pkt2/partn(9,lp1)/2.0d0
        partn(9,l1)  = partn(9,lp1) - qmns
        partn(8,l1)  = partn(8,lp1) - qpls
c
        partn(7,lp1)  = (partn(8,lp1) + partn(9,lp1))/sqrt(2.0d0)
        partn(6,lp1)  = (partn(8,lp1) - partn(9,lp1))/sqrt(2.0d0)
        partn(7,l1)   = (partn(8,l1)  + partn(9,l1) )/sqrt(2.0d0)
        partn(6,l1)   = (partn(8,l1)  - partn(9,l1) )/sqrt(2.0d0)
        pkt = sqrt(pkt2)
        partn(4,lp1)  = -pkt*sn
        partn(5,lp1)  =  pkt*cs
        partn(4,l1)   = -pkt*sn
        partn(5,l1)   =  pkt*cs
        listp(5,lp1)=  2
        listp(5,l1) =  2

       
c
c       momentu assign of scattered parton(tentative)
c
        p2pls = partn(8,lp2)+qpls
        p2mns = partn(9,lp2)+qmns
        partn(7,l2)=(p2pls+p2mns)/sqrt(2.0d0)
        p2z=(p2pls-p2mns)/sqrt(2.0d0)
        p2x = pt2 + p2z**2
        partn(8,l2) = (partn(7,l2) + sqrt(p2x))/sqrt(2.0d0)
        partn(9,l2) = (partn(7,l2) - sqrt(p2x))/sqrt(2.0d0)
        partn(6,l2) = sqrt(p2x)
        partn(5,l2) = 0.0d0
        partn(4,l2) = 0.0d0
c
        listp(5,l2) =  2
c
c        write(6,*)'pt2,p2z,ptx,pty=',pt2,p2z,ptx,pty
        cos2 = p2z/sqrt(p2z**2 + pt2 )
        if( ptx .eq. 0)then
           phi2 = pi/2.0d0
           if(pty .le. 0.0d0)phi2=-phi2
                      else
           phi2 = atan(pty/ptx)
              if(ptx.lt. 0.0d0)then
             phi2 = phi2 + pi
                              endif
                      endif 
          phi2x = phi2
        
          return
          end

       SUBROUTINE CMSINT(NAD)
         IMPLICIT REAL*8(A-H,O-Z)
      common/qcevnt/partn(10,1000),listp(10,1000),NTOP
      common/color/cf,ca,tr,pi,beta0
      common/qcpar/roots,q2max,q02,q2min,flamd,qmass(10),nflav


        NPR = NAD - listp(2,NAD)
        L1 = NAD
        L2 = L1 + 1
        X= partn(1,L1) / partn(1,NPR)
        partn(8,L1)= partn(8,NPR)*X
        partn(8,L2)= partn(8,NPR)*(1.0D0-X)
        PT2=X*(1.0D0-X)
     &  *(partn(2,NPR)-partn(2,L1)/X-partn(2,L2)/(1.0D0-X))
        IF( PT2 .LT. 0.0)THEN
              WRITE(6,*)' ERROR AT CMSINT, PT2<0, PT2=',PT2
              CALL qclist
                         ENDIF

         THETA = 2.0d0*PI*RANDOM(1)
         PT = SQRT( PT2 )
         PX = PT*COS(THETA)
         PY = PT*SIN(THETA)

         partn(3,L1) = PT2
         partn(3,L2) = PT2
         partn(4,L1) = partn(4,NPR)*X         + PX
         partn(4,L2) = partn(4,NPR)*(1.0D0-X) - PX
         partn(5,L1) = partn(5,NPR)*X         + PY
         partn(5,L2) = partn(5,NPR)*(1.0D0-X) - PY
         partn(9,L1)=(partn(2,L1)
     &   +partn(4,L1)**2+partn(5,L1)**2)/partn(8,L1)/2.0d0
         partn(9,L2)=(partn(2,L2)
     &   +partn(4,L2)**2+partn(5,L2)**2)/partn(8,L2)/2.0d0
         partn(6,L1)=( partn(8,L1)-partn(9,L1) )/SQRT(2.0d0)
         partn(6,L2)=( partn(8,L2)-partn(9,L2) )/SQRT(2.0d0)

         partn(7,L1)=( partn(8,L1)+partn(9,L1) )/SQRT(2.0d0)
         partn(7,L2)=( partn(8,L2)+partn(9,L2) )/SQRT(2.0d0)

         listp(5,L1) = 2
         listp(5,L2) = 2

          RETURN
          END


          SUBROUTINE CMSFNL(NAD)
         IMPLICIT REAL*8(A-H,O-Z)
      common/qcevnt/partn(10,1000),listp(10,1000),NTOP
      common/color/cf,ca,tr,pi,beta0
      common/qcpar/roots,q2max,q02,q2min,flamd,qmass(10),nflav
C
C
        NPR = NAD - listp(2,NAD)
        L1 = NAD
        L2 = L1 + 1
        X= partn(1,L1) / partn(1,npr)
      IF( X .LE. 0.0d0  .OR.  X .GE. 1.0d0)THEN
        WRITE(6,*)'npr,L1,L2,X at cmsfnl =',npr,L1,L2,X
          CALL qclist
          stop
                                       ENDIF
        partn(8,L1)= partn(8,npr)*X
        partn(8,L2)= partn(8,npr)*(1.0D0-X)
        PT2 = X*(1.0D0-X)
     &  *(partn(2,npr)-partn(2,L1)/X-partn(2,L2)/(1.0D0-X))
        IF( PT2 .LT. 0.0)THEN
              WRITE(6,*)' ERROR AT CMSFNL, PT2<0, PT2=',PT2
              WRITE(6,*)' npr,L1,L2,X =',npr,L1,L2,X
              CALL qclist
                         ENDIF

         THETA = 2.0d0*PI*RANDOM(1)
         PT = SQRT( PT2 )
         PX = PT*COS(THETA)
         PY = PT*SIN(THETA)

         partn(3,L1) = PT2
         partn(3,L2) = PT2
         partn(4,L1) = partn(4,npr)*X         + PX
         partn(4,L2) = partn(4,npr)*(1.0D0-X) - PX
         partn(5,L1) = partn(5,npr)*X         + PY
         partn(5,L2) = partn(5,npr)*(1.0D0-X) - PY
         partn(9,L1)=(partn(2,L1)
     &    +partn(4,L1)**2+partn(5,L1)**2)/partn(8,L1)/2.0d0
         partn(9,L2)=(partn(2,L2)
     &    +partn(4,L2)**2+partn(5,L2)**2)/partn(8,L2)/2.0d0
         partn(6,L1)=( partn(8,L1)-partn(9,L1) )/SQRT(2.0d0)
         partn(6,L2)=( partn(8,L2)-partn(9,L2) )/SQRT(2.0d0)

         partn(7,L1)=( partn(8,L1)+partn(9,L1) )/SQRT(2.0d0)
         partn(7,L2)=( partn(8,L2)+partn(9,L2) )/SQRT(2.0d0)

C
C        WRITE(6,*)'npr,P(I)',npr,(partn(I,npr),I=1,9)
C        WRITE(6,*)'L1 ,P(I)',L1 ,(partn(I,L1 ),I=1,9)

C        WRITE(6,*)'L2 ,P(I)',L2 ,(partn(I,L2 ),I=1,9)
C
         listp(5,L1) = 2
         listp(5,L2) = 2

          RETURN
          END
 
C**********
C MODULE NO.=034
C.....................MAINTENANCE RECORD AREA........................
C   VERSION 84/06/08
C   UPDATE  00/00/00 : PERSON( XXXXXXX )
C
C....................................................................
C
      SUBROUTINE BOOST( PPPP , PREF , NNN1 , NNN2 )
C
      IMPLICIT REAL*8(A-H,O-Z)
C
      DIMENSION PPPP(4,1), PREF(4)
CCC   DATA EPS/1.0E-5/
C
C     BOOST A SET OF 4-VECTOR PPPP(4,I) I=NNN1 TO NNN2
C     BY 4-VECTOR PREF(4)
C
C     BE CAREFUL FOR SIZE OF ARRAY PPPP
C
      SUM=0.0d0
      DO 10 K=1,3
   10 SUM = SUM+PREF(K)**2
      P = SQRT(SUM)
      E0 = PREF(4)
      IF(E0.LE.0.0) THEN
          WRITE(6,601) PREF
  601     FORMAT(' #EPOCS3401(W,BOOST) ENERGY<=0. PREF=',4G15.5)
          E0=1.0d-30
          ENDIF
      BET = P/E0
C                           NO NEED OF BOOST
CCC   IF(BET.LE.EPS) RETURN
      IF(BET.GE.1.0d0) THEN
          WRITE(6,602) PREF
  602     FORMAT(' #EPOCS3402(W,BOOST) VELOCITY>=1. PREF=',4G15.5)
          BET=0.99999d0
          ENDIF
      BE2 =    SQRT((1.0d0-BET)*(1.0d0+BET))
      GAM = 1.0d0/BE2
      GA1 = 1.0d0/(BE2*(1.0d0+BE2))
C
      DO 100 N = NNN1 , NNN2
        S = 0.0d0
        DO 20 K = 1,3
        S = S + PREF(K)*PPPP(K,N)
   20   CONTINUE
        E = PPPP(4,N)
CCC     PX = (GAM-1.0d0)*S/P**2+GAM*E/E0
        PX =  GA1    *S/E0**2+GAM*E/E0
        DO 30 K = 1,3
        PPPP(K,N) = PPPP(K,N) + PX*PREF(K)
   30   CONTINUE
        PPPP(4,N) = GAM*( E + S/E0 )
  100 CONTINUE
C
      RETURN
      END
C**********
C* ROTAT *
C**********
C MODULE NO.=039
C.....................MAINTENANCE RECORD AREA........................
C   VERSION 84/06/08
C   UPDATE  00/00/00 : PERSON( XXXXXXX )
C
C....................................................................
C
        SUBROUTINE ROTAT( PPPP , PHI1, COST , PHI2 , nnn1,nnn2 , NDIM )
C        EULER ROTATION OF PPPP(I,J)
C            I=1 TO NDIM, (I=VECTOR INDEX , NDIM= 3 OR 4 )
Cxxxxx            J=1 TO NUMP, (J=NO. OF PARTICLE)
C            J=nnn1,nnnn2, (J=NO. OF PARTICLE)
C
C         ........................
          IMPLICIT REAL*8(A-H,O-Z)
C         ........................
       DIMENSION PPPP( NDIM , nnn2)
C
       COP=COS(PHI1)
       SIP=SIN(PHI1)
       COT=COST
       SIT=SQRT(1.0d0-COST**2)
       COR=COS(PHI2)
       SIR=SIN(PHI2)
C (1) ROTATE PHI1 AROUND Z-AXIS
      DO 10 J=nnn1,nnn2
        PX=PPPP(1,J)
        PY=PPPP(2,J)
        PPPP(1,J)=PX*COP-PY*SIP
        PPPP(2,J)=PX*SIP+PY*COP
   10 CONTINUE
C (2) ROTATE THETA AROUND Y-AXIS
      DO 20 J=nnn1,nnn2
        PX=PPPP(1,J)
        PZ=PPPP(3,J)
        PPPP(3,J)=PZ*COT-PX*SIT
        PPPP(1,J)=PZ*SIT+PX*COT
   20 CONTINUE
C (3) ROTATE PHI2 AROUND Z-AXIS
      DO 30 J=nnn1,nnn2
        PX=PPPP(1,J)
        PY=PPPP(2,J)
       PPPP(1,J)=PX*COR-PY*SIR
        PPPP(2,J)=PX*SIR+PY*COR
   30 CONTINUE
      RETURN
      END
