c$Id: make22.f,v 1.36 2003/05/02 11:06:54 weber Exp $
C####C##1#########2#########3#########4#########5#########6#########7##
      subroutine make22(iio,e,ii1,iiz1,mm1,xfac1,ii2,iiz2,mm2,xfac2)
c
cinput iio   : label for exit-channel
cinput e     : $\sqrt{s}$ of process
cinput ii1   : ID of incoming particle 1
cinput iiz1  : $2\cdot I_3$ of incoming particle 1
cinput mm1   : mass of incoming particle 1
cinput xfac1 : scaling factor for preformed hadron
cinput ii2   : ID of incoming particle 2
cinput iiz2  : $2\cdot I_3$ of incoming particle 2
cinput mm2   : mass of incoming particle 2
cinput xfac2 : scaling factor for preformed hadron
c
c  output:   exit channel via common-blocks in {\tt newpart.f}
c
c {\tt make22} generates the final state for all scatterings and
c decays. Due to the diverse nature of the interactions handled
c many special cases have to be taken care of. The label {\tt iio}
c matches in most cases the respective label in subroutine {\tt crossx},
c which returns the respective partial cross sections.
c
C####C##1#########2#########3#########4#########5#########6#########7##
      implicit none

      include 'coms.f'
      include 'options.f'
      include 'newpart.f'
      include 'comres.f'



      integer io,iio,i1,i2,i3,i4,i,k,i1p,i1m,i2p,i2m
      integer ifqrk1,ifqrk2,ifdiq1,ifdiq2,ifqrk3,ifqrk4,ifdiq3,ifdiq4
      integer iq1(2),iq2(3),kq1,kq2,kqq1,kqq2,iii,iff(3)
      real*8 sig,e,m1,m2,m3,m4,gam,mm1,mm2,m1m,m2m,xfac1,xfac2

      integer iz1,iz2,iz3,iz4,errflg,icnt,ntry,ib1,ib2,i1old,i2old

      logical bit


c...functions
      integer isoit,whichres
      real*8 getmass,fmsr,ranf,pcms,massit,widit,mminit

