c-----------------------------------------------------------------------
      subroutine conaa(iret)
c-----------------------------------------------------------------------
c  determines interaction configuration
c-----------------------------------------------------------------------

      include 'epos.inc'
      include 'epos.incems'
      include 'epos.incsem'
      include 'epos.incpar'

      common/geom/rmproj,rmtarg,bmax,bkmx
      common/nucl3/phi,bimp
      common /cncl/xproj(mamx),yproj(mamx),zproj(mamx)
     *            ,xtarg(mamx),ytarg(mamx),ztarg(mamx)
      common/cfacmss/facmss
      double precision yyrmax
      common/scrangle/ phik3(kollmx),thetak3(kollmx)

      call utpri('conaa ',ish,ishini,4)

      iret=0

c     initialisations
c     ---------------

      vel=tanh(ypjtl-yhaha)+tanh(yhaha)

c     determine phi, bimp, coll, iproj, itarg, x/y/zproj, x/y/ztarg
c     ---------------------------------------------------------------

           if(iokoll.eq.1)then

      koll=matarg
      do k=1,koll
        do n=1,4
          coord(n,k)=0.
        enddo
      enddo
      bimp=0
      phi=0
      xproj(1)=0
      yproj(1)=0
      zproj(1)=0
      lproj(1)=koll
      lproj3(1)=0
      do k=1,koll
        bij=bkmx*sqrt(rangen())
        bk(k)=bij
        iproj(k)=1
        itarg(k)=k
        phi=2.*pi*rangen()
        xtarg(k)=bij*cos(phi)
        ytarg(k)=bij*sin(phi)
        ztarg(k)=0
        ltarg(k)=1
        kproj(1,k)=k
        ktarg(k,1)=k
        ltarg3(k)=0
        ktarg3(k,1)=0
        kproj3(k,1)=0
        if(iscreen.ne.0.and.bij.le.bkmxndif)then
          if(zbrmax.le.0..or. bij.lt.zbcutx+zbrmax*rangen())then
            lproj3(1)=lproj3(1)+1
            ltarg3(k)=1
            kproj3(1,lproj3(1))=k
            ktarg3(k,1)=k
          endif
        endif
      enddo

           elseif(maproj.eq.1.and.matarg.eq.1)then

      b1=bminim
      b2=amin1(bkmx,bmaxim)
      if(b1.gt.b2)call utstop('conaa: bmin > bmax&')
      bimp=sqrt(b1*b1+(b2*b2-b1*b1)*rangen())
      koll=1
      do n=1,4
        coord(n,1)=0.
      enddo
      bk(1)=bimp
      iproj(1)=1
      itarg(1)=1
      phi=2.*pi*rangen()
      xproj(1)=bimp*cos(phi)
      yproj(1)=bimp*sin(phi)
      zproj(1)=0
      xtarg(1)=0
      ytarg(1)=0
      ztarg(1)=0
      lproj(1)=1
      ltarg(1)=1
      lproj3(1)=1
      ltarg3(1)=1
      kproj3(1,1)=1
      ktarg3(1,1)=1
      kproj(1,1)=1
      ktarg(1,1)=1

           else

      call conxyz('p',mamx,xproj,yproj,zproj,ypjtl-yhaha)
      call conxyz('t',mamx,xtarg,ytarg,ztarg,yhaha)

      bx=0
      by=0
      if(maproj.gt.0)then
      if(bimevt.lt.0)then
        b1=bminim
        b2=amin1(rmproj+rmtarg,bmaxim)
        if(b1.gt.b2)call utstop('conaa: bmin > bmax&')
        bimp=sqrt(b1**2+(b2**2-b1**2)*rangen())
        if(nbarray.gt.0)bimp=barray(mod(nrevt,nbarray)+1)
        if(jpsi.gt.0)then
          bimp=b1+(b2-b1)*(float(mod(nrevt,12))+rangen())/12.
          bimevt=bimp
        endif
        phi=phimin+rangen()*(phimax-phimin)
      else
        phi=phievt
        bimp=bimevt
      endif
      bx=cos(phi)*bimp
      by=sin(phi)*bimp
      endif
      if(jpsi.lt.0)then !modify b
        bx=xtarg(1)
        by=ytarg(1)
      endif
      if(maproj.eq.0)goto1000
      koll=0
      do i=1,maproj
      lproj(i)=0
      lproj3(i)=0
      enddo
      do j=1,matarg
      ltarg(j)=0
      ltarg3(j)=0
      enddo
      do 12 i=1,maproj
      do 11 j=1,matarg
      if(jpsi.lt.0.and.ztarg(j).le.ztarg(1))goto11
      bij=sqrt((xproj(i)+bx-xtarg(j))**2+(yproj(i)+by-ytarg(j))**2)
      if(ish.ge.7)write(ifch,*)'i_p:',i,' i_t:',j,' b_ij:',bij
      if(bij.gt.bkmx)goto 11

      koll=koll+1
      if(koll.gt.kollmx)call utstop('conaa: kollmx too small&')
      bk(koll)=bij
      bkx(koll)=xproj(i)+bx-xtarg(j)
      bky(koll)=yproj(i)+by-ytarg(j)
      iproj(koll)=i
      itarg(koll)=j
      lproj(i)=lproj(i)+1
      ltarg(j)=ltarg(j)+1
      kproj(i,lproj(i))=koll
      ktarg(j,ltarg(j))=koll
      phik3(koll)=0.
      thetak3(koll)=0.
      if(iscreen.ne.0.and.bij.le.bkmxndif)then
        if(zbrmax.gt.0..and.bij.gt.zbcutx+zbrmax*rangen())goto 11
        lproj3(i)=lproj3(i)+1
        ltarg3(j)=ltarg3(j)+1
        kproj3(i,lproj3(i))=koll
        ktarg3(j,ltarg3(j))=koll
