EPOS/epos167/epos-uti-165.f

   1 c-----------------------------------------------------------------------
   2       subroutine utresc(iret)
   3 c-----------------------------------------------------------------------
   4 c  if irescl=1 rescaling is done, otherwise the purpose of going through
   5 c  this routine is to not change the seed in case of no rescaling
   6 c-----------------------------------------------------------------------
   7       include 'epos.inc'
   8       include 'epos.incems'
   9       double precision p1,esoll,ppp,seedp
  10       dimension p1(5),p0(5,mamx+mamx),pini(mxptl)
  11       logical force,nolead(mxptl),lspec(mxptl)
  12       data scalmean/0./
  13       save scalmean
  14       call utpri('utresc',ish,ishini,4)
  15
  16       errlim=0.001 !max(0.001,1./engy)
  17
  18       iret=0
  19       scal=1.
  20       nptlpt=iabs(maproj)+iabs(matarg)
  21       call ranfgt(seedp)        !not to change the seed ...
  22       if(nptl.le.nptlpt) goto 9999
  23
  24       esoll=0.d0
  25       p1(1)=0.d0
  26       p1(2)=0.d0
  27       p1(3)=0.d0
  28       p1(4)=0.d0
  29       ipmax=4
  30       imin=nptlpt+1
  31       if(iappl.eq.1)then
  32         imin=1
  33         ipmax=2
  34       endif
  35       do i=1,nptlpt
  36          nolead(i)=.false.
  37          do j=1,5
  38            p0(j,i)=pptl(j,i)
  39          enddo
  40 c         if(mod(istptl(i),10).eq.1)then
  41            do j=ipmax+1,4
  42              p1(j)=p1(j)+dble(pptl(j,i))
  43            enddo
  44 c         endif
  45       enddo
  46       do i=nptlpt+1,nptl
  47        if(mod(istptl(i),10).eq.0)then
  48          do j=1,ipmax
  49            p1(j)=p1(j)+dble(pptl(j,i))
  50          enddo
  51          lspec(i)=.false.
  52          if((ityptl(i).eq.45.and.maproj.ge.100)
  53      &  .or.(ityptl(i).eq.55.and.matarg.ge.100))lspec(i)=.true.
  54          if((abs(pptl(3,i)/pnullx).le.0.9
  55      & .and.abs(pptl(3,i)).gt.pptl(5,i)).or.lspec(i))then
  56            nolead(i)=.true.
  57 c           write(ifch,*)'nolead',i
  58          else
  59            nolead(i)=.false.
  60 c           write(ifch,*)'lead',i
  61          endif
  62        endif
  63       enddo
  64       ppp=(p1(4)+p1(3))*(p1(4)-p1(3))-p1(2)*p1(2)-p1(1)*p1(1)
  65       if(ppp.gt.0.d0)then
  66         p1(5)=dsqrt(ppp)
  67       else
  68         iret=1
  69         write(ifch,*)'p1',p1(1),p1(2),p1(3),p1(4),ppp
  70         if(ish.ge.2)write (ifch,*) 'Problem in utresc (0): redo event'
  71         write (ifmt,*) 'Problem in utresc (0): redo event'
  72 c       call utstop('utresc&')
  73         goto 9999
  74       endif
  75       esoll=p1(5)
  76       if(ish.ge.4) write (ifch,'(a,5g13.6)') 'boost-vector: ',p1
  77
  78 c     trafo
  79 c     -----
  80       pmax=0.d0
  81       do i=imin,nptl
  82         if(mod(istptl(i),10).le.1)then
  83           call utlob4(1,p1(1),p1(2),p1(3),p1(4),p1(5)
  84      $         ,pptl(1,i),pptl(2,i),pptl(3,i),pptl(4,i))
  85         endif
  86         if(i.gt.nptlpt.and.nolead(i))pmax=max(pmax,abs(pptl(3,i)))
  87       enddo
  88
  89       if(ish.ge.6)then
  90         call alistf('list after boost&')
  91       endif
  92
  93       if(ish.ge.5)write(ifch,'(a)')'--------rescale momenta----------'
  94
  95 c     rescale momenta in rest frame
  96 c     -----------------------------
  97       scal=1.
  98       diff0=0.
  99 c      ndif0=1
 100       ferr=0.05
 101       force=.false.
 102       npart=nptl-imin+1
 103       do ipass=1,300
 104         sum=0.
 105         sum3=0.
 106         difft=diff0
 107         ndif=0
 108         nfirst=int(rangen()*float(npart))    !start list randomly
 109         do  i=0,npart-1
 110           j=imin+i+nfirst
 111           if(j.gt.nptl)j=imin+i+nfirst-npart
 112           if(mod(istptl(j),10).eq.0)then
 113 c modified particles
 114             if(nolead(j))then
 115 c            if(j.gt.nptlpt)then
 116 c            if(abs(pptl(3,j))/pnullx.lt.0.9)then  !not spectator or diffraction
 117               if(scal.eq.1..and.diff0.eq.0.)then
 118                 ndif=ndif+1
 119                 pini(j)=pptl(3,j)
 120               else
 121                 pptl3new=0.
 122                 if( force .or.(
 123      &            ityptl(j)/10.eq.4.or.ityptl(j)/10.eq.5
 124      &                      ))then !try just remnant first
 125                   ndif=ndif+1
 126                   diff=sign(min(0.3*abs(pini(j)),
 127      &                      rangen()*abs(difft)),difft)
 128                   pptl3new=scal*(pptl(3,j)-diff)
 129 c                  write(ifch,*)'par',j,pptl3new,pptl(3,j),diff,difft
 130 c     &                 ,ndif,pmax,scal
 131                   if(abs(pptl3new).lt.pmax)then
 132                     if(abs(pptl3new-pini(j)).lt.ferr*abs(pini(j))
 133      &         .or.(lspec(j).and.abs(pptl3new).lt.abs(0.8*pini(j))))then !particle should not be too fast or too modified
 134 c                  write(ifch,*)'used'
 135                       difft=difft-diff
 136                       pptl(3,j)=scal*(pptl(3,j)-diff)
 137                       pptl(4,j)=sqrt(pptl(1,j)**2+pptl(2,j)**2
 138      *                     +pptl(3,j)**2+pptl(5,j)**2)
 139                     endif
 140                   endif
 141                 endif
 142               endif
 143             endif
 144 c sum over all particles
 145             sum=sum+pptl(4,j)
 146             sum3=sum3+pptl(3,j)
 147           endif
 148         enddo
 149
 150         diff=sum3
 151         scal=real(esoll)/(sum-diff)
 152         if(ish.ge.6)write(ifch,*)
 153      $       'ipass,scal,diff/pnullx,e,esoll,pz,pmax,ndif,f:'
 154      $   ,ipass,scal,diff/pnullx,sum,esoll,sum3,pnullx,ndif,force
 155         if(abs(scal-1.).le.errlim.and.abs(diff/pnullx).lt.errlim)
 156      $  goto 300
 157         if(ndif.gt.0.and.(force.or.ipass.lt.150))then
 158 c          ndif0=ndif
 159           diff0=diff
 160         elseif(abs(scal-1.).le.1e-2.and.abs(diff/pnullx).lt.5e-2)then
 161           goto 300
 162         elseif(force)then
 163           if(ish.ge.2)
 164      $    write(ifmt,*)'Warning in utresc: no more allowed particle'
 165           goto 302
 166         else
 167           force=.true.
 168           ferr=0.1
 169           diff=0.
 170         endif
 171  301  enddo
 172  302  if(ish.ge.1)then
 173         call utmsg('utrescl')
 174         write(ifch,*)'*****  scal=',scal,diff
 175         call utmsgf
 176       endif
 177       if(abs(scal)+abs(diff/pnullx).gt.2.5)then
 178         write(ifmt,*)'Warning in utresc !'
 179         write(ifch,*)'redo EPOS event ...'
 180         iret=1
 181         goto 9999
 182       endif
 183
 184 c     trafo
 185 c     -----
 186  300  continue
 187       do i=1,nptl
 188         if(i.le.nptlpt)then
 189           do j=1,5
 190             pptl(j,i)=p0(j,i)
 191           enddo
 192         else
 193           if(mod(istptl(i),10).le.1)then
 194             call utlob4(-1,p1(1),p1(2),p1(3),p1(4),p1(5)
 195      $           ,pptl(1,i),pptl(2,i),pptl(3,i),pptl(4,i))
 196           endif
 197         endif
 198       enddo
 199
 200       if(ish.ge.4)call alist('list after rescaling&',1,nptl)
 201
 202  9999 continue
 203       if(ish.ge.2)then
 204         scalmean=scalmean+scal
 205         write(ifch,*)' average rescaling factor: ',scalmean
 206      &                                            /float(nrevt+1)
 207       endif
 208       call ranfst(seedp)        ! ... after this subroutine
 209       call utprix('utresc',ish,ishini,4)
 210
 211       end
 212
 213 c-----------------------------------------------------------------------
 214       subroutine utghost(iret)
 215 c-----------------------------------------------------------------------
 216 c  if irescl=1 make particle on-shell if not
 217 c-----------------------------------------------------------------------
 218       include 'epos.inc'
 219       include 'epos.incems'
 220       double precision seedp
 221       call utpri('ughost',ish,ishini,4)
 222
 223       iret=0
 224       nptlpt=iabs(maproj)+iabs(matarg)
 225       call ranfgt(seedp)        !not to change the seed ...
 226       if(nptl.le.nptlpt) goto 9999
 227
 228       if(ish.ge.5)write(ifch,'(a)')'---------mark ghosts---------'
 229
 230 c     mark ghosts
 231 c     -----------
 232       do  j=nptlpt+1,nptl
 233         if(mod(istptl(j),10).eq.0.and.pptl(4,j).gt.0.d0)then
 234           amass=pptl(5,j)
 235           call idmass(idptl(j),amass)
 236           if(abs(pptl(5,j)-amass).gt.0.01*amass)then
 237             if(ish.ge.5)write(ifch,*)'wrong particle mass',j,idptl(j)
 238      &                                           ,pptl(5,j),amass
 239             amass=pptl(5,j)
 240             call idres(idptl(j),amass,idr,iadj,0)
 241             if(idr.ne.0)then
 242               pptl(5,j)=amass
 243               idptl(j)=idr
 244             else
 245               call idmass(idptl(j),amass)
 246               pptl(5,j)=amass
 247             endif
 248             call idtau(idptl(j),pptl(4,j),pptl(5,j),taugm)
 249             tivptl(2,j)=tivptl(1,j)+taugm*(-alog(rangen()))
 250           endif
 251           if(abs((pptl(4,j)+pptl(3,j))*(pptl(4,j)-pptl(3,j))
 252      $         -pptl(2,j)**2-pptl(1,j)**2-pptl(5,j)**2).gt.0.3
 253      $       .and.abs(1.-abs(pptl(3,j))/pptl(4,j)).gt.0.01)then
 254         !print*,'ghost',ityptl(j),idptl(j)
 255            if(ish.ge.1)write(ifmt,*)'ghost:',j,idptl(j),ityptl(j)
 256            if(ish.ge.5)then
 257               write(ifch,'(a,$)')'ghost:'
 258               call alistc("&",j,j)
 259             endif
 260             ityptl(j)=100+ityptl(j)/10
 261           endif
 262         elseif(mod(istptl(j),10).eq.0)then
 263            if(ish.ge.1)then
 264              write(ifmt,*)'Lost particle (E=0)'
 265              write(ifch,*)'Lost particle (E=0) :'
 266              call alistc("utghost&",j,j)
 267            endif
 268            istptl(j)=istptl(j)+2
 269         endif
 270       enddo
 271
 272       if(ish.ge.5)write(ifch,'(a)')'---------treat ghosts---------'
 273
 274 c     treat ghosts
 275 c     ------------
 276       ifirst=1
 277       scal=1.
 278       pfif=0.
 279       efif=0.
 280       ntry=0
 281  132  nfif=0
 282       psum=0
 283       esum=0.
 284       ntry=ntry+1
 285       do  j=nptlpt+1,nptl
 286         if(mod(istptl(j),10).eq.0)then
 287           if(ityptl(j).ge.101.and.ityptl(j).le.105)then
 288             nfif=nfif+1
 289             if(ifirst.eq.1)then
 290               pfif=pfif+pptl(3,j)
 291               if(pptl(4,j).gt.0.)efif=efif+pptl(4,j)
 292             endif
 293             if(irescl.ge.1) then
 294               if(ifirst.gt.1)then
 295                 if(pptl(4,j).gt.0.)then
 296                   Einv=1./pptl(4,j)
 297                   amt=1.-(pptl(5,j)*Einv)**2+(pptl(1,j)*Einv)**2
 298      $                +(pptl(2,j)*Einv)**2
 299                 else
 300                   amt=-1.
 301                 endif
 302                 if(amt.gt.0.)then
 303                   pptl(3,j)=sign(pptl(4,j),pptl(3,j))*sqrt(amt)
 304                 else
 305                   y=(rangen()+rangen()+rangen()+rangen()-2.)/2.*yhaha
 306                   y=sign(abs(y),pptl(3,j))
 307                   pptl(3,j)
 308      $                 =sqrt(pptl(5,j)**2+pptl(1,j)**2
 309      $                 +pptl(2,j)**2)*sinh(y)
 310                   pptl(4,j)
 311      $                 =sqrt(pptl(5,j)**2+pptl(1,j)**2
 312      $                 +pptl(2,j)**2)*cosh(y)
 313                   efif=efif+pptl(4,j)
 314                 endif
 315                 ifirst=0
 316               else
 317 c                do k=1,3
 318                 do k=3,3
 319                   pptl(k,j)=pptl(k,j)*scal
 320                 enddo
 321                 pptl(4,j)=sqrt(pptl(1,j)**2+pptl(2,j)**2+pptl(3,j)**2
 322      *                 +pptl(5,j)**2)
 323               endif
 324             endif
 325             psum=psum+pptl(3,j)
 326             esum=esum+pptl(4,j)
 327             if(ish.ge.5)
 328      $           write (ifch,*) 'nrevt,psum,esum,pfif,efif,nfif,scal'
 329      $           ,nrevt,psum,esum,pfif,efif,nfif,scal
 330           endif
 331         endif
 332       enddo
 333       if ( ish.ge.5 )  write (ifch,*) 'tot',nfif,efif,pfif,esum,psum
 334
 335
 336       if(nfif.gt.5.or.(esum.gt.0.05*engy.and.nfif.ne.1))then
 337         if(ifirst.eq.0)then
 338           do  j=nptlpt+1,nptl
 339             if ( ityptl(j).ge.101 .and. ityptl(j).le.105 )then
 340               if((psum-pfif)*(1.-scal).ge.0)
 341      &             pptl(3,j)=pptl(3,j)-(psum-pfif)/nfif
 342             endif
 343           enddo
 344         else
 345           ifirst=2
 346           goto 132
 347         endif
 348         scal=efif/esum
 349         if ( ish.ge.5 )  write (ifch,*) 'scal',scal
 350         if ( abs(scal-1.) .gt. 0.05 ) then
 351           if(ntry.le.1000)then
 352             goto 132
 353           else
 354             iret=1
 355             if(ish.ge.2)write (ifch,*) 'Problem in utghost : redo event'
 356             if(ish.ge.1)write (ifmt,*) 'Problem in utghost : redo event'
 357             goto 9999
 358          endif
 359         endif
 360       else
 361         do  j=nptlpt+1,nptl
 362           if ( ityptl(j).ge.101 .and. ityptl(j).le.105 )then
 363             pptl(4,j)=sqrt(pptl(1,j)**2+pptl(2,j)**2+pptl(3,j)**2
 364      *                 +pptl(5,j)**2)
 365           endif
 366         enddo
 367       endif
 368
 369       if(ish.ge.5)write(ifch,'(a)')'---------Check Ghost list---------'
 370
 371 c Check Ghost list
 372
 373       if(ish.ge.5)then
 374         do  j=nptlpt+1,nptl
 375           if(mod(istptl(j),10).eq.0)then
 376             if(ityptl(j).le.105.and.ityptl(j).ge.101)then
 377               write(ifch,'(a,$)')'ghost:'
 378               call alistc("&",j,j)
 379             endif
 380           endif
 381         enddo
 382       endif
 383
 384
 385  9999 continue
 386       call ranfst(seedp)        ! ... after this subroutine
 387       call utprix('ughost',ish,ishini,4)
 388
 389       end
 390
 391 c-----------------------------------------------------------------------
 392       subroutine utrsph(iret)
 393 c-----------------------------------------------------------------------
 394 c  if irescl=1 and ispherio=1 rescaling is done for particle used by
 395 c  spherio as initial condition.
 396 c-----------------------------------------------------------------------
 397       include 'epos.inc'
 398       include 'epos.incems'
 399       double precision p1,esoll
 400       dimension p1(5),p0(5,mamx+mamx)
 401       call utpri('utrsph',ish,ishini,4)
 402
 403       errlim=0.0001
 404
 405       iret=0
 406       nptlpt=iabs(maproj)+iabs(matarg)
 407       if(nptl.le.nptlpt) goto 9999
 408
 409       esoll=0.d0
 410       p1(1)=0.d0
 411       p1(2)=0.d0
 412       p1(3)=0.d0
 413       p1(4)=0.d0
 414       do i=nptlpt+1,nptl
 415         if((istptl(i).le.11
 416      $   .and.(iorptl(i).ge.1.and.istptl(iorptl(i)).eq.41))
 417      $   .or.istptl(i).eq.20.or.istptl(i).eq.21)then
 418          do j=1,2
 419            p1(j)=p1(j)+dble(pptl(j,i))
 420          enddo
 421        endif
 422       enddo
 423       do i=1,nptlpt
 424          do j=1,5
 425            p0(j,i)=pptl(j,i)
 426          enddo
 427          do j=3,4
 428            p1(j)=p1(j)+dble(pptl(j,i))
 429          enddo
 430       enddo
 431       p1(5)=dsqrt((p1(4)+p1(3))*(p1(4)-p1(3))-p1(2)**2.d0-p1(1)**2.d0)
 432       esoll=p1(5)
 433       if(ish.ge.4) write (ifch,'(a,5g13.6)') 'boost-vector',p1
 434
 435 c     trafo
 436 c     -----
 437       do i=1,nptl
 438         if((istptl(i).le.11
 439      $   .and.(iorptl(i).ge.1.and.istptl(iorptl(i)).eq.41))
 440      $   .or.istptl(i).eq.20.or.istptl(i).eq.21
 441      $   .or.(istptl(i).eq.0.and.i.le.nptlpt))then
 442           call utlob4(1,p1(1),p1(2),p1(3),p1(4),p1(5)
 443      $         ,pptl(1,i),pptl(2,i),pptl(3,i),pptl(4,i))
 444         endif
 445       enddo
 446
 447
 448       if(ish.ge.5)write(ifch,'(a)')'------------------'
 449
 450 c     rescale momenta in rest frame
 451 c     -----------------------------
 452
 453       scal=1.
 454       diff=0.
 455       do ipass=1,1000
 456         sum=0.
 457         sum3=0.
 458         ndif=0
 459         do  j=1,nptl
 460         if((istptl(j).le.11
 461      $   .and.(iorptl(j).ge.1.and.istptl(iorptl(j)).eq.41))
 462      $   .or.istptl(j).eq.20.or.istptl(j).eq.21
 463      $   .or.(istptl(j).eq.0.and.j.le.nptlpt))then
 464             if(j.gt.nptlpt)then
 465               ndif=ndif+1
 466               pptl(3,j)=scal*(pptl(3,j)-diff)
 467               pptl(4,j)=sqrt(pptl(1,j)**2+pptl(2,j)**2+pptl(3,j)**2
 468      *           +pptl(5,j)**2)
 469             endif
 470             sum=sum+pptl(4,j)
 471             sum3=sum3+pptl(3,j)
 472           endif
 473         enddo
 474
 475         diff=sum3/real(ndif)
 476         scal=real(esoll)/sum
 477         if(ish.ge.6)write(ifch,*)'ipass,scal,diff,e,esoll,pz,ndif:'
 478      $       ,ipass,scal,diff,sum,esoll,sum3,ndif
 479         if(abs(scal-1.).le.errlim.and.abs(diff).lt.10.*errlim) goto300
 480  301  enddo
 481       if(ish.ge.1)then
 482         call utmsg('hresph')
 483         write(ifch,*)'*****  scal=',scal,diff
 484         call utmsgf
 485       endif
 486
 487
 488 c     trafo
 489 c     -----
 490  300  continue
 491 c      do i=nptlpt+1,nptl
 492       do i=1,nptl
 493         if((istptl(i).le.11
 494      $   .and.(iorptl(i).ge.1.and.istptl(iorptl(i)).eq.41))
 495      $   .or.istptl(i).eq.20.or.istptl(i).eq.21
 496      $   .or.(istptl(i).eq.0.and.i.le.nptlpt))then
 497           call utlob4(-1,p1(1),p1(2),p1(3),p1(4),p1(5)
 498      $         ,pptl(1,i),pptl(2,i),pptl(3,i),pptl(4,i))
 499         endif
 500         if(i.le.nptlpt)then
 501           do j=1,5
 502             pptl(j,i)=p0(j,i)
 503           enddo
 504         endif
 505       enddo
 506
 507  9999 call utprix('utrsph',ish,ishini,4)
 508
 509       end
 510
 511 cc-----------------------------------------------------------------------
 512 c      double precision function dddlog(xxx)
 513 cc-----------------------------------------------------------------------
 514 c      double precision xxx
 515 c      dddlog=-1d50
 516 c      if(xxx.gt.0d0)dddlog=dlog(xxx)
 517 c      end
 518 c
 519 ccc-----------------------------------------------------------------------
 520 c      subroutine randfl(jc,iqa0,iflav,ic,isame)
 521 cc-----------------------------------------------------------------------
 522 cc     returns random flavour ic(2) (iqa0=1:quark,2:antiquark,11:diquark)
 523 cc-----------------------------------------------------------------------
 524 c      include 'epos.inc'
 525 c      real probab(nflav),probsu(nflav+1)
 526 c      integer jc(nflav,2),jc0(nflav,2),ic(2)
 527 c      if(ish.ge.6)then
 528 c      write(ifch,*)('-',i=1,10)
 529 c     *,' entry sr randfl ',('-',i=1,30)
 530 c      write(ifch,*)'iqa0:',iqa0
 531 c      write(ifch,*)'jc:'
 532 c      write(ifch,*)jc
 533 c      endif
 534 c      iflav=0
 535 c      ic(1)=0
 536 c      ic(2)=0
 537 c      do 10 n=1,nflav
 538 c      do 10 i=1,2
 539 c10    jc0(n,i)=0
 540 c      iqa1=iqa0*10
 541 c9999  iqa1=iqa1/10
 542 c      if(iqa1.eq.0)goto9998
 543 c      iqa=mod(iqa1,10)
 544 c      su=0
 545 c      do 20 i=1,nflav
 546 c      probab(i)=jc(i,iqa)-jc0(i,iqa)
 547 c      if(isame.eq.1)probab(i)=probab(i)*(jc(i,3-iqa)-jc0(i,3-iqa))
 548 c20    su=su+probab(i)
 549 c      if(su.lt..5)then
 550 c      iflav=0
 551 c      ic(1)=0
 552 c      ic(2)=0
 553 c      goto9998
 554 c      endif
 555 c      probsu(1)=0.
 556 c      do 30 i=1,nflav
 557 c      probsu(i+1)=probsu(i)+probab(i)/su
 558 c      if(probsu(i+1)-probsu(i).lt.1e-5)probsu(i+1)=probsu(i)
 559 c30    continue
 560 c      r=rangen()*probsu(nflav+1)
 561 c      do 50 i=1,nflav
 562 c      if(probsu(i).le.r.and.r.lt.probsu(i+1))iflav=i
 563 c50    continue
 564 c      jc0(iflav,iqa)=jc0(iflav,iqa)+1
 565 c      if(isame.eq.1)jc0(iflav,3-iqa)=jc0(iflav,3-iqa)+1
 566 c      call idenco(jc0,ic,ireten)
 567 c      if(ireten.eq.1)call utstop('randfl: idenco ret code = 1&')
 568 c      if(ish.ge.6)then
 569 c      write(ifch,*)'probab:'
 570 c      write(ifch,*)probab
 571 c      write(ifch,*)'probsu:'
 572 c      write(ifch,*)probsu
 573 c      write(ifch,*)'ran#:',r,'   flav:',iflav
 574 c      endif
 575 c      goto9999
 576 c9998  continue
 577 c      if(ish.ge.6)write(ifch,*)('-',i=1,30)
 578 c     *,' exit sr randfl ',('-',i=1,10)
 579 c      return
 580 c      end
 581 c
 582 c
 583 cc-----------------------------------------------------------------------
 584 c      subroutine ranhvy(x,eps)
 585 cc-----------------------------------------------------------------------
 586 cc     generates x for heavy particle fragmentation according to
 587 cc     the peterson form
 588 cc          d(x)=1/(x*(1-1/x-eps/(1-x))**2)
 589 cc              =d0(x)*d1(x)*d2(x)
 590 cc          d0(x)=(1-x)**2/((1-x)**2+eps)**2
 591 cc          d1(x)=x
 592 cc          d2(x)=(((1-x)**2+eps)/((1-x)**2+eps*x))**2
 593 cc     using x=1-y**pow
 594 cc     generates flat in x if eps>1.
 595 cc-----------------------------------------------------------------------
 596 c      data aln4/1.3863/
 597 c      if(eps.lt.1.) then
 598 c        pow=alog((3.+eps)/eps)/aln4
 599 c        ymx=(eps*(3.*pow-1.)/(pow+1.))**(.5/pow)
 600 c        zmx=1-ymx**pow
 601 c        d0mx=(1-zmx)**2/((1.-zmx)**2+eps)**2*pow*ymx**(pow-1.)
 602 c        d2mx=2./(2.-sqrt(eps))
 603 c      else
 604 c        pow=1.
 605 c        zmx=0.
 606 c        d0mx=(1.-zmx)**2/((1.-zmx)**2+eps)**2
 607 c        d2mx=1.+eps
 608 c      endif
 609 cc
 610 cc          generate z according to (1-z)**2/((1-z)**2+eps*z)**2
 611 c1     continue
 612 c      y=rangen()
 613 c      z=1.-y**pow
 614 cc
 615 c      d0z=(1.-z)**2/((1.-z)**2+eps)**2*pow*y**(pow-1.)
 616 c      if(d0z.lt.rangen()*d0mx) goto1
 617 cc
 618 cc          check remaining factors
 619 c      d1=z
 620 c      d2=(((1.-z)**2+eps)/((1.-z)**2+eps*z))**2
 621 c      if(d1*d2.lt.rangen()*d2mx) goto1
 622 cc
 623 cc          good x
 624 c      x=z
 625 c      return
 626 c      end
 627 c
 628 c-----------------------------------------------------------------------
 629       function ransig()
 630 c-----------------------------------------------------------------------
 631 c     returns randomly +1 or -1
 632 c-----------------------------------------------------------------------
 633       ransig=1
 634       if(rangen().gt.0.5)ransig=-1
 635       return
 636       end
 637
 638 cc-----------------------------------------------------------------------
 639 c      function ranxq(n,x,q,xmin)
 640 cc-----------------------------------------------------------------------
 641 cc     returns random number according to x(i) q(i) with x>=xmin
 642 cc-----------------------------------------------------------------------
 643 c      include 'epos.inc'
 644 c      real x(n),q(n)
 645 c      imin=1
 646 c      if(xmin.eq.0.)goto3
 647 c      i1=1
 648 c      i2=n
 649 c1     i=i1+(i2-i1)/2
 650 c      if(x(i).lt.xmin)then
 651 c      i1=i
 652 c      elseif(x(i).gt.xmin)then
 653 c      i2=i
 654 c      else
 655 c      imin=i
 656 c      goto3
 657 c      endif
 658 c      if(i2-i1.gt.1)goto1
 659 c      imin=i2
 660 c3     continue
 661 c      if(q(imin).gt.q(n)*.9999)then
 662 c      ranxq=xmin
 663 c      goto4
 664 c      endif
 665 c      qran=q(imin)+rangen()*(q(n)-q(imin))
 666 c      ranxq=utinvt(n,x,q,qran)
 667 c4     continue
 668 c
 669 c      if(ranxq.lt.xmin)then
 670 c      call utmsg('ranxq ')
 671 c      write(ifch,*)'*****  ranxq=',ranxq,' <       xmin=',xmin
 672 c      write(ifch,*)'q(imin) q q(n):',q(imin),qran,q(n)
 673 c      write(ifch,*)'x(imin) x x(n):',x(imin),ranxq,x(n)
 674 c      call utmsgf
 675 c      ranxq=xmin
 676 c      endif
 677 c
 678 c      return
 679 c      end
 680 c
 681 cc  ***** end r-routines
 682 cc  ***** beg s-routines
 683 c
 684 cc-----------------------------------------------------------------------
 685 c      function sbet(z,w)
 686 cc-----------------------------------------------------------------------
 687 c      sbet=utgam1(z)*utgam1(w)/utgam1(z+w)
 688 c      return
 689 c      end
 690 c
 691 cc-----------------------------------------------------------------------
 692 c      function smass(a,y,z)
 693 cc-----------------------------------------------------------------------
 694 cc     returns droplet mass (in gev) (per droplet, not (!) per nucleon)
 695 cc     according to berger/jaffe mass formula, prc35(1987)213 eq.2.31,
 696 cc     see also c. dover, BNL-46322, intersections-meeting, tucson, 91.
 697 cc     a: massnr, y: hypercharge, z: charge,
 698 cc-----------------------------------------------------------------------
 699 c      common/cmass/thet,epsi,as,ac,dy,dz,ym,cz,zm,sigma,rzero
 700 c      ymin=ym*a
 701 c      zmin=cz/(dz/a+zm/a**(1./3.))
 702 c      smass=epsi*a+as*a**(2./3.)+(ac/a**(1./3.)+dz/a/2.)*(z-zmin)**2
 703 c     *+dy/a/2.*(y-ymin)**2
 704 c      return
 705 c      end
 706 c
 707 cc-----------------------------------------------------------------------
 708 c      subroutine smassi(theta)
 709 cc-----------------------------------------------------------------------
 710 cc     initialization for smass.
 711 cc     calculates parameters for berger/jaffe mass formula
 712 cc     (prc35(1987)213 eq.2.31, see also c. dover, BNL-46322).
 713 cc     theta: parameter that determines all parameters in mass formula.
 714 cc-----------------------------------------------------------------------
 715 c      common/cmass/thet,epsi,as,ac,dy,dz,ym,cz,zm,sigma,rzero
 716 c      thet=theta
 717 c
 718 c      astr=.150
 719 c      pi=3.14159
 720 c      alp=1./137.
 721 c
 722 c      co=cos(theta)
 723 c      si=sin(theta)
 724 c      bet=(1+co**3)/2.
 725 c      rzero=si/astr/(  2./3./pi*(1+co**3)  )**(1./3.)
 726 cctp060829      cs=astr/si
 727 c      cz=-astr/si*(  (  .5*(1+co**3)  )**(1./3.)-1  )
 728 c      sigma=6./8./pi*(astr/si)**3*(co**2/6.-si**2*(1-si)/3.-
 729 c     *1./3./pi*(pi/2.-theta-sin(2*theta)+si**3*alog((1+co)/si)))
 730 c
 731 c      epsi=astr*((.5*(1+co**3))**(1./3.)+2)/si
 732 c      as=4*pi*sigma*rzero**2
 733 c      ac=3./5.*alp/rzero
 734 c      dz=astr/si*bet**(1./3.)*co**2*
 735 c     *(co**4*(1+bet**(2./3.))+(1+bet)**2)/
 736 c     *(  (2*co**2+bet**(1./3.))*(co**4*(1+bet**(2./3.))+(1+bet)**2)-
 737 c     *(co**4+bet**(1./3.)*(1+bet))*((2*bet**(2./3.)-1)*co**2+1+bet)  )
 738 c      dy=astr/6.*(1+co**3)**3/si*
 739 c     *(  1+(1+co)/(4*(1+co**3))**(2./3.)  )/
 740 c     *(co**6+co+co*(.5*(1+co**3))**(4./3.))
 741 c      zm=6*alp/(5*rzero)
 742 c      ym=(1-co**3)/(1+co**3)
 743 c
 744 c      return
 745 c      end
 746 c
 747 cc-----------------------------------------------------------------------
 748 c      subroutine smassp
 749 cc-----------------------------------------------------------------------
 750 cc     prints smass.
 751 cc-----------------------------------------------------------------------
 752 c      include 'epos.inc'
 753 c      common/cmass/thet,epsi,as,ac,dy,dz,ym,cz,zm,sigma,rzero
 754 c      real eng(14),ymi(14),zmi(14)
 755 c      pi=3.14159
 756 c      write(ifch,*)'parameters of mass formula:'
 757 c      write(ifch,*)'theta=',thet,'   epsi=',epsi
 758 c      write(ifch,*)'as=',as,'   ac=',ac
 759 c      write(ifch,*)'dy=',dy,'   dz=',dz
 760 c      write(ifch,*)'ym=',ym
 761 c      write(ifch,*)'cz dz zm=',cz,dz,zm
 762 c      write(ifch,*)'sigma**1/3=',sigma**(1./3.),'   rzero=',rzero
 763 c      write(ifch,*)'mass:'
 764 c      write(ifch,5000)(j,j=1,14)
 765 c5000  format(5x,'a:',14i5)
 766 c      do 4 j=1,14
 767 c      a=j
 768 c      ymi(j)=ym*a
 769 c4     zmi(j)=cz/(dz/a+zm/a**(1./3.))
 770 c      write(ifch,5002)(ymi(j),j=1,14)
 771 c5002  format(1x,'ymin: ',14f5.2)
 772 c      write(ifch,5003)(zmi(j),j=1,14)
 773 c5003  format(1x,'zmin: ',14f5.2)
 774 c      do 2 i=1,15
 775 c      ns=11-i
 776 c      do 3 j=1,14
 777 c      a=j
 778 c      y=a-ns
 779 c      z=0.
 780 c3     eng(j)=smass(a,y,z)/a
 781 c      write(ifch,5001)ns,(eng(j),j=1,14)
 782 c5001  format(1x,'s=',i2,2x,14f5.2)
 783 c2     continue
 784 c      write(ifch,*)'mass-mass(free):'
 785 c      write(ifch,5000)(j,j=1,14)
 786 c      do 5 i=1,15
 787 c      ns=11-i
 788 c      do 6 j=1,14
 789 c      a=j
 790 c      y=a-ns
 791 c      z=0.
 792 c      call smassu(a,y,z,ku,kd,ks,kc)
 793 c6     eng(j)=(smass(a,y,z)-utamnu(ku,kd,ks,kc,0,0,3))/a
 794 c      write(ifch,5001)ns,(eng(j),j=1,14)
 795 c5     continue
 796 c
 797 c      stop
 798 c      end
 799 c
 800 cc-----------------------------------------------------------------------
 801 c      subroutine smasst(kux,kdx,ksx,kcx,a,y,z)
 802 cc-----------------------------------------------------------------------
 803 cc     input: kux,kdx,ksx,kcx = net quark numbers (for u,d,s,c quarks).
 804 cc     output: massnr a, hypercharge y and charge z.
 805 cc-----------------------------------------------------------------------
 806 c      sg=1
 807 c      if(kux+kdx+ksx+kcx.lt.0.)sg=-1
 808 c      ku=sg*kux
 809 c      kd=sg*kdx
 810 c      ks=sg*ksx
 811 c      kc=sg*kcx
 812 c      k=ku+kd+ks+kc
 813 c      if(mod(k,3).ne.0)stop'noninteger baryon number'
 814 c      a=k/3
 815 c      y=a-ks
 816 c      nz=2*ku-kd-ks+2*kc
 817 c      if(mod(nz,3).ne.0)stop'noninteger charge'
 818 c      z=nz/3
 819 c      return
 820 c      end
 821 c
 822 cc-----------------------------------------------------------------------
 823 c      subroutine smassu(ax,yx,zx,ku,kd,ks,kc)
 824 cc-----------------------------------------------------------------------
 825 cc     input: massnr ax, hypercharge yx and charge zx.
 826 cc     output: ku,kd,ks,kc = net quark numbers (for u,d,s,c quarks).
 827 cc-----------------------------------------------------------------------
 828 c      sg=1
 829 c      if(ax.lt.0.)sg=-1
 830 c      a=sg*ax
 831 c      y=sg*yx
 832 c      z=sg*zx
 833 c      ku=nint(a+z)
 834 c      kd=nint(a-z+y)
 835 c      ks=nint(a-y)
 836 c      kc=0
 837 c      return
 838 c      end
 839 c
 840 cc-----------------------------------------------------------------------
 841 c      function spoc(a,b,c,d,x)
 842 cc-----------------------------------------------------------------------
 843 cc     power fctn with cutoff
 844 cc-----------------------------------------------------------------------
 845 c      spoc=0
 846 c      if(a.eq.0..and.b.eq.0.)return
 847 c      spoc =a+b*x**c
 848 c      spoc0=a+b*d**c
 849 c      spoc=amin1(spoc,spoc0)
 850 c      spoc=amax1(0.,spoc)
 851 c      return
 852 c      end
 853 c
 854 c-----------------------------------------------------------------------
 855       function utacos(x)
 856 c-----------------------------------------------------------------------
 857 c     returns acos(x) for -1 <= x <= 1 , acos(+-1) else
 858 c-----------------------------------------------------------------------
 859       include 'epos.inc'
 860       argum=x
 861       if(x.lt.-1.)then
 862       if(ish.ge.1)then
 863       call utmsg('utacos')
 864       write(ifch,*)'*****  argum = ',argum,' set -1'
 865       call utmsgf
 866       endif
 867       argum=-1.
 868       elseif(x.gt.1.)then
 869       if(ish.ge.1)then
 870       call utmsg('utacos')
 871       write(ifch,*)'*****  argum = ',argum,' set 1'
 872       call utmsgf
 873       endif
 874       argum=1.
 875       endif
 876       utacos=acos(argum)
 877       return
 878       end
 879
 880 c----------------------------------------------------------------------
 881       function utamnu(keux,kedx,kesx,kecx,kebx,ketx,modus)
 882 c----------------------------------------------------------------------
 883 c     returns min mass of droplet with given u,d,s,c content
 884 c     keux: net u quark number
 885 c     kedx: net d quark number
 886 c     kesx: net s quark number
 887 c     kecx: net c quark number
 888 c     kebx: net b quark number
 889 c     ketx: net t quark number
 890 c     modus: 4=two lowest multiplets; 5=lowest multiplet
 891 c----------------------------------------------------------------------
 892       common/files/ifop,ifmt,ifch,ifcx,ifhi,ifdt,ifcp,ifdr
 893       common/csjcga/amnull,asuha(7)
 894       common/drop4/asuhax(7),asuhay(7)
 895
 896       if(modus.lt.4.or.modus.gt.5)stop'UTAMNU: not supported'
 897  1    format(' flavours:',6i5 )
 898  100  format(' flavours+mass:',6i5,f8.2 )
 899 c      write(ifch,1)keux,kedx,kesx,kecx,kebx,ketx      c          write(ifmt,*)'wrong mass in gethad ',damss
 900
 901       amnull=0.
 902
 903       do i=1,7
 904       if(modus.eq.4)asuha(i)=asuhax(i)    !two lowest multiplets
 905       if(modus.eq.5)asuha(i)=asuhay(i)    !lowest multiplet
 906       enddo
 907
 908       ke=iabs(keux+kedx+kesx+kecx+kebx+ketx)
 909
 910       if(keux+kedx+kesx+kecx+kebx+ketx.ge.0)then
 911       keu=keux
 912       ked=kedx
 913       kes=kesx
 914       kec=kecx
 915       keb=kebx
 916       ket=ketx
 917       else
 918       keu=-keux
 919       ked=-kedx
 920       kes=-kesx
 921       kec=-kecx
 922       keb=-kebx
 923       ket=-ketx
 924       endif
 925
 926 c      write(ifch,*)keu,ked,kes,kec,keb,ket
 927
 928 c   removing top mesons  to remove t quarks or antiquarks
 929       if(ket.ne.0)then
 930 12    continue
 931       ii=sign(1,ket)
 932       ket=ket-ii
 933       if(ii*keu.le.ii*ked)then
 934       keu=keu+ii
 935       else
 936       ked=ked+ii
 937       endif
 938       amnull=amnull+200.    ! ???????
 939       if(ket.ne.0)goto12
 940       endif
 941
 942 c   removing bottom mesons  to remove b quarks or antiquarks
 943       if(keb.ne.0)then
 944 11    continue
 945       ii=sign(1,keb)
 946       keb=keb-ii
 947       if(ii*keu.le.ii*ked)then
 948       keu=keu+ii
 949       else
 950       ked=ked+ii
 951       endif
 952       amnull=amnull+5.28   ! (B-meson)
 953       if(keb.ne.0)goto11
 954       endif
 955
 956 c   removing charm mesons  to remove c quarks or antiquarks
 957       if(kec.ne.0)then
 958 10    continue
 959       ii=sign(1,kec)
 960       kec=kec-ii
 961       if(keu*ii.le.ked*ii)then
 962       keu=keu+ii
 963       else
 964       ked=ked+ii
 965       endif
 966       amnull=amnull+1.87  ! (D-meson)
 967       if(kec.ne.0)goto10
 968       endif
 969
 970 c      write(ifch,100)keu,ked,kes,kec,keb,ket,amnull
 971
 972 c   removing mesons to remove s antiquarks
 973 5     continue
 974       if(kes.lt.0)then
 975       amnull=amnull+asuha(6)
 976       if(keu.ge.ked)then
 977       keu=keu-1
 978       else
 979       ked=ked-1
 980       endif
 981       kes=kes+1
 982       goto5
 983       endif
 984
 985 c   removing mesons to remove d antiquarks
 986 6     continue
 987       if(ked.lt.0)then
 988       if(keu.ge.kes)then
 989       amnull=amnull+asuha(5)
 990       keu=keu-1
 991       else
 992       amnull=amnull+asuha(6)
 993       kes=kes-1
 994       endif
 995       ked=ked+1
 996       goto6
 997       endif
 998
 999 c   removing mesons to remove u antiquarks