c...vacuum quantumnumber(s) for special string decay (don't touch)
      real*8 valint(1)
      common /values/ valint

c...string
      real*8 b1,b2,ba1(3),ba2(3)
      integer j,l,ii1,ii2,iiz1,iiz2,iexopt,iddum,ibar,jbar
      logical fboost,switips



      integer ident(2,mprt)
      real*8 part(9,mprt),ms1,ms2,msmin1,msmin2,tau,esum


      bit=.true.
      switips=.false.
      io=mod(iio,200)
      i1=ii1
      i2=ii2
      iz1=iiz1
      iz2=iiz2
      m1=mm1
      m2=mm2
      icnt=0
      ntry=0
      ibar=0

      if(iabs(i1).lt.minmes)ibar=ibar+isign(1,i1)
      if(iabs(i2).lt.minmes)ibar=ibar+isign(1,i2)

c in case of a MB-reaction, sort particles (but keep track of
c any id-switch with the 'switips'-flag)
      if(iabs(i1).ge.minmes.and.iabs(i1).le.maxmes.and.
     &   iabs(i2).ge.minbar.and.iabs(i2).le.maxbar)then
        call swpizm(i1,iz1,m1,i2,iz2,m2)
        switips=.true.
      endif

      if(i1+i2.eq.0.and.iz1+iz2.eq.0.and.CTOption(20).ne.0.and.
     .    io.gt.20)goto 27 !e+e-

      if(io.lt.0)goto(100,100,100,100,100,100,100,29)-io

c      if(i1+i2.gt.2)write(6,*)'make22:',i1,i2
ctp060202 1007 continue
      goto(1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,9,17,9,17,20,
     ,9,13,23,14,12,26,27,15,14,100,29,100,14,15,14,36,36,13)io

      write(6,*)'make22: unknown channel requested io:',io
      write(6,*)'  ',e,i1,iz1,m1,i2,iz2,m2
      stop

 1    continue
c...pp->ND
      i3=minnuc
      i4=mindel
      m3=massit(i3)

      call getmas(massit(i4),widit(i4),i4,isoit(i4),mminit(i4),
     .     e-m3,m3,m4)

      if(ranf(0).gt.0.5d0)call swpizm(i3,iz3,m3,i4,iz4,m4)
      goto 1008

 2    continue
c...pp->pp*
      i3=minnuc
      m3=massit(i3)
      if(bit)call getres(io,e,minnuc+1,maxnuc,i4)
      call getmas(massit(i4),widit(i4),i4,isoit(i4),mminit(i4),
     .     e-m3,m3,m4)

      if(ranf(0).gt.0.5d0)call swpizm(i3,iz3,m3,i4,iz4,m4)
      goto 1008

 3    continue
c...pp->ND*
      i3=minnuc
      m3=massit(i3)
      if(bit)call getres(io,e,mindel+1,maxdel,i4)
      call getmas(massit(i4),widit(i4),i4,isoit(i4),mminit(i4),
     .     e-m3,m3,m4)

      if(ranf(0).gt.0.5d0)call swpizm(i3,iz3,m3,i4,iz4,m4)
      goto 1008

 4    continue
c...pp->DD
      i3=mindel
      i4=mindel
      call getmas(massit(i3),widit(i3),i3,isoit(i3),
     .     mminit(i4),e-mminit(i4),mminit(i4),m3)
      call getmas(massit(i4),widit(i4),i4,isoit(i4),mminit(i4),
     .     e-m3,m3,m4)

      if(ranf(0).gt.0.5d0)call swpizm(i3,iz3,m3,i4,iz4,m4)
      goto 1008

 5    continue
c...pp->DN*, DN* with factor 4/3
      i3=mindel
      if(bit)call getres(io,e,minnuc+1,maxnuc,i4)
      call getmas(massit(i3),widit(i3),i3,isoit(i3),
     .     mminit(i3),e-mminit(i4),mminit(i4),m3)
      call getmas(massit(i4),widit(i4),i4,isoit(i4),mminit(i4),
     .     e-m3,m3,m4)

      if(ranf(0).gt.0.5d0)call swpizm(i3,iz3,m3,i4,iz4,m4)
      goto 1008

 6    continue
c...pp->DD*, DN* with factor 4/3
      i3=mindel
      if(bit)call getres(io,e,mindel+1,maxdel,i4)
      call getmas(massit(i3),widit(i3),i3,isoit(i3),
     .     mminit(i3),e-mminit(i4),mminit(i4),m3)
      call getmas(massit(i4),widit(i4),i4,isoit(i4),mminit(i4),
     .     e-m3,m3,m4)

      if(ranf(0).gt.0.5d0)call swpizm(i3,iz3,m3,i4,iz4,m4)
      goto 1008

 7    continue
c...pp -> generate N*N*,N*D*,D*D*
      m3=fmsr(mminit(mindel),e-mresmin)
      m4=fmsr(mminit(mindel),e-m3)
      if(bit)i3=whichres(m3,3)
      if(bit)i4=whichres(m4,3)
      if(ranf(0).gt.0.5d0)call swpizm(i3,iz3,m3,i4,iz4,m4)
      goto 1008

 8    continue
c...ND->DD
      i3=mindel
      i4=mindel
      call getmas(massit(i3),widit(i3),i3,isoit(i3),
     .     mminit(i4),e-mminit(i4),mminit(i4),m3)
      call getmas(massit(i4),widit(i4),i4,isoit(i4),mminit(i4),
     .     e-m3,m3,m4)

      if(ranf(0).gt.0.5d0)call swpizm(i3,iz3,m3,i4,iz4,m4)
      goto 1008

 9    continue
      write(6,*)'make22: channel no.',io,'not implemented.'
      stop

 10   continue
c...MB->B',MM->M* annihilations
      call anndec(0,m1,i1,iz1,m2,i2,iz2,e,sig,gam)

      goto 2002
c      return

 11   continue
c...MM->M'
      call anndec(0,m1,i1,iz1,m2,i2,iz2,e,sig,gam)

      goto 2002
c      return


 12   continue
      write(6,*)'make22: channel no.',io,'should correspond to',
     ,'total cross section, i.e. make22 should not be called.'
      stop

 13   continue
c...elastic scattering
c
      if(switips)then
         call swpizm(i1,iz1,m1,i2,iz2,m2)
         switips=.false.
      endif

      call setizm(i1,iz1,m1,i2,iz2,m2,i3,iz3,m3,i4,iz4,m4)

      if(mminit(i4)+mminit(i3).gt.e)then
ctp060926         write(6,*)'make22(el):threshold violated'
ctp060926         write(6,*)'m3:',m3,mminit(i3)
ctp060926         write(6,*)'m4:',m4,mminit(i4)
      else
c
c
        if(m3.lt.mminit(i3))m3=mminit(i3)
        if(m4.lt.mminit(i4))m4=mminit(i4)
      end if

      goto 2001

 14   continue
c...inelastic scattering (aqm for nonstrange resonances)
c   for arbitrary resonances getinw should be modified
c   arbitrary particles are excited (whithout charge exchange)

c...get minimal masses
      call getirg(i1,i1m,i1p)
      m1m=mminit(i1m)
      call getirg(i2,i2m,i2p)
      m2m=mminit(i2m)

      if(e-m1m-m2m.le.0d0)then
c...elastic scattering
        call setizm(i1,iz1,m1,i2,iz2,m2,i3,iz3,m3,i4,iz4,m4)
        goto 2001
c        return
      end if

c...masses
      if(ranf(0).lt.5d-1)then
         m3=fmsr(m1m,e-m2m)
         m4=fmsr(m2m,e-m3)
      else
         m4=fmsr(m2m,e-m1m)
         m3=fmsr(m1m,e-m4)
      end if

      if(e-m3-m4.lt.0d0)goto 13

c...itype3
      if(i1m.lt.i1p.and.widit(i1m).lt.1d-3.and.
     .    m3.le.massit(i1m+1)-0.5*widit(i1m+1))then
c...lowest itype of this kind of particles is stable
c & mass .lt. approximate minimal mass of lowest itype + 1
        m3=massit(i1m)
        i3=i1m
      else if(i1m.eq.i1p.and.widit(i1m).lt.1d-3) then
c...class with only one narrow particle
        i3=i1
        m3=massit(i3)
      else
c...itypes of this kind are all unstable
        call whichi(i3,i1m,i1p,m3)
      end if

c...itype4
      if(i2m.lt.i2p.and.widit(i2m).lt.1d-3.and.
     .    m4.le.massit(i2m+1)-0.5*widit(i2m+1))then
c...lowest itype of this kind of particles is stable
        m4=massit(i2m)
        i4=i2m
      else if(i2m.eq.i2p.and.widit(i2m).lt.1d-3) then
c...class with only one narrow particle
        i4=i2
        m4=massit(i4)
      else
c...itypes of this kind are all unstable
        call whichi(i4,i2m,i2p,m4)
      end if

c...no charge transfer
      iz3=iz1
      iz4=iz2
      i3=isign(i3,ii1)
      i4=isign(i4,ii2)
      goto 2001
c      return

 15    continue
c...XX -> 2 strings
      if(CTOption(12).ne.0)then
        write(6,*)' *** error(make22): string section is called ',
     .            'while strings are switched off:CTOption(12).ne.0'
        stop
      end if

c
      if(switips)then
         call swpizm(i1,iz1,m1,i2,iz2,m2)
         switips=.false.
      endif


c store old itypes
      i1old=i1
      i2old=i2

 155  continue

c allow for deexcitation: 'excitation' starts from groundstate

      if(iabs(i1).ge.minmes)then
c.. meson resonances may also be deexcitated: assign the particle
c id's of the lowest multiplet (with the same quark content)
        call ityp2id(i1,iz1,iq1(1),iq1(2))
        if(iabs(iq1(1)).gt.iabs(iq1(2)))then
          iddum=iq1(1)
          iq1(1)=iq1(2)
          iq1(2)=iddum
        endif
        iddum=isign(100*iabs(iq1(1))+10*iabs(iq1(2)),iq1(1))
        call id2ityp(iddum,0d0,i1,iz1)
      else
         call getirg(i1,i1m,i1p)
      endif

c same for second particle
      if(iabs(i2).ge.minmes)then
        call ityp2id(i2,iz2,iq1(1),iq1(2))
        if(iabs(iq1(1)).gt.iabs(iq1(2)))then
          iddum=iq1(1)
          iq1(1)=iq1(2)
          iq1(2)=iddum
        endif
        iddum=isign(100*iabs(iq1(1))+10*iabs(iq1(2)),iq1(1))
        call id2ityp(iddum,0d0,i2,iz2)
      else
         call getirg(i2,i2m,i2p)
      endif

c
c BEWARE: now i1 and i2 do not contain anymore the ityps of the ingoing
c particles, but the ityps of the lowest possible states (groundstates).
c This is needed in order for the string-excitation to be able to
c excite ALL states and not only those above the state of the ingoing
c particle.
c If you need the old ityps (i.e. for elastic scattering) then reset them
c via i1old and i2old


c if the incoming masses are changed due to excitation, then the
c new masses should not be lower than:
      m1m=mminit(i1)
      m2m=mminit(i2)

c..discriminate between BB and MB collisions:
      ib1=1
      ib2=1
      if(iabs(i1).ge.minmes)ib1=0
      if(iabs(i2).ge.minmes)ib2=0


c set minimum energy for the two strings
      msmin1=m1m+CTParam(2)
      msmin2=m2m+CTParam(2)

c no (meson) string below 1 gev:
      if(msmin1.lt.1d0)msmin1=1d0
      if(msmin2.lt.1d0)msmin2=1d0

c if energy too low: do elastic collision
      if(m1m+m2m+CTParam(34).ge.e)then
ctp060926        write(*,*)'make22: not enough energy for string exc. ->elastic/
ctp060926     &deexcitation'
ctp060926        write(*,*)' i1, i2, m1, m2, e: ',i1,i2,m1,m2,e
         i1=i1old
         i2=i2old
        goto 13
      endif


c convert to quark-IDs
         call ityp2id(i1,iz1,ifdiq1,ifqrk1)
         call ityp2id(i2,iz2,ifdiq2,ifqrk2)

         iexopt=CTOption(22)
 81      continue
c 100 tries for excitation, otherwise elastic scattering
         ntry=ntry+1
         if(ntry.gt.100)then
ctp060926          write(*,*)'make22: too many tries for string exc. ->elastic/
ctp060926     &           deexcitation'
ctp060926          write(*,*)' i1, i2, m1, m2, e: ',i1,i2,m1,m2,e
         i1=i1old
         i2=i2old
          goto 13
         endif

c string-excitation:
c get string masses ms1,ms2 and the leading quarks
         call STREXCT(IFdiq1,IFqrk1,ib1,M1m,
     &           ifdiq2,ifqrk2,ib2,M2m,E,
     &           iexopt,
     &           ba1,ms1,ba2,ms2,
     &           ifdiq3,ifqrk3,ifdiq4,ifqrk4)
c the boost parameters are now fixed for the masses ms1, ms2. If the
c masses will be changed, set the parameter fboost to "false":
      fboost=.true.

c accept deexcitation of one of the hadrons:
      if(ms1.le.msmin1.and.ms2.ge.msmin2)then
        ms1=massit(i1)
        fboost=.false.
      else if(ms2.le.msmin2.and.ms1.ge.msmin1)then
        ms2=massit(i2)
        fboost=.false.
c don't accept elastic-like (both masses too low):
      else if(ms1.lt.msmin1.and.ms2.lt.msmin2)then
        goto 81
      endif

c in case of deexcitation new masses are necessary
c single diffractive, mass excitation according 1/m
      if(ms1.le.msmin1)then
        ms2=fmsr(msmin2,e-ms1)
        fboost=.false.
      elseif(ms2.le.msmin2)then
        ms1=fmsr(msmin1,e-ms2)
        fboost=.false.
      endif

c quark-quark scattering -> only elastic !
      if(xfac1.lt..999d0.and.xfac2.lt..999d0
     &     .and.ranf(0).lt..25) then
         i1=i1old
         i2=i2old
         goto 13
      endif

c single diffractive if one particle is a quark state,
c mass excitation according 1/m
      if(xfac1.lt..999d0.and.ranf(0).lt..5)then
        ms1=massit(i1)
        ms2=fmsr(msmin2,e-ms1)
        fboost=.false.
      elseif(xfac2.lt..999d0.and.ranf(0).lt..5)then
        ms2=massit(i2)
        ms1=fmsr(msmin1,e-ms2)
        fboost=.false.
      endif

c take care that the particles will be able to decay lateron:
      if(ms1.lt.m1m)then
        ms1=m1m
        fboost=.false.
      endif
      if(ms2.lt.m2m)then
        ms2=m2m
        fboost=.false.
      endif

c avoid energy conservation violation
      if(ms1+ms2.gt.e)goto 155

       if(CTOption(22).ne.1.or..not.fboost)then
c the boost parameters have to be calculated:
         b1=2.*e*pcms(e,ms1,ms2)/(e**2+ms1**2-ms2**2)
         b2=2.*e*pcms(e,ms1,ms2)/(e**2-ms1**2+ms2**2)
         ba1(3)=+b1
         ba2(3)=-b2
         do 151 j=1,2
            ba1(j)=0d0
 151        ba2(j)=0d0
      end if

c now write information to newpart common-blocks
      mstring(1)=ms1
      mstring(2)=ms2

c string#1
      if(ms1.le.msmin1)then
        nstring1=1
        l=1
        do j=1,3
          part(j,l)=0d0
          part(j+5,l)=0d0
        enddo
        part(4,l)=ms1
        part(4+5,l)=0d0
        part(5,l)=ms1
        ident(1,l)=i1
        ident(2,l)=iz1
      else
       call qstring(ifdiq3,ifqrk3,ms1,part,ident,nstring1)
       if(nstring1.eq.0) goto 155
      end if

      esum=0d0

      jbar=0

      do l=1,nstring1
         pnew(5,l)=part(5,l)
         itypnew(l)=ident(1,l)

         if(iabs(ident(1,l)).lt.minmes)jbar=jbar+isign(1,ident(1,l))

         i3new(l)=ident(2,l)
         do j=1,4
           pnew(j,l)=part(j,l)
           xnew(j,l)=part(j+5,l)
         enddo
         pnew(3,l)=part(3,l)
         xnew(3,l)=part(3+5,l)
         call rotbos(0d0,0d0,ba1(1),ba1(2),ba1(3),
     ,    pnew(1,l),pnew(2,l),pnew(3,l),pnew(4,l))
         call rotbos(0d0,0d0,ba1(1),ba1(2),ba1(3),
     ,    xnew(1,l),xnew(2,l),xnew(3,l),xnew(4,l))
         esum=esum+pnew(4,l)
      enddo
      call leadhad(1,nstring1,1)


c string #2
      if(ms2.le.msmin2) then
        nstring2=1
        l=1
        do j=1,3
          part(j,l)=0d0
          part(j+5,l)=0d0
        enddo
        part(4,l)=ms2
        part(4+5,l)=0d0
        part(5,l)=ms2
        ident(1,l)=i2
        ident(2,l)=iz2
      else
       call qstring(ifdiq4,ifqrk4,ms2,part,ident,nstring2)
       if(nstring2.eq.0) goto 155
      end if

      esum=0d0
       do l=1,nstring2
        pnew(5,nstring1+l)=part(5,l)
        itypnew(nstring1+l)=ident(1,l)

        if(iabs(ident(1,l)).lt.minmes)jbar=jbar+isign(1,ident(1,l))

        i3new(nstring1+l)=ident(2,l)
        do j=1,4
          pnew(j,nstring1+l)=part(j,l)
          xnew(j,nstring1+l)=part(j+5,l)
        enddo
        pnew(3,nstring1+l)=-pnew(3,nstring1+l)
        xnew(3,nstring1+l)=-xnew(3,nstring1+l)
        call rotbos(0d0,0d0,ba2(1),ba2(2),ba2(3),
     ,    pnew(1,nstring1+l),pnew(2,nstring1+l),pnew(3,nstring1+l),
     ,     pnew(4,nstring1+l))
        call rotbos(0d0,0d0,ba2(1),ba2(2),ba2(3),
     ,    xnew(1,nstring1+l),xnew(2,nstring1+l),xnew(3,nstring1+l),
     ,     xnew(4,nstring1+l))
        esum=esum+pnew(4,nstring1+l)
      enddo
      call leadhad(nstring1+1,nstring1+nstring2,1)

c error check
ctp060926      if(ibar.ne.jbar)then
ctp060926         write(6,*)' *** (E) no baryon number conservation', ibar,jbar
ctp060926         write(6,*)'     ',i1,i2,ms1,ms2
ctp060926         write(6,'(5i4)')(itypnew(l),l=1,nstring1+nstring2)
ctp060926         end if


      return


ctp060202 718  format(i2,i4,i3,1x,10(f10.4,1x))



 17   continue

      iz3=iz1
      iz4=iz2
      i3=i1
      i4=i2


      m3=m1
      m4=m2

c the following lines MUST be there in order to set nucleons on-shell
c after their first collision
      if(m3.lt.mminit(i3))m3=mminit(i3)
      if(m4.lt.mminit(i4))m4=mminit(i4)


      goto 2001

 20   continue
c...decays
      if(ityptd(1,pslot(1)).eq.0) then
c normal decay
c note: m4,i4 and iz4 are dummies in this call
         call anndec(1,m1,i1,iz1,m4,i4,iz4,e,sig,gam)
      else
c forward time-delay
         pnew(5,1)=fmasstd(1,pslot(1))
         itypnew(1)=ityptd(1,pslot(1))
         i3new(1)=iso3td(1,pslot(1))
         pnew(5,2)=fmasstd(2,pslot(1))
         itypnew(2)=ityptd(2,pslot(1))
         i3new(2)=iso3td(2,pslot(1))
      endif

      if(nexit.eq.2) then
         i3=itypnew(1)
         iz3=i3new(1)
         m3=pnew(5,1)
         i4=itypnew(2)
         iz4=i3new(2)
         m4=pnew(5,2)
         goto 2001
      else
c three or four body decay
         nstring1=1
         nstring2=nexit-1
         do i=1,4
            do j=1,nexit
               pnew(i,j)=0d0
               xnew(i,j)=0d0
            end do
         end do
         mstring(1)=pnew(5,1)
c
         mstring(2)=pnew(5,2)
         do 91 j=3,nexit
            mstring(2)=mstring(2)+pnew(5,j)
 91      continue
c
c now call routine for momentum phase space...
         call nbodydec(e)
         return
      endif

 23   continue
c...annihilation -> string
      if(CTOption(12).ne.0)then
        write(6,*)' *** error(make22): string section is called ',
     .            'while strings are switched off:CTOption(12).ne.0'
        stop
      end if

      ms1=e/2.
      ms2=e/2.
      mstring(1)=ms1
      mstring(2)=ms2

c determine flavour content of b-bbar-system
      call ityp2id(i1,iz1,ifdiq1,ifqrk1)
      call ityp2id(i2,iz2,ifdiq2,ifqrk2)
c...create string 1 out of quark-antiquark pair
      call qstring(ifqrk1,ifqrk2,ms1,part,ident,nstring1)
      esum=0d0

      do k=1,nstring1
        l=k
        do j=1,4
          pnew(j,l)=part(j,l)
          xnew(j,l)=part(j+5,l)
        enddo
        esum=esum+pnew(4,l)
        pnew(5,l)=part(5,l)
        itypnew(l)=ident(1,l)
        i3new(l)=ident(2,l)
        tau=part(9,l)/ (part(4,l)/part(5,l))
      enddo
      call leadhad(1,nstring1,0)

c...create string 2 out of diquark-antidiquark-pair
c   use one quark and one antiquark for the string-ends:
       ifqrk1=int(ifdiq1/1000)
       ifqrk2=int(ifdiq2/1000)
c...store remaining flavour quantum numbers in ctp(26), they will
c   be passed to the 'clustr'-routine:
       ifdiq1=mod(ifdiq1/100,10)
       ifdiq2=mod(ifdiq2/100,10)
       valint(1)=dble(((abs(ifdiq1*10.d0)+abs(ifdiq2*1.d0))
     &            *isign(1,ifdiq1)))
       valint(1)=sign(valint(1),dble(ifdiq1))
       call qstring(ifqrk1,ifqrk2,ms2,part,ident,nstring2)
       valint(1)=0.d0
      esum=0d0

      do l=1,nstring2
        do j=1,4
          pnew(j,nstring1+l)=part(j,l)
          xnew(j,nstring1+l)=part(j+5,l)
        enddo
        esum=esum+pnew(4,nstring1+l)
        pnew(5,nstring1+l)=part(5,l)
        itypnew(nstring1+l)=ident(1,l)
        i3new(nstring1+l)=ident(2,l)
        tau=part(9,nstring1+l)/ (part(4,l)/part(5,l))
      enddo
      call leadhad(nstring1+1,nstring1+nstring2,0)

      return

 26   continue
c...elastic MB scattering (the outgoing particle id's must not be
c   assigned randomly like at label 2001)
c
      if(switips)then
         call swpizm(i1,iz1,m1,i2,iz2,m2)
         switips=.false.
      endif

      call setizm(i1,iz1,m1,i2,iz2,m2,i3,iz3,m3,i4,iz4,m4)

      if(mminit(i4)+mminit(i3).gt.e)then
ctp060926         write(6,*)'make22(el):threshold violated'
ctp060926         write(6,*)'m3:',m3,mminit(i3)
ctp060926         write(6,*)'m4:',m4,mminit(i4)
      else
        if(m3.lt.mminit(i3))m3=mminit(i3)
        if(m4.lt.mminit(i4))m4=mminit(i4)
      end if

c... get momenta & fill newpart, 2 particle exit-channel

      nstring1=1
      nstring2=1
      nexit=2
      do i=1,4
         do j=1,2
            pnew(i,j)=0d0
            xnew(i,j)=0d0
         end do
      end do

c...boost to 2-particle cms

      pnew(3,1)=pcms(e,m3,m4)
      pnew(3,2)=-pcms(e,m3,m4)

      pnew(4,1)=sqrt(m3**2+pnew(3,1)**2)
      pnew(4,2)=sqrt(m4**2+pnew(3,2)**2)

      pnew(5,1)=m3
      mstring(1)=m3
      itypnew(1)=i3
      i3new(1)=iz3

      pnew(5,2)=m4
      mstring(2)=m4
      itypnew(2)=i4
      i3new(2)=iz4

      return

 27   continue
c XX-> 1 string : e+e- , MB
      if(CTOption(12).ne.0)then
        write(6,*)' *** error(make22): string section is called ',
     .            'while strings are switched off:CTOption(12).ne.0'
        stop
      end if

c quark-quark scattering -> only elastic !
      if(xfac1.lt..999d0.and.xfac2.lt..999d0
     &   .and.ranf(0).lt.0.25) goto 26

      ms1=e
      mstring(1)=ms1
      mstring(2)=0d0

c determine flavour content of string
      if(CTOption(20).eq.1)then
c..e+e- annihilation
      if(ranf(0).lt.CTParam(6))then
        ifqrk1=3  ! ssbar
        ifqrk2=-3
      else
        call ityp2id(104,0,ifqrk1,ifqrk2)  ! qqbar
      end if
      else
c...MB annihilation. the quark content must be known:
      call ityp2id(i2,iz2,iq1(1),iq1(2))
      call ityp2id(i1,iz1,ifdiq2,iq2(3))

      if(abs(i1).ge.minmes) then
         iq2(1)=ifdiq2
         iq2(2)=iq2(3)
         iq2(3)=0
      else
         iq2(1)=mod(ifdiq2/100,10)
         iq2(2)=int(ifdiq2/1000)
      endif

      do 312 kq1=1,2
       do 412 kq2=1,3
c..two of the quarks must be able to annihilate:
        if(iq1(kq1)+iq2(kq2).eq.0) then
         kqq1=kq1
         kqq2=kq2
         goto 414
        endif
 412   continue
 312  continue
      goto 26 ! could not create double charged strange baryon string
 414  continue
c.. the 'iff'-quarks constitute the produced (anti-)baron
      iff(1)=iq1(3-kqq1)
      iii=1
      do 512 kq2=1,3
       if (kq2.ne.kqq2) then
        iii=iii+1
        iff(iii)=iq2(kq2)
       endif
 512  continue

      if(abs(i1).lt.minmes) then
         call mquarks(iff,ifqrk1,ifqrk2)
      else
         ifqrk1=iff(1)
         ifqrk2=iff(2)
      endif

      endif

c...create string 1 out of quark-antiquark pair
      call qstring(ifqrk1,ifqrk2,ms1,part,ident,nstring1)
      esum=0d0

c primitive bug-fix:
      if(nstring1.eq.0)then
ctp060926        write(6,*)'make22: iline 27 not completed. ->elastic'
        goto 26
      endif

      do 101 k=1,nstring1
        l=k
c no leading hadron in e+e-
        if(CTOption(20).eq.1)then
          leadfac(l)=0.0d0
        endif
        do 102 j=1,4
          pnew(j,l)=part(j,l)
          xnew(j,l)=part(j+5,l)
 102    continue
        esum=esum+pnew(4,l)
        pnew(5,l)=part(5,l)
        itypnew(l)=ident(1,l)
        i3new(l)=ident(2,l)
        tau=part(9,l)/ (part(4,l)/part(5,l))
 101    continue
      if(CTOption(20).eq.1)then
        call leadhad(1,nstring1,3)
      else
        call leadhad(1,nstring1,1)
      endif

      nstring2=0
      return


 29   continue
c...DD->ND detailed balance
      i3=minnuc
      i4=mindel
      m3=massit(i3)
      m4=getmass(e-m3,0)
      if(iabs(iz1+iz2).gt.isoit(i3)+isoit(i4))then
       iz3=-9
       iz4=-9
      else
         nexit=2
         itot(1)=isoit(i3)
         itot(2)=isoit(i4)
         call isocgk4(isoit(i1),iz1,isoit(i2),iz2,itot,i3new,errflg)
         i3=isign(i3,ii1)
         i4=isign(i4,ii2)
         iz3=i3new(1)
         iz4=i3new(2)
      end if

      if(ranf(0).gt.0.5d0)call swpizm(i3,iz3,m3,i4,iz4,m4)
      goto 2001


 100  continue
c...??->NN detailed balance inverse channels
c iso3 are assigned in detbal
      nexit=2
      i3=minnuc
      i4=minnuc
      m3=massit(minnuc)
      m4=massit(minnuc)
      itot(1)=isoit(i3)
      itot(2)=isoit(i4)
      call isocgk4(isoit(i1),iz1,isoit(i2),iz2,itot,i3new,errflg)
      i3=isign(i3,ii1)
      i4=isign(i4,ii2)
      iz3=i3new(1)
      iz4=i3new(2)

      if(ranf(0).gt.0.5d0)call swpizm(i3,iz3,m3,i4,iz4,m4)
      goto 2001


 36   continue
c...MB->B',MM->M* annihilations (forward time delay)
      call anndec(0,m1,i1,iz1,m2,i2,iz2,e,sig,gam)
      goto 2003

 1008 continue
c...get isospin-3 components

      nexit=2
      itot(1)=isoit(i3)
      itot(2)=isoit(i4)
      call isocgk4(isoit(i1),iz1,isoit(i2),iz2,itot,i3new,errflg)
      iz3=i3new(1)
      iz4=i3new(2)

ctp060926      if(errflg.ne.0)then
ctp060926        write(6,*)'make22: iso-spin conservation ',
ctp060926     ,            'not possible in isocgk: error-flag=',errflg
ctp060926        write(6,*)'      ',isoit(i1),iz1,isoit(i2),iz2,'>',
ctp060926     ,            isoit(i3),isoit(i4),iz3,iz4,
ctp060926     ,            ' process:',e,i1,m1,i2,m2,'>',i3,m3,i4,m4,'io=',io
ctp060926      end if
      i3=isign(i3,ii1)
      i4=isign(i4,ii2)

 2001 continue


c... get momenta & fill newpart, 2 particle exit-channel

      nstring1=1
      nstring2=1
      nexit=2
      do i=1,4
         do j=1,2
            pnew(i,j)=0d0
            xnew(i,j)=0d0
         end do
      end do

c...boost to 2-particle cms


      if(.not.(ityptd(1,pslot(1)).ne.0.and.CTOption(34).eq.2.and.
     &     iline.eq.20)) then
c normal decay

         pnew(3,1)=pcms(e,m3,m4)
         pnew(3,2)=-pcms(e,m3,m4)

         pnew(4,1)=sqrt(m3**2+pnew(3,1)**2)
         pnew(4,2)=sqrt(m4**2+pnew(3,2)**2)
      else
c forward time delay
         pnew(1,1)=pxtd(1,pslot(1))
         pnew(1,2)=pxtd(2,pslot(1))
         pnew(2,1)=pytd(1,pslot(1))
         pnew(2,2)=pytd(2,pslot(1))
         pnew(3,1)=pztd(1,pslot(1))
         pnew(3,2)=pztd(2,pslot(1))
         pnew(4,1)=p0td(1,pslot(1))
         pnew(4,2)=p0td(2,pslot(1))
      endif

      pnew(5,1)=m3
      mstring(1)=m3
      itypnew(1)=i3
      i3new(1)=iz3

      pnew(5,2)=m4
      mstring(2)=m4
      itypnew(2)=i4
      i3new(2)=iz4

      return

 2002 continue
c...  fill newpart, one particle exit channel
      nstring1=1
      nstring2=0
      nexit=1
      do i=1,4
        pnew(i,1)=0d0
        xnew(i,1)=0d0
      end do
      pnew(4,1)=e
      mstring(1)=e
c     the rest of the relevant new particle data have been
c     filled into the newpart arrays by anndex

      return

 2003 continue
c bookkeeping for forward time-delay
      do 204 j=1,2
         if(pslot(j).lt.1) goto 204
         pold(1,j)=px(pslot(j))
         pold(2,j)=py(pslot(j))
         pold(3,j)=pz(pslot(j))
         pold(4,j)=p0(pslot(j))
         pold(5,j)=fmass(pslot(j))
         itypold(j)=ityp(pslot(j))
         iso3old(j)=iso3(pslot(j))
 204  continue

c...  fill newpart, one particle exit channel
      nstring1=1
      nstring2=0
      nexit=1
      do i=1,4
        pnew(i,1)=0d0
        xnew(i,1)=0d0
      end do
      pnew(4,1)=e
      mstring(1)=e

      return
      end


C####C##1#########2#########3#########4#########5#########6#########7##
      real*8 function rnfxf(b,xm,xp)
c
cinput b  : parameter
cinput xm : lower interval boundary
cinput xp : upper interval boundary
c
c {\tt rnfxf} yields a value $x$ between {\tt xm} and {\tt xp}
c distributetd like $(1-x)^b/x$.
c
cccccCcc1ccccccccc2ccccccccc3ccccccccc4ccccccccc5ccccccccc6ccccccccc7cc
      implicit none
      real*8 xm,xp,i,inv,x,fmin,fmax,f,ranf,b

      i(x)=log(x)
      inv(x)=exp(x)

      fmin=i(xm)
      fmax=i(xp)
 3    f=fmin+(fmax-fmin)*ranf(0)
      x=inv(f)
      if(ranf(0).gt.(1d0-x/xp)**b)goto 3
      rnfxf=x

      return
      end

C####C##1#########2#########3#########4#########5#########6#########7##
      subroutine setizm(i1,iz1,m1,i2,iz2,m2,i3,iz3,m3,i4,iz4,m4)
c
c input  : {\tt i1,iz1,m1,i2,iz2,m2}
c output : {\tt i3,iz3,m3,i4,iz4,m4}
c
c This subroutine simply maps {\tt iz1} $\to$ {\tt iz3} $\ldots$
c {\tt m2} $\to$ {\tt m4}
c
cccccCcc1ccccccccc2ccccccccc3ccccccccc4ccccccccc5ccccccccc6ccccccccc7cc
      real*8 m1,m2,m3,m4
      integer i1,iz1,i2,iz2,i3,iz3,i4,iz4
      i3=i1
      i4=i2
      m3=m1
      m4=m2
      iz3=iz1
      iz4=iz2
      return
      end

C####C##1#########2#########3#########4#########5#########6#########7##
      subroutine swpizm(i1,iz1,m1,i2,iz2,m2)
c
c
c input   : {\tt i1,iz1,m1,i2,iz2,m2}
c output  : {\tt i1,iz1,m1,i2,iz2,m2}
c
c This subroutine simply swaps {\tt 1} $\to$ {\tt 2} and vice versa.
c
c
cccccCcc1ccccccccc2ccccccccc3ccccccccc4ccccccccc5ccccccccc6ccccccccc7cc
      real*8 m1,m2,ms
      integer i1,iz1,i2,iz2,is,izs
      is=i1
      i1=i2
      ms=m1
      m1=m2
      izs=iz1
      iz1=iz2
      i2=is
      iz2=izs
      m2=ms
      return
      end


C####C##1#########2#########3#########4#########5#########6#########7##
      real*8 function rnfpow(n,mmin,mmax)
c
cinput n    : parameter
cinput mmin : lower interval boundary
cinput mmax : upper interval boundary
c
c {\tt rnfpow} yields a value $x$ between {\tt mmin} and {\tt mmax}
c distributetd like $x^n$.
c
cccccCcc1ccccccccc2ccccccccc3ccccccccc4ccccccccc5ccccccccc6ccccccccc7cc
      implicit none
      real*8 mmin,mmax,i,inv,ii,iinv,x,fmin,fmax,f,ranf
      integer n

      i(x)=x**(n+1)/(n+1)
      inv(x)=((n+1)*x)**(1/(n+1))
      ii(x)=log(x)
      iinv(x)=exp(x)

      if(n.eq.-1)then
        fmin=ii(mmin)
        fmax=ii(mmax)
        f=fmin+(fmax-fmin)*ranf(0)
        rnfpow=iinv(f)
      else
        fmin=i(mmin)
        fmax=i(mmax)
        f=fmin+(fmax-fmin)*ranf(0)
        rnfpow=inv(f)
      end if
      return
      end

C####C##1#########2#########3#########4#########5#########6#########7##
      real*8 function fmsr(mmin,mmax)
c
cinput mmin : minimum mass
cinput mmax : maximum mass
c
c  {\tt fmsr} yields a mass according to the finite mass sum rule (FMSR)
c  A.I.Sanda Phys.~Rev.~{\bf D6}~(1973)~231 and
c  M.B.Einhorn  et al., Phys. Rev. {\bf D5}~(1972)~2063
c
cccccCcc1ccccccccc2ccccccccc3ccccccccc4ccccccccc5ccccccccc6ccccccccc7cc
      implicit none
      real*8 mmin,mmax,rnfpow!,x
      ! i(x)=log(x)
      ! inv(x)=exp(x)

      fmsr=rnfpow(-1,mmin,mmax)
      return
      end

C####C##1#########2#########3#########4#########5#########6#########7##
      subroutine getfmsr(io,mmin,mmax,ir,mr)
c
cinput io   : class of resonance
cinput mmin : minimum mass
cinput mmax : maximum mass
coutput ir  : ID of resonance
coutput mr  : mass of resonance
c
c  {\tt getfmsr} is an extension of {\tt fmsr} to simulate
c  a resonance structure at low masses
c
cccccCcc1ccccccccc2ccccccccc3ccccccccc4ccccccccc5ccccccccc6ccccccccc7cc

      implicit none
      integer io,ir,whichres
      real*8 mmin,mmax,mr,massit,widit
      include 'comres.f'
      real*8 fmsr,mmean

      mr=fmsr(mmin,mmax)
c...if strings come into play > modify whichres!
      ir=whichres(mr,io)
      mr=mmean(1,massit(ir),widit(ir),mmin,mmax)
      return
      end

C####C##1#########2#########3#########4#########5#########6#########7##
      SUBROUTINE STREXCT(IFL11,IFL12,IB1,AM1,IFL21,IFL22,IB2,AM2,ECM,
     *IOPT,ba1,ams1,ba2,ams2,ifl31,ifl32,ifl41,ifl42)
c
cinput ifl11  : ID of (di)quark of projectile hadron
cinput ifl12  : ID of remaining quark of projectile hadron
cinput ib1    : baryon number of projectile hadron
cinput am1    : mass of projectile hadron
cinput ifl21  : ID of (di)quark of target hadron
cinput ifl22  : ID of remaining quark of target hadron
cinput ib2    : baryon number of target hadron
cinput am2    : mass of target hadron
cinput ecm    : $\sqrt{s}$ of excitation
cinput iopt   : flag for excitation ansatz (1: Fritiof, 2:QGSM)
coutput ba1   : velocity of 1st string
coutput ams1  : mass of 1st string
coutput ba2   : velocity vector of 2nd string
coutput ams2  : mass of 2nd string
coutput ifl31 : (di)quark content of 1st string
coutput ifl32 : remainig quark content of 1st string
coutput ifl41 : (di)quark content of 2nd string
coutput ifl42 : remaining quark content of  2nd string
c
c output : particle IDs, masses and  momenta via common block {\tt newpart}
c
c This routine performs the excitation of two strings according
c to different ansatzes defined via parameter {\tt iopt}
c
cccccCcc1ccccccccc2ccccccccc3ccccccccc4ccccccccc5ccccccccc6ccccccccc7cc
      implicit real*8 (a-h,o-z)
      implicit integer (i-n)
      include 'options.f'

      real*8 ba1(3),ba2(3),P(2,5),ranf
      integer IFLS(2,2)

      PARAMETER(PI=3.14159)! value of Pi
      PARAMETER(XMAX=1.0)
      LOGICAL SPINT,bt

      sigma=CTParam(31)
      alpha=CTParam(32)
      betav=CTParam(33)
      ampit=CTParam(34)
      betas=CTParam(40)
c
C   CHECK INITIAL ENERGY
      IF(AM1+AM2+ampit.gt.ECM) THEN
         WRITE(6,*) '....STOP! Initial energy is too low for string
     * excitation, ECM=', ECM,'<',am1,'+',am2
         STOP
      ENDIF

      bt=.false.
      if(ranf(0).lt.5d-1)bt=.true.

C   COMPUTE XMIN
      XMIN=AMPIT/ECM

C   COMPUTE C.M.S. HADRON MOMENTUM AND ENERGIES
      PCMS=SQRT(ALAMB(ECM**2,AM1**2,AM2**2))/(2.0*ECM)
      ECMS1=SQRT(PCMS**2+AM1**2)
      ECMS2=SQRT(PCMS**2+AM2**2)

C   COMPUTE LIGHT-CONE VARIABLES FOR HADRON-PROJECTILE :
C                                           (PPLS01,PMINS01,PTX01,PTY01)
      PPLS01=ECMS1+PCMS
      PMINS01=ECMS1-PCMS
      PTX01=0.
      PTY01=0.

C   COMPUTE LIGHT-CONE VARIABLES FOR HADRON-TARGET :
C                                           (PPLS02,PMINS02,PTX02,PTY02)
      PPLS02=ECMS2-PCMS
      PMINS02=ECMS2+PCMS
      PTX02=0.
      PTY02=0.

C   COMPUTE TRANSFERRED TRANSVERSE MOMENTUM :
C                                           (QX,QY)
100   CALL GAUSPT(QT,SIGMA)
      PHI=2.*PI*ranf(0)
      QX=QT*COS(PHI)
      QY=QT*SIN(PHI)

C   COMPUTE XPLUS AND XMINUS VALUES FOR INTERACTING PARTONS
c      IF(IOPT.EQ.1) THEN
       if(IOPT.eq.1.or.dabs(alpha).lt.1d-4)then
         CALL XSDIS(XPLUS,XMIN,XMAX,BETAS)
         CALL XSDIS(XMINUS,XMIN,XMAX,BETAS)
       else ! IF(IOPT.EQ.2) THEN
         XPLUS=XVDIS(XMIN,ALPHA,BETAV)
         XMINUS=XVDIS(XMIN,ALPHA,BETAV)
c      else
c        write(6,*)'strexct: undefined excitation model'
c        stop
      ENDIF
C  COMPUTE LIGHT-CONE PARAMETERS OF INTERACTING PARTONS :
C                                      (PPRPLS1,PPRMINS1,PPRX1,PPRY1)
      PPRPLS1=XPLUS*PPLS01
      PPRMINS1=0.
      PPRX1=0.
      PPRY1=0.
C
C                                      (PPRPLS2,PPRMINS2,PPRX2,PPRY2)
      PPRPLS2=0.
      PPRMINS2=XMINUS*PMINS02
      PPRX2=0.
      PPRY2=0.
C
C  COMPUTE LIGHT-CONE COMPONENT OF TRANSFERRED MOMENTUM :
C                                                       (QPLS,QMINS)
      QPLS=-(QT**2/PPRMINS2)
      QMINS= QT**2/PPRPLS1
C
C  COMPUTE LIGHT-CONE PARAMETERS AND QUARK CONTENTS FOR EXCITED STRINGS
      IF(IOPT.EQ.1) THEN
C  FRITIOF ANSATZ :
C  P1 ---> P1 + Q
C  P2 ---> P2 - Q
         PPLS1=PPLS01+QPLS
         PMINS1=PMINS01+QMINS
         PTX1=PTX01+QX
         PTY1=PTY01+QY
C
         PPLS2=PPLS02-QPLS
         PMINS2=PMINS02-QMINS
         PTX2=PTX02-QX
         PTY2=PTY02-QY
C
C  NO QUARK REARRANGEMENT
         IFLS(1,1)=IFL11
         IFLS(1,2)=IFL12
         IFLS(2,1)=IFL21
         IFLS(2,2)=IFL22
      ELSEIF(IOPT.EQ.2) THEN
C  QGSM ANSATZ :
C  P1 ---> P1 - PPR1 + PPR2 - Q
C  P2 ---> P2 - PPR2 + PPR1 + Q
C
         PPLS1=PPLS01-PPRPLS1+PPRPLS2-QPLS
         PMINS1=PMINS01-PPRMINS1+PPRMINS2-QMINS
         PTX1=PTX01-PPRX1+PPRX2-QX
         PTY1=PTY01-PPRY1+PPRY2-QY
C
         PPLS2=PPLS02-PPRPLS2+PPRPLS1+QPLS
         PMINS2=PMINS02-PPRMINS2+PPRMINS1+QMINS
         PTX2=PTX02-PPRX2+PPRX1+QX
         PTY2=PTY02-PPRY2+PPRY1+QY

C  QUARK REARRANGEMENT IS NEEDED TO CREATE COLOUR NEUTRAL STRINGS
C  (In case of Baryon-Baryon or Meson-Meson or
C  Antibaryon-Antibaryon or Meson-Antibaryon or Antibaryon-Meson
C  interaction)
         IF((IB1.EQ.1.AND.IB2.EQ.1).OR.
     .        (IB1.EQ.-1.AND.IB2.EQ.-1).OR.
     .        (IB1.EQ.0.AND.IB2.EQ.0).OR.
     .        (IB1.EQ.0.AND.IB2.EQ.-1).OR.
     .        (IB1.EQ.-1.AND.IB2.EQ.0)) THEN
            IFLS(1,1)=IFL11
            IFLS(1,2)=IFL22
            IFLS(2,1)=IFL21
            IFLS(2,2)=IFL12
         ENDIF

C  QUARK REARRANGEMENT  (In case of Meson-Baryon or Antibaryon-Baryon or v.v.
C  interaction)
         IF((IB1.EQ.0.AND.IB2.EQ.1).OR.
     .        (IB1.EQ.1.AND.IB2.EQ.0).OR.
     .        (IB1.EQ.-1.AND.IB2.EQ.1).OR.
     .        (IB1.EQ.1.AND.IB2.EQ.-1)) THEN
            IFLS(1,1)=IFL11
            IFLS(1,2)=IFL21
            IFLS(2,1)=IFL22
            IFLS(2,2)=IFL12
         ENDIF
      ELSE
         PPLS1=0
         PMINS1=0
         PTX1=0
         PTY1=0
C
         PPLS2=0
         PMINS2=0
         PTX2=0
         PTY2=0
C
         WRITE(6,*) 'ERROR IN MAKE22: WRONG IOPT'
         STOP
      ENDIF

C  COMPUTE OUTGOING STRING MASSES. They should be more than stable hadron
C  (with the same quark content) masses. In the contrary case generation
C  should be repeated
      AMS1S=PPLS1*PMINS1-QT**2
      SPINT=.TRUE.
      AMS0=0
C  COMPUTE AMS0
      IF(MOD(IFLS(1,1),100).EQ.0.AND.MOD(IFLS(1,2),100).EQ.0) THEN
C  qq-qqbar string
       IFLU=1  ! add u- quark to construct hadron
C  AMS0 is sum of masses of lowest baryon states
       IFLH1=IDPARS(IFLS(1,1), -ISIGN(IFLU,IFLS(1,1)),SPINT,2)
       IFLH2=IDPARS(IFLS(1,2), -ISIGN(IFLU,IFLS(1,2)),SPINT,2)
       AMS0=amass(IFLH1)+amass(IFLH2)
      ENDIF
      IF(.NOT.(MOD(IFLS(1,1),100).EQ.0
     .    .AND.MOD(IFLS(1,2),100).EQ.0))THEN
C  AMS0 is mass of lowest hadron state
        IKH=IDPARS(IFLS(1,1),IFLS(1,2),SPINT,2)
        AMS0=amass(IKH)
      ENDIF
C
      IF(AMS1S.LT.AMS0**2) GO TO 100
C
      AMS2S=PPLS2*PMINS2-QT**2
      SPINT=.TRUE.
C  COMPUTE AMS0
      IF(MOD(IFLS(2,1),100).EQ.0.AND.MOD(IFLS(2,2),100).EQ.0) THEN
C  qq-qqbar string
       IFLU=1  ! add u- quark to construct hadron
C  AMS0 is sum of masses of lowest baryon states
       IFLH1=IDPARS(IFLS(2,1),-ISIGN(IFLU,IFLS(2,1)),SPINT,2)
       IFLH2=IDPARS(IFLS(2,2),-ISIGN(IFLU,IFLS(2,2)),SPINT,2)
       AMS0=amass(IFLH1)+amass(IFLH2)
      ENDIF
      IF(.NOT.(MOD(IFLS(2,1),100).EQ.0
     .    .AND.MOD(IFLS(2,2),100).EQ.0)) THEN
C  AMS0 is mass of lowest hadron state
        IKH=IDPARS(IFLS(2,1),IFLS(2,2),SPINT,2)
        AMS0=amass(IKH)
      ENDIF
C
      IF(AMS2S.LT.AMS0**2) GO TO 100
C
      AMS1=SQRT(AMS1S)
      AMS2=SQRT(AMS2S)
C
C SUM OF MASSES OF EXCITED STRINGS SHOULD BE LESS THAN INITIAL ENERGY
      IF(AMS1+AMS2.GT.ECM) GO TO 100
C
      P(1,1)=PTX1
      P(1,2)=PTY1
      P(1,3)=0.5*(PPLS1-PMINS1)
      P(1,4)=0.5*(PPLS1+PMINS1)
      P(1,5)=AMS1
      P(2,1)=PTX2
      P(2,2)=PTY2
      P(2,3)=0.5*(PPLS2-PMINS2)
      P(2,4)=0.5*(PPLS2+PMINS2)
      P(2,5)=AMS2
c
      do 152 j=1,3
      ba1(j)=p(1,j)/p(1,4)
 152  ba2(j)=p(2,j)/p(2,4)

      ifl31=ifls(1,1)
      ifl32=ifls(1,2)
      ifl41=ifls(2,1)
      ifl42=ifls(2,2)

C
C  CHECK ENERGY-MOMENTUM CONSERVATION
      ESTR=P(1,4)+P(2,4)
      PSTRX=PTX1+PTX2
      PSTRY=PTY1+PTY2
      PSTRZ=P(1,3)+P(2,3)
C
ctp060926      if(abs(estr-ecm).gt.1d-12)then
ctp060926           WRITE(6,*) 'ECM=',ECM
ctp060926           WRITE(6,*) 'ESTR=',ESTR
ctp060926      end if
ctp060926      if(pstrx.gt.1d-12)WRITE(6,*) 'PSTRX=',PSTRX
ctp060926      if(pstry.gt.1d-12)WRITE(6,*) 'PSTRY=',PSTRY
ctp060926      if(pstrz.gt.1d-12)WRITE(6,*) 'PSTRZ=',PSTRZ
C Note! One should also check flavour conservation!
C
      RETURN
      END
C####C##1#########2#########3#########4#########5#########6#########7##
      real*8 FUNCTION ALAMB(X,Y,Z)
c
c input :  {\tt x,y,z}
c
C   THIS ROUTINE COMPUTES KINEMATICAL FUNCTION
cccccCcc1ccccccccc2ccccccccc3ccccccccc4ccccccccc5ccccccccc6ccccccccc7cc
      implicit none
      real*8 x,y,z

      ALAMB=max(0d0,(X-Y-Z)*(X-Y-Z) - 4.D0*Y*Z)

      RETURN
      END

C####C##1#########2#########3#########4#########5#########6#########7##
      real*8 FUNCTION XVDIS(XMIN,ALFA,BETA)
c
cinput xmin : lower boundary
cinput alfa : parameter
cinput beta : parameter
c
C This function returns {\tt xmin}$ < x < 1$
c values distributed according to the Beta-function
C
cccccCcc1ccccccccc2ccccccccc3ccccccccc4ccccccccc5ccccccccc6ccccccccc7cc
      implicit none
      real*8 x,xmin,alfa,beta

 100  CALL SBETA(X,ALFA,BETA)
      IF(X.LE.XMIN) GO TO 100
      XVDIS=X
      RETURN
      END

C####C##1#########2#########3#########4#########5#########6#########7##
      SUBROUTINE SBETA(X,ALFA,BETA)
c
coutput x   : Beta-distributed value
cinput alfa : parameter
cinput beta : parameter

C THIS ROUTINE GENERATES X ACCORDING TO BETA DISTRIBUTION
C   $U(X)=C*X**(ALFA-1)*(1-X)**(BETA-1)$
C IONK,S METHOD IS USED
c
cccccCcc1ccccccccc2ccccccccc3ccccccccc4ccccccccc5ccccccccc6ccccccccc7cc
      implicit none
      real*8 ran1,ran2,alfa,beta,r1a,r2b,r12,x,ranf

 1    RAN1=ranf(0)
      RAN2=ranf(0)
      R1A=RAN1**(1./ALFA)
      R2B=RAN2**(1./BETA)
      R12=R1A+R2B
      IF(R12.GE.1.) GO TO 1
      X=R1A/R12
      RETURN
      END



C####C##1#########2#########3#########4#########5#########6#########7##
      SUBROUTINE XSDIS(X,XMIN,XMAX,BETA)
c
cinput xmin : lower boundary
cinput xmax : upper boundary
cinput beta : parameter
c
C  THIS FUNCTION GENERATES $XMIN < X < XMAX$ ACCORDING TO
C DISTRIBUTION $U(X)= 1./X*(1.-X)**(BETA+1)$ DISTRIBUTION
C PARAMETER $BETA > 0$
c
cccccCcc1ccccccccc2ccccccccc3ccccccccc4ccccccccc5ccccccccc6ccccccccc7cc
      implicit real*8 (a-h,o-z)
      implicit integer (i-n)

      PE1=(1.-XMIN)**(BETA+1.)/(BETA+1.)-(1.-XMAX)**(BETA+1.)/
     *(BETA+1.)
      PE=PE1+DLOG(XMAX/XMIN)
      PE1=PE1/PE
 108  RND=ranf(0)
      RNDMPE1=ranf(0)
      IF(RNDMPE1.GT.PE1) GO TO 200
      X=1.-((1.-XMIN)**(BETA+1.)*(1.-RND)+(1.-XMAX)**(BETA+1.)*
     *RND)**(1./(BETA+1.))
      GO TO 300
200   X=XMIN*(XMAX/XMIN)**RND
300   PPE1=(1.-X)**BETA
      PPE2=1./X
      PPE1=PPE1+PPE2
      PPE2=PPE1*PPE2
      IF(PPE1*ranf(0).GT.PPE2) GO TO 108
      RETURN
      END


C####C##1#########2#########3#########4#########5#########6#########7##
      subroutine getres(io,e,im,ip,ir)
c
cinput io : resonance class
cinput e  : $\sqrt{s}$ of excitation process
cinput im : lower boundary for resonance ID's
cinput ip : upper boundary for resonance ID's
coutput ir : ID of resonance
c
c This subroutine randomly selects a resonance ID for a given
c resonance class ({\tt io}), collision $\sqrt{s}$ and ID possible range.
c The selected resonance choice is weighted according to the respective
c partial cross section for the excitation of such a resonance.
c
C####C##1#########2#########3#########4#########5#########6#########7##
      implicit none
      include 'comres.f'
      include 'options.f'
      integer io,i,im,ip,ir
      real*8 e,x(minnuc:maxmes),xmax

      do i=im,ip
         call crossf(io,e,i,x(i))
      end do
      call getbran(x,minnuc,maxmes,xmax,im,ip,ir)

ctp060926      if(ir.gt.ip.or.ir.lt.im) then
ctp060926         write(6,*)'***(E) getres: no final state selected...'
ctp060926         write(6,*)io,e,xmax,im,ip,ir
ctp060926      endif

      return
      end

C####C##1#########2#########3#########4#########5#########6#########7##
      subroutine aqm(i1,i2,xt,xel)
c
cinput :  {\tt i1,i2} : ID's of particle 1 and 2
coutput : {\tt xt,xel}: total and elastic cross section
c
c This subroutine returns cross sections according
c to the additive quark model.
c
C####C##1#########2#########3#########4#########5#########6#########7##
      integer i(2),ij,j,i1,i2,strit
      real*8 s(2),xt,xel,mn(2)
      integer ifa, ifb, isoit
      include 'comres.f'

      i(1)=i1
      i(2)=i2
      do 108 j=1,2
        ij=i(j)
        if(iabs(ij).ge.minmes)then
          s(j)=1d0*abs(strit(ij))
          call ityp2id(ij, isoit(ij), ifa, ifb)
c take care of hidden strangeness
          if (abs(ifa).eq.3 .and. abs(ifb).eq.3) s(j)=2
          mn(j)=1.d0
        else
          s(j)=1d0*abs(strit(ij))
          mn(j)=0.d0
        end if
 108  continue

      xt=max(0d0,40d0*0.666667**(mn(1)+mn(2))
     *  *( 1d0-0.4*s(1)/(3d0-mn(1)) )
     *  *( 1d0-0.4*s(2)/(3d0-mn(2)) ) )
      xel=0.039*xt**1.5d0
      return
      end

C####C##1#########2#########3#########4#########5#########6#########7##
      real* 8 function fppfit(iio,e,i3,i4)
c
c     Version: 1.0
c
cinput iio    : tag for cross section class
cinput e     : $\sqrt{s}$ of collision
cinput i3    : ID of first outgoing particle
cinput i4    : ID of second outgoing particle
c
c
c     {\tt fppfit} returns the isospin-independent part
c     of the production cross section
c     for one or two outgoing resonances ({\tt i3} and {\tt i4})
c     in a proton-proton
c     collision. {\tt io} sets the class of cross section
c     which is returned (i.e. $p p \rightarrow N \Delta$). If
c        {\tt io} is set 99 the class will be determined according
c        to {\tt i3} and {\tt i4}.
c
C####C##1#########2#########3#########4#########5#########6#########7##
      implicit none
        include 'comres.f'
      include 'comwid.f'
      include 'options.f'

      real*8 e,tmp
      integer i,iio,io,im,jm,id1,id2,id3,id4,i3,i4
        integer i1m,i1p,i2m,i2p,class1,class2,iim
c funct.
      real*8 pmean,widit,massit,splint
      integer jit,isoit

c...parameters for resonance production in NN
      integer nfit
      parameter (nfit=7)
      real*8 ar(nfit)
      integer rr(2,nfit)
c io =         1      2     3     4       5      6      7
c              nd    nn*    nd*   dd      dn*    dd*    B*B*
      data ar/ 4d4,  6.3d0, 12d0, 2.8d0,   3.5d0, 3.5d0, 0.d0 /
      data rr/ 4,0,   4,1,   4,2,  0,0,   0,1,  0,2,    3,3 /
c rr tells of which particle class the out particles are
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc

c hard cut: no resonant cross section at highest energies
        if(e.gt.maxtab2)then
          fppfit=0d0
          return
        endif
c

      id1=0
      id2=0
      id3=0
      id4=0
      im=i3
      jm=i4

        if(iio.eq.99)then
c determine iline

c get class of 1 particle
          do 101 i=0,4
            call getrange(i,i1m,i1p)
            class1=0
            if(im.ge.i1m.and.im.le.i1p)then
              class1=i
              goto 102
            endif
101          continue
          write(*,*)'(E) fppfit: no x-section parametrized for ',im
          stop
102          continue
c get class of 2 particle
          do 103 i=0,4
            call getrange(i,i2m,i2p)
            class2=0
            if(jm.ge.i2m.and.jm.le.i2p)then
              class2=i
              goto 104
            endif
103          continue
          write(*,*)'(E) fppfit: no x-section parametrized for ',jm
          stop
104          continue
c get iline corrsponding to pp->i3,i4
          do 105 i=1,nfit
            if((rr(1,i).eq.class1.and.rr(2,i).eq.class2).or.
     v       (rr(1,i).eq.class2.and.rr(2,i).eq.class1))then
              io=i
              goto 106
            endif
105          continue
c maybe we have an R*R* reaction?
          if(im.ge.minres.and.im.le.maxres
     &     .and.jm.ge.minres.and.jm.le.maxres)then
                  io=7
                  goto 106
               endif
          write(*,*)'(E) fppfit: no iline found for particles',im,jm
          stop
106          continue
        else
          io=iio
        endif

c sort particles acc. to itypes
        if(io.ne.5.and.im.gt.jm)then
          im=i4
          jm=i3
        endif
        if(io.eq.5.and.im.lt.jm)then
          im=i4
          jm=i3
        endif

c consistency checks
        call getrange(rr(1,io),i1m,i1p)
        if(im.lt.i1m.or.im.gt.i1p)then
          write(*,*)'(E) fppfit: wrong iline for outgoing particles',
     &      io,im,jm
          stop
        endif
        call getrange(rr(2,io),i2m,i2p)
        if(jm.lt.i2m.or.jm.gt.i2p)then
          write(*,*)'(E) fppfit: wrong iline for outgoing particles',
     &      io,im,jm
          stop
        endif


      goto(1,2,2,2,2,2,4) io
      write(6,*) '****(E) wrong x-section ID in fppfit *****',io
      stop

c pp ->ND
 1    continue
        if(wtabflg.ge.3.and.CTOption(9).eq.0)then
c table lookup
           iim=0
           if(im.eq.minnuc)then
             iim=1
           elseif(im.eq.mindel)then
             iim=2
           else
             write(*,*)'(E) fppfit: First particle should be N or Delta'
             stop
           endif
           fppfit=max(0d0,splint(tabxnd,frrtaby(1,1,iim,jm),
     &       frrtaby(1,2,iim,jm),widnsp,e))
      else
c calculate cross section
         fppfit=pmean(e,im,isoit(im),jm,isoit(jm),id1,id2,id3,id4,1)
     &       /(pmean(e,1,1,1,1,id1,id2,id3,id4,1)
     &       *e*e)*ar(io)
     &       *(massit(jm)**2*widit(jm)**2/((e**2-massit(jm)**2)**2
     &       +widit(jm)**2*massit(jm)**2))
     &       *dble((jit(im)+1)*(jit(jm)+1))
        endif

      return

c pp->NN* pp->ND* pp->DD pp->DN* pp->DD*
 2        continue
        if(wtabflg.ge.3.and.CTOption(9).eq.0)then
c table lookup
          iim=0
           if(im.eq.minnuc)then
             iim=1
           elseif(im.eq.mindel)then
             iim=2
           else
             write(*,*)'(E) fppfit: First particle should be N or Delta'
             stop
           endif
           fppfit=max(0d0,splint(tabxnd,frrtaby(1,1,iim,jm),
     .           frrtaby(1,2,iim,jm),widnsp,e))
      else
c calculate cross section
         tmp=pmean(e,im,isoit(im),jm,isoit(jm),id1,id2,id3,id4,1)
     &         /(pmean(e,1,1,1,1,id1,id2,id3,id4,1)
     &         *e*e)*ar(io)
     &         *dble((jit(im)+1)*(jit(jm)+1))
         if(im.ne.jm) then
            tmp=tmp
     &         /((massit(jm)-massit(im))**2*(massit(jm)+massit(im))**2)
         endif
         fppfit=tmp
      endif
      return

c pp->B*B*
 4    continue
c sofar set to zero
      fppfit=0d0
      return

      end



C####C##1#########2#########3#########4#########5#########6#########7##
      subroutine crossx(iio,e,ii1,iiz1,mm1,ii2,iiz2,mm2,sig)
c
cinput iio  : pointer to cross section
cinput e    : $\sqrt{s}$ of collision
cinput ii1  : ID of particle 1
cinput iiz1 : $2\cdot I_3$ of particle 1
cinput mm1  : mass of particle 1
cinput ii2  : ID of particle 2
cinput iiz2 : $2\cdot I_3$ of particle 2
cinput mm2  : mass of particle 2
coutput sig : cross section
c
c This routine returns cross sections which are accessed via tags (pointers)
c {\tt iio}. The cross sections can either be partial ones for specific
c exit channels or total cross sections for the two incoming hadrons
c
C####C##1#########2#########3#########4#########5#########6#########7##
      implicit none

      include 'comres.f'
      include 'options.f'
      include 'newpart.f'

      integer io,iio,i1,i2,iz1,iz2
      integer i,im,ip,j,jm,jp,ii1,iiz1,ii2,iiz2
c      integer itmp1,itmp2,itag
      real*8 sig,e,e0,sum,m1,m2,m3,gam,dum,s1,s2,s3
      real*8 c1,c2,mm1,mm2,aaqm,dbfact,cgkcor,qflg,factor,ggam
      real*8 sighelp,meltpoint2,pfrome,pmelt,pmelt2,plab
      logical kplflag
c...functions
      real*8 clebsch,siglookup,bcms,massit,widit,mminit,xmelt
      real*8 fppfit,fwidth,fbwnorm,ppiso
c...additional functions and variables for ppbar&pp scattering:
      real*8 bbphi,nnphi,sighera,sapptot,
     &       sappela,sappann,sappdiff,dgcgkfct
      integer isoit,iq1(2),iq2(3),ifdiq2,kq1,kq2

c
      real*8 sigheramb,sigmb,meltpoint

      logical excl

      real*8 ef,sigf,ee
      integer jf,if

c...parameters for resonance production in NN
      integer nfit
      parameter (nfit=7)
c      real*8 er(nfit)
      integer rr(2,nfit),in(nfit)
      real*8 x,x0,xa,xb,xc,bf
      real*8 f

      save kplflag

      f(x,x0,xa,xb,xc)=max(0d0,2.*xb*xa*(x-x0)/(xb**2+(x-x0)**2))*
     *     ((x0+xb)/(x))**xc


c io =         1      2     3     4       5      6      7
c              nd    nn*    nd*   dd      dn*    dd*    B*B*
c      data er/ 4d-2, 0d0,   0d0,  0d0,    0d0,   0d0,   0d0/
      data rr/ 4,0,   4,1,   4,2,  0,0, 0,1,  0,2,    3,3 /
c rr tells of which particle class the out particles are
      data in/ 1,     1,     1,    1,  -1,    1,      1  /
c in tells which out particle has the higher itype


      excl=.false. !flag for exclusive cross section(detailed balance)

      kplflag=.false. ! flag for k+ p total x-section

      goto 107

      entry crossf(if,ef,jf,sigf)
      io=if
      ee=ef
      j=jf
      call getrange(rr(1,io),i,ip) ! call getrange(rr(2,io),jm,jp)

      sigf=fppfit(io,ee,i,j)

      return

      entry crossz(iio,e,ii1,iiz1,mm1,ii2,iiz2,mm2,sig)

      excl=.true.  !flag for exclusive cross section(detailed balance)

 107  continue

c...some settings for all channels enter here
      sig=0d0
      io=iabs(iio)
      aaqm=CTParam(3)      !exponent for bcm scaled AQM-cross sections

      call setizm(ii1,iiz1,mm1,ii2,iiz2,mm2,
     ,i1,iz1,m1,i2,iz2,m2)

c in case of a MB-reaction, sort particles: meson must be second (i2)
      if(iabs(i1).ge.minmes.and.iabs(i1).le.maxmes.and.
     &   iabs(i2).ge.minbar.and.iabs(i2).le.maxbar)then
        call swpizm(i1,iz1,m1,i2,iz2,m2)
      endif


      goto(1,1,1,1,1,1,1,8,9,10,
     ,       11,12,13,14,15,16,17,18,19,9,21,22,23,24,
     ,       25,26,27,28,29,30,31,32,33,34,35,36,37,38)io


      write(6,*)'cross[x,z]: ',
     ,' unknown channel requested, wrong io:',io
      stop


 1    continue
c...pp-> ND pp->NN* pp->ND* pp->DD pp->DN* pp->DD*
      call setizm(iabs(ii1),iiz1,mm1,iabs(ii2),iiz2,mm2,
     ,i1,iz1,m1,i2,iz2,m2)
c     fist outgoing particle has lower ID
      if(i1.gt.i2)call swpizm(i1,iz1,m1,i2,iz2,m2)

c     determine range for first outgoing particle
      call getrange(rr(1,io),im,ip)


      if(excl)then
cut down the loop to one resonance for inverse exclusive processes
        jm=in(io)*max0(in(io)*i1,in(io)*i2)
        jp=jm
      else
        call getrange(rr(2,io),jm,jp)

      end if
      sum=0.

      do i=im,ip
c...loop over i3 (1st outgoing particle)
        do j=jm,jp
c...loop over i4 (2nd outgoing particle)
c   sum over x-section

c explicit isospin dependence
           cgkcor=ppiso(iio,i1,iz1,i2,iz2,i,j)

           sum=sum+fppfit(io,e,i,j)*cgkcor
        end do
      end do
      sig=sum
c
c     detailed balance
      if(iio.lt.0.and.sig.gt.1.d-12) then


            call detbal(e,i1,i2,iz1,iz2,max(mminit(i1),m1),
     &           max(mminit(i2),m2),1,1,dbfact)
            sig=sig*dbfact
      endif

c     check energy conservation (cut-off is one 1 MeV)
      if(max(m1,mminit(i1))+max(m2,mminit(i2))+1.d-3.gt.e) sig=0.d0

      return


 8    continue
c...ND->DD  iso-spin summed value
      call setizm(iabs(ii1),iiz1,mm1,iabs(ii2),iiz2,mm2,
     ,i1,iz1,m1,i2,iz2,m2)

      e0=2.*massit(mindel)-widit(mindel)

      c1=clebsch(isoit(i1),isoit(i2),iz1,iz2,1)
      c2=clebsch(isoit(i1),isoit(i2),iz1,iz2,2)

ctp060926      if(io.lt.0)write(6,*)'crossz(DD->ND):c1,c2=',c1,c2,
ctp060926     ,           isoit(i1),isoit(i2),iz1,iz2,'itypes:',i1,i2
      sig=f(e,e0,12.0d0,0.02d0,2d0)*(0.66667*c1+4d0/dsqrt(20d0)*c2)

      if(iio.lt.0.and.sig.gt.1.d-12) then
            call detbal(e,i1,i2,iz1,iz2,max(mminit(i1),m1),
     &           max(mminit(i2),m2),minnuc,mindel,dbfact)
            sig=sig*dbfact
      endif
      return

 9    continue
      write(6,*)'crossx: channel no.',io,'not implemented.'
      stop

 10   continue
c...MB->B'
      call anndec(0,m1,i1,iz1,m2,i2,iz2,e,sig,gam)
      return

 11   continue
c...MM->M'
      call anndec(0,m1,i1,iz1,m2,i2,iz2,e,sig,gam)
      return

 12   continue
c...??->X => additive quark model total cross section
      call aqm(i1,i2,sig,dum)
      if(sig.gt.1d3)then
        write(6,*)'sig=',sig
      end if

c     check energy conservation (cut-off is one 1 MeV)
      if(max(m1,mminit(i1))+max(m2,mminit(i2))+1.d-3.gt.e) sig=0.d0

      return

 13   continue
c...??->X => additive quark model elastic cross section
      call aqm(i1,i2,dum,sig)

c     check energy conservation (cut-off is one 1 MeV)
      if(max(m1,mminit(i1))+max(m2,mminit(i2))+1.d-3.gt.e) sig=0.d0

      return

 14   continue
c...??->X => additive quark model inelastic cross section
      call aqm(i1,i2,s1,s2)
      s3=s1-s2
      s1=s3
c DD->DD is possible only for e<4gev, then DD->Strings (15) is used
      e0=mminit(i1)+mminit(i2)+2d0*CTParam(4)
      bf=2d0*dabs(CTParam(4)-CTParam(2))
      s2=f(e,e0,s3,bf,1d0)
      sig=xmelt(e,s1,s2,max(mminit(i1),m1)+max(mminit(i2),m2),
     @     e0+2d0*CTparam(2))

c     check energy conservation (cut-off is one 1 MeV)
      if(max(m1,mminit(i1))+max(m2,mminit(i2))+1.d-3.gt.e) sig=0.d0

      return

 15   continue
c...??->X => additive quark model string cross section
      if(CTOption(12).ne.0)return
      if(CTOption(12).eq.0)then
        call aqm(i1,i2,sig,dum)
        sig=sig-dum
c dum is the inelastic xsec
c string threshold is at 3.2 GeV by default
        if(e.gt.2d0*(1.08d0+ctparam(2)))then
          sig=sig*bcms(e,1.08d0+ctparam(2),1.08d0+ctparam(2))**(aaqm)
        else
          sig=0.d0
        end if
      end if

c     check energy conservation (cut-off is one 1 MeV)
      if(max(m1,mminit(i1))+max(m2,mminit(i2))+1.d-3.gt.e) sig=0.d0

      return

 16   continue
c...pn-total...low energy tables & high energy fit
      sig=xmelt(e,siglookup(1,e),sighera(.938d0,.938d0,e,3),3.d0,5.d0)
      return

 17   continue
c...pn-elastic...low energy tables & high energy fit (then pn=pp)
      sig=xmelt(e,siglookup(2,e),sighera(.938d0,.938d0,e,2),3.d0,5.d0)
      return

 18   continue
c...pp-total...low energy tables & high energy fit
      sig=xmelt(e,siglookup(3,e),sighera(.938d0,.938d0,e,1),3.d0,5.d0)
      return

 19   continue
c...pp-elastic...low energy tables & high energy fit
      sig=xmelt(e,siglookup(4,e),sighera(.938d0,.938d0,e,2),3.d0,5.d0)
      return

ctp060202 20   continue
      write(6,*)'crossx: cross section for decay(io=20) requested'
      stop
      return

 21   continue
c...bbar total...from ppbar via AQM & phase-space
      call aqm(1,-1,nnphi,dum)
      call aqm(i1,i2,bbphi,dum)
      sig=sapptot(e,m1,m2)*bbphi/nnphi
      return


 22   continue
c...bbar elastic...from ppbar via AQM & phase-space
      call aqm(1,-1,nnphi,dum)
      call aqm(i1,i2,bbphi,dum)
      sig=sappela(e,m1,m2)*bbphi/nnphi
      return

 23   continue
c...bbar annihilation from ppbar via AQM & phase-space
      call aqm(1,-1,nnphi,dum)
      call aqm(i1,i2,bbphi,dum)
      sig=sappann(e,m1,m2)*bbphi/nnphi
      return

 24   continue
c...bbar diffractive...from ppbar via AQM & phase-space
      call aqm(1,-1,nnphi,dum)
      call aqm(i1,i2,bbphi,dum)
      sig=sappdiff(e,m1,m2)*bbphi/nnphi
      return

 25   continue
c...pi^+ +p
      if ((i2.eq.minmes+1).and.(iz2.eq.2)
     &   .and.(i1.eq.1).and.(iz1.eq.1)) then
         sigheramb=sighera(m2,m1,e,7)
         meltpoint=2.00d0
c...pi^- +p
      elseif ((i2.eq.minmes+1).and.(iz2.eq.-2)
     &   .and.(i1.eq.1).and.(iz1.eq.1)) then
         sigheramb=sighera(m2,m1,e,9)
         meltpoint=2.18d0
c...K^+ +p
      elseif ((i2.eq.minmes+6).and.(iz2.eq.1)
     &   .and.(i1.eq.1).and.(iz1.eq.1)) then
         sigheramb=sighera(m2,m1,e ,11)
         meltpoint=1.84d0
C        here second meltpoint to enable linear interpolation between
C        meltpoint2 and meltpoint (in plab)
         meltpoint2=1.7d0
         kplflag=.true.
c...K^- +p
      elseif ((i2.eq.-(minmes+6)).and.(iz2.eq.-1)
     &   .and.(i1.eq.1).and.(iz1.eq.1)) then
         sigheramb=sighera(m2,m1,e,14)
         meltpoint=2.12d0
c...K^+ +n
      elseif ((i2.eq.minmes+6).and.(iz2.eq.1)
     &   .and.(i1.eq.1).and.(iz1.eq.-1)) then
         sigheramb=sighera(m2,m1,e,13)
         meltpoint=1.75d0
c...K^- +n
      elseif ((i2.eq.-(minmes+6)).and.(iz2.eq.-1)
     &   .and.(i1.eq.1).and.(iz1.eq.-1)) then
         sigheramb=sighera(m2,m1,e,16)
         meltpoint=1.6d0
c...gamma + p
      elseif ((i2.eq.minmes).and.(iz2.eq.0)
     &   .and.iabs(i1).le.maxbar) then
         sigheramb=sighera(m2,m1,e,6)
         meltpoint=1.75d0
      else
c...MB total (->B*/->Strings/el.)
        call anndec(0,m1,i1,iz1,m2,i2,iz2,e,sigmb,gam)
        if (sigmb.gt.1d-8) then
c meltpoint for resonant meson absorption moved to higher energies
          meltpoint=max(1.7d0, m1+m2+.2d0)
        else
          meltpoint=1.7d0
        endif
        call aqm(pimeson,nucleon,dum,nnphi)
        call aqm(i1,i2,dum,bbphi)
        sigheramb=sighera(m2,m1,e,7)*bbphi/nnphi
C
      endif
C
C
c...breit-wigners...
      call anndec(0,m1,i1,iz1,m2,i2,iz2,e,sigmb,gam)
c now melt the low energy resonance x-sec and the high energy hera-fits:
      if (e.lt.meltpoint) then
C
c cross section for Danielewicz forward delay
c     here the DP-cross section has the same form as the normal one
C
         if(CTOption(34).eq.2) then
            sig=sigmb+CTParam(58)*sigmb
         elseif(CTOption(34).eq.3) then
            m3=e
            call anndec(0,m1,i1,iz1,m2,i2,iz2,e,sigmb,gam)
            if(sig.gt.1d-5) sig=sigmb+CTParam(58)
         elseif(CTOption(34).eq.4) then
            m3=e
            call anndec(0,m1,i1,iz1,m2,i2,iz2,e,sigmb,gam)
            if(sigmb.gt.1d-5)then
               ggam=fwidth(itypnew(1),i3new(1),m3)
               sig=sigmb+CTParam(58)*
     &              sig*dsqrt(2.d0/(ggam*3.1415d0*fbwnorm(m3,
     &              itypnew(1),i3new(1))))
            endif
         else
C
C actually default for e.lt.meltpoint
C    take result from anndec as sigma
C
            sig=sigmb
C
         endif
C
c-- nonres adds s-channel strings for strange meson baryon reactions
c-- resonances are not sufficient to fit data - sig gets modified in nonres !!!
C-- so that k- p cross section describes data reasonable.
C-- nonres call only below meltpoint
C
      call nonres(e,i1,iz1,i2,iz2,sig)
C
      else
C
C        for e.ge.meltpoint
C
         sig=sigheramb
C
      endif
C


c.. is annihilation possible ? (quark content?)
      call ityp2id(i2,iz2,iq1(1),iq1(2))
      call ityp2id(i1,iz1,ifdiq2,iq2(3))
      if(abs(i1).ge.minmes) then
        iq2(1)=ifdiq2
        iq2(2)=iq2(3)
        iq2(3)=0
      else
        iq2(1)=mod(ifdiq2/100,10)
        iq2(2)=int(ifdiq2/1000)
      endif
      qflg=0.d0
      do 7312 kq1=1,2
       do 7412 kq2=1,3
c.. annihilation ?
        if(iq1(kq1)+iq2(kq2).eq.0) qflg=1.d0
 7412   continue
 7312  continue
c..--> elastic X-section gets minimum val. of 12.5 mb (K+ P)
c..    if annihilation is not possible
c..      (exclude gamma baryon elastic)
      if (i2.ne.minmes) then
          sig=max(sig,12.5d0*(1d0-qflg))

C       for k+ p linear interpolation between meltpoint2 and meltpoint
          if (kplflag.and.
     &        qflg.eq.0.d0.and.e.lt.meltpoint
     &       .and.e.gt.meltpoint2) then
             sighelp=sighera(m2,m1,meltpoint,11)
             pfrome=plab(m2,m1,e)
             pmelt2=plab(m2,m1,meltpoint2)
             pmelt=plab(m2,m1,meltpoint)
             sig=xmelt(pfrome,1.25d1,sighelp,pmelt2,pmelt)
          endif
      endif

c     check energy conservation (cut-off is one 1 MeV)
      if(max(m1,mminit(i1))+max(m2,mminit(i2))+1.d-3.gt.e) sig=0.d0


C
C  end of MB total
C
      return

 26   continue
      if(i2.eq.minmes)then
c..no elastic gamma B scattering:
       sig=0.d0
       return
      endif
c...MB->MB elastic (pi+ p scaled with aqm)
      if(mminit(i1)+mminit(i2).gt.e) then
c     no collision if sqrts is insufficient to put particles on-shell
         sig=0.d0
      else
         call aqm(pimeson,nucleon,dum,nnphi)
         call aqm(i1,i2,dum,bbphi)
         sig=sighera(m2,m1,e,8)*bbphi/nnphi
      endif
      return

 27   continue
c..MB-> 1 String (s-channel) (pi+ p scaled with aqm)
      if(CTOption(12).ne.0)return
      call aqm(pimeson,nucleon,nnphi,dum)
      call aqm(i1,i2,bbphi,dum)
      sig=sighera(m2,m1,e,7)*bbphi/nnphi
c-- minus elastic x-sect. (scaled with el. aqm)
      call aqm(pimeson,nucleon,dum,nnphi)
      call aqm(i1,i2,dum,bbphi)
      sig=sig-sighera(m2,m1,e,8)*bbphi/nnphi
c.. is s-channel possible ? (quark content?)
      call ityp2id(i2,iz2,iq1(1),iq1(2))
      call ityp2id(i1,iz1,ifdiq2,iq2(3))
      if(abs(i1).ge.minmes) then
        iq2(1)=ifdiq2
        iq2(2)=iq2(3)
        iq2(3)=0
      else
        iq2(1)=mod(ifdiq2/100,10)
        iq2(2)=int(ifdiq2/1000)
      endif
      qflg=0.d0
      do 312 kq1=1,2
       do 412 kq2=1,3
c..two of the quarks must be able to annihilate:
        if(iq1(kq1)+iq2(kq2).eq.0) qflg=1.d0
 412   continue
 312  continue
      sig=xmelt(e,sig,0.d0,3d0,6.d0)*qflg
      if(e.le.1.7d0)sig=0.d0
      call nonres(e,i1,iz1,i2,iz2,sig)

      return

 28   continue
c..MB-> 2 Strings (t-channel) (pi+ p scaled with aqm)
      if(CTOption(12).ne.0)return
      call aqm(pimeson,nucleon,nnphi,dum)
      call aqm(i1,i2,bbphi,dum)
      sig=sighera(m2,m1,e,7)*bbphi/nnphi
c-- minus elastic x-sect. (scaled with el. aqm)
      call aqm(pimeson,nucleon,dum,nnphi)
      call aqm(i1,i2,dum,bbphi)
      sig=sig-sighera(m2,m1,e,8)*bbphi/nnphi
c-- melting
      sig=xmelt(e,0d0,sig,3d0,6.d0)
      return

 29   continue
c...??->X => additive quark model inelastic cross section
c specially reduced to avoid double-counting in case of DN and DD
c    cross sections
      call aqm(i1,i2,sig,dum)
      sig=sig-dum
      sig=sig*bcms(e,max(mminit(i1),m1),max(mminit(i2),m2))**aaqm
      sig=xmelt(e,sig,0d0,2d0*(1.08d0+ctparam(2)),5.0d0)
c here comes the reduction
      sig=xmelt(e,0d0,sig,2.75d0,3.4d0)
      return

 30   continue
c parameterized detailed balance cross sections for ND->NN
      factor=dgcgkfct(i1,i2,iz1,iz2,nucleon,nucleon)
      if(factor.le.1.d-8) then
         sig=0.d0
         return
      endif
      cgkcor=ppiso(-1,i1,iz1,i2,iz2,nucleon,nucleon)

         sig=factor*cgkcor*(1.3d8*(e**(-17.5d0))+3.6d4*(e**(-7d0)))

c         write(6,*)factor,cgkcor,sig

      return

 31   continue
c parameterized detailed balance cross sections for DD->DN
      e0=massit(i1)-0.5*widit(i1)+massit(i2)-0.5*widit(i2)

      factor=dgcgkfct(i1,i2,iz1,iz2,nucleon,mindel)
      if(factor.le.1.d-8.or.e.lt.e0)then
         sig=0.d0
      else
         c1=clebsch(isoit(i1),isoit(i2),iz1,iz2,1)
         c2=clebsch(isoit(i1),isoit(i2),iz1,iz2,2)

c param
         sig=
     &        (2.5d56*exp(-(50.0d0*e))
     &        +4.9d14*exp(-(12.d0*e))
     &        +1.1d6*exp(-(4.50d0*e)))
     &        *factor
c the following factors are from Heinz Sorge's Habilitation
     &        *(0.66667*c1+4d0/dsqrt(20d0)*c2)

      endif
      return

 32   continue
C parameterized detailed balance cross sections for DD->NN
      factor=dgcgkfct(i1,i2,iz1,iz2,nucleon,nucleon)
      if(factor.le.1.d-8.or.e.lt.2.15d0) then
         sig=0.d0
         return
      endif
      cgkcor=ppiso(-4,i1,iz1,i2,iz2,nucleon,nucleon)
c param
      sig=(7.27d0*(e-2.14d0)**(-1.2176d0)+0.05d0*(e-2.14d0)**(-3.257d0))
     &     *factor*cgkcor

      return

 33   continue
c...??->X => additive quark model resonance cross section
       call aqm(i1,i2,sig,dum)
       sig=sig-dum  ! dum is the elastic xsec
      if(e.gt.mminit(i1)+mminit(i2))then
c x_string+x_resonances
          sig=sig*bcms(e,mminit(i1),mminit(i2))**(aaqm*3)
c get x_string
          sig=xmelt(e,sig,0d0,mminit(i1)+mminit(i2)+1.5d0,
     @         mminit(i1)+mminit(i2)+3d0)

      end if

      return

 34   continue
c...??->X => additive quark model string cross section
      if(CTOption(12).ne.0)return
        call aqm(i1,i2,sig,dum)
        sig=sig-dum !  dum is the elastic xsec
       if(e.gt.mminit(i1)+mminit(i2))then
c x_string+x_resonances
          sig=sig*bcms(e,mminit(i1),mminit(i2))**(aaqm*3)
c get x_resonances
          sig=xmelt(e,0d0,sig,mminit(i1)+mminit(i2)+1.5d0,
     @         mminit(i1)+mminit(i2)+3d0)
      end if
      return

 35   continue
      e0=mminit(i1)+mminit(i2)+CTParam(4)*2d0
      call aqm(i1,i2,sig,dum)
      sig=f(e,e0,sig-dum,1d0,1d0)
      return

 36   continue
c cross section for Danielewicz forward delay
c...MB->B'
      if(CTOption(34).eq.2) then
         call anndec(0,m1,i1,iz1,m2,i2,iz2,e,sig,gam)
         sig=CTParam(58)*sig
      elseif(CTOption(34).eq.3) then
         m3=e
         call anndec(0,m1,i1,iz1,m2,i2,iz2,e,sig,gam)
         if(sig.gt.1d-5) sig=CTParam(58)
      elseif(CTOption(34).eq.4) then
         m3=e
         call anndec(0,m1,i1,iz1,m2,i2,iz2,e,sig,gam)
         if(sig.gt.1d-5)then
            ggam=fwidth(itypnew(1),i3new(1),m3)
            sig=CTParam(58)*
     &           sig*dsqrt(2.d0/(ggam*3.1415d0*
     &           fbwnorm(m3,itypnew(1),i3new(1))))
         endif
      else
         sig=0.d0
      endif
      return

 37   continue
c cross section for Danielewicz forward delay
c...MM->M'
      if(CTOption(34).eq.2) then
         call anndec(0,m1,i1,iz1,m2,i2,iz2,e,sig,gam)
         sig=CTParam(58)*sig
       elseif(CTOption(34).eq.3) then
          m3=e
          call anndec(0,m1,i1,iz1,m2,i2,iz2,e,sig,gam)
          if(sig.gt.1d-5) sig=CTParam(58)
       elseif(CTOption(34).eq.4) then
          m3=e
          call anndec(0,m1,i1,iz1,m2,i2,iz2,e,sig,gam)
          if(sig.gt.1d-5)then
             ggam=fwidth(itypnew(1),i3new(1),m3)
             sig=CTParam(58)*
     &            sig*dsqrt(2.d0/(ggam*3.1415d0*
     &           fbwnorm(m3,itypnew(1),i3new(1))))
          endif
      else
         sig=0.d0
      endif
      return
 38   continue
c elastic meson-meson cross section
      sig=5.d0

c     check energy conservation (cut-off is one 1 MeV)
      if(max(m1,mminit(i1))+max(m2,mminit(i2))+1.d-3.gt.e) sig=0.d0

      return

ctp060202 99   continue
c...single diffr. pp
      sig=.68*(1.+36/e**2)*log(0.6+0.1*e**2)
      return

      end

C####C##1#########2#########3#########4#########5#########6#########7##
      subroutine nonres(e,ii1,iiz1,ii2,iiz2,sig)
c
cinput e    : $\sqrt{s}$ of collision
cinput ii1  : ID of particle 1
cinput iiz1 : $2\cdot I_3$ of particle 1
cinput ii2  : ID of particle 2
cinput iiz2 : $2\cdot I_3$ of particle 2
coutput sig : cross section
c
c {\tt nonres} adds s-channel strings for strange meson baryon reactions,
c since
c resonances are not sufficient to fit data -
c {\tt sig} gets modified in {\tt nonres},
C so that $k^- p$ cross section describes data reasonable.
C {\tt nonres} should be only called below {\tt meltpoint}
c
C####C##1#########2#########3#########4#########5#########6#########7##
      implicit none
      logical b
      real*8 m1,m2,e,sig
      integer strit,i1,i2,iz1,iz2,ii1,ii2,iiz1,iiz2

      b=.false.
      call setizm(ii1,iiz1,0d0,ii2,iiz2,0d0,
     @      i1,iz1,m1,i2,iz2,m2)
      if(iabs(i1).gt.iabs(i2))
     @     call swpizm(i1,iz1,m1,i2,iz2,m2)

C    check wether sigma should be modified
C        1) strange meson required
C        2) nonstrange baryon required
C        3) combination of strange meson + baryon (k- p)
C         or antistrange meson + antibaryon (k+ pbar) required
C
      if(strit(i2).ne.0.and.i1*i2.lt.0.and.strit(i1).eq.0)b=.true.
