EPOS/epos167/epos-fra-163.f

   1 c-----------------------------------------------------------------------
   2       subroutine gakfra(iret)
   3 c-----------------------------------------------------------------------
   4
   5       include 'epos.inc'
   6       include 'epos.incems'
   7       include 'epos.incsem'
   8       parameter (eta=0.4,etap=0.1)
   9       parameter (mxnbr=500,mxnba=5000)
  10       common /gag/nob,pst(4,0:mxnba),ipst(0:mxnba)
  11      $     ,nbr,ijb(2,0:mxnbr),xyb(2,0:mxnbr)
  12      &     ,ptb(4,0:mxnbr),iflb(0:mxnbr),ip(0:mxnbr),co,x,y
  13 c      common /cpptr/ pptr(4,mxptl),ystr(mxptl)
  14       double precision p1,p12,p2
  15       dimension co(12),ind(0:mxnbr),p1(5),p2(5)
  16       dimension ic2(2),ic1(2),ic(2)
  17       common/pb/pb /cnsbp/nsbp  /cn8ptl/n8ptl
  18       common/czz/kky,krm,zzod,zzop,kdrm,zzoq,zzos,zzrm,zzipp,zzipt
  19      &,zzid,zzip,zzzzr,itp
  20       double precision psg
  21       common/zpsg/psg(5)
  22       logical leadstr,go
  23
  24       pf(i,j)=(pst(4,i)-pst(3,i))*(pst(4,j)+pst(3,j))
  25      &     -pst(2,i)*pst(2,j)-pst(1,i)*pst(1,j)
  26
  27       nob=0
  28       call utpri('gakfra',ish,ishini,4)
  29       iret=0
  30 c.....search string in pptl-list...................................
  31       nkmax=nptl
  32       nk1=maproj+matarg+1
  33  2    continue
  34       if(nclean.gt.0.and.nptl.gt.mxptl/2)then
  35         nnnpt=0
  36         do iii=maproj+matarg+1,nptl
  37           go=.true.
  38           if(nclean.eq.1.and.istptl(iii).le.istmax)go=.false.
  39           if(go.and.mod(istptl(iii),10).ne.0)then
  40             istptl(iii)=99
  41             nnnpt=nnnpt+1
  42           endif
  43         enddo
  44         if(nnnpt.gt.mxptl-nptl)then
  45           nptl0=nptl
  46           call utclea(maproj+matarg+1,nptl0)
  47           nkmax=nptl
  48           nk1=maproj+matarg+1
  49 c          write(ifch,'(a)')' ---------after utclea------------'
  50 c          call blist('&',nk1,nkmax)
  51         endif
  52       endif
  53
  54       do while(nk1.le.nkmax)
  55         if(istptl(nk1).eq.20)then
  56           nk2=nk1+1
  57           do while(idptl(nk2).eq.9)
  58             nk2=nk2+1
  59           enddo
  60           goto 3                !ok, string from nk1 to nk2
  61         else
  62           nk1=nk1+1
  63         endif
  64       enddo
  65
  66  431  continue
  67
  68
  69       goto 9999                 !no more string
  70
  71
  72 c.......................................................................
  73 c                       decay string  nk1 - nk2
  74 c.......................................................................
  75
  76  3    nob=-1
  77       if(ish.ge.3)then
  78         write(ifch,'(/1x,25a1/5x,a/1x,25a1/)')
  79      *  ('-',k=1,25),'string decay',('-',k=1,25)
  80         write(ifch,'(a)')' ---------string------------'
  81         call blist('&',nk1,nk2)
  82       endif
  83
  84       kky=0
  85       if(nk2-nk1.gt.1.or.iappl.eq.6)kky=1
  86       itp=ityptl(nk1)
  87       krm=0
  88       if(itp.ge.50.and.itp.le.59.or.itp.ge.40.and.itp.le.49)krm=1
  89       kdrm=0
  90       if(itp.eq.43.or.itp.eq.53)kdrm=1
  91 c special behaviour of remnant of reactions with lost Pomerons
  92       if(itp.eq.47.or.itp.eq.57)then
  93         kdrm=1
  94         krm=0
  95       endif
  96       zzzz=1.
  97       zzzzp=1.
  98       zzzzt=1.
  99       zzod=1.
 100       zzos=1.
 101       zzop=1.
 102       zzoq=1.
 103       zzrm=1.
 104       zzipp=1.
 105       zzipt=1.
 106       zzid=1.
 107       zzip=1.
 108       if(isplit.eq.1)then
 109         zzzz=zpaptl(nk1)+zpaptl(nk2)
 110         zzzzp=zpaptl(nk1)
 111         zzzzt=zpaptl(nk2)
 112         zzid=min( 1+zidmax , 1+zidinc  *zzzz ) !<-----
 113         zzip=min( 1+zopmax , 1+zopinc  *zzzz ) !<-----
 114 c        zzrm=min( zrmmax   , zrminc    *zzzz ) !<-----
 115         if(itp.ge.40.and.itp.le.49)then
 116           zzzzr=zzzzp
 117         elseif(itp.ge.50.and.itp.le.59)then
 118           zzzzr=zzzzt
 119         endif
 120         zzod=min( 1+zodmax , 1+zodinc  *zzzzr ) !<-----
 121         zzos=min( 1+zosmax , 1+zosinc  *zzzzr ) !<-----
 122         zzop=min( 1+zopmax , 1+zopinc  *zzzzr ) !<-----
 123         zzoq=min( 1+2.*zopmax , 1+ 2.*zopinc  *zzzzr ) !<-----
 124         zzipp=min( 1+zipmax , 1+zipinc  *zzzzp ) !<-----
 125         zzipt=min( 1+zipmax , 1+zipinc  *zzzzt ) !<-----
 126       endif
 127       do k=1,5
 128         p1(k)=0d0
 129       enddo
 130       do i=nk1,nk2
 131         do k=1,5
 132           p2(k)=dble(pptl(k,i))
 133         enddo
 134         p2(4)=sqrt(p2(3)**2+p2(2)**2+p2(1)**2+p2(5)**2)
 135         do k=1,4
 136           p1(k)=p1(k)+p2(k)
 137         enddo
 138       enddo
 139       p1(5)=(p1(4)-p1(3))*(p1(4)+p1(3))-p1(2)**2-p1(1)**2
 140       if(p1(5).gt.0.d0)then
 141         p1(5)=sqrt(p1(5))
 142       else
 143         iorrem=iorptl(nk1)
 144         write(*,*)'Precision problem in gakfra, p:',
 145      &             p1(5),p1(4),p1(3),p1(2),p1(1)
 146         write(*,*)'origin : ',iorrem
 147         p1(5)=0.d0
 148 c        if(iorrem.ne.0.and.
 149 c     &    (ityptl(iorrem).eq.40.or.ityptl(iorrem).eq.10))
 150 c     &    p1(5)=dble(pptl(5,iorrem))
 151         if(iorrem.ne.0)then
 152           write(*,*)'string mass : ',ityptl(iorrem),pptl(5,iorrem)
 153           p1(5)=dble(pptl(5,iorrem))
 154         endif
 155       endif
 156       do k=1,5
 157         psg(k)=p1(k)
 158       enddo
 159
 160       if(iappl.ne.6)qsqptl(nk1)=p1(5)*p1(5) !hadronic string energy
 161                                                         !(already defined for lepton)
 162
 163 c.....light like ...............................................
 164
 165       if(ish.ge.6) call blist("before light like&",nk1,nk2)
 166       if(rangen().lt.0.5)then
 167         i1=nk1
 168         i2=nk1+1
 169       else
 170         i1=nk2-1
 171         i2=nk2
 172       endif
 173       ii=1
 174       do while(ii.ge.0)
 175         do i=1,4
 176           p2(i)=dble(pptl(i,i1))+dble(pptl(i,i2))
 177         enddo
 178 c       p2(5)=sqrt(p2(4)**2-p2(3)**2-p2(2)**2-p2(1)**2)
 179         am1=pptl(5,i1)**2
 180         am2=pptl(5,i2)**2
 181         am12=max(0d0,p2(4)**2-p2(3)**2-p2(2)**2-p2(1)**2-am1-am2)/2
 182         if(am12**2.gt.am1*am2)then
 183           ak1=.5*((am2+am12)/sqrt(am12**2-am1*am2))-.5
 184           ak2=.5*((am1+am12)/sqrt(am12**2-am1*am2))-.5
 185         else
 186           ak1=0.
 187           ak2=0.
 188         endif
 189         if(ish.ge.6)write(ifch,'(a,2i4,9f12.6)') 'overlaps'
 190      $       ,i1,i2,ak1,ak2,sqrt(2.*am12)
 191      &       ,am1,am2
 192         do j=1,4
 193           a1=(1.+ak1)*pptl(j,i1)-ak2*pptl(j,i2)
 194           a2=(1.+ak2)*pptl(j,i2)-ak1*pptl(j,i1)
 195           pptl(j,i1)=a1
 196           pptl(j,i2)=a2
 197         enddo
 198         pptl(5,i1)=0.
 199         pptl(5,i2)=0.
 200         if(nk2-nk1.eq.1) ii=-1
 201         ii=ii-1
 202         if (nk1.eq.i1)then
 203           i1=nk2-1
 204           i2=nk2
 205         else
 206           i1=nk1
 207           i2=nk1+1
 208         endif
 209       enddo
 210       if(ish.ge.6) call blist("after  light like&",nk1,nk2)
 211
 212 c.....on-shell and copy ...............................................
 213
 214       if(ish.ge.6) write (ifch,*) 'on-shell check'
 215       do i=nk1,nk2
 216         do k=1,5
 217           p2(k)=dble(pptl(k,i))
 218         enddo
 219         p12=p2(4)
 220         p2(4)=sqrt(p2(3)**2+p2(2)**2+p2(1)**2+p2(5)**2)
 221         if(abs(p12-p2(4))/max(.1d0,p12,p2(4)).ge.1e-1.and.ish.ge.2)then
 222           write(ifch,*)'warning: on-shell problem:'
 223      &           ,i, idptl(i),(pptl(k,i),k=1,4),' (',sngl(p2(4)),') '
 224      &           ,pptl(5,i), istptl(i)
 225         endif
 226         if(ish.ge.6)  write (ifch,*)  p12 ,'->',p2(4)
 227         call utlob2(1,p1(1),p1(2),p1(3),p1(4),p1(5)
 228      $       ,p2(1),p2(2),p2(3),p2(4),60)
 229         f=1.0
 230         if(i.ne.nk1.and.i.ne.nk2)f=0.5
 231         nmax=1
 232         ff=1.
 233         aemax=1000.              ! here max band mass
 234         if(f*p2(4).gt.aemax) then
 235           nmax = int(f*p2(4)/aemax)+1
 236           ff=1./real(nmax)
 237 c          print *,"nmax",nmax
 238         endif
 239         nn=1
 240         do while(nn.le.nmax)
 241           f=.5
 242           if(i.eq.nk1.and.nn.eq. 1  ) f=1.
 243           if(i.eq.nk2.and.nn.eq.nmax) f=1.
 244           nob=nob+1
 245           if(nob.gt.mxnba-7)stop'gakfra: increase parameter mxnba '
 246 c         if(nmax.ge.2) print *, nob,ff,f,i
 247           do k=1,4
 248             pst(k,nob)=ff*f*p2(k)
 249             ipst(nob)=i
 250           enddo
 251           nn=nn+1
 252         enddo
 253       enddo                     !i
 254
 255       do i1=nob-1,1,-1          ! copy again gluons
 256         nob=nob+1
 257         if(nob.gt.mxnba-7)stop'gakfra: increase parameter mxnba '
 258         do k=1,4
 259           pst(k,nob)=pst(k,i1)
 260           ipst(nob)=ipst(i1)
 261         enddo
 262       enddo
 263       nob=nob+1                 ! nob is number of kinks
 264       if(nob.gt.mxnba-7)stop'gakfra: increase parameter mxnba '
 265       if(ish.ge.6) then
 266         write (ifch,*) 'bands:'
 267         do i=0,nob-1
 268           write (ifch,'(i4,4g20.13)') i,(pst(k,i),k=1,4)
 269         enddo
 270       endif
 271
 272 c.....total momentum...............................................
 273       do k=1,4
 274         pst(k,nob)=0.
 275       enddo
 276       do i=0,nob-1
 277         do k=1,4
 278           pst(k,nob)=pst(k,nob)+pst(k,i)
 279         enddo
 280       enddo
 281
 282
 283 c.....pbreak ................................
 284
 285       x=sqrt(qsqptl(nk1))
 286       if(krm.eq.1)then
 287         if(kdrm.eq.0)then
 288           pbi=pbreakr
 289         else
 290           pbi=pbreakd
 291         endif
 292       elseif(kky.eq.0)then
 293         pbi=pbreak
 294       elseif(kky.eq.1)then
 295         pbi=pbreakk
 296       else
 297         stop'in gakfra: unknown kky'
 298       endif
 299       if(pbi.gt.0.d0)then
 300         pb=pbi
 301       else
 302         if(x.lt.14)then
 303 c          pb=0.38
 304          pb=0.4
 305         else
 306          pb=0.14369+ 854.965/x**3 - 111.1144/x**2 + 7.163/x
 307 c         !! {{14,0.4},{22,0.32},{34,0.28},{91.2,0.21}}
 308 c         pb=0.143+ 900/x**3 - 77/x**2 + 4.5/x
 309 c         pb=0.165+ 1020/x**3 - 80/x**2 + 3.8/x
 310 c         pb=0.17+ 950/x**3 - 76/x**2 + 3.8/x
 311 c         pb=0.165+ 600/x**3 - 50/x**2 + 3.8/x
 312 c          pb=0.135+ 855./x**3 - 111./x**2 + 7./x
 313 c          pb=0.05+ 1100./x**3 - 80./x**2 + 5./x
 314         endif
 315         pb=pb*abs(pbi)
 316       endif
 317 c      write(*,*)'pb',kky,krm,kdrm,x,pb,zzrm,zzzz
 318
 319
 320 c.....write total string on pptl - list.............................
 321       nptl=nptl+1
 322       if(nptl.gt.mxptl)call utstop('gakfra: mxptl too small ...&')
 323       call idtr5(idptl(nk1),ic1)
 324       call idtr5(idptl(nk2),ic2)
 325       ic(1)=ic1(1)+ic2(1)
 326       ic(2)=ic1(2)+ic2(2)
 327       idptl(nptl) = 8*10**8+ ic(1)*100 + ic(2)/100 !nbr+1
 328       istptl(nptl)=10
 329       n8ptl=nptl
 330       do i=1,5
 331         pptl(i,nptl)=p1(i)
 332       enddo
 333       if(ish.ge.3)then
 334         write(ifch,'(a)')' ---------total string------------'
 335         call blist('&',nptl,nptl)
 336       endif
 337
 338 c....................................................................
 339       ijb(1,0)=-1
 340       ijb(2,0)=-1
 341       ijb(1,1)=-1
 342       ijb(2,1)=-1
 343       iflb(0)=-idptl(nk1)
 344       iflb(1)= idptl(nk2)
 345       do i=1,4
 346         ptb(i,0)=0.
 347         ptb(i,1)=0.
 348       enddo
 349 c calculate isospin to conserve isospin symmetry in popcorn
 350 c p -> pi+ <=> n -> pi-
 351       if(abs(iflb(0)).gt.6)then
 352         iq01=mod(abs(iflb(0))/100,10)
 353         iq02=abs(iflb(0))/1000
 354         iso=(isospin(iq01)+isospin(iq02))*sign(1,-iflb(0))
 355       else
 356         iso=isospin(abs(iflb(0)))*sign(1,-iflb(0))
 357       endif
 358       if(abs(iflb(1)).gt.6)then
 359         iq11=mod(abs(iflb(1))/100,10)
 360         iq12=abs(iflb(1))/1000
 361         iso=iso+(isospin(iq11)+isospin(iq12))*sign(1,iflb(1))
 362       else
 363         iso=iso+isospin(abs(iflb(1)))*sign(1,iflb(1))
 364       endif
 365
 366 c...................light string....................................
 367       amm=ammin(idptl(nk1),idptl(nk2))
 368
 369       if(sqrt(max(0.,pf(nob,nob))).lt.amm)then
 370         id=idsp(  idptl(nk1),idptl(nk2) )
 371           if(ish.ge.1)then
 372        write (ifch,*)
 373      .       '-------string too light to decay--------'
 374           write (ifch,*) id,p1(5),amm
 375           write (ifch,*) id, sqrt(max(0.,pf(nob,nob))) ,amm
 376        endif
 377         am=sqrt(max(0.,pf(nob,nob)))
 378         call idress(id,am,idr,iadj)
 379        if(ish.ge.1)write (ifch,*) id,am,idr,iadj
 380         nptl=nptl+1
 381         if(nptl.gt.mxptl)
 382      &       call utstop('gakfra (light): mxptl too small ...&')
 383         do i=1,5
 384           pptl(i,nptl)=p1(i)
 385         enddo
 386         do i=nk1,nk2
 387           istptl(i)=21
 388           ifrptl(1,i)=n8ptl
 389           ifrptl(2,i)=0
 390         enddo
 391         istptl(n8ptl)=29        !code for fragmented string
 392         ityptl(n8ptl)=ityptl(nk1)
 393         itsptl(n8ptl)=itsptl(nk1)
 394         rinptl(n8ptl)=-9999
 395         iorptl(n8ptl)=nk1
 396         jorptl(n8ptl)=nk2
 397         ifrptl(1,n8ptl)=n8ptl+1
 398         ifrptl(2,n8ptl)=nptl
 399         ityptl(nptl)=ityptl(nk1)
 400         itsptl(nptl)=itsptl(nk1)
 401         rinptl(nptl)=-9999
 402         istptl(nptl)=0
 403         iorptl(nptl)=n8ptl
 404         jorptl(nptl)=0
 405         ifrptl(1,nptl)=0
 406         ifrptl(2,nptl)=0
 407         if(idr.ne.0)then
 408           idptl(nptl)=idr
 409         else
 410           idptl(nptl)=id
 411         endif
 412         do j=1,4
 413           xorptl(j,nptl)=xorptl(j,nk1)
 414         enddo
 415         tivptl(1,nptl)=xorptl(4,nptl)
 416         call idtau(idptl(nptl),pptl(4,nptl),pptl(5,nptl),taugm)
 417         tivptl(2,nptl)=tivptl(1,nptl)+taugm*(-alog(rangen()))
 418         if(abs(p1(5)-am).gt.0.01)then
 419 c          write (*,*) 'light string  ---  particle off-shell!!!!'
 420 c          write (*,*) idr,'  mass:',p1(5),'  should be:',am
 421         endif
 422         goto 98
 423       endif
 424
 425 c................search breakpoints...................................
 426
 427       n=0
 428       nsbp=1
 429       nbr=0
 430       iok=0
 431        !----------------------!
 432  11      call gaksbp(0,1,iok)
 433        !----------------------!
 434       nbrtry=0
 435       goto 10
 436
 437 c..............delete breakpoint....................................
 438  9    if(ish.ge.4)write(ifch,*) 'delete breakpoint',n,' -> '
 439      &,nbr-1,' breakpoints left'
 440       call gakdel(n)              !hier loeschen
 441
 442 c..............no more breakpoints -> redo..........................
 443       if(nbr.eq.0)then
 444         nsbp=nsbp+1
 445         if(ish.ge.3)write (ifch,*)
 446      &    'no breakpoints left ... try again (',nsbp,')'
 447         if(ish.ge.3)write (ifch,*)' '
 448         goto11
 449       endif
 450
 451 c................make index list of breakpoints to adjust...........
 452  10   continue
 453       nbrtry=nbrtry+1
 454       nind=0
 455       xlastdiq=0.
 456       do i=1,nbr
 457         if(ip(i).eq.0.or.ip(i+1).eq.0)then
 458           nind=nind+1
 459           ind(nind)=i
 460         endif
 461       enddo
 462       if(nbrtry.gt.1000000)then
 463         if(ish.ge.1)then
 464         write(*,*)'Gakfra : string can not be fragmented, redo event !'
 465         write(*,*)nk1,nk2,nbr,nind,pb
 466         endif
 467         iret=1
 468         goto 9999
 469       endif
 470
 471 c.....no more breakpoint to adjust...............................
 472       if(nind.eq.0) goto 100
 473
 474 c................................................................
 475       if(ish.ge.5)then
 476         write(ifch,*) 'breakpoints:'
 477         write(ifch,'(i3,i5)') 0,-iflb(0)
 478         do i=1,nbr
 479           write(ifch,'(i3,2i5,1x,4e12.4,2x,2i3,2f6.3)')
 480      &      i,iflb(i),-iflb(i),(ptb(j,i),j=1,4)
 481      &         ,ijb(1,i),ijb(2,i),xyb(1,i),xyb(2,i)
 482         enddo
 483         write(ifch,'(i3,i5)') nbr+1,iflb(nbr+1)
 484       endif
 485
 486 c.....choose breakpoint, calculate masses, lambda..................
 487       r=rangen()
 488       nn=1+int(r*real(nind))
 489 c      nn=1+(nind-1)*int(r+0.5)
 490       n=ind(nn)
 491       if(ish.ge.4)write(ifch,*) 'choose breakpoint',n
 492       nl=n-1
 493       nr=n+1
 494       do while (nl.gt.0.and.ip(nl).eq.0)
 495          nl=nl-1
 496       enddo
 497       do while (nr.lt.nbr+1.and.ip(nr+1).eq.0)
 498          nr=nr+1
 499       enddo
 500       if(ish.ge.6)write (ifch,*) 'nl,n,nr:',nl,n,nr
 501 c      print *,'------------------------------------',1
 502       call gaksco(n-1,n,n+1,ijb(1,n),ijb(2,n),x1,x2,y1,y2)
 503       if(x2.le.x1.or.y2.le.y1)goto 9
 504 cc      x=(xyb(1,n)-x1)/(x2-x1)
 505 cc      y=(xyb(2,n)-y1)/(y2-y1)
 506       x=xyb(1,n)
 507       y=xyb(2,n)
 508       aml2=co(1)+x*co(2)+y*co(3)+x*y*co(4)
 509       amr2=co(5)+x*co(6)+y*co(7)+x*y*co(8)
 510       ala2=co(9)+x*co(10)+y*co(11)+x*y*co(12)
 511
 512 c.....determine id of both particles..............................
 513       aml=sqrt(max(0.,aml2))
 514       idl=idsp( -iflb(n-1) , iflb(n)   )
 515       amr=sqrt(max(0.,amr2))
 516       idr=idsp( -iflb(n)   , iflb(n+1) )
 517
 518 c.....if mass to small (because of spacelike pt)  reject..........
 519 c      amin=0.0
 520 c      if (aml2.le.amin.or.amr2.le.amin) goto 9
 521
 522 c.....if no left or right ptl id -> reject.........
 523       if(idl.eq.0.or.idr.eq.0)then
 524         if(ish.ge.5)write(ifch,*)'no left or right ptl id'
 525         goto 9
 526       endif
 527
 528       iadjl=0
 529       iadjr=0
 530       idrl=0
 531       idrr=0
 532
 533 c.....if no adjusted mass on left side, check for resonance...........
 534       if(ip(n)  .eq.0) then
 535         aml=sqrt(max(0.,aml2+0.*min(0.,amr2))) !!!????
 536         amlxx=aml
 537         call idress(idl,aml,idrl,iadjl)
 538         r=rangen()
 539         if(idrl.eq.110.and.r.lt.0.5)then
 540           if (r.gt.eta+etap) goto 9      !rare numerical errors
 541           idl=220
 542           aml=.549
 543           if(r.lt.etap)aml=.958
 544           amlxx=aml
 545           call idress(idl,aml,idrl,iadjl)
 546           iadjl=1
 547         endif
 548         if(ish.ge.5)write(ifch,'(a,i5,2f12.6,1i10,i5)')
 549      &    ' l:  ',idl,amlxx,aml,idrl,iadjl
 550       else
 551         if(ish.ge.5)write(ifch,'(a,i5,2f12.6,i10,a5)')
 552      &    ' l:  ',idl,aml,aml,ip(n),'ok'
 553       endif
 554
 555 c.....if no adjusted mass on right side, check for resonance...........
 556       if(ip(n+1).eq.0) then
 557         amr=sqrt(max(0.,amr2+0.*min(0.,aml2))) !!!????
 558         amrxx=amr
 559         call idress(idr,amr,idrr,iadjr)
 560         r=rangen()
 561         if(idrr.eq.110.and.r.lt.0.5)then
 562           if (r.gt.eta+etap) goto 9    !rare numerical errors
 563           idr=220
 564           amr=.549
 565           if(r.lt.etap)amr=.958
 566           amrxx=amr
 567           call idress(idr,amr,idrr,iadjr)
 568           iadjr=1
 569         endif
 570         if(ish.ge.5)write(ifch,'(a,i5,2f12.6,1i10,i5)')
 571      &    ' r:  ',idr,amrxx,amr,idrr,iadjr
 572       else
 573         if(ish.ge.5)write(ifch,'(a,i5,2f12.6,i10,a5)')
 574      &    ' r:  ',idr,amr,amr,ip(n+1),'ok'
 575       endif
 576
 577 c.....mass adjustments.........................................
 578       iok=0
 579       if(ip(n+1).ne.0)then  !.........adjusted mass on right side
 580         if(idrl.eq.0)then
 581           call gaksbp(n-1,n,iok)
 582           if(iok.eq.1)goto 9
 583           goto 10
 584         endif
 585         if(iadjl.eq.1)then
 586            if(ish.ge.5)write(ifch,*)'mass adjustment 1'
 587            n2=n+1
 588            call gakanp(n-1,n,n2,aml**2,amr**2,0.,iok)
 589         endif
 590         if(iok.eq.1)goto 9
 591         ip(n)=idrl
 592       elseif(ip(n).ne.0)then !.........adjusted mass on left side
 593         if(idrr.eq.0)then
 594           call gaksbp(n,n+1,iok)
 595           if(iok.eq.1)goto 9
 596           goto 10
 597         endif
 598         if(iadjr.eq.1)then
 599            if(ish.ge.5)write(ifch,*)'mass adjustment 2'
 600            n2=n+1
 601            call gakanp(n-1,n,n2,aml**2,amr**2,0.,iok)
 602         endif
 603         if(iok.eq.1)goto 9
 604         ip(n+1)=idrr
 605       else       !.........adjusted mass neither left nor right
 606         if(idrr.eq.0.and.idrl.eq.0)then
 607           call gaksbp(n,n+1,iok)
 608           if(iok.eq.1)goto 9
 609           call gaksbp(n-1,n,iok)
 610           if(iok.eq.1)goto 9
 611           goto 10
 612         elseif(idrl.ne.0.and.idrr.eq.0)then
 613           if(iadjl.eq.1) then
 614            if(ish.ge.5)write(ifch,*)'mass adjustment 3'
 615            call gakanp(n-1,n,nr,aml**2,0.,ala2,iok)
 616           endif
 617         elseif(idrl.eq.0.and.idrr.ne.0)then
 618           if(iadjr.eq.1) then
 619            if(ish.ge.5)write(ifch,*)'mass adjustment 4'
 620            n2=n+1
 621            call gakanp(nl,n,n2,0.,amr**2,ala2,iok)
 622           endif
 623         else  !if(idrl.ne.0.and.idrr.ne.0)then
 624           if(iadjl.eq.1.or.iadjr.eq.1) then
 625            if(ish.ge.5)write(ifch,*)'mass adjustment 5'
 626            n2=n+1
 627            call gakanp(n-1,n,n2,aml**2,amr**2,0.,iok)
 628           endif
 629         endif
 630         if(iok.eq.1)goto 9
 631         ip(n)=idrl
 632         ip(n+1)=idrr
 633       endif
 634       if(ish.ge.4)then
 635         write(ifch,*) 'left/right particles:'
 636      &         ,ip(n),aml,'  ',ip(n+1),amr
 637       endif
 638       goto 10
 639
 640 c................................................................
 641 c                         end of string decay
 642 c................................................................
 643
 644 c.....final list...............................................
 645  100  if(ish.ge.4)then
 646         write(ifch,*) '   ************ OK **************'
 647         write(ifch,*) 'final breakpoints:'
 648         write(ifch,'(i3,i5)') 0,-iflb(0)
 649         do i=1,nbr
 650           write(ifch,'(i3,2i5,1x,4e12.4,2x,2i3,2f6.3)')
 651      &      i,iflb(i),-iflb(i),(ptb(j,i),j=1,4)
 652      &         ,ijb(1,i),ijb(2,i),xyb(1,i),xyb(2,i)
 653         enddo
 654         write(ifch,'(i3,i5)') nbr+1,iflb(nbr+1)
 655       endif
 656
 657 c.....write particles in pptl-list................................
 658       if(ish.ge.3)then
 659         write(ifch,'(a)')' ---------produced particles---------'
 660       endif
 661 c to avoid leading strange baryon for diffractive string, rotate string
 662
 663       zz=0.
 664       if(kdrm.eq.1.and.krm.ne.0.and.nbr.gt.0)then
 665         iq1=abs(iflb(1))
 666         iq2=abs(iflb(2))
 667         if(iq1.gt.3.or.iq2.gt.2)then
 668           am=0.
 669         elseif(iq1.eq.3)then
 670           am=qmass(iq1)
 671         elseif(iso.gt.0)then   !for u or d, use proper mass if isospin > 0 (proton)
 672           am=qmass(iq1)
 673         elseif(iso.lt.0)then                  !us d mass for u and u mass for u if isospin < 0
 674           am=qmass(3-iq1)
 675         else
 676           am=0.5*(qmass(1)+qmass(2))
 677         endif
 678         zz=zrminc*pf(nob,nob)**2.
 679         zz=zz*(1.+zrmmax*log(float(matarg+maproj-1)))
 680 c        if(itp/10.eq.4)then
 681 c          zz=zz*(1.+zrmmax*log(max(1.,float(matarg)/float(maproj))))
 682 c        elseif(itp/10.eq.5)then
 683 c          zz=zz*(1.+zrmmax*log(max(1.,float(maproj)/float(matarg))))
 684 c        endif
 685         zz=am*zz
 686 c        print *,iso,itp,pf(nob,nob),am,zz,iq1,iq2,exp(-min(100.,zz))
 687       endif
 688 c      if(kdrm.eq.1)print *,'reminv',itp,am,pf(nob,nob),zz,iq1,iq2
 689 c     &                       ,reminv+exp(-min(100.,zz))
 690       if(rangen().gt.reminv+exp(-min(100.,zz)))then
 691         leadstr=.true.
 692       else
 693         leadstr=.false.
 694       endif
 695
 696       nptlini=nptl
 697       if(nptlini+nbr+1.gt.mxptl)
 698      &       call utstop('gakfra (end): mxptl too small ...&')
 699       do i=1,nbr+1
 700          nptl=nptl+1
 701          if(i.lt.nbr+1)then     !particle = left side of breakpoints
 702            call gaksco(i-1,i,i+1,ijb(1,i),ijb(2,i),x1,x2,y1,y2)
 703            taubrr=pst(4,nob+7)+x*pst(4,nob+8)+y*pst(4,nob+9)
 704            x=xyb(1,i)
 705            y=xyb(2,i)
 706            aml2=co(1)+x*co(2)+y*co(3)+x*y*co(4)
 707            amr2=co(5)+x*co(6)+y*co(7)+x*y*co(8)
 708
 709            ala2=co(9)+x*co(10)+y*co(11)+x*y*co(12)
 710            aml=sign(sqrt(abs(aml2)),aml2)
 711            amr=sign(sqrt(abs(amr2)),amr2)
 712            if(aml.le.0.d0.or.amr.le.0.d0)then
 713              if(ish.ge.4)write(ifch,*)
 714      & 'Negative mass, fragmentation not OK -> continue ...'
 715              n=i
 716              nptl=nptlini
 717              goto 9
 718            endif
 719            qsqptl(nptl)=ala2
 720            pptl(5,nptl)=aml
 721            do j=1,4
 722              pptl(j,nptl)=pst(j,nob+1)-
 723      &            x*pst(j,nob+2)+y*pst(j,nob+3)
 724 c             pptr(j,nptl)=ptb(j,i)
 725            enddo
 726          else                   !last particle = right side of last breakpoint
 727            pptl(5,nptl)=amr
 728            qsqptl(nptl)=0.
 729            do j=1,4
 730              pptl(j,nptl)=pst(j,nob+4)+
 731      &            x*pst(j,nob+5)-y*pst(j,nob+6)
 732 c             pptr(j,nptl)=ptb(j,i) !should be zero
 733            enddo
 734            taubrr=0.
 735          endif
 736          idptl(nptl)=ip(i)
 737
 738       enddo
 739
 740       if(leadstr)then
 741         iptmp=idptl(nptlini+2)
 742         idptl(nptlini+2)=idptl(nptlini+1)
 743         idptl(nptlini+1)=iptmp
 744         pptltmp=pptl(5,nptlini+2)
 745         pptl(5,nptlini+2)=pptl(5,nptlini+1)
 746         pptl(5,nptlini+1)=pptltmp
 747         pptl(4,nptlini+1)=pptl(5,nptlini+1)*pptl(5,nptlini+1)
 748         pptl(4,nptlini+2)=pptl(5,nptlini+2)*pptl(5,nptlini+2)
 749         do k=1,3
 750           pptl(4,nptlini+1)=pptl(4,nptlini+1)
 751      &                     +pptl(k,nptlini+1)*pptl(k,nptlini+1)
 752           pptl(4,nptlini+2)=pptl(4,nptlini+2)
 753      &                     +pptl(k,nptlini+2)*pptl(k,nptlini+2)
 754         enddo
 755         pptl(4,nptlini+1)=sqrt(pptl(4,nptlini+1))
 756         pptl(4,nptlini+2)=sqrt(pptl(4,nptlini+2))
 757       endif
 758       nptl=nptlini
 759       do i=1,nbr+1
 760          nptl=nptl+1
 761
 762          call utlob4(-1,p1(1),p1(2),p1(3),p1(4),p1(5)
 763      $        ,pptl(1,nptl),pptl(2,nptl),pptl(3,nptl),pptl(4,nptl))
 764 c         call utlob4(-1,p1(1),p1(2),p1(3),p1(4),p1(5)
 765 c     $        ,pptr(1,nptl),pptr(2,nptl),pptr(3,nptl),pptr(4,nptl))
 766
 767 c........Origin..................................................
 768          istptl(nptl)=0
 769          iorptl(nptl)=n8ptl
 770          jorptl(nptl)=0
 771          ifrptl(1,nptl)=0
 772          ifrptl(2,nptl)=0
 773
 774          r=rangen()
 775          tauran=-taurea*alog(r)*pptl(4,nptl)/pptl(5,nptl)
 776 c         if(jpsi.ge.1)tauran=0.
 777 c         tauran=max(taubrl,taubrr) !take formation time from string-theory
 778          do j=1,4
 779            xorptl(j,nptl)=xorptl(j,nk1)+pptl(j,nptl)
 780      &       /pptl(4,nptl)*tauran
 781          enddo
 782          tivptl(1,nptl)=xorptl(4,nptl)
 783          call idtau(idptl(nptl),pptl(4,nptl),pptl(5,nptl),taugm)
 784          tivptl(2,nptl)=tivptl(1,nptl)+taugm*(-alog(rangen()))
 785          ityptl(nptl)=ityptl(nk1)
 786          itsptl(nptl)=itsptl(nk1)
 787          rinptl(nptl)=-9999
 788          if(ish.ge.3)call blist('&',nptl,nptl)
 789 ctp060829          taubrl=taubrr
 790       enddo
 791       iorptl(n8ptl)=nk1
 792       jorptl(n8ptl)=nk2
 793       ityptl(n8ptl)=ityptl(nk1)
 794       do i=nk1,nk2
 795          istptl(i)=21
 796          ifrptl(1,i)=n8ptl
 797          ifrptl(2,i)=0
 798       enddo
 799       istptl(n8ptl)=29          !code for fragmented string
 800       ifrptl(1,n8ptl)=n8ptl+1
 801       ifrptl(2,n8ptl)=nptl
 802       if(ish.ge.5)then
 803         write(ifch,*)'string momentum sum:'
 804         do kk=1,5
 805           pptl(kk,nptl+1)=0
 806           do ii=n8ptl+1,nptl
 807             pptl(kk,nptl+1)=pptl(kk,nptl+1)+pptl(kk,ii)
 808           enddo
 809         enddo
 810         call alist2('&',n8ptl,n8ptl,nptl+1,nptl+1)
 811       endif
 812
 813
 814
 815 c.....another string?.........................................
 816  98   nk1=nk2+1
 817       goto 2                    !next string
 818
 819 c.....end of fragmentation.....................................
 820  9999 continue
 821       call utprix('gakfra',ish,ishini,4)
 822       return
 823       end
 824
 825 c---------------------------------------------------------------------
 826       subroutine gaksbp(ibr1,ibr2,iok)
 827 c---------------------------------------------------------------------
 828       ! search break points
 829       !-----------------------------------------
 830       ! nbr ... number of break points
 831       !
 832       !
 833       !
 834       !--------------------------------------------------------------
 835       include 'epos.inc'
 836
 837       parameter (mxnbr=500,mxnba=5000)
 838       common /gag/nob,pst(4,0:mxnba),ipst(0:mxnba)
 839      $     ,nbr,ijb(2,0:mxnbr),xyb(2,0:mxnbr)
 840      &     ,ptb(4,0:mxnbr),iflb(0:mxnbr),ip(0:mxnbr),co,x,y
 841       common/pb/pb /cnsbp/nsbp /cn8ptl/n8ptl
 842       double precision ax,ay,az,ae,am,bx,by,bz,be
 843       dimension co(12)
 844       logical weiter
 845       common/czz/kky,krm,zzod,zzop,kdrm,zzoq,zzos,zzrm,zzipp,zzipt
 846      &,zzid,zzip,zzzzr,itp
 847       double precision psg
 848       common/zpsg/psg(5)
 849       dimension pleft(5),pright(5)
 850
 851       pf(i,j)=pst(4,i)*pst(4,j)-pst(3,i)*pst(3,j)
 852      &     -pst(2,i)*pst(2,j)-pst(1,i)*pst(1,j)
 853       mmod(i)=mod(mod(i,nob)+nob,nob)
 854
 855       call utpri('gaksbp',ish,ishini,5)
 856
 857
 858       ib1=ibr1
 859       ib2=ibr2
 860       iside=1
 861       if(rangen().lt.0.5)iside=2
 862       if(ish.ge.6)write(ifch,*)'iside:',iside
 863       Amxx=80./pb
 864       if(ish.ge.4)write(ifch,*)
 865      &'search brk between ib1=',ib1,' and ib2=',ib2
 866       ntry=0
 867       nbrold=nbr
 868  24   Amin=0.
 869  25   Amax=Amxx
 870  26   ntry=ntry+1
 871       A0=-log(exp(-pb*Amin)+rangen()*(exp(-pb*Amax)-exp(-pb*Amin)))/pb
 872       if(ish.ge.6)write(ifch,*)'pb, Amin, Amax, A0:',pb, Amin, Amax, A0
 873       A=A0
 874       Amaxn=0.
 875       if(iside.eq.2)goto 51
 876 c.....................................................................
 877       if(ib1.eq.0)then          !startwert der j-schleife
 878          jj=nob-1
 879       else
 880          jj=ijb(2,ib1)
 881       endif
 882       do while (jj.gt.0)
 883          if(jj.eq.ijb(2,ib1) )then       !linker brk im Streifen jj?
 884             y1=xyb(2,ib1)
 885          else                   !nein
 886             y1=0.
 887          endif
 888          if(jj.eq.ijb(2,ib2))then       !rechter brk im Streifen jj?
 889             y2=xyb(2,ib2)
 890          else                   !nein
 891             y2=1.
 892          endif
 893          if(y1.eq.y2) goto 9999
 894          if(abs(y1-y2)/abs(y1+y2).le.1e-7) goto 9999
 895          if(ish.ge.6)write(ifch,*)'y1,y2,A',y1,y2,A
 896                                 !startwert der i-schleife
 897          if(ib1.eq.0)then       !ohne linken brkpt
 898             ii=mmod(jj+1)
 899             if(jj.lt.nob/2)ii=mmod(nob-jj+1)
 900             if(ib2.le.nbr.and.mmod(ii+nob/2)
 901      &        .gt.mmod(ijb(1,ib2)+nob/2))then
 902               if(ish.ge.6)write(ifch,*) 'very special case',ii,jj
 903               goto12
 904             endif
 905          else
 906             ii=ijb(1,ib1)
 907             if(jj.lt.nob/2 .and.
 908      $           mmod(nob-jj+1+nob/2).gt. mmod(ii+nob/2)
 909      $           )ii=mmod(nob-jj+1)
 910          endif
 911          weiter=.true.
 912          aa=0.
 913          if(ii.eq.jj) then
 914             if(ish.ge.6)write(ifch,*) 'Rand erreicht'
 915             goto 15
 916          endif
 917          do while(weiter)
 918             if(ii.eq.ijb(1,ib1))then    !linker brk im Feld (ii,jj)
 919                x2=xyb(1,ib1)
 920             else
 921                x2=1.
 922             endif
 923
 924             if(ii.eq.ijb(1,ib2))then    !rechter brk im Feld (ii,jj)
 925                x1=xyb(1,ib2)
 926             else
 927                x1=0.
 928             endif
 929             if(x1.eq.x2) goto 9999
 930             if(abs(x1-x2)/abs(x1+x2).le.1e-7) goto 9999
 931             f=1.0
 932             if(ipst(ii).ne.ipst(jj))aa=aa+2.*(x2-x1)*(y2-y1)*f*pf(ii,jj)
 933             if(ish.ge.6)write(ifch,*)'ii,jj,x1,x2,aa:',ii,jj
 934      $           ,ipst(ii).ne.ipst(jj),x1,x2,aa
 935      &           ,pf(ii,jj)
 936             if(ii.eq.ijb(1,ib2))then
 937                weiter=.false.
 938             else
 939                ii=mmod(ii+1)
 940                if(ii.eq.jj.or.mmod(ii+jj).eq.0)weiter=.false.
 941             endif
 942          enddo
 943          Amaxn=Amaxn+aa
 944          if(aa.gt.A)goto 10
 945          A=A-aa
 946          if(jj.eq.ijb(2,ib2)) then
 947             if(ish.ge.6)write(ifch,*) 'brk erreicht'
 948             goto 15
 949          endif
 950  12      jj=mmod(jj-1)
 951       enddo
 952       goto 15
 953
 954  10   continue
 955       yb=A/aa*(y2-y1)+y1
 956       if(ish.ge.6)write(ifch,*)'jj,yb ok:',jj,yb
 957       r=rangen()
 958       ra=aa*r
 959       if(ish.ge.6)write(ifch,*)'r,ra,aa',r,ra,aa
 960                                 !nochmal die letzte ii-Schleife
 961       if(ib1.eq.0)then          !ohne linken brkpt
 962          ii=mmod(jj+1)
 963          if(jj.lt.nob/2)ii=mmod(nob-jj+1)
 964       else
 965         ii=ijb(1,ib1)
 966         if(jj.lt.nob/2 .and.
 967      $       mmod(nob-jj+1+nob/2).gt. mmod(ii+nob/2)
 968      $       )ii=mmod(nob-jj+1)
 969       endif
 970       do while(ra.gt.0.)
 971          if(ii.eq.ijb(1,ib1))then       !linker brk im Feld (ii,jj)
 972             x2=xyb(1,ib1)
 973          else
 974             x2=1.
 975          endif
 976
 977          if(ii.eq.ijb(1,ib2))then       !rechter brk im Feld (ii,jj)
 978             x1=xyb(1,ib2)
 979          else
 980             x1=0.
 981          endif
 982          f=1.0
 983          ab=0.
 984          if(ipst(ii).ne.ipst(jj)) ab=2.*(x2-x1)*(y2-y1)*f*pf(ii,jj)
 985          if(ish.ge.6)write(ifch,*)'ii,jj,x1,x2,ab,ra:',ii,jj
 986      $        ,ipst(ii).ne.ipst(jj),x1,x2,ab,ra
 987          if(ab.gt.ra)then
 988             xb=ra/ab*(x2-x1)+x1
 989          else
 990             ii=mmod(ii+1)
 991          endif
 992          ra=ra-ab
 993       enddo
 994       if(ish.ge.5)write(ifch,*) 'breakpoint in field ',ii,jj
 995      & ,' at position ',xb,yb
 996       goto 95
 997
 998 c......................................................................
 999       !von rechts
