Back to home page

Project CMSSW displayed by LXR

 
 

    


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

0001 c-----------------------------------------------------------------------
0002       subroutine utresc(iret)
0003 c-----------------------------------------------------------------------
0004 c  if irescl=1 rescaling is done, otherwise the purpose of going through
0005 c  this routine is to not change the seed in case of no rescaling
0006 c-----------------------------------------------------------------------
0007       include 'epos.inc'
0008       include 'epos.incems'
0009       double precision p1,esoll,ppp,seedp,psoll,pt1soll,pt2soll
0010       dimension p1(5),p2(4),p0(5,mamx+mamx),pini(mxptl)
0011       logical force,nolead(mxptl),lspec(mxptl),lim
0012       data scalmean/0./scalevt/0./
0013       save scalmean,scalevt
0014       call utpri('utresc',ish,ishini,4)
0015 
0016       errlim=0.005 !max(0.001,1./engy)
0017       if(iLHC.eq.1)errlim=max(0.00005,0.5/engy)
0018 
0019       iret=0
0020       nptlpt=iabs(maproj)+iabs(matarg)
0021       call ranfgt(seedp)        !not to change the seed ...
0022       if(nptl.le.nptlpt) goto 9999
0023 
0024       if(ish.ge.8)then
0025         call alistf('list before boost&')
0026       endif
0027       esoll=0.d0
0028       psoll=0.d0
0029       p1(1)=0.d0
0030       p1(2)=0.d0
0031       p1(3)=0.d0
0032       p1(4)=0.d0
0033       p2(3)=0.d0
0034       p2(4)=0.d0
0035       ipmax=4
0036       imin=nptlpt+1
0037       if(iappl.eq.1)then
0038         imin=1
0039         ipmax=2
0040         if(iLHC.eq.1)ipmax=0
0041       endif
0042 c store projectile and target in p0 and sum pz an E in p1(3) and p1(4)
0043       do i=1,nptlpt
0044         nolead(i)=.false.
0045         do j=1,5
0046           p0(j,i)=dble(pptl(j,i))
0047         enddo
0048 c calculate total energy of primaries
0049 c       if(mod(istptl(i),10).eq.1)then
0050           do j=ipmax+1,4
0051             p1(j)=p1(j)+dble(pptl(j,i))
0052           enddo
0053 c       endif
0054 c calculate total energy of secondaries
0055         if(mod(istptl(i),10).eq.0)then
0056           do j=ipmax+1,4
0057             p2(j)=p2(j)+dble(pptl(j,i))
0058           enddo
0059         endif
0060       enddo
0061 c fix secondaries counted in the system
0062       do i=nptlpt+1,nptl
0063        if(mod(istptl(i),10).eq.0)then
0064 c check maximum energy
0065          if(iLHC.eq.1.and.pptl(4,i).gt.engy*0.51)then
0066            if(ish.ge.1)write(ifch,*)'Energy of particle too high !'
0067      &                                     ,i,ityptl(i),pptl(4,i)
0068 c           call utstop('utresc&')
0069            if(ityptl(i).eq.48.or.ityptl(i).eq.58  !diffractive resonance
0070      &    .or.ityptl(i).eq.47.or.ityptl(i).eq.57)then   !active spectators
0071              pptl(4,i)=0.5*engy
0072              amt=sqrt(pptl(1,i)**2+pptl(2,i)**2+pptl(5,i)**2)
0073              pptl(3,i)=(pptl(4,i)+amt)*(pptl(4,i)-amt)
0074              if(pptl(3,i).gt.0.)then
0075                pptl(3,i)=sqrt(pptl(3,i))
0076              else
0077                iret=1
0078              endif
0079            else
0080              iret=1
0081            endif
0082            if(iret.eq.1)then
0083              if(ish.ge.1)write(ifch,*)'Warning in utresc: redo event...'
0084 c             call utstop('Energy of particle too high !&')
0085              goto 9999
0086            endif
0087          endif
0088 c fix pt (p1(1) and p2(1)) from secondaries
0089          do j=1,ipmax
0090            p1(j)=p1(j)+dble(pptl(j,i))
0091          enddo
0092 c calculate total energy of secondaries
0093          do j=ipmax+1,4
0094            p2(j)=p2(j)+dble(pptl(j,i))
0095          enddo
0096          lspec(i)=.false.
0097          if(((ityptl(i).eq.45.or.ityptl(i).eq.48).and.maproj.ge.100)
0098      &  .or.((ityptl(i).eq.55.or.ityptl(i).eq.58).and.matarg.ge.100))
0099      &   lspec(i)=.true.
0100          if((abs(pptl(3,i)/pnullx).le.0.9
0101      & .and.abs(pptl(3,i)).gt.pptl(5,i)).or.lspec(i))then
0102            nolead(i)=.true.
0103 c           write(ifch,*)'nolead',i
0104          else
0105            nolead(i)=.false.
0106 c           write(ifch,*)'lead',i
0107          endif
0108        endif
0109       enddo
0110       psoll=max(dble(pnullx),abs(p1(3)))
0111 c check if energy is conserved before boost
0112       if(iappl.eq.1)then
0113         diff=abs(p2(3)-p1(3))
0114         scal=p2(4)/p1(4)
0115         if(abs(scal-1.).le.errlim.and.abs(diff/psoll).lt.errlim
0116      &     .and.(iLHC.eq.0.or.
0117      &    (abs(p2(1)).lt.errlim.and.abs(p2(2)).lt.errlim)))then
0118           if(ish.ge.4) 
0119      & write (ifch,'(2x,a,2g14.6)') 'Energy OK: ',scal,abs(diff/psoll)
0120           goto 9999
0121         else
0122          diff=0.
0123          scal=1.
0124         endif
0125       endif
0126 c calculate boost vector to have sum of pt=0 and sum of pz=0
0127       ppp=(p1(4)+p1(3))*(p1(4)-p1(3))-p1(2)*p1(2)-p1(1)*p1(1)
0128       if(ppp.gt.0.d0)then
0129         p1(5)=sqrt(ppp)
0130       else
0131         iret=1
0132         write(ifch,*)'p1',p1(1),p1(2),p1(3),p1(4),ppp
0133         if(ish.ge.2)write (ifch,*) 'Problem in utresc (1): redo event'
0134         write (ifmt,*) 'Problem in utresc (1): redo event'
0135 c       call utstop('utresc&')
0136         goto 9999
0137       endif
0138       esoll=p1(5)
0139       if(ish.ge.4) write (ifch,'(a,5g14.6)') 'boost-vector: ',p1
0140 
0141 c     trafo
0142 c     -----
0143       pmax=0.d0
0144       npart=0
0145       pt1soll=0.d0
0146       pt2soll=0.d0
0147       do i=imin,nptl
0148         if(mod(istptl(i),10).le.1)then
0149           call utlob4(1,p1(1),p1(2),p1(3),p1(4),p1(5)
0150      $         ,pptl(1,i),pptl(2,i),pptl(3,i),pptl(4,i))
0151           if(mod(istptl(i),10).eq.0.and.i.gt.nptlpt)then
0152             npart=npart+1
0153             pt1soll=pt1soll+dble(pptl(1,i))
0154             pt2soll=pt2soll+dble(pptl(2,i))
0155           endif
0156         endif
0157         if(i.gt.nptlpt.and.nolead(i))pmax=max(pmax,abs(pptl(3,i)))
0158       enddo
0159 
0160       if(ish.ge.6)then
0161         call alistf('list after boost&')
0162       endif
0163 
0164       if(ish.ge.5)write(ifch,'(a)')'--------rescale momenta----------'
0165 
0166       if(iLHC.eq.1)then
0167 
0168         if(ish.ge.6)write(ifch,*)'ptsoll:',pt1soll,pt2soll,npart
0169         pt1soll=pt1soll/dble(npart)
0170         pt2soll=pt2soll/dble(npart)
0171         do i=nptlpt+1,nptl
0172           if(mod(istptl(i),10).eq.0)then
0173             pptl(1,i)=pptl(1,i)-sngl(pt1soll)
0174             pptl(2,i)=pptl(2,i)-sngl(pt2soll)
0175             pptl(4,i)=sqrt(pptl(1,i)**2+pptl(2,i)**2
0176      &                    +pptl(3,i)**2+pptl(5,i)**2)
0177           endif
0178         enddo
0179 
0180       endif
0181 
0182       if(ish.ge.6)write(ifch,*)'esoll,psoll,pmax:',esoll,psoll,pmax
0183 
0184 c     rescale momenta in rest frame
0185 c     -----------------------------
0186       scal=1.
0187       diff0=0.
0188 c      ndif0=1
0189       ferr=0.05
0190       force=.false.
0191       npart=nptl-imin+1
0192       do ipass=1,300
0193         sum=0.
0194         sum3=0.
0195         difft=diff0
0196         ndif=0
0197         nfirst=int(rangen()*float(npart))    !start list randomly
0198         do  i=0,npart-1
0199           j=imin+i+nfirst
0200           if(j.gt.nptl)j=imin+i+nfirst-npart
0201           if(mod(istptl(j),10).eq.0)then
0202 c modified particles
0203             if(nolead(j))then
0204 c            if(j.gt.nptlpt)then
0205 c            if(abs(pptl(3,j))/pnullx.lt.0.9)then  !not spectator or diffraction
0206               if(scal.eq.1..and.abs(diff0).lt.1.e-6)then
0207                 ndif=ndif+1
0208                 pini(j)=pptl(3,j)
0209               else
0210                 pptl3new=0.
0211                 if( force .or.(
0212      &            ityptl(j)/10.eq.4.or.ityptl(j)/10.eq.5
0213      &                      ))then !try just remnant first
0214                   ndif=ndif+1
0215                   diff=sign(min(0.3*abs(pini(j)),
0216      &                      rangen()*abs(difft)),difft)
0217                   pptl3new=scal*(pptl(3,j)-diff)
0218 c                  write(ifch,*)'par',j,pptl3new,pptl(3,j),diff,difft
0219 c     &                 ,ndif,pmax,scal
0220                   if(abs(pptl3new).lt.pmax)then
0221 c particle should not be too fast or too modified
0222                     if(abs(pptl3new-pini(j)).lt.ferr*abs(pini(j))
0223      &         .or.(lspec(j).and.abs(pptl3new).lt.abs(0.8*pini(j))))then
0224 c                  write(ifch,*)'used'
0225                       difft=difft-diff
0226                       pptl(3,j)=scal*(pptl(3,j)-diff)
0227                       pptl(4,j)=sqrt(pptl(1,j)**2+pptl(2,j)**2
0228      *                     +pptl(3,j)**2+pptl(5,j)**2)
0229                     endif
0230                   endif
0231                 endif
0232               endif
0233             endif
0234 c sum over all particles
0235             sum=sum+pptl(4,j)
0236             sum3=sum3+pptl(3,j)
0237           endif
0238         enddo
0239 
0240         diff=sum3
0241         scal=sngl(esoll)/(sum-diff)
0242         if(ish.ge.6)write(ifch,*)
0243      $       'ipass,scal,diff/psoll,e,pz,ndif,f:'
0244      $   ,ipass,scal,diff/psoll,sum,sum3,ndif,force
0245         if(abs(scal-1.).le.errlim.and.abs(diff/psoll).lt.errlim)
0246      $  goto 300
0247         if(ndif.gt.0.and.(force.or.ipass.lt.150))then
0248 c          ndif0=ndif
0249           diff0=diff
0250         elseif(abs(scal-1.).le.1e-2.and.abs(diff/psoll).lt.5e-2
0251      &                                           .and.iLHC.eq.0)then
0252           goto 300
0253         elseif(force)then
0254           if(ish.ge.2)
0255      $    write(ifmt,*)'Warning in utresc: no more allowed particle'
0256           goto 302
0257         else
0258           force=.true.
0259           ferr=0.1
0260           diff=0.
0261         endif
0262       enddo
0263  302  if(iLHC.eq.1)then
0264         lim=.not.(abs(scal-1.).le.errlim.and.abs(diff/psoll).lt.errlim)
0265       else
0266         lim=abs(scal)+abs(diff/psoll).gt.2.5
0267       endif
0268  
0269       if(ish.ge.1)then
0270         call utmsg('utrescl')
0271         write(ifch,*)'*****  scal=',scal,diff/psoll,lim
0272         call utmsgf
0273       endif
0274 
0275       if(lim)then
0276         if(ish.ge.1)then
0277           write(ifmt,*)'Warning in utresc !'
0278           write(ifch,'(a,i10,d25.15)')'redo EPOS event ...'
0279      &                                ,nint(seedj),seedc
0280         endif
0281         iret=1
0282         goto 9999
0283       endif
0284 
0285 c     trafo
0286 c     -----
0287  300  continue
0288 
0289       do i=1,nptl
0290         if(i.le.nptlpt)then
0291           do j=1,5
0292             pptl(j,i)=p0(j,i)
0293           enddo
0294         else
0295           if(mod(istptl(i),10).le.1)then
0296             call utlob4(-1,p1(1),p1(2),p1(3),p1(4),p1(5)
0297      $           ,pptl(1,i),pptl(2,i),pptl(3,i),pptl(4,i))
0298           endif
0299         endif
0300       enddo
0301 
0302       if(ish.ge.4)call alist('list after rescaling&',1,nptl)
0303 
0304  9999 continue
0305       if(ish.ge.2)then
0306         scalevt=scalevt+1.
0307         scalmean=scalmean+scal
0308         write(ifch,*)' average rescaling factor: ',scalmean
0309      &                                            /scalevt
0310       endif
0311       call ranfst(seedp)        ! ... after this subroutine
0312       call utprix('utresc',ish,ishini,4)
0313 
0314       end
0315 
0316 c-----------------------------------------------------------------------
0317       subroutine utghost(iret)
0318 c-----------------------------------------------------------------------
0319 c  if irescl=1 make particle on-shell if not
0320 c-----------------------------------------------------------------------
0321       include 'epos.inc'
0322       include 'epos.incems'
0323       double precision seedp
0324       call utpri('ughost',ish,ishini,4)
0325 
0326       iret=0
0327       nptlpt=iabs(maproj)+iabs(matarg)
0328       if(iappl.eq.6.or.iappl.eq.8)nptlpt=3   ! ee or DIS
0329       call ranfgt(seedp)        !not to change the seed ...
0330       if(nptl.le.nptlpt) goto 9999
0331 
0332       if(ish.ge.5)write(ifch,'(a)')'---------mark ghosts---------'
0333 
0334 c     mark ghosts
0335 c     -----------
0336       do  j=nptlpt+1,nptl
0337         if(istptl(j).le.1.and.pptl(4,j).gt.0.d0)then
0338           if(iLHC.eq.0.or.mod(abs(idptl(j)),10).le.1)then !for LHC tune don't fix mass of resonnances (to keep width)
0339             amass=pptl(5,j)
0340             call idmass(idptl(j),amass)
0341             if(abs(idptl(j)).gt.100.and.
0342      &       abs(pptl(5,j)-amass).gt.0.01*amass)then
0343               if(ish.ge.5)write(ifch,*)'wrong particle mass',j,idptl(j)
0344      &                                           ,pptl(5,j),amass
0345               amass=pptl(5,j)
0346               call idres(idptl(j),amass,idr,iadj)
0347               if(idr.ne.0)then
0348                 pptl(5,j)=amass
0349                 idptl(j)=idr
0350               else
0351                 call idmass(idptl(j),amass)
0352                 pptl(5,j)=amass
0353               endif
0354               call idtau(idptl(j),pptl(4,j),pptl(5,j),taugm)
0355               tivptl(2,j)=tivptl(1,j)+taugm*(-alog(rangen()))
0356             else
0357               pptl(5,j)=amass
0358             endif
0359           endif
0360           if(abs((pptl(4,j)+pptl(3,j))*(pptl(4,j)-pptl(3,j))
0361      $         -pptl(2,j)**2-pptl(1,j)**2-pptl(5,j)**2).gt.0.3
0362      $       .and.abs(1.-abs(pptl(3,j))/pptl(4,j)).gt.0.01)then
0363         !print*,'ghost',ityptl(j),idptl(j)
0364            if(ish.ge.1)write(ifmt,*)'ghost:',j,idptl(j),ityptl(j)
0365            if(ish.ge.5)then
0366               write(ifch,'(a,$)')'ghost:'
0367               call alistc("&",j,j)
0368             endif
0369             ityptl(j)=100+ityptl(j)/10
0370           elseif(irescl.ge.1)then
0371 c ensure that all particles are really on-shell
0372             pptl(4,j)=sqrt(pptl(1,j)**2+pptl(2,j)**2
0373      *                    +pptl(3,j)**2+pptl(5,j)**2)
0374           endif
0375         elseif(mod(istptl(j),10).eq.0)then
0376 c if not droplet with fusion
0377           if(istptl(j).ne.10.or.iorsdf.ne.3)then
0378             if(ish.ge.1)then
0379               write(ifmt,*)'Lost particle (E=0)'
0380               write(ifch,*)'Lost particle (E=0) :'
0381               call alistc("utghost&",j,j)
0382             endif
0383             istptl(j)=istptl(j)+2
0384           endif
0385         endif
0386       enddo
0387 
0388       if(ish.ge.5)write(ifch,'(a)')'---------treat ghosts---------'
0389 
0390 c     treat ghosts
0391 c     ------------
0392       ifirst=1
0393       scal=1.
0394       pfif=0.
0395       efif=0.
0396       ntry=0
0397  132  nfif=0
0398       psum=0
0399       esum=0.
0400       ntry=ntry+1
0401       do  j=nptlpt+1,nptl
0402         if(mod(istptl(j),10).eq.0)then
0403           if(ityptl(j).gt.100)then
0404             nfif=nfif+1
0405             if(ifirst.eq.1)then
0406               pfif=pfif+pptl(3,j)
0407               if(pptl(4,j).gt.0.)efif=efif+pptl(4,j)
0408             endif
0409             if(irescl.ge.1) then
0410               if(ifirst.gt.1)then
0411                 if(pptl(4,j).gt.0.)then
0412                   Einv=1./pptl(4,j)
0413                   amt=1.-(pptl(5,j)*Einv)**2+(pptl(1,j)*Einv)**2
0414      $                +(pptl(2,j)*Einv)**2
0415                 else
0416                   amt=-1.
0417                 endif
0418                 if(amt.gt.0.)then
0419                   pptl(3,j)=sign(pptl(4,j),pptl(3,j))*sqrt(amt)
0420                 else
0421                   y=(rangen()+rangen()+rangen()+rangen()-2.)/2.*yhaha
0422                   y=sign(abs(y),pptl(3,j))
0423                   pptl(3,j)
0424      $                 =sqrt(pptl(5,j)**2+pptl(1,j)**2
0425      $                 +pptl(2,j)**2)*sinh(y)
0426                   pptl(4,j)
0427      $                 =sqrt(pptl(5,j)**2+pptl(1,j)**2
0428      $                 +pptl(2,j)**2)*cosh(y)
0429                   efif=efif+pptl(4,j)
0430                 endif
0431                 ifirst=0
0432               else
0433 c                do k=1,3
0434                 do k=3,3
0435                   pptl(k,j)=pptl(k,j)*scal
0436                 enddo
0437                 pptl(4,j)=sqrt(pptl(1,j)**2+pptl(2,j)**2+pptl(3,j)**2
0438      *                 +pptl(5,j)**2)
0439               endif
0440             endif
0441             psum=psum+pptl(3,j)
0442             esum=esum+pptl(4,j)
0443             if(ish.ge.5)
0444      $           write (ifch,*) 'nrevt,psum,esum,pfif,efif,nfif,scal'
0445      $           ,nrevt,psum,esum,pfif,efif,nfif,scal
0446           endif
0447         endif
0448       enddo
0449       if ( ish.ge.5 )  write (ifch,*) 'tot',nfif,efif,pfif,esum,psum
0450 
0451 
0452       if(nfif.gt.5.or.(esum.gt.0.05*engy.and.nfif.ne.1))then
0453         if(ifirst.eq.0)then
0454           do  j=nptlpt+1,nptl
0455             if ( ityptl(j).ge.101 .and. ityptl(j).le.105 )then
0456               if((psum-pfif)*(1.-scal).ge.0)
0457      &             pptl(3,j)=pptl(3,j)-(psum-pfif)/nfif
0458             endif
0459           enddo
0460         else
0461           ifirst=2
0462           goto 132
0463         endif
0464         scal=efif/esum
0465         if ( ish.ge.5 )  write (ifch,*) 'scal',scal
0466         if ( abs(scal-1.) .gt. 0.05 ) then
0467           if(ntry.le.1000)then
0468             goto 132
0469           else
0470             iret=1
0471             if(ish.ge.2)write (ifch,*) 'Problem in utghost : redo event'
0472             if(ish.ge.1)write (ifmt,*) 'Problem in utghost : redo event'
0473             goto 9999
0474          endif
0475         endif
0476       else
0477         do  j=nptlpt+1,nptl
0478           if ( ityptl(j).ge.101 .and. ityptl(j).le.105 )then
0479             pptl(4,j)=sqrt(pptl(1,j)**2+pptl(2,j)**2+pptl(3,j)**2
0480      *                 +pptl(5,j)**2)
0481           endif
0482         enddo
0483       endif
0484 
0485       if(ish.ge.5)write(ifch,'(a)')'---------Check Ghost list---------'
0486 
0487 c Check Ghost list
0488 
0489       if(ish.ge.5)then
0490         do  j=nptlpt+1,nptl
0491           if(mod(istptl(j),10).eq.0)then
0492             if(ityptl(j).le.105.and.ityptl(j).ge.101)then
0493               write(ifch,'(a,$)')'ghost:'
0494               call alistc("&",j,j)
0495             endif
0496           endif
0497         enddo
0498       endif
0499 
0500 
0501  9999 continue
0502       call ranfst(seedp)        ! ... after this subroutine
0503       call utprix('ughost',ish,ishini,4)
0504 
0505       end
0506 
0507 c-----------------------------------------------------------------------
0508       subroutine utrsph(iret)
0509 c-----------------------------------------------------------------------
0510 c  if irescl=1 and ispherio=1 rescaling is done for particle used by
0511 c  spherio as initial condition.
0512 c-----------------------------------------------------------------------
0513       include 'epos.inc'
0514       include 'epos.incems'
0515       double precision p1,esoll
0516       dimension p1(5),p0(5,mamx+mamx)
0517       call utpri('utrsph',ish,ishini,4)
0518 
0519       errlim=0.0001
0520 
0521       iret=0
0522       nptlpt=iabs(maproj)+iabs(matarg)
0523       if(nptl.le.nptlpt) goto 9999
0524 
0525       esoll=0.d0
0526       p1(1)=0.d0
0527       p1(2)=0.d0
0528       p1(3)=0.d0
0529       p1(4)=0.d0
0530       do i=nptlpt+1,nptl
0531         if((istptl(i).le.11
0532      $   .and.(iorptl(i).ge.1.and.istptl(iorptl(i)).eq.41))
0533      $   .or.istptl(i).eq.20.or.istptl(i).eq.21)then
0534          do j=1,2
0535            p1(j)=p1(j)+dble(pptl(j,i))
0536          enddo
0537        endif
0538       enddo
0539       do i=1,nptlpt
0540          do j=1,5
0541            p0(j,i)=pptl(j,i)
0542          enddo
0543          do j=3,4
0544            p1(j)=p1(j)+dble(pptl(j,i))
0545          enddo
0546       enddo
0547       p1(5)=dsqrt((p1(4)+p1(3))*(p1(4)-p1(3))-p1(2)**2.d0-p1(1)**2.d0)
0548       esoll=p1(5)
0549       if(ish.ge.4) write (ifch,'(a,5g13.6)') 'boost-vector',p1
0550 
0551 c     trafo
0552 c     -----
0553       do i=1,nptl
0554         if((istptl(i).le.11
0555      $   .and.(iorptl(i).ge.1.and.istptl(iorptl(i)).eq.41))
0556      $   .or.istptl(i).eq.20.or.istptl(i).eq.21
0557      $   .or.(istptl(i).eq.0.and.i.le.nptlpt))then
0558           call utlob4(1,p1(1),p1(2),p1(3),p1(4),p1(5)
0559      $         ,pptl(1,i),pptl(2,i),pptl(3,i),pptl(4,i))
0560         endif
0561       enddo
0562 
0563 
0564       if(ish.ge.5)write(ifch,'(a)')'------------------'
0565 
0566 c     rescale momenta in rest frame
0567 c     -----------------------------
0568 
0569       scal=1.
0570       diff=0.
0571       do ipass=1,1000
0572         sum=0.
0573         sum3=0.
0574         ndif=0
0575         do  j=1,nptl
0576         if((istptl(j).le.11
0577      $   .and.(iorptl(j).ge.1.and.istptl(iorptl(j)).eq.41))
0578      $   .or.istptl(j).eq.20.or.istptl(j).eq.21
0579      $   .or.(istptl(j).eq.0.and.j.le.nptlpt))then
0580             if(j.gt.nptlpt)then
0581               ndif=ndif+1
0582               pptl(3,j)=scal*(pptl(3,j)-diff)
0583               pptl(4,j)=sqrt(pptl(1,j)**2+pptl(2,j)**2+pptl(3,j)**2
0584      *           +pptl(5,j)**2)
0585             endif
0586             sum=sum+pptl(4,j)
0587             sum3=sum3+pptl(3,j)
0588           endif
0589         enddo
0590 
0591         diff=sum3/real(ndif)
0592         scal=real(esoll)/sum
0593         if(ish.ge.6)write(ifch,*)'ipass,scal,diff,e,esoll,pz,ndif:'
0594      $       ,ipass,scal,diff,sum,esoll,sum3,ndif
0595         if(abs(scal-1.).le.errlim.and.abs(diff).lt.10.*errlim) goto300
0596       enddo
0597       if(ish.ge.1)then
0598         call utmsg('hresph')
0599         write(ifch,*)'*****  scal=',scal,diff
0600         call utmsgf
0601       endif
0602 
0603 
0604 c     trafo
0605 c     -----
0606  300  continue
0607 c      do i=nptlpt+1,nptl
0608       do i=1,nptl
0609         if((istptl(i).le.11
0610      $   .and.(iorptl(i).ge.1.and.istptl(iorptl(i)).eq.41))
0611      $   .or.istptl(i).eq.20.or.istptl(i).eq.21
0612      $   .or.(istptl(i).eq.0.and.i.le.nptlpt))then
0613           call utlob4(-1,p1(1),p1(2),p1(3),p1(4),p1(5)
0614      $         ,pptl(1,i),pptl(2,i),pptl(3,i),pptl(4,i))
0615         endif
0616         if(i.le.nptlpt)then
0617           do j=1,5
0618             pptl(j,i)=p0(j,i)
0619           enddo
0620         endif
0621       enddo
0622 
0623  9999 call utprix('utrsph',ish,ishini,4)
0624 
0625       end
0626 
0627 cc-----------------------------------------------------------------------
0628 c      double precision function dddlog(xxx)
0629 cc-----------------------------------------------------------------------
0630 c      double precision xxx
0631 c      dddlog=-1d50
0632 c      if(xxx.gt.0d0)dddlog=dlog(xxx)
0633 c      end
0634 c
0635 ccc-----------------------------------------------------------------------
0636 c      subroutine randfl(jc,iqa0,iflav,ic,isame)
0637 cc-----------------------------------------------------------------------
0638 cc     returns random flavour ic(2) (iqa0=1:quark,2:antiquark,11:diquark)
0639 cc-----------------------------------------------------------------------
0640 c      include 'epos.inc'
0641 c      real probab(nflav),probsu(nflav+1)
0642 c      integer jc(nflav,2),jc0(nflav,2),ic(2)
0643 c      if(ish.ge.6)then
0644 c      write(ifch,*)('-',i=1,10)
0645 c     *,' entry sr randfl ',('-',i=1,30)
0646 c      write(ifch,*)'iqa0:',iqa0
0647 c      write(ifch,*)'jc:'
0648 c      write(ifch,*)jc
0649 c      endif
0650 c      iflav=0
0651 c      ic(1)=0
0652 c      ic(2)=0
0653 c      do 10 n=1,nflav
0654 c      do 10 i=1,2
0655 c10    jc0(n,i)=0
0656 c      iqa1=iqa0*10
0657 c9999  iqa1=iqa1/10
0658 c      if(iqa1.eq.0)goto9998
0659 c      iqa=mod(iqa1,10)
0660 c      su=0
0661 c      do 20 i=1,nflav
0662 c      probab(i)=jc(i,iqa)-jc0(i,iqa)
0663 c      if(isame.eq.1)probab(i)=probab(i)*(jc(i,3-iqa)-jc0(i,3-iqa))
0664 c20    su=su+probab(i)
0665 c      if(su.lt..5)then
0666 c      iflav=0
0667 c      ic(1)=0
0668 c      ic(2)=0
0669 c      goto9998
0670 c      endif
0671 c      probsu(1)=0.
0672 c      do 30 i=1,nflav
0673 c      probsu(i+1)=probsu(i)+probab(i)/su
0674 c      if(probsu(i+1)-probsu(i).lt.1e-5)probsu(i+1)=probsu(i)
0675 c30    continue
0676 c      r=rangen()*probsu(nflav+1)
0677 c      do 50 i=1,nflav
0678 c      if(probsu(i).le.r.and.r.lt.probsu(i+1))iflav=i
0679 c50    continue
0680 c      jc0(iflav,iqa)=jc0(iflav,iqa)+1
0681 c      if(isame.eq.1)jc0(iflav,3-iqa)=jc0(iflav,3-iqa)+1
0682 c      call idenco(jc0,ic,ireten)
0683 c      if(ireten.eq.1)call utstop('randfl: idenco ret code = 1&')
0684 c      if(ish.ge.6)then
0685 c      write(ifch,*)'probab:'
0686 c      write(ifch,*)probab
0687 c      write(ifch,*)'probsu:'
0688 c      write(ifch,*)probsu
0689 c      write(ifch,*)'ran#:',r,'   flav:',iflav
0690 c      endif
0691 c      goto9999
0692 c9998  continue
0693 c      if(ish.ge.6)write(ifch,*)('-',i=1,30)
0694 c     *,' exit sr randfl ',('-',i=1,10)
0695 c      return
0696 c      end
0697 c
0698 c
0699 cc-----------------------------------------------------------------------
0700 c      subroutine ranhvy(x,eps)
0701 cc-----------------------------------------------------------------------
0702 cc     generates x for heavy particle fragmentation according to
0703 cc     the peterson form
0704 cc          d(x)=1/(x*(1-1/x-eps/(1-x))**2)
0705 cc              =d0(x)*d1(x)*d2(x)
0706 cc          d0(x)=(1-x)**2/((1-x)**2+eps)**2
0707 cc          d1(x)=x
0708 cc          d2(x)=(((1-x)**2+eps)/((1-x)**2+eps*x))**2
0709 cc     using x=1-y**pow
0710 cc     generates flat in x if eps>1.
0711 cc-----------------------------------------------------------------------
0712 c      data aln4/1.3863/
0713 c      if(eps.lt.1.) then
0714 c        pow=alog((3.+eps)/eps)/aln4
0715 c        ymx=(eps*(3.*pow-1.)/(pow+1.))**(.5/pow)
0716 c        zmx=1-ymx**pow
0717 c        d0mx=(1-zmx)**2/((1.-zmx)**2+eps)**2*pow*ymx**(pow-1.)
0718 c        d2mx=2./(2.-sqrt(eps))
0719 c      else
0720 c        pow=1.
0721 c        zmx=0.
0722 c        d0mx=(1.-zmx)**2/((1.-zmx)**2+eps)**2
0723 c        d2mx=1.+eps
0724 c      endif
0725 cc
0726 cc          generate z according to (1-z)**2/((1-z)**2+eps*z)**2
0727 c1     continue
0728 c      y=rangen()
0729 c      z=1.-y**pow
0730 cc
0731 c      d0z=(1.-z)**2/((1.-z)**2+eps)**2*pow*y**(pow-1.)
0732 c      if(d0z.lt.rangen()*d0mx) goto1
0733 cc
0734 cc          check remaining factors
0735 c      d1=z
0736 c      d2=(((1.-z)**2+eps)/((1.-z)**2+eps*z))**2
0737 c      if(d1*d2.lt.rangen()*d2mx) goto1
0738 cc
0739 cc          good x
0740 c      x=z
0741 c      return
0742 c      end
0743 c
0744 *-- Author :    D. HECK IK FZK KARLSRUHE       27/04/1994
0745 C=======================================================================
0746 
0747       SUBROUTINE EPOVAPOR( MAPRO,INEW,JFIN,ITYP,PFRX,PFRY,PFRZ )
0748 
0749 C-----------------------------------------------------------------------
0750 C  (E)VAPOR(ATION OF NUCLEONS AND ALPHA PARTICLES FROM FRAGMENT)
0751 C
0752 C  TREATES THE REMAINING UNFRAGMENTED NUCLEUS
0753 C  EVAPORATION FOLLOWING CAMPI APPROXIMATION.
0754 C  SEE: X. CAMPI AND J. HUEFNER, PHYS.REV. C24 (1981) 2199
0755 C  AND  J.J. GAIMARD, THESE UNIVERSITE PARIS 7, (1990)
0756 C  THIS SUBROUTINE IS CALLED FROM SDPM, DPMJST, NSTORE, AND VSTORE.
0757 C  ARGUMENTS INPUT:
0758 C   MAPRO        = NUMBER OF NUCLEONS OF PROJECTILE
0759 C   INEW         = PARTICLE TYPE OF SPECTATOR FRAGMENT
0760 C  ARGUMENTS OUTPUT:
0761 C   JFIN         = NUMBER OF FRAGMENTS
0762 C   ITYP(1:JFIN) = NATURE (PARTICLE CODE) OF FRAGMENTS
0763 C   PFRX(1:JFIN) = TRANSVERSE MOMENTUM OF FRAGMENTS IN X-DIRECTION
0764 C   PFRY(1:JFIN) = TRANSVERSE MOMENTUM OF FRAGMENTS IN Y-DIRECTION
0765 C   PFRZ(1:JFIN) = LONGITUDINAL MOMENTUM OF FRAGMENTS IN Z-DIRECTION
0766 C
0767 C  FROM CORSIKA AND ADAPTED BY T. PIEROG TO INCLUDE LONG MOMENTUM AND
0768 C  MORE REALISTIC FRAGMENTS
0769 C-----------------------------------------------------------------------
0770 
0771       IMPLICIT NONE
0772       include 'epos.inc'
0773       common/eporansto2/irndmseq
0774       integer irndmseq
0775 
0776       DOUBLE PRECISION PFR(mamxx),PFRX(mamxx),PFRY(mamxx),PFRZ(mamxx)
0777      *                ,RD(2*mamxx),SPFRY,SPFRZ,drangen
0778       DOUBLE PRECISION AFIN,AGLH,APRF,BGLH,EEX,PHIFR,RANNORM,SPFRX
0779       INTEGER          ITYP(mamxx),IARM,INEW,ITYPRM,INRM,IS,IZRM,JC,lseq
0780      *                ,JFIN,K,L,LS,MAPRO,MF,NFIN,NINTA,NNUC,NPRF,NNSTEP
0781       SAVE
0782       EXTERNAL         RANNORM
0783 C-----------------------------------------------------------------------
0784 
0785       IF(ish.ge.7)WRITE(ifch,*)'EPOVAPOR : MAPRO,INEW=',MAPRO,INEW
0786 
0787       lseq   = irndmseq
0788       ITYPRM = INEW
0789       NPRF   = INEW/100
0790       NINTA  = MAPRO - NPRF
0791       IF ( NINTA .EQ. 0 ) THEN
0792 C  NO NUCLEON HAS INTERACTED
0793         JFIN    = 1
0794         PFRX(1)  = 0.D0
0795         PFRY(1)  = 0.D0
0796         PFRZ(1)  = 0.D0
0797         ITYP(1) = -ITYPRM
0798         IF(ish.ge.7)WRITE(ifch,*) 'EPOVAPOR : JFIN,NINTA=',JFIN,NINTA
0799         RETURN
0800       ENDIF
0801 
0802 C  EXCITATION ENERGY EEX OF PREFRAGMENT
0803 C  SEE: J.J. GAIMARD, THESE UNIVERSITE PARIS 7, (1990), CHPT. 4.2
0804       EEX = 0.D0
0805       CALL RMMARD( RD,2*NINTA,lseq )
0806       DO  L = 1, NINTA
0807         IF ( RD(NINTA+L) .LT. RD(L) ) RD(L) = 1.D0 - RD(L)
0808         EEX = EEX + RD(L)
0809       ENDDO
0810 C  DEPTH OF WOODS-SAXON POTENTIAL TO FERMI SURFACE IS 0.040 GEV
0811       IF(ish.ge.7)WRITE(ifch,*)'EPOVAPOR : EEX=',SNGL(EEX*0.04D0),
0812      &                                            ' GEV'
0813 C  EVAPORATION: EACH EVAPORATION STEP NEEDS ABOUT 0.020 GEV, THEREFORE
0814 C  NNSTEP IS EEX * 0.04/0.02 = EEX * 2.
0815       NNSTEP = INT( EEX*2.D0 )
0816 
0817       IF ( NNSTEP .LE. 0 ) THEN
0818 C  EXCITATION ENERGY TOO SMALL, NO EVAPORATION
0819         JFIN    = 1
0820         PFRX(1)  = 0.D0
0821         PFRY(1)  = 0.D0
0822         PFRZ(1)  = 0.D0
0823         ITYP(1) = -ITYPRM
0824         IF(ish.ge.7)WRITE(ifch,*) 'EPOVAPOR : JFIN,EEX=',JFIN,SNGL(EEX)
0825         RETURN
0826       ENDIF
0827 
0828 C  AFIN IS ATOMIC NUMBER OF FINAL NUCLEUS
0829       APRF = DBLE(NPRF)
0830       AFIN = APRF - 1.6D0 * DBLE(NNSTEP)
0831       NFIN = MAX( 0, INT( AFIN+0.5D0 ) )
0832 C  CORRESPONDS TO DEFINITION; FRAGMENTATION-EVAPORATION
0833 C  CONVOLUTION EMU07 /MODEL ABRASION EVAPORATION (JNC FZK APRIL 94)
0834 C  NNUC IS NUMBER OF EVAPORATING NUCLEONS
0835       NNUC = NPRF - NFIN
0836       IF(ish.ge.7)WRITE(ifch,*) 'EPOVAPOR : NFIN,NNUC=',NFIN,NNUC
0837       JC   = 0
0838 
0839       IF     ( NNUC .LE. 0 ) THEN
0840 C  NO EVAPORATION
0841         JFIN    = 1
0842         PFRX(1)  = 0.D0
0843         PFRY(1)  = 0.D0
0844         PFRZ(1)  = 0.D0
0845         ITYP(1) = -ITYPRM
0846         RETURN
0847 
0848       ELSEIF ( NNUC .GE. 4 ) THEN
0849 C  EVAPORATION WITH FORMATION OF ALPHA PARTICLES POSSIBLE
0850 C  IARM, IZRM, INRM ARE NUMBER OF NUCLEONS, PROTONS, NEUTRONS OF
0851 C  REMAINDER
0852         DO  LS = 1, NNSTEP
0853           IARM = ITYPRM/100
0854           IF ( IARM .LE. 0 ) GOTO 100
0855           IZRM = MOD(ITYPRM,100)
0856           INRM = IARM - IZRM
0857           JC   = JC + 1
0858           CALL RMMARD( RD,2,lseq )
0859           IF ( RD(1) .LT. 0.2D0  .AND.  IZRM .GE. 2
0860      *                           .AND.  INRM .GE. 2 ) THEN
0861             ITYP(JC) = -402          !alpha
0862             NNUC     = NNUC - 4
0863             ITYPRM   = ITYPRM - 402
0864           ELSE
0865             IF ( IZRM .EQ. 1 .AND. INRM .GT. IZRM ) THEN
0866               ITYP(JC) = 1220
0867               ITYPRM   = ITYPRM - 100              
0868             ELSEIF ( INRM .EQ. 1 .AND. IZRM .GT. INRM ) THEN
0869               ITYP(JC) = 1120
0870               ITYPRM   = ITYPRM - 101              
0871             ELSEIF(RD(2)*(IZRM+INRM).LT.IZRM.AND.IZRM.GE.INRM)THEN
0872               ITYP(JC) = 1120
0873               ITYPRM   = ITYPRM - 101
0874             ELSE
0875               ITYP(JC) = 1220
0876               ITYPRM   = ITYPRM - 100
0877             ENDIF
0878             NNUC = NNUC - 1
0879           ENDIF
0880           IF ( NNUC .LE. 0 ) GOTO 50
0881         ENDDO
0882       ENDIF
0883 
0884       IF ( NNUC .LT. 4 ) THEN
0885 C  EVAPORATION WITHOUT FORMATION OF ALPHA PARTICLES
0886         CALL RMMARD( RD,NNUC,lseq )
0887         DO  IS = 1, NNUC
0888           IARM = ITYPRM/100
0889           IF ( IARM .LE. 0 ) GOTO 100
0890           IZRM = MOD(ITYPRM,100)
0891           INRM = IARM - IZRM
0892           JC   = JC + 1
0893           IF ( IZRM .EQ. 1 .AND. INRM .GT. IZRM ) THEN
0894             ITYP(JC) = 1220
0895             ITYPRM   = ITYPRM - 100              
0896           ELSEIF ( INRM .EQ. 1 .AND. IZRM .GT. IZRM ) THEN
0897             ITYP(JC) = 1120
0898             ITYPRM   = ITYPRM - 101              
0899           ELSEIF ( RD(IS)*IARM .LT. IZRM .AND. IZRM .GE. INRM ) THEN
0900             ITYP(JC) = 1120
0901             ITYPRM   = ITYPRM - 101
0902           ELSE
0903             ITYP(JC) = 1220
0904             ITYPRM   = ITYPRM - 100
0905           ENDIF
0906         ENDDO
0907       ENDIF
0908 
0909  50   CONTINUE
0910       IARM = ITYPRM/100
0911       IF ( IARM .LE. 0 ) GOTO 100
0912       IZRM = MOD(ITYPRM,100)
0913       INRM = IARM - IZRM
0914       JC = JC + 1
0915       IF     ( IARM .EQ. 5 ) THEN     !EXCLUDED
0916         IF ( IZRM .GE. INRM ) THEN
0917           ITYP(JC) = 1120
0918           ITYPRM   = ITYPRM - 101              
0919         ELSE
0920           ITYP(JC) = 1220
0921           ITYPRM   = ITYPRM - 100
0922         ENDIF
0923         JC = JC + 1
0924         ITYP(JC) = -ITYPRM
0925       ELSEIF     ( ITYPRM .GT. 200 ) THEN
0926         ITYP(JC) = -ITYPRM
0927       ELSEIF ( ITYPRM .EQ. 200 ) THEN
0928         ITYP(JC) = 1220
0929         JC = JC + 1
0930         ITYP(JC) = 1220
0931       ELSEIF ( ITYPRM .EQ. 101 ) THEN
0932         ITYP(JC) = 1120
0933       ELSEIF ( ITYPRM .EQ. 100 ) THEN
0934         ITYP(JC) = 1220
0935       ELSE
0936         JC = JC - 1
0937         IF ( ITYPRM .NE. 0 ) WRITE(*,*)
0938      *                  'EPOVAPOR : ILLEGAL PARTICLE ITYPRM =',ITYPRM
0939       ENDIF
0940 
0941   100 CONTINUE
0942       JFIN = JC
0943       IF(ish.ge.7)WRITE(ifch,*) 
0944      *   'EPOVAPOR :  NO        ITYP     PFR       PFL'
0945       IF     ( infragm .EQ. 2 ) THEN
0946 C  EVAPORATION WITH PT AFTER PARAMETRIZED JACEE DATA
0947         DO  MF = 1, JFIN
0948           IF(ITYP(MF).LT.0)THEN
0949             IARM=-ITYP(MF)/100
0950           ELSE
0951             IARM=1
0952           ENDIF
0953           PFR(MF) = RANNORM(0.088D0,0.044D0)
0954           PFRZ(MF)= (2*int(drangen(PFR(MF))+0.5d0)-1)
0955      &   *RANNORM(0.300D0/DBLE(IARM),0.100D0/SQRT(DBLE(IARM)))    !Fermi motion about 300 MeV
0956           IF(ish.ge.7)WRITE(ifch,*) MF,ITYP(MF),SNGL(PFR(MF))
0957      &                                         ,SNGL(PFRZ(MF))
0958         ENDDO
0959       ELSEIF ( infragm .EQ. 3 ) THEN
0960 C  EVAPORATION WITH PT AFTER GOLDHABER''S MODEL (PHYS.LETT.53B(1974)306)
0961         DO  MF = 1, JFIN
0962           K    = MAX( 1, -ITYP(MF)/100 )
0963           BGLH = K * (MAPRO - K) / DBLE(MAPRO-1)
0964 C  THE VALUE 0.103 [GEV] IS SIGMA(0)=P(FERMI)/SQRT(5.)
0965 *         AGLH = 0.103D0 * SQRT( BGLH )
0966 C  THE VALUE 0.090 [GEV] IS EXPERIMENTALLY DETERMINED SIGMA(0)
0967           AGLH = 0.090D0 * SQRT( BGLH )
0968           PFR(MF) = RANNORM(0.D0,AGLH)
0969           PFRZ(MF)= RANNORM(0.000D0,0.500D0)    !from pAg at 100 GeV
0970           IF(ish.ge.7)WRITE(ifch,*) MF,ITYP(MF),SNGL(PFR(MF))
0971      &                                         ,SNGL(PFRZ(MF))
0972         ENDDO
0973       ELSE
0974 C  EVAPORATION WITHOUT TRANSVERSE MOMENTUM
0975         DO  MF = 1, JFIN
0976           PFR(MF) = 0.D0
0977           PFRZ(MF)= 0.D0
0978           IF(ish.ge.7)WRITE(ifch,*) MF,ITYP(MF),SNGL(PFR(MF))
0979      &                                         ,SNGL(PFRZ(MF))
0980         ENDDO
0981       ENDIF
0982 C  CALCULATE RESIDUAL TRANSVERSE MOMENTUM
0983       SPFRX = 0.D0
0984       SPFRY = 0.D0
0985       SPFRZ = 0.D0
0986       CALL RMMARD( RD,JFIN,lseq )
0987       DO  MF = 1, JFIN
0988         PHIFR = PI * RD(MF)
0989         PFRX(MF) = PFR(MF) * COS( PHIFR )
0990         PFRY(MF) = PFR(MF) * SIN( PHIFR )
0991         SPFRY = SPFRY + PFRY(MF)
0992         SPFRX = SPFRX + PFRX(MF)
0993         SPFRZ = SPFRZ + PFRZ(MF)
0994       ENDDO
0995 C  CORRECT ALL TRANSVERSE MOMENTA FOR MOMENTUM CONSERVATION
0996       SPFRX = SPFRX / JFIN
0997       SPFRY = SPFRY / JFIN
0998       SPFRZ = SPFRZ / JFIN
0999       DO  MF = 1, JFIN
1000         PFRX(MF) = PFRX(MF) - SPFRX
1001         PFRY(MF) = PFRY(MF) - SPFRY
1002         PFRZ(MF) = PFRZ(MF) - SPFRZ
1003       ENDDO
1004 
1005       IF(ish.ge.7)WRITE(ifch,*) 'EPOVAPOR : NINTA,JFIN=',NINTA,JFIN
1006 
1007       RETURN
1008       END
1009 
1010 *-- Author :    The CORSIKA development group   21/04/1994
1011 C=======================================================================
1012 
1013       DOUBLE PRECISION FUNCTION RANNORM( A,B )
1014 
1015 C-----------------------------------------------------------------------
1016 C  RAN(DOM NUMBER) NOR(MALLY DISTRIBUTED)
1017 C
1018 C  GENERATES NORMAL DISTRIBUTED RANDOM NUMBER
1019 C  DELIVERS 2 UNCORRELATED RANDOM NUMBERS,
1020 C  THEREFORE RANDOM CALLS ARE ONLY NECESSARY EVERY SECOND TIME.
1021 C  REFERENCE : NUMERICAL RECIPES, W.H. PRESS ET AL.,
1022 C              CAMBRIDGE UNIVERSITY PRESS, 1992  ISBN 0 521 43064 X
1023 C  ARGUMENTS:
1024 C   A      = MEAN VALUE
1025 C   B      = STANDARD DEVIATION
1026 C
1027 C  FROM CORSIKA
1028 C-----------------------------------------------------------------------
1029 
1030       IMPLICIT NONE
1031       double precision facrdm,u1rdm,u2rdm,drangen
1032       logical knordm
1033       data knordm/.true./
1034 
1035       DOUBLE PRECISION A,B,RR
1036       SAVE facrdm,u1rdm,u2rdm,knordm
1037 C-----------------------------------------------------------------------
1038 
1039       IF ( KNORdm ) THEN
1040   1     CONTINUE
1041         U1rdm = 2.D0*drangen(a) - 1.D0
1042         U2rdm = 2.D0*drangen(b) - 1.D0
1043         RR = U1rdm**2 + U2rdm**2
1044         IF ( RR .GE. 1.D0  .OR.  RR .EQ. 0.D0 ) GOTO 1
1045         FACrdm = SQRT( (-2.D0) * LOG(RR) / RR )
1046 
1047         RANNORM = FACrdm * U1rdm * B + A
1048         KNORdm   = .FALSE.
1049       ELSE
1050         RANNORM = FACrdm * U2rdm * B + A
1051         KNORdm   = .TRUE.
1052       ENDIF
1053 
1054       RETURN
1055       END
1056 
1057 c-----------------------------------------------------------------------
1058       function ransig()
1059 c-----------------------------------------------------------------------
1060 c     returns randomly +1 or -1
1061 c-----------------------------------------------------------------------
1062       ransig=1
1063       if(rangen().gt.0.5)ransig=-1
1064       return
1065       end
1066 
1067 cc-----------------------------------------------------------------------
1068 c      function ranxq(n,x,q,xmin)
1069 cc-----------------------------------------------------------------------
1070 cc     returns random number according to x(i) q(i) with x>=xmin
1071 cc-----------------------------------------------------------------------
1072 c      include 'epos.inc'
1073 c      real x(n),q(n)
1074 c      imin=1
1075 c      if(xmin.eq.0.)goto3
1076 c      i1=1
1077 c      i2=n
1078 c1     i=i1+(i2-i1)/2
1079 c      if(x(i).lt.xmin)then
1080 c      i1=i
1081 c      elseif(x(i).gt.xmin)then
1082 c      i2=i
1083 c      else
1084 c      imin=i
1085 c      goto3
1086 c      endif
1087 c      if(i2-i1.gt.1)goto1
1088 c      imin=i2
1089 c3     continue
1090 c      if(q(imin).gt.q(n)*.9999)then
1091 c      ranxq=xmin
1092 c      goto4
1093 c      endif
1094 c      qran=q(imin)+rangen()*(q(n)-q(imin))
1095 c      ranxq=utinvt(n,x,q,qran)
1096 c4     continue
1097 c
1098 c      if(ranxq.lt.xmin)then
1099 c      call utmsg('ranxq ')
1100 c      write(ifch,*)'*****  ranxq=',ranxq,' <       xmin=',xmin
1101 c      write(ifch,*)'q(imin) q q(n):',q(imin),qran,q(n)
1102 c      write(ifch,*)'x(imin) x x(n):',x(imin),ranxq,x(n)
1103 c      call utmsgf
1104 c      ranxq=xmin
1105 c      endif
1106 c
1107 c      return
1108 c      end
1109 c
1110 cc  ***** end r-routines
1111 cc  ***** beg s-routines
1112 c
1113 cc-----------------------------------------------------------------------
1114 c      function sbet(z,w)
1115 cc-----------------------------------------------------------------------
1116 c      sbet=utgam1(z)*utgam1(w)/utgam1(z+w)
1117 c      return
1118 c      end
1119 c
1120 cc-----------------------------------------------------------------------
1121 c      function smass(a,y,z)
1122 cc-----------------------------------------------------------------------
1123 cc     returns droplet mass (in gev) (per droplet, not (!) per nucleon)
1124 cc     according to berger/jaffe mass formula, prc35(1987)213 eq.2.31,
1125 cc     see also c. dover, BNL-46322, intersections-meeting, tucson, 91.
1126 cc     a: massnr, y: hypercharge, z: charge,
1127 cc-----------------------------------------------------------------------
1128 c      common/cmass/thet,epsi,as,ac,dy,dz,ym,cz,zm,sigma,rzero
1129 c      ymin=ym*a
1130 c      zmin=cz/(dz/a+zm/a**(1./3.))
1131 c      smass=epsi*a+as*a**(2./3.)+(ac/a**(1./3.)+dz/a/2.)*(z-zmin)**2
1132 c     *+dy/a/2.*(y-ymin)**2
1133 c      return
1134 c      end
1135 c
1136 cc-----------------------------------------------------------------------
1137 c      subroutine smassi(theta)
1138 cc-----------------------------------------------------------------------
1139 cc     initialization for smass.
1140 cc     calculates parameters for berger/jaffe mass formula
1141 cc     (prc35(1987)213 eq.2.31, see also c. dover, BNL-46322).
1142 cc     theta: parameter that determines all parameters in mass formula.
1143 cc-----------------------------------------------------------------------
1144 c      common/cmass/thet,epsi,as,ac,dy,dz,ym,cz,zm,sigma,rzero
1145 c      thet=theta
1146 c
1147 c      astr=.150
1148 c      pi=3.14159
1149 c      alp=1./137.
1150 c
1151 c      co=cos(theta)
1152 c      si=sin(theta)
1153 c      bet=(1+co**3)/2.
1154 c      rzero=si/astr/(  2./3./pi*(1+co**3)  )**(1./3.)
1155 cctp060829      cs=astr/si
1156 c      cz=-astr/si*(  (  .5*(1+co**3)  )**(1./3.)-1  )
1157 c      sigma=6./8./pi*(astr/si)**3*(co**2/6.-si**2*(1-si)/3.-
1158 c     *1./3./pi*(pi/2.-theta-sin(2*theta)+si**3*alog((1+co)/si)))
1159 c
1160 c      epsi=astr*((.5*(1+co**3))**(1./3.)+2)/si
1161 c      as=4*pi*sigma*rzero**2
1162 c      ac=3./5.*alp/rzero
1163 c      dz=astr/si*bet**(1./3.)*co**2*
1164 c     *(co**4*(1+bet**(2./3.))+(1+bet)**2)/
1165 c     *(  (2*co**2+bet**(1./3.))*(co**4*(1+bet**(2./3.))+(1+bet)**2)-
1166 c     *(co**4+bet**(1./3.)*(1+bet))*((2*bet**(2./3.)-1)*co**2+1+bet)  )
1167 c      dy=astr/6.*(1+co**3)**3/si*
1168 c     *(  1+(1+co)/(4*(1+co**3))**(2./3.)  )/
1169 c     *(co**6+co+co*(.5*(1+co**3))**(4./3.))
1170 c      zm=6*alp/(5*rzero)
1171 c      ym=(1-co**3)/(1+co**3)
1172 c
1173 c      return
1174 c      end
1175 c
1176 cc-----------------------------------------------------------------------
1177 c      subroutine smassp
1178 cc-----------------------------------------------------------------------
1179 cc     prints smass.
1180 cc-----------------------------------------------------------------------
1181 c      include 'epos.inc'
1182 c      common/cmass/thet,epsi,as,ac,dy,dz,ym,cz,zm,sigma,rzero
1183 c      real eng(14),ymi(14),zmi(14)
1184 c      pi=3.14159
1185 c      write(ifch,*)'parameters of mass formula:'
1186 c      write(ifch,*)'theta=',thet,'   epsi=',epsi
1187 c      write(ifch,*)'as=',as,'   ac=',ac
1188 c      write(ifch,*)'dy=',dy,'   dz=',dz
1189 c      write(ifch,*)'ym=',ym
1190 c      write(ifch,*)'cz dz zm=',cz,dz,zm
1191 c      write(ifch,*)'sigma**1/3=',sigma**(1./3.),'   rzero=',rzero
1192 c      write(ifch,*)'mass:'
1193 c      write(ifch,5000)(j,j=1,14)
1194 c5000  format(5x,'a:',14i5)
1195 c      do 4 j=1,14
1196 c      a=j
1197 c      ymi(j)=ym*a
1198 c4     zmi(j)=cz/(dz/a+zm/a**(1./3.))
1199 c      write(ifch,5002)(ymi(j),j=1,14)
1200 c5002  format(1x,'ymin: ',14f5.2)
1201 c      write(ifch,5003)(zmi(j),j=1,14)
1202 c5003  format(1x,'zmin: ',14f5.2)
1203 c      do 2 i=1,15
1204 c      ns=11-i
1205 c      do 3 j=1,14
1206 c      a=j
1207 c      y=a-ns
1208 c      z=0.
1209 c3     eng(j)=smass(a,y,z)/a
1210 c      write(ifch,5001)ns,(eng(j),j=1,14)
1211 c5001  format(1x,'s=',i2,2x,14f5.2)
1212 c2     continue
1213 c      write(ifch,*)'mass-mass(free):'
1214 c      write(ifch,5000)(j,j=1,14)
1215 c      do 5 i=1,15
1216 c      ns=11-i
1217 c      do 6 j=1,14
1218 c      a=j
1219 c      y=a-ns
1220 c      z=0.
1221 c      call smassu(a,y,z,ku,kd,ks,kc)
1222 c6     eng(j)=(smass(a,y,z)-utamnu(ku,kd,ks,kc,0,0,3))/a
1223 c      write(ifch,5001)ns,(eng(j),j=1,14)
1224 c5     continue
1225 c
1226 c      stop
1227 c      end
1228 c
1229 cc-----------------------------------------------------------------------
1230 c      subroutine smasst(kux,kdx,ksx,kcx,a,y,z)
1231 cc-----------------------------------------------------------------------
1232 cc     input: kux,kdx,ksx,kcx = net quark numbers (for u,d,s,c quarks).
1233 cc     output: massnr a, hypercharge y and charge z.
1234 cc-----------------------------------------------------------------------
1235 c      sg=1
1236 c      if(kux+kdx+ksx+kcx.lt.0.)sg=-1
1237 c      ku=sg*kux
1238 c      kd=sg*kdx
1239 c      ks=sg*ksx
1240 c      kc=sg*kcx
1241 c      k=ku+kd+ks+kc
1242 c      if(mod(k,3).ne.0)stop'noninteger baryon number'
1243 c      a=k/3
1244 c      y=a-ks
1245 c      nz=2*ku-kd-ks+2*kc
1246 c      if(mod(nz,3).ne.0)stop'noninteger charge'
1247 c      z=nz/3
1248 c      return
1249 c      end
1250 c
1251 cc-----------------------------------------------------------------------
1252 c      subroutine smassu(ax,yx,zx,ku,kd,ks,kc)
1253 cc-----------------------------------------------------------------------
1254 cc     input: massnr ax, hypercharge yx and charge zx.
1255 cc     output: ku,kd,ks,kc = net quark numbers (for u,d,s,c quarks).
1256 cc-----------------------------------------------------------------------
1257 c      sg=1
1258 c      if(ax.lt.0.)sg=-1
1259 c      a=sg*ax
1260 c      y=sg*yx
1261 c      z=sg*zx
1262 c      ku=nint(a+z)
1263 c      kd=nint(a-z+y)
1264 c      ks=nint(a-y)
1265 c      kc=0
1266 c      return
1267 c      end
1268 c
1269 cc-----------------------------------------------------------------------
1270 c      function spoc(a,b,c,d,x)
1271 cc-----------------------------------------------------------------------
1272 cc     power fctn with cutoff
1273 cc-----------------------------------------------------------------------
1274 c      spoc=0
1275 c      if(a.eq.0..and.b.eq.0.)return
1276 c      spoc =a+b*x**c
1277 c      spoc0=a+b*d**c
1278 c      spoc=amin1(spoc,spoc0)
1279 c      spoc=amax1(0.,spoc)
1280 c      return
1281 c      end
1282 c
1283 c-----------------------------------------------------------------------
1284       function utacos(x)
1285 c-----------------------------------------------------------------------
1286 c     returns acos(x) for -1 <= x <= 1 , acos(+-1) else
1287 c-----------------------------------------------------------------------
1288       include 'epos.inc'
1289       argum=x
1290       if(x.lt.-1.)then
1291       if(ish.ge.1)then
1292       call utmsg('utacos')
1293       write(ifch,*)'*****  argum = ',argum,' set -1'
1294       call utmsgf
1295       endif
1296       argum=-1.
1297       elseif(x.gt.1.)then
1298       if(ish.ge.1)then
1299       call utmsg('utacos')
1300       write(ifch,*)'*****  argum = ',argum,' set 1'
1301       call utmsgf
1302       endif
1303       argum=1.
1304       endif
1305       utacos=acos(argum)
1306       return
1307       end
1308 
1309 c----------------------------------------------------------------------
1310       function utamnu(keux,kedx,kesx,kecx,kebx,ketx,modus)
1311 c----------------------------------------------------------------------
1312 c     returns min mass of droplet with given u,d,s,c content
1313 c     keux: net u quark number
1314 c     kedx: net d quark number
1315 c     kesx: net s quark number
1316 c     kecx: net c quark number
1317 c     kebx: net b quark number
1318 c     ketx: net t quark number
1319 c     modus: 4=two lowest multiplets; 5=lowest multiplet
1320 c----------------------------------------------------------------------
1321       common/files/ifop,ifmt,ifch,ifcx,ifhi,ifdt,ifcp,ifdr
1322       common/csjcga/amnull,asuha(7)
1323       common/drop4/asuhax(7),asuhay(7)
1324 
1325       if(modus.lt.4.or.modus.gt.5)stop'UTAMNU: not supported'
1326 c 1    format(' flavours:',6i5 )
1327 c 100  format(' flavours+mass:',6i5,f8.2 )
1328 c      write(ifch,1)keux,kedx,kesx,kecx,kebx,ketx
1329 
1330       amnull=0.
1331 
1332       do i=1,7
1333       if(modus.eq.4)asuha(i)=asuhax(i)    !two lowest multiplets
1334       if(modus.eq.5)asuha(i)=asuhay(i)    !lowest multiplet
1335       enddo
1336 
1337       ke=iabs(keux+kedx+kesx+kecx+kebx+ketx)
1338 
1339       if(keux+kedx+kesx+kecx+kebx+ketx.ge.0)then
1340       keu=keux
1341       ked=kedx
1342       kes=kesx
1343       kec=kecx
1344       keb=kebx
1345       ket=ketx
1346       else
1347       keu=-keux
1348       ked=-kedx
1349       kes=-kesx
1350       kec=-kecx
1351       keb=-kebx
1352       ket=-ketx
1353       endif
1354 
1355 c      write(ifch,*)keu,ked,kes,kec,keb,ket
1356 
1357 c   removing top mesons  to remove t quarks or antiquarks
1358       if(ket.ne.0)then
1359 12    continue
1360       ii=sign(1,ket)
1361       ket=ket-ii
1362       if(ii*keu.le.ii*ked)then
1363       keu=keu+ii
1364       else
1365       ked=ked+ii
1366       endif
1367       amnull=amnull+200.    ! ???????
1368       if(ket.ne.0)goto12
1369       endif
1370 
1371 c   removing bottom mesons  to remove b quarks or antiquarks
1372       if(keb.ne.0)then
1373 11    continue
1374       ii=sign(1,keb)
1375       keb=keb-ii
1376       if(ii*keu.le.ii*ked)then
1377       keu=keu+ii
1378       else
1379       ked=ked+ii
1380       endif
1381       amnull=amnull+6. !5.28   ! (more than B-meson)
1382       if(keb.ne.0)goto11
1383       endif
1384 
1385 c   removing charm mesons  to remove c quarks or antiquarks
1386       if(kec.ne.0)then
1387 10    continue
1388       ii=sign(1,kec)
1389       kec=kec-ii
1390       if(keu*ii.le.ked*ii)then
1391       keu=keu+ii
1392       else
1393       ked=ked+ii
1394       endif
1395       amnull=amnull+2.2 !1.87  ! (more than D-meson)
1396       if(kec.ne.0)goto10
1397       endif
1398 
1399 c      write(ifch,100)keu,ked,kes,kec,keb,ket,amnull
1400 
1401 c   removing mesons to remove s antiquarks
1402 5     continue
1403       if(kes.lt.0)then
1404       amnull=amnull+asuha(6)
1405       if(keu.ge.ked)then
1406       keu=keu-1
1407       else
1408       ked=ked-1
1409       endif
1410       kes=kes+1
1411       goto5
1412       endif
1413 
1414 c   removing mesons to remove d antiquarks
1415 6     continue
1416       if(ked.lt.0)then
1417       if(keu.ge.kes)then
1418       amnull=amnull+asuha(5)
1419       keu=keu-1
1420       else
1421       amnull=amnull+asuha(6)
1422       kes=kes-1
1423       endif
1424       ked=ked+1
1425       goto6
1426       endif
1427 
1428 c   removing mesons to remove u antiquarks
1429 7     continue
1430       if(keu.lt.0)then
1431       if(ked.ge.kes)then
1432       amnull=amnull+asuha(5)
1433       ked=ked-1
1434       else
1435       amnull=amnull+asuha(6)
1436       kes=kes-1
1437       endif
1438       keu=keu+1
1439       goto7
1440       endif
1441 
1442 c      write(ifch,100)keu,ked,kes,kec,keb,ket,amnull
1443 c      print*,keu,ked,kes,kec,keb,ket,amnull
1444 
1445       if(keu+ked+kes+kec+keb+ket.ne.ke)
1446      *call utstop('utamnu: sum_kei /= ke&')
1447       keq=keu+ked
1448       keqx=keq
1449       amnux=0
1450 
1451 c   removing strange baryons
1452       i=4
1453 2     i=i-1
1454 3     continue
1455       if((4-i)*kes.gt.(i-1)*keq)then
1456       amnux=amnux+asuha(1+i)
1457       kes=kes-i
1458       keq=keq-3+i
1459       if(kes.lt.0)call utstop('utamnu: negative kes&')
1460       if(keq.lt.0)call utstop('utamnu: negative keq&')
1461       goto3
1462       endif
1463       if(i.gt.1)goto2
1464       if(keqx.gt.keq)then
1465       do 8 k=1,keqx-keq
1466       if(keu.ge.ked)then
1467       keu=keu-1
1468       else
1469       ked=ked-1
1470       endif
1471 8     continue
1472       endif
1473 
1474       if(keu+ked.ne.keq)call utstop('utamnu: keu+ked /= keq&')
1475 c      write(ifch,100)keu,ked,kes,kec,keb,ket,amnull+amnux
1476 c      print*,keu,ked,kes,kec,keb,ket,amnull+amnux
1477 
1478 c   removing nonstrange baryons
1479 9     continue
1480       if(keu.gt.2*ked)then
1481       amnux=amnux+asuha(7)
1482       keu=keu-3
1483       if(keu.lt.0)call utstop('utamnu: negative keu&')
1484       goto9
1485       endif
1486       if(ked.gt.2*keu)then
1487       amnux=amnux+asuha(7)
1488       ked=ked-3
1489       if(ked.lt.0)call utstop('utamnu: negative ked&')
1490       goto9
1491       endif
1492       keq=keu+ked
1493 
1494 c      write(ifch,100)keu,ked,kes,kec,keb,ket,amnull+amnux
1495 c      print*,keu,ked,kes,kec,keb,ket,amnull+amnux
1496 
1497       if(mod(keq,3).ne.0)call utstop('utamnu: mod(keq,3) /= 0&')
1498       amnux=amnux+asuha(1)*keq/3
1499 
1500 c      write(ifch,100)keu,ked,kes,kec,keb,ket,amnull+amnux
1501 c      print*,keu,ked,kes,kec,keb,ket,amnull+amnux
1502 
1503       amnull=amnull+amnux
1504 
1505       if(amnull.eq.0)amnull=asuha(5)
1506 
1507       utamnu=amnull
1508       return
1509       end
1510 
1511 c-----------------------------------------------------------------------
1512       function utamnx(jcp,jcm)
1513 c-----------------------------------------------------------------------
1514 c returns minimum mass for the decay of jcp---jcm (by calling utamnu).
1515 c-----------------------------------------------------------------------
1516       parameter (nflav=6)
1517       integer jcp(nflav,2),jcm(nflav,2)
1518 
1519       do i=1,nflav
1520       do j=1,2
1521       if(jcp(i,j).ne.0)goto1
1522       enddo
1523       enddo
1524       keu=jcm(1,1)-jcm(1,2)
1525       ked=jcm(2,1)-jcm(2,2)
1526       kes=jcm(3,1)-jcm(3,2)
1527       kec=jcm(4,1)-jcm(4,2)
1528       keb=jcm(5,1)-jcm(5,2)
1529       ket=jcm(6,1)-jcm(6,2)
1530       utamnx=utamnu(keu,ked,kes,kec,keb,ket,5)
1531       return
1532 1     continue
1533 
1534       do i=1,nflav
1535       do j=1,2
1536       if(jcm(i,j).ne.0)goto2
1537       enddo
1538       enddo
1539       keu=jcp(1,1)-jcp(1,2)
1540       ked=jcp(2,1)-jcp(2,2)
1541       kes=jcp(3,1)-jcp(3,2)
1542       kec=jcp(4,1)-jcp(4,2)
1543       keb=jcp(5,1)-jcp(5,2)
1544       ket=jcp(6,1)-jcp(6,2)
1545       utamnx=utamnu(keu,ked,kes,kec,keb,ket,5)
1546       return
1547 2     continue
1548 
1549       keu=jcp(1,1)-jcp(1,2)
1550       ked=jcp(2,1)-jcp(2,2)
1551       kes=jcp(3,1)-jcp(3,2)
1552       kec=jcp(4,1)-jcp(4,2)
1553       keb=jcp(5,1)-jcp(5,2)
1554       ket=jcp(6,1)-jcp(6,2)
1555       ke=keu+ked+kes+kec+keb+ket
1556       if(mod(ke+1,3).eq.0)then
1557         keu=keu+1
1558         amms1=utamnu(keu,ked,kes,kec,keb,ket,5)
1559         keu=keu-1
1560         ked=ked+1
1561         amms2=utamnu(keu,ked,kes,kec,keb,ket,5)
1562       elseif(mod(ke-1,3).eq.0)then
1563         keu=keu-1
1564         amms1=utamnu(keu,ked,kes,kec,keb,ket,5)
1565         keu=keu+1
1566         ked=ked-1
1567         amms2=utamnu(keu,ked,kes,kec,keb,ket,5)
1568       else
1569         amms1=0
1570         amms2=0
1571         amms3=0
1572         amms4=0
1573         call utstop('utamnx: no singlet possible (1)&')
1574       endif
1575       keu=jcm(1,1)-jcm(1,2)
1576       ked=jcm(2,1)-jcm(2,2)
1577       kes=jcm(3,1)-jcm(3,2)
1578       kec=jcm(4,1)-jcm(4,2)
1579       keb=jcm(5,1)-jcm(5,2)
1580       ket=jcm(6,1)-jcm(6,2)
1581       ke=keu+ked+kes+kec+keb+ket
1582       if(mod(ke+1,3).eq.0)then
1583         keu=keu+1
1584         amms3=utamnu(keu,ked,kes,kec,keb,ket,5)
1585         keu=keu-1
1586         ked=ked+1
1587         amms4=utamnu(keu,ked,kes,kec,keb,ket,5)
1588       elseif(mod(ke-1,3).eq.0)then
1589         keu=keu-1
1590         amms3=utamnu(keu,ked,kes,kec,keb,ket,5)
1591         keu=keu+1
1592         ked=ked-1
1593         amms4=utamnu(keu,ked,kes,kec,keb,ket,5)
1594       else
1595         call utstop('utamnx: no singlet possible (2)&')
1596       endif
1597       utamnx=min(amms1+amms3,amms2+amms4)
1598 c       print *,amms1,amms3,amms2,amms4,jcp,jcm
1599       return
1600       end
1601 
1602 
1603 
1604 cc-----------------------------------------------------------------------
1605 c      function utamny(jcp,jcm)
1606 cc-----------------------------------------------------------------------
1607 cc returns minimum mass of jcp+jcm (by calling utamnu).
1608 cc-----------------------------------------------------------------------
1609 c      parameter (nflav=6)
1610 c      integer jcp(nflav,2),jcm(nflav,2),jc(nflav,2)
1611 c      do 7 nf=1,nflav
1612 c      jc(nf,1)=jcp(nf,1)+jcm(nf,1)
1613 c7     jc(nf,2)=jcp(nf,2)+jcm(nf,2)
1614 c      keu=jc(1,1)-jc(1,2)
1615 c      ked=jc(2,1)-jc(2,2)
1616 c      kes=jc(3,1)-jc(3,2)
1617 c      kec=jc(4,1)-jc(4,2)
1618 c      keb=jc(5,1)-jc(5,2)
1619 c      ket=jc(6,1)-jc(6,2)
1620 c      utamny=utamnu(keu,ked,kes,kec,keb,ket,5)
1621 c      return
1622 c      end
1623 c
1624 c-----------------------------------------------------------------------
1625       function utamnz(jc,modus)
1626 c-----------------------------------------------------------------------
1627 c returns minimum mass of jc (by calling utamnu).
1628 c-----------------------------------------------------------------------
1629       parameter (nflav=6)
1630       integer jc(nflav,2)
1631       keu=jc(1,1)-jc(1,2)
1632       ked=jc(2,1)-jc(2,2)
1633       kes=jc(3,1)-jc(3,2)
1634       kec=jc(4,1)-jc(4,2)
1635       keb=jc(5,1)-jc(5,2)
1636       ket=jc(6,1)-jc(6,2)
1637       utamnz=utamnu(keu,ked,kes,kec,keb,ket,modus)
1638       return
1639       end
1640 
1641 c-----------------------------------------------------------------------
1642       subroutine utar(i1,i2,i3,x0,x1,x2,x3,xx)
1643 c-----------------------------------------------------------------------
1644 c     returns the array xx with xx(1)=x0 <= xx(i) <= xx(i3)=x3
1645 c-----------------------------------------------------------------------
1646       real xx(i3)
1647       do 1 i=1,i1-1
1648   1   xx(i)=x0+(i-1.)/(i1-1.)*(x1-x0)
1649       do 2 i=i1,i2-1
1650   2   xx(i)=x1+(i-i1*1.)/(i2-i1*1.)*(x2-x1)
1651       do 3 i=i2,i3
1652   3   xx(i)=x2+(i-i2*1.)/(i3-i2*1.)*(x3-x2)
1653       return
1654       end
1655 
1656 cc---------------------------------------------------------------------
1657 c      subroutine utaxis(i,j,a1,a2,a3)
1658 cc-----------------------------------------------------------------------
1659 cc     calculates the axis defined by the ptls i,j in the i,j cm system
1660 cc---------------------------------------------------------------------
1661 c      include 'epos.inc'
1662 c      double precision pi1,pi2,pi3,pi4,pj1,pj2,pj3,pj4,p1,p2,p3,p4,p5
1663 c     *,err,a
1664 c      a1=0
1665 c      a2=0
1666 c      a3=1
1667 c      pi1=dble(pptl(1,i))
1668 c      pi2=dble(pptl(2,i))
1669 c      pi3=dble(pptl(3,i))
1670 c      pi4=dble(pptl(4,i))
1671 c      pj1=dble(pptl(1,j))
1672 c      pj2=dble(pptl(2,j))
1673 c      pj3=dble(pptl(3,j))
1674 c      pj4=dble(pptl(4,j))
1675 c      p1=pi1+pj1
1676 c      p2=pi2+pj2
1677 c      p3=pi3+pj3
1678 c      p4=pi4+pj4
1679 c      p5=dsqrt(p4**2-p3**2-p2**2-p1**2)
1680 c      call utlob2(1,p1,p2,p3,p4,p5,pi1,pi2,pi3,pi4,50)
1681 c      call utlob2(1,p1,p2,p3,p4,p5,pj1,pj2,pj3,pj4,51)
1682 c           err=(pi1+pj1)**2+(pi2+pj2)**2+(pi3+pj3)**2
1683 c           if(err.gt.1d-3)then
1684 c      call utmsg('utaxis')
1685 c      write(ifch,*)'*****  err=',err
1686 c      write(ifch,*)'pi:',pi1,pi2,pi3,pi4
1687 c      write(ifch,*)'pj:',pj1,pj2,pj3,pj4
1688 c      call utmsgf
1689 c           endif
1690 c      a=dsqrt( (pj1-pi1)**2 + (pj2-pi2)**2 + (pj3-pi3)**2 )
1691 c      if(a.eq.0.d0)return
1692 c      a1=sngl((pi1-pj1)/a)
1693 c      a2=sngl((pi2-pj2)/a)
1694 c      a3=sngl((pi3-pj3)/a)
1695 c      return
1696 c      end
1697 c
1698 cc---------------------------------------------------------------------
1699 c      subroutine uthalf(i,j,zz1,zz2,iret)
1700 cc-----------------------------------------------------------------------
1701 cc     give equal energy (E_i+E_j)/2 to particle i+j in their cm system
1702 cc---------------------------------------------------------------------
1703 c      include 'epos.inc'
1704 c      double precision pi1,pi2,pi3,pi4,pi5,pj1,pj2,pj3,pj4,pj5
1705 c     *,p1,p2,p3,p4,p5,err,pt,pti,sinp,cosp,pmax,phi,drangen!,rrr
1706 c      iret=0
1707 c      pi1=dble(pptl(1,i))
1708 c      pi2=dble(pptl(2,i))
1709 c      pi3=dble(pptl(3,i))
1710 c      pi5=dble(pptl(5,i))
1711 c      pi4=sqrt(pi1**2+pi2**2+pi3**2+pi5**2)
1712 c      pj1=dble(pptl(1,j))
1713 c      pj2=dble(pptl(2,j))
1714 c      pj3=dble(pptl(3,j))
1715 c      pj5=dble(pptl(5,j))
1716 c      pj4=sqrt(pj1**2+pj2**2+pj3**2+pj5**2)
1717 c      if(ish.ge.6)then
1718 c        write(ifch,*)'uthalf for ',i,' and ',j
1719 c        write(ifch,*)'in ',idptl(i),pi1,pi2,pi3,pi4,pi5
1720 c        write(ifch,*)'<> ',idptl(j),pj1,pj2,pj3,pj4,pj5
1721 c      endif
1722 c      p1=pi1+pj1
1723 c      p2=pi2+pj2
1724 c      p3=pi3+pj3
1725 c      p4=pi4+pj4
1726 c      p5=(p4-p3)*(p4+p3)-p2**2-p1**2
1727 c      if(p5.lt.0d0.or.(pi3.lt.0.99d0*pj3.and.mod(ityptl(j)/10,10).eq.4)
1728 c     &       .or.(pi3.gt.0.99d0*pj3.and.mod(ityptl(j)/10,10).eq.5))then
1729 cc      if(p5.lt.0d0)then
1730 c        if(ish.ge.7)write(ifch,*)'Inversion not possible (1)',p5,pi3,pj3
1731 c        iret=1
1732 c        return
1733 c      else
1734 c        p5=sqrt(p5)
1735 c      endif
1736 c      call utlob2(1,p1,p2,p3,p4,p5,pi1,pi2,pi3,pi4,50)
1737 c      call utlob2(1,p1,p2,p3,p4,p5,pj1,pj2,pj3,pj4,51)
1738 c           err=(pi1+pj1)**2+(pi2+pj2)**2+(pi3+pj3)**2
1739 c           if(err.gt.1d-3)then
1740 c      call utmsg('uthalf')
1741 c      write(ifch,*)'*****  err=',err
1742 c      write(ifch,*)'pi:',pi1,pi2,pi3,pi4
1743 c      write(ifch,*)'pj:',pj1,pj2,pj3,pj4
1744 c      call utmsgf
1745 c           endif
1746 c      if(ish.ge.8)then
1747 c        write(ifch,*)'pi 1:',pi1,pi2,pi3,pi4,pi5
1748 c     &       ,sqrt(pi1**2+pi2**2+pi3**2+pi5**2)
1749 c        write(ifch,*)'pj 1:',pj1,pj2,pj3,pj4,pj5
1750 c     &       ,sqrt(pj1**2+pj2**2+pj3**2+pj5**2)
1751 c      endif
1752 c
1753 c      phi=drangen(p5)*2d0*dble(pi)
1754 c      pti=sqrt(pi1*pi1+pi2*pi2)
1755 c
1756 cc      sinp=abs(asin(pi2/pti))
1757 c
1758 c      cosp=cos(phi)
1759 c      sinp=sin(phi)
1760 cc      cosp=pi1/pti
1761 cc      sinp=pi2/pti
1762 c      pini=abs(pi3)
1763 c      pmax=pini
1764 c
1765 c      ntry=0
1766 c 10   ntry=ntry+1
1767 c      pi1=-pj1
1768 c      pi2=-pj2
1769 c      pi3=-pj3
1770 c
1771 c
1772 cc      rrr=dble(ranptcut(4.))*max(0.1d0,dble(zz))+0.01d0
1773 cc      rrr=dble(ranptcut(zz*max(1.,float(ntry/10)**3)))        !flip if pt too large
1774 cc      rrr=dble(min(1.,ranptcut(zz/float(max(1,ntry/10)))))        !return if pt too large
1775 c
1776 c      rrr=rangen()
1777 c      phi=dble(pi*zz1*(1.-rrr*zz2/float(max(1,ntry/10))))
1778 c      call utroa2(phi,cosp,sinp,0d0,pi1,pi2,pi3)    !rotation around an axis perpendicular to p3
1779 c      pt=pi1*pi1+pi2*pi2
1780 c
1781 ccc      pi3=pini*(1d0-(pmax/pini+1d0)*min(1d0,rrr))
1782 ccc     &   *sign(1.d0,pi3)
1783 cc      pi3=pini*(1d0-(pmax/pini+1d0)*(1d0-min(1d0,rrr)))
1784 cc     &   *sign(1.d0,pi3)
1785 cccc      pi3=(1d0-exp(-drangen(pti)))*pi3
1786 cccc      pi3=min(1d0,(drangen()+exp(-min(100d0,dble(zz)))))*pi3
1787 cccc      pi3=(dble(reminv)+exp(-min(100d0,dble(zz))))*pi3
1788 cc      pt=((pi4+pi3)*(pi4-pi3)-pi5*pi5)
1789 c      if(ish.ge.8)write(ifch,*)'ut',ntry,zz,rrr,-pi3/pj3,pt,pti*pti
1790 c     &,idptl(i),idptl(j)
1791 c      if((pt.lt.0d0.or.pt.gt.max(2.*pti*pti,1.d0))
1792 cc      if(pt.lt.0d0
1793 c     &                               .and.ntry.lt.1000)then
1794 c        goto 10
1795 c      elseif(ntry.ge.1000)then
1796 c        if(ish.ge.7)write(ifch,*)'Inversion not possible (2)',pt,pti*pti
1797 c        iret=1
1798 c        return         !pion distribution fall at xf=1
1799 cc        pi3=pj3       !if flip all particle with large pt
1800 cc        pt=pti*pti    !then very very hard pi+ spectra (flat at xf=1 like in pp400 !)
1801 c                       !but to hard for K, rho, eta, etc ..
1802 c      endif
1803 cc      print *,'ut',ntry,phi/(pi-sinp),zz,rrr,-pi3/pj3,pt/pti/pti
1804 cc      pt=sqrt(pt)
1805 cc      pi1=cosp*pt
1806 cc      pi2=sinp*pt
1807 c
1808 c      pj1=-pi1
1809 c      pj2=-pi2
1810 c      pj3=-pi3
1811 cc      pi4=sqrt(pi3*pi3+pt*pt+pi5*pi5)
1812 cc      pj4=sqrt(pj3*pj3+pt*pt+pj5*pj5)
1813 c      pi4=sqrt(pi3*pi3+pt+pi5*pi5)
1814 c      pj4=sqrt(pj3*pj3+pt+pj5*pj5)
1815 c      if(ish.ge.8)then
1816 c        write(ifch,*)'pi 2:',pi1,pi2,pi3,pi4,pi5
1817 c     &       ,sqrt(pi1**2+pi2**2+pi3**2+pi5**2)
1818 c        write(ifch,*)'pj 2:',pj1,pj2,pj3,pj4,pj5
1819 c     &       ,sqrt(pj1**2+pj2**2+pj3**2+pj5**2)
1820 c      endif
1821 c      call utlob2(-1,p1,p2,p3,p4,p5,pi1,pi2,pi3,pi4,-50)
1822 c      call utlob2(-1,p1,p2,p3,p4,p5,pj1,pj2,pj3,pj4,-51)
1823 c      if(pi3/dble(pptl(3,i)).gt.1.00001d0
1824 c     &.or.(pi3.gt.1.001d0*pj3.and.mod(ityptl(j)/10,10).eq.4)
1825 c     &.or.(pi3.lt.1.001d0*pj3.and.mod(ityptl(j)/10,10).eq.5))then
1826 c        if(ish.ge.7)write(ifch,*)'Inversion not possible (3)',pi3,pj3
1827 c        iret=1
1828 c        return
1829 c      endif
1830 c      id=idptl(i)
1831 c      pptl(1,i)=sngl(pi1)
1832 c      pptl(2,i)=sngl(pi2)
1833 c      pptl(3,i)=sngl(pi3)
1834 c      pptl(4,i)=sngl(pi4)
1835 c      pptl(5,i)=sngl(pi5)
1836 cc      pptl(5,i)=sngl(pj5)
1837 cc      idptl(i)=idptl(j)
1838 c      pptl(1,j)=sngl(pj1)
1839 c      pptl(2,j)=sngl(pj2)
1840 c      pptl(3,j)=sngl(pj3)
1841 c      pptl(4,j)=sngl(pj4)
1842 c      pptl(5,j)=sngl(pj5)
1843 cc      pptl(5,j)=sngl(pi5)
1844 cc      idptl(j)=id
1845 c      if(ish.ge.6)then
1846 c        write(ifch,*)'out ',idptl(i),pi1,pi2,pi3,pi4
1847 c        write(ifch,*)' <> ',idptl(j),pj1,pj2,pj3,pj4
1848 c      endif
1849 c
1850 c      return
1851 c      end
1852 c
1853 cc-----------------------------------------------------------------------
1854 c      subroutine utchm(arp,arm,ii)
1855 cc-----------------------------------------------------------------------
1856 cc     checks whether arp**2=0 and arm**2=0.
1857 cc-----------------------------------------------------------------------
1858 c      include 'epos.inc'
1859 c      double precision arp(4),arm(4),difp,difm
1860 c      difp=arp(4)**2-arp(1)**2-arp(2)**2-arp(3)**2
1861 c      difm=arm(4)**2-arm(1)**2-arm(2)**2-arm(3)**2
1862 c      if(dabs(difp).gt.1e-3*arp(4)**2
1863 c     *.or.dabs(difm).gt.1e-3*arm(4)**2)then
1864 c      call utmsg('utchm ')
1865 c      write(ifch,*)'*****  mass non zero  -  ',ii
1866 c      write(ifch,*)'jet-mass**2`s:    ',difp,difm
1867 c      write(ifch,*)'energy**2`s:      ',arp(4)**2,arm(4)**2
1868 c      write(ifch,*)(sngl(arp(i)),i=1,4)
1869 c      write(ifch,*)(sngl(arm(i)),i=1,4)
1870 c      call utmsgf
1871 c      endif
1872 c      return
1873 c      end
1874 c
1875 c-----------------------------------------------------------------------
1876       subroutine utclea(nptlii,nptl0)
1877 c-----------------------------------------------------------------------
1878 c     starting from nptlii
1879 c     overwrites istptl=99 particles in /cptl/, reduces so nptl
1880 c     and update minfra and maxfra
1881 c-----------------------------------------------------------------------
1882       include 'epos.inc'
1883       integer newptl(mxptl)!,oldptl(mxptl),ii(mxptl)
1884 
1885       ish0=ish
1886       if(ishsub/100.eq.18)ish=mod(ishsub,100)
1887 
1888       call utpri('utclea',ish,ishini,2)
1889 
1890       nptli=max(maproj+matarg+1,nptlii)
1891       minfra0=minfra
1892       maxfra0=maxfra
1893       minfra1=maxfra
1894       maxfra1=minfra
1895       if(ish.ge.2)write(ifch,*)'entering subr utclea:',nptl
1896      &                                                ,minfra,maxfra
1897       if(ish.ge.7)then
1898       write(ifch,*)('-',l=1,68)
1899       write(ifch,*)'sr utclea. initial.'
1900       write(ifch,*)('-',l=1,68)
1901       do 34 n=nptli,nptl
1902       write(ifch,116)iorptl(n),jorptl(n),n,ifrptl(1,n),ifrptl(2,n)
1903      *,idptl(n),sqrt(pptl(1,n)**2+pptl(2,n)**2),pptl(3,n),pptl(5,n)
1904      *,istptl(n),ityptl(n)
1905 34    continue
1906 116   format(1x,i6,i6,4x,i6,4x,i6,i6,i12,3x,3(e8.2,1x),i3,i3)
1907       endif
1908 
1909 c      ish=ish0
1910 c      ish0=ish
1911 c      if(ishsub/100.eq.18)ish=mod(ishsub,100)
1912 
1913       i=nptli-1
1914 1     i=i+1
1915       if(i.gt.nptl)goto 1000
1916       if(istptl(i).eq.99)goto 2
1917       newptl(i)=i
1918 c      oldptl(i)=i
1919       goto 1
1920 
1921 2     i=i-1
1922       j=i
1923 3     i=i+1
1924 4     j=j+1
1925       if(j.gt.nptl)goto 5
1926       newptl(j)=0
1927       if(istptl(j).eq.99)goto 4
1928       newptl(j)=i
1929 c      oldptl(i)=j
1930 c      write(ifch,*)'move',j,' to ',i
1931 c       write(ifch,*)idptl(i),ityptl(i),idptl(j),ityptl(j),minfra,maxfra
1932       call utrepl(i,j)
1933       if(j.ge.minfra0.and.j.le.maxfra0)then
1934         minfra1=min(minfra1,i)
1935         maxfra1=max(maxfra1,i)
1936       endif
1937       goto 3
1938 
1939 5     nptl=i-1
1940       if(nptl.eq.0)then
1941         nptl0=0
1942         goto 1000
1943       endif
1944 
1945 20    n0=newptl(nptl0)
1946       if(n0.gt.0)then
1947       nptl0=n0
1948       else
1949       nptl0=nptl0-1
1950       if(nptl0.gt.0)goto 20
1951       endif
1952 
1953 
1954 c      do 11 k=1,nptl
1955 c      io=iorptl(k)
1956 c      if(io.le.0)ii(k)=io
1957 c      if(io.gt.0)ii(k)=newptl(io)
1958 c11    continue
1959 c      do 12 k=1,nptl
1960 c12    iorptl(k)=ii(k)
1961 c
1962 c      do 13 k=1,nptl
1963 c      jo=jorptl(k)
1964 c      if(jo.le.0)ii(k)=jo
1965 c      if(jo.gt.0)ii(k)=newptl(jo)
1966 c13    continue
1967 c      do 14 k=1,nptl
1968 c14    jorptl(k)=ii(k)
1969 c
1970 c      do 15 k=1,nptl
1971 c      if1=ifrptl(1,k)
1972 c      if(if1.le.0)ii(k)=if1
1973 c      if(if1.gt.0)ii(k)=newptl(if1)
1974 c15    continue
1975 c      do 16 k=1,nptl
1976 c16    ifrptl(1,k)=ii(k)
1977 c
1978 c      do 17 k=1,nptl
1979 c      if2=ifrptl(2,k)
1980 c      if(if2.le.0)ii(k)=if2
1981 c      if(if2.gt.0)ii(k)=newptl(if2)
1982 c17    continue
1983 c      do 18 k=1,nptl
1984 c18    ifrptl(2,k)=ii(k)
1985 c
1986 c      do 19 k=1,nptl
1987 c      if(ifrptl(1,k).eq.0.and.ifrptl(2,k).gt.0)ifrptl(1,k)=ifrptl(2,k)
1988 c      if(ifrptl(2,k).eq.0.and.ifrptl(1,k).gt.0)ifrptl(2,k)=ifrptl(1,k)
1989 c19    continue
1990 
1991 1000  continue
1992 
1993       if(minfra1.lt.minfra0)minfra=minfra1
1994       if(maxfra1.ge.minfra1)maxfra=maxfra1
1995 
1996       if(ish.ge.2)then
1997       write(ifch,*)'exiting subr utclea:'
1998       do 35 n=1,nptl
1999       write(ifch,116)iorptl(n),jorptl(n),n,ifrptl(1,n),ifrptl(2,n)
2000      *,idptl(n),sqrt(pptl(1,n)**2+pptl(2,n)**2),pptl(3,n),pptl(5,n)
2001      *,istptl(n),ityptl(n)
2002 35    continue
2003       endif
2004 
2005       if(ish.ge.2)write(ifch,*)'exiting subr utclea:',nptl
2006      &                                                ,minfra,maxfra
2007 
2008       call utprix('utclea',ish,ishini,2)
2009       ish=ish0
2010       return
2011       end
2012 
2013 c---------------------------------------------------------------------
2014       subroutine utfit(x,y,ndata,sig,mwt,a,b,siga,sigb,chi2,q)
2015 c---------------------------------------------------------------------
2016 c linear fit to data
2017 c input:
2018 c    ndata: nr of data points
2019 c    x(),y(),sig(): data
2020 c    mwt: unweighted (0) or weighted (else) data points
2021 c output:
2022 c    a,b: parameters of linear fit a+b*x
2023 c---------------------------------------------------------------------
2024       INTEGER mwt,ndata
2025       REAL a,b,chi2,q,siga,sigb,sig(ndata),x(ndata),y(ndata)
2026 CU    USES utgmq
2027       INTEGER i
2028       REAL sigdat,ss,st2,sx,sxoss,sy,t,wt,utgmq
2029       sx=0.
2030       sy=0.
2031       st2=0.
2032       b=0.
2033       if(mwt.ne.0) then
2034         ss=0.
2035         do 11 i=1,ndata
2036           wt=1./(sig(i)**2)
2037           ss=ss+wt
2038           sx=sx+x(i)*wt
2039           sy=sy+y(i)*wt
2040 11      continue
2041       else
2042         do 12 i=1,ndata
2043           sx=sx+x(i)
2044           sy=sy+y(i)
2045 12      continue
2046         ss=float(ndata)
2047       endif
2048       sxoss=sx/ss
2049       if(mwt.ne.0) then
2050         do 13 i=1,ndata
2051           t=(x(i)-sxoss)/sig(i)
2052           st2=st2+t*t
2053           b=b+t*y(i)/sig(i)
2054 13      continue
2055       else
2056         do 14 i=1,ndata
2057           t=x(i)-sxoss
2058           st2=st2+t*t
2059           b=b+t*y(i)
2060 14      continue
2061       endif
2062       b=b/st2
2063       a=(sy-sx*b)/ss
2064       siga=sqrt((1.+sx*sx/(ss*st2))/ss)
2065       sigb=sqrt(1./st2)
2066       chi2=0.
2067       if(mwt.eq.0) then
2068         do 15 i=1,ndata
2069           chi2=chi2+(y(i)-a-b*x(i))**2
2070 15      continue
2071         q=1.
2072         sigdat=sqrt(chi2/(ndata-2))
2073         siga=siga*sigdat
2074         sigb=sigb*sigdat
2075       else
2076         do 16 i=1,ndata
2077           chi2=chi2+((y(i)-a-b*x(i))/sig(i))**2
2078 16      continue
2079         q=utgmq(0.5*(ndata-2),0.5*chi2)
2080       endif
2081       return
2082       END
2083 
2084 c-----------------------------------------------------------------------
2085       function utgam1(x)
2086 c-----------------------------------------------------------------------
2087 c  gamma fctn tabulated
2088 c  single precision
2089 c-----------------------------------------------------------------------
2090       double precision utgamtab,utgam,al,dl
2091       common/gamtab/utgamtab(10000)
2092 
2093       if(x.gt.0.01.and.x.lt.99.99)then
2094         al=100.d0*dble(x)
2095         k1=int(al)
2096         k2=k1+1
2097         dl =al-dble(k1)
2098         utgam1=real(utgamtab(k2)*dl+utgamtab(k1)*(1.d0-dl))
2099       elseif(x.eq.0.)then
2100         utgam1=0.
2101       else
2102         utgam1=real(utgam(dble(x)))
2103       endif
2104 
2105       end
2106 
2107 c-----------------------------------------------------------------------
2108       double precision function utgam2(x)
2109 c-----------------------------------------------------------------------
2110 c  gamma fctn tabulated
2111 c  double precision
2112 c-----------------------------------------------------------------------
2113       double precision utgamtab,x,al,dl,utgam
2114       common/gamtab/utgamtab(10000)
2115 
2116       if(x.gt.0.01d0.and.x.le.99.99d0)then
2117         al=100.d0*x
2118         k1=int(al)
2119         k2=k1+1
2120         dl =al-dble(k1)
2121         utgam2=utgamtab(k2)*dl+utgamtab(k1)*(1.d0-dl)
2122       elseif(x.eq.0.d0)then
2123         utgam2=0.d0
2124       else
2125         utgam2=utgam(x)
2126       endif
2127 
2128       end
2129 
2130 c-----------------------------------------------------------------------
2131       double precision function utgam(x)
2132 c-----------------------------------------------------------------------
2133 c  gamma fctn
2134 c  double precision
2135 c-----------------------------------------------------------------------
2136       include 'epos.inc'
2137       double precision c(13),x,z,f
2138       data c
2139      1/ 0.00053 96989 58808, 0.00261 93072 82746, 0.02044 96308 23590,
2140      2  0.07309 48364 14370, 0.27964 36915 78538, 0.55338 76923 85769,
2141      3  0.99999 99999 99998,-0.00083 27247 08684, 0.00469 86580 79622,
2142      4  0.02252 38347 47260,-0.17044 79328 74746,-0.05681 03350 86194,
2143      5  1.13060 33572 86556/
2144       utgam=0d0
2145       z=x
2146       if(x .gt. 170.d0) goto6
2147       if(x .gt. 0.0d0) goto1
2148       if(x .eq. int(x)) goto5
2149       z=1.0d0-z
2150     1 f=1.0d0/z
2151       if(z .le. 1.0d0) goto4
2152       f=1.0d0
2153     2 continue
2154       if(z .lt. 2.0d0) goto3
2155       z=z-1.0d0
2156       f=f*z
2157       goto2
2158     3 z=z-1.0d0
2159     4 utgam=
2160      1 f*((((((c(1)*z+c(2))*z+c(3))*z+c(4))*z+c(5))*z+c(6))*z+c(7))/
2161      2 ((((((c(8)*z+c(9))*z+c(10))*z+c(11))*z+c(12))*z+c(13))*z+1.0d0)
2162       if(x .gt. 0.0d0) return
2163       utgam=3.141592653589793d0/(sin(3.141592653589793d0*x)*utgam)
2164       return
2165     5 write(ifch,10)sngl(x)
2166    10 format(1x,'argument of gamma fctn = ',e20.5)
2167       call utstop('utgam : negative integer argument&')
2168     6 write(ifch,11)sngl(x)
2169    11 format(1x,'argument of gamma fctn = ',e20.5)
2170       call utstop('utgam : argument too large&')
2171       end
2172 
2173 c---------------------------------------------------------------------
2174       subroutine utgcf(gammcf,a,x,gln)
2175 c---------------------------------------------------------------------
2176       INTEGER ITMAX
2177       REAL a,gammcf,gln,x,EPS,FPMIN
2178       PARAMETER (ITMAX=100,EPS=3.e-7,FPMIN=1.e-30)
2179 CU    USES utgmln
2180       INTEGER i
2181       REAL an,b,c,d,del,h,utgmln
2182       gln=utgmln(a)
2183       b=x+1.-a
2184       c=1./FPMIN
2185       d=1./b
2186       h=d
2187       do 11 i=1,ITMAX
2188         an=-i*(i-a)
2189         b=b+2.
2190         d=an*d+b
2191         if(abs(d).lt.FPMIN)d=FPMIN
2192         c=b+an/c
2193         if(abs(c).lt.FPMIN)c=FPMIN
2194         d=1./d
2195         del=d*c
2196         h=h*del
2197         if(abs(del-1.).lt.EPS)goto 1
2198 11    continue
2199       call utstop("a too large, ITMAX too small in utgcf&")
2200 1     gammcf=exp(-x+a*log(x)-gln)*h
2201       return
2202       END
2203 
2204 c---------------------------------------------------------------------
2205       function utgmln(xx)
2206 c---------------------------------------------------------------------
2207       REAL utgmln,xx
2208       INTEGER j
2209       DOUBLE PRECISION ser,stp,tmp,x,y,cof(6)
2210       SAVE cof,stp
2211       DATA cof,stp/76.18009172947146d0,-86.50532032941677d0,
2212      *24.01409824083091d0,-1.231739572450155d0,.1208650973866179d-2,
2213      *-.5395239384953d-5,2.5066282746310005d0/
2214       x=xx
2215       y=x
2216       tmp=x+5.5d0
2217       tmp=(x+0.5d0)*log(tmp)-tmp
2218       ser=1.000000000190015d0
2219       do 11 j=1,6
2220         y=y+1.d0
2221         ser=ser+cof(j)/y
2222 11    continue
2223       utgmln=tmp+log(stp*ser/x)
2224       return
2225       END
2226 
2227 c---------------------------------------------------------------------
2228       function utgmq(a,x)
2229 c---------------------------------------------------------------------
2230       REAL a,utgmq,x
2231 CU    USES utgcf,utgser
2232       REAL gammcf,gamser,gln
2233       if(x.lt.0..or.a.le.0.) call utstop("bad arguments in utgmq&")
2234       if(x.lt.a+1.)then
2235         call utgser(gamser,a,x,gln)
2236         utgmq=1.-gamser
2237       else
2238         call utgcf(gammcf,a,x,gln)
2239         utgmq=gammcf
2240       endif
2241       return
2242       END
2243 
2244 c---------------------------------------------------------------------
2245       subroutine utgser(gamser,a,x,gln)
2246 c---------------------------------------------------------------------
2247       INTEGER ITMAX
2248       REAL a,gamser,gln,x,EPS
2249       PARAMETER (ITMAX=100,EPS=3.e-7)
2250 CU    USES utgmln
2251       INTEGER n
2252       REAL ap,del,sum,utgmln
2253       gln=utgmln(a)
2254       if(x.le.0.)then
2255         if(x.lt.0.)call utstop("x < 0 in utgser&")
2256         gamser=0.
2257         return
2258       endif
2259       ap=a
2260       sum=1./a
2261       del=sum
2262       do 11 n=1,ITMAX
2263         ap=ap+1.
2264         del=del*x/ap
2265         sum=sum+del
2266         if(abs(del).lt.abs(sum)*EPS)goto 1
2267 11    continue
2268       call utstop("a too large, ITMAX too small in utgser&")
2269 1     gamser=sum*exp(-x+a*log(x)-gln)
2270       return
2271       END
2272 
2273 c-------------------------------------------------------------------------
2274       subroutine uticpl(ic,ifla,iqaq,iret)
2275 c-------------------------------------------------------------------------
2276 c  adds a quark (iqaq=1) or antiquark (iqaq=2) of flavour ifla
2277 c  to 2-id ic
2278 c-------------------------------------------------------------------------
2279       include 'epos.inc'
2280       integer jc(nflav,2),ic(2)
2281       iret=0
2282       if(ifla.eq.0)return
2283       call iddeco(ic,jc)
2284       if(ish.ge.8)write(ifch,'(2i8,12i3)')ic,jc
2285       jqaq=3-iqaq
2286       if(jc(ifla,jqaq).gt.0)then
2287       jc(ifla,jqaq)=jc(ifla,jqaq)-1
2288       else
2289       jc(ifla,iqaq)=jc(ifla,iqaq)+1
2290       endif
2291       call idcomj(jc)
2292       call idenco(jc,ic,ireten)
2293       if(ish.ge.8)write(ifch,'(2i8,12i3)')ic,jc
2294       if(ireten.eq.1)iret=1
2295       if(ic(1).eq.0.and.ic(2).eq.0.and.ireten.eq.0)then
2296       ic(1)=100000
2297       ic(2)=100000
2298       endif
2299       return
2300       end
2301 
2302 cc-----------------------------------------------------------------------
2303 c      subroutine utindx(n,xar,x,i)
2304 cc-----------------------------------------------------------------------
2305 cc  input:  dimension n
2306 cc          array xar(n) with xar(i) > xar(i-1)
2307 cc          some number x between xar(1) and xar(n)
2308 cc  output: the index i such that x is between xar(i)  and xar(i+1)
2309 cc-----------------------------------------------------------------------
2310 c      include 'epos.inc'
2311 c      real xar(n)
2312 c           if(x.lt.xar(1))then
2313 c      if(ish.ge.5)then
2314 c      call utmsg('utindx')
2315 c      write(ifch,*)'*****  x=',x,' < xar(1)=',xar(1)
2316 c      call utmsgf
2317 c      endif
2318 c      i=1
2319 c      return
2320 c           elseif(x.gt.xar(n))then
2321 c      if(ish.ge.5)then
2322 c      call utmsg('utindx')
2323 c      write(ifch,*)'*****  x=',x,' > xar(n)=',xar(n)
2324 c      call utmsgf
2325 c      endif
2326 c      i=n
2327 c      return
2328 c           endif
2329 c      lu=1
2330 c      lo=n
2331 c1     lz=(lo+lu)/2
2332 c      if((xar(lu).le.x).and.(x.le.xar(lz)))then
2333 c      lo=lz
2334 c      elseif((xar(lz).lt.x).and.(x.le.xar(lo)))then
2335 c      lu=lz
2336 c      else
2337 c      call utstop('utindx: no interval found&')
2338 c      endif
2339 c      if((lo-lu).ge.2) goto1
2340 c      if(lo.le.lu)call utstop('utinvt: lo.le.lu&')
2341 c      i=lu
2342 c      return
2343 c      end
2344 c
2345 c-----------------------------------------------------------------------
2346       function utinvt(n,x,q,y)
2347 c-----------------------------------------------------------------------
2348 c     returns x with y=q(x)
2349 c-----------------------------------------------------------------------
2350       include 'epos.inc'
2351       real x(n),q(n)
2352       if(q(n).eq.0.)call utstop('utinvt: q(n)=0&')
2353            if(y.lt.0.)then
2354       if(ish.ge.1)then
2355       call utmsg('utinvt')
2356       write(ifch,*)'*****  y=',y,' < 0'
2357       call utmsgf
2358       endif
2359       y=0.
2360            elseif(y.gt.q(n))then
2361       if(ish.ge.1)then
2362       call utmsg('utinvt')
2363       write(ifch,*)'*****  y=',y,' > ',q(n)
2364       call utmsgf
2365       endif
2366       y=q(n)
2367            endif
2368       lu=1
2369       lo=n
2370 1     lz=(lo+lu)/2
2371       if((q(lu).le.y).and.(y.le.q(lz)))then
2372       lo=lz
2373       elseif((q(lz).lt.y).and.(y.le.q(lo)))then
2374       lu=lz
2375       else
2376       write(ifch,*)'q(1),y,q(n):',q(1),y,q(n)
2377       write(ifch,*)'lu,lz,lo:',lu,lz,lo
2378       write(ifch,*)'q(lu),q(lz),q(lo):',q(lu),q(lz),q(lo)
2379       call utstop('utinvt: no interval found&')
2380       endif
2381       if((lo-lu).ge.2) goto1
2382       if(lo.le.lu)call utstop('utinvt: lo.le.lu&')
2383       utinvt=x(lu)+(y-q(lu))*(x(lo)-x(lu))/(q(lo)-q(lu))
2384       return
2385       end
2386 
2387 c-----------------------------------------------------------------------
2388       subroutine utlob2(isig,p1,p2,p3,p4,p5,x1,x2,x3,x4,idi)
2389 c-----------------------------------------------------------------------
2390 c  performs a lorentz boost, double prec.
2391 c  isig=+1 is to boost the four vector x1,x2,x3,x4 such as to obtain it
2392 c  in the frame specified by the 5-vector p1...p5 (5-vector=4-vector+mass).
2393 c  isig=-1: the other way round, that means,
2394 c  if the 4-vector x1...x4 is given in some frame characterized by
2395 c  p1...p5 with respect to to some lab-frame, utlob2 returns the 4-vector
2396 c  x1...x4  in the lab frame.
2397 c  idi is a call identifyer (integer) to identify the call in case of problem
2398 c-----------------------------------------------------------------------
2399       include 'epos.inc'
2400       double precision beta(4),z(4),p1,p2,p3,p4,p5,pp,bp,x1,x2,x3,x4
2401      *,xx0,x10,x20,x30,x40,x4x,x0123
2402            if(ish.ge.2)then
2403       if(ish.ge.9)then
2404       write(ifch,101)x1,x2,x3,x4,(x4-x3)*(x4+x3)-x2*x2-x1*x1
2405       write(ifch,301)p1,p2,p3,p4,p5,(p4-p3)*(p4+p3)-p2*p2-p1*p1
2406 101   format(' utlob2: x =  ',5e13.5)
2407 301   format('         p =  ',6e13.5)
2408       endif
2409       pp=(p4-p3)*(p4+p3)-p2*p2-p1*p1
2410       if(dabs(pp-p5*p5).gt.1e-3*p4*p4.and.dabs(pp-p5*p5).gt.1e-3)then
2411       call utmsg('utlob2')
2412       write(ifch,*)'*****  p**2 .ne. p5**2'
2413       write(ifch,*)'call identifyer:',idi
2414       write(ifch,*)'p**2,p5**2: ',pp,p5*p5
2415       write(ifch,*)'p: ',p1,p2,p3,p4,p5
2416       call utmsgf
2417       endif
2418       x10=x1
2419       x20=x2
2420       x30=x3
2421       x40=x4
2422            endif
2423       xx0=(x4-x3)*(x4+x3)-x2*x2-x1*x1
2424       if(p5.le.0.)then
2425       call utmsg('utlob2')
2426       write(ifch,*)'*****  p5 negative.'
2427       write(ifch,*)'call identifyer:',idi
2428       write(ifch,*)'p(5): ',p1,p2,p3,p4,p5
2429       write(ifmt,*)'call identifyer:',idi
2430       write(ifmt,*)'p(5): ',p1,p2,p3,p4,p5
2431       call utmsgf
2432       call utstop('utlob2: p5 negative.&')
2433       endif
2434       z(1)=x1
2435       z(2)=x2
2436       z(3)=x3
2437       z(4)=x4
2438       beta(1)=-p1/p5
2439       beta(2)=-p2/p5
2440       beta(3)=-p3/p5
2441       beta(4)= p4/p5
2442       bp=0.
2443       do 220 k=1,3
2444 220   bp=bp+z(k)*isig*beta(k)
2445       do 230 k=1,3
2446 230   z(k)=z(k)+isig*beta(k)*z(4)
2447      *+isig*beta(k)*bp/(beta(4)+1.)
2448       z(4)=beta(4)*z(4)+bp
2449       x1=z(1)
2450       x2=z(2)
2451       x3=z(3)
2452       x4=z(4)
2453       if(ish.ge.9)
2454      *write(ifch,101)x1,x2,x3,x4,(x4-x3)*(x4+x3)-x2*x2-x1*x1
2455       x4x=x4
2456       x0123=xx0+x1*x1+x2*x2+x3*x3
2457       if(x0123.gt.0.)then
2458       x4=sign( dsqrt(x0123) , x4x )
2459       else
2460       x4=0
2461       endif
2462       if(ish.ge.9)then
2463       write(ifch,101)x1,x2,x3,x4,(x4-x3)*(x4+x3)-x2*x2-x1*x1
2464       endif
2465            if(ish.ge.2)then
2466       if(ish.ge.9)write(ifch,*)'check x**2_ini -- x**2_fin'
2467       if(dabs(x4-x4x).gt.1d-2*dabs(x4).and.dabs(x4-x4x).gt.1d-2)then
2468       call utmsg('utlob2')
2469       write(ifch,*)'*****  x**2_ini .ne. x**2_fin.'
2470       write(ifch,*)'call identifyer:',idi
2471       write(ifch,*)'x1 x2 x3 x4 x**2 (initial/final/corrected):'
2472 102   format(5e13.5)
2473       write(ifch,102)x10,x20,x30,x40,(x40-x30)*(x40+x30)-x20*x20-x10*x10
2474       write(ifch,102)x1,x2,x3,x4x,(x4x-x3)*(x4x+x3)-x2*x2-x1*x1
2475       write(ifch,102)x1,x2,x3,x4,(x4-x3)*(x4+x3)-x2*x2-x1*x1
2476       call utmsgf
2477       endif
2478            endif
2479       if(ish.ge.9)write(ifch,*)'return from utlob2'
2480       return
2481       end
2482 
2483 c-----------------------------------------------------------------------
2484       subroutine utlob3(isig,p1,p2,p3,p4,p5,x1,x2,x3,x4)
2485 c-----------------------------------------------------------------------
2486 c  performs a lorentz boost, double prec.
2487 c  but arguments are single precision
2488 c-----------------------------------------------------------------------
2489       double precision xx1,xx2,xx3,xx4
2490       xx1=dble(x1)
2491       xx2=dble(x2)
2492       xx3=dble(x3)
2493       xx4=dble(x4)
2494       call utlob2(isig
2495      *,dble(p1),dble(p2),dble(p3),dble(p4),dble(p5)
2496      *,xx1,xx2,xx3,xx4,52)
2497       x1=sngl(xx1)
2498       x2=sngl(xx2)
2499       x3=sngl(xx3)
2500       x4=sngl(xx4)
2501       return
2502       end
2503 
2504 c-----------------------------------------------------------------------
2505       subroutine utlob5(yboost,x1,x2,x3,x4,x5)
2506 c-----------------------------------------------------------------------
2507       amt=sqrt(x5**2+x1**2+x2**2)
2508       y=sign(1.,x3)*alog((x4+abs(x3))/amt)
2509       y=y-yboost
2510       x4=amt*cosh(y)
2511       x3=amt*sinh(y)
2512       return
2513       end
2514 
2515 c-----------------------------------------------------------------------
2516       subroutine utlob4(isig,pp1,pp2,pp3,pp4,pp5,x1,x2,x3,x4)
2517 c-----------------------------------------------------------------------
2518 c  performs a lorentz boost, double prec.
2519 c  but arguments are partly single precision
2520 c-----------------------------------------------------------------------
2521       double precision xx1,xx2,xx3,xx4,pp1,pp2,pp3,pp4,pp5
2522       xx1=dble(x1)
2523       xx2=dble(x2)
2524       xx3=dble(x3)
2525       xx4=dble(x4)
2526       call utlob2(isig,pp1,pp2,pp3,pp4,pp5,xx1,xx2,xx3,xx4,53)
2527       x1=sngl(xx1)
2528       x2=sngl(xx2)
2529       x3=sngl(xx3)
2530       x4=sngl(xx4)
2531       return
2532       end
2533 
2534 
2535 c-----------------------------------------------------------------------
2536       subroutine utlobo(isig,p1,p2,p3,p4,p5,x1,x2,x3,x4)
2537 c-----------------------------------------------------------------------
2538 c     performs a lorentz boost
2539 c-----------------------------------------------------------------------
2540       include 'epos.inc'
2541       real beta(4),z(4)
2542       if(p5.le.0.)then
2543       call utmsg('utlobo')
2544       write(ifch,*)'*****  mass <= 0.'
2545       write(ifch,*)'p(5): ',p1,p2,p3,p4,p5
2546       call utmsgf
2547       call utstop('utlobo: mass <= 0.&')
2548       endif
2549       z(1)=x1
2550       z(2)=x2
2551       z(3)=x3
2552       z(4)=x4
2553       beta(1)=-p1/p5
2554       beta(2)=-p2/p5
2555       beta(3)=-p3/p5
2556       beta(4)= p4/p5
2557       bp=0.
2558       do 220 k=1,3
2559 220   bp=bp+z(k)*isig*beta(k)
2560       do 230 k=1,3
2561 230   z(k)=z(k)+isig*beta(k)*z(4)
2562      *+isig*beta(k)*bp/(beta(4)+1.)
2563       z(4)=beta(4)*z(4)+bp
2564       x1=z(1)
2565       x2=z(2)
2566       x3=z(3)
2567       x4=z(4)
2568       return
2569       end
2570 
2571 c-----------------------------------------------------------------------
2572       subroutine utloc(ar,n,a,l)
2573 c-----------------------------------------------------------------------
2574       real ar(n)
2575       do 1 i=1,n
2576       l=i-1
2577       if(a.lt.ar(i))return
2578 1     continue
2579       l=n
2580       return
2581       end
2582 
2583 cc-----------------------------------------------------------------------
2584 c      subroutine utlow(cone)
2585 cc-----------------------------------------------------------------------
2586 c      character*1 cone
2587 c      if(cone.eq.'A')cone='a'
2588 c      if(cone.eq.'B')cone='b'
2589 c      if(cone.eq.'C')cone='c'
2590 c      if(cone.eq.'D')cone='d'
2591 c      if(cone.eq.'E')cone='e'
2592 c      if(cone.eq.'F')cone='f'
2593 c      if(cone.eq.'G')cone='g'
2594 c      if(cone.eq.'H')cone='h'
2595 c      if(cone.eq.'I')cone='i'
2596 c      if(cone.eq.'J')cone='j'
2597 c      if(cone.eq.'K')cone='k'
2598 c      if(cone.eq.'L')cone='l'
2599 c      if(cone.eq.'M')cone='m'
2600 c      if(cone.eq.'N')cone='n'
2601 c      if(cone.eq.'O')cone='o'
2602 c      if(cone.eq.'P')cone='p'
2603 c      if(cone.eq.'Q')cone='q'
2604 c      if(cone.eq.'R')cone='r'
2605 c      if(cone.eq.'S')cone='s'
2606 c      if(cone.eq.'T')cone='t'
2607 c      if(cone.eq.'U')cone='u'
2608 c      if(cone.eq.'V')cone='v'
2609 c      if(cone.eq.'W')cone='w'
2610 c      if(cone.eq.'X')cone='x'
2611 c      if(cone.eq.'Y')cone='y'
2612 c      if(cone.eq.'Z')cone='z'
2613 c      return
2614 c      end
2615 c
2616 cc-----------------------------------------------------------------------
2617 c      subroutine utlow3(cthree)
2618 cc-----------------------------------------------------------------------
2619 c      character cthree*3
2620 c      do 1 i=1,3
2621 c1     call utlow(cthree(i:i))
2622 c      return
2623 c      end
2624 c
2625 cc-----------------------------------------------------------------------
2626 c      subroutine utlow6(csix)
2627 cc-----------------------------------------------------------------------
2628 c      character csix*6
2629 c      do 1 i=1,6
2630 c1     call utlow(csix(i:i))
2631 c      return
2632 c      end
2633 c
2634 cc-----------------------------------------------------------------------
2635 c      function utmom(k,n,x,q)
2636 cc-----------------------------------------------------------------------
2637 cc     calculates kth moment for f(x) with q(i)=int[0,x(i)]f(z)dz
2638 cc-----------------------------------------------------------------------
2639 c      real x(n),q(n)
2640 c      if(n.lt.2)call utstop('utmom : dimension too small&')
2641 c      utmom=0
2642 c      do 1 i=2,n
2643 c1     utmom=utmom+((x(i)+x(i-1))/2)**k*(q(i)-q(i-1))
2644 c      utmom=utmom/q(n)
2645 c      return
2646 c      end
2647 c
2648 c-----------------------------------------------------------------------
2649       function utpcm(a,b,c)
2650 c-----------------------------------------------------------------------
2651 c     calculates cm momentum for a-->b+c
2652 c-----------------------------------------------------------------------
2653       val=(a*a-b*b-c*c)*(a*a-b*b-c*c)-(2.*b*c)*(2.*b*c)
2654       if(val.lt.0..and.val.gt.-1e-4)then
2655       utpcm=0
2656       return
2657       endif
2658       utpcm=sqrt(val)/(2.*a)
2659       return
2660       end
2661 
2662 c-----------------------------------------------------------------------
2663       double precision function utpcmd(a,b,c,iret)
2664 c-----------------------------------------------------------------------
2665 c     calculates cm momentum for a-->b+c
2666 c-----------------------------------------------------------------------
2667       double precision a,b,c,val
2668       iret=0
2669       val=(a*a-b*b-c*c)*(a*a-b*b-c*c)-(2.*b*c)*(2.*b*c)
2670       utpcmd=0d0
2671       if(val.lt.0d0.and.val.gt.-1d-4)then
2672         return
2673       elseif(val.lt.0d0)then
2674         iret=1
2675         return
2676       endif
2677       utpcmd=sqrt(val)/(2.d0*a)
2678       return
2679       end
2680 
2681 c-----------------------------------------------------------------------
2682       subroutine utpri(text,ishi,ishini,ishx)
2683 c-----------------------------------------------------------------------
2684       include 'epos.inc'
2685       character*6 text
2686 c      double precision seedx                               !!!
2687       ishini=ishi
2688       if(ishevt.ne.0.and.nrevt+1.ne.ishevt)return
2689       if(nrpri.gt.0)then
2690       do nr=1,nrpri
2691       if(subpri(nr)(1:6).eq.text)then
2692       ishi=ishpri(nr)
2693       endif
2694       enddo
2695       endif
2696       if(ish.ge.ishx)then
2697         write(ifch,'(1x,43a)')
2698      *  ('-',i=1,10),' entry ',text,' ',('-',i=1,30)
2699 c       call ranfgt(seedx)                                   !!!
2700 c       if(ish.ge.ishx)write(ifch,*)'seed:',seedx            !!!
2701       endif
2702       return
2703       end
2704 
2705 c-----------------------------------------------------------------------
2706       subroutine utprix(text,ishi,ishini,ishx)
2707 c-----------------------------------------------------------------------
2708       include 'epos.inc'
2709       character*6 text
2710       if(ishevt.ne.0.and.nrevt+1.ne.ishevt)return
2711       if(ish.ge.ishx)write(ifch,'(1x,44a)')
2712      *('-',i=1,30),' exit ',text,' ',('-',i=1,11)
2713       ishi=ishini
2714       return
2715       end
2716 
2717 c-----------------------------------------------------------------------
2718       subroutine utprj(text,ishi,ishini,ishx)
2719 c-----------------------------------------------------------------------
2720       include 'epos.inc'
2721       character*20 text
2722 c      double precision seedx                               !!!
2723       idx=index(text,' ')-1
2724       ishini=ishi
2725       if(ishevt.ne.0.and.nrevt+1.ne.ishevt)return
2726       if(nrpri.gt.0)then
2727       do nr=1,nrpri
2728       if(subpri(nr)(1:idx).eq.text(1:idx))then
2729       ishi=ishpri(nr)
2730       endif
2731       enddo
2732       endif
2733       if(ish.ge.ishx)then
2734         write(ifch,'(1x,43a)')
2735      *  ('-',i=1,10),' entry ',text(1:idx),' ',('-',i=1,30)
2736 c       call ranfgt(seedx)                                   !!!
2737 c       if(ish.ge.ishx)write(ifch,*)'seed:',seedx            !!!
2738       endif
2739       return
2740       end
2741 
2742 c-----------------------------------------------------------------------
2743       subroutine utprjx(text,ishi,ishini,ishx)
2744 c-----------------------------------------------------------------------
2745       include 'epos.inc'
2746       character*20 text
2747       idx=index(text,' ')-1
2748       if(ishevt.ne.0.and.nrevt+1.ne.ishevt)return
2749       if(ish.ge.ishx)write(ifch,'(1x,44a)')
2750      *('-',i=1,30),' exit ',text(1:idx),' ',('-',i=1,11)
2751       ishi=ishini
2752       return
2753       end
2754 
2755 c-----------------------------------------------------------------------
2756       function utquad(m,x,f,k)
2757 c-----------------------------------------------------------------------
2758 c     performs an integration according to simpson
2759 c-----------------------------------------------------------------------
2760       real x(m),f(m)
2761       utquad=0
2762       do 1 i=1,k-1
2763   1   utquad=utquad+(f(i)+f(i+1))/2*(x(i+1)-x(i))
2764       return
2765       end
2766 
2767 c-----------------------------------------------------------------------
2768       subroutine utquaf(fu,n,x,q,x0,x1,x2,x3)
2769 c-----------------------------------------------------------------------
2770 c     returns q(i) = integral [x(1)->x(i)] fu(x) dx
2771 c-----------------------------------------------------------------------
2772       include 'epos.inc'
2773       real x(n),q(n)
2774       parameter (m=10)
2775       real xa(m),fa(m)
2776       external fu
2777       if(x1.lt.x0.or.x2.lt.x1.or.x3.lt.x2)then
2778       if(ish.ge.1)then
2779       call utmsg('utquaf')
2780       write(ifch,*)'   xi=',x0,x1,x2,x3
2781       call utmsgf
2782       endif
2783       endif
2784       call utar(n/3,n*2/3,n,x0,x1,x2,x3,x)
2785       q(1)=0
2786       do 2 i=2,n
2787       do 3 k=1,m
2788       z=x(i-1)+(k-1.)/(m-1.)*(x(i)-x(i-1))
2789       xa(k)=z
2790 3     fa(k)=fu(z)
2791       q(i)=q(i-1)+utquad(m,xa,fa,m)
2792 2     continue
2793       return
2794       end
2795 
2796 c-----------------------------------------------------------------------
2797       subroutine utrepl(i,j)
2798 c-----------------------------------------------------------------------
2799 c     i is replaced by j in /cptl/
2800 c-----------------------------------------------------------------------
2801       include 'epos.inc'
2802       do 1 k=1,5
2803 1     pptl(k,i)  =pptl(k,j)
2804       iorptl(i)  = 0 !iorptl(j)
2805       idptl(i)   =idptl(j)
2806       istptl(i)  =istptl(j)
2807       do 2 k=1,2
2808 2     tivptl(k,i)=tivptl(k,j)
2809       do 3 k=1,2
2810 3     ifrptl(k,i)= 0 !ifrptl(k,j)
2811       jorptl(i)  = 0 !jorptl(j)
2812       do 4 k=1,4
2813 4     xorptl(k,i)=xorptl(k,j)
2814       do 5 k=1,4
2815 5     ibptl(k,i) =ibptl(k,j)
2816       ityptl(i)  =ityptl(j)
2817       iaaptl(i)  =iaaptl(j)
2818       radptl(i)  =radptl(j)
2819       desptl(i)  =desptl(j)
2820       dezptl(i)  =dezptl(j)
2821       qsqptl(i)  =qsqptl(j)
2822       zpaptl(1,i)=zpaptl(1,j)
2823       zpaptl(2,i)=zpaptl(2,j)
2824       itsptl(i)  =itsptl(j)
2825       rinptl(i)  =rinptl(j)
2826       return
2827       end
2828 
2829 c-----------------------------------------------------------------------
2830       subroutine utrepla(i,j)
2831 c-----------------------------------------------------------------------
2832 c     i is replaced by j in /cptl/
2833 c-----------------------------------------------------------------------
2834       include 'epos.inc'
2835       do 1 k=1,5
2836 1     pptl(k,i)  =pptl(k,j)
2837       iorptl(i)  = iorptl(j)
2838       idptl(i)   =idptl(j)
2839       istptl(i)  =istptl(j)
2840       do 2 k=1,2
2841 2     tivptl(k,i)=tivptl(k,j)
2842       do 3 k=1,2
2843 3     ifrptl(k,i)= ifrptl(k,j)
2844       jorptl(i)  = jorptl(j)
2845       do 4 k=1,4
2846 4     xorptl(k,i)=xorptl(k,j)
2847       do 5 k=1,4
2848 5     ibptl(k,i) =ibptl(k,j)
2849       ityptl(i)  =ityptl(j)
2850       iaaptl(i)  =iaaptl(j)
2851       radptl(i)  =radptl(j)
2852       desptl(i)  =desptl(j)
2853       dezptl(i)  =dezptl(j)
2854       qsqptl(i)  =qsqptl(j)
2855       zpaptl(1,i)=zpaptl(1,j)
2856       zpaptl(2,i)=zpaptl(2,j)
2857       itsptl(i)  =itsptl(j)
2858       rinptl(i)  =rinptl(j)
2859       return
2860       end
2861 
2862 cc-----------------------------------------------------------------------
2863 c      subroutine utresm(icp1,icp2,icm1,icm2,amp,idpr,iadj,ireten)
2864 cc-----------------------------------------------------------------------
2865 c      parameter (nflav=6)
2866 c      integer icm(2),icp(2),jcm(nflav,2),jcp(nflav,2)
2867 c      icm(1)=icm1
2868 c      icm(2)=icm2
2869 c      icp(1)=icp1
2870 c      icp(2)=icp2
2871 c      CALL IDDECO(ICM,JCM)
2872 c      CALL IDDECO(ICP,JCP)
2873 c      do 37 nf=1,nflav
2874 c      do 37 k=1,2
2875 c37    jcP(nf,k)=jcp(nf,k)+jcm(nf,k)
2876 c      CALL IDENCO(JCP,ICP,IRETEN)
2877 c      IDP=IDTRA(ICP,0,0,3)
2878 c      call idres(idp,amp,idpr,iadj)
2879 c      return
2880 c      end
2881 c
2882 cc-----------------------------------------------------------------------
2883 c      subroutine utroa1(phi,a1,a2,a3,x1,x2,x3)
2884 cc-----------------------------------------------------------------------
2885 cc  rotates x by angle phi around axis a.
2886 cc  normalization of a is irrelevant.
2887 cc-----------------------------------------------------------------------
2888 c      double precision aaa,aa(3),xxx,xx(3),e1(3),e2(3),e3(3),xp,xt,dphi
2889 c      dphi=phi
2890 c      xx(1)=x1
2891 c      xx(2)=x2
2892 c      xx(3)=x3
2893 c      aa(1)=a1
2894 c      aa(2)=a2
2895 c      aa(3)=a3
2896 c      aaa=0
2897 c      xxx=0
2898 c      do i=1,3
2899 c      aaa=aaa+aa(i)**2
2900 c      xxx=xxx+xx(i)**2
2901 c      enddo
2902 c      if(xxx.eq.0d0)return
2903 c      if(aaa.eq.0d0)call utstop('utroa1: zero rotation axis&')
2904 c      aaa=dsqrt(aaa)
2905 c      xxx=dsqrt(xxx)
2906 cc e3 = a / !a!
2907 c      do i=1,3
2908 c      e3(i)=aa(i)/aaa
2909 c      enddo
2910 cc x_parallel
2911 c      xp=0
2912 c      do i=1,3
2913 c      xp=xp+xx(i)*e3(i)
2914 c      enddo
2915 cc x_transverse
2916 c      if(xxx**2-xp**2.le.0.)return
2917 c      xt=dsqrt(xxx**2-xp**2)
2918 cc e1 = vector x_transverse / absolute value x_transverse
2919 c      do i=1,3
2920 c      e1(i)=(xx(i)-e3(i)*xp)/xt
2921 c      enddo
2922 cc e2 orthogonal e3,e1
2923 c      call utvec2(e3,e1,e2)
2924 cc rotate x
2925 c      do i=1,3
2926 c      xx(i)=xp*e3(i)+xt*dcos(dphi)*e1(i)+xt*dsin(dphi)*e2(i)
2927 c      enddo
2928 cc back to single precision
2929 c      x1=xx(1)
2930 c      x2=xx(2)
2931 c      x3=xx(3)
2932 c      return
2933 c      end
2934 c
2935 c-----------------------------------------------------------------------
2936       subroutine utroa1(phi,a1,a2,a3,x1,x2,x3)
2937 c-----------------------------------------------------------------------
2938 c  rotates x by angle phi around axis a (argument single precision)
2939 c  normalization of a is irrelevant.
2940 c-----------------------------------------------------------------------
2941       double precision aa(3),xx(3),dphi
2942       dphi=phi
2943       xx(1)=x1
2944       xx(2)=x2
2945       xx(3)=x3
2946       aa(1)=a1
2947       aa(2)=a2
2948       aa(3)=a3
2949       call utroa2(dphi,aa(1),aa(2),aa(3),xx(1),xx(2),xx(3))
2950 c back to single precision
2951       x1=sngl(xx(1))
2952       x2=sngl(xx(2))
2953       x3=sngl(xx(3))
2954       return
2955       end
2956 
2957 c-----------------------------------------------------------------------
2958       subroutine utroa2(phi,a1,a2,a3,x1,x2,x3)
2959 c-----------------------------------------------------------------------
2960 c  rotates x by angle phi around axis a.
2961 c  normalization of a is irrelevant.
2962 c  double precision phi,a1,a2,a3,x1,x2,x3
2963 c-----------------------------------------------------------------------
2964       double precision phi,a1,a2,a3,x1,x2,x3
2965       double precision aaa,aa(3),xxx,xx(3),e1(3),e2(3),e3(3),xp,xt,dphi
2966       dphi=phi
2967       xx(1)=x1
2968       xx(2)=x2
2969       xx(3)=x3
2970       aa(1)=a1
2971       aa(2)=a2
2972       aa(3)=a3
2973       aaa=0d0
2974       xxx=0d0
2975       do i=1,3
2976       aaa=aaa+aa(i)**2
2977       xxx=xxx+xx(i)**2
2978       enddo
2979       if(xxx.eq.0d0)return
2980       if(aaa.eq.0d0)call utstop('utroa1: zero rotation axis&')
2981       aaa=1.0/dsqrt(aaa)
2982 c e3 = a / !a!
2983       do i=1,3
2984       e3(i)=aa(i)*aaa
2985       enddo
2986 c x_parallel
2987       xp=0
2988       do i=1,3
2989       xp=xp+xx(i)*e3(i)
2990       enddo
2991 c x_transverse
2992       if(xxx-xp**2.le.0.)return
2993       xt=dsqrt(xxx-xp**2)
2994 c e1 = vector x_transverse / absolute value x_transverse
2995       do i=1,3
2996       e1(i)=(xx(i)-e3(i)*xp)/xt
2997       enddo
2998 c e2 orthogonal e3,e1
2999       call utvec2(e3,e1,e2)
3000 c rotate x
3001       do i=1,3
3002       xx(i)=xp*e3(i)+xt*cos(dphi)*e1(i)+xt*sin(dphi)*e2(i)
3003       enddo
3004       xxx=0d0
3005       do i=1,3
3006       xxx=xxx+xx(i)**2
3007       enddo
3008 c back to single precision
3009       x1=xx(1)
3010       x2=xx(2)
3011       x3=xx(3)
3012       return
3013       end
3014 
3015 cc--------------------------------------------------------------------
3016 c      function utroot(funcd,x1,x2,xacc)
3017 cc--------------------------------------------------------------------
3018 cc combination of newton-raphson and bisection method for root finding
3019 cc input:
3020 cc   funcd: subr returning fctn value and first derivative
3021 cc   x1,x2: x-interval
3022 cc   xacc:  accuracy
3023 cc output:
3024 cc   utroot: root
3025 cc--------------------------------------------------------------------
3026 c      include 'epos.inc'
3027 c      INTEGER MAXIT
3028 c      REAL utroot,x1,x2,xacc
3029 c      EXTERNAL funcd
3030 c      PARAMETER (MAXIT=100)
3031 c      INTEGER j
3032 c      REAL df,dx,dxold,f,fh,fl,temp,xh,xl
3033 c      call funcd(x1,fl,df)
3034 c      call funcd(x2,fh,df)
3035 c      if((fl.gt.0..and.fh.gt.0.).or.(fl.lt.0..and.fh.lt.0.))
3036 c     *call utstop('utroot: root must be bracketed&')
3037 c      if(fl.eq.0.)then
3038 c        utroot=x1
3039 c        return
3040 c      else if(fh.eq.0.)then
3041 c        utroot=x2
3042 c        return
3043 c      else if(fl.lt.0.)then
3044 c        xl=x1
3045 c        xh=x2
3046 c      else
3047 c        xh=x1
3048 c        xl=x2
3049 c      endif
3050 c      utroot=.5*(x1+x2)
3051 c      dxold=abs(x2-x1)
3052 c      dx=dxold
3053 c      call funcd(utroot,f,df)
3054 c      do 11 j=1,MAXIT
3055 c        if(((utroot-xh)*df-f)*((utroot-xl)*df-f).ge.0..or. abs(2.*
3056 c     *f).gt.abs(dxold*df) ) then
3057 c          dxold=dx
3058 c          dx=0.5*(xh-xl)
3059 c          utroot=xl+dx
3060 c          if(xl.eq.utroot)return
3061 c        else
3062 c          dxold=dx
3063 c          dx=f/df
3064 c          temp=utroot
3065 c          utroot=utroot-dx
3066 c          if(temp.eq.utroot)return
3067 c        endif
3068 c        if(abs(dx).lt.xacc) return
3069 c        call funcd(utroot,f,df)
3070 c        if(f.lt.0.) then
3071 c          xl=utroot
3072 c        else
3073 c          xh=utroot
3074 c        endif
3075 c11    continue
3076 c      call utmsg('utroot')
3077 c      write(ifch,*)'*****  exceeding maximum iterations'
3078 c      write(ifch,*)'dx:',dx
3079 c      call utmsgf
3080 c      return
3081 c      END
3082 c
3083 c-----------------------------------------------------------------------
3084       subroutine utrot2(isig,ax,ay,az,x,y,z)
3085 c-----------------------------------------------------------------------
3086 c     performs a rotation, double prec.
3087 c-----------------------------------------------------------------------
3088       include 'epos.inc'
3089       double precision ax,ay,az,x,y,z,rx,ry,rz
3090      *,alp,bet,cosa,sina,cosb,sinb,xs,ys,zs
3091          if(ax**2.eq.0.and.ay**2.eq.0.and.az**2.eq.0.)then
3092       write(ifch,*)'ax**2,ay**2,az**2:',ax**2,ay**2,az**2
3093       write(ifch,*)'ax,ay,az:',ax,ay,az
3094       call utstop('utrot2: zero vector.&')
3095          endif
3096          if(az.ge.0.)then
3097       rx=ax
3098       ry=ay
3099       rz=az
3100          else
3101       rx=-ax
3102       ry=-ay
3103       rz=-az
3104          endif
3105       if(rz**2+ry**2.ne.0.)then
3106       alp=dabs(dacos(rz/dsqrt(rz**2+ry**2)))*sign(1.,sngl(ry))
3107       bet=
3108      *dabs(dacos(dsqrt(rz**2+ry**2)/dsqrt(rz**2+ry**2+rx**2)))*
3109      *sign(1.,sngl(rx))
3110       else
3111       alp=3.1415927d0/2d0
3112       bet=3.1415927d0/2d0
3113       endif
3114       cosa=dcos(alp)
3115       sina=dsin(alp)
3116       cosb=dcos(bet)
3117       sinb=dsin(bet)
3118            if(isig.ge.0)then
3119       xs=x*cosb-y*sina*sinb-z*cosa*sinb
3120       ys=       y*cosa     -z*sina
3121       zs=x*sinb+y*sina*cosb+z*cosa*cosb
3122            else     !if(isig.lt.0)then
3123       xs= x*cosb            +z*sinb
3124       ys=-x*sinb*sina+y*cosa+z*cosb*sina
3125       zs=-x*sinb*cosa-y*sina+z*cosb*cosa
3126            endif
3127       x=xs
3128       y=ys
3129       z=zs
3130       return
3131       end
3132 
3133 
3134 c-----------------------------------------------------------------------
3135       subroutine utrot4(isig,ax,ay,az,x,y,z)
3136 c-----------------------------------------------------------------------
3137 c     performs a rotation, double prec.
3138 c     arguments partly single
3139 c-----------------------------------------------------------------------
3140       double precision ax,ay,az,xx,yy,zz
3141       xx=dble(x)
3142       yy=dble(y)
3143       zz=dble(z)
3144       call utrot2(isig,ax,ay,az,xx,yy,zz)
3145       x=sngl(xx)
3146       y=sngl(yy)
3147       z=sngl(zz)
3148       return
3149       end
3150 
3151 c-----------------------------------------------------------------------
3152       subroutine utrota(isig,ax,ay,az,x,y,z)
3153 c-----------------------------------------------------------------------
3154 c     performs a rotation
3155 c-----------------------------------------------------------------------
3156          if(az.ge.0.)then
3157       rx=ax
3158       ry=ay
3159       rz=az
3160          else
3161       rx=-ax
3162       ry=-ay
3163       rz=-az
3164          endif
3165       if(rz.eq.0..and.ry.eq.0.)then
3166         alp=0.
3167         stop
3168       else
3169         alp=abs(utacos(rz/sqrt(rz**2+ry**2)))*sign(1.,ry)
3170       endif
3171       bet=
3172      *abs(utacos(sqrt(rz**2+ry**2)/sqrt(rz**2+ry**2+rx**2)))*sign(1.,rx)
3173       cosa=cos(alp)
3174       sina=sin(alp)
3175       cosb=cos(bet)
3176       sinb=sin(bet)
3177            if(isig.ge.0)then
3178       xs=x*cosb-y*sina*sinb-z*cosa*sinb
3179       ys=       y*cosa     -z*sina
3180       zs=x*sinb+y*sina*cosb+z*cosa*cosb
3181            else        !if(isig.lt.0)then
3182       xs= x*cosb            +z*sinb
3183       ys=-x*sinb*sina+y*cosa+z*cosb*sina
3184       zs=-x*sinb*cosa-y*sina+z*cosb*cosa
3185            endif
3186       x=xs
3187       y=ys
3188       z=zs
3189       return
3190       end
3191 
3192 c-----------------------------------------------------------------------
3193       subroutine utstop(text)
3194 c-----------------------------------------------------------------------
3195 c  returns error message and stops execution.
3196 c  text is an optonal text to appear in the error message.
3197 c  text is a character string of length 40;
3198 c     for shorter text, it has to be terminated by &;
3199 c        example: call utstop('error in subr xyz&')
3200 c-----------------------------------------------------------------------
3201       include 'epos.inc'
3202 c      parameter(itext=40)
3203       character  text*(*)  ,txt*6
3204       imax=index(text,'&')
3205       do 1 j=1,2
3206       if(j.eq.1)then
3207         ifi=ifch
3208       else        !if(j.eq.2)
3209         ifi=ifmt
3210       endif
3211       if(imax.gt.1)then
3212       write(ifi,101)('*',k=1,72),text(1:imax-1)
3213      *,nrevt+1,nint(seedj),seedc,('*',k=1,72)
3214       else
3215       write(ifi,101)('*',k=1,72),' '
3216      *,nrevt+1,nint(seedj),seedc,('*',k=1,72)
3217       endif
3218 101   format(
3219      *1x,72a1
3220      */1x,'***** stop in ',a
3221      */1x,'***** current event number: ',i12
3222      */1x,'***** initial seed for current run:',i10
3223      */1x,'***** initial seed for current event:',d25.15
3224      */1x,72a1)
3225 1     continue
3226 c      c=0.
3227 c      b=a/c
3228       stop
3229       entry utmsg(txt)
3230       imsg=imsg+1
3231       write(ifch,'(1x,74a1)')('*',j=1,72)
3232       write(ifch,100)txt,nrevt+1,nint(seedj),seedc
3233 100   format(1x,'***** msg from ',a6,'.   es:',i7,2x,i9,2x,d23.17)
3234       return
3235       entry utmsgf
3236       if(ish.eq.1)return
3237       write(ifch,'(1x,74a1)')('*',j=1,72)
3238       end
3239 
3240 c-----------------------------------------------------------------
3241       subroutine uttrap(func,a,b,s)
3242 c-----------------------------------------------------------------
3243 c trapezoidal method for integration.
3244 c input: fctn func and limits a,b
3245 c output: value s of the integral
3246 c-----------------------------------------------------------------
3247       include 'epos.inc'
3248 
3249       INTEGER JMAX
3250       REAL a,b,func,s
3251       EXTERNAL func
3252       PARAMETER (JMAX=10)
3253 CU    USES uttras
3254       INTEGER j
3255       REAL olds
3256       olds=-1.e30
3257       do 11 j=1,JMAX
3258         if(ish.ge.9)write(ifch,*)'sr uttrap:   j:',j
3259         call uttras(func,a,b,s,j)
3260         ds=abs(s-olds)
3261         if (ds.lt.epsr*abs(olds)) return
3262         olds=s
3263 11    continue
3264 
3265 c-c   nepsr=nepsr+1
3266       if(ish.ge.9)then
3267       call utmsg('uttrap')
3268       write(ifch,*)
3269      *'*****  requested accuracy could not be achieved'
3270       write(ifch,*)'achieved accuracy: ',ds/abs(olds)
3271       write(ifch,*)'requested accuracy:',epsr
3272       call utmsgf
3273       endif
3274 
3275       END
3276 
3277 c-----------------------------------------------------------------
3278       subroutine uttraq(func,a,b,s)
3279 c-----------------------------------------------------------------
3280 c trapezoidal method for integration.
3281 c input: function func and limits a,b
3282 c output: value s of the integral
3283 c-----------------------------------------------------------------
3284 
3285       REAL a,b,func,s
3286       EXTERNAL func
3287       PARAMETER (eps=1.e-6)
3288 CU    USES uttras
3289       INTEGER j
3290       REAL olds
3291       olds=-1.e30
3292       j=1
3293 10      call uttras(func,a,b,s,j)
3294         ds=abs(s-olds)
3295         if (ds.le.eps*abs(olds)) return
3296         olds=s
3297         if(j.ge.15)return
3298         j=j+1
3299       goto10
3300       END
3301 
3302 c-----------------------------------------------------------------
3303       subroutine uttras(func,a,b,s,n)
3304 c-----------------------------------------------------------------
3305 c performs one iteration of the trapezoidal method for integration
3306 c-----------------------------------------------------------------
3307       INTEGER n
3308       REAL a,b,s,func
3309       EXTERNAL func
3310       INTEGER it,j
3311       REAL del,sum,tnm,x
3312       if (n.eq.1) then
3313         s=0.5*(b-a)*(func(a)+func(b))
3314       else
3315         it=2**(n-2)
3316         tnm=it
3317         del=(b-a)/tnm
3318         x=a+0.5*del
3319         sum=0.
3320         do 11 j=1,it
3321           sum=sum+func(x)
3322           x=x+del
3323 11      continue
3324         s=0.5*(s+(b-a)*sum/tnm)
3325       endif
3326       return
3327       END
3328 
3329 cc-----------------------------------------------------------------------
3330 c      subroutine utvec1(a,b,c)
3331 cc-----------------------------------------------------------------------
3332 cc  returns vector product c = a x b .
3333 cc-----------------------------------------------------------------------
3334 c      real a(3),b(3),c(3)
3335 c      c(1)=a(2)*b(3)-a(3)*b(2)
3336 c      c(2)=a(3)*b(1)-a(1)*b(3)
3337 c      c(3)=a(1)*b(2)-a(2)*b(1)
3338 c      return
3339 c      end
3340 c
3341 c-----------------------------------------------------------------------
3342       subroutine utvec2(a,b,c)
3343 c-----------------------------------------------------------------------
3344 c  returns vector product c = a x b .
3345 c  a,b,c double precision.
3346 c-----------------------------------------------------------------------
3347       double precision a(3),b(3),c(3)
3348       c(1)=a(2)*b(3)-a(3)*b(2)
3349       c(2)=a(3)*b(1)-a(1)*b(3)
3350       c(3)=a(1)*b(2)-a(2)*b(1)
3351       return
3352       end
3353 
3354 c-------------------------------------------------------------------
3355       subroutine utword(line,i,j,iqu)
3356 c-------------------------------------------------------------------
3357 c  finds the first word of the character string line(j+1:1000).
3358 c  the word is line(i:j) (with new i and j).
3359 c  if j<0 or if no word found --> new line read.
3360 c  a text between quotes "..." is considered one word;
3361 c  stronger: a text between double quotes ""..."" is consid one word
3362 c  stronger: a text between "{ and }" is considered one word
3363 c-------------------------------------------------------------------
3364 c  input:
3365 c    line: character string (*1000)
3366 c    i: integer between 1 and 1000
3367 c    iqu: for iqu=1 a "?" is written to output before reading a line,
3368 c         otherwise (iqu/=1) nothing is typed
3369 c  output:
3370 c    i,j: left and right end of word (word=line(i:j))
3371 c-------------------------------------------------------------------
3372       include 'epos.inc'
3373       parameter(mempty=2)
3374       character*1 empty(mempty),mk
3375       character line*1000
3376       character*2 mrk
3377       data empty/' ',','/
3378       parameter(mxdefine=40)
3379       character w1define*100,w2define*100
3380       common/cdefine/ndefine,l1define(mxdefine),l2define(mxdefine)
3381      &               ,w1define(mxdefine),w2define(mxdefine)
3382 
3383       j0=0
3384       if(j.ge.0)then
3385       i=j
3386       goto 1
3387       endif
3388 
3389     5 continue
3390       if(iqu.eq.1.and.iprmpt.gt.0)write(ifmt,'(a)')'?'
3391       if(nopen.eq.0)then
3392         ifopx=ifop
3393       elseif(nopen.gt.0)then
3394         ifopx=20+nopen
3395       else !if(nopen.lt.0)
3396         ifopx=ifcp
3397       endif
3398       read(ifopx,'(a1000)',end=9999)line
3399       if(iecho.eq.1.or.(nopen.ge.0.and.kcpopen.eq.1))then
3400        kmax=2
3401        do k=3,1000
3402        if(line(k:k).ne.' ')kmax=k
3403        enddo
3404       else
3405         kmax=0
3406       endif
3407       if(nopen.ge.0.and.kcpopen.eq.1)
3408      &  write(ifcp,'(a)')line(1:kmax)
3409       if(iecho.eq.1)
3410      &  write(ifmt,'(a)')line(1:kmax)
3411       i=0
3412 
3413     1 i=i+1
3414       if(i.gt.1000)goto 5
3415       if(line(i:i).eq.'!')goto 5
3416       do ne=1,mempty
3417       if(line(i:i).eq.empty(ne))goto 1
3418       enddo
3419 
3420       nbla=1
3421       mrk='  '
3422       mk=' '
3423       if(line(i:i).eq.'~')mk='~'
3424       if(line(i:i+1).eq.'"{')mrk='}"'
3425       if(line(i:i+1).eq.'""')mrk='""'
3426       if(mrk.ne.'  ')goto 10
3427       if(line(i:i).eq.'"')mk='"'
3428       if(mk.ne.' ')goto 8
3429       j=i-1
3430     6 j=j+1
3431       if(j.gt.1000)goto 7
3432       if(line(j:j).eq.'!')goto 7
3433       do ne=1,mempty
3434       if(line(j:j).eq.empty(ne))goto 7
3435       enddo
3436       goto 6
3437 
3438     8 continue
3439       if(i.ge.1000-1)stop'utword: make line shorter!!!         '
3440       i=i+1
3441       j=i
3442       if(line(j:j).eq.mk)stop'utword: empty string!!!           '
3443     9 j=j+1
3444       if(j.gt.1000)then                 !reach the end of the line
3445         j=j-nbla+2
3446         goto 7
3447       endif
3448       if(line(j:j).eq.' ')then
3449         nbla=nbla+1
3450       else
3451         nbla=2
3452       endif
3453       if(line(j:j).eq.mk)then
3454       line(i-1:i-1)=' '
3455       line(j:j)=' '
3456       goto 7
3457       endif
3458       goto 9
3459 
3460    10 continue
3461       if(i.ge.1000-3)stop'utword: make line shorter!!!!          '
3462       i=i+2
3463       j=i
3464       if(line(j:j+1).eq.mrk)stop'utword: empty string!!!!        '
3465    11 j=j+1
3466       if(j.gt.1000-1)then                 !reach the end of the line
3467         j=j-nbla+2
3468         goto 7
3469       endif
3470       if(line(j:j+1).eq.mrk)then
3471       line(i-2:i-1)='  '
3472       line(j:j+1)='  '
3473       goto 7
3474       endif
3475       if(line(j:j).eq.' ')then
3476         nbla=nbla+1
3477       else
3478         nbla=2
3479       endif
3480       goto 11
3481 
3482     7 j=j-1
3483       !--------#define---------------
3484       if(ndefine.gt.0)then
3485         do ndf=1,ndefine
3486           l1=l1define(ndf)
3487           l2=l2define(ndf)
3488           do i0=i,j+1-l1
3489             if(line(i0:i0-1+l1).eq.w1define(ndf)(1:l1))then
3490               if(l2.eq.l1)then
3491                 line(i0:i0-1+l1)=w2define(ndf)(1:l2)
3492               elseif(l2.lt.l1)then
3493                 line(i0:i0+l2-1)=w2define(ndf)(1:l2)
3494                 do k=i0+l2,i0-1+l1
3495                   line(k:k)=' '
3496                 enddo
3497               elseif(l2.gt.l1)then
3498                 do k=i0+l1,i0+l2-1
3499                   if(line(k:k).ne.' ')
3500      &        stop'utword: no space for `define` replacement.   '
3501                 enddo
3502                 line(i0:i0+l2-1)=w2define(ndf)(1:l2)
3503                 j=i0+l2-1
3504               endif
3505             endif
3506           enddo
3507         enddo
3508        do k=i,j
3509          if(line(k:k).ne.' ')j0=j
3510        enddo
3511        j=j0
3512       endif
3513       !--------
3514       return
3515 
3516 9999  close(ifopx)
3517       nopen=nopen-1
3518       if(nopen.eq.0.and.iprmpt.eq.-1)iprmpt=1
3519       goto 5
3520       end
3521 
3522 c--------------------------------------------------------------------
3523       subroutine utworn(line,j,ne)
3524 c--------------------------------------------------------------------
3525 c  returns number ne of nonempty characters of line(j+1:1000)
3526 c--------------------------------------------------------------------
3527       character line*1000
3528       ne=0
3529       do l=j+1,1000
3530       if(line(l:l).ne.' ')ne=ne+1
3531       enddo
3532       return
3533       end
3534 
3535 c-----------------------------------------------------------------------
3536       subroutine getairmol(iz,ia)
3537 c-----------------------------------------------------------------------
3538       include 'epos.inc'
3539       i=0
3540       r=rangen()
3541       do while(r.gt.0.)  ! choose air-molecule
3542         i=i+1
3543         r=r-airwnxs(i)
3544       enddo
3545       iz = nint(airznxs(i))
3546       ia = nint(airanxs(i))
3547       end
3548 
3549 c----------------------------------------------------------------------
3550 
3551       subroutine factoriel
3552 
3553 c----------------------------------------------------------------------
3554 c tabulation of fctrl(n)=n!, facto(n)=1/n! and utgamtab(x) for x=0 to 50
3555 c----------------------------------------------------------------------
3556       include 'epos.incems'
3557       double precision utgamtab,utgam,x
3558       common/gamtab/utgamtab(10000)
3559 
3560       nfctrl=100
3561       fctrl(0)=1.D0
3562       facto(0)=1.D0
3563       do i=1,min(npommx,nfctrl)
3564         fctrl(i)=fctrl(i-1)*dble(i)
3565         facto(i)=1.d0/fctrl(i)
3566       enddo
3567 
3568       do k=1,10000
3569         x=dble(k)/100.d0
3570         utgamtab(k)=utgam(x)
3571       enddo
3572 
3573       return
3574       end
3575 
3576 c-----------------------------------------------------------------------
3577       subroutine fremnu(amin,ca,cb,ca0,cb0,ic1,ic2,ic3,ic4)
3578 c-----------------------------------------------------------------------
3579       common/hadr2/iomodl,idproj,idtarg,wexcit
3580       real pnll,ptq
3581       common/hadr1/pnll,ptq,exmass,cutmss,wproj,wtarg
3582       parameter (nflav=6)
3583       integer ic(2),jc(nflav,2)
3584       ic(1)=ca
3585       ic(2)=cb
3586       call iddeco(ic,jc)
3587       keu=jc(1,1)-jc(1,2)
3588       ked=jc(2,1)-jc(2,2)
3589       kes=jc(3,1)-jc(3,2)
3590       kec=jc(4,1)-jc(4,2)
3591       keb=jc(5,1)-jc(5,2)
3592       ket=jc(6,1)-jc(6,2)
3593       amin=utamnu(keu,ked,kes,kec,keb,ket,5)  !???4=2mults, 5=1mult
3594         if(ca-ca0.eq.0..and.cb-cb0.eq.0..and.rangen().gt.wexcit)then
3595       ic3=0
3596       ic4=0
3597       ic1=ca
3598       ic2=cb
3599        else
3600       amin=amin+exmass
3601       n=0
3602       do i=1,4
3603       do j=1,2
3604       n=n+jc(i,j)
3605       enddo
3606       enddo
3607       k=1+rangen()*n
3608       do i=1,4
3609       do j=1,2
3610       k=k-jc(i,j)
3611       if(k.le.0)goto 1
3612       enddo
3613       enddo
3614 1     if(j.eq.1)then
3615       ic3=10**(6-i)
3616       ic4=0
3617       else
3618       ic3=0
3619       ic4=10**(6-i)
3620       endif
3621       ic1=int(ca)-ic3
3622       ic2=int(cb)-ic4
3623         endif
3624       return
3625       end
3626 
3627 
3628 c-----------------------------------------------------------------------
3629       function fremnux(jc)
3630 c-----------------------------------------------------------------------
3631       real pnll,ptq
3632       common/hadr1/pnll,ptq,exmass,cutmss,wproj,wtarg
3633       parameter (nflav=6)
3634       integer jc(nflav,2)!,ic(2)
3635 c      ic(1)=210000
3636 c      ic(2)=0
3637 c      call iddeco(ic,jc)
3638       keu=jc(1,1)-jc(1,2)
3639       ked=jc(2,1)-jc(2,2)
3640       kes=jc(3,1)-jc(3,2)
3641       kec=jc(4,1)-jc(4,2)
3642       keb=jc(5,1)-jc(5,2)
3643       ket=jc(6,1)-jc(6,2)
3644       fremnux=utamnu(keu,ked,kes,kec,keb,ket,4) !+exmass  !???4=2mults, 5=1mult
3645       return
3646       end
3647 
3648 c-----------------------------------------------------------------------
3649       function fremnux2(jc)
3650 c-----------------------------------------------------------------------
3651       real pnll,ptq
3652       common/hadr1/pnll,ptq,exmass,cutmss,wproj,wtarg
3653       parameter (nflav=6)
3654       integer jc(nflav,2)!,ic(2)
3655 c      ic(1)=210000
3656 c      ic(2)=0
3657 c      call iddeco(ic,jc)
3658       keu=jc(1,1)-jc(1,2)
3659       ked=jc(2,1)-jc(2,2)
3660       kes=jc(3,1)-jc(3,2)
3661       kec=jc(4,1)-jc(4,2)
3662       keb=jc(5,1)-jc(5,2)
3663       ket=jc(6,1)-jc(6,2)
3664       fremnux2=utamnu(keu,ked,kes,kec,keb,ket,5) !+exmass  !???4=2mults, 5=1mult
3665       return
3666       end
3667 
3668 c-----------------------------------------------------------------------
3669       function fremnux3(jci)
3670 c-----------------------------------------------------------------------
3671 c minimum mass from ic counting all quarks
3672 c-----------------------------------------------------------------------
3673       include 'epos.inc'
3674       integer jc(nflav,2),jci(nflav,2)!,ic(2)
3675 c      ic(1)=210000
3676 c      ic(2)=0
3677 c      print *,'start',ic
3678       fremnux3=0.
3679       do j=1,2
3680         do i=1,nflav
3681           jc(i,j)=jci(i,j)
3682         enddo
3683       enddo
3684 c      call iddeco(ic,jc)
3685       call idquacjc(jc,nqua,naqu)
3686         do ii=1,2
3687       if(ii.eq.1)then
3688         nqu=nqua
3689       else
3690         nqu=naqu
3691       endif
3692       if(nqu.ge.3)then
3693         do while(jc(3,ii).ne.0.and.nqu.ge.3)  !count baryons with s quark
3694           jc(3,ii)=jc(3,ii)-1
3695           if(jc(3,ii).gt.0)then
3696             jc(3,ii)=jc(3,ii)-1
3697             if(jc(3,ii).gt.0)then
3698               jc(3,ii)=jc(3,ii)-1
3699               fremnux3=fremnux3+asuhax(4)
3700             elseif(jc(2,ii).gt.0)then
3701               jc(2,ii)=jc(2,ii)-1
3702               fremnux3=fremnux3+asuhax(4)
3703             elseif(jc(1,ii).gt.0)then
3704               jc(1,ii)=jc(1,ii)-1
3705               fremnux3=fremnux3+asuhax(4)
3706             endif
3707           elseif(jc(2,ii).gt.0)then
3708             jc(2,ii)=jc(2,ii)-1
3709             if(jc(1,ii).gt.0)then
3710               jc(1,ii)=jc(1,ii)-1
3711               fremnux3=fremnux3+asuhax(3)
3712             elseif(jc(2,ii).gt.0)then
3713               jc(2,ii)=jc(2,ii)-1
3714               fremnux3=fremnux3+asuhax(3)
3715             endif
3716           elseif(jc(1,ii).gt.0)then
3717             jc(1,ii)=jc(1,ii)-2
3718             fremnux3=fremnux3+asuhay(3)
3719           endif
3720           nqu=nqu-3
3721         enddo
3722         do while(jc(2,ii).ne.0.and.nqu.ge.3)  !count baryons with d quark
3723           jc(2,ii)=jc(2,ii)-1
3724           if(jc(1,ii).gt.0)then
3725             jc(1,ii)=jc(1,ii)-1
3726             if(jc(2,ii).gt.0)then
3727               jc(2,ii)=jc(2,ii)-1
3728               fremnux3=fremnux3+asuhay(2)
3729             elseif(jc(1,ii).gt.0)then
3730               jc(1,ii)=jc(1,ii)-1
3731               fremnux3=fremnux3+asuhay(2)
3732             endif
3733           elseif(jc(2,ii).gt.0)then
3734             jc(2,ii)=jc(2,ii)-2
3735             fremnux3=fremnux3+asuhay(3)
3736           endif
3737           nqu=nqu-3
3738         enddo
3739         do while(jc(1,ii).ne.0.and.nqu.ge.3)  !count baryons with s quark
3740           jc(1,ii)=jc(1,ii)-3
3741           fremnux3=fremnux3+asuhay(3)
3742           nqu=nqu-3
3743         enddo
3744         if(ii.eq.1)then
3745           nqua=nqu
3746         else
3747           naqu=nqu
3748         endif
3749       endif
3750 c      print *,ii,nqua,naqu,jc,fremnux3
3751       enddo
3752       if(nqua+naqu.ne.0)then
3753       do while(jc(3,1).ne.0)    !count mesons with s quark
3754         jc(3,1)=jc(3,1)-1
3755         if(jc(3,2).gt.0)then
3756           jc(3,2)=jc(3,2)-1
3757           fremnux3=fremnux3+asuhax(6)
3758         elseif(jc(2,2).gt.0)then
3759           jc(2,2)=jc(2,2)-1
3760           fremnux3=fremnux3+asuhay(6)
3761         elseif(jc(1,2).gt.0)then
3762           jc(1,2)=jc(1,2)-1
3763           fremnux3=fremnux3+asuhay(6)
3764         endif
3765       enddo
3766       do while(jc(2,1).ne.0)    !count mesons with d quark
3767         jc(2,1)=jc(2,1)-1
3768         if(jc(2,2).gt.0)then
3769           jc(2,2)=jc(2,2)-1
3770           fremnux3=fremnux3+asuhay(5)
3771         elseif(jc(1,2).gt.0)then
3772           jc(1,2)=jc(1,2)-1
3773           fremnux3=fremnux3+asuhay(5)
3774         endif
3775       enddo
3776       do while(jc(1,1).ne.0)    !count mesons with s quark
3777         jc(1,1)=jc(1,1)-1
3778         if(jc(1,2).gt.0)then
3779           jc(1,2)=jc(1,2)-1
3780           fremnux3=fremnux3+asuhay(5)
3781         endif
3782       enddo
3783       endif
3784 c      fremnux3=fremnux3+0.5
3785 c      print *,'stop',nqua,naqu,fremnux3
3786 
3787       return
3788       end
3789 
3790 c-----------------------------------------------------------------------
3791       subroutine fremnx(ammax,amin,sm,ic3,ic4,iret)
3792 c-----------------------------------------------------------------------
3793       common/psar9/ alpr
3794       include 'epos.inc'
3795       iret=0
3796       if(ic3.eq.0.and.ic4.eq.0)then
3797         if(ammax.lt.amin**2)then
3798           iret=1
3799           return
3800         endif
3801         sm=amin**2
3802       else
3803 c       ammax1=min(ammax,(engy/4.)**2)
3804         ammax1=ammax
3805         if(ammax1.lt.amin**2)then
3806           iret=1
3807           return
3808         endif
3809         if(alpr.eq.-1.)then
3810           sm=amin**2*(ammax1/amin**2)**rangen()
3811         else
3812           sm=amin**2*(1.+((ammax1/amin**2)**(1.+alpr)-1.)
3813      *    *rangen())**(1./(1.+alpr))
3814         endif
3815       endif
3816       return
3817       end
3818 
3819       SUBROUTINE gaulag(x,w,n,alf)
3820       INTEGER n,MAXIT
3821       REAL alf,w(n),x(n)
3822       DOUBLE PRECISION EPS
3823       PARAMETER (EPS=3.D-14,MAXIT=10)
3824 CU    USES gammln
3825       INTEGER i,its,j
3826       REAL ai,gammln
3827       DOUBLE PRECISION p1,p2,p3,pp,z,z1
3828       z=0.
3829       do 13 i=1,n
3830         if(i.eq.1)then
3831           z=(1.+alf)*(3.+.92*alf)/(1.+2.4*n+1.8*alf)
3832         else if(i.eq.2)then
3833           z=z+(15.+6.25*alf)/(1.+.9*alf+2.5*n)
3834         else
3835           ai=i-2
3836           z=z+((1.+2.55*ai)/(1.9*ai)+1.26*ai*alf/(1.+3.5*ai))*
3837      *(z-x(i-2))/(1.+.3*alf)
3838         endif
3839         do 12 its=1,MAXIT
3840           p1=1.d0
3841           p2=0.d0
3842           do 11 j=1,n
3843             p3=p2
3844             p2=p1
3845             p1=((2*j-1+alf-z)*p2-(j-1+alf)*p3)/j
3846 11        continue
3847           pp=(n*p1-(n+alf)*p2)/z
3848           z1=z
3849           z=z1-p1/pp
3850           if(abs(z-z1).le.EPS)goto 1
3851 12      continue
3852         call utstop("too many iterations in gaulag")
3853 1       x(i)=z
3854         w(i)=-exp(gammln(alf+n)-gammln(float(n)))/(pp*n*p2)
3855 13    continue
3856       return
3857       END
3858 C  (C) Copr. 1986-92 Numerical Recipes Software 4+1$!].
3859 
3860       FUNCTION gammln(xx)
3861       REAL gammln,xx
3862       INTEGER j
3863       DOUBLE PRECISION ser,stp,tmp,x,y,cof(6)
3864       SAVE cof,stp
3865       DATA cof,stp/76.18009172947146d0,-86.50532032941677d0,
3866      *24.01409824083091d0,-1.231739572450155d0,.1208650973866179d-2,
3867      *-.5395239384953d-5,2.5066282746310005d0/
3868       x=xx
3869       y=x
3870       tmp=x+5.5d0
3871       tmp=(x+0.5d0)*log(tmp)-tmp
3872       ser=1.000000000190015d0
3873       do 11 j=1,6
3874         y=y+1.d0
3875         ser=ser+cof(j)/y
3876 11    continue
3877       gammln=tmp+log(stp*ser/x)
3878       return
3879       END
3880 C  (C) Copr. 1986-92 Numerical Recipes Software 4+1$!].
3881 
3882       function polar(x,y)
3883 
3884       pi=3.1415927
3885       if(abs(x).gt.1.e-6)then
3886         phi=atan(y/x)
3887         if(x.lt.0.)phi=pi+phi
3888         if(phi.lt.0)phi=2*pi+phi
3889       else
3890         phi=0.5*pi
3891         if(y.lt.0)phi=phi+pi
3892       endif
3893       polar=phi
3894 
3895       end
3896 
3897       subroutine getJKNcentr
3898       end