C   if this condition is fullfilled modify sigma
C   add exponential underground and corrections
C   (two gaussians,constant,low energy cut off)
C       to describe experimental k- p cross section
C    four  cases for different corrections
C
      if ((b).and.e.gt.1.433.and.e.lt.1.4738188) then
        sig=sig+120.
      elseif ((b).and.e.ge.1.4738188.and.e.lt.1.485215) then
        sig=sig-1.296457765d7*(e-1.433)**4+
     &        2.160975431d4*(e-1.433)**2+120.
      elseif ((b).and.e.ge.1.485215.and.e.lt.1.977) then
        sig=sig+1.07769d+06*exp(-(6.44463d0*e))-
     &         10.*exp(-(((e-1.644)**2)/0.004))+
     &         10.*exp(-(((e-1.977)**2)/0.004))
      elseif ((b).and.e.ge.1.977.and.e.lt.2.12) then
C keep maximum value of gaussian above e=1.977 GeV
C                         (e=2.12 GeV is meltpoint)
        sig=sig+1.07769d+06*exp(-(6.44463d0*e))+10.
      endif

C
      return
      end

C####C##1#########2#########3#########4#########5#########6#########7##
      real*8 function sappann(sroot,m1,m2)
c
cinput sroot : $\sqrt{s}$ of collision
cinput m1    : mass of 1st (anti-)baryon
cinput m2    : mass of 2nd (anti-)baryon
c
c  OUTPUT: {\ttsappann} = annihilation cross section for $\bar N N$
c
c  Taken from: P. Koch, C.B. Dover, Phys. Rev. {\bf C40} (1989) 145
c
C####C##1#########2#########3#########4#########5#########6#########7##

      implicit none

      include 'options.f'

      real*8 sroot,m1,m2,s0,a,b,sig0,s,srootnn,snn
      parameter(a=0.05d0,b=0.6d0,sig0=120.d0,s0=3.52)

      if (CTOption(38).eq.1) then