1000  51   if(ib2.eq.nbr+1)then      !startwert der i-schleife
1001          ii=nob/2-1
1002       else
1003          ii=ijb(1,ib2)
1004       endif
1005       do while (ii.ne.nob/2)
1006          if(ii.eq.ijb(1,ib1))then !linker brk im Streifen (ii)
1007             x2=xyb(1,ib1)
1008          else
1009             x2=1.
1010          endif
1011          if(ii.eq.ijb(1,ib2))then !rechter brk im Streifen (ii)
1012             x1=xyb(1,ib2)
1013          else
1014             x1=0.
1015          endif
1016          if(x1.eq.x2) goto 9999
1017          if(abs(x1-x2)/abs(x1+x2).le.1e-7) goto 9999
1018          if(ish.ge.6)write(ifch,*)'x1,x2 A',x1,x2,A
1019
1020          if(ib2.eq.nbr+1)then
1021             jj=mmod(ii+1)
1022             if(ii.gt.nob/2)jj=mmod(nob-ii+1)
1023             if(ib1.ge.1.and.jj.gt.ijb(2,ib1))then
1024                if(ish.ge.6)write(ifch,*) 'very special case',ii,jj
1025                goto 13
1026             endif
1027          else
1028             jj=ijb(2,ib2)
1029             if(ii.gt.nob/2 .and. mmod(nob-ii+1).gt.jj)jj=mmod(nob-ii+1)
1030          endif
1031          weiter=.true.
1032          aa=0.
1033          if(ii.eq.jj) then
1034             if(ish.ge.6)write(ifch,*) 'Rand erreicht'
1035             goto 15
1036          endif
1037          do while(weiter)
1038             if(jj.eq.ijb(2,ib1))then    !linker brk im Feld (ii,jj)
1039                y1=xyb(2,ib1)
1040             else
1041                y1=0.
1042             endif
1043             if(jj.eq.ijb(2,ib2))then !rechter brk im Feld (ii,jj)
1044                y2=xyb(2,ib2)
1045             else
1046                y2=1.
1047             endif
1048             if(y1.eq.y2) goto 9999
1049             if(abs(y1-y2)/abs(y1+y2).le.1e-7) goto 9999
1050             f=1.0
1051             if(ipst(ii).ne.ipst(jj))aa=aa+2.*(x2-x1)*(y2-y1)*f*pf(ii,jj)
1052             if(ish.ge.6)write(ifch,*)'ii,jj,x1,x2,aa:',ii,jj
1053      $           ,ipst(ii).ne.ipst(jj),x1,x2,aa
1054      &           ,pf(ii,jj)
1055             if(jj.eq.ijb(2,ib1))then
1056                weiter=.false.
1057             else
1058                jj=mmod(jj+1)
1059                if(jj.eq.ii.or.mmod(ii+jj).eq.0)weiter=.false.
1060             endif
1061          enddo
1062          Amaxn=Amaxn+aa
1063          if(aa.gt.A)goto 14
1064          A=A-aa
1065          if(ii.eq.ijb(1,ib1)) then
1066             if(ish.ge.6)write(ifch,*) 'brk erreicht'
1067             goto 15
1068          endif
1069  13      ii=mmod(ii-1)
1070       enddo
1071       goto 15
1072
1073
1074
1075  14   continue
1076       xb=A/aa*(x2-x1)+x1
1077       if(ish.ge.6)write(ifch,*)'ii,xb ok:',ii,xb
1078       r=rangen()
1079       ra=aa*r
1080       if(ish.ge.6)write(ifch,*)'r,ra,aa',r,ra,aa
1081                                 !nochmal die letzte jj-Schleife
1082       if(ib2.eq.nbr+1)then
1083          jj=mmod(ii+1)
1084          if(ii.gt.nob/2)jj=mmod(nob-ii+1)
1085       else
1086          jj=ijb(2,ib2)
1087          if(ii.gt.nob/2 .and. mmod(nob-ii+1).gt.jj)jj=mmod(nob-ii+1)
1088       endif
1089
1090       do while(ra.gt.0.)
1091          if(jj.eq.ijb(2,ib1))then       !linker brk im Feld (ii,jj)
1092             y1=xyb(2,ib1)
1093          else
1094             y1=0.
1095          endif
1096
1097          if(jj.eq.ijb(2,ib2))then       !rechter brk im Feld (ii,jj)
1098             y2=xyb(2,ib2)
1099          else
1100             y2=1.
1101          endif
1102          f=1.0
1103          ab=0.
1104          if(ipst(ii).ne.ipst(jj)) ab=2.*(x2-x1)*(y2-y1)*f*pf(ii,jj)
1105          if(ish.ge.6)write(ifch,*)'ii,jj,x1,x2,ab,ra:',ii,jj
1106      $        ,ipst(ii).ne.ipst(jj),x1,x2,ab,ra
1107          if(ab.gt.ra)then
1108             yb=ra/ab*(y2-y1)+y1
1109          else
1110             jj=mmod(jj+1)
1111          endif
1112          ra=ra-ab
1113       enddo
1114       if(ish.ge.5)write(ifch,*) 'breakpoint in field ',ii,jj
1115      & ,' at position ',xb,yb
1116
1117  95   continue
1118
1119 c.....breakpoint accepted......................
1120       nbr=nbr+1
1121       if(nbr.gt.mxnbr-2) stop 'gaksbp: increase parameter mxnbr'
1122       do i=nbr+1,ib1+1,-1
1123          do j=1,2
1124             ijb(j,i)=ijb(j,i-1)
1125             xyb(j,i)=xyb(j,i-1)
1126          enddo
1127          do k=1,4
1128             ptb(k,i)=ptb(k,i-1)
1129          enddo
1130          iflb(i)=iflb(i-1)
1131          ip(i)=ip(i-1)
1132       enddo
1133       ip(ib1+1)=0
1134       ip(ib1+2)=0
1135       ijb(1,ib1+1)=ii
1136       ijb(2,ib1+1)=jj
1137       xyb(1,ib1+1)=xb
1138       xyb(2,ib1+1)=yb
1139
1140       if(kdrm.eq.1)then
1141         f=2.
1142         if(ib1.eq.0)f=1.
1143         if(abs(iflb(ib1)).gt.3)then
1144           iq1=mod(abs(iflb(ib1))/1000,10)
1145           iq2=mod(abs(iflb(ib1))/100,10)
1146           amlast=qmass(0)+f*(qmass(iq1)+qmass(iq2))
1147         else
1148           amlast=f*qmass(abs(iflb(ib1)))
1149         endif
1150         if(itp.ge.50)then
1151           ampt=amtarg
1152         elseif(itp.ge.40)then
1153           ampt=amproj
1154         endif
1155         xlastdiq=0.
1156         do i=0,ib1-1
1157           f=2.
1158           if(i.eq.0)f=1.
1159           if(abs(iflb(i)).gt.3)then
1160             iq1=mod(abs(iflb(i))/1000,10)
1161             iq2=mod(abs(iflb(i))/100,10)
1162             xlastdiq=xlastdiq+qmass(0)+f*(qmass(iq1)+qmass(iq2))
1163           else
1164             xlastdiq=xlastdiq+f*qmass(abs(iflb(i)))
1165           endif
1166         enddo
1167         if(pf(nob,nob).le.5.)amlast=1.e-10 !if very low energy, strangeness should not be totally killed
1168 c        if(pf(nob,nob).le.5.)then
1169 c          xlastdiq=-(ampt+2.*qmass(3)) !if very low energy, strangeness should not be totally killed
1170 c          amlast=strcut**2
1171 c        endif
1172        if(ish.ge.5)write(ifch,*) 'diff cut:',pf(nob,nob),xlastdiq,amlast
1173      &  ,float(ib1+1)-ampt-xlastdiq+diqcut*sqrt(pf(nob,nob))
1174      &  ,strcut/amlast/ampt*pf(nob,nob)
1175       endif
1176
1177 c ..........diquark...............................................
1178
1179
1180       if(krm.eq.1)then
1181         if(kdrm.eq.0)then
1182           pdiqu=pdiquar*zzod
1183          else
1184           pdiqu=pdiquad
1185         endif
1186       elseif(kky.eq.1)then
1187         pdiqu=pdiquak*zzid
1188       else
1189         pdiqu=pdiqua!*zzid
1190       endif
1191
1192       if(kdrm.eq.1)pdiqu=pdiqu*(1.-exp(-min(100.,
1193      & strcut/amlast/qmass(0)/ampt
1194      & *max(0.,float(ib1+1)-ampt-xlastdiq-2.*qmass(0)+diqcut*
1195      &                      sqrt(pf(nob,nob)))*pf(nob,nob))))
1196 c      if(kdrm.eq.1.and.sqrt(pf(nob,nob)).lt.diqcut)pdiqu=0.
1197
1198       if(rangen().lt.pdiqu.and.iabs(iflb(ib1)).lt.6
1199      &  .and.iabs(iflb(ib2)).lt.6)then
1200         jqu=2
1201       else
1202         jqu=1
1203       endif
1204 c ..........flavor...............................................
1205   77  continue
1206       if(krm.eq.1)then
1207         if(kdrm.eq.0)then
1208           pud1= ( 1 - min(0.33,zzos*(1-2*pudr)) ) /2.
1209 c          pud2= pud1
1210         else
1211           pud1= pudd
1212 c          pud2= pud1
1213         endif
1214       elseif(kky.eq.1)then
1215       pud1=pudk
1216 c     pud2=pud1
1217       else
1218       pud1= pud
1219 c     pud2= pud1
1220       endif
1221       if(kdrm.eq.1)pud1=0.5*(1.-(1.-2.*pud1)
1222      &  *(1.-exp(-min(100.,strcut/amlast/qmass(3)/ampt
1223      &  *max(0.,float(ib1+1)-ampt-xlastdiq-2.*qmass(3)+diqcut
1224      &                      *sqrt(pf(nob,nob)))*pf(nob,nob)))))
1225 c      if(kdrm.eq.1.and.sqrt(pf(nob,nob)).lt.strcut)pud1=0.5
1226
1227       pud2=(1.-(1.-2.*pud1)*puds)*0.5   !from e+e- data
1228
1229       if(iabs(iflb(ib1)).gt.1000)then
1230         ifl1=int(rangen()**0.625/pud1)+1
1231       else
1232         ifl1=int(rangen()**1.15/pud1)+1     !assymmetric u and d relevant because of assymmetry Kch / K0 in e+e-
1233       endif
1234
1235 c      amk0=2.*qmass(ifl1)                 !mass for mt distribution
1236       if(ifl1.eq.4)ifl1=3
1237       ifl2=0
1238       if(jqu.eq.2)then  !diquark ------
1239         ifl2=int(rangen()**0.625/pud2)+1
1240 c        amk0=amk0+2.*qmass(ifl2)+qmass(0)  !mass for mt distribution with bounding energy for diquark (qmass(0))
1241         if(ifl2.eq.4)ifl2=3
1242         ifl=-min(ifl1,ifl2)*1000-max(ifl1,ifl2)*100
1243       else !quark ------
1244         ifl=ifl1
1245       endif
1246       if(iflb(ib1).lt.0.and.iflb(ib1).ge.-6)ifl=-ifl
1247       if(iflb(ib1).gt.1000)ifl=-ifl
1248       iflb(ib1+1)=ifl
1249       if(ish.ge.5)write(ifch,*) 'flavor,pud1/2:',ifl,pud1,pud2,pdiqu
1250 c..............................pt.......................................
1251       if(krm.ne.1.and.kdrm.ne.1)then
1252         icub=ib1+1
1253         !---------------------------------------------------------------------------
1254         !  ib1+1 is the current break point index
1255         !             (between 1 and nbr)
1256         !  ijb(1,ib1) and ijb(2,ib1) are band indices
1257         !         each index from 0 to nob-1 (nob= number of bands)
1258         !         0 is left, then come the gluons, then antiquark, then agin gluons
1259         !         like q - g - g - g - ~q - g - g - g
1260         !--------------------------------------------------------------------------
1261         call gaksco(icub-1,icub,icub+1
1262      &    ,ijb(1,icub),ijb(2,icub),x1,x2,y1,y2)
1263         x=xyb(1,icub)
1264         y=xyb(2,icub)
1265         aml2=co(1)+x*co(2)+y*co(3)+x*y*co(4)
1266         amr2=co(5)+x*co(6)+y*co(7)+x*y*co(8)
1267         aml=sign(sqrt(abs(aml2)),aml2)
1268         amr=sign(sqrt(abs(amr2)),amr2)
1269         !------segment left of current breakpoint icub -----
1270         pleft(5)=aml
1271         do j=1,4
1272           pleft(j)=pst(j,nob+1)-x*pst(j,nob+2)+y*pst(j,nob+3)
1273         enddo
1274         call utlob4(-1,psg(1),psg(2),psg(3),psg(4),psg(5)
1275      $     ,pleft(1),pleft(2),pleft(3),pleft(4))
1276         !------segment right of current breakpoint icub-----
1277         pright(5)=amr
1278         do j=1,4
1279           pright(j)=pst(j,nob+4)+x*pst(j,nob+5)-y*pst(j,nob+6)
1280         enddo
1281         call utlob4(-1,psg(1),psg(2),psg(3),psg(4),psg(5)
1282      $     ,pright(1),pright(2),pright(3),pright(4))
1283         !-------------------------
1284         amt=pleft(5)**2+pleft(1)**2+pleft(2)**2
1285         if(amt.gt.0..and.pleft(4)+abs(pleft(3)).gt.0.d0)then
1286           amt=sqrt(amt)
1287           yleft=sign(1.,pleft(3))*alog((pleft(4)+abs(pleft(3)))/amt)
1288         else
1289           yleft=0.                  !
1290         endif
1291         amt=pright(5)**2+pright(1)**2+pright(2)**2
1292         if(amt.gt.0..and.pright(4)+abs(pright(3)).gt.0.d0)then
1293           amt=sqrt(amt)
1294           yright=sign(1.,pright(3))*alog((pright(4)+abs(pright(3)))/amt)
1295         else
1296           yright=0.                  !
1297         endif
1298         ybrk=(yleft+yright)/2.
1299         yhax=2                  !0.5*yhaha
1300         zzipx=1
1301         if(ybrk.gt.yhax)then
1302           zzipx=zzipx+(zzipp-1)
1303         elseif(ybrk.gt.-yhax)then
1304           zzipx=zzipx+(zzipp-1)*(ybrk+yhax)/(2*yhax)
1305         endif
1306         if(ybrk.lt.-yhax)then
1307           zzipx=zzipx+(zzipt-1)
1308         elseif(ybrk.lt.yhax)then
1309           zzipx=zzipx+(zzipt-1)*(yhax-ybrk)/(2*yhax)
1310         endif
1311       endif
1312       delptfra=0
1313       if(ifl1.eq.3.and.ifl2.eq.0)delptfra=delptfra+ptfrasr
1314 c      if(ifl1.eq.3.or.ifl2.eq.3)delptfra=delptfra+ptfrasr
1315       fnsbp=1
1316       if(nsbp.gt.9)fnsbp=0
1317
1318 c for a better pt distribution, we generate mt=sqrt(pt2+m02) (with m0=amk0) instead of pt (do not give good results !)
1319       if(krm.eq.1)then
1320         if(kdrm.eq.0)then
1321 c          ptadd=ranpt()*(ptfrair+delptfra)*(zzop-1)*fnsbp  !*jqu
1322 c          ptadd=ranpt()*ptfrair*(zzop-1)*fnsbp  !*jqu
1323           ptadd=(ptfrair+delptfra)*(zzop-1)*fnsbp  !*jqu
1324           if(jqu.eq.1)then
1325 c            pt=ranpt()*(ptfrair +delptfra )  + ptadd
1326 c            pt=ranpt()*ptfrair
1327             pt=ranpt()*(ptfrair              + ptadd )
1328           else
1329             pt=ranpt()*(ptfrair   + ptadd          )
1330 c            pt=ranpt()*ptfraqq
1331 c            pt=ranpt()*(ptfraqq    + ptadd )
1332           endif
1333         else
1334 c          ptadd=0.
1335            pt=ranpt()*ptfradr
1336           if(jqu.eq.1)then
1337 c            pt=ranpt()*(ptfradr +delptfra   )
1338             pt=ranpt()*ptfradr
1339           else
1340             pt=ranpt()*ptfraqq
1341 c            pt=ranpt()*(ptfradr       )
1342           endif
1343         endif
1344       elseif(kky.eq.1)then
1345 c        ptadd=ranptk()*(ptfrak+delptfra)*(zzipx-1)*fnsbp  !*jqu
1346         pttra=(ptfrak+delptfra)*(zzipx-1)*fnsbp    !/jqu
1347         if(jqu.eq.1)then
1348             !pt=ranptk()*(ptfrak+delptfra    +  pttra )
1349             !pt=ranptk()*(ptfrak   +  pttra)
1350          pt=ranptk()*ptfrak  +  ranptk()*pttra
1351         else
1352             !pt=ranptk()*(ptfrakd   +  pttra)
1353          pt=ranptk()*ptfrakd  +  ranptk()*pttra
1354         endif
1355       else
1356 c        ptadd=ranpt()*(ptfra+delptfra)*(zzip-1)*fnsbp     !*jqu
1357 c        pttra=ranpt()*ptfra*(zzip-1)*fnsbp     !*jqu
1358         pttra=(ptfra+delptfra)*(zzip-1)*fnsbp     !*jqu
1359         if(jqu.eq.1)then
1360 c          pt=ranpt()*(ptfra+delptfra )   +  pttra
1361 c          pt=ranpt()*ptfra
1362           pt=ranpt()*(ptfra          +  pttra )
1363         else
1364 c          pt=ranpt()*(ptfra+delptfra)    +  pttra
1365 c          pt=ranpt()*ptfrakd
1366           pt=ranpt()*(ptfra         +  pttra )
1367         endif
1368       endif
1369 c      pt=sqrt(pt*pt+2.*pt*amk0)+ptadd                  !sample mt-m0 instead of pt ... doesn't work at RHIC
1370       beta=2.*pi*rangen()
1371  60   continue
1372       if(ish.ge.5)then
1373         write(ifch,*) 'pt:',pt
1374       endif
1375       ptb(1,ib1+1)=pt*cos(beta)
1376       ptb(2,ib1+1)=pt*sin(beta)
1377       ptb(3,ib1+1)=0.
1378       ptb(4,ib1+1)=0.
1379       if(ish.ge.8)then
1380         write(ifch,*) 'the bands'
1381         write(ifch,'(4g12.6)') (pst(i,ii),i=1,4)
1382         write(ifch,'(4g12.6)') (pst(i,jj),i=1,4)
1383         write(ifch,*) 'pt before rotation'
1384         write(ifch,'(4f12.8)') (ptb(i,ib1+1),i=1,4)
1385       endif
1386       ax=dble(pst(1,ii))+dble(pst(1,jj))
1387       ay=dble(pst(2,ii))+dble(pst(2,jj))
1388       az=dble(pst(3,ii))+dble(pst(3,jj))
1389       ae=dble(pst(4,ii))+dble(pst(4,jj))
1390       am=sqrt(max(1d-10,(ae+az)*(ae-az)-ax**2-ay**2)) !?????????
1391       if(ish.ge.8)then
1392         write(ifch,*) 'boost vector ( region ) '
1393         write(ifch,'(5g12.6)') ax,ay,az,ae,am,pf(ii,jj)
1394       endif
1395       bx=pst(1,ii)
1396       by=pst(2,ii)
1397       bz=pst(3,ii)
1398       be=pst(4,ii)
1399       call utlob2(1,ax,ay,az,ae,am,bx,by,bz,be,60)
1400       if(ish.ge.8) then
1401         write (ifch,*) 'boost of b'
1402         write (ifch,*) 'bx,by,bz,be',bx,by,bz,be
1403       endif
1404       call utrot4(-1,bx,by,bz,ptb(1,ib1+1),ptb(2,ib1+1),ptb(3,ib1+1))
1405       if(ish.ge.8) then
1406         write (ifch,*) 'rot of pt'
1407         write(ifch,'(4f12.8)') (ptb(i,ib1+1),i=1,4)
1408       endif
1409
1410       call utlob4(-1,ax,ay,az,ae,am
1411      &  ,ptb(1,ib1+1),ptb(2,ib1+1),ptb(3,ib1+1),ptb(4,ib1+1))
1412
1413       if(ish.ge.8) then
1414         write (ifch,*) 'backboost of pt'
1415         write(ifch,'(4f12.8)') (ptb(i,ib1+1),i=1,4)
1416       endif
1417       if(az.eq.0..and.ay.eq.0.)az=1e-8
1418       if(ish.ge.8)write(ifch,*)'rot vector:',ax,ay,az,ae,am
1419
1420 c.....call utrota(-1,ax,ay,az,ptb(1,ib1+1),ptb(2,ib1+1),ptb(3,ib1+1))
1421
1422       if(ish.ge.6)then
1423         write(ifch,*) 'pt'
1424         write(ifch,'(4g12.6)') (ptb(i,ib1+1),i=1,4)
1425         write (ifch,*) ijb(1,ib1+1),ijb(2,ib1+1),xyb(1,ib1+1)
1426      $       ,xyb(2,ib1+1)
1427       endif
1428
1429       if(iside.eq.1)then
1430          ib1=ib1+1
1431          ib2=ib2+1
1432       endif
1433
1434
1435 c      Amin=0.
1436 c      if(Amax.lt.Amxx) goto 15
1437 c      goto 25
1438
1439  15   continue
1440       if(ish.ge.6)write(ifch,*) 'Amax:',Amax,Amaxn
1441       if (nbr.eq.nbrold) then
1442          Amax=Amaxn
1443          Amin=0.
1444          if(pb*Amax.le.0..or.ntry.ge.1000)then
1445            goto 9999
1446          endif
1447          goto 26
1448       endif
1449
1450       if(ish.ge.6)then
1451          write(ifch,*) 0,iflb(0)
1452          do i=1,nbr
1453            if(i.eq.ib2) write(ifch,*) '.................'
1454             write(ifch,'(i3,2x,2(i3),2(" ",g14.7),3x,i5,4(" ",g12.6))')
1455      &           i,ijb(1,i),ijb(2,i),xyb(1,i),xyb(2,i)
1456      &           ,iflb(i),(ptb(j,i),j=1,4)
1457            if(i.eq.ibr1) write(ifch,*) '.................'
1458          enddo
1459          write(ifch,*) nbr+1,iflb(nbr+1)
1460       endif
1461
1462  9999 if(nbr.eq.nbrold)iok=1
1463       call utprix('gaksbp',ish,ishini,5)
1464       return
1465       end
1466
1467 c----------------------------------------------------------------------
1468       function ranptcut(xcut)
1469 c----------------------------------------------------------------------
1470 c .........exp(-x**2)
1471  12   x=sqrt(-(log(rangen()))/(3.1415927/4.)/xcut) !gauss
1472 C 12   x=sqrt(-(log(rangen()))/(3.1415927/4.)/xcut) !gauss
1473
1474
1475 c      if(xcut.gt.0.)then
1476 c        if(rangen().lt.x/xcut)goto 12
1477 c      endif
1478
1479       ranptcut=x
1480
1481       return
1482
1483 c .........exp(-x)
1484 c  12  xmx=50
1485 c      r=2.
1486 c      do while (r.gt.1.)
1487 c  11    x=sqrt(exp(rangen()*log(1+xmx**2))-1)
1488 c        if(x.eq.0.)goto11
1489 c        r=rangen()  /  ( exp(-x)*(1+x**2) )
1490 c      enddo
1491 c      x=x/2.
1492
1493       end
1494
1495 cc----------------------------------------------------------------------
1496 c      function ranpticut(xcut)
1497 cc----------------------------------------------------------------------
1498 c
1499 cc .........exp(-x)
1500 c
1501 c  12  xmx=50
1502 c      r=2.
1503 c      do while (r.gt.1.)
1504 c  11    x=sqrt(exp(rangen()*log(1+xmx**2))-1)
1505 c        if(x.eq.0.)goto11
1506 c        r=rangen()  /  ( exp(-x)*(1+x**2) )
1507 c      enddo
1508 c      x=x/2.
1509 c      if(rangen().gt.xcut/x)goto 12
1510 c
1511 c      ranpticut=x
1512 c
1513 c      end
1514
1515 c----------------------------------------------------------------------
1516       function ranpt()
1517 c----------------------------------------------------------------------
1518
1519 c .........exp(-x)
1520       xmx=50
1521       r=2.
1522       do while (r.gt.1.)
1523   11    x=sqrt(exp(rangen()*log(1+xmx**2))-1)
1524         if(x.eq.0.)goto11
1525         r=rangen()  /  ( exp(-x)*(1+x**2) )
1526       enddo
1527       ranpte=x/2.
1528
1529 c     -------------
1530       ranpt=ranpte
1531       return
1532 c     -------------
1533 c .........exp(-x**2)
1534       ranptg=sqrt(-log(rangen())/(3.1415927/4.)) !gauss
1535
1536 c .........exp(-sqrt(x))
1537       xmx=500
1538       r=2.
1539       do while (r.gt.1.)
1540         x=sqrt(exp(rangen()*log(1+xmx**2))-1)
1541         r=rangen()  /  ( exp(-sqrt(x))*(1+x**2)/5. )
1542       enddo
1543       ranpts=x/20.
1544
1545
1546       end
1547 c----------------------------------------------------------------------
1548       function ranptk()
1549 c----------------------------------------------------------------------
1550
1551 c .........exp(-x)
1552       xmx=50
1553       r=2.
1554       do while (r.gt.1.)
1555   11    x=sqrt(exp(rangen()*log(1+xmx**2))-1)
1556         if(x.eq.0.)goto11
1557         r=rangen()  /  ( exp(-x)*(1+x**2) )
1558       enddo
1559       ranpte=x/2.
1560
1561 c     -------------
1562       ranptk=ranpte
1563       return
1564 c     -------------
1565
1566 c .........exp(-x**2)
1567       ranptg=sqrt(-log(rangen())/(3.1415927/4.)) !gauss
1568
1569 c .........exp(-sqrt(x))
1570       xmx=500
1571       r=2.
1572       do while (r.gt.1.)
1573         x=sqrt(exp(rangen()*log(1+xmx**2))-1)
1574         r=rangen()  /  ( exp(-sqrt(x))*(1+x**2)/5. )
1575       enddo
1576       ranpts=x/20.
1577
1578       end
1579
1580 c----------------------------------------------------------------------
1581       function ranptd()
1582 c----------------------------------------------------------------------
1583
1584 c .........exp(-x**2)
1585       ranptg=sqrt(-log(rangen())/(3.1415927/4.)) !gauss
1586
1587
1588 c     -------------
1589       ranptd=ranptg
1590       return
1591 c     -------------
1592
1593 c .........exp(-x)
1594       xmx=50
1595       r=2.
1596       do while (r.gt.1.)
1597   11    x=sqrt(exp(rangen()*log(1+xmx**2))-1)
1598         if(x.eq.0.)goto11
1599         r=rangen()  /  ( exp(-x)*(1+x**2) )
1600       enddo
1601       ranpte=x/2.
1602
1603 c .........exp(-sqrt(x))
1604       xmx=500
1605       r=2.
1606       do while (r.gt.1.)
1607         x=sqrt(exp(rangen()*log(1+xmx**2))-1)
1608         r=rangen()  /  ( exp(-sqrt(x))*(1+x**2)/5. )
1609       enddo
1610       ranpts=x/20.
1611
1612
1613
1614
1615
1616       end
1617 c----------------------------------------------------------------------
1618       subroutine gakdel(ibr)
1619 c----------------------------------------------------------------------
1620       parameter (mxnbr=500,mxnba=5000)
1621       common /gag/nob,pst(4,0:mxnba),ipst(0:mxnba)
1622      $     ,nbr,ijb(2,0:mxnbr),xyb(2,0:mxnbr)
1623      &     ,ptb(4,0:mxnbr),iflb(0:mxnbr),ip(0:mxnbr),co,x,y
1624       dimension co(12)
1625
1626       do i=ibr,nbr+1
1627          do j=1,2
1628             ijb(j,i)=ijb(j,i+1)
1629             xyb(j,i)=xyb(j,i+1)
1630          enddo
1631          do k=1,4
1632             ptb(k,i)=ptb(k,i+1)
1633          enddo
1634          iflb(i)=iflb(i+1)
1635          ip(i)=ip(i+1)
1636       enddo
1637       ip(ibr)=0
1638       nbr=nbr-1
1639       end
1640
1641 c----------------------------------------------------------------------
1642       subroutine gaksco(ibr1,ibr,ibr2,ib,jb,x1,x2,y1,y2)
1643 c----------------------------------------------------------------------
1644       include 'epos.inc'
1645
1646       parameter (mxnbr=500,mxnba=5000)
1647       common /gag/nob,pst(4,0:mxnba),ipst(0:mxnba)
1648      $     ,nbr,ijb(2,0:mxnbr),xyb(2,0:mxnbr)
1649      &     ,ptb(4,0:mxnbr),iflb(0:mxnbr),ip(0:mxnbr),co,x,y
1650       dimension co(12)
1651       logical weiter
1652
1653       pf(i,j)=pst(4,i)*pst(4,j)-pst(3,i)*pst(3,j)
1654      &     -pst(2,i)*pst(2,j)-pst(1,i)*pst(1,j)
1655       mmod(i)=mod(mod(i,nob)+nob,nob)
1656
1657       call utpri('gaksco',ish,ishini,8)
1658
1659       if(ish.ge.8)then
1660          write(ifch,*) 'zwischen brk:',ibr1,ibr,ibr2,'(',nob,')',nbr
1661          write(ifch,*) 'region:',ib,jb
1662       endif
1663       if(ib.eq.ijb(1,ibr1))then
1664          x2=xyb(1,ibr1)
1665       else
1666          x2=1.
1667       endif
1668       if(ib.eq.ijb(1,ibr2))then
1669          x1=xyb(1,ibr2)
1670       else
1671          x1=0.
1672       endif
1673       if(jb.eq.ijb(2,ibr1))then
1674          y1=xyb(2,ibr1)
1675       else
1676          y1=0.
1677       endif
1678       if(jb.eq.ijb(2,ibr2))then
1679          y2=xyb(2,ibr2)
1680       else
1681          y2=1.
1682       endif
1683
1684 c.....left side...................................................
1685       n=nob+1
1686       if(ish.ge.8)write(ifch,*)'x1,x2',x1,x2
1687       do i=1,4
1688 cc        pst(i,n)=(x2-x1)*pst(i,ib)+ptb(i,ibr)-ptb(i,ibr1)
1689         pst(i,n)=     x2*pst(i,ib)+ptb(i,ibr)-ptb(i,ibr1)-y1*pst(i,jb)
1690       enddo
1691       if(ish.ge.8)write(ifch,*) 'add a  1 ii',1.,ib
1692       ii=mmod(ib-1)
1693       weiter=.true.
1694       if(ib.eq.ijb(1,ibr1))weiter=.false.
1695       do while(ii.ne.jb.and.weiter) !linker Rand??
1696          f1=1.
1697          if(ii.eq.ijb(1,ibr1))f1=xyb(1,ibr1)
1698          do i=1,4
1699             pst(i,n)=pst(i,n)+f1*pst(i,ii)
1700          enddo
1701          if(ish.ge.8)write(ifch,*) 'add a f1 ii',f1,ii
1702          if(ii.eq.ijb(1,ibr1))weiter=.false.
1703          ii=mmod(ii-1)
1704       enddo
1705       jj=mmod(jb+1)
1706       weiter=.not.weiter
1707       if(jb.eq.ijb(2,ibr1))weiter=.false.
1708       do while(weiter)
1709          f1=1.
1710          if(jj.eq.ijb(2,ibr1))f1=1.-xyb(2,ibr1)
1711          do i=1,4
1712             pst(i,n)=pst(i,n)+f1*pst(i,jj)
1713          enddo
1714          if(ish.ge.8)write(ifch,*) 'add b f1 ii',f1,jj
1715          if(jj.eq.ijb(2,ibr1))weiter=.false.
1716          jj=mmod(jj+1)
1717       enddo
1718       do i=1,4
1719 cc        pst(i,n+1)=(x2-x1)*pst(i,ib)
1720         pst(i,n+1)=        pst(i,ib)
1721 cc        pst(i,n+2)=(y2-y1)*pst(i,jb)
1722         pst(i,n+2)=        pst(i,jb)
1723       enddo
1724       co(1)= pf(n,n)
1725       co(2)=-2.*pf(n,n+1)
1726       co(3)= 2.*pf(n,n+2)
1727       co(4)=-2.*pf(n+1,n+2)
1728       if(ish.ge.8) then
1729         do i=n,n+2
1730           write (ifch,'(4g12.5)') (pst(j,i),j=1,4)
1731         enddo
1732       endif
1733       if(ish.ge.8)write(ifch,*) 'co left:',co(1),co(2),co(3),co(4)
1734
1735 c.....right side...................................................
1736       n=nob+4
1737       if(ish.ge.8)write(ifch,*)'y1,y2',y1,y2
1738       do i=1,4
1739 cc         pst(i,n)=(y2-y1)*pst(i,jb)-ptb(i,ibr)+ptb(i,ibr2)
1740          pst(i,n)=    y2*pst(i,jb)-ptb(i,ibr)+ptb(i,ibr2)-x1*pst(i,ib)
1741       enddo
1742       if(ish.ge.8)write(ifch,*) 'add a  1 jj',1.,jb
1743       ii=mmod(ib+1)
1744       weiter=.true.
1745       if(ib.eq.ijb(1,ibr2))weiter=.false.
1746       do while(ii.ne.jb.and.weiter)
1747          f1=1.
1748          if(ii.eq.ijb(1,ibr2))f1=1.-xyb(1,ibr2)
1749          do i=1,4
1750             pst(i,n)=pst(i,n)+f1*pst(i,ii)
1751          enddo
1752          if(ish.ge.8)write(ifch,*) 'add a f1 ii',f1,ii
1753          if(ii.eq.ijb(1,ibr2))weiter=.false.
1754          ii=mmod(ii+1)
1755       enddo
1756       jj=mmod(jb-1)
1757       weiter=.not.weiter
1758       if(jb.eq.ijb(2,ibr2))weiter=.false.
1759       do while(weiter)
1760          f1=1.
1761          if(jj.eq.ijb(2,ibr2))f1=xyb(2,ibr2)
1762          do i=1,4
1763             pst(i,n)=pst(i,n)+f1*pst(i,jj)
1764          enddo
1765          if(ish.ge.8)write(ifch,*) 'add b f1 ii',f1,jj
1766          if(jj.eq.ijb(2,ibr2))weiter=.false.
1767          jj=mmod(jj-1)
1768       enddo
1769       do i=1,4
1770 cc         pst(i,n+1)=(x2-x1)*pst(i,ib)
1771          pst(i,n+1)=        pst(i,ib)
1772 cc         pst(i,n+2)=(y2-y1)*pst(i,jb)
1773          pst(i,n+2)=         pst(i,jb)
1774       enddo
1775       co(5)=pf(n,n)
1776       co(6)= 2.*pf(n,n+1)
1777       co(7)=-2.*pf(n,n+2)
1778       co(8)=-2.*pf(n+1,n+2)
1779       if(ish.ge.8) then
1780         do i=n,n+2
1781           write (ifch,'(4g12.5)') (pst(j,i),j=1,4)
1782         enddo
1783       endif
1784       if(ish.ge.8)write(ifch,*) 'co right:',co(5),co(6),co(7),co(8)
1785
1786 c.....lambda (absolute past).....................................
1787       n=nob+7
1788       do i=1,4
1789 cc         pst(i,n)= x1*pst(i,ib)+y1*pst(i,jb)
1790          pst(i,n)= 0.
1791       enddo
1792       ii=mmod(ib+1)
1793       do while (mmod(ii+jb).ne.0)
1794          if(ish.ge.8)write(ifch,*) 'add lambda',ii
1795          do i=1,4
1796             pst(i,n)=pst(i,n)+pst(i,ii)
1797          enddo
1798          ii=mmod(ii+1)
1799       enddo
1800       do i=1,4
1801 cc         pst(i,n+1)=(x2-x1)*pst(i,ib)
1802 cc         pst(i,n+2)=(y2-y1)*pst(i,jb)
1803          pst(i,n+1)= pst(i,ib)
1804          pst(i,n+2)= pst(i,jb)
1805       enddo
1806       co(9)=     pf(n,n)
1807       co(10)= 2.*pf(n,n+1)
1808       co(11)= 2.*pf(n,n+2)
1809       co(12)= 2.*pf(n+1,n+2)
1810       if(ish.ge.8)write(ifch,*)'co lambda:',co(9),co(10),co(11),co(12)
1811       call utprix('gaksco',ish,ishini,8)
1812       end
1813
1814 c---------------------------------------------------------------------
1815       subroutine gakanp(ibr1,ibr,ibrr2,aml2,amr2,ala2,iok)
1816 c---------------------------------------------------------------------
1817 c   mass adjustment of fragments
1818 c        ibr1-ibr   ibr-ibrr2
1819 c   where ibr1,ibr,ibrr2 are break point indices
1820 c   aml2,amr2 are the reqired squared masses (if zero -> any mass)
1821 c   iok=0 (ok) or ok=1 (error)
1822 c---------------------------------------------------------------------
1823       include 'epos.inc'
1824       parameter (mxnbr=500,mxnba=5000,mxnin=2000)
1825       common /gag/nob,pst(4,0:mxnba),ipst(0:mxnba)
1826      $     ,nbr,ijb(2,0:mxnbr),xyb(2,0:mxnbr)
1827      &     ,ptb(4,0:mxnbr),iflb(0:mxnbr),ip(0:mxnbr),co,x,y
1828       double precision ax,ay,az,ae,am,bx,by,bz,be,A,B,C
1829       dimension co(12),am2(0:2),nin(2,0:mxnin)
1830       logical weiter
1831 c      pf(i,j)=pst(4,i)*pst(4,j)-pst(3,i)*pst(3,j)
1832 c     &     -pst(2,i)*pst(2,j)-pst(1,i)*pst(1,j)
1833       mmod(i)=mod(mod(i,nob)+nob,nob)
1834       call utpri('gakanp',ish,ishini,6)
1835
1836       ibr2=ibrr2
1837       if(ish.ge.6)write(ifch,*) ibr1,ibr,ibr2,aml2,amr2,ala2,iok
1838       iok=0
1839       ib=ijb(1,ibr)
1840       jb=ijb(2,ibr)
1841       ni=0
1842  10   do i=1,ni
1843         if((nin(1,i).eq.ib.and.nin(2,i).eq.jb)
1844      $       .or.(ipst(ib).eq.ipst(jb)))then
1845           iok=1
1846           if(ish.ge.4)then
1847             write(ifch,*) 'error ... endless loop'
1848             if(ib.eq.ipst(jb))  write(ifch,*) ' in zero mass region'
1849           endif
1850           goto 9999
1851         endif
1852       enddo
1853       ni=ni+1
1854       if(ni.gt.mxnin)stop'gakanp: increase parameter mxnin  '
1855       nin(1,ni)=ib
1856       nin(2,ni)=jb
1857       if(ish.ge.6)write(ifch,*)
1858       if(ish.ge.6)write(ifch,*) 'ib,jb=',ib,jb
1859       if(ni.ge.2)then
1860          if(ish.ge.6)write(ifch,*)'rotate pt to new band'
1861          if(ish.ge.6)write(ifch,*)'from',nin(1,ni-1),nin(2,ni-1)
1862          if(ish.ge.6)write(ifch,*)'  to',ib,jb
1863          if(ish.ge.6)write(ifch,'(4f12.8)') (ptb(i,ibr),i=1,4)
1864          ax=pst(1,nin(1,ni-1))+pst(1,nin(2,ni-1))
1865          ay=pst(2,nin(1,ni-1))+pst(2,nin(2,ni-1))
1866          az=pst(3,nin(1,ni-1))+pst(3,nin(2,ni-1))
1867          ae=pst(4,nin(1,ni-1))+pst(4,nin(2,ni-1))
1868          am=sngl(sqrt(max(1d-8,dble(ae)**2-dble(ax)**2
1869      $        -dble(ay)**2-dble(az)**2))) !???????????????????????
1870          bx=pst(1,nin(1,ni-1))
1871          by=pst(2,nin(1,ni-1))
1872          bz=pst(3,nin(1,ni-1))
1873          be=pst(4,nin(1,ni-1))
1874          if(ish.ge.6)write (ifch,*) 'bx,by,bz,be',bx,by,bz,be
1875          call utlob2(1,ax,ay,az,ae,am,bx,by,bz,be,60)
1876          if(ish.ge.6)write (ifch,*) 'bx,by,bz,be',bx,by,bz,be
1877          call utlob4( 1,ax,ay,az,ae,am
1878      &        ,ptb(1,ibr),ptb(2,ibr),ptb(3,ibr),ptb(4,ibr))
1879          if(ish.ge.6)write(ifch,'(4f12.8)') (ptb(i,ibr),i=1,4)
1880          call utrot4( 1,bx,by,bz,ptb(1,ibr),ptb(2,ibr),ptb(3,ibr))
1881          if(ish.ge.6)write(ifch,'(4f12.8)') (ptb(i,ibr),i=1,4)
1882          ax=pst(1,ib)+pst(1,jb)
1883          ay=pst(2,ib)+pst(2,jb)
1884          az=pst(3,ib)+pst(3,jb)
1885          ae=pst(4,ib)+pst(4,jb)
1886          am=sngl(sqrt(max(1d-8,dble(ae)**2-dble(ax)**2
1887      $        -dble(ay)**2-dble(az)**2))) !???????????????????????
1888          if(am.le.1.1e-4)then
1889            if(ish.ge.5)write(ifch,*)'error ... am<1.1e-4'
1890            iok=1
1891            goto 9999
1892          endif
1893          bx=pst(1,ib)
1894          by=pst(2,ib)
1895          bz=pst(3,ib)
1896          be=pst(4,ib)
1897          if(ish.ge.6)write (ifch,*) 'bx,by,bz,be',bx,by,bz,be
1898          if(ish.ge.6)write (ifch,*) 'ax,ay,az,ae',ax,ay,az,ae,am
1899          call utlob2(1,ax,ay,az,ae,am,bx,by,bz,be,60)
1900          if(ish.ge.6)write (ifch,*) 'bx,by,bz,be',bx,by,bz,be
1901          call utrot4(-1,bx,by,bz,ptb(1,ibr),ptb(2,ibr),ptb(3,ibr))
1902          if(ish.ge.6)write(ifch,'(4f12.8)') (ptb(i,ibr),i=1,4)
1903          call utlob4(-1,ax,ay,az,ae,am
1904      &        ,ptb(1,ibr),ptb(2,ibr),ptb(3,ibr),ptb(4,ibr))
1905          if(ish.ge.6)write(ifch,'(4f12.8)') (ptb(i,ibr),i=1,4)
1906       endif
1907       call gaksco(ibr1,ibr,ibr2,ib,jb,x1,x2,y1,y2)
1908 c      if(ni.eq.1)print *,'------------------------------------',2
1909 cc      x=(xyb(1,ibr)-x1)/(x2-x1)
1910 cc      y=(xyb(2,ibr)-y1)/(y2-y1)
1911       x=xyb(1,ibr)
1912       y=xyb(2,ibr)
1913       am2(0)=aml2
1914       am2(1)=amr2
1915       am2(2)=ala2
1916       if(ish.ge.6)write(ifch,*) ibr1,ibr,ibr2,aml2,amr2,ala2,iok
1917       if(amr2.le.0.)then
1918          l1=2
1919          l2=0
1920       elseif(aml2.le.0.)then
1921          l1=1
1922          l2=2
1923       elseif(ala2.le.0.)then
1924          l1=1
1925          l2=0
1926       else
1927          stop' not like this , please...'
1928       endif
1929       if(ish.ge.6.and.amr2.le.0)write(ifch,*) 'adjust: 1',l1,l2
1930       if(ish.ge.6.and.aml2.le.0)write(ifch,*) 'adjust: 2' ,l1,l2
1931       if(ish.ge.6.and.ala2.le.0)write(ifch,*) 'adjust: 3',l1,l2
1932       i=4*l1
1933       j=4*l2
1934       A = dble(co(i+4))*dble(co(j+3)) - dble(co(j+4))*dble(co(i+3))
1935       B = dble(co(i+4))*dble(co(j+1)) - dble(co(i+3))*dble(co(j+2))
1936      &  + dble(co(i+2))*dble(co(j+3)) - dble(co(i+1))*dble(co(j+4))
1937      &  - dble(am2(l2))*dble(co(i+4)) + dble(am2(l1))*dble(co(j+4))
1938       C = dble(co(i+2))*dble(co(j+1)) - dble(co(i+1))*dble(co(j+2))
1939      &  + dble(am2(l1))*dble(co(j+2)) - dble(am2(l2))*dble(co(i+2))
1940       if (ish.ge.7) then
1941          write(ifch,*) 'ABC,q ',A,B,C,B**2-4.*A*C
1942          if(abs(A).gt.0.d0)then
1943             write(ifch,*) sqrt(max(0d0,B**2-4d0*A*C))/2d0/A-B/A/2d0
1944             write(ifch,*) -sqrt(max(0d0,B**2-4d0*A*C))/2d0/A-B/A/2d0
1945          endif
1946       endif
1947       x=0.
1948       y=0.
1949       xx=0.
1950       yy=0.
1951       if(abs(A).gt.1.d-20.and.B*B-4.*A*C.ge.0.d0)then
1952         y=sngl(sqrt(max(0.d0,B**2-4.*A*C))/2.d0/A-B/A/2.d0)
1953         if(abs(co(i+2)+y*co(i+4)).gt.0.)
1954      &       x=(am2(l1)-co(i+1)-y*co(i+3))/(co(i+2)+y*co(i+4))
1955       elseif(abs(A).le.1.d-20.and.abs(B).gt.0.d0)then
1956         y=-sngl(C/B)
1957         if(abs(co(i+2)+y*co(i+4)).gt.0.)
1958      &       x=(am2(l1)-co(i+1)-y*co(i+3))/(co(i+2)+y*co(i+4))
1959       else
1960         if(ish.ge.5)write(ifch,*)'error ... no solution of quadr equ'
1961         iok=1
1962         goto 9999
1963       endif
1964       if(abs(A).gt.1.d-20.and.B**2-4.*A*C.ge.0.d0)then
1965         yy=sngl(-sqrt(max(0.d0,B**2-4.*A*C))/2.d0/A-B/A/2.d0)
1966         if(abs(co(i+2)+yy*co(i+4)).gt.0.)
1967      &       xx=(am2(l1)-co(i+1)-yy*co(i+3))/(co(i+2)+yy*co(i+4))
1968       elseif(abs(A).le.1.d-20.and.abs(B).gt.0.d0)then
1969         yy=-sngl(C/B)
1970         if(abs(co(i+2)+yy*co(i+4)).gt.0.)
1971      &       xx=(am2(l1)-co(i+1)-yy*co(i+3))/(co(i+2)+yy*co(i+4))
1972       else
1973          if(ish.ge.5)write(ifch,*)'error ... no solution (2) '
1974          iok=1
1975          goto 9999
1976       endif
1977       if(ish.ge.6)then
1978          write(ifch,*) x ,y ,(co(i+2)+ y*co(i+4)),' OK '
1979          write(ifch,*) xx,yy,(co(i+2)+yy*co(i+4)),' OK '
1980       endif
1981       weiter=.true.
1982  50   if(x.gt.x1.and.x.lt.x2.and.y.gt.y1.and.y.lt.y2)then
1983 cc         xyb(1,ibr)=x1+(x2-x1)*x
1984 cc         xyb(2,ibr)=y1+(y2-y1)*y
1985          xyb(1,ibr)=x
1986          xyb(2,ibr)=y
1987          ijb(1,ibr)=ib
1988          ijb(2,ibr)=jb
1989          e1=pst(4,nob+1)-x*pst(4,nob+2)+y*pst(4,nob+3)
1990          e2=pst(4,nob+4)+x*pst(4,nob+5)-y*pst(4,nob+6)
1991          if( e1.lt.0. .or. e2.lt.0. ) then
1992            if(ish.ge.5)write(ifch,*)'error ... e1<0 or e2<0'
1993            iok=1
1994            goto 9999
1995          endif
1996          !amal2=co(1)+co(2)*x+co(3)*y+co(4)*x*y
1997          !amar2=co(5)+co(6)*x+co(7)*y+co(8)*x*y
1998          if(ish.ge.6)then
1999            amal2=co(1)+co(2)*x+co(3)*y+co(4)*x*y
2000            amar2=co(5)+co(6)*x+co(7)*y+co(8)*x*y
2001            write(ifch,*) 'brkshift:',xyb(1,ibr),xyb(2,ibr),ib,jb
2002            write (ifch,*)'E:',e1
2003            write (ifch,*)'E:',e2
2004            write(ifch,'(2(a6,1g12.6))') 'aml:'
2005      &          ,sqrt(max(0.,amal2)),'amr:',sqrt(max(0.,amar2))
2006          endif
2007          i=ibr1+1
2008          do while(i.le.ibr-1)
2009             if((mmod(ijb(1,i)+nob/2) .lt. mmod(ijb(1,ibr)+nob/2)
2010      &           .or.(mmod(ijb(1,i)+nob/2) .eq. mmod(ijb(1,ibr)+nob/2)
2011      &           .and.(xyb(1,i).gt.xyb(1,ibr))))
2012      &           .and.
2013      &           (ijb(2,i) .gt. ijb(2,ibr)
2014      &           .or.(ijb(2,i) .eq. ijb(2,ibr)
2015      &           .and.xyb(2,i).lt.xyb(2,ibr)))) goto 150
2016             if(ish.ge.6) then
2017                write(ifch,*) 'away:'
2018      &          ,i,xyb(1,i),xyb(2,i),ijb(1,i),ijb(2,i)
2019             endif
2020             call gakdel(i)
2021             i=i-1
2022             ibr=ibr-1
2023             ibr2=ibr2-1
2024  150        i=i+1
2025          enddo
2026          i=ibr+1
2027          do while (i.le.ibr2-1)
2028             if((mmod(ijb(1,i)+nob/2) .gt. mmod(ijb(1,ibr)+nob/2)
2029      &           .or.(mmod(ijb(1,i)+nob/2) .eq. mmod(ijb(1,ibr)+nob/2)
2030      &           .and.(xyb(1,i).lt.xyb(1,ibr))))
2031      &           .and.
2032      &           (ijb(2,i) .lt. ijb(2,ibr)
2033      &           .or.(ijb(2,i) .eq. ijb(2,ibr)
2034      &           .and.xyb(2,i).gt.xyb(2,ibr)))) goto 160
2035             if(ish.ge.6) then
2036                write(ifch,*) 'away:'
2037      &          ,i,xyb(1,i),xyb(2,i),ijb(1,i),ijb(2,i)
2038             endif
2039             call gakdel(i)
2040             ibr2=ibr2-1
2041             i=i-1
2042  160        i=i+1
2043          enddo
2044          goto 9999
2045       else
2046         if(x.gt.x2
2047      &    .and.ib.ne.ijb(1,ibr1) !brk-begrenzung
2048      &    .and.mmod(ib-1).ne.jb !linker oder rechter Rand
2049      &    .and.mmod(ib-1+jb).ne.0)then !oben oder unten
2050           ib=mmod(ib-1)
2051           goto 10
2052         endif
2053         if(x.lt.x1
2054      &    .and.ib.ne.ijb(1,ibr2) !brk-begrenzung
2055      &    .and.mmod(ib+1).ne.jb !linker oder rechter Rand
2056      &    .and.mmod(ib+1+jb).ne.0)then !oben oder unten
2057           ib=mmod(ib+1)
2058           goto 10
2059         endif
2060         if(y.gt.y2
2061      &    .and.jb.ne.ijb(2,ibr2) !brk-begrenzung
2062      &    .and.mmod(jb-1).ne.ib !linker oder rechter Rand
2063      &    .and.mmod(jb-1+ib).ne.0)then !oben oder unten
2064           jb=mmod(jb-1)
2065           goto 10
2066         endif
2067         if(y.lt.y1
2068      &    .and.jb.ne.ijb(2,ibr1) !brk-begrenzung
2069      &    .and.mmod(jb+1).ne.ib !linker oder rechter Rand
2070      &    .and.mmod(jb+1+ib).ne.0)then !oben oder unten
2071           jb=mmod(jb+1)
2072           goto 10
2073         endif
2074         if(weiter)then
2075           weiter=.false.
2076           x=xx
2077           y=yy
2078           goto 50
2079         endif
2080         if(ish.ge.5)write(ifch,*)'error ... x,y not in allowed range'
2081         iok=1
2082       endif
2083  9999 if(amr2.eq.0.) ibrr2=ibr2
2084       call utprix('gakanp',ish,ishini,6)
2085       end
2086
2087 cc----------------------------------------------------------------------
2088 c      subroutine gakstr(ifl)
2089 cc----------------------------------------------------------------------
2090 cc
2091 cc     calculates string-fragments by taking off pt of breakup
2092 cc     do with ifl=1   undo with ifl=-1
2093 cc
2094 cc----------------------------------------------------------------------
2095 c      include 'epos.inc'
2096 c      common /cpptr/ pptr(4,mxptl),ystr(mxptl)
2097 c
2098 c      do i=1,nptl
2099 c        if(istptl(i).eq.29)then
2100 c          nk1=ifrptl(1,i)
2101 c          nk2=ifrptl(2,i)
2102 c          do j=nk1,nk2
2103 c            if ((istptl(j).eq.0.or.istptl(j-1).eq.0).and.j.ne.nk1) then
2104 c              do k=1,4
2105 c                pptl(k,j)=pptl(k,j)+pptr(k,j-1)*ifl
2106 c              enddo
2107 c              !write(ifch,*)"left side back  to ",j,(pptr(k,j-1),k=1,4)
2108 c            endif
2109 c            if ((istptl(j).eq.0.or.istptl(j+1).eq.0).and.j.ne.nk2) then
2110 c              do k=1,4
2111 c                pptl(k,j)=pptl(k,j)-pptr(k,j)*ifl
2112 c              enddo
2113 c              !write(ifch,*)"right side back to ",j,(-pptr(k,j),k=1,4)
2114 c            endif
2115 c            if(ifl.eq.-1.and.istptl(j).eq.0)then
2116 c           e=pptl(1,j)**2+pptl(2,j)**2+pptl(3,j)**2
2117 c     &           +pptl(5,j)**2
2118 c              e=sqrt(e)
2119 c           !dif=abs(e-pptl(4,j))
2120 c           !if(dif.gt.0.01.and.dif/e.gt.0.1)print*,j,e,pptl(4,j)
2121 c              pptl(4,j)=e
2122 c            endif
2123 c          enddo
2124 c        endif
2125 c        if(istptl(i).eq.0)then
2126 c          if ( ifl.eq.1 ) then
2127 c            ystr(i)=sign(1.,pptl(3,i))*alog((pptl(4,i)+abs(pptl(3,i)))
2128 c     *           /sqrt(pptl(5,i)**2+pptl(1,i)**2+pptl(2,i)**2) )
2129 c          endif
2130 c        endif
2131 c      enddo
2132 c      end
2133
2134 c----------------------------------------------------------------------
2135       subroutine gakli2(nn1,nn2)
2136 c----------------------------------------------------------------------
2137
2138       include 'epos.inc'
2139       double precision pgampr,rgampr
2140       common/cgampr/pgampr(5),rgampr(4)
2141 c      double precision db1,db2,db3,db4,db5
2142       character label*8,idlabl*8
2143 c      db1=0d0
2144 c      db2=0d0
2145 c      db3=rgampr(4)
2146 c      db4=sqrt(rgampr(1)**2+rgampr(2)**2+rgampr(3)**2)
2147 c      db5=sqrt( db4**2-rgampr(4)**2)
2148       n1=nn1
2149       n2=nn2
2150       if (n1.eq.0) n1=1
2151       if (n2.eq.0) n2=nptl
2152       write (ifch,'(1a4,5a12,4a4,a10,2a4)')
2153      &'no.','px','py','pz','E','m','ior','jor','if1','if2'
2154      &,'id','ist','ity'
2155       do i=n1,n2
2156          if (idptl(i).lt.10000)then
2157             label='        '
2158             label=idlabl(idptl(i))
2159          else
2160             label='        '
2161          endif
2162          chrg=0.
2163          if(iabs(idptl(i)).le.9999)call idchrg(idptl(i),chrg)
2164          write (ifch,125) i,(pptl(j,i),j=1,5),iorptl(i),jorptl(i)
2165      &        ,ifrptl(1,i),ifrptl(2,i),idptl(i)
2166      $        ,chrg             !charge
2167      &        ,istptl(i),ityptl(i),label
2168       enddo
2169  125  format (1i4,5g18.10,4i6,1i10
2170      $     ,1f5.2               !charge
2171      $     ,2i4,'  ',A8
2172 c     $     ,7g12.4,i5
2173      $     )
2174  126  format (A4,7g12.4,i5)
2175       return
2176       end
2177
2178 c----------------------------------------------------------------------
2179 c      subroutine gakli4
2180 cc----------------------------------------------------------------------
2181 c
2182 c      include 'epos.inc'
2183 c      parameter (mxnbr=500,mxnba=5000)
2184 c      common /gag/nob,pst(4,0:mxnba),ipst(0:mxnba)
2185 c     $     ,nbr,ijb(2,0:mxnbr),xyb(2,0:mxnbr)
2186 c     &     ,ptb(4,0:mxnbr),iflb(0:mxnbr),ip(0:mxnbr),co,x,y
2187 c      dimension co(12)
2188 c      do i=0,nob-1
2189 c        write (ifch,10) i,(pst(j,i),j=1,4)
2190 c      enddo
2191 c 10   format(1i4,5g18.10)
2192 c      return
2193 c      end
2194 c
2195 cc----------------------------------------------------------------------
2196 c      subroutine gakli3
2197 cc----------------------------------------------------------------------
2198 c
2199 c      include 'epos.inc'
2200 c      parameter (mxnbr=500,mxnba=5000)
2201 c      common /gag/nob,pst(4,0:mxnba),ipst(0:mxnba)
2202 c     $     ,nbr,ijb(2,0:mxnbr),xyb(2,0:mxnbr)
2203 c     &     ,ptb(4,0:mxnbr),iflb(0:mxnbr),ip(0:mxnbr),co,x,y
2204 c      dimension co(12),p1(5),p2(5)
2205 c
2206 c      write(ifch,*) 'particle list of string decay'
2207 c      do i=1,5
2208 c        p1(i)=0.
2209 c      enddo
2210 c      do i=1,nbr+1
2211 c        if(i.lt.nbr+1)then
2212 c          call gaksco(i-1,i,i+1,ijb(1,i),ijb(2,i),x1,x2,y1,y2)
2213 c          if(x2.gt.x1)then
2214 ccc            x=(xyb(1,i)-x1)/(x2-x1)
2215 c            x=xyb(1,i)
2216 c          else
2217 c            x=0.
2218 c          endif
2219 c          if(y2.gt.y1)then
2220 ccc            y=(xyb(2,i)-y1)/(y2-y1)
2221 c            y=xyb(2,i)
2222 c          else
2223 c            y=0.
2224 c          endif
2225 c          aml2=co(1)+x*co(2)+y*co(3)+x*y*co(4)
2226 c          amr2=co(5)+x*co(6)+y*co(7)+x*y*co(8)
2227 c          aml=sign(sqrt(abs(aml2)),aml2)
2228 c          amr=sign(sqrt(abs(amr2)),amr2)
2229 c          do j=1,4
2230 c            p2(j)=pst(j,nob+1)-x*pst(j,nob+2)+y*pst(j,nob+3)
2231 c            p1(j)=p1(j)+p2(j)
2232 c          enddo
2233 c          p2(5)=aml
2234 c        else
2235 c          do j=1,4
2236 c            p2(j)=pst(j,nob+4)+x*pst(j,nob+5)-y*pst(j,nob+6)
2237 c            p1(j)=p1(j)+p2(j)
2238 c          enddo
2239 c          p2(5)=amr
2240 c        endif
2241 c        write(ifch,'(2i4,i6,a,i5,i10,5g14.6)') i-1,i
2242 c     &    ,-iflb(i-1),'==',iflb(i),ip(i)
2243 c     &    ,(p2(j),j=1,5)
2244 c      enddo
2245 c      am2=p1(4)**2-p1(3)**2-p1(2)**2-p1(1)**2
2246 c      p1(5)=sign(sqrt(abs(am2)),am2)
2247 c      write(ifch,'(12x,a60)')
2248 c     &  '------------------------------------------------------------'
2249 c      write(ifch,'(14x,5g14.6)') (p1(j),j=1,5)
2250 c      write(ifch,*)
2251 c
2252 c      end
2253 c
2254
2255 c---------------------------------------------------------------------
2256       subroutine idress(id,am,idr,iadj)
2257 c---------------------------------------------------------------------
2258       include 'epos.inc'
2259       call idres(id,am,idr,iadj,1)
2260       if(idr.eq.0)then
2261         return
2262       endif
2263       ids=max(mod(iabs(id)/100,10),mod(iabs(id)/10,10))
2264       if(iabs(idr).le.999) then
2265 c        write (*,*) '  ',id,idr
2266         if(ids.le.3)return
2267         if(ids.le.2)then
2268           idr=sign(iabs(id)+int(rangen()+0.5),id)
2269         elseif(ids.eq.3)then
2270           idr=sign(iabs(id)+int(rangen()+0.6),id)
2271         else
2272           idr=sign(iabs(id)+int(rangen()+0.75),id)
2273         endif
2274 c        write (*,*) '->',id,idr
2275         call idmass(idr,am)
2276       elseif(iabs(idr).le.9999)then
2277         if(ids.le.3)return
2278         if(mod(iabs(idr),10).gt.1)then
2279           if(iabs(id).ne.1111.and.iabs(id).ne.2221.and.iabs(id).ne.3331)
2280      $         then
2281             idr=sign(iabs(id)+1,id)
2282             call idmass(idr,am)
2283           else
2284             idr=id
2285             call idmass(idr,am)
2286           endif
2287         endif
2288       endif
2289
2290       end
2291
2292
2293 c---------------------------------------------------------------------
2294       SUBROUTINE gaksphe(sphe,r,mstu41)
2295
2296 C...Purpose: to perform sphericity tensor analysis to give sphericity,
2297 C...aplanarity and the related event axes. stolen from jetset ;-)
2298       include 'epos.inc'
2299       dimension sphe(4,3)
2300       DIMENSION SM(3,3),SV(3,3)
2301
2302 C...Calculate matrix to be diagonalized.
2303       NP=0
2304 c      MSTU41=2
2305       DO 110 J1=1,3
2306          DO 100 J2=J1,3
2307             SM(J1,J2)=0.
2308  100     CONTINUE
2309  110  CONTINUE
2310       PS=0.
2311       DO 140 I=1,nptl
2312       IF(istptl(i).ne.0)  GOTO 140
2313       IF(MSTU41.GE.2) THEN
2314          ida=iabs(idptl(i))
2315          IF(ida.EQ.0.OR.ida.EQ.11.OR.ida.EQ.13.OR.ida.EQ.15) GOTO 140
2316          IF(MSTU41.GE.3) then
2317             call idchrg(idptl(i),chrg)
2318             if (abs(chrg).le.0.1) goto 140
2319          endif
2320       ENDIF
2321       NP=NP+1
2322       PA=SQRT(pptl(1,i)**2+pptl(2,I)**2+pptl(3,i)**2)
2323       PWT=1.
2324       IF(ABS(r-2.).GT.0.001) PWT=MAX(1E-10,PA)**(r-2.)
2325       DO 130 J1=1,3
2326          DO 120 J2=J1,3
2327             SM(J1,J2)=SM(J1,J2)+PWT*pptl(j1,i)*pptl(j2,i)
2328  120     CONTINUE
2329  130  CONTINUE
2330       PS=PS+PWT*PA**2
2331  140  CONTINUE
2332
2333 C...Very low multiplicities (0 or 1) not considered.
2334       IF(NP.LE.1) THEN
2335         if(ish.ge.1)then
2336           CALL utmsg('sphe  ')
2337           write(ifch,*) 'too few particles for analysis'
2338           call utmsgf
2339         endif
2340         sphe(4,1)=-1.
2341         RETURN
2342       ENDIF
2343       DO 160 J1=1,3
2344          DO 150 J2=J1,3
2345             SM(J1,J2)=SM(J1,J2)/PS
2346  150     CONTINUE
2347  160  CONTINUE
2348
2349 C...Find eigenvalues to matrix (third degree equation).
2350       SQ=(SM(1,1)*SM(2,2)+SM(1,1)*SM(3,3)+SM(2,2)*SM(3,3)-SM(1,2)**2-
2351      &     SM(1,3)**2-SM(2,3)**2)/3.-1./9.
2352       SR=-0.5*(SQ+1./9.+SM(1,1)*SM(2,3)**2+SM(2,2)*SM(1,3)**2+SM(3,3)*
2353      &     SM(1,2)**2-SM(1,1)*SM(2,2)*SM(3,3))
2354      &     +SM(1,2)*SM(1,3)*SM(2,3)+1./27.
2355       SP=COS(ACOS(MAX(MIN(SR/SQRT(-SQ**3),1.),-1.))/3.)
2356       sphe(4,1)=1./3.+SQRT(-SQ)*MAX(2.*SP,SQRT(3.*(1.-SP**2))-SP)
2357       sphe(4,3)=1./3.+SQRT(-SQ)*MIN(2.*SP,-SQRT(3.*(1.-SP**2))-SP)
2358       sphe(4,2)=1.-sphe(4,1)-sphe(4,3)
2359       IF(sphe(4,2).LT.1E-5) THEN
2360         if(ish.ge.1)then
2361           CALL utmsg('gaksphe')
2362           write(ifch,*) 'all particles back-to-back'
2363           call utmsgf
2364         endif
2365         sphe(4,1)=-1.
2366         RETURN
2367       ENDIF
2368
2369 C...Find first and last eigenvector by solving equation system.
2370       DO 240 I=1,3,2
2371          DO 180 J1=1,3
2372             SV(J1,J1)=SM(J1,J1)-sphe(4,I)
2373             DO 170 J2=J1+1,3
2374                SV(J1,J2)=SM(J1,J2)
2375                SV(J2,J1)=SM(J1,J2)
2376  170        CONTINUE
2377  180     CONTINUE
2378          SMAX=0.
2379          DO 200 J1=1,3
2380             DO 190 J2=1,3
2381                IF(ABS(SV(J1,J2)).LE.SMAX) GOTO 190
2382                JA=J1
2383                JB=J2
2384                SMAX=ABS(SV(J1,J2))
2385  190        CONTINUE
2386  200     CONTINUE
2387          SMAX=0.
2388          DO 220 J3=JA+1,JA+2
2389             J1=J3-3*((J3-1)/3)
2390             RL=SV(J1,JB)/SV(JA,JB)
2391             DO 210 J2=1,3
2392                SV(J1,J2)=SV(J1,J2)-RL*SV(JA,J2)
2393                IF(ABS(SV(J1,J2)).LE.SMAX) GOTO 210
2394                JC=J1
2395                SMAX=ABS(SV(J1,J2))
2396  210        CONTINUE
2397  220     CONTINUE
2398          JB1=JB+1-3*(JB/3)
2399          JB2=JB+2-3*((JB+1)/3)
2400          sphe(JB1,I)=-SV(JC,JB2)
2401          sphe(jb2,I)=SV(JC,JB1)
2402          sphe(jb,I)=-(SV(JA,JB1)*sphe(jb1,I)+SV(JA,JB2)
2403      &        *sphe(jb2,I))/SV(JA,JB)
2404          PA=SQRT(sphe(1,I)**2+sphe(2,I)**2+sphe(3,I)**2)
2405          SGN=(-1.)**INT(rangen()+0.5)
2406          DO 230 J=1,3
2407             sphe(j,I)=SGN*sphe(j,I)/PA
2408  230     CONTINUE
2409  240  CONTINUE
2410
2411 C...Middle axis orthogonal to other two. Fill other codes.
2412       SGN=(-1.)**INT(rangen()+0.5)
2413       sphe(1,2)=SGN*(sphe(2,1)*sphe(3,3)-sphe(3,1)*sphe(2,3))
2414       sphe(2,2)=SGN*(sphe(3,1)*sphe(1,3)-sphe(1,1)*sphe(3,3))
2415       sphe(3,2)=SGN*(sphe(1,1)*sphe(2,3)-sphe(2,1)*sphe(1,3))
2416
2417       do i=1,3
2418          do j=1,4
2419             pptl(j,nptl+i)=sphe(j,i)
2420          enddo
2421       enddo
2422
2423
2424 C...Calculate sphericity and aplanarity. Select storing option.
2425 ccc      SPH=1.5*(sphe(4,2)+sphe(4,3))
2426 ccc      APL=1.5*sphe(4,3)
2427
2428       RETURN
2429       END
2430
2431 C*********************************************************************
2432
2433       SUBROUTINE gakthru(thru,mstu41)
2434
2435 C...Purpose: to perform thrust analysis to give thrust, oblateness
2436 C...and the related event axes. stolen from jetset ;-)
2437       include 'epos.inc'
2438       DIMENSION TDI(3),TPR(3)
2439       dimension thru(4,3)
2440
2441 C...Take copy of particles that are to be considered in thrust analysis.
2442       NP=0
2443       PS=0.
2444 c      MSTU41=2
2445       MSTU44=4
2446       MSTU45=2
2447       PARU42=1.0
2448       PARU48=0.0000001
2449       DO 100 I=1,nptl
2450          IF(istptl(i).ne.0)goto 100
2451          ida=iabs(idptl(i))
2452          IF(ida.EQ.0.OR.ida.EQ.11.OR.ida.EQ.13.OR.ida.EQ.15)GOTO 100
2453          IF(MSTU41.GE.3) then
2454             call idchrg(idptl(i),chrg)
2455             if (abs(chrg).le.0.1) goto 100
2456          endif
2457
2458          IF(nptl+NP.GE.mxptl) THEN
2459             CALL utstop('gakthru: no more memory left in cptl&')
2460             thru(4,1)=-1.
2461             RETURN
2462          ENDIF
2463          NP=NP+1
2464 c         K(nptl+NP,1)=23
2465          pptl(1,nptl+NP)=pptl(1,I)
2466          pptl(2,nptl+NP)=pptl(2,I)
2467          pptl(3,nptl+NP)=pptl(3,I)
2468          pptl(4,nptl+NP)=SQRT(pptl(1,I)**2+pptl(2,I)**2+pptl(3,I)**2)
2469          pptl(5,nptl+NP)=1.
2470          IF(ABS(PARU42-1.).GT.0.001)
2471      &        pptl(5,nptl+NP)=pptl(4,nptl+NP)**(PARU42-1.)
2472          PS=PS+pptl(4,nptl+NP)*pptl(5,nptl+NP)
2473  100  CONTINUE
2474
2475 C...Very low multiplicities (0 or 1) not considered.
2476       IF(NP.LE.1) THEN
2477          CALL utmsg('thru  ')
2478          write(ifch,*) 'too few particles for analysis'
2479          call utmsgf
2480          thru(4,1)=-1
2481          RETURN
2482       ENDIF
2483
2484 C...Loop over thrust and major. T axis along z direction in latter case.
2485       DO 320 ILD=1,2
2486       IF(ILD.EQ.2) THEN
2487 c        PHI=GAKANG(pptl(1,nptl+NP+1),pptl(2,nptl+NP+1))
2488 c        CALL lurot(nptl+1,nptl+NP+1,0.,-PHI)
2489 c        THE=GAKANG(pptl(3,nptl+NP+1),pptl(1,nptl+NP+1))
2490 c        CALL lurot(nptl+1,nptl+NP+1,-THE,0.)
2491          ax=pptl(1,nptl+NP+1)
2492          ay=pptl(2,nptl+NP+1)
2493          az=pptl(3,nptl+NP+1)
2494          do irot=nptl+1,nptl+NP+1
2495             call utrota(1,ax,ay,az
2496      &           ,pptl(1,irot),pptl(2,irot),pptl(3,irot))
2497          enddo
2498          if(np.eq.2)then
2499            pptl(1,nptl+NP+2)=1.
2500            pptl(2,nptl+NP+2)=0.
2501            pptl(3,nptl+NP+2)=0.
2502            pptl(4,nptl+NP+2)=0.
2503            goto 325
2504          endif
2505       ENDIF
2506
2507 C...Find and order particles with highest p (pT for major).
2508       DO 110 ILF=nptl+NP+4,nptl+NP+MSTU44+4
2509          pptl(4,ILF)=0.
2510  110  CONTINUE
2511       DO 160 I=nptl+1,nptl+NP
2512          IF(ILD.EQ.2) pptl(4,I)=SQRT(pptl(1,I)**2+pptl(2,I)**2)
2513          DO 130 ILF=nptl+NP+MSTU44+3,nptl+NP+4,-1
2514             IF(pptl(4,I).LE.pptl(4,ILF)) GOTO 140
2515             DO 120 J=1,5
2516                pptl(j,ILF+1)=pptl(j,ILF)
2517  120        CONTINUE
2518  130     CONTINUE
2519          ILF=nptl+NP+3
2520  140     DO 150 J=1,5
2521             pptl(j,ILF+1)=pptl(j,I)
2522  150     CONTINUE
2523  160  CONTINUE
2524
2525 C...Find and order initial axes with highest thrust (major).
2526       DO 170 ILG=nptl+NP+MSTU44+5,nptl+NP+MSTU44+15
2527          pptl(4,ILG)=0.
2528  170  CONTINUE
2529       NC=2**(MIN(MSTU44,NP)-1)
2530       DO 250 ILC=1,NC
2531          DO 180 J=1,3
2532             TDI(J)=0.
2533  180     CONTINUE
2534          DO 200 ILF=1,MIN(MSTU44,NP)
2535             SGN=pptl(5,nptl+NP+ILF+3)
2536             IF(2**ILF*((ILC+2**(ILF-1)-1)/2**ILF).GE.ILC) SGN=-SGN
2537             DO 190 J=1,4-ILD
2538                TDI(J)=TDI(J)+SGN*pptl(j,nptl+NP+ILF+3)
2539  190        CONTINUE
2540  200     CONTINUE
2541          TDS=TDI(1)**2+TDI(2)**2+TDI(3)**2
2542          DO 220 ILG=nptl+NP+MSTU44+MIN(ILC,10)+4,nptl+NP+MSTU44+5,-1
2543             IF(TDS.LE.pptl(4,ILG)) GOTO 230
2544             DO 210 J=1,4
2545                pptl(j,ILG+1)=pptl(j,ILG)
2546  210        CONTINUE
2547  220     CONTINUE
2548          ILG=nptl+NP+MSTU44+4
2549  230     DO 240 J=1,3
2550             pptl(j,ILG+1)=TDI(J)
2551  240     CONTINUE
2552          pptl(4,ILG+1)=TDS
2553  250  CONTINUE
2554
2555 C...  Iterate direction of axis until stable maximum.
2556       pptl(4,nptl+NP+ILD)=0.
2557       ILG=0
2558  260  ILG=ILG+1
2559       THP=0.
2560  270  THPS=THP
2561       DO 280 J=1,3
2562          IF(THP.LE.1E-10) TDI(J)=pptl(j,nptl+NP+MSTU44+4+ILG)
2563          IF(THP.GT.1E-10) TDI(J)=TPR(J)
2564          TPR(J)=0.
2565  280  CONTINUE
2566       DO 300 I=nptl+1,nptl+NP
2567          SGN=SIGN(pptl(5,I),TDI(1)*pptl(1,I)
2568      &        +TDI(2)*pptl(2,I)+TDI(3)*pptl(3,I))
2569          DO 290 J=1,4-ILD
2570             TPR(J)=TPR(J)+SGN*pptl(j,I)
2571  290     CONTINUE
2572  300  CONTINUE
2573       THP=SQRT(TPR(1)**2+TPR(2)**2+TPR(3)**2)/PS
2574       IF(THP.GE.THPS+PARU48) GOTO 270
2575
2576 C...  Save good axis. Try new initial axis until a number of tries agree.
2577       IF(THP.LT.pptl(4,nptl+NP+ILD)-PARU48.AND.ILG.LT.MIN(10,NC))
2578      &     GOTO 260
2579       IF(THP.GT.pptl(4,nptl+NP+ILD)+PARU48)
2580      $     THEN
2581          IAGR=0
2582          SGN=(-1.)**INT(rangen()+0.5)
2583          DO 310 J=1,3
2584             pptl(j,nptl+NP+ILD)=SGN*TPR(J)/(PS*THP)
2585  310     CONTINUE
2586          pptl(4,nptl+NP+ILD)=THP
2587          pptl(5,nptl+NP+ILD)=0.
2588       ENDIF
2589       IAGR=IAGR+1
2590       IF(IAGR.LT.MSTU45.AND.ILG.LT.MIN(10,NC)) GOTO 260
2591  320  CONTINUE
2592
2593 C...  Find minor axis and value by orthogonality.
2594  325  SGN=(-1.)**INT(rangen()+0.5)
2595       pptl(1,nptl+NP+3)=-SGN*pptl(2,nptl+NP+2)
2596       pptl(2,nptl+NP+3)=SGN*pptl(1,nptl+NP+2)
2597       pptl(3,nptl+NP+3)=0.
2598       THP=0.
2599       DO 330 I=nptl+1,nptl+NP
2600          THP=THP+pptl(5,I)*ABS(pptl(1,nptl+NP+3)*pptl(1,I)
2601      &        +pptl(2,nptl+NP+3)*pptl(2,I))
2602  330  CONTINUE
2603       pptl(4,nptl+NP+3)=THP/PS
2604       pptl(5,nptl+NP+3)=0.
2605
2606
2607 C...  Fill axis information. Rotate back to original coordinate system.
2608       do irot=nptl+NP+1,nptl+NP+3
2609          call utrota(-1,ax,ay,az
2610      &        ,pptl(1,irot),pptl(2,irot),pptl(3,irot))
2611       enddo
2612
2613       do ild=1,3
2614          do j=1,4
2615             thru(j,ild)=pptl(j,nptl+NP+ild)
2616          enddo
2617       enddo
2618
2619 C...Calculate thrust and oblateness. Select storing option.
2620 ccc      THR=thru(4,1)
2621 ccc      OBL=thru(4,2)-thru(4,3)
2622
2623       RETURN
2624       END
2625
2626       subroutine gakjet(ijadu)
2627       include 'epos.inc'
2628       common/njjjj/njets(5,2,mxbins)
2629 c      nmin=xpar1
2630 c      nmax=xpar2
2631       if(nrevt.eq.1)then
2632         do i=1,5
2633           do j=1,mxbins
2634             njets(i,ijadu,j)=0
2635           enddo
2636         enddo
2637       endif
2638       do i=1,nrbins
2639         ycut=xminim*(xmaxim/xminim)**((real(i)-0.5)/nrbins)
2640 c        if(iologe.ne.1)ycut=xminim+(xmaxim-xminim)/nrbins*(nrbins)
2641         nj=gaknjt(ycut,ijadu)
2642         if(nj.ge.1.and.nj.le.5)then
2643           njets(nj,ijadu,i)=njets(nj,ijadu,i)+1
2644         endif
2645       enddo
2646       end
2647
2648       subroutine gakjto
2649       include 'epos.inc'
2650       common/njjjj/njets(5,2,mxbins)
2651       n=xpar4
2652       ijadu=xpar3
2653       do i=1,nrbins
2654         ycut=xminim*(xmaxim/xminim)**((real(i)-0.5)/nrbins)
2655         write (ifhi,*) ycut,real(njets(n,ijadu,i))/nrevt
2656      $       ,sqrt(1.*njets(n,ijadu,i))/nrevt
2657       enddo
2658       end
2659 C*********************************************************************
2660       function gaknjt(ycut,ijadu)
2661 c
2662 c     ijadu   1 =  JADE    2=DURHAM
2663 c     ycut - max. distance
2664 c
2665       include 'epos.inc'
2666
2667       a2j(i,j)=2.*pptl(4,i)*pptl(4,j)*(1.-(pptl(1,i)*pptl(1,j)
2668      &     +pptl(2,i)*pptl(2,j)+pptl(3,i)*pptl(3,j))
2669      &     /(sqrt(pptl(1,i)**2+pptl(2,i)**2+pptl(3,i)**2)
2670      &       *sqrt(pptl(1,j)**2+pptl(2,j)**2+pptl(3,j)**2)))/evis**2
2671
2672       a2d(i,j)=2.*min(pptl(4,i)**2,pptl(4,j)**2)
2673      &     *(1.-(pptl(1,i)*pptl(1,j)
2674      &     +pptl(2,i)*pptl(2,j)+pptl(3,i)*pptl(3,j))
2675      &     /(sqrt(pptl(1,i)**2+pptl(2,i)**2+pptl(3,i)**2)
2676      &       *sqrt(pptl(1,j)**2+pptl(2,j)**2+pptl(3,j)**2)))/evis**2
2677
2678       bet(i)=sqrt(pptl(1,i)**2+pptl(2,i)**2+pptl(3,i)**2)
2679
2680       ska(i,j)=pptl(1,i)*pptl(1,j)+pptl(2,i)*pptl(2,j)
2681      &     +pptl(3,i)*pptl(3,j)
2682
2683       a2c(i,j)= ((bet(i)*bet(j)-ska(i,j))
2684      &     *2.*bet(i)*bet(j)  /  (0.00001+bet(i)+bet(j))**2 )
2685
2686       evis=0.
2687       nn=0
2688       do i=1,nptl
2689          if (istptl(i).eq.0) then
2690             nn=nn+1
2691             do j=1,5
2692                pptl(j,nptl+nn)=pptl(j,i)
2693             enddo
2694             evis=evis+pptl(4,i)
2695             jorptl(i)=nn
2696          endif
2697       enddo
2698       iflag=0
2699       do while (iflag.eq.0.and.nn.ge.2)
2700          a2min=ycut
2701          iflag=1
2702          do i=nptl+1,nptl+nn-1
2703             do j=i+1,nptl+nn
2704                if(ijadu.eq.1)then
2705                   a2=a2j(i,j)
2706                elseif(ijadu.eq.2) then
2707                   a2=a2d(i,j)
2708                else
2709                   a2=a2c(i,j)
2710                endif
2711                if (a2.lt.a2min) then
2712                   iflag=0
2713                   i1=i
2714                   j1=j
2715                   a2min=a2
2716                endif
2717             enddo
2718          enddo
2719          if (iflag.eq.0) then
2720             do i=1,4
2721                pptl(i,i1)=pptl(i,i1)+pptl(i,j1)
2722             enddo
2723             do i=1,nptl
2724                if(istptl(i).eq.0.and.jorptl(i).eq.j1-nptl)
2725      &              jorptl(i)=i1-nptl
2726                if(istptl(i).eq.0.and.jorptl(i)+nptl.gt.j1)
2727      &              jorptl(i)=jorptl(i)-1
2728             enddo
2729             do i=j1,nptl+nn
2730                do j=1,5
2731                   pptl(j,i)=pptl(j,i+1)
2732                enddo
2733             enddo
2734             nn=nn-1
2735          endif
2736       enddo
2737       do i=nptl+1,nptl+nn
2738          istptl(i)=-1
2739          jorptl(i)=i-nptl
2740          pptl(5,i)=sqrt(max(0.,(pptl(4,i)+pptl(3,i))
2741      &   *(pptl(4,i)-pptl(3,i))-pptl(2,i)**2-pptl(1,i)**2))
2742       enddo
2743       do i=nptl+1,nptl+nn-1
2744          do j=i+1,nptl+nn
2745             if(pptl(4,jorptl(i)+nptl).lt.pptl(4,jorptl(j)+nptl))then
2746                k=jorptl(i)
2747                jorptl(i)=jorptl(j)
2748                jorptl(j)=k
2749             endif
2750          enddo
2751       enddo
2752       do i=nptl+1,nptl+nn
2753          idptl(nptl+jorptl(i))=9910+i-nptl
2754       enddo
2755       do i=1,nptl
2756         jorptl(i)=0             !jorptl(nptl+jorptl(i))
2757       enddo
2758 c      nptl=nptl+nn
2759
2760       gaknjt=nn
2761       return
2762       end
2763
2764       subroutine idtr5(id,ic)
2765       integer ic(2)
2766       ic(1)=0
2767       ic(2)=0
2768       ii=1
2769       if(id.lt.0)ii=2
2770       i1=1
2771       if(iabs(id).gt.999)i1=3
2772       do i=i1,int(log(abs(real(id)))/log(10.))+1
2773         j=mod(iabs(id)/10**(i-1),10)
2774         if(j.gt.0)then
2775           ic(ii)=ic(ii)+10**(6-j)
2776         endif
2777       enddo
2778       return
2779       end
2780
2781       function ammin(id1,id2)
2782       dimension ic1(2),ic2(2),jc2(6,2),jc1(6,2)
2783       call idtr5(id1,ic1)
2784       call idtr5(id2,ic2)
2785       call idcomk(ic1)
2786       call idcomk(ic2)
2787       call iddeco(ic1,jc1)
2788       call iddeco(ic2,jc2)
2789       ammin=utamnx(jc1,jc2)
2790       end
2791
2792
2793 c      function idtr(id1,id2)
2794 c      integer ic(2),id(2)
2795 c      id(1)=id1
2796 c      id(2)=id2
2797 c      do i=1,2
2798 c        ic(i)=0
2799 c      enddo
2800 c      do j=1,2
2801 c        ii=1
2802 c        if(id(j).lt.0)ii=2
2803 c        i1=1
2804 c        if(iabs(id(j)).gt.999)i1=3
2805 c        do i=i1,int(log(abs(real(id(j))))/log(10.))+1
2806 c          jj=mod(iabs(id(j))/10**(i-1),10)
2807 c          if(jj.gt.0)then
2808 c            ic(ii)=ic(ii)+10**(6-jj)
2809 c          endif
2810 c        enddo
2811 c      enddo
2812 c      idtr=idtra(ic,0,0,3)
2813 c      if(idtr.ne.idsp(id1,id2))then
2814 c        write (*,'(4i6)') idtr,idsp(id1,id2),id1,id2
2815 c      endif
2816 c      return
2817 c      end
2818
2819       function idsp(id1,id2)
2820       ia1=iabs(id1)
2821       ia2=iabs(id2)
2822       if(ia1.ge.1000.and.ia2.ge.1000)then
2823         idsp=0
2824         isign=0
2825       elseif(ia1.le.1000.and.ia2.le.1000)then
2826         idsp=min(ia1,ia2)*100+max(ia1,ia2)*10
2827         isign=1
2828         if(max(ia1,ia2).ne.-min(id1,id2)) isign = -1
2829         if(idsp.eq.220)idsp=110
2830         if(idsp.eq.330)idsp=220
2831       else
2832         isign=1
2833         if(id1.lt.0.and.id2.lt.0)isign=-1
2834         idb=min(ia1,ia2)
2835         if(idb.eq.5)then
2836           idsp=0
2837           return
2838         endif
2839         ida=max(ia1,ia2)
2840         ida1=ida/1000
2841         ida2=mod(ida/100,10)
2842         if(idb.le.ida1)then
2843           idsp=idb*1000+ida/10
2844         elseif(idb.le.ida2)then
2845           idsp=ida1*1000+idb*100+ida2*10
2846         else
2847           idsp=ida+idb*10
2848         endif
2849         if(ida1.eq.ida2.and.ida2.eq.idb)idsp=idsp+1
2850       endif
2851       idsp=idsp*isign
2852       return
2853       end
2854