1000 7     continue
1001       if(keu.lt.0)then
1002       if(ked.ge.kes)then
1003       amnull=amnull+asuha(5)
1004       ked=ked-1
1005       else
1006       amnull=amnull+asuha(6)
1007       kes=kes-1
1008       endif
1009       keu=keu+1
1010       goto7
1011       endif
1012
1013 c      write(ifch,100)keu,ked,kes,kec,keb,ket,amnull
1014 c      print*,keu,ked,kes,kec,keb,ket,amnull
1015
1016       if(keu+ked+kes+kec+keb+ket.ne.ke)
1017      *call utstop('utamnu: sum_kei /= ke&')
1018       keq=keu+ked
1019       keqx=keq
1020       amnux=0
1021
1022 c   removing strange baryons
1023       i=4
1024 2     i=i-1
1025 3     continue
1026       if((4-i)*kes.gt.(i-1)*keq)then
1027       amnux=amnux+asuha(1+i)
1028       kes=kes-i
1029       keq=keq-3+i
1030       if(kes.lt.0)call utstop('utamnu: negative kes&')
1031       if(keq.lt.0)call utstop('utamnu: negative keq&')
1032       goto3
1033       endif
1034       if(i.gt.1)goto2
1035       if(keqx.gt.keq)then
1036       do 8 k=1,keqx-keq
1037       if(keu.ge.ked)then
1038       keu=keu-1
1039       else
1040       ked=ked-1
1041       endif
1042 8     continue
1043       endif
1044
1045       if(keu+ked.ne.keq)call utstop('utamnu: keu+ked /= keq&')
1046 c      write(ifch,100)keu,ked,kes,kec,keb,ket,amnull+amnux
1047 c      print*,keu,ked,kes,kec,keb,ket,amnull+amnux
1048
1049 c   removing nonstrange baryons
1050 9     continue
1051       if(keu.gt.2*ked)then
1052       amnux=amnux+asuha(7)
1053       keu=keu-3
1054       if(keu.lt.0)call utstop('utamnu: negative keu&')
1055       goto9
1056       endif
1057       if(ked.gt.2*keu)then
1058       amnux=amnux+asuha(7)
1059       ked=ked-3
1060       if(ked.lt.0)call utstop('utamnu: negative ked&')
1061       goto9
1062       endif
1063       keq=keu+ked
1064
1065 c      write(ifch,100)keu,ked,kes,kec,keb,ket,amnull+amnux
1066 c      print*,keu,ked,kes,kec,keb,ket,amnull+amnux
1067
1068       if(mod(keq,3).ne.0)call utstop('utamnu: mod(keq,3) /= 0&')
1069       amnux=amnux+asuha(1)*keq/3
1070
1071 c      write(ifch,100)keu,ked,kes,kec,keb,ket,amnull+amnux
1072 c      print*,keu,ked,kes,kec,keb,ket,amnull+amnux
1073
1074       amnull=amnull+amnux
1075
1076       if(amnull.eq.0)amnull=asuha(5)
1077
1078 1000  utamnu=amnull
1079       return
1080       end
1081
1082 c-----------------------------------------------------------------------
1083       function utamnx(jcp,jcm)
1084 c-----------------------------------------------------------------------
1085 c returns minimum mass for the decay of jcp---jcm (by calling utamnu).
1086 c-----------------------------------------------------------------------
1087       parameter (nflav=6)
1088       integer jcp(nflav,2),jcm(nflav,2)
1089
1090       do i=1,nflav
1091       do j=1,2
1092       if(jcp(i,j).ne.0)goto1
1093       enddo
1094       enddo
1095       keu=jcm(1,1)-jcm(1,2)
1096       ked=jcm(2,1)-jcm(2,2)
1097       kes=jcm(3,1)-jcm(3,2)
1098       kec=jcm(4,1)-jcm(4,2)
1099       keb=jcm(5,1)-jcm(5,2)
1100       ket=jcm(6,1)-jcm(6,2)
1101       utamnx=utamnu(keu,ked,kes,kec,keb,ket,5)
1102       return
1103 1     continue
1104
1105       do i=1,nflav
1106       do j=1,2
1107       if(jcm(i,j).ne.0)goto2
1108       enddo
1109       enddo
1110       keu=jcp(1,1)-jcp(1,2)
1111       ked=jcp(2,1)-jcp(2,2)
1112       kes=jcp(3,1)-jcp(3,2)
1113       kec=jcp(4,1)-jcp(4,2)
1114       keb=jcp(5,1)-jcp(5,2)
1115       ket=jcp(6,1)-jcp(6,2)
1116       utamnx=utamnu(keu,ked,kes,kec,keb,ket,5)
1117       return
1118 2     continue
1119
1120       keu=jcp(1,1)-jcp(1,2)
1121       ked=jcp(2,1)-jcp(2,2)
1122       kes=jcp(3,1)-jcp(3,2)
1123       kec=jcp(4,1)-jcp(4,2)
1124       keb=jcp(5,1)-jcp(5,2)
1125       ket=jcp(6,1)-jcp(6,2)
1126       ke=keu+ked+kes+kec+keb+ket
1127       if(mod(ke+1,3).eq.0)then
1128         keu=keu+1
1129         amms1=utamnu(keu,ked,kes,kec,keb,ket,5)
1130         keu=keu-1
1131         ked=ked+1
1132         amms2=utamnu(keu,ked,kes,kec,keb,ket,5)
1133       elseif(mod(ke-1,3).eq.0)then
1134         keu=keu-1
1135         amms1=utamnu(keu,ked,kes,kec,keb,ket,5)
1136         keu=keu+1
1137         ked=ked-1
1138         amms2=utamnu(keu,ked,kes,kec,keb,ket,5)
1139       else
1140         call utstop('utamnx: no singlet possible (1)&')
1141       endif
1142       keu=jcm(1,1)-jcm(1,2)
1143       ked=jcm(2,1)-jcm(2,2)
1144       kes=jcm(3,1)-jcm(3,2)
1145       kec=jcm(4,1)-jcm(4,2)
1146       keb=jcm(5,1)-jcm(5,2)
1147       ket=jcm(6,1)-jcm(6,2)
1148       ke=keu+ked+kes+kec+keb+ket
1149       if(mod(ke+1,3).eq.0)then
1150         keu=keu+1
1151         amms3=utamnu(keu,ked,kes,kec,keb,ket,5)
1152         keu=keu-1
1153         ked=ked+1
1154         amms4=utamnu(keu,ked,kes,kec,keb,ket,5)
1155       elseif(mod(ke-1,3).eq.0)then
1156         keu=keu-1
1157         amms3=utamnu(keu,ked,kes,kec,keb,ket,5)
1158         keu=keu+1
1159         ked=ked-1
1160         amms4=utamnu(keu,ked,kes,kec,keb,ket,5)
1161       else
1162         call utstop('utamnx: no singlet possible (2)&')
1163       endif
1164       utamnx=min(amms1+amms3,amms2+amms4)
1165 c       print *,amms1,amms3,amms2,amms4,jcp,jcm
1166       return
1167       end
1168
1169
1170
1171 cc-----------------------------------------------------------------------
1172 c      function utamny(jcp,jcm)
1173 cc-----------------------------------------------------------------------
1174 cc returns minimum mass of jcp+jcm (by calling utamnu).
1175 cc-----------------------------------------------------------------------
1176 c      parameter (nflav=6)
1177 c      integer jcp(nflav,2),jcm(nflav,2),jc(nflav,2)
1178 c      do 7 nf=1,nflav
1179 c      jc(nf,1)=jcp(nf,1)+jcm(nf,1)
1180 c7     jc(nf,2)=jcp(nf,2)+jcm(nf,2)
1181 c      keu=jc(1,1)-jc(1,2)
1182 c      ked=jc(2,1)-jc(2,2)
1183 c      kes=jc(3,1)-jc(3,2)
1184 c      kec=jc(4,1)-jc(4,2)
1185 c      keb=jc(5,1)-jc(5,2)
1186 c      ket=jc(6,1)-jc(6,2)
1187 c      utamny=utamnu(keu,ked,kes,kec,keb,ket,5)
1188 c      return
1189 c      end
1190 c
1191 c-----------------------------------------------------------------------
1192       function utamnz(jc,modus)
1193 c-----------------------------------------------------------------------
1194 c returns minimum mass of jc (by calling utamnu).
1195 c-----------------------------------------------------------------------
1196       parameter (nflav=6)
1197       integer jc(nflav,2)
1198       keu=jc(1,1)-jc(1,2)
1199       ked=jc(2,1)-jc(2,2)
1200       kes=jc(3,1)-jc(3,2)
1201       kec=jc(4,1)-jc(4,2)
1202       keb=jc(5,1)-jc(5,2)
1203       ket=jc(6,1)-jc(6,2)
1204       utamnz=utamnu(keu,ked,kes,kec,keb,ket,modus)
1205       return
1206       end
1207
1208 c-----------------------------------------------------------------------
1209       subroutine utar(i1,i2,i3,x0,x1,x2,x3,xx)
1210 c-----------------------------------------------------------------------
1211 c     returns the array xx with xx(1)=x0 <= xx(i) <= xx(i3)=x3
1212 c-----------------------------------------------------------------------
1213       real xx(i3)
1214       do 1 i=1,i1-1
1215   1   xx(i)=x0+(i-1.)/(i1-1.)*(x1-x0)
1216       do 2 i=i1,i2-1
1217   2   xx(i)=x1+(i-i1*1.)/(i2-i1*1.)*(x2-x1)
1218       do 3 i=i2,i3
1219   3   xx(i)=x2+(i-i2*1.)/(i3-i2*1.)*(x3-x2)
1220       return
1221       end
1222
1223 cc---------------------------------------------------------------------
1224 c      subroutine utaxis(i,j,a1,a2,a3)
1225 cc-----------------------------------------------------------------------
1226 cc     calculates the axis defined by the ptls i,j in the i,j cm system
1227 cc---------------------------------------------------------------------
1228 c      include 'epos.inc'
1229 c      double precision pi1,pi2,pi3,pi4,pj1,pj2,pj3,pj4,p1,p2,p3,p4,p5
1230 c     *,err,a
1231 c      a1=0
1232 c      a2=0
1233 c      a3=1
1234 c      pi1=dble(pptl(1,i))
1235 c      pi2=dble(pptl(2,i))
1236 c      pi3=dble(pptl(3,i))
1237 c      pi4=dble(pptl(4,i))
1238 c      pj1=dble(pptl(1,j))
1239 c      pj2=dble(pptl(2,j))
1240 c      pj3=dble(pptl(3,j))
1241 c      pj4=dble(pptl(4,j))
1242 c      p1=pi1+pj1
1243 c      p2=pi2+pj2
1244 c      p3=pi3+pj3
1245 c      p4=pi4+pj4
1246 c      p5=dsqrt(p4**2-p3**2-p2**2-p1**2)
1247 c      call utlob2(1,p1,p2,p3,p4,p5,pi1,pi2,pi3,pi4,50)
1248 c      call utlob2(1,p1,p2,p3,p4,p5,pj1,pj2,pj3,pj4,51)
1249 c           err=(pi1+pj1)**2+(pi2+pj2)**2+(pi3+pj3)**2
1250 c           if(err.gt.1d-3)then
1251 c      call utmsg('utaxis')
1252 c      write(ifch,*)'*****  err=',err
1253 c      write(ifch,*)'pi:',pi1,pi2,pi3,pi4
1254 c      write(ifch,*)'pj:',pj1,pj2,pj3,pj4
1255 c      call utmsgf
1256 c           endif
1257 c      a=dsqrt( (pj1-pi1)**2 + (pj2-pi2)**2 + (pj3-pi3)**2 )
1258 c      if(a.eq.0.d0)return
1259 c      a1=sngl((pi1-pj1)/a)
1260 c      a2=sngl((pi2-pj2)/a)
1261 c      a3=sngl((pi3-pj3)/a)
1262 c      return
1263 c      end
1264 c
1265 cc-----------------------------------------------------------------------
1266 c      subroutine utchm(arp,arm,ii)
1267 cc-----------------------------------------------------------------------
1268 cc     checks whether arp**2=0 and arm**2=0.
1269 cc-----------------------------------------------------------------------
1270 c      include 'epos.inc'
1271 c      double precision arp(4),arm(4),difp,difm
1272 c      difp=arp(4)**2-arp(1)**2-arp(2)**2-arp(3)**2
1273 c      difm=arm(4)**2-arm(1)**2-arm(2)**2-arm(3)**2
1274 c      if(dabs(difp).gt.1e-3*arp(4)**2
1275 c     *.or.dabs(difm).gt.1e-3*arm(4)**2)then
1276 c      call utmsg('utchm ')
1277 c      write(ifch,*)'*****  mass non zero  -  ',ii
1278 c      write(ifch,*)'jet-mass**2`s:    ',difp,difm
1279 c      write(ifch,*)'energy**2`s:      ',arp(4)**2,arm(4)**2
1280 c      write(ifch,*)(sngl(arp(i)),i=1,4)
1281 c      write(ifch,*)(sngl(arm(i)),i=1,4)
1282 c      call utmsgf
1283 c      endif
1284 c      return
1285 c      end
1286 c
1287 c-----------------------------------------------------------------------
1288       subroutine utclea(nptlii,nptl0)
1289 c-----------------------------------------------------------------------
1290 c     starting from nptlii
1291 c     overwrites istptl=99 particles in /cptl/, reduces so nptl
1292 c     and update minfra and maxfra
1293 c-----------------------------------------------------------------------
1294       include 'epos.inc'
1295       integer newptl(mxptl)!,oldptl(mxptl),ii(mxptl)
1296
1297       ish0=ish
1298       if(ishsub/100.eq.18)ish=mod(ishsub,100)
1299
1300       call utpri('utclea',ish,ishini,2)
1301
1302       nptli=max(maproj+matarg+1,nptlii)
1303       minfra0=minfra
1304       maxfra0=maxfra
1305       minfra1=maxfra
1306       maxfra1=minfra
1307       if(ish.ge.2)write(ifch,*)'entering subr utclea:',nptl
1308      &                                                ,minfra,maxfra
1309       if(ish.ge.7)then
1310       write(ifch,*)('-',l=1,68)
1311       write(ifch,*)'sr utclea. initial.'
1312       write(ifch,*)('-',l=1,68)
1313       do 34 n=nptli,nptl
1314       write(ifch,116)iorptl(n),jorptl(n),n,ifrptl(1,n),ifrptl(2,n)
1315      *,idptl(n),sqrt(pptl(1,n)**2+pptl(2,n)**2),pptl(3,n),pptl(5,n)
1316      *,istptl(n),ityptl(n)
1317 34    continue
1318 116   format(1x,i6,i6,4x,i6,4x,i6,i6,i12,3x,3(e8.2,1x),i3,i3)
1319       endif
1320
1321 c      ish=ish0
1322 c      ish0=ish
1323 c      if(ishsub/100.eq.18)ish=mod(ishsub,100)
1324
1325       i=nptli-1
1326 1     i=i+1
1327       if(i.gt.nptl)goto 1000
1328       if(istptl(i).eq.99)goto 2
1329       newptl(i)=i
1330 c      oldptl(i)=i
1331       goto 1
1332
1333 2     i=i-1
1334       j=i
1335 3     i=i+1
1336 4     j=j+1
1337       if(j.gt.nptl)goto 5
1338       newptl(j)=0
1339       if(istptl(j).eq.99)goto 4
1340       newptl(j)=i
1341 c      oldptl(i)=j
1342 c      write(ifch,*)'move',j,' to ',i
1343 c       write(ifch,*)idptl(i),ityptl(i),idptl(j),ityptl(j),minfra,maxfra
1344       call utrepl(i,j)
1345       if(j.ge.minfra0.and.j.le.maxfra0)then
1346         minfra1=min(minfra1,i)
1347         maxfra1=max(maxfra1,i)
1348       endif
1349       goto 3
1350
1351 5     nptl=i-1
1352       if(nptl.eq.0)then
1353         nptl0=0
1354         goto 1000
1355       endif
1356
1357 20    n0=newptl(nptl0)
1358       if(n0.gt.0)then
1359       nptl0=n0
1360       else
1361       nptl0=nptl0-1
1362       if(nptl0.gt.0)goto 20
1363       endif
1364
1365
1366 c      do 11 k=1,nptl
1367 c      io=iorptl(k)
1368 c      if(io.le.0)ii(k)=io
1369 c      if(io.gt.0)ii(k)=newptl(io)
1370 c11    continue
1371 c      do 12 k=1,nptl
1372 c12    iorptl(k)=ii(k)
1373 c
1374 c      do 13 k=1,nptl
1375 c      jo=jorptl(k)
1376 c      if(jo.le.0)ii(k)=jo
1377 c      if(jo.gt.0)ii(k)=newptl(jo)
1378 c13    continue
1379 c      do 14 k=1,nptl
1380 c14    jorptl(k)=ii(k)
1381 c
1382 c      do 15 k=1,nptl
1383 c      if1=ifrptl(1,k)
1384 c      if(if1.le.0)ii(k)=if1
1385 c      if(if1.gt.0)ii(k)=newptl(if1)
1386 c15    continue
1387 c      do 16 k=1,nptl
1388 c16    ifrptl(1,k)=ii(k)
1389 c
1390 c      do 17 k=1,nptl
1391 c      if2=ifrptl(2,k)
1392 c      if(if2.le.0)ii(k)=if2
1393 c      if(if2.gt.0)ii(k)=newptl(if2)
1394 c17    continue
1395 c      do 18 k=1,nptl
1396 c18    ifrptl(2,k)=ii(k)
1397 c
1398 c      do 19 k=1,nptl
1399 c      if(ifrptl(1,k).eq.0.and.ifrptl(2,k).gt.0)ifrptl(1,k)=ifrptl(2,k)
1400 c      if(ifrptl(2,k).eq.0.and.ifrptl(1,k).gt.0)ifrptl(2,k)=ifrptl(1,k)
1401 c19    continue
1402
1403 1000  continue
1404
1405       if(minfra1.lt.minfra0)minfra=minfra1
1406       if(maxfra1.ge.minfra1)maxfra=maxfra1
1407
1408       if(ish.ge.2)then
1409       write(ifch,*)'before exiting subr utclea:'
1410       do 35 n=1,nptl
1411       write(ifch,116)iorptl(n),jorptl(n),n,ifrptl(1,n),ifrptl(2,n)
1412      *,idptl(n),sqrt(pptl(1,n)**2+pptl(2,n)**2),pptl(3,n),pptl(5,n)
1413      *,istptl(n),ityptl(n)
1414 35    continue
1415       endif
1416
1417       if(ish.ge.2)write(ifch,*)'exiting subr utclea:',nptl
1418      &                                                ,minfra,maxfra
1419
1420       call utprix('utclea',ish,ishini,2)
1421       ish=ish0
1422       return
1423       end
1424
1425 c---------------------------------------------------------------------
1426       subroutine utfit(x,y,ndata,sig,mwt,a,b,siga,sigb,chi2,q)
1427 c---------------------------------------------------------------------
1428 c linear fit to data
1429 c input:
1430 c    ndata: nr of data points
1431 c    x(),y(),sig(): data
1432 c    mwt: unweighted (0) or weighted (else) data points
1433 c output:
1434 c    a,b: parameters of linear fit a+b*x
1435 c---------------------------------------------------------------------
1436       INTEGER mwt,ndata
1437       REAL a,b,chi2,q,siga,sigb,sig(ndata),x(ndata),y(ndata)
1438 CU    USES utgmq
1439       INTEGER i
1440       REAL sigdat,ss,st2,sx,sxoss,sy,t,wt,utgmq
1441       sx=0.
1442       sy=0.
1443       st2=0.
1444       b=0.
1445       if(mwt.ne.0) then
1446         ss=0.
1447         do 11 i=1,ndata
1448           wt=1./(sig(i)**2)
1449           ss=ss+wt
1450           sx=sx+x(i)*wt
1451           sy=sy+y(i)*wt
1452 11      continue
1453       else
1454         do 12 i=1,ndata
1455           sx=sx+x(i)
1456           sy=sy+y(i)
1457 12      continue
1458         ss=float(ndata)
1459       endif
1460       sxoss=sx/ss
1461       if(mwt.ne.0) then
1462         do 13 i=1,ndata
1463           t=(x(i)-sxoss)/sig(i)
1464           st2=st2+t*t
1465           b=b+t*y(i)/sig(i)
1466 13      continue
1467       else
1468         do 14 i=1,ndata
1469           t=x(i)-sxoss
1470           st2=st2+t*t
1471           b=b+t*y(i)
1472 14      continue
1473       endif
1474       b=b/st2
1475       a=(sy-sx*b)/ss
1476       siga=sqrt((1.+sx*sx/(ss*st2))/ss)
1477       sigb=sqrt(1./st2)
1478       chi2=0.
1479       if(mwt.eq.0) then
1480         do 15 i=1,ndata
1481           chi2=chi2+(y(i)-a-b*x(i))**2
1482 15      continue
1483         q=1.
1484         sigdat=sqrt(chi2/(ndata-2))
1485         siga=siga*sigdat
1486         sigb=sigb*sigdat
1487       else
1488         do 16 i=1,ndata
1489           chi2=chi2+((y(i)-a-b*x(i))/sig(i))**2
1490 16      continue
1491         q=utgmq(0.5*(ndata-2),0.5*chi2)
1492       endif
1493       return
1494       END
1495
1496 c-----------------------------------------------------------------------
1497       function utgam1(x)
1498 c-----------------------------------------------------------------------
1499 c  gamma fctn tabulated
1500 c  single precision
1501 c-----------------------------------------------------------------------
1502       double precision utgamtab,utgam,al,dl
1503       common/gamtab/utgamtab(10000)
1504
1505       if(x.gt.0.01.and.x.lt.99.99)then
1506         al=100.d0*dble(x)
1507         k1=int(al)
1508         k2=k1+1
1509         dl =al-dble(k1)
1510         utgam1=real(utgamtab(k2)*dl+utgamtab(k1)*(1.d0-dl))
1511       elseif(x.eq.0.)then
1512         utgam1=0.
1513       else
1514         utgam1=real(utgam(dble(x)))
1515       endif
1516
1517       end
1518
1519 c-----------------------------------------------------------------------
1520       double precision function utgam2(x)
1521 c-----------------------------------------------------------------------
1522 c  gamma fctn tabulated
1523 c  double precision
1524 c-----------------------------------------------------------------------
1525       double precision utgamtab,x,al,dl,utgam
1526       common/gamtab/utgamtab(10000)
1527
1528       if(x.gt.0.01d0.and.x.le.99.99d0)then
1529         al=100.d0*x
1530         k1=int(al)
1531         k2=k1+1
1532         dl =al-dble(k1)
1533         utgam2=utgamtab(k2)*dl+utgamtab(k1)*(1.d0-dl)
1534       elseif(x.eq.0.d0)then
1535         utgam2=0.d0
1536       else
1537         utgam2=utgam(x)
1538       endif
1539
1540       end
1541
1542 c-----------------------------------------------------------------------
1543       double precision function utgam(x)
1544 c-----------------------------------------------------------------------
1545 c  gamma fctn
1546 c  double precision
1547 c-----------------------------------------------------------------------
1548       include 'epos.inc'
1549       double precision c(13),x,z,f
1550       data c
1551      1/ 0.00053 96989 58808, 0.00261 93072 82746, 0.02044 96308 23590,
1552      2  0.07309 48364 14370, 0.27964 36915 78538, 0.55338 76923 85769,
1553      3  0.99999 99999 99998,-0.00083 27247 08684, 0.00469 86580 79622,
1554      4  0.02252 38347 47260,-0.17044 79328 74746,-0.05681 03350 86194,
1555      5  1.13060 33572 86556/
1556       utgam=0d0
1557       z=x
1558       if(x .gt. 170.d0) goto6
1559       if(x .gt. 0.0d0) goto1
1560       if(x .eq. int(x)) goto5
1561       z=1.0d0-z
1562     1 f=1.0d0/z
1563       if(z .le. 1.0d0) goto4
1564       f=1.0d0
1565     2 continue
1566       if(z .lt. 2.0d0) goto3
1567       z=z-1.0d0
1568       f=f*z
1569       goto2
1570     3 z=z-1.0d0
1571     4 utgam=
1572      1 f*((((((c(1)*z+c(2))*z+c(3))*z+c(4))*z+c(5))*z+c(6))*z+c(7))/
1573      2 ((((((c(8)*z+c(9))*z+c(10))*z+c(11))*z+c(12))*z+c(13))*z+1.0d0)
1574       if(x .gt. 0.0d0) return
1575       utgam=3.141592653589793d0/(sin(3.141592653589793d0*x)*utgam)
1576       return
1577     5 write(ifch,10)sngl(x)
1578    10 format(1x,'argument of gamma fctn = ',e20.5)
1579       call utstop('utgam : negative integer argument&')
1580     6 write(ifch,11)sngl(x)
1581    11 format(1x,'argument of gamma fctn = ',e20.5)
1582       call utstop('utgam : argument too large&')
1583       end
1584
1585 c---------------------------------------------------------------------
1586       subroutine utgcf(gammcf,a,x,gln)
1587 c---------------------------------------------------------------------
1588       INTEGER ITMAX
1589       REAL a,gammcf,gln,x,EPS,FPMIN
1590       PARAMETER (ITMAX=100,EPS=3.e-7,FPMIN=1.e-30)
1591 CU    USES utgmln
1592       INTEGER i
1593       REAL an,b,c,d,del,h,utgmln
1594       gln=utgmln(a)
1595       b=x+1.-a
1596       c=1./FPMIN
1597       d=1./b
1598       h=d
1599       do 11 i=1,ITMAX
1600         an=-i*(i-a)
1601         b=b+2.
1602         d=an*d+b
1603         if(abs(d).lt.FPMIN)d=FPMIN
1604         c=b+an/c
1605         if(abs(c).lt.FPMIN)c=FPMIN
1606         d=1./d
1607         del=d*c
1608         h=h*del
1609         if(abs(del-1.).lt.EPS)goto 1
1610 11    continue
1611       pause 'a too large, ITMAX too small in utgcf'
1612 1     gammcf=exp(-x+a*log(x)-gln)*h
1613       return
1614       END
1615
1616 c---------------------------------------------------------------------
1617       function utgmln(xx)
1618 c---------------------------------------------------------------------
1619       REAL utgmln,xx
1620       INTEGER j
1621       DOUBLE PRECISION ser,stp,tmp,x,y,cof(6)
1622       SAVE cof,stp
1623       DATA cof,stp/76.18009172947146d0,-86.50532032941677d0,
1624      *24.01409824083091d0,-1.231739572450155d0,.1208650973866179d-2,
1625      *-.5395239384953d-5,2.5066282746310005d0/
1626       x=xx
1627       y=x
1628       tmp=x+5.5d0
1629       tmp=(x+0.5d0)*log(tmp)-tmp
1630       ser=1.000000000190015d0
1631       do 11 j=1,6
1632         y=y+1.d0
1633         ser=ser+cof(j)/y
1634 11    continue
1635       utgmln=tmp+log(stp*ser/x)
1636       return
1637       END
1638
1639 c---------------------------------------------------------------------
1640       function utgmq(a,x)
1641 c---------------------------------------------------------------------
1642       REAL a,utgmq,x
1643 CU    USES utgcf,utgser
1644       REAL gammcf,gamser,gln
1645       if(x.lt.0..or.a.le.0.)pause 'bad arguments in utgmq'
1646       if(x.lt.a+1.)then
1647         call utgser(gamser,a,x,gln)
1648         utgmq=1.-gamser
1649       else
1650         call utgcf(gammcf,a,x,gln)
1651         utgmq=gammcf
1652       endif
1653       return
1654       END
1655
1656 c---------------------------------------------------------------------
1657       subroutine utgser(gamser,a,x,gln)
1658 c---------------------------------------------------------------------
1659       INTEGER ITMAX
1660       REAL a,gamser,gln,x,EPS
1661       PARAMETER (ITMAX=100,EPS=3.e-7)
1662 CU    USES utgmln
1663       INTEGER n
1664       REAL ap,del,sum,utgmln
1665       gln=utgmln(a)
1666       if(x.le.0.)then
1667         if(x.lt.0.)pause 'x < 0 in utgser'
1668         gamser=0.
1669         return
1670       endif
1671       ap=a
1672       sum=1./a
1673       del=sum
1674       do 11 n=1,ITMAX
1675         ap=ap+1.
1676         del=del*x/ap
1677         sum=sum+del
1678         if(abs(del).lt.abs(sum)*EPS)goto 1
1679 11    continue
1680       pause 'a too large, ITMAX too small in utgser'
1681 1     gamser=sum*exp(-x+a*log(x)-gln)
1682       return
1683       END
1684
1685 c-------------------------------------------------------------------------
1686       subroutine uticpl(ic,ifla,iqaq,iret)
1687 c-------------------------------------------------------------------------
1688 c  adds a quark (iqaq=1) or antiquark (iqaq=2) of flavour ifla
1689 c  to 2-id ic
1690 c-------------------------------------------------------------------------
1691       include 'epos.inc'
1692       integer jc(nflav,2),ic(2)
1693       iret=0
1694       if(ifla.eq.0)return
1695       call iddeco(ic,jc)
1696       if(ish.ge.8)write(ifch,'(2i8,12i3)')ic,jc
1697       jqaq=3-iqaq
1698       if(jc(ifla,jqaq).gt.0)then
1699       jc(ifla,jqaq)=jc(ifla,jqaq)-1
1700       else
1701       jc(ifla,iqaq)=jc(ifla,iqaq)+1
1702       endif
1703       call idcomj(jc)
1704       call idenco(jc,ic,ireten)
1705       if(ish.ge.8)write(ifch,'(2i8,12i3)')ic,jc
1706       if(ireten.eq.1)iret=1
1707       if(ic(1).eq.0.and.ic(2).eq.0.and.ireten.eq.0)then
1708       ic(1)=100000
1709       ic(2)=100000
1710       endif
1711       return
1712       end
1713
1714 cc-----------------------------------------------------------------------
1715 c      subroutine utindx(n,xar,x,i)
1716 cc-----------------------------------------------------------------------
1717 cc  input:  dimension n
1718 cc          array xar(n) with xar(i) > xar(i-1)
1719 cc          some number x between xar(1) and xar(n)
1720 cc  output: the index i such that x is between xar(i)  and xar(i+1)
1721 cc-----------------------------------------------------------------------
1722 c      include 'epos.inc'
1723 c      real xar(n)
1724 c           if(x.lt.xar(1))then
1725 c      if(ish.ge.5)then
1726 c      call utmsg('utindx')
1727 c      write(ifch,*)'*****  x=',x,' < xar(1)=',xar(1)
1728 c      call utmsgf
1729 c      endif
1730 c      i=1
1731 c      return
1732 c           elseif(x.gt.xar(n))then
1733 c      if(ish.ge.5)then
1734 c      call utmsg('utindx')
1735 c      write(ifch,*)'*****  x=',x,' > xar(n)=',xar(n)
1736 c      call utmsgf
1737 c      endif
1738 c      i=n
1739 c      return
1740 c           endif
1741 c      lu=1
1742 c      lo=n
1743 c1     lz=(lo+lu)/2
1744 c      if((xar(lu).le.x).and.(x.le.xar(lz)))then
1745 c      lo=lz
1746 c      elseif((xar(lz).lt.x).and.(x.le.xar(lo)))then
1747 c      lu=lz
1748 c      else
1749 c      call utstop('utindx: no interval found&')
1750 c      endif
1751 c      if((lo-lu).ge.2) goto1
1752 c      if(lo.le.lu)call utstop('utinvt: lo.le.lu&')
1753 c      i=lu
1754 c      return
1755 c      end
1756 c
1757 c-----------------------------------------------------------------------
1758       function utinvt(n,x,q,y)
1759 c-----------------------------------------------------------------------
1760 c     returns x with y=q(x)
1761 c-----------------------------------------------------------------------
1762       include 'epos.inc'
1763       real x(n),q(n)
1764       if(q(n).eq.0.)call utstop('utinvt: q(n)=0&')
1765            if(y.lt.0.)then
1766       if(ish.ge.1)then
1767       call utmsg('utinvt')
1768       write(ifch,*)'*****  y=',y,' < 0'
1769       call utmsgf
1770       endif
1771       y=0.
1772            elseif(y.gt.q(n))then
1773       if(ish.ge.1)then
1774       call utmsg('utinvt')
1775       write(ifch,*)'*****  y=',y,' > ',q(n)
1776       call utmsgf
1777       endif
1778       y=q(n)
1779            endif
1780       lu=1
1781       lo=n
1782 1     lz=(lo+lu)/2
1783       if((q(lu).le.y).and.(y.le.q(lz)))then
1784       lo=lz
1785       elseif((q(lz).lt.y).and.(y.le.q(lo)))then
1786       lu=lz
1787       else
1788       write(ifch,*)'q(1),y,q(n):',q(1),y,q(n)
1789       write(ifch,*)'lu,lz,lo:',lu,lz,lo
1790       write(ifch,*)'q(lu),q(lz),q(lo):',q(lu),q(lz),q(lo)
1791       call utstop('utinvt: no interval found&')
1792       endif
1793       if((lo-lu).ge.2) goto1
1794       if(lo.le.lu)call utstop('utinvt: lo.le.lu&')
1795       utinvt=x(lu)+(y-q(lu))*(x(lo)-x(lu))/(q(lo)-q(lu))
1796       return
1797       end
1798
1799 c-----------------------------------------------------------------------
1800       subroutine utlob2(isig,p1,p2,p3,p4,p5,x1,x2,x3,x4,idi)
1801 c-----------------------------------------------------------------------
1802 c  performs a lorentz boost, double prec.
1803 c  isig=+1 is to boost the four vector x1,x2,x3,x4 such as to obtain it
1804 c  in the frame specified by the 5-vector p1...p5 (5-vector=4-vector+mass).
1805 c  isig=-1: the other way round, that means,
1806 c  if the 4-vector x1...x4 is given in some frame characterized by
1807 c  p1...p5 with respect to to some lab-frame, utlob2 returns the 4-vector
1808 c  x1...x4  in the lab frame.
1809 c  idi is a call identifyer (integer) to identify the call in case of problem
1810 c-----------------------------------------------------------------------
1811       include 'epos.inc'
1812       double precision beta(4),z(4),p1,p2,p3,p4,p5,pp,bp,x1,x2,x3,x4
1813      *,xx0,x10,x20,x30,x40,x4x,x0123
1814            if(ish.ge.2)then
1815       if(ish.ge.9)then
1816       write(ifch,101)x1,x2,x3,x4,(x4-x3)*(x4+x3)-x2*x2-x1*x1
1817       write(ifch,301)p1,p2,p3,p4,p5,(p4-p3)*(p4+p3)-p2*p2-p1*p1
1818 101   format(' utlob2: x =  ',5e13.5)
1819 301   format('         p =  ',6e13.5)
1820       endif
1821       pp=(p4-p3)*(p4+p3)-p2*p2-p1*p1
1822       if(dabs(pp-p5*p5).gt.1e-3*p4*p4.and.dabs(pp-p5*p5).gt.1e-3)then
1823       call utmsg('utlob2')
1824       write(ifch,*)'*****  p**2 .ne. p5**2'
1825       write(ifch,*)'call identifyer:',idi
1826       write(ifch,*)'p**2,p5**2: ',pp,p5*p5
1827       write(ifch,*)'p: ',p1,p2,p3,p4,p5
1828       call utmsgf
1829       endif
1830       x10=x1
1831       x20=x2
1832       x30=x3
1833       x40=x4
1834            endif
1835       xx0=(x4-x3)*(x4+x3)-x2*x2-x1*x1
1836       if(p5.le.0.)then
1837       call utmsg('utlob2')
1838       write(ifch,*)'*****  p5 negative.'
1839       write(ifch,*)'call identifyer:',idi
1840       write(ifch,*)'p(5): ',p1,p2,p3,p4,p5
1841       write(ifmt,*)'call identifyer:',idi
1842       write(ifmt,*)'p(5): ',p1,p2,p3,p4,p5
1843       call utmsgf
1844       call utstop('utlob2: p5 negative.&')
1845       endif
1846       z(1)=x1
1847       z(2)=x2
1848       z(3)=x3
1849       z(4)=x4
1850       beta(1)=-p1/p5
1851       beta(2)=-p2/p5
1852       beta(3)=-p3/p5
1853       beta(4)= p4/p5
1854       bp=0.
1855       do 220 k=1,3
1856 220   bp=bp+z(k)*isig*beta(k)
1857       do 230 k=1,3
1858 230   z(k)=z(k)+isig*beta(k)*z(4)
1859      *+isig*beta(k)*bp/(beta(4)+1.)
1860       z(4)=beta(4)*z(4)+bp
1861       x1=z(1)
1862       x2=z(2)
1863       x3=z(3)
1864       x4=z(4)
1865       if(ish.ge.9)
1866      *write(ifch,101)x1,x2,x3,x4,(x4-x3)*(x4+x3)-x2*x2-x1*x1
1867       x4x=x4
1868       x0123=xx0+x1*x1+x2*x2+x3*x3
1869       if(x0123.gt.0.)then
1870       x4=sign( dsqrt(x0123) , x4x )
1871       else
1872       x4=0
1873       endif
1874       if(ish.ge.9)then
1875       write(ifch,101)x1,x2,x3,x4,(x4-x3)*(x4+x3)-x2*x2-x1*x1
1876       endif
1877            if(ish.ge.2)then
1878       if(ish.ge.9)write(ifch,*)'check x**2_ini -- x**2_fin'
1879       if(dabs(x4-x4x).gt.1d-2*dabs(x4).and.dabs(x4-x4x).gt.1d-2)then
1880       call utmsg('utlob2')
1881       write(ifch,*)'*****  x**2_ini .ne. x**2_fin.'
1882       write(ifch,*)'call identifyer:',idi
1883       write(ifch,*)'x1 x2 x3 x4 x**2 (initial/final/corrected):'
1884 102   format(5e13.5)
1885       write(ifch,102)x10,x20,x30,x40,(x40-x30)*(x40+x30)-x20*x20-x10*x10
1886       write(ifch,102)x1,x2,x3,x4x,(x4x-x3)*(x4x+x3)-x2*x2-x1*x1
1887       write(ifch,102)x1,x2,x3,x4,(x4-x3)*(x4+x3)-x2*x2-x1*x1
1888       call utmsgf
1889       endif
1890            endif
1891       if(ish.ge.9)write(ifch,*)'return from utlob2'
1892       return
1893       end
1894
1895 c-----------------------------------------------------------------------
1896       subroutine utlob3(isig,p1,p2,p3,p4,p5,x1,x2,x3,x4)
1897 c-----------------------------------------------------------------------
1898 c  performs a lorentz boost, double prec.
1899 c  but arguments are single precision
1900 c-----------------------------------------------------------------------
1901       double precision xx1,xx2,xx3,xx4
1902       xx1=dble(x1)
1903       xx2=dble(x2)
1904       xx3=dble(x3)
1905       xx4=dble(x4)
1906       call utlob2(isig
1907      *,dble(p1),dble(p2),dble(p3),dble(p4),dble(p5)
1908      *,xx1,xx2,xx3,xx4,52)
1909       x1=sngl(xx1)
1910       x2=sngl(xx2)
1911       x3=sngl(xx3)
1912       x4=sngl(xx4)
1913       return
1914       end
1915
1916 c-----------------------------------------------------------------------
1917       subroutine utlob5(yboost,x1,x2,x3,x4,x5)
1918 c-----------------------------------------------------------------------
1919       amt=sqrt(x5**2+x1**2+x2**2)
1920       y=sign(1.,x3)*alog((x4+abs(x3))/amt)
1921       y=y-yboost
1922       x4=amt*cosh(y)
1923       x3=amt*sinh(y)
1924       return
1925       end
1926
1927 c-----------------------------------------------------------------------
1928       subroutine utlob4(isig,pp1,pp2,pp3,pp4,pp5,x1,x2,x3,x4)
1929 c-----------------------------------------------------------------------
1930 c  performs a lorentz boost, double prec.
1931 c  but arguments are partly single precision
1932 c-----------------------------------------------------------------------
1933       double precision xx1,xx2,xx3,xx4,pp1,pp2,pp3,pp4,pp5
1934       xx1=dble(x1)
1935       xx2=dble(x2)
1936       xx3=dble(x3)
1937       xx4=dble(x4)
1938       call utlob2(isig,pp1,pp2,pp3,pp4,pp5,xx1,xx2,xx3,xx4,53)
1939       x1=sngl(xx1)
1940       x2=sngl(xx2)
1941       x3=sngl(xx3)
1942       x4=sngl(xx4)
1943       return
1944       end
1945
1946
1947 c-----------------------------------------------------------------------
1948       subroutine utlobo(isig,p1,p2,p3,p4,p5,x1,x2,x3,x4)
1949 c-----------------------------------------------------------------------
1950 c     performs a lorentz boost
1951 c-----------------------------------------------------------------------
1952       include 'epos.inc'
1953       real beta(4),z(4)
1954       if(p5.le.0.)then
1955       call utmsg('utlobo')
1956       write(ifch,*)'*****  mass <= 0.'
1957       write(ifch,*)'p(5): ',p1,p2,p3,p4,p5
1958       call utmsgf
1959       call utstop('utlobo: mass <= 0.&')
1960       endif
1961       z(1)=x1
1962       z(2)=x2
1963       z(3)=x3
1964       z(4)=x4
1965       beta(1)=-p1/p5
1966       beta(2)=-p2/p5
1967       beta(3)=-p3/p5
1968       beta(4)= p4/p5
1969       bp=0.
1970       do 220 k=1,3
1971 220   bp=bp+z(k)*isig*beta(k)
1972       do 230 k=1,3
1973 230   z(k)=z(k)+isig*beta(k)*z(4)
1974      *+isig*beta(k)*bp/(beta(4)+1.)
1975       z(4)=beta(4)*z(4)+bp
1976       x1=z(1)
1977       x2=z(2)
1978       x3=z(3)
1979       x4=z(4)
1980       return
1981       end
1982
1983 c-----------------------------------------------------------------------
1984       subroutine utloc(ar,n,a,l)
1985 c-----------------------------------------------------------------------
1986       real ar(n)
1987       do 1 i=1,n
1988       l=i-1
1989       if(a.lt.ar(i))return
1990 1     continue
1991       l=n
1992       return
1993       end
1994
1995 cc-----------------------------------------------------------------------
1996 c      subroutine utlow(cone)
1997 cc-----------------------------------------------------------------------
1998 c      character*1 cone
1999 c      if(cone.eq.'A')cone='a'
2000 c      if(cone.eq.'B')cone='b'
2001 c      if(cone.eq.'C')cone='c'
2002 c      if(cone.eq.'D')cone='d'
2003 c      if(cone.eq.'E')cone='e'
2004 c      if(cone.eq.'F')cone='f'
2005 c      if(cone.eq.'G')cone='g'
2006 c      if(cone.eq.'H')cone='h'
2007 c      if(cone.eq.'I')cone='i'
2008 c      if(cone.eq.'J')cone='j'
2009 c      if(cone.eq.'K')cone='k'
2010 c      if(cone.eq.'L')cone='l'
2011 c      if(cone.eq.'M')cone='m'
2012 c      if(cone.eq.'N')cone='n'
2013 c      if(cone.eq.'O')cone='o'
2014 c      if(cone.eq.'P')cone='p'
2015 c      if(cone.eq.'Q')cone='q'
2016 c      if(cone.eq.'R')cone='r'
2017 c      if(cone.eq.'S')cone='s'
2018 c      if(cone.eq.'T')cone='t'
2019 c      if(cone.eq.'U')cone='u'
2020 c      if(cone.eq.'V')cone='v'
2021 c      if(cone.eq.'W')cone='w'
2022 c      if(cone.eq.'X')cone='x'
2023 c      if(cone.eq.'Y')cone='y'
2024 c      if(cone.eq.'Z')cone='z'
2025 c      return
2026 c      end
2027 c
2028 cc-----------------------------------------------------------------------
2029 c      subroutine utlow3(cthree)
2030 cc-----------------------------------------------------------------------
2031 c      character cthree*3
2032 c      do 1 i=1,3
2033 c1     call utlow(cthree(i:i))
2034 c      return
2035 c      end
2036 c
2037 cc-----------------------------------------------------------------------
2038 c      subroutine utlow6(csix)
2039 cc-----------------------------------------------------------------------
2040 c      character csix*6
2041 c      do 1 i=1,6
2042 c1     call utlow(csix(i:i))
2043 c      return
2044 c      end
2045 c
2046 cc-----------------------------------------------------------------------
2047 c      function utmom(k,n,x,q)
2048 cc-----------------------------------------------------------------------
2049 cc     calculates kth moment for f(x) with q(i)=int[0,x(i)]f(z)dz
2050 cc-----------------------------------------------------------------------
2051 c      real x(n),q(n)
2052 c      if(n.lt.2)call utstop('utmom : dimension too small&')
2053 c      utmom=0
2054 c      do 1 i=2,n
2055 c1     utmom=utmom+((x(i)+x(i-1))/2)**k*(q(i)-q(i-1))
2056 c      utmom=utmom/q(n)
2057 c      return
2058 c      end
2059 c
2060 c-----------------------------------------------------------------------
2061       function utpcm(a,b,c)
2062 c-----------------------------------------------------------------------
2063 c     calculates cm momentum for a-->b+c
2064 c-----------------------------------------------------------------------
2065       val=(a*a-b*b-c*c)*(a*a-b*b-c*c)-(2.*b*c)*(2.*b*c)
2066       if(val.lt.0..and.val.gt.-1e-4)then
2067       utpcm=0
2068       return
2069       endif
2070       utpcm=sqrt(val)/(2.*a)
2071       return
2072       end
2073
2074 c-----------------------------------------------------------------------
2075       double precision function utpcmd(a,b,c)
2076 c-----------------------------------------------------------------------
2077 c     calculates cm momentum for a-->b+c
2078 c-----------------------------------------------------------------------
2079       double precision a,b,c,val
2080       val=(a*a-b*b-c*c)*(a*a-b*b-c*c)-(2.*b*c)*(2.*b*c)
2081       if(val.lt.0d0.and.val.gt.-1d-4)then
2082       utpcmd=0d0
2083       return
2084       elseif(val.lt.0d0)then
2085         call utstop('utpcmd: problem with masses&')
2086       endif
2087       utpcmd=sqrt(val)/(2.d0*a)
2088       return
2089       end
2090
2091 c-----------------------------------------------------------------------
2092       subroutine utpri(text,idmmy,ishini,ishx)
2093 c-----------------------------------------------------------------------
2094       include 'epos.inc'
2095       character*6 text
2096 c      double precision seedx                               !!!
2097       ishini=ish
2098       if(ishevt.ne.0.and.nrevt+1.ne.ishevt)return
2099       if(nrpri.gt.0)then
2100       do nr=1,nrpri
2101       if(subpri(nr)(1:6).eq.text)then
2102       ish=ishpri(nr)
2103       endif
2104       enddo
2105       endif
2106       if(ish.ge.ishx)then
2107         write(ifch,'(1x,43a)')
2108      *  ('-',i=1,10),' entry ',text,' ',('-',i=1,30)
2109 c       call ranfgt(seedx)                                   !!!
2110 c       if(ish.ge.ishx)write(ifch,*)'seed:',seedx            !!!
2111       endif
2112       return
2113       end
2114
2115 c-----------------------------------------------------------------------
2116       subroutine utprix(text,idmmy,ishini,ishx)
2117 c-----------------------------------------------------------------------
2118       include 'epos.inc'
2119       character*6 text
2120       if(ishevt.ne.0.and.nrevt+1.ne.ishevt)return
2121       if(ish.ge.ishx)write(ifch,'(1x,44a)')
2122      *('-',i=1,30),' exit ',text,' ',('-',i=1,11)
2123       ish=ishini
2124       return
2125       end
2126
2127 c-----------------------------------------------------------------------
2128       subroutine utprj(text,idmmy,ishini,ishx)
2129 c-----------------------------------------------------------------------
2130       include 'epos.inc'
2131       character*20 text
2132 c      double precision seedx                               !!!
2133       idx=index(text,' ')-1
2134       ishini=ish
2135       if(ishevt.ne.0.and.nrevt+1.ne.ishevt)return
2136       if(nrpri.gt.0)then
2137       do nr=1,nrpri
2138       if(subpri(nr)(1:idx).eq.text(1:idx))then
2139       ish=ishpri(nr)
2140       endif
2141       enddo
2142       endif
2143       if(ish.ge.ishx)then
2144         write(ifch,'(1x,43a)')
2145      *  ('-',i=1,10),' entry ',text(1:idx),' ',('-',i=1,30)
2146 c       call ranfgt(seedx)                                   !!!
2147 c       if(ish.ge.ishx)write(ifch,*)'seed:',seedx            !!!
2148       endif
2149       return
2150       end
2151
2152 c-----------------------------------------------------------------------
2153       subroutine utprjx(text,idmmy,ishini,ishx)
2154 c-----------------------------------------------------------------------
2155       include 'epos.inc'
2156       character*20 text
2157       idx=index(text,' ')-1
2158       if(ishevt.ne.0.and.nrevt+1.ne.ishevt)return
2159       if(ish.ge.ishx)write(ifch,'(1x,44a)')
2160      *('-',i=1,30),' exit ',text(1:idx),' ',('-',i=1,11)
2161       ish=ishini
2162       return
2163       end
2164
2165 c-----------------------------------------------------------------------
2166       function utquad(m,x,f,k)
2167 c-----------------------------------------------------------------------
2168 c     performs an integration according to simpson
2169 c-----------------------------------------------------------------------
2170       real x(m),f(m)
2171       utquad=0
2172       do 1 i=1,k-1
2173   1   utquad=utquad+(f(i)+f(i+1))/2*(x(i+1)-x(i))
2174       return
2175       end
2176
2177 c-----------------------------------------------------------------------
2178       subroutine utquaf(fu,n,x,q,wo,x0,x1,x2,x3)
2179 c-----------------------------------------------------------------------
2180 c     returns q(i) = integral [x(1)->x(i)] fu(x) dx
2181 c-----------------------------------------------------------------------
2182       include 'epos.inc'
2183       real x(n),q(n),wo(n)
2184       parameter (m=10)
2185       real xa(m),fa(m)
2186       external fu
2187       if(x1.lt.x0.or.x2.lt.x1.or.x3.lt.x2)then
2188       if(ish.ge.1)then
2189       call utmsg('utquaf')
2190       write(ifch,*)'   xi=',x0,x1,x2,x3
2191       call utmsgf
2192       endif
2193       endif
2194       call utar(n/3,n*2/3,n,x0,x1,x2,x3,x)
2195       q(1)=0
2196       do 2 i=2,n
2197       do 3 k=1,m
2198       z=x(i-1)+(k-1.)/(m-1.)*(x(i)-x(i-1))
2199       xa(k)=z
2200 3     fa(k)=fu(z)
2201       q(i)=q(i-1)+utquad(m,xa,fa,m)
2202 2     continue
2203       return
2204       end
2205
2206 c-----------------------------------------------------------------------
2207       subroutine utrepl(i,j)
2208 c-----------------------------------------------------------------------
2209 c     i is replaced by j in /cptl/
2210 c-----------------------------------------------------------------------
2211       include 'epos.inc'
2212       do 1 k=1,5
2213 1     pptl(k,i)  =pptl(k,j)
2214       iorptl(i)  = 0 !iorptl(j)
2215       idptl(i)   =idptl(j)
2216       istptl(i)  =istptl(j)
2217       do 2 k=1,2
2218 2     tivptl(k,i)=tivptl(k,j)
2219       do 3 k=1,2
2220 3     ifrptl(k,i)= 0 !ifrptl(k,j)
2221       jorptl(i)  = 0 !jorptl(j)
2222       do 4 k=1,4
2223 4     xorptl(k,i)=xorptl(k,j)
2224       do 5 k=1,4
2225 5     ibptl(k,i) =ibptl(k,j)
2226       ityptl(i)  =ityptl(j)
2227       iaaptl(i)  =iaaptl(j)
2228       radptl(i)  =radptl(j)
2229       desptl(i)  =desptl(j)
2230       dezptl(i)  =dezptl(j)
2231       qsqptl(i)  =qsqptl(j)
2232       zpaptl(i)  =zpaptl(j)
2233       itsptl(i)  =itsptl(j)
2234       rinptl(i)  =rinptl(j)
2235       return
2236       end
2237
2238 cc-----------------------------------------------------------------------
2239 c      subroutine utresm(icp1,icp2,icm1,icm2,amp,idpr,iadj,ireten)
2240 cc-----------------------------------------------------------------------
2241 c      parameter (nflav=6)
2242 c      integer icm(2),icp(2),jcm(nflav,2),jcp(nflav,2)
2243 c      icm(1)=icm1
2244 c      icm(2)=icm2
2245 c      icp(1)=icp1
2246 c      icp(2)=icp2
2247 c      CALL IDDECO(ICM,JCM)
2248 c      CALL IDDECO(ICP,JCP)
2249 c      do 37 nf=1,nflav
2250 c      do 37 k=1,2
2251 c37    jcP(nf,k)=jcp(nf,k)+jcm(nf,k)
2252 c      CALL IDENCO(JCP,ICP,IRETEN)
2253 c      IDP=IDTRA(ICP,0,0,3)
2254 c      call idres(idp,amp,idpr,iadj,0)
2255 c      return
2256 c      end
2257 c
2258 cc-----------------------------------------------------------------------
2259 c      subroutine utroa1(phi,a1,a2,a3,x1,x2,x3)
2260 cc-----------------------------------------------------------------------
2261 cc  rotates x by angle phi around axis a.
2262 cc  normalization of a is irrelevant.
2263 cc-----------------------------------------------------------------------
2264 c      double precision aaa,aa(3),xxx,xx(3),e1(3),e2(3),e3(3),xp,xt,dphi
2265 c      dphi=phi
2266 c      xx(1)=x1
2267 c      xx(2)=x2
2268 c      xx(3)=x3
2269 c      aa(1)=a1
2270 c      aa(2)=a2
2271 c      aa(3)=a3
2272 c      aaa=0
2273 c      xxx=0
2274 c      do i=1,3
2275 c      aaa=aaa+aa(i)**2
2276 c      xxx=xxx+xx(i)**2
2277 c      enddo
2278 c      if(xxx.eq.0d0)return
2279 c      if(aaa.eq.0d0)call utstop('utroa1: zero rotation axis&')
2280 c      aaa=dsqrt(aaa)
2281 c      xxx=dsqrt(xxx)
2282 cc e3 = a / !a!
2283 c      do i=1,3
2284 c      e3(i)=aa(i)/aaa
2285 c      enddo
2286 cc x_parallel
2287 c      xp=0
2288 c      do i=1,3
2289 c      xp=xp+xx(i)*e3(i)
2290 c      enddo
2291 cc x_transverse
2292 c      if(xxx**2-xp**2.le.0.)return
2293 c      xt=dsqrt(xxx**2-xp**2)
2294 cc e1 = vector x_transverse / absolute value x_transverse
2295 c      do i=1,3
2296 c      e1(i)=(xx(i)-e3(i)*xp)/xt
2297 c      enddo
2298 cc e2 orthogonal e3,e1
2299 c      call utvec2(e3,e1,e2)
2300 cc rotate x
2301 c      do i=1,3
2302 c      xx(i)=xp*e3(i)+xt*dcos(dphi)*e1(i)+xt*dsin(dphi)*e2(i)
2303 c      enddo
2304 cc back to single precision
2305 c      x1=xx(1)
2306 c      x2=xx(2)
2307 c      x3=xx(3)
2308 c      return
2309 c      end
2310 c
2311 cc-----------------------------------------------------------------------
2312 c      subroutine utroa2(phi,a1,a2,a3,x1,x2,x3)
2313 cc-----------------------------------------------------------------------
2314 cc  rotates x by angle phi around axis a.
2315 cc  normalization of a is irrelevant.
2316 cc  double precision phi,a1,a2,a3,x1,x2,x3
2317 cc-----------------------------------------------------------------------
2318 c      double precision phi,a1,a2,a3,x1,x2,x3
2319 c      double precision aaa,aa(3),xxx,xx(3),e1(3),e2(3),e3(3),xp,xt,dphi
2320 c      dphi=phi
2321 c      xx(1)=x1
2322 c      xx(2)=x2
2323 c      xx(3)=x3
2324 c      aa(1)=a1
2325 c      aa(2)=a2
2326 c      aa(3)=a3
2327 c      aaa=0
2328 c      xxx=0
2329 c      do i=1,3
2330 c      aaa=aaa+aa(i)**2
2331 c      xxx=xxx+xx(i)**2
2332 c      enddo
2333 c      if(xxx.eq.0d0)return
2334 c      if(aaa.eq.0d0)call utstop('utroa1: zero rotation axis&')
2335 c      aaa=1.0/dsqrt(aaa)
2336 cc e3 = a / !a!
2337 c      do i=1,3
2338 c      e3(i)=aa(i)*aaa
2339 c      enddo
2340 cc x_parallel
2341 c      xp=0
2342 c      do i=1,3
2343 c      xp=xp+xx(i)*e3(i)
2344 c      enddo
2345 cc x_transverse
2346 c      if(xxx-xp**2.le.0.)return
2347 c      xt=dsqrt(xxx-xp**2)
2348 cc e1 = vector x_transverse / absolute value x_transverse
2349 c      do i=1,3
2350 c      e1(i)=(xx(i)-e3(i)*xp)/xt
2351 c      enddo
2352 cc e2 orthogonal e3,e1
2353 c      call utvec2(e3,e1,e2)
2354 cc rotate x
2355 c      do i=1,3
2356 c      xx(i)=xp*e3(i)+xt*dcos(dphi)*e1(i)+xt*dsin(dphi)*e2(i)
2357 c      enddo
2358 cc back to single precision
2359 c      x1=xx(1)
2360 c      x2=xx(2)
2361 c      x3=xx(3)
2362 c      return
2363 c      end
2364 c
2365 cc--------------------------------------------------------------------
2366 c      function utroot(funcd,x1,x2,xacc)
2367 cc--------------------------------------------------------------------
2368 cc combination of newton-raphson and bisection method for root finding
2369 cc input:
2370 cc   funcd: subr returning fctn value and first derivative
2371 cc   x1,x2: x-interval
2372 cc   xacc:  accuracy
2373 cc output:
2374 cc   utroot: root
2375 cc--------------------------------------------------------------------
2376 c      include 'epos.inc'
2377 c      INTEGER MAXIT
2378 c      REAL utroot,x1,x2,xacc
2379 c      EXTERNAL funcd
2380 c      PARAMETER (MAXIT=100)
2381 c      INTEGER j
2382 c      REAL df,dx,dxold,f,fh,fl,temp,xh,xl
2383 c      call funcd(x1,fl,df)
2384 c      call funcd(x2,fh,df)
2385 c      if((fl.gt.0..and.fh.gt.0.).or.(fl.lt.0..and.fh.lt.0.))
2386 c     *call utstop('utroot: root must be bracketed&')
2387 c      if(fl.eq.0.)then
2388 c        utroot=x1
2389 c        return
2390 c      else if(fh.eq.0.)then
2391 c        utroot=x2
2392 c        return
2393 c      else if(fl.lt.0.)then
2394 c        xl=x1
2395 c        xh=x2
2396 c      else
2397 c        xh=x1
2398 c        xl=x2
2399 c      endif
2400 c      utroot=.5*(x1+x2)
2401 c      dxold=abs(x2-x1)
2402 c      dx=dxold
2403 c      call funcd(utroot,f,df)
2404 c      do 11 j=1,MAXIT
2405 c        if(((utroot-xh)*df-f)*((utroot-xl)*df-f).ge.0..or. abs(2.*
2406 c     *f).gt.abs(dxold*df) ) then
2407 c          dxold=dx
2408 c          dx=0.5*(xh-xl)
2409 c          utroot=xl+dx
2410 c          if(xl.eq.utroot)return
2411 c        else
2412 c          dxold=dx
2413 c          dx=f/df
2414 c          temp=utroot
2415 c          utroot=utroot-dx
2416 c          if(temp.eq.utroot)return
2417 c        endif
2418 c        if(abs(dx).lt.xacc) return
2419 c        call funcd(utroot,f,df)
2420 c        if(f.lt.0.) then
2421 c          xl=utroot
2422 c        else
2423 c          xh=utroot
2424 c        endif
2425 c11    continue
2426 c      call utmsg('utroot')
2427 c      write(ifch,*)'*****  exceeding maximum iterations'
2428 c      write(ifch,*)'dx:',dx
2429 c      call utmsgf
2430 c      return
2431 c      END
2432 c
2433 c-----------------------------------------------------------------------
2434       subroutine utrot2(isig,ax,ay,az,x,y,z)
2435 c-----------------------------------------------------------------------
2436 c     performs a rotation, double prec.
2437 c-----------------------------------------------------------------------
2438       include 'epos.inc'
2439       double precision ax,ay,az,x,y,z,rx,ry,rz
2440      *,alp,bet,cosa,sina,cosb,sinb,xs,ys,zs
2441          if(ax**2.eq.0.and.ay**2.eq.0.and.az**2.eq.0.)then
2442       write(ifch,*)'ax**2,ay**2,az**2:',ax**2,ay**2,az**2
2443       write(ifch,*)'ax,ay,az:',ax,ay,az
2444       call utstop('utrot2: zero vector.&')
2445          endif
2446          if(az.ge.0.)then
2447       rx=ax
2448       ry=ay
2449       rz=az
2450          else
2451       rx=-ax
2452       ry=-ay
2453       rz=-az
2454          endif
2455       if(rz**2+ry**2.ne.0.)then
2456       alp=dabs(dacos(rz/dsqrt(rz**2+ry**2)))*sign(1.,sngl(ry))
2457       bet=
2458      *dabs(dacos(dsqrt(rz**2+ry**2)/dsqrt(rz**2+ry**2+rx**2)))*
2459      *sign(1.,sngl(rx))
2460       else
2461       alp=3.1415927d0/2d0
2462       bet=3.1415927d0/2d0
2463       endif
2464       cosa=dcos(alp)
2465       sina=dsin(alp)
2466       cosb=dcos(bet)
2467       sinb=dsin(bet)
2468            if(isig.gt.0)then
2469       xs=x*cosb-y*sina*sinb-z*cosa*sinb
2470       ys=       y*cosa     -z*sina
2471       zs=x*sinb+y*sina*cosb+z*cosa*cosb
2472            elseif(isig.lt.0)then
2473       xs= x*cosb            +z*sinb
2474       ys=-x*sinb*sina+y*cosa+z*cosb*sina
2475       zs=-x*sinb*cosa-y*sina+z*cosb*cosa
2476            endif
2477       x=xs
2478       y=ys
2479       z=zs
2480       return
2481       end
2482
2483
2484 c-----------------------------------------------------------------------
2485       subroutine utrot4(isig,ax,ay,az,x,y,z)
2486 c-----------------------------------------------------------------------
2487 c     performs a rotation, double prec.
2488 c     arguments partly single
2489 c-----------------------------------------------------------------------
2490       double precision ax,ay,az,xx,yy,zz
2491       xx=dble(x)
2492       yy=dble(y)
2493       zz=dble(z)
2494       call utrot2(isig,ax,ay,az,xx,yy,zz)
2495       x=sngl(xx)
2496       y=sngl(yy)
2497       z=sngl(zz)
2498       return
2499       end
2500
2501 c-----------------------------------------------------------------------
2502       subroutine utrota(isig,ax,ay,az,x,y,z)
2503 c-----------------------------------------------------------------------
2504 c     performs a rotation
2505 c-----------------------------------------------------------------------
2506          if(az.ge.0.)then
2507       rx=ax
2508       ry=ay
2509       rz=az
2510          else
2511       rx=-ax
2512       ry=-ay
2513       rz=-az
2514          endif
2515       if(rz.eq.0..and.ry.eq.0.)then
2516         alp=0.
2517         stop
2518       else
2519         alp=abs(utacos(rz/sqrt(rz**2+ry**2)))*sign(1.,ry)
2520       endif
2521       bet=
2522      *abs(utacos(sqrt(rz**2+ry**2)/sqrt(rz**2+ry**2+rx**2)))*sign(1.,rx)
2523       cosa=cos(alp)
2524       sina=sin(alp)
2525       cosb=cos(bet)
2526       sinb=sin(bet)
2527            if(isig.gt.0)then
2528       xs=x*cosb-y*sina*sinb-z*cosa*sinb
2529       ys=       y*cosa     -z*sina
2530       zs=x*sinb+y*sina*cosb+z*cosa*cosb
2531            elseif(isig.lt.0)then
2532       xs= x*cosb            +z*sinb
2533       ys=-x*sinb*sina+y*cosa+z*cosb*sina
2534       zs=-x*sinb*cosa-y*sina+z*cosb*cosa
2535            endif
2536       x=xs
2537       y=ys
2538       z=zs
2539       return
2540       end
2541
2542 c-----------------------------------------------------------------------
2543       subroutine utstop(text)
2544 c-----------------------------------------------------------------------
2545 c  returns error message and stops execution.
2546 c  text is an optonal text to appear in the error message.
2547 c  text is a character string of length 40;
2548 c     for shorter text, it has to be terminated by &;
2549 c        example: call utstop('error in subr xyz&')
2550 c-----------------------------------------------------------------------
2551       include 'epos.inc'
2552       parameter(itext=40)
2553       character  text*40  ,txt*6
2554       imax=itext+1
2555       do i=itext,1,-1
2556       if(text(i:i).eq.'&')imax=i
2557       enddo
2558       do 1 j=1,2
2559       if(j.eq.1)ifi=ifch
2560       if(j.eq.2)ifi=ifmt
2561       write(ifi,101)('*',k=1,72),text(1:imax-1)
2562      *,nrevt+1,seedc,('*',k=1,72)
2563 101   format(
2564      *1x,72a1
2565      */1x,'***** stop in ',a
2566      */1x,'***** current event number: ',i12
2567      */1x,'***** initial seed for current event:',d25.15
2568      */1x,72a1)
2569 1     continue
2570       stop
2571       entry utmsg(txt)
2572       imsg=imsg+1
2573       write(ifch,'(1x,74a1)')('*',j=1,72)
2574       write(ifch,100)txt,nrevt+1,seedc
2575 100   format(1x,'***** msg from ',a6,'.   es:',i7,2x,d23.17)
2576       return
2577       entry utmsgf
2578       if(ish.eq.1)return
2579       write(ifch,'(1x,74a1)')('*',j=1,72)
2580       end
2581
2582 c-----------------------------------------------------------------
2583       subroutine uttrap(func,a,b,s)
2584 c-----------------------------------------------------------------
2585 c trapezoidal method for integration.
2586 c input: fctn func and limits a,b
2587 c output: value s of the integral
2588 c-----------------------------------------------------------------
2589       include 'epos.inc'
2590
2591       INTEGER JMAX
2592       REAL a,b,func,s,epsr
2593       EXTERNAL func
2594       PARAMETER (JMAX=10)
2595 CU    USES uttras
2596       INTEGER j
2597       REAL olds
2598       olds=-1.e30
2599       do 11 j=1,JMAX
2600         if(ish.ge.9)write(ifch,*)'sr uttrap:   j:',j
2601         call uttras(func,a,b,s,j)
2602         ds=abs(s-olds)
2603         if (ds.lt.epsr*abs(olds)) return
2604         olds=s
2605 11    continue
2606
2607 c-c   nepsr=nepsr+1
2608       if(ish.ge.9)then
2609       call utmsg('uttrap')
2610       write(ifch,*)
2611      *'*****  requested accuracy could not be achieved'
2612       write(ifch,*)'achieved accuracy: ',ds/abs(olds)
2613       write(ifch,*)'requested accuracy:',epsr
2614       call utmsgf
2615       endif
2616
2617       END
2618
2619 c-----------------------------------------------------------------
2620       subroutine uttraq(func,a,b,s)
2621 c-----------------------------------------------------------------
2622 c trapezoidal method for integration.
2623 c input: function func and limits a,b
2624 c output: value s of the integral
2625 c-----------------------------------------------------------------
2626
2627       REAL a,b,func,s
2628       EXTERNAL func
2629       PARAMETER (eps=1.e-6)
2630 CU    USES uttras
2631       INTEGER j
2632       REAL olds
2633       olds=-1.e30
2634       j=1
2635 10      call uttras(func,a,b,s,j)
2636         ds=abs(s-olds)
2637         if (ds.le.eps*abs(olds)) return
2638         olds=s
2639         if(j.ge.15)return
2640         j=j+1
2641       goto10
2642       END
2643
2644 c-----------------------------------------------------------------
2645       subroutine uttras(func,a,b,s,n)
2646 c-----------------------------------------------------------------
2647 c performs one iteration of the trapezoidal method for integration
2648 c-----------------------------------------------------------------
2649       INTEGER n
2650       REAL a,b,s,func
2651       EXTERNAL func
2652       INTEGER it,j
2653       REAL del,sum,tnm,x
2654       if (n.eq.1) then
2655         s=0.5*(b-a)*(func(a)+func(b))
2656       else
2657         it=2**(n-2)
2658         tnm=it
2659         del=(b-a)/tnm
2660         x=a+0.5*del
2661         sum=0.
2662         do 11 j=1,it
2663           sum=sum+func(x)
2664           x=x+del
2665 11      continue
2666         s=0.5*(s+(b-a)*sum/tnm)
2667       endif
2668       return
2669       END
2670
2671 cc-----------------------------------------------------------------------
2672 c      subroutine utvec1(a,b,c)
2673 cc-----------------------------------------------------------------------
2674 cc  returns vector product c = a x b .
2675 cc-----------------------------------------------------------------------
2676 c      real a(3),b(3),c(3)
2677 c      c(1)=a(2)*b(3)-a(3)*b(2)
2678 c      c(2)=a(3)*b(1)-a(1)*b(3)
2679 c      c(3)=a(1)*b(2)-a(2)*b(1)
2680 c      return
2681 c      end
2682 c
2683 cc-----------------------------------------------------------------------
2684 c      subroutine utvec2(a,b,c)
2685 cc-----------------------------------------------------------------------
2686 cc  returns vector product c = a x b .
2687 cc  a,b,c double precision.
2688 cc-----------------------------------------------------------------------
2689 c      double precision a(3),b(3),c(3)
2690 c      c(1)=a(2)*b(3)-a(3)*b(2)
2691 c      c(2)=a(3)*b(1)-a(1)*b(3)
2692 c      c(3)=a(1)*b(2)-a(2)*b(1)
2693 c      return
2694 c      end
2695 c
2696 c-------------------------------------------------------------------
2697       subroutine utword(line,i,j,iqu)
2698 c-------------------------------------------------------------------
2699 c  finds the first word of the character string line(j+1:160).
2700 c  the word is line(i:j) (with new i and j).
2701 c  if j<0 or if no word found --> new line read.
2702 c  a text between quotes "..." is considered one word;
2703 c  stronger: a text between double quotes ""..."" is consid one word
2704 c  stronger: a text between "{ and }" is considered one word
2705 c-------------------------------------------------------------------
2706 c  input:
2707 c    line: character string (*160)
2708 c    i: integer between 1 and 160
2709 c    iqu: for iqu=1 a "?" is written to output before reading a line,
2710 c         otherwise (iqu/=1) nothing is typed
2711 c  output:
2712 c    i,j: left and right end of word (word=line(i:j))
2713 c-------------------------------------------------------------------
2714       include 'epos.inc'
2715       parameter(mempty=2)
2716       character*1 empty(mempty),mk
2717       character line*160
2718       character*2 mrk
2719       data empty/' ',','/
2720       parameter(mxdefine=40)
2721       character w1define*100,w2define*100
2722       common/cdefine/ndefine,l1define(mxdefine),l2define(mxdefine)
2723      &               ,w1define(mxdefine),w2define(mxdefine)
2724
2725       if(j.ge.0)then
2726       i=j
2727       goto1
2728       endif
2729
2730     5 continue
2731       if(iqu.eq.1.and.iprmpt.gt.0)write(ifmt,'(a)')'?'
2732       if(nopen.eq.0)ifopx=ifop
2733       if(nopen.gt.0)ifopx=20+nopen
2734       if(nopen.lt.0)ifopx=ifcp
2735       read(ifopx,'(a160)',end=9999)line
2736       if(iecho.eq.1.or.(nopen.ge.0.and.kcpopen.eq.1))then
2737        kmax=2
2738        do k=3,160
2739        if(line(k:k).ne.' ')kmax=k
2740        enddo
2741       endif
2742       if(nopen.ge.0.and.kcpopen.eq.1)
2743      &  write(ifcp,'(a)')line(1:kmax)
2744       if(iecho.eq.1)
2745      &  write(ifmt,'(a)')line(1:kmax)
2746       i=0
2747
2748     1 i=i+1
2749       if(i.gt.160)goto5
2750       if(line(i:i).eq.'!')goto5
2751       do ne=1,mempty
2752       if(line(i:i).eq.empty(ne))goto1
2753       enddo
2754
2755       mrk='  '
2756       mk=' '
2757       if(line(i:i).eq.'~')mk='~'
2758       if(line(i:i+1).eq.'"{')mrk='}"'
2759       if(line(i:i+1).eq.'""')mrk='""'
2760       if(mrk.ne.'  ')goto10
2761       if(line(i:i).eq.'"')mk='"'
2762       if(mk.ne.' ')goto8
2763       j=i-1
2764     6 j=j+1
2765       if(j.gt.160)goto7
2766       if(line(j:j).eq.'!')goto7
2767       do ne=1,mempty
2768       if(line(j:j).eq.empty(ne))goto7
2769       enddo
2770       goto6
2771
2772     8 continue
2773       if(i.ge.160-1)stop'utword: make line shorter!!!         '
2774       i=i+1
2775       j=i
2776       if(line(j:j).eq.mk)stop'utword: empty string!!!           '
2777     9 j=j+1
2778       if(j.gt.160)goto7
2779       if(line(j:j).eq.mk)then
2780       line(i-1:i-1)=' '
2781       line(j:j)=' '
2782       goto7
2783       endif
2784       goto9
2785
2786    10 continue
2787       if(i.ge.160-3)stop'utword: make line shorter!!!!          '
2788       i=i+2
2789       j=i
2790       if(line(j:j+1).eq.mrk)stop'utword: empty string!!!!        '
2791    11 j=j+1
2792       if(j+1.gt.160)goto7
2793       if(line(j:j+1).eq.mrk)then
2794       line(i-2:i-1)='  '
2795       line(j:j+1)='  '
2796       goto7
2797       endif
2798       goto11
2799
2800     7 j=j-1
2801       !--------#define---------------
2802       if(ndefine.gt.0)then
2803         do ndf=1,ndefine
2804           l1=l1define(ndf)
2805           l2=l2define(ndf)
2806           do i0=i,j+1-l1
2807             if(line(i0:i0-1+l1).eq.w1define(ndf)(1:l1))then
2808               if(l2.eq.l1)then
2809                 line(i0:i0-1+l1)=w2define(ndf)(1:l2)
2810               elseif(l2.lt.l1)then
2811                 line(i0:i0+l2-1)=w2define(ndf)(1:l2)
2812                 do k=i0+l2,i0-1+l1
2813                   line(k:k)=' '
2814                 enddo
2815               elseif(l2.gt.l1)then
2816                 do k=i0+l1,i0+l2-1
2817                   if(line(k:k).ne.' ')
2818      &        stop'utword: no space for `define` replacement.   '
2819                 enddo
2820                 line(i0:i0+l2-1)=w2define(ndf)(1:l2)
2821                 j=i0+l2-1
2822               endif
2823             endif
2824           enddo
2825         enddo
2826        do k=i,j
2827          if(line(k:k).ne.' ')j0=j
2828        enddo
2829        j=j0
2830       endif
2831       !--------
2832       return
2833
2834 9999  close(ifopx)
2835       nopen=nopen-1
2836       if(nopen.eq.0.and.iprmpt.eq.-1)iprmpt=1
2837       goto5
2838       end
2839
2840 c--------------------------------------------------------------------
2841       subroutine utworn(line,j,ne)
2842 c--------------------------------------------------------------------
2843 c  returns number ne of nonempty characters of line(j+1:160)
2844 c--------------------------------------------------------------------
2845       character line*160
2846       ne=0
2847       do l=j+1,160
2848       if(line(l:l).ne.' ')ne=ne+1
2849       enddo
2850       return
2851       end
2852
2853 c-----------------------------------------------------------------------
2854       subroutine getairmol(iz,ia)
2855 c-----------------------------------------------------------------------
2856       include 'epos.inc'
2857       i=0
2858       r=rangen()
2859       do while(r.gt.0.)  ! choose air-molecule
2860         i=i+1
2861         r=r-airwnxs(i)
2862       enddo
2863       iz = nint(airznxs(i))
2864       ia = nint(airanxs(i))
2865       end
2866
2867 c----------------------------------------------------------------------
2868
2869       subroutine factoriel
2870
2871 c----------------------------------------------------------------------
2872 c tabulation of fctrl(n)=n!, facto(n)=1/n! and utgamtab(x) for x=0 to 50
2873 c----------------------------------------------------------------------
2874       include 'epos.incems'
2875       double precision utgamtab,utgam,x
2876       common/gamtab/utgamtab(10000)
2877
2878       nfctrl=100
2879       fctrl(0)=1.D0
2880       facto(0)=1.D0
2881       do i=1,min(npommx,nfctrl)
2882         fctrl(i)=fctrl(i-1)*dble(i)
2883         facto(i)=1.d0/fctrl(i)
2884       enddo
2885
2886       do k=1,10000
2887         x=dble(k)/100.d0
2888         utgamtab(k)=utgam(x)
2889       enddo
2890
2891       return
2892       end
2893
2894 c-----------------------------------------------------------------------
2895       subroutine fremnu(amin,ca,cb,ca0,cb0,ic1,ic2,ic3,ic4)
2896 c-----------------------------------------------------------------------
2897       common/hadr2/iomodl,idproj,idtarg,wexcit
2898       real pnll,ptq
2899       common/hadr1/pnll,ptq,exmass,cutmss,wproj,wtarg
2900       parameter (nflav=6)
2901       integer ic(2),jc(nflav,2)
2902       ic(1)=ca
2903       ic(2)=cb
2904       call iddeco(ic,jc)
2905       keu=jc(1,1)-jc(1,2)
2906       ked=jc(2,1)-jc(2,2)
2907       kes=jc(3,1)-jc(3,2)
2908       kec=jc(4,1)-jc(4,2)
2909       keb=jc(5,1)-jc(5,2)
2910       ket=jc(6,1)-jc(6,2)
2911       amin=utamnu(keu,ked,kes,kec,keb,ket,5)  !???4=2mults, 5=1mult
2912         if(ca-ca0.eq.0..and.cb-cb0.eq.0..and.rangen().gt.wexcit)then
2913       ic3=0
2914       ic4=0
2915       ic1=ca
2916       ic2=cb
2917        else
2918       amin=amin+exmass
2919       n=0
2920       do i=1,4
2921       do j=1,2
2922       n=n+jc(i,j)
2923       enddo
2924       enddo
2925       k=1+rangen()*n
2926       do i=1,4
2927       do j=1,2
2928       k=k-jc(i,j)
2929       if(k.le.0)goto 1
2930       enddo
2931       enddo
2932 1     if(j.eq.1)then
2933       ic3=10**(6-i)
2934       ic4=0
2935       else
2936       ic3=0
2937       ic4=10**(6-i)
2938       endif
2939       ic1=int(ca)-ic3
2940       ic2=int(cb)-ic4
2941         endif
2942       return
2943       end
2944
2945
2946 c-----------------------------------------------------------------------
2947       function fremnux(ic)
2948 c-----------------------------------------------------------------------
2949       real pnll,ptq
2950       common/hadr1/pnll,ptq,exmass,cutmss,wproj,wtarg
2951       parameter (nflav=6)
2952       integer ic(2),jc(nflav,2)
2953 c      ic(1)=210000
2954 c      ic(2)=0
2955       call iddeco(ic,jc)
2956       keu=jc(1,1)-jc(1,2)
2957       ked=jc(2,1)-jc(2,2)
2958       kes=jc(3,1)-jc(3,2)
2959       kec=jc(4,1)-jc(4,2)
2960       keb=jc(5,1)-jc(5,2)
2961       ket=jc(6,1)-jc(6,2)
2962       fremnux=utamnu(keu,ked,kes,kec,keb,ket,4) !+exmass  !???4=2mults, 5=1mult
2963       return
2964       end
2965
2966 c-----------------------------------------------------------------------
2967       function fremnux2(ic)
2968 c-----------------------------------------------------------------------
2969       real pnll,ptq
2970       common/hadr1/pnll,ptq,exmass,cutmss,wproj,wtarg
2971       parameter (nflav=6)
2972       integer ic(2),jc(nflav,2)
2973 c      ic(1)=210000
2974 c      ic(2)=0
2975       call iddeco(ic,jc)
2976       keu=jc(1,1)-jc(1,2)
2977       ked=jc(2,1)-jc(2,2)
2978       kes=jc(3,1)-jc(3,2)
2979       kec=jc(4,1)-jc(4,2)
2980       keb=jc(5,1)-jc(5,2)
2981       ket=jc(6,1)-jc(6,2)
2982       fremnux2=utamnu(keu,ked,kes,kec,keb,ket,5) !+exmass  !???4=2mults, 5=1mult
2983       return
2984       end
2985
2986 c-----------------------------------------------------------------------
2987       subroutine fremnx(ammax,amin,sm,ic3,ic4,iret)
2988 c-----------------------------------------------------------------------
2989       common/psar9/ alpr
2990       include 'epos.inc'
2991       iret=0
2992       if(ic3.eq.0.and.ic4.eq.0)then
2993         if(ammax.lt.amin**2)then
2994           iret=1
2995           return
2996         endif
2997         sm=amin**2
2998       else
2999 c       ammax1=min(ammax,(engy/4.)**2)
3000         ammax1=ammax
3001         if(ammax1.lt.amin**2)then
3002           iret=1
3003           return
3004         endif
3005         if(alpr.eq.-1.)then
3006           sm=amin**2*(ammax1/amin**2)**rangen()
3007         else
3008           sm=amin**2*(1.+((ammax1/amin**2)**(1.+alpr)-1.)
3009      *    *rangen())**(1./(1.+alpr))
3010         endif
3011       endif
3012       return
3013       end
3014
3015