Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2021-02-14 13:07:29

0001 c  reshuffled from sem, sto, sha
0002 
0003 c    contains psahot and related stuff
0004 c             ------
0005 
0006 
0007 
0008 c-----------------------------------------------------------------------
0009       subroutine psahot(kcol,ncolp,iret)
0010 c-----------------------------------------------------------------------
0011 c psahot - showering (semihard parton-parton interaction)
0012 c-----------------------------------------------------------------------
0013       include 'epos.inc'
0014       include 'epos.incsem'
0015       include 'epos.incems'
0016       include 'epos.incpar'
0017       double precision ept(4),ept1(4),xx,wpt(2),eprt,pl,plprt,wplc
0018      *,wp0,wm0,s,si,smin,xmin,xp1,wpi,wmi,xp,xm,wp1,wm1,wp2,wm2
0019      *,wpq,wmq,wpi1,wmi1,pxh,pyh,pth,xmm,xp2,xg,zg,smax,xmax
0020      *,xmax1,xmin1,zp0,psutz,xpmax0,xpmin0,gb0,tmax0,tmin0,zpm,gb
0021      *,gbyj,tmax,tmin,t,gb7,x,s2,discr,qt,qt2,x1min,x1max,t1min,t1max
0022      *,t1,xq1,qq,qqmax,qqmin,pt2,xmin2,ptnew,xpmin,xpmax,xm0,psuds
0023 c     *,xprh,xmrh
0024       dimension ep3(4),ey(3),bx(6)
0025      *,qmin(2),iqc(2),nqc(2),ncc(2,2),amqt(4)
0026       parameter (mjstr=20000)
0027       common /psar7/  delx,alam3p,gam3p
0028       common /psar29/eqj(4,mjstr),iqj(mjstr),ncj(2,mjstr),ioj(mjstr),nj
0029       common /psar30/iorj(mjstr),ityj(mjstr),bxj(6,mjstr),q2j(mjstr)
0030       common /testj/  ajeth(4),ajete(5),ajet0(7)
0031       parameter (ntim=1000)
0032       common/cprt/pprt(5,ntim),q2prt(ntim),idaprt(2,ntim),idprt(ntim)
0033      &,iorprt(ntim),jorprt(ntim),nprtj
0034       common/emsptl/nppr(npommx,kollmx),npproj(mamx),nptarg(mamx)
0035       integer icp(2),ict(2),nemis(2)
0036       integer icp1(2),icp2(2),icm1(2),icm2(2)
0037       integer jcp(nflav,2),jct(nflav,2),jcpr(nflav,2),jctr(nflav,2)
0038       common/cprtx/nprtjx,pprtx(5,2)/ciptl/iptl
0039 
0040       call utpri('psahot',ish,ishini,3)
0041 
0042       iret=0
0043       alpq=-(alppar+1.)/2.
0044       qqcut=q2min  !????????????pt2cut
0045 
0046 
0047       nptl1=nptl
0048       iptl=nppr(ncolp,kcol)
0049       ip=iproj(kcol)
0050       it=itarg(kcol)
0051       do i=1,2
0052         icp(i)=icproj(i,ip)
0053         ict(i)=ictarg(i,it)
0054       enddo
0055       idpomr=idhpr(ncolp,kcol)
0056       bpomr=bhpr(ncolp,kcol)
0057 
0058       q2finsave=q2fin
0059       zzzz=zparpro(kcol)+zpartar(kcol)
0060       nnn=ltarg3(it)+lproj3(ip)
0061       zz=1.+zoeinc*zzzz**2               !<-----
0062       q2fin=q2fin*zz
0063 
0064 c      print *,kcol,zzzz,zz,q2fin
0065 
0066       ajeth(idpomr+1)=ajeth(idpomr+1)+1.
0067 c      write(*,*)ajeth
0068       idfpomr=idfpr(ncolp,kcol)
0069       if(ish.ge.3)write(ifch,*)'Start psahot (icp,ict):',ip,icp,it,ict
0070      *,ncolp,kcol,iptl,idpomr,idfpomr,bpomr
0071 
0072 
0073       if(idfpomr.eq.0)stop'idfpomr??????'
0074       if(ish.ge.3)then
0075         write(ifch,20)iptl
0076      *  ,sqrt(pptl(1,iptl)**2+pptl(2,iptl)**2),pptl(3,iptl)
0077      *  ,pptl(4,iptl),pptl(5,iptl)
0078      *  ,istptl(iptl),ityptl(iptl)
0079 20      format(1x,i4,3x,4(e11.5,1x),i2,1x,i3)
0080       endif
0081       istptl(iptl)=31
0082       
0083 csp  initialise to 0
0084 1        do i=1,nj
0085       ncj(1,i)=0
0086       ncj(2,i)=0
0087       enddo  
0088 csp end initialisation
0089 
0090       nj=0
0091       nptl=nptl1
0092       if(iremn.ge.2)then
0093         call iddeco(icp,jcp)
0094         call iddeco(ict,jct)
0095       endif
0096 
0097       wp0=dsqrt(xpr(ncolp,kcol))*dexp(ypr(ncolp,kcol))*dble(engy)     !???? new
0098       wm0=dsqrt(xpr(ncolp,kcol))*dexp(-ypr(ncolp,kcol))*dble(engy)    !double
0099 
0100 
0101       amqt(1)=sqrt(sngl(xxp1pr(ncolp,kcol)**2+xyp1pr(ncolp,kcol)**2))
0102       amqt(2)=sqrt(sngl(xxp2pr(ncolp,kcol)**2+xyp2pr(ncolp,kcol)**2))
0103       amqt(3)=sqrt(sngl(xxm2pr(ncolp,kcol)**2+xym2pr(ncolp,kcol)**2))
0104       amqt(4)=sqrt(sngl(xxm1pr(ncolp,kcol)**2+xym1pr(ncolp,kcol)**2))
0105       amqpt=amqt(1)+amqt(2)+amqt(3)+amqt(4)
0106 
0107       s2min=4.*q2min
0108       if(sngl(wp0*wm0).le.(sqrt(s2min)+amqpt)**2)then
0109         if(ish.ge.1)then
0110           call utmsg('psahot: insufficient pomeron mass&')
0111           write (ifch,*)'mass:',dsqrt(wp0*wm0),amqpt+sqrt(s2min)
0112           call utmsgf
0113         endif
0114         iret=1
0115         goto 16
0116       endif
0117 
0118       ih=iproj(kcol)
0119       jh=itarg(kcol)
0120 c      xprh=xpp(ih)
0121 c      xmrh=xmt(jh)
0122       rp=r2had(iclpro)+r2had(icltar)+slopom*log(engy**2)
0123       z=exp(-bpomr**2/(4.*.0389*rp)) !coef for rejection
0124 
0125       if(z.eq.0)then
0126        write(ifch,*)'psahot : z,ih,jh ! -> ',z,ih,jh
0127        call gakli2(ih,ih)
0128        call gakli2(jh,jh)
0129        call gakli2(iptl,iptl)
0130        stop
0131       endif
0132 
0133       do l=1,4
0134         bx(l)=xorptl(l,iptl)
0135       enddo
0136       bx(5)=tivptl(1,iptl)
0137       bx(6)=tivptl(2,iptl)
0138       ity=ityptl(iptl)
0139 
0140       if(idpomr.eq.0)then        !gg-pomeron
0141         iqq=0   !type of the hard interaction: 0 - gg, 1 - qg, 2 - gq, 3 - qq
0142         pxh=0.d0  !p_x for sh pomeron
0143         pyh=0.d0  !p_y for sh pomeron
0144       elseif(idpomr.eq.1)then    !qg-pomeron
0145         iqq=1
0146         pxh=xxp1pr(ncolp,kcol)
0147         pyh=xyp1pr(ncolp,kcol)
0148         amqpt=amqpt-amqt(1)
0149       elseif(idpomr.eq.2)then    !gq-pomeron
0150         iqq=2
0151         pxh=xxm2pr(ncolp,kcol)
0152         pyh=xym2pr(ncolp,kcol)
0153         amqpt=amqpt-amqt(3)
0154       elseif(idpomr.eq.3)then    !qq-pomeron
0155         iqq=3
0156         pxh=xxp1pr(ncolp,kcol)+xxm2pr(ncolp,kcol)
0157         pyh=xyp1pr(ncolp,kcol)+xym2pr(ncolp,kcol)
0158         amqpt=amqpt-amqt(1)-amqt(3)
0159       else
0160         stop'unknown pomeron'
0161       endif
0162       pth=pxh**2+pyh**2
0163 
0164       nj0=nj
0165       if(ish.ge.6)then
0166         write(ifch,*)'iptl,nptl,wp0,wm0,z,iqq,bx:'
0167         write(ifch,*) iptl,nptl,wp0,wm0,z,iqq,bx
0168       endif
0169 
0170       s=wp0*wm0                !lc+*lc- for the semihard interaction
0171       smin=dble(s2min)+pth          !mass cutoff for the hard pomeron
0172       xmin=smin/s
0173       smax=(dsqrt(s)-dble(amqpt))**2  !max mass for the hard pomeron
0174       xmax=smax/s
0175       !wwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwww
0176       if(iLHC.eq.1.and.iqq.ne.3)then     !not for val-val
0177 c xmin is the minimum value of zpm which is the fraction of the Pomeron
0178 c energy going into the hard part. If zpm=1 nothing left is the soft
0179 c preevolution part (momentum of soft string ends)
0180 c take into account z (impact parameter) dependence to avoid hard Pom
0181 c at large b (lot of rejection and not realistic)
0182 c        print *,xmin,zopinc*fegypp*exp(-bpomr**2/b2xscr),xmax
0183 c     *,smin,max(xmin,min(0.99*xmax,zzzz*zzsoft))*ss
0184       xmin=max(xmin,min((sqrt(s)-dble(amqpt+ammsdd))**2/s
0185      &  ,dble(max(nnn*fegypp*exp(-bpomr**2/b2xscr),zzzz)*zopinc*z)))
0186 c      xmin=max(xmin,1d0-exp(-dble(zzsoft*zzzz**2)))  !???????????????
0187       smin=xmin*s         
0188       endif
0189       !wwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwww
0190        if(smax.le.smin)then
0191         write (ifmt,*)'smax,smin',smax,smin
0192         iret=1
0193         goto 16
0194       endif
0195       xpmax0=psutz(s,smin,dble(amqpt**2))   !max x+ for the hard pomeron
0196       xpmin0=xmin/xpmax0              !min x+ for the hard pomeron
0197       xp1=wp0/dble(engy)              !lc+ share for the semihard interaction
0198       xp2=wm0/dble(engy)              !lc- share for the semihard interaction
0199 
0200 
0201 c-----------------------------------------------------------------
0202 c  determine LC momenta wpi,wmi for hard Pomeron
0203 c-----------------------------------------------------------------
0204 
0205       if(ish.ge.4)write(ifch,*)
0206      & 'determine LC momenta wpi,wmi for hard Pomeron'
0207 
0208       if(iremn.ge.2)then
0209         if(iclpro.eq.2)then
0210           if(iabs(idptl(ip)).eq.1120)then !proj=proton
0211             nquu1=jcp(1,1)+jcp(1,2)
0212             nqud1=jcp(2,1)+jcp(2,2)
0213           elseif(iabs(idptl(ip)).eq.1220)then !proj=neutron
0214             nquu1=jcp(2,1)+jcp(2,2)
0215             nqud1=jcp(1,1)+jcp(1,2)
0216           else    !to avoid flavor problem with exotic projectile (but should not happen (only gg)
0217             nquu1=0
0218             nqud1=0
0219           endif
0220         elseif(iclpro.eq.1)then
0221           if(iabs(idptl(ip)).eq.120)then
0222             nquu1=jcp(1,1)+jcp(1,2)
0223             nqud1=jcp(2,1)+jcp(2,2)
0224           else    !to avoid flavor problem with exotic projectile (but should not happen (only gg)
0225             nquu1=0
0226             nqud1=0
0227           endif
0228         elseif(iclpro.eq.3)then
0229           if(iabs(idptl(ip)).eq.130)then !proj=Kch
0230             nquu1=jcp(1,1)+jcp(1,2)
0231             nqud1=jcp(3,1)+jcp(3,2)
0232           elseif(iabs(idptl(ip)).eq.230)then  !proj=K0
0233             nquu1=jcp(2,1)+jcp(2,2)
0234             nqud1=jcp(3,1)+jcp(3,2)
0235           else    !to avoid flavor problem with exotic projectile (but should not happen (only gg)
0236             nquu1=0
0237             nqud1=0
0238           endif
0239         else                    !charm
0240           if(iabs(idptl(ip)).eq.140)then
0241             nquu1=jcp(1,1)+jcp(1,2)
0242             nqud1=jcp(4,1)+jcp(4,2)
0243           elseif(iabs(idptl(ip)).eq.240)then
0244             nquu1=jcp(2,1)+jcp(2,2)
0245             nqud1=jcp(4,1)+jcp(4,2)
0246           elseif(iabs(idptl(ip)).eq.340)then
0247             nquu1=jcp(3,1)+jcp(3,2)
0248             nqud1=jcp(4,1)+jcp(4,2)
0249           else
0250             nquu1=jcp(4,1)
0251             nqud1=jcp(4,2)
0252           endif
0253         endif
0254         if(iabs(idptl(maproj+it)).eq.1220)then !targ=neutron
0255           nquu2=jct(2,1)+jct(2,2)
0256           nqud2=jct(1,1)+jct(1,2)
0257         else
0258           nquu2=jct(1,1)+jct(1,2)
0259           nqud2=jct(2,1)+jct(2,2)
0260         endif
0261       endif
0262 
0263       iq1=0
0264       iq2=0
0265       wwgg=0.
0266       wwgq=0.
0267       wwqg=0.
0268       wwqq=0.
0269       wwdd=0.
0270 
0271       !------------------------------------------
0272       if(iqq.eq.3)then     !     val-val
0273       !------------------------------------------
0274 
0275        if(ish.ge.4)write(ifch,*)'val-val'
0276         xmin1=xmin**dble(delh+.4)
0277         xmax1=xmax**dble(delh+.4)
0278         zp0=dsqrt(xmin)
0279         if(zp0.ge.1.d0)call utstop('zp0 in sem&')
0280       !........ kinematical bounds
0281         tmax0=dlog((1.d0+dsqrt(1.d0-zp0))/(1.d0-dsqrt(1.d0-zp0)))
0282         tmin0=dlog((1.d0+dsqrt(1.d0-xpmax0))
0283      *                   /(1.d0-dsqrt(1.d0-xpmax0)))
0284         if(iclpro.ne.4)then
0285           call psjti0(sngl(smax-pth),sqq,sqqb,1,1)
0286           call psjti0(sngl(smax-pth),sqqp,sqqpb,1,2)
0287           call psjti0(sngl(smax-pth),sqaq,sqaqb,-1,1)
0288         else
0289           call psjti0(sngl(smax-pth),sqqp,sqqpb,4,2)
0290           sqq=0.
0291           sqaq=0.
0292         endif
0293         if(iremn.ge.2)then
0294           if(nquu1.gt.nqud1.or.iclpro.ne.2)then
0295             uv1=psdfh4(sngl(zp0*xp1),q2min,0.,iclpro,1)
0296             dv1=psdfh4(sngl(zp0*xp1),q2min,0.,iclpro,2)
0297           else         !if nquu<nqud => no u or no d
0298             uv1=psdfh4(sngl(zp0*xp1),q2min,0.,1,1)
0299             dv1=uv1
0300           endif
0301           if(nquu1.eq.0)uv1=0d0
0302           if(nqud1.eq.0)dv1=0d0
0303           if(nquu2.gt.nqud2)then
0304             uv2=psdfh4(sngl(zp0*xp2),q2min,0.,icltar,1)
0305             dv2=psdfh4(sngl(zp0*xp2),q2min,0.,icltar,2)
0306           else                  !if nquu<nqud => no u or no d
0307             uv2=psdfh4(sngl(zp0*xp2),q2min,0.,1,1)
0308             dv2=uv2
0309           endif
0310           if(nquu2.eq.0)uv2=0d0
0311           if(nqud2.eq.0)dv2=0d0
0312         else
0313           uv1=psdfh4(sngl(zp0*xp1),q2min,0.,iclpro,1)
0314           dv1=psdfh4(sngl(zp0*xp1),q2min,0.,iclpro,2)
0315           uv2=psdfh4(sngl(zp0*xp2),q2min,0.,icltar,1)
0316           dv2=psdfh4(sngl(zp0*xp2),q2min,0.,icltar,2)
0317         endif
0318         wwuu=uv1*uv2*sqq
0319         if(iclpro.eq.2)then
0320           wwdd=dv1*dv2*sqq
0321         elseif(iclpro.eq.1)then
0322           wwdd=dv1*dv2*sqaq
0323         elseif(iclpro.eq.3)then
0324           wwdd=dv1*dv2*sqqp
0325         elseif(iclpro.eq.4)then
0326           wwuu=uv1*uv2*sqqp
0327           wwdd=0.
0328         endif
0329         wwud=uv1*dv2*sqqp
0330         wwdu=dv1*uv2*sqqp
0331         wudt=wwuu+wwdd+wwud+wwdu
0332         gb0=dble(wudt)/xmax**dble(delh)/xmin**0.4*
0333      *    (1.d0-zp0*xp1)**dble(-1.-alpq-alplea(iclpro))*
0334      *    (1.d0-zp0*xp2)**dble(-1.-alpq-alplea(icltar))*(tmax0-tmin0)
0335      *    *(1.d0-zp0)**dble(.5+alpq)*(1.d0-zp0)**dble(alpq)  *5.d0
0336 
0337 3       zpm=(xmin1+dble(rangen())*(xmax1-xmin1))**dble(1./(delh+.4))   !zpm
0338         if(iclpro.ne.4)then
0339           call psjti0(sngl(zpm*s-pth),sqq,sqqb,1,1)
0340           call psjti0(sngl(zpm*s-pth),sqqp,sqqpb,1,2)
0341           call psjti0(sngl(zpm*s-pth),sqaq,sqaqb,-1,1)
0342         else
0343           call psjti0(sngl(zpm*s-pth),sqqp,sqqpb,4,2)
0344           sqq=0.
0345           sqaq=0.
0346         endif
0347         xpmax=psutz(s,zpm*s,dble(amqpt**2))  !max x+ for sh pomeron
0348         tmax=dlog((1.d0+dsqrt(1.d0-dsqrt(zpm)))
0349      *             /(1.d0-dsqrt(1.d0-dsqrt(zpm))))
0350         tmin=dlog((1.d0+dsqrt(1.d0-xpmax))/(1.d0-dsqrt(1.d0-xpmax)))
0351 
0352         t=(tmin+dble(rangen())*(tmax-tmin))
0353         xp=1.d0-((1.d0-dexp(-t))/(1.d0+dexp(-t)))**2  !x+_v
0354         xm=zpm/xp                             !x-_v
0355         if(xm.gt.xp.and.ish.ge.1)write(ifmt,*)'xm,xp',xm,xp
0356         gb=(1.d0-xm)**alpq*(1.d0-xp)**(.5+alpq)*(tmax-tmin)
0357         if(rangen().lt..5)then
0358           xp=xm
0359           xm=zpm/xp
0360         endif
0361 
0362         if(iremn.ge.2)then
0363           if(nquu1.gt.nqud1.or.iclpro.ne.2)then
0364             uv1=psdfh4(sngl(xp*xp1),q2min,0.,iclpro,1)
0365             dv1=psdfh4(sngl(xp*xp1),q2min,0.,iclpro,2)
0366           else         !if nquu<nqud => no u or no d
0367             uv1=psdfh4(sngl(xp*xp1),q2min,0.,1,1)
0368             dv1=uv1
0369           endif
0370           if(nquu1.eq.0)uv1=0d0
0371           if(nqud1.eq.0)dv1=0d0
0372           if(nquu2.gt.nqud2)then
0373             uv2=psdfh4(sngl(xm*xp2),q2min,0.,icltar,1)
0374             dv2=psdfh4(sngl(xm*xp2),q2min,0.,icltar,2)
0375           else                  !if nquu<nqud => no u or no d
0376             uv2=psdfh4(sngl(xm*xp2),q2min,0.,1,1)
0377             dv2=uv2
0378           endif
0379           if(nquu2.eq.0)uv2=0d0
0380           if(nqud2.eq.0)dv2=0d0
0381         else
0382           uv1=psdfh4(sngl(xp*xp1),q2min,0.,iclpro,1)
0383           dv1=psdfh4(sngl(xp*xp1),q2min,0.,iclpro,2)
0384           uv2=psdfh4(sngl(xm*xp2),q2min,0.,icltar,1)
0385           dv2=psdfh4(sngl(xm*xp2),q2min,0.,icltar,2)
0386         endif
0387 
0388         wwuu=uv1*uv2*sqq
0389         if(iclpro.eq.2)then
0390           wwdd=dv1*dv2*sqq
0391         elseif(iclpro.eq.1)then
0392           wwdd=dv1*dv2*sqaq
0393         elseif(iclpro.eq.3)then
0394           wwdd=dv1*dv2*sqqp
0395         elseif(iclpro.eq.4)then
0396           wwuu=uv1*uv2*sqqp
0397           wwdd=0.
0398         endif
0399         wwud=uv1*dv2*sqqp
0400         wwdu=dv1*uv2*sqqp
0401         wudt=wwuu+wwdd+wwud+wwdu
0402         if(wudt.lt.1d-16)then
0403           if(ish.ge.1)write(ifmt,*)'No more valence quark for psahot !'
0404           write(ifch,*)'No more valence quark for psahot !'
0405      &                     ,ip,it,nquu1,nqud1,nquu2,nqud2
0406           iret=1
0407           goto 16
0408         endif
0409 
0410         gb=gb*dble(wudt)/zpm**dble(delh+0.4)
0411      *    *(1.d0-xp*xp1)**dble(-1.-alpq-alplea(iclpro))
0412      *    *(1.d0-xm*xp2)**dble(-1.-alpq-alplea(icltar))/gb0
0413 c          if(ish.ge.4)then
0414           if(gb.gt.1.d0.and.ish.ge.1)write (ifch,*)
0415      *      'gb-qq,iclpro,zpm,xp,tmax,tmin,xpmax',
0416      *      gb,iclpro,zpm,xp,tmax,tmin,xpmax
0417 c          endif
0418         if(dble(rangen()).gt.gb)goto 3
0419 
0420         aks=rangen()*wudt
0421         if(aks.le.wwuu)then
0422           if(iclpro.le.2)then
0423             iq1=1
0424           elseif(iclpro.eq.3)then
0425             if(iabs(idptl(ip)).eq.130)then !proj=Kch
0426               iq1=1
0427             else !proj=K0
0428               iq1=2
0429             endif
0430           else   !charm
0431             if(iabs(idptl(ip)).eq.140)then
0432               iq1=1
0433             elseif(iabs(idptl(ip)).eq.240)then
0434               iq1=2
0435             elseif(iabs(idptl(ip)).eq.340)then
0436               iq1=3
0437             else
0438               iq1=4
0439             endif
0440           endif
0441           iq2=1
0442         elseif(aks.le.wwuu+wwdd)then
0443           if(iclpro.eq.2)then
0444             iq1=2
0445           elseif(iclpro.eq.1)then
0446             iq1=-2
0447           elseif(iclpro.eq.3)then
0448             iq1=-3
0449           else
0450             iq1=-4
0451           endif
0452           iq2=2
0453         elseif(aks.le.wwuu+wwdd+wwud)then
0454           if(iclpro.le.2)then
0455             iq1=1
0456           elseif(iclpro.eq.3)then
0457             if(iabs(idptl(ip)).eq.130)then !proj=Kch
0458               iq1=1
0459             else !proj=K0
0460               iq1=2
0461             endif
0462           else   !charm
0463             if(iabs(idptl(ip)).eq.140)then
0464               iq1=1
0465             elseif(iabs(idptl(ip)).eq.240)then
0466               iq1=2
0467             elseif(iabs(idptl(ip)).eq.340)then
0468               iq1=3
0469             else
0470               iq1=4
0471             endif
0472           endif
0473           iq2=2
0474         else
0475           if(iclpro.eq.2)then
0476             iq1=2
0477           elseif(iclpro.eq.1)then
0478             iq1=-2
0479           elseif(iclpro.eq.3)then
0480             iq1=-3
0481           else
0482             iq1=-4
0483           endif
0484           iq2=1
0485         endif
0486 
0487         wpi=xp*wp0       !lc+ for the semihard interaction
0488         wmi=xm*wm0       !lc- for the semihard interaction
0489         wp1=(wp0-wpi)
0490         wm1=(wm0-wmi)
0491         wp1=wp1*psutz(wp1*wm1,dble(amqt(2)**2),dble(amqt(4)**2))
0492         wm1=wm1-amqt(2)**2/wp1
0493 
0494       !-------------------------------------
0495       else  ! sea-sea  val-sea  sea-val
0496       !-------------------------------------
0497 
0498        if(ish.ge.4)write(ifch,*)'sea-sea  val-sea  sea-val'
0499         xmin1=xmin**(delh-dels)
0500         xmax1=xmax**(delh-dels)
0501 
0502         if(iqq.eq.0)then    !rejection function normalization
0503           gb0=dlog(xpmax0/xpmin0)*(1.d0-dsqrt(xmin))**(2.*betpom) !y_soft =
0504         else
0505           tmax0=acos(dsqrt(xpmin0))  !kinematical limits for t=cos(x+-)**2
0506           tmin0=acos(dsqrt(xpmax0))
0507           if(iqq.eq.1)then
0508             uv1=psdfh4(sngl(xpmin0*xp1),q2min,0.,iclpro,1)
0509      *        *sngl(1.d0-xpmin0*xp1)**(-1.-alpq-alplea(iclpro))
0510             dv1=psdfh4(sngl(xpmin0*xp1),q2min,0.,iclpro,2)
0511      *        *sngl(1.d0-xpmin0*xp1)**(-1.-alpq-alplea(iclpro))
0512           else
0513             uv1=psdfh4(sngl(xpmin0*xp2),q2min,0.,icltar,1)
0514      *        *sngl(1.d0-xpmin0*xp2)**(-1.-alpq-alplea(icltar))
0515             dv1=psdfh4(sngl(xpmin0*xp2),q2min,0.,icltar,2)
0516      *        *sngl(1.d0-xpmin0*xp2)**(-1.-alpq-alplea(icltar))
0517           endif
0518           gb0=(1.d0-xmin/xpmax0)**betpom*dble(uv1+dv1)
0519      *      *xpmin0**(-0.5+dels)
0520      *      *(1.d0-xpmin0)**(0.5+alpq)*(tmax0-tmin0)
0521           if(ish.ge.6)write (ifch,*)
0522      *      'gb0,tmax0,tmin0,xpmax0,xpmin0,xmin,xp1,xp2',
0523      *      gb0,tmax0,tmin0,xpmax0,xpmin0,xmin,xp1,xp2
0524         endif
0525         if(iclpro.ne.4.or.iqq.ne.1)then
0526           call psjti0(sngl(smax-pth),sj,sjb,iqq,0) !inclusive (sj) and born
0527         else
0528           call psjti0(sngl(smax-pth),sj,sjb,4,0)
0529         endif
0530         gb0=gb0*dble(sj)/xmax**delh      *1.5d0    !rejection function norm.
0531         if(ish.ge.6)write (ifch,*)'gb0,sj,z',gb0,sj,z
0532 
0533         if(gb0.le.0.)then
0534           write (ifmt,*)'gb0<0, smax,pth',smax,pth
0535           iret=1
0536           goto 16
0537         endif
0538 
0539         ! sharing of light cone momenta between soft preevolution and
0540         ! hard interaction: ( first energy-momentum is shared according to
0541         ! f_hard(yj)~zpm**(delh-dels-1) and then rejected as
0542         ! w_rej ~sigma_hard_tot(yj) / exp(delh*yj)
0543 
0544  4        continue
0545         zpm=(xmin1+dble(rangen())*(xmax1-xmin1))**dble(1./(delh-dels)) !zpm
0546         if(iclpro.ne.4.or.iqq.ne.1)then
0547           call psjti0(sngl(zpm*s-pth),sgg,sggb,0,0)!inclusive (sgg) and born
0548           call psjti0(sngl(zpm*s-pth),sgq,sgqb,0,1)
0549           call psjti0(sngl(zpm*s-pth),sqq,sqqb,1,1)
0550           call psjti0(sngl(zpm*s-pth),sqaq,sqaqb,-1,1)
0551           call psjti0(sngl(zpm*s-pth),sqqp,sqqpb,1,2)
0552           sqq=(sqq+sqaq+2.*(naflav-1)*sqqp)/naflav/2.
0553         else
0554           call psjti0(sngl(zpm*s-pth),sgq,sgqb,4,0)
0555           call psjti0(sngl(zpm*s-pth),sqq,sqqb,4,1)
0556         endif
0557         xpmax=psutz(s,zpm*s,dble(amqpt**2))  !max x+ for sh pomeron
0558         xpmin=zpm/xpmax
0559         if(ish.ge.8)write (ifch,*)'zpm,xpmax,xpmin',zpm,xpmax,xpmin
0560 
0561         if(iqq.eq.0)then
0562           xp=xpmin*(xpmax/xpmin)**rangen()  !lc+ share for the hard interaction
0563           xm=zpm/xp           !lc- share for the hard interaction
0564         else
0565           tmax=acos(dsqrt(xpmin))  !kinematical limits for t=cos(x+-)**2
0566           tmin=acos(dsqrt(xpmax))
0567           t=tmin+dble(rangen())*(tmax-tmin)
0568           xp=cos(t)**2
0569           xm=zpm/xp
0570         endif
0571 
0572         if(ish.ge.8)write(ifch,*)'zpm,xp,xm,xpmax,xpmin:',
0573      *      zpm,xp,xm,xpmax,xpmin
0574 
0575 
0576         if(iqq.eq.0)then  ! --------- sea-sea -----------
0577 
0578           dwwgg1=0.
0579           dwwgq1=0.
0580           dwwqg1=0.
0581           dwwqq1=0.
0582 
0583           dwwgg2=0.
0584           dwwgq2=0.
0585           dwwqg2=0.
0586           dwwqq2=0.
0587 
0588           wwgg=sngl((1.d0-xp)*(1.d0-xm))**betpom
0589           wwgq=sngl(1.d0-xp)**betpom*EsoftQZero(sngl(xm))
0590           wwqg=sngl(1.d0-xm)**betpom*EsoftQZero(sngl(xp))
0591           wwqq=EsoftQZero(sngl(xp))*EsoftQZero(sngl(xm))
0592 
0593           if(idfpomr.eq.1)then
0594             rh=r2had(iclpro)+r2had(icltar)-slopom*sngl(dlog(zpm))
0595             wwgg=(wwgg*z**(rp/rh)/rh+dwwgg1+dwwgg2)
0596      *        *(r2had(iclpro)+r2had(icltar))/z
0597             wwgq=(wwgq*z**(rp/rh)/rh+dwwgq1+dwwgq2)
0598      *        *(r2had(iclpro)+r2had(icltar))/z
0599             wwqg=(wwqg*z**(rp/rh)/rh+dwwqg1+dwwqg2)
0600      *        *(r2had(iclpro)+r2had(icltar))/z
0601             wwqq=(wwqq*z**(rp/rh)/rh+dwwqq1+dwwqq2)
0602      *        *(r2had(iclpro)+r2had(icltar))/z
0603           elseif(idfpomr.eq.2)then
0604             rh=r2had(iclpro)+alam3p/2.-slopom*sngl(dlog(zpm))
0605             wwgg=wwgg/rh*(r2had(iclpro)+alam3p/2.)*z**(rp/rh-1.)
0606             wwqg=wwqg/rh*(r2had(iclpro)+alam3p/2.)*z**(rp/rh-1.)
0607             wwgq=wwgq/rh*(r2had(iclpro)+alam3p/2.)*z**(rp/rh-1.)
0608             wwqq=wwqq/rh*(r2had(iclpro)+alam3p/2.)*z**(rp/rh-1.)
0609           elseif(idfpomr.eq.3)then
0610             rh=r2had(icltar)+alam3p/2.-slopom*sngl(dlog(zpm))
0611             wwgg=wwgg/rh*(r2had(icltar)+alam3p/2.)*z**(rp/rh-1.)
0612             wwqg=wwqg/rh*(r2had(icltar)+alam3p/2.)*z**(rp/rh-1.)
0613             wwgq=wwgq/rh*(r2had(icltar)+alam3p/2.)*z**(rp/rh-1.)
0614             wwqq=wwqq/rh*(r2had(icltar)+alam3p/2.)*z**(rp/rh-1.)
0615           elseif(idfpomr.eq.4)then
0616             rh=alam3p-slopom*sngl(dlog(zpm))
0617             wwgg=wwgg/rh*alam3p*z**(rp/rh-1.)
0618             wwqg=wwqg/rh*alam3p*z**(rp/rh-1.)
0619             wwgq=wwgq/rh*alam3p*z**(rp/rh-1.)
0620             wwqq=wwqq/rh*alam3p*z**(rp/rh-1.)
0621           else
0622             stop'psahot-idfpomr????'
0623           endif
0624 
0625           wwgg=wwgg*sgg*(1.-glusea)**2
0626           wwgq=wwgq*sgq*(1.-glusea)*glusea
0627           wwqg=wwqg*sgq*(1.-glusea)*glusea
0628           wwqq=wwqq*sqq*glusea**2
0629           gbyj=dlog(xpmax/xpmin)*dble(wwgg+wwgq+wwqg+wwqq)
0630           wpi=wp0*xp               !lc+ for the hard interaction
0631           wmi=wm0*xm               !lc+-for the hard interaction
0632           gbyj=gbyj/zpm**dble(delh)/gb0     !rejection fu
0633           if(gbyj.ge.1.d0.and.ish.ge.2)write(ifmt,*)'gbyj',gbyj
0634           if(dble(rangen()).gt.gbyj)goto 4                     !rejection
0635           wp1=wp0-wpi
0636           wm1=wm0-wmi
0637           call pslcsh(wp1,wm1,wp2,wm2,amqt,dble(amqpt))
0638 
0639         else  ! --------- val-sea  sea-val -----------
0640 
0641           dwwg=0.
0642           dwwq=0.
0643 
0644           wwgq=sngl(1.d0-xm)**betpom
0645           wwqq=EsoftQZero(sngl(xm))
0646           if(idfpomr.eq.1)then
0647             rh=r2had(iclpro)+r2had(icltar)-slopom*sngl(dlog(xm))
0648             wwgq=(wwgq*z**(rp/rh)/rh+dwwg)
0649      *        *(r2had(iclpro)+r2had(icltar))/z
0650             wwqq=(wwqq*z**(rp/rh)/rh+dwwq)
0651      *        *(r2had(iclpro)+r2had(icltar))/z
0652           else              !tp071031 not used anymore
0653             if(iqq.eq.1)then
0654               rh=r2had(iclpro)+alam3p/2.-slopom*sngl(dlog(xm))
0655               wwgq=wwgq/rh*(r2had(iclpro)+alam3p/2.)*z**(rp/rh-1.)
0656               wwqq=wwqq/rh*(r2had(iclpro)+alam3p/2.)*z**(rp/rh-1.)
0657             else
0658               rh=r2had(icltar)+alam3p/2.-slopom*sngl(dlog(xm))
0659               wwgq=wwgq/rh*(r2had(icltar)+alam3p/2.)*z**(rp/rh-1.)
0660               wwqq=wwqq/rh*(r2had(icltar)+alam3p/2.)*z**(rp/rh-1.)
0661             endif
0662           endif
0663 
0664           wwgq=wwgq*sgq*(1.-glusea)*sngl(xp)**(-0.5+dels)
0665      *      *sngl(1.d0-xp)**(0.5+alpq)
0666           wwqq=wwqq*sqq*glusea*sngl(xp)**(-0.5+dels)
0667      *      *sngl(1.d0-xp)**(0.5+alpq)
0668 
0669           if(iqq.eq.1)then         !valence quark-gluon hard interaction
0670             if(iremn.ge.2)then
0671               if(nquu1.gt.nqud1.or.iclpro.ne.2)then
0672                 uv1=psdfh4(sngl(xp*xp1),q2min,0.,iclpro,1)
0673                 dv1=psdfh4(sngl(xp*xp1),q2min,0.,iclpro,2)
0674               else       !if nquu<nqud => no u or no d
0675                 uv1=psdfh4(sngl(xp*xp1),q2min,0.,1,1)
0676                 dv1=uv1
0677               endif
0678               if(nquu1.eq.0)uv1=0d0
0679               if(nqud1.eq.0)dv1=0d0
0680             else
0681               uv1=psdfh4(sngl(xp*xp1),q2min,0.,iclpro,1)
0682               dv1=psdfh4(sngl(xp*xp1),q2min,0.,iclpro,2)
0683             endif
0684             if(uv1+dv1.lt.1d-16)then
0685            if(ish.ge.1)write(ifmt,*)'No more valence quark for psahot !'
0686          write(ifch,*)'No more valence quark in projectile for psahot !'
0687      &                     ,ip,nquu1,nqud1
0688               iret=1
0689               goto 16
0690             endif
0691             wpi=wp0*xp             !lc+ for the hard interaction
0692             wmi=wm0*xm             !lc+-for the hard interaction
0693             aks=rangen()
0694             if(aks.le.uv1/(uv1+dv1))then
0695               if(iclpro.le.2)then
0696                 iq1=1
0697               elseif(iclpro.eq.3)then
0698                 if(iabs(idptl(ip)).eq.130)then !proj=Kch
0699                   iq1=1
0700                 else !proj=K0
0701                   iq1=2
0702                 endif
0703               else              !charm
0704                 if(iabs(idptl(ip)).eq.140)then
0705                   iq1=1
0706                 elseif(iabs(idptl(ip)).eq.240)then
0707                   iq1=2
0708                 elseif(iabs(idptl(ip)).eq.340)then
0709                   iq1=3
0710                 else
0711                   iq1=4
0712                 endif
0713               endif
0714             else
0715               if(iclpro.eq.2)then
0716                 iq1=2
0717               elseif(iclpro.eq.1)then
0718                 iq1=-2
0719               elseif(iclpro.eq.3)then
0720                 iq1=-3
0721               else
0722                 iq1=-4
0723               endif
0724             endif
0725             gbyj=dble((wwgq+wwqq)*(uv1+dv1))*(1.d0-xp*xp1)**
0726      *        (-1.-alpq-alplea(iclpro))
0727           else                !gluon-valence quark hard interaction
0728             xm=xp
0729             xp=zpm/xm
0730             if(iremn.ge.2)then
0731               if(nquu2.gt.nqud2)then
0732                 uv1=psdfh4(sngl(xm*xp2),q2min,0.,icltar,1)
0733                 dv1=psdfh4(sngl(xm*xp2),q2min,0.,icltar,2)
0734               else              !if nquu<nqud => no u or no d
0735                 uv1=psdfh4(sngl(xm*xp2),q2min,0.,1,1)
0736                 dv1=uv1
0737               endif
0738               if(nquu2.eq.0)uv1=0d0
0739               if(nqud2.eq.0)dv1=0d0
0740             else
0741               uv1=psdfh4(sngl(xm*xp2),q2min,0.,icltar,1)
0742               dv1=psdfh4(sngl(xm*xp2),q2min,0.,icltar,2)
0743             endif
0744             if(uv1+dv1.lt.1d-16)then
0745            if(ish.ge.1)write(ifmt,*)'No more valence quark for psahot !'
0746              write(ifch,*)'No more valence quark in target for psahot !'
0747      &                     ,it,nquu2,nqud2
0748               iret=1
0749               goto 16
0750             endif
0751             wpi=wp0*xp             !lc+ for the hard interaction
0752             wmi=wm0*xm             !lc+-for the hard interaction
0753             aks=rangen()
0754             if(aks.le.uv1/(uv1+dv1))then
0755               iq2=1
0756             else
0757               iq2=2
0758             endif
0759             gbyj=dble(wwgq+wwqq)*dble(uv1+dv1)*
0760      *        (1.d0-xm*xp2)**(-1.-alpq-alplea(icltar))
0761           endif
0762 
0763           gbyj=gbyj*(tmax-tmin)/zpm**delh /gb0 /2.1d0 !rejection
0764           if(ish.ge.6)write (ifch,*)
0765      *      'gbyj,zpm,xp,tmax,tmin,xpmax,xpmin',
0766      *      gbyj,zpm,xp,tmax,tmin,xpmax,xpmin
0767           if(dble(rangen()).gt.gbyj)goto 4
0768 
0769           wp1=wp0-wpi
0770           wm1=wm0-wmi
0771           if(ish.ge.8)write (ifch,*)'q_sea mass check',wp1*wm1,amqpt
0772           if(iqq.eq.1)then         !valence quark-gluon hard interaction
0773             amq1=amqt(3)**2
0774             s24=(amqt(2)+amqt(4))**2
0775           else
0776             amq1=amqt(1)**2
0777             s24=(amqt(2)+amqt(4))**2
0778             xp=xm
0779             xm=zpm/xp
0780           endif
0781           x1max=psutz(wp1*wm1,dble(amq1),dble(s24))
0782           x1min=dble(amq1)/x1max/wp1/wm1
0783           t1min=(1.d0/x1max-1.d0)
0784           t1max=(1.d0/x1min-1.d0)
0785 5           t1=t1min*(t1max/t1min)**dble(rangen())
0786           if(ish.ge.8)write (ifch,*)'t1,t1min,t1max',t1,t1min,t1max
0787           xq1=1.d0/(1.d0+t1)
0788           if(dble(rangen()).gt.(xq1*(1.d0-xq1))**(1.+(-alpqua)))goto 5
0789           if(iqq.eq.1)then         !valence quark-gluon hard interacti
0790             wm2=wm1*(1.d0-xq1)
0791             wm1=wm1*xq1
0792             wp1=wp1-dble(amq1)/wm1
0793             if(ish.ge.8)write (ifch,*)'q_sea+ mass check',
0794      *        wp1*wm2,s24
0795             wp1=wp1*psutz(wp1*wm2,dble(amqt(2)**2),dble(amqt(4)**2))
0796             wm2=wm2-dble(amqt(2))**2/wp1
0797           else
0798             wp2=wp1*(1.d0-xq1)
0799             wp1=wp1*xq1
0800             wm1=wm1-dble(amq1)/wp1
0801             if(ish.ge.8)write (ifch,*)'q_sea- mass check',
0802      *        wp2*wm1,s24
0803             wm1=wm1*psutz(wm1*wp2,dble(amqt(4)**2),dble(amqt(2)**2))
0804             wp2=wp2-amqt(4)**2/wm1
0805           endif
0806 
0807         endif  ! -------------------
0808 
0809       !-------------------------------
0810       endif
0811       !-------------------------------
0812 
0813 c-------------------------------------------------------------------------
0814 c  flavor and momenta of end partons of the hard Pomeron
0815 c-------------------------------------------------------------------------
0816 
0817        if(ish.ge.4)write(ifch,*)
0818      &  'flavor and momenta of end partons of the hard Pomeron'
0819 
0820   6   continue
0821       wpi1=wpi
0822       wmi1=wmi
0823       nj=nj0          !initialization for the number of jets
0824       if(ish.ge.8)write (ifch,*)'5-ww,smin',wpi*wmi,smin
0825 
0826       rrr=rangen()
0827       jqq=0
0828 
0829       call iddeco(icp,jcp)     !save valence quarks in jcp and jct
0830       call iddeco(ict,jct)
0831       if(iremn.ge.2)then
0832         do j=1,2
0833           do n=1,nrflav
0834             jcpr(n,j)=jcpref(n,j,ip)   !remnant sea quarks in jcpr and jctr
0835             jctr(n,j)=jctref(n,j,it)
0836           enddo
0837           do n=nrflav+1,nflav
0838             jcpr(n,j)=0
0839             jctr(n,j)=0
0840           enddo
0841         enddo
0842       endif
0843       do i=1,2
0844         icp1(i)=0
0845         icp2(i)=0
0846         icm1(i)=0
0847         icm2(i)=0
0848       enddo
0849 
0850       iret1=0
0851       iret2=0
0852 
0853 c if one or 2 of the 2 initial pomeron ends is sea, choose whether the first
0854 c emitted parton is a gluon or a quark to define the hard interaction
0855 c  jqq=0 : gg interaction after soft emission
0856 c  jqq=1 : qg interaction after soft emission
0857 c  jqq=2 : gq interaction after soft emission
0858 c  jqq=3 : qq interaction after soft emission
0859 c If the initial parton is a quark, the flavor has to be defined and the
0860 c antiflavor remains as string end.
0861 
0862       if(iqq.eq.1.or.iqq.eq.2)then
0863         if(rrr.lt.wwqq/(wwgq+wwqq))jqq=1
0864       elseif(iqq.eq.0)then
0865         if(rrr.lt.wwqg/(wwgg+wwqg+wwgq+wwqq))then
0866           jqq=1
0867         elseif(rrr.lt.(wwqg+wwgq)/(wwgg+wwqg+wwgq+wwqq))then
0868           jqq=2
0869         elseif(rrr.lt.(wwqg+wwgq+wwqq)/(wwgg+wwqg+wwgq+wwqq))then
0870           jqq=3
0871         endif
0872       endif
0873 
0874       if(iLHC.eq.1)then
0875 
0876       if((iqq-1)*(iqq-3).eq.0)then    !valence quark on projectile side
0877 c valence quark used to start space like cascade
0878         iqc(1)=iq1                                !proj=particle
0879         if(iabs(idptl(ip)).eq.1220)iqc(1)=3-iq1      !proj=neutron
0880         if(idptl(ip).lt.0)iqc(1)=-iqc(1)               !proj=antiparticle
0881         ifl1=iabs(iqc(1))
0882 
0883 c check string ends for mesons : 
0884 c switch q<->aq if needed (done randomly in ProSeTyp)
0885         if(iabs(idptl(ip)).lt.1000)then
0886           if(iqc(1).gt.0.and.idp1pr(ncolp,kcol).ne.2)then
0887             idp2pr(ncolp,kcol)=idp1pr(ncolp,kcol)
0888             idp1pr(ncolp,kcol)=2
0889           elseif(iqc(1).lt.0.and.idp2pr(ncolp,kcol).ne.2)then
0890             idp1pr(ncolp,kcol)=idp2pr(ncolp,kcol)
0891             idp2pr(ncolp,kcol)=2
0892           endif
0893         endif
0894 
0895 c store flavor of used valence quark in icp1(1) (quark) or icp2(2) (antiquark)
0896         idsp=100
0897         if(iqc(1).gt.0)then
0898           icp1(1)=ifl1
0899           if(idp1pr(ncolp,kcol).ne.2)
0900      &    call utstop("psahot: Problem with SE (1)!&")
0901         else
0902           icp2(2)=ifl1
0903           if(idp2pr(ncolp,kcol).ne.2)
0904      &    call utstop("psahot: Problem with SE (2)!&")
0905         endif
0906       else
0907         idsp=idsppr(ncolp,kcol)
0908       endif
0909 
0910       if((iqq-2)*(iqq-3).eq.0)then    !valence quark on target side
0911 c valence quark used to start space like cascade
0912         iqc(2)=iq2             !tar=particle (can not be pion or kaon !)
0913         if(iabs(idptl(maproj+it)).eq.1220)iqc(2)=3-iq2     !targ=neutron
0914         if(idptl(maproj+it).lt.0)iqc(2)=-iqc(2)  !targ=antiparticle
0915         ifl2=iabs(iqc(2))
0916 
0917 c store flavor of used valence quark in icm1(1) (quark) or icm2(2) (antiquark)
0918         idst=100
0919         if(iqc(2).gt.0)then
0920           icm1(1)=ifl2
0921           if(idm1pr(ncolp,kcol).ne.2)
0922      &    call utstop("psahot: Problem with SE (3)!&")
0923         else
0924           icm2(2)=ifl2
0925           if(idm2pr(ncolp,kcol).ne.2)
0926      &    call utstop("psahot: Problem with SE (4)!&")
0927         endif
0928       else
0929         idst=idstpr(ncolp,kcol)
0930       endif
0931 
0932 c determine soft string end flavors
0933 
0934       iret=0
0935 
0936       call fstrfl(jcpr,jctr,jcp,jct,icp1,icp2,icm1,icm2
0937      *           ,idp1pr(ncolp,kcol),idp2pr(ncolp,kcol)
0938      *           ,idm1pr(ncolp,kcol),idm2pr(ncolp,kcol)
0939      *                                  ,idsp,idst,iret)
0940 
0941       if(iret.ne.0)goto 16
0942 
0943 c put partons in stack for cascades
0944 
0945       if((iqq-1)*(iqq-3).eq.0)then    !valence quark on projectile side
0946 
0947 
0948 c define the second string end from sea
0949         if(iqc(1).gt.0)then
0950           ifl=idtra(icp2,1,0,2)
0951         elseif(iqc(1).lt.0)then
0952           ifl=idtra(icp1,1,0,2)
0953         else
0954           call utstop('No quark for hard Pomeron+ in psahot!&')
0955         endif
0956 
0957         if(ish.ge.5)write(ifch,*)'flavor vq+,sqb+',iqc(1)
0958      &                                            ,ifl
0959 
0960         nj=nj+1
0961 c put the sea quark (or diquark) into the parton list as string end
0962         if(abs(ifl).eq.3)then
0963           iqj(nj)=ifl*4/3
0964         elseif(abs(ifl).ge.4)then
0965           iqj(nj)=ifl*10
0966         else
0967           iqj(nj)=ifl
0968         endif
0969         ioj(nj)=7
0970 
0971         eqj(1,nj)=.5*sngl(wp1+dble(amqt(2))**2/wp1)
0972         eqj(2,nj)=wp1-eqj(1,nj)
0973         eqj(3,nj)=xxp2pr(ncolp,kcol)
0974         eqj(4,nj)=xyp2pr(ncolp,kcol)
0975         if(ish.ge.8)write (ifch,*)'q_v+ mass check',(eqj(1,nj)-
0976      *    eqj(2,nj))*(eqj(1,nj)+eqj(2,nj))-eqj(3,nj)**2-eqj(4,nj)**2
0977         eqj(1,nj)=sqrt(eqj(2,nj)**2+eqj(3,nj)**2+eqj(4,nj)**2)
0978         ncc(1,1)=nj
0979         ncc(2,1)=0
0980 
0981       else               !gluon (soft sea) parton on projectile side
0982 
0983 
0984 c define the 2 string ends from sea
0985 
0986         iflq=idtra(icp1,1,0,2)
0987         iflqb=idtra(icp2,1,0,2)
0988 
0989 
0990         if(ish.ge.5)write(ifch,*)'flavor sq+,sqb+',iflq,iflqb
0991 
0992         nj=nj+1
0993 c put the sea quarks (or diquarks) into the parton list as string end
0994         if(abs(iflqb).eq.3)then
0995           iqj(nj)=iflqb*4/3
0996         elseif(abs(iflqb).ge.4)then
0997           iqj(nj)=iflqb*10
0998         else
0999           iqj(nj)=iflqb
1000         endif
1001         ioj(nj)=7
1002         if(abs(iflq).eq.3)then
1003           iqj(nj+1)=iflq*4/3
1004         elseif(abs(iflq).ge.4)then
1005           iqj(nj+1)=iflq*10
1006         else
1007           iqj(nj+1)=iflq
1008         endif
1009         ioj(nj+1)=7
1010 
1011         eqj(1,nj)=.5*sngl(wp1+dble(amqt(1))**2/wp1)
1012         eqj(2,nj)=wp1-eqj(1,nj)
1013         eqj(3,nj)=xxp1pr(ncolp,kcol)
1014         eqj(4,nj)=xyp1pr(ncolp,kcol)
1015         if(ish.ge.8)write (ifch,*)'q_s1+ mass check',(eqj(1,nj)-
1016      *    eqj(2,nj))*(eqj(1,nj)+eqj(2,nj))-eqj(3,nj)**2-eqj(4,nj)**2
1017         eqj(1,nj)=sqrt(eqj(2,nj)**2+eqj(3,nj)**2+eqj(4,nj)**2)
1018         eqj(1,nj+1)=.5*sngl(wp2+dble(amqt(2))**2/wp2)
1019         eqj(2,nj+1)=wp2-eqj(1,nj+1)
1020         eqj(3,nj+1)=xxp2pr(ncolp,kcol)
1021         eqj(4,nj+1)=xyp2pr(ncolp,kcol)
1022         nj=nj+1
1023         if(ish.ge.8)write (ifch,*)'q_s2+ mass check',(eqj(1,nj)-
1024      *    eqj(2,nj))*(eqj(1,nj)+eqj(2,nj))-eqj(3,nj)**2-eqj(4,nj)**2
1025         eqj(1,nj)=sqrt(eqj(2,nj)**2+eqj(3,nj)**2+eqj(4,nj)**2)
1026 
1027 c gluon initiate space like cascade
1028         if(jqq.eq.0.or.iqq.eq.0.and.jqq.eq.2)then   
1029           iqc(1)=0
1030           ncc(1,1)=nj-1
1031           ncc(2,1)=nj
1032 c quark initiate space like cascade
1033         else
1034 c choose the first parton for the space like cascade (projectile)
1035           iqc(1)=int(3*rangen()+1.)*(2.*int(.5+rangen())-1.) 
1036           if(1.-1./(wp0*(wmi1-smin/wpi1)).le.xp)goto 1
1037 7         zg=1d0-dble(rangen())*(1.d0-xp)
1038           if(ish.ge.8)write (ifch,*)'6-zg,xp',zg,xp
1039           if(dble(rangen()).gt.zg**dels*((1.d0-xp/zg)/(1.d0-xp))
1040      *                                              **betpom)goto 7
1041           xg=xp/zg
1042           wpq=wp0*(xg-xp)
1043           wmq=1.d0/wpq
1044           wmi1=wmi1-wmq
1045           if(wmi1*wpi1.le.smin)goto 1
1046 
1047 
1048 c add the corresponding anti-parton in the list to compensate the emitted one
1049           nj=nj+1
1050           if(iabs(iqc(1)).eq.3)then
1051             iqj(nj)=-iqc(1)*4/3
1052           elseif(iabs(iqc(1)).ge.3)then
1053             iqj(nj)=-iqc(1)*10
1054           else
1055             iqj(nj)=-iqc(1)
1056           endif
1057           ioj(nj)=-7
1058           eqj(1,nj)=.5*wmq
1059           eqj(2,nj)=-eqj(1,nj)
1060           eqj(3,nj)=0.
1061           eqj(4,nj)=0.
1062           if(ish.ge.8)write (ifch,*)'q_s3+ mass check',eqj(1,nj)**2-
1063      *      eqj(2,nj)**2-eqj(3,nj)**2-eqj(4,nj)**2
1064           if(iqc(1).gt.0)then
1065             ncj(1,nj)=nj-1
1066             ncj(1,nj-1)=nj
1067             ncj(2,nj)=0
1068             ncj(2,nj-1)=0
1069             ncc(1,1)=nj-2
1070             ncc(2,1)=0
1071           else
1072             ncj(1,nj)=nj-2
1073             ncj(1,nj-2)=nj
1074             ncj(2,nj)=0
1075             ncj(2,nj-2)=0
1076             ncc(1,1)=nj-1
1077             ncc(2,1)=0
1078           endif
1079         endif
1080       endif
1081 
1082 
1083 
1084       if((iqq-2)*(iqq-3).eq.0)then    !valence quark on target side
1085 
1086 c define the second string end from sea
1087         if(iqc(2).gt.0)then
1088           ifl=idtra(icm2,1,0,2)
1089         elseif(iqc(2).lt.0)then
1090           ifl=idtra(icm1,1,0,2)
1091         else
1092           call utstop('No quark for hard Pomeron- in psahot!&')
1093         endif
1094 
1095         if(ish.ge.5)write(ifch,*)'flavor vq-,sqb-',iqc(2)
1096      &                                            ,ifl
1097 
1098         nj=nj+1
1099 c put the sea quark (or diquark) into the parton list as string end
1100         if(abs(ifl).eq.3)then
1101           iqj(nj)=ifl*4/3
1102         elseif(abs(ifl).ge.4)then
1103           iqj(nj)=ifl*10
1104         else
1105           iqj(nj)=ifl
1106         endif
1107         ioj(nj)=8
1108 
1109         eqj(1,nj)=.5*sngl(wm1+dble(amqt(4))**2/wm1)
1110         eqj(2,nj)=eqj(1,nj)-sngl(wm1)
1111         eqj(3,nj)=xxm1pr(ncolp,kcol)
1112         eqj(4,nj)=xym1pr(ncolp,kcol)
1113         if(ish.ge.8)write (ifch,*)'q_v- mass check',(eqj(1,nj)
1114      *    +eqj(2,nj))*(eqj(1,nj)-eqj(2,nj))-eqj(3,nj)**2-eqj(4,nj)**2
1115         eqj(1,nj)=sqrt(eqj(2,nj)**2+eqj(3,nj)**2+eqj(4,nj)**2)
1116         ncc(1,2)=nj
1117         ncc(2,2)=0
1118 
1119       else               !gluon (soft sea) parton on target side
1120 
1121 c define the 2 string ends from sea
1122 
1123         iflq=idtra(icm1,1,0,2)
1124         iflqb=idtra(icm2,1,0,2)
1125 
1126 
1127         if(ish.ge.5)write(ifch,*)'flavor sq-,sqb-',iflq,iflqb
1128 
1129         nj=nj+1
1130 c put the sea quarks (or diquarks) into the parton list as string end
1131         if(abs(iflqb).eq.3)then
1132           iqj(nj)=iflqb*4/3
1133         elseif(abs(iflqb).ge.4)then
1134           iqj(nj)=iflqb*10
1135         else
1136           iqj(nj)=iflqb
1137         endif
1138         ioj(nj)=8
1139         if(abs(iflq).eq.3)then
1140           iqj(nj+1)=iflq*4/3
1141         elseif(abs(iflq).ge.4)then
1142           iqj(nj+1)=iflq*10
1143         else
1144           iqj(nj+1)=iflq
1145         endif
1146         ioj(nj+1)=8
1147 
1148         eqj(1,nj)=.5*sngl(wm1+dble(amqt(3))**2/wm1)
1149         eqj(2,nj)=eqj(1,nj)-sngl(wm1)
1150         eqj(3,nj)=xxm2pr(ncolp,kcol)
1151         eqj(4,nj)=xym2pr(ncolp,kcol)
1152         if(ish.ge.8)write (ifch,*)'q_s1- mass check',(eqj(1,nj)-
1153      *    eqj(2,nj))*(eqj(1,nj)+eqj(2,nj))-eqj(3,nj)**2-eqj(4,nj)**2
1154         eqj(1,nj)=sqrt(eqj(2,nj)**2+eqj(3,nj)**2+eqj(4,nj)**2)
1155         eqj(1,nj+1)=.5*sngl(wm2+dble(amqt(4))**2/wm2)
1156         eqj(2,nj+1)=eqj(1,nj+1)-wm2
1157         eqj(3,nj+1)=xxm1pr(ncolp,kcol)
1158         eqj(4,nj+1)=xym1pr(ncolp,kcol)
1159         nj=nj+1
1160         if(ish.ge.8)write (ifch,*)'q_s2- mass check',(eqj(1,nj)-
1161      *    eqj(2,nj))*(eqj(1,nj)+eqj(2,nj))-eqj(3,nj)**2-eqj(4,nj)**2
1162         eqj(1,nj)=sqrt(eqj(2,nj)**2+eqj(3,nj)**2+eqj(4,nj)**2)
1163 
1164 c gluon initiate space like cascade
1165         if(jqq.eq.0.or.iqq.eq.0.and.jqq.eq.1)then
1166           iqc(2)=0
1167           ncc(1,2)=nj-1
1168           ncc(2,2)=nj
1169 c quark initiate space like cascade
1170         else
1171 c choose the first parton for the space like cascade (target)
1172           iqc(2)=int(3*rangen()+1.)*(2.*int(.5+rangen())-1.)
1173           if(1.-1./(wm0*(wpi1-smin/wmi1)).le.xm)goto 1
1174 8           zg=1.d0-dble(rangen())*(1.d0-xm)
1175           if(ish.ge.8)write (ifch,*)'7-zg,xm',zg,xm
1176           if(rangen().gt.zg**dels*((1.d0-xm/zg)/(1.d0-xm))**betpom)
1177      *      goto 8
1178           xg=xm/zg
1179           wmq=wm0*(xg-xm)
1180           wpq=1.d0/wmq
1181           wpi1=wpi1-wpq
1182           if(wmi1*wpi1.le.smin)goto 1
1183 
1184 c add the corresponding anti-parton in the list to compensate the emitted one
1185           nj=nj+1
1186           if(iabs(iqc(2)).eq.3)then
1187             iqj(nj)=-iqc(2)*4/3
1188           elseif(iabs(iqc(2)).ge.4)then
1189             iqj(nj)=-iqc(2)*10
1190           else
1191             iqj(nj)=-iqc(2)
1192           endif
1193           ioj(nj)=-8
1194 
1195           eqj(1,nj)=.5*sngl(wpq)
1196           eqj(2,nj)=eqj(1,nj)
1197           eqj(3,nj)=0.
1198           eqj(4,nj)=0.
1199           if(ish.ge.8)write (ifch,*)'q_s3- mass check',(eqj(1,nj)-
1200      *      eqj(2,nj))*(eqj(1,nj)+eqj(2,nj))-eqj(3,nj)**2-eqj(4,nj)**2
1201           if(iqc(2).gt.0)then
1202             ncj(1,nj)=nj-1
1203             ncj(1,nj-1)=nj
1204             ncj(2,nj)=0
1205             ncj(2,nj-1)=0
1206             ncc(1,2)=nj-2
1207             ncc(2,2)=0
1208           else
1209             ncj(1,nj)=nj-2
1210             ncj(1,nj-2)=nj
1211             ncj(2,nj)=0
1212             ncj(2,nj-2)=0
1213             ncc(1,2)=nj-1
1214             ncc(2,2)=0
1215           endif
1216         endif
1217       endif
1218 
1219       else
1220 
1221 c put partons in stack for cascades
1222 
1223       if((iqq-1)*(iqq-3).eq.0)then    !valence quark on projectile side
1224 c valence quark used to start space like cascade
1225         iqc(1)=iq1              !proj=particle
1226         if(iabs(idptl(ip)).eq.1220)iqc(1)=3-iq1 !proj=neutron
1227         if(idptl(ip).lt.0)iqc(1)=-iqc(1) !proj=antiparticle
1228         ifl1=iabs(iqc(1))
1229         
1230         nj=nj+1
1231 
1232         if(iqc(1).gt.0)then
1233           if(iremn.ne.0)then
1234             if(iremn.eq.3)then
1235               call idsufl3(ifl1,1,jcp) !remove valence quark from remnant
1236               ifl=idrafl(iclpro,jcpr,2,'v',1,iret3) !take sea antiquark from jcpr
1237             elseif(iremn.eq.2)then
1238               call idsufl3(ifl1,1,jcp) !remove valence quark from remnant
1239               ifl=idrafl(iclpro,jcpr,2,'s',0,iret3) !take sea antiquark freely
1240               call idsufl3(ifl,2,jcpr) !remove sea antiquark from remnant
1241             else
1242               call idsufl(ifl1,1,jcp,iret1) !remove valence quark from remnant
1243               if(iret1.ne.1)then
1244                 ifl=idrafl(iclpro,jcp,2,'s',1,iret2) !take sea antiquark
1245               endif
1246             endif
1247           else
1248             ifl=ifl1
1249           endif
1250         elseif(iqc(1).lt.0)then
1251           if(iremn.ne.0)then
1252             if(iremn.eq.3)then
1253               call idsufl3(ifl1,2,jcp)
1254               ifl=idrafl(iclpro,jcpr,1,'v',1,iret3) !take sea quark
1255             elseif(iremn.eq.2)then
1256               call idsufl3(ifl1,2,jcp)
1257               ifl=idrafl(iclpro,jcpr,1,'s',0,iret3) !take sea quark
1258               call idsufl3(ifl,1,jcpr) !remove sea quark from remnant
1259             else
1260               call idsufl(ifl1,2,jcp,iret1)
1261               if(iret1.ne.1)then
1262                 ifl=idrafl(iclpro,jcp,1,'s',1,iret2)
1263               endif
1264             endif
1265           else
1266             ifl=ifl1
1267           endif
1268         else
1269           call utstop('No quark for hard Pomeron+ in psahot!&')
1270         endif
1271 
1272         if(iret1.eq.1.or.iret2.eq.1)then
1273           ifl=ifl1
1274          if(ish.ge.3)write(ifmt,*)'Not enough space in rem (psahot 1)'
1275          if(ish.ge.5)write(ifch,*)'Not enough space in rem (psahot 1)'
1276           call iddeco(icp,jcp)
1277           iret1=0
1278           iret2=0
1279         endif
1280 
1281         if(ish.ge.5)write(ifch,*)'flavor vq+,sqb+',iqc(1)
1282      &                                            ,-isign(ifl,iqc(1))
1283 
1284         if(ifl.eq.3)then
1285           iqj(nj)=-isign(ifl,iqc(1))*4/3
1286         elseif(ifl.eq.4)then
1287           iqj(nj)=-isign(ifl,iqc(1))/4
1288         else
1289           iqj(nj)=-isign(ifl,iqc(1))
1290         endif
1291         eqj(1,nj)=.5*sngl(wp1+dble(amqt(2))**2/wp1)
1292         eqj(2,nj)=wp1-eqj(1,nj)
1293         eqj(3,nj)=xxp2pr(ncolp,kcol)
1294         eqj(4,nj)=xyp2pr(ncolp,kcol)
1295         if(ish.ge.8)write (ifch,*)'q_v+ mass check',(eqj(1,nj)-
1296      *    eqj(2,nj))*(eqj(1,nj)+eqj(2,nj))-eqj(3,nj)**2-eqj(4,nj)**2
1297         eqj(1,nj)=sqrt(eqj(2,nj)**2+eqj(3,nj)**2+eqj(4,nj)**2)
1298         ncc(1,1)=nj
1299         ncc(2,1)=0
1300 
1301       else
1302 
1303         nj=nj+1
1304         if(idfpomr.lt.3.and.iremn.ne.0)then
1305           if(iremn.eq.3)then
1306             iflq=idrafl(iclpro,jcpr,1,'v',1,iret3) !take sea quark
1307             iflqb=idrafl(iclpro,jcpr,2,'v',1,iret3) !take sea antiquark
1308           elseif(iremn.eq.2)then
1309             iflq=idrafl(iclpro,jcpr,1,'s',0,iret3) !take sea quark
1310             call idsufl3(iflq,1,jcpr) !remove sea quark from remnant
1311             iflqb=idrafl(iclpro,jcpr,2,'s',0,iret3) !take sea antiquark
1312             call idsufl3(iflqb,2,jcpr) !remove sea antiquark from remnant
1313           else
1314             iflq=idrafl(iclpro,jcp,1,'s',1,iret1) !Choose flavor of sea quark.
1315             if(iret1.ne.1)iflqb=idrafl(iclpro,jcp,2,'s',1,iret2) !Choose antiquark.
1316           endif
1317         else
1318           iflq=idrafl(iclpro,jcp,0,'s',0,iret1)   !Choose flavor of sea quark.
1319           iflqb=iflq                   !antiquark=quark (vertex end)
1320         endif
1321 
1322         if(iret1.eq.1.or.iret2.eq.1)then
1323           iflqb=iflq
1324          if(ish.ge.3)write(ifmt,*)'Not enough space in rem (psahot 2)'
1325          if(ish.ge.5)write(ifch,*)'Not enough space in rem (psahot 2)'
1326           call iddeco(icp,jcp)
1327           iret1=0
1328           iret2=0
1329         endif
1330 
1331         if(ish.ge.5)write(ifch,*)'flavor sq+,sqb+',iflq,-iflqb
1332 
1333         if(iflqb.eq.3)then
1334           iqj(nj)=-iflqb*4/3
1335         elseif(iflqb.eq.4)then
1336           iqj(nj)=-iflqb/4
1337         else
1338           iqj(nj)=-iflqb
1339         endif
1340         ioj(nj)=7
1341         if(iflq.eq.3)then
1342           iqj(nj+1)=iflq*4/3
1343         elseif(iflq.eq.4)then
1344           iqj(nj+1)=iflq/4
1345         else
1346           iqj(nj+1)=iflq
1347         endif
1348         ioj(nj+1)=7
1349 
1350         eqj(1,nj)=.5*sngl(wp1+dble(amqt(1))**2/wp1)
1351         eqj(2,nj)=wp1-eqj(1,nj)
1352         eqj(3,nj)=xxp1pr(ncolp,kcol)
1353         eqj(4,nj)=xyp1pr(ncolp,kcol)
1354         if(ish.ge.8)write (ifch,*)'q_s1+ mass check',(eqj(1,nj)-
1355      *    eqj(2,nj))*(eqj(1,nj)+eqj(2,nj))-eqj(3,nj)**2-eqj(4,nj)**2
1356         eqj(1,nj)=sqrt(eqj(2,nj)**2+eqj(3,nj)**2+eqj(4,nj)**2)
1357         eqj(1,nj+1)=.5*sngl(wp2+dble(amqt(2))**2/wp2)
1358         eqj(2,nj+1)=wp2-eqj(1,nj+1)
1359         eqj(3,nj+1)=xxp2pr(ncolp,kcol)
1360         eqj(4,nj+1)=xyp2pr(ncolp,kcol)
1361         nj=nj+1
1362         if(ish.ge.8)write (ifch,*)'q_s2+ mass check',(eqj(1,nj)-
1363      *    eqj(2,nj))*(eqj(1,nj)+eqj(2,nj))-eqj(3,nj)**2-eqj(4,nj)**2
1364         eqj(1,nj)=sqrt(eqj(2,nj)**2+eqj(3,nj)**2+eqj(4,nj)**2)
1365 
1366 c gluon initiate space like cascade
1367         if(jqq.eq.0.or.iqq.eq.0.and.jqq.eq.2)then
1368           iqc(1)=0
1369           ncc(1,1)=nj-1
1370           ncc(2,1)=nj
1371 c quark initiate space like cascade
1372         else
1373 c choose the first parton for the space like cascade (projectile)
1374           iqc(1)=int(3*rangen()+1.)*(2.*int(.5+rangen())-1.)
1375           if(1.-1./(wp0*(wmi1-smin/wpi1)).le.xp)goto 1
1376  17       zg=1d0-dble(rangen())*(1.d0-xp)
1377           if(ish.ge.8)write (ifch,*)'6-zg,xp',zg,xp
1378           if(dble(rangen()).gt.zg**dels*((1.d0-xp/zg)/(1.d0-xp))
1379      *                                              **betpom)goto 17
1380           xg=xp/zg
1381           wpq=wp0*(xg-xp)
1382           wmq=1.d0/wpq
1383           wmi1=wmi1-wmq
1384           if(wmi1*wpi1.le.smin)goto 1
1385 
1386 
1387 c add the corresponding anti-parton in the list to compensate the emitted one
1388           nj=nj+1
1389           if(iabs(iqc(1)).eq.3)then
1390             iqj(nj)=-iqc(1)*4/3
1391           else
1392             iqj(nj)=-iqc(1)
1393           endif
1394           ioj(nj)=-7
1395           eqj(1,nj)=.5*wmq
1396           eqj(2,nj)=-eqj(1,nj)
1397           eqj(3,nj)=0.
1398           eqj(4,nj)=0.
1399           if(ish.ge.8)write (ifch,*)'q_s3+ mass check',eqj(1,nj)**2-
1400      *      eqj(2,nj)**2-eqj(3,nj)**2-eqj(4,nj)**2
1401           if(iqc(1).gt.0)then
1402             ncj(1,nj)=nj-1
1403             ncj(1,nj-1)=nj
1404             ncj(2,nj)=0
1405             ncj(2,nj-1)=0
1406             ncc(1,1)=nj-2
1407             ncc(2,1)=0
1408           else
1409             ncj(1,nj)=nj-2
1410             ncj(1,nj-2)=nj
1411             ncj(2,nj)=0
1412             ncj(2,nj-2)=0
1413             ncc(1,1)=nj-1
1414             ncc(2,1)=0
1415           endif
1416         endif
1417       endif
1418 
1419       if((iqq-2)*(iqq-3).eq.0)then
1420         iqc(2)=iq2              !tar=particle (can not be pion or kaon !)
1421         if(iabs(idptl(maproj+it)).eq.1220)iqc(2)=3-iq2 !targ=neutron
1422         if(idptl(maproj+it).lt.0)iqc(2)=-iqc(2) !targ=antiparticle
1423         ifl2=iabs(iqc(2))
1424         
1425         nj=nj+1
1426         if(iqc(2).gt.0)then
1427           if(iremn.ne.0)then
1428             if(iremn.eq.3)then
1429               call idsufl3(ifl2,1,jct) !remove valence quark from remnant
1430               ifl=idrafl(icltar,jctr,2,'v',1,iret3) !take sea antiquark from jctr
1431             elseif(iremn.eq.2)then
1432               call idsufl3(ifl2,1,jct) !remove valence quark from remnant
1433               ifl=idrafl(icltar,jctr,2,'s',0,iret3) !take sea antiquark freely
1434               call idsufl3(ifl,2,jctr) !remove sea antiquark from remnant
1435             else
1436               call idsufl(ifl2,1,jct,iret1) !remove valence quark from remnant
1437               if(iret1.ne.1)then
1438                 ifl=idrafl(icltar,jct,2,'s',1,iret2) !take sea antiquark
1439               endif
1440             endif
1441           else
1442             ifl=ifl2
1443           endif
1444         elseif(iqc(2).lt.0)then
1445           if(iremn.ne.0)then
1446             if(iremn.eq.3)then
1447               call idsufl3(ifl2,2,jct)
1448               ifl=idrafl(icltar,jctr,1,'v',1,iret3) !take sea quark
1449             elseif(iremn.eq.2)then
1450               call idsufl3(ifl2,2,jct)
1451               ifl=idrafl(icltar,jctr,1,'s',0,iret3) !take sea quark
1452               call idsufl3(ifl,1,jctr) !remove sea quark from remnant
1453             else
1454               call idsufl(ifl2,2,jct,iret1)
1455               if(iret1.ne.1)then
1456                 ifl=idrafl(icltar,jct,1,'s',1,iret2)
1457               endif
1458             endif
1459           else
1460             ifl=ifl2
1461           endif
1462         else
1463           call utstop('No quark for hard Pomeron- in psahot!&')
1464         endif
1465 
1466         if(iret1.eq.1.or.iret2.eq.1)then
1467           ifl=ifl2
1468          if(ish.ge.3)write(ifmt,*)'Not enough space in rem (psahot 3)'
1469          if(ish.ge.5)write(ifch,*)'Not enough space in rem (psahot 3)'
1470           call iddeco(ict,jct)
1471           iret1=0
1472           iret2=0
1473         endif
1474 
1475         if(ish.ge.5)write(ifch,*)'flavor vq-,sqb-',iqc(2)
1476      &                                            ,-isign(ifl,iqc(2))
1477 
1478         if(ifl.eq.3)then
1479           iqj(nj)=-isign(ifl,iqc(2))*4/3
1480         elseif(ifl.eq.4)then
1481           iqj(nj)=-isign(ifl,iqc(2))/4
1482         else
1483           iqj(nj)=-isign(ifl,iqc(2))
1484         endif
1485 
1486         eqj(1,nj)=.5*sngl(wm1+dble(amqt(4))**2/wm1)
1487         eqj(2,nj)=eqj(1,nj)-sngl(wm1)
1488         eqj(3,nj)=xxm1pr(ncolp,kcol)
1489         eqj(4,nj)=xym1pr(ncolp,kcol)
1490         if(ish.ge.8)write (ifch,*)'q_v- mass check',(eqj(1,nj)
1491      *    +eqj(2,nj))*(eqj(1,nj)-eqj(2,nj))-eqj(3,nj)**2-eqj(4,nj)**2
1492         eqj(1,nj)=sqrt(eqj(2,nj)**2+eqj(3,nj)**2+eqj(4,nj)**2)
1493         ncc(1,2)=nj
1494         ncc(2,2)=0
1495       else
1496         nj=nj+1
1497         if(mod(idfpomr,2).ne.0.and.iremn.ne.0)then
1498           if(iremn.eq.3)then
1499             iflq=idrafl(icltar,jctr,1,'v',1,iret3) !take sea quark
1500             iflqb=idrafl(icltar,jctr,2,'v',1,iret3) !take sea antiquark
1501           elseif(iremn.eq.2)then
1502             iflq=idrafl(icltar,jctr,1,'s',0,iret3) !take sea quark
1503             call idsufl3(iflq,1,jctr) !remove sea antiquark from remnant
1504             iflqb=idrafl(icltar,jctr,2,'s',0,iret3) !take sea antiquark
1505             call idsufl3(iflqb,2,jctr) !remove sea antiquark from remnant
1506           else
1507             iflq=idrafl(icltar,jct,1,'s',1,iret2) !Choose flavor of sea quark.
1508             if(iret1.ne.1)iflqb=idrafl(icltar,jct,2,'s',1,iret2) !Choose  antiquark.
1509           endif
1510         else
1511           iflq=idrafl(iclpro,jcp,0,'s',0,iret1)   !Choose flavor of sea quark.
1512           iflqb=iflq                   !antiquark=quark (vertex end)
1513         endif
1514 
1515         if(iret1.eq.1.or.iret2.eq.1)then
1516           iflqb=iflq
1517          if(ish.ge.3)write(ifmt,*)'Not enough space in rem (psahot 4)'
1518          if(ish.ge.5)write(ifch,*)'Not enough space in rem (psahot 4)'
1519           call iddeco(ict,jct)
1520           iret1=0
1521           iret2=0
1522         endif
1523 
1524         if(ish.ge.5)write(ifch,*)'flavor sq-,sqb-',iflq,-iflqb
1525 
1526         if(iflqb.eq.3)then
1527           iqj(nj)=-iflqb*4/3
1528         elseif(iflqb.eq.4)then
1529           iqj(nj)=-iflqb/4
1530         else
1531           iqj(nj)=-iflqb
1532         endif
1533         if(iflq.eq.3)then
1534           iqj(nj+1)=iflq*4/3
1535         elseif(iflq.eq.4)then
1536           iqj(nj+1)=iflq/4
1537         else
1538           iqj(nj+1)=iflq
1539         endif
1540         ioj(nj+1)=8
1541 
1542         eqj(1,nj)=.5*sngl(wm1+dble(amqt(3))**2/wm1)
1543         eqj(2,nj)=eqj(1,nj)-sngl(wm1)
1544         eqj(3,nj)=xxm2pr(ncolp,kcol)
1545         eqj(4,nj)=xym2pr(ncolp,kcol)
1546         if(ish.ge.8)write (ifch,*)'q_s1- mass check',(eqj(1,nj)-
1547      *    eqj(2,nj))*(eqj(1,nj)+eqj(2,nj))-eqj(3,nj)**2-eqj(4,nj)**2
1548         eqj(1,nj)=sqrt(eqj(2,nj)**2+eqj(3,nj)**2+eqj(4,nj)**2)
1549         eqj(1,nj+1)=.5*sngl(wm2+dble(amqt(4))**2/wm2)
1550         eqj(2,nj+1)=eqj(1,nj+1)-wm2
1551         eqj(3,nj+1)=xxm1pr(ncolp,kcol)
1552         eqj(4,nj+1)=xym1pr(ncolp,kcol)
1553         nj=nj+1
1554         if(ish.ge.8)write (ifch,*)'q_s2- mass check',(eqj(1,nj)-
1555      *    eqj(2,nj))*(eqj(1,nj)+eqj(2,nj))-eqj(3,nj)**2-eqj(4,nj)**2
1556         eqj(1,nj)=sqrt(eqj(2,nj)**2+eqj(3,nj)**2+eqj(4,nj)**2)
1557 
1558 c gluon initiate space like cascade
1559         if(jqq.eq.0.or.iqq.eq.0.and.jqq.eq.1)then
1560           iqc(2)=0
1561           ncc(1,2)=nj-1
1562           ncc(2,2)=nj
1563 c quark initiate space like cascade
1564         else
1565 c choose the first parton for the space like cascade (target)
1566           iqc(2)=int(3*rangen()+1.)*(2.*int(.5+rangen())-1.)
1567           if(1.-1./(wm0*(wpi1-smin/wmi1)).le.xm)goto 1
1568  18       zg=1.d0-dble(rangen())*(1.d0-xm)
1569           if(ish.ge.8)write (ifch,*)'7-zg,xm',zg,xm
1570           if(rangen().gt.zg**dels*((1.d0-xm/zg)/(1.d0-xm))**betpom)
1571      *      goto 18
1572           xg=xm/zg
1573           wmq=wm0*(xg-xm)
1574           wpq=1.d0/wmq
1575           wpi1=wpi1-wpq
1576           if(wmi1*wpi1.le.smin)goto 1
1577 
1578 c add the corresponding anti-parton in the list to compensate the emitted one
1579           nj=nj+1
1580           if(iabs(iqc(2)).eq.3)then
1581             iqj(nj)=-iqc(2)*4/3
1582           else
1583             iqj(nj)=-iqc(2)
1584           endif
1585           ioj(nj)=-8
1586 
1587           eqj(1,nj)=.5*sngl(wpq)
1588           eqj(2,nj)=eqj(1,nj)
1589           eqj(3,nj)=0.
1590           eqj(4,nj)=0.
1591           if(ish.ge.8)write (ifch,*)'q_s3- mass check',(eqj(1,nj)-
1592      *      eqj(2,nj))*(eqj(1,nj)+eqj(2,nj))-eqj(3,nj)**2-eqj(4,nj)**2
1593           if(iqc(2).gt.0)then
1594             ncj(1,nj)=nj-1
1595             ncj(1,nj-1)=nj
1596             ncj(2,nj)=0
1597             ncj(2,nj-1)=0
1598             ncc(1,2)=nj-2
1599             ncc(2,2)=0
1600           else
1601             ncj(1,nj)=nj-2
1602             ncj(1,nj-2)=nj
1603             ncj(2,nj)=0
1604             ncj(2,nj-2)=0
1605             ncc(1,2)=nj-1
1606             ncc(2,2)=0
1607           endif
1608         endif
1609       endif
1610 
1611       endif        !iLHC
1612 
1613       if(jqq.ne.0)then
1614         if(iqq.ne.0.or.iqq.eq.0.and.jqq.eq.3)then
1615           if(iclpro.ne.4.or.iqq.ne.1)then
1616             call psjti0(sngl(wpi1*wmi1-pth),sqq1,sqqb1,1,1)
1617             call psjti0(sngl(wpi1*wmi1-pth),sqaq1,sqaqb1,-1,1)
1618             call psjti0(sngl(wpi1*wmi1-pth),sqqp1,sqqpb1,1,2)
1619             sqq1=(sqq1+sqaq1+2.*(naflav-1)*sqqp1)/naflav/2.
1620           else
1621             call psjti0(sngl(wpi1*wmi1-pth),sqq1,sqqb1,4,1)
1622           endif
1623           gbs=sqq1/sqq
1624         else
1625           call psjti0(sngl(wpi1*wmi1-pth),sgq1,sgqb1,0,1)
1626           gbs=sgq1/sgq
1627         endif
1628         if(ish.ge.8)write (ifch,*)'gbs',gbs
1629         if(rangen().gt.gbs)goto 6
1630       endif
1631 
1632 
1633 c---------------------------------------------------------------
1634 c        inner partons of the hard Pomeron
1635 c---------------------------------------------------------------
1636 
1637        if(ish.ge.4)write(ifch,*)
1638      &  'inner partons of the hard Pomeron'
1639       nj00=nj
1640       wpi=wpi1
1641       wmi=wmi1
1642       si=wpi*wmi-pxh**2-pyh**2     !total energy squared for the hard
1643       if(ish.ge.7)write (ifch,*)'si,wpi,wmi',si,wpi,wmi
1644 
1645       ept(1)=.5d0*(wpi+wmi)
1646       ept(2)=.5d0*(wpi-wmi)
1647       ept(3)=pxh
1648       ept(4)=pyh
1649       qmin(1)=q2min          !effective momentum cutoff above current la
1650       qmin(2)=q2min          !effective momentum cutoff below current la
1651       qminn=max(qmin(1),qmin(2)) !effective momentum cutoff for the bo
1652 c        si=psnorm(ept)       !4-momentum squared for the hard pomeron
1653       jfirst=1
1654       jj=int(1.5+rangen())
1655       nemis(1)=0
1656       nemis(2)=0
1657 
1658   9   continue ! <<<<----------- ladder rung ---------------------------
1659 
1660       if(ish.ge.4)write(ifch,*)'ladder rung'
1661       pt2=ept(3)**2+ept(4)**2
1662       pt=sqrt(sngl(pt2))
1663       if(iabs(iqc(1)).ne.4)then
1664         q2mass=0.
1665       else
1666         q2mass=qcmass**2
1667         si=si-dble(q2mass)
1668       endif
1669       s2min=4.*qminn+q2mass         !mass cutoff for born scattering
1670       wwmin=5.*qminn+q2mass-2.*pt*sqrt(q2ini)
1671       wwmin=(wwmin+sqrt(wwmin**2+4.*(q2mass+pt2)*(qminn-q2ini)))
1672      */(1.-q2ini/qminn)/2.
1673       if(ish.ge.5)write(ifch,*)'qminn,q2mass,pt,wwmin,iqc(1):'
1674       if(ish.ge.5)write(ifch,*)qminn,q2mass,pt,wwmin
1675 
1676       wpt(1)=ept(1)+ept(2)            !lc+ for the current jet emissi
1677       wpt(2)=ept(1)-ept(2)            !lc- for the current jet emissi
1678       
1679       sjord=0.
1680       if(iabs(iqc(1)).ne.4)then
1681         if(jfirst.eq.1)then
1682           sj=psjti(qmin(jj),qqcut,sngl(si),iqc(jj),iqc(3-jj),0)
1683           sj2=psjti1(qmin(3-jj),qmin(jj),qqcut
1684      *              ,sngl(si),iqc(3-jj),iqc(jj),0)
1685           if(ish.ge.5)write(ifch,*)'si,sj,sj2:',si,sj,sj2
1686           if(rangen().gt.sj2/sj.and.si.gt.1.1d0*dble(wwmin))goto 112
1687           jfirst=0
1688           jj=3-jj
1689           sj=sj2
1690           goto 111
1691         endif
1692        sj=psjti1(qmin(jj),qmin(3-jj),qqcut,sngl(si),iqc(jj),iqc(3-jj),0)
1693 111     sjb=psbint(qmin(1),qmin(2),qqcut,sngl(si),iqc(1),iqc(2),0) !born parton-parton
1694       else
1695         sjord=psjci(qmin(2),sngl(si),iqc(2))
1696         sj=sjord
1697       sjb=psbint(qmin(1),qmin(2),qqcut,sngl(si)+q2mass,iqc(1),iqc(2),0)
1698         if(qmin(2).eq.q2min)then
1699           wwmins=2.5*q2min*(1.+sqrt(1.+
1700      *    4.*(pt2+q2mass)/q2min))
1701        if(si.gt.dble(wwmins))call psjti0(sngl(si)+q2mass,sj,sjb
1702      *                                                  ,iqc(1),iqc(2))
1703         endif
1704       endif
1705       if(ish.ge.5)write(ifch,*)'si,pt2,wwmin,s2min,sj,sjb,iqc:'
1706       if(ish.ge.5)write(ifch,*)si,pt2,wwmin,s2min,sj,sjb,iqc
1707 
1708       if(si.lt.1.1d0*dble(wwmin))goto 12 !------>>>>>>>
1709       if(rangen().lt.sjb/sj)goto 12 !------>>>>>>>
1710 
1711       if(iabs(iqc(1)).eq.4)jj=min(2,int(sjord/sj+rangen())+1)  !?????????
1712 
1713 112   continue
1714 
1715       if(iabs(iqc(jj)).ne.4)then
1716 
1717         discr=dble((sngl(si)+2.*pt*sqrt(q2ini)-q2mass)**2
1718      *  -4.*q2ini*(5.*sngl(si)+q2mass+sngl(pt2)))
1719         if(discr.lt.0.d0.and.ish.ge.1)write(ifmt,*)'discr,si,pt,wwmin',
1720      *  discr,si,pt,wwmin
1721         discr=dsqrt(discr)
1722         qqmax=(si+2.d0*dble(pt*sqrt(q2ini))-dble(q2mass)+discr)/2.d0
1723      *  /(5.d0+(dble(q2mass)+pt2)/si)
1724         qqmin=qqmax-discr/(5.d0+(dble(q2mass)+pt2)/si)
1725         if(s2min.gt.4.d0*qqmin+dble(q2mass))then
1726           xmm=.5d0*(si-s2min+2.d0*dble(pt*sqrt(q2ini)))
1727           discr=xmm**2-si*dble(q2ini)*(1.d0+(dble(q2mass)+pt2)/si)
1728           if(discr.lt.0..and.ish.ge.1)write(ifmt,*)'discr1,si,pt,wwmin',
1729      *    discr,si,pt,wwmin
1730           qqmin=(xmm-dsqrt(discr))/(1.d0+(dble(q2mass)+pt2)/si)
1731         endif
1732 
1733         xmax=min(1.d0-dble(q2ini)/qqmax,.9999d0)
1734         if(qqmin.lt.dble(qmin(jj)))then
1735           qqmin=dble(qmin(jj))
1736           xmin=max(1.d0-((dble(pt)*dsqrt(qqmin)+dsqrt(pt2*qqmin+
1737      *    si*(si-s2min-qqmin*(1.d0+(dble(q2mass)+pt2)/si))))/si)**2,
1738      *    (s2min+qqmin*(1.d0+(dble(q2mass)+pt2)/si)-2.d0*dble(pt)
1739      *    *dsqrt(qqmin))/si)
1740           if(xmin.le.0.d0)xmin=s2min/si
1741         else
1742           xmin=1.d0-dble(q2ini)/qqmin
1743         endif
1744 
1745         qm0=qmin(jj)
1746         xm0=1.0-dble(q2ini/qm0)
1747         if(xm0.gt..999*xmax.or.xm0.lt.1.001*xmin)then
1748           xm0=.5d0*(xmax+xmin)
1749         endif
1750         s2max=sngl(xm0*si-qm0*(dble(q2mass)+pt2)/si)
1751      *       +2.*pt*sqrt(q2ini)
1752         xx=xm0
1753 
1754         if(iabs(iqc(1)).ne.4)then
1755           if(jfirst.eq.1)then
1756             sj0=psjti(qm0,qqcut,s2max,0,iqc(3-jj),0)*
1757      *      psfap(xx,iqc(jj),0)+
1758      *      psjti(qm0,qqcut,s2max,7,iqc(3-jj),0)
1759      *      *psfap(xx,iqc(jj),1)
1760           else
1761             sj0=psjti1(qm0,qmin(3-jj),qqcut,s2max,0,iqc(3-jj),0)*
1762      *      psfap(xx,iqc(jj),0)+
1763      *      psjti1(qm0,qmin(3-jj),qqcut,s2max,7,iqc(3-jj),0)
1764      *      *psfap(xx,iqc(jj),1)
1765           endif
1766         else
1767           sj0=psjci(qm0,s2max,0)*psfap(xx,iqc(jj),0)+
1768      *    psjci(qm0,s2max,1)*psfap(xx,iqc(jj),1)
1769         endif
1770 
1771         gb0=dble(sj0/log(q2ini/qcdlam)*qm0*5.)*psuds(qm0,iqc(jj))
1772         if(ish.ge.5)write(ifch,*)'gb0,qm0,xm0,s2max:',
1773      *  gb0,qm0,xm0,s2max
1774 c        if(gb0.le.0.)stop'gb0<=0'
1775         if(gb0.le.0.d0)then
1776           write(ifmt,*)'gb0.le.0.  si,pt2:',si,pt2
1777           iret=1
1778           goto 16
1779         endif
1780 
1781         if(xm0.le..5d0)then
1782           gb0=gb0*xm0**(1.-delh)
1783         else
1784           gb0=gb0*(1.d0-xm0)*2.d0**delh
1785         endif
1786 
1787         xmin2=max(.5d0,xmin)
1788         xmin1=xmin**delh
1789         xmax1=min(xmax,.5d0)**delh
1790         if(xmin.ge..5d0)then
1791           djl=1.
1792         elseif(xmax.le..5d0)then
1793           djl=0.
1794         else
1795           djl=1./(1.+((2.*sngl(xmin))**delh-1.)/delh/
1796      *    log(2.*(1.-sngl(xmax))))
1797         endif
1798 
1799       else
1800 
1801         xmin=5.d0*dble(q2min)/si
1802         xmax=min(si/(si+5.0*(pt2+dble(q2mass))),.9999d0)
1803         qqmax=xmax*si/5.d0
1804         qqmin=dble(q2min)
1805 
1806         qm0=sngl(qqmin)
1807         xm0=2.d0/(1.d0+sqrt(1.d0+4.d0*(pt2+dble(q2mass))/qm0))
1808         if(xm0.gt..999d0*xmax.or.xm0.lt.1.001d0*xmin)then
1809           xm0=.5d0*(xmax+xmin)
1810         endif
1811         s2max=sngl(xm0*si)
1812         xx=xm0
1813 
1814         sj0=psjti(qm0,qmin(3-jj),s2max,0,iqc(3-jj),0)*
1815      *  psfap(xx,iqc(jj),0)
1816         gb0=dble(sj0/log(qm0/qcdlam)*qm0   *5.)
1817         gb0=gb0*xm0**(1.-delh)
1818         if(gb0.le.0.d0)then
1819           if(ish.ge.2)write(ifch,*)'gb0.le.0. (charm)  si,pt2:',si,pt2
1820           iret=1
1821           goto 16
1822         endif
1823         djl=0.
1824         xmin2=0d0
1825         xmin1=xmin**delh
1826         xmax1=xmax**delh
1827 
1828       endif
1829 
1830       if(ish.ge.5)write(ifch,*)'xmin,xmax,qqmin,qqmax:',
1831      *xmin,xmax,qqmin,qqmax
1832 
1833       ntry=0
1834 10    continue
1835       ntry=ntry+1
1836       if(ntry.gt.5000000)then
1837         if(ish.ge.1)write(*,*)'Reject hard string (too many rejection)'
1838      &,kcol,ncolp,nptl,gb7
1839         iret=1
1840         goto 16
1841       endif
1842       if(rangen().gt.djl)then        !lc momentum share in the cur
1843         x=(xmin1+dble(rangen())*(xmax1-xmin1))**(1./delh)
1844       else
1845         x=1.d0-(1.d0-xmin2)*((1.d0-xmax)/(1.d0-xmin2))**rangen()
1846       endif
1847       qq=qqmin/(1.d0+dble(rangen())*(qqmin/qqmax-1.d0))
1848       if(ish.ge.7)write(ifch,*)'x,qq:',x,qq,ntry
1849 
1850       if(iabs(iqc(jj)).ne.4)then
1851 
1852         qt2=qq*(1.d0-x)
1853         if(qt2.lt.dble(q2ini))then
1854           if(ish.ge.7)write(ifch,*)'qt2:',qt2
1855           goto 10
1856         endif
1857 
1858         qmin2=max(qminn,sngl(qq))
1859         qt=dsqrt(qt2)
1860         call pscs(bcos,bsin)
1861         ep3(3)=sngl(qt)*bcos   !new parton pt-s
1862         ep3(4)=sngl(qt)*bsin
1863         ptnew=(ept(3)-dble(ep3(3)))**2+(ept(4)-dble(ep3(4)))**2
1864         s2min2=4.*qmin2+q2mass
1865 
1866         s2=x*(si-qq)-ptnew-qq*(dble(q2mass)+pt2)/si+pt2  !new ladder mass
1867         if(s2.lt.dble(s2min2))then
1868           if(ish.ge.7)write(ifch,*)'s2,s2min2:',s2,s2min2
1869           goto 10  !rejection in case of too low mass
1870         endif
1871 
1872         xx=x
1873         if(iabs(iqc(1)).ne.4)then
1874           if(jfirst.eq.1)then
1875             sj1=psjti(sngl(qq),qqcut,sngl(s2),0,iqc(3-jj),0)
1876             if(iqc(jj).ne.0)then                             !f1 - f2
1877               sj2=psjti(sngl(qq),qqcut,sngl(s2),iqc(jj),iqc(3-jj),0)
1878             elseif(iqc(3-jj).eq.0)then                       !f1 - g
1879               sj2=psjti(sngl(qq),qqcut,sngl(s2),1,0,0)
1880             else                                             !g  - f2
1881               sj2=psjti(sngl(qq),qqcut,sngl(s2),1,1,0)/naflav/2.      !q  - q
1882      *        +psjti(sngl(qq),qqcut,sngl(s2),-1,1,0)/naflav/2.        !q~ - q
1883      *        +psjti(sngl(qq),qqcut,sngl(s2),1,2,0)*(1.-1./naflav)    !q' - q
1884             endif
1885           else
1886             sj1=psjti1(sngl(qq),qmin(3-jj),qqcut,sngl(s2),0,iqc(3-jj),0)
1887             if(iqc(jj).ne.0)then                             !f1 - f2
1888               sj2=psjti1(sngl(qq),qmin(3-jj),qqcut
1889      *                   ,sngl(s2),iqc(jj),iqc(3-jj),0)
1890             elseif(iqc(3-jj).eq.0)then                       !f1 - g
1891               sj2=psjti1(sngl(qq),qmin(3-jj),qqcut,sngl(s2),1,0,0)
1892             else                                             !g  - f2
1893               sj2=psjti1(sngl(qq),qmin(3-jj),qqcut
1894      *                         ,sngl(s2),1,1,0)/naflav/2.         !q - q
1895      *           +psjti1(sngl(qq),qmin(3-jj),qqcut
1896      *                         ,sngl(s2),-1,1,0)/naflav/2.   !q~ - q
1897      *           +psjti1(sngl(qq),qmin(3-jj),qqcut
1898      *                      ,sngl(s2),1,2,0)*(1.-1./naflav)!q' - q
1899             endif
1900           endif
1901         else
1902           sj1=psjci(sngl(qq),sngl(s2),0)
1903           sj2=psjci(sngl(qq),sngl(s2),1)
1904         endif
1905 
1906         !...gb7 - rejection function for x and q**2 simulation
1907         gb7=dble((sj1*psfap(xx,iqc(jj),0)+sj2*psfap(xx,iqc(jj),1))/
1908      *  log(sngl(qt2)/qcdlam))*psuds(sngl(qq),iqc(jj))*qq/gb0
1909 
1910         if(x.le..5d0)then
1911           gb7=gb7*x**(1.-delh)
1912         else
1913           gb7=gb7*(1.d0-x)*2.d0**delh
1914         endif
1915         if(gb7.gt.1..and.ish.ge.2)write(ifmt,*)'gb7,qq,x,qm0,xm0',
1916      *  gb7,qq,x,qm0,xm0
1917         if(ish.ge.7)write(ifch,*)'gb7:',gb7
1918         if(dble(rangen()).gt.gb7)goto 10
1919       else
1920         qmin2=max(qminn,sngl(qq))
1921         s2min2=4.*qmin2
1922         s2=x*si-qq               !new ladder mass
1923         if(s2.lt.dble(s2min2))goto 10  !rejection in case of too low mass
1924 
1925         call pscs(bcos,bsin)
1926         xmm=x*(ept(3)*dble(bcos)+ept(4)*dble(bsin))
1927         discr=xmm**2+qq*(1.d0-x)-x**2*(pt2+dble(q2mass))
1928         if(discr.lt.0.d0)goto 10
1929         qt=xmm+dsqrt(discr)
1930         ep3(3)=sngl(ept(3)-qt*dble(bcos))   !new parton pt-s
1931         ep3(4)=sngl(ept(4)-qt*dble(bsin))
1932         qt2=dble(ep3(3)**2+ep3(4)**2)
1933 
1934         xx=x
1935         sj1=psjti(sngl(qq),qqcut,sngl(s2),0,iqc(3-jj),0)
1936         sj2=0.
1937         !.... gb7 - rejection function for x and q**2 simulation
1938         gb7=dble(sj1*psfap(xx,iqc(jj),0)/
1939      *  log(sngl(qq)/qcdlam))*qq/gb0
1940         gb7=gb7*x**(1.-delh)
1941         if(gb7.gt.1..and.ish.ge.2)write(ifmt,*)'gb7,qq,x,qm0,xm0',
1942      *  gb7,qq,x,qm0,xm0
1943         if(ish.ge.7)write(ifch,*)'gb7:',gb7
1944         if(dble(rangen()).gt.gb7)goto 10
1945       endif
1946 
1947 ! parton type selection iqc -> iqnew (space like) + iq1 (time like)
1948       nqc(2)=0
1949       iqnew=iqc(jj)
1950       if(rangen().lt.sj1/(sj1+sj2))then
1951         if(iqc(jj).eq.0)then    !g -> g + g  (jt=1)
1952           jt=1
1953           jq=int(1.5+rangen())
1954           nqc(1)=ncc(jq,jj)
1955         else                    !q -> g + q  (jt=2)
1956           jt=2
1957           if(iqc(jj).gt.0)then
1958             jq=1
1959           else
1960             jq=2
1961           endif
1962           nqc(1)=0
1963           iqnew=0
1964         endif
1965         iq1=iqc(jj)
1966         jo=ncc(jq,jj)           !parton origin
1967         if(jo.ne.0)then
1968           jo=ioj(jo)
1969         else
1970           jo=ioj(ncc(3-jq,jj))
1971         endif
1972       else
1973         if(iqc(jj).ne.0)then    !q -> q + g  (jt=3)
1974           iq1=0
1975           jt=3
1976           if(iqc(jj).gt.0)then
1977             jq=1
1978           else
1979             jq=2
1980           endif
1981           nqc(1)=ncc(1,jj)
1982         else                    !g -> q + q  (jt=4)
1983           jt=4
1984           jq=int(1.5+rangen())
1985           iq1=int(naflav*rangen()+1.)*(3-2*jq)
1986           iqnew=-iq1
1987           nqc(1)=ncc(jq,jj)
1988         endif
1989         jo=ncc(1,jj)             !parton origin
1990         if(jo.ne.0)then
1991           jo=ioj(jo)
1992         else
1993           jo=ioj(ncc(2,jj))
1994         endif
1995       endif
1996       if(jo.ge.70.and.jo.ne.99)jo=jo/10
1997       if(ish.ge.5)write(ifch,'(a,i3,a4,i3,a4,i3)')' Process :'
1998      *                       ,iqc(jj),' -> ',iqnew,'  + ',iq1
1999 
2000       if(iabs(iqc(jj)).ne.4)then
2001         q2part=0.
2002 c        write(*,*)'sem:',wpt(jj),pt2,q2mass,wpt(3-jj),jj
2003         if(dabs(wpt(3-jj)).gt.1.d-20)then
2004           wplc=wpt(jj)-(pt2+dble(q2mass))/wpt(3-jj)
2005         else
2006           if(ish.gt.1)write(ifmt,*)'Problem with wpt in sem',wpt(3-jj)
2007           wplc=wpt(jj)-(pt2+dble(q2mass))
2008         endif
2009         qp2max=max(0.,sngl(qt2))
2010       else
2011         q2part=q2mass
2012         wplc=wpt(jj)
2013         qp2max=max(0.,sngl(qq))
2014       endif
2015       eprt=max(dsqrt(qt2+dble(q2part)),.5d0*((1.d0-x)*wplc+
2016      *(qt2+dble(q2part))/(1.d0-x)/wplc))
2017       pl=((1.d0-x)*wplc-eprt)*dble(3-2*jj)
2018       zeta=sqrt(qp2max/si)/sqrt(x*(1.-x))
2019       if(iq1.eq.0)then
2020         iq2ini=9
2021         jo=iq2ini
2022         if(zeta.gt.zetacut)jo=-jo
2023       else
2024         iq2ini=iq1
2025         jo=iq2ini
2026       endif
2027       if(ish.ge.5)write(ifch,*)'qq,eprt,iq2ini,jo,E2-Q2',qp2max,eprt
2028      *,iq2ini,jo,eprt**2-qt2
2029       ntest=0
2030 11    ntest=ntest+1
2031       call timsh1(qp2max,sngl(eprt),iq2ini,jo)
2032       amprt=pprt(5,1)**2
2033       plprt=max(0.d0,eprt**2-qt2)-dble(amprt)
2034       if(plprt.lt.0.d0)then
2035         if(ntest.lt.10000)then
2036           goto 11
2037         else
2038           iret=1
2039           goto 16
2040         endif
2041       endif
2042       ep3(1)=sngl(eprt)
2043       ep3(2)=sngl(dsqrt(plprt))
2044       if(pl.lt.0.d0)ep3(2)=-ep3(2)
2045       ey(1)=1.
2046       ey(2)=1.
2047       ey(3)=1.
2048       do i=1,4
2049         ept1(i)=ept(i)-dble(ep3(i))
2050       enddo
2051       call psdefrot(ep3,s0xh,c0xh,s0h,c0h)
2052       s2new=psnorm(ept1)
2053 
2054       if(iabs(iqc(jj)).ne.4.and.s2new-q2mass.gt.s2min2.or.
2055      *iabs(iqc(jj)).eq.4.and.s2new.gt.s2min2)then
2056         if(iabs(iqc(1)).ne.4.or.jj.eq.1)then
2057           if(jfirst.eq.1)then
2058             gb=dble(psjti(sngl(qq),qqcut,s2new,iqnew,iqc(3-jj),0))
2059           else
2060             gb=dble(psjti1(sngl(qq),qmin(3-jj),qqcut,s2new,iqnew,
2061      *                                               iqc(3-jj),0))
2062           endif
2063         else
2064           gb=dble(psjci(sngl(qq),s2new-q2mass,iqnew))
2065         endif
2066         if(iqnew.eq.0)then
2067           gb=gb/dble(sj1)
2068         else
2069           gb=gb/dble(sj2)
2070         endif
2071         if(dble(rangen()).gt.gb)goto 10
2072       else
2073         goto 10
2074       endif
2075 
2076       if(ish.ge.6)write(ifch,*)'jt,jj,jq,nqc:',jt,jj,jq,nqc
2077       nprtjx=0
2078       call psreti(nqc,jq,1,ey,s0xh,c0xh,s0h,c0h)
2079 
2080       if(jt.eq.1)then
2081         ncc(jq,jj)=nqc(2)
2082       elseif(jt.eq.2)then
2083         ncc(jq,jj)=ncc(1,jj)
2084         ncc(3-jq,jj)=nqc(1)
2085       elseif(jt.eq.3)then
2086         ncc(1,jj)=nqc(2)
2087       elseif(jt.eq.4)then
2088         ncc(1,jj)=ncc(3-jq,jj)
2089       endif
2090       iqc(jj)=iqnew
2091 
2092       do i=1,4
2093         ept(i)=ept1(i)
2094       enddo
2095       ! c.m. energy squared, minimal  4-momentum transfer square and gluon 4-v
2096       ! for the next ladder run
2097       qmin(jj)=sngl(qq)
2098       qminn=qmin2
2099       si=dble(s2new)
2100       nemis(jj)=nemis(jj)+1
2101 
2102       goto 9  !  ---------------next ladder rung ------------>>>>>>>>>>>
2103 
2104 
2105  12   continue !------------------- Born process------------------------
2106 
2107       if(ish.ge.4)write(ifch,*)'Born process'
2108       if(ish.ge.6)write(ifch,*)'iqc,si:',iqc,si
2109       qq=dble(qminn)
2110 
2111    !kkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkk
2112       xpprbor(ncolp,kcol)=(ept(1)+ept(2))/dble(engy)
2113       xmprbor(ncolp,kcol)=(ept(1)-ept(2))/dble(engy)
2114       nemispr(1,ncolp,kcol)=nemis(1)
2115       nemispr(2,ncolp,kcol)=nemis(2)
2116    !   write(*,'(a,2f8.3,i3,3x,2i3)')'------------'
2117    !  *          ,xpprbor(ncolp,kcol),xmprbor(ncolp,kcol),nj-nj00,nemis
2118    !kkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkk
2119 
2120       call psabor(sngl(si),sngl(qq),iqc,ncc,ept,0,iptl,bx)
2121 
2122    !kkkkkkkkkkkkk out Born partons without timelike cascade
2123       if(nprtjx.eq.2)then
2124        do ii=1,2
2125         ptprboo(ii,ncolp,kcol)=sqrt(pprtx(1,ii)**2+pprtx(2,ii)**2)
2126        enddo
2127       else
2128         stop'psahot: should not happen!!!!                '
2129       endif
2130    !kkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkk
2131    !     ptxxx=max(  ptprboo(1,ncolp,kcol) , ptprboo(2,ncolp,kcol)  )
2132    !            print*,sqrt(wm0*wp0),sqrt(wmi*wpi),ptxxx  !++++++++++++++++++++++++++
2133    !      if(ptxxx.lt.10)goto1
2134    !            print*,'  ++++++++++++++++++++++++++++++++++++ '
2135    !kkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkk
2136       if(nj.ne.0.)then
2137         do i=1,nj
2138           do l=1,6
2139             bxj(l,i)=bx(l)
2140           enddo
2141           ityj(i)=ity
2142           iorj(i)=iptl
2143         enddo
2144       endif
2145 
2146       call psjarr(jfl)       !kinky strings formation
2147 
2148       if(jfl.eq.0.and.ish.ge.4)write(ifch,*)
2149      *'jfl,nj,nptl',jfl,nj,nptl
2150       if(jfl.eq.0)goto 1
2151 
2152 c --- update remnant flavour ---
2153 
2154       iret1=0
2155       call idenco(jcp,icp,iret1)
2156       if(iret1.eq.1)call utstop('Problem with proj rem in psahot !&')
2157       iret2=0
2158       call idenco(jct,ict,iret2)
2159       if(iret2.eq.1)call utstop('Problem with targ rem in psahot !&')
2160 
2161       do i=1,2
2162         icproj(i,ip)=icp(i)
2163         ictarg(i,it)=ict(i)
2164       enddo
2165 
2166       if(ish.ge.3)write(ifch,*)'End psahot (icp,ict):',icp,ict
2167 
2168       if(iremn.ge.2)then        !uses precalculated flavors
2169 
2170         do j=1,2
2171           do n=1,nrflav
2172             jcpref(n,j,ip)=jcpr(n,j)
2173             jctref(n,j,it)=jctr(n,j)
2174           enddo
2175         enddo
2176         if(ish.ge.3)then
2177           write(ifch,'(a,6i3,2x,6i3)')' proj:  ',jcpr
2178           write(ifch,'(a,6i3,2x,6i3)')' targ:  ',jctr
2179         endif
2180 
2181       endif
2182 
2183 
2184 c ------------------------------
2185 
2186 16    continue
2187       call utprix('psahot',ish,ishini,3)
2188 
2189       q2fin=q2finsave
2190       return
2191       end
2192 
2193 c------------------------------------------------------------------------
2194       subroutine psjarr(jfl)
2195 c-----------------------------------------------------------------------
2196 c
2197 c   final jets rearrangement according to their colour connection
2198 c   and write to /cptl/
2199 c
2200 c jfl - flag for the rejection (in case of jfl=0)
2201 c-----------------------------------------------------------------------
2202 c Input:
2203 
2204       parameter (mjstr=20000)
2205       common /psar29/ eqj(4,mjstr),iqj(mjstr),ncj(2,mjstr),ioj(mjstr),nj
2206 
2207 c eqj(1:4,k) - 4-momentum (qgs) for k-th parton;
2208 c bxj(1:4,k) - coordinates for k-th parton formation point;
2209 c iqj(k) - ID (qgs) for k-th parton;
2210 c ncj(1:2,k) - colour connected partons indexes for k-th parton;
2211 c nj - number of partons
2212 c-----------------------------------------------------------------------
2213       dimension mark(mjstr)
2214       double precision  ept(4),eptot(2)
2215       include 'epos.inc'
2216       include 'epos.incsem'
2217 
2218       if(nj.eq.0)then
2219         jfl=1
2220         return
2221       endif
2222 csp : TEST if string are full  
2223      
2224       do k=1,nj      
2225         if(iqj(k).eq.0)then   !gluon must have two neighbours        
2226           if(ncj(1,k).eq.0)then !first neigbour missing
2227 csp          write(*,*)'correction'
2228             do kk=1,nj          !look to which parton he is connected
2229               if(ncj(1,kk).eq.k)then
2230                 if(ncj(2,k).ne.kk)ncj(1,k)=kk !if not already connected : connection
2231               elseif(ncj(2,kk).eq.k)then
2232                 if(ncj(1,k).ne.kk)ncj(1,k)=kk
2233                 endif
2234             enddo  
2235            endif
2236            if(ncj(2,k).eq.0)then !second neigbour missing
2237 csp           write(*,*)'correction'
2238             do kk=1,nj
2239               if(ncj(2,kk).eq.k)then
2240                 if(ncj(1,k).ne.kk)ncj(2,k)=kk
2241                 elseif(ncj(1,kk).eq.k)then
2242                 if(ncj(2,k).ne.kk)ncj(2,k)=kk
2243                 endif
2244             enddo  
2245            endif          
2246         endif
2247       enddo
2248       
2249       
2250 csp END OF TEST
2251 c      do k=1,nj  !???????????????
2252 c       eqj(1,k)=dsqrt(0d0+eqj(2,k)**2+eqj(3,k)**2+eqj(4,k)**2)
2253 c      enddo
2254       if(ish.ge.3)then
2255         write (ifch,*)'psjarr: nj',nj
2256         do k=1,nj
2257           if(iqj(k).ne.0)ncj(2,k)=0
2258           write(ifch,'(a,i4)')' parton',k
2259           write(ifch,'(i6,2x,4e10.3,2x,2i3)')iqj(k)
2260      *    ,(eqj(j,k),j=1,4),(ncj(j,k),j=1,2)
2261         enddo
2262       endif
2263 
2264       jfl=0
2265       do i=1,nj
2266         mark(i)=1
2267       enddo
2268       nptl0=nptl
2269 c total energy of the two half of the Pomeron
2270       eptot(1)=0d0
2271       eptot(2)=0d0
2272 
2273 1     continue
2274       do ij=1,nj
2275         if(mark(ij).ne.0.and.iqj(ij).ne.0)goto 2
2276       enddo
2277 2     continue
2278       jfirst=1
2279 c to calculate the total energy of the 2 big strings produce by a pomeron
2280 c we first fix to which string (1 or 2) the sub-string belong to
2281       ij0=ncj(1,ij)
2282       kkk=1
2283       if(ij0.gt.0)then
2284         if(ncj(1,ij0).eq.ij)kkk=1
2285         if(ncj(2,ij0).eq.ij)kkk=2
2286       endif
2287       if(iabs(iqj(ij)).le.2)then
2288         am1=amhadr(1)
2289       elseif(iabs(iqj(ij)).eq.4)then
2290         am1=amhadr(3)
2291       elseif(iabs(iqj(ij)).eq.40)then
2292         am1=qcmass
2293       else
2294         am1=amhadr(2)
2295       endif
2296 
2297       do i=1,4
2298         ept(i)=0.
2299       enddo
2300 
2301 3     continue
2302       call pspawr(ij,kkk)
2303 
2304       mark(ij)=0
2305 
2306       do i=1,4
2307         ept(i)=ept(i)+eqj(i,ij)
2308       enddo
2309       eptot(kkk)=eptot(kkk)+eqj(1,ij)
2310 
2311       if(iqj(ij).ne.0)then
2312         if(jfirst.ne.1)then
2313           if(iabs(iqj(ij)).le.2)then
2314             am2=amhadr(1)
2315           elseif(iabs(iqj(ij)).eq.4)then
2316             am2=amhadr(3)
2317           elseif(iabs(iqj(ij)).eq.40)then
2318             am2=qcmass
2319           else
2320             am2=amhadr(2)
2321           endif
2322           amj=(am1+am2+stmass)**2
2323           sm=psnorm(ept)
2324           if(sm.lt.amj)then
2325             nptl=nptl0
2326             goto 999
2327           endif
2328 
2329           if(nptl-nptl0.lt.nj)then
2330             goto 1
2331           else
2332             if(iLHC.ne.1)then
2333 c at the end of the process, save eptot(kkk) is qsqptl of each particle
2334               do k=nptl0+1,nptl
2335                 qsqptl(k)=sngl(eptot(nint(qsqptl(k)))**2)
2336               enddo
2337             endif
2338             
2339             if(ish.ge.3)then
2340               write (ifch,*)'psjarr: nptl',nptl
2341               do k=nptl0+1,nptl
2342                 write(ifch,'(a,i4)')' particle',k
2343                 write(ifch,'(i5,2x,6e10.3)')idptl(k)
2344      *          ,(pptl(j,k),j=1,5),sqrt(qsqptl(k))
2345               enddo
2346             endif
2347             jfl=1
2348             goto 999
2349           endif
2350         else
2351           jfirst=0
2352           njpar=ij
2353           ij=ncj(1,ij)
2354           if(ij.eq.0)write(ifch,*)
2355      &'IN PSJARR ij=0 :parton,',njpar,'with no connection'
2356           goto 3
2357         endif
2358 
2359       else
2360         if(ncj(1,ij).eq.njpar)then
2361           njdau=ncj(2,ij)
2362         else
2363           njdau=ncj(1,ij)
2364         endif
2365         njpar=ij
2366         ij=njdau
2367         if(ij.eq.0)write(ifch,*)
2368      &'IN PSJARR ij=0 :parton,',njpar, 'with no connection'
2369         goto 3
2370       endif
2371       
2372       
2373   999 continue
2374 
2375       end
2376 
2377 c------------------------------------------------------------------------
2378       subroutine pspawr(kj,kkk)
2379 c-----------------------------------------------------------------------
2380 c pspawr - writing final parton kj into particle list
2381 c------------------------------------------------------------------------
2382 c Input:
2383       parameter (mjstr=20000)
2384       common /psar29/ eqj(4,mjstr),iqj(mjstr),ncj(2,mjstr),ioj(mjstr),nj
2385       common /psar30/ iorj(mjstr),ityj(mjstr),bxj(6,mjstr),q2j(mjstr)
2386 
2387 c eqj(1:4,k) - 4-momentum (qgs) for k-th parton;
2388 c bxj(1:4,k) - coordinates for k-th parton formation point;
2389 c iqj(k) - ID (qgs) for k-th parton;
2390 c ncj(1:2,k) - colour connected partons indexes for k-th parton;
2391 c nj - number of partons
2392 c kkk is the indice of the original half of the pomeron
2393 c-----------------------------------------------------------------------
2394       include 'epos.inc'
2395       include 'epos.incsem'
2396       common/ciptl/iptl
2397 
2398       nptl=nptl+1
2399       if(ish.ge.9)write (ifch,*)'nptl,kj (sto)',nptl,kj
2400       if(nptl.ge.mxptl.or.kj.le.0)then
2401        write (ifmt,*)'nptl,kj',nptl,kj
2402        call alist('Error in pspawr: nptl or kj out of bounds &',1,nptl)
2403        call utstop('nptl or kj out of bounds&')
2404       endif
2405 
2406       if(ifrptl(1,iptl).eq.0)ifrptl(1,iptl)=nptl
2407       ifrptl(2,iptl)=nptl
2408 
2409       pptl(1,nptl)=eqj(3,kj)
2410       pptl(2,nptl)=eqj(4,kj)
2411       pptl(3,nptl)=eqj(2,kj)
2412       pptl(4,nptl)=eqj(1,kj)
2413       pptl(5,nptl)=0.
2414       idptl(nptl)=psidd(iqj(kj))
2415       iorptl(nptl)=iorj(kj)
2416       jorptl(nptl)=ioj(kj)
2417       istptl(nptl)=20
2418       do i=1,4
2419         xorptl(i,nptl)=bxj(i,kj)
2420       enddo
2421       tivptl(1,nptl)=bxj(5,kj)
2422       tivptl(2,nptl)=bxj(6,kj)
2423       ityptl(nptl)=ityj(kj)
2424 c register to which big string the particle belongs to
2425       if(iLHC.eq.1)then
2426         qsqptl(nptl)=q2j(kj)
2427       else
2428         qsqptl(nptl)=float(kkk)
2429       endif
2430           !kkkkkkkkkkkkkkkkkkkkkkk
2431        !    write(*,'(a,2i4,i6,f8.3)')'.... ',kj,nptl,idptl(nptl)
2432        !*     ,sqrt(pptl(1,nptl)**2+pptl(2,nptl)**2)
2433       return
2434       end
2435 
2436 c------------------------------------------------------------------------
2437       subroutine psabor(si,qq,iqc,ncc,ept,jdis,iptl,coordo)
2438 c------------------------------------------------------------------------
2439 c psabor - highest virtuality subprocess in the ladder
2440 c si - c.m. energy squared for the process,
2441 c qq - p_t-cutoff for the process due to the evolution,
2442 c iqc(i), i=1,2 - incoming parton types(0-g,(+-)1-u(u~),(+-)2-d(d~)etc.),
2443 c ncc(i,j), i,j=1,2 - incoming partons color connections,
2444 c ept(4) - total 4-momentum for the system of the 2 partons
2445 c jdis=0 - hadronic process in pp; 1 - resolved photon process
2446 c------------------------------------------------------------------------
2447       double precision ept(4),psutz,psuds
2448       dimension ep3(4),ey(3),wsub(5),iqc(2),ncc(2,2),nqc(2),coordo(6)
2449       parameter (mjstr=20000)
2450       common /psar29/eqj(4,mjstr),iqj(mjstr),ncj(2,mjstr),ioj(mjstr),nj
2451       parameter (ntim=1000)
2452       common/cprt/pprt(5,ntim),q2prt(ntim),idaprt(2,ntim),idprt(ntim)
2453      &,iorprt(ntim),jorprt(ntim),nprtj
2454       common/cprtx/nprtjx,pprtx(5,2)
2455       include 'epos.incsem'
2456       include 'epos.inc'
2457       call utpri('psabor',ish,ishini,5)
2458 
2459       if(iabs(iqc(1)).ne.4)then   !gluon or light quark
2460         q2mass=0.
2461       else                        !c-quark
2462         q2mass=qcmass**2
2463       endif
2464       p1=si/(1.+q2mass/si)
2465       if(p1.gt.4.*qq)then                 !|t|-cutoff (p^2_t>qq)
2466         tmin=2.*qq/(1.+sqrt(1.-4.*qq/p1))
2467       else
2468         tmin=2.*qq
2469       endif
2470       tmax=.5*p1
2471 
2472       fborn=0.
2473       qt2=tmin*(1.d0-tmin/p1)
2474       if(qt2.lt..999d0*qq.and.ish.ge.2)write(ifmt,*)'qt20,qq',qt2,qq
2475       if(iqc(1).ne.0.or.iqc(2).ne.0)then
2476         do l=1,5       !sum over different subprocesses
2477           wsub(l)=psbori(si,tmin,iqc(1),iqc(2),l)  !matrix element
2478         if(l.le.3)then
2479           wsub(l)=wsub(l)*pssalf(qt2/qcdlam)**2
2480         elseif(l.le.4)then
2481           wsub(l)=wsub(l)*pssalf(qt2/qcdlam)*alfe/2/pi
2482         else
2483           wsub(l)=wsub(l)*(alfe/2/pi)**2
2484         endif
2485           fborn=fborn+wsub(l)
2486         enddo
2487         fborn=tmin**2*fborn
2488       else
2489         do l=1,5
2490           wsub(l)=psbori(si,.5*si,iqc(1),iqc(2),l)
2491         if(l.le.3)then
2492           wsub(l)=wsub(l)*pssalf(qt2/qcdlam)**2
2493         elseif(l.le.4)then
2494           wsub(l)=wsub(l)*pssalf(qt2/qcdlam)*alfe/2/pi
2495         else
2496           wsub(l)=wsub(l)*(alfe/2/pi)**2
2497         endif
2498           fborn=fborn+wsub(l)
2499         enddo
2500         fborn=.25*si**2*fborn
2501       endif
2502        if(jdis.eq.0)then
2503          scale1=qt2
2504        else
2505          scale1=4.*qt2
2506        endif
2507       gb0=dble(fborn)*psuds(scale1,iqc(1))*psuds(qt2,iqc(2))
2508       if(ish.ge.7)write(ifch,*)'tmin,gb0:',tmin,gb0
2509 
2510 c------------------------------------------------
2511 c 4-momentum transfer squared is simulated first as dq_t**2/q_t**4 from
2512 c tmin to s/2
2513 14    q2=tmin/(1.-rangen()*(1.-tmin/tmax))    !q2=min(|t|,|u|)
2514       qt2=q2*(1.-q2/p1)                       !qt2=p_t^2 for the process
2515       if(qt2.lt.qq.and.ish.ge.2)write(ifmt,*)'qt2,qq',qt2,qq
2516       if(rangen().lt..5)then  !|u|=q2, |t|=p1-q2
2517         jm=2                  !first parton to be considered
2518         tq=p1-q2
2519       else                    !|t|=q2, |u|=p1-q2
2520         jm=1                  !first parton to be considered
2521         tq=q2
2522       endif
2523 
2524       fborn=0.
2525       do l=1,5                       !sum over different subprocesses
2526         wsub(l)=psbori(si,tq,iqc(1),iqc(2),l)
2527         if(l.le.3)then
2528           wsub(l)=wsub(l)*pssalf(qt2/qcdlam)**2
2529         elseif(l.le.4)then
2530           wsub(l)=wsub(l)*pssalf(qt2/qcdlam)*alfe/2/pi
2531         else
2532           wsub(l)=wsub(l)*(alfe/2/pi)**2
2533         endif
2534         fborn=fborn+wsub(l)
2535       enddo
2536        if(jdis.eq.0)then
2537          scale1=qt2
2538        else
2539          scale1=4.*qt2
2540        endif
2541       gb=dble(q2**2*fborn)
2542      &*psuds(scale1,iqc(1))*psuds(qt2,iqc(2))/gb0 !rejection function
2543       if(ish.ge.7)write(ifch,*)'q2,qt2,gb:',q2,qt2,gb
2544 
2545       if(dble(rangen()).gt.gb)goto 14                   !rejection
2546 
2547 c determination of the color configuration
2548       nqc(2)=0
2549       if(iqc(1).eq.0.and.iqc(2).eq.0)then      !g+g
2550         jq=int(1.5+rangen())    !jq=1(2) - transfer of color (anticolor)
2551         nqc(1)=ncc(jq,jm)
2552 
2553         if(rangen().lt.wsub(1)/fborn)then      !gg->gg
2554           if(rangen().lt..5)then
2555             jt=1                !anticolor-color annihilation
2556             nqc(2)=0
2557             njc1=ncc(3-jq,jm)
2558             njc2=ncc(jq,3-jm)
2559             if(iqj(njc1).ne.0)then
2560               ncj(1,njc1)=njc2
2561             else
2562               ncj(jq,njc1)=njc2
2563             endif
2564             if(iqj(njc2).ne.0)then
2565               ncj(1,njc2)=njc1
2566             else
2567               ncj(3-jq,njc2)=njc1
2568             endif
2569           else
2570             jt=2                    !produced gluons get color and
2571             nqc(2)=ncc(3-jq,3-jm)   !anticolor from the 2 parents
2572           endif
2573         else                                   !gg->qq~
2574           jt=9                  !anticolor-color annihilation
2575           iqc(jm)=int(naflav*rangen()+1)*(3-2*jq) !(anti)quark flavor
2576           iqc(3-jm)=-iqc(jm)
2577           njc1=ncc(3-jq,jm)
2578           njc2=ncc(jq,3-jm)
2579           if(iqj(njc1).ne.0)then
2580             ncj(1,njc1)=njc2
2581           else
2582             ncj(jq,njc1)=njc2
2583           endif
2584           if(iqj(njc2).ne.0)then
2585             ncj(1,njc2)=njc1
2586           else
2587             ncj(3-jq,njc2)=njc1
2588           endif
2589         endif
2590 
2591       elseif(iqc(1)*iqc(2).eq.0)then       !q(q~)+g
2592         if(iqc(1)+iqc(2).gt.0)then
2593           jq=1                             !q
2594         else
2595           jq=2                             !q~
2596         endif
2597        if(rangen().lt.wsub(1)/fborn)then  !q(q~)g->q(q~)g
2598         if(rangen().lt..5)then      !anticolor-color annihilation
2599           if(iqc(jm).eq.0)then
2600             jt=3                    !first parton=g
2601             nqc(1)=ncc(jq,jm)
2602             njc1=ncc(3-jq,jm)
2603             njc2=ncc(1,3-jm)
2604             if(iqj(njc1).ne.0)then
2605               ncj(1,njc1)=njc2
2606             else
2607               ncj(jq,njc1)=njc2
2608             endif
2609             if(iqj(njc2).ne.0)then
2610               ncj(1,njc2)=njc1
2611             else
2612               ncj(3-jq,njc2)=njc1
2613             endif
2614 
2615           else
2616             jt=4                    !first parton=q(q~)
2617             nqc(1)=0
2618             njc1=ncc(1,jm)
2619             njc2=ncc(3-jq,3-jm)
2620             if(iqj(njc1).ne.0)then
2621               ncj(1,njc1)=njc2
2622             else
2623               ncj(3-jq,njc1)=njc2
2624             endif
2625             if(iqj(njc2).ne.0)then
2626               ncj(1,njc2)=njc1
2627             else
2628               ncj(jq,njc2)=njc1
2629             endif
2630           endif
2631 
2632         else                        !color transfer
2633           if(iqc(jm).eq.0)then
2634             jt=5                    !first parton=g
2635             nqc(2)=ncc(3-jq,jm)
2636             nqc(1)=ncc(1,3-jm)
2637           else                      !first parton=q(q~)
2638             jt=6
2639             nqc(1)=ncc(jq,3-jm)
2640           endif
2641         endif
2642 
2643        else          !q(q~)g->q(q~)-gamma (+-color annihilation)
2644           if(iqc(jm).eq.0)then
2645             jt=11                    !first parton=g
2646             nqc(1)=ncc(jq,jm)
2647             njc1=ncc(3-jq,jm)
2648             njc2=ncc(1,3-jm)
2649             if(iqj(njc1).ne.0)then
2650               ncj(1,njc1)=njc2
2651             else
2652               ncj(jq,njc1)=njc2
2653             endif
2654             if(iqj(njc2).ne.0)then
2655               ncj(1,njc2)=njc1
2656             else
2657               ncj(3-jq,njc2)=njc1
2658             endif
2659             iqc(jm)=iqc(3-jm)
2660             iqc(3-jm)=10              !make the second output is gamma.
2661 
2662           else
2663             jt=12                    !first parton=q(q~)
2664             nqc(1)=ncc(jq,3-jm)                 !here nqc(1) is gluon.
2665             njc1=ncc(1,jm)
2666             njc2=ncc(3-jq,3-jm)
2667             if(iqj(njc1).ne.0)then
2668               ncj(1,njc1)=njc2
2669             else
2670               ncj(3-jq,njc1)=njc2
2671             endif
2672             if(iqj(njc2).ne.0)then
2673               ncj(1,njc2)=njc1
2674             else
2675               ncj(jq,njc2)=njc1
2676             endif
2677             iqc(3-jm)=10
2678           endif
2679        endif
2680 
2681        elseif(iqc(1)*iqc(2).gt.0)then
2682         jt=7                        !qq->qq (q~q~->q~q~)
2683         if(iqc(1).gt.0)then
2684           jq=1
2685         else
2686           jq=2
2687         endif
2688         nqc(1)=ncc(1,3-jm)
2689 
2690       else                          ! qq~ ->
2691         if(iqc(jm).gt.0)then
2692           jq=1
2693         else
2694           jq=2
2695         endif
2696         aks=rangen()
2697         if(aks.lt.(wsub(1)+wsub(2))/fborn)then
2698           jt=8                     ! qq~->qq~ (anticolor-color annihilation)
2699           if(aks.gt.wsub(1)/fborn)then
2700             iqa=iabs(iqc(jm))
2701             iq=int((naflav-1)*rangen())+1
2702             if(iq.eq.iqa)iq=naflav
2703             iqc(jm)=iq*iqc(jm)/iqa
2704             iqc(3-jm)=-iqc(jm)
2705           endif
2706           nqc(1)=0
2707           njc1=ncc(1,jm)
2708           njc2=ncc(1,3-jm)
2709           if(iqj(njc1).ne.0)then
2710             ncj(1,njc1)=njc2
2711           else
2712             ncj(3-jq,njc1)=njc2
2713           endif
2714           if(iqj(njc2).ne.0)then
2715             ncj(1,njc2)=njc1
2716           else
2717             ncj(jq,njc2)=njc1
2718           endif
2719         elseif(aks.lt.(wsub(1)+wsub(2)+wsub(3))/fborn)then
2720           jt=10                    !color transfer  qq~->gg
2721           iqc(1)=0
2722           iqc(2)=0
2723           nqc(1)=ncc(1,jm)
2724           nqc(2)=0
2725         elseif(aks.lt.(wsub(1)+wsub(2)+wsub(3)+wsub(4))/fborn)then
2726           jt=13                   ! qq~->g+gamma
2727           nqc(1)=ncc(1,jm)
2728           nqc(2)=ncc(1,3-jm)
2729           iqc(jm)=0
2730           iqc(3-jm)=10
2731         else
2732           jt=14                  ! qq~->gamma+gamma
2733           njc1=ncc(1,jm)
2734           njc2=ncc(1,3-jm)
2735           if(iqj(njc1).ne.0)then
2736             ncj(jq,njc1)=njc2
2737           else
2738             ncj(3-jq,njc1)=njc2
2739           endif
2740           if(iqj(njc2).ne.0)then
2741             ncj(3-jq,njc2)=njc1
2742           else
2743             ncj(jq,njc2)=njc1
2744           endif
2745           iqc(jm)=10
2746           iqc(3-jm)=10
2747         endif
2748       endif
2749 
2750       if(jt.ne.8.and.jt.ne.9)then
2751         jq2=jq
2752       else
2753         jq2=3-jq
2754       endif
2755 
2756       call psdeftr(si+q2mass,ept,ey)    !lorentz boost to c.m. frame
2757 
2758       qt=sqrt(qt2)                      !p_t
2759       call pscs(bcos,bsin)              !cos and sin of the polar angle
2760       if(iabs(iqc(1)).ne.4)then
2761 clight cone momentum share for the first parton
2762         z=sngl(psutz(dble(si),dble(qt2),dble(qt2)))
2763         if((jt.eq.11.and.jm.eq.1).or.(jt.eq.12.and.jm.eq.2)
2764      $   .or.(jt.eq.13.and.jm.eq.2))z=1-z
2765         wp3=z*sqrt(si)
2766         wm3=qt2/wp3
2767       elseif(jm.eq.1)then
2768         z=sngl(psutz(dble(si),dble(qt2+q2mass),dble(qt2)))
2769         wp3=z*sqrt(si)
2770         wm3=(qt2+q2mass)/wp3
2771       else
2772         z=sngl(psutz(dble(si),dble(qt2),dble(qt2+dble(q2mass))))
2773         wp3=z*sqrt(si)
2774         wm3=qt2/wp3
2775       endif
2776       ep3(1)=.5*(wp3+wm3)               !parton 4-momentum
2777       ep3(2)=.5*(wp3-wm3)
2778       ep3(3)=qt*bcos
2779       ep3(4)=qt*bsin
2780       call psdefrot(ep3,s0xh,c0xh,s0h,c0h)   !spacial rotation to z-axis
2781 
2782       zeta=2.                        !2=back-to-back emission (angle=pi)
2783       if(iqc(jm).eq.0)then
2784         iq2ini1=9
2785         jo1=iq2ini1
2786         if(zeta.gt.zetacut)jo1=-jo1
2787 c        q2fin=q2fin+zoeinc
2788       else
2789         iq2ini1=iqc(jm)
2790         jo1=iq2ini1
2791       endif
2792       if(iqc(3-jm).eq.0)then
2793         iq2ini2=9
2794         jo2=iq2ini2
2795         if(zeta.gt.zetacut)jo2=-jo2
2796 c        q2fin=q2fin+zoeinc
2797       else
2798         iq2ini2=iqc(3-jm)
2799         jo2=iq2ini2
2800       endif
2801       if(jt.le.10)then
2802         qq1=qt2
2803         qq2=qt2
2804       elseif(jt.le.13)then
2805         qq1=qt2
2806         qq2=0
2807       else
2808         qq1=0
2809         qq2=0
2810       endif
2811 
2812       call timsh2(qq1,qq2,sqrt(si),iq2ini1,iq2ini2,jo1,jo2)  !final state cascade
2813 
2814       if(jt.le.10)then      !color connection for the 2nd parton
2815         if(ish.ge.6)write(ifch,*)'jt,jq,nqc:',jt,jq,nqc
2816         call psreti(nqc,jq,1,ey,s0xh,c0xh,s0h,c0h) !color conn. reconstruction
2817         if(jt.eq.1)then
2818           nqc(1)=nqc(2)
2819           nqc(2)=ncc(3-jq,3-jm)
2820         elseif(jt.eq.2)then
2821           nqc(2)=ncc(3-jq,jm)
2822           nqc(1)=ncc(jq,3-jm)
2823         elseif(jt.eq.3)then
2824           nqc(1)=nqc(2)
2825         elseif(jt.eq.4)then
2826           nqc(2)=nqc(1)
2827           nqc(1)=ncc(jq,3-jm)
2828         elseif(jt.eq.5)then
2829           nqc(1)=ncc(jq,jm)
2830         elseif(jt.eq.6)then
2831           nqc(2)=ncc(3-jq,3-jm)
2832           nqc(1)=ncc(1,jm)
2833         elseif(jt.eq.7)then
2834           nqc(1)=ncc(1,jm)
2835         elseif(jt.eq.9)then
2836           nqc(1)=ncc(3-jq,3-jm)
2837         elseif(jt.eq.10)then
2838           nqc(1)=nqc(2)
2839           nqc(2)=ncc(1,3-jm)
2840         endif
2841         if(ish.ge.6)write(ifch,*)'jt,jq2,nqc:',jt,jq2,nqc
2842         call psreti(nqc,jq2,2,ey,s0xh,c0xh,s0h,c0h) !color conn. reconstr.
2843       elseif(jt.le.13)then
2844         if(ish.ge.6)write(ifch,*)'jt,jq,nqc:',jt,jq,nqc
2845         call psreti(nqc,jq,1,ey,s0xh,c0xh,s0h,c0h) !color conn. reconstruction
2846         ep3(1)=pprt(4,2)
2847         ep3(2)=pprt(3,2)
2848         ep3(3)=pprt(1,2)
2849         ep3(4)=pprt(2,2)
2850         call psrotat(ep3,s0xh,c0xh,s0h,c0h)  !special rotation for photon.
2851         call pstrans(ep3,ey,1)
2852         nptl=nptl+1
2853         pptl(1,nptl)=ep3(3)
2854         pptl(2,nptl)=ep3(4)
2855         pptl(3,nptl)=ep3(2)
2856         pptl(4,nptl)=ep3(1)
2857         pptl(5,nptl)=0
2858         idptl(nptl)=10
2859         iorptl(nptl)=iptl
2860         istptl(nptl)=0
2861         jorptl(nptl)=0
2862         do i=1,4
2863           xorptl(i,nptl)=coordo(i)
2864         enddo
2865         tivptl(1,nptl)=coordo(5)
2866         tivptl(2,nptl)=coordo(6)
2867         ityptl(nptl)=71
2868         ifrptl(1,nptl)=0
2869         ifrptl(2,nptl)=0
2870       else
2871         if(ish.ge.6)write(ifch,*)'jt,iqc:',jt,iqc
2872         do j=1,2
2873           ep3(1)=pprt(4,j)
2874           ep3(2)=pprt(3,j)
2875           ep3(3)=pprt(1,j)
2876           ep3(4)=pprt(2,j)
2877           call psrotat(ep3,s0xh,c0xh,s0h,c0h)  !special rotation for photon.
2878           call pstrans(ep3,ey,1)
2879           nptl=nptl+1
2880           pptl(1,nptl)=ep3(3)
2881           pptl(2,nptl)=ep3(4)
2882           pptl(3,nptl)=ep3(2)
2883           pptl(4,nptl)=ep3(1)
2884           pptl(5,nptl)=0
2885           idptl(nptl)=10
2886           iorptl(nptl)=iptl
2887           istptl(nptl)=0
2888           jorptl(nptl)=0
2889           do i=1,4
2890             xorptl(i,nptl)=coordo(i)
2891           enddo
2892           tivptl(1,nptl)=coordo(5)
2893           tivptl(2,nptl)=coordo(6)
2894           ityptl(nptl)=72
2895           ifrptl(1,nptl)=0
2896           ifrptl(2,nptl)=0
2897         enddo
2898       endif
2899       call utprix('psabor',ish,ishini,5)
2900       return
2901       end
2902 
2903 c------------------------------------------------------------------------
2904       subroutine psreti(nqc,jort,nfprt,ey,s0xh,c0xh,s0h,c0h)
2905 c-----------------------------------------------------------------------
2906 c jet reconstructuring procedure - 4-momenta for all final jets
2907 c nqc(i) - colour connections for the jet
2908 c jort - color orientation for gluons (=1 if +color goes first, =-1 otherwise)
2909 
2910       parameter (ntim=1000)
2911       common/cprt/pprt(5,ntim),q2prt(ntim),idaprt(2,ntim),idprt(ntim)
2912      &,iorprt(ntim),jorprt(ntim),nprtj
2913 c nprtj - number of partons in the jet (including virtual ones)
2914 c pprt - 5-momenta for the partons
2915 c idprt - parton id
2916 c iorprt - parent parton position in the list
2917 c idaprt - daughter partons positions in the list
2918 
2919 c output:
2920       parameter (mjstr=20000)
2921       common /psar29/ eqj(4,mjstr),iqj(mjstr),ncj(2,mjstr),ioj(mjstr),nj
2922 c nj - number of final jets
2923 c eqj(i,j) - 4-momentum for the final jet j
2924 c iqj(j) - flavour for the final jet j
2925 c ncj(m,j) - colour connections for the final jet j
2926       common /psar30/ iorj(mjstr),ityj(mjstr),bxj(6,mjstr),q2j(mjstr)
2927 c-----------------------------------------------------------------------
2928       dimension ep3(4),nqc(2),ncc(2,ntim),ey(3)
2929       include 'epos.inc'
2930       include 'epos.incsem'
2931       common/cprtx/nprtjx,pprtx(5,2)
2932 
2933       if(ish.ge.6)then
2934         write (ifch,*)'nprtj',nprtj
2935         do i=1,nprtj
2936           write (ifch,*)'i,ic,np,ndd',i,idprt(i),iorprt(i),
2937      *    idaprt(1,i),idaprt(2,i)
2938         enddo
2939       endif
2940 
2941       ncc(1,nfprt)=nqc(1)
2942       if(idprt(nfprt).eq.9)ncc(2,nfprt)=nqc(2)
2943       iprt=nfprt
2944 
2945       if(nprtjx.eq.2)then !out Born before timelike cascade
2946        ep3(1)=pprtx(4,iprt)
2947        ep3(2)=pprtx(3,iprt)
2948        ep3(3)=pprtx(1,iprt)
2949        ep3(4)=pprtx(2,iprt)
2950        call psrotat(ep3,s0xh,c0xh,s0h,c0h)
2951        call pstrans(ep3,ey,1)
2952        pprtx(4,iprt)=ep3(1)
2953        pprtx(3,iprt)=ep3(2)
2954        pprtx(1,iprt)=ep3(3)
2955        pprtx(2,iprt)=ep3(4)
2956       endif
2957 
2958 1     continue
2959 
2960       idau1=idaprt(1,iprt)
2961       idau2=idaprt(2,iprt)
2962       icp=idprt(iprt)
2963 
2964       if(ish.ge.6)then
2965         write (ifch,*)'1-iprt,icp,idau1,idau2',iprt,icp,idau1,idau2,
2966      *  ncc(1,iprt)
2967         if(icp.eq.9)write (ifch,*)'ncc2',ncc(2,iprt)
2968       endif
2969 
2970       if(idau1.ne.0.)then         !virtual parton
2971         icd1=idprt(idau1)
2972 
2973         if(icp.eq.9)then
2974           if(icd1.ne.9)then      !g -> qq~
2975             ncc(1,idau1)=ncc(jort,iprt)
2976             ncc(1,idau2)=ncc(3-jort,iprt)
2977           else                    !g -> gg
2978             ncc(1,idau1)=ncc(1,iprt)
2979             ncc(2,idau1)=0
2980             ncc(2,idau2)=ncc(2,iprt)
2981             ncc(1,idau2)=0
2982           endif
2983         else                      !q -> qg
2984           ncc(1,idau1)=0
2985           if(icp*(3-2*jort).gt.0)then
2986             ncc(1,idau2)=ncc(1,iprt)
2987             ncc(2,idau2)=0
2988           else
2989             ncc(1,idau2)=0
2990             ncc(2,idau2)=ncc(1,iprt)
2991           endif
2992         endif
2993         iprt=idau1
2994         goto 1
2995       else
2996 
2997         nj=nj+1
2998         ep3(1)=pprt(4,iprt)
2999         ep3(2)=pprt(3,iprt)
3000         ep3(3)=pprt(1,iprt)
3001         ep3(4)=pprt(2,iprt)
3002         call psrotat(ep3,s0xh,c0xh,s0h,c0h)
3003         call pstrans(ep3,ey,1)
3004         do i=1,4
3005           eqj(i,nj)=ep3(i)
3006         enddo
3007 
3008         if(icp.eq.9)then
3009           iqj(nj)=0
3010         elseif(iabs(icp).lt.3)then
3011           iqj(nj)=icp
3012         elseif(iabs(icp).eq.3)then
3013           iqj(nj)=icp*4/3
3014         else
3015           iqj(nj)=icp*10
3016         endif
3017 
3018         ioj(nj)=jorprt(iprt) !flavor of mother parton
3019         q2j(nj)=q2prt(iprt)
3020 
3021         if(iqj(nj).ne.0)then
3022           njc=ncc(1,iprt)
3023           if(njc.ne.0)then
3024             ncj(1,nj)=njc
3025             iqc=iqj(njc)
3026             if(iqc.ne.0)then
3027               ncj(1,njc)=nj
3028             else
3029               if(iqj(nj).gt.0)then
3030                 ncj(2,njc)=nj
3031               else
3032                 ncj(1,njc)=nj
3033               endif
3034             endif
3035           else
3036             ncc(1,iprt)=nj
3037           endif
3038         else
3039 
3040           do m=1,2
3041             if(jort.eq.1)then
3042               m1=m
3043             else
3044               m1=3-m
3045             endif
3046             njc=ncc(m1,iprt)
3047             if(njc.ne.0)then
3048               ncj(m,nj)=njc
3049               iqc=iqj(njc)
3050               if(iqc.ne.0)then
3051                 ncj(1,njc)=nj
3052               else
3053                 ncj(3-m,njc)=nj
3054               endif
3055             else
3056               ncc(m1,iprt)=nj
3057             endif
3058           enddo
3059         endif
3060         if(ish.ge.6)then
3061           write (ifch,*)'jet-nj,iprt,icp,iqj(nj),ioj(nj),ncj',
3062      *    nj,iprt,icp,iqj(nj),ioj(nj),ncj(1,nj)
3063           if(icp.eq.9)write (ifch,*)'ncj2',ncj(2,nj)
3064         endif
3065       endif
3066 
3067 2     continue
3068       if(iprt.ne.nfprt)then
3069         icp=idprt(iprt)
3070         ipar=iorprt(iprt)
3071         idau1=idaprt(1,ipar)
3072         idau2=idaprt(2,ipar)
3073         if(ish.ge.6)then
3074           write (ifch,*)'2-iprt,icp,idau1,idau2,ipar',
3075      *    iprt,icp,idau1,idau2,ipar,ncc(1,iprt)
3076           if(icp.eq.9)write (ifch,*)ncc(2,iprt)
3077         endif
3078 
3079         if(idau1.eq.iprt)then
3080           if(icp.eq.9)then                   !g -> gg
3081             ncc(1,ipar)=ncc(1,iprt)
3082             ncc(1,idau2)=ncc(2,iprt)
3083           else
3084             icpar=idprt(ipar)
3085             if(icpar.eq.9)then               !g -> qq~
3086               ncc(jort,ipar)=ncc(1,iprt)
3087             else                             !q -> qg
3088               if(icp*(3-2*jort).gt.0)then
3089                 ncc(2,idau2)=ncc(1,iprt)
3090               else
3091                 ncc(1,idau2)=ncc(1,iprt)
3092               endif
3093             endif
3094           endif
3095           iprt=idau2
3096           goto 1
3097 
3098         else
3099           if(icp.eq.9)then
3100             icpar=idprt(ipar)
3101             if(icpar.eq.9)then                !g -> gg
3102               ncc(2,ipar)=ncc(2,iprt)
3103               ncc(2,idau1)=ncc(1,iprt)
3104             else                              !q -> qg
3105               if(icpar*(3-2*jort).gt.0)then
3106                 ncc(1,ipar)=ncc(1,iprt)
3107                 ncc(1,idau1)=ncc(2,iprt)
3108               else
3109                 ncc(1,ipar)=ncc(2,iprt)
3110                 ncc(1,idau1)=ncc(1,iprt)
3111               endif
3112             endif
3113           else
3114             ncc(3-jort,ipar)=ncc(1,iprt)
3115           endif
3116           iprt=ipar
3117           goto 2
3118         endif
3119       else
3120         if(ish.ge.6)write (ifch,*)'3-iprt,ncc',iprt,ncc(1,iprt)
3121         nqc(1)=ncc(1,nfprt)
3122         if(idprt(nfprt).eq.9)nqc(2)=ncc(2,nfprt)
3123       endif
3124       return
3125       end
3126