c evaluate the parametrization at the same relative momentum as in nbar-n
         srootnn=snn(sroot,m1,m2)
         s=srootnn**2
      else
c evaluate parametrization now at the same sqrts
         s=sroot**2
      endif
      sappann=sig0*(s0/s)*(a**2*s0/((s-s0)**2+a**2*s0)+b)

      return
      end

C####C##1#########2#########3#########4#########5#########6#########7##
      real*8 function sapptot(sroot,m1,m2)
c
cinput sroot : $\sqrt{s}$ of collision
cinput m1    : mass of 1st (anti-)baryon
cinput m2    : mass of 2nd (anti-)baryon
c
c  OUTPUT: {\ttsapptot} = total cross section for $\bar N N$
c  high energy paramametrization and data taken from PRD 50 (1994)\\
c  $p_{lab} > 5 $GeV:  CERN/HERA parametrization\\
c  0.3 GeV $< p_{lab} <$ 5 GeV: polynomial fit to the data (by C.S.)\\
c  $p_{lab} <0.3$ GeV: another fit, only constrained by sigma annihilation
c
C####C##1#########2#########3#########4#########5#########6#########7##

      implicit none

      include 'options.f'

      real*8 p,sroot,plab,sighera,m1,m2,srootnn,snn

      if (CTOption(38).eq.1) then