c define angle for anti-shadowing
        if(abs(bky(koll)).gt.1.e-6)then
          if(abs(bkx(koll)).gt.1.e-6)then
            phik3(koll)=atan(bky(koll)/bkx(koll))
          else
            phik3(koll)=sign(0.5*pi,bky(koll))
          endif
        elseif(bkx(koll).lt.0.)then
          phik3(koll)=pi
        endif
        if(bk(koll).gt.0.)then
          thetak3(koll)=atan(bglaubx/bk(koll))
        else
          thetak3(koll)=0.5*pi
        endif
      endif

11    continue
12    continue

      do k=1,koll
        do n=1,4
          coord(n,k)=0.
        enddo
      enddo


          endif

      if(ish.ge.3)write(ifch,*)'koll=',koll
      if(koll.eq.0)goto 1001


c     determine coord
c     ---------------
      do kl=1,koll
      i=iproj(kl)
      j=itarg(kl)
      dist=ztarg(j)-zproj(i)
      coord(1,kl)=(xproj(i)+xtarg(j))*0.5
      coord(2,kl)=(yproj(i)+ytarg(j))*0.5
      coord(3,kl)=(zproj(i)+ztarg(j))*0.5
      coord(4,kl)=dist/vel
      enddo

      if(iscreen.ne.0)call CalcScrPair(bimp)

    !~~~~~redefine energy in case of imposed radial flow~~~~~~~~~~~~~~~~
      yrmaxi=max(0.,yradmx+yradmi*log(1.+engy/sqrt(float(koll))))
      if(yrmaxi.gt.1e-5)then
        yyrmax=dble(yrmaxi)
        fradflii=sngl(1d0/
     &  ((sinh(yyrmax)*yyrmax-cosh(yyrmax)+1d0)/(yyrmax**2/2d0)))
      else
        fradflii=1.
      endif
      if(ish.ge.3)write(ifch,*)'yrmaxi=',yrmaxi


c     exit
c     ----
1000  continue
      if(ish.ge.5)then
      write(ifch,*)'ia,x/y/zproj:'
      do mm=1,maproj
      write(ifch,*)mm,xproj(mm),yproj(mm),zproj(mm)
      enddo
      write(ifch,*)'ia,x/y/ztarg:'
      do mm=1,matarg
      write(ifch,*)mm,xtarg(mm),ytarg(mm),ztarg(mm)
      enddo
      write(ifch,*)'iret',iret
      endif
      call utprix('conaa ',ish,ishini,4)
      return

1001  continue !iret=1 causes redo of whole collision
      iret=1
      if(ish.ge.3)then
      write(ifch,*)
      write(ifch,*)'***** subroutine conaa:'
      write(ifch,*)'***** no nucleon pair found --> no interaction'
      write(ifch,*)
      endif
      goto 1000

      end

c-----------------------------------------------------------------------
      function frdmzcut(b)
c-----------------------------------------------------------------------
c  determines interaction configuration
c-----------------------------------------------------------------------

      include 'epos.inc'
      include 'epos.incpar'

      frdmzcut=zbcutx*exp(-(b*rangen()))
c      frdmzcut=zbcutx*rangen()
c      bdum=b
      return
      end

c-----------------------------------------------------------------------
      subroutine CalcScrPair(b)
c-----------------------------------------------------------------------
c  determines interaction configuration
c-----------------------------------------------------------------------

      include 'epos.inc'
      include 'epos.incems'
      include 'epos.incpar'

      common/scrangle/ phik3(kollmx),thetak3(kollmx)
      logical lascr(kollmx),cont


c     order pairs for splitting as a function of b
c     ---------------
      do i=1,maproj
        if(lproj3(i).gt.1)then
          do l=2,lproj3(i)
            m=l
 100        continue
              kp=kproj3(i,m-1)
              kl=kproj3(i,m)
              if(bk(kl).lt.bk(kp))then
                kproj3(i,m-1)=kl
                kproj3(i,m)=kp
                m=m-1
            if(m.gt.1)goto 100
              endif
          enddo
        endif
      enddo
      do j=1,matarg
        if(ltarg3(j).gt.1)then
          do l=2,ltarg3(j)
            m=l
 200        continue
              kp=ktarg3(j,m-1)
              kl=ktarg3(j,m)
              if(bk(kl).lt.bk(kp))then
                ktarg3(j,m-1)=kl
                ktarg3(j,m)=kp
                m=m-1
            if(m.gt.1)goto 200
              endif
          enddo
        endif
      enddo

      if(koll.gt.1)then
c Define anti-shadowing as a consequence of geometrical screening of
c nucleons at large b by the one at small b

