Project CMSSW displayed by LXR CMSSW_14_1_X_2024-07-25-2300/​GeneratorInterface/​ReggeGribovPartonMCInterface/​src/​epos-fra.f	Source navigation Diff markup Identifier search General search
	Version: CMSSW_14_1_X_2024-07-25-2300 CMSSW_14_1_X_2024-07-24-2300 CMSSW_14_1_X_2024-07-23-2300 CMSSW_14_1_X_2024-07-22-2300 CMSSW_14_1_X_2024-07-21-2300 CMSSW_14_1_X_2024-07-19-2300 Revert   Change

File indexing completed on 2024-04-06 12:14:02

 c-----------------------------------------------------------------------
       subroutine gakfra(iclu,iret)
 c-----------------------------------------------------------------------
 
       include 'epos.inc'
       include 'epos.incems'
       include 'epos.incsem'
       include 'epos.incpar'
       parameter (eta=0.4,etap=0.1)
       parameter (mxnbr=500,mxnba=5000)
       common /gag/nob,pst(4,0:mxnba),ipst(0:mxnba)
      $     ,nbr,ijb(2,0:mxnbr),xyb(2,0:mxnbr)
      &     ,ptb(4,0:mxnbr),iflb(0:mxnbr),ip(0:mxnbr),co,x,y
 c      common /cpptr/ pptr(4,mxptl),ystr(mxptl)
       double precision p1,p12,p2
       dimension co(12),ind(0:mxnbr),p1(5),p2(5)
       dimension ic2(2),ic1(2),ic(2),fkap(3)
       common/pb/pb /cnsbp/nsbp  /cn8ptl/n8ptl
       common/czz/kky,krm,kdrm,kzrm,pui(3),pdi(3),psi(3),pci(3),zzzz,itp
      &,pduui(3),pdudi(3),pdusi(3),pduci(3),pdddi(3),pddsi(3),pddci(3)
      &,pdssi(3),pdsci(3),pdcci(3)
       double precision psg
       common/zpsg/psg(5)
       logical go
 
       pf(i,j)=(pst(4,i)-pst(3,i))*(pst(4,j)+pst(3,j))
      &     -pst(2,i)*pst(2,j)-pst(1,i)*pst(1,j)
 
       nob=0
       call utpri('gakfra',ish,ishini,4)
       iret=0
 c.....search string in pptl-list...................................
       nkmax=nptl
       nk1=maproj+matarg+1
  2    continue
       if(nclean.gt.0.and.nptl.gt.mxptl/2)then
         nnnpt=0
         do iii=maproj+matarg+1,nptl
           go=.true.
           if(nclean.eq.1.and.istptl(iii).le.istmax)go=.false.
           if(go.and.mod(istptl(iii),10).ne.0)then
             istptl(iii)=99
             nnnpt=nnnpt+1
           endif
         enddo
         if(nnnpt.gt.mxptl-nptl)then
           nptl0=nptl
           call utclea(maproj+matarg+1,nptl0)
           nkmax=nptl
           nk1=maproj+matarg+1
 c          write(ifch,'(a)')' ---------after utclea------------'
 c          call blist('&',nk1,nkmax)
         endif
       endif
 
       do while(nk1.le.nkmax)
         if(istptl(nk1).eq.20)then
           nk2=nk1+1
           do while(idptl(nk2).eq.9)
             nk2=nk2+1
           enddo
           goto 3                !ok, string from nk1 to nk2
         else
           nk1=nk1+1
         endif
       enddo
 
       goto 9999                 !no more string
 
 
 c.......................................................................
 c                       decay string  nk1 - nk2
 c.......................................................................
 
  3    nob=-1
       if(ish.ge.3)then
         write(ifch,'(/1x,25a1/5x,a/1x,25a1/)')
      *  ('-',k=1,25),'string decay',('-',k=1,25)
         write(ifch,'(a)')' ---------string------------'
         call blist('&',nk1,nk2)
       endif
 
       itp=ityptl(nk1)
       if(iLHC.eq.1)then
         if(jorptl(nk1).ne.0.or.jorptl(nk2).ne.0)then
           kky=1
         else
           kky=0
         endif
       else
         kky=0
 c      if(nk2-nk1.gt.1.or.iappl.eq.6)kky=1
         if((itp.ge.30.and.itp.le.39).or.iappl.eq.6)kky=1
       endif
       krm=0
       if((itp.ge.50.and.itp.le.59).or.(itp.ge.40.and.itp.le.49))krm=1
       kdrm=0
       if(itp.eq.43.or.itp.eq.53)kdrm=1
 c Excited remnant due to split
       kzrm=0
       if(itp.eq.44.or.itp.eq.54.or.itp.eq.46.or.itp.eq.56)kzrm=1
 
       do k=1,5
         p1(k)=0d0
       enddo
       do i=nk1,nk2
         do k=1,5
           p2(k)=dble(pptl(k,i))
         enddo
         p2(4)=sqrt(p2(3)**2+p2(2)**2+p2(1)**2+p2(5)**2)
         do k=1,4
           p1(k)=p1(k)+p2(k)
         enddo
       enddo
       p1(5)=(p1(4)-p1(3))*(p1(4)+p1(3))-p1(2)**2-p1(1)**2
       if(p1(5).gt.0.d0)then
         p1(5)=sqrt(p1(5))
       else
         iorrem=iorptl(nk1)
         write(*,*)'Precision problem in gakfra, p:',
      &             p1(5),p1(4),p1(3),p1(2),p1(1)
         write(*,*)'origin : ',iorrem
         p1(5)=0.d0
 c        if(iorrem.ne.0.and.
 c     &    (ityptl(iorrem).eq.40.or.ityptl(iorrem).eq.10))
 c     &    p1(5)=dble(pptl(5,iorrem))
         if(iorrem.ne.0)then
           write(*,*)'string mass : ',ityptl(iorrem),pptl(5,iorrem)
           p1(5)=dble(pptl(5,iorrem))
         endif
       endif
       do k=1,5
         psg(k)=p1(k)
       enddo
 
 c.....pbreak ................................
 
       pbi=pbreak
       zz=0.
       if(iLHC.eq.1)then
         x=p1(5)                 !sub string mass (energy in cms) for e+e- and pp
         if(pbi.gt.0.d0.or.(pbi.gt.-0.99999d0.and.kky.eq.0))then
 c around 0.4 for soft because of pi0 spectra for CR and pp multiplicity 
 c and pt spectra at RHIC
           pb=abs(pbi)
         else
 c fit is important not only in e+e- but for all particle production
           pb0=abs(pbi)
 
           pb=pb0+(1.-exp(-x/40.))*(pbreakg-pb0)
 
 c       {{14,0.25},{22,0.22},{34,0.22},{91.2,0.15}} !new
         endif
       else
         x0=1.
         x=sqrt(qsqptl(nk1))     !string energy
         if(iappl.ne.6)then  
 c in hadronic collisions, strings are substring of big string from a Pomeron
 c the ratio Energy_tot / Energy_sub is used to quantify the modification of fkappa 
           x0=x                  !pomeron energy
           x=p1(4)               !sub string energy
         endif
 c       if(pbi.gt.0.d0.or.(pbi.gt.-0.99999d0.and.krm.eq.1))then 
         if(pbi.gt.0.d0.or.(pbi.gt.-0.99999d0.and.kky.eq.0))then
 c around 0.4 for soft because of pi0 spectra for CR and pp multiplicity 
 c and pt spectra at RHIC
           pb=abs(pbi)
         else
 c fit is important not only in e+e- but for all particle production
           pb0=0.33
 
 c       pb=0.065 + 3600./x**3 - 400./x**2 + 17./x   !<--------
           pb=pb0*(1.+exp(-x**2/1000))*0.5
 
 c       {{14,0.25},{22,0.22},{34,0.22},{91.2,0.15}} !new
           if(iappl.ne.6)then
             zz=x0/x             !energy dependence from ratio
             if(iLHC.eq.1)zz=log(max(1.,zz))
           endif
         endif
       endif
 
       if(iLHC.eq.1)then
         fkap(1)=fkappa
         fkap(3)=fkappag
         fkap(2)=0.5*(fkap(1)+fkap(3))
       else
         fkap(1)=fkappa
         if(zz.gt.0.)then
           zzz=min(fkamax,fkainc*zz)
           zzz=sqrt(1.+zzz)
           if(iLHC.eq.0)pb=pb/zzz
           fkap(1)=fkap(1)*zzz         !fix increase of ap and K at Tevatron
         endif
 c      write(*,*)'pb',kky,krm,kdrm,x,x0,pb,zz,fkap
       endif
 
 c      write(*,*)'pb',kky,krm,kdrm,x,pb,zpaptl(1,nk1)+zpaptl(1,nk2)
 
       zzzz=1.
       if(isplit.eq.1.and.iappl.ne.6)then       !for pt
         zzzz=0.
         zzz=0.
 c        if(krm.eq.0)then  !both hard and soft (not good for mesons for HERA-B or pp RHIC, but good for baryons)
         if(kky.eq.1)then   !only hard
 c        if((kky.eq.1.or.(iLHC.eq.1.and.krm.eq.0)))then   
 c          zzzz=zipinc*(zpaptl(1,nk1)+zpaptl(1,nk2))
 c          zzzz=sqrt(1.+zzzz)
           zz=(zpaptl(1,nk1)*(1.+zodinc*log(zpaptl(2,nk2)))     !zpaptl(2,nk1)=# of connected pom
      &       +zpaptl(1,nk2)*(1.+zodinc*log(zpaptl(2,nk2))))
 c          print *,zpaptl(2,nk1),zpaptl(2,nk2)
 c     &           ,1./cosh(0.007*max(0.,zpaptl(2,nk1)-50.))
 c reduce zz before fusion (make it proportionnal to b and E because at low energy not so many fragments are used and then if only one fragment is used all others particles are fragmented without particular effect 
           if(iclu.eq.1)zz=zz/cosh(0.0085*max(0.,zpaptl(2,nk1)-30.))   
           zzzz=zipinc*zz
           if(iLHC.eq.1)then
             zzz=min(fkamax,fkainc*zz)
             zzz=sqrt(1.+zzz)
             fkap(1)=fkap(1)*zzz !fix increase of ap and K at Tevatron
             fkap(2)=fkap(2)*zzz 
             fkap(3)=fkap(3)*zzz
           endif
 c          print *,bimevt,zzzz,zz,0.75*(zpaptl(1,nk1)+zpaptl(1,nk2))
 c          zzz=fkainc*zz
         elseif(kzrm.eq.1)then
           zz=zpaptl(1,nk1)
           zzzz=0. !zopinc*zz
           zzz=0.!zodinc*zz
           zzz=sqrt(1.+zzz)
           fkap(1)=fkap(1)*zzz
           pb=pb/zzz
         endif
         zzzz=sqrt(1.+zzzz)
       endif
       if(ish.ge.4)write(ifch,*)"String parameters:"
      &                      ,iclu,itp,krm,kdrm,kzrm,fkap,zzzz,pb,x
       if(fkap(1).le.0..or.zzzz.le.0.)
      &   call utstop("fkappa<=0. not allowed !&")
 
       igmx=1
       if(iLHC.eq.1)igmx=3
       do ig=1,igmx
       pui(ig)=1.
       pdi(ig)=pui(ig)*exp(-pi*difud/fkap(ig))    !assymmetric u and d relevant because of assymmetry Kch / K0 in e+e-
       psi(ig)=pui(ig)*exp(-pi*difus/fkap(ig))
       pduui(ig)=pui(ig)*exp(-pi*difuuu/fkap(ig))
       pdudi(ig)=pui(ig)*exp(-pi*difuud/fkap(ig))
       pdusi(ig)=pui(ig)*exp(-pi*difuus/fkap(ig))
       pdddi(ig)=pui(ig)*exp(-pi*difudd/fkap(ig))
       pddsi(ig)=pui(ig)*exp(-pi*difuds/fkap(ig))
       pdssi(ig)=pui(ig)*exp(-pi*difuss/fkap(ig))
       if(nrflav.gt.3)then
         pci(ig)=pui(ig)*exp(-pi*difuc/fkap(ig))
         pduci(ig)=pui(ig)*exp(-pi*difuuc/fkap(ig))
         pddci(ig)=pui(ig)*exp(-pi*difudc/fkap(ig))
         pdsci(ig)=pui(ig)*exp(-pi*difusc/fkap(ig))
         pdcci(ig)=pui(ig)*exp(-pi*difucc/fkap(ig))
       else
         pci(ig)=0.
         pduci(ig)=0.
         pddci(ig)=0.
         pdsci(ig)=0.
         pdcci(ig)=0.
       endif
       enddo
 
 c.....light like ...............................................
 
       if(ish.ge.6) call blist("before light like&",nk1,nk2)
       if(rangen().lt.0.5)then
         i1=nk1
         i2=nk1+1
       else
         i1=nk2-1
         i2=nk2
       endif
       ii=1
       do while(ii.ge.0)
         do i=1,4
           p2(i)=dble(pptl(i,i1))+dble(pptl(i,i2))
         enddo
 c       p2(5)=sqrt(p2(4)**2-p2(3)**2-p2(2)**2-p2(1)**2)
         am1=pptl(5,i1)**2
         am2=pptl(5,i2)**2
         am12=max(0d0,p2(4)**2-p2(3)**2-p2(2)**2-p2(1)**2-am1-am2)/2
         if(am12**2.gt.am1*am2)then
           ak1=.5*((am2+am12)/sqrt(am12**2-am1*am2))-.5
           ak2=.5*((am1+am12)/sqrt(am12**2-am1*am2))-.5
         else
           ak1=0.
           ak2=0.
         endif
         if(ish.ge.6)write(ifch,'(a,2i4,9f12.6)') 'overlaps'
      $       ,i1,i2,ak1,ak2,sqrt(2.*am12)
      &       ,am1,am2
         do j=1,4
           a1=(1.+ak1)*pptl(j,i1)-ak2*pptl(j,i2)
           a2=(1.+ak2)*pptl(j,i2)-ak1*pptl(j,i1)
           pptl(j,i1)=a1
           pptl(j,i2)=a2
         enddo
         pptl(5,i1)=0.
         pptl(5,i2)=0.
         if(nk2-nk1.eq.1) ii=-1
         ii=ii-1
         if (nk1.eq.i1)then
           i1=nk2-1
           i2=nk2
         else
           i1=nk1
           i2=nk1+1
         endif
       enddo
 
 c.....inverse  ...............................................
 
 c      if(nk2.eq.nk1+1.and.nptl.lt.mxptl)then
 c        call utrepla(nptl+1,nk2)
 c        call utrepla(nk2,nk1)
 c        call utrepla(nk1,nptl+1)
 c      endif
 
       if(ish.ge.6) call blist("after  light like&",nk1,nk2)
 
 c.....on-shell and copy ...............................................
 
       if(ish.ge.6) write (ifch,*) 'on-shell check'
       do i=nk1,nk2
         do k=1,5
           p2(k)=dble(pptl(k,i))
         enddo
         p12=p2(4)
         p2(4)=sqrt(p2(3)**2+p2(2)**2+p2(1)**2+p2(5)**2)
         if(abs(p12-p2(4))/max(.1d0,p12,p2(4)).ge.1e-1.and.ish.ge.2)then
           write(ifch,*)'warning: on-shell problem:'
      &           ,i, idptl(i),(pptl(k,i),k=1,4),' (',sngl(p2(4)),') '
      &           ,pptl(5,i), istptl(i)
         endif
         if(ish.ge.6)  write (ifch,*)  p12 ,'->',p2(4)
         call utlob2(1,p1(1),p1(2),p1(3),p1(4),p1(5)
      $       ,p2(1),p2(2),p2(3),p2(4),60)
         f=1.0
         if(i.ne.nk1.and.i.ne.nk2)f=0.5
         nmax=1
         ff=1.
         aemax=1000.              ! here max band mass
 c          write(ifch,*)"test",f,p2(4),aemax,f*p2(4).gt.aemax
         if(f*p2(4).gt.aemax) then
           nmax = int(f*p2(4)/aemax)+1
           ff=1./float(nmax)
 c          print *,"nmax",nmax
 c          write(ifch,*)"nmax",nmax
         endif
         nn=1
         do while(nn.le.nmax)
           f=.5
           if(i.eq.nk1.and.nn.eq. 1  ) f=1.
           if(i.eq.nk2.and.nn.eq.nmax) f=1.
           nob=nob+1
           if(nob.gt.mxnba-7)stop'gakfra: increase parameter mxnba '
 c         if(nmax.ge.2) print *, nob,ff,f,i
 c         if(nmax.ge.2) write(ifch,*)'nob',nob,ff,f,i
           do k=1,4
             pst(k,nob)=ff*f*p2(k)
             ipst(nob)=i
           enddo
           nn=nn+1
         enddo
       enddo                     !i
 
       do i1=nob-1,1,-1          ! copy again gluons
         nob=nob+1
         if(nob.gt.mxnba-7)stop'gakfra: increase parameter mxnba '
         do k=1,4
           pst(k,nob)=pst(k,i1)
           ipst(nob)=ipst(i1)
         enddo
       enddo
       nob=nob+1                 ! nob is number of kinks
       if(nob.gt.mxnba-7)stop'gakfra: increase parameter mxnba '
       if(ish.ge.6) then
         write (ifch,*) 'bands:'
         do i=0,nob-1
           write (ifch,'(i4,4g20.13)') i,(pst(k,i),k=1,4)
         enddo
       endif
 
 c.....total momentum...............................................
       do k=1,4
         pst(k,nob)=0.
       enddo
       do i=0,nob-1
         do k=1,4
           pst(k,nob)=pst(k,nob)+pst(k,i)
         enddo
       enddo
 
 
 
 c.....write total string on pptl - list.............................
       nptl=nptl+1
       if(nptl.gt.mxptl)call utstop('gakfra: mxptl too small ...&')
       call idtr5(idptl(nk1),ic1)
       call idtr5(idptl(nk2),ic2)
       ic(1)=ic1(1)+ic2(1)
       ic(2)=ic1(2)+ic2(2)
       idptl(nptl) = 8*10**8+ ic(1)*100 + ic(2)/100 !nbr+1
       istptl(nptl)=10
       n8ptl=nptl
       do i=1,5
         pptl(i,nptl)=p1(i)
       enddo
       if(ish.ge.3)then
         write(ifch,'(a)')' ---------total string------------'
         call blist('&',nptl,nptl)
       endif
 
 c....................................................................
       ijb(1,0)=-1
       ijb(2,0)=-1
       ijb(1,1)=-1
       ijb(2,1)=-1
       iflb(0)=-idptl(nk1)
       iflb(1)= idptl(nk2)
       do i=1,4
         ptb(i,0)=0.
         ptb(i,1)=0.
       enddo
 
 c...................light string....................................
       amm=ammin(idptl(nk1),idptl(nk2))
 
       if(sqrt(max(0.,pf(nob,nob))).lt.amm)then
         id=idsp(  idptl(nk1),idptl(nk2) )
           if(ish.ge.1)then
        write (ifch,*)
      .       '-------string too light to decay--------'
           write (ifch,*) id,p1(5),amm
           write (ifch,*) id, sqrt(max(0.,pf(nob,nob))) ,amm
        endif
       if(id.ne.0.or.iLHC.ne.1)then
         am=sqrt(max(0.,pf(nob,nob)))
         call idress(id,am,idr,iadj)
        if(ish.ge.1)write (ifch,*) id,am,idr,iadj
         nptl=nptl+1
         if(nptl.gt.mxptl)
      &       call utstop('gakfra (light): mxptl too small ...&')
         do i=1,5
           pptl(i,nptl)=p1(i)
         enddo
         do i=nk1,nk2
           istptl(i)=21
           ifrptl(1,i)=n8ptl
           ifrptl(2,i)=0
         enddo
         istptl(n8ptl)=29        !code for fragmented string
         ityptl(n8ptl)=ityptl(nk1)
         itsptl(n8ptl)=itsptl(nk1)
         rinptl(n8ptl)=-9999
         iorptl(n8ptl)=nk1
         jorptl(n8ptl)=nk2
         ifrptl(1,n8ptl)=n8ptl+1
         ifrptl(2,n8ptl)=nptl
         ityptl(nptl)=ityptl(nk1)
         itsptl(nptl)=itsptl(nk1)
         rinptl(nptl)=-9999
         istptl(nptl)=0
         iorptl(nptl)=n8ptl
         jorptl(nptl)=0
         ifrptl(1,nptl)=0
         ifrptl(2,nptl)=0
         if(idr.ne.0)then
           idptl(nptl)=idr
         else
           idptl(nptl)=id
         endif
         do j=1,4
           xorptl(j,nptl)=xorptl(j,nk1)
         enddo
         tivptl(1,nptl)=xorptl(4,nptl)
         call idtau(idptl(nptl),pptl(4,nptl),pptl(5,nptl),taugm)
         tivptl(2,nptl)=tivptl(1,nptl)+taugm*(-alog(rangen()))
         if(abs(p1(5)-am).gt.0.01)then
 c          write (*,*) 'light string  ---  particle off-shell!!!!'
 c          write (*,*) idr,'  mass:',p1(5),'  should be:',am
         endif
         goto 98
       endif
       endif
 
 c................search breakpoints...................................
 
       n=0
       nsbp=1
       nbr=0
       iok=0
 
  11   continue
       if(nsbp.gt.10000)then
         if(ish.ge.1)then
         write(*,*)'Gakfra : string can not be fragmented, redo event !'
         write(*,*)nk1,idptl(nk1),nk2,idptl(nk2),istptl(nk1)
         write(ifch,*)'Gakfra : redo event !'
         endif
         iret=1
         goto 9999
       endif
 
       !----------------------!
       call gaksbp(0,1,iok)
       !----------------------!
       nbrtry=0
       goto 10
 
 c..............delete breakpoint....................................
  9    if(ish.ge.4)write(ifch,*) 'delete breakpoint',n,' -> '
      &,nbr-1,' breakpoints left'
       call gakdel(n)              !hier loeschen
 
 c..............no more breakpoints -> redo..........................
       if(nbr.eq.0)then
         nsbp=nsbp+1
         if(ish.ge.3)write (ifch,*)
      &    'no breakpoints left ... try again (',nsbp,')'
         if(ish.ge.3)write (ifch,*)' '
         goto 11
       endif
 
 c................make index list of breakpoints to adjust...........
  10   continue
       nbrtry=nbrtry+1
       nind=0
       do i=1,nbr
         if(ip(i).eq.0.or.ip(i+1).eq.0)then
           nind=nind+1
           ind(nind)=i
         endif
       enddo
       if(nbrtry.gt.1000000)then
         if(ish.ge.1)then
         write(*,*)'Gakfra : string can not be fragmented, redo event !'
         write(*,*)nk1,nk2,nbr,nind,pb
         endif
         iret=1
         goto 9999
       endif
 
 c.....no more breakpoint to adjust...............................
       if(nind.eq.0) goto 100
 
 c................................................................
       if(ish.ge.5)then
         write(ifch,*) 'breakpoints:'
         write(ifch,'(i3,i5)') 0,-iflb(0)
         do i=1,nbr
           write(ifch,'(i3,2i5,1x,4e12.4,2x,2i3,2f6.3)')
      &      i,iflb(i),-iflb(i),(ptb(j,i),j=1,4)
      &         ,ijb(1,i),ijb(2,i),xyb(1,i),xyb(2,i)
         enddo
         write(ifch,'(i3,i5)') nbr+1,iflb(nbr+1)
       endif
 
 c.....choose breakpoint, calculate masses, lambda..................
       r=rangen()
       nn=1+int(r*float(nind))
 c      nn=1+(nind-1)*int(r+0.5)
       n=ind(nn)
       if(ish.ge.4)write(ifch,*) 'choose breakpoint',n
       nl=n-1
       nr=n+1
       do while (nl.gt.0.and.ip(nl).eq.0)
          nl=nl-1
       enddo
       do while (nr.lt.nbr+1.and.ip(nr+1).eq.0)
          nr=nr+1
       enddo
       if(ish.ge.6)write (ifch,*) 'nl,n,nr:',nl,n,nr
 c      print *,'------------------------------------',1
       call gaksco(n-1,n,n+1,ijb(1,n),ijb(2,n),x1,x2,y1,y2)
       if(x2.le.x1.or.y2.le.y1)goto 9
 cc      x=(xyb(1,n)-x1)/(x2-x1)
 cc      y=(xyb(2,n)-y1)/(y2-y1)
       x=xyb(1,n)
       y=xyb(2,n)
       aml2=co(1)+x*co(2)+y*co(3)+x*y*co(4)
       amr2=co(5)+x*co(6)+y*co(7)+x*y*co(8)
       ala2=co(9)+x*co(10)+y*co(11)+x*y*co(12)
 
 c.....determine id of both particles..............................
       aml=sqrt(max(0.,aml2))
       idl=idsp( -iflb(n-1) , iflb(n)   )
       amr=sqrt(max(0.,amr2))
       idr=idsp( -iflb(n)   , iflb(n+1) )
 
 c.....if mass to small (because of spacelike pt)  reject..........
 c      amin=0.0
 c      if (aml2.le.amin.or.amr2.le.amin) goto 9
 
 c.....if no left or right ptl id -> reject.........
       if(idl.eq.0.or.idr.eq.0)then
         if(ish.ge.5)write(ifch,*)'no left or right ptl id'
         goto 9
       endif
 
       iadjl=0
       iadjr=0
       idrl=0
       idrr=0
 
 c.....if no adjusted mass on left side, check for resonance...........
       if(ip(n)  .eq.0) then
         aml=sqrt(max(0.,aml2+0.*min(0.,amr2))) !!!????
         amlxx=aml
         call idress(idl,aml,idrl,iadjl)
         r=rangen()
         if(idrl.eq.110.and.r.lt.0.5)then
           if (r.gt.eta+etap) goto 9      !rare numerical errors
           idl=220
           aml=.549
           if(r.lt.etap)aml=.958
           amlxx=aml
           call idress(idl,aml,idrl,iadjl)
           iadjl=1
         endif
         if(ish.ge.5)write(ifch,'(a,i5,2f12.6,1i10,i5)')
      &    ' l:  ',idl,amlxx,aml,idrl,iadjl
       else
         if(ish.ge.5)write(ifch,'(a,i5,2f12.6,i10,a5)')
      &    ' l:  ',idl,aml,aml,ip(n),'ok'
       endif
 
 c.....if no adjusted mass on right side, check for resonance...........
       if(ip(n+1).eq.0) then
         amr=sqrt(max(0.,amr2+0.*min(0.,aml2))) !!!????
         amrxx=amr
         call idress(idr,amr,idrr,iadjr)
         r=rangen()
         if(idrr.eq.110.and.r.lt.0.5)then
           if (r.gt.eta+etap) goto 9    !rare numerical errors
           idr=220
           amr=.549
           if(r.lt.etap)amr=.958
           amrxx=amr
           call idress(idr,amr,idrr,iadjr)
           iadjr=1
         endif
         if(ish.ge.5)write(ifch,'(a,i5,2f12.6,1i10,i5)')
      &    ' r:  ',idr,amrxx,amr,idrr,iadjr
       else
         if(ish.ge.5)write(ifch,'(a,i5,2f12.6,i10,a5)')
      &    ' r:  ',idr,amr,amr,ip(n+1),'ok'
       endif
 
 c.....mass adjustments.........................................
       iok=0
       if(ip(n+1).ne.0)then  !.........adjusted mass on right side
         if(idrl.eq.0)then
           call gaksbp(n-1,n,iok)
           if(iok.eq.1)goto 9
           goto 10
         endif
         if(iadjl.eq.1)then
            if(ish.ge.5)write(ifch,*)'mass adjustment 1'
            n2=n+1
            call gakanp(n-1,n,n2,aml**2,amr**2,0.,iok)
         endif
         if(iok.eq.1)goto 9
         ip(n)=idrl
       elseif(ip(n).ne.0)then !.........adjusted mass on left side
         if(idrr.eq.0)then
           call gaksbp(n,n+1,iok)
           if(iok.eq.1)goto 9
           goto 10
         endif
         if(iadjr.eq.1)then
            if(ish.ge.5)write(ifch,*)'mass adjustment 2'
            n2=n+1
            call gakanp(n-1,n,n2,aml**2,amr**2,0.,iok)
         endif
         if(iok.eq.1)goto 9
         ip(n+1)=idrr
       else       !.........adjusted mass neither left nor right
         if(idrr.eq.0.and.idrl.eq.0)then
           call gaksbp(n,n+1,iok)
           if(iok.eq.1)goto 9
           call gaksbp(n-1,n,iok)
           if(iok.eq.1)goto 9
           goto 10
         elseif(idrl.ne.0.and.idrr.eq.0)then
           if(iadjl.eq.1) then
            if(ish.ge.5)write(ifch,*)'mass adjustment 3'
            call gakanp(n-1,n,nr,aml**2,0.,ala2,iok)
           endif
         elseif(idrl.eq.0.and.idrr.ne.0)then
           if(iadjr.eq.1) then
            if(ish.ge.5)write(ifch,*)'mass adjustment 4'
            n2=n+1
            call gakanp(nl,n,n2,0.,amr**2,ala2,iok)
           endif
         else  !if(idrl.ne.0.and.idrr.ne.0)then
           if(iadjl.eq.1.or.iadjr.eq.1) then
            if(ish.ge.5)write(ifch,*)'mass adjustment 5'
            n2=n+1
            call gakanp(n-1,n,n2,aml**2,amr**2,0.,iok)
           endif
         endif
         if(iok.eq.1)goto 9
         ip(n)=idrl
         ip(n+1)=idrr
       endif
       if(ish.ge.4)then
         write(ifch,*) 'left/right particles:'
      &         ,ip(n),aml,'  ',ip(n+1),amr
       endif
       goto 10
 
 c................................................................
 c                         end of string decay
 c................................................................
 
 c.....final list...............................................
  100  if(ish.ge.4)then
         write(ifch,*) '   ************ OK **************'
         write(ifch,*) 'final breakpoints:'
         write(ifch,'(i3,i5)') 0,-iflb(0)
         do i=1,nbr
           write(ifch,'(i3,2i5,1x,4e12.4,2x,2i3,2f6.3)')
      &      i,iflb(i),-iflb(i),(ptb(j,i),j=1,4)
      &         ,ijb(1,i),ijb(2,i),xyb(1,i),xyb(2,i)
         enddo
         write(ifch,'(i3,i5)') nbr+1,iflb(nbr+1)
       endif
 
 c.....write particles in pptl-list................................
       if(ish.ge.3)then
         write(ifch,'(a)')' ---------produced particles---------'
       endif
 
       nptlini=nptl
       if(nptlini+nbr+1.gt.mxptl)
      &       call utstop('gakfra (end): mxptl too small ...&')
       do i=1,nbr+1
          nptl=nptl+1
          if(i.lt.nbr+1)then     !particle = left side of breakpoints
            call gaksco(i-1,i,i+1,ijb(1,i),ijb(2,i),x1,x2,y1,y2)
 c           taubrr=pst(4,nob+7)+x*pst(4,nob+8)+y*pst(4,nob+9)
            x=xyb(1,i)
            y=xyb(2,i)
            aml2=co(1)+x*co(2)+y*co(3)+x*y*co(4)
            amr2=co(5)+x*co(6)+y*co(7)+x*y*co(8)
 
            ala2=co(9)+x*co(10)+y*co(11)+x*y*co(12)
            aml=sign(sqrt(abs(aml2)),aml2)
            amr=sign(sqrt(abs(amr2)),amr2)
            if(aml.le.0.d0.or.amr.le.0.d0)then
              if(ish.ge.4)write(ifch,*)
      & 'Negative mass, fragmentation not OK -> redo ...'
              n=i
              nptl=nptlini
              goto 9
            endif
            qsqptl(nptl)=ala2
            pptl(5,nptl)=aml
            do j=1,4
              pptl(j,nptl)=pst(j,nob+1)-
      &            x*pst(j,nob+2)+y*pst(j,nob+3)
 c             pptr(j,nptl)=ptb(j,i)
            enddo
          else                   !last particle = right side of last breakpoint
            pptl(5,nptl)=amr
            qsqptl(nptl)=0.
            do j=1,4
              pptl(j,nptl)=pst(j,nob+4)+
      &            x*pst(j,nob+5)-y*pst(j,nob+6)
 c             pptr(j,nptl)=ptb(j,i) !should be zero
            enddo
 c           taubrr=0.
          endif
          idptl(nptl)=ip(i)
          if(ish.ge.7)call blist('&',nptl,nptl)
          if(pptl(4,nptl).le.0.d0)then
            if(ish.ge.4)write(ifch,*)
      & 'Negative energy, fragmentation not OK -> redo ...'
            n=i
            nptl=nptlini
            goto 9
          endif
 
       enddo
 
       nptl=nptlini
 
       if(ish.ge.7)then
         write(ifch,'(a)')' ---------produced particles---------'
       endif
 
 
       do i=1,nbr+1
          nptl=nptl+1
 
          call utlob4(-1,p1(1),p1(2),p1(3),p1(4),p1(5)
      $        ,pptl(1,nptl),pptl(2,nptl),pptl(3,nptl),pptl(4,nptl))
 c         call utlob4(-1,p1(1),p1(2),p1(3),p1(4),p1(5)
 c     $        ,pptr(1,nptl),pptr(2,nptl),pptr(3,nptl),pptr(4,nptl))
 
 
 c........Origin..................................................
          istptl(nptl)=0
          iorptl(nptl)=n8ptl
          jorptl(nptl)=0
          ifrptl(1,nptl)=0
          ifrptl(2,nptl)=0
 
          r=rangen()
          tauran=-taurea*alog(r)*pptl(4,nptl)/pptl(5,nptl)
 c         if(jpsi.ge.1)tauran=0.
 c         tauran=max(taubrl,taubrr) !take formation time from string-theory
          do j=1,4
            xorptl(j,nptl)=xorptl(j,nk1)+pptl(j,nptl)
      &       /pptl(4,nptl)*tauran
          enddo
          tivptl(1,nptl)=xorptl(4,nptl)
          call idtau(idptl(nptl),pptl(4,nptl),pptl(5,nptl),taugm)
          tivptl(2,nptl)=tivptl(1,nptl)+taugm*(-alog(rangen()))
          ityptl(nptl)=ityptl(nk1)
          itsptl(nptl)=itsptl(nk1)
          rinptl(nptl)=-9999
          if(ish.ge.3)call blist('&',nptl,nptl)
 ctp060829          taubrl=taubrr
       enddo
       iorptl(n8ptl)=nk1
       jorptl(n8ptl)=nk2
       ityptl(n8ptl)=ityptl(nk1)
       do i=nk1,nk2
          istptl(i)=21
          ifrptl(1,i)=n8ptl
          ifrptl(2,i)=0
       enddo
       istptl(n8ptl)=29          !code for fragmented string
       ifrptl(1,n8ptl)=n8ptl+1
       ifrptl(2,n8ptl)=nptl
       if(ish.ge.5)then
         write(ifch,*)'string momentum sum:'
         do kk=1,5
           pptl(kk,nptl+1)=0
           do ii=n8ptl+1,nptl
             pptl(kk,nptl+1)=pptl(kk,nptl+1)+pptl(kk,ii)
           enddo
         enddo
         call alist2('&',n8ptl,n8ptl,nptl+1,nptl+1)
       endif
 
 
 
 c.....another string?.........................................
  98   nk1=nk2+1
       goto 2                    !next string
 
 c.....end of fragmentation.....................................
  9999 continue
       call utprix('gakfra',ish,ishini,4)
       return
       end
 
 c---------------------------------------------------------------------
       subroutine gaksbp(ibr1,ibr2,iok)
 c---------------------------------------------------------------------
       ! search break points
       !-----------------------------------------
       ! nbr ... number of break points
       !
       !
       !
       !--------------------------------------------------------------
       include 'epos.inc'
 
       parameter (mxnbr=500,mxnba=5000)
       common /gag/nob,pst(4,0:mxnba),ipst(0:mxnba)
      $     ,nbr,ijb(2,0:mxnbr),xyb(2,0:mxnbr)
      &     ,ptb(4,0:mxnbr),iflb(0:mxnbr),ip(0:mxnbr),co,x,y
       common/pb/pb /cnsbp/nsbp /cn8ptl/n8ptl
       double precision ax,ay,az,ae,am,bx,by,bz,be
       dimension co(12)
       logical weiter
       common/czz/kky,krm,kdrm,kzrm,pui(3),pdi(3),psi(3),pci(3),zzzz,itp
      &,pduui(3),pdudi(3),pdusi(3),pduci(3),pdddi(3),pddsi(3),pddci(3)
      &,pdssi(3),pdsci(3),pdcci(3)
       double precision psg
       common/zpsg/psg(5)
 c      dimension pleft(5),pright(5)
 
       pf(i,j)=pst(4,i)*pst(4,j)-pst(3,i)*pst(3,j)
      &     -pst(2,i)*pst(2,j)-pst(1,i)*pst(1,j)
       mmod(i)=mod(mod(i,nob)+nob,nob)
 
       call utpri('gaksbp',ish,ishini,5)
 
 
 
       ib1=ibr1
       ib2=ibr2
       iside=1
       if(rangen().lt.0.5)iside=2
       if(ish.ge.6)write(ifch,*)'iside:',iside
       Amxx=80./pb
       if(ish.ge.4)write(ifch,*)
      &'search brk between ib1=',ib1,' and ib2=',ib2
       ntry=0
       nbrold=nbr
       Amin=0.
       Amax=Amxx
  26   ntry=ntry+1
       A0=-log(exp(-pb*Amin)+rangen()*(exp(-pb*Amax)-exp(-pb*Amin)))/pb
       if(ish.ge.6)write(ifch,*)'pb, Amin, Amax, A0:',pb, Amin, Amax, A0
       A=A0
       Amaxn=0.
       if(iside.eq.2)goto 51
 c.....................................................................
       if(ib1.eq.0)then          !startwert der j-schleife
          jj=nob-1
       else
          jj=ijb(2,ib1)
       endif
       do while (jj.gt.0)
          if(jj.eq.ijb(2,ib1) )then       !linker brk im Streifen jj?
             y1=xyb(2,ib1)
          else                   !nein
             y1=0.
          endif
          if(jj.eq.ijb(2,ib2))then       !rechter brk im Streifen jj?
             y2=xyb(2,ib2)
          else                   !nein
             y2=1.
          endif
          if(y1.eq.y2) goto 9999
          if(abs(y1-y2)/abs(y1+y2).le.1e-7) goto 9999
          if(ish.ge.6)write(ifch,*)'y1,y2,A',y1,y2,A
                                 !startwert der i-schleife
          if(ib1.eq.0)then       !ohne linken brkpt
             ii=mmod(jj+1)
             if(jj.lt.nob/2)ii=mmod(nob-jj+1)
             if(ib2.le.nbr.and.mmod(ii+nob/2)
      &        .gt.mmod(ijb(1,ib2)+nob/2))then
               if(ish.ge.6)write(ifch,*) 'very special case',ii,jj
               goto12
             endif
          else
             ii=ijb(1,ib1)
             if(jj.lt.nob/2 .and.
      $           mmod(nob-jj+1+nob/2).gt. mmod(ii+nob/2)
      $           )ii=mmod(nob-jj+1)
          endif
          weiter=.true.
          aa=0.
          if(ii.eq.jj) then
             if(ish.ge.6)write(ifch,*) 'Rand erreicht'
             goto 15
          endif
          do while(weiter)
             if(ii.eq.ijb(1,ib1))then    !linker brk im Feld (ii,jj)
                x2=xyb(1,ib1)
             else
                x2=1.
             endif
 
             if(ii.eq.ijb(1,ib2))then    !rechter brk im Feld (ii,jj)
                x1=xyb(1,ib2)
             else
                x1=0.
             endif
             if(x1.eq.x2) goto 9999
             if(abs(x1-x2)/abs(x1+x2).le.1e-7) goto 9999
             f=1.0
             if(ipst(ii).ne.ipst(jj))aa=aa+2.*(x2-x1)*(y2-y1)*f*pf(ii,jj)
             if(ish.ge.6)write(ifch,*)'ii,jj,x1,x2,aa:',ii,jj
      $           ,ipst(ii).ne.ipst(jj),x1,x2,aa
      &           ,pf(ii,jj)
             if(ii.eq.ijb(1,ib2))then
                weiter=.false.
             else
                ii=mmod(ii+1)
                if(ii.eq.jj.or.mmod(ii+jj).eq.0)weiter=.false.
             endif
          enddo
          Amaxn=Amaxn+aa
          if(aa.gt.A)goto 10
          A=A-aa
          if(jj.eq.ijb(2,ib2)) then
             if(ish.ge.6)write(ifch,*) 'brk erreicht'
             goto 15
          endif
  12      jj=mmod(jj-1)
       enddo
       goto 15
 
  10   continue
       yb=A/aa*(y2-y1)+y1
       if(ish.ge.6)write(ifch,*)'jj,yb ok:',jj,yb
       r=rangen()
       ra=aa*r
       if(ish.ge.6)write(ifch,*)'r,ra,aa',r,ra,aa
                                 !nochmal die letzte ii-Schleife
       if(ib1.eq.0)then          !ohne linken brkpt
          ii=mmod(jj+1)
          if(jj.lt.nob/2)ii=mmod(nob-jj+1)
       else
         ii=ijb(1,ib1)
         if(jj.lt.nob/2 .and.
      $       mmod(nob-jj+1+nob/2).gt. mmod(ii+nob/2)
      $       )ii=mmod(nob-jj+1)
       endif
       do while(ra.gt.0.)
          if(ii.eq.ijb(1,ib1))then       !linker brk im Feld (ii,jj)
             x2=xyb(1,ib1)
          else
             x2=1.
          endif
 
          if(ii.eq.ijb(1,ib2))then       !rechter brk im Feld (ii,jj)
             x1=xyb(1,ib2)
          else
             x1=0.
          endif
          f=1.0
          ab=0.
          if(ipst(ii).ne.ipst(jj)) ab=2.*(x2-x1)*(y2-y1)*f*pf(ii,jj)
          if(ish.ge.6)write(ifch,*)'ii,jj,x1,x2,ab,ra:',ii,jj
      $        ,ipst(ii).ne.ipst(jj),x1,x2,ab,ra
          if(ab.gt.ra)then
             xb=ra/ab*(x2-x1)+x1
          else
             ii=mmod(ii+1)
          endif
          ra=ra-ab
       enddo
       if(ish.ge.5)write(ifch,*) 'breakpoint in field ',ii,jj
      & ,' at position ',xb,yb
       goto 95
 
 c......................................................................
       !von rechts
  51   if(ib2.eq.nbr+1)then      !startwert der i-schleife
          ii=nob/2-1
       else
          ii=ijb(1,ib2)
       endif
       do while (ii.ne.nob/2)
          if(ii.eq.ijb(1,ib1))then !linker brk im Streifen (ii)
             x2=xyb(1,ib1)
          else
             x2=1.
          endif
          if(ii.eq.ijb(1,ib2))then !rechter brk im Streifen (ii)
             x1=xyb(1,ib2)
          else
             x1=0.
          endif
          if(x1.eq.x2) goto 9999
          if(abs(x1-x2)/abs(x1+x2).le.1e-7) goto 9999
          if(ish.ge.6)write(ifch,*)'x1,x2 A',x1,x2,A
 
          if(ib2.eq.nbr+1)then
             jj=mmod(ii+1)
             if(ii.gt.nob/2)jj=mmod(nob-ii+1)
             if(ib1.ge.1.and.jj.gt.ijb(2,ib1))then
                if(ish.ge.6)write(ifch,*) 'very special case',ii,jj
                goto 13
             endif
          else
             jj=ijb(2,ib2)
             if(ii.gt.nob/2 .and. mmod(nob-ii+1).gt.jj)jj=mmod(nob-ii+1)
          endif
          weiter=.true.
          aa=0.
          if(ii.eq.jj) then
             if(ish.ge.6)write(ifch,*) 'Rand erreicht'
             goto 15
          endif
          do while(weiter)
             if(jj.eq.ijb(2,ib1))then    !linker brk im Feld (ii,jj)
                y1=xyb(2,ib1)
             else
                y1=0.
             endif
             if(jj.eq.ijb(2,ib2))then !rechter brk im Feld (ii,jj)
                y2=xyb(2,ib2)
             else
                y2=1.
             endif
             if(y1.eq.y2) goto 9999
             if(abs(y1-y2)/abs(y1+y2).le.1e-7) goto 9999
             f=1.0
             if(ipst(ii).ne.ipst(jj))aa=aa+2.*(x2-x1)*(y2-y1)*f*pf(ii,jj)
             if(ish.ge.6)write(ifch,*)'ii,jj,x1,x2,aa:',ii,jj
      $           ,ipst(ii).ne.ipst(jj),x1,x2,aa
      &           ,pf(ii,jj)
             if(jj.eq.ijb(2,ib1))then
                weiter=.false.
             else
                jj=mmod(jj+1)
                if(jj.eq.ii.or.mmod(ii+jj).eq.0)weiter=.false.
             endif
          enddo
          Amaxn=Amaxn+aa
          if(aa.gt.A)goto 14
          A=A-aa
          if(ii.eq.ijb(1,ib1)) then
             if(ish.ge.6)write(ifch,*) 'brk erreicht'
             goto 15
          endif
  13      ii=mmod(ii-1)
       enddo
       goto 15
 
 
 
  14   continue
       xb=A/aa*(x2-x1)+x1
       if(ish.ge.6)write(ifch,*)'ii,xb ok:',ii,xb
       r=rangen()
       ra=aa*r
       if(ish.ge.6)write(ifch,*)'r,ra,aa',r,ra,aa
                                 !nochmal die letzte jj-Schleife
       if(ib2.eq.nbr+1)then
          jj=mmod(ii+1)
          if(ii.gt.nob/2)jj=mmod(nob-ii+1)
       else
          jj=ijb(2,ib2)
          if(ii.gt.nob/2 .and. mmod(nob-ii+1).gt.jj)jj=mmod(nob-ii+1)
       endif
 
       do while(ra.gt.0.)
          if(jj.eq.ijb(2,ib1))then       !linker brk im Feld (ii,jj)
             y1=xyb(2,ib1)
          else
             y1=0.
          endif
 
          if(jj.eq.ijb(2,ib2))then       !rechter brk im Feld (ii,jj)
             y2=xyb(2,ib2)
          else
             y2=1.
          endif
          f=1.0
          ab=0.
          if(ipst(ii).ne.ipst(jj)) ab=2.*(x2-x1)*(y2-y1)*f*pf(ii,jj)
          if(ish.ge.6)write(ifch,*)'ii,jj,x1,x2,ab,ra:',ii,jj
      $        ,ipst(ii).ne.ipst(jj),x1,x2,ab,ra
          if(ab.gt.ra)then
             yb=ra/ab*(y2-y1)+y1
          else
             jj=mmod(jj+1)
          endif
          ra=ra-ab
       enddo
       if(ish.ge.5)write(ifch,*) 'breakpoint in field ',ii,jj
      & ,' at position ',xb,yb
 
  95   continue
 
 c.....breakpoint accepted......................
       nbr=nbr+1
       if(nbr.gt.mxnbr-2) stop 'gaksbp: increase parameter mxnbr'
       do i=nbr+1,ib1+1,-1
          do j=1,2
             ijb(j,i)=ijb(j,i-1)
             xyb(j,i)=xyb(j,i-1)
          enddo
          do k=1,4
             ptb(k,i)=ptb(k,i-1)
          enddo
          iflb(i)=iflb(i-1)
          ip(i)=ip(i-1)
       enddo
       ip(ib1+1)=0
       ip(ib1+2)=0
       ijb(1,ib1+1)=ii
       ijb(2,ib1+1)=jj
       xyb(1,ib1+1)=xb
       xyb(2,ib1+1)=yb
 
 
 cc ..........diquark...............................................
 c
 c
 c      pdiqu=pdiqua
 c      if(kzrm.eq.1.or.krm.eq.0)then
 cc        pdiqu=exp(log(pdiqua)/sqrt(1.+fkainc*(zzzz-1.)))
 cc        pdiqu=exp(log(pdiqua)/(1.+fkainc*(zzzz-1.)))
 c      endif
 cc      if(iappl.eq.1.and.krm.eq.0)then
 cc        pdiqu=pdiquak
 cc      endif
 c      
 c
 c      if(nbr.le.2.and.abs(psg(3)/psg(4)).gt.diqcut)pdiqu=0.
 c
 c      
 c      
 c      if(rangen().lt.pdiqu.and.iabs(iflb(ib1)).lt.6
 c     &  .and.iabs(iflb(ib2)).lt.6)then
 c        jqu=2
 c      else
 c        jqu=1
 c      endif
 c
 
 c ..........flavor...............................................
 
       ig=1
       if(iLHC.eq.1)then
         if(jorptl(ipst(ii)).eq.-9)ig=ig+1
         if(jorptl(ipst(jj)).eq.-9)ig=ig+1
       endif
 
       pu=pui(ig)
       pd=pdi(ig)
       ps=psi(ig)
       pc=pci(ig)
       pduu=pduui(ig)
       pdud=pdudi(ig)
       pdus=pdusi(ig)
       pduc=pduci(ig)
       pddd=pdddi(ig)
       pdds=pddsi(ig)
       pddc=pddci(ig)
       pdss=pdssi(ig)
       pdsc=pdsci(ig)
       pdcc=pdcci(ig)
 
 
 !      print *,pu,pd,ps,pc,difud,difus,difuc
 c      print *,krm,kdrm,nbr,abs(psg(3)/psg(4))
 
 c suppress forward (proj) or backward (targ) diquark production 
 c in non central strings (baryon spectra in pion collisions or NA49 data)
 c ----> this should not be used for meson projectile because of kaon spectra
 c       at 100 GeV (Barton, etc) <--------
 c      if(kdrm.eq.1.and.ib1.le.1.and.abs(psg(3)/psg(4)).gt.strcut)then
 c        ps=0.
 c        pc=0.
 c      endif
 
       if(iLHC.eq.1)then
 c kill baryon production in case of high pt or forward production
         if((ib1.le.1.or.ib2.ge.nbr-1)
      &   .and.sqrt(psg(1)**2+psg(2)**2+psg(3)**2)/psg(4).gt.diqcut)then
           pduu=0.
           pdud=0.
           pdus=0.
           pduc=0.
           pddd=0.
           pdds=0.
           pddc=0.
           pdss=0.
           pdsc=0.
           pdcc=0.
         endif
       else
 c suppress forward (proj) or backward (targ) diquark production 
 c in non central strings (baryon spectra in pion collisions or NA49 data)
         if(kdrm.eq.1.and.ib1.le.1.and.abs(psg(3)/psg(4)).gt.diqcut)then
           pduu=0.
           pdud=0.
           pdus=0.
           pduc=0.
           pddd=0.
           pdds=0.
           pddc=0.
           pdss=0.
           pdsc=0.
           pdcc=0.
         endif
       endif
 
 
 
       pdiqu=pduu+pdud+pdus+pduc+pddd+pdds+pddc+pdss+pdsc+pdcc
       psum=pu+pd+ps+pc+pdiqu
 c      print *,krm,kdrm,kky,zzzz,ps/(pu+pd+ps),pdiqu/psum
 
       r=rangen()*psum
       jqu=1
       ifl2=0
       if(r.gt.psum-pdiqu)then
         jqu=2
         if(r.gt.pu+pd+ps+pc+pduu+pdud+pdus+pduc+pddd+pdds+pddc+pdss+pdsc
      &       .and.pdcc.gt.0.)then
           ifl1=4
           ifl2=4
         elseif(r.gt.pu+pd+ps+pc+pduu+pdud+pdus+pduc+pddd+pdds+pddc+pdss
      &         .and.pdsc.gt.0.)then
           if(rangen().gt.0.5)then
             ifl1=3
             ifl2=4
           else
             ifl1=4
             ifl2=3
           endif
         elseif(r.gt.pu+pd+ps+pc+pduu+pdud+pdus+pduc+pddd+pdds+pddc
      &         .and.pdss.gt.0.)then
           ifl1=3
           ifl2=3
         elseif(r.gt.pu+pd+ps+pc+pduu+pdud+pdus+pduc+pddd+pdds
      &         .and.pddc.gt.0.)then
           if(rangen().gt.0.5)then
             ifl1=2
             ifl2=4
           else
             ifl1=4
             ifl2=2
           endif
         elseif(r.gt.pu+pd+ps+pc+pduu+pdud+pdus+pduc+pddd
      &         .and.pdds.gt.0.)then
           if(rangen().gt.0.5)then
             ifl1=2
             ifl2=3
           else
             ifl1=3
             ifl2=2
           endif
         elseif(r.gt.pu+pd+ps+pc+pduu+pdud+pdus+pduc.and.pddd.gt.0.)then
           ifl1=2
           ifl2=2
           jqu=2
         elseif(r.gt.pu+pd+ps+pc+pduu+pdud+pdus.and.pduc.gt.0.)then
           if(rangen().gt.0.5)then
             ifl1=1
             ifl2=4
           else
             ifl1=4
             ifl2=1
           endif
         elseif(r.gt.pu+pd+ps+pc+pduu+pdud.and.pdus.gt.0.)then
           if(rangen().gt.0.5)then
             ifl1=1
             ifl2=3
           else
             ifl1=3
             ifl2=1
           endif
         elseif(r.gt.pu+pd+ps+pc+pduu.and.pdud.gt.0.)then
           if(rangen().gt.0.5)then
             ifl1=1
             ifl2=2
           else
             ifl1=2
             ifl2=1
           endif
         else
           ifl1=1
           ifl2=1
         endif
       elseif(r.gt.pu+pd+ps.and.pc.gt.0.)then
         ifl1=4
       elseif(r.gt.pu+pd.and.ps.gt.0.)then
         ifl1=3
       elseif(r.gt.pu.and.pd.gt.0.)then
         ifl1=2
       else
         ifl1=1
       endif
 
       if(jqu.eq.2)then  !diquark ------
         ifl=-min(ifl1,ifl2)*1000-max(ifl1,ifl2)*100
       else              !quark ------
         ifl=ifl1
       endif
 
       if(iflb(ib1).lt.0.and.iflb(ib1).ge.-6)ifl=-ifl
       if(iflb(ib1).gt.1000)ifl=-ifl
       iflb(ib1+1)=ifl
       if(ish.ge.5)write(ifch,20) ig,ifl,pu,pd,ps,pc,pdiqu/psum
      &,sqrt(psg(1)**2+psg(2)**2+psg(3)**2)/psg(4).gt.diqcut
      &,pduu,pdud,pdus,pduc,pddd,pdds,pddc,pdss,pdsc,pdcc
  20   format('ig,flavor,pu,pd,ps,pc,pdiqu,dcut:',i2,i6,5g13.5,l2/
      &,'  diquark u:',4g15.5,/,'  diquark d:',3g15.5,/
      &,'  diquark s:',2g15.5,10x,'diquark c:',g15.5)
 c..............................pt.......................................
 c      if(krm.ne.1.and.kdrm.ne.1)then
 c        icub=ib1+1
 c        !---------------------------------------------------------------------------
 c        !  ib1+1 is the current break point index
 c        !             (between 1 and nbr)
 c        !  ijb(1,ib1) and ijb(2,ib1) are band indices
 c        !         each index from 0 to nob-1 (nob= number of bands)
 c        !         0 is left, then come the gluons, then antiquark, then agin gluons
 c        !         like q - g - g - g - ~q - g - g - g
 c        !--------------------------------------------------------------------------
 c        call gaksco(icub-1,icub,icub+1
 c     &    ,ijb(1,icub),ijb(2,icub),x1,x2,y1,y2)
 c        x=xyb(1,icub)
 c        y=xyb(2,icub)
 c        aml2=co(1)+x*co(2)+y*co(3)+x*y*co(4)
 c        amr2=co(5)+x*co(6)+y*co(7)+x*y*co(8)
 c        aml=sign(sqrt(abs(aml2)),aml2)
 c        amr=sign(sqrt(abs(amr2)),amr2)
 c        !------segment left of current breakpoint icub -----
 c        pleft(5)=aml
 c        do j=1,4
 c          pleft(j)=pst(j,nob+1)-x*pst(j,nob+2)+y*pst(j,nob+3)
 c        enddo
 c        call utlob4(-1,psg(1),psg(2),psg(3),psg(4),psg(5)
 c     $     ,pleft(1),pleft(2),pleft(3),pleft(4))
 c        !------segment right of current breakpoint icub-----
 c        pright(5)=amr
 c        do j=1,4
 c          pright(j)=pst(j,nob+4)+x*pst(j,nob+5)-y*pst(j,nob+6)
 c        enddo
 c        call utlob4(-1,psg(1),psg(2),psg(3),psg(4),psg(5)
 c     $     ,pright(1),pright(2),pright(3),pright(4))
 c        !-------------------------
 c        amt=pleft(5)**2+pleft(1)**2+pleft(2)**2
 c        if(amt.gt.0..and.pleft(4)+abs(pleft(3)).gt.0.d0)then
 c          amt=sqrt(amt)
 c          yleft=sign(1.,pleft(3))*alog((pleft(4)+abs(pleft(3)))/amt)
 c        else
 c          yleft=0.                  !
 c        endif
 c        amt=pright(5)**2+pright(1)**2+pright(2)**2
 c        if(amt.gt.0..and.pright(4)+abs(pright(3)).gt.0.d0)then
 c          amt=sqrt(amt)
 c          yright=sign(1.,pright(3))*alog((pright(4)+abs(pright(3)))/amt)
 c        else
 c          yright=0.                  !
 c        endif
 c        ybrk=(yleft+yright)/2.
 c        yhax=2                  !0.5*yhaha
 c        zzipx=1
 c        if(ybrk.gt.yhax)then
 c          zzipx=zzipx+(zzipp-1)
 c        elseif(ybrk.gt.-yhax)then
 c          zzipx=zzipx+(zzipp-1)*(ybrk+yhax)/(2*yhax)
 c        endif
 c        if(ybrk.lt.-yhax)then
 c          zzipx=zzipx+(zzipt-1)
 c        elseif(ybrk.lt.yhax)then
 c          zzipx=zzipx+(zzipt-1)*(yhax-ybrk)/(2*yhax)
 c        endif
 c      endif
       delptfra=0.
 c      if(ifl1.eq.3.and.ifl2.eq.0)delptfra=delptfra+ptfrasr
 c      if(ifl1.eq.3.or.ifl2.eq.3)delptfra=delptfra+ptfrasr
       fnsbp=1
       if(nsbp.gt.9)fnsbp=0
 c pt kink due to split elastic scattering (if zzzz.ne.1.)
       if(iLHC.eq.1)then
         pttra=(ptfra+delptfra)*(zzzz-1.)*fnsbp
         if(ig.ge.0)then         !quark and soft strings
 c        if(ig.eq.1)then         !quark and soft strings
 c        if(krm.eq.0)then         !quark and soft strings
           pt=ranpt()*(ptfra +  pttra)!/float(jqu)
         else
           pt=ranpt()*(ptfraqq   +  pttra)!/float(jqu)
         endif
       else
         pttra=(ptfra+delptfra)*(zzzz-1.)*fnsbp
         if(jqu.eq.1)then
           pt=ranpt()*(ptfra                +  pttra)
         else
           pt=ranpt()*(ptfraqq              +  pttra)
           pt=pt*0.5
         endif
       endif
 
       beta=2.*pi*rangen()
 
       if(ish.ge.5)then
         write(ifch,*) 'pt:',pt
       endif
       ptb(1,ib1+1)=pt*cos(beta)
       ptb(2,ib1+1)=pt*sin(beta)
       ptb(3,ib1+1)=0.
       ptb(4,ib1+1)=0.
       if(ish.ge.8)then
         write(ifch,*) 'the bands'
         write(ifch,'(4g12.6)') (pst(i,ii),i=1,4)
         write(ifch,'(4g12.6)') (pst(i,jj),i=1,4)
         write(ifch,*) 'pt before rotation'
         write(ifch,'(4f12.8)') (ptb(i,ib1+1),i=1,4)
       endif
       ax=dble(pst(1,ii))+dble(pst(1,jj))
       ay=dble(pst(2,ii))+dble(pst(2,jj))
       az=dble(pst(3,ii))+dble(pst(3,jj))
       ae=dble(pst(4,ii))+dble(pst(4,jj))
       am=sqrt(max(1d-10,(ae+az)*(ae-az)-ax**2-ay**2)) !?????????
       if(ish.ge.8)then
         write(ifch,*) 'boost vector ( region ) '
         write(ifch,'(5g12.6)') ax,ay,az,ae,am,pf(ii,jj)
       endif
       bx=pst(1,ii)
       by=pst(2,ii)
       bz=pst(3,ii)
       be=pst(4,ii)
       call utlob2(1,ax,ay,az,ae,am,bx,by,bz,be,60)
       if(ish.ge.8) then
         write (ifch,*) 'boost of b'
         write (ifch,*) 'bx,by,bz,be',bx,by,bz,be
       endif
       if(abs(bx)+abs(by)+abs(bz).gt.0.)then
         call utrot4(-1,bx,by,bz,ptb(1,ib1+1),ptb(2,ib1+1),ptb(3,ib1+1))
       else
         write(ifmt,*) 'null rot of pt',bx,by,bz
         write(ifmt,'(4f12.8)') (ptb(i,ib1+1),i=1,4)
       endif
       if(ish.ge.8) then
         write (ifch,*) 'rot of pt'
         write(ifch,'(4f12.8)') (ptb(i,ib1+1),i=1,4)
       endif
 
       call utlob4(-1,ax,ay,az,ae,am
      &  ,ptb(1,ib1+1),ptb(2,ib1+1),ptb(3,ib1+1),ptb(4,ib1+1))
 
       if(ish.ge.8) then
         write (ifch,*) 'backboost of pt'
         write(ifch,'(4f12.8)') (ptb(i,ib1+1),i=1,4)
       endif
 c      if(az.eq.0..and.ay.eq.0.)az=1e-8    !not needed if following line commented
 c      if(ish.ge.8)write(ifch,*)'rot vector:',ax,ay,az,ae,am
 
 c.....call utrota(-1,sngl(ax),sngl(ay),sngl(az)  !already commented in nexus2
 c....&              ,ptb(1,ib1+1),ptb(2,ib1+1),ptb(3,ib1+1))
 
       if(ish.ge.6)then
         write(ifch,*) 'pt'
         write(ifch,'(4g12.6)') (ptb(i,ib1+1),i=1,4)
         write (ifch,*) ijb(1,ib1+1),ijb(2,ib1+1),xyb(1,ib1+1)
      $       ,xyb(2,ib1+1)
       endif
 
       if(iside.eq.1)then
          ib1=ib1+1
          ib2=ib2+1
       endif
 
 
 c      Amin=0.
 c      if(Amax.lt.Amxx) goto 15
 c      goto 25
 
  15   continue
       if(ish.ge.6)write(ifch,*) 'Amax:',Amax,Amaxn
       if (nbr.eq.nbrold) then
          Amax=Amaxn
          Amin=0.
          if(pb*Amax.le.0..or.ntry.ge.1000)then
            goto 9999
          endif
          goto 26
       endif
 
       if(ish.ge.6)then
          write(ifch,*) 0,iflb(0)
          do i=1,nbr
            if(i.eq.ib2) write(ifch,*) '.................'
             write(ifch,'(i3,2x,2(i3),2(" ",g14.7),3x,i5,4(" ",g12.6))')
      &           i,ijb(1,i),ijb(2,i),xyb(1,i),xyb(2,i)
      &           ,iflb(i),(ptb(j,i),j=1,4)
            if(i.eq.ibr1) write(ifch,*) '.................'
          enddo
          write(ifch,*) nbr+1,iflb(nbr+1)
       endif
 
  9999 if(nbr.eq.nbrold)iok=1
       call utprix('gaksbp',ish,ishini,5)
       return
       end
 
 c----------------------------------------------------------------------
       function ranptcut(xcut)
 c----------------------------------------------------------------------
 c .........exp(-x**2)
 c      if(xcut.gt.0.)then
 c        x=sqrt(-(log(rangen()))/(3.1415927/4.)/xcut) !gauss
 c      else
 c        x=sqrt(-(log(rangen()))/(3.1415927/4.)) !gauss
 c      endif
  12   x=sqrt(-(log(rangen()))/(3.1415927/4.)) !gauss
 
 
       if(xcut.gt.0.)then
         if(rangen().lt.x/xcut)goto 12
       endif
 
       ranptcut=x
 
       return
 
 c .........exp(-x)
 c  12  xmx=50
 c      r=2.
 c      do while (r.gt.1.)
 c  11    x=sqrt(exp(rangen()*log(1+xmx**2))-1)
 c        if(x.eq.0.)goto11
 c        r=rangen()  /  ( exp(-x)*(1+x**2) )
 c      enddo
 c      x=x/2.
 
       end
 
 cc----------------------------------------------------------------------
 c      function ranpticut(xcut)
 cc----------------------------------------------------------------------
 c
 cc .........exp(-x)
 c
 c  12  xmx=50
 c      r=2.
 c      do while (r.gt.1.)
 c  11    x=sqrt(exp(rangen()*log(1+xmx**2))-1)
 c        if(x.eq.0.)goto11
 c        r=rangen()  /  ( exp(-x)*(1+x**2) )
 c      enddo
 c      x=x/2.
 c      if(rangen().gt.xcut/x)goto 12
 c
 c      ranpticut=x
 c
 c      end
 
 c----------------------------------------------------------------------
       function ranpt()
 c----------------------------------------------------------------------
 
 c .........exp(-x)
       xmx=50
       r=2.
       do while (r.gt.1.)
   11    x=sqrt(exp(rangen()*log(1+xmx**2))-1)
         if(x.eq.0.)goto11
         r=rangen()  /  ( exp(-x)*(1+x**2) )
       enddo
       ranpte=x/2.
 
 c     -------------
       ranpt=ranpte
       return
 c     -------------
 cc .........exp(-x**2)
 c      ranptg=sqrt(-log(rangen())/(3.1415927/4.)) !gauss
 c
 cc .........exp(-sqrt(x))
 c      xmx=500
 c      r=2.
 c      do while (r.gt.1.)
 c        x=sqrt(exp(rangen()*log(1+xmx**2))-1)
 c        r=rangen()  /  ( exp(-sqrt(x))*(1+x**2)/5. )
 c      enddo
 c      ranpts=x/20.
 c
 
       end
 c----------------------------------------------------------------------
       function ranptk()
 c----------------------------------------------------------------------
 
 c .........exp(-x)
       xmx=50
       r=2.
       do while (r.gt.1.)
   11    x=sqrt(exp(rangen()*log(1+xmx**2))-1)
         if(x.eq.0.)goto11
         r=rangen()  /  ( exp(-x)*(1+x**2) )
       enddo
       ranpte=x/2.
 
 c     -------------
       ranptk=ranpte
       return
 c     -------------
 c
 cc .........exp(-x**2)
 c      ranptg=sqrt(-log(rangen())/(3.1415927/4.)) !gauss
 c
 cc .........exp(-sqrt(x))
 c      xmx=500
 c      r=2.
 c      do while (r.gt.1.)
 c        x=sqrt(exp(rangen()*log(1+xmx**2))-1)
 c        r=rangen()  /  ( exp(-sqrt(x))*(1+x**2)/5. )
 c      enddo
 c      ranpts=x/20.
 
       end
 
 c----------------------------------------------------------------------
       function ranptd()
 c----------------------------------------------------------------------
 
 c .........exp(-x**2)
       ranptg=sqrt(-log(rangen())/(3.1415927/4.)) !gauss
 
 
 c     -------------
       ranptd=ranptg
       return
 c     -------------
 c
 cc .........exp(-x)
 c      xmx=50
 c      r=2.
 c      do while (r.gt.1.)
 c  11    x=sqrt(exp(rangen()*log(1+xmx**2))-1)
 c        if(x.eq.0.)goto11
 c        r=rangen()  /  ( exp(-x)*(1+x**2) )
 c      enddo
 c      ranpte=x/2.
 c
 cc .........exp(-sqrt(x))
 c      xmx=500
 c      r=2.
 c      do while (r.gt.1.)
 c        x=sqrt(exp(rangen()*log(1+xmx**2))-1)
 c        r=rangen()  /  ( exp(-sqrt(x))*(1+x**2)/5. )
 c      enddo
 c      ranpts=x/20.
 
 
 
 
 
       end
 c----------------------------------------------------------------------
       subroutine gakdel(ibr)
 c----------------------------------------------------------------------
       parameter (mxnbr=500,mxnba=5000)
       common /gag/nob,pst(4,0:mxnba),ipst(0:mxnba)
      $     ,nbr,ijb(2,0:mxnbr),xyb(2,0:mxnbr)
      &     ,ptb(4,0:mxnbr),iflb(0:mxnbr),ip(0:mxnbr),co,x,y
       dimension co(12)
 
       do i=ibr,nbr+1
          do j=1,2
             ijb(j,i)=ijb(j,i+1)
             xyb(j,i)=xyb(j,i+1)
          enddo
          do k=1,4
             ptb(k,i)=ptb(k,i+1)
          enddo
          iflb(i)=iflb(i+1)
          ip(i)=ip(i+1)
       enddo
       ip(ibr)=0
       nbr=nbr-1
       end
 
 c----------------------------------------------------------------------
       subroutine gaksco(ibr1,ibr,ibr2,ib,jb,x1,x2,y1,y2)
 c----------------------------------------------------------------------
       include 'epos.inc'
 
       parameter (mxnbr=500,mxnba=5000)
       common /gag/nob,pst(4,0:mxnba),ipst(0:mxnba)
      $     ,nbr,ijb(2,0:mxnbr),xyb(2,0:mxnbr)
      &     ,ptb(4,0:mxnbr),iflb(0:mxnbr),ip(0:mxnbr),co,x,y
       dimension co(12)
       logical weiter
 
       pf(i,j)=pst(4,i)*pst(4,j)-pst(3,i)*pst(3,j)
      &     -pst(2,i)*pst(2,j)-pst(1,i)*pst(1,j)
       mmod(i)=mod(mod(i,nob)+nob,nob)
 
       call utpri('gaksco',ish,ishini,8)
 
       if(ish.ge.8)then
          write(ifch,*) 'zwischen brk:',ibr1,ibr,ibr2,'(',nob,')',nbr
          write(ifch,*) 'region:',ib,jb
       endif
       if(ib.eq.ijb(1,ibr1))then
          x2=xyb(1,ibr1)
       else
          x2=1.
       endif
       if(ib.eq.ijb(1,ibr2))then
          x1=xyb(1,ibr2)
       else
          x1=0.
       endif
       if(jb.eq.ijb(2,ibr1))then
          y1=xyb(2,ibr1)
       else
          y1=0.
       endif
       if(jb.eq.ijb(2,ibr2))then
          y2=xyb(2,ibr2)
       else
          y2=1.
       endif
 
 c.....left side...................................................
       n=nob+1
       if(ish.ge.8)write(ifch,*)'x1,x2',x1,x2
       do i=1,4
 cc        pst(i,n)=(x2-x1)*pst(i,ib)+ptb(i,ibr)-ptb(i,ibr1)
         pst(i,n)=     x2*pst(i,ib)+ptb(i,ibr)-ptb(i,ibr1)-y1*pst(i,jb)
       enddo
       if(ish.ge.8)write(ifch,*) 'add a  1 ii',1.,ib
       ii=mmod(ib-1)
       weiter=.true.
       if(ib.eq.ijb(1,ibr1))weiter=.false.
       do while(ii.ne.jb.and.weiter) !linker Rand??
          f1=1.
          if(ii.eq.ijb(1,ibr1))f1=xyb(1,ibr1)
          do i=1,4
             pst(i,n)=pst(i,n)+f1*pst(i,ii)
          enddo
          if(ish.ge.8)write(ifch,*) 'add a f1 ii',f1,ii
          if(ii.eq.ijb(1,ibr1))weiter=.false.
          ii=mmod(ii-1)
       enddo
       jj=mmod(jb+1)
       weiter=.not.weiter
       if(jb.eq.ijb(2,ibr1))weiter=.false.
       do while(weiter)
          f1=1.
          if(jj.eq.ijb(2,ibr1))f1=1.-xyb(2,ibr1)
          do i=1,4
             pst(i,n)=pst(i,n)+f1*pst(i,jj)
          enddo
          if(ish.ge.8)write(ifch,*) 'add b f1 ii',f1,jj
          if(jj.eq.ijb(2,ibr1))weiter=.false.
          jj=mmod(jj+1)
       enddo
       do i=1,4
 cc        pst(i,n+1)=(x2-x1)*pst(i,ib)
         pst(i,n+1)=        pst(i,ib)
 cc        pst(i,n+2)=(y2-y1)*pst(i,jb)
         pst(i,n+2)=        pst(i,jb)
       enddo
       co(1)= pf(n,n)
       co(2)=-2.*pf(n,n+1)
       co(3)= 2.*pf(n,n+2)
       co(4)=-2.*pf(n+1,n+2)
       if(ish.ge.8) then
         do i=n,n+2
           write (ifch,'(4g12.5)') (pst(j,i),j=1,4)
         enddo
       endif
       if(ish.ge.8)write(ifch,*) 'co left:',co(1),co(2),co(3),co(4)
 
 c.....right side...................................................
       n=nob+4
       if(ish.ge.8)write(ifch,*)'y1,y2',y1,y2
       do i=1,4
 cc         pst(i,n)=(y2-y1)*pst(i,jb)-ptb(i,ibr)+ptb(i,ibr2)
          pst(i,n)=    y2*pst(i,jb)-ptb(i,ibr)+ptb(i,ibr2)-x1*pst(i,ib)
       enddo
       if(ish.ge.8)write(ifch,*) 'add a  1 jj',1.,jb
       ii=mmod(ib+1)
       weiter=.true.
       if(ib.eq.ijb(1,ibr2))weiter=.false.
       do while(ii.ne.jb.and.weiter)
          f1=1.
          if(ii.eq.ijb(1,ibr2))f1=1.-xyb(1,ibr2)
          do i=1,4
             pst(i,n)=pst(i,n)+f1*pst(i,ii)
          enddo
          if(ish.ge.8)write(ifch,*) 'add a f1 ii',f1,ii
          if(ii.eq.ijb(1,ibr2))weiter=.false.
          ii=mmod(ii+1)
       enddo
       jj=mmod(jb-1)
       weiter=.not.weiter
       if(jb.eq.ijb(2,ibr2))weiter=.false.
       do while(weiter)
          f1=1.
          if(jj.eq.ijb(2,ibr2))f1=xyb(2,ibr2)
          do i=1,4
             pst(i,n)=pst(i,n)+f1*pst(i,jj)
          enddo
          if(ish.ge.8)write(ifch,*) 'add b f1 ii',f1,jj
          if(jj.eq.ijb(2,ibr2))weiter=.false.
          jj=mmod(jj-1)
       enddo
       do i=1,4
 cc         pst(i,n+1)=(x2-x1)*pst(i,ib)
          pst(i,n+1)=        pst(i,ib)
 cc         pst(i,n+2)=(y2-y1)*pst(i,jb)
          pst(i,n+2)=         pst(i,jb)
       enddo
       co(5)=pf(n,n)
       co(6)= 2.*pf(n,n+1)
       co(7)=-2.*pf(n,n+2)
       co(8)=-2.*pf(n+1,n+2)
       if(ish.ge.8) then
         do i=n,n+2
           write (ifch,'(4g12.5)') (pst(j,i),j=1,4)
         enddo
       endif
       if(ish.ge.8)write(ifch,*) 'co right:',co(5),co(6),co(7),co(8)
 
 c.....lambda (absolute past).....................................
       n=nob+7
       do i=1,4
 cc         pst(i,n)= x1*pst(i,ib)+y1*pst(i,jb)
          pst(i,n)= 0.
       enddo
       ii=mmod(ib+1)
       do while (mmod(ii+jb).ne.0)
          if(ish.ge.8)write(ifch,*) 'add lambda',ii
          do i=1,4
             pst(i,n)=pst(i,n)+pst(i,ii)
          enddo
          ii=mmod(ii+1)
       enddo
       do i=1,4
 cc         pst(i,n+1)=(x2-x1)*pst(i,ib)
 cc         pst(i,n+2)=(y2-y1)*pst(i,jb)
          pst(i,n+1)= pst(i,ib)
          pst(i,n+2)= pst(i,jb)
       enddo
       co(9)=     pf(n,n)
       co(10)= 2.*pf(n,n+1)
       co(11)= 2.*pf(n,n+2)
       co(12)= 2.*pf(n+1,n+2)
       if(ish.ge.8)write(ifch,*)'co lambda:',co(9),co(10),co(11),co(12)
       call utprix('gaksco',ish,ishini,8)
       end
 
 c---------------------------------------------------------------------
       subroutine gakanp(ibr1,ibr,ibrr2,aml2,amr2,ala2,iok)
 c---------------------------------------------------------------------
 c   mass adjustment of fragments
 c        ibr1-ibr   ibr-ibrr2
 c   where ibr1,ibr,ibrr2 are break point indices
 c   aml2,amr2 are the reqired squared masses (if zero -> any mass)
 c   iok=0 (ok) or ok=1 (error)
 c---------------------------------------------------------------------
       include 'epos.inc'
       parameter (mxnbr=500,mxnba=5000,mxnin=2000)
       common /gag/nob,pst(4,0:mxnba),ipst(0:mxnba)
      $     ,nbr,ijb(2,0:mxnbr),xyb(2,0:mxnbr)
      &     ,ptb(4,0:mxnbr),iflb(0:mxnbr),ip(0:mxnbr),co,x,y
       double precision ax,ay,az,ae,am,bx,by,bz,be,A,B,C
       dimension co(12),am2(0:2),nin(2,0:mxnin)
       logical weiter
 c      pf(i,j)=pst(4,i)*pst(4,j)-pst(3,i)*pst(3,j)
 c     &     -pst(2,i)*pst(2,j)-pst(1,i)*pst(1,j)
       mmod(i)=mod(mod(i,nob)+nob,nob)
       call utpri('gakanp',ish,ishini,6)
 
       ibr2=ibrr2
       if(ish.ge.6)write(ifch,*) ibr1,ibr,ibr2,aml2,amr2,ala2,iok
       iok=0
       ib=ijb(1,ibr)
       jb=ijb(2,ibr)
       ni=0
  10   do i=1,ni
         if((nin(1,i).eq.ib.and.nin(2,i).eq.jb)
      $       .or.(ipst(ib).eq.ipst(jb)))then
           iok=1
           if(ish.ge.4)then
             write(ifch,*) 'error ... endless loop'
             if(ib.eq.ipst(jb))  write(ifch,*) ' in zero mass region'
           endif
           goto 9999
         endif
       enddo
       ni=ni+1
       if(ni.gt.mxnin)stop'gakanp: increase parameter mxnin  '
       nin(1,ni)=ib
       nin(2,ni)=jb
       if(ish.ge.6)write(ifch,*)
       if(ish.ge.6)write(ifch,*) 'ib,jb=',ib,jb
       if(ni.ge.2)then
          if(ish.ge.6)write(ifch,*)'rotate pt to new band'
          if(ish.ge.6)write(ifch,*)'from',nin(1,ni-1),nin(2,ni-1)
          if(ish.ge.6)write(ifch,*)'  to',ib,jb
          if(ish.ge.6)write(ifch,'(4f12.8)') (ptb(i,ibr),i=1,4)
          ax=pst(1,nin(1,ni-1))+pst(1,nin(2,ni-1))
          ay=pst(2,nin(1,ni-1))+pst(2,nin(2,ni-1))
          az=pst(3,nin(1,ni-1))+pst(3,nin(2,ni-1))
          ae=pst(4,nin(1,ni-1))+pst(4,nin(2,ni-1))
          am=sngl(sqrt(max(1d-8,dble(ae)**2-dble(ax)**2
      $        -dble(ay)**2-dble(az)**2))) !???????????????????????
          bx=pst(1,nin(1,ni-1))
          by=pst(2,nin(1,ni-1))
          bz=pst(3,nin(1,ni-1))
          be=pst(4,nin(1,ni-1))
          if(ish.ge.6)write (ifch,*) 'bx,by,bz,be',bx,by,bz,be
          call utlob2(1,ax,ay,az,ae,am,bx,by,bz,be,60)
          if(ish.ge.6)write (ifch,*) 'bx,by,bz,be',bx,by,bz,be
          call utlob4( 1,ax,ay,az,ae,am
      &        ,ptb(1,ibr),ptb(2,ibr),ptb(3,ibr),ptb(4,ibr))
          if(ish.ge.6)write(ifch,'(4f12.8)') (ptb(i,ibr),i=1,4)
          if(abs(bx)+abs(by)+abs(bz).gt.0.)then
            call utrot4( 1,bx,by,bz,ptb(1,ibr),ptb(2,ibr),ptb(3,ibr))
          else
            write(ifmt,*) 'null rot of pt (2)',bx,by,bz
            write(ifmt,'(4f12.8)') (ptb(i,ibr),i=1,4)
          endif
          if(ish.ge.6)write(ifch,'(4f12.8)') (ptb(i,ibr),i=1,4)
          ax=pst(1,ib)+pst(1,jb)
          ay=pst(2,ib)+pst(2,jb)
          az=pst(3,ib)+pst(3,jb)
          ae=pst(4,ib)+pst(4,jb)
          am=sngl(sqrt(max(1d-8,dble(ae)**2-dble(ax)**2
      $        -dble(ay)**2-dble(az)**2))) !???????????????????????
          if(am.le.1.1e-4)then
            if(ish.ge.5)write(ifch,*)'error ... am<1.1e-4'
            iok=1
            goto 9999
          endif
          bx=pst(1,ib)
          by=pst(2,ib)
          bz=pst(3,ib)
          be=pst(4,ib)
          if(ish.ge.6)write (ifch,*) 'bx,by,bz,be',bx,by,bz,be
          if(ish.ge.6)write (ifch,*) 'ax,ay,az,ae',ax,ay,az,ae,am
          call utlob2(1,ax,ay,az,ae,am,bx,by,bz,be,60)
          if(ish.ge.6)write (ifch,*) 'bx,by,bz,be',bx,by,bz,be
          if(abs(bx)+abs(by)+abs(bz).gt.0.)then
            call utrot4(-1,bx,by,bz,ptb(1,ibr),ptb(2,ibr),ptb(3,ibr))
          else
            write(ifmt,*) 'null rot of pt (3)',bx,by,bz
            write(ifmt,'(4f12.8)') (ptb(i,ibr),i=1,4)
          endif
          if(ish.ge.6)write(ifch,'(4f12.8)') (ptb(i,ibr),i=1,4)
          call utlob4(-1,ax,ay,az,ae,am
      &        ,ptb(1,ibr),ptb(2,ibr),ptb(3,ibr),ptb(4,ibr))
          if(ish.ge.6)write(ifch,'(4f12.8)') (ptb(i,ibr),i=1,4)
       endif
       call gaksco(ibr1,ibr,ibr2,ib,jb,x1,x2,y1,y2)
 c      if(ni.eq.1)print *,'------------------------------------',2
 cc      x=(xyb(1,ibr)-x1)/(x2-x1)
 cc      y=(xyb(2,ibr)-y1)/(y2-y1)
       x=xyb(1,ibr)
       y=xyb(2,ibr)
       am2(0)=aml2
       am2(1)=amr2
       am2(2)=ala2
       if(ish.ge.6)write(ifch,*) ibr1,ibr,ibr2,aml2,amr2,ala2,iok
       if(amr2.le.0.)then
          l1=2
          l2=0
       elseif(aml2.le.0.)then
          l1=1
          l2=2
       elseif(ala2.le.0.)then
          l1=1
          l2=0
       else
          stop' not like this , please...'
       endif
       if(ish.ge.6.and.amr2.le.0)write(ifch,*) 'adjust: 1',l1,l2
       if(ish.ge.6.and.aml2.le.0)write(ifch,*) 'adjust: 2' ,l1,l2
       if(ish.ge.6.and.ala2.le.0)write(ifch,*) 'adjust: 3',l1,l2
       i=4*l1
       j=4*l2
       A = dble(co(i+4))*dble(co(j+3)) - dble(co(j+4))*dble(co(i+3))
       B = dble(co(i+4))*dble(co(j+1)) - dble(co(i+3))*dble(co(j+2))
      &  + dble(co(i+2))*dble(co(j+3)) - dble(co(i+1))*dble(co(j+4))
      &  - dble(am2(l2))*dble(co(i+4)) + dble(am2(l1))*dble(co(j+4))
       C = dble(co(i+2))*dble(co(j+1)) - dble(co(i+1))*dble(co(j+2))
      &  + dble(am2(l1))*dble(co(j+2)) - dble(am2(l2))*dble(co(i+2))
       if (ish.ge.7) then
          write(ifch,*) 'ABC,q ',A,B,C,B**2-4.*A*C
          if(abs(A).gt.0.d0)then
             write(ifch,*) sqrt(max(0d0,B**2-4d0*A*C))/2d0/A-B/A/2d0
             write(ifch,*) -sqrt(max(0d0,B**2-4d0*A*C))/2d0/A-B/A/2d0
          endif
       endif
       x=0.
       y=0.
       xx=0.
       yy=0.
       if(abs(A).gt.1.d-20.and.B*B-4.*A*C.ge.0.d0)then
         y=sngl(sqrt(max(0.d0,B**2-4.*A*C))/2.d0/A-B/A/2.d0)
         if(abs(co(i+2)+y*co(i+4)).gt.0.)
      &       x=(am2(l1)-co(i+1)-y*co(i+3))/(co(i+2)+y*co(i+4))
       elseif(abs(A).le.1.d-20.and.abs(B).gt.0.d0)then
         y=-sngl(C/B)
         if(abs(co(i+2)+y*co(i+4)).gt.0.)
      &       x=(am2(l1)-co(i+1)-y*co(i+3))/(co(i+2)+y*co(i+4))
       else
         if(ish.ge.5)write(ifch,*)'error ... no solution of quadr equ'
         iok=1
         goto 9999
       endif
       if(abs(A).gt.1.d-20.and.B**2-4.*A*C.ge.0.d0)then
         yy=sngl(-sqrt(max(0.d0,B**2-4.*A*C))/2.d0/A-B/A/2.d0)
         if(abs(co(i+2)+yy*co(i+4)).gt.0.)
      &       xx=(am2(l1)-co(i+1)-yy*co(i+3))/(co(i+2)+yy*co(i+4))
       elseif(abs(A).le.1.d-20.and.abs(B).gt.0.d0)then
         yy=-sngl(C/B)
         if(abs(co(i+2)+yy*co(i+4)).gt.0.)
      &       xx=(am2(l1)-co(i+1)-yy*co(i+3))/(co(i+2)+yy*co(i+4))
       else
          if(ish.ge.5)write(ifch,*)'error ... no solution (2) '
          iok=1
          goto 9999
       endif
       if(ish.ge.6)then
          write(ifch,*) x ,y ,(co(i+2)+ y*co(i+4)),' OK '
          write(ifch,*) xx,yy,(co(i+2)+yy*co(i+4)),' OK '
       endif
       weiter=.true.
  50   if(x.gt.x1.and.x.lt.x2.and.y.gt.y1.and.y.lt.y2)then
 cc         xyb(1,ibr)=x1+(x2-x1)*x
 cc         xyb(2,ibr)=y1+(y2-y1)*y
          xyb(1,ibr)=x
          xyb(2,ibr)=y
          ijb(1,ibr)=ib
          ijb(2,ibr)=jb
          e1=pst(4,nob+1)-x*pst(4,nob+2)+y*pst(4,nob+3)
          e2=pst(4,nob+4)+x*pst(4,nob+5)-y*pst(4,nob+6)
          if( e1.lt.0. .or. e2.lt.0. ) then
            if(ish.ge.5)write(ifch,*)'error ... e1<0 or e2<0'
            iok=1
            goto 9999
          endif
          !amal2=co(1)+co(2)*x+co(3)*y+co(4)*x*y
          !amar2=co(5)+co(6)*x+co(7)*y+co(8)*x*y
          if(ish.ge.6)then
            amal2=co(1)+co(2)*x+co(3)*y+co(4)*x*y
            amar2=co(5)+co(6)*x+co(7)*y+co(8)*x*y
            write(ifch,*) 'brkshift:',xyb(1,ibr),xyb(2,ibr),ib,jb
            write (ifch,*)'E:',e1
            write (ifch,*)'E:',e2
            write(ifch,'(2(a6,1g12.6))') 'aml:'
      &          ,sqrt(max(0.,amal2)),'amr:',sqrt(max(0.,amar2))
          endif
          i=ibr1+1
          do while(i.le.ibr-1)
             if((mmod(ijb(1,i)+nob/2) .lt. mmod(ijb(1,ibr)+nob/2)
      &           .or.(mmod(ijb(1,i)+nob/2) .eq. mmod(ijb(1,ibr)+nob/2)
      &           .and.(xyb(1,i).gt.xyb(1,ibr))))
      &           .and.
      &           (ijb(2,i) .gt. ijb(2,ibr)
      &           .or.(ijb(2,i) .eq. ijb(2,ibr)
      &           .and.xyb(2,i).lt.xyb(2,ibr)))) goto 150
             if(ish.ge.6) then
                write(ifch,*) 'away:'
      &          ,i,xyb(1,i),xyb(2,i),ijb(1,i),ijb(2,i)
             endif
             call gakdel(i)
             i=i-1
             ibr=ibr-1
             ibr2=ibr2-1
  150        i=i+1
          enddo
          i=ibr+1
          do while (i.le.ibr2-1)
             if((mmod(ijb(1,i)+nob/2) .gt. mmod(ijb(1,ibr)+nob/2)
      &           .or.(mmod(ijb(1,i)+nob/2) .eq. mmod(ijb(1,ibr)+nob/2)
      &           .and.(xyb(1,i).lt.xyb(1,ibr))))
      &           .and.
      &           (ijb(2,i) .lt. ijb(2,ibr)
      &           .or.(ijb(2,i) .eq. ijb(2,ibr)
      &           .and.xyb(2,i).gt.xyb(2,ibr)))) goto 160
             if(ish.ge.6) then
                write(ifch,*) 'away:'
      &          ,i,xyb(1,i),xyb(2,i),ijb(1,i),ijb(2,i)
             endif
             call gakdel(i)
             ibr2=ibr2-1
             i=i-1
  160        i=i+1
          enddo
          goto 9999
       else
         if(x.gt.x2
      &    .and.ib.ne.ijb(1,ibr1) !brk-begrenzung
      &    .and.mmod(ib-1).ne.jb !linker oder rechter Rand
      &    .and.mmod(ib-1+jb).ne.0)then !oben oder unten
           ib=mmod(ib-1)
           goto 10
         endif
         if(x.lt.x1
      &    .and.ib.ne.ijb(1,ibr2) !brk-begrenzung
      &    .and.mmod(ib+1).ne.jb !linker oder rechter Rand
      &    .and.mmod(ib+1+jb).ne.0)then !oben oder unten
           ib=mmod(ib+1)
           goto 10
         endif
         if(y.gt.y2
      &    .and.jb.ne.ijb(2,ibr2) !brk-begrenzung
      &    .and.mmod(jb-1).ne.ib !linker oder rechter Rand
      &    .and.mmod(jb-1+ib).ne.0)then !oben oder unten
           jb=mmod(jb-1)
           goto 10
         endif
         if(y.lt.y1
      &    .and.jb.ne.ijb(2,ibr1) !brk-begrenzung
      &    .and.mmod(jb+1).ne.ib !linker oder rechter Rand
      &    .and.mmod(jb+1+ib).ne.0)then !oben oder unten
           jb=mmod(jb+1)
           goto 10
         endif
         if(weiter)then
           weiter=.false.
           x=xx
           y=yy
           goto 50
         endif
         if(ish.ge.5)write(ifch,*)'error ... x,y not in allowed range'
         iok=1
       endif
  9999 if(amr2.eq.0.) ibrr2=ibr2
       call utprix('gakanp',ish,ishini,6)
       end
 
 cc----------------------------------------------------------------------
 c      subroutine gakstr(ifl)
 cc----------------------------------------------------------------------
 cc
 cc     calculates string-fragments by taking off pt of breakup
 cc     do with ifl=1   undo with ifl=-1
 cc
 cc----------------------------------------------------------------------
 c      include 'epos.inc'
 c      common /cpptr/ pptr(4,mxptl),ystr(mxptl)
 c
 c      do i=1,nptl
 c        if(istptl(i).eq.29)then
 c          nk1=ifrptl(1,i)
 c          nk2=ifrptl(2,i)
 c          do j=nk1,nk2
 c            if ((istptl(j).eq.0.or.istptl(j-1).eq.0).and.j.ne.nk1) then
 c              do k=1,4
 c                pptl(k,j)=pptl(k,j)+pptr(k,j-1)*ifl
 c              enddo
 c              !write(ifch,*)"left side back  to ",j,(pptr(k,j-1),k=1,4)
 c            endif
 c            if ((istptl(j).eq.0.or.istptl(j+1).eq.0).and.j.ne.nk2) then
 c              do k=1,4
 c                pptl(k,j)=pptl(k,j)-pptr(k,j)*ifl
 c              enddo
 c              !write(ifch,*)"right side back to ",j,(-pptr(k,j),k=1,4)
 c            endif
 c            if(ifl.eq.-1.and.istptl(j).eq.0)then
 c           e=pptl(1,j)**2+pptl(2,j)**2+pptl(3,j)**2
 c     &           +pptl(5,j)**2
 c              e=sqrt(e)
 c           !dif=abs(e-pptl(4,j))
 c           !if(dif.gt.0.01.and.dif/e.gt.0.1)print*,j,e,pptl(4,j)
 c              pptl(4,j)=e
 c            endif
 c          enddo
 c        endif
 c        if(istptl(i).eq.0)then
 c          if ( ifl.eq.1 ) then
 c            ystr(i)=sign(1.,pptl(3,i))*alog((pptl(4,i)+abs(pptl(3,i)))
 c     *           /sqrt(pptl(5,i)**2+pptl(1,i)**2+pptl(2,i)**2) )
 c          endif
 c        endif
 c      enddo
 c      end
 
 c----------------------------------------------------------------------
       subroutine gakli2(nn1,nn2)
 c----------------------------------------------------------------------
 
       include 'epos.inc'
       double precision pgampr,rgampr
       common/cgampr/pgampr(5),rgampr(4)
 c      double precision db1,db2,db3,db4,db5
       character label*8,idlabl*8
 c      db1=0d0
 c      db2=0d0
 c      db3=rgampr(4)
 c      db4=sqrt(rgampr(1)**2+rgampr(2)**2+rgampr(3)**2)
 c      db5=sqrt( db4**2-rgampr(4)**2)
       n1=nn1
       n2=nn2
       if (n1.eq.0) n1=1
       if (n2.eq.0) n2=nptl
       write (ifch,'(1a4,5a12,4a4,a10,2a4)')
      &'no.','px','py','pz','E','m','ior','jor','if1','if2'
      &,'id','ist','ity'
       do i=n1,n2
          if (idptl(i).lt.10000)then
             label='        '
             label=idlabl(idptl(i))
          else
             label='        '
          endif
          chrg=0.
          if(iabs(idptl(i)).le.9999)call idchrg(idptl(i),chrg)
          write (ifch,125) i,(pptl(j,i),j=1,5),iorptl(i),jorptl(i)
      &        ,ifrptl(1,i),ifrptl(2,i),idptl(i)
      $        ,chrg             !charge
      &        ,istptl(i),ityptl(i),label
       enddo
  125  format (1i4,5g18.10,4i6,1i10
      $     ,1f5.2               !charge
      $     ,2i4,'  ',A8
 c     $     ,7g12.4,i5
      $     )
       return
       end
 
 c----------------------------------------------------------------------
 c      subroutine gakli4
 cc----------------------------------------------------------------------
 c
 c      include 'epos.inc'
 c      parameter (mxnbr=500,mxnba=5000)
 c      common /gag/nob,pst(4,0:mxnba),ipst(0:mxnba)
 c     $     ,nbr,ijb(2,0:mxnbr),xyb(2,0:mxnbr)
 c     &     ,ptb(4,0:mxnbr),iflb(0:mxnbr),ip(0:mxnbr),co,x,y
 c      dimension co(12)
 c      do i=0,nob-1
 c        write (ifch,10) i,(pst(j,i),j=1,4)
 c      enddo
 c 10   format(1i4,5g18.10)
 c      return
 c      end
 c
 cc----------------------------------------------------------------------
 c      subroutine gakli3
 cc----------------------------------------------------------------------
 c
 c      include 'epos.inc'
 c      parameter (mxnbr=500,mxnba=5000)
 c      common /gag/nob,pst(4,0:mxnba),ipst(0:mxnba)
 c     $     ,nbr,ijb(2,0:mxnbr),xyb(2,0:mxnbr)
 c     &     ,ptb(4,0:mxnbr),iflb(0:mxnbr),ip(0:mxnbr),co,x,y
 c      dimension co(12),p1(5),p2(5)
 c
 c      write(ifch,*) 'particle list of string decay'
 c      do i=1,5
 c        p1(i)=0.
 c      enddo
 c      do i=1,nbr+1
 c        if(i.lt.nbr+1)then
 c          call gaksco(i-1,i,i+1,ijb(1,i),ijb(2,i),x1,x2,y1,y2)
 c          if(x2.gt.x1)then
 ccc            x=(xyb(1,i)-x1)/(x2-x1)
 c            x=xyb(1,i)
 c          else
 c            x=0.
 c          endif
 c          if(y2.gt.y1)then
 ccc            y=(xyb(2,i)-y1)/(y2-y1)
 c            y=xyb(2,i)
 c          else
 c            y=0.
 c          endif
 c          aml2=co(1)+x*co(2)+y*co(3)+x*y*co(4)
 c          amr2=co(5)+x*co(6)+y*co(7)+x*y*co(8)
 c          aml=sign(sqrt(abs(aml2)),aml2)
 c          amr=sign(sqrt(abs(amr2)),amr2)
 c          do j=1,4
 c            p2(j)=pst(j,nob+1)-x*pst(j,nob+2)+y*pst(j,nob+3)
 c            p1(j)=p1(j)+p2(j)
 c          enddo
 c          p2(5)=aml
 c        else
 c          do j=1,4
 c            p2(j)=pst(j,nob+4)+x*pst(j,nob+5)-y*pst(j,nob+6)
 c            p1(j)=p1(j)+p2(j)
 c          enddo
 c          p2(5)=amr
 c        endif
 c        write(ifch,'(2i4,i6,a,i5,i10,5g14.6)') i-1,i
 c     &    ,-iflb(i-1),'==',iflb(i),ip(i)
 c     &    ,(p2(j),j=1,5)
 c      enddo
 c      am2=p1(4)**2-p1(3)**2-p1(2)**2-p1(1)**2
 c      p1(5)=sign(sqrt(abs(am2)),am2)
 c      write(ifch,'(12x,a60)')
 c     &  '------------------------------------------------------------'
 c      write(ifch,'(14x,5g14.6)') (p1(j),j=1,5)
 c      write(ifch,*)
 c
 c      end
 c
 
 c---------------------------------------------------------------------
       subroutine idress(id,am,idr,iadj)
 c---------------------------------------------------------------------
       include 'epos.inc'
       call idres(id,am,idr,iadj)
       if(idr.eq.0)then
         return
       endif
       ids=max(mod(iabs(id)/100,10),mod(iabs(id)/10,10))
       if(iabs(idr).le.999) then
 c        write (ifch,*) '  ',id,idr,ids
         if(ids.le.4)return  !???? if the following is used, bad result
         if(ids.le.2)then
           idr=sign(iabs(id)+int(rangen()+0.5),id)
         elseif(ids.eq.3)then
           idr=sign(iabs(id)+int(rangen()+0.6),id)
         else
           idr=sign(iabs(id)+int(rangen()+0.75),id)
         endif
 c        write (ifch,*) '->',id,idr
         call idmass(idr,am)
       elseif(iabs(idr).le.9999)then
         if(ids.le.3)return
         if(mod(iabs(idr),10).gt.1)then
           if(iabs(id).ne.1111.and.iabs(id).ne.2221.and.iabs(id).ne.3331)
      $         then
             idr=sign(iabs(id)+1,id)
             call idmass(idr,am)
           else
             idr=id
             call idmass(idr,am)
           endif
         endif
       endif
 
       end
 
 
 c---------------------------------------------------------------------
       SUBROUTINE gaksphe(sphe,r,mstu41)
 
 C...Purpose: to perform sphericity tensor analysis to give sphericity,
 C...aplanarity and the related event axes. stolen from jetset ;-)
       include 'epos.inc'
       dimension sphe(4,3)
       DIMENSION SM(3,3),SV(3,3)
 
 C...Calculate matrix to be diagonalized.
       NP=0
       JA=0
       JB=0
       JC=0
 c      MSTU41=2
       DO 110 J1=1,3
          DO 100 J2=J1,3
             SM(J1,J2)=0.
  100     CONTINUE
  110  CONTINUE
       PS=0.
       DO 140 I=1,nptl
       IF(istptl(i).ne.0)  GOTO 140
       IF(MSTU41.GE.2) THEN
          ida=iabs(idptl(i))
          IF(ida.EQ.0.OR.ida.EQ.11.OR.ida.EQ.13.OR.ida.EQ.15) GOTO 140
          IF(MSTU41.GE.3) then
             call idchrg(idptl(i),chrg)
             if (abs(chrg).le.0.1) goto 140
          endif
       ENDIF
       NP=NP+1
       PA=SQRT(pptl(1,i)**2+pptl(2,I)**2+pptl(3,i)**2)
       PWT=1.
       IF(ABS(r-2.).GT.0.001) PWT=MAX(1E-10,PA)**(r-2.)
       DO 130 J1=1,3
          DO 120 J2=J1,3
             SM(J1,J2)=SM(J1,J2)+PWT*pptl(j1,i)*pptl(j2,i)
  120     CONTINUE
  130  CONTINUE
       PS=PS+PWT*PA**2
  140  CONTINUE
 
 C...Very low multiplicities (0 or 1) not considered.
       IF(NP.LE.1) THEN
         if(ish.ge.1)then
           CALL utmsg('sphe  ')
           write(ifch,*) 'too few particles for analysis'
           call utmsgf
         endif
         sphe(4,1)=-1.
         RETURN
       ENDIF
       DO 160 J1=1,3
          DO 150 J2=J1,3
             SM(J1,J2)=SM(J1,J2)/PS
  150     CONTINUE
  160  CONTINUE
 
 C...Find eigenvalues to matrix (third degree equation).
       SQ=(SM(1,1)*SM(2,2)+SM(1,1)*SM(3,3)+SM(2,2)*SM(3,3)-SM(1,2)**2-
      &     SM(1,3)**2-SM(2,3)**2)/3.-1./9.
       SR=-0.5*(SQ+1./9.+SM(1,1)*SM(2,3)**2+SM(2,2)*SM(1,3)**2+SM(3,3)*
      &     SM(1,2)**2-SM(1,1)*SM(2,2)*SM(3,3))
      &     +SM(1,2)*SM(1,3)*SM(2,3)+1./27.
       SP=COS(ACOS(MAX(MIN(SR/SQRT(-SQ**3),1.),-1.))/3.)
       sphe(4,1)=1./3.+SQRT(-SQ)*MAX(2.*SP,SQRT(3.*(1.-SP**2))-SP)
       sphe(4,3)=1./3.+SQRT(-SQ)*MIN(2.*SP,-SQRT(3.*(1.-SP**2))-SP)
       sphe(4,2)=1.-sphe(4,1)-sphe(4,3)
       IF(sphe(4,2).LT.1E-5) THEN
         if(ish.ge.1)then
           CALL utmsg('gaksphe')
           write(ifch,*) 'all particles back-to-back'
           call utmsgf
         endif
         sphe(4,1)=-1.
         RETURN
       ENDIF
 
 C...Find first and last eigenvector by solving equation system.
       DO 240 I=1,3,2
          DO 180 J1=1,3
             SV(J1,J1)=SM(J1,J1)-sphe(4,I)
             DO 170 J2=J1+1,3
                SV(J1,J2)=SM(J1,J2)
                SV(J2,J1)=SM(J1,J2)
  170        CONTINUE
  180     CONTINUE
          SMAX=0.
          DO 200 J1=1,3
             DO 190 J2=1,3
                IF(ABS(SV(J1,J2)).LE.SMAX) GOTO 190
                JA=J1
                JB=J2
                SMAX=ABS(SV(J1,J2))
  190        CONTINUE
  200     CONTINUE
          SMAX=0.
          DO 220 J3=JA+1,JA+2
             J1=J3-3*((J3-1)/3)
             RL=SV(J1,JB)/SV(JA,JB)
             DO 210 J2=1,3
                SV(J1,J2)=SV(J1,J2)-RL*SV(JA,J2)
                IF(ABS(SV(J1,J2)).LE.SMAX) GOTO 210
                JC=J1
                SMAX=ABS(SV(J1,J2))
  210        CONTINUE
  220     CONTINUE
          JB1=JB+1-3*(JB/3)
          JB2=JB+2-3*((JB+1)/3)
          sphe(JB1,I)=-SV(JC,JB2)
          sphe(jb2,I)=SV(JC,JB1)
          sphe(jb,I)=-(SV(JA,JB1)*sphe(jb1,I)+SV(JA,JB2)
      &        *sphe(jb2,I))/SV(JA,JB)
          PA=SQRT(sphe(1,I)**2+sphe(2,I)**2+sphe(3,I)**2)
          SGN=(-1.)**INT(rangen()+0.5)
          DO 230 J=1,3
             sphe(j,I)=SGN*sphe(j,I)/PA
  230     CONTINUE
  240  CONTINUE
 
 C...Middle axis orthogonal to other two. Fill other codes.
       SGN=(-1.)**INT(rangen()+0.5)
       sphe(1,2)=SGN*(sphe(2,1)*sphe(3,3)-sphe(3,1)*sphe(2,3))
       sphe(2,2)=SGN*(sphe(3,1)*sphe(1,3)-sphe(1,1)*sphe(3,3))
       sphe(3,2)=SGN*(sphe(1,1)*sphe(2,3)-sphe(2,1)*sphe(1,3))
 
       do i=1,3
          do j=1,4
             pptl(j,nptl+i)=sphe(j,i)
          enddo
       enddo
 
 
 C...Calculate sphericity and aplanarity. Select storing option.
 ccc      SPH=1.5*(sphe(4,2)+sphe(4,3))
 ccc      APL=1.5*sphe(4,3)
 
       RETURN
       END
 
 C*********************************************************************
 
       SUBROUTINE gakthru(thru,mstu41)
 
 C...Purpose: to perform thrust analysis to give thrust, oblateness
 C...and the related event axes. stolen from jetset ;-)
       include 'epos.inc'
       DIMENSION TDI(3),TPR(3)
       dimension thru(4,3)
 
 C...Take copy of particles that are to be considered in thrust analysis.
       IAGR=0
       NP=0
       PS=0.
 c      MSTU41=2
       MSTU44=4
       MSTU45=2
       PARU42=1.0
       PARU48=0.0000001
       DO 100 I=1,nptl
          IF(istptl(i).ne.0)goto 100
          ida=iabs(idptl(i))
          IF(ida.EQ.0.OR.ida.EQ.11.OR.ida.EQ.13.OR.ida.EQ.15)GOTO 100
          IF(MSTU41.GE.3) then
             call idchrg(idptl(i),chrg)
             if (abs(chrg).le.0.1) goto 100
          endif
 
          IF(nptl+NP.GE.mxptl) THEN
             CALL utstop('gakthru: no more memory left in cptl&')
             thru(4,1)=-1.
             RETURN
          ENDIF
          NP=NP+1
 c         K(nptl+NP,1)=23
          pptl(1,nptl+NP)=pptl(1,I)
          pptl(2,nptl+NP)=pptl(2,I)
          pptl(3,nptl+NP)=pptl(3,I)
          pptl(4,nptl+NP)=SQRT(pptl(1,I)**2+pptl(2,I)**2+pptl(3,I)**2)
          pptl(5,nptl+NP)=1.
          IF(ABS(PARU42-1.).GT.0.001)
      &        pptl(5,nptl+NP)=pptl(4,nptl+NP)**(PARU42-1.)
          PS=PS+pptl(4,nptl+NP)*pptl(5,nptl+NP)
  100  CONTINUE
 
 C...Very low multiplicities (0 or 1) not considered.
       IF(NP.LE.1) THEN
          CALL utmsg('thru  ')
          write(ifch,*) 'too few particles for analysis'
          call utmsgf
          thru(4,1)=-1
          RETURN
       ENDIF
 
 C...Loop over thrust and major. T axis along z direction in latter case.
       DO 320 ILD=1,2
       IF(ILD.EQ.2) THEN
 c        PHI=GAKANG(pptl(1,nptl+NP+1),pptl(2,nptl+NP+1))
 c        CALL lurot(nptl+1,nptl+NP+1,0.,-PHI)
 c        THE=GAKANG(pptl(3,nptl+NP+1),pptl(1,nptl+NP+1))
 c        CALL lurot(nptl+1,nptl+NP+1,-THE,0.)
          ax=pptl(1,nptl+NP+1)
          ay=pptl(2,nptl+NP+1)
          az=pptl(3,nptl+NP+1)
          do irot=nptl+1,nptl+NP+1
             call utrota(1,ax,ay,az
      &           ,pptl(1,irot),pptl(2,irot),pptl(3,irot))
          enddo
          if(np.eq.2)then
            pptl(1,nptl+NP+2)=1.
            pptl(2,nptl+NP+2)=0.
            pptl(3,nptl+NP+2)=0.
            pptl(4,nptl+NP+2)=0.
            goto 325
          endif
       ENDIF
 
 C...Find and order particles with highest p (pT for major).
       DO 110 ILF=nptl+NP+4,nptl+NP+MSTU44+4
          pptl(4,ILF)=0.
  110  CONTINUE
       DO 160 I=nptl+1,nptl+NP
          IF(ILD.EQ.2) pptl(4,I)=SQRT(pptl(1,I)**2+pptl(2,I)**2)
          DO 130 ILF=nptl+NP+MSTU44+3,nptl+NP+4,-1
             IF(pptl(4,I).LE.pptl(4,ILF)) GOTO 140
             DO 120 J=1,5
                pptl(j,ILF+1)=pptl(j,ILF)
  120        CONTINUE
  130     CONTINUE
          ILF=nptl+NP+3
  140     DO 150 J=1,5
             pptl(j,ILF+1)=pptl(j,I)
  150     CONTINUE
  160  CONTINUE
 
 C...Find and order initial axes with highest thrust (major).
       DO 170 ILG=nptl+NP+MSTU44+5,nptl+NP+MSTU44+15
          pptl(4,ILG)=0.
  170  CONTINUE
       NC=2**(MIN(MSTU44,NP)-1)
       DO 250 ILC=1,NC
          DO 180 J=1,3
             TDI(J)=0.
  180     CONTINUE
          DO 200 ILF=1,MIN(MSTU44,NP)
             SGN=pptl(5,nptl+NP+ILF+3)
             IF(2**ILF*((ILC+2**(ILF-1)-1)/2**ILF).GE.ILC) SGN=-SGN
             DO 190 J=1,4-ILD
                TDI(J)=TDI(J)+SGN*pptl(j,nptl+NP+ILF+3)
  190        CONTINUE
  200     CONTINUE
          TDS=TDI(1)**2+TDI(2)**2+TDI(3)**2
          DO 220 ILG=nptl+NP+MSTU44+MIN(ILC,10)+4,nptl+NP+MSTU44+5,-1
             IF(TDS.LE.pptl(4,ILG)) GOTO 230
             DO 210 J=1,4
                pptl(j,ILG+1)=pptl(j,ILG)
  210        CONTINUE
  220     CONTINUE
          ILG=nptl+NP+MSTU44+4
  230     DO 240 J=1,3
             pptl(j,ILG+1)=TDI(J)
  240     CONTINUE
          pptl(4,ILG+1)=TDS
  250  CONTINUE
 
 C...  Iterate direction of axis until stable maximum.
       pptl(4,nptl+NP+ILD)=0.
       ILG=0
  260  ILG=ILG+1
       THP=0.
  270  THPS=THP
       DO 280 J=1,3
          IF(THP.LE.1E-10) TDI(J)=pptl(j,nptl+NP+MSTU44+4+ILG)
          IF(THP.GT.1E-10) TDI(J)=TPR(J)
          TPR(J)=0.
  280  CONTINUE
       DO 300 I=nptl+1,nptl+NP
          SGN=SIGN(pptl(5,I),TDI(1)*pptl(1,I)
      &        +TDI(2)*pptl(2,I)+TDI(3)*pptl(3,I))
          DO 290 J=1,4-ILD
             TPR(J)=TPR(J)+SGN*pptl(j,I)
  290     CONTINUE
  300  CONTINUE
       THP=SQRT(TPR(1)**2+TPR(2)**2+TPR(3)**2)/PS
       IF(THP.GE.THPS+PARU48) GOTO 270
 
 C...  Save good axis. Try new initial axis until a number of tries agree.
       IF(THP.LT.pptl(4,nptl+NP+ILD)-PARU48.AND.ILG.LT.MIN(10,NC))
      &     GOTO 260
       IF(THP.GT.pptl(4,nptl+NP+ILD)+PARU48)
      $     THEN
          IAGR=0
          SGN=(-1.)**INT(rangen()+0.5)
          DO 310 J=1,3
             pptl(j,nptl+NP+ILD)=SGN*TPR(J)/(PS*THP)
  310     CONTINUE
          pptl(4,nptl+NP+ILD)=THP
          pptl(5,nptl+NP+ILD)=0.
       ENDIF
       IAGR=IAGR+1
       IF(IAGR.LT.MSTU45.AND.ILG.LT.MIN(10,NC)) GOTO 260
  320  CONTINUE
 
 C...  Find minor axis and value by orthogonality.
  325  SGN=(-1.)**INT(rangen()+0.5)
       pptl(1,nptl+NP+3)=-SGN*pptl(2,nptl+NP+2)
       pptl(2,nptl+NP+3)=SGN*pptl(1,nptl+NP+2)
       pptl(3,nptl+NP+3)=0.
       THP=0.
       DO 330 I=nptl+1,nptl+NP
          THP=THP+pptl(5,I)*ABS(pptl(1,nptl+NP+3)*pptl(1,I)
      &        +pptl(2,nptl+NP+3)*pptl(2,I))
  330  CONTINUE
       pptl(4,nptl+NP+3)=THP/PS
       pptl(5,nptl+NP+3)=0.
 
 
 C...  Fill axis information. Rotate back to original coordinate system.
       do irot=nptl+NP+1,nptl+NP+3
          call utrota(-1,ax,ay,az
      &        ,pptl(1,irot),pptl(2,irot),pptl(3,irot))
       enddo
 
       do ild=1,3
          do j=1,4
             thru(j,ild)=pptl(j,nptl+NP+ild)
          enddo
       enddo
 
 C...Calculate thrust and oblateness. Select storing option.
 ccc      THR=thru(4,1)
 ccc      OBL=thru(4,2)-thru(4,3)
 
       RETURN
       END
 
       subroutine gakjet(ijadu)
       include 'epos.inc'
       common/njjjj/njets(5,2,mxbins)
 c      nmin=xpar1
 c      nmax=xpar2
       if(nrevt.eq.1)then
         do i=1,5
           do j=1,mxbins
             njets(i,ijadu,j)=0
           enddo
         enddo
       endif
       do i=1,nrbins
         ycut=xminim*(xmaxim/xminim)**((float(i)-0.5)/nrbins)
 c        if(iologe.ne.1)ycut=xminim+(xmaxim-xminim)/nrbins*(nrbins)
         nj=gaknjt(ycut,ijadu)
         if(nj.ge.1.and.nj.le.5)then
           njets(nj,ijadu,i)=njets(nj,ijadu,i)+1
         endif
       enddo
       end
 
       subroutine gakjto
       include 'epos.inc'
       common/njjjj/njets(5,2,mxbins)
       n=xpar4
       ijadu=xpar3
       do i=1,nrbins
         ycut=xminim*(xmaxim/xminim)**((float(i)-0.5)/nrbins)
         write (ifhi,*) ycut,float(njets(n,ijadu,i))/nrevt
      $       ,sqrt(1.*njets(n,ijadu,i))/nrevt
       enddo
       end
 C*********************************************************************
       function gaknjt(ycut,ijadu)
 c
 c     ijadu   1 =  JADE    2=DURHAM
 c     ycut - max. distance
 c
       include 'epos.inc'
 
       a2j(i,j)=2.*pptl(4,i)*pptl(4,j)*(1.-(pptl(1,i)*pptl(1,j)
      &     +pptl(2,i)*pptl(2,j)+pptl(3,i)*pptl(3,j))
      &     /(sqrt(pptl(1,i)**2+pptl(2,i)**2+pptl(3,i)**2)
      &       *sqrt(pptl(1,j)**2+pptl(2,j)**2+pptl(3,j)**2)))/evis**2
 
       a2d(i,j)=2.*min(pptl(4,i)**2,pptl(4,j)**2)
      &     *(1.-(pptl(1,i)*pptl(1,j)
      &     +pptl(2,i)*pptl(2,j)+pptl(3,i)*pptl(3,j))
      &     /(sqrt(pptl(1,i)**2+pptl(2,i)**2+pptl(3,i)**2)
      &       *sqrt(pptl(1,j)**2+pptl(2,j)**2+pptl(3,j)**2)))/evis**2
 
       bet(i)=sqrt(pptl(1,i)**2+pptl(2,i)**2+pptl(3,i)**2)
 
       ska(i,j)=pptl(1,i)*pptl(1,j)+pptl(2,i)*pptl(2,j)
      &     +pptl(3,i)*pptl(3,j)
 
       a2c(i,j)= ((bet(i)*bet(j)-ska(i,j))
      &     *2.*bet(i)*bet(j)  /  (0.00001+bet(i)+bet(j))**2 )
 
       evis=0.
       nn=0
       do i=1,nptl
          if (istptl(i).eq.0) then
             nn=nn+1
             do j=1,5
                pptl(j,nptl+nn)=pptl(j,i)
             enddo
             evis=evis+pptl(4,i)
             jorptl(i)=nn
          endif
       enddo
       iflag=0
       i1=0
       j1=0
       do while (iflag.eq.0.and.nn.ge.2)
          a2min=ycut
          iflag=1
          do i=nptl+1,nptl+nn-1
             do j=i+1,nptl+nn
                if(ijadu.eq.1)then
                   a2=a2j(i,j)
                elseif(ijadu.eq.2) then
                   a2=a2d(i,j)
                else
                   a2=a2c(i,j)
                endif
                if (a2.lt.a2min) then
                   iflag=0
                   i1=i
                   j1=j
                   a2min=a2
                endif
             enddo
          enddo
          if (iflag.eq.0) then
             do i=1,4
                pptl(i,i1)=pptl(i,i1)+pptl(i,j1)
             enddo
             do i=1,nptl
                if(istptl(i).eq.0.and.jorptl(i).eq.j1-nptl)
      &              jorptl(i)=i1-nptl
                if(istptl(i).eq.0.and.jorptl(i)+nptl.gt.j1)
      &              jorptl(i)=jorptl(i)-1
             enddo
             do i=j1,nptl+nn
                do j=1,5
                   pptl(j,i)=pptl(j,i+1)
                enddo
             enddo
             nn=nn-1
          endif
       enddo
       do i=nptl+1,nptl+nn
          istptl(i)=-1
          jorptl(i)=i-nptl
          pptl(5,i)=sqrt(max(0.,(pptl(4,i)+pptl(3,i))
      &   *(pptl(4,i)-pptl(3,i))-pptl(2,i)**2-pptl(1,i)**2))
       enddo
       do i=nptl+1,nptl+nn-1
          do j=i+1,nptl+nn
             if(pptl(4,jorptl(i)+nptl).lt.pptl(4,jorptl(j)+nptl))then
                k=jorptl(i)
                jorptl(i)=jorptl(j)
                jorptl(j)=k
             endif
          enddo
       enddo
       do i=nptl+1,nptl+nn
          idptl(nptl+jorptl(i))=9910+i-nptl
       enddo
       do i=1,nptl
         jorptl(i)=0             !jorptl(nptl+jorptl(i))
       enddo
 c      nptl=nptl+nn
 
       gaknjt=nn
       return
       end
 
       subroutine idtr5(id,ic)
       integer ic(2)
       ic(1)=0
       ic(2)=0
       ii=1
       if(id.lt.0)ii=2
       i1=1
       if(iabs(id).gt.999)i1=3
       do i=i1,int(log(abs(float(id)))/log(10.))+1
         j=mod(iabs(id)/10**(i-1),10)
         if(j.gt.0)then
           ic(ii)=ic(ii)+10**(6-j)
         endif
       enddo
       return
       end
 
       function ammin(id1,id2)
       dimension ic1(2),ic2(2),jc2(6,2),jc1(6,2)
       call idtr5(id1,ic1)
       call idtr5(id2,ic2)
       call idcomk(ic1)
       call idcomk(ic2)
       call iddeco(ic1,jc1)
       call iddeco(ic2,jc2)
       ammin=utamnx(jc1,jc2)
       end
 
 
 c      function idtr(id1,id2)
 c      integer ic(2),id(2)
 c      id(1)=id1
 c      id(2)=id2
 c      do i=1,2
 c        ic(i)=0
 c      enddo
 c      do j=1,2
 c        ii=1
 c        if(id(j).lt.0)ii=2
 c        i1=1
 c        if(iabs(id(j)).gt.999)i1=3
 c        do i=i1,int(log(abs(float(id(j))))/log(10.))+1
 c          jj=mod(iabs(id(j))/10**(i-1),10)
 c          if(jj.gt.0)then
 c            ic(ii)=ic(ii)+10**(6-jj)
 c          endif
 c        enddo
 c      enddo
 c      idtr=idtra(ic,0,0,3)
 c      if(idtr.ne.idsp(id1,id2))then
 c        write (*,'(4i6)') idtr,idsp(id1,id2),id1,id2
 c      endif
 c      return
 c      end
 
       function idsp(id1,id2)
       ia1=iabs(id1)
       ia2=iabs(id2)
       if(ia1.ge.1000.and.ia2.ge.1000)then
         idsp=0
         isign=0
       elseif(ia1.le.1000.and.ia2.le.1000)then
         idsp=min(ia1,ia2)*100+max(ia1,ia2)*10
         isign=1
         if(max(ia1,ia2).ne.-min(id1,id2)) isign = -1
         if(idsp.eq.220)idsp=110
         if(idsp.eq.330)idsp=220
       else
         isign=1
         if(id1.lt.0.and.id2.lt.0)isign=-1
         idb=min(ia1,ia2)
         if(idb.eq.5)then
           idsp=0
           return
         endif
         ida=max(ia1,ia2)
         ida1=ida/1000
         ida2=mod(ida/100,10)
         if(idb.le.ida1)then
           idsp=idb*1000+ida/10
         elseif(idb.le.ida2)then
           idsp=ida1*1000+idb*100+ida2*10
         else
           idsp=ida+idb*10
         endif
         if(ida1.eq.ida2.and.ida2.eq.idb)idsp=idsp+1
       endif
       idsp=idsp*isign
       return
       end