c evaluate the parametrization at the same relative momentum as in nbar-n
         srootnn=snn(sroot,m1,m2)
      else
c evaluate parametrization now at the same sqrts
         srootnn=sroot
      endif

      p=plab(0.938d0,0.938d0,srootnn)

      if(p.ge.5.d0)then
         sapptot=sighera(0.938d0,0.938d0,srootnn,4)
        return
      else if(p.ge.0.3d0)then
        sapptot=75.0146d0+43.1276d0/p+2.58298d0/p**2-3.90783d0*p
        return
      else
        sapptot=271.6d0*exp(-(1.1d0*p**2))
        return
      endif

      end


C####C##1#########2#########3#########4#########5#########6#########7##
      real*8 function sappela(sroot,m1,m2)
c
cinput sroot : $\sqrt{s}$ of collision
cinput m1    : mass of 1st (anti-)baryon
cinput m2    : mass of 2nd (anti-)baryon
c
c  OUTPUT: {\tt sappela} = elastic cross section for $\bar N N$
c  high energy paramametrization and data taken from PRD 50 (1994)\\
c  $p_{lab} > 5 $GeV:  CERN/HERA parametrization\\
c  0.3 GeV $< p_{lab} <$ 5 GeV: polynomial fit to the data (by C.S.)\\
c  $p_{lab} <0.3$ GeV: no data, set constant
c
C####C##1#########2#########3#########4#########5#########6#########7##

      implicit none

      include 'options.f'

      real*8 p,sroot,plab,sighera,m1,m2,srootnn,snn

      if (CTOption(38).eq.1) then
