EPOS/epos167/epos-qsh-165.f

   1 c  reshuffled from sem, sto, sha
   2
   3 c    contains DIS, and unused 3P stuff
   4 c             ---      ---------------
   5
   6
   7
   8
   9 c###########################################################################
  10 c###########################################################################
  11 c###########################################################################
  12 c###########################################################################
  13 c
  14 c                              DIS
  15 c
  16 c###########################################################################
  17 c###########################################################################
  18 c###########################################################################
  19 c###########################################################################
  20
  21
  22
  23
  24
  25
  26
  27
  28
  29
  30 c-----------------------------------------------------------------------
  31       subroutine lepexp(rxbj,rqsq)
  32 c-----------------------------------------------------------------------
  33 c     generates x_bjorken and q**2 according to an experimental
  34 c     distribution ( given in array xq(nxbj,nqsq) ).
  35 c-----------------------------------------------------------------------
  36       parameter (nxbj=10,nqsq=10)
  37       parameter (xbjmin=0.,qsqmin=4.)
  38       parameter (xbjwid=0.025, qsqwid=4.)
  39       dimension xq(nxbj,nqsq),vxq(nxbj*nqsq)
  40       equivalence (xq(1,1),vxq(1))
  41
  42       data (vxq(i),i=1,50)/
  43      &         1304.02,   366.40,    19.84,    10.79,     6.42,
  44      &            4.54,     4.15,     3.38,     2.03,     1.56,
  45      &          241.63,  1637.26,   427.36,   164.51,    73.72,
  46      &           43.07,    20.73,    12.78,     9.34,     5.83,
  47      &            0.01,   724.66,   563.79,   275.08,   176.13,
  48      &          106.44,    85.82,    54.52,    37.12,    28.65,
  49      &            0.01,   202.40,   491.10,   245.13,   157.07,
  50      &          104.43,    61.05,    49.42,    37.84,    26.79,
  51      &            0.01,     3.77,   316.38,   226.92,   133.45,
  52      &           90.30,    63.67,    48.42,    35.73,    28.04/
  53       data (vxq(i),i=51,100)/
  54      &            0.01,     0.01,   153.74,   213.09,   114.14,
  55      &           76.26,    60.02,    43.15,    43.47,    25.60,
  56      &            0.01,     0.01,    39.31,   185.74,   108.56,
  57      &           88.40,    47.29,    39.35,    31.80,    22.91,
  58      &            0.01,     0.01,     0.01,   104.61,   107.01,
  59      &           66.24,    45.34,    37.45,    33.44,    23.78,
  60      &            0.01,     0.01,     0.01,    56.58,    99.39,
  61      &           67.78,    43.28,    35.98,    34.63,    18.31,
  62      &            0.01,     0.01,     0.01,    13.56,    76.25,
  63      &           64.30,    42.80,    28.56,    21.19,    20.75 /
  64
  65       data init/0/
  66       init=init+1
  67       if(init.eq.1) then
  68       n=nxbj*nqsq
  69       sum=0.
  70       do 1 i=1,n
  71       sum=sum+vxq(i)
  72 1     continue
  73       do 2 i=2,n
  74 2     vxq(i)=vxq(i)+vxq(i-1)
  75       do 3 i=1,n
  76 3     vxq(i)=vxq(i)/sum
  77       endif
  78
  79       n=nxbj*nqsq
  80       r=rangen()
  81       call utloc(vxq,n,r,iloc)
  82       if(iloc.ge.n) iloc=iloc-1
  83       i=mod(iloc,nxbj)+1
  84       if(i.eq.0) i=nxbj
  85       j=iloc/nxbj + 1
  86       dxint=vxq(1)
  87       if(iloc.gt.0) dxint=vxq(iloc+1)-vxq(iloc)
  88       dxbj=xbjwid*abs(r-vxq(iloc+1))/dxint
  89       dy=qsqwid*rangen()
  90       rxbj=xbjmin+xbjwid*float(i-1)+dxbj
  91       rqsq=qsqmin+qsqwid*float(j-1)+dy
  92       return
  93       end
  94
  95 c-----------------------------------------------------------------------
  96       subroutine fremny(wp1,wm1,pnx,pny,sm,ic1,ic2,ic3,ic4,coord,ey0)
  97 c-----------------------------------------------------------------------
  98 c  treats remnant from deep inelastic process;
  99 c-----------------------------------------------------------------------
 100       include 'epos.inc'
 101       include 'epos.incsem'
 102       dimension coord(6),ic(2),ep(4),ey(3),ey0(3),ep3(4)
 103       double precision  ept(4),ept1(4)
 104
 105       call utpri('fremny',ish,ishini,5)
 106       if(ish.ge.5)write(ifch,*)'writing remnant'
 107       if(ish.ge.5)write(ifch,*)
 108      *'wp1,wm1,pnx,pny,sm,ic1,ic2,ic3,ic4,coord,ey0:'
 109       if(ish.ge.5)write(ifch,*)
 110      *wp1,wm1,pnx,pny,sm,ic1,ic2,ic3,ic4,coord,ey0
 111
 112         if(ic3.eq.0.and.ic4.eq.0)then
 113
 114       ic(1)=ic1
 115       ic(2)=ic2
 116       nptl=nptl+1
 117       ep3(3)=pnx
 118       ep3(4)=pny
 119       ep3(2)=(wp1-wm1)/2
 120       ep3(1)=(wp1+wm1)/2
 121       call pstrans(ep3,ey0,-1)
 122       pptl(1,nptl)=ep3(3)
 123       pptl(2,nptl)=ep3(4)
 124       pptl(3,nptl)=ep3(2)
 125       pptl(4,nptl)=ep3(1)
 126       pptl(5,nptl)=sqrt(sm)
 127       idptl(nptl)=idtra(ic,0,0,3)
 128       iorptl(nptl)=1
 129       istptl(nptl)=0
 130       jorptl(nptl)=0
 131       do i=1,4
 132       xorptl(i,nptl)=coord(i)
 133       enddo
 134       tivptl(1,nptl)=coord(5)
 135       tivptl(2,nptl)=coord(6)
 136       ityptl(nptl)=40
 137
 138        if(ish.ge.7)then
 139       write(ifch,*)'proj: nptl, mass**2,id',nptl,sm,idptl(nptl)
 140       write(ifch,*)'ept',
 141      *pptl(1,nptl),pptl(2,nptl),pptl(3,nptl),pptl(4,nptl)
 142        endif
 143
 144         else
 145
 146       ic(1)=ic1
 147       ic(2)=ic2
 148       nptl=nptl+1
 149       idptl(nptl)=idtra(ic,0,0,3)
 150       istptl(nptl)=20
 151       iorptl(nptl)=1
 152       jorptl(nptl)=0
 153       do i=1,4
 154       xorptl(i,nptl)=coord(i)
 155       enddo
 156       tivptl(1,nptl)=coord(5)
 157       tivptl(2,nptl)=coord(6)
 158       ityptl(nptl)=40
 159
 160       ic(1)=ic3
 161       ic(2)=ic4
 162       idptl(nptl+1)=idtra(ic,0,0,3)
 163       istptl(nptl+1)=20
 164       iorptl(nptl+1)=1
 165       jorptl(nptl+1)=0
 166       do i=1,4
 167       xorptl(i,nptl+1)=coord(i)
 168       enddo
 169       tivptl(1,nptl+1)=coord(5)
 170       tivptl(2,nptl+1)=coord(6)
 171       ityptl(nptl+1)=40
 172
 173       ep3(3)=pnx
 174       ep3(4)=pny
 175       ep3(2)=(wp1-wm1)/2
 176       ep3(1)=(wp1+wm1)/2
 177       call pstrans(ep3,ey0,-1)   !boost to hadronic c.m.s.
 178       ept(1)=ep3(3)
 179       ept(2)=ep3(4)
 180       ept(3)=ep3(2)
 181       ept(4)=ep3(1)
 182       do i=1,4
 183         ept1(i)=ep3(i)
 184       enddo
 185
 186       sww=sqrt(sm)
 187       call psdeftr(sm,ept1,ey)
 188       ep(1)=.5*sww
 189       ep(2)=.5*sww
 190       ep(3)=0.
 191       ep(4)=0.
 192       call pstrans(ep,ey,1)
 193       pptl(1,nptl)=ep(3)
 194       pptl(2,nptl)=ep(4)
 195       pptl(3,nptl)=ep(2)
 196       pptl(4,nptl)=ep(1)
 197       do i=1,4
 198         pptl(i,nptl+1)=ept(i)-pptl(i,nptl)
 199       enddo
 200       nptl=nptl+1
 201         endif
 202
 203  1000 continue
 204       if(ish.ge.5)write(ifch,*)'fremny: final nptl',nptl
 205       call utprix('fremny',ish,ishini,5)
 206       return
 207       end
 208
 209 c-----------------------------------------------------------------------
 210       subroutine psadis(iret)
 211 c-----------------------------------------------------------------------
 212 c psadis - DIS interaction
 213 c-----------------------------------------------------------------------
 214       double precision ept(4),ept1(4),xx,wpt(2),eprt,pl,plprt,psutz
 215      *,psuds
 216       dimension ep3(4),ey(3),ey0(3),bx(6),
 217      *qmin(2),iqc(2),nqc(2),ncc(2,2),gdv(2),gds(2),dfp(4)
 218       parameter (mjstr=20000)
 219       common /psar29/ eqj(4,mjstr),iqj(mjstr),ncj(2,mjstr),nj
 220       common /psar30/ iorj(mjstr),ityj(mjstr),bxj(6,mjstr)
 221       double precision pgampr,rgampr
 222       common/cgampr/pgampr(5),rgampr(4)
 223       parameter (ntim=1000)
 224       common/cprt/nprtj,pprt(5,ntim),idprt(ntim),iorprt(ntim)
 225      &,idaprt(2,ntim)
 226       common/ciptl/iptl
 227       include 'epos.inc'
 228       include 'epos.incsem'
 229
 230       call utpri('psadis',ish,ishini,3)
 231       if(ish.ge.3)write (ifch,*)'engy,elepti,iolept:'
 232       if(ish.ge.3)write (ifch,*)engy,elepti,iolept
 233       nptl=nptl+1
 234       idptl(nptl)=1220
 235       istptl(nptl)=1
 236       nptlh=nptl
 237       iptl=nptl
 238       s00=1.
 239
 240       pptl(1,nptl)=0.
 241       pptl(2,nptl)=0.
 242       pptl(3,nptl)=-engy/2
 243       pptl(4,nptl)=engy/2
 244       pptl(5,nptl)=0
 245
 246 1     continue
 247       if(iolept.eq.1)then
 248         wtot=engy**2
 249         engypr=wtot/4./elepti
 250         gdv01=psdh(ydmax*wtot,qdmin,iclpro,0)
 251         gdv02=psdh(ydmax*wtot,qdmin,iclpro,1)
 252         gds01=psdsh(ydmax*wtot,qdmin,iclpro,dqsh,0)
 253         gds02=psdsh(ydmax*wtot,qdmin,iclpro,dqsh1,1)
 254         gb0=(1.+(1.-ydmax)**2)*(gdv01+gds01)
 255      *  +2.*(1.-ydmax)*(gdv02+gds02)
 256
 257 2       continue
 258         qq=qdmin*(qdmax/qdmin)**rangen()
 259         yd=ydmin*(ydmax/ydmin)**rangen()
 260         wd=yd*wtot
 261         if(ish.ge.4)write (ifch,*)'qq,wd,yd,ydmin,ydmax:'
 262         if(ish.ge.4)write (ifch,*)qq,wd,yd,ydmin,ydmax
 263         if(wd.lt.qq)goto 2
 264
 265         gdv(1)=psdh(wd,qq,iclpro,0)
 266         gdv(2)=psdh(wd,qq,iclpro,1)
 267         gds(1)=psdsh(wd,qq,iclpro,dqsh,0)
 268         gds(2)=psdsh(wd,qq,iclpro,dqsh1,1)
 269         gbtr=(1.+(1.-yd)**2)*(gdv(1)+gds(1))
 270         gblong=2.*(1.-yd)*(gdv(2)+gds(2))
 271 c        gblong=0.   !???????
 272         gb=(gbtr+gblong)/gb0*.7
 273         if(ish.ge.4)then
 274           if(gb.gt.1.)write(ifmt,*)'gb,qq,yd,wd',gb,qq,yd,wd
 275           write (ifch,*)'gb,gdv,gds,gdv0,gds0,yd:'
 276           write (ifch,*)gb,gdv,gds,gdv01,gds01,
 277      *    gdv02,gds02,yd
 278         endif
 279         if(rangen().gt.gb)goto 2
 280
 281         long=int(rangen()+gblong/(gbtr+gblong))
 282         elepto=qq/elepti/4.+elepti*(1.-yd)
 283         costhet=1.-qq/elepti/elepto/2.
 284         theta=acos(costhet)
 285         if(theta/pi*180..lt.themin)goto 2
 286         if(theta/pi*180..gt.themax)goto 2
 287         if(elepto.lt.elomin)goto 2
 288         if(ish.ge.3)write (ifch,*)'theta,elepto,elepti,iclpro:'
 289         if(ish.ge.3)write (ifch,*)theta/pi*180.,elepto,elepti,iclpro
 290         xbjevt=qq/wd
 291         qsqevt=qq
 292
 293         call pscs(bcos,bsin)
 294         rgampr(1)=-elepto*sin(theta)*bcos
 295         rgampr(2)=-elepto*sin(theta)*bsin
 296         rgampr(3)=elepti-elepto*costhet
 297         rgampr(4)=elepti-elepto
 298
 299         pgampr(1)=rgampr(1)
 300         pgampr(2)=rgampr(2)
 301         pgampr(3)=rgampr(3)-engypr
 302         pgampr(4)=rgampr(4)+engypr
 303         sm2=pgampr(4)*pgampr(4)
 304      *  -pgampr(1)*pgampr(1)-pgampr(2)*pgampr(2)-pgampr(3)*pgampr(3)
 305         pgampr(5)=sqrt(sm2)
 306         call utlob2(1,pgampr(1),pgampr(2),pgampr(3),pgampr(4),pgampr(5)
 307      *  ,rgampr(1),rgampr(2),rgampr(3),rgampr(4),40)
 308         if(ish.ge.4)write (ifch,*)'rgampr:',rgampr
 309
 310       elseif(iolept.lt.0)then
 311 21      call lepexp(xbjevt,qsq)
 312         qq=qsq
 313         wd=qq/xbjevt
 314         if(qq.lt.qdmin.or.qq.gt.qdmax)goto21
 315
 316         gdv(1)=psdh(wd,qq,iclpro,0)
 317         gdv(2)=psdh(wd,qq,iclpro,1)
 318         gds(1)=psdsh(wd,qq,iclpro,dqsh,0)
 319         gds(2)=psdsh(wd,qq,iclpro,dqsh1,1)
 320         yd=wd/engy**2
 321         gbtr=(1.+(1.-yd)**2)*(gdv(1)+gds(1))
 322         gblong=2.*(1.-yd)*(gdv(2)+gds(2))
 323         gblong=0. !????????????
 324         long=int(rangen()+gblong/(gbtr+gblong))
 325       else
 326         stop'wrong iolept'
 327       endif
 328       if(ish.ge.3)write (ifch,*)'qq,xbj,wd,gdv,gds,dqsh:'
 329       if(ish.ge.3)write (ifch,*)qq,xbjevt,wd,gdv,gds,dqsh
 330
 331       egyevt=sqrt(wd-qq)
 332       pmxevt=.5*egyevt
 333
 334       wp0=sqrt(qq)       !breit frame
 335       wm0=(wd-qq)/wp0
 336       ey0(1)=egyevt/wp0  !boost to the hadronic c.m.s.
 337       ey0(2)=1.
 338       ey0(3)=1.
 339       do i=1,6
 340         bx(i)=0.
 341       enddo
 342
 343       if(long.eq.0)then
 344         sdmin=qq/(1.-sqrt(q2ini/qq))
 345         sqmin=sdmin
 346       else
 347         sdmin=4.*max(q2min,qcmass**2)+qq  !minimal mass for born
 348         xmm=(5.*sdmin-qq)/4.
 349         sqmin=1.1*(xmm+sqrt(xmm**2-qq*(sdmin-qq-4.*q2ini)))
 350      *  /2./(1.-4.*q2ini/(sdmin-qq))
 351       endif
 352       if(long.eq.1.and.wd.lt.1.001*sdmin)goto 1
 353
 354       proja=210000.
 355       projb=0.
 356       call fremnu(ammin,proja,projb,proja,projb,
 357      *icp1,icp2,icp3,icp4)
 358
 359       nj=0
 360
 361       if((rangen().lt.gdv(long+1)/(gdv(long+1)+gds(long+1)).or.
 362      *egyevt.lt.1.0001*(ammin+sqrt(sdmin-qq))).and.
 363      *(long.eq.0.or.wd.gt.sqmin))then
 364         if(long.eq.0)then
 365           xd=qq/wd
 366           tu=psdfh4(xd,q2min,0.,iclpro,1)/2.25
 367           td=psdfh4(xd,q2min,0.,iclpro,2)/9.
 368           gdv0=(tu+td)*4.*pi**2*alfe/qq
 369      *    *sngl(psuds(qq,1)/psuds(q2min,1))
 370           if(ish.ge.4)write (ifch,*)'gdv0:',gdv0,sdmin
 371
 372           if(rangen().lt.gdv0/gdv(1).or.wd.le.1.0001*sdmin)then    !?????
 373             if(ish.ge.3)write (ifch,*)'no cascade,gdv0,gdv',gdv0,gdv
 374             if(rangen().lt.tu/(tu+td))then
 375               iq1=1
 376               izh=3
 377             else
 378               iq1=2
 379               izh=6
 380             endif
 381             jq=1
 382             if(ish.ge.8)write (ifch,*)'before call timsh2: ',
 383      *      'qq,egyevt,iq1',qq,egyevt,iq1
 384             call timsh2(qq,0.,egyevt,iq1,-iq1)
 385
 386             nj=nj+1
 387             iqj(nj)=izh
 388             nqc(1)=nj
 389             nqc(2)=0
 390
 391             ep3(1)=pprt(4,2)
 392             ep3(2)=pprt(3,2)
 393             ep3(3)=0.
 394             ep3(4)=0.
 395             call pstrans(ep3,ey0,1)
 396             do i=1,4
 397               eqj(i,nj)=ep3(i)
 398             enddo
 399             s0h=0.
 400             c0h=1.
 401             s0xh=0.
 402             c0xh=1.
 403             call psreti(nqc,jq,1,ey0,s0xh,c0xh,s0h,c0h)
 404             goto 17
 405           endif
 406         endif
 407
 408         call psdint(wd,qq,sds0,sdn0,sdb0,sdt0,sdr0,1,long)
 409         if(ish.ge.3)write (ifch,*)'wd,qq,sds0,sdn0,sdb0,sdt0,sdr0:'
 410         if(ish.ge.3)write (ifch,*)wd,qq,sds0,sdn0,sdb0,sdt0
 411         gb10=(sdn0+sdt0)*(1.-qq/wd)
 412         xdmin=sqmin/wd
 413
 414 3       continue
 415         xd=(xdmin-qq/wd)/((xdmin-qq/wd)/(1.-qq/wd))
 416      *  **rangen()+qq/wd
 417         call psdint(xd*wd,qq,sds,sdn,sdb,sdt,sdr,1,long)
 418         if(ish.ge.3)write (ifch,*)'wdhard,qq,sds,sdn,sdb,sdt:'
 419         if(ish.ge.3)write (ifch,*)xd*wd,qq,sds,sdn,sdb,sdt
 420         tu=psdfh4(xd,q2min,0.,iclpro,1)
 421         td=psdfh4(xd,q2min,0.,iclpro,2)
 422         gb1=(sdn*(tu/2.25+td/9.)+sdt*(tu+td)/4.5)
 423      *  *(1.-qq/wd/xd)/gb10
 424         if(gb1.gt.1..and.ish.ge.1)write(ifmt,*)'gb1,xd,wd,qq,sdt0,sdt',
 425      *   gb1,xd,wd,qq,sdt0,sdt
 426         if(ish.ge.6)write (ifch,*)'gb1,xd,wd,qq,sdt0,sdt:'
 427         if(ish.ge.6)write (ifch,*)gb1,xd,wd,qq,sdt0,sdt
 428         if(rangen().gt.gb1)goto 3
 429
 430         gdres=(sdt-sds)/4.5
 431         gdrga=sdr/4.5
 432         gdsin=sds/4.5
 433         dtu=tu*(sdn/2.25+sdt/4.5)
 434         dtd=td*(sdn/9.+sdt/4.5)
 435         if(rangen().lt.dtu/(dtu+dtd))then
 436           iq1=1
 437           izh=3
 438           gdbor=sdb/2.25
 439           gdnon=sdn/2.25
 440         else
 441           iq1=2
 442           izh=6
 443           gdbor=sdb/9.
 444           gdnon=sdn/9.
 445         endif
 446
 447         wpi=wp0
 448         wmi=(xd*wd-qq)/wpi
 449         iqc(2)=iq1
 450         nj=nj+1
 451         iqj(nj)=izh
 452         eqj(1,nj)=.5*(wm0-wmi)
 453         eqj(2,nj)=-eqj(1,nj)
 454         eqj(3,nj)=0.
 455         eqj(4,nj)=0.
 456         ncc(1,2)=nj
 457         ncc(2,2)=0
 458         if(ish.ge.3)write (ifch,*)'wp0,wm0,wpi,wmi,iqc(2),eqj'
 459         if(ish.ge.3)write (ifch,*)wp0,wm0,wpi,wmi,iqc(2),eqj(2,nj)
 460
 461       else
 462         xdmin=sdmin/wd
 463         xpmax=((egyevt-ammin)**2+qq)/wd
 464         iq1=int(3.*rangen()+1.)*(2.*int(.5+rangen())-1.)
 465
 466         aks=rangen()
 467         if(long.eq.0.and.aks.lt.dqsh/gds(1).and.
 468      *  egyevt.gt.ammin+sqrt(s00))then
 469           if(ish.ge.3)write (ifch,*)'no cascade for q_s',
 470      *    aks,dqsh/gds(1)
 471           xd=qq/wd
 472           xpmin=xd+s00/wd
 473           jcasc=0
 474           if(iq1.gt.0)then
 475             jq=1
 476           else
 477             jq=2
 478           endif
 479         else
 480           jcasc=1
 481           call psdint(xpmax*wd,qq,sds0,sdn0,sdb0,sdt0,sdr0,0,long)
 482           call psdint(xpmax*wd,qq,sdsq0,sdnq0,sdbq0,sdtq0,sdrq0,1,long)
 483           if(ish.ge.3)write (ifch,*)
 484      *    'xpmax*wd,qq,sds0,sdn0,sdb0,sdt0,sdr0:'
 485           if(ish.ge.3)write (ifch,*)
 486      *    xpmax*wd,qq,sds0,sdn0,sdb0,sdt0,sdr0
 487         gb10=sdt0*fzeroGluZZ(0.,iclpro)+(sdnq0+sdtq0)
 488      *        *fzeroSeaZZ(0.,iclpro)
 489           gb10=gb10*15.
 490
 491 4         xd=xdmin*(xpmax/xdmin)**rangen()
 492           xpmin=xd
 493           call psdint(xd*wd,qq,sds,sdn,sdb,sdt,sdr,0,long)
 494           call psdint(xd*wd,qq,sdsq,sdnq,sdbq,sdtq,sdrq,1,long)
 495           if(ish.ge.3)write (ifch,*)'xd*wd,qq,sds,sdn,sdb,sdt,sdr:'
 496           if(ish.ge.3)write (ifch,*)xd*wd,qq,sds,sdn,sdb,sdt,sdr
 497           wwg=sdt*fzeroGluZZ(xd,iclpro)
 498           wwq=(sdnq+sdtq)*fzeroSeaZZ(xd,iclpro)
 499           gb12=(wwq+wwg)/gb10*(xpmax/xd)**dels
 500           if(gb12.gt.1..and.ish.ge.1)write(ifmt,*)
 501      *    'gb12,xpmax*wd,xd*wd,sdt0,sdnq0+sdtq0,sdt,sdnq+sdtq',
 502      *    gb12,xpmax*wd,xd*wd,sdt0,sdnq0+sdtq0,sdt,sdnq+sdtq,
 503      *    wwq,wwg,(xpmax/xd)**dels,gb10
 504           if(ish.ge.5)write (ifch,*)'gb12,xd,xpmax,wwq,wwg:'
 505           if(ish.ge.5)write (ifch,*)gb12,xd,xpmax,wwq,wwg
 506           if(rangen().gt.gb12)goto 4
 507         endif
 508
 509         if(jcasc.ne.0)then
 510           gb20=(1.-xd/xpmax)**betpom*sdt*(1.-glusea)+
 511      *    EsoftQZero(xd/xpmax)*(sdnq+sdtq)*glusea
 512         else
 513           gb20=EsoftQZero(xd/xpmax)
 514         endif
 515         if(1.+2.*(-alpqua)+dels.ge.0.)then
 516           xpminl=(1.-xpmax)**(alplea(iclpro)+1.)
 517           xpmaxl=(1.-xpmin)**(alplea(iclpro)+1.)
 518
 519 5         xp=1.-(xpminl+(xpmaxl-xpminl)*rangen())**
 520      *    (1./(alplea(iclpro)+1.))
 521           if(jcasc.ne.0)then
 522             gb2=((1.-xd/xp)**betpom*sdt*(1.-glusea)+
 523      *      EsoftQZero(xd/xp)*(sdnq+sdtq)*glusea)*(xp/xpmax)**
 524      *      (1.+2.*(-alpqua)+dels)/gb20
 525           else
 526            gb2=EsoftQZero(xd/xp)*(xp/xpmax)**(1.+2.*(-alpqua)+dels)/gb20
 527           endif
 528           if(gb2.gt.1..and.ish.ge.1)then
 529             write(ifmt,*)'gb2,xp:',gb2,xp
 530 c            read (*,*)
 531           endif
 532           if(rangen().gt.gb2)goto 5
 533         else
 534           xpmaxl=xpmax**(2.+2.*(-alpqua)+dels)
 535           xpminl=xpmin**(2.+2.*(-alpqua)+dels)
 536
 537 6         xp=(xpminl+(xpmaxl-xpminl)*rangen())**
 538      *    (1./(2.+2.*(-alpqua)+dels))
 539           if(jcasc.ne.0)then
 540             gb21=((1.-xd/xp)**betpom*sdt*(1.-glusea)+
 541      *      EsoftQZero(xd/xp)*(sdnq+sdtq)*glusea)*
 542      *      ((1.-xp)/(1.-xd))**alplea(iclpro)/gb20
 543           else
 544           gb21=EsoftQZero(xd/xp)*((1.-xp)/(1.-xd))**alplea(iclpro)/gb20
 545           endif
 546           if(gb21.gt.1..and.ish.ge.1)then
 547             write(ifmt,*)'gb21,xp:',gb21,xp
 548 c            read (*,*)
 549           endif
 550           if(rangen().gt.gb21)goto 6
 551         endif
 552
 553         wwh=xd*wd-qq
 554         wwsh=xp*wd-qq
 555         ammax=(egyevt-sqrt(wwsh))**2
 556 22      call fremnx(ammax,ammin,sm,icp3,icp4,iret)
 557         if(iret.ne.0.and.ish.ge.1)write(ifmt,*)'iret.ne.0!'
 558      *                                         ,ammax,ammin**2
 559         wmn=(1.-xp)*wd/wp0
 560         wpn=sm/wmn
 561         pnx=0.
 562         pny=0.
 563         wpp=wp0-wpn
 564         wmp=wm0-wmn
 565         if(ish.ge.5)write(ifch,*)'wp0,wm0,wpn,wmn,wpp,wmp:'
 566         if(ish.ge.5)write(ifch,*)wp0,wm0,wpn,wmn,wpp,wmp
 567
 568         if(jcasc.eq.0.or.rangen().lt.wwq/(wwg+wwq).
 569      *  and.xd*wd.gt.sqmin.and.wwsh.gt.
 570      *  (sqrt(wwh)+sqrt(s00))**2)then
 571           zgmin=xd/xp
 572           zgmax=1./(1.+wp0/xd/wd/(wpp-wwh/wmp))
 573           if(zgmin.gt.zgmax)goto 22
 574 23        zg=zgmin-rangen()*(zgmin-zgmax)
 575           if(rangen().gt.zg**dels*((1.-xd/xp/zg)/ (1.-xd/xp))**betpom)
 576      *    goto 23
 577           xg=xd/zg             !w- share for the struck quark
 578           wmq=wd/wp0*(xg-xd)   !w- for its counterpart
 579           wpq=s00/wmq          !1. gev^2 / wmq
 580           wmq=0.
 581           wpp=wpp-wpq
 582           wmp=wmp-wmq
 583           sxx=wpp*wmp
 584           if(ish.ge.5)write (ifch,*)'wpq,wmq,wpp,wmp,sxx:'
 585           if(ish.ge.5)write (ifch,*)wpq,wmq,wpp,wmp,sxx
 586
 587           if(jcasc.eq.0)then
 588             if(ish.ge.6)write (ifch,*)'before call timsh2: qq,sxx,iq1',
 589      *      qq,sxx,iq1
 590             call timsh2(qq,0.,sqrt(sxx),iq1,-iq1)
 591             ept(1)=.5*(wpp+wmp)
 592             ept(2)=.5*(wpp-wmp)
 593             ept(3)=0.
 594             ept(4)=0.
 595             call psdeftr(sxx,ept,ey)
 596             ep3(1)=pprt(4,2)
 597             ep3(2)=pprt(3,2)
 598             ep3(3)=0.
 599             ep3(4)=0.
 600
 601             call pstrans(ep3,ey,1)
 602             wmp=ep3(1)-ep3(2)
 603             goto 24
 604           endif
 605         else
 606           iq1=0
 607           sxx=wpp*wmp
 608         endif
 609
 610         if(ish.ge.3)write (ifch,*)'wwh,wwsh,sxx,wpp,wmp:',
 611      *  wwh,wwsh,sxx,wpp,wmp
 612
 613         wpi=wpp
 614         wmi=wwh/wpp
 615         wmp=wmp-wmi
 616 24      call fremny(wpn,wmn,pnx,pny,sm,icp1,icp2,icp3,icp4,bx,ey0)
 617
 618         if((-alpqua).eq.-1.)stop'dis does not work for 1/x'
 619 25      aks=rangen()
 620         z=.5*aks**(1./(1.+(-alpqua)))
 621         if(z.lt.1.e-5.or.rangen().gt.(2.*(1.-z))**(-alpqua))goto 25
 622         if(rangen().gt..5)z=1.-z
 623         wm2=wmp*z
 624         wm1=wmp-wm2
 625
 626         iqc(2)=iq1
 627         nj=nj+1
 628         iqj(nj)=-int(2.*rangen()+1.)
 629         iqj(nj+1)=-iqj(nj)
 630         eqj(1,nj)=.5*wm1
 631         eqj(2,nj)=-eqj(1,nj)
 632         eqj(3,nj)=0.
 633         eqj(4,nj)=0.
 634         eqj(1,nj+1)=.5*wm2
 635         eqj(2,nj+1)=-eqj(1,nj+1)
 636         eqj(3,nj+1)=0.
 637         eqj(4,nj+1)=0.
 638         nj=nj+1
 639
 640         if(iq1.eq.0)then
 641           ncc(1,2)=nj-1
 642           ncc(2,2)=nj
 643           gdres=sdt-sds
 644           gdrga=sdr
 645           gdsin=sds
 646           gdbor=sdb
 647           gdnon=sdn
 648         else
 649           nj=nj+1
 650           if(iabs(iq1).eq.3)then
 651             iqj(nj)=-iq1*4/3
 652           else
 653             iqj(nj)=-iq1
 654           endif
 655           eqj(1,nj)=.5*(wpq+wmq)
 656           eqj(2,nj)=.5*(wpq-wmq)
 657           eqj(3,nj)=0.
 658           eqj(4,nj)=0.
 659           if(iq1.gt.0)then
 660             ncj(1,nj)=nj-1
 661             ncj(1,nj-1)=nj
 662             ncj(2,nj)=0
 663             ncj(2,nj-1)=0
 664           else
 665             ncj(1,nj)=nj-2
 666             ncj(1,nj-2)=nj
 667             ncj(2,nj)=0
 668             ncj(2,nj-2)=0
 669           endif
 670
 671           if(jcasc.eq.0)then
 672             if(iq1.gt.0)then
 673               nqc(1)=nj-2
 674               nqc(2)=0
 675             else
 676               nqc(1)=nj-1
 677               nqc(2)=0
 678             endif
 679             s0h=0.
 680             c0h=1.
 681             s0xh=0.
 682             c0xh=1.
 683             call psreti(nqc,jq,1,ey,s0xh,c0xh,s0h,c0h)
 684             goto 17
 685           else
 686             gdres=(sdtq-sdsq)/4.5
 687             gdrga=sdrq/4.5
 688             gdsin=sdsq/4.5
 689             gdbor=sdbq/4.5
 690             gdnon=sdnq/4.5
 691             if(iq1.gt.0)then
 692               ncc(1,2)=nj-2
 693               ncc(2,2)=0
 694             else
 695               ncc(1,2)=nj-1
 696               ncc(2,2)=0
 697             endif
 698           endif
 699         endif
 700
 701         if(ish.ge.3)write (ifch,*)'wpn,wmn,wpi,wmi,wm1,wm2,nj'
 702         if(ish.ge.3)write (ifch,*)wpn,wmn,wpi,wmi,wm1,wm2,nj
 703       endif
 704
 705       si=wpi*wmi+qq
 706       qmin(2)=q2min                 !effective momentum cutoff below
 707       s2min=max(4.*qq,16.*q2min)    !mass cutoff for born scattering
 708
 709       if(rangen().gt.gdres/(gdres+gdsin+gdnon).or.
 710      *si.lt.(s2min+qq))goto 12
 711
 712 c---------------------------------------
 713 c hard pomeron (resolved photon)
 714 c---------------------------------------
 715       if(ish.ge.3)write(ifmt,*)'resolved,gdrga,gdres',gdrga,gdres
 716
 717       jj=1
 718       if(rangen().gt.gdrga/gdres.and.si.gt.1.1*s2min+qq)then
 719         if(ish.ge.3)write(ifmt,*)'dir-res,si,qq',si,qq
 720         pt=0.
 721         pt2=0.
 722         iqc(1)=0
 723         ept(1)=.5*(wpi+wmi)
 724         ept(2)=.5*(wpi-wmi)
 725         ept(3)=0.
 726         ept(4)=0.
 727         wpt(1)=wpi           !lc+ for the current jet emission
 728         wpt(2)=si/wpi        !lc- for the current jet emission
 729
 730         qqmin=max(q2min,s2min/(si/qq-1.))
 731         qqmax=min(si/2.,si-s2min)
 732         qmin(1)=qqmin
 733         xmax=1.
 734         xmin=(s2min+qq)/si
 735         if(qqmin.ge.qqmax.or.xmin.ge.xmax)stop'min>max'
 736         gb0=psjti(qmin(1),qq,si-qq,7,iqc(2),1)*psfap(1.d0,0,1)
 737
 738         ncc(1,1)=0
 739         ncc(2,1)=0
 740         jgamma=1
 741         ntry=0
 742         goto 9
 743       else
 744         if(ish.ge.3)write(ifmt,*)'res,si,qq',si,qq
 745         qmin(1)=q2min                 !effective momentum cutoff above
 746         si=si-qq
 747         zmin=s2min/si
 748         dft0=psdfh4(zmin,q2min,qq,0,0)*psjti(q2min,qq,si,0,iqc(2),1)
 749      *  +(psdfh4(zmin,q2min,qq,0,1)+psdfh4(zmin,q2min,qq,0,2)+
 750      *  psdfh4(zmin,q2min,qq,0,3))*psjti(q2min,qq,si,7,iqc(2),1)
 751
 752 7       continue
 753         z=zmin**rangen()
 754         do i=1,4
 755           dfp(i)=psdfh4(z,q2min,qq,0,i-1)
 756         enddo
 757         dfp(1)=dfp(1)*psjti(q2min,qq,z*si,0,iqc(2),1)
 758         dfptot=dfp(1)
 759         if(iqc(2).eq.0)then
 760           sjq=psjti(q2min,qq,z*si,1,0,1)
 761           do i=2,4
 762             dfp(i)=dfp(i)*sjq
 763             dfptot=dfptot+dfp(i)
 764           enddo
 765         else
 766           sjqqp=psjti(q2min,qq,z*si,1,2,1)
 767           do i=2,4
 768             if(iabs(iqc(2)).eq.i-1)then
 769               dfp(i)=dfp(i)*(psjti(q2min,qq,z*si,1,1,1)+
 770      *        psjti(q2min,qq,z*si,1,-1,1))/2.
 771             else
 772               dfp(i)=dfp(i)*sjqqp
 773             endif
 774             dfptot=dfptot+dfp(i)
 775           enddo
 776         endif
 777
 778         if(rangen().gt.dfptot/dft0)goto 7
 779
 780         wpq=wpi*(1.-z)
 781         wpi=wpi*z
 782         aks=dfptot*rangen()
 783         if(aks.lt.dfp(1))then
 784           iqc(1)=0
 785           nj=nj+1
 786           ncc(1,1)=nj
 787           ncc(2,1)=nj+1
 788
 789           iqj(nj)=-int(2.*rangen()+1.)
 790           iqj(nj+1)=-iqj(nj)
 791           wpq1=wpq*rangen()
 792           eqj(1,nj)=.5*wpq1
 793           eqj(2,nj)=eqj(1,nj)
 794           eqj(3,nj)=0.
 795           eqj(4,nj)=0.
 796           eqj(1,nj+1)=.5*(wpq-wpq1)
 797           eqj(2,nj+1)=eqj(1,nj+1)
 798           eqj(3,nj+1)=0.
 799           eqj(4,nj+1)=0.
 800           nj=nj+1
 801
 802         else
 803           if(aks.lt.dfp(1)+dfp(2))then
 804             iqc(1)=1
 805           elseif(aks.lt.dfp(1)+dfp(2)+dfp(3))then
 806             iqc(1)=2
 807           else
 808             iqc(1)=3
 809           endif
 810           iqc(1)=iqc(1)*(2*int(2.*rangen())-1)
 811           nj=nj+1
 812           ncc(1,1)=nj
 813           ncc(2,1)=0
 814
 815           iqj(nj)=-iqc(1)
 816           eqj(1,nj)=.5*wpq
 817           eqj(2,nj)=eqj(1,nj)
 818           eqj(3,nj)=0.
 819           eqj(4,nj)=0.
 820         endif
 821
 822         ept(1)=.5*(wpi+wmi)
 823         ept(2)=.5*(wpi-wmi)
 824         ept(3)=0.
 825         ept(4)=0.
 826         jgamma=0
 827         ntry=0
 828       endif
 829
 830 8     continue
 831
 832 c ladder rung
 833 c---------------------------------------
 834       pt2=ept(3)**2+ept(4)**2
 835       pt=sqrt(pt2)
 836
 837       wpt(1)=ept(1)+ept(2)              !lc+ for the current jet emissi
 838       wpt(2)=ept(1)-ept(2)              !lc- for the current jet emissi
 839
 840       s2min=max(qmin(1),16.*qmin(2))    !mass cutoff for born
 841       s2min=max(s2min,4.*qq)
 842
 843       if(jj.eq.1)then
 844         wwmin=2.*s2min-2.*pt*sqrt(q2ini)
 845         wwmin=(wwmin+sqrt(wwmin**2+4.*pt2*(s2min-q2ini)))
 846      *  /(1.-q2ini/s2min)/2.
 847         sj=psjti(qmin(1),qq,si,iqc(1),iqc(2),1)   !total jet
 848         sj2=psjti1(q2min,qmin(1),qq,si,iqc(2),iqc(1),1)
 849         if(ish.ge.3)write(ifch,*)'resolved - si,wwmin,s2min,sj,sj2:'
 850         if(ish.ge.3)write(ifch,*)si,wwmin,s2min,sj,sj2
 851         if(sj.eq.0.)stop'sj=0'
 852         if(rangen().gt.sj2/sj.and.si.gt.1.1*wwmin)goto 26
 853         jj=2
 854       endif
 855       sj=psjti1(qmin(2),qmin(1),qq,si,iqc(2),iqc(1),1)
 856       sjb=psbint(qmin(1),qmin(2),qq,si,iqc(1),iqc(2),1) !born parton-parton
 857       wwmin=17./16*s2min-2.*pt*sqrt(q2ini)
 858       wwmin=(wwmin+sqrt(wwmin**2+pt2*(s2min/4.-4.*q2ini)))
 859      */(1.-16.*q2ini/s2min)/2.
 860       if(rangen().lt.sjb/sj.or.si.lt.1.1*wwmin)goto 10
 861
 862 26    continue
 863       wpt(jj)=wpt(jj)-pt2/wpt(3-jj)
 864
 865       if(jj.eq.1)then
 866         discr=(si+2.*pt*sqrt(q2ini))**2-4.*q2ini*(2.*si+pt2)
 867         if(discr.lt.0..and.ish.ge.1)write(ifmt,*)'discr,si,pt,wwmin',
 868      *  discr,si,pt,wwmin
 869         discr=sqrt(discr)
 870         qqmax=(si+2.*pt*sqrt(q2ini)+discr)/2./(2.+pt2/si)
 871       else
 872         discr=(si+2.*pt*sqrt(q2ini))**2-4.*q2ini*(17.*si+pt2)
 873         if(discr.lt.0..and.ish.ge.1)write(ifmt,*)'discr,si,pt,wwmin',
 874      *  discr,si,pt,wwmin
 875         discr=sqrt(discr)
 876         qqmax=(si+2.*pt*sqrt(q2ini)+discr)/2./(17.+pt2/si)
 877       endif
 878       qqmin=2.*q2ini*si/(si+2.*pt*sqrt(q2ini)+discr)
 879       if(jj.eq.1.and.s2min.gt.qqmin.or.
 880      *jj.eq.2.and.s2min.gt.16.*qqmin)then
 881         xmm=.5*(si-s2min+2.*pt*sqrt(q2ini))
 882         discr=xmm**2-q2ini*(si+pt2)
 883         if(discr.lt.0..and.ish.ge.1)write(ifmt,*)'discr1,si,pt,wwmin',
 884      *  discr,si,pt,wwmin
 885         qqmin=q2ini*si/(xmm+sqrt(discr))
 886       endif
 887
 888       xmin=1.-q2ini/qqmin
 889       xmax=1.-q2ini/qqmax
 890       if(ish.ge.6)write(ifch,*)'qqmin,qqmax,xmin,xmax',
 891      *qqmin,qqmax,xmin,xmax
 892       if(qqmin.lt.qmin(jj))then
 893         qqmin=qmin(jj)
 894         xmi=max(1.-((pt*sqrt(qqmin)+sqrt(pt2*qqmin+
 895      *  si*(si-s2min-qqmin*(1.+pt2/si))))/si)**2,
 896      *  (s2min+qqmin*(1.+pt2/si)-2.*pt*sqrt(qqmin))/si)
 897         xmin=max(xmin,xmi)
 898         if(xmin.le.0.)xmin=(s2min+qqmin*(1.+pt2/si))/si
 899         if(ish.ge.6)write(ifch,*)'qqmin,qmin(jj),xmin,s2min',
 900      *  qqmin,qmin(jj),xmin,s2min
 901       endif
 902
 903       qm0=qmin(jj)
 904       xm0=1.-q2ini/qm0
 905       if(xm0.gt.xmax.or.xm0.lt.xmin)then
 906         xm0=.5*(xmax+xmin)
 907       endif
 908 c      s2max=xm0*si
 909       s2max=xm0*si-qm0*(1.+pt2/si)+2.*pt*sqrt(q2ini)  !new ladder mass squared
 910       xx=xm0
 911
 912       if(jj.eq.1)then
 913         sj0=psjti(qm0,qq,s2max,0,iqc(2),1)*psfap(xx,iqc(1),0)+
 914      *  psjti(qm0,qq,s2max,7,iqc(2),1)*psfap(xx,iqc(1),1)
 915         gb0=sj0/log(q2ini/qcdlam)*sngl(psuds(qm0,iqc(1)))*qm0*2.
 916       else
 917         sj0=psjti1(qm0,qmin(1),qq,s2max,0,iqc(1),1)*psfap(xx,iqc(2),0)
 918      *  +psjti1(qm0,qmin(1),qq,s2max,7,iqc(1),1)*psfap(xx,iqc(2),1)
 919         gb0=sj0/log(q2ini/qcdlam)*sngl(psuds(qm0,iqc(2)))*qm0*2.
 920       endif
 921       if(gb0.le.0.)then
 922          write(ifmt,*)'gb0.le.0.  si,qq,pt2:',si,qq,pt2
 923          iret=1
 924          goto 9999
 925       endif
 926       if(xm0.le..5)then
 927         gb0=gb0*xm0**(1.-delh)
 928       else
 929         gb0=gb0*(1.-xm0)*2.**delh
 930       endif
 931
 932       xmin2=max(.5,xmin)
 933       xmin1=xmin**delh                 !xmin, xmax are put into powe
 934       xmax1=min(xmax,.5)**delh       !to simulate x value below
 935       if(xmin.ge..5)then
 936         djl=1.
 937       elseif(xmax.lt..5)then
 938         djl=0.
 939       else
 940         djl=1./(1.+((2.*xmin)**delh-1.)/delh/
 941      *  log(2.*(1.-xmax)))
 942       endif
 943
 944       ntry=0
 945 9     continue
 946       ntry=ntry+1
 947       if(ntry.ge.10000)then
 948         print *,"ntry.ge.10000"
 949         iret=1
 950         goto 9999
 951       endif
 952       if(jgamma.ne.1)then
 953         if(rangen().gt.djl)then        !lc momentum share in the cur
 954           x=(xmin1+rangen()*(xmax1-xmin1))**(1./delh)
 955         else
 956           x=1.-(1.-xmin2)*((1.-xmax)/(1.-xmin2))**rangen()
 957         endif
 958         q2=qqmin/(1.+rangen()*(qqmin/qqmax-1.))
 959         qt2=q2*(1.-x)
 960         if(ish.ge.6)write(ifch,*)'jj,q2,x,qt2',jj,q2,x,qt2
 961         if(qt2.lt.q2ini)goto 9
 962       else
 963         x=xmin+rangen()*(xmax-xmin)
 964         q2=qqmin*(qqmax/qqmin)**rangen()
 965         qt2=(q2-x*qq)*(1.-x)
 966         if(ish.ge.6)write(ifch,*)'jj,q2,x,qt2',jj,q2,x,qt2
 967         if(qt2.lt.0.)goto 9
 968       endif
 969
 970       qt=sqrt(qt2)
 971       call pscs(bcos,bsin)
 972 c ep3 is now 4-vector for s-channel gluon produced in current ladder run
 973       ep3(3)=qt*bcos
 974       ep3(4)=qt*bsin
 975       ptnew=(ept(3)-ep3(3))**2+(ept(4)-ep3(4))**2
 976       if(jj.eq.1)then
 977         s2min2=max(q2,s2min)
 978       else
 979         s2min2=max(s2min,16.*q2)
 980       endif
 981
 982       if(jgamma.ne.1)then
 983         s2=x*si-q2*(1.+pt2/si)-ptnew+pt2+qt2  !new ladder mass squared
 984         if(s2.lt.s2min2)goto 9      !rejection in case of too low mass
 985         xx=x
 986
 987         if(jj.eq.1)then
 988           sj1=psjti(q2,qq,s2,0,iqc(2),1)
 989           if(iqc(1).ne.0)then
 990             sj2=psjti(q2,qq,s2,iqc(1),iqc(2),1)
 991           elseif(iqc(2).eq.0)then
 992             sj2=psjti(q2,qq,s2,1,0,1)
 993           else
 994             sj2=psjti(q2,qq,s2,1,1,1)/6.+
 995      *      psjti(q2,qq,s2,-1,1,1)/6.+
 996      *      psjti(q2,qq,s2,2,1,1)/1.5
 997           endif
 998         else
 999           sj1=psjti1(q2,qmin(1),qq,s2,0,iqc(1),1)
