EPOS/epos167/epos-rsh-165.f

   1 c  reshuffled from sem, sto, sha
   2
   3 c    contains psahot and related stuff
   4 c             ------
   5
   6
   7
   8 c-----------------------------------------------------------------------
   9       subroutine psahot(kcol,ncolp,iret)
  10 c-----------------------------------------------------------------------
  11 c psahot - showering (semihard parton-parton interaction)
  12 c-----------------------------------------------------------------------
  13       include 'epos.inc'
  14       include 'epos.incsem'
  15       include 'epos.incems'
  16       double precision ept(4),ept1(4),xx,wpt(2),eprt,pl,plprt,wplc
  17      *,wp0,wm0,s,si,smin,xmin,xp1,wpi,wmi,xp,xm,wp1,wm1,wp2,wm2
  18      *,wpq,wmq,wpi1,wmi1,pxh,pyh,pth,xmm,xp2,xg,zg,smax,xmax
  19      *,xmax1,xmin1,zp0,psutz,xpmax0,xpmin0,gb0,tmax0,tmin0,zpm,gb
  20      *,gbyj,tmax,tmin,t,gb7,x,s2,discr,qt,qt2,x1min,x1max,t1min,t1max
  21      *,t1,xq1,qq,qqmax,qqmin,pt2,xmin2,ptnew,xpmin,xpmax,xm0,psuds
  22 c     *,xprh,xmrh
  23       dimension ep3(4),ey(3),bx(6)
  24      *,qmin(2),iqc(2),nqc(2),ncc(2,2),amqt(4)
  25       parameter (mjstr=20000)
  26       common /psar7/  delx,alam3p,gam3p
  27       common /psar29/eqj(4,mjstr),iqj(mjstr),ncj(2,mjstr),nj
  28       common /psar30/iorj(mjstr),ityj(mjstr),bxj(6,mjstr)
  29       common /testj/  ajeth(4),ajete(5),ajet0(7)
  30       parameter (ntim=1000)
  31       common/cprt/nprtj,pprt(5,ntim),idprt(ntim),iorprt(ntim)
  32      &,idaprt(2,ntim)
  33       common/emsptl/nppr(npommx,kollmx),npproj(mamx),nptarg(mamx)
  34       integer icp(2),ict(2),icini(2)
  35       integer jcp(nflav,2),jct(nflav,2),nemis(2)
  36       common/cprtx/nprtjx,pprtx(5,2)/ciptl/iptl
  37
  38       call utpri('psahot',ish,ishini,3)
  39
  40       iret=0
  41       alpq=-(alppar+1.)/2.
  42       qqcut=q2min  !????????????pt2cut
  43
  44
  45       nptl1=nptl
  46       iptl=nppr(ncolp,kcol)
  47       ip=iproj(kcol)
  48       it=itarg(kcol)
  49       do i=1,2
  50         icp(i)=icproj(i,ip)
  51         ict(i)=ictarg(i,it)
  52       enddo
  53       idpomr=idhpr(ncolp,kcol)
  54       bpomr=bhpr(ncolp,kcol)
  55
  56       q2finsave=q2fin
  57       zzzz=zzremn(ip,1)+zzremn(it,2)
  58       zz=min( 1+zoemax , 1+zoeinc  *zzzz ) !<-----
  59       q2fin=q2fin*zz
  60
  61       ajeth(idpomr+1)=ajeth(idpomr+1)+1.
  62 c      write(*,*)ajeth
  63       idfpomr=idfpr(ncolp,kcol)
  64       if(ish.ge.3)write(ifch,*)'Start psahot (icp,ict):',ip,icp,it,ict
  65      *,ncolp,kcol,iptl,idpomr,idfpomr,bpomr
  66
  67
  68       if(idfpomr.eq.0)stop'idfpomr??????'
  69       if(ish.ge.3)then
  70         write(ifch,20)iptl
  71      *  ,sqrt(pptl(1,iptl)**2+pptl(2,iptl)**2),pptl(3,iptl)
  72      *  ,pptl(4,iptl),pptl(5,iptl)
  73      *  ,istptl(iptl),ityptl(iptl)
  74 20      format(1x,i4,3x,4(e11.5,1x),i2,1x,i3)
  75       endif
  76       istptl(iptl)=31
  77
  78 1     nj=0
  79       nptl=nptl1
  80
  81       wp0=dsqrt(xpr(ncolp,kcol))*dexp(ypr(ncolp,kcol))*dble(engy)     !???? new
  82       wm0=dsqrt(xpr(ncolp,kcol))*dexp(-ypr(ncolp,kcol))*dble(engy)    !double
  83
  84
  85       amqt(1)=sqrt(sngl(xxp1pr(ncolp,kcol)**2+xyp1pr(ncolp,kcol)**2))
  86       amqt(2)=sqrt(sngl(xxp2pr(ncolp,kcol)**2+xyp2pr(ncolp,kcol)**2))
  87       amqt(3)=sqrt(sngl(xxm2pr(ncolp,kcol)**2+xym2pr(ncolp,kcol)**2))
  88       amqt(4)=sqrt(sngl(xxm1pr(ncolp,kcol)**2+xym1pr(ncolp,kcol)**2))
  89       amqpt=amqt(1)+amqt(2)+amqt(3)+amqt(4)
  90
  91       s2min=4.*q2min
  92       if(sngl(wp0*wm0).le.(sqrt(s2min)+amqpt)**2)then
  93         if(ish.ge.1)then
  94           call utmsg('psahot: insufficient pomeron mass&')
  95           write (ifch,*)'mass:',dsqrt(wp0*wm0),amqpt+sqrt(s2min)
  96           call utmsgf
  97         endif
  98         iret=1
  99         goto 16
 100       endif
 101
 102       ih=iproj(kcol)
 103       jh=itarg(kcol)
 104 c      xprh=xpp(ih)
 105 c      xmrh=xmt(jh)
 106       rp=r2had(iclpro)+r2had(icltar)+slopom*log(engy**2)
 107       z=exp(-bpomr**2/(4.*.0389*rp))
 108
 109       if(z.eq.0)then
 110        write(ifch,*)'psahot : z,ih,jh ! -> ',z,ih,jh
 111        call gakli2(ih,ih)
 112        call gakli2(jh,jh)
 113        call gakli2(iptl,iptl)
 114        stop
 115       endif
 116
 117       do l=1,4
 118         bx(l)=xorptl(l,iptl)
 119       enddo
 120       bx(5)=tivptl(1,iptl)
 121       bx(6)=tivptl(2,iptl)
 122       ity=ityptl(iptl)
 123
 124       if(idpomr.eq.0)then        !gg-pomeron
 125         iqq=0   !type of the hard interaction: 0 - gg, 1 - qg, 2 - gq, 3 - qq
 126         pxh=0.d0  !p_x for sh pomeron
 127         pyh=0.d0  !p_y for sh pomeron
 128       elseif(idpomr.eq.1)then    !qg-pomeron
 129         iqq=1
 130         pxh=xxp1pr(ncolp,kcol)
 131         pyh=xyp1pr(ncolp,kcol)
 132         amqpt=amqpt-amqt(1)
 133       elseif(idpomr.eq.2)then    !gq-pomeron
 134         iqq=2
 135         pxh=xxm2pr(ncolp,kcol)
 136         pyh=xym2pr(ncolp,kcol)
 137         amqpt=amqpt-amqt(3)
 138       elseif(idpomr.eq.3)then    !qq-pomeron
 139         iqq=3
 140         pxh=xxp1pr(ncolp,kcol)+xxm2pr(ncolp,kcol)
 141         pyh=xyp1pr(ncolp,kcol)+xym2pr(ncolp,kcol)
 142         amqpt=amqpt-amqt(1)-amqt(3)
 143       else
 144         stop'unknown pomeron'
 145       endif
 146       pth=pxh**2+pyh**2
 147
 148       nj0=nj
 149       if(ish.ge.6)write(ifch,*)'iptl,nptl,wp0,wm0,z,iqq,bx:'
 150       if(ish.ge.6)write(ifch,*) iptl,nptl,wp0,wm0,z,iqq,bx
 151
 152       s=wp0*wm0                !lc+*lc- for the semihard interaction
 153       smin=dble(s2min)+pth          !mass cutoff for the hard pomeron
 154       xmin=smin/s
 155       smax=(dsqrt(s)-dble(amqpt))**2  !max mass for the hard pomeron
 156       xmax=smax/s
 157       if(smax.le.smin)then
 158         write (ifmt,*)'smax,smin',smax,smin
 159         iret=1
 160         goto 16
 161       endif
 162       xpmax0=psutz(s,smin,dble(amqpt**2))   !max x+ for the hard pomeron
 163       xpmin0=xmin/xpmax0              !min x+ for the hard pomeron
 164       xp1=wp0/dble(engy)              !lc+ share for the semihard interaction
 165       xp2=wm0/dble(engy)              !lc- share for the semihard interaction
 166
 167
 168 c-----------------------------------------------------------------
 169 c  determine LC momenta wpi,wmi for hard Pomeron
 170 c-----------------------------------------------------------------
 171
 172       if(ish.ge.4)write(ifch,*)
 173      & 'determine LC momenta wpi,wmi for hard Pomeron'
 174
 175       !------------------------------------------
 176       if(iqq.eq.3)then     !     val-val
 177       !------------------------------------------
 178
 179        if(ish.ge.4)write(ifch,*)'val-val'
 180         xmin1=xmin**dble(delh+.4)
 181         xmax1=xmax**dble(delh+.4)
 182         zp0=dsqrt(xmin)
 183         if(zp0.ge.1.d0)call utstop('zp0 in sem&')
 184       !........ kinematical bounds
 185         tmax0=dlog((1.d0+dsqrt(1.d0-zp0))/(1.d0-dsqrt(1.d0-zp0)))
 186         tmin0=dlog((1.d0+dsqrt(1.d0-xpmax0))
 187      *                   /(1.d0-dsqrt(1.d0-xpmax0)))
 188         if(iclpro.ne.4)then
 189           call psjti0(sngl(smax-pth),sqq,sqqb,1,1)
 190           call psjti0(sngl(smax-pth),sqqp,sqqpb,1,2)
 191           call psjti0(sngl(smax-pth),sqaq,sqaqb,-1,1)
 192         else
 193           call psjti0(sngl(smax-pth),sqqp,sqqpb,4,2)
 194           sqq=0.
 195           sqaq=0.
 196         endif
 197         uv1=psdfh4(sngl(zp0*xp1),q2min,0.,iclpro,1)
 198         dv1=psdfh4(sngl(zp0*xp1),q2min,0.,iclpro,2)
 199         uv2=psdfh4(sngl(zp0*xp2),q2min,0.,icltar,1)
 200         dv2=psdfh4(sngl(zp0*xp2),q2min,0.,icltar,2)
 201         wwuu=uv1*uv2*sqq
 202         if(iclpro.eq.2)then
 203           wwdd=dv1*dv2*sqq
 204         elseif(iclpro.eq.1)then
 205           wwdd=dv1*dv2*sqaq
 206         elseif(iclpro.eq.3)then
 207           wwdd=dv1*dv2*sqqp
 208         elseif(iclpro.eq.4)then
 209           wwuu=uv1*uv2*sqqp
 210           wwdd=0.
 211         endif
 212         wwud=uv1*dv2*sqqp
 213         wwdu=dv1*uv2*sqqp
 214         wudt=wwuu+wwdd+wwud+wwdu
 215         gb0=dble(wudt)/xmax**dble(delh)/xmin**0.4*
 216      *    (1.d0-zp0*xp1)**dble(-1.-alpq-alplea(iclpro))*
 217      *    (1.d0-zp0*xp2)**dble(-1.-alpq-alplea(icltar))*(tmax0-tmin0)
 218      *    *(1.d0-zp0)**dble(.5+alpq)*(1.d0-zp0)**dble(alpq)  *5.d0
 219
 220 3       zpm=(xmin1+dble(rangen())*(xmax1-xmin1))**dble(1./(delh+.4))   !zpm
 221         if(iclpro.ne.4)then
 222           call psjti0(sngl(zpm*s-pth),sqq,sqqb,1,1)
 223           call psjti0(sngl(zpm*s-pth),sqqp,sqqpb,1,2)
 224           call psjti0(sngl(zpm*s-pth),sqaq,sqaqb,-1,1)
 225         else
 226           call psjti0(sngl(zpm*s-pth),sqqp,sqqpb,4,2)
 227           sqq=0.
 228           sqaq=0.
 229         endif
 230         xpmax=psutz(s,zpm*s,dble(amqpt**2))  !max x+ for sh pomeron
 231         tmax=dlog((1.d0+dsqrt(1.d0-dsqrt(zpm)))
 232      *             /(1.d0-dsqrt(1.d0-dsqrt(zpm))))
 233         tmin=dlog((1.d0+dsqrt(1.d0-xpmax))/(1.d0-dsqrt(1.d0-xpmax)))
 234
 235         t=(tmin+dble(rangen())*(tmax-tmin))
 236         xp=1.d0-((1.d0-dexp(-t))/(1.d0+dexp(-t)))**2  !x+_v
 237         xm=zpm/xp                             !x-_v
 238         if(xm.gt.xp.and.ish.ge.1)write(ifmt,*)'xm,xp',xm,xp
 239         gb=(1.d0-xm)**alpq*(1.d0-xp)**(.5+alpq)*(tmax-tmin)
 240         if(rangen().lt..5)then
 241           xp=xm
 242           xm=zpm/xp
 243         endif
 244
 245         uv1=psdfh4(sngl(xp*xp1),q2min,0.,iclpro,1)
 246         dv1=psdfh4(sngl(xp*xp1),q2min,0.,iclpro,2)
 247         uv2=psdfh4(sngl(xm*xp2),q2min,0.,icltar,1)
 248         dv2=psdfh4(sngl(xm*xp2),q2min,0.,icltar,2)
 249         wwuu=uv1*uv2*sqq
 250         if(iclpro.eq.2)then
 251           wwdd=dv1*dv2*sqq
 252         elseif(iclpro.eq.1)then
 253           wwdd=dv1*dv2*sqaq
 254         elseif(iclpro.eq.3)then
 255           wwdd=dv1*dv2*sqqp
 256         elseif(iclpro.eq.4)then
 257           wwuu=uv1*uv2*sqqp
 258           wwdd=0.
 259         endif
 260         wwud=uv1*dv2*sqqp
 261         wwdu=dv1*uv2*sqqp
 262         wudt=wwuu+wwdd+wwud+wwdu
 263
 264         gb=gb*dble(wudt)/zpm**dble(delh+0.4)
 265      *    *(1.d0-xp*xp1)**dble(-1.-alpq-alplea(iclpro))
 266      *    *(1.d0-xm*xp2)**dble(-1.-alpq-alplea(icltar))/gb0
 267 c          if(ish.ge.4)then
 268           if(gb.gt.1.d0.and.ish.ge.1)write (ifch,*)
 269      *      'gb-qq,iclpro,zpm,xp,tmax,tmin,xpmax',
 270      *      gb,iclpro,zpm,xp,tmax,tmin,xpmax
 271 c          endif
 272         if(dble(rangen()).gt.gb)goto 3
 273
 274         aks=rangen()*wudt
 275         if(aks.lt.wwuu)then
 276           if(iclpro.ne.4)then
 277             iq1=1
 278           else
 279             iq1=4*(2.*int(.5+rangen())-1.)
 280           endif
 281           iq2=1
 282         elseif(aks.lt.wwuu+wwdd)then
 283           if(iclpro.eq.2)then
 284             iq1=2
 285           elseif(iclpro.eq.1)then
 286             iq1=-2
 287           elseif(iclpro.eq.3)then
 288             iq1=-3
 289           else
 290             stop'?????'
 291           endif
 292           iq2=2
 293         elseif(aks.lt.wwuu+wwdd+wwud)then
 294           if(iclpro.ne.4)then
 295             iq1=1
 296           else
 297             iq1=4*(2.*int(.5+rangen())-1.)
 298           endif
 299           iq2=2
 300         else
 301           if(iclpro.eq.2)then
 302             iq1=2
 303           elseif(iclpro.eq.1)then
 304             iq1=-2
 305           elseif(iclpro.eq.3)then
 306             iq1=-3
 307           else
 308             stop'?????'
 309           endif
 310           iq2=1
 311         endif
 312
 313         wpi=xp*wp0       !lc+ for the semihard interaction
 314         wmi=xm*wm0       !lc- for the semihard interaction
 315         wp1=(wp0-wpi)
 316         wm1=(wm0-wmi)
 317         wp1=wp1*psutz(wp1*wm1,dble(amqt(2)**2),dble(amqt(4)**2))
 318         wm1=wm1-amqt(2)**2/wp1
 319
 320       !-------------------------------------
 321       else  ! sea-sea  val-sea  sea-val
 322       !-------------------------------------
 323
 324        if(ish.ge.4)write(ifch,*)'sea-sea  val-sea  sea-val'
 325         xmin1=xmin**(delh-dels)
 326         xmax1=xmax**(delh-dels)
 327
 328         if(iqq.eq.0)then    !rejection function normalization
 329           gb0=dlog(xpmax0/xpmin0)*(1.d0-dsqrt(xmin))**(2.*betpom) !y_soft =
 330         else
 331           tmax0=acos(dsqrt(xpmin0))  !kinematical limits for t=cos(x+-)**2
 332           tmin0=acos(dsqrt(xpmax0))
 333           if(iqq.eq.1)then
 334             uv1=psdfh4(sngl(xpmin0*xp1),q2min,0.,iclpro,1)
 335      *        *sngl(1.d0-xpmin0*xp1)**(-1.-alpq-alplea(iclpro))
 336             dv1=psdfh4(sngl(xpmin0*xp1),q2min,0.,iclpro,2)
 337      *        *sngl(1.d0-xpmin0*xp1)**(-1.-alpq-alplea(iclpro))
 338           else
 339             uv1=psdfh4(sngl(xpmin0*xp2),q2min,0.,icltar,1)
 340      *        *sngl(1.d0-xpmin0*xp2)**(-1.-alpq-alplea(icltar))
 341             dv1=psdfh4(sngl(xpmin0*xp2),q2min,0.,icltar,2)
 342      *        *sngl(1.d0-xpmin0*xp2)**(-1.-alpq-alplea(icltar))
 343           endif
 344           gb0=(1.d0-xmin/xpmax0)**betpom*dble(uv1+dv1)
 345      *      *xpmin0**(-0.5+dels)
 346      *      *(1.d0-xpmin0)**(0.5+alpq)*(tmax0-tmin0)
 347           if(ish.ge.6)write (ifch,*)
 348      *      'gb0,tmax0,tmin0,xpmax0,xpmin0,xmin,xp1,xp2',
 349      *      gb0,tmax0,tmin0,xpmax0,xpmin0,xmin,xp1,xp2
 350         endif
 351         if(iclpro.ne.4.or.iqq.ne.1)then
 352           call psjti0(sngl(smax-pth),sj,sjb,iqq,0) !inclusive (sj) and born
 353         else
 354           call psjti0(sngl(smax-pth),sj,sjb,4,0)
 355         endif
 356         gb0=gb0*dble(sj)/xmax**delh      *1.5d0    !rejection function norm.
 357         if(ish.ge.6)write (ifch,*)'gb0,sj,z',gb0,sj,z
 358
 359         if(gb0.le.0.)then
 360           write (ifmt,*)'gb0<0, smax,pth',smax,pth
 361           iret=1
 362           goto 16
 363         endif
 364
 365         ! sharing of light cone momenta between soft preevolution and
 366         ! hard interaction: ( first energy-momentum is shared according to
 367         ! f_hard(yj)~zpm**(delh-dels-1) and then rejected as
 368         ! w_rej ~sigma_hard_tot(yj) / exp(delh*yj)
 369
 370  4        continue
 371         zpm=(xmin1+dble(rangen())*(xmax1-xmin1))**dble(1./(delh-dels)) !zpm
 372         if(iclpro.ne.4.or.iqq.ne.1)then
 373           call psjti0(sngl(zpm*s-pth),sgg,sggb,0,0)!inclusive (sgg) and born
 374           call psjti0(sngl(zpm*s-pth),sgq,sgqb,0,1)
 375           call psjti0(sngl(zpm*s-pth),sqq,sqqb,1,1)
 376           call psjti0(sngl(zpm*s-pth),sqaq,sqaqb,-1,1)
 377           call psjti0(sngl(zpm*s-pth),sqqp,sqqpb,1,2)
 378           sqq=(sqq+sqaq+2.*(naflav-1)*sqqp)/naflav/2.
 379         else
 380           call psjti0(sngl(zpm*s-pth),sgq,sgqb,4,0)
 381           call psjti0(sngl(zpm*s-pth),sqq,sqqb,4,1)
 382         endif
 383         xpmax=psutz(s,zpm*s,dble(amqpt**2))  !max x+ for sh pomeron
 384         xpmin=zpm/xpmax
 385         if(ish.ge.8)write (ifch,*)'zpm,xpmax,xpmin',zpm,xpmax,xpmin
 386
 387         if(iqq.eq.0)then
 388           xp=xpmin*(xpmax/xpmin)**rangen()  !lc+ share for the hard interaction
 389           xm=zpm/xp           !lc- share for the hard interaction
 390         else
 391           tmax=acos(dsqrt(xpmin))  !kinematical limits for t=cos(x+-)**2
 392           tmin=acos(dsqrt(xpmax))
 393           t=tmin+dble(rangen())*(tmax-tmin)
 394           xp=cos(t)**2
 395           xm=zpm/xp
 396         endif
 397
 398         if(ish.ge.8)write(ifch,*)'zpm,xp,xm,xpmax,xpmin:',
 399      *      zpm,xp,xm,xpmax,xpmin
 400
 401         if(iqq.eq.0)then  ! --------- sea-sea -----------
 402
 403           dwwgg1=0.
 404           dwwgq1=0.
 405           dwwqg1=0.
 406           dwwqq1=0.
 407
 408           dwwgg2=0.
 409           dwwgq2=0.
 410           dwwqg2=0.
 411           dwwqq2=0.
 412
 413           wwgg=sngl((1.d0-xp)*(1.d0-xm))**betpom
 414           wwgq=sngl(1.d0-xp)**betpom*EsoftQZero(sngl(xm))
 415           wwqg=sngl(1.d0-xm)**betpom*EsoftQZero(sngl(xp))
 416           wwqq=EsoftQZero(sngl(xp))*EsoftQZero(sngl(xm))
 417
 418           if(idfpomr.eq.1)then
 419             rh=r2had(iclpro)+r2had(icltar)-slopom*sngl(dlog(zpm))
 420             wwgg=(wwgg*z**(rp/rh)/rh+dwwgg1+dwwgg2)
 421      *        *(r2had(iclpro)+r2had(icltar))/z
 422             wwgq=(wwgq*z**(rp/rh)/rh+dwwgq1+dwwgq2)
 423      *        *(r2had(iclpro)+r2had(icltar))/z
 424             wwqg=(wwqg*z**(rp/rh)/rh+dwwqg1+dwwqg2)
 425      *        *(r2had(iclpro)+r2had(icltar))/z
 426             wwqq=(wwqq*z**(rp/rh)/rh+dwwqq1+dwwqq2)
 427      *        *(r2had(iclpro)+r2had(icltar))/z
 428           elseif(idfpomr.eq.2)then
 429             rh=r2had(iclpro)+alam3p/2.-slopom*sngl(dlog(zpm))
 430             wwgg=wwgg/rh*(r2had(iclpro)+alam3p/2.)*z**(rp/rh-1.)
 431             wwqg=wwqg/rh*(r2had(iclpro)+alam3p/2.)*z**(rp/rh-1.)
 432             wwgq=wwgq/rh*(r2had(iclpro)+alam3p/2.)*z**(rp/rh-1.)
 433             wwqq=wwqq/rh*(r2had(iclpro)+alam3p/2.)*z**(rp/rh-1.)
 434           elseif(idfpomr.eq.3)then
 435             rh=r2had(icltar)+alam3p/2.-slopom*sngl(dlog(zpm))
 436             wwgg=wwgg/rh*(r2had(icltar)+alam3p/2.)*z**(rp/rh-1.)
 437             wwqg=wwqg/rh*(r2had(icltar)+alam3p/2.)*z**(rp/rh-1.)
 438             wwgq=wwgq/rh*(r2had(icltar)+alam3p/2.)*z**(rp/rh-1.)
 439             wwqq=wwqq/rh*(r2had(icltar)+alam3p/2.)*z**(rp/rh-1.)
 440           else
 441             stop'psahot-idfpomr????'
 442           endif
 443
 444           wwgg=wwgg*sgg*(1.-glusea)**2
 445           wwgq=wwgq*sgq*(1.-glusea)*glusea
 446           wwqg=wwqg*sgq*(1.-glusea)*glusea
 447           wwqq=wwqq*sqq*glusea**2
 448           gbyj=dlog(xpmax/xpmin)*dble(wwgg+wwgq+wwqg+wwqq)
 449           wpi=wp0*xp               !lc+ for the hard interaction
 450           wmi=wm0*xm               !lc+-for the hard interaction
 451           gbyj=gbyj/zpm**dble(delh)/gb0     !rejection fu
 452           if(gbyj.gt.1.d0.and.ish.ge.2)write(ifmt,*)'gbyj',gbyj
 453           if(dble(rangen()).gt.gbyj)goto 4                     !rejection
 454           wp1=wp0-wpi
 455           wm1=wm0-wmi
 456           call pslcsh(wp1,wm1,wp2,wm2,amqt,dble(amqpt))
 457
 458         else  ! --------- val-sea  sea-val -----------
 459
 460           dwwg=0.
 461           dwwq=0.
 462
 463           wwgq=sngl(1.d0-xm)**betpom
 464           wwqq=EsoftQZero(sngl(xm))
 465           if(idfpomr.eq.1)then
 466             rh=r2had(iclpro)+r2had(icltar)-slopom*sngl(dlog(xm))
 467             wwgq=(wwgq*z**(rp/rh)/rh+dwwg)
 468      *        *(r2had(iclpro)+r2had(icltar))/z
 469             wwqq=(wwqq*z**(rp/rh)/rh+dwwq)
 470      *        *(r2had(iclpro)+r2had(icltar))/z
 471           else
 472             if(iqq.eq.1)then
 473               rh=r2had(iclpro)+alam3p/2.-slopom*sngl(dlog(xm))
 474               wwgq=wwgq/rh*(r2had(iclpro)+alam3p/2.)*z**(rp/rh-1.)
 475               wwqq=wwqq/rh*(r2had(iclpro)+alam3p/2.)*z**(rp/rh-1.)
 476             else
 477               rh=r2had(icltar)+alam3p/2.-slopom*sngl(dlog(xm))
 478               wwgq=wwgq/rh*(r2had(icltar)+alam3p/2.)*z**(rp/rh-1.)
 479               wwqq=wwqq/rh*(r2had(icltar)+alam3p/2.)*z**(rp/rh-1.)
 480             endif
 481           endif
 482
 483           wwgq=wwgq*sgq*(1.-glusea)*sngl(xp)**(-0.5+dels)
 484      *      *sngl(1.d0-xp)**(0.5+alpq)
 485           wwqq=wwqq*sqq*glusea*sngl(xp)**(-0.5+dels)
 486      *      *sngl(1.d0-xp)**(0.5+alpq)
 487
 488           if(iqq.eq.1)then         !valence quark-gluon hard interaction
 489             uv1=psdfh4(sngl(xp*xp1),q2min,0.,iclpro,1)
 490             dv1=psdfh4(sngl(xp*xp1),q2min,0.,iclpro,2)
 491             wpi=wp0*xp             !lc+ for the hard interaction
 492             wmi=wm0*xm             !lc+-for the hard interaction
 493             aks=rangen()
 494             if(aks.lt.uv1/(uv1+dv1))then
 495               if(iclpro.ne.4)then
 496                 iq1=1
 497               else
 498                 iq1=4*(2.*int(.5+rangen())-1.)
 499               endif
 500             else
 501               iq1=2
 502             endif
 503             gbyj=dble((wwgq+wwqq)*(uv1+dv1))*(1.d0-xp*xp1)**
 504      *        (-1.-alpq-alplea(iclpro))
 505           else                !gluon-valence quark hard interaction
 506             xm=xp
 507             xp=zpm/xm
 508             uv1=psdfh4(sngl(xm*xp2),q2min,0.,icltar,1)
 509             dv1=psdfh4(sngl(xm*xp2),q2min,0.,icltar,2)
 510             wpi=wp0*xp             !lc+ for the hard interaction
 511             wmi=wm0*xm             !lc+-for the hard interaction
 512             aks=rangen()
 513             if(aks.lt.uv1/(uv1+dv1))then
 514               iq2=1
 515             else
 516               iq2=2
 517             endif
 518             gbyj=dble(wwgq+wwqq)*dble(uv1+dv1)*
 519      *        (1.d0-xm*xp2)**(-1.-alpq-alplea(icltar))
 520           endif
 521
 522           gbyj=gbyj*(tmax-tmin)/zpm**delh /gb0 /2.1d0 !rejection
 523           if(ish.ge.6)write (ifch,*)
 524      *      'gbyj,zpm,xp,tmax,tmin,xpmax,xpmin',
 525      *      gbyj,zpm,xp,tmax,tmin,xpmax,xpmin
 526           if(dble(rangen()).gt.gbyj)goto 4
 527
 528           wp1=wp0-wpi
 529           wm1=wm0-wmi
 530           if(ish.ge.8)write (ifch,*)'q_sea mass check',wp1*wm1,amqpt
 531           if(iqq.eq.1)then         !valence quark-gluon hard interaction
 532             amq1=amqt(3)**2
 533             s24=(amqt(2)+amqt(4))**2
 534           else
 535             amq1=amqt(1)**2
 536             s24=(amqt(2)+amqt(4))**2
 537             xp=xm
 538             xm=zpm/xp
 539           endif
 540           x1max=psutz(wp1*wm1,dble(amq1),dble(s24))
 541           x1min=dble(amq1)/x1max/wp1/wm1
 542           t1min=(1.d0/x1max-1.d0)
 543           t1max=(1.d0/x1min-1.d0)
 544 5           t1=t1min*(t1max/t1min)**dble(rangen())
 545           if(ish.ge.8)write (ifch,*)'t1,t1min,t1max',t1,t1min,t1max
 546           xq1=1.d0/(1.d0+t1)
 547           if(dble(rangen()).gt.(xq1*(1.d0-xq1))**(1.+(-alpqua)))goto 5
 548           if(iqq.eq.1)then         !valence quark-gluon hard interacti
 549             wm2=wm1*(1.d0-xq1)
 550             wm1=wm1*xq1
 551             wp1=wp1-dble(amq1)/wm1
 552             if(ish.ge.8)write (ifch,*)'q_sea+ mass check',
 553      *        wp1*wm2,s24
 554             wp1=wp1*psutz(wp1*wm2,dble(amqt(2)**2),dble(amqt(4)**2))
 555             wm2=wm2-dble(amqt(2))**2/wp1
 556           else
 557             wp2=wp1*(1.d0-xq1)
 558             wp1=wp1*xq1
 559             wm1=wm1-dble(amq1)/wp1
 560             if(ish.ge.8)write (ifch,*)'q_sea- mass check',
 561      *        wp2*wm1,s24
 562             wm1=wm1*psutz(wm1*wp2,dble(amqt(4)**2),dble(amqt(2)**2))
 563             wp2=wp2-amqt(4)**2/wm1
 564           endif
 565
 566         endif  ! -------------------
 567
 568       !-------------------------------
 569       endif
 570       !-------------------------------
 571
 572 c-------------------------------------------------------------------------
 573 c  flavor and momenta of end partons of the hard Pomeron
 574 c-------------------------------------------------------------------------
 575
 576        if(ish.ge.4)write(ifch,*)
 577      &  'flavor and momenta of end partons of the hard Pomeron'
 578
 579   6   continue
 580       wpi1=wpi
 581       wmi1=wmi
 582       nj=nj0          !initialization for the number of jets
 583       if(ish.ge.8)write (ifch,*)'5-ww,smin',wpi*wmi,smin
 584
 585       rrr=rangen()
 586       jqq=0
 587
 588       call iddeco(icp,jcp)
 589       call iddeco(ict,jct)
 590
 591       iret1=0
 592       iret2=0
 593
 594       if(iqq.eq.1.or.iqq.eq.2)then
 595         if(rrr.lt.wwqq/(wwgq+wwqq))jqq=1
 596       elseif(iqq.eq.0)then
 597         if(rrr.lt.wwqg/(wwgg+wwqg+wwgq+wwqq))then
 598           jqq=1
 599         elseif(rrr.lt.(wwqg+wwgq)/(wwgg+wwqg+wwgq+wwqq))then
 600           jqq=2
 601         elseif(rrr.lt.(wwqg+wwgq+wwqq)/(wwgg+wwqg+wwgq+wwqq))then
 602           jqq=3
 603         endif
 604       endif
 605
 606       if((iqq-1)*(iqq-3).eq.0)then
 607
 608         if(iabs(idptl(ip)).eq.1220)iq1=3-iq1      !proj=neutron
 609         iqc(1)=iq1                                !proj=particle
 610         if(idptl(ip).lt.-1000)iqc(1)=-iq1         !proj=antiparticle
 611         nj=nj+1
 612         ifl1=iabs(iq1)
 613
 614         if(iqc(1).gt.0)then
 615           call idsufl(ifl1,1,jcp,iret1)   !remove valence quark from remnant
 616           if(iret1.ne.1)ifl=idrafl(iclpro,jcp,2,'s',iret2)     !take sea antiquark
 617         elseif(iqc(1).lt.0)then
 618           call idsufl(ifl1,2,jcp,iret1)
 619           if(iret1.ne.1)ifl=idrafl(iclpro,jcp,1,'s',iret2)
 620         else
 621           call utstop('No quark for hard Pomeron+ in psahot!&')
 622         endif
 623
 624         if(iret1.eq.1.or.iret2.eq.1)then
 625           ifl=ifl1
 626          if(ish.ge.3)write(ifmt,*)'Not enough place in rem (psahot 1)'
 627          if(ish.ge.5)write(ifch,*)'Not enough place in rem (psahot 1)'
 628           call iddeco(icp,jcp)
 629           iret1=0
 630           iret2=0
 631         endif
 632
 633         if(ish.ge.5)write(ifch,*)'flavor vq+,sqb+',iq1,-ifl
 634
 635         if(ifl.eq.3)then
 636           iqj(nj)=-isign(ifl,iqc(1))*4/3
 637         elseif(ifl.eq.4)then
 638           iqj(nj)=-isign(ifl,iqc(1))/4
 639         else
 640           iqj(nj)=-isign(ifl,iqc(1))
 641         endif
 642         eqj(1,nj)=.5*sngl(wp1+dble(amqt(2))**2/wp1)
 643         eqj(2,nj)=wp1-eqj(1,nj)
 644         eqj(3,nj)=xxp2pr(ncolp,kcol)
 645         eqj(4,nj)=xyp2pr(ncolp,kcol)
 646         if(ish.ge.8)write (ifch,*)'q_v+ mass check',(eqj(1,nj)-
 647      *    eqj(2,nj))*(eqj(1,nj)+eqj(2,nj))-eqj(3,nj)**2-eqj(4,nj)**2
 648         eqj(1,nj)=sqrt(eqj(2,nj)**2+eqj(3,nj)**2+eqj(4,nj)**2)
 649         ncc(1,1)=nj
 650         ncc(2,1)=0
 651
 652       else
 653
 654         nj=nj+1
 655         if(idfpomr.ne.3)then
 656           iflq=idrafl(iclpro,jcp,1,'s',iret1)       !Choose flavor of sea quark.
 657           if(iret1.ne.1)iflqb=idrafl(iclpro,jcp,2,'s',iret2)  !Choose antiquark.
 658         else
 659           iflq=jdrafl(iclpro,jcp,1,iret1)       !Choose flavor of sea quark.
 660           iflqb=iflq                   !antiquark=quark (vertex end)
 661         endif
 662
 663         if(iret1.eq.1.or.iret2.eq.1)then
 664           iflqb=iflq
 665          if(ish.ge.3)write(ifmt,*)'Not enough place in rem (psahot 2)'
 666          if(ish.ge.5)write(ifch,*)'Not enough place in rem (psahot 2)'
 667           call iddeco(icp,jcp)
 668           iret1=0
 669           iret2=0
 670         endif
 671
 672         if(ish.ge.5)write(ifch,*)'flavor sq+,sqb+',iflq,-iflqb
 673
 674         if(iflqb.eq.3)then
 675           iqj(nj)=-iflqb*4/3
 676         elseif(iflqb.eq.4)then
 677           iqj(nj)=-iflqb/4
 678         else
 679           iqj(nj)=-iflqb
 680         endif
 681         if(iflq.eq.3)then
 682           iqj(nj+1)=iflq*4/3
 683         elseif(iflq.eq.4)then
 684           iqj(nj+1)=iflq/4
 685         else
 686           iqj(nj+1)=iflq
 687         endif
 688         eqj(1,nj)=.5*sngl(wp1+dble(amqt(1))**2/wp1)
 689         eqj(2,nj)=wp1-eqj(1,nj)
 690         eqj(3,nj)=xxp1pr(ncolp,kcol)
 691         eqj(4,nj)=xyp1pr(ncolp,kcol)
 692         if(ish.ge.8)write (ifch,*)'q_s1+ mass check',(eqj(1,nj)-
 693      *    eqj(2,nj))*(eqj(1,nj)+eqj(2,nj))-eqj(3,nj)**2-eqj(4,nj)**2
 694         eqj(1,nj)=sqrt(eqj(2,nj)**2+eqj(3,nj)**2+eqj(4,nj)**2)
 695         eqj(1,nj+1)=.5*sngl(wp2+dble(amqt(2))**2/wp2)
 696         eqj(2,nj+1)=wp2-eqj(1,nj+1)
 697         eqj(3,nj+1)=xxp2pr(ncolp,kcol)
 698         eqj(4,nj+1)=xyp2pr(ncolp,kcol)
 699         nj=nj+1
 700         if(ish.ge.8)write (ifch,*)'q_s2+ mass check',(eqj(1,nj)-
 701      *    eqj(2,nj))*(eqj(1,nj)+eqj(2,nj))-eqj(3,nj)**2-eqj(4,nj)**2
 702         eqj(1,nj)=sqrt(eqj(2,nj)**2+eqj(3,nj)**2+eqj(4,nj)**2)
 703
 704         if(jqq.eq.0.or.iqq.eq.0.and.jqq.eq.2)then
 705           iqc(1)=0
 706           ncc(1,1)=nj-1
 707           ncc(2,1)=nj
 708         else
 709           iqc(1)=int(3*rangen()+1.)*(2.*int(.5+rangen())-1.)
 710 7           zg=1.d0-dble(rangen())*(1.d0-xp)
 711           if(ish.ge.8)write (ifch,*)'6-zg,xp',zg,xp
 712           if(dble(rangen()).gt.zg**dels*((1.d0-xp/zg)/(1.d0-xp))
 713      *                                              **betpom)goto 7
 714           xg=xp/zg
 715           wpq=wp0*(xg-xp)
 716           wmq=1.d0/wpq
 717           wmi1=wmi1-wmq
 718           if(wmi1*wpi1.le.smin)goto 6
 719
 720           nj=nj+1
 721           if(iabs(iqc(1)).eq.3)then
 722             iqj(nj)=-iqc(1)*4/3
 723           else
 724             iqj(nj)=-iqc(1)
 725           endif
 726           eqj(1,nj)=.5*wmq
 727           eqj(2,nj)=-eqj(1,nj)
 728           eqj(3,nj)=0.
 729           eqj(4,nj)=0.
 730           if(ish.ge.8)write (ifch,*)'q_s3+ mass check',eqj(1,nj)**2-
 731      *      eqj(2,nj)**2-eqj(3,nj)**2-eqj(4,nj)**2
 732           if(iqc(1).gt.0)then
 733             ncj(1,nj)=nj-1
 734             ncj(1,nj-1)=nj
 735             ncj(2,nj)=0
 736             ncj(2,nj-1)=0
 737             ncc(1,1)=nj-2
 738             ncc(2,1)=0
 739           else
 740             ncj(1,nj)=nj-2
 741             ncj(1,nj-2)=nj
 742             ncj(2,nj)=0
 743             ncj(2,nj-2)=0
 744             ncc(1,1)=nj-1
 745             ncc(2,1)=0
 746           endif
 747         endif
 748       endif
 749
 750       if((iqq-2)*(iqq-3).eq.0)then
 751         if(iabs(idptl(maproj+it)).eq.1220)iq2=3-iq2     !targ=neutron
 752         iqc(2)=iq2             !tar=particle (can not be pion or kaon !)
 753         if(idptl(maproj+it).lt.-1000)iqc(2)=-iq2  !targ=antiparticle
 754         ifl2=iabs(iq2)
 755
 756         nj=nj+1
 757         if(iqc(2).gt.0)then
 758           call idsufl(ifl2,1,jct,iret1)   !remove valence quark from remnant
 759           if(iret1.ne.1)ifl=idrafl(icltar,jct,2,'s',iret2)    !take sea antiquark
 760         elseif(iqc(2).lt.0)then
 761           call idsufl(ifl2,2,jct,iret1)
 762           if(iret1.ne.1)ifl=idrafl(icltar,jct,1,'s',iret2)
 763         else
 764           call utstop('No quark for hard Pomeron- in psahot!&')
 765         endif
 766
 767         if(iret1.eq.1.or.iret2.eq.1)then
 768           ifl=ifl2
 769          if(ish.ge.3)write(ifmt,*)'Not enough place in rem (psahot 3)'
 770          if(ish.ge.5)write(ifch,*)'Not enough place in rem (psahot 3)'
 771           call iddeco(ict,jct)
 772           iret1=0
 773           iret2=0
 774         endif
 775
 776         if(ish.ge.5)write(ifch,*)'flavor vq-,sqb-',iq2,-ifl
 777
 778         if(ifl.eq.3)then
 779           iqj(nj)=-isign(ifl,iqc(2))*4/3
 780         elseif(ifl.eq.4)then
 781           iqj(nj)=-isign(ifl,iqc(2))/4
 782         else
 783           iqj(nj)=-isign(ifl,iqc(2))
 784         endif
 785
 786         eqj(1,nj)=.5*sngl(wm1+dble(amqt(4))**2/wm1)
 787         eqj(2,nj)=eqj(1,nj)-sngl(wm1)
 788         eqj(3,nj)=xxm1pr(ncolp,kcol)
 789         eqj(4,nj)=xym1pr(ncolp,kcol)
 790         if(ish.ge.8)write (ifch,*)'q_v- mass check',(eqj(1,nj)
 791      *    +eqj(2,nj))*(eqj(1,nj)-eqj(2,nj))-eqj(3,nj)**2-eqj(4,nj)**2
 792         eqj(1,nj)=sqrt(eqj(2,nj)**2+eqj(3,nj)**2+eqj(4,nj)**2)
 793         ncc(1,2)=nj
 794         ncc(2,2)=0
 795       else
 796         nj=nj+1
 797         if(idfpomr.ne.2)then
 798           iflq=idrafl(icltar,jct,1,'s',iret1)       !Choose flavor of sea quark.
 799           if(iret1.ne.1)iflqb=idrafl(icltar,jct,2,'s',iret2) !Choose  antiquark.
 800         else
 801           iflq=jdrafl(icltar,jct,1,iret1)       !Choose flavor of sea quark.
 802           iflqb=iflq                   !antiquark=quark (vertex end)
 803         endif
 804
 805         if(iret1.eq.1.or.iret2.eq.1)then
 806           iflqb=iflq
 807          if(ish.ge.3)write(ifmt,*)'Not enough place in rem (psahot 4)'
 808          if(ish.ge.5)write(ifch,*)'Not enough place in rem (psahot 4)'
 809           call iddeco(ict,jct)
 810           iret1=0
 811           iret2=0
 812         endif
 813
 814         if(ish.ge.5)write(ifch,*)'flavor sq-,sqb-',iflq,-iflqb
 815
 816         if(iflqb.eq.3)then
 817           iqj(nj)=-iflqb*4/3
 818         elseif(iflqb.eq.4)then
 819           iqj(nj)=-iflqb/4
 820         else
 821           iqj(nj)=-iflqb
 822         endif
 823         if(iflq.eq.3)then
 824           iqj(nj+1)=iflq*4/3
 825         elseif(iflq.eq.4)then
 826           iqj(nj+1)=iflq/4
 827         else
 828           iqj(nj+1)=iflq
 829         endif
 830         eqj(1,nj)=.5*sngl(wm1+dble(amqt(3))**2/wm1)
 831         eqj(2,nj)=eqj(1,nj)-sngl(wm1)
 832         eqj(3,nj)=xxm2pr(ncolp,kcol)
 833         eqj(4,nj)=xym2pr(ncolp,kcol)
 834         if(ish.ge.8)write (ifch,*)'q_s1- mass check',(eqj(1,nj)-
 835      *    eqj(2,nj))*(eqj(1,nj)+eqj(2,nj))-eqj(3,nj)**2-eqj(4,nj)**2
 836         eqj(1,nj)=sqrt(eqj(2,nj)**2+eqj(3,nj)**2+eqj(4,nj)**2)
 837         eqj(1,nj+1)=.5*sngl(wm2+dble(amqt(4))**2/wm2)
 838         eqj(2,nj+1)=eqj(1,nj+1)-wm2
 839         eqj(3,nj+1)=xxm1pr(ncolp,kcol)
 840         eqj(4,nj+1)=xym1pr(ncolp,kcol)
 841         nj=nj+1
 842         if(ish.ge.8)write (ifch,*)'q_s2- mass check',(eqj(1,nj)-
 843      *    eqj(2,nj))*(eqj(1,nj)+eqj(2,nj))-eqj(3,nj)**2-eqj(4,nj)**2
 844         eqj(1,nj)=sqrt(eqj(2,nj)**2+eqj(3,nj)**2+eqj(4,nj)**2)
 845
 846         if(jqq.eq.0.or.iqq.eq.0.and.jqq.eq.1)then
 847           iqc(2)=0
 848           ncc(1,2)=nj-1
 849           ncc(2,2)=nj
 850         else
 851           iqc(2)=int(3*rangen()+1.)*(2.*int(.5+rangen())-1.)
 852 8           zg=1.d0-dble(rangen())*(1.d0-xm)
 853           if(ish.ge.8)write (ifch,*)'7-zg,xm',zg,xm
 854           if(rangen().gt.zg**dels*((1.d0-xm/zg)/(1.d0-xm))**betpom)
 855      *      goto 8
 856           xg=xm/zg
 857           wmq=wm0*(xg-xm)
 858           wpq=1.d0/wmq
 859           wpi1=wpi1-wpq
 860           if(wmi1*wpi1.le.smin)goto 6
 861
 862           nj=nj+1
 863           if(iabs(iqc(2)).eq.3)then
 864             iqj(nj)=-iqc(2)*4/3
 865           else
 866             iqj(nj)=-iqc(2)
 867           endif
 868           eqj(1,nj)=.5*sngl(wpq)
 869           eqj(2,nj)=eqj(1,nj)
 870           eqj(3,nj)=0.
 871           eqj(4,nj)=0.
 872           if(ish.ge.8)write (ifch,*)'q_s3- mass check',(eqj(1,nj)-
 873      *      eqj(2,nj))*(eqj(1,nj)+eqj(2,nj))-eqj(3,nj)**2-eqj(4,nj)**2
 874           if(iqc(2).gt.0)then
 875             ncj(1,nj)=nj-1
 876             ncj(1,nj-1)=nj
 877             ncj(2,nj)=0
 878             ncj(2,nj-1)=0
 879             ncc(1,2)=nj-2
 880             ncc(2,2)=0
 881           else
 882             ncj(1,nj)=nj-2
 883             ncj(1,nj-2)=nj
 884             ncj(2,nj)=0
 885             ncj(2,nj-2)=0
 886             ncc(1,2)=nj-1
 887             ncc(2,2)=0
 888           endif
 889         endif
 890       endif
 891       if(jqq.ne.0)then
 892         if(iqq.ne.0.or.iqq.eq.0.and.jqq.eq.3)then
 893           if(iclpro.ne.4.or.iqq.ne.1)then
 894             call psjti0(sngl(wpi1*wmi1-pth),sqq1,sqqb1,1,1)
 895             call psjti0(sngl(wpi1*wmi1-pth),sqaq1,sqaqb1,-1,1)
 896             call psjti0(sngl(wpi1*wmi1-pth),sqqp1,sqqpb1,1,2)
 897             sqq1=(sqq1+sqaq1+2.*(naflav-1)*sqqp1)/naflav/2.
 898           else
 899             call psjti0(sngl(wpi1*wmi1-pth),sqq1,sqqb1,4,1)
 900           endif
 901           gbs=sqq1/sqq
 902         else
 903           call psjti0(sngl(wpi1*wmi1-pth),sgq1,sgqb1,0,1)
 904           gbs=sgq1/sgq
 905         endif
 906         if(ish.ge.8)write (ifch,*)'gbs',gbs
 907         if(rangen().gt.gbs)goto 6
 908       endif
 909
 910
 911 c---------------------------------------------------------------
 912 c        inner partons of the hard Pomeron
 913 c---------------------------------------------------------------
 914
 915        if(ish.ge.4)write(ifch,*)
 916      &  'inner partons of the hard Pomeron'
 917       nj00=nj
 918       wpi=wpi1
 919       wmi=wmi1
 920       si=wpi*wmi-pxh**2-pyh**2     !total energy squared for the hard
 921       if(ish.ge.7)write (ifch,*)'si,wpi,wmi',si,wpi,wmi
 922
 923       ept(1)=.5d0*(wpi+wmi)
 924       ept(2)=.5d0*(wpi-wmi)
 925       ept(3)=pxh
 926       ept(4)=pyh
 927       qmin(1)=q2min          !effective momentum cutoff above current la
 928       qmin(2)=q2min          !effective momentum cutoff below current la
 929       qminn=max(qmin(1),qmin(2)) !effective momentum cutoff for the bo
 930 c        si=psnorm(ept)       !4-momentum squared for the hard pomeron
 931       jfirst=1
 932       jj=int(1.5+rangen())
 933       nemis(1)=0
 934       nemis(2)=0
 935
 936   9   continue ! <<<<----------- ladder rung ---------------------------
 937
 938       if(ish.ge.4)write(ifch,*)'ladder rung'
 939       pt2=ept(3)**2+ept(4)**2
 940       pt=sqrt(sngl(pt2))
 941       if(iabs(iqc(1)).ne.4)then
 942         q2mass=0.
 943       else
 944         q2mass=qcmass**2
 945         si=si-dble(q2mass)
 946       endif
 947       s2min=4.*qminn+q2mass         !mass cutoff for born scattering
 948       wwmin=5.*qminn+q2mass-2.*pt*sqrt(q2ini)
 949       wwmin=(wwmin+sqrt(wwmin**2+4.*(q2mass+pt2)*(qminn-q2ini)))
 950      */(1.-q2ini/qminn)/2.
 951       if(ish.ge.5)write(ifch,*)'qminn,q2mass,pt,wwmin:'
 952       if(ish.ge.5)write(ifch,*)qminn,q2mass,pt,wwmin
 953
 954       wpt(1)=ept(1)+ept(2)            !lc+ for the current jet emissi
 955       wpt(2)=ept(1)-ept(2)            !lc- for the current jet emissi
 956
 957       if(iabs(iqc(1)).ne.4)then
 958         if(jfirst.eq.1)then
 959           sj=psjti(qmin(jj),qqcut,sngl(si),iqc(jj),iqc(3-jj),0)
 960           sj2=psjti1(qmin(3-jj),qmin(jj),qqcut
 961      *              ,sngl(si),iqc(3-jj),iqc(jj),0)
 962           if(ish.ge.5)write(ifch,*)'si,sj,sj2:',si,sj,sj2
 963           if(rangen().gt.sj2/sj.and.si.gt.1.1d0*dble(wwmin))goto 112
 964           jfirst=0
 965           jj=3-jj
 966           sj=sj2
 967           goto 111
 968         endif
 969        sj=psjti1(qmin(jj),qmin(3-jj),qqcut,sngl(si),iqc(jj),iqc(3-jj),0)
 970 111     sjb=psbint(qmin(1),qmin(2),qqcut,sngl(si),iqc(1),iqc(2),0) !born parton-parton
 971       else
 972         sjord=psjci(qmin(2),sngl(si),iqc(2))
 973         sj=sjord
 974       sjb=psbint(qmin(1),qmin(2),qqcut,sngl(si)+q2mass,iqc(1),iqc(2),0)
 975         if(qmin(2).eq.q2min)then
 976           wwmins=2.5*q2min*(1.+sqrt(1.+
 977      *    4.*(pt2+q2mass)/q2min))
 978        if(si.gt.dble(wwmins))call psjti0(sngl(si)+q2mass,sj,sjb
 979      *                                                  ,iqc(1),iqc(2))
 980         endif
 981       endif
 982       if(ish.ge.5)write(ifch,*)'si,pt2,wwmin,s2min,sj,sjb,iqc:'
 983       if(ish.ge.5)write(ifch,*)si,pt2,wwmin,s2min,sj,sjb,iqc
 984
 985       if(si.lt.1.1d0*dble(wwmin))goto 12 !------>>>>>>>
 986       if(rangen().lt.sjb/sj)goto 12 !------>>>>>>>
 987
 988       if(iabs(iqc(1)).eq.4)jj=min(2,int(sjord/sj+rangen())+1)  !?????????
 989
 990 112   continue
 991
 992       if(iabs(iqc(jj)).ne.4)then
 993
 994         discr=dble((sngl(si)+2.*pt*sqrt(q2ini)-q2mass)**2
 995      *  -4.*q2ini*(5.*sngl(si)+q2mass+sngl(pt2)))
 996         if(discr.lt.0.d0.and.ish.ge.1)write(ifmt,*)'discr,si,pt,wwmin',
 997      *  discr,si,pt,wwmin
 998         discr=dsqrt(discr)
 999         qqmax=(si+2.d0*dble(pt*sqrt(q2ini))-dble(q2mass)+discr)/2.d0
1000      *  /(5.d0+(dble(q2mass)+pt2)/si)
1001         qqmin=qqmax-discr/(5.d0+(dble(q2mass)+pt2)/si)
1002         if(s2min.gt.4.d0*qqmin+dble(q2mass))then
1003           xmm=.5d0*(si-s2min+2.d0*dble(pt*sqrt(q2ini)))
1004           discr=xmm**2-si*dble(q2ini)*(1.d0+(dble(q2mass)+pt2)/si)
1005           if(discr.lt.0..and.ish.ge.1)write(ifmt,*)'discr1,si,pt,wwmin',
1006      *    discr,si,pt,wwmin
1007           qqmin=(xmm-dsqrt(discr))/(1.d0+(dble(q2mass)+pt2)/si)
1008         endif
1009
1010         xmax=min(1.d0-dble(q2ini)/qqmax,.9999d0)
1011         if(qqmin.lt.dble(qmin(jj)))then
1012           qqmin=dble(qmin(jj))
1013           xmin=max(1.d0-((dble(pt)*dsqrt(qqmin)+dsqrt(pt2*qqmin+
1014      *    si*(si-s2min-qqmin*(1.d0+(dble(q2mass)+pt2)/si))))/si)**2,
1015      *    (s2min+qqmin*(1.d0+(dble(q2mass)+pt2)/si)-2.d0*dble(pt)
1016      *    *dsqrt(qqmin))/si)
1017           if(xmin.le.0.d0)xmin=s2min/si
1018         else
1019           xmin=1.d0-dble(q2ini)/qqmin
1020         endif
1021
1022         qm0=qmin(jj)
1023         xm0=1.0-dble(q2ini/qm0)
1024         if(xm0.gt..999*xmax.or.xm0.lt.1.001*xmin)then
1025           xm0=.5d0*(xmax+xmin)
1026         endif
1027         s2max=sngl(xm0*si-qm0*(dble(q2mass)+pt2)/si)
1028      *       +2.*pt*sqrt(q2ini)
1029         xx=xm0
1030
1031         if(iabs(iqc(1)).ne.4)then
1032           if(jfirst.eq.1)then
1033             sj0=psjti(qm0,qqcut,s2max,0,iqc(3-jj),0)*
1034      *      psfap(xx,iqc(jj),0)+
1035      *      psjti(qm0,qqcut,s2max,7,iqc(3-jj),0)
1036      *      *psfap(xx,iqc(jj),1)
1037           else
1038             sj0=psjti1(qm0,qmin(3-jj),qqcut,s2max,0,iqc(3-jj),0)*
1039      *      psfap(xx,iqc(jj),0)+
1040      *      psjti1(qm0,qmin(3-jj),qqcut,s2max,7,iqc(3-jj),0)
1041      *      *psfap(xx,iqc(jj),1)
1042           endif
1043         else
1044           sj0=psjci(qm0,s2max,0)*psfap(xx,iqc(jj),0)+
1045      *    psjci(qm0,s2max,1)*psfap(xx,iqc(jj),1)
1046         endif
1047
1048         gb0=dble(sj0/log(q2ini/qcdlam)*qm0*5.)*psuds(qm0,iqc(jj))
1049         if(ish.ge.5)write(ifch,*)'gb0,qm0,xm0,s2max:',
1050      *  gb0,qm0,xm0,s2max
1051 c        if(gb0.le.0.)stop'gb0<=0'
1052         if(gb0.le.0.d0)then
1053           write(ifmt,*)'gb0.le.0.  si,pt2:',si,pt2
1054           iret=1
1055           goto 16
1056         endif
1057
1058         if(xm0.le..5d0)then
1059           gb0=gb0*xm0**(1.-delh)
1060         else
1061           gb0=gb0*(1.d0-xm0)*2.d0**delh
1062         endif
1063
1064         xmin2=max(.5d0,xmin)
1065         xmin1=xmin**delh
1066         xmax1=min(xmax,.5d0)**delh
1067         if(xmin.ge..5d0)then
1068           djl=1.
1069         elseif(xmax.le..5d0)then
1070           djl=0.
1071         else
1072           djl=1./(1.+((2.*sngl(xmin))**delh-1.)/delh/
1073      *    log(2.*(1.-sngl(xmax))))
1074         endif
1075
1076       else
1077
1078         xmin=5.d0*dble(q2min)/si
1079         xmax=min(si/(si+5.0*(pt2+dble(q2mass))),.9999d0)
1080         qqmax=xmax*si/5.d0
1081         qqmin=dble(q2min)
1082
1083         qm0=sngl(qqmin)
1084         xm0=2.d0/(1.d0+sqrt(1.d0+4.d0*(pt2+dble(q2mass))/qm0))
1085         if(xm0.gt..999d0*xmax.or.xm0.lt.1.001d0*xmin)then
1086           xm0=.5d0*(xmax+xmin)
1087         endif
1088         s2max=sngl(xm0*si)
1089         xx=xm0
1090
1091         sj0=psjti(qm0,qmin(3-jj),s2max,0,iqc(3-jj),0)*
1092      *  psfap(xx,iqc(jj),0)
1093         gb0=dble(sj0/log(qm0/qcdlam)*qm0   *5.)
1094         gb0=gb0*xm0**(1.-delh)
1095         if(gb0.le.0.d0)then
1096           if(ish.ge.2)write(ifch,*)'gb0.le.0. (charm)  si,pt2:',si,pt2
1097           iret=1
1098           goto 16
1099         endif
1100         djl=0.
1101         xmin1=xmin**delh
1102         xmax1=xmax**delh
1103
1104       endif
1105
1106       if(ish.ge.5)write(ifch,*)'xmin,xmax,qqmin,qqmax:',
1107      *xmin,xmax,qqmin,qqmax
1108
1109       ntry=0
1110 10    continue
1111       ntry=ntry+1
1112       if(ntry.gt.5000000)then
1113         if(ish.ge.1)write(*,*)'Reject hard string (too many rejection)'
1114      &,kcol,ncolp,nptl,gb7
1115         iret=1
1116         goto 16
1117       endif
1118       if(rangen().gt.djl)then        !lc momentum share in the cur
1119         x=(xmin1+dble(rangen())*(xmax1-xmin1))**(1./delh)
1120       else
1121         x=1.d0-(1.d0-xmin2)*((1.d0-xmax)/(1.d0-xmin2))**rangen()
1122       endif
1123       qq=qqmin/(1.d0+dble(rangen())*(qqmin/qqmax-1.d0))
1124       if(ish.ge.7)write(ifch,*)'x,qq:',x,qq,ntry
1125
1126       if(iabs(iqc(jj)).ne.4)then
1127
1128         qt2=qq*(1.d0-x)
1129         if(qt2.lt.dble(q2ini))then
1130           if(ish.ge.7)write(ifch,*)'qt2:',qt2
1131           goto 10
1132         endif
1133
1134         qmin2=max(qminn,sngl(qq))
1135         qt=dsqrt(qt2)
1136         call pscs(bcos,bsin)
1137         ep3(3)=sngl(qt)*bcos   !new parton pt-s
1138         ep3(4)=sngl(qt)*bsin
1139         ptnew=(ept(3)-dble(ep3(3)))**2+(ept(4)-dble(ep3(4)))**2
1140         s2min2=4.*qmin2+q2mass
1141
1142         s2=x*(si-qq)-ptnew-qq*(dble(q2mass)+pt2)/si+pt2  !new ladder mass
1143         if(s2.lt.dble(s2min2))then
1144           if(ish.ge.7)write(ifch,*)'s2,s2min2:',s2,s2min2
1145           goto 10  !rejection in case of too low mass
1146         endif
1147
1148         xx=x
1149         if(iabs(iqc(1)).ne.4)then
1150           if(jfirst.eq.1)then
1151             sj1=psjti(sngl(qq),qqcut,sngl(s2),0,iqc(3-jj),0)
1152             if(iqc(jj).ne.0)then                             !f1 - f2
1153               sj2=psjti(sngl(qq),qqcut,sngl(s2),iqc(jj),iqc(3-jj),0)
1154             elseif(iqc(3-jj).eq.0)then                       !f1 - g
1155               sj2=psjti(sngl(qq),qqcut,sngl(s2),1,0,0)
1156             else                                             !g  - f2
1157               sj2=psjti(sngl(qq),qqcut,sngl(s2),1,1,0)/naflav/2.      !q  - q
1158      *        +psjti(sngl(qq),qqcut,sngl(s2),-1,1,0)/naflav/2.        !q~ - q
1159      *        +psjti(sngl(qq),qqcut,sngl(s2),1,2,0)*(1.-1./naflav)    !q' - q
1160             endif
1161           else
1162             sj1=psjti1(sngl(qq),qmin(3-jj),qqcut,sngl(s2),0,iqc(3-jj),0)
1163             if(iqc(jj).ne.0)then                             !f1 - f2
1164               sj2=psjti1(sngl(qq),qmin(3-jj),qqcut
1165      *                   ,sngl(s2),iqc(jj),iqc(3-jj),0)
1166             elseif(iqc(3-jj).eq.0)then                       !f1 - g
1167               sj2=psjti1(sngl(qq),qmin(3-jj),qqcut,sngl(s2),1,0,0)
1168             else                                             !g  - f2
1169               sj2=psjti1(sngl(qq),qmin(3-jj),qqcut
1170      *                         ,sngl(s2),1,1,0)/naflav/2.         !q - q
1171      *           +psjti1(sngl(qq),qmin(3-jj),qqcut
1172      *                         ,sngl(s2),-1,1,0)/naflav/2.   !q~ - q
1173      *           +psjti1(sngl(qq),qmin(3-jj),qqcut
1174      *                      ,sngl(s2),1,2,0)*(1.-1./naflav)!q' - q
1175             endif
1176           endif
1177         else
1178           sj1=psjci(sngl(qq),sngl(s2),0)
1179           sj2=psjci(sngl(qq),sngl(s2),1)
1180         endif
1181
1182         !...gb7 - rejection function for x and q**2 simulation
1183         gb7=dble((sj1*psfap(xx,iqc(jj),0)+sj2*psfap(xx,iqc(jj),1))/
1184      *  log(sngl(qt2)/qcdlam))*psuds(sngl(qq),iqc(jj))*qq/gb0
1185
1186         if(x.le..5d0)then
1187           gb7=gb7*x**(1.-delh)
1188         else
1189           gb7=gb7*(1.d0-x)*2.d0**delh
1190         endif
1191         if(gb7.gt.1..and.ish.ge.2)write(ifmt,*)'gb7,qq,x,qm0,xm0',
1192      *  gb7,qq,x,qm0,xm0
1193         if(ish.ge.7)write(ifch,*)'gb7:',gb7
1194         if(dble(rangen()).gt.gb7)goto 10
1195       else
1196         qmin2=max(qminn,sngl(qq))
1197         s2min2=4.*qmin2
1198         s2=x*si-qq               !new ladder mass
1199         if(s2.lt.dble(s2min2))goto 10  !rejection in case of too low mass
1200
1201         call pscs(bcos,bsin)
1202         xmm=x*(ept(3)*dble(bcos)+ept(4)*dble(bsin))
1203         discr=xmm**2+qq*(1.d0-x)-x**2*(pt2+dble(q2mass))
1204         if(discr.lt.0.d0)goto 10
1205         qt=xmm+dsqrt(discr)
1206         ep3(3)=sngl(ept(3)-qt*dble(bcos))   !new parton pt-s
1207         ep3(4)=sngl(ept(4)-qt*dble(bsin))
1208         qt2=dble(ep3(3)**2+ep3(4)**2)
1209
1210         xx=x
1211         sj1=psjti(sngl(qq),qqcut,sngl(s2),0,iqc(3-jj),0)
1212         sj2=0.
1213         !.... gb7 - rejection function for x and q**2 simulation
1214         gb7=dble(sj1*psfap(xx,iqc(jj),0)/
1215      *  log(sngl(qq)/qcdlam))*qq/gb0
1216         gb7=gb7*x**(1.-delh)
1217         if(gb7.gt.1..and.ish.ge.2)write(ifmt,*)'gb7,qq,x,qm0,xm0',
1218      *  gb7,qq,x,qm0,xm0
1219         if(ish.ge.7)write(ifch,*)'gb7:',gb7
1220         if(dble(rangen()).gt.gb7)goto 10
1221       endif
1222
1223       nqc(2)=0
1224       iqnew=iqc(jj)
1225       if(rangen().lt.sj1/(sj1+sj2))then
1226         if(iqc(jj).eq.0)then
1227           jt=1
1228           jq=int(1.5+rangen())
1229           nqc(1)=ncc(jq,jj)
1230         else
1231           jt=2
1232           if(iqc(jj).gt.0)then
1233             jq=1
1234           else
1235             jq=2
1236           endif
1237           nqc(1)=0
1238           iqnew=0
1239         endif
1240         iq1=iqc(jj)
1241
1242       else
1243         if(iqc(jj).ne.0)then
1244           iq1=0
1245           jt=3
1246           if(iqc(jj).gt.0)then
1247             jq=1
1248           else
1249             jq=2
1250           endif
1251           nqc(1)=ncc(1,jj)
1252         else
1253           jt=4
1254           jq=int(1.5+rangen())
1255           iq1=int(naflav*rangen()+1.)*(3-2*jq)
1256           iqnew=-iq1
1257           nqc(1)=ncc(jq,jj)
1258         endif
1259       endif
1260       if(iabs(iqc(jj)).ne.4)then
1261         q2part=0.
1262 c        write(*,*)'sem:',wpt(jj),pt2,q2mass,wpt(3-jj),jj
1263         if(dabs(wpt(3-jj)).gt.1.d-20)then
1264           wplc=wpt(jj)-(pt2+dble(q2mass))/wpt(3-jj)
1265         else
1266           if(ish.gt.1)write(ifmt,*)'Problem with wpt in sem',wpt(3-jj)
1267           wplc=wpt(jj)-(pt2+dble(q2mass))
1268         endif
1269         qp2max=max(0.,sngl(qt2))
1270       else
1271         q2part=q2mass
1272         wplc=wpt(jj)
1273         qp2max=max(0.,sngl(qq))
1274       endif
1275       eprt=max(dsqrt(qt2+dble(q2part)),.5d0*((1.d0-x)*wplc+
1276      *(qt2+dble(q2part))/(1.d0-x)/wplc))
1277       pl=((1.d0-x)*wplc-eprt)*dble(3-2*jj)
1278       if(iq1.eq.0)then
1279         iq2ini=9
1280       else
1281         iq2ini=iq1
1282       endif
1283       if(ish.ge.5)write(ifch,*)'qq,eprt,iq2ini',qq,eprt,iq2ini,
1284      *eprt**2-qt2
1285       ntest=0
1286 11    ntest=ntest+1
1287       call timsh1(qp2max,sngl(eprt),iq2ini)
1288       amprt=pprt(5,1)**2
1289       plprt=max(0.d0,eprt**2-qt2)-dble(amprt)
1290       if(plprt.lt.0.d0)then
1291         if(ntest.lt.10000)then
1292           goto 11
1293         else
1294           iret=1
1295           goto 16
1296         endif
1297       endif
1298       ep3(1)=sngl(eprt)
1299       ep3(2)=sngl(dsqrt(plprt))
1300       if(pl.lt.0.d0)ep3(2)=-ep3(2)
1301       ey(1)=1.
1302       ey(2)=1.
1303       ey(3)=1.
1304       do i=1,4
1305         ept1(i)=ept(i)-dble(ep3(i))
1306       enddo
1307       call psdefrot(ep3,s0xh,c0xh,s0h,c0h)
1308       s2new=psnorm(ept1)
1309
1310       if(iabs(iqc(jj)).ne.4.and.s2new-q2mass.gt.s2min2.or.
1311      *iabs(iqc(jj)).eq.4.and.s2new.gt.s2min2)then
1312         if(iabs(iqc(1)).ne.4.or.jj.eq.1)then
1313           if(jfirst.eq.1)then
1314             gb=dble(psjti(sngl(qq),qqcut,s2new,iqnew,iqc(3-jj),0))
1315           else
1316             gb=dble(psjti1(sngl(qq),qmin(3-jj),qqcut,s2new,iqnew,
1317      *                                               iqc(3-jj),0))
1318           endif
1319         else
1320           gb=dble(psjci(sngl(qq),s2new-q2mass,iqnew))
1321         endif
1322         if(iqnew.eq.0)then
1323           gb=gb/dble(sj1)
1324         else
1325           gb=gb/dble(sj2)
1326         endif
1327         if(dble(rangen()).gt.gb)goto 10
1328       else
1329         goto 10
1330       endif
1331
1332       if(ish.ge.6)write(ifch,*)'jt,jj,jq,nqc:',jt,jj,jq,nqc
1333       nprtjx=0
1334       call psreti(nqc,jq,1,ey,s0xh,c0xh,s0h,c0h)
1335
1336       if(jt.eq.1)then
1337         ncc(jq,jj)=nqc(2)
1338       elseif(jt.eq.2)then
1339         ncc(jq,jj)=ncc(1,jj)
1340         ncc(3-jq,jj)=nqc(1)
1341       elseif(jt.eq.3)then
1342         ncc(1,jj)=nqc(2)
1343       elseif(jt.eq.4)then
1344         ncc(1,jj)=ncc(3-jq,jj)
1345       endif
1346       iqc(jj)=iqnew
1347
1348       do i=1,4
1349         ept(i)=ept1(i)
1350       enddo
1351       ! c.m. energy squared, minimal  4-momentum transfer square and gluon 4-v
1352       ! for the next ladder run
1353       qmin(jj)=sngl(qq)
1354       qminn=qmin2
1355       si=dble(s2new)
1356       nemis(jj)=nemis(jj)+1
1357
1358       goto 9  !  ---------------next ladder rung ------------>>>>>>>>>>>
1359
1360
1361  12   continue !------------------- Born process------------------------
1362
1363       if(ish.ge.4)write(ifch,*)'Born process'
1364       if(ish.ge.6)write(ifch,*)'iqc,si:',iqc,si
1365       qq=dble(qminn)
1366
1367    !kkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkk
1368       xpprbor(ncolp,kcol)=(ept(1)+ept(2))/dble(engy)
1369       xmprbor(ncolp,kcol)=(ept(1)-ept(2))/dble(engy)
1370       nemispr(1,ncolp,kcol)=nemis(1)
1371       nemispr(2,ncolp,kcol)=nemis(2)
1372    !   write(*,'(a,2f8.3,i3,3x,2i3)')'------------'
1373    !  *          ,xpprbor(ncolp,kcol),xmprbor(ncolp,kcol),nj-nj00,nemis
1374    !kkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkk
1375
1376       call psabor(sngl(si),sngl(qq),iqc,ncc,ept,0,iptl,bx)
1377
1378    !kkkkkkkkkkkkk out Born partons without timelike cascade
1379       if(nprtjx.eq.2)then
1380        do ii=1,2
1381         ptprboo(ii,ncolp,kcol)=sqrt(pprtx(1,ii)**2+pprtx(2,ii)**2)
1382        enddo
1383       else
1384         stop'psahot: should not happen!!!!                '
1385       endif
1386    !kkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkk
1387    !     ptxxx=max(  ptprboo(1,ncolp,kcol) , ptprboo(2,ncolp,kcol)  )
1388    !            print*,sqrt(wm0*wp0),sqrt(wmi*wpi),ptxxx  !++++++++++++++++++++++++++
1389    !      if(ptxxx.lt.10)goto1
1390    !            print*,'  ++++++++++++++++++++++++++++++++++++ '
1391    !kkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkk
1392       if(nj.ne.0.)then
1393         do i=1,nj
1394           do l=1,6
1395             bxj(l,i)=bx(l)
1396           enddo
1397           ityj(i)=ity
1398           iorj(i)=iptl
1399         enddo
1400       endif
1401
1402       call psjarr(jfl)       !kinky strings formation
1403
1404       if(jfl.eq.0.and.ish.ge.4)write(ifch,*)
1405      *'jfl,nj,nptl',jfl,nj,nptl
1406       if(jfl.eq.0)goto 1
1407
1408 c --- update remnant flavour ---
1409
1410       iret1=0
1411       call idenco(jcp,icp,iret1)
1412       if(iret1.eq.1)call utstop('Problem with proj rem in psahot!&')
1413       iret2=0
1414       call idenco(jct,ict,iret2)
1415       if(iret2.eq.1)call utstop('Problem with targ rem in psahot!&')
1416
1417       do i=1,2
1418         icproj(i,ip)=icp(i)
1419         ictarg(i,it)=ict(i)
1420       enddo
1421
1422       if(ish.ge.3)write(ifch,*)'End psahot (icp,ict):',icp,ict
1423
1424       if(iep(ip).ne.1)then
1425       call idtr4(idptl(ip),icini) !excited remnant ?
1426       if(ish.ge.5)write(ifch,*)'icini proj',icini
1427      &     ,(icp(1)-icini(1)),(icp(2)-icini(2))
1428       if((icp(1)-icini(1))+(icp(2)-icini(2)).ne.0)iep(ip)=1
1429       endif
1430       if(iet(it).ne.1)then
1431       call idtr4(idptl(maproj+it),icini)
1432       if(ish.ge.5)write(ifch,*)'icini targ',icini
1433      &     ,(ict(1)-icini(1)),(ict(2)-icini(2))
1434       if((ict(1)-icini(1))+(ict(2)-icini(2)).ne.0)iet(it)=1
1435       endif
1436       if(ish.ge.5)write(ifch,*)'iep,iet ',iep(ip),iet(it)
1437
1438 c ------------------------------
1439
1440 16    continue
1441       call utprix('psahot',ish,ishini,3)
1442
1443       q2fin=q2finsave
1444       return
1445       end
1446
1447 c------------------------------------------------------------------------
1448       subroutine psjarr(jfl)
1449 c-----------------------------------------------------------------------
1450 c
1451 c   final jets rearrangement according to their colour connection
1452 c   and write to /cptl/
1453 c
1454 c jfl - flag for the rejection (in case of jfl=0)
1455 c-----------------------------------------------------------------------
1456 c Input:
1457
1458       parameter (mjstr=20000)
1459       common /psar29/ eqj(4,mjstr),iqj(mjstr),ncj(2,mjstr),nj
1460
1461 c eqj(1:4,k) - 4-momentum (qgs) for k-th parton;
1462 c bxj(1:4,k) - coordinates for k-th parton formation point;
1463 c iqj(k) - ID (qgs) for k-th parton;
1464 c ncj(1:2,k) - colour connected partons indexes for k-th parton;
1465 c nj - number of partons
1466 c-----------------------------------------------------------------------
1467       dimension mark(mjstr)
1468       double precision  ept(4)
1469       include 'epos.inc'
1470       include 'epos.incsem'
1471
1472       if(nj.eq.0)then
1473         jfl=1
1474         return
1475       endif
1476 c      do k=1,nj  !???????????????
1477 c       eqj(1,k)=dsqrt(0d0+eqj(2,k)**2+eqj(3,k)**2+eqj(4,k)**2)
1478 c      enddo
1479       if(ish.ge.3)then
1480         write (ifch,*)'psjarr: nj',nj
1481         do k=1,nj
1482           if(iqj(k).ne.0)ncj(2,k)=0
1483           write(ifch,'(a,i4)')' parton',k
1484           write(ifch,'(i3,2x,4e10.3,2x,2i3)')iqj(k)
1485      *    ,(eqj(j,k),j=1,4),(ncj(j,k),j=1,2)
1486         enddo
1487       endif
1488
1489       jfl=0
1490       do i=1,nj
1491         mark(i)=1
1492       enddo
1493       nptl0=nptl
1494       if(ish.ge.5)write (ifch,*)'nptl,nj',nptl,nj
1495
1496 1     continue
1497       do ij=1,nj
1498         if(mark(ij).ne.0.and.iqj(ij).ne.0)goto 2
1499       enddo
1500 2     continue
1501       jfirst=1
1502       if(iabs(iqj(ij)).le.2)then
1503         am1=amhadr(1)
1504       elseif(iabs(iqj(ij)).eq.4)then
1505         am1=amhadr(3)
1506       elseif(iabs(iqj(ij)).eq.40)then
1507         am1=qcmass
1508       else
1509         am1=amhadr(2)
1510       endif
1511
1512       do i=1,4
1513         ept(i)=0.
1514       enddo
1515
1516 3     continue
1517       call pspawr(ij)
1518       mark(ij)=0
1519
1520       do i=1,4
1521         ept(i)=ept(i)+eqj(i,ij)
1522       enddo
1523
1524       if(iqj(ij).ne.0)then
1525         if(jfirst.ne.1)then
1526           if(iabs(iqj(ij)).le.2)then
1527             am2=amhadr(1)
1528           elseif(iabs(iqj(ij)).eq.4)then
1529             am2=amhadr(3)
1530           elseif(iabs(iqj(ij)).eq.40)then
1531             am2=qcmass
1532           else
1533             am2=amhadr(2)
1534           endif
1535           amj=(am1+am2+stmass)**2
1536           sm=psnorm(ept)
1537           if(sm.lt.amj)then
1538             nptl=nptl0
1539             goto999
1540           endif
1541
1542           if(nptl-nptl0.lt.nj)then
1543             goto 1
1544           else
1545             if(ish.ge.3)then
1546               write (ifch,*)'psjarr: nptl',nptl
1547               do k=nptl0,nptl
1548                 write(ifch,'(a,i4)')' particle',k
1549                 write(ifch,'(i5,2x,5e10.3)')idptl(k)
1550      *          ,(pptl(j,k),j=1,5)
1551               enddo
1552             endif
1553             jfl=1
1554             goto999
1555           endif
1556         else
1557           jfirst=0
1558           njpar=ij
1559           ij=ncj(1,ij)
1560           goto 3
1561         endif
1562
1563       else
1564         if(ncj(1,ij).eq.njpar)then
1565           njdau=ncj(2,ij)
1566         else
1567           njdau=ncj(1,ij)
1568         endif
1569         njpar=ij
1570         ij=njdau
1571         goto 3
1572       endif
1573   999 continue
1574       end
1575
1576 c------------------------------------------------------------------------
1577       subroutine pspawr(kj)
1578 c-----------------------------------------------------------------------
1579 c pspawr - writing final parton kj into particle list
1580 c------------------------------------------------------------------------
1581 c Input:
1582
1583       parameter (mjstr=20000)
1584       common /psar29/ eqj(4,mjstr),iqj(mjstr),ncj(2,mjstr),nj
1585       common /psar30/ iorj(mjstr),ityj(mjstr),bxj(6,mjstr)
1586
1587 c eqj(1:4,k) - 4-momentum (qgs) for k-th parton;
1588 c bxj(1:4,k) - coordinates for k-th parton formation point;
1589 c iqj(k) - ID (qgs) for k-th parton;
1590 c ncj(1:2,k) - colour connected partons indexes for k-th parton;
1591 c nj - number of partons
1592 c-----------------------------------------------------------------------
1593       include 'epos.inc'
1594       include 'epos.incsem'
1595       common/ciptl/iptl
1596
1597       if(ish.ge.9)write (ifch,*)'nptl,kj (sto)',nptl,kj
1598       if(nptl.ge.mxptl)then
1599        write (ifmt,*)'nptl,kj',nptl,kj
1600        call alist('Error in pspawr: nptl out of bounds &',1,nptl)
1601        stop'nptl out of bounds'
1602       endif
1603
1604       if(ifrptl(1,iptl).eq.0)ifrptl(1,iptl)=nptl+1
1605       ifrptl(2,iptl)=nptl+1
1606
1607       nptl=nptl+1
1608       pptl(1,nptl)=eqj(3,kj)
1609       pptl(2,nptl)=eqj(4,kj)
1610       pptl(3,nptl)=eqj(2,kj)
1611       pptl(4,nptl)=eqj(1,kj)
1612       pptl(5,nptl)=0.
1613       idptl(nptl)=psidd(iqj(kj))
1614       iorptl(nptl)=iorj(kj)
1615       jorptl(nptl)=0
1616       istptl(nptl)=20
1617       do i=1,4
1618         xorptl(i,nptl)=bxj(i,kj)
1619       enddo
1620       tivptl(1,nptl)=bxj(5,kj)
1621       tivptl(2,nptl)=bxj(6,kj)
1622       ityptl(nptl)=ityj(kj)
1623           !kkkkkkkkkkkkkkkkkkkkkkk
1624        !    write(*,'(a,2i4,i6,f8.3)')'.... ',kj,nptl,idptl(nptl)
1625        !*     ,sqrt(pptl(1,nptl)**2+pptl(2,nptl)**2)
1626       return
1627       end
1628
1629 c------------------------------------------------------------------------
1630       subroutine psabor(si,qq,iqc,ncc,ept,jdis,iptl,coordo)
1631 c------------------------------------------------------------------------
1632 c psabor - highest virtuality subprocess in the ladder
1633 c si - c.m. energy squared for the process,
1634 c qq - p_t-cutoff for the process due to the evolution,
1635 c iqc(i), i=1,2 - incoming parton types(0-g,(+-)1-u(u~),(+-)2-d(d~)etc.),
1636 c ncc(i,j), i,j=1,2 - incoming partons color connections,
1637 c ept(4) - total 4-momentum for the system of the 2 partons
1638 c jdis=0 - hadronic process in pp; 1 - resolved photon process
1639 c------------------------------------------------------------------------
1640       double precision ept(4),psutz,psuds
1641       dimension ep3(4),ey(3),wsub(5),iqc(2),ncc(2,2),nqc(2),coordo(6)
1642       parameter (mjstr=20000)
1643       common /psar29/eqj(4,mjstr),iqj(mjstr),ncj(2,mjstr),nj
1644       parameter (ntim=1000)
1645       common/cprt/nprtj,pprt(5,ntim),idprt(ntim),iorprt(ntim)
1646      &,idaprt(2,ntim)
1647       common/cprtx/nprtjx,pprtx(5,2)
1648       include 'epos.incsem'
1649       include 'epos.inc'
1650       call utpri('psabor',ish,ishini,5)
1651
1652       if(iabs(iqc(1)).ne.4)then   !gluon or light quark
1653         q2mass=0.
1654       else                        !c-quark
1655         q2mass=qcmass**2
1656       endif
1657       p1=si/(1.+q2mass/si)
1658       if(p1.gt.4.*qq)then                 !|t|-cutoff (p^2_t>qq)
1659         tmin=2.*qq/(1.+sqrt(1.-4.*qq/p1))
1660       else
1661         tmin=2.*qq
1662       endif
1663       tmax=.5*p1
1664
1665       fborn=0.
1666       qt2=tmin*(1.d0-tmin/p1)
1667       if(qt2.lt..999d0*qq.and.ish.ge.2)write(ifmt,*)'qt20,qq',qt2,qq
1668       if(iqc(1).ne.0.or.iqc(2).ne.0)then
1669         do l=1,5       !sum over different subprocesses
1670           wsub(l)=psbori(si,tmin,iqc(1),iqc(2),l)  !matrix element
1671         if(l.le.3)then
1672           wsub(l)=wsub(l)*pssalf(qt2/qcdlam)**2
1673         elseif(l.le.4)then
1674           wsub(l)=wsub(l)*pssalf(qt2/qcdlam)*alfe/2/pi
1675         else
1676           wsub(l)=wsub(l)*(alfe/2/pi)**2
1677         endif
1678           fborn=fborn+wsub(l)
1679         enddo
1680         fborn=tmin**2*fborn
1681       else
1682         do l=1,5
1683           wsub(l)=psbori(si,.5*si,iqc(1),iqc(2),l)
1684         if(l.le.3)then
1685           wsub(l)=wsub(l)*pssalf(qt2/qcdlam)**2
1686         elseif(l.le.4)then
1687           wsub(l)=wsub(l)*pssalf(qt2/qcdlam)*alfe/2/pi
1688         else
1689           wsub(l)=wsub(l)*(alfe/2/pi)**2
1690         endif
1691           fborn=fborn+wsub(l)
1692         enddo
1693         fborn=.25*si**2*fborn
1694       endif
1695        if(jdis.eq.0)then
1696          scale1=qt2
1697        else
1698          scale1=4.*qt2
1699        endif
1700       gb0=dble(fborn)*psuds(scale1,iqc(1))*psuds(qt2,iqc(2))
1701       if(ish.ge.7)write(ifch,*)'tmin,gb0:',tmin,gb0
1702
1703 c------------------------------------------------
1704 c 4-momentum transfer squared is simulated first as dq_t**2/q_t**4 from
1705 c tmin to s/2
1706 14    q2=tmin/(1.-rangen()*(1.-tmin/tmax))    !q2=min(|t|,|u|)
1707       qt2=q2*(1.-q2/p1)                       !qt2=p_t^2 for the process
1708       if(qt2.lt.qq.and.ish.ge.2)write(ifmt,*)'qt2,qq',qt2,qq
1709       if(rangen().lt..5)then  !|u|=q2, |t|=p1-q2
1710         jm=2                  !first parton to be considered
1711         tq=p1-q2
1712       else                    !|t|=q2, |u|=p1-q2
1713         jm=1                  !first parton to be considered
1714         tq=q2
1715       endif
1716
1717       fborn=0.
1718       do l=1,5                       !sum over different subprocesses
1719         wsub(l)=psbori(si,tq,iqc(1),iqc(2),l)
1720         if(l.le.3)then
1721           wsub(l)=wsub(l)*pssalf(qt2/qcdlam)**2
1722         elseif(l.le.4)then
1723           wsub(l)=wsub(l)*pssalf(qt2/qcdlam)*alfe/2/pi
1724         else
1725           wsub(l)=wsub(l)*(alfe/2/pi)**2
1726         endif
1727         fborn=fborn+wsub(l)
1728       enddo
1729        if(jdis.eq.0)then
1730          scale1=qt2
1731        else
1732          scale1=4.*qt2
1733        endif
1734       gb=dble(q2**2*fborn)
1735      &*psuds(scale1,iqc(1))*psuds(qt2,iqc(2))/gb0 !rejection function
1736       if(ish.ge.7)write(ifch,*)'q2,qt2,gb:',q2,qt2,gb
1737
1738       if(dble(rangen()).gt.gb)goto 14                   !rejection
1739
1740 c determination of the color configuration
1741       nqc(2)=0
1742       if(iqc(1).eq.0.and.iqc(2).eq.0)then      !g+g
1743         jq=int(1.5+rangen())    !jq=1(2) - transfer of color (anticolor)
1744         nqc(1)=ncc(jq,jm)
1745
1746         if(rangen().lt.wsub(1)/fborn)then      !gg->gg
1747           if(rangen().lt..5)then
1748             jt=1                !anticolor-color annihilation
1749             nqc(2)=0
1750             njc1=ncc(3-jq,jm)
1751             njc2=ncc(jq,3-jm)
1752             if(iqj(njc1).ne.0)then
1753               ncj(1,njc1)=njc2
1754             else
1755               ncj(jq,njc1)=njc2
1756             endif
1757             if(iqj(njc2).ne.0)then
1758               ncj(1,njc2)=njc1
1759             else
1760               ncj(3-jq,njc2)=njc1
1761             endif
1762           else
1763             jt=2                    !produced gluons get color and
1764             nqc(2)=ncc(3-jq,3-jm)   !anticolor from the 2 parents
1765           endif
1766         else                                   !gg->qq~
1767           jt=9                  !anticolor-color annihilation
1768           iqc(jm)=int(naflav*rangen()+1)*(3-2*jq) !(anti)quark flavor
1769           iqc(3-jm)=-iqc(jm)
1770           njc1=ncc(3-jq,jm)
1771           njc2=ncc(jq,3-jm)
1772           if(iqj(njc1).ne.0)then
1773             ncj(1,njc1)=njc2
1774           else
1775             ncj(jq,njc1)=njc2
1776           endif
1777           if(iqj(njc2).ne.0)then
1778             ncj(1,njc2)=njc1
1779           else
1780             ncj(3-jq,njc2)=njc1
1781           endif
1782         endif
1783
1784       elseif(iqc(1)*iqc(2).eq.0)then       !q(q~)+g
1785         if(iqc(1)+iqc(2).gt.0)then
1786           jq=1                             !q
1787         else
1788           jq=2                             !q~
1789         endif
1790        if(rangen().lt.wsub(1)/fborn)then  !q(q~)g->q(q~)g
1791         if(rangen().lt..5)then      !anticolor-color annihilation
1792           if(iqc(jm).eq.0)then
1793             jt=3                    !first parton=g
1794             nqc(1)=ncc(jq,jm)
1795             njc1=ncc(3-jq,jm)
1796             njc2=ncc(1,3-jm)
1797             if(iqj(njc1).ne.0)then
1798               ncj(1,njc1)=njc2
1799             else
1800               ncj(jq,njc1)=njc2
1801             endif
1802             if(iqj(njc2).ne.0)then
1803               ncj(1,njc2)=njc1
1804             else
1805               ncj(3-jq,njc2)=njc1
1806             endif
1807
1808           else
1809             jt=4                    !first parton=q(q~)
1810             nqc(1)=0
1811             njc1=ncc(1,jm)
1812             njc2=ncc(3-jq,3-jm)
1813             if(iqj(njc1).ne.0)then
1814               ncj(1,njc1)=njc2
1815             else
1816               ncj(3-jq,njc1)=njc2
1817             endif
1818             if(iqj(njc2).ne.0)then
1819               ncj(1,njc2)=njc1
1820             else
1821               ncj(jq,njc2)=njc1
1822             endif
1823           endif
1824
1825         else                        !color transfer
1826           if(iqc(jm).eq.0)then
1827             jt=5                    !first parton=g
1828             nqc(2)=ncc(3-jq,jm)
1829             nqc(1)=ncc(1,3-jm)
1830           else                      !first parton=q(q~)
1831             jt=6
1832             nqc(1)=ncc(jq,3-jm)
1833           endif
1834         endif
1835
1836        else          !q(q~)g->q(q~)-gamma (+-color annihilation)
1837           if(iqc(jm).eq.0)then
1838             jt=11                    !first parton=g
1839             nqc(1)=ncc(jq,jm)
1840             njc1=ncc(3-jq,jm)
1841             njc2=ncc(1,3-jm)
1842             if(iqj(njc1).ne.0)then
1843               ncj(1,njc1)=njc2
1844             else
1845               ncj(jq,njc1)=njc2
1846             endif
1847             if(iqj(njc2).ne.0)then
1848               ncj(1,njc2)=njc1
1849             else
1850               ncj(3-jq,njc2)=njc1
1851             endif
1852             iqc(jm)=iqc(3-jm)
1853             iqc(3-jm)=10              !make the second output is gamma.
1854
1855           else
1856             jt=12                    !first parton=q(q~)
1857             nqc(1)=ncc(jq,3-jm)                 !here nqc(1) is gluon.
1858             njc1=ncc(1,jm)
1859             njc2=ncc(3-jq,3-jm)
1860             if(iqj(njc1).ne.0)then
1861               ncj(1,njc1)=njc2
1862             else
1863               ncj(3-jq,njc1)=njc2
1864             endif
1865             if(iqj(njc2).ne.0)then
1866               ncj(1,njc2)=njc1
1867             else
1868               ncj(jq,njc2)=njc1
1869             endif
1870             iqc(3-jm)=10
1871           endif
1872        endif
1873
1874        elseif(iqc(1)*iqc(2).gt.0)then
1875         jt=7                        !qq->qq (q~q~->q~q~)
1876         if(iqc(1).gt.0)then
1877           jq=1
1878         else
1879           jq=2
1880         endif
1881         nqc(1)=ncc(1,3-jm)
1882
1883       else                          ! qq~ ->
1884         if(iqc(jm).gt.0)then
1885           jq=1
1886         else
1887           jq=2
1888         endif
1889         aks=rangen()
1890         if(aks.lt.(wsub(1)+wsub(2))/fborn)then
1891           jt=8                     ! qq~->qq~ (anticolor-color annihilation)
1892           if(aks.gt.wsub(1)/fborn)then
1893             iqa=iabs(iqc(jm))
1894             iq=int((naflav-1)*rangen())+1
1895             if(iq.eq.iqa)iq=naflav
1896             iqc(jm)=iq*iqc(jm)/iqa
1897             iqc(3-jm)=-iqc(jm)
1898           endif
1899           nqc(1)=0
1900           njc1=ncc(1,jm)
1901           njc2=ncc(1,3-jm)
1902           if(iqj(njc1).ne.0)then
1903             ncj(1,njc1)=njc2
1904           else
1905             ncj(3-jq,njc1)=njc2
1906           endif
1907           if(iqj(njc2).ne.0)then
1908             ncj(1,njc2)=njc1
1909           else
1910             ncj(jq,njc2)=njc1
1911           endif
1912         elseif(aks.lt.(wsub(1)+wsub(2)+wsub(3))/fborn)then
1913           jt=10                    !color transfer  qq~->gg
1914           iqc(1)=0
1915           iqc(2)=0
1916           nqc(1)=ncc(1,jm)
1917           nqc(2)=0
1918         elseif(aks.lt.(wsub(1)+wsub(2)+wsub(3)+wsub(4))/fborn)then
1919           jt=13                   ! qq~->g+gamma
1920           nqc(1)=ncc(1,jm)
1921           nqc(2)=ncc(1,3-jm)
1922           iqc(jm)=0
1923           iqc(3-jm)=10
1924         else
1925           jt=14                  ! qq~->gamma+gamma
1926           njc1=ncc(1,jm)
1927           njc2=ncc(1,3-jm)
1928           if(iqj(njc1).ne.0)then
1929             ncj(1,njc1)=njc2
1930           else
1931             ncj(jq,njc1)=njc2
1932           endif
1933           if(iqj(njc2).ne.0)then
1934             ncj(1,njc2)=njc1
1935           else
1936             ncj(3-jq,njc2)=njc1
1937           endif
1938           iqc(jm)=10
1939           iqc(3-jm)=10
1940         endif
1941       endif
1942
1943       if(jt.ne.8.and.jt.ne.9)then
1944         jq2=jq
1945       else
1946         jq2=3-jq
1947       endif
1948
1949       call psdeftr(si+q2mass,ept,ey)    !lorentz boost to c.m. frame
1950
1951       qt=sqrt(qt2)                      !p_t
1952       call pscs(bcos,bsin)              !cos and sin of the polar angle
1953       if(iabs(iqc(1)).ne.4)then
1954 clight cone momentum share for the first parton
1955         z=sngl(psutz(dble(si),dble(qt2),dble(qt2)))
1956         if((jt.eq.11.and.jm.eq.1).or.(jt.eq.12.and.jm.eq.2)
1957      $   .or.(jt.eq.13.and.jm.eq.2))z=1-z
1958         wp3=z*sqrt(si)
1959         wm3=qt2/wp3
1960       elseif(jm.eq.1)then
1961         z=sngl(psutz(dble(si),dble(qt2+q2mass),dble(qt2)))
1962         wp3=z*sqrt(si)
1963         wm3=(qt2+q2mass)/wp3
1964       else
1965         z=sngl(psutz(dble(si),dble(qt2),dble(qt2+dble(q2mass))))
1966         wp3=z*sqrt(si)
1967         wm3=qt2/wp3
1968       endif
1969       ep3(1)=.5*(wp3+wm3)               !parton 4-momentum
1970       ep3(2)=.5*(wp3-wm3)
1971       ep3(3)=qt*bcos
1972       ep3(4)=qt*bsin
1973       call psdefrot(ep3,s0xh,c0xh,s0h,c0h)   !spacial rotation to z-axis
1974       if(iqc(jm).eq.0)then
1975         iq2ini1=9
1976       else
1977         iq2ini1=iqc(jm)
1978       endif
1979       if(iqc(3-jm).eq.0)then
1980         iq2ini2=9
1981       else
1982         iq2ini2=iqc(3-jm)
1983       endif
1984       if(jt.le.10)then
1985         qq1=qt2
1986         qq2=qt2
1987       elseif(jt.le.13)then
1988         qq1=qt2
1989         qq2=0
1990       else
1991         qq1=0
1992         qq2=0
1993       endif
1994
1995       call timsh2(qq1,qq2,sqrt(si),iq2ini1,iq2ini2)  !final state cascade
1996
1997       if(jt.le.10)then      !color connection for the 2nd parton
1998         if(ish.ge.6)write(ifch,*)'jt,jq,nqc:',jt,jq,nqc
1999         call psreti(nqc,jq,1,ey,s0xh,c0xh,s0h,c0h) !color conn. reconstruction
2000         if(jt.eq.1)then
2001           nqc(1)=nqc(2)
2002           nqc(2)=ncc(3-jq,3-jm)
2003         elseif(jt.eq.2)then
2004           nqc(2)=ncc(3-jq,jm)
2005           nqc(1)=ncc(jq,3-jm)
2006         elseif(jt.eq.3)then
2007           nqc(1)=nqc(2)
2008         elseif(jt.eq.4)then
2009           nqc(2)=nqc(1)
2010           nqc(1)=ncc(jq,3-jm)
2011         elseif(jt.eq.5)then
2012           nqc(1)=ncc(jq,jm)
2013         elseif(jt.eq.6)then
2014           nqc(2)=ncc(3-jq,3-jm)
2015           nqc(1)=ncc(1,jm)
2016         elseif(jt.eq.7)then
2017           nqc(1)=ncc(1,jm)
2018         elseif(jt.eq.9)then
2019           nqc(1)=ncc(3-jq,3-jm)
2020         elseif(jt.eq.10)then
2021           nqc(1)=nqc(2)
2022           nqc(2)=ncc(1,3-jm)
2023         endif
2024         if(ish.ge.6)write(ifch,*)'jt,jq2,nqc:',jt,jq2,nqc
2025         call psreti(nqc,jq2,2,ey,s0xh,c0xh,s0h,c0h) !color conn. reconstr.
2026       elseif(jt.le.13)then
2027         if(ish.ge.6)write(ifch,*)'jt,jq,nqc:',jt,jq,nqc
2028         call psreti(nqc,jq,1,ey,s0xh,c0xh,s0h,c0h) !color conn. reconstruction
2029         ep3(1)=pprt(4,2)
2030         ep3(2)=pprt(3,2)
2031         ep3(3)=pprt(1,2)
2032         ep3(4)=pprt(2,2)
2033         call psrotat(ep3,s0xh,c0xh,s0h,c0h)  !special rotation for photon.
2034         call pstrans(ep3,ey,1)
2035         nptl=nptl+1
2036         pptl(1,nptl)=ep3(3)
2037         pptl(2,nptl)=ep3(4)
2038         pptl(3,nptl)=ep3(2)
2039         pptl(4,nptl)=ep3(1)
2040         pptl(5,nptl)=0
2041         idptl(nptl)=10
2042         iorptl(nptl)=iptl
2043         istptl(nptl)=0
2044         jorptl(nptl)=0
2045         do i=1,4
2046           xorptl(i,nptl)=coordo(i)
2047         enddo
2048         tivptl(1,nptl)=coordo(5)
2049         tivptl(2,nptl)=coordo(6)
2050         ityptl(nptl)=71
2051         ifrptl(1,nptl)=0
2052         ifrptl(2,nptl)=0
2053       else
2054         do j=1,2
2055           ep3(1)=pprt(4,j)
2056           ep3(2)=pprt(3,j)
2057           ep3(3)=pprt(1,j)
2058           ep3(4)=pprt(2,j)
2059           call psrotat(ep3,s0xh,c0xh,s0h,c0h)  !special rotation for photon.
2060           call pstrans(ep3,ey,1)
2061           nptl=nptl+1
2062           pptl(1,nptl)=ep3(3)
2063           pptl(2,nptl)=ep3(4)
2064           pptl(3,nptl)=ep3(2)
2065           pptl(4,nptl)=ep3(1)
2066           pptl(5,nptl)=0
2067           idptl(nptl)=10
2068           iorptl(nptl)=iptl
2069           istptl(nptl)=0
2070           jorptl(nptl)=0
2071           do i=1,4
2072             xorptl(i,nptl)=coordo(i)
2073           enddo
2074           tivptl(1,nptl)=coordo(5)
2075           tivptl(2,nptl)=coordo(6)
2076           ityptl(nptl)=72
2077           ifrptl(1,nptl)=0
2078           ifrptl(2,nptl)=0
2079         enddo
2080       endif
2081       call utprix('psabor',ish,ishini,5)
2082       return
2083       end
2084
2085 c------------------------------------------------------------------------
2086       subroutine psreti(nqc,jort,nfprt,ey,s0xh,c0xh,s0h,c0h)
2087 c-----------------------------------------------------------------------
2088 c jet reconstructuring procedure - 4-momenta for all final jets
2089 c nqc(i) - colour connections for the jet
2090 c jort - color orientation for gluons (=1 if +color goes first, =-1 otherwise)
2091
2092       parameter (ntim=1000)
2093       common/cprt/nprtj,pprt(5,ntim),idprt(ntim),iorprt(ntim)
2094      &,idaprt(2,ntim)
2095 c nprtj - number of partons in the jet (including virtual ones)
2096 c pprt - 5-momenta for the partons
2097 c idprt - parton id
2098 c iorprt - parent parton position in the list
2099 c idaprt - daughter partons positions in the list
2100
2101 c output:
2102       parameter (mjstr=20000)
2103       common /psar29/ eqj(4,mjstr),iqj(mjstr),ncj(2,mjstr),nj
2104 c nj - number of final jets
2105 c eqj(i,j) - 4-momentum for the final jet j
2106 c iqj(j) - flavour for the final jet j
2107 c ncj(m,j) - colour connections for the final jet j
2108 c-----------------------------------------------------------------------
2109       dimension ep3(4),nqc(2),ncc(2,ntim),ey(3)
2110       include 'epos.inc'
2111       include 'epos.incsem'
2112       common/cprtx/nprtjx,pprtx(5,2)
2113
2114       if(ish.ge.6)then
2115         write (ifch,*)'nprtj',nprtj
2116         do i=1,nprtj
2117           write (ifch,*)'i,ic,np,ndd',i,idprt(i),iorprt(i),
2118      *    idaprt(1,i),idaprt(2,i)
2119         enddo
2120       endif
2121
2122       ncc(1,nfprt)=nqc(1)
2123       if(idprt(nfprt).eq.9)ncc(2,nfprt)=nqc(2)
2124       iprt=nfprt
2125
2126       if(nprtjx.eq.2)then !out Born before timelike cascade
2127        ep3(1)=pprtx(4,iprt)
2128        ep3(2)=pprtx(3,iprt)
2129        ep3(3)=pprtx(1,iprt)
2130        ep3(4)=pprtx(2,iprt)
2131        call psrotat(ep3,s0xh,c0xh,s0h,c0h)
2132        call pstrans(ep3,ey,1)
2133        pprtx(4,iprt)=ep3(1)
2134        pprtx(3,iprt)=ep3(2)
2135        pprtx(1,iprt)=ep3(3)
2136        pprtx(2,iprt)=ep3(4)
2137       endif
2138
2139 1     continue
2140
2141       idau1=idaprt(1,iprt)
2142       idau2=idaprt(2,iprt)
2143       icp=idprt(iprt)
2144
2145       if(ish.ge.6)then
2146         write (ifch,*)'1-iprt,icp,idau1,idau2',iprt,icp,idau1,idau2,
2147      *  ncc(1,iprt)
2148         if(icp.eq.9)write (ifch,*)'ncc2',ncc(2,iprt)
2149       endif
2150
2151       if(idau1.ne.0.)then         !virtual parton
2152         icd1=idprt(idau1)
2153
2154         if(icp.eq.9)then
2155           if(icd1.ne.9)then      !g -> qq~
2156             ncc(1,idau1)=ncc(jort,iprt)
2157             ncc(1,idau2)=ncc(3-jort,iprt)
2158           else                    !g -> gg
2159             ncc(1,idau1)=ncc(1,iprt)
2160             ncc(2,idau1)=0
2161             ncc(2,idau2)=ncc(2,iprt)
2162             ncc(1,idau2)=0
2163           endif
2164         else                      !q -> qg
2165           ncc(1,idau1)=0
2166           if(icp*(3-2*jort).gt.0)then
2167             ncc(1,idau2)=ncc(1,iprt)
2168             ncc(2,idau2)=0
2169           else
2170             ncc(1,idau2)=0
2171             ncc(2,idau2)=ncc(1,iprt)
2172           endif
2173         endif
2174         iprt=idau1
2175         goto 1
2176       else
2177
2178         nj=nj+1
2179         ep3(1)=pprt(4,iprt)
2180         ep3(2)=pprt(3,iprt)
2181         ep3(3)=pprt(1,iprt)
2182         ep3(4)=pprt(2,iprt)
2183         call psrotat(ep3,s0xh,c0xh,s0h,c0h)
2184         call pstrans(ep3,ey,1)
2185         do i=1,4
2186           eqj(i,nj)=ep3(i)
2187         enddo
2188
2189         if(icp.eq.9)then
2190           iqj(nj)=0
2191         elseif(iabs(icp).lt.3)then
2192           iqj(nj)=icp
2193         elseif(iabs(icp).eq.3)then
2194           iqj(nj)=icp*4/3
2195         else
2196           iqj(nj)=icp*10
2197         endif
2198
2199         if(iqj(nj).ne.0)then
2200           njc=ncc(1,iprt)
2201           if(njc.ne.0)then
2202             ncj(1,nj)=njc
2203             iqc=iqj(njc)
2204             if(iqc.ne.0)then
2205               ncj(1,njc)=nj
2206             else
2207               if(iqj(nj).gt.0)then
2208                 ncj(2,njc)=nj
2209               else
2210                 ncj(1,njc)=nj
2211               endif
2212             endif
2213           else
2214             ncc(1,iprt)=nj
2215           endif
2216         else
2217
2218           do m=1,2
2219             if(jort.eq.1)then
2220               m1=m
2221             else
2222               m1=3-m
2223             endif
2224             njc=ncc(m1,iprt)
2225             if(njc.ne.0)then
2226               ncj(m,nj)=njc
2227               iqc=iqj(njc)
2228               if(iqc.ne.0)then
2229                 ncj(1,njc)=nj
2230               else
2231                 ncj(3-m,njc)=nj
2232               endif
2233             else
2234               ncc(m1,iprt)=nj
2235             endif
2236           enddo
2237         endif
2238         if(ish.ge.6)then
2239           write (ifch,*)'jet-nj,iprt,icp,iqj(nj),ncj',
2240      *    nj,iprt,icp,iqj(nj),ncj(1,nj)
2241           if(icp.eq.9)write (ifch,*)'ncj2',ncj(2,nj)
2242         endif
2243       endif
2244
2245 2     continue
2246       if(iprt.ne.nfprt)then
2247         icp=idprt(iprt)
2248         ipar=iorprt(iprt)
2249         idau1=idaprt(1,ipar)
2250         idau2=idaprt(2,ipar)
2251         if(ish.ge.6)then
2252           write (ifch,*)'2-iprt,icp,idau1,idau2,ipar',
2253      *    iprt,icp,idau1,idau2,ipar,ncc(1,iprt)
2254           if(icp.eq.9)write (ifch,*)ncc(2,iprt)
2255         endif
2256
2257         if(idau1.eq.iprt)then
2258           if(icp.eq.9)then                   !g -> gg
2259             ncc(1,ipar)=ncc(1,iprt)
2260             ncc(1,idau2)=ncc(2,iprt)
2261           else
2262             icpar=idprt(ipar)
2263             if(icpar.eq.9)then               !g -> qq~
2264               ncc(jort,ipar)=ncc(1,iprt)
2265             else                             !q -> qg
2266               if(icp*(3-2*jort).gt.0)then
2267                 ncc(2,idau2)=ncc(1,iprt)
2268               else
2269                 ncc(1,idau2)=ncc(1,iprt)
2270               endif
2271             endif
2272           endif
2273           iprt=idau2
2274           goto 1
2275
2276         else
2277           if(icp.eq.9)then
2278             icpar=idprt(ipar)
2279             if(icpar.eq.9)then                !g -> gg
2280               ncc(2,ipar)=ncc(2,iprt)
2281               ncc(2,idau1)=ncc(1,iprt)
2282             else                              !q -> qg
2283               if(icpar*(3-2*jort).gt.0)then
2284                 ncc(1,ipar)=ncc(1,iprt)
2285                 ncc(1,idau1)=ncc(2,iprt)
2286               else
2287                 ncc(1,ipar)=ncc(2,iprt)
2288                 ncc(1,idau1)=ncc(1,iprt)
2289               endif
2290             endif
2291           else
2292             ncc(3-jort,ipar)=ncc(1,iprt)
2293           endif
2294           iprt=ipar
2295           goto 2
2296         endif
2297       else
2298         if(ish.ge.6)write (ifch,*)'3-iprt,ncc',iprt,ncc(1,iprt)
2299         nqc(1)=ncc(1,nfprt)
2300         if(idprt(nfprt).eq.9)nqc(2)=ncc(2,nfprt)
2301       endif
2302       return
2303       end
2304