c evaluate the parametrization at the same relative momentum as in nbar-n
         srootnn=snn(sroot,m1,m2)
      else
c evaluate parametrization now at the same sqrts
         srootnn=sroot
      endif

      p=plab(0.938d0,0.938d0,srootnn)
      if(p.ge.5.d0)then
         sappela=sighera(0.938d0,0.938d0,srootnn,5)
        return
      else if(p.ge.0.3d0)then
        sappela=31.6166d0+18.2842d0/p-1.14896d0/p**2-3.79508d0*p
        return
      else
        sappela=78.6d0
        return
      endif

      end

C####C##1#########2#########3#########4#########5#########6#########7##
      real*8 function sappdiff(sroot,m1,m2)
c
cinput sroot : $\sqrt{s}$ of collision
cinput m1    : mass of 1st (anti-)baryon
cinput m2    : mass of 2nd (anti-)baryon
c
c  OUTPUT: {\tt sappdiff} = diffractive cross section for $\bar N N$
c
c  This cross section is totally determined by
c $ \sigma_{diff}=\sigma_{tot}-\sigma_{elast.}-\sigma_{annihil.}$
c
C####C##1#########2#########3#########4#########5#########6#########7##
      implicit none
      real*8 p,sroot,m1,m2,plab,sappann,sapptot,sappela

      p=plab(0.938d0,0.938d0,sroot)
      if(p.le.0.1d0)then
        sappdiff=0.d0
        return
      else
        sappdiff=max(0.d0,sapptot(sroot,m1,m2)
     &           -sappela(sroot,m1,m2)-sappann(sroot,m1,m2))
        return
      endif
      end