1000           if(iqc(2).ne.0)then
1001             sj2=psjti1(q2,qmin(1),qq,s2,iqc(2),iqc(1),1)
1002           elseif(iqc(1).eq.0)then
1003             sj2=psjti1(q2,qmin(1),qq,s2,1,0,1)
1004           else
1005             sj2=psjti1(q2,qmin(1),qq,s2,1,1,1)/6.+
1006      *      psjti1(q2,qmin(1),qq,s2,-1,1,1)/6.+
1007      *      psjti1(q2,qmin(1),qq,s2,2,1,1)/1.5
1008           endif
1009         endif
1010 c gb7 is the rejection function for x and q**2 simulation
1011         gb7=(sj1*psfap(xx,iqc(jj),0)+sj2*psfap(xx,iqc(jj),1))
1012      *  /log(qt2/qcdlam)*sngl(psuds(q2,iqc(jj)))*q2/gb0
1013
1014         if(x.le..5)then
1015           gb7=gb7*x**(1.-delh)
1016         else
1017           gb7=gb7*(1.-x)*2.**delh
1018         endif
1019       else
1020         s2=x*si-q2               !new ladder mass squared
1021         if(s2.lt.s2min2)goto 9   !rejection in case of too low mass
1022
1023         sj1=0.
1024         xx=x
1025         if(iqc(2).eq.0)then
1026           sj2=psjti(q2,qq,s2,1,0,1)
1027         else
1028           sj2=psjti(q2,qq,s2,1,1,1)/naflav/2.+
1029      *    psjti(q2,qq,s2,-1,1,1)/naflav/2.+
1030      *    psjti(q2,qq,s2,2,1,1)*(1.-1./naflav)
1031         endif
1032         gb7=sj2*psfap(xx,0,1)/gb0  !????*(1.-x*qq/q2)
1033       endif
1034       if(gb7.gt.1..or.gb7.lt.0..and.ish.ge.1)write(ifmt,*)'gb7,q2,x,gb0'
1035      *,gb7,q2,x,gb0
1036       if(rangen().gt.gb7)goto 9
1037
1038       if(ish.ge.6)write(ifch,*)'res: jj,iqc,ncc:',
1039      *jj,iqc(jj),ncc(1,jj),ncc(2,jj)
1040
1041       nqc(2)=0
1042       iqnew=iqc(jj)
1043       if(jgamma.ne.1)then
1044         if(rangen().lt.sj1/(sj1+sj2))then
1045           if(iqc(jj).eq.0)then
1046             jt=1
1047             jq=int(1.5+rangen())
1048             nqc(1)=ncc(jq,jj)
1049           else
1050             jt=2
1051             if(iqc(jj).gt.0)then
1052               jq=1
1053             else
1054               jq=2
1055             endif
1056             nqc(1)=0
1057             iqnew=0
1058           endif
1059           iq1=iqc(jj)
1060         else
1061           if(iqc(jj).ne.0)then
1062             iq1=0
1063             jt=3
1064             if(iqc(jj).gt.0)then
1065               jq=1
1066             else
1067               jq=2
1068             endif
1069             nqc(1)=ncc(1,jj)
1070           else
1071             jt=4
1072             jq=int(1.5+rangen())
1073             iq1=int(naflav*rangen()+1.)*(3-2*jq)
1074             nqc(1)=ncc(jq,jj)
1075             iqnew=-iq1
1076           endif
1077         endif
1078       else
1079         jt=5
1080         jq=int(1.5+rangen())
1081         iq1=int(naflav*rangen()+1.)*(3-2*jq)
1082         iqnew=-iq1
1083         nqc(1)=0
1084       endif
1085       eprt=max(1.d0*qt,
1086      *.5d0*((1.d0-x)*wpt(jj)+qt2/(1.d0-x)/wpt(jj)))
1087       pl=((1.d0-x)*wpt(jj)-eprt)*(3-2*jj)
1088       if(iq1.eq.0)then
1089         iq2ini=9
1090       else
1091         iq2ini=iq1
1092       endif
1093 27    call timsh1(q2,sngl(eprt),iq2ini)
1094       amprt=pprt(5,1)**2
1095       plprt=eprt**2-amprt-qt2
1096       if(plprt.lt.-1d-6)goto 27
1097       ep3(1)=eprt
1098       ep3(2)=dsqrt(max(0.d0,plprt))
1099       if(pl.lt.0.d0)ep3(2)=-ep3(2)
1100       ey(1)=1.
1101       ey(2)=1.
1102       ey(3)=1.
1103       do i=1,4
1104         ept1(i)=ept(i)-ep3(i)
1105       enddo
1106       call psdefrot(ep3,s0xh,c0xh,s0h,c0h)
1107       if(ish.ge.6)then
1108       write(ifch,*)'q2,amprt,qt2',q2,amprt,qt2
1109       write(ifch,*)'eprt,plprt',eprt,plprt
1110       write(ifch,*)'ep3',ep3
1111       write(ifch,*)'ept',ept
1112       write(ifch,*)'ept1',ept1
1113       endif
1114       s2new=psnorm(ept1)
1115
1116       if(s2new.gt.s2min2)then
1117         if(jj.eq.1)then
1118           gb=psjti(q2,qq,s2new,iqnew,iqc(2),1)
1119         else
1120           gb=psjti1(q2,qmin(1),qq,s2new,iqnew,iqc(1),1)
1121         endif
1122         if(iqnew.eq.0)then
1123           gb=gb/sj1
1124         else
1125           gb=gb/sj2
1126         endif
1127         if(ish.ge.1)then
1128           if(gb.gt.1.)write (ifch,*)'gb,s2new,s2,q2,iqnew',
1129      *    gb,s2new,s2,q2,iqnew
1130         endif
1131         if(rangen().gt.gb)goto 9
1132       else
1133         goto 9
1134       endif
1135       jgamma=0
1136
1137       call psreti(nqc,jq,1,ey,s0xh,c0xh,s0h,c0h)
1138
1139       if(jt.eq.1)then
1140         ncc(jq,jj)=nqc(2)
1141       elseif(jt.eq.2)then
1142         ncc(jq,jj)=ncc(1,jj)
1143         ncc(3-jq,jj)=nqc(1)
1144       elseif(jt.eq.3)then
1145         ncc(1,jj)=nqc(2)
1146       elseif(jt.eq.4)then
1147         ncc(1,jj)=ncc(3-jq,jj)
1148       elseif(jt.eq.5)then
1149         ncc(1,jj)=nqc(1)
1150         ncc(2,jj)=0
1151       endif
1152       iqc(jj)=iqnew
1153       if(ish.ge.6)write(ifch,*)'qt2,amprt,ncc:',
1154      *qt2,amprt,ncc(1,jj),ncc(2,jj)
1155
1156       do i=1,4
1157         ept(i)=ept1(i)
1158       enddo
1159 c c.m. energy squared, minimal  4-momentum transfer square and gluon 4-v
1160 c for the next ladder run
1161       qmin(jj)=q2
1162       si=s2new
1163       if(ish.ge.3)write (ifch,*)'res: new jet - iqj,ncj,ep3,ept',
1164      *iqj(nj),ncj(1,nj),ncj(2,nj),ep3,ept
1165
1166       goto 8            !next simulation step will be considered
1167
1168 10    continue
1169       if(ish.ge.3)write(ifch,*)'res: iqc,si,ept:',iqc,si,ept
1170
1171 c highest virtuality subprocess in the ladder
1172 c---------------------------------------
1173       qqs=max(qmin(1)/4.,4.*qmin(2))
1174       qqs=max(qqs,qq)
1175       call psabor(si,qqs,iqc,ncc,ept,1,nptlh,bx)
1176       goto 17
1177
1178 12    continue
1179 c---------------------------------------
1180 c hard pomeron (direct photon)
1181 c---------------------------------------
1182       ept(1)=.5*(wpi+wmi)
1183       ept(2)=.5*(wpi-wmi)
1184       ept(3)=0.
1185       ept(4)=0.
1186       if(ish.ge.3)write (ifch,*)'direct photon - ept,si,qq:',ept,si,qq
1187
1188 13    continue
1189
1190 c ladder rung
1191 c---------------------------------------
1192       pt2=ept(3)**2+ept(4)**2
1193       pt=sqrt(pt2)
1194       wpt(1)=ept(1)+ept(2)
1195       wpt(2)=si/wpt(1)
1196
1197       gdbor=psdbin(qmin(2),qq,si,iqc(2),long)
1198       gdtot=psdsin(qmin(2),qq,si,iqc(2),long)
1199       if(iqc(2).ne.0)then
1200         if(ish.ge.8)write (ifch,*)'qmin(2),qq,si',qmin(2),qq,si
1201         gdnon=psdnsi(qmin(2),qq,si,long)
1202         if(iabs(iqc(2)).eq.1.or.iabs(iqc(2)).eq.4)then
1203           gdbor=gdbor/2.25
1204           gdtot=gdnon/2.25+gdtot/4.5
1205         else
1206           gdbor=gdbor/9.
1207           gdtot=gdnon/9.+gdtot/4.5
1208         endif
1209       else
1210         gdnon=0.
1211       endif
1212
1213       if(long.ne.0.or.qmin(2).ge.qq)then
1214         s2min=qq+4.*max(qmin(2),qcmass**2)
1215         wwmin=(5.*s2min-qq)/4.-2.*pt*sqrt(q2ini)
1216         wwmin=(wwmin+sqrt(wwmin**2-(qq-pt2)*(s2min-qq-4.*q2ini)))
1217      *  /2./(1.-4.*q2ini/(s2min-qq))
1218       else
1219         s2min=qq/(1.-sqrt(q2ini/qq))
1220         wwmin=s2min+qq-2.*pt*sqrt(q2ini)
1221         wwmin=(wwmin+sqrt(wwmin**2-4.*(qq-pt2)*(qq-q2ini)))
1222      *  /2./(1.-q2ini/qq)
1223       endif
1224
1225       if(ish.ge.3)write(ifch,*)'si,s2min,wwmin,qmin(2),gdtot,gdbor:'
1226       if(ish.ge.3)write(ifch,*)si,s2min,wwmin,qmin(2),gdtot,gdbor
1227
1228       if((rangen().lt.gdbor/gdtot.or.si.lt.1.1*wwmin).and.
1229      *(long.eq.0.and.qmin(2).lt.qq.or.iqc(2).eq.0))goto 15
1230       if(si.lt.1.1*wwmin)stop'si<1.1*wwmin'
1231
1232       qqmax=0.
1233       qqmin=0.
1234
1235       xmm=si+2.*sqrt(q2ini)*pt-qq
1236       discr=xmm**2-4.*q2ini*(5.*si-qq+pt2)
1237       if(discr.lt.0.)goto 29
1238       discr=sqrt(discr)
1239       qqmax=(xmm+discr)/2./(5.-(qq-pt2)/si)
1240       qqmin=2.*q2ini*si/(xmm+discr)
1241
1242 29    continue
1243       if(4.*qqmin.lt.s2min-qq.or.long.eq.0.and.
1244      *qmin(2).lt.qq)then
1245         xmm=si-s2min+2.*sqrt(q2ini)*pt
1246         qqmin=2.*q2ini*si/(xmm+sqrt(xmm**2-4.*q2ini*(si-qq+pt2)))
1247       endif
1248       xmin=1.-q2ini/qqmin
1249
1250       if(qqmin.lt.qmin(2))then
1251         qqmin=qmin(2)
1252         xmi=max(1.-((pt*sqrt(qqmin)+sqrt(pt2*qqmin+
1253      *  si*(si-s2min-qqmin*(1.-(qq-pt2)/si))))/si)**2,
1254      *  (s2min+qqmin*(1.-(qq-pt2)/si)-2.*pt*sqrt(qqmin))/si)
1255         xmin=max(xmin,xmi)
1256       endif
1257       if(xmin.le.qq/si)xmin=1.001*qq/si
1258
1259       if(long.eq.0.and.qmin(2).lt.qq)qqmax=max(qqmax,qq)
1260       xmax=1.-q2ini/qqmax
1261
1262       if(ish.ge.6)write(ifch,*)'qqmax,qqmin,xmax,xmin:',
1263      *qqmax,qqmin,xmax,xmin
1264       if(qqmax.lt.qqmin)stop'qqmax<qqmin'
1265
1266       qm0=qqmin
1267       xm0=1.-q2ini/qm0
1268       s2max=si*xm0-qm0*(1.-qq/si)
1269
1270       sds=psdsin(qm0,qq,s2max,0,long)/4.5
1271       sdv=psdsin(qm0,qq,s2max,1,long)/4.5
1272
1273       sdn=psdnsi(qm0,qq,s2max,long)
1274       if(iqc(2).eq.0)then
1275         sdn=sdn/4.5
1276       elseif(iabs(iqc(2)).eq.1.or.iabs(iqc(2)).eq.4)then
1277         sdn=sdn/2.25
1278       else
1279         sdn=sdn/9.
1280       endif
1281       sdv=sdv+sdn
1282       xx=xm0
1283
1284       sj0=sds*psfap(xx,iqc(2),0)+sdv*psfap(xx,iqc(2),1)
1285       gb0=sj0/log(q2ini/qcdlam)*sngl(psuds(qm0,iqc(2)))*qm0*5.
1286       if(gb0.le.0.)then
1287          write(ifmt,*)'gb0.le.0.  si,qq,pt2:',si,qq,pt2
1288          iret=1
1289          goto 9999
1290       endif
1291
1292       if(xm0.le..5)then
1293         gb0=gb0*(xm0-qq/si)/(1.-2.*qq/si)
1294       else
1295         gb0=gb0*(1.-xm0)
1296       endif
1297
1298       xmin2=max(.5,xmin)
1299       xmax1=min(xmax,.5)
1300       if(xmin.ge..5)then
1301         djl=1.
1302       elseif(xmax.lt..5)then
1303         djl=0.
1304       else
1305         djl=1./(1.-(1.-2.*qq/si)*log((.5-qq/si)/(xmin-qq/si))/
1306      *  log(2.*(1.-xmax)))
1307       endif
1308
1309 14    continue
1310       if(rangen().gt.djl)then        !lc momentum share in the cur
1311         x=(xmin-qq/si)*((xmax1-qq/si)/(xmin-qq/si))**rangen()+qq/si
1312       else
1313         x=1.-(1.-xmin2)*((1.-xmax)/(1.-xmin2))**rangen()
1314       endif
1315
1316       q2=qqmin/(1.+rangen()*(qqmin/qqmax-1.))
1317
1318       qt2=q2*(1.-x)
1319       if(ish.ge.9)write(ifch,*)'q2,x,qt2,qq,qqmin,qqmax:',
1320      *q2,x,qt2,qq,qqmin,qqmax
1321       if(qt2.lt.q2ini)goto 14   !p_t check
1322
1323       if(long.ne.0.or.q2.ge.qq)then
1324         s2min2=max(4.*q2+qq,s2min)
1325       else
1326         s2min2=s2min
1327       endif
1328       qt=sqrt(qt2)
1329       call pscs(bcos,bsin)
1330 c ep3 is now 4-vector for s-channel gluon produced in current ladder run
1331       ep3(3)=qt*bcos
1332       ep3(4)=qt*bsin
1333       ptnew=(ept(3)-ep3(3))**2+(ept(4)-ep3(4))**2
1334
1335       s2=x*si-ptnew+pt2-q2*(x-(qq-pt2)/si)
1336       if(s2.lt.s2min2)goto 14   !check of the kinematics
1337       sds=psdsin(q2,qq,s2,0,long)/4.5
1338       sdv0=psdsin(q2,qq,s2,1,long)/4.5
1339       if(ish.ge.8)write (ifch,*)'q2,qq,s2',q2,qq,s2
1340       sdn0=psdnsi(q2,qq,s2,long)
1341
1342       if(iqc(2).eq.0)then
1343         sdn=sdn0/4.5
1344       else
1345         if(iabs(iqc(2)).eq.1.or.iabs(iqc(2)).eq.4)then
1346           sdn=sdn0/2.25
1347         else
1348           sdn=sdn0/9.
1349         endif
1350       endif
1351       sdv=sdv0+sdn
1352
1353       xx=x
1354       sj1=sds*psfap(xx,iqc(2),0)
1355       sj2=sdv*psfap(xx,iqc(2),1)
1356
1357 c gb7 is the rejection function for x and q**2 simulation.
1358       gb7=(sj1+sj2)/log(qt2/qcdlam)*sngl(psuds(q2,iqc(2)))/gb0*q2
1359       if(x.le..5)then
1360         gb7=gb7*(x-qq/si)/(1.-2.*qq/si)
1361       else
1362         gb7=gb7*(1.-x)
1363       endif
1364
1365       if(gb7.gt.1..and.ish.ge.1)write(ifmt,*)'gb7,q2,x,qt2,iqc(2),'
1366      * ,'gb0,sj1,sj2',gb7,q2,x,qt2,iqc(2),gb0,sj1,sj2
1367       if(ish.ge.3)write (ifch,*)'gb7,q2,x,qt2,iqc(2),gb0,sj1,sj2,long',
1368      * gb7,q2,x,qt2,iqc(2),gb0,sj1,sj2,long
1369       if(rangen().gt.gb7)goto 14
1370
1371
1372       if(ish.ge.6)write(ifch,*)'iqc,ncc:',iqc(2),ncc(1,2),ncc(2,2)
1373       iqcnew=iqc(2)
1374       nqc(2)=0         !emitted parton color connections
1375       if(rangen().lt.sj1/(sj1+sj2).or.(long.ne.0.or.q2.ge.qq).and.
1376      *s2.lt.1.5*s2min2)then
1377         if(iqc(2).eq.0)then
1378           jt=1
1379           jq=int(1.5+rangen())
1380           nqc(1)=ncc(jq,2)
1381         else
1382           jt=2
1383           if(iqc(2).gt.0)then
1384             jq=1
1385           else
1386             jq=2
1387           endif
1388           nqc(1)=0
1389         endif
1390         iq1=iqc(2)
1391         iqcnew=0
1392
1393       else
1394         if(iqc(2).ne.0)then
1395           jt=3
1396           iq1=0
1397           if(iqc(2).gt.0)then
1398             jq=1
1399           else
1400             jq=2
1401           endif
1402           nqc(1)=ncc(1,2)
1403
1404         else
1405           tu=sdn0/2.25+sdv0
1406           if(naflav.eq.4)tu=tu*2.
1407           td=sdn0/9.+sdv0
1408           if(rangen().lt.tu/(tu+2.*td))then
1409             if(naflav.eq.3)then
1410               iq1=1
1411             else
1412               iq1=1+3*int(.5+rangen())
1413             endif
1414           else
1415             iq1=int(2.5+rangen())
1416           endif
1417           jq=int(1.5+rangen())
1418           iq1=iq1*(3-2*jq)
1419           iqcnew=-iq1
1420           jt=4
1421           nqc(1)=ncc(jq,2)
1422         endif
1423       endif
1424
1425       eprt=max(1.d0*qt,
1426      *.5d0*((1.d0-x)*wpt(2)+qt2/(1.d0-x)/wpt(2)))
1427       pl=eprt-(1.d0-x)*wpt(2)
1428       if(iq1.eq.0)then
1429         iq2ini=9
1430       else
1431         iq2ini=iq1
1432       endif
1433 28    call timsh1(q2,sngl(eprt),iq2ini)
1434       amprt=pprt(5,1)**2
1435       plprt=eprt**2-amprt-qt2
1436       if(plprt.lt.-1d-6)goto 28
1437       ep3(1)=eprt
1438       ep3(2)=dsqrt(max(0.d0,plprt))
1439       if(pl.lt.0.d0)ep3(2)=-ep3(2)
1440       ey(1)=1.
1441       ey(2)=1.
1442       ey(3)=1.
1443       do i=1,4
1444         ept1(i)=ept(i)-ep3(i)
1445       enddo
1446       call psdefrot(ep3,s0xh,c0xh,s0h,c0h)
1447       call psrotat(ep3,s0xh,c0xh,s0h,c0h)
1448       s2new=psnorm(ept1)+qq
1449
1450       if((long.ne.0.or.q2.ge.qq).and.iqcnew.ne.0)then
1451         xmm=(5.*s2min2-qq)/4.-2.*sqrt(ptnew*q2ini)
1452         s2min2=1.1*(xmm+sqrt(xmm**2-(qq-ptnew)*
1453      *  (s2min2-qq-4.*q2ini)))/2./(1.-4.*q2ini/(s2min2-qq))
1454       endif
1455       if(s2new.gt.s2min2)then
1456         sds1=psdsin(q2,qq,s2new,iqcnew,long)/4.5
1457         if(iqcnew.eq.0)then
1458           gb=sds1/sds
1459         else
1460           if(ish.ge.8)write (ifch,*)'q2,qq,s2new',q2,qq,s2new
1461           sdn1=psdnsi(q2,qq,s2new,long)
1462           if(iabs(iqcnew).eq.1.or.iabs(iqcnew).eq.4)then
1463             sdn1=sdn1/2.25
1464             sdv=sdv0+sdn0/2.25
1465           else
1466             sdn1=sdn1/9.
1467             sdv=sdv0+sdn0/9.
1468           endif
1469           gb=.9999*(sds1+sdn1)/sdv
1470         endif
1471         if(ish.ge.3.and.gb.gt.1..and.ish.ge.1)write(ifmt,*)'gbs2',gb
1472         if(rangen().gt.gb)goto 14
1473       else
1474         goto 14
1475       endif
1476
1477       call psreti(nqc,jq,1,ey,s0xh,c0xh,s0h,c0h)
1478
1479       iqc(2)=iqcnew
1480       if(jt.eq.1)then      !current parton color connections
1481         ncc(jq,2)=nqc(2)
1482       elseif(jt.eq.2)then
1483         ncc(jq,2)=ncc(1,2)
1484         ncc(3-jq,2)=nqc(1)
1485       elseif(jt.eq.3)then
1486         ncc(1,2)=nqc(2)
1487       elseif(jt.eq.4)then
1488         ncc(1,2)=ncc(3-jq,2)
1489         ncc(2,2)=0
1490       endif
1491
1492       do i=1,4
1493         ept(i)=ept1(i)
1494       enddo
1495       if(ish.ge.3)write (ifch,*)'new jet - iqj,ncj,ep3,ept',
1496      *iqj(nj),ncj(1,nj),ncj(2,nj),ep3,ept
1497 c c.m. energy squared, minimal  4-momentum transfer square and gluon 4-v
1498 c for the next ladder run
1499       qmin(2)=q2
1500       si=s2new
1501       goto 13            !next simulation step will be considered
1502
1503 15    continue
1504       if(ish.ge.3)write(ifch,*)'iqc,si,qmin(2),nj:',
1505      *iqc(2),si,qmin(2),nj
1506 c highest virtuality subprocess in the ladder
1507 c---------------------------------------
1508       gb01=0.
1509       tmax=0.
1510       tmin=si
1511       if(iqc(2).eq.0.and.si.gt.qq+4.*max(qcmass**2,qmin(2)))then
1512         qminn=max(qcmass**2,qmin(2))
1513         tmin1=2.*qminn/(1.-qq/si)/(1.+sqrt(1.-4.*qminn/(si-qq)))
1514         tmin=tmin1
1515         tmax=si/2.
1516         fb01=psdbom(si,si/2.,si/2.,qq,long)
1517         if(long.eq.0)fb01=fb01*si/2.
1518         gb01=fb01/log(qminn/qcdlam)*sngl(psuds(qminn,iqc(2)))/si**2
1519         gb0=gb01
1520       endif
1521
1522       if(long.eq.0.and.qmin(2).lt.qq)then
1523         tmax=max(tmax,qq)
1524         tmin=max(qmin(2),
1525      *  2.*q2ini/(1.-qq/si)/(1.+sqrt(1.-4.*q2ini/(si-qq))))
1526         ze=qq/si+tmin/si*(1.-qq/si)
1527         xx=ze
1528         qt2=tmin*(1.-ze)
1529         if(qt2.lt..999*q2ini.and.ish.ge.1)write(ifmt,*)'bor-dir:qt20'
1530      *                                                 ,qt2
1531         gb0=gb01+psfap(xx,iqc(2),1)/log(qt2/qcdlam)
1532      *  *sngl(psuds(tmin,iqc(2))/psuds(tmin,1)*psuds(qq,1))
1533      *  /si*(1.-tmin*qq/si**2/ze)
1534       endif
1535       gb0=gb0*2.
1536
1537       call psdeftr(si-qq,ept,ey)
1538
1539       if(ish.ge.6)write(ifch,*)'tmin,tmax,qq,si-qq,gb0:'
1540       if(ish.ge.6)write(ifch,*)tmin,tmax,qq,si-qq,psnorm(ept),gb0
1541
1542 c------------------------------------------------
1543 16    continue
1544       if(long.eq.0)then
1545         t=tmin*(tmax/tmin)**rangen()
1546       else
1547         t=tmin+(tmax-tmin)*rangen()
1548       endif
1549
1550       u=si-t
1551       ze=qq/si+t/si*(1.-qq/si)
1552       qt2=t*(1.-ze)
1553       if(t.le.qq.and.long.eq.0)then
1554         xx=ze
1555         gb=psfap(xx,iqc(2),1)/log(qt2/qcdlam)*sngl(psuds(t,iqc(2))
1556      *  /psuds(t,1)*psuds(qq,1))/si*(1.-t*qq/si**2/ze)/gb0
1557       else
1558         gb=0.
1559       endif
1560
1561       gb1=0.
1562       if(iqc(2).eq.0..and.si.gt.qq+4.*max(qcmass**2,qmin(2)).
1563      *and.qt2.gt.qcmass**2.and.t.le.si/2..and.t.ge.tmin1)then
1564         fb1=psdbom(si,t,u,qq,long)
1565         if(long.eq.0)fb1=fb1*t
1566         gb1=fb1/log(qt2/qcdlam)*sngl(psuds(qt2,iqc(2)))/si**2/gb0
1567 c        gb1=0.  !???????????????????????
1568         gb=gb+gb1
1569       endif
1570
1571
1572       if(ish.ge.6)write(ifch,*)'gb,t,iqc(2),si,qq,qmin(2),long:',
1573      *gb,t,iqc(2),si,qq,qmin(2),long
1574       if (ish.ge.1) then
1575         if(gb.gt.1.)write(*,*)'gb,gb1,gb0,gb01',
1576      *  ',t,iqc(2),si,qq,qmin(2),long:',
1577      *  gb,gb1,gb0,gb01,fb1,fb01,t,iqc(2),si,qq,qmin(2),long
1578       endif
1579
1580       if(rangen().gt.gb)goto 16
1581       if(ish.ge.3)write(ifch,*)'born:t,qt2:',t,qt2
1582
1583       nqc(2)=0
1584       if(iqc(2).eq.0)then
1585         jq=int(1.5+rangen())
1586         jq2=3-jq
1587         if(rangen().gt.gb1/gb)then
1588           iq1=(1+int(3.*rangen()))*(3-2*jq)
1589         else
1590           iq1=4*(3-2*jq)
1591         endif
1592         iq2=-iq1                             !quark flavors
1593         nqc(1)=ncc(jq,2)
1594       else
1595         if(iqc(2).gt.0)then
1596           jq=1
1597         else
1598           jq=2
1599         endif
1600         jq2=jq
1601         iq1=0
1602         iq2=iqc(2)
1603         nqc(1)=ncc(1,2)
1604       endif
1605
1606       call pscs(bcos,bsin)
1607       z=sngl(psutz(dble(si-qq),dble(qt2),dble(qt2)))
1608       if(t.lt..5*si)z=1.-z
1609       wp3=z*sqrt(si-qq)
1610       wm3=qt2/wp3
1611       if(iabs(iq1).eq.4)qt2=qt2-qcmass**2
1612       qt=sqrt(qt2)
1613       ep3(1)=.5*(wp3+wm3)
1614       ep3(2)=.5*(wp3-wm3)
1615       ep3(3)=qt*bcos
1616       ep3(4)=qt*bsin
1617       call psdefrot(ep3,s0xh,c0xh,s0h,c0h)
1618       if(iq1.eq.0)then
1619         iq2ini1=9
1620       else
1621         iq2ini1=iq1
1622       endif
1623       if(iq2.eq.0)then
1624         iq2ini2=9
1625       else
1626         iq2ini2=iq2
1627       endif
1628       if(ish.ge.5)write (ifch,*)'jq,jt,iq2ini1,iq2ini2',
1629      *jq,jt,iq2ini1,iq2ini2
1630
1631       if(t.lt.qq.and.iabs(iq1).ne.4)then
1632         qq1=t*(1.-ze)
1633         qq2=qq
1634       else
1635         qq1=qt2
1636         qq2=qt2
1637       endif
1638       call timsh2(qq1,qq2,sqrt(si-qq),iq2ini1,iq2ini2)
1639       nfprt=1
1640       call psreti(nqc,jq,nfprt,ey,s0xh,c0xh,s0h,c0h)
1641
1642       if(iqc(2).eq.0)then
1643         nqc(1)=ncc(3-jq,2)
1644         nqc(2)=0
1645       else
1646         nqc(1)=nqc(2)
1647         nqc(2)=0
1648       endif
1649
1650       nfprt=2
1651       call psreti(nqc,jq2,nfprt,ey,s0xh,c0xh,s0h,c0h)
1652
1653 17    continue
1654       if(ish.ge.3)write (ifch,*)'nj',nj
1655       if(nj.gt.0)then
1656
1657           ityj(i)=30
1658           iorj(i)=nptlh
1659         do n=1,nj
1660           do i=1,4
1661             ep3(i)=eqj(i,n)
1662           enddo
1663           call pstrans(ep3,ey0,-1)         !boost to the c.m. system
1664           do i=1,4
1665             eqj(i,n)=ep3(i)
1666           enddo
1667           do l=1,6
1668             bxj(l,n)=bx(l)
1669           enddo
1670           ityj(n)=0
1671           iorj(n)=1
1672         enddo
1673       endif
1674       call psjarr(jfl)       !kinky strings formation
1675       if(ish.ge.3)write (ifch,*)'jfl',jfl
1676       if(jfl.eq.0)then
1677         iret=1
1678       else
1679         iret=0
1680         ep3(4)=egyevt
1681         ep3(2)=0.
1682         ep3(3)=0.
1683         ep3(1)=0.
1684         do i=2,nptl
1685           do l=1,4
1686             ep3(l)=ep3(l)-pptl(l,i)
1687           enddo
1688         enddo
1689         if(ish.ge.3)write(ifch,*)'energy-momentum balance:'
1690         if(ish.ge.3)write(ifch,*)ep3
1691         if(abs(ep3(4)).gt.3.e-2)write(*,*)'energy-momentum balance:',ep3
1692       endif
1693  9999 call utprix('psadis',ish,ishini,3)
1694       return
1695       end
1696
1697 c-----------------------------------------------------------------------
1698       subroutine psaevc
1699 c-----------------------------------------------------------------------
1700       include 'epos.inc'
1701 c structure functions calculation
1702       logical lcalc
1703       double precision xx,xmax
1704       dimension evs(21,21,135)
1705       common /psar2/  edmax,epmax
1706       common /psar31/ evk0(21,21,54)
1707       common /psar32/ evk(21,21,135)
1708       include 'epos.incsem'
1709
1710       inquire(file=fnie(1:nfnie),exist=lcalc)
1711       if(lcalc)then
1712        if(inicnt.eq.1)then
1713         write(ifmt,'(3a)')'read from ',fnie(1:nfnie),' ...'
1714         open(1,file=fnie(1:nfnie),status='old')
1715         read (1,*)qcdlam0,q2min0,q2ini0,naflav0,epmax0
1716         if(qcdlam0.ne.qcdlam)write(ifmt,'(a)')'iniev: wrong qcdlam'
1717         if(q2min0 .ne.q2min )write(ifmt,'(a)')'iniev: wrong q2min'
1718         if(q2ini0 .ne.q2ini )write(ifmt,'(a)')'iniev: wrong q2ini'
1719         if(naflav0.ne.naflav)write(ifmt,'(a)')'iniev: wrong naflav'
1720         if(epmax0 .ne.epmax )write(ifmt,'(a)')'iniev: wrong epmax'
1721         if(qcdlam0.ne.qcdlam.or.q2min0.ne.q2min.or.q2ini0.ne.q2ini
1722      *  .or.naflav0.ne.naflav.or.epmax0.ne.epmax)then
1723            write(6,'(//a//)')'   iniev has to be reinitialized!!!'
1724            stop
1725         endif
1726         read (1,*)evk0,evk
1727         close(1)
1728        endif
1729        goto 101
1730       endif
1731
1732       write(ifmt,'(a)')'iniev does not exist -> calculate tables  ...'
1733       xmax=1.d0-2.d0*q2ini/epmax
1734       do l=1,27
1735         if(l.le.12)then
1736           xx=.1d0*exp(l-13.d0)
1737         elseif(l.le.21)then
1738           xx=.1d0*(l-12.d0)
1739         else
1740           xx=1.d0-.1d0*(10.d0*(1.d0-xmax))**((l-21)/6.)
1741         endif
1742
1743         qmin=max(1.d0*q2min,q2ini/(1.-xx))
1744       do i=1,21
1745         qq=qmin*(.5*epmax/qmin)**((i-1)/20.)
1746       do j=1,21
1747         qj=qmin*(qq/qmin)**((j-1)/20.)
1748         if(l.eq.27.or.i.eq.1.or.j.eq.21)then
1749           evk0(i,j,l)=0.
1750           evk0(i,j,l+27)=0.
1751           do k=1,5
1752             evk(i,j,l+27*(k-1))=0.
1753           enddo
1754         else
1755           do k=1,2
1756             evk0(i,j,l+27*(k-1))=log(psev0(qj,qq,xx,k))
1757           enddo
1758         endif
1759       enddo
1760       enddo
1761       enddo
1762
1763       n=1
1764
1765 1     n=n+1
1766       write(ifmt,2)n
1767 2     format(5x,i2,'-th order contribution')
1768
1769       do l=1,26
1770         write(ifmt,*)'l',l
1771         if(l.le.12)then
1772           xx=.1d0*exp(l-13.d0)
1773         elseif(l.le.21)then
1774           xx=.1d0*(l-12.d0)
1775         else
1776           xx=1.d0-.1d0*(10.d0*(1.d0-xmax))**((l-21)/6.)
1777         endif
1778
1779         qmin=max(1.d0*q2min,q2ini/(1.d0-xx))
1780       do i=2,21
1781         qq=qmin*(.5*epmax/qmin)**((i-1)/20.)
1782       do j=1,20
1783         qj=qmin*(qq/qmin)**((j-1)/20.)
1784         do m=1,3
1785         do k=1,2
1786           if(m.ne.3)then
1787             ev=psev(qj,qq,xx,m,k,n)
1788             ev0=psevi0(qj,qq,xx,m,k)
1789             evs(i,j,l+27*(m-1)+54*(k-1))=log((ev+ev0)/psfap(xx,m-1,k-1)
1790      *      /log(log(qq*(1.d0-xx)/qcdlam)/log(qj*(1.d0-xx)/qcdlam))*4.5)
1791           elseif(k.ne.1)then
1792             evs(i,j,l+108)=log((psev(qj,qq,xx,m,k,n)+
1793      *      psevi0(qj,qq,xx,2,2))/psfap(xx,2,2)
1794      *      /log(log(qq*(1.d0-xx)/qcdlam)/log(qj*(1.d0-xx)/qcdlam))*4.5)
1795           endif
1796         enddo
1797         enddo
1798       enddo
1799       enddo
1800       enddo
1801
1802       jec=0
1803       do i=2,21
1804       do j=1,20
1805       do l=1,26
1806       do k=1,5
1807         if(n.eq.2.or.evs(i,j,l+27*(k-1)).ne.0..and.
1808      *  abs(1.-evk(i,j,l+27*(k-1))/evs(i,j,l+27*(k-1))).gt.1.e-2)then
1809           jec=1
1810           evk(i,j,l+27*(k-1))=evs(i,j,l+27*(k-1))
1811         endif
1812       enddo
1813       enddo
1814       enddo
1815       enddo
1816
1817       if(jec.ne.0)goto 1
1818
1819       write(ifmt,'(a)')'write to iniev ...'
1820       open(1,file=fnie(1:nfnie),status='unknown')
1821       write (1,*)qcdlam,q2min,q2ini,naflav,epmax
1822       write (1,*)evk0,evk
1823       close(1)
1824
1825 101   continue
1826       return
1827       end
1828
1829 c------------------------------------------------------------------------
1830       function psdbom(s,t,u,qq,long)
1831 c-----------------------------------------------------------------------
1832 c psdbom - integrand for DIS c-quark cross-sections (matrix element squared)
1833 c s  - total c.m. energy squared for the scattering (for n=2: s+qq),
1834 c t  - invariant variable for the scattering |(p1-p3)**2|
1835 c u  - invariant variable for the scattering |(p1-p4)**2|
1836 c qq - photon virtuality
1837 c long: 0 - contr. to (F2-F_L), 1 - contr. to F_L
1838 c-----------------------------------------------------------------------
1839       include 'epos.incsem'
1840       if(long.eq.0)then       !F2-F_L
1841         psdbom=(2.*(t/u+u/t)*(qq**2+(s-qq)**2)/s**2+
1842      *  4.*(qcmass*s/t/u)**2*(qq-2.*qcmass**2)+
1843      *  8.*qcmass**2/t/u*(s-2.*qq))   *2.    !=4.5/2.25
1844       else                    !F_L_C
1845         psdbom=16.*qq*((s-qq)/s**2-qcmass**2/t/u) *2.  !=4.5/2.25
1846       endif
1847       return
1848       end
1849
1850 c------------------------------------------------------------------------
1851       function psdbin(q1,qq,s,m1,long)
1852 c-----------------------------------------------------------------------
1853 c psdbin - DIS born cross-section
1854 c q1      - virtuality cutoff for current end of the ladder
1855 c qq      - photon virtuality
1856 c s=2(pq) - s_true + qq
1857 c s2min   - mass cutoff for born scattering
1858 c m1       - incoming parton type (0 - g, 1,2 - q)
1859 c-----------------------------------------------------------------------
1860       double precision xx
1861       include 'epos.incsem'
1862       include 'epos.inc'
1863
1864       psdbin=0.
1865       q2mass=qcmass**2
1866       s2min=4.*max(q1,q2mass)+qq
1867       if(m1.eq.0.and.s.gt.s2min.and.(idisco.eq.0.or.idisco.eq.2))then
1868         tmax=s/2.
1869         qtq=4.*max(q2mass,q1)/(s-qq)
1870         if(qtq.lt.1.)then
1871           tmin=.5*s*qtq/(1.+sqrt(1.-qtq))
1872         else
1873           tmin=.5*s
1874         endif
1875         psdbin=psdbin+psdbor(q1,qq,s,long)*(1./tmin-1./tmax)
1876       endif
1877
1878       if(long.eq.0.and.q1.lt.qq.and.s.gt.qq/(1.-q2ini/qq)
1879      *.and.(idisco.eq.0.or.idisco.eq.1))then
1880         m=min(1,iabs(m1))+1
1881         xx=qq/s
1882         psdbin=psdbin+psevi0(q1,qq,xx,m,2)*4.*pi**2*alfe/s
1883       endif
1884       return
1885       end
1886
1887 c------------------------------------------------------------------------
1888       function psdbor(q1,qq,s,long)
1889 c-----------------------------------------------------------------------
1890 c psdbor - DIS born cross-section
1891 c q1      - virtuality cutoff for current end of the ladder
1892 c qq      - photon virtuality
1893 c s=2(pq) - s_true + qq
1894 c s2min   - mass cutoff for born scattering
1895 c-----------------------------------------------------------------------
1896       common /ar3/    x1(7),a1(7)
1897       include 'epos.inc'
1898       include 'epos.incsem'
1899       double precision psuds
1900
1901       psdbor=0.
1902       q2mass=qcmass**2
1903       qtq=4.*max(q2mass,q1)/(s-qq)
1904       j=0   !Gluon
1905
1906       tmax=s/2.
1907       if(qtq.lt.1.)then
1908         tmin=.5*s*qtq/(1.+sqrt(1.-qtq))
1909       else
1910         tmin=.5*s
1911       endif
1912       if(tmax.lt.tmin.and.ish.ge.1)write(ifmt,*)'s,q1,qq,tmin,tmax',
1913      *s,q1,qq,tmin,tmax
1914
1915       ft=0.
1916       do i=1,7
1917       do m=1,2
1918         t=2.*tmin/(1.+tmin/tmax+(2*m-3)*x1(i)*(1.-tmin/tmax))
1919         u=s-t
1920
1921         qt=t*u/s*(1.-qq/s)
1922         if(qt.lt..999*max(q2mass,q1).and.ish.ge.1)
1923      &  write(ifmt,*)'psdbor:qt,q1',qt,q1
1924         fb=psdbom(s,t,u,qq,long)*t**2
1925         ft=ft+a1(i)*fb*pssalf(qt/qcdlam)*sngl(psuds(qt,j))
1926       enddo
1927       enddo
1928       psdbor=ft/s**2*pi**2*alfe/sngl(psuds(q1,j))
1929       return
1930       end
1931
1932 c------------------------------------------------------------------------
1933       subroutine psdint(s,qq,sds,sdn,sdb,sdt,sdr,m1,long)
1934 c-----------------------------------------------------------------------
1935 c psdint - dis cross-sections interpolation - for minimal
1936 c effective momentum cutoff in the ladder
1937 c s   - total c.m. energy squared for the ladder,
1938 c qq  - photon virtuality,
1939 c sds - dis singlet cross-section,
1940 c sdn - dis nonsinglet cross-section,
1941 c sdb - dis born cross-section,
1942 c sdt - dis singlet+resolved cross-section,
1943 c m1  - parton type at current end of the ladder (0 - g, 1,2 - q)
1944 c-----------------------------------------------------------------------
1945       double precision xx
1946       dimension wk(3),wj(3)
1947       common /psar2/  edmax,epmax
1948       common /psar27/ csds(21,26,4),csdt(21,26,2),csdr(21,26,2)
1949       include 'epos.incsem'
1950       include 'epos.inc'
1951
1952       sds=0.
1953       sdn=0.
1954       sdt=0.
1955       sdr=0.
1956       sdb=psdbin(q2min,qq,s,m1,long)
1957
1958       m=min(1,iabs(m1))+1
1959       qlj=log(qq/q2min)*2.+1.
1960       j=int(qlj)
1961       if(j.lt.1)j=1
1962       if(j.gt.19)j=19
1963       wj(2)=qlj-j
1964       wj(3)=wj(2)*(wj(2)-1.)*.5
1965       wj(1)=1.-wj(2)+wj(3)
1966       wj(2)=wj(2)-2.*wj(3)
1967
1968       s2min=4.*max(q2min,qcmass**2)+qq
1969       if(m1.ne.0)s2min=s2min/(1.-4.*q2ini/(s2min-qq))
1970       if(s.le.s2min.or.idisco.ne.0.and.idisco.ne.2)goto 1
1971
1972       qtq=4.*max(q2min,qcmass**2)/(s-qq)
1973       if(qtq.lt.1.)then
1974         tmin=.5*s*qtq/(1.+sqrt(1.-qtq))
1975       else
1976         tmin=.5*s
1977       endif
1978       tmax=s/2.
1979
1980       sl=log(s/s2min)/log(edmax/s2min)*25.+1.
1981       k=int(sl)
1982       if(k.lt.1)k=1
1983       if(k.gt.24)k=24
1984       wk(2)=sl-k
1985       wk(3)=wk(2)*(wk(2)-1.)*.5
1986       wk(1)=1.-wk(2)+wk(3)
1987       wk(2)=wk(2)-2.*wk(3)
1988
1989       do k1=1,3
1990         k2=k+k1-1
1991       do j1=1,3
1992         sds=sds+csds(j+j1-1,k2,m+2*long)*wj(j1)*wk(k1)
1993       enddo
1994       enddo
1995       if(m.eq.1)then
1996         sds=exp(sds)*(1./tmin-1./tmax)
1997       else
1998         sds=max(sds,0.)
1999       endif
2000
2001 1     continue
2002       s2min=max(4.*qq,16.*q2min)+qq
2003       if(s.le.s2min.or.long.ne.0.or.idisco.ne.0.and.idisco.ne.3)then
2004         sdt=sds
2005         goto 2
2006       endif
2007
2008       sl=log(s/s2min)/log(edmax/s2min)*25.+1.
2009       k=int(sl)
2010       if(k.lt.1)k=1
2011       if(k.gt.24)k=24
2012       wk(2)=sl-k
2013       wk(3)=wk(2)*(wk(2)-1.)*.5
2014       wk(1)=1.-wk(2)+wk(3)
2015       wk(2)=wk(2)-2.*wk(3)
2016
2017       do k1=1,3
2018         k2=k+k1-1
2019       do j1=1,3
2020         sdr=sdr+csdr(j+j1-1,k2,m)*wj(j1)*wk(k1)
2021         sdt=sdt+csdt(j+j1-1,k2,m)*wj(j1)*wk(k1)
2022       enddo
2023       enddo
2024
2025       sdr=max(sdr,0.)
2026       sdt=max(sds,sds+sdt)
2027       sdt=sdt+sdr
2028
2029 2     continue
2030       if(long.eq.0.and.q2min.lt.qq.and.s.gt.qq/(1.-q2ini/qq)
2031      *.and.(idisco.eq.0.or.idisco.eq.1))then
2032         xx=qq/s
2033         dsi=psevi(q2min,qq,xx,m,2)*4.*pi**2*alfe/s
2034         if(m1.eq.0)then
2035           sds=sds+dsi
2036           sdt=sdt+dsi
2037         else
2038           dnsi=psevi(q2min,qq,xx,3,2)*4.*pi**2*alfe/s
2039           sdn=sdn+dnsi
2040           sds=sds+max(dsi-dnsi,0.)
2041           sdt=sdt+max(dsi-dnsi,0.)
2042         endif
2043       endif
2044
2045       if(m1.eq.0)then
2046         sds=max(sds,sdb)
2047         sdt=max(sdt,sdb)
2048       else
2049         sdn=max(sdn,sdb)
2050       endif
2051       return
2052       end
2053
2054 c-----------------------------------------------------------------------
2055       function psdnsi(q1,qq,s,long)
2056 c-----------------------------------------------------------------------
2057 c psdnsi - DIS nonsinglet cross-section interpolation
2058 c q1 - effective momentum cutoff for current end of the ladder,
2059 c qq - photon virtuality,
2060 c s - total c.m. energy squared for the ladder,
2061 c-----------------------------------------------------------------------
2062       double precision xx
2063       include 'epos.incsem'
2064       include 'epos.inc'
2065
2066       psdnsi=0.
2067       if(long.eq.0.and.q1.lt.qq.and.s.gt.qq/(1.-q2ini/qq))then
2068         xx=qq/s
2069         psdnsi=psdnsi+max(0.,psevi(q1,qq,xx,3,2)*4.*pi**2*alfe/s)
2070       endif
2071       return
2072       end
2073
2074 c-----------------------------------------------------------------------
2075       function psdrga(qq,s,s2min,j)
2076 c-----------------------------------------------------------------------
2077 c psdrga - DIS resolved cross-section (photon sf)
2078 c qq    - photon virtuality
2079 c s     - total c.m. energy squared for the process
2080 c s2min - mass cutoff for born scattering
2081 c j     - parton type at current end of the ladder (0 - g, 1,2 etc. - q)
2082 c-----------------------------------------------------------------------
2083       common /ar3/   x1(7),a1(7)
2084       include 'epos.incsem'
2085
2086       psdrga=0.
2087       if(s.le.s2min)return
2088
2089       xmin=s2min/s
2090       do i=1,7
2091       do m=1,2
2092         z=xmin**(.5+(m-1.5)*x1(i))
2093         tu=psdfh4(z,q2min,qq,0,1)
2094         td=psdfh4(z,q2min,qq,0,2)
2095         ts=psdfh4(z,q2min,qq,0,3)
2096         tg=psdfh4(z,q2min,qq,0,0)
2097         if(j.eq.0)then
2098           sj=tg*psjti(q2min,qq,z*s,0,j,1)+
2099      *    (tu+td+ts)*psjti(q2min,qq,z*s,1,j,1)
2100         else
2101           sj=tg*psjti(q2min,qq,z*s,0,j,1)+
2102      *    (tu+td)*(psjti(q2min,qq,z*s,1,1,1)/4.+
2103      *    psjti(q2min,qq,z*s,-1,1,1)/4.+
2104      *    psjti(q2min,qq,z*s,2,1,1)/2.)+
2105      *    ts*psjti(q2min,qq,z*s,2,1,1)
2106         endif
2107         psdrga=psdrga+a1(i)*sj
2108       enddo
2109       enddo
2110       psdrga=-psdrga*log(xmin)*alfe/2.  *4.5 !mean e^2 is taken out
2111       return
2112       end
2113
2114 c-----------------------------------------------------------------------
2115       function psdres(qq,s,s2min,j)
2116 c-----------------------------------------------------------------------
2117 c psdres - DIS resolved photon cross-section
2118 c qq    - photon virtuality
2119 c s     - total w squared for the ladder (s+qq)
2120 c s2min - mass cutoff for born scattering
2121 c j     - parton type at current end of the ladder (0 - g, 1,2 etc. - q)
2122 c-----------------------------------------------------------------------
2123       double precision xx
2124       common /ar3/   x1(7),a1(7)
2125       include 'epos.inc'
2126       include 'epos.incsem'
2127
2128       psdres=0.
2129       if(s.le.s2min+qq)return
2130
2131       qmin=max(q2min,s2min/(s/qq-1.))
2132       qmax=min(s-s2min,s/2.)
2133
2134 c numerical integration over transverse momentum squared;
2135 c gaussian integration is used
2136       do i=1,7
2137       do m=1,2
2138         qi=2.*qmin/(1.+qmin/qmax+(2*m-3)*x1(i)*(1.-qmin/qmax))
2139
2140         zmax=min(1.,qi/qq)
2141         zmin=(max(qi,s2min)+qi)/s
2142
2143         fsj=0.
2144         if(zmax.gt.zmin)then
2145           do i1=1,7
2146           do m1=1,2
2147             z=.5*(zmax+zmin+(2*m1-3)*x1(i1)*(zmax-zmin))
2148             s2=z*s-qi
2149             xx=z
2150             if(j.eq.0)then
2151               sj=psfap(xx,0,1)*psjti(qi,qq,s2,1,j,1)
2152             else
2153               sj=psfap(xx,0,1)*(psjti(qi,qq,s2,1,1,1)/6.+
2154      *        psjti(qi,qq,s2,-1,1,1)/6.+
2155      *        psjti(qi,qq,s2,2,1,1)/1.5)
2156             endif
2157             fsj=fsj+a1(i1)*sj*qi    !????????(qi-z*qq)
2158           enddo
2159           enddo
2160           fsj=fsj*(zmax-zmin)
2161         elseif(ish.ge.1)then
2162           write(ifmt,*)'psdres:zmax,zmin',zmax,zmin
2163         endif
2164         psdres=psdres+a1(i)*fsj
2165       enddo
2166       enddo
2167       psdres=psdres*(1./qmin-1./qmax)*alfe*.75/pi  !alpha_s -> 6 alpha_e
2168       return
2169       end
2170
2171 c------------------------------------------------------------------------
2172       function psds(q1,qq,s,j,long)
2173 c-----------------------------------------------------------------------
2174 c psds - DIS singlet cross-section
2175 c q1      - virtuality cutoff for current end of the ladder
2176 c qq      - photon virtuality
2177 c s=2(pq) - s_true + qq
2178 c s2min   - mass cutoff for born scattering
2179 c-----------------------------------------------------------------------
2180       double precision xxe,xmax,xmin,xmax1,xmin1
2181       common /ar3/    x1(7),a1(7)
2182       include 'epos.inc'
2183       include 'epos.incsem'
2184
2185       psds=0.
2186       q2mass=qcmass**2
2187       s2min=4.*max(q1,q2mass)
2188       smin=(s2min+qq)/(1.-4.*q2ini/s2min)
2189       if(s.le.1.001*smin)return
2190
2191       xmax=.5d0*(1.d0+qq/s+dsqrt((1.d0-qq/s)**2-16.d0*q2ini/s))
2192       xmin=max(1.d0+qq/s-xmax,1.d0*(s2min+qq)/s)
2193       if(xmin.gt.xmax.and.ish.ge.1)write(ifmt,*)'xmin,xmax,q1,qq,s,smin'
2194      *,xmin,xmax,q1,qq,s,smin
2195
2196       fx1=0.
2197       fx2=0.
2198       if(xmax.gt..9d0)then
2199         xmin1=max(xmin,.9d0)
2200         do i=1,7
2201         do m=1,2
2202           xxe=1.d0-(1.d0-xmax)*((1.d0-xmin1)/(1.d0-xmax))**
2203      *    (.5d0-x1(i)*(m-1.5))
2204           xx=xxe
2205
2206           sh=xx*s
2207           qtmin=max(1.d0*max(q2mass,q1),q2ini/(1.d0-xxe))
2208           qtq=4.*qtmin/(sh-qq)
2209
2210           tmin=.5*sh*qtq/(1.+sqrt(1.-qtq))
2211           tmax=.5*sh
2212           if(tmin.gt.tmax.and.ish.ge.1)write(ifmt,*)'psds:tmin,tmax'
2213      &                                              ,tmin,tmax
2214
2215           ft=0.
2216           do i1=1,7
2217           do m1=1,2
2218             t=.5*(tmin+tmax+(2*m1-3)*x1(i1)*(tmin-tmax))
2219             u=sh-t
2220             qt=t*u/sh*(1.-qq/sh)
2221             if(qt.lt.qtmin.and.ish.ge.1)write(ifmt,*)'psds:qt,qtmin'
2222      &                                               ,qt,qtmin
2223
2224             fb=psdsj(q1,xxe,sh,qt,t,u,qq,j,long)
2225             ft=ft+a1(i1)*fb*pssalf(qt/qcdlam)
2226           enddo
2227           enddo
2228           ft=ft*(tmax-tmin)
2229           fx1=fx1+a1(i)*ft*(1.-xx)/sh**2
2230         enddo
2231         enddo
2232         fx1=fx1*log((1.d0-xmin1)/(1.d0-xmax))
2233       endif
2234
2235       if(xmin.lt..9d0)then
2236         xmax1=min(xmax,.9d0)
2237         do i=1,7
2238         do m=1,2
2239           xxe=xmin*(xmax1/xmin)**(.5-x1(i)*(m-1.5))
2240           xx=xxe
2241
2242           sh=xx*s
2243           qtmin=max(1.d0*max(q2mass,q1),q2ini/(1.d0-xxe))
2244           qtq=4.*qtmin/(sh-qq)
2245
2246           tmin=.5*sh*qtq/(1.+sqrt(1.-qtq))
2247           tmax=.5*sh
2248           if(tmin.gt.tmax.and.ish.ge.1)write(ifmt,*)'psds:tmin,tmax'
2249      *                                              ,tmin,tmax
2250
2251           ft=0.
2252           do i1=1,7
2253           do m1=1,2
2254             t=(.5*(tmin+tmax+(2*m1-3)*x1(i1)*
2255      *      (tmin-tmax)))
2256             u=sh-t
2257             qt=t*u/sh*(1.-qq/sh)
2258             if(qt.lt.qtmin.and.ish.ge.1)write(ifmt,*)'psds:qt,qtmin'
2259      *                                               ,qt,qtmin
2260
2261             fb=psdsj(q1,xxe,sh,qt,t,u,qq,j,long)
2262             ft=ft+a1(i1)*fb*pssalf(qt/qcdlam)
2263           enddo
2264           enddo
2265           ft=ft*(tmax-tmin)
2266           fx2=fx2+a1(i)*ft*xx/sh**2
2267         enddo
2268         enddo
2269         fx2=fx2*log(xmax1/xmin)
2270       endif
2271       psds=(fx1+fx2)*pi**2*alfe
2272       return
2273       end
2274
2275 c-----------------------------------------------------------------------
2276       function psdsj(q1,xx,s,qt,t,u,qq,j,long)
2277 c-----------------------------------------------------------------------
2278 c psdsj - integrand for dis singlet cross-section
2279 c q1 - virtuality cutoff for current end of the ladder
2280 c xx - lc momentum ratio between initial (j) and final (l) partons
2281 c s  - c.m. energy squared for the born scattering,
2282 c t  - invariant variable for the born scattering |(p1-p3)**2|
2283 c u  - invariant variable for the born scattering |(p1-p4)**2|
2284 c qq - photon virtuality
2285 c j  - initial parton at the end of the ladder (0 - g, 1,2 - q)
2286 c-----------------------------------------------------------------------
2287       double precision xx
2288       include 'epos.incsem'
2289
2290       fb=psdbom(s,t,u,qq,long)
2291       psdsj=psevi(q1,qt,xx,min(1,iabs(j))+1,1)*fb
2292       return
2293       end
2294
2295 c------------------------------------------------------------------------
2296       function psdsin(q1,qq,s,m1,long)
2297 c-----------------------------------------------------------------------
2298 c psdsin - DIS singlet cross-section interpolation
2299 c q1 - effective momentum cutoff for current end of the ladder,
2300 c qq - photon virtuality,
2301 c s -  total c.m. energy squared for the ladder,
2302 c m1 - parton type at current end of the ladder (0 - g, 1,2 - q)
2303 c-----------------------------------------------------------------------
2304       double precision xx
2305       dimension wi(3),wj(3),wk(3)
2306       common /psar2/  edmax,epmax
2307       common /psar25/ csdsi(21,21,104)
2308       include 'epos.incsem'
2309       include 'epos.inc'
2310
2311       psdsin=0.
2312       m=min(1,iabs(m1))+1
2313
2314       q2mass=qcmass**2
2315       s2min=4.*max(q2min,q2mass)+qq
2316       sdmin=4.*max(q1,q2mass)+qq
2317       if(m1.ne.0)then
2318         s2min=s2min/(1.-4.*q2ini/(s2min-qq))
2319         sdmin=sdmin/(1.-4.*q2ini/(sdmin-qq))
2320       endif
2321 c      if(s.le.1.e8*sdmin)goto 2  !????????????????
2322       if(s.le.sdmin)goto 2
2323
2324       qmin=q2min
2325       qmax=(s-qq)/4.
2326       if(m1.ne.0)qmax=(s-qq+sqrt((s-qq)**2-16.*s*q2ini))/8.
2327       qtq=4.*max(q2mass,q1)/(s-qq)
2328       if(qtq.lt.1.)then
2329         tmin=.5*s*qtq/(1.+sqrt(1.-qtq))
2330       else
2331         tmin=.5*s
2332       endif
2333       tmax=s/2.
2334
2335       qlj=log(qq/q2min)*2.+1.
2336       j=int(qlj)
2337       if(j.lt.1)j=1
2338       if(j.gt.19)j=19
2339       wj(2)=qlj-j
2340       wj(3)=wj(2)*(wj(2)-1.)*.5
2341       wj(1)=1.-wj(2)+wj(3)
2342       wj(2)=wj(2)-2.*wj(3)
2343
2344       qli=log(q1/qmin)/log(qmax/qmin)*20.+1.
2345       i=int(qli)
2346       if(i.lt.1)i=1
2347       if(i.gt.19)i=19
2348       wi(2)=qli-i
2349       wi(3)=wi(2)*(wi(2)-1.)*.5
2350       wi(1)=1.-wi(2)+wi(3)
2351       wi(2)=wi(2)-2.*wi(3)
2352
2353       sl=log(s/s2min)/log(edmax/s2min)*25.+1.
2354       k=int(sl)
2355       if(k.lt.1)k=1
2356       if(k.gt.24)k=24
2357       wk(2)=sl-k
2358       wk(3)=wk(2)*(wk(2)-1.)*.5
2359       wk(1)=1.-wk(2)+wk(3)
2360       wk(2)=wk(2)-2.*wk(3)
2361
2362       dsin1=0.
2363       do k1=1,3
2364         k2=k+k1-1+26*(m-1)+52*long
2365       do i1=1,3
2366       do j1=1,3
2367         dsin1=dsin1+csdsi(i+i1-1,j+j1-1,k2)*wi(i1)*wj(j1)*wk(k1)
2368       enddo
2369       enddo
2370       enddo
2371       if(m1.eq.0)then
2372         psdsin=psdsin+exp(dsin1)*(1./tmin-1./tmax)
2373       else
2374         psdsin=psdsin+max(0.,dsin1)
2375       endif
2376
2377 2     continue
2378       if(long.eq.0.and.q1.lt.qq.and.s.gt.qq/(1.-q2ini/qq))then
2379         xx=qq/s
2380         dsi=psevi(q1,qq,xx,m,2)*4.*pi**2*alfe/s
2381         if(m1.eq.0)then
2382           psdsin=psdsin+max(dsi,0.)
2383         else
2384           dnsi=psevi(q1,qq,xx,3,2)*4.*pi**2*alfe/s
2385           psdsin=psdsin+max(dsi-dnsi,0.)
2386         endif
2387       endif
2388       return
2389       end
2390
2391
2392
2393
2394
2395
2396
2397
2398
2399
2400
2401
2402
2403
2404 c###########################################################################
2405 c###########################################################################
2406 c###########################################################################
2407 c###########################################################################
2408 c
2409 c                        unused 3P
2410 c
2411 c###########################################################################
2412 c###########################################################################
2413 c###########################################################################
2414 c###########################################################################
2415
2416
2417 c------------------------------------------------------------------------
2418       function psvy(xpp0,xpr0,xpm0,xmr0,b,iqq)
2419 c-----------------------------------------------------------------------
2420 c psvy - 3p-contributions to the interaction eikonal
2421 c xpp - lc+ for the pomeron,
2422 c xpr - lc+ for the remnant,
2423 c xpm - lc- for the pomeron,
2424 c xpr - lc- for the remnant,
2425 c b   - impact parameter,
2426 c iqq=1  - Y-proj-uncut
2427 c iqq=2  - Y-proj-1-cut
2428 c iqq=3  - Y-proj-2-cut
2429 c iqq=4  - Y-proj-soft-cut
2430 c iqq=5  - Y-proj-gss-cut
2431 c iqq=6  - Y-proj-qss-cut
2432 c iqq=7  - Y-proj-ssg-cut
2433 c iqq=8  - Y-proj-ssq-cut
2434 c iqq=9  - Y-proj-difr
2435 c iqq=-1 - Y-targ-uncut
2436 c iqq=-2 - Y-targ-1-cut
2437 c iqq=-3 - Y-targ-2-cut
2438 c iqq=-4 - Y-targ-soft-cut
2439 c iqq=-5 - Y-targ-gss-cut
2440 c iqq=-6 - Y-targ-qss-cut
2441 c iqq=-7 - Y-targ-ssg-cut
2442 c iqq=-8 - Y-targ-ssq-cut
2443 c iqq=-9 - Y-targ-difr
2444 c------------------------------------------------------------------------
2445
2446       psvy=0.
2447       return
2448       end
2449
2450 c------------------------------------------------------------------------
2451       function psvx(xpp,xpr,xpm,xmr,b,iqq)
2452 c-----------------------------------------------------------------------
2453 c psvx - 4p-contributions to the interaction eikonal
2454 c xpp - lc+ for the pomeron,
2455 c xpr - lc+ for the remnant,
2456 c xpm - lc- for the pomeron,
2457 c xpr - lc- for the remnant,
2458 c b   - impact parameter,
2459 c iqq=0  - X-uncut
2460 c iqq=1  - X-1-cut
2461 c iqq=2  - X-Y+cut
2462 c iqq=-2 - X-Y-cut
2463 c iqq=3  - X-1-cut-soft
2464 c iqq=4  - X-1-cut-gss
2465 c iqq=-4 - X-1-cut-gss
2466 c iqq=5  - X-1-cut-qss
2467 c iqq=-5 - X-1-cut-qss
2468 c iqq=6  - X-difr+
2469 c iqq=-6 - X-difr-
2470 c------------------------------------------------------------------------
2471
2472       psvx=0.
2473       return
2474       end
2475
2476
2477
2478
2479
2480 c
2481 c
2482 c
2483 c
2484 c
2485 cc------------------------------------------------------------------------
2486 c      function psftig(x,xpomr,jj)
2487 cc-----------------------------------------------------------------------
2488 c
2489 c        psftig=0.
2490 c      end
2491 c
2492 c
2493 cc------------------------------------------------------------------------
2494 c      function psftih(zz,ddd)
2495 cc-----------------------------------------------------------------------
2496 c
2497 c      psftih=0.
2498 c      return
2499 c      end
2500 c
2501 cc------------------------------------------------------------------------
2502 c      function psftij(zz,ddd)
2503 cc------------------------------------------------------------------------
2504 c
2505 c      psftij=0.
2506 c      return
2507 c      end
2508 c
2509 cc------------------------------------------------------------------------
2510 c      function psftik(xp,del1,rh1,rh2,rh3,rp,z)
2511 cc-----------------------------------------------------------------------
2512 c
2513 c      psftik=0.
2514 c      return
2515 c      end
2516 c
2517 cc------------------------------------------------------------------------
2518 c      function psftim(xp1,xp,del1,rh1,rh2,rh3,rp,z)
2519 cc-----------------------------------------------------------------------
2520 c
2521 c      psftim=0.
2522 c      return
2523 c      end
2524 c
2525 cc------------------------------------------------------------------------
2526 c      function psftil(xp,del1,rh1,rh2,rh3,rp,z)
2527 cc------------------------------------------------------------------------
2528 c
2529 c      psftil=0.
2530 c      return
2531 c      end
2532 c
2533 cc------------------------------------------------------------------------
2534 c      function psftigt(x)
2535 cc-----------------------------------------------------------------------
2536 c
2537 c      psftigt=0.
2538 c       return
2539 c       end
2540 c
2541 cc------------------------------------------------------------------------
2542 c      function psftist(x)
2543 cc-----------------------------------------------------------------------
2544 c      psftist=0.
2545 c       return
2546 c       end
2547 c
2548 cc------------------------------------------------------------------------
2549 c      function psftig1(x,xpomr,jj)
2550 cc-----------------------------------------------------------------------
2551 c      psftig1=0.
2552 c      return
2553 c      end
2554 c
2555 cc------------------------------------------------------------------------
2556 c      function psftis1(x,xpomr,jj)
2557 cc-----------------------------------------------------------------------
2558 c      psftis1=0.
2559 c      return
2560 c      end
2561 c
2562 cc------------------------------------------------------------------------
2563 c      function psd3p1(xd,xpomr,qq,jj)
2564 cc-----------------------------------------------------------------------
2565 cc psd3p1 - df2difr/dx_pomr
2566 cc xd - bjorken x,
2567 cc xpomr - pomeron x,
2568 cc qq - photon virtuality
2569 cc jj=1 - 1st order
2570 cc-----------------------------------------------------------------------
2571 c
2572 c      psd3p1=0.
2573 c
2574 c      return
2575 c      end
2576 c
2577 cc------------------------------------------------------------------------
2578 c      function psv3p(sy,xpp,xpm,zb)
2579 cc-----------------------------------------------------------------------
2580 cc psv3p - 3p-contributions to the interaction eikonal
2581 cc sy  - energy squared for the hard interaction,
2582 cc xpp - lc+ for the sh pomeron,
2583 cc xpm - lc- for the sh pomeron,
2584 cc z   - impact parameter factor, z=exp(-b**2/4*rp),
2585 cc------------------------------------------------------------------------
2586 c
2587 c      psv3p=0.
2588 c      return
2589 c      end
2590 c
2591 cc------------------------------------------------------------------------
2592 c      function psfro(xpomr,zb,iclp,icdpro)
2593 cc-----------------------------------------------------------------------
2594 cc psfro - generalized froissaron between proj. and 3p-vertex
2595 cc xpp   - lc+ for the proj. side,
2596 cc xpomr - lc+ for the vertex,
2597 cc zb    - impact parameter factor, z=exp(-b**2/4*rp),
2598 cc------------------------------------------------------------------------
2599 c      psfro=0.
2600 c      return
2601 c      end
2602 c
2603 cc------------------------------------------------------------------------
2604 c      function psv2(xpomr,zb,iclp,iqq)
2605 cc-----------------------------------------------------------------------
2606 cc psv2 - 2-pom contribution to the froissaron
2607 cc xpomr - lc+ for the vertex,
2608 cc zb    - impact parameter factor, z=exp(-b**2/4*rp),
2609 cc------------------------------------------------------------------------
2610 c
2611 c      psv2=0.
2612 c      return
2613 c      end
2614 c
2615 cc------------------------------------------------------------------------
2616 c      function psvfro(xpp,xpr,xpomr,zb,iclp,icdpro,iqq)
2617 cc-----------------------------------------------------------------------
2618 cc psvfro - effective froissaron contributions
2619 cc xpomr - lc+ for the vertex,
2620 cc zb    - impact parameter factor, z=exp(-b**2/4*rp),
2621 cc iqq=0 - total uncut
2622 cc iqq=1 - total 1-cut
2623 cc iqq=2 - total 2-cut
2624 cc iqq=3 - soft 1-cut
2625 cc iqq=4 - gg 1-cut
2626 cc iqq=5 - qg 1-cut
2627 cc iqq=6 - difr
2628 cc------------------------------------------------------------------------
2629 c
2630 c      psvfro=0.
2631 c      return
2632 c      end
2633 c
2634 cc------------------------------------------------------------------------
2635 c      function psfroin(xpomr,z,iclp,icdpro)
2636 cc-----------------------------------------------------------------------
2637 cc psfroin - interpolation of effective froissaron corrections
2638 cc xpomr - lc+ for the 3p-vertex,
2639 cc z   - impact parameter factor, z=exp(-b**2/4*rp),
2640 cc-----------------------------------------------------------------------
2641 c
2642 c      psfroin=0.
2643 c      return
2644 c      end
2645 c
2646 cc------------------------------------------------------------------------
2647 c      function psvnorm(b)
2648 cc-----------------------------------------------------------------------
2649 cc psvnorm - X-contribution normalization
2650 cc b   - impact parameter
2651 cc-----------------------------------------------------------------------
2652 c
2653 c      psvnorm=0.
2654 c      return
2655 c      end
2656 c
2657 cc------------------------------------------------------------------------
2658 c      function psvxb(dxpp,dxpr,dxpm,dxmr,dxpomr,bb1,bb2
2659 c     *,iclp,iclt,icdpro,icdtar,iqq)
2660 cc-----------------------------------------------------------------------
2661 cc psvxb   - integrand for X-contributions
2662 cc dxpomr  - lc+ for the vertex,
2663 cc bb1,bb2 - impact parameters to proj(targ),
2664 cc iqq=0   - uncut
2665 cc iqq=1   - 1-cut
2666 cc iqq=2   - Y-cut
2667 cc iqq=3   - soft 1-cut
2668 cc iqq=4   - gg 1-cut
2669 cc iqq=5   - qg 1-cut
2670 cc iqq=6   - difr
2671 cc------------------------------------------------------------------------
2672 c      double precision dxpp,dxpr,dxpm,dxmr,dxpomr
2673 c
2674 c      psvxb=0.
2675 c      return
2676 c      end
2677 c
2678 cc------------------------------------------------------------------------
2679 c      function pscoef(zz,alp1,alp2,iclp,iqq)
2680 cc-----------------------------------------------------------------------
2681 cc pscoef - integrated vertexes
2682 cc zz=xpomr/xpr(xpp)
2683 cc iqq=0 - 2-cut
2684 cc iqq=1 - 1-uncut
2685 cc iqq=2 - 2-uncut
2686 cc------------------------------------------------------------------------
2687 c
2688 c      pscoef=0.
2689 c      return
2690 c      end
2691 c
2692 cc------------------------------------------------------------------------
2693 c      function pscoefi(z,i1,i2,iclp,iqq)
2694 cc-----------------------------------------------------------------------
2695 cc pscoefi - interpolation of integrated vertexes
2696 cc z=xpomr/xpr(xpp)
2697 cc iqq=0 - 2-cut
2698 cc iqq=1 - 1-uncut
2699 cc iqq=2 - 2-uncut
2700 cc------------------------------------------------------------------------
2701 c      pscoefi=0.
2702 c      return
2703 c      end
2704 c
2705 c