c Projectile
c Check phi not to be in the range of a nucleon with smaller b
        do i=1,maproj
          if(lproj3(i).gt.1)then
            kl=kproj3(i,1)
            lascr(kl)=.true.
            do 300 l=lproj3(i),2,-1
              kl=kproj3(i,l)
              lascr(kl)=.true.
              if(bk(kl).ge.frdmzcut(b))then
                do m=1,l-1
                  km=kproj3(i,m)
                  if(kl.ne.km.and.bk(km).ge.frdmzcut(b)
     &          .and.phik3(kl).ge.phik3(km)-thetak3(km)
     &          .and.phik3(kl).le.phik3(km)+thetak3(km))then
                    lascr(kl)=.false.
                    goto 300
                  endif
                enddo
              endif
 300        continue
          endif
        enddo
c suppress screened pair from the list
        do i=1,maproj
          if(lproj3(i).gt.1)then
            do l=1,lproj3(i)
              kl=kproj3(i,l)
            enddo
          endif
        enddo
        do i=1,maproj
          if(lproj3(i).gt.1)then
            n=2
            cont=lproj3(i).gt.1
            do while(cont)
              l=lproj3(i)
              kl=kproj3(i,l)
c suppress end of the list in order to have the last pair active
              do while(.not.lascr(kl).and.l.gt.1)
                kproj3(i,l)=0
                lproj3(i)=lproj3(i)-1
                l=lproj3(i)
                kl=kproj3(i,l)
              enddo
              cont=lproj3(i).gt.n
              if(cont)then
c compress list
                kn=kproj3(i,n)
                if(.not.lascr(kn))then
                  m=lproj3(i)   !last pair always active
                  km=kproj3(i,m)
                  kproj3(i,n)=km
                  kproj3(i,m)=0
                  lproj3(i)=lproj3(i)-1
                endif
              endif
              n=min(lproj3(i)-1,n)+1
              cont=lproj3(i).ne.n
            enddo
          endif
        enddo

c Target
c Check phi not to be in the range of a nucleon with smaller b
        do j=1,matarg
          if(ltarg3(j).gt.1)then
            kl=ktarg3(j,1)
            lascr(kl)=.true.
            do 400 l=ltarg3(j),2,-1
              kl=ktarg3(j,l)
              lascr(kl)=.true.
              if(bk(kl).ge.frdmzcut(b))then
                do m=1,l-1
                  km=ktarg3(j,m)
                  if(km.ne.kl.and.bk(km).ge.frdmzcut(b)
     &          .and.phik3(kl).ge.phik3(km)-thetak3(km)
     &          .and.phik3(kl).le.phik3(km)+thetak3(km))then
                    lascr(kl)=.false.
                    goto 400
                  endif
                enddo
              endif
 400        continue
          endif
        enddo
c suppress screened pair from the list
        do j=1,matarg
          if(ltarg3(j).gt.1)then
            n=2
            cont=ltarg3(j).gt.1
            do while(cont)
              l=ltarg3(j)
              kl=ktarg3(j,l)
c suppress end of the list in order to have the last pair active
              do while(.not.lascr(kl).and.l.gt.1)
                ktarg3(j,l)=0
                ltarg3(j)=ltarg3(j)-1
                l=ltarg3(j)
                kl=ktarg3(j,l)
              enddo
              cont=ltarg3(j).gt.n
              if(cont)then
c compress list
                kn=ktarg3(j,n)
                if(.not.lascr(kn))then
                  m=ltarg3(j)   !last pair always active
                  km=ktarg3(j,m)
                  ktarg3(j,n)=km
                  ktarg3(j,m)=0
                  ltarg3(j)=ltarg3(j)-1
                endif
              endif
              n=min(ltarg3(j)-1,n)+1
              cont=ltarg3(j).ne.n
            enddo
          endif
        enddo
      endif


      end

c-----------------------------------------------------------------------
      subroutine xGeometry(iii)