C####C##1#########2#########3#########4#########5#########6#########7##
      blockdata herafits
c  cross section parameters for specific collsion type
c  CERN/HERA fits, taken from PRD 50 (1994)
c  see: function 'sighera'
C####C##1#########2#########3#########4#########5#########6#########7##
      implicit none
      integer nfit
c     integer io
      parameter (nfit=16)
      real*8 a(nfit),b(nfit),n(nfit),c(nfit),d(nfit)
     &       ,p1(nfit),p2(nfit)
      common /HERA/ a,b,n,c,d,p1,p2

c io =          1       2       3       4        5        6
c            pp(tot) pp(ela) pn(tot) app(tot) app(ela) gammap(tot)
c               7       8       9      10       11       12
c         pi+p(tot) pi+p(el) pi-p(tot) pi-p(el) k+p(tot) k+p(el)
c              13      14      15      16
c          k+n(tot) k-p(tot) k-p(ela) k-n(tot)
c
c p1 (p2) give the momentum range (lab momentum) of the fit

      data a/  48.0,   11.9,  47.3,   38.4,     10.2,   0.147,
     @         16.4,     0.,  33.0,   1.76,     18.1,     5.0,
     @         18.7,   32.1,   7.3,   25.2/
      data b/    0.,   26.9,    0.,   77.6,     52.7,      0.,
     @         19.3,   11.4,   14.0,  11.2,       0.,     8.1,
     @           0.,     0.,     0.,    0./
      data n/    0.,  -1.21,    0.,  -0.64,    -1.16,      0.,
     @         -.42,    -.4, -1.36,  -0.64,       0.,    -1.8,
     @           0.,     0.,    0.,     0./
      data c/ 0.522,  0.169, 0.513,   0.26,    0.125,   .0022,
     @          .19,   .079,  .456,   .043,      .26,     .16,
     @          .21,    .66,   .29,    .38/
      data d/ -4.51,  -1.85, -4.27,   -1.2,    -1.28,   -.017,
     @           0.,     0., -4.03,     0.,      -1.,    -1.3,
     @         -.89,   -5.6,  -2.4,   -2.9/
      data p1/   3.,     2.,    3.,     5.,       5.,      3.,
     @           3.,     2.,   1.8,     1.8,       2.,      2.,
     @           2.,     1.75,    3.,    1.8/
c   for (k- p) p1=1.75 GeV/c works reasonable
      data p2/2100.,  2100.,  370.,  1.7d6,    1.7d6,    183.,
     @         340.,   200.,  370.,   360.,     310.,    175.,
     @         310.,   310.,  175.,   310./

      end

C####C##1#########2#########3#########4#########5#########6#########7##
      real*8 function sighera(m1,m2,sroot,io)
c
c
cinput sroot : $\sqrt{s}$ of collision
cinput m1    : mass of 1st (anti-)baryon
cinput m2    : mass of 2nd (anti-)baryon
cinput io    : flag for respective cross section
c
c  OUTPUT: {\tt sighera} = cross section for specific collsion type
c
c {\tt io} can have the following values:
c \begin{tabular}{rl}
c               1 & $pp$ total \\
c               2 & $pp$ elastic \\
c               3 & $pn$ total \\
c               4 & $\bar p p$ total\\
c               5 & $\bar p p$ elastic\\
c               6 & $\gamma p$ (tot) \\
c               7 & $\pi^+ p$ (tot) \\
c               8 & $\pi^+ p$ (el) \\
c               9 & $\pi^- p$ (tot)\\
c              10 & $\ pi- p$ (el) \\
c              11 & $k^+ p$ (tot)\\
c              12 & $k^+ p$ (el)\\
c              13 & $k^+ n$ (tot)\\
c              14 & $k^- p$ (tot)\\
c              15 & $k^- p$ (ela)\\
c              16 & $k^- n$ (tot)\\
c \end{tabular}
c
c This subroutine returns CERN/HERA parametrizations for cross sections
c {\tt p1} and {\tt p2} in the {\tt blockdata} routine
c give the momentum range (lab momentum) of the fit.
c The fits have been taken from
c Phys. Rev. {\bf D50} (1994).
c
C####C##1#########2#########3#########4#########5#########6#########7##
      implicit none
      include 'coms.f'
      include 'comres.f'
      integer nfit,io
      parameter (nfit=16)
      real*8 a(nfit),b(nfit),n(nfit),c(nfit),d(nfit),sroot,p,plab
     &       ,p1(nfit),p2(nfit),m1,m2,massit
      common /HERA/ a,b,n,c,d,p1,p2

c      p=plab(m1,m2,sroot)
      p=0.d0
      if(io.ge.1.and.io.le.5)then
       p=plab(massit(minnuc),massit(minnuc),sroot)
      elseif(io.eq.6)then
       p=plab(massit(minmes),massit(minnuc),sroot)
      elseif(io.ge.7.and.io.le.10)then
       p=plab(massit(pimeson),massit(minnuc),sroot)
      elseif(io.ge.11.and.io.le.16)then
       p=plab(massit(itkaon),massit(minnuc),sroot)
      else
       write(*,*)'#make22 error. sighera called with io=',io
       stop
      endif

         if(p.lt.1d-15) then
c if energy conservation is not possible (p.eq.0) then return zero
            sighera=0.d0
            return
         endif
C
      if(p.lt.p1(io))then
        p=p1(io)
C
      elseif (p.gt.p2(io).and.(warn)) then
        write(6,*)'sighera: sroot=',sroot,' high!, io=',io
        write(6,*)'         m1,m2,plab,p2(io)=',m1,m2,p,p2(io)
        write(6,*)'sighera fit used above upper limit (extrapolation)'
      endif
        sighera=a(io)+b(io)*p**n(io)+c(io)*log(p)**2+d(io)*log(p)
      return
      end

C####C##1#########2#########3#########4#########5#########6#########7##
      real*8 function xmelt(e,x1,x2,em,ep)
c
cinput e  :  $\sqrt{s}$ of process
cinput x1 : value at {\tt em}
cinput x2 : value at {\tt ep}
cinput em : lower boundary in $\sqrt{s}$
cinput ep : upper boundary in $\sqrt{s}$
c
c  {\tt xmelt } yields an interpolation beween {\tt x1} and {\tt x2}
c  in the $\sqrt{s}$ range beween {\tt em} and {\tt ep}.
c  For {\tt e < em}  it equals {\tt x2},  if {\tt e > em} it yields
c  a combination of both {\tt x1} and {\tt x2} such that there is a
c  continuous  transition from {\tt x1} to {\tt x2}.
c
c  For parameter b=.false. a linear combination is used,
c  for parameter b=.true. a sin-form is used.
c
C####C##1#########2#########3#########4#########5#########6#########7##
      implicit none
      real*8 e,x1,x2,em,ep,lam,pi
      logical b
      parameter(b=.false.,pi=.31415927d-1)


      lam=max(0d0,min(1d0,(e-min(em,ep))/abs(ep-em)))
      if(b)lam=5d-1*sin((lam-5d-1)*pi)+5d-1

      xmelt=(1d0-lam)*x1+lam*x2

      return
      end

C####C##1#########2#########3#########4#########5#########6#########7##
      real*8 function plab(m1,m2,sroot)
c
cinput sroot : $\sqrt{s}$ of collision
cinput m1    : mass of 1st (anti-)baryon
cinput m2    : mass of 2nd (anti-)baryon
c
c {\tt plab} returns the lab-momentum of particle 1
c
C####C##1#########2#########3#########4#########5#########6#########7##
      implicit none
      real*8 sroot,m1,m2
      if(sroot-m1-m2.lt.0d0)then
        plab=0d0
      else
        plab=sqrt((sroot**2-(m1+m2)**2)*(sroot**2-(m1-m2)**2))/(2*m2)
      end if
      return
      end
C####C##1#########2#########3#########4#########5#########6#########7##
      real*8 function snn(sbb,m1,m2)
c
cinput sbb   : $\sqrt{s}$ of a $BB$ collision
cinput m1    : mass of 1st (anti-)baryon
cinput m2    : mass of 2nd (anti-)baryon
c
c  {\tt snn} returns the  equivalent c.o.m. energy of a $NN$-collision
c                with the same relative momentum
c
C####C##1#########2#########3#########4#########5#########6#########7##
      implicit none
      real*8 sbb,prel,m1,m2
      prel=sqrt((sbb**2-(m1+m2)**2)*(sbb**2-(m1-m2)**2))/sbb
      snn=sqrt(prel**2+3.52d0)
      return
      end
