Back to home page

Project CMSSW displayed by LXR

 
 

    


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

0001 c-----------------------------------------------------------------------
0002       subroutine gakfra(iclu,iret)
0003 c-----------------------------------------------------------------------
0004 
0005       include 'epos.inc'
0006       include 'epos.incems'
0007       include 'epos.incsem'
0008       include 'epos.incpar'
0009       parameter (eta=0.4,etap=0.1)
0010       parameter (mxnbr=500,mxnba=5000)
0011       common /gag/nob,pst(4,0:mxnba),ipst(0:mxnba)
0012      $     ,nbr,ijb(2,0:mxnbr),xyb(2,0:mxnbr)
0013      &     ,ptb(4,0:mxnbr),iflb(0:mxnbr),ip(0:mxnbr),co,x,y
0014 c      common /cpptr/ pptr(4,mxptl),ystr(mxptl)
0015       double precision p1,p12,p2
0016       dimension co(12),ind(0:mxnbr),p1(5),p2(5)
0017       dimension ic2(2),ic1(2),ic(2),fkap(3)
0018       common/pb/pb /cnsbp/nsbp  /cn8ptl/n8ptl
0019       common/czz/kky,krm,kdrm,kzrm,pui(3),pdi(3),psi(3),pci(3),zzzz,itp
0020      &,pduui(3),pdudi(3),pdusi(3),pduci(3),pdddi(3),pddsi(3),pddci(3)
0021      &,pdssi(3),pdsci(3),pdcci(3)
0022       double precision psg
0023       common/zpsg/psg(5)
0024       logical go
0025 
0026       pf(i,j)=(pst(4,i)-pst(3,i))*(pst(4,j)+pst(3,j))
0027      &     -pst(2,i)*pst(2,j)-pst(1,i)*pst(1,j)
0028 
0029       nob=0
0030       call utpri('gakfra',ish,ishini,4)
0031       iret=0
0032 c.....search string in pptl-list...................................
0033       nkmax=nptl
0034       nk1=maproj+matarg+1
0035  2    continue
0036       if(nclean.gt.0.and.nptl.gt.mxptl/2)then
0037         nnnpt=0
0038         do iii=maproj+matarg+1,nptl
0039           go=.true.
0040           if(nclean.eq.1.and.istptl(iii).le.istmax)go=.false.
0041           if(go.and.mod(istptl(iii),10).ne.0)then
0042             istptl(iii)=99
0043             nnnpt=nnnpt+1
0044           endif
0045         enddo
0046         if(nnnpt.gt.mxptl-nptl)then
0047           nptl0=nptl
0048           call utclea(maproj+matarg+1,nptl0)
0049           nkmax=nptl
0050           nk1=maproj+matarg+1
0051 c          write(ifch,'(a)')' ---------after utclea------------'
0052 c          call blist('&',nk1,nkmax)
0053         endif
0054       endif
0055 
0056       do while(nk1.le.nkmax)
0057         if(istptl(nk1).eq.20)then
0058           nk2=nk1+1
0059           do while(idptl(nk2).eq.9)
0060             nk2=nk2+1
0061           enddo
0062           goto 3                !ok, string from nk1 to nk2
0063         else
0064           nk1=nk1+1
0065         endif
0066       enddo
0067 
0068       goto 9999                 !no more string
0069 
0070 
0071 c.......................................................................
0072 c                       decay string  nk1 - nk2
0073 c.......................................................................
0074 
0075  3    nob=-1
0076       if(ish.ge.3)then
0077         write(ifch,'(/1x,25a1/5x,a/1x,25a1/)')
0078      *  ('-',k=1,25),'string decay',('-',k=1,25)
0079         write(ifch,'(a)')' ---------string------------'
0080         call blist('&',nk1,nk2)
0081       endif
0082 
0083       itp=ityptl(nk1)
0084       if(iLHC.eq.1)then
0085         if(jorptl(nk1).ne.0.or.jorptl(nk2).ne.0)then
0086           kky=1
0087         else
0088           kky=0
0089         endif
0090       else
0091         kky=0
0092 c      if(nk2-nk1.gt.1.or.iappl.eq.6)kky=1
0093         if((itp.ge.30.and.itp.le.39).or.iappl.eq.6)kky=1
0094       endif
0095       krm=0
0096       if((itp.ge.50.and.itp.le.59).or.(itp.ge.40.and.itp.le.49))krm=1
0097       kdrm=0
0098       if(itp.eq.43.or.itp.eq.53)kdrm=1
0099 c Excited remnant due to split
0100       kzrm=0
0101       if(itp.eq.44.or.itp.eq.54.or.itp.eq.46.or.itp.eq.56)kzrm=1
0102 
0103       do k=1,5
0104         p1(k)=0d0
0105       enddo
0106       do i=nk1,nk2
0107         do k=1,5
0108           p2(k)=dble(pptl(k,i))
0109         enddo
0110         p2(4)=sqrt(p2(3)**2+p2(2)**2+p2(1)**2+p2(5)**2)
0111         do k=1,4
0112           p1(k)=p1(k)+p2(k)
0113         enddo
0114       enddo
0115       p1(5)=(p1(4)-p1(3))*(p1(4)+p1(3))-p1(2)**2-p1(1)**2
0116       if(p1(5).gt.0.d0)then
0117         p1(5)=sqrt(p1(5))
0118       else
0119         iorrem=iorptl(nk1)
0120         write(*,*)'Precision problem in gakfra, p:',
0121      &             p1(5),p1(4),p1(3),p1(2),p1(1)
0122         write(*,*)'origin : ',iorrem
0123         p1(5)=0.d0
0124 c        if(iorrem.ne.0.and.
0125 c     &    (ityptl(iorrem).eq.40.or.ityptl(iorrem).eq.10))
0126 c     &    p1(5)=dble(pptl(5,iorrem))
0127         if(iorrem.ne.0)then
0128           write(*,*)'string mass : ',ityptl(iorrem),pptl(5,iorrem)
0129           p1(5)=dble(pptl(5,iorrem))
0130         endif
0131       endif
0132       do k=1,5
0133         psg(k)=p1(k)
0134       enddo
0135 
0136 c.....pbreak ................................
0137 
0138       pbi=pbreak
0139       zz=0.
0140       if(iLHC.eq.1)then
0141         x=p1(5)                 !sub string mass (energy in cms) for e+e- and pp
0142         if(pbi.gt.0.d0.or.(pbi.gt.-0.99999d0.and.kky.eq.0))then
0143 c around 0.4 for soft because of pi0 spectra for CR and pp multiplicity 
0144 c and pt spectra at RHIC
0145           pb=abs(pbi)
0146         else
0147 c fit is important not only in e+e- but for all particle production
0148           pb0=abs(pbi)
0149 
0150           pb=pb0+(1.-exp(-x/40.))*(pbreakg-pb0)
0151 
0152 c       {{14,0.25},{22,0.22},{34,0.22},{91.2,0.15}} !new
0153         endif
0154       else
0155         x0=1.
0156         x=sqrt(qsqptl(nk1))     !string energy
0157         if(iappl.ne.6)then  
0158 c in hadronic collisions, strings are substring of big string from a Pomeron
0159 c the ratio Energy_tot / Energy_sub is used to quantify the modification of fkappa 
0160           x0=x                  !pomeron energy
0161           x=p1(4)               !sub string energy
0162         endif
0163 c       if(pbi.gt.0.d0.or.(pbi.gt.-0.99999d0.and.krm.eq.1))then 
0164         if(pbi.gt.0.d0.or.(pbi.gt.-0.99999d0.and.kky.eq.0))then
0165 c around 0.4 for soft because of pi0 spectra for CR and pp multiplicity 
0166 c and pt spectra at RHIC
0167           pb=abs(pbi)
0168         else
0169 c fit is important not only in e+e- but for all particle production
0170           pb0=0.33
0171 
0172 c       pb=0.065 + 3600./x**3 - 400./x**2 + 17./x   !<--------
0173           pb=pb0*(1.+exp(-x**2/1000))*0.5
0174 
0175 c       {{14,0.25},{22,0.22},{34,0.22},{91.2,0.15}} !new
0176           if(iappl.ne.6)then
0177             zz=x0/x             !energy dependence from ratio
0178             if(iLHC.eq.1)zz=log(max(1.,zz))
0179           endif
0180         endif
0181       endif
0182 
0183       if(iLHC.eq.1)then
0184         fkap(1)=fkappa
0185         fkap(3)=fkappag
0186         fkap(2)=0.5*(fkap(1)+fkap(3))
0187       else
0188         fkap(1)=fkappa
0189         if(zz.gt.0.)then
0190           zzz=min(fkamax,fkainc*zz)
0191           zzz=sqrt(1.+zzz)
0192           if(iLHC.eq.0)pb=pb/zzz
0193           fkap(1)=fkap(1)*zzz         !fix increase of ap and K at Tevatron
0194         endif
0195 c      write(*,*)'pb',kky,krm,kdrm,x,x0,pb,zz,fkap
0196       endif
0197 
0198 c      write(*,*)'pb',kky,krm,kdrm,x,pb,zpaptl(1,nk1)+zpaptl(1,nk2)
0199 
0200       zzzz=1.
0201       if(isplit.eq.1.and.iappl.ne.6)then       !for pt
0202         zzzz=0.
0203         zzz=0.
0204 c        if(krm.eq.0)then  !both hard and soft (not good for mesons for HERA-B or pp RHIC, but good for baryons)
0205         if(kky.eq.1)then   !only hard
0206 c        if((kky.eq.1.or.(iLHC.eq.1.and.krm.eq.0)))then   
0207 c          zzzz=zipinc*(zpaptl(1,nk1)+zpaptl(1,nk2))
0208 c          zzzz=sqrt(1.+zzzz)
0209           zz=(zpaptl(1,nk1)*(1.+zodinc*log(zpaptl(2,nk2)))     !zpaptl(2,nk1)=# of connected pom
0210      &       +zpaptl(1,nk2)*(1.+zodinc*log(zpaptl(2,nk2))))
0211 c          print *,zpaptl(2,nk1),zpaptl(2,nk2)
0212 c     &           ,1./cosh(0.007*max(0.,zpaptl(2,nk1)-50.))
0213 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 
0214           if(iclu.eq.1)zz=zz/cosh(0.0085*max(0.,zpaptl(2,nk1)-30.))   
0215           zzzz=zipinc*zz
0216           if(iLHC.eq.1)then
0217             zzz=min(fkamax,fkainc*zz)
0218             zzz=sqrt(1.+zzz)
0219             fkap(1)=fkap(1)*zzz !fix increase of ap and K at Tevatron
0220             fkap(2)=fkap(2)*zzz 
0221             fkap(3)=fkap(3)*zzz
0222           endif
0223 c          print *,bimevt,zzzz,zz,0.75*(zpaptl(1,nk1)+zpaptl(1,nk2))
0224 c          zzz=fkainc*zz
0225         elseif(kzrm.eq.1)then
0226           zz=zpaptl(1,nk1)
0227           zzzz=0. !zopinc*zz
0228           zzz=0.!zodinc*zz
0229           zzz=sqrt(1.+zzz)
0230           fkap(1)=fkap(1)*zzz
0231           pb=pb/zzz
0232         endif
0233         zzzz=sqrt(1.+zzzz)
0234       endif
0235       if(ish.ge.4)write(ifch,*)"String parameters:"
0236      &                      ,iclu,itp,krm,kdrm,kzrm,fkap,zzzz,pb,x
0237       if(fkap(1).le.0..or.zzzz.le.0.)
0238      &   call utstop("fkappa<=0. not allowed !&")
0239 
0240       igmx=1
0241       if(iLHC.eq.1)igmx=3
0242       do ig=1,igmx
0243       pui(ig)=1.
0244       pdi(ig)=pui(ig)*exp(-pi*difud/fkap(ig))    !assymmetric u and d relevant because of assymmetry Kch / K0 in e+e-
0245       psi(ig)=pui(ig)*exp(-pi*difus/fkap(ig))
0246       pduui(ig)=pui(ig)*exp(-pi*difuuu/fkap(ig))
0247       pdudi(ig)=pui(ig)*exp(-pi*difuud/fkap(ig))
0248       pdusi(ig)=pui(ig)*exp(-pi*difuus/fkap(ig))
0249       pdddi(ig)=pui(ig)*exp(-pi*difudd/fkap(ig))
0250       pddsi(ig)=pui(ig)*exp(-pi*difuds/fkap(ig))
0251       pdssi(ig)=pui(ig)*exp(-pi*difuss/fkap(ig))
0252       if(nrflav.gt.3)then
0253         pci(ig)=pui(ig)*exp(-pi*difuc/fkap(ig))
0254         pduci(ig)=pui(ig)*exp(-pi*difuuc/fkap(ig))
0255         pddci(ig)=pui(ig)*exp(-pi*difudc/fkap(ig))
0256         pdsci(ig)=pui(ig)*exp(-pi*difusc/fkap(ig))
0257         pdcci(ig)=pui(ig)*exp(-pi*difucc/fkap(ig))
0258       else
0259         pci(ig)=0.
0260         pduci(ig)=0.
0261         pddci(ig)=0.
0262         pdsci(ig)=0.
0263         pdcci(ig)=0.
0264       endif
0265       enddo
0266 
0267 c.....light like ...............................................
0268 
0269       if(ish.ge.6) call blist("before light like&",nk1,nk2)
0270       if(rangen().lt.0.5)then
0271         i1=nk1
0272         i2=nk1+1
0273       else
0274         i1=nk2-1
0275         i2=nk2
0276       endif
0277       ii=1
0278       do while(ii.ge.0)
0279         do i=1,4
0280           p2(i)=dble(pptl(i,i1))+dble(pptl(i,i2))
0281         enddo
0282 c       p2(5)=sqrt(p2(4)**2-p2(3)**2-p2(2)**2-p2(1)**2)
0283         am1=pptl(5,i1)**2
0284         am2=pptl(5,i2)**2
0285         am12=max(0d0,p2(4)**2-p2(3)**2-p2(2)**2-p2(1)**2-am1-am2)/2
0286         if(am12**2.gt.am1*am2)then
0287           ak1=.5*((am2+am12)/sqrt(am12**2-am1*am2))-.5
0288           ak2=.5*((am1+am12)/sqrt(am12**2-am1*am2))-.5
0289         else
0290           ak1=0.
0291           ak2=0.
0292         endif
0293         if(ish.ge.6)write(ifch,'(a,2i4,9f12.6)') 'overlaps'
0294      $       ,i1,i2,ak1,ak2,sqrt(2.*am12)
0295      &       ,am1,am2
0296         do j=1,4
0297           a1=(1.+ak1)*pptl(j,i1)-ak2*pptl(j,i2)
0298           a2=(1.+ak2)*pptl(j,i2)-ak1*pptl(j,i1)
0299           pptl(j,i1)=a1
0300           pptl(j,i2)=a2
0301         enddo
0302         pptl(5,i1)=0.
0303         pptl(5,i2)=0.
0304         if(nk2-nk1.eq.1) ii=-1
0305         ii=ii-1
0306         if (nk1.eq.i1)then
0307           i1=nk2-1
0308           i2=nk2
0309         else
0310           i1=nk1
0311           i2=nk1+1
0312         endif
0313       enddo
0314 
0315 c.....inverse  ...............................................
0316 
0317 c      if(nk2.eq.nk1+1.and.nptl.lt.mxptl)then
0318 c        call utrepla(nptl+1,nk2)
0319 c        call utrepla(nk2,nk1)
0320 c        call utrepla(nk1,nptl+1)
0321 c      endif
0322 
0323       if(ish.ge.6) call blist("after  light like&",nk1,nk2)
0324 
0325 c.....on-shell and copy ...............................................
0326 
0327       if(ish.ge.6) write (ifch,*) 'on-shell check'
0328       do i=nk1,nk2
0329         do k=1,5
0330           p2(k)=dble(pptl(k,i))
0331         enddo
0332         p12=p2(4)
0333         p2(4)=sqrt(p2(3)**2+p2(2)**2+p2(1)**2+p2(5)**2)
0334         if(abs(p12-p2(4))/max(.1d0,p12,p2(4)).ge.1e-1.and.ish.ge.2)then
0335           write(ifch,*)'warning: on-shell problem:'
0336      &           ,i, idptl(i),(pptl(k,i),k=1,4),' (',sngl(p2(4)),') '
0337      &           ,pptl(5,i), istptl(i)
0338         endif
0339         if(ish.ge.6)  write (ifch,*)  p12 ,'->',p2(4)
0340         call utlob2(1,p1(1),p1(2),p1(3),p1(4),p1(5)
0341      $       ,p2(1),p2(2),p2(3),p2(4),60)
0342         f=1.0
0343         if(i.ne.nk1.and.i.ne.nk2)f=0.5
0344         nmax=1
0345         ff=1.
0346         aemax=1000.              ! here max band mass
0347 c          write(ifch,*)"test",f,p2(4),aemax,f*p2(4).gt.aemax
0348         if(f*p2(4).gt.aemax) then
0349           nmax = int(f*p2(4)/aemax)+1
0350           ff=1./float(nmax)
0351 c          print *,"nmax",nmax
0352 c          write(ifch,*)"nmax",nmax
0353         endif
0354         nn=1
0355         do while(nn.le.nmax)
0356           f=.5
0357           if(i.eq.nk1.and.nn.eq. 1  ) f=1.
0358           if(i.eq.nk2.and.nn.eq.nmax) f=1.
0359           nob=nob+1
0360           if(nob.gt.mxnba-7)stop'gakfra: increase parameter mxnba '
0361 c         if(nmax.ge.2) print *, nob,ff,f,i
0362 c         if(nmax.ge.2) write(ifch,*)'nob',nob,ff,f,i
0363           do k=1,4
0364             pst(k,nob)=ff*f*p2(k)
0365             ipst(nob)=i
0366           enddo
0367           nn=nn+1
0368         enddo
0369       enddo                     !i
0370 
0371       do i1=nob-1,1,-1          ! copy again gluons
0372         nob=nob+1
0373         if(nob.gt.mxnba-7)stop'gakfra: increase parameter mxnba '
0374         do k=1,4
0375           pst(k,nob)=pst(k,i1)
0376           ipst(nob)=ipst(i1)
0377         enddo
0378       enddo
0379       nob=nob+1                 ! nob is number of kinks
0380       if(nob.gt.mxnba-7)stop'gakfra: increase parameter mxnba '
0381       if(ish.ge.6) then
0382         write (ifch,*) 'bands:'
0383         do i=0,nob-1
0384           write (ifch,'(i4,4g20.13)') i,(pst(k,i),k=1,4)
0385         enddo
0386       endif
0387 
0388 c.....total momentum...............................................
0389       do k=1,4
0390         pst(k,nob)=0.
0391       enddo
0392       do i=0,nob-1
0393         do k=1,4
0394           pst(k,nob)=pst(k,nob)+pst(k,i)
0395         enddo
0396       enddo
0397 
0398 
0399 
0400 c.....write total string on pptl - list.............................
0401       nptl=nptl+1
0402       if(nptl.gt.mxptl)call utstop('gakfra: mxptl too small ...&')
0403       call idtr5(idptl(nk1),ic1)
0404       call idtr5(idptl(nk2),ic2)
0405       ic(1)=ic1(1)+ic2(1)
0406       ic(2)=ic1(2)+ic2(2)
0407       idptl(nptl) = 8*10**8+ ic(1)*100 + ic(2)/100 !nbr+1
0408       istptl(nptl)=10
0409       n8ptl=nptl
0410       do i=1,5
0411         pptl(i,nptl)=p1(i)
0412       enddo
0413       if(ish.ge.3)then
0414         write(ifch,'(a)')' ---------total string------------'
0415         call blist('&',nptl,nptl)
0416       endif
0417 
0418 c....................................................................
0419       ijb(1,0)=-1
0420       ijb(2,0)=-1
0421       ijb(1,1)=-1
0422       ijb(2,1)=-1
0423       iflb(0)=-idptl(nk1)
0424       iflb(1)= idptl(nk2)
0425       do i=1,4
0426         ptb(i,0)=0.
0427         ptb(i,1)=0.
0428       enddo
0429 
0430 c...................light string....................................
0431       amm=ammin(idptl(nk1),idptl(nk2))
0432 
0433       if(sqrt(max(0.,pf(nob,nob))).lt.amm)then
0434         id=idsp(  idptl(nk1),idptl(nk2) )
0435           if(ish.ge.1)then
0436        write (ifch,*)
0437      .       '-------string too light to decay--------'
0438           write (ifch,*) id,p1(5),amm
0439           write (ifch,*) id, sqrt(max(0.,pf(nob,nob))) ,amm
0440        endif
0441       if(id.ne.0.or.iLHC.ne.1)then
0442         am=sqrt(max(0.,pf(nob,nob)))
0443         call idress(id,am,idr,iadj)
0444        if(ish.ge.1)write (ifch,*) id,am,idr,iadj
0445         nptl=nptl+1
0446         if(nptl.gt.mxptl)
0447      &       call utstop('gakfra (light): mxptl too small ...&')
0448         do i=1,5
0449           pptl(i,nptl)=p1(i)
0450         enddo
0451         do i=nk1,nk2
0452           istptl(i)=21
0453           ifrptl(1,i)=n8ptl
0454           ifrptl(2,i)=0
0455         enddo
0456         istptl(n8ptl)=29        !code for fragmented string
0457         ityptl(n8ptl)=ityptl(nk1)
0458         itsptl(n8ptl)=itsptl(nk1)
0459         rinptl(n8ptl)=-9999
0460         iorptl(n8ptl)=nk1
0461         jorptl(n8ptl)=nk2
0462         ifrptl(1,n8ptl)=n8ptl+1
0463         ifrptl(2,n8ptl)=nptl
0464         ityptl(nptl)=ityptl(nk1)
0465         itsptl(nptl)=itsptl(nk1)
0466         rinptl(nptl)=-9999
0467         istptl(nptl)=0
0468         iorptl(nptl)=n8ptl
0469         jorptl(nptl)=0
0470         ifrptl(1,nptl)=0
0471         ifrptl(2,nptl)=0
0472         if(idr.ne.0)then
0473           idptl(nptl)=idr
0474         else
0475           idptl(nptl)=id
0476         endif
0477         do j=1,4
0478           xorptl(j,nptl)=xorptl(j,nk1)
0479         enddo
0480         tivptl(1,nptl)=xorptl(4,nptl)
0481         call idtau(idptl(nptl),pptl(4,nptl),pptl(5,nptl),taugm)
0482         tivptl(2,nptl)=tivptl(1,nptl)+taugm*(-alog(rangen()))
0483         if(abs(p1(5)-am).gt.0.01)then
0484 c          write (*,*) 'light string  ---  particle off-shell!!!!'
0485 c          write (*,*) idr,'  mass:',p1(5),'  should be:',am
0486         endif
0487         goto 98
0488       endif
0489       endif
0490 
0491 c................search breakpoints...................................
0492 
0493       n=0
0494       nsbp=1
0495       nbr=0
0496       iok=0
0497 
0498  11   continue
0499       if(nsbp.gt.10000)then
0500         if(ish.ge.1)then
0501         write(*,*)'Gakfra : string can not be fragmented, redo event !'
0502         write(*,*)nk1,idptl(nk1),nk2,idptl(nk2),istptl(nk1)
0503         write(ifch,*)'Gakfra : redo event !'
0504         endif
0505         iret=1
0506         goto 9999
0507       endif
0508 
0509       !----------------------!
0510       call gaksbp(0,1,iok)
0511       !----------------------!
0512       nbrtry=0
0513       goto 10
0514 
0515 c..............delete breakpoint....................................
0516  9    if(ish.ge.4)write(ifch,*) 'delete breakpoint',n,' -> '
0517      &,nbr-1,' breakpoints left'
0518       call gakdel(n)              !hier loeschen
0519 
0520 c..............no more breakpoints -> redo..........................
0521       if(nbr.eq.0)then
0522         nsbp=nsbp+1
0523         if(ish.ge.3)write (ifch,*)
0524      &    'no breakpoints left ... try again (',nsbp,')'
0525         if(ish.ge.3)write (ifch,*)' '
0526         goto 11
0527       endif
0528 
0529 c................make index list of breakpoints to adjust...........
0530  10   continue
0531       nbrtry=nbrtry+1
0532       nind=0
0533       do i=1,nbr
0534         if(ip(i).eq.0.or.ip(i+1).eq.0)then
0535           nind=nind+1
0536           ind(nind)=i
0537         endif
0538       enddo
0539       if(nbrtry.gt.1000000)then
0540         if(ish.ge.1)then
0541         write(*,*)'Gakfra : string can not be fragmented, redo event !'
0542         write(*,*)nk1,nk2,nbr,nind,pb
0543         endif
0544         iret=1
0545         goto 9999
0546       endif
0547 
0548 c.....no more breakpoint to adjust...............................
0549       if(nind.eq.0) goto 100
0550 
0551 c................................................................
0552       if(ish.ge.5)then
0553         write(ifch,*) 'breakpoints:'
0554         write(ifch,'(i3,i5)') 0,-iflb(0)
0555         do i=1,nbr
0556           write(ifch,'(i3,2i5,1x,4e12.4,2x,2i3,2f6.3)')
0557      &      i,iflb(i),-iflb(i),(ptb(j,i),j=1,4)
0558      &         ,ijb(1,i),ijb(2,i),xyb(1,i),xyb(2,i)
0559         enddo
0560         write(ifch,'(i3,i5)') nbr+1,iflb(nbr+1)
0561       endif
0562 
0563 c.....choose breakpoint, calculate masses, lambda..................
0564       r=rangen()
0565       nn=1+int(r*float(nind))
0566 c      nn=1+(nind-1)*int(r+0.5)
0567       n=ind(nn)
0568       if(ish.ge.4)write(ifch,*) 'choose breakpoint',n
0569       nl=n-1
0570       nr=n+1
0571       do while (nl.gt.0.and.ip(nl).eq.0)
0572          nl=nl-1
0573       enddo
0574       do while (nr.lt.nbr+1.and.ip(nr+1).eq.0)
0575          nr=nr+1
0576       enddo
0577       if(ish.ge.6)write (ifch,*) 'nl,n,nr:',nl,n,nr
0578 c      print *,'------------------------------------',1
0579       call gaksco(n-1,n,n+1,ijb(1,n),ijb(2,n),x1,x2,y1,y2)
0580       if(x2.le.x1.or.y2.le.y1)goto 9
0581 cc      x=(xyb(1,n)-x1)/(x2-x1)
0582 cc      y=(xyb(2,n)-y1)/(y2-y1)
0583       x=xyb(1,n)
0584       y=xyb(2,n)
0585       aml2=co(1)+x*co(2)+y*co(3)+x*y*co(4)
0586       amr2=co(5)+x*co(6)+y*co(7)+x*y*co(8)
0587       ala2=co(9)+x*co(10)+y*co(11)+x*y*co(12)
0588 
0589 c.....determine id of both particles..............................
0590       aml=sqrt(max(0.,aml2))
0591       idl=idsp( -iflb(n-1) , iflb(n)   )
0592       amr=sqrt(max(0.,amr2))
0593       idr=idsp( -iflb(n)   , iflb(n+1) )
0594 
0595 c.....if mass to small (because of spacelike pt)  reject..........
0596 c      amin=0.0
0597 c      if (aml2.le.amin.or.amr2.le.amin) goto 9
0598 
0599 c.....if no left or right ptl id -> reject.........
0600       if(idl.eq.0.or.idr.eq.0)then
0601         if(ish.ge.5)write(ifch,*)'no left or right ptl id'
0602         goto 9
0603       endif
0604 
0605       iadjl=0
0606       iadjr=0
0607       idrl=0
0608       idrr=0
0609 
0610 c.....if no adjusted mass on left side, check for resonance...........
0611       if(ip(n)  .eq.0) then
0612         aml=sqrt(max(0.,aml2+0.*min(0.,amr2))) !!!????
0613         amlxx=aml
0614         call idress(idl,aml,idrl,iadjl)
0615         r=rangen()
0616         if(idrl.eq.110.and.r.lt.0.5)then
0617           if (r.gt.eta+etap) goto 9      !rare numerical errors
0618           idl=220
0619           aml=.549
0620           if(r.lt.etap)aml=.958
0621           amlxx=aml
0622           call idress(idl,aml,idrl,iadjl)
0623           iadjl=1
0624         endif
0625         if(ish.ge.5)write(ifch,'(a,i5,2f12.6,1i10,i5)')
0626      &    ' l:  ',idl,amlxx,aml,idrl,iadjl
0627       else
0628         if(ish.ge.5)write(ifch,'(a,i5,2f12.6,i10,a5)')
0629      &    ' l:  ',idl,aml,aml,ip(n),'ok'
0630       endif
0631 
0632 c.....if no adjusted mass on right side, check for resonance...........
0633       if(ip(n+1).eq.0) then
0634         amr=sqrt(max(0.,amr2+0.*min(0.,aml2))) !!!????
0635         amrxx=amr
0636         call idress(idr,amr,idrr,iadjr)
0637         r=rangen()
0638         if(idrr.eq.110.and.r.lt.0.5)then
0639           if (r.gt.eta+etap) goto 9    !rare numerical errors
0640           idr=220
0641           amr=.549
0642           if(r.lt.etap)amr=.958
0643           amrxx=amr
0644           call idress(idr,amr,idrr,iadjr)
0645           iadjr=1
0646         endif
0647         if(ish.ge.5)write(ifch,'(a,i5,2f12.6,1i10,i5)')
0648      &    ' r:  ',idr,amrxx,amr,idrr,iadjr
0649       else
0650         if(ish.ge.5)write(ifch,'(a,i5,2f12.6,i10,a5)')
0651      &    ' r:  ',idr,amr,amr,ip(n+1),'ok'
0652       endif
0653 
0654 c.....mass adjustments.........................................
0655       iok=0
0656       if(ip(n+1).ne.0)then  !.........adjusted mass on right side
0657         if(idrl.eq.0)then
0658           call gaksbp(n-1,n,iok)
0659           if(iok.eq.1)goto 9
0660           goto 10
0661         endif
0662         if(iadjl.eq.1)then
0663            if(ish.ge.5)write(ifch,*)'mass adjustment 1'
0664            n2=n+1
0665            call gakanp(n-1,n,n2,aml**2,amr**2,0.,iok)
0666         endif
0667         if(iok.eq.1)goto 9
0668         ip(n)=idrl
0669       elseif(ip(n).ne.0)then !.........adjusted mass on left side
0670         if(idrr.eq.0)then
0671           call gaksbp(n,n+1,iok)
0672           if(iok.eq.1)goto 9
0673           goto 10
0674         endif
0675         if(iadjr.eq.1)then
0676            if(ish.ge.5)write(ifch,*)'mass adjustment 2'
0677            n2=n+1
0678            call gakanp(n-1,n,n2,aml**2,amr**2,0.,iok)
0679         endif
0680         if(iok.eq.1)goto 9
0681         ip(n+1)=idrr
0682       else       !.........adjusted mass neither left nor right
0683         if(idrr.eq.0.and.idrl.eq.0)then
0684           call gaksbp(n,n+1,iok)
0685           if(iok.eq.1)goto 9
0686           call gaksbp(n-1,n,iok)
0687           if(iok.eq.1)goto 9
0688           goto 10
0689         elseif(idrl.ne.0.and.idrr.eq.0)then
0690           if(iadjl.eq.1) then
0691            if(ish.ge.5)write(ifch,*)'mass adjustment 3'
0692            call gakanp(n-1,n,nr,aml**2,0.,ala2,iok)
0693           endif
0694         elseif(idrl.eq.0.and.idrr.ne.0)then
0695           if(iadjr.eq.1) then
0696            if(ish.ge.5)write(ifch,*)'mass adjustment 4'
0697            n2=n+1
0698            call gakanp(nl,n,n2,0.,amr**2,ala2,iok)
0699           endif
0700         else  !if(idrl.ne.0.and.idrr.ne.0)then
0701           if(iadjl.eq.1.or.iadjr.eq.1) then
0702            if(ish.ge.5)write(ifch,*)'mass adjustment 5'
0703            n2=n+1
0704            call gakanp(n-1,n,n2,aml**2,amr**2,0.,iok)
0705           endif
0706         endif
0707         if(iok.eq.1)goto 9
0708         ip(n)=idrl
0709         ip(n+1)=idrr
0710       endif
0711       if(ish.ge.4)then
0712         write(ifch,*) 'left/right particles:'
0713      &         ,ip(n),aml,'  ',ip(n+1),amr
0714       endif
0715       goto 10
0716 
0717 c................................................................
0718 c                         end of string decay
0719 c................................................................
0720 
0721 c.....final list...............................................
0722  100  if(ish.ge.4)then
0723         write(ifch,*) '   ************ OK **************'
0724         write(ifch,*) 'final breakpoints:'
0725         write(ifch,'(i3,i5)') 0,-iflb(0)
0726         do i=1,nbr
0727           write(ifch,'(i3,2i5,1x,4e12.4,2x,2i3,2f6.3)')
0728      &      i,iflb(i),-iflb(i),(ptb(j,i),j=1,4)
0729      &         ,ijb(1,i),ijb(2,i),xyb(1,i),xyb(2,i)
0730         enddo
0731         write(ifch,'(i3,i5)') nbr+1,iflb(nbr+1)
0732       endif
0733 
0734 c.....write particles in pptl-list................................
0735       if(ish.ge.3)then
0736         write(ifch,'(a)')' ---------produced particles---------'
0737       endif
0738 
0739       nptlini=nptl
0740       if(nptlini+nbr+1.gt.mxptl)
0741      &       call utstop('gakfra (end): mxptl too small ...&')
0742       do i=1,nbr+1
0743          nptl=nptl+1
0744          if(i.lt.nbr+1)then     !particle = left side of breakpoints
0745            call gaksco(i-1,i,i+1,ijb(1,i),ijb(2,i),x1,x2,y1,y2)
0746 c           taubrr=pst(4,nob+7)+x*pst(4,nob+8)+y*pst(4,nob+9)
0747            x=xyb(1,i)
0748            y=xyb(2,i)
0749            aml2=co(1)+x*co(2)+y*co(3)+x*y*co(4)
0750            amr2=co(5)+x*co(6)+y*co(7)+x*y*co(8)
0751 
0752            ala2=co(9)+x*co(10)+y*co(11)+x*y*co(12)
0753            aml=sign(sqrt(abs(aml2)),aml2)
0754            amr=sign(sqrt(abs(amr2)),amr2)
0755            if(aml.le.0.d0.or.amr.le.0.d0)then
0756              if(ish.ge.4)write(ifch,*)
0757      & 'Negative mass, fragmentation not OK -> redo ...'
0758              n=i
0759              nptl=nptlini
0760              goto 9
0761            endif
0762            qsqptl(nptl)=ala2
0763            pptl(5,nptl)=aml
0764            do j=1,4
0765              pptl(j,nptl)=pst(j,nob+1)-
0766      &            x*pst(j,nob+2)+y*pst(j,nob+3)
0767 c             pptr(j,nptl)=ptb(j,i)
0768            enddo
0769          else                   !last particle = right side of last breakpoint
0770            pptl(5,nptl)=amr
0771            qsqptl(nptl)=0.
0772            do j=1,4
0773              pptl(j,nptl)=pst(j,nob+4)+
0774      &            x*pst(j,nob+5)-y*pst(j,nob+6)
0775 c             pptr(j,nptl)=ptb(j,i) !should be zero
0776            enddo
0777 c           taubrr=0.
0778          endif
0779          idptl(nptl)=ip(i)
0780          if(ish.ge.7)call blist('&',nptl,nptl)
0781          if(pptl(4,nptl).le.0.d0)then
0782            if(ish.ge.4)write(ifch,*)
0783      & 'Negative energy, fragmentation not OK -> redo ...'
0784            n=i
0785            nptl=nptlini
0786            goto 9
0787          endif
0788 
0789       enddo
0790 
0791       nptl=nptlini
0792 
0793       if(ish.ge.7)then
0794         write(ifch,'(a)')' ---------produced particles---------'
0795       endif
0796 
0797 
0798       do i=1,nbr+1
0799          nptl=nptl+1
0800 
0801          call utlob4(-1,p1(1),p1(2),p1(3),p1(4),p1(5)
0802      $        ,pptl(1,nptl),pptl(2,nptl),pptl(3,nptl),pptl(4,nptl))
0803 c         call utlob4(-1,p1(1),p1(2),p1(3),p1(4),p1(5)
0804 c     $        ,pptr(1,nptl),pptr(2,nptl),pptr(3,nptl),pptr(4,nptl))
0805 
0806 
0807 c........Origin..................................................
0808          istptl(nptl)=0
0809          iorptl(nptl)=n8ptl
0810          jorptl(nptl)=0
0811          ifrptl(1,nptl)=0
0812          ifrptl(2,nptl)=0
0813 
0814          r=rangen()
0815          tauran=-taurea*alog(r)*pptl(4,nptl)/pptl(5,nptl)
0816 c         if(jpsi.ge.1)tauran=0.
0817 c         tauran=max(taubrl,taubrr) !take formation time from string-theory
0818          do j=1,4
0819            xorptl(j,nptl)=xorptl(j,nk1)+pptl(j,nptl)
0820      &       /pptl(4,nptl)*tauran
0821          enddo
0822          tivptl(1,nptl)=xorptl(4,nptl)
0823          call idtau(idptl(nptl),pptl(4,nptl),pptl(5,nptl),taugm)
0824          tivptl(2,nptl)=tivptl(1,nptl)+taugm*(-alog(rangen()))
0825          ityptl(nptl)=ityptl(nk1)
0826          itsptl(nptl)=itsptl(nk1)
0827          rinptl(nptl)=-9999
0828          if(ish.ge.3)call blist('&',nptl,nptl)
0829 ctp060829          taubrl=taubrr
0830       enddo
0831       iorptl(n8ptl)=nk1
0832       jorptl(n8ptl)=nk2
0833       ityptl(n8ptl)=ityptl(nk1)
0834       do i=nk1,nk2
0835          istptl(i)=21
0836          ifrptl(1,i)=n8ptl
0837          ifrptl(2,i)=0
0838       enddo
0839       istptl(n8ptl)=29          !code for fragmented string
0840       ifrptl(1,n8ptl)=n8ptl+1
0841       ifrptl(2,n8ptl)=nptl
0842       if(ish.ge.5)then
0843         write(ifch,*)'string momentum sum:'
0844         do kk=1,5
0845           pptl(kk,nptl+1)=0
0846           do ii=n8ptl+1,nptl
0847             pptl(kk,nptl+1)=pptl(kk,nptl+1)+pptl(kk,ii)
0848           enddo
0849         enddo
0850         call alist2('&',n8ptl,n8ptl,nptl+1,nptl+1)
0851       endif
0852 
0853 
0854 
0855 c.....another string?.........................................
0856  98   nk1=nk2+1
0857       goto 2                    !next string
0858 
0859 c.....end of fragmentation.....................................
0860  9999 continue
0861       call utprix('gakfra',ish,ishini,4)
0862       return
0863       end
0864 
0865 c---------------------------------------------------------------------
0866       subroutine gaksbp(ibr1,ibr2,iok)
0867 c---------------------------------------------------------------------
0868       ! search break points
0869       !-----------------------------------------
0870       ! nbr ... number of break points
0871       !
0872       !
0873       !
0874       !--------------------------------------------------------------
0875       include 'epos.inc'
0876 
0877       parameter (mxnbr=500,mxnba=5000)
0878       common /gag/nob,pst(4,0:mxnba),ipst(0:mxnba)
0879      $     ,nbr,ijb(2,0:mxnbr),xyb(2,0:mxnbr)
0880      &     ,ptb(4,0:mxnbr),iflb(0:mxnbr),ip(0:mxnbr),co,x,y
0881       common/pb/pb /cnsbp/nsbp /cn8ptl/n8ptl
0882       double precision ax,ay,az,ae,am,bx,by,bz,be
0883       dimension co(12)
0884       logical weiter
0885       common/czz/kky,krm,kdrm,kzrm,pui(3),pdi(3),psi(3),pci(3),zzzz,itp
0886      &,pduui(3),pdudi(3),pdusi(3),pduci(3),pdddi(3),pddsi(3),pddci(3)
0887      &,pdssi(3),pdsci(3),pdcci(3)
0888       double precision psg
0889       common/zpsg/psg(5)
0890 c      dimension pleft(5),pright(5)
0891 
0892       pf(i,j)=pst(4,i)*pst(4,j)-pst(3,i)*pst(3,j)
0893      &     -pst(2,i)*pst(2,j)-pst(1,i)*pst(1,j)
0894       mmod(i)=mod(mod(i,nob)+nob,nob)
0895 
0896       call utpri('gaksbp',ish,ishini,5)
0897 
0898 
0899 
0900       ib1=ibr1
0901       ib2=ibr2
0902       iside=1
0903       if(rangen().lt.0.5)iside=2
0904       if(ish.ge.6)write(ifch,*)'iside:',iside
0905       Amxx=80./pb
0906       if(ish.ge.4)write(ifch,*)
0907      &'search brk between ib1=',ib1,' and ib2=',ib2
0908       ntry=0
0909       nbrold=nbr
0910       Amin=0.
0911       Amax=Amxx
0912  26   ntry=ntry+1
0913       A0=-log(exp(-pb*Amin)+rangen()*(exp(-pb*Amax)-exp(-pb*Amin)))/pb
0914       if(ish.ge.6)write(ifch,*)'pb, Amin, Amax, A0:',pb, Amin, Amax, A0
0915       A=A0
0916       Amaxn=0.
0917       if(iside.eq.2)goto 51
0918 c.....................................................................
0919       if(ib1.eq.0)then          !startwert der j-schleife
0920          jj=nob-1
0921       else
0922          jj=ijb(2,ib1)
0923       endif
0924       do while (jj.gt.0)
0925          if(jj.eq.ijb(2,ib1) )then       !linker brk im Streifen jj?
0926             y1=xyb(2,ib1)
0927          else                   !nein
0928             y1=0.
0929          endif
0930          if(jj.eq.ijb(2,ib2))then       !rechter brk im Streifen jj?
0931             y2=xyb(2,ib2)
0932          else                   !nein
0933             y2=1.
0934          endif
0935          if(y1.eq.y2) goto 9999
0936          if(abs(y1-y2)/abs(y1+y2).le.1e-7) goto 9999
0937          if(ish.ge.6)write(ifch,*)'y1,y2,A',y1,y2,A
0938                                 !startwert der i-schleife
0939          if(ib1.eq.0)then       !ohne linken brkpt
0940             ii=mmod(jj+1)
0941             if(jj.lt.nob/2)ii=mmod(nob-jj+1)
0942             if(ib2.le.nbr.and.mmod(ii+nob/2)
0943      &        .gt.mmod(ijb(1,ib2)+nob/2))then
0944               if(ish.ge.6)write(ifch,*) 'very special case',ii,jj
0945               goto12
0946             endif
0947          else
0948             ii=ijb(1,ib1)
0949             if(jj.lt.nob/2 .and.
0950      $           mmod(nob-jj+1+nob/2).gt. mmod(ii+nob/2)
0951      $           )ii=mmod(nob-jj+1)
0952          endif
0953          weiter=.true.
0954          aa=0.
0955          if(ii.eq.jj) then
0956             if(ish.ge.6)write(ifch,*) 'Rand erreicht'
0957             goto 15
0958          endif
0959          do while(weiter)
0960             if(ii.eq.ijb(1,ib1))then    !linker brk im Feld (ii,jj)
0961                x2=xyb(1,ib1)
0962             else
0963                x2=1.
0964             endif
0965 
0966             if(ii.eq.ijb(1,ib2))then    !rechter brk im Feld (ii,jj)
0967                x1=xyb(1,ib2)
0968             else
0969                x1=0.
0970             endif
0971             if(x1.eq.x2) goto 9999
0972             if(abs(x1-x2)/abs(x1+x2).le.1e-7) goto 9999
0973             f=1.0
0974             if(ipst(ii).ne.ipst(jj))aa=aa+2.*(x2-x1)*(y2-y1)*f*pf(ii,jj)
0975             if(ish.ge.6)write(ifch,*)'ii,jj,x1,x2,aa:',ii,jj
0976      $           ,ipst(ii).ne.ipst(jj),x1,x2,aa
0977      &           ,pf(ii,jj)
0978             if(ii.eq.ijb(1,ib2))then
0979                weiter=.false.
0980             else
0981                ii=mmod(ii+1)
0982                if(ii.eq.jj.or.mmod(ii+jj).eq.0)weiter=.false.
0983             endif
0984          enddo
0985          Amaxn=Amaxn+aa
0986          if(aa.gt.A)goto 10
0987          A=A-aa
0988          if(jj.eq.ijb(2,ib2)) then
0989             if(ish.ge.6)write(ifch,*) 'brk erreicht'
0990             goto 15
0991          endif
0992  12      jj=mmod(jj-1)
0993       enddo
0994       goto 15
0995 
0996  10   continue
0997       yb=A/aa*(y2-y1)+y1
0998       if(ish.ge.6)write(ifch,*)'jj,yb ok:',jj,yb
0999       r=rangen()
1000       ra=aa*r
1001       if(ish.ge.6)write(ifch,*)'r,ra,aa',r,ra,aa
1002                                 !nochmal die letzte ii-Schleife
1003       if(ib1.eq.0)then          !ohne linken brkpt
1004          ii=mmod(jj+1)
1005          if(jj.lt.nob/2)ii=mmod(nob-jj+1)
1006       else
1007         ii=ijb(1,ib1)
1008         if(jj.lt.nob/2 .and.
1009      $       mmod(nob-jj+1+nob/2).gt. mmod(ii+nob/2)
1010      $       )ii=mmod(nob-jj+1)
1011       endif
1012       do while(ra.gt.0.)
1013          if(ii.eq.ijb(1,ib1))then       !linker brk im Feld (ii,jj)
1014             x2=xyb(1,ib1)
1015          else
1016             x2=1.
1017          endif
1018 
1019          if(ii.eq.ijb(1,ib2))then       !rechter brk im Feld (ii,jj)
1020             x1=xyb(1,ib2)
1021          else
1022             x1=0.
1023          endif
1024          f=1.0
1025          ab=0.
1026          if(ipst(ii).ne.ipst(jj)) ab=2.*(x2-x1)*(y2-y1)*f*pf(ii,jj)
1027          if(ish.ge.6)write(ifch,*)'ii,jj,x1,x2,ab,ra:',ii,jj
1028      $        ,ipst(ii).ne.ipst(jj),x1,x2,ab,ra
1029          if(ab.gt.ra)then
1030             xb=ra/ab*(x2-x1)+x1
1031          else
1032             ii=mmod(ii+1)
1033          endif
1034          ra=ra-ab
1035       enddo
1036       if(ish.ge.5)write(ifch,*) 'breakpoint in field ',ii,jj
1037      & ,' at position ',xb,yb
1038       goto 95
1039 
1040 c......................................................................
1041       !von rechts
1042  51   if(ib2.eq.nbr+1)then      !startwert der i-schleife
1043          ii=nob/2-1
1044       else
1045          ii=ijb(1,ib2)
1046       endif
1047       do while (ii.ne.nob/2)
1048          if(ii.eq.ijb(1,ib1))then !linker brk im Streifen (ii)
1049             x2=xyb(1,ib1)
1050          else
1051             x2=1.
1052          endif
1053          if(ii.eq.ijb(1,ib2))then !rechter brk im Streifen (ii)
1054             x1=xyb(1,ib2)
1055          else
1056             x1=0.
1057          endif
1058          if(x1.eq.x2) goto 9999
1059          if(abs(x1-x2)/abs(x1+x2).le.1e-7) goto 9999
1060          if(ish.ge.6)write(ifch,*)'x1,x2 A',x1,x2,A
1061 
1062          if(ib2.eq.nbr+1)then
1063             jj=mmod(ii+1)
1064             if(ii.gt.nob/2)jj=mmod(nob-ii+1)
1065             if(ib1.ge.1.and.jj.gt.ijb(2,ib1))then
1066                if(ish.ge.6)write(ifch,*) 'very special case',ii,jj
1067                goto 13
1068             endif
1069          else
1070             jj=ijb(2,ib2)
1071             if(ii.gt.nob/2 .and. mmod(nob-ii+1).gt.jj)jj=mmod(nob-ii+1)
1072          endif
1073          weiter=.true.
1074          aa=0.
1075          if(ii.eq.jj) then
1076             if(ish.ge.6)write(ifch,*) 'Rand erreicht'
1077             goto 15
1078          endif
1079          do while(weiter)
1080             if(jj.eq.ijb(2,ib1))then    !linker brk im Feld (ii,jj)
1081                y1=xyb(2,ib1)
1082             else
1083                y1=0.
1084             endif
1085             if(jj.eq.ijb(2,ib2))then !rechter brk im Feld (ii,jj)
1086                y2=xyb(2,ib2)
1087             else
1088                y2=1.
1089             endif
1090             if(y1.eq.y2) goto 9999
1091             if(abs(y1-y2)/abs(y1+y2).le.1e-7) goto 9999
1092             f=1.0
1093             if(ipst(ii).ne.ipst(jj))aa=aa+2.*(x2-x1)*(y2-y1)*f*pf(ii,jj)
1094             if(ish.ge.6)write(ifch,*)'ii,jj,x1,x2,aa:',ii,jj
1095      $           ,ipst(ii).ne.ipst(jj),x1,x2,aa
1096      &           ,pf(ii,jj)
1097             if(jj.eq.ijb(2,ib1))then
1098                weiter=.false.
1099             else
1100                jj=mmod(jj+1)
1101                if(jj.eq.ii.or.mmod(ii+jj).eq.0)weiter=.false.
1102             endif
1103          enddo
1104          Amaxn=Amaxn+aa
1105          if(aa.gt.A)goto 14
1106          A=A-aa
1107          if(ii.eq.ijb(1,ib1)) then
1108             if(ish.ge.6)write(ifch,*) 'brk erreicht'
1109             goto 15
1110          endif
1111  13      ii=mmod(ii-1)
1112       enddo
1113       goto 15
1114 
1115 
1116 
1117  14   continue
1118       xb=A/aa*(x2-x1)+x1
1119       if(ish.ge.6)write(ifch,*)'ii,xb ok:',ii,xb
1120       r=rangen()
1121       ra=aa*r
1122       if(ish.ge.6)write(ifch,*)'r,ra,aa',r,ra,aa
1123                                 !nochmal die letzte jj-Schleife
1124       if(ib2.eq.nbr+1)then
1125          jj=mmod(ii+1)
1126          if(ii.gt.nob/2)jj=mmod(nob-ii+1)
1127       else
1128          jj=ijb(2,ib2)
1129          if(ii.gt.nob/2 .and. mmod(nob-ii+1).gt.jj)jj=mmod(nob-ii+1)
1130       endif
1131 
1132       do while(ra.gt.0.)
1133          if(jj.eq.ijb(2,ib1))then       !linker brk im Feld (ii,jj)
1134             y1=xyb(2,ib1)
1135          else
1136             y1=0.
1137          endif
1138 
1139          if(jj.eq.ijb(2,ib2))then       !rechter brk im Feld (ii,jj)
1140             y2=xyb(2,ib2)
1141          else
1142             y2=1.
1143          endif
1144          f=1.0
1145          ab=0.
1146          if(ipst(ii).ne.ipst(jj)) ab=2.*(x2-x1)*(y2-y1)*f*pf(ii,jj)
1147          if(ish.ge.6)write(ifch,*)'ii,jj,x1,x2,ab,ra:',ii,jj
1148      $        ,ipst(ii).ne.ipst(jj),x1,x2,ab,ra
1149          if(ab.gt.ra)then
1150             yb=ra/ab*(y2-y1)+y1
1151          else
1152             jj=mmod(jj+1)
1153          endif
1154          ra=ra-ab
1155       enddo
1156       if(ish.ge.5)write(ifch,*) 'breakpoint in field ',ii,jj
1157      & ,' at position ',xb,yb
1158 
1159  95   continue
1160 
1161 c.....breakpoint accepted......................
1162       nbr=nbr+1
1163       if(nbr.gt.mxnbr-2) stop 'gaksbp: increase parameter mxnbr'
1164       do i=nbr+1,ib1+1,-1
1165          do j=1,2
1166             ijb(j,i)=ijb(j,i-1)
1167             xyb(j,i)=xyb(j,i-1)
1168          enddo
1169          do k=1,4
1170             ptb(k,i)=ptb(k,i-1)
1171          enddo
1172          iflb(i)=iflb(i-1)
1173          ip(i)=ip(i-1)
1174       enddo
1175       ip(ib1+1)=0
1176       ip(ib1+2)=0
1177       ijb(1,ib1+1)=ii
1178       ijb(2,ib1+1)=jj
1179       xyb(1,ib1+1)=xb
1180       xyb(2,ib1+1)=yb
1181 
1182 
1183 cc ..........diquark...............................................
1184 c
1185 c
1186 c      pdiqu=pdiqua
1187 c      if(kzrm.eq.1.or.krm.eq.0)then
1188 cc        pdiqu=exp(log(pdiqua)/sqrt(1.+fkainc*(zzzz-1.)))
1189 cc        pdiqu=exp(log(pdiqua)/(1.+fkainc*(zzzz-1.)))
1190 c      endif
1191 cc      if(iappl.eq.1.and.krm.eq.0)then
1192 cc        pdiqu=pdiquak
1193 cc      endif
1194 c      
1195 c
1196 c      if(nbr.le.2.and.abs(psg(3)/psg(4)).gt.diqcut)pdiqu=0.
1197 c
1198 c      
1199 c      
1200 c      if(rangen().lt.pdiqu.and.iabs(iflb(ib1)).lt.6
1201 c     &  .and.iabs(iflb(ib2)).lt.6)then
1202 c        jqu=2
1203 c      else
1204 c        jqu=1
1205 c      endif
1206 c
1207 
1208 c ..........flavor...............................................
1209 
1210       ig=1
1211       if(iLHC.eq.1)then
1212         if(jorptl(ipst(ii)).eq.-9)ig=ig+1
1213         if(jorptl(ipst(jj)).eq.-9)ig=ig+1
1214       endif
1215 
1216       pu=pui(ig)
1217       pd=pdi(ig)
1218       ps=psi(ig)
1219       pc=pci(ig)
1220       pduu=pduui(ig)
1221       pdud=pdudi(ig)
1222       pdus=pdusi(ig)
1223       pduc=pduci(ig)
1224       pddd=pdddi(ig)
1225       pdds=pddsi(ig)
1226       pddc=pddci(ig)
1227       pdss=pdssi(ig)
1228       pdsc=pdsci(ig)
1229       pdcc=pdcci(ig)
1230 
1231 
1232 !      print *,pu,pd,ps,pc,difud,difus,difuc
1233 c      print *,krm,kdrm,nbr,abs(psg(3)/psg(4))
1234 
1235 c suppress forward (proj) or backward (targ) diquark production 
1236 c in non central strings (baryon spectra in pion collisions or NA49 data)
1237 c ----> this should not be used for meson projectile because of kaon spectra
1238 c       at 100 GeV (Barton, etc) <--------
1239 c      if(kdrm.eq.1.and.ib1.le.1.and.abs(psg(3)/psg(4)).gt.strcut)then
1240 c        ps=0.
1241 c        pc=0.
1242 c      endif
1243 
1244       if(iLHC.eq.1)then
1245 c kill baryon production in case of high pt or forward production
1246         if((ib1.le.1.or.ib2.ge.nbr-1)
1247      &   .and.sqrt(psg(1)**2+psg(2)**2+psg(3)**2)/psg(4).gt.diqcut)then
1248           pduu=0.
1249           pdud=0.
1250           pdus=0.
1251           pduc=0.
1252           pddd=0.
1253           pdds=0.
1254           pddc=0.
1255           pdss=0.
1256           pdsc=0.
1257           pdcc=0.
1258         endif
1259       else
1260 c suppress forward (proj) or backward (targ) diquark production 
1261 c in non central strings (baryon spectra in pion collisions or NA49 data)
1262         if(kdrm.eq.1.and.ib1.le.1.and.abs(psg(3)/psg(4)).gt.diqcut)then
1263           pduu=0.
1264           pdud=0.
1265           pdus=0.
1266           pduc=0.
1267           pddd=0.
1268           pdds=0.
1269           pddc=0.
1270           pdss=0.
1271           pdsc=0.
1272           pdcc=0.
1273         endif
1274       endif
1275 
1276 
1277 
1278       pdiqu=pduu+pdud+pdus+pduc+pddd+pdds+pddc+pdss+pdsc+pdcc
1279       psum=pu+pd+ps+pc+pdiqu
1280 c      print *,krm,kdrm,kky,zzzz,ps/(pu+pd+ps),pdiqu/psum
1281 
1282       r=rangen()*psum
1283       jqu=1
1284       ifl2=0
1285       if(r.gt.psum-pdiqu)then
1286         jqu=2
1287         if(r.gt.pu+pd+ps+pc+pduu+pdud+pdus+pduc+pddd+pdds+pddc+pdss+pdsc
1288      &       .and.pdcc.gt.0.)then
1289           ifl1=4
1290           ifl2=4
1291         elseif(r.gt.pu+pd+ps+pc+pduu+pdud+pdus+pduc+pddd+pdds+pddc+pdss
1292      &         .and.pdsc.gt.0.)then
1293           if(rangen().gt.0.5)then
1294             ifl1=3
1295             ifl2=4
1296           else
1297             ifl1=4
1298             ifl2=3
1299           endif
1300         elseif(r.gt.pu+pd+ps+pc+pduu+pdud+pdus+pduc+pddd+pdds+pddc
1301      &         .and.pdss.gt.0.)then
1302           ifl1=3
1303           ifl2=3
1304         elseif(r.gt.pu+pd+ps+pc+pduu+pdud+pdus+pduc+pddd+pdds
1305      &         .and.pddc.gt.0.)then
1306           if(rangen().gt.0.5)then
1307             ifl1=2
1308             ifl2=4
1309           else
1310             ifl1=4
1311             ifl2=2
1312           endif
1313         elseif(r.gt.pu+pd+ps+pc+pduu+pdud+pdus+pduc+pddd
1314      &         .and.pdds.gt.0.)then
1315           if(rangen().gt.0.5)then
1316             ifl1=2
1317             ifl2=3
1318           else
1319             ifl1=3
1320             ifl2=2
1321           endif
1322         elseif(r.gt.pu+pd+ps+pc+pduu+pdud+pdus+pduc.and.pddd.gt.0.)then
1323           ifl1=2
1324           ifl2=2
1325           jqu=2
1326         elseif(r.gt.pu+pd+ps+pc+pduu+pdud+pdus.and.pduc.gt.0.)then
1327           if(rangen().gt.0.5)then
1328             ifl1=1
1329             ifl2=4
1330           else
1331             ifl1=4
1332             ifl2=1
1333           endif
1334         elseif(r.gt.pu+pd+ps+pc+pduu+pdud.and.pdus.gt.0.)then
1335           if(rangen().gt.0.5)then
1336             ifl1=1
1337             ifl2=3
1338           else
1339             ifl1=3
1340             ifl2=1
1341           endif
1342         elseif(r.gt.pu+pd+ps+pc+pduu.and.pdud.gt.0.)then
1343           if(rangen().gt.0.5)then
1344             ifl1=1
1345             ifl2=2
1346           else
1347             ifl1=2
1348             ifl2=1
1349           endif
1350         else
1351           ifl1=1
1352           ifl2=1
1353         endif
1354       elseif(r.gt.pu+pd+ps.and.pc.gt.0.)then
1355         ifl1=4
1356       elseif(r.gt.pu+pd.and.ps.gt.0.)then
1357         ifl1=3
1358       elseif(r.gt.pu.and.pd.gt.0.)then
1359         ifl1=2
1360       else
1361         ifl1=1
1362       endif
1363 
1364       if(jqu.eq.2)then  !diquark ------
1365         ifl=-min(ifl1,ifl2)*1000-max(ifl1,ifl2)*100
1366       else              !quark ------
1367         ifl=ifl1
1368       endif
1369 
1370       if(iflb(ib1).lt.0.and.iflb(ib1).ge.-6)ifl=-ifl
1371       if(iflb(ib1).gt.1000)ifl=-ifl
1372       iflb(ib1+1)=ifl
1373       if(ish.ge.5)write(ifch,20) ig,ifl,pu,pd,ps,pc,pdiqu/psum
1374      &,sqrt(psg(1)**2+psg(2)**2+psg(3)**2)/psg(4).gt.diqcut
1375      &,pduu,pdud,pdus,pduc,pddd,pdds,pddc,pdss,pdsc,pdcc
1376  20   format('ig,flavor,pu,pd,ps,pc,pdiqu,dcut:',i2,i6,5g13.5,l2/
1377      &,'  diquark u:',4g15.5,/,'  diquark d:',3g15.5,/
1378      &,'  diquark s:',2g15.5,10x,'diquark c:',g15.5)
1379 c..............................pt.......................................
1380 c      if(krm.ne.1.and.kdrm.ne.1)then
1381 c        icub=ib1+1
1382 c        !---------------------------------------------------------------------------
1383 c        !  ib1+1 is the current break point index
1384 c        !             (between 1 and nbr)
1385 c        !  ijb(1,ib1) and ijb(2,ib1) are band indices
1386 c        !         each index from 0 to nob-1 (nob= number of bands)
1387 c        !         0 is left, then come the gluons, then antiquark, then agin gluons
1388 c        !         like q - g - g - g - ~q - g - g - g
1389 c        !--------------------------------------------------------------------------
1390 c        call gaksco(icub-1,icub,icub+1
1391 c     &    ,ijb(1,icub),ijb(2,icub),x1,x2,y1,y2)
1392 c        x=xyb(1,icub)
1393 c        y=xyb(2,icub)
1394 c        aml2=co(1)+x*co(2)+y*co(3)+x*y*co(4)
1395 c        amr2=co(5)+x*co(6)+y*co(7)+x*y*co(8)
1396 c        aml=sign(sqrt(abs(aml2)),aml2)
1397 c        amr=sign(sqrt(abs(amr2)),amr2)
1398 c        !------segment left of current breakpoint icub -----
1399 c        pleft(5)=aml
1400 c        do j=1,4
1401 c          pleft(j)=pst(j,nob+1)-x*pst(j,nob+2)+y*pst(j,nob+3)
1402 c        enddo
1403 c        call utlob4(-1,psg(1),psg(2),psg(3),psg(4),psg(5)
1404 c     $     ,pleft(1),pleft(2),pleft(3),pleft(4))
1405 c        !------segment right of current breakpoint icub-----
1406 c        pright(5)=amr
1407 c        do j=1,4
1408 c          pright(j)=pst(j,nob+4)+x*pst(j,nob+5)-y*pst(j,nob+6)
1409 c        enddo
1410 c        call utlob4(-1,psg(1),psg(2),psg(3),psg(4),psg(5)
1411 c     $     ,pright(1),pright(2),pright(3),pright(4))
1412 c        !-------------------------
1413 c        amt=pleft(5)**2+pleft(1)**2+pleft(2)**2
1414 c        if(amt.gt.0..and.pleft(4)+abs(pleft(3)).gt.0.d0)then
1415 c          amt=sqrt(amt)
1416 c          yleft=sign(1.,pleft(3))*alog((pleft(4)+abs(pleft(3)))/amt)
1417 c        else
1418 c          yleft=0.                  !
1419 c        endif
1420 c        amt=pright(5)**2+pright(1)**2+pright(2)**2
1421 c        if(amt.gt.0..and.pright(4)+abs(pright(3)).gt.0.d0)then
1422 c          amt=sqrt(amt)
1423 c          yright=sign(1.,pright(3))*alog((pright(4)+abs(pright(3)))/amt)
1424 c        else
1425 c          yright=0.                  !
1426 c        endif
1427 c        ybrk=(yleft+yright)/2.
1428 c        yhax=2                  !0.5*yhaha
1429 c        zzipx=1
1430 c        if(ybrk.gt.yhax)then
1431 c          zzipx=zzipx+(zzipp-1)
1432 c        elseif(ybrk.gt.-yhax)then
1433 c          zzipx=zzipx+(zzipp-1)*(ybrk+yhax)/(2*yhax)
1434 c        endif
1435 c        if(ybrk.lt.-yhax)then
1436 c          zzipx=zzipx+(zzipt-1)
1437 c        elseif(ybrk.lt.yhax)then
1438 c          zzipx=zzipx+(zzipt-1)*(yhax-ybrk)/(2*yhax)
1439 c        endif
1440 c      endif
1441       delptfra=0.
1442 c      if(ifl1.eq.3.and.ifl2.eq.0)delptfra=delptfra+ptfrasr
1443 c      if(ifl1.eq.3.or.ifl2.eq.3)delptfra=delptfra+ptfrasr
1444       fnsbp=1
1445       if(nsbp.gt.9)fnsbp=0
1446 c pt kink due to split elastic scattering (if zzzz.ne.1.)
1447       if(iLHC.eq.1)then
1448         pttra=(ptfra+delptfra)*(zzzz-1.)*fnsbp
1449         if(ig.ge.0)then         !quark and soft strings
1450 c        if(ig.eq.1)then         !quark and soft strings
1451 c        if(krm.eq.0)then         !quark and soft strings
1452           pt=ranpt()*(ptfra +  pttra)!/float(jqu)
1453         else
1454           pt=ranpt()*(ptfraqq   +  pttra)!/float(jqu)
1455         endif
1456       else
1457         pttra=(ptfra+delptfra)*(zzzz-1.)*fnsbp
1458         if(jqu.eq.1)then
1459           pt=ranpt()*(ptfra                +  pttra)
1460         else
1461           pt=ranpt()*(ptfraqq              +  pttra)
1462           pt=pt*0.5
1463         endif
1464       endif
1465 
1466       beta=2.*pi*rangen()
1467 
1468       if(ish.ge.5)then
1469         write(ifch,*) 'pt:',pt
1470       endif
1471       ptb(1,ib1+1)=pt*cos(beta)
1472       ptb(2,ib1+1)=pt*sin(beta)
1473       ptb(3,ib1+1)=0.
1474       ptb(4,ib1+1)=0.
1475       if(ish.ge.8)then
1476         write(ifch,*) 'the bands'
1477         write(ifch,'(4g12.6)') (pst(i,ii),i=1,4)
1478         write(ifch,'(4g12.6)') (pst(i,jj),i=1,4)
1479         write(ifch,*) 'pt before rotation'
1480         write(ifch,'(4f12.8)') (ptb(i,ib1+1),i=1,4)
1481       endif
1482       ax=dble(pst(1,ii))+dble(pst(1,jj))
1483       ay=dble(pst(2,ii))+dble(pst(2,jj))
1484       az=dble(pst(3,ii))+dble(pst(3,jj))
1485       ae=dble(pst(4,ii))+dble(pst(4,jj))
1486       am=sqrt(max(1d-10,(ae+az)*(ae-az)-ax**2-ay**2)) !?????????
1487       if(ish.ge.8)then
1488         write(ifch,*) 'boost vector ( region ) '
1489         write(ifch,'(5g12.6)') ax,ay,az,ae,am,pf(ii,jj)
1490       endif
1491       bx=pst(1,ii)
1492       by=pst(2,ii)
1493       bz=pst(3,ii)
1494       be=pst(4,ii)
1495       call utlob2(1,ax,ay,az,ae,am,bx,by,bz,be,60)
1496       if(ish.ge.8) then
1497         write (ifch,*) 'boost of b'
1498         write (ifch,*) 'bx,by,bz,be',bx,by,bz,be
1499       endif
1500       if(abs(bx)+abs(by)+abs(bz).gt.0.)then
1501         call utrot4(-1,bx,by,bz,ptb(1,ib1+1),ptb(2,ib1+1),ptb(3,ib1+1))
1502       else
1503         write(ifmt,*) 'null rot of pt',bx,by,bz
1504         write(ifmt,'(4f12.8)') (ptb(i,ib1+1),i=1,4)
1505       endif
1506       if(ish.ge.8) then
1507         write (ifch,*) 'rot of pt'
1508         write(ifch,'(4f12.8)') (ptb(i,ib1+1),i=1,4)
1509       endif
1510 
1511       call utlob4(-1,ax,ay,az,ae,am
1512      &  ,ptb(1,ib1+1),ptb(2,ib1+1),ptb(3,ib1+1),ptb(4,ib1+1))
1513 
1514       if(ish.ge.8) then
1515         write (ifch,*) 'backboost of pt'
1516         write(ifch,'(4f12.8)') (ptb(i,ib1+1),i=1,4)
1517       endif
1518 c      if(az.eq.0..and.ay.eq.0.)az=1e-8    !not needed if following line commented
1519 c      if(ish.ge.8)write(ifch,*)'rot vector:',ax,ay,az,ae,am
1520 
1521 c.....call utrota(-1,sngl(ax),sngl(ay),sngl(az)  !already commented in nexus2
1522 c....&              ,ptb(1,ib1+1),ptb(2,ib1+1),ptb(3,ib1+1))
1523 
1524       if(ish.ge.6)then
1525         write(ifch,*) 'pt'
1526         write(ifch,'(4g12.6)') (ptb(i,ib1+1),i=1,4)
1527         write (ifch,*) ijb(1,ib1+1),ijb(2,ib1+1),xyb(1,ib1+1)
1528      $       ,xyb(2,ib1+1)
1529       endif
1530 
1531       if(iside.eq.1)then
1532          ib1=ib1+1
1533          ib2=ib2+1
1534       endif
1535 
1536 
1537 c      Amin=0.
1538 c      if(Amax.lt.Amxx) goto 15
1539 c      goto 25
1540 
1541  15   continue
1542       if(ish.ge.6)write(ifch,*) 'Amax:',Amax,Amaxn
1543       if (nbr.eq.nbrold) then
1544          Amax=Amaxn
1545          Amin=0.
1546          if(pb*Amax.le.0..or.ntry.ge.1000)then
1547            goto 9999
1548          endif
1549          goto 26
1550       endif
1551 
1552       if(ish.ge.6)then
1553          write(ifch,*) 0,iflb(0)
1554          do i=1,nbr
1555            if(i.eq.ib2) write(ifch,*) '.................'
1556             write(ifch,'(i3,2x,2(i3),2(" ",g14.7),3x,i5,4(" ",g12.6))')
1557      &           i,ijb(1,i),ijb(2,i),xyb(1,i),xyb(2,i)
1558      &           ,iflb(i),(ptb(j,i),j=1,4)
1559            if(i.eq.ibr1) write(ifch,*) '.................'
1560          enddo
1561          write(ifch,*) nbr+1,iflb(nbr+1)
1562       endif
1563 
1564  9999 if(nbr.eq.nbrold)iok=1
1565       call utprix('gaksbp',ish,ishini,5)
1566       return
1567       end
1568 
1569 c----------------------------------------------------------------------
1570       function ranptcut(xcut)
1571 c----------------------------------------------------------------------
1572 c .........exp(-x**2)
1573 c      if(xcut.gt.0.)then
1574 c        x=sqrt(-(log(rangen()))/(3.1415927/4.)/xcut) !gauss
1575 c      else
1576 c        x=sqrt(-(log(rangen()))/(3.1415927/4.)) !gauss
1577 c      endif
1578  12   x=sqrt(-(log(rangen()))/(3.1415927/4.)) !gauss
1579 
1580 
1581       if(xcut.gt.0.)then
1582         if(rangen().lt.x/xcut)goto 12
1583       endif
1584 
1585       ranptcut=x
1586 
1587       return
1588 
1589 c .........exp(-x)
1590 c  12  xmx=50
1591 c      r=2.
1592 c      do while (r.gt.1.)
1593 c  11    x=sqrt(exp(rangen()*log(1+xmx**2))-1)
1594 c        if(x.eq.0.)goto11
1595 c        r=rangen()  /  ( exp(-x)*(1+x**2) )
1596 c      enddo
1597 c      x=x/2.
1598 
1599       end
1600 
1601 cc----------------------------------------------------------------------
1602 c      function ranpticut(xcut)
1603 cc----------------------------------------------------------------------
1604 c
1605 cc .........exp(-x)
1606 c
1607 c  12  xmx=50
1608 c      r=2.
1609 c      do while (r.gt.1.)
1610 c  11    x=sqrt(exp(rangen()*log(1+xmx**2))-1)
1611 c        if(x.eq.0.)goto11
1612 c        r=rangen()  /  ( exp(-x)*(1+x**2) )
1613 c      enddo
1614 c      x=x/2.
1615 c      if(rangen().gt.xcut/x)goto 12
1616 c
1617 c      ranpticut=x
1618 c
1619 c      end
1620 
1621 c----------------------------------------------------------------------
1622       function ranpt()
1623 c----------------------------------------------------------------------
1624 
1625 c .........exp(-x)
1626       xmx=50
1627       r=2.
1628       do while (r.gt.1.)
1629   11    x=sqrt(exp(rangen()*log(1+xmx**2))-1)
1630         if(x.eq.0.)goto11
1631         r=rangen()  /  ( exp(-x)*(1+x**2) )
1632       enddo
1633       ranpte=x/2.
1634 
1635 c     -------------
1636       ranpt=ranpte
1637       return
1638 c     -------------
1639 cc .........exp(-x**2)
1640 c      ranptg=sqrt(-log(rangen())/(3.1415927/4.)) !gauss
1641 c
1642 cc .........exp(-sqrt(x))
1643 c      xmx=500
1644 c      r=2.
1645 c      do while (r.gt.1.)
1646 c        x=sqrt(exp(rangen()*log(1+xmx**2))-1)
1647 c        r=rangen()  /  ( exp(-sqrt(x))*(1+x**2)/5. )
1648 c      enddo
1649 c      ranpts=x/20.
1650 c
1651 
1652       end
1653 c----------------------------------------------------------------------
1654       function ranptk()
1655 c----------------------------------------------------------------------
1656 
1657 c .........exp(-x)
1658       xmx=50
1659       r=2.
1660       do while (r.gt.1.)
1661   11    x=sqrt(exp(rangen()*log(1+xmx**2))-1)
1662         if(x.eq.0.)goto11
1663         r=rangen()  /  ( exp(-x)*(1+x**2) )
1664       enddo
1665       ranpte=x/2.
1666 
1667 c     -------------
1668       ranptk=ranpte
1669       return
1670 c     -------------
1671 c
1672 cc .........exp(-x**2)
1673 c      ranptg=sqrt(-log(rangen())/(3.1415927/4.)) !gauss
1674 c
1675 cc .........exp(-sqrt(x))
1676 c      xmx=500
1677 c      r=2.
1678 c      do while (r.gt.1.)
1679 c        x=sqrt(exp(rangen()*log(1+xmx**2))-1)
1680 c        r=rangen()  /  ( exp(-sqrt(x))*(1+x**2)/5. )
1681 c      enddo
1682 c      ranpts=x/20.
1683 
1684       end
1685 
1686 c----------------------------------------------------------------------
1687       function ranptd()
1688 c----------------------------------------------------------------------
1689 
1690 c .........exp(-x**2)
1691       ranptg=sqrt(-log(rangen())/(3.1415927/4.)) !gauss
1692 
1693 
1694 c     -------------
1695       ranptd=ranptg
1696       return
1697 c     -------------
1698 c
1699 cc .........exp(-x)
1700 c      xmx=50
1701 c      r=2.
1702 c      do while (r.gt.1.)
1703 c  11    x=sqrt(exp(rangen()*log(1+xmx**2))-1)
1704 c        if(x.eq.0.)goto11
1705 c        r=rangen()  /  ( exp(-x)*(1+x**2) )
1706 c      enddo
1707 c      ranpte=x/2.
1708 c
1709 cc .........exp(-sqrt(x))
1710 c      xmx=500
1711 c      r=2.
1712 c      do while (r.gt.1.)
1713 c        x=sqrt(exp(rangen()*log(1+xmx**2))-1)
1714 c        r=rangen()  /  ( exp(-sqrt(x))*(1+x**2)/5. )
1715 c      enddo
1716 c      ranpts=x/20.
1717 
1718 
1719 
1720 
1721 
1722       end
1723 c----------------------------------------------------------------------
1724       subroutine gakdel(ibr)
1725 c----------------------------------------------------------------------
1726       parameter (mxnbr=500,mxnba=5000)
1727       common /gag/nob,pst(4,0:mxnba),ipst(0:mxnba)
1728      $     ,nbr,ijb(2,0:mxnbr),xyb(2,0:mxnbr)
1729      &     ,ptb(4,0:mxnbr),iflb(0:mxnbr),ip(0:mxnbr),co,x,y
1730       dimension co(12)
1731 
1732       do i=ibr,nbr+1
1733          do j=1,2
1734             ijb(j,i)=ijb(j,i+1)
1735             xyb(j,i)=xyb(j,i+1)
1736          enddo
1737          do k=1,4
1738             ptb(k,i)=ptb(k,i+1)
1739          enddo
1740          iflb(i)=iflb(i+1)
1741          ip(i)=ip(i+1)
1742       enddo
1743       ip(ibr)=0
1744       nbr=nbr-1
1745       end
1746 
1747 c----------------------------------------------------------------------
1748       subroutine gaksco(ibr1,ibr,ibr2,ib,jb,x1,x2,y1,y2)
1749 c----------------------------------------------------------------------
1750       include 'epos.inc'
1751 
1752       parameter (mxnbr=500,mxnba=5000)
1753       common /gag/nob,pst(4,0:mxnba),ipst(0:mxnba)
1754      $     ,nbr,ijb(2,0:mxnbr),xyb(2,0:mxnbr)
1755      &     ,ptb(4,0:mxnbr),iflb(0:mxnbr),ip(0:mxnbr),co,x,y
1756       dimension co(12)
1757       logical weiter
1758 
1759       pf(i,j)=pst(4,i)*pst(4,j)-pst(3,i)*pst(3,j)
1760      &     -pst(2,i)*pst(2,j)-pst(1,i)*pst(1,j)
1761       mmod(i)=mod(mod(i,nob)+nob,nob)
1762 
1763       call utpri('gaksco',ish,ishini,8)
1764 
1765       if(ish.ge.8)then
1766          write(ifch,*) 'zwischen brk:',ibr1,ibr,ibr2,'(',nob,')',nbr
1767          write(ifch,*) 'region:',ib,jb
1768       endif
1769       if(ib.eq.ijb(1,ibr1))then
1770          x2=xyb(1,ibr1)
1771       else
1772          x2=1.
1773       endif
1774       if(ib.eq.ijb(1,ibr2))then
1775          x1=xyb(1,ibr2)
1776       else
1777          x1=0.
1778       endif
1779       if(jb.eq.ijb(2,ibr1))then
1780          y1=xyb(2,ibr1)
1781       else
1782          y1=0.
1783       endif
1784       if(jb.eq.ijb(2,ibr2))then
1785          y2=xyb(2,ibr2)
1786       else
1787          y2=1.
1788       endif
1789 
1790 c.....left side...................................................
1791       n=nob+1
1792       if(ish.ge.8)write(ifch,*)'x1,x2',x1,x2
1793       do i=1,4
1794 cc        pst(i,n)=(x2-x1)*pst(i,ib)+ptb(i,ibr)-ptb(i,ibr1)
1795         pst(i,n)=     x2*pst(i,ib)+ptb(i,ibr)-ptb(i,ibr1)-y1*pst(i,jb)
1796       enddo
1797       if(ish.ge.8)write(ifch,*) 'add a  1 ii',1.,ib
1798       ii=mmod(ib-1)
1799       weiter=.true.
1800       if(ib.eq.ijb(1,ibr1))weiter=.false.
1801       do while(ii.ne.jb.and.weiter) !linker Rand??
1802          f1=1.
1803          if(ii.eq.ijb(1,ibr1))f1=xyb(1,ibr1)
1804          do i=1,4
1805             pst(i,n)=pst(i,n)+f1*pst(i,ii)
1806          enddo
1807          if(ish.ge.8)write(ifch,*) 'add a f1 ii',f1,ii
1808          if(ii.eq.ijb(1,ibr1))weiter=.false.
1809          ii=mmod(ii-1)
1810       enddo
1811       jj=mmod(jb+1)
1812       weiter=.not.weiter
1813       if(jb.eq.ijb(2,ibr1))weiter=.false.
1814       do while(weiter)
1815          f1=1.
1816          if(jj.eq.ijb(2,ibr1))f1=1.-xyb(2,ibr1)
1817          do i=1,4
1818             pst(i,n)=pst(i,n)+f1*pst(i,jj)
1819          enddo
1820          if(ish.ge.8)write(ifch,*) 'add b f1 ii',f1,jj
1821          if(jj.eq.ijb(2,ibr1))weiter=.false.
1822          jj=mmod(jj+1)
1823       enddo
1824       do i=1,4
1825 cc        pst(i,n+1)=(x2-x1)*pst(i,ib)
1826         pst(i,n+1)=        pst(i,ib)
1827 cc        pst(i,n+2)=(y2-y1)*pst(i,jb)
1828         pst(i,n+2)=        pst(i,jb)
1829       enddo
1830       co(1)= pf(n,n)
1831       co(2)=-2.*pf(n,n+1)
1832       co(3)= 2.*pf(n,n+2)
1833       co(4)=-2.*pf(n+1,n+2)
1834       if(ish.ge.8) then
1835         do i=n,n+2
1836           write (ifch,'(4g12.5)') (pst(j,i),j=1,4)
1837         enddo
1838       endif
1839       if(ish.ge.8)write(ifch,*) 'co left:',co(1),co(2),co(3),co(4)
1840 
1841 c.....right side...................................................
1842       n=nob+4
1843       if(ish.ge.8)write(ifch,*)'y1,y2',y1,y2
1844       do i=1,4
1845 cc         pst(i,n)=(y2-y1)*pst(i,jb)-ptb(i,ibr)+ptb(i,ibr2)
1846          pst(i,n)=    y2*pst(i,jb)-ptb(i,ibr)+ptb(i,ibr2)-x1*pst(i,ib)
1847       enddo
1848       if(ish.ge.8)write(ifch,*) 'add a  1 jj',1.,jb
1849       ii=mmod(ib+1)
1850       weiter=.true.
1851       if(ib.eq.ijb(1,ibr2))weiter=.false.
1852       do while(ii.ne.jb.and.weiter)
1853          f1=1.
1854          if(ii.eq.ijb(1,ibr2))f1=1.-xyb(1,ibr2)
1855          do i=1,4
1856             pst(i,n)=pst(i,n)+f1*pst(i,ii)
1857          enddo
1858          if(ish.ge.8)write(ifch,*) 'add a f1 ii',f1,ii
1859          if(ii.eq.ijb(1,ibr2))weiter=.false.
1860          ii=mmod(ii+1)
1861       enddo
1862       jj=mmod(jb-1)
1863       weiter=.not.weiter
1864       if(jb.eq.ijb(2,ibr2))weiter=.false.
1865       do while(weiter)
1866          f1=1.
1867          if(jj.eq.ijb(2,ibr2))f1=xyb(2,ibr2)
1868          do i=1,4
1869             pst(i,n)=pst(i,n)+f1*pst(i,jj)
1870          enddo
1871          if(ish.ge.8)write(ifch,*) 'add b f1 ii',f1,jj
1872          if(jj.eq.ijb(2,ibr2))weiter=.false.
1873          jj=mmod(jj-1)
1874       enddo
1875       do i=1,4
1876 cc         pst(i,n+1)=(x2-x1)*pst(i,ib)
1877          pst(i,n+1)=        pst(i,ib)
1878 cc         pst(i,n+2)=(y2-y1)*pst(i,jb)
1879          pst(i,n+2)=         pst(i,jb)
1880       enddo
1881       co(5)=pf(n,n)
1882       co(6)= 2.*pf(n,n+1)
1883       co(7)=-2.*pf(n,n+2)
1884       co(8)=-2.*pf(n+1,n+2)
1885       if(ish.ge.8) then
1886         do i=n,n+2
1887           write (ifch,'(4g12.5)') (pst(j,i),j=1,4)
1888         enddo
1889       endif
1890       if(ish.ge.8)write(ifch,*) 'co right:',co(5),co(6),co(7),co(8)
1891 
1892 c.....lambda (absolute past).....................................
1893       n=nob+7
1894       do i=1,4
1895 cc         pst(i,n)= x1*pst(i,ib)+y1*pst(i,jb)
1896          pst(i,n)= 0.
1897       enddo
1898       ii=mmod(ib+1)
1899       do while (mmod(ii+jb).ne.0)
1900          if(ish.ge.8)write(ifch,*) 'add lambda',ii
1901          do i=1,4
1902             pst(i,n)=pst(i,n)+pst(i,ii)
1903          enddo
1904          ii=mmod(ii+1)
1905       enddo
1906       do i=1,4
1907 cc         pst(i,n+1)=(x2-x1)*pst(i,ib)
1908 cc         pst(i,n+2)=(y2-y1)*pst(i,jb)
1909          pst(i,n+1)= pst(i,ib)
1910          pst(i,n+2)= pst(i,jb)
1911       enddo
1912       co(9)=     pf(n,n)
1913       co(10)= 2.*pf(n,n+1)
1914       co(11)= 2.*pf(n,n+2)
1915       co(12)= 2.*pf(n+1,n+2)
1916       if(ish.ge.8)write(ifch,*)'co lambda:',co(9),co(10),co(11),co(12)
1917       call utprix('gaksco',ish,ishini,8)
1918       end
1919 
1920 c---------------------------------------------------------------------
1921       subroutine gakanp(ibr1,ibr,ibrr2,aml2,amr2,ala2,iok)
1922 c---------------------------------------------------------------------
1923 c   mass adjustment of fragments
1924 c        ibr1-ibr   ibr-ibrr2
1925 c   where ibr1,ibr,ibrr2 are break point indices
1926 c   aml2,amr2 are the reqired squared masses (if zero -> any mass)
1927 c   iok=0 (ok) or ok=1 (error)
1928 c---------------------------------------------------------------------
1929       include 'epos.inc'
1930       parameter (mxnbr=500,mxnba=5000,mxnin=2000)
1931       common /gag/nob,pst(4,0:mxnba),ipst(0:mxnba)
1932      $     ,nbr,ijb(2,0:mxnbr),xyb(2,0:mxnbr)
1933      &     ,ptb(4,0:mxnbr),iflb(0:mxnbr),ip(0:mxnbr),co,x,y
1934       double precision ax,ay,az,ae,am,bx,by,bz,be,A,B,C
1935       dimension co(12),am2(0:2),nin(2,0:mxnin)
1936       logical weiter
1937 c      pf(i,j)=pst(4,i)*pst(4,j)-pst(3,i)*pst(3,j)
1938 c     &     -pst(2,i)*pst(2,j)-pst(1,i)*pst(1,j)
1939       mmod(i)=mod(mod(i,nob)+nob,nob)
1940       call utpri('gakanp',ish,ishini,6)
1941 
1942       ibr2=ibrr2
1943       if(ish.ge.6)write(ifch,*) ibr1,ibr,ibr2,aml2,amr2,ala2,iok
1944       iok=0
1945       ib=ijb(1,ibr)
1946       jb=ijb(2,ibr)
1947       ni=0
1948  10   do i=1,ni
1949         if((nin(1,i).eq.ib.and.nin(2,i).eq.jb)
1950      $       .or.(ipst(ib).eq.ipst(jb)))then
1951           iok=1
1952           if(ish.ge.4)then
1953             write(ifch,*) 'error ... endless loop'
1954             if(ib.eq.ipst(jb))  write(ifch,*) ' in zero mass region'
1955           endif
1956           goto 9999
1957         endif
1958       enddo
1959       ni=ni+1
1960       if(ni.gt.mxnin)stop'gakanp: increase parameter mxnin  '
1961       nin(1,ni)=ib
1962       nin(2,ni)=jb
1963       if(ish.ge.6)write(ifch,*)
1964       if(ish.ge.6)write(ifch,*) 'ib,jb=',ib,jb
1965       if(ni.ge.2)then
1966          if(ish.ge.6)write(ifch,*)'rotate pt to new band'
1967          if(ish.ge.6)write(ifch,*)'from',nin(1,ni-1),nin(2,ni-1)
1968          if(ish.ge.6)write(ifch,*)'  to',ib,jb
1969          if(ish.ge.6)write(ifch,'(4f12.8)') (ptb(i,ibr),i=1,4)
1970          ax=pst(1,nin(1,ni-1))+pst(1,nin(2,ni-1))
1971          ay=pst(2,nin(1,ni-1))+pst(2,nin(2,ni-1))
1972          az=pst(3,nin(1,ni-1))+pst(3,nin(2,ni-1))
1973          ae=pst(4,nin(1,ni-1))+pst(4,nin(2,ni-1))
1974          am=sngl(sqrt(max(1d-8,dble(ae)**2-dble(ax)**2
1975      $        -dble(ay)**2-dble(az)**2))) !???????????????????????
1976          bx=pst(1,nin(1,ni-1))
1977          by=pst(2,nin(1,ni-1))
1978          bz=pst(3,nin(1,ni-1))
1979          be=pst(4,nin(1,ni-1))
1980          if(ish.ge.6)write (ifch,*) 'bx,by,bz,be',bx,by,bz,be
1981          call utlob2(1,ax,ay,az,ae,am,bx,by,bz,be,60)
1982          if(ish.ge.6)write (ifch,*) 'bx,by,bz,be',bx,by,bz,be
1983          call utlob4( 1,ax,ay,az,ae,am
1984      &        ,ptb(1,ibr),ptb(2,ibr),ptb(3,ibr),ptb(4,ibr))
1985          if(ish.ge.6)write(ifch,'(4f12.8)') (ptb(i,ibr),i=1,4)
1986          if(abs(bx)+abs(by)+abs(bz).gt.0.)then
1987            call utrot4( 1,bx,by,bz,ptb(1,ibr),ptb(2,ibr),ptb(3,ibr))
1988          else
1989            write(ifmt,*) 'null rot of pt (2)',bx,by,bz
1990            write(ifmt,'(4f12.8)') (ptb(i,ibr),i=1,4)
1991          endif
1992          if(ish.ge.6)write(ifch,'(4f12.8)') (ptb(i,ibr),i=1,4)
1993          ax=pst(1,ib)+pst(1,jb)
1994          ay=pst(2,ib)+pst(2,jb)
1995          az=pst(3,ib)+pst(3,jb)
1996          ae=pst(4,ib)+pst(4,jb)
1997          am=sngl(sqrt(max(1d-8,dble(ae)**2-dble(ax)**2
1998      $        -dble(ay)**2-dble(az)**2))) !???????????????????????
1999          if(am.le.1.1e-4)then
2000            if(ish.ge.5)write(ifch,*)'error ... am<1.1e-4'
2001            iok=1
2002            goto 9999
2003          endif
2004          bx=pst(1,ib)
2005          by=pst(2,ib)
2006          bz=pst(3,ib)
2007          be=pst(4,ib)
2008          if(ish.ge.6)write (ifch,*) 'bx,by,bz,be',bx,by,bz,be
2009          if(ish.ge.6)write (ifch,*) 'ax,ay,az,ae',ax,ay,az,ae,am
2010          call utlob2(1,ax,ay,az,ae,am,bx,by,bz,be,60)
2011          if(ish.ge.6)write (ifch,*) 'bx,by,bz,be',bx,by,bz,be
2012          if(abs(bx)+abs(by)+abs(bz).gt.0.)then
2013            call utrot4(-1,bx,by,bz,ptb(1,ibr),ptb(2,ibr),ptb(3,ibr))
2014          else
2015            write(ifmt,*) 'null rot of pt (3)',bx,by,bz
2016            write(ifmt,'(4f12.8)') (ptb(i,ibr),i=1,4)
2017          endif
2018          if(ish.ge.6)write(ifch,'(4f12.8)') (ptb(i,ibr),i=1,4)
2019          call utlob4(-1,ax,ay,az,ae,am
2020      &        ,ptb(1,ibr),ptb(2,ibr),ptb(3,ibr),ptb(4,ibr))
2021          if(ish.ge.6)write(ifch,'(4f12.8)') (ptb(i,ibr),i=1,4)
2022       endif
2023       call gaksco(ibr1,ibr,ibr2,ib,jb,x1,x2,y1,y2)
2024 c      if(ni.eq.1)print *,'------------------------------------',2
2025 cc      x=(xyb(1,ibr)-x1)/(x2-x1)
2026 cc      y=(xyb(2,ibr)-y1)/(y2-y1)
2027       x=xyb(1,ibr)
2028       y=xyb(2,ibr)
2029       am2(0)=aml2
2030       am2(1)=amr2
2031       am2(2)=ala2
2032       if(ish.ge.6)write(ifch,*) ibr1,ibr,ibr2,aml2,amr2,ala2,iok
2033       if(amr2.le.0.)then
2034          l1=2
2035          l2=0
2036       elseif(aml2.le.0.)then
2037          l1=1
2038          l2=2
2039       elseif(ala2.le.0.)then
2040          l1=1
2041          l2=0
2042       else
2043          stop' not like this , please...'
2044       endif
2045       if(ish.ge.6.and.amr2.le.0)write(ifch,*) 'adjust: 1',l1,l2
2046       if(ish.ge.6.and.aml2.le.0)write(ifch,*) 'adjust: 2' ,l1,l2
2047       if(ish.ge.6.and.ala2.le.0)write(ifch,*) 'adjust: 3',l1,l2
2048       i=4*l1
2049       j=4*l2
2050       A = dble(co(i+4))*dble(co(j+3)) - dble(co(j+4))*dble(co(i+3))
2051       B = dble(co(i+4))*dble(co(j+1)) - dble(co(i+3))*dble(co(j+2))
2052      &  + dble(co(i+2))*dble(co(j+3)) - dble(co(i+1))*dble(co(j+4))
2053      &  - dble(am2(l2))*dble(co(i+4)) + dble(am2(l1))*dble(co(j+4))
2054       C = dble(co(i+2))*dble(co(j+1)) - dble(co(i+1))*dble(co(j+2))
2055      &  + dble(am2(l1))*dble(co(j+2)) - dble(am2(l2))*dble(co(i+2))
2056       if (ish.ge.7) then
2057          write(ifch,*) 'ABC,q ',A,B,C,B**2-4.*A*C
2058          if(abs(A).gt.0.d0)then
2059             write(ifch,*) sqrt(max(0d0,B**2-4d0*A*C))/2d0/A-B/A/2d0
2060             write(ifch,*) -sqrt(max(0d0,B**2-4d0*A*C))/2d0/A-B/A/2d0
2061          endif
2062       endif
2063       x=0.
2064       y=0.
2065       xx=0.
2066       yy=0.
2067       if(abs(A).gt.1.d-20.and.B*B-4.*A*C.ge.0.d0)then
2068         y=sngl(sqrt(max(0.d0,B**2-4.*A*C))/2.d0/A-B/A/2.d0)
2069         if(abs(co(i+2)+y*co(i+4)).gt.0.)
2070      &       x=(am2(l1)-co(i+1)-y*co(i+3))/(co(i+2)+y*co(i+4))
2071       elseif(abs(A).le.1.d-20.and.abs(B).gt.0.d0)then
2072         y=-sngl(C/B)
2073         if(abs(co(i+2)+y*co(i+4)).gt.0.)
2074      &       x=(am2(l1)-co(i+1)-y*co(i+3))/(co(i+2)+y*co(i+4))
2075       else
2076         if(ish.ge.5)write(ifch,*)'error ... no solution of quadr equ'
2077         iok=1
2078         goto 9999
2079       endif
2080       if(abs(A).gt.1.d-20.and.B**2-4.*A*C.ge.0.d0)then
2081         yy=sngl(-sqrt(max(0.d0,B**2-4.*A*C))/2.d0/A-B/A/2.d0)
2082         if(abs(co(i+2)+yy*co(i+4)).gt.0.)
2083      &       xx=(am2(l1)-co(i+1)-yy*co(i+3))/(co(i+2)+yy*co(i+4))
2084       elseif(abs(A).le.1.d-20.and.abs(B).gt.0.d0)then
2085         yy=-sngl(C/B)
2086         if(abs(co(i+2)+yy*co(i+4)).gt.0.)
2087      &       xx=(am2(l1)-co(i+1)-yy*co(i+3))/(co(i+2)+yy*co(i+4))
2088       else
2089          if(ish.ge.5)write(ifch,*)'error ... no solution (2) '
2090          iok=1
2091          goto 9999
2092       endif
2093       if(ish.ge.6)then
2094          write(ifch,*) x ,y ,(co(i+2)+ y*co(i+4)),' OK '
2095          write(ifch,*) xx,yy,(co(i+2)+yy*co(i+4)),' OK '
2096       endif
2097       weiter=.true.
2098  50   if(x.gt.x1.and.x.lt.x2.and.y.gt.y1.and.y.lt.y2)then
2099 cc         xyb(1,ibr)=x1+(x2-x1)*x
2100 cc         xyb(2,ibr)=y1+(y2-y1)*y
2101          xyb(1,ibr)=x
2102          xyb(2,ibr)=y
2103          ijb(1,ibr)=ib
2104          ijb(2,ibr)=jb
2105          e1=pst(4,nob+1)-x*pst(4,nob+2)+y*pst(4,nob+3)
2106          e2=pst(4,nob+4)+x*pst(4,nob+5)-y*pst(4,nob+6)
2107          if( e1.lt.0. .or. e2.lt.0. ) then
2108            if(ish.ge.5)write(ifch,*)'error ... e1<0 or e2<0'
2109            iok=1
2110            goto 9999
2111          endif
2112          !amal2=co(1)+co(2)*x+co(3)*y+co(4)*x*y
2113          !amar2=co(5)+co(6)*x+co(7)*y+co(8)*x*y
2114          if(ish.ge.6)then
2115            amal2=co(1)+co(2)*x+co(3)*y+co(4)*x*y
2116            amar2=co(5)+co(6)*x+co(7)*y+co(8)*x*y
2117            write(ifch,*) 'brkshift:',xyb(1,ibr),xyb(2,ibr),ib,jb
2118            write (ifch,*)'E:',e1
2119            write (ifch,*)'E:',e2
2120            write(ifch,'(2(a6,1g12.6))') 'aml:'
2121      &          ,sqrt(max(0.,amal2)),'amr:',sqrt(max(0.,amar2))
2122          endif
2123          i=ibr1+1
2124          do while(i.le.ibr-1)
2125             if((mmod(ijb(1,i)+nob/2) .lt. mmod(ijb(1,ibr)+nob/2)
2126      &           .or.(mmod(ijb(1,i)+nob/2) .eq. mmod(ijb(1,ibr)+nob/2)
2127      &           .and.(xyb(1,i).gt.xyb(1,ibr))))
2128      &           .and.
2129      &           (ijb(2,i) .gt. ijb(2,ibr)
2130      &           .or.(ijb(2,i) .eq. ijb(2,ibr)
2131      &           .and.xyb(2,i).lt.xyb(2,ibr)))) goto 150
2132             if(ish.ge.6) then
2133                write(ifch,*) 'away:'
2134      &          ,i,xyb(1,i),xyb(2,i),ijb(1,i),ijb(2,i)
2135             endif
2136             call gakdel(i)
2137             i=i-1
2138             ibr=ibr-1
2139             ibr2=ibr2-1
2140  150        i=i+1
2141          enddo
2142          i=ibr+1
2143          do while (i.le.ibr2-1)
2144             if((mmod(ijb(1,i)+nob/2) .gt. mmod(ijb(1,ibr)+nob/2)
2145      &           .or.(mmod(ijb(1,i)+nob/2) .eq. mmod(ijb(1,ibr)+nob/2)
2146      &           .and.(xyb(1,i).lt.xyb(1,ibr))))
2147      &           .and.
2148      &           (ijb(2,i) .lt. ijb(2,ibr)
2149      &           .or.(ijb(2,i) .eq. ijb(2,ibr)
2150      &           .and.xyb(2,i).gt.xyb(2,ibr)))) goto 160
2151             if(ish.ge.6) then
2152                write(ifch,*) 'away:'
2153      &          ,i,xyb(1,i),xyb(2,i),ijb(1,i),ijb(2,i)
2154             endif
2155             call gakdel(i)
2156             ibr2=ibr2-1
2157             i=i-1
2158  160        i=i+1
2159          enddo
2160          goto 9999
2161       else
2162         if(x.gt.x2
2163      &    .and.ib.ne.ijb(1,ibr1) !brk-begrenzung
2164      &    .and.mmod(ib-1).ne.jb !linker oder rechter Rand
2165      &    .and.mmod(ib-1+jb).ne.0)then !oben oder unten
2166           ib=mmod(ib-1)
2167           goto 10
2168         endif
2169         if(x.lt.x1
2170      &    .and.ib.ne.ijb(1,ibr2) !brk-begrenzung
2171      &    .and.mmod(ib+1).ne.jb !linker oder rechter Rand
2172      &    .and.mmod(ib+1+jb).ne.0)then !oben oder unten
2173           ib=mmod(ib+1)
2174           goto 10
2175         endif
2176         if(y.gt.y2
2177      &    .and.jb.ne.ijb(2,ibr2) !brk-begrenzung
2178      &    .and.mmod(jb-1).ne.ib !linker oder rechter Rand
2179      &    .and.mmod(jb-1+ib).ne.0)then !oben oder unten
2180           jb=mmod(jb-1)
2181           goto 10
2182         endif
2183         if(y.lt.y1
2184      &    .and.jb.ne.ijb(2,ibr1) !brk-begrenzung
2185      &    .and.mmod(jb+1).ne.ib !linker oder rechter Rand
2186      &    .and.mmod(jb+1+ib).ne.0)then !oben oder unten
2187           jb=mmod(jb+1)
2188           goto 10
2189         endif
2190         if(weiter)then
2191           weiter=.false.
2192           x=xx
2193           y=yy
2194           goto 50
2195         endif
2196         if(ish.ge.5)write(ifch,*)'error ... x,y not in allowed range'
2197         iok=1
2198       endif
2199  9999 if(amr2.eq.0.) ibrr2=ibr2
2200       call utprix('gakanp',ish,ishini,6)
2201       end
2202 
2203 cc----------------------------------------------------------------------
2204 c      subroutine gakstr(ifl)
2205 cc----------------------------------------------------------------------
2206 cc
2207 cc     calculates string-fragments by taking off pt of breakup
2208 cc     do with ifl=1   undo with ifl=-1
2209 cc
2210 cc----------------------------------------------------------------------
2211 c      include 'epos.inc'
2212 c      common /cpptr/ pptr(4,mxptl),ystr(mxptl)
2213 c
2214 c      do i=1,nptl
2215 c        if(istptl(i).eq.29)then
2216 c          nk1=ifrptl(1,i)
2217 c          nk2=ifrptl(2,i)
2218 c          do j=nk1,nk2
2219 c            if ((istptl(j).eq.0.or.istptl(j-1).eq.0).and.j.ne.nk1) then
2220 c              do k=1,4
2221 c                pptl(k,j)=pptl(k,j)+pptr(k,j-1)*ifl
2222 c              enddo
2223 c              !write(ifch,*)"left side back  to ",j,(pptr(k,j-1),k=1,4)
2224 c            endif
2225 c            if ((istptl(j).eq.0.or.istptl(j+1).eq.0).and.j.ne.nk2) then
2226 c              do k=1,4
2227 c                pptl(k,j)=pptl(k,j)-pptr(k,j)*ifl
2228 c              enddo
2229 c              !write(ifch,*)"right side back to ",j,(-pptr(k,j),k=1,4)
2230 c            endif
2231 c            if(ifl.eq.-1.and.istptl(j).eq.0)then
2232 c           e=pptl(1,j)**2+pptl(2,j)**2+pptl(3,j)**2
2233 c     &           +pptl(5,j)**2
2234 c              e=sqrt(e)
2235 c           !dif=abs(e-pptl(4,j))
2236 c           !if(dif.gt.0.01.and.dif/e.gt.0.1)print*,j,e,pptl(4,j)
2237 c              pptl(4,j)=e
2238 c            endif
2239 c          enddo
2240 c        endif
2241 c        if(istptl(i).eq.0)then
2242 c          if ( ifl.eq.1 ) then
2243 c            ystr(i)=sign(1.,pptl(3,i))*alog((pptl(4,i)+abs(pptl(3,i)))
2244 c     *           /sqrt(pptl(5,i)**2+pptl(1,i)**2+pptl(2,i)**2) )
2245 c          endif
2246 c        endif
2247 c      enddo
2248 c      end
2249 
2250 c----------------------------------------------------------------------
2251       subroutine gakli2(nn1,nn2)
2252 c----------------------------------------------------------------------
2253 
2254       include 'epos.inc'
2255       double precision pgampr,rgampr
2256       common/cgampr/pgampr(5),rgampr(4)
2257 c      double precision db1,db2,db3,db4,db5
2258       character label*8,idlabl*8
2259 c      db1=0d0
2260 c      db2=0d0
2261 c      db3=rgampr(4)
2262 c      db4=sqrt(rgampr(1)**2+rgampr(2)**2+rgampr(3)**2)
2263 c      db5=sqrt( db4**2-rgampr(4)**2)
2264       n1=nn1
2265       n2=nn2
2266       if (n1.eq.0) n1=1
2267       if (n2.eq.0) n2=nptl
2268       write (ifch,'(1a4,5a12,4a4,a10,2a4)')
2269      &'no.','px','py','pz','E','m','ior','jor','if1','if2'
2270      &,'id','ist','ity'
2271       do i=n1,n2
2272          if (idptl(i).lt.10000)then
2273             label='        '
2274             label=idlabl(idptl(i))
2275          else
2276             label='        '
2277          endif
2278          chrg=0.
2279          if(iabs(idptl(i)).le.9999)call idchrg(idptl(i),chrg)
2280          write (ifch,125) i,(pptl(j,i),j=1,5),iorptl(i),jorptl(i)
2281      &        ,ifrptl(1,i),ifrptl(2,i),idptl(i)
2282      $        ,chrg             !charge
2283      &        ,istptl(i),ityptl(i),label
2284       enddo
2285  125  format (1i4,5g18.10,4i6,1i10
2286      $     ,1f5.2               !charge
2287      $     ,2i4,'  ',A8
2288 c     $     ,7g12.4,i5
2289      $     )
2290       return
2291       end
2292 
2293 c----------------------------------------------------------------------
2294 c      subroutine gakli4
2295 cc----------------------------------------------------------------------
2296 c
2297 c      include 'epos.inc'
2298 c      parameter (mxnbr=500,mxnba=5000)
2299 c      common /gag/nob,pst(4,0:mxnba),ipst(0:mxnba)
2300 c     $     ,nbr,ijb(2,0:mxnbr),xyb(2,0:mxnbr)
2301 c     &     ,ptb(4,0:mxnbr),iflb(0:mxnbr),ip(0:mxnbr),co,x,y
2302 c      dimension co(12)
2303 c      do i=0,nob-1
2304 c        write (ifch,10) i,(pst(j,i),j=1,4)
2305 c      enddo
2306 c 10   format(1i4,5g18.10)
2307 c      return
2308 c      end
2309 c
2310 cc----------------------------------------------------------------------
2311 c      subroutine gakli3
2312 cc----------------------------------------------------------------------
2313 c
2314 c      include 'epos.inc'
2315 c      parameter (mxnbr=500,mxnba=5000)
2316 c      common /gag/nob,pst(4,0:mxnba),ipst(0:mxnba)
2317 c     $     ,nbr,ijb(2,0:mxnbr),xyb(2,0:mxnbr)
2318 c     &     ,ptb(4,0:mxnbr),iflb(0:mxnbr),ip(0:mxnbr),co,x,y
2319 c      dimension co(12),p1(5),p2(5)
2320 c
2321 c      write(ifch,*) 'particle list of string decay'
2322 c      do i=1,5
2323 c        p1(i)=0.
2324 c      enddo
2325 c      do i=1,nbr+1
2326 c        if(i.lt.nbr+1)then
2327 c          call gaksco(i-1,i,i+1,ijb(1,i),ijb(2,i),x1,x2,y1,y2)
2328 c          if(x2.gt.x1)then
2329 ccc            x=(xyb(1,i)-x1)/(x2-x1)
2330 c            x=xyb(1,i)
2331 c          else
2332 c            x=0.
2333 c          endif
2334 c          if(y2.gt.y1)then
2335 ccc            y=(xyb(2,i)-y1)/(y2-y1)
2336 c            y=xyb(2,i)
2337 c          else
2338 c            y=0.
2339 c          endif
2340 c          aml2=co(1)+x*co(2)+y*co(3)+x*y*co(4)
2341 c          amr2=co(5)+x*co(6)+y*co(7)+x*y*co(8)
2342 c          aml=sign(sqrt(abs(aml2)),aml2)
2343 c          amr=sign(sqrt(abs(amr2)),amr2)
2344 c          do j=1,4
2345 c            p2(j)=pst(j,nob+1)-x*pst(j,nob+2)+y*pst(j,nob+3)
2346 c            p1(j)=p1(j)+p2(j)
2347 c          enddo
2348 c          p2(5)=aml
2349 c        else
2350 c          do j=1,4
2351 c            p2(j)=pst(j,nob+4)+x*pst(j,nob+5)-y*pst(j,nob+6)
2352 c            p1(j)=p1(j)+p2(j)
2353 c          enddo
2354 c          p2(5)=amr
2355 c        endif
2356 c        write(ifch,'(2i4,i6,a,i5,i10,5g14.6)') i-1,i
2357 c     &    ,-iflb(i-1),'==',iflb(i),ip(i)
2358 c     &    ,(p2(j),j=1,5)
2359 c      enddo
2360 c      am2=p1(4)**2-p1(3)**2-p1(2)**2-p1(1)**2
2361 c      p1(5)=sign(sqrt(abs(am2)),am2)
2362 c      write(ifch,'(12x,a60)')
2363 c     &  '------------------------------------------------------------'
2364 c      write(ifch,'(14x,5g14.6)') (p1(j),j=1,5)
2365 c      write(ifch,*)
2366 c
2367 c      end
2368 c
2369 
2370 c---------------------------------------------------------------------
2371       subroutine idress(id,am,idr,iadj)
2372 c---------------------------------------------------------------------
2373       include 'epos.inc'
2374       call idres(id,am,idr,iadj)
2375       if(idr.eq.0)then
2376         return
2377       endif
2378       ids=max(mod(iabs(id)/100,10),mod(iabs(id)/10,10))
2379       if(iabs(idr).le.999) then
2380 c        write (ifch,*) '  ',id,idr,ids
2381         if(ids.le.4)return  !???? if the following is used, bad result
2382         if(ids.le.2)then
2383           idr=sign(iabs(id)+int(rangen()+0.5),id)
2384         elseif(ids.eq.3)then
2385           idr=sign(iabs(id)+int(rangen()+0.6),id)
2386         else
2387           idr=sign(iabs(id)+int(rangen()+0.75),id)
2388         endif
2389 c        write (ifch,*) '->',id,idr
2390         call idmass(idr,am)
2391       elseif(iabs(idr).le.9999)then
2392         if(ids.le.3)return
2393         if(mod(iabs(idr),10).gt.1)then
2394           if(iabs(id).ne.1111.and.iabs(id).ne.2221.and.iabs(id).ne.3331)
2395      $         then
2396             idr=sign(iabs(id)+1,id)
2397             call idmass(idr,am)
2398           else
2399             idr=id
2400             call idmass(idr,am)
2401           endif
2402         endif
2403       endif
2404 
2405       end
2406 
2407 
2408 c---------------------------------------------------------------------
2409       SUBROUTINE gaksphe(sphe,r,mstu41)
2410 
2411 C...Purpose: to perform sphericity tensor analysis to give sphericity,
2412 C...aplanarity and the related event axes. stolen from jetset ;-)
2413       include 'epos.inc'
2414       dimension sphe(4,3)
2415       DIMENSION SM(3,3),SV(3,3)
2416 
2417 C...Calculate matrix to be diagonalized.
2418       NP=0
2419       JA=0
2420       JB=0
2421       JC=0
2422 c      MSTU41=2
2423       DO 110 J1=1,3
2424          DO 100 J2=J1,3
2425             SM(J1,J2)=0.
2426  100     CONTINUE
2427  110  CONTINUE
2428       PS=0.
2429       DO 140 I=1,nptl
2430       IF(istptl(i).ne.0)  GOTO 140
2431       IF(MSTU41.GE.2) THEN
2432          ida=iabs(idptl(i))
2433          IF(ida.EQ.0.OR.ida.EQ.11.OR.ida.EQ.13.OR.ida.EQ.15) GOTO 140
2434          IF(MSTU41.GE.3) then
2435             call idchrg(idptl(i),chrg)
2436             if (abs(chrg).le.0.1) goto 140
2437          endif
2438       ENDIF
2439       NP=NP+1
2440       PA=SQRT(pptl(1,i)**2+pptl(2,I)**2+pptl(3,i)**2)
2441       PWT=1.
2442       IF(ABS(r-2.).GT.0.001) PWT=MAX(1E-10,PA)**(r-2.)
2443       DO 130 J1=1,3
2444          DO 120 J2=J1,3
2445             SM(J1,J2)=SM(J1,J2)+PWT*pptl(j1,i)*pptl(j2,i)
2446  120     CONTINUE
2447  130  CONTINUE
2448       PS=PS+PWT*PA**2
2449  140  CONTINUE
2450 
2451 C...Very low multiplicities (0 or 1) not considered.
2452       IF(NP.LE.1) THEN
2453         if(ish.ge.1)then
2454           CALL utmsg('sphe  ')
2455           write(ifch,*) 'too few particles for analysis'
2456           call utmsgf
2457         endif
2458         sphe(4,1)=-1.
2459         RETURN
2460       ENDIF
2461       DO 160 J1=1,3
2462          DO 150 J2=J1,3
2463             SM(J1,J2)=SM(J1,J2)/PS
2464  150     CONTINUE
2465  160  CONTINUE
2466 
2467 C...Find eigenvalues to matrix (third degree equation).
2468       SQ=(SM(1,1)*SM(2,2)+SM(1,1)*SM(3,3)+SM(2,2)*SM(3,3)-SM(1,2)**2-
2469      &     SM(1,3)**2-SM(2,3)**2)/3.-1./9.
2470       SR=-0.5*(SQ+1./9.+SM(1,1)*SM(2,3)**2+SM(2,2)*SM(1,3)**2+SM(3,3)*
2471      &     SM(1,2)**2-SM(1,1)*SM(2,2)*SM(3,3))
2472      &     +SM(1,2)*SM(1,3)*SM(2,3)+1./27.
2473       SP=COS(ACOS(MAX(MIN(SR/SQRT(-SQ**3),1.),-1.))/3.)
2474       sphe(4,1)=1./3.+SQRT(-SQ)*MAX(2.*SP,SQRT(3.*(1.-SP**2))-SP)
2475       sphe(4,3)=1./3.+SQRT(-SQ)*MIN(2.*SP,-SQRT(3.*(1.-SP**2))-SP)
2476       sphe(4,2)=1.-sphe(4,1)-sphe(4,3)
2477       IF(sphe(4,2).LT.1E-5) THEN
2478         if(ish.ge.1)then
2479           CALL utmsg('gaksphe')
2480           write(ifch,*) 'all particles back-to-back'
2481           call utmsgf
2482         endif
2483         sphe(4,1)=-1.
2484         RETURN
2485       ENDIF
2486 
2487 C...Find first and last eigenvector by solving equation system.
2488       DO 240 I=1,3,2
2489          DO 180 J1=1,3
2490             SV(J1,J1)=SM(J1,J1)-sphe(4,I)
2491             DO 170 J2=J1+1,3
2492                SV(J1,J2)=SM(J1,J2)
2493                SV(J2,J1)=SM(J1,J2)
2494  170        CONTINUE
2495  180     CONTINUE
2496          SMAX=0.
2497          DO 200 J1=1,3
2498             DO 190 J2=1,3
2499                IF(ABS(SV(J1,J2)).LE.SMAX) GOTO 190
2500                JA=J1
2501                JB=J2
2502                SMAX=ABS(SV(J1,J2))
2503  190        CONTINUE
2504  200     CONTINUE
2505          SMAX=0.
2506          DO 220 J3=JA+1,JA+2
2507             J1=J3-3*((J3-1)/3)
2508             RL=SV(J1,JB)/SV(JA,JB)
2509             DO 210 J2=1,3
2510                SV(J1,J2)=SV(J1,J2)-RL*SV(JA,J2)
2511                IF(ABS(SV(J1,J2)).LE.SMAX) GOTO 210
2512                JC=J1
2513                SMAX=ABS(SV(J1,J2))
2514  210        CONTINUE
2515  220     CONTINUE
2516          JB1=JB+1-3*(JB/3)
2517          JB2=JB+2-3*((JB+1)/3)
2518          sphe(JB1,I)=-SV(JC,JB2)
2519          sphe(jb2,I)=SV(JC,JB1)
2520          sphe(jb,I)=-(SV(JA,JB1)*sphe(jb1,I)+SV(JA,JB2)
2521      &        *sphe(jb2,I))/SV(JA,JB)
2522          PA=SQRT(sphe(1,I)**2+sphe(2,I)**2+sphe(3,I)**2)
2523          SGN=(-1.)**INT(rangen()+0.5)
2524          DO 230 J=1,3
2525             sphe(j,I)=SGN*sphe(j,I)/PA
2526  230     CONTINUE
2527  240  CONTINUE
2528 
2529 C...Middle axis orthogonal to other two. Fill other codes.
2530       SGN=(-1.)**INT(rangen()+0.5)
2531       sphe(1,2)=SGN*(sphe(2,1)*sphe(3,3)-sphe(3,1)*sphe(2,3))
2532       sphe(2,2)=SGN*(sphe(3,1)*sphe(1,3)-sphe(1,1)*sphe(3,3))
2533       sphe(3,2)=SGN*(sphe(1,1)*sphe(2,3)-sphe(2,1)*sphe(1,3))
2534 
2535       do i=1,3
2536          do j=1,4
2537             pptl(j,nptl+i)=sphe(j,i)
2538          enddo
2539       enddo
2540 
2541 
2542 C...Calculate sphericity and aplanarity. Select storing option.
2543 ccc      SPH=1.5*(sphe(4,2)+sphe(4,3))
2544 ccc      APL=1.5*sphe(4,3)
2545 
2546       RETURN
2547       END
2548 
2549 C*********************************************************************
2550 
2551       SUBROUTINE gakthru(thru,mstu41)
2552 
2553 C...Purpose: to perform thrust analysis to give thrust, oblateness
2554 C...and the related event axes. stolen from jetset ;-)
2555       include 'epos.inc'
2556       DIMENSION TDI(3),TPR(3)
2557       dimension thru(4,3)
2558 
2559 C...Take copy of particles that are to be considered in thrust analysis.
2560       IAGR=0
2561       NP=0
2562       PS=0.
2563 c      MSTU41=2
2564       MSTU44=4
2565       MSTU45=2
2566       PARU42=1.0
2567       PARU48=0.0000001
2568       DO 100 I=1,nptl
2569          IF(istptl(i).ne.0)goto 100
2570          ida=iabs(idptl(i))
2571          IF(ida.EQ.0.OR.ida.EQ.11.OR.ida.EQ.13.OR.ida.EQ.15)GOTO 100
2572          IF(MSTU41.GE.3) then
2573             call idchrg(idptl(i),chrg)
2574             if (abs(chrg).le.0.1) goto 100
2575          endif
2576 
2577          IF(nptl+NP.GE.mxptl) THEN
2578             CALL utstop('gakthru: no more memory left in cptl&')
2579             thru(4,1)=-1.
2580             RETURN
2581          ENDIF
2582          NP=NP+1
2583 c         K(nptl+NP,1)=23
2584          pptl(1,nptl+NP)=pptl(1,I)
2585          pptl(2,nptl+NP)=pptl(2,I)
2586          pptl(3,nptl+NP)=pptl(3,I)
2587          pptl(4,nptl+NP)=SQRT(pptl(1,I)**2+pptl(2,I)**2+pptl(3,I)**2)
2588          pptl(5,nptl+NP)=1.
2589          IF(ABS(PARU42-1.).GT.0.001)
2590      &        pptl(5,nptl+NP)=pptl(4,nptl+NP)**(PARU42-1.)
2591          PS=PS+pptl(4,nptl+NP)*pptl(5,nptl+NP)
2592  100  CONTINUE
2593 
2594 C...Very low multiplicities (0 or 1) not considered.
2595       IF(NP.LE.1) THEN
2596          CALL utmsg('thru  ')
2597          write(ifch,*) 'too few particles for analysis'
2598          call utmsgf
2599          thru(4,1)=-1
2600          RETURN
2601       ENDIF
2602 
2603 C...Loop over thrust and major. T axis along z direction in latter case.
2604       DO 320 ILD=1,2
2605       IF(ILD.EQ.2) THEN
2606 c        PHI=GAKANG(pptl(1,nptl+NP+1),pptl(2,nptl+NP+1))
2607 c        CALL lurot(nptl+1,nptl+NP+1,0.,-PHI)
2608 c        THE=GAKANG(pptl(3,nptl+NP+1),pptl(1,nptl+NP+1))
2609 c        CALL lurot(nptl+1,nptl+NP+1,-THE,0.)
2610          ax=pptl(1,nptl+NP+1)
2611          ay=pptl(2,nptl+NP+1)
2612          az=pptl(3,nptl+NP+1)
2613          do irot=nptl+1,nptl+NP+1
2614             call utrota(1,ax,ay,az
2615      &           ,pptl(1,irot),pptl(2,irot),pptl(3,irot))
2616          enddo
2617          if(np.eq.2)then
2618            pptl(1,nptl+NP+2)=1.
2619            pptl(2,nptl+NP+2)=0.
2620            pptl(3,nptl+NP+2)=0.
2621            pptl(4,nptl+NP+2)=0.
2622            goto 325
2623          endif
2624       ENDIF
2625 
2626 C...Find and order particles with highest p (pT for major).
2627       DO 110 ILF=nptl+NP+4,nptl+NP+MSTU44+4
2628          pptl(4,ILF)=0.
2629  110  CONTINUE
2630       DO 160 I=nptl+1,nptl+NP
2631          IF(ILD.EQ.2) pptl(4,I)=SQRT(pptl(1,I)**2+pptl(2,I)**2)
2632          DO 130 ILF=nptl+NP+MSTU44+3,nptl+NP+4,-1
2633             IF(pptl(4,I).LE.pptl(4,ILF)) GOTO 140
2634             DO 120 J=1,5
2635                pptl(j,ILF+1)=pptl(j,ILF)
2636  120        CONTINUE
2637  130     CONTINUE
2638          ILF=nptl+NP+3
2639  140     DO 150 J=1,5
2640             pptl(j,ILF+1)=pptl(j,I)
2641  150     CONTINUE
2642  160  CONTINUE
2643 
2644 C...Find and order initial axes with highest thrust (major).
2645       DO 170 ILG=nptl+NP+MSTU44+5,nptl+NP+MSTU44+15
2646          pptl(4,ILG)=0.
2647  170  CONTINUE
2648       NC=2**(MIN(MSTU44,NP)-1)
2649       DO 250 ILC=1,NC
2650          DO 180 J=1,3
2651             TDI(J)=0.
2652  180     CONTINUE
2653          DO 200 ILF=1,MIN(MSTU44,NP)
2654             SGN=pptl(5,nptl+NP+ILF+3)
2655             IF(2**ILF*((ILC+2**(ILF-1)-1)/2**ILF).GE.ILC) SGN=-SGN
2656             DO 190 J=1,4-ILD
2657                TDI(J)=TDI(J)+SGN*pptl(j,nptl+NP+ILF+3)
2658  190        CONTINUE
2659  200     CONTINUE
2660          TDS=TDI(1)**2+TDI(2)**2+TDI(3)**2
2661          DO 220 ILG=nptl+NP+MSTU44+MIN(ILC,10)+4,nptl+NP+MSTU44+5,-1
2662             IF(TDS.LE.pptl(4,ILG)) GOTO 230
2663             DO 210 J=1,4
2664                pptl(j,ILG+1)=pptl(j,ILG)
2665  210        CONTINUE
2666  220     CONTINUE
2667          ILG=nptl+NP+MSTU44+4
2668  230     DO 240 J=1,3
2669             pptl(j,ILG+1)=TDI(J)
2670  240     CONTINUE
2671          pptl(4,ILG+1)=TDS
2672  250  CONTINUE
2673 
2674 C...  Iterate direction of axis until stable maximum.
2675       pptl(4,nptl+NP+ILD)=0.
2676       ILG=0
2677  260  ILG=ILG+1
2678       THP=0.
2679  270  THPS=THP
2680       DO 280 J=1,3
2681          IF(THP.LE.1E-10) TDI(J)=pptl(j,nptl+NP+MSTU44+4+ILG)
2682          IF(THP.GT.1E-10) TDI(J)=TPR(J)
2683          TPR(J)=0.
2684  280  CONTINUE
2685       DO 300 I=nptl+1,nptl+NP
2686          SGN=SIGN(pptl(5,I),TDI(1)*pptl(1,I)
2687      &        +TDI(2)*pptl(2,I)+TDI(3)*pptl(3,I))
2688          DO 290 J=1,4-ILD
2689             TPR(J)=TPR(J)+SGN*pptl(j,I)
2690  290     CONTINUE
2691  300  CONTINUE
2692       THP=SQRT(TPR(1)**2+TPR(2)**2+TPR(3)**2)/PS
2693       IF(THP.GE.THPS+PARU48) GOTO 270
2694 
2695 C...  Save good axis. Try new initial axis until a number of tries agree.
2696       IF(THP.LT.pptl(4,nptl+NP+ILD)-PARU48.AND.ILG.LT.MIN(10,NC))
2697      &     GOTO 260
2698       IF(THP.GT.pptl(4,nptl+NP+ILD)+PARU48)
2699      $     THEN
2700          IAGR=0
2701          SGN=(-1.)**INT(rangen()+0.5)
2702          DO 310 J=1,3
2703             pptl(j,nptl+NP+ILD)=SGN*TPR(J)/(PS*THP)
2704  310     CONTINUE
2705          pptl(4,nptl+NP+ILD)=THP
2706          pptl(5,nptl+NP+ILD)=0.
2707       ENDIF
2708       IAGR=IAGR+1
2709       IF(IAGR.LT.MSTU45.AND.ILG.LT.MIN(10,NC)) GOTO 260
2710  320  CONTINUE
2711 
2712 C...  Find minor axis and value by orthogonality.
2713  325  SGN=(-1.)**INT(rangen()+0.5)
2714       pptl(1,nptl+NP+3)=-SGN*pptl(2,nptl+NP+2)
2715       pptl(2,nptl+NP+3)=SGN*pptl(1,nptl+NP+2)
2716       pptl(3,nptl+NP+3)=0.
2717       THP=0.
2718       DO 330 I=nptl+1,nptl+NP
2719          THP=THP+pptl(5,I)*ABS(pptl(1,nptl+NP+3)*pptl(1,I)
2720      &        +pptl(2,nptl+NP+3)*pptl(2,I))
2721  330  CONTINUE
2722       pptl(4,nptl+NP+3)=THP/PS
2723       pptl(5,nptl+NP+3)=0.
2724 
2725 
2726 C...  Fill axis information. Rotate back to original coordinate system.
2727       do irot=nptl+NP+1,nptl+NP+3
2728          call utrota(-1,ax,ay,az
2729      &        ,pptl(1,irot),pptl(2,irot),pptl(3,irot))
2730       enddo
2731 
2732       do ild=1,3
2733          do j=1,4
2734             thru(j,ild)=pptl(j,nptl+NP+ild)
2735          enddo
2736       enddo
2737 
2738 C...Calculate thrust and oblateness. Select storing option.
2739 ccc      THR=thru(4,1)
2740 ccc      OBL=thru(4,2)-thru(4,3)
2741 
2742       RETURN
2743       END
2744 
2745       subroutine gakjet(ijadu)
2746       include 'epos.inc'
2747       common/njjjj/njets(5,2,mxbins)
2748 c      nmin=xpar1
2749 c      nmax=xpar2
2750       if(nrevt.eq.1)then
2751         do i=1,5
2752           do j=1,mxbins
2753             njets(i,ijadu,j)=0
2754           enddo
2755         enddo
2756       endif
2757       do i=1,nrbins
2758         ycut=xminim*(xmaxim/xminim)**((float(i)-0.5)/nrbins)
2759 c        if(iologe.ne.1)ycut=xminim+(xmaxim-xminim)/nrbins*(nrbins)
2760         nj=gaknjt(ycut,ijadu)
2761         if(nj.ge.1.and.nj.le.5)then
2762           njets(nj,ijadu,i)=njets(nj,ijadu,i)+1
2763         endif
2764       enddo
2765       end
2766 
2767       subroutine gakjto
2768       include 'epos.inc'
2769       common/njjjj/njets(5,2,mxbins)
2770       n=xpar4
2771       ijadu=xpar3
2772       do i=1,nrbins
2773         ycut=xminim*(xmaxim/xminim)**((float(i)-0.5)/nrbins)
2774         write (ifhi,*) ycut,float(njets(n,ijadu,i))/nrevt
2775      $       ,sqrt(1.*njets(n,ijadu,i))/nrevt
2776       enddo
2777       end
2778 C*********************************************************************
2779       function gaknjt(ycut,ijadu)
2780 c
2781 c     ijadu   1 =  JADE    2=DURHAM
2782 c     ycut - max. distance
2783 c
2784       include 'epos.inc'
2785 
2786       a2j(i,j)=2.*pptl(4,i)*pptl(4,j)*(1.-(pptl(1,i)*pptl(1,j)
2787      &     +pptl(2,i)*pptl(2,j)+pptl(3,i)*pptl(3,j))
2788      &     /(sqrt(pptl(1,i)**2+pptl(2,i)**2+pptl(3,i)**2)
2789      &       *sqrt(pptl(1,j)**2+pptl(2,j)**2+pptl(3,j)**2)))/evis**2
2790 
2791       a2d(i,j)=2.*min(pptl(4,i)**2,pptl(4,j)**2)
2792      &     *(1.-(pptl(1,i)*pptl(1,j)
2793      &     +pptl(2,i)*pptl(2,j)+pptl(3,i)*pptl(3,j))
2794      &     /(sqrt(pptl(1,i)**2+pptl(2,i)**2+pptl(3,i)**2)
2795      &       *sqrt(pptl(1,j)**2+pptl(2,j)**2+pptl(3,j)**2)))/evis**2
2796 
2797       bet(i)=sqrt(pptl(1,i)**2+pptl(2,i)**2+pptl(3,i)**2)
2798 
2799       ska(i,j)=pptl(1,i)*pptl(1,j)+pptl(2,i)*pptl(2,j)
2800      &     +pptl(3,i)*pptl(3,j)
2801 
2802       a2c(i,j)= ((bet(i)*bet(j)-ska(i,j))
2803      &     *2.*bet(i)*bet(j)  /  (0.00001+bet(i)+bet(j))**2 )
2804 
2805       evis=0.
2806       nn=0
2807       do i=1,nptl
2808          if (istptl(i).eq.0) then
2809             nn=nn+1
2810             do j=1,5
2811                pptl(j,nptl+nn)=pptl(j,i)
2812             enddo
2813             evis=evis+pptl(4,i)
2814             jorptl(i)=nn
2815          endif
2816       enddo
2817       iflag=0
2818       i1=0
2819       j1=0
2820       do while (iflag.eq.0.and.nn.ge.2)
2821          a2min=ycut
2822          iflag=1
2823          do i=nptl+1,nptl+nn-1
2824             do j=i+1,nptl+nn
2825                if(ijadu.eq.1)then
2826                   a2=a2j(i,j)
2827                elseif(ijadu.eq.2) then
2828                   a2=a2d(i,j)
2829                else
2830                   a2=a2c(i,j)
2831                endif
2832                if (a2.lt.a2min) then
2833                   iflag=0
2834                   i1=i
2835                   j1=j
2836                   a2min=a2
2837                endif
2838             enddo
2839          enddo
2840          if (iflag.eq.0) then
2841             do i=1,4
2842                pptl(i,i1)=pptl(i,i1)+pptl(i,j1)
2843             enddo
2844             do i=1,nptl
2845                if(istptl(i).eq.0.and.jorptl(i).eq.j1-nptl)
2846      &              jorptl(i)=i1-nptl
2847                if(istptl(i).eq.0.and.jorptl(i)+nptl.gt.j1)
2848      &              jorptl(i)=jorptl(i)-1
2849             enddo
2850             do i=j1,nptl+nn
2851                do j=1,5
2852                   pptl(j,i)=pptl(j,i+1)
2853                enddo
2854             enddo
2855             nn=nn-1
2856          endif
2857       enddo
2858       do i=nptl+1,nptl+nn
2859          istptl(i)=-1
2860          jorptl(i)=i-nptl
2861          pptl(5,i)=sqrt(max(0.,(pptl(4,i)+pptl(3,i))
2862      &   *(pptl(4,i)-pptl(3,i))-pptl(2,i)**2-pptl(1,i)**2))
2863       enddo
2864       do i=nptl+1,nptl+nn-1
2865          do j=i+1,nptl+nn
2866             if(pptl(4,jorptl(i)+nptl).lt.pptl(4,jorptl(j)+nptl))then
2867                k=jorptl(i)
2868                jorptl(i)=jorptl(j)
2869                jorptl(j)=k
2870             endif
2871          enddo
2872       enddo
2873       do i=nptl+1,nptl+nn
2874          idptl(nptl+jorptl(i))=9910+i-nptl
2875       enddo
2876       do i=1,nptl
2877         jorptl(i)=0             !jorptl(nptl+jorptl(i))
2878       enddo
2879 c      nptl=nptl+nn
2880 
2881       gaknjt=nn
2882       return
2883       end
2884 
2885       subroutine idtr5(id,ic)
2886       integer ic(2)
2887       ic(1)=0
2888       ic(2)=0
2889       ii=1
2890       if(id.lt.0)ii=2
2891       i1=1
2892       if(iabs(id).gt.999)i1=3
2893       do i=i1,int(log(abs(float(id)))/log(10.))+1
2894         j=mod(iabs(id)/10**(i-1),10)
2895         if(j.gt.0)then
2896           ic(ii)=ic(ii)+10**(6-j)
2897         endif
2898       enddo
2899       return
2900       end
2901 
2902       function ammin(id1,id2)
2903       dimension ic1(2),ic2(2),jc2(6,2),jc1(6,2)
2904       call idtr5(id1,ic1)
2905       call idtr5(id2,ic2)
2906       call idcomk(ic1)
2907       call idcomk(ic2)
2908       call iddeco(ic1,jc1)
2909       call iddeco(ic2,jc2)
2910       ammin=utamnx(jc1,jc2)
2911       end
2912 
2913 
2914 c      function idtr(id1,id2)
2915 c      integer ic(2),id(2)
2916 c      id(1)=id1
2917 c      id(2)=id2
2918 c      do i=1,2
2919 c        ic(i)=0
2920 c      enddo
2921 c      do j=1,2
2922 c        ii=1
2923 c        if(id(j).lt.0)ii=2
2924 c        i1=1
2925 c        if(iabs(id(j)).gt.999)i1=3
2926 c        do i=i1,int(log(abs(float(id(j))))/log(10.))+1
2927 c          jj=mod(iabs(id(j))/10**(i-1),10)
2928 c          if(jj.gt.0)then
2929 c            ic(ii)=ic(ii)+10**(6-jj)
2930 c          endif
2931 c        enddo
2932 c      enddo
2933 c      idtr=idtra(ic,0,0,3)
2934 c      if(idtr.ne.idsp(id1,id2))then
2935 c        write (*,'(4i6)') idtr,idsp(id1,id2),id1,id2
2936 c      endif
2937 c      return
2938 c      end
2939 
2940       function idsp(id1,id2)
2941       ia1=iabs(id1)
2942       ia2=iabs(id2)
2943       if(ia1.ge.1000.and.ia2.ge.1000)then
2944         idsp=0
2945         isign=0
2946       elseif(ia1.le.1000.and.ia2.le.1000)then
2947         idsp=min(ia1,ia2)*100+max(ia1,ia2)*10
2948         isign=1
2949         if(max(ia1,ia2).ne.-min(id1,id2)) isign = -1
2950         if(idsp.eq.220)idsp=110
2951         if(idsp.eq.330)idsp=220
2952       else
2953         isign=1
2954         if(id1.lt.0.and.id2.lt.0)isign=-1
2955         idb=min(ia1,ia2)
2956         if(idb.eq.5)then
2957           idsp=0
2958           return
2959         endif
2960         ida=max(ia1,ia2)
2961         ida1=ida/1000
2962         ida2=mod(ida/100,10)
2963         if(idb.le.ida1)then
2964           idsp=idb*1000+ida/10
2965         elseif(idb.le.ida2)then
2966           idsp=ida1*1000+idb*100+ida2*10
2967         else
2968           idsp=ida+idb*10
2969         endif
2970         if(ida1.eq.ida2.and.ida2.eq.idb)idsp=idsp+1
2971       endif
2972       idsp=idsp*isign
2973       return
2974       end
2975