c-----------------------------------------------------------------------
      include 'epos.inc'
      include 'epos.incems'
      include 'epos.incsem'
      common/xgeom/nnn,naa(kollmx),nbb(kollmx)
      character*5 fmt1,fmt2

      if(iii.eq.0)then

      do k=1,kollmx
        naa(k)=0
      enddo
      nnn=0


      elseif(iii.eq.1)then

      ngl=0
      do k=1,koll
        r=bk(k)
        if(r.le.sqrt(sigine/10./pi))ngl=ngl+1
      enddo
      if(ngl.ne.0)then
        nnn=nnn+1
        naa(ngl)=naa(ngl)+1
      endif

      elseif(iii.eq.2)then

      if(xpar1.eq.0..and.xpar2.eq.0.)then
       print*,'---------- xGeometry -----------------------'
       return
      endif
      x1=1-0.01*xpar2
      x2=1-0.01*xpar1
      kmx=0
      nbb(1)=naa(1)
      do k=2,kollmx
        if(naa(k).ne.0.)kmx=k
        nbb(k)=nbb(k-1)+naa(k)
      enddo
      k1=0
      k2=0
      do k=1,kmx
        x=nbb(k)/float(nnn)
        if(x.lt.x1)k1=k
        if(x.lt.x2)k2=k
      enddo
      k1=k1+1
      k2=k2+1
      x=0
      av=0
      su=0
      p1=0.
      p2=0.
      do k=1,kmx
        xb=x
        x=nbb(k)/float(nnn)
        y=naa(k)/float(nnn)
        dx=x-xb
        p=0.
        if(k.eq.k1)then
          p=(x-x1)/dx
          p1=p
        elseif(k.eq.k2)then
          p=(x2-xb)/dx
          p2=p
        elseif(k.gt.k1.and.k.lt.k2)then
          p=1
        endif
      av=av+y*p*k
      su=su+y*p
      enddo
      av=av/su
      n1=nint(100*p1)
      n2=nint(100*p2)
      if(n1.eq.0)then
        k1=k1+1
        n1=100
      endif
      if(k1.le.9)fmt1='i1,4x'
      if(k1.gt.9.and.k1.le.99)fmt1='i2,3x'
      if(k1.gt.99.and.k1.le.999)fmt1='i3,2x'
      if(k1.gt.999.and.k1.le.9999)fmt1='i4,1x'
      if(k2.le.9)fmt2='i1,4x'
      if(k2.gt.9.and.k2.le.99)fmt2='i2,3x'
      if(k2.gt.99.and.k2.le.999)fmt2='i3,2x'
      if(k2.gt.999.and.k2.le.9999)fmt2='i4,1x'
      write(6,'(i4,a,i5,a,'//fmt1//',i6,a,i5,a,'//fmt2//',5x,a,f8.2)')
     &       nint(xpar2),':MIN',n1,'%',k1
     &      ,nint(xpar1),':MAX',n2,'%',k2   ,'av:',av
      endif

      end

c-----------------------------------------------------------------------
      function conbmx()
c-----------------------------------------------------------------------
      double precision om1intbc,p,eps
      include 'epos.inc'
      include 'epos.incsem'

      conbmx=0.
      b1=0.
      b2=7.
      eps=5.0d-3
      p=1.d0-dexp(-om1intbc(b2))
      if(p.gt.2.d0*eps)return

      ntry=0

10    ntry=ntry+1
      b=b1+.5*(b2-b1)

      p=(1.d0-dexp(-om1intbc(b)))

      if(p.gt.eps)then
       if(p.gt.2.d0*eps)then
        b1=b
       else
        conbmx=b
       return
       endif
      else
       if(p.lt.eps/5.d0)then
        b2=b
       else
        conbmx=b
        return
       endif
      endif

      if(ntry.le.1000)goto 10
      write(ifmt,*)'Too much try in conbmx ... bmax=',b
      conbmx=b
      return

      end

c-----------------------------------------------------------------------
      function conbmxndif()
c-----------------------------------------------------------------------
      double precision om1intbc,p,eps
      include 'epos.inc'
      include 'epos.incsem'


      iomegasave=iomega
      iomega=2
      conbmxndif=0.
      b1=0.
      b2=7.
      conbmxndif=b2
      eps=1d-10
      p=1.d0-dexp(-om1intbc(b2))
      if(p.gt.2.d0*eps)goto 100

      ntry=0

10    ntry=ntry+1
      b=b1+.5*(b2-b1)

      p=(1.d0-dexp(-om1intbc(b)))

      if(p.gt.eps)then
       if(p.gt.2.d0*eps)then
        b1=b
       else
        conbmxndif=b
       goto 100
       endif
      else
       if(p.lt.eps/5.d0)then
        b2=b
       else
        conbmxndif=b
        goto 100
       endif
      endif

      if(ntry.le.1000)goto 10
      write(ifmt,*)'Too much try in conbmxndif ... bkmxndif=',b
      conbmxndif=b
 100  iomega=iomegasave
      return

      end

c-----------------------------------------------------------------------
      function conbmxdif()
c-----------------------------------------------------------------------
c find b to have (1-exp(-om))pmax=pdiff
c-----------------------------------------------------------------------
      double precision om1intbc,pmax,drootom,pdiff
      include 'epos.inc'
      include 'epos.incsem'

      conbmxdif=0.
      b1=0.
      bmax=7.
      iomegasave=iomega
      iomega=2

      eps=1.e-5
      pmax=1.d0-dexp(-om1intbc(b1))
      pdiff=facdif
      if(pmax.gt.eps)then
        conbmxdif=drootom(pdiff,pmax,bmax,eps)
      endif
      iomega=iomegasave

      return

      end

c-----------------------------------------------------------------------
      subroutine conre
c-----------------------------------------------------------------------
c  initializes remnants
c-----------------------------------------------------------------------
      include "epos.incems"
      include 'epos.inc'

      call utpri('conre ',ish,ishini,6)

c     proj
c     ----
      la=laproj
      ma=iabs(maproj)
      las=0
      mas=0
      do l=1,ma
      if(la.lt.0)then
        if(iabs(idproj).lt.20)then
          id=idproj
        else
          ia=iabs(idproj/10)
          is=idproj/iabs(idproj)
          if(ia.ne.111.and.ia.ne.222.and.ia.ne.333)id=idproj/10*10
          if(ia.eq.111.or. ia.eq.222.or. ia.eq.333)id=idproj/10*10+is
          if(ia.eq.213)id=1230*is
        endif
      else
        id=1220
        if(rangen().le.(la-las)*1./(ma-mas))id=1120
        if(id.eq.1120)las=las+1
        mas=mas+1
      endif
      ic1=idtrai(1,id,1)
      ic2=idtrai(2,id,1)
      icproj(1,l)=ic1
      icproj(2,l)=ic2
      enddo

c     targ
c     ----
      la=latarg
      ma=iabs(matarg)
      las=0
      mas=0
      do l=1,ma
      if(la.lt.0)then
        if(iabs(idtarg).lt.20)then
          id=idtarg
        else
          ia=iabs(idtarg/10)
          is=idtarg/iabs(idtarg)
          if(ia.ne.111.and.ia.ne.222.and.ia.ne.333)id=idtarg/10*10
          if(ia.eq.111.or. ia.eq.222.or. ia.eq.333)id=idtarg/10*10+is
          if(ia.eq.213)id=1230*is
        endif
      else
        id=1220
        if(rangen().le.(la-las)*1./(ma-mas))id=1120
        if(id.eq.1120)las=las+1
        mas=mas+1
      endif
      ic1=idtrai(1,id,1)
      ic2=idtrai(2,id,1)
      ictarg(1,l)=ic1
      ictarg(2,l)=ic2
      enddo

      call utprix('conre ',ish,ishini,6)
      return
      end

c-----------------------------------------------------------------------
      subroutine conrl
c-----------------------------------------------------------------------
c  initializes target remnant in case of appl lepton
c-----------------------------------------------------------------------
      include "epos.incems"
      common/nucl1/laproj,maproj,latarg,matarg,core,fctrmx
      common/hadr2/iomodl,idproj,idtarg,wexcit

c     targ
c     ----
      la=latarg
      ma=iabs(matarg)
      las=0
      mas=0
      do l=1,ma
      if(la.lt.0)then
      id=idtarg
      else
      id=1220
      if(rangen().le.(la-las)*1./(ma-mas))id=1120
      if(id.eq.1120)las=las+1
      mas=mas+1
      endif
      ic1=idtrai(1,id,1)
      ic2=idtrai(2,id,1)
      ictarg(1,l)=ic1
      ictarg(2,l)=ic2
      enddo

      return
      end

c-----------------------------------------------------------------------
      subroutine conwr
c-----------------------------------------------------------------------
c     writes /cptl/
c-----------------------------------------------------------------------
      include "epos.inc"
      include "epos.incems"
      double precision XA(64,3),XB(64,3),BQGS,BMAXQGS,BMAXNEX,BMINNEX
      COMMON /Q_QGSNEX1/ XA,XB,BQGS,BMAXQGS,BMAXNEX,BMINNEX
      common/nucl3/phi,bimp
      common /cncl/xproj(mamx),yproj(mamx),zproj(mamx)
     *,xtarg(mamx),ytarg(mamx),ztarg(mamx)
      parameter(iapmax=208)
      double precision bqgs2,bmaxqgs2,bmaxnex2,bminnex2,xan,xbn
      common /qgsIInex1/xan(iapmax,3),xbn(iapmax,3)
     *,bqgs2,bmaxqgs2,bmaxnex2,bminnex2
      common/photrans/phoele(4),ebeam
      integer ic(2)

      call utpri('conwr ',ish,ishini,6)

      bx=cos(phi)*bimp
      by=sin(phi)*bimp

c     write /cptl/
c     ------------
      nptl=0

           if(iokoll.eq.1)then   ! precisely matarg collisions

      nptl=nptl+1
      do 3 i=1,4
3     xorptl(i,nptl)=0
      tivptl(1,nptl)=-ainfin
      tivptl(2,nptl)=0
      istptl(nptl)=1
      iorptl(nptl)=-1
      jorptl(nptl)=0
      do 1 k=1,koll
      nptl=nptl+1
      do 4 i=1,4
4     xorptl(i,nptl)=0
      tivptl(1,nptl)=-ainfin
      tivptl(2,nptl)=0
      istptl(nptl)=1
      iorptl(nptl)=-1
      jorptl(nptl)=0
1     continue

           elseif(iappl.ne.7)then

c             print *,'proj'

      do 6 i=1,maproj
      nptl=nptl+1
      istptl(nptl)=0
      iorptl(nptl)=0
      jorptl(nptl)=0
      if(model.eq.2)then       !QGSJet
      xproj(i)=XA(i,1)
      yproj(i)=XA(i,2)
      zproj(i)=XA(i,3)
      istptl(nptl)=1
      iorptl(nptl)=-1
      elseif(model.eq.7.or.model.eq.11)then       !QGSJetII
      xproj(i)=xan(i,1)
      yproj(i)=xan(i,2)
      zproj(i)=xan(i,3)
      istptl(nptl)=1
      iorptl(nptl)=-1
      elseif(model.ge.3)then       !Gheisha, ...
      istptl(nptl)=1
      iorptl(nptl)=-1
      endif
      xorptl(1,nptl)=xproj(i)+bx/2
      xorptl(2,nptl)=yproj(i)+by/2
      xorptl(3,nptl)=zproj(i)
      xorptl(4,nptl)=0
      tivptl(1,nptl)=-ainfin
c for visualisation uncomment
c-c   tivptl(1,nptl)=-100
      tivptl(2,nptl)= ainfin
c             print *,i,xorptl(1,nptl),xorptl(2,nptl),xorptl(3,nptl)
6     continue
c             print *,'targ'
      do 7 i=1,matarg
      nptl=nptl+1
      istptl(nptl)=0
      iorptl(nptl)=0
      jorptl(nptl)=0
      if(model.eq.2)then       !QGSJet
      xtarg(i)=XB(i,1)
      ytarg(i)=XB(i,2)
      ztarg(i)=XB(i,3)
      istptl(nptl)=1
      iorptl(nptl)=-1
      elseif(model.eq.7.or.model.eq.11)then       !QGSJetII
      xtarg(i)=xbn(i,1)
      ytarg(i)=xbn(i,2)
      ztarg(i)=xbn(i,3)
      istptl(nptl)=1
      iorptl(nptl)=-1
      elseif(model.ge.3)then       !Gheisha, ...
      istptl(nptl)=1
      iorptl(nptl)=-1
      endif
      xorptl(1,nptl)=xtarg(i)-bx/2
      xorptl(2,nptl)=ytarg(i)-by/2
      xorptl(3,nptl)=ztarg(i)
      xorptl(4,nptl)=0
      tivptl(1,nptl)=-ainfin
c for visualisation uncomment
c-c   tivptl(1,nptl)=-100
      tivptl(2,nptl)= ainfin
c             print *,i,xorptl(1,nptl),xorptl(2,nptl),xorptl(3,nptl)
7     continue
      if(abs(idprojin).eq.12)then   !electron for fake DIS
c electron projectile
        nptl=nptl+1
        istptl(nptl)=41
        iorptl(nptl)=-1
        jorptl(nptl)=-1
        iorptl(1)=nptl         !pi0 (porjectile) coming from lepton
        xorptl(1,nptl)=bx/2
        xorptl(2,nptl)=by/2
        xorptl(3,nptl)=0.
        xorptl(4,nptl)=0.
        tivptl(1,nptl)=-ainfin
        tivptl(2,nptl)=0.
c target nucleons (in lab frame)
        do i=1,matarg
          nptl=nptl+1
          istptl(nptl)=41
          iorptl(nptl)=-1
          jorptl(nptl)=-1
          xorptl(1,nptl)=xtarg(i)-bx/2
          xorptl(2,nptl)=ytarg(i)-by/2
          xorptl(3,nptl)=ztarg(i)
          xorptl(4,nptl)=0
          tivptl(1,nptl)=-ainfin
c         for visualisation uncomment
c         -c   tivptl(1,nptl)=-100
          tivptl(2,nptl)= ainfin
        enddo
c electron remnant
        nptl=nptl+1
        istptl(nptl)=0
        iorptl(nptl)=maproj+matarg+1
        jorptl(nptl)=-1
        xorptl(1,nptl)=bx/2
        xorptl(2,nptl)=by/2
        xorptl(3,nptl)=0.
        xorptl(4,nptl)=0.
        tivptl(1,nptl)=0.
        tivptl(2,nptl)= ainfin
      endif

          endif

      nptl=0
      if(iappl.le.2)then
      do i=1,maproj
      nptl=nptl+1
      ic(1)=icproj(1,i)
      ic(2)=icproj(2,i)
      id=idtra(ic,0,0,0)
c      id=idtra(ic,0,0,3)      !tp071107 imix=3 ??????????
      call idmass(id,ams)
      idptl(nptl)=id
      pptl(1,nptl)=0.
      pptl(2,nptl)=0.
      pptl(3,nptl)=pnullx
      pptl(4,nptl)=sqrt(pnullx**2+ams**2)
      pptl(5,nptl)=ams
      ifrptl(1,nptl)=0
      ifrptl(2,nptl)=0
      ityptl(nptl)=0
      enddo
      endif
      if(iappl.ne.7)then
      do i=1,matarg
      nptl=nptl+1
      ic(1)=ictarg(1,i)
      ic(2)=ictarg(2,i)
      id=idtra(ic,0,0,0)
c      id=idtra(ic,0,0,3)      !tp071107 imix=3 ??????????
      call idmass(id,ams)
      idptl(nptl)=id
      pptl(1,nptl)=0.
      pptl(2,nptl)=0.
      pptl(3,nptl)=-pnullx
      pptl(4,nptl)=sqrt(pnullx**2+ams**2)
      pptl(5,nptl)=ams
      ifrptl(1,nptl)=0
      ifrptl(2,nptl)=0
      ityptl(nptl)=0
      enddo
      if(abs(idprojin).eq.12)then   !electron for fake DIS
c electron projectile
        nptl=nptl+1
        id=idprojin
        call idmass(id,ams)
        idptl(nptl)=id
        pptl(1,nptl)=0.
        pptl(2,nptl)=0.
        pptl(3,nptl)=sqrt(max(0.,(elepti+ams)*(elepti-ams)))
        pptl(4,nptl)=elepti
        pptl(5,nptl)=ams
        ifrptl(1,nptl)=1
        ifrptl(2,nptl)=1
        ityptl(nptl)=40
c target nucleons (in lab frame)
        do i=1,matarg
          nptl=nptl+1
          idptl(nptl)=idptl(maproj+i)
          pptl(1,nptl)=0.
          pptl(2,nptl)=0.
          pptl(3,nptl)=-pnll
          pptl(4,nptl)=ebeam
          pptl(5,nptl)=pptl(5,maproj+i)
          ifrptl(1,nptl)=maproj+i
          ifrptl(2,nptl)=maproj+i
          ityptl(nptl)=50
        enddo
c electron remnant
        nptl=nptl+1
        idptl(nptl)=id
        pptl(1,nptl)=phoele(1)
        pptl(2,nptl)=phoele(2)
        pptl(3,nptl)=phoele(3)
        pptl(4,nptl)=phoele(4)
        pptl(5,nptl)=ams
        ifrptl(1,nptl)=0
        ifrptl(2,nptl)=0
        ityptl(nptl)=40
      endif

      else

      nptl=nptl+1
      id=idproj
      call idmass(id,ams)
      idptl(nptl)=id
      pptl(1,nptl)=0.
      pptl(2,nptl)=0.
      pptl(3,nptl)=pnullx
      pptl(4,nptl)=sqrt(pnullx**2+ams**2)
      pptl(5,nptl)=ams
      ifrptl(1,nptl)=0
      ifrptl(2,nptl)=0
      ityptl(nptl)=0
      iorptl(nptl)=-1
      jorptl(nptl)=0
      istptl(nptl)=0
      do 5 i=1,4
 5      xorptl(i,nptl)=0
      tivptl(1,nptl)=0
      tivptl(2,nptl)=0
      endif

c     exit
c     ----

      call utprix('conwr ',ish,ishini,6)
      return
      end

c------------------------------------------------------------------------
      subroutine conxyz(ch,n,x,y,z,ynuc)
c-----------------------------------------------------------------------
      include 'epos.inc'

      real x(n),y(n),z(n)
      character ch*1

      massnr=0
      iii=0
      if(ch.eq.'p')then
      massnr=maproj
      iii=1
      elseif(ch.eq.'t')then
      massnr=matarg
      iii=2
      else
      call utstop('conxyz: nucleus neither proj nor targ&')
      endif

      if(massnr.eq.0)return
      if(massnr.gt.n)call utstop('conxyz: massnr.gt.n&')
      if(massnr.eq.1)then
      x(1)=0
      y(1)=0
      z(1)=0
      return
      endif

      rad=radnuc(massnr)

      if(massnr.ge.10)then !---wood-saxon density---

        rad=rad/difnuc(massnr)
        cr1=1.+3./rad+6./rad**2+6./rad**3
        cr2=3./rad
        cr3=3./rad+6./rad**2
        do i=1,massnr
   1      zuk=rangen()*cr1-1.
          if(zuk.le.0.)then
            tt=rad*(rangen()**.3333-1.)
          elseif(zuk.le.cr2 )then
            tt=-log(rangen())
          elseif(zuk.lt.cr3 )then
            tt=-log(rangen())-log(rangen())
          else
            tt=-log(rangen())-log(rangen())-log(rangen())
          endif
          if(rangen().gt.1./(1.+exp(-abs(tt))))goto 1
          rim=tt+rad
          zz=rim*(2.*rangen()-1.)
          rim=sqrt(rim*rim-zz*zz)
          z(i)=zz*difnuc(massnr)
          call pscs(c,s)
          x(i)=rim*c*difnuc(massnr)
          y(i)=rim*s*difnuc(massnr)
        enddo

      elseif(massnr.ge.3)then  ! ---gaussian density---

        rad=rad*sqrt(2.*massnr/(massnr-1.))   !van hove simulation
        do l=1,3
          summ=0.
          do i=1,massnr-1
            j=massnr-i
            aks=rad *(rangen()+rangen()+rangen()-1.5)
            k=j+1
            if(l.eq.1)x(k)=summ-aks*sqrt(float(j)/k)
            if(l.eq.2)y(k)=summ-aks*sqrt(float(j)/k)
            if(l.eq.3)z(k)=summ-aks*sqrt(float(j)/k)
            summ=summ+aks/sqrt(float(j*k))
          enddo
          if(l.eq.1)x(1)=summ
          if(l.eq.2)y(1)=summ
          if(l.eq.3)z(1)=summ
        enddo

      elseif(massnr.eq.2)then  ! ---deuteron---

        !.........r=t*difnuc(massnr), t~exp(-2*t)*(1-exp(-a*t))
        a=radnuc(massnr)
  2     t=-0.5*alog(rangen())  !~exp(-2*t)
        if(rangen().gt.(1-exp(-a*t))**2)goto2
        r=t*difnuc(massnr)
        zz=r*(2.*rangen()-1.)
        call pscs(c,s)
        rxy=sqrt(r*r-zz*zz)
        z(1)=0.5*zz
        x(1)=0.5*rxy*c
        y(1)=0.5*rxy*s
        z(2)=-z(1)
        x(2)=-x(1)
        y(2)=-y(1)

      else

        stop'conxyz: wrong massnr.     '

      endif

c...plot preparation

      rmax=(radnuc(massnr)+3)
      drnucl(iii)=rmax/mxnucl
      nrnucl(iii)=nrnucl(iii)+1
      do i=1,massnr
        r=sqrt(x(i)**2+y(i)**2+z(i)**2)
        k=1+int(r/drnucl(iii))
        if(k.le.mxnucl)rnucl(k,iii)=rnucl(k,iii)+1
      enddo

c...lorentz trafo

      do i=1,massnr
      z(i)=z(i)/cosh(ynuc)
      enddo

      return
      end

c-----------------------------------------------------------------------
      subroutine conini
c-----------------------------------------------------------------------
      include 'epos.inc'

      imax=max(maproj,matarg)
      if(idtargin.eq.0)imax=max(imax,40)
      do massnr=1,mamxx
        dif=0.54
        rad=0.
        if(massnr.gt.imax.or.massnr.eq.1)then
          dif=0
          rad=0
        elseif(massnr.eq.197)then
          dif=0.562
          rad=6.5
        elseif(massnr.ge.10)then
          rad=1.12*massnr**0.33333-0.86*massnr**(-0.33333)
        elseif(massnr.ge.3)then
          rad=.9*float(massnr)**.3333
        elseif(massnr.eq.2)then
          dif=4.316
          rad=4.68
        endif
        difnuc(massnr)=dif
        radnuc(massnr)=rad
      enddo

      end

c-----------------------------------------------------------------------
      subroutine xConNuclDens(iii)
c-----------------------------------------------------------------------
c plots distribution of nucleons in nuclei
c  iii = 1 (proj) or 2 (targ)
c-----------------------------------------------------------------------
      include 'epos.inc'
      if(model.ne.1)return
      massnr=1
      if(iii.eq.1)then
        massnr=maproj
      elseif(iii.eq.2)then
        massnr=matarg
      endif
      if(massnr.eq.1)return
      a=1./4.316
      b=4.68
      write(ifhi,'(a)') '!-----------------------------------------'
      write(ifhi,'(a)') '!          nuclear density          '
      write(ifhi,'(a)') '!-----------------------------------------'
      write(ifhi,'(a)')       'openhisto'
      if(massnr.ge.10)write(ifhi,'(a)')'htyp lin xmod lin ymod lin'
      if(massnr.lt.10)write(ifhi,'(a)')'htyp lin xmod lin ymod log'
      write(ifhi,'(a)')       'text 0 0 "title nuclear density"'
      write(ifhi,'(a)')       'text 0.99 -0.15 " r (fm) "'
      write(ifhi,'(a)')       'text 0 0 "yaxis rho(r)"'
      write(ifhi,'(a,2e11.3)')'xrange',0.,mxnucl*drnucl(iii)
      write(ifhi,'(3a)')'yrange',' 0 ',' auto'
      write(ifhi,'(a)')       'array 2'
      do j=1,mxnucl
      r=(j-0.5)*drnucl(iii)
      d=0.5*drnucl(iii)
      write(ifhi,'(2e12.4)')  r,rnucl(j,iii)/nrnucl(iii)/
     *                     (4./3.*pi*((r+d)**3-(r-d)**3))
      enddo
      write(ifhi,'(a)')       '  endarray'
      write(ifhi,'(a)')       'closehisto plot 0-'
      write(ifhi,'(a)')       'openhisto'
      write(ifhi,'(a)')       'htyp lbo '
      write(ifhi,'(a)')       'array 2'
      do j=1,mxnucl
      r=(j-0.5)*drnucl(iii)
      rr=2*r
      rho=1
      if(massnr.eq.2)then
        rho=1.00*((1-exp(-b*a*rr))*exp(-a*rr)/rr)**2
      elseif(massnr.eq.197)then
        rho=0.16/(1+exp((r-6.5)/0.562))
      elseif(massnr.ge.10)then
        rho=0.16/(1+exp((r-radnuc(massnr))/difnuc(massnr)))
      endif
      write(ifhi,'(2e12.4)')  r,rho
      enddo
      write(ifhi,'(a)')       '  endarray'
      write(ifhi,'(a)')       'closehisto plot 0'
      end

c-----------------------------------------------------------------------
      subroutine xConThick(iii)
c-----------------------------------------------------------------------
      ! plots sigma_pp *T_A (b)  (=average number of collisions)
      ! T_A = thickness function
      !  iii = 1 (proj) or 2 (targ)
      !----------------------------------------------------------------
      include 'epos.inc'
      parameter(iconimax=20,iconkmax=100)
      real thick(2,0:iconimax)
      imax=iconimax
      kmax=iconkmax
      if(model.ne.1)return
      ramx=mxnucl*drnucl(iii)
      do i=0,imax
        x=i/float(imax)*ramx
        sum=0
        rho0=conrho(iii,0.)
        h=ramx/kmax
        do k=1,kmax
          z=k*h
          r=sqrt(x*x+z*z)
          rho2=conrho(iii,r)
          z=(k-0.5)*h
          r=sqrt(x*x+z*z)
          rho1=conrho(iii,r)
          sum=sum+h/6.*(rho0+4*rho1+rho2)
          rho0=rho2
        enddo
        sum=sum*2 ! integral fro -infty to +infty
        thick(1,i)=x
        thick(2,i)=sum
      enddo
      write(ifhi,'(a)')       'openhisto'
      write(ifhi,'(a)')       'htyp lbo '
      write(ifhi,'(a)')       'txt "xaxis b (fm)" '
      write(ifhi,'(a)')       'txt "yaxis [s]?pp! T?A! (b) " '
      write(ifhi,'(a)')       'array 2'
      do i=0,imax
        write(ifhi,'(2e12.4)') thick(1,i),sigine/10.*thick(2,i)
      enddo
      write(ifhi,'(a)')       '  endarray'
      write(ifhi,'(a)')       'closehisto plot 0'

      end

c-----------------------------------------------------------------------
      function conrho(iii,r)
c-----------------------------------------------------------------------
      ! nuclear density
      !  iii = 1 (proj) or 2 (targ)
      !----------------------------------------------------------------
      include 'epos.inc'
      conrho=1.
      if(model.ne.1)return
      massnr=1
      if(iii.eq.1)then
        massnr=maproj
      elseif(iii.eq.2)then
        massnr=matarg
      endif
      if(massnr.eq.1)return
      a=1./4.316
      b=4.68
      rr=2*r
      rho=1
      if(massnr.eq.2.and.rr.gt.0.)then
        rho=1.00*((1-exp(-b*a*rr))*exp(-a*rr)/rr)**2
      elseif(massnr.eq.197)then
        rho=0.16/(1+exp((r-6.5)/0.562))
      elseif(massnr.ge.10)then
        rho=0.16/(1+exp((r-radnuc(massnr))/difnuc(massnr)))
      endif
      conrho=rho
      end







