Back to home page

Project CMSSW displayed by LXR

 
 

    


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

0001 c-----------------------------------------------------------------------
0002       subroutine conaa(iret)
0003 c-----------------------------------------------------------------------
0004 c  determines interaction configuration
0005 c-----------------------------------------------------------------------
0006 
0007       include 'epos.inc'
0008       include 'epos.incems'
0009       include 'epos.incsem'
0010       include 'epos.incpar'
0011 
0012       common/geom/rmproj,rmtarg,bmax,bkmx
0013       common/nucl3/phi,bimp
0014       common /cncl/xproj(mamx),yproj(mamx),zproj(mamx)
0015      *            ,xtarg(mamx),ytarg(mamx),ztarg(mamx)
0016       common/cfacmss/facmss
0017       double precision yyrmax
0018       common/scrangle/ phik3(kollmx),thetak3(kollmx)
0019 
0020       call utpri('conaa ',ish,ishini,4)
0021 
0022       iret=0
0023 
0024 c     initialisations
0025 c     ---------------
0026 
0027       vel=tanh(ypjtl-yhaha)+tanh(yhaha)
0028 
0029 c     determine phi, bimp, coll, iproj, itarg, x/y/zproj, x/y/ztarg
0030 c     ---------------------------------------------------------------
0031 
0032            if(iokoll.eq.1)then
0033 
0034       koll=matarg
0035       do k=1,koll
0036         do n=1,4
0037           coord(n,k)=0.
0038         enddo
0039       enddo
0040       bimp=0
0041       phi=0
0042       xproj(1)=0
0043       yproj(1)=0
0044       zproj(1)=0
0045       lproj(1)=koll
0046       lproj3(1)=0
0047       do k=1,koll
0048         bij=bkmx*sqrt(rangen())
0049         bk(k)=bij
0050         iproj(k)=1
0051         itarg(k)=k
0052         phi=2.*pi*rangen()
0053         xtarg(k)=bij*cos(phi)
0054         ytarg(k)=bij*sin(phi)
0055         ztarg(k)=0
0056         ltarg(k)=1
0057         kproj(1,k)=k
0058         ktarg(k,1)=k
0059         ltarg3(k)=0
0060         ktarg3(k,1)=0
0061         kproj3(k,1)=0
0062         if(iscreen.ne.0.and.bij.le.bkmxndif)then
0063           if(zbrmax.le.0..or. bij.lt.zbcutx+zbrmax*rangen())then
0064             lproj3(1)=lproj3(1)+1
0065             ltarg3(k)=1
0066             kproj3(1,lproj3(1))=k
0067             ktarg3(k,1)=k
0068           endif
0069         endif
0070       enddo
0071 
0072            elseif(maproj.eq.1.and.matarg.eq.1)then
0073 
0074       b1=bminim
0075       b2=amin1(bkmx,bmaxim)
0076       if(b1.gt.b2)call utstop('conaa: bmin > bmax&')
0077       bimp=sqrt(b1*b1+(b2*b2-b1*b1)*rangen())
0078       koll=1
0079       do n=1,4
0080         coord(n,1)=0.
0081       enddo
0082       bk(1)=bimp
0083       iproj(1)=1
0084       itarg(1)=1
0085       phi=2.*pi*rangen()
0086       xproj(1)=bimp*cos(phi)
0087       yproj(1)=bimp*sin(phi)
0088       zproj(1)=0
0089       xtarg(1)=0
0090       ytarg(1)=0
0091       ztarg(1)=0
0092       lproj(1)=1
0093       ltarg(1)=1
0094       lproj3(1)=1
0095       ltarg3(1)=1
0096       kproj3(1,1)=1
0097       ktarg3(1,1)=1
0098       kproj(1,1)=1
0099       ktarg(1,1)=1
0100 
0101            else
0102 
0103       call conxyz('p',mamx,xproj,yproj,zproj,ypjtl-yhaha)
0104       call conxyz('t',mamx,xtarg,ytarg,ztarg,yhaha)
0105 
0106       bx=0
0107       by=0
0108       if(maproj.gt.0)then
0109       if(bimevt.lt.0)then
0110         b1=bminim
0111         b2=amin1(rmproj+rmtarg,bmaxim)
0112         if(b1.gt.b2)call utstop('conaa: bmin > bmax&')
0113         bimp=sqrt(b1**2+(b2**2-b1**2)*rangen())
0114         if(nbarray.gt.0)bimp=barray(mod(nrevt,nbarray)+1)
0115         if(jpsi.gt.0)then
0116           bimp=b1+(b2-b1)*(float(mod(nrevt,12))+rangen())/12.
0117           bimevt=bimp
0118         endif
0119         phi=phimin+rangen()*(phimax-phimin)
0120       else
0121         phi=phievt
0122         bimp=bimevt
0123       endif
0124       bx=cos(phi)*bimp
0125       by=sin(phi)*bimp
0126       endif
0127       if(jpsi.lt.0)then !modify b
0128         bx=xtarg(1)
0129         by=ytarg(1)
0130       endif
0131       if(maproj.eq.0)goto1000
0132       koll=0
0133       do i=1,maproj
0134       lproj(i)=0
0135       lproj3(i)=0
0136       enddo
0137       do j=1,matarg
0138       ltarg(j)=0
0139       ltarg3(j)=0
0140       enddo
0141       do 12 i=1,maproj
0142       do 11 j=1,matarg
0143       if(jpsi.lt.0.and.ztarg(j).le.ztarg(1))goto11
0144       bij=sqrt((xproj(i)+bx-xtarg(j))**2+(yproj(i)+by-ytarg(j))**2)
0145       if(ish.ge.7)write(ifch,*)'i_p:',i,' i_t:',j,' b_ij:',bij
0146       if(bij.gt.bkmx)goto 11
0147 
0148       koll=koll+1
0149       if(koll.gt.kollmx)call utstop('conaa: kollmx too small&')
0150       bk(koll)=bij
0151       bkx(koll)=xproj(i)+bx-xtarg(j)
0152       bky(koll)=yproj(i)+by-ytarg(j)
0153       iproj(koll)=i
0154       itarg(koll)=j
0155       lproj(i)=lproj(i)+1
0156       ltarg(j)=ltarg(j)+1
0157       kproj(i,lproj(i))=koll
0158       ktarg(j,ltarg(j))=koll
0159       phik3(koll)=0.
0160       thetak3(koll)=0.
0161       if(iscreen.ne.0.and.bij.le.bkmxndif)then
0162         if(zbrmax.gt.0..and.bij.gt.zbcutx+zbrmax*rangen())goto 11
0163         lproj3(i)=lproj3(i)+1
0164         ltarg3(j)=ltarg3(j)+1
0165         kproj3(i,lproj3(i))=koll
0166         ktarg3(j,ltarg3(j))=koll
0167 c define angle for anti-shadowing
0168         if(abs(bky(koll)).gt.1.e-6)then
0169           if(abs(bkx(koll)).gt.1.e-6)then
0170             phik3(koll)=atan(bky(koll)/bkx(koll))
0171           else
0172             phik3(koll)=sign(0.5*pi,bky(koll))
0173           endif
0174         elseif(bkx(koll).lt.0.)then
0175           phik3(koll)=pi
0176         endif
0177         if(bk(koll).gt.0.)then
0178           thetak3(koll)=atan(bglaubx/bk(koll))
0179         else
0180           thetak3(koll)=0.5*pi
0181         endif
0182       endif
0183 
0184 11    continue
0185 12    continue
0186 
0187       do k=1,koll
0188         do n=1,4
0189           coord(n,k)=0.
0190         enddo
0191       enddo
0192 
0193 
0194           endif
0195 
0196       if(ish.ge.3)write(ifch,*)'koll=',koll
0197       if(koll.eq.0)goto 1001
0198 
0199 
0200 c     determine coord
0201 c     ---------------
0202       do kl=1,koll
0203       i=iproj(kl)
0204       j=itarg(kl)
0205       dist=ztarg(j)-zproj(i)
0206       coord(1,kl)=(xproj(i)+xtarg(j))*0.5
0207       coord(2,kl)=(yproj(i)+ytarg(j))*0.5
0208       coord(3,kl)=(zproj(i)+ztarg(j))*0.5
0209       coord(4,kl)=dist/vel
0210       enddo
0211 
0212       if(iscreen.ne.0)call CalcScrPair(bimp)
0213 
0214     !~~~~~redefine energy in case of imposed radial flow~~~~~~~~~~~~~~~~
0215       yrmaxi=max(0.,yradmx+yradmi*log(1.+engy/sqrt(float(koll))))
0216       if(yrmaxi.gt.1e-5)then
0217         yyrmax=dble(yrmaxi)
0218         fradflii=sngl(1d0/
0219      &  ((sinh(yyrmax)*yyrmax-cosh(yyrmax)+1d0)/(yyrmax**2/2d0)))
0220       else
0221         fradflii=1.
0222       endif
0223       if(ish.ge.3)write(ifch,*)'yrmaxi=',yrmaxi
0224 
0225 
0226 c     exit
0227 c     ----
0228 1000  continue
0229       if(ish.ge.5)then
0230       write(ifch,*)'ia,x/y/zproj:'
0231       do mm=1,maproj
0232       write(ifch,*)mm,xproj(mm),yproj(mm),zproj(mm)
0233       enddo
0234       write(ifch,*)'ia,x/y/ztarg:'
0235       do mm=1,matarg
0236       write(ifch,*)mm,xtarg(mm),ytarg(mm),ztarg(mm)
0237       enddo
0238       write(ifch,*)'iret',iret
0239       endif
0240       call utprix('conaa ',ish,ishini,4)
0241       return
0242 
0243 1001  continue !iret=1 causes redo of whole collision
0244       iret=1
0245       if(ish.ge.3)then
0246       write(ifch,*)
0247       write(ifch,*)'***** subroutine conaa:'
0248       write(ifch,*)'***** no nucleon pair found --> no interaction'
0249       write(ifch,*)
0250       endif
0251       goto 1000
0252 
0253       end
0254 
0255 c-----------------------------------------------------------------------
0256       function frdmzcut(b)
0257 c-----------------------------------------------------------------------
0258 c  determines interaction configuration
0259 c-----------------------------------------------------------------------
0260 
0261       include 'epos.inc'
0262       include 'epos.incpar'
0263 
0264       frdmzcut=zbcutx*exp(-(b*rangen()))
0265 c      frdmzcut=zbcutx*rangen()
0266 c      bdum=b
0267       return
0268       end
0269 
0270 c-----------------------------------------------------------------------
0271       subroutine CalcScrPair(b)
0272 c-----------------------------------------------------------------------
0273 c  determines interaction configuration
0274 c-----------------------------------------------------------------------
0275 
0276       include 'epos.inc'
0277       include 'epos.incems'
0278       include 'epos.incpar'
0279 
0280       common/scrangle/ phik3(kollmx),thetak3(kollmx)
0281       logical lascr(kollmx),cont
0282 
0283 
0284 c     order pairs for splitting as a function of b
0285 c     ---------------
0286       do i=1,maproj
0287         if(lproj3(i).gt.1)then
0288           do l=2,lproj3(i)
0289             m=l
0290  100        continue
0291               kp=kproj3(i,m-1)
0292               kl=kproj3(i,m)
0293               if(bk(kl).lt.bk(kp))then
0294                 kproj3(i,m-1)=kl
0295                 kproj3(i,m)=kp
0296                 m=m-1
0297             if(m.gt.1)goto 100
0298               endif
0299           enddo
0300         endif
0301       enddo
0302       do j=1,matarg
0303         if(ltarg3(j).gt.1)then
0304           do l=2,ltarg3(j)
0305             m=l
0306  200        continue
0307               kp=ktarg3(j,m-1)
0308               kl=ktarg3(j,m)
0309               if(bk(kl).lt.bk(kp))then
0310                 ktarg3(j,m-1)=kl
0311                 ktarg3(j,m)=kp
0312                 m=m-1
0313             if(m.gt.1)goto 200
0314               endif
0315           enddo
0316         endif
0317       enddo
0318 
0319       if(koll.gt.1)then
0320 c Define anti-shadowing as a consequence of geometrical screening of
0321 c nucleons at large b by the one at small b
0322 
0323 c Projectile
0324 c Check phi not to be in the range of a nucleon with smaller b
0325         do i=1,maproj
0326           if(lproj3(i).gt.1)then
0327             kl=kproj3(i,1)
0328             lascr(kl)=.true.
0329             do 300 l=lproj3(i),2,-1
0330               kl=kproj3(i,l)
0331               lascr(kl)=.true.
0332               if(bk(kl).ge.frdmzcut(b))then
0333                 do m=1,l-1
0334                   km=kproj3(i,m)
0335                   if(kl.ne.km.and.bk(km).ge.frdmzcut(b)
0336      &          .and.phik3(kl).ge.phik3(km)-thetak3(km)
0337      &          .and.phik3(kl).le.phik3(km)+thetak3(km))then
0338                     lascr(kl)=.false.
0339                     goto 300
0340                   endif
0341                 enddo
0342               endif
0343  300        continue
0344           endif
0345         enddo
0346 c suppress screened pair from the list
0347         do i=1,maproj
0348           if(lproj3(i).gt.1)then
0349             do l=1,lproj3(i)
0350               kl=kproj3(i,l)
0351             enddo
0352           endif
0353         enddo
0354         do i=1,maproj
0355           if(lproj3(i).gt.1)then
0356             n=2
0357             cont=lproj3(i).gt.1
0358             do while(cont)
0359               l=lproj3(i)
0360               kl=kproj3(i,l)
0361 c suppress end of the list in order to have the last pair active
0362               do while(.not.lascr(kl).and.l.gt.1)
0363                 kproj3(i,l)=0
0364                 lproj3(i)=lproj3(i)-1
0365                 l=lproj3(i)
0366                 kl=kproj3(i,l)
0367               enddo
0368               cont=lproj3(i).gt.n
0369               if(cont)then
0370 c compress list
0371                 kn=kproj3(i,n)
0372                 if(.not.lascr(kn))then
0373                   m=lproj3(i)   !last pair always active
0374                   km=kproj3(i,m)
0375                   kproj3(i,n)=km
0376                   kproj3(i,m)=0
0377                   lproj3(i)=lproj3(i)-1
0378                 endif
0379               endif
0380               n=min(lproj3(i)-1,n)+1
0381               cont=lproj3(i).ne.n
0382             enddo
0383           endif
0384         enddo
0385 
0386 c Target
0387 c Check phi not to be in the range of a nucleon with smaller b
0388         do j=1,matarg
0389           if(ltarg3(j).gt.1)then
0390             kl=ktarg3(j,1)
0391             lascr(kl)=.true.
0392             do 400 l=ltarg3(j),2,-1
0393               kl=ktarg3(j,l)
0394               lascr(kl)=.true.
0395               if(bk(kl).ge.frdmzcut(b))then
0396                 do m=1,l-1
0397                   km=ktarg3(j,m)
0398                   if(km.ne.kl.and.bk(km).ge.frdmzcut(b)
0399      &          .and.phik3(kl).ge.phik3(km)-thetak3(km)
0400      &          .and.phik3(kl).le.phik3(km)+thetak3(km))then
0401                     lascr(kl)=.false.
0402                     goto 400
0403                   endif
0404                 enddo
0405               endif
0406  400        continue
0407           endif
0408         enddo
0409 c suppress screened pair from the list
0410         do j=1,matarg
0411           if(ltarg3(j).gt.1)then
0412             n=2
0413             cont=ltarg3(j).gt.1
0414             do while(cont)
0415               l=ltarg3(j)
0416               kl=ktarg3(j,l)
0417 c suppress end of the list in order to have the last pair active
0418               do while(.not.lascr(kl).and.l.gt.1)
0419                 ktarg3(j,l)=0
0420                 ltarg3(j)=ltarg3(j)-1
0421                 l=ltarg3(j)
0422                 kl=ktarg3(j,l)
0423               enddo
0424               cont=ltarg3(j).gt.n
0425               if(cont)then
0426 c compress list
0427                 kn=ktarg3(j,n)
0428                 if(.not.lascr(kn))then
0429                   m=ltarg3(j)   !last pair always active
0430                   km=ktarg3(j,m)
0431                   ktarg3(j,n)=km
0432                   ktarg3(j,m)=0
0433                   ltarg3(j)=ltarg3(j)-1
0434                 endif
0435               endif
0436               n=min(ltarg3(j)-1,n)+1
0437               cont=ltarg3(j).ne.n
0438             enddo
0439           endif
0440         enddo
0441       endif
0442 
0443 
0444       end
0445 
0446 c-----------------------------------------------------------------------
0447       subroutine xGeometry(iii)
0448 c-----------------------------------------------------------------------
0449       include 'epos.inc'
0450       include 'epos.incems'
0451       include 'epos.incsem'
0452       common/xgeom/nnn,naa(kollmx),nbb(kollmx)
0453       character*5 fmt1,fmt2
0454 
0455       if(iii.eq.0)then
0456 
0457       do k=1,kollmx
0458         naa(k)=0
0459       enddo
0460       nnn=0
0461 
0462 
0463       elseif(iii.eq.1)then
0464 
0465       ngl=0
0466       do k=1,koll
0467         r=bk(k)
0468         if(r.le.sqrt(sigine/10./pi))ngl=ngl+1
0469       enddo
0470       if(ngl.ne.0)then
0471         nnn=nnn+1
0472         naa(ngl)=naa(ngl)+1
0473       endif
0474 
0475       elseif(iii.eq.2)then
0476 
0477       if(xpar1.eq.0..and.xpar2.eq.0.)then
0478        print*,'---------- xGeometry -----------------------'
0479        return
0480       endif
0481       x1=1-0.01*xpar2
0482       x2=1-0.01*xpar1
0483       kmx=0
0484       nbb(1)=naa(1)
0485       do k=2,kollmx
0486         if(naa(k).ne.0.)kmx=k
0487         nbb(k)=nbb(k-1)+naa(k)
0488       enddo
0489       k1=0
0490       k2=0
0491       do k=1,kmx
0492         x=nbb(k)/float(nnn)
0493         if(x.lt.x1)k1=k
0494         if(x.lt.x2)k2=k
0495       enddo
0496       k1=k1+1
0497       k2=k2+1
0498       x=0
0499       av=0
0500       su=0
0501       p1=0.
0502       p2=0.
0503       do k=1,kmx
0504         xb=x
0505         x=nbb(k)/float(nnn)
0506         y=naa(k)/float(nnn)
0507         dx=x-xb
0508         p=0.
0509         if(k.eq.k1)then
0510           p=(x-x1)/dx
0511           p1=p
0512         elseif(k.eq.k2)then
0513           p=(x2-xb)/dx
0514           p2=p
0515         elseif(k.gt.k1.and.k.lt.k2)then
0516           p=1
0517         endif
0518       av=av+y*p*k
0519       su=su+y*p
0520       enddo
0521       av=av/su
0522       n1=nint(100*p1)
0523       n2=nint(100*p2)
0524       if(n1.eq.0)then
0525         k1=k1+1
0526         n1=100
0527       endif
0528       if(k1.le.9)fmt1='i1,4x'
0529       if(k1.gt.9.and.k1.le.99)fmt1='i2,3x'
0530       if(k1.gt.99.and.k1.le.999)fmt1='i3,2x'
0531       if(k1.gt.999.and.k1.le.9999)fmt1='i4,1x'
0532       if(k2.le.9)fmt2='i1,4x'
0533       if(k2.gt.9.and.k2.le.99)fmt2='i2,3x'
0534       if(k2.gt.99.and.k2.le.999)fmt2='i3,2x'
0535       if(k2.gt.999.and.k2.le.9999)fmt2='i4,1x'
0536       write(6,'(i4,a,i5,a,'//fmt1//',i6,a,i5,a,'//fmt2//',5x,a,f8.2)')
0537      &       nint(xpar2),':MIN',n1,'%',k1
0538      &      ,nint(xpar1),':MAX',n2,'%',k2   ,'av:',av
0539       endif
0540 
0541       end
0542 
0543 c-----------------------------------------------------------------------
0544       function conbmx()
0545 c-----------------------------------------------------------------------
0546       double precision om1intbc,p,eps
0547       include 'epos.inc'
0548       include 'epos.incsem'
0549 
0550       conbmx=0.
0551       b1=0.
0552       b2=7.
0553       eps=5.0d-3
0554       p=1.d0-dexp(-om1intbc(b2))
0555       if(p.gt.2.d0*eps)return
0556 
0557       ntry=0
0558 
0559 10    ntry=ntry+1
0560       b=b1+.5*(b2-b1)
0561 
0562       p=(1.d0-dexp(-om1intbc(b)))
0563 
0564       if(p.gt.eps)then
0565        if(p.gt.2.d0*eps)then
0566         b1=b
0567        else
0568         conbmx=b
0569        return
0570        endif
0571       else
0572        if(p.lt.eps/5.d0)then
0573         b2=b
0574        else
0575         conbmx=b
0576         return
0577        endif
0578       endif
0579 
0580       if(ntry.le.1000)goto 10
0581       write(ifmt,*)'Too much try in conbmx ... bmax=',b
0582       conbmx=b
0583       return
0584 
0585       end
0586 
0587 c-----------------------------------------------------------------------
0588       function conbmxndif()
0589 c-----------------------------------------------------------------------
0590       double precision om1intbc,p,eps
0591       include 'epos.inc'
0592       include 'epos.incsem'
0593 
0594 
0595       iomegasave=iomega
0596       iomega=2
0597       conbmxndif=0.
0598       b1=0.
0599       b2=7.
0600       conbmxndif=b2
0601       eps=1d-10
0602       p=1.d0-dexp(-om1intbc(b2))
0603       if(p.gt.2.d0*eps)goto 100
0604 
0605       ntry=0
0606 
0607 10    ntry=ntry+1
0608       b=b1+.5*(b2-b1)
0609 
0610       p=(1.d0-dexp(-om1intbc(b)))
0611 
0612       if(p.gt.eps)then
0613        if(p.gt.2.d0*eps)then
0614         b1=b
0615        else
0616         conbmxndif=b
0617        goto 100
0618        endif
0619       else
0620        if(p.lt.eps/5.d0)then
0621         b2=b
0622        else
0623         conbmxndif=b
0624         goto 100
0625        endif
0626       endif
0627 
0628       if(ntry.le.1000)goto 10
0629       write(ifmt,*)'Too much try in conbmxndif ... bkmxndif=',b
0630       conbmxndif=b
0631  100  iomega=iomegasave
0632       return
0633 
0634       end
0635 
0636 c-----------------------------------------------------------------------
0637       function conbmxdif()
0638 c-----------------------------------------------------------------------
0639 c find b to have (1-exp(-om))pmax=pdiff
0640 c-----------------------------------------------------------------------
0641       double precision om1intbc,pmax,drootom,pdiff
0642       include 'epos.inc'
0643       include 'epos.incsem'
0644 
0645       conbmxdif=0.
0646       b1=0.
0647       bmax=7.
0648       iomegasave=iomega
0649       iomega=2
0650 
0651       eps=1.e-5
0652       pmax=1.d0-dexp(-om1intbc(b1))
0653       pdiff=facdif
0654       if(pmax.gt.eps)then
0655         conbmxdif=drootom(pdiff,pmax,bmax,eps)
0656       endif
0657       iomega=iomegasave
0658 
0659       return
0660 
0661       end
0662 
0663 c-----------------------------------------------------------------------
0664       subroutine conre
0665 c-----------------------------------------------------------------------
0666 c  initializes remnants
0667 c-----------------------------------------------------------------------
0668       include "epos.incems"
0669       include 'epos.inc'
0670 
0671       call utpri('conre ',ish,ishini,6)
0672 
0673 c     proj
0674 c     ----
0675       la=laproj
0676       ma=iabs(maproj)
0677       las=0
0678       mas=0
0679       do l=1,ma
0680       if(la.lt.0)then
0681         if(iabs(idproj).lt.20)then
0682           id=idproj
0683         else
0684           ia=iabs(idproj/10)
0685           is=idproj/iabs(idproj)
0686           if(ia.ne.111.and.ia.ne.222.and.ia.ne.333)id=idproj/10*10
0687           if(ia.eq.111.or. ia.eq.222.or. ia.eq.333)id=idproj/10*10+is
0688           if(ia.eq.213)id=1230*is
0689         endif
0690       else
0691         id=1220
0692         if(rangen().le.(la-las)*1./(ma-mas))id=1120
0693         if(id.eq.1120)las=las+1
0694         mas=mas+1
0695       endif
0696       ic1=idtrai(1,id,1)
0697       ic2=idtrai(2,id,1)
0698       icproj(1,l)=ic1
0699       icproj(2,l)=ic2
0700       enddo
0701 
0702 c     targ
0703 c     ----
0704       la=latarg
0705       ma=iabs(matarg)
0706       las=0
0707       mas=0
0708       do l=1,ma
0709       if(la.lt.0)then
0710         if(iabs(idtarg).lt.20)then
0711           id=idtarg
0712         else
0713           ia=iabs(idtarg/10)
0714           is=idtarg/iabs(idtarg)
0715           if(ia.ne.111.and.ia.ne.222.and.ia.ne.333)id=idtarg/10*10
0716           if(ia.eq.111.or. ia.eq.222.or. ia.eq.333)id=idtarg/10*10+is
0717           if(ia.eq.213)id=1230*is
0718         endif
0719       else
0720         id=1220
0721         if(rangen().le.(la-las)*1./(ma-mas))id=1120
0722         if(id.eq.1120)las=las+1
0723         mas=mas+1
0724       endif
0725       ic1=idtrai(1,id,1)
0726       ic2=idtrai(2,id,1)
0727       ictarg(1,l)=ic1
0728       ictarg(2,l)=ic2
0729       enddo
0730 
0731       call utprix('conre ',ish,ishini,6)
0732       return
0733       end
0734 
0735 c-----------------------------------------------------------------------
0736       subroutine conrl
0737 c-----------------------------------------------------------------------
0738 c  initializes target remnant in case of appl lepton
0739 c-----------------------------------------------------------------------
0740       include "epos.incems"
0741       common/nucl1/laproj,maproj,latarg,matarg,core,fctrmx
0742       common/hadr2/iomodl,idproj,idtarg,wexcit
0743 
0744 c     targ
0745 c     ----
0746       la=latarg
0747       ma=iabs(matarg)
0748       las=0
0749       mas=0
0750       do l=1,ma
0751       if(la.lt.0)then
0752       id=idtarg
0753       else
0754       id=1220
0755       if(rangen().le.(la-las)*1./(ma-mas))id=1120
0756       if(id.eq.1120)las=las+1
0757       mas=mas+1
0758       endif
0759       ic1=idtrai(1,id,1)
0760       ic2=idtrai(2,id,1)
0761       ictarg(1,l)=ic1
0762       ictarg(2,l)=ic2
0763       enddo
0764 
0765       return
0766       end
0767 
0768 c-----------------------------------------------------------------------
0769       subroutine conwr
0770 c-----------------------------------------------------------------------
0771 c     writes /cptl/
0772 c-----------------------------------------------------------------------
0773       include "epos.inc"
0774       include "epos.incems"
0775       double precision XA(210,3),XB(210,3),BQGS,BMAXQGS,BMAXNEX,BMINNEX
0776       COMMON /Q_QGSNEX1/ XA,XB,BQGS,BMAXQGS,BMAXNEX,BMINNEX
0777       common/nucl3/phi,bimp
0778       common /cncl/xproj(mamx),yproj(mamx),zproj(mamx)
0779      *,xtarg(mamx),ytarg(mamx),ztarg(mamx)
0780       parameter(iapmax=208)
0781       double precision bqgs2,bmaxqgs2,bmaxnex2,bminnex2,xan,xbn
0782       common /qgsIInex1/xan(iapmax,3),xbn(iapmax,3)
0783      *,bqgs2,bmaxqgs2,bmaxnex2,bminnex2
0784       common/photrans/phoele(4),ebeam
0785       integer ic(2)
0786 
0787       call utpri('conwr ',ish,ishini,6)
0788 
0789       bx=cos(phi)*bimp
0790       by=sin(phi)*bimp
0791 
0792 c     write /cptl/
0793 c     ------------
0794       nptl=0
0795 
0796            if(iokoll.eq.1)then   ! precisely matarg collisions
0797 
0798       nptl=nptl+1
0799       do 3 i=1,4
0800 3     xorptl(i,nptl)=0
0801       tivptl(1,nptl)=-ainfin
0802       tivptl(2,nptl)=0
0803       istptl(nptl)=1
0804       iorptl(nptl)=-1
0805       jorptl(nptl)=0
0806       do 1 k=1,koll
0807       nptl=nptl+1
0808       do 4 i=1,4
0809 4     xorptl(i,nptl)=0
0810       tivptl(1,nptl)=-ainfin
0811       tivptl(2,nptl)=0
0812       istptl(nptl)=1
0813       iorptl(nptl)=-1
0814       jorptl(nptl)=0
0815 1     continue
0816 
0817            elseif(iappl.ne.7)then
0818 
0819 c             print *,'proj'
0820 
0821       do 6 i=1,maproj
0822       nptl=nptl+1
0823       istptl(nptl)=0
0824       iorptl(nptl)=0
0825       jorptl(nptl)=0
0826       if(model.eq.2)then       !QGSJet
0827       xproj(i)=XA(i,1)
0828       yproj(i)=XA(i,2)
0829       zproj(i)=XA(i,3)
0830       istptl(nptl)=1
0831       iorptl(nptl)=-1
0832       elseif(model.eq.7.or.model.eq.11)then       !QGSJetII
0833       xproj(i)=xan(i,1)
0834       yproj(i)=xan(i,2)
0835       zproj(i)=xan(i,3)
0836       istptl(nptl)=1
0837       iorptl(nptl)=-1
0838       elseif(model.ge.3)then       !Gheisha, ...
0839       istptl(nptl)=1
0840       iorptl(nptl)=-1
0841       endif
0842       xorptl(1,nptl)=xproj(i)+bx/2
0843       xorptl(2,nptl)=yproj(i)+by/2
0844       xorptl(3,nptl)=zproj(i)
0845       xorptl(4,nptl)=0
0846       tivptl(1,nptl)=-ainfin
0847 c for visualisation uncomment
0848 c-c   tivptl(1,nptl)=-100
0849       tivptl(2,nptl)= ainfin
0850 c             print *,i,xorptl(1,nptl),xorptl(2,nptl),xorptl(3,nptl)
0851 6     continue
0852 c             print *,'targ'
0853       do 7 i=1,matarg
0854       nptl=nptl+1
0855       istptl(nptl)=0
0856       iorptl(nptl)=0
0857       jorptl(nptl)=0
0858       if(model.eq.2)then       !QGSJet
0859       xtarg(i)=XB(i,1)
0860       ytarg(i)=XB(i,2)
0861       ztarg(i)=XB(i,3)
0862       istptl(nptl)=1
0863       iorptl(nptl)=-1
0864       elseif(model.eq.7.or.model.eq.11)then       !QGSJetII
0865       xtarg(i)=xbn(i,1)
0866       ytarg(i)=xbn(i,2)
0867       ztarg(i)=xbn(i,3)
0868       istptl(nptl)=1
0869       iorptl(nptl)=-1
0870       elseif(model.ge.3)then       !Gheisha, ...
0871       istptl(nptl)=1
0872       iorptl(nptl)=-1
0873       endif
0874       xorptl(1,nptl)=xtarg(i)-bx/2
0875       xorptl(2,nptl)=ytarg(i)-by/2
0876       xorptl(3,nptl)=ztarg(i)
0877       xorptl(4,nptl)=0
0878       tivptl(1,nptl)=-ainfin
0879 c for visualisation uncomment
0880 c-c   tivptl(1,nptl)=-100
0881       tivptl(2,nptl)= ainfin
0882 c             print *,i,xorptl(1,nptl),xorptl(2,nptl),xorptl(3,nptl)
0883 7     continue
0884       if(abs(idprojin).eq.12)then   !electron for fake DIS
0885 c electron projectile
0886         nptl=nptl+1
0887         istptl(nptl)=41
0888         iorptl(nptl)=-1
0889         jorptl(nptl)=-1
0890         iorptl(1)=nptl         !pi0 (porjectile) coming from lepton
0891         xorptl(1,nptl)=bx/2
0892         xorptl(2,nptl)=by/2
0893         xorptl(3,nptl)=0.
0894         xorptl(4,nptl)=0.
0895         tivptl(1,nptl)=-ainfin
0896         tivptl(2,nptl)=0.
0897 c target nucleons (in lab frame)
0898         do i=1,matarg
0899           nptl=nptl+1
0900           istptl(nptl)=41
0901           iorptl(nptl)=-1
0902           jorptl(nptl)=-1
0903           xorptl(1,nptl)=xtarg(i)-bx/2
0904           xorptl(2,nptl)=ytarg(i)-by/2
0905           xorptl(3,nptl)=ztarg(i)
0906           xorptl(4,nptl)=0
0907           tivptl(1,nptl)=-ainfin
0908 c         for visualisation uncomment
0909 c         -c   tivptl(1,nptl)=-100
0910           tivptl(2,nptl)= ainfin
0911         enddo
0912 c electron remnant
0913         nptl=nptl+1
0914         istptl(nptl)=0
0915         iorptl(nptl)=maproj+matarg+1
0916         jorptl(nptl)=-1
0917         xorptl(1,nptl)=bx/2
0918         xorptl(2,nptl)=by/2
0919         xorptl(3,nptl)=0.
0920         xorptl(4,nptl)=0.
0921         tivptl(1,nptl)=0.
0922         tivptl(2,nptl)= ainfin
0923       endif
0924 
0925           endif
0926 
0927       nptl=0
0928       if(iappl.le.2)then
0929       do i=1,maproj
0930       nptl=nptl+1
0931       ic(1)=icproj(1,i)
0932       ic(2)=icproj(2,i)
0933       id=idtra(ic,0,0,0)
0934 c      id=idtra(ic,0,0,3)      !tp071107 imix=3 ??????????
0935       call idmass(id,ams)
0936       idptl(nptl)=id
0937       pptl(1,nptl)=0.
0938       pptl(2,nptl)=0.
0939       pptl(3,nptl)=pnullx
0940       pptl(4,nptl)=sqrt(pnullx**2+ams**2)
0941       pptl(5,nptl)=ams
0942       ifrptl(1,nptl)=0
0943       ifrptl(2,nptl)=0
0944       ityptl(nptl)=0
0945       enddo
0946       endif
0947       if(iappl.ne.7)then
0948       do i=1,matarg
0949       nptl=nptl+1
0950       ic(1)=ictarg(1,i)
0951       ic(2)=ictarg(2,i)
0952       id=idtra(ic,0,0,0)
0953 c      id=idtra(ic,0,0,3)      !tp071107 imix=3 ??????????
0954       call idmass(id,ams)
0955       idptl(nptl)=id
0956       pptl(1,nptl)=0.
0957       pptl(2,nptl)=0.
0958       pptl(3,nptl)=-pnullx
0959       pptl(4,nptl)=sqrt(pnullx**2+ams**2)
0960       pptl(5,nptl)=ams
0961       ifrptl(1,nptl)=0
0962       ifrptl(2,nptl)=0
0963       ityptl(nptl)=0
0964       enddo
0965       if(abs(idprojin).eq.12)then   !electron for fake DIS
0966 c electron projectile
0967         nptl=nptl+1
0968         id=idprojin
0969         call idmass(id,ams)
0970         idptl(nptl)=id
0971         pptl(1,nptl)=0.
0972         pptl(2,nptl)=0.
0973         pptl(3,nptl)=sqrt(max(0.,(elepti+ams)*(elepti-ams)))
0974         pptl(4,nptl)=elepti
0975         pptl(5,nptl)=ams
0976         ifrptl(1,nptl)=1
0977         ifrptl(2,nptl)=1
0978         ityptl(nptl)=40
0979 c target nucleons (in lab frame)
0980         do i=1,matarg
0981           nptl=nptl+1
0982           idptl(nptl)=idptl(maproj+i)
0983           pptl(1,nptl)=0.
0984           pptl(2,nptl)=0.
0985           pptl(3,nptl)=-pnll
0986           pptl(4,nptl)=ebeam
0987           pptl(5,nptl)=pptl(5,maproj+i)
0988           ifrptl(1,nptl)=maproj+i
0989           ifrptl(2,nptl)=maproj+i
0990           ityptl(nptl)=50
0991         enddo
0992 c electron remnant
0993         nptl=nptl+1
0994         idptl(nptl)=id
0995         pptl(1,nptl)=phoele(1)
0996         pptl(2,nptl)=phoele(2)
0997         pptl(3,nptl)=phoele(3)
0998         pptl(4,nptl)=phoele(4)
0999         pptl(5,nptl)=ams
1000         ifrptl(1,nptl)=0
1001         ifrptl(2,nptl)=0
1002         ityptl(nptl)=40
1003       endif
1004 
1005       else
1006 
1007       nptl=nptl+1
1008       id=idproj
1009       call idmass(id,ams)
1010       idptl(nptl)=id
1011       pptl(1,nptl)=0.
1012       pptl(2,nptl)=0.
1013       pptl(3,nptl)=pnullx
1014       pptl(4,nptl)=sqrt(pnullx**2+ams**2)
1015       pptl(5,nptl)=ams
1016       ifrptl(1,nptl)=0
1017       ifrptl(2,nptl)=0
1018       ityptl(nptl)=0
1019       iorptl(nptl)=-1
1020       jorptl(nptl)=0
1021       istptl(nptl)=0
1022       do 5 i=1,4
1023  5      xorptl(i,nptl)=0
1024       tivptl(1,nptl)=0
1025       tivptl(2,nptl)=0
1026       endif
1027 
1028 c     exit
1029 c     ----
1030 
1031       call utprix('conwr ',ish,ishini,6)
1032       return
1033       end
1034 
1035 c------------------------------------------------------------------------
1036       subroutine conxyz(ch,n,x,y,z,ynuc)
1037 c-----------------------------------------------------------------------
1038       include 'epos.inc'
1039 
1040       real x(n),y(n),z(n)
1041       character ch*1
1042 
1043       massnr=0
1044       iii=0
1045       if(ch.eq.'p')then
1046       massnr=maproj
1047       iii=1
1048       elseif(ch.eq.'t')then
1049       massnr=matarg
1050       iii=2
1051       else
1052       call utstop('conxyz: nucleus neither proj nor targ&')
1053       endif
1054 
1055       if(massnr.eq.0)return
1056       if(massnr.gt.n)call utstop('conxyz: massnr.gt.n&')
1057       if(massnr.eq.1)then
1058       x(1)=0
1059       y(1)=0
1060       z(1)=0
1061       return
1062       endif
1063 
1064       rad=radnuc(massnr)
1065 
1066       if(massnr.ge.10)then !---wood-saxon density---
1067 
1068         rad=rad/difnuc(massnr)
1069         cr1=1.+3./rad+6./rad**2+6./rad**3
1070         cr2=3./rad
1071         cr3=3./rad+6./rad**2
1072         do i=1,massnr
1073    1      zuk=rangen()*cr1-1.
1074           if(zuk.le.0.)then
1075             tt=rad*(rangen()**.3333-1.)
1076           elseif(zuk.le.cr2 )then
1077             tt=-log(rangen())
1078           elseif(zuk.lt.cr3 )then
1079             tt=-log(rangen())-log(rangen())
1080           else
1081             tt=-log(rangen())-log(rangen())-log(rangen())
1082           endif
1083           if(rangen().gt.1./(1.+exp(-abs(tt))))goto 1
1084           rim=tt+rad
1085           zz=rim*(2.*rangen()-1.)
1086           rim=sqrt(rim*rim-zz*zz)
1087           z(i)=zz*difnuc(massnr)
1088           call pscs(c,s)
1089           x(i)=rim*c*difnuc(massnr)
1090           y(i)=rim*s*difnuc(massnr)
1091         enddo
1092 
1093       elseif(massnr.ge.3)then  ! ---gaussian density---
1094 
1095         rad=rad*sqrt(2.*massnr/(massnr-1.))   !van hove simulation
1096         do l=1,3
1097           summ=0.
1098           do i=1,massnr-1
1099             j=massnr-i
1100             aks=rad *(rangen()+rangen()+rangen()-1.5)
1101             k=j+1
1102             if(l.eq.1)x(k)=summ-aks*sqrt(float(j)/k)
1103             if(l.eq.2)y(k)=summ-aks*sqrt(float(j)/k)
1104             if(l.eq.3)z(k)=summ-aks*sqrt(float(j)/k)
1105             summ=summ+aks/sqrt(float(j*k))
1106           enddo
1107           if(l.eq.1)x(1)=summ
1108           if(l.eq.2)y(1)=summ
1109           if(l.eq.3)z(1)=summ
1110         enddo
1111 
1112       elseif(massnr.eq.2)then  ! ---deuteron---
1113 
1114         !.........r=t*difnuc(massnr), t~exp(-2*t)*(1-exp(-a*t))
1115         a=radnuc(massnr)
1116   2     t=-0.5*alog(rangen())  !~exp(-2*t)
1117         if(rangen().gt.(1-exp(-a*t))**2)goto2
1118         r=t*difnuc(massnr)
1119         zz=r*(2.*rangen()-1.)
1120         call pscs(c,s)
1121         rxy=sqrt(r*r-zz*zz)
1122         z(1)=0.5*zz
1123         x(1)=0.5*rxy*c
1124         y(1)=0.5*rxy*s
1125         z(2)=-z(1)
1126         x(2)=-x(1)
1127         y(2)=-y(1)
1128 
1129       else
1130 
1131         stop'conxyz: wrong massnr.     '
1132 
1133       endif
1134 
1135 c...plot preparation
1136 
1137       rmax=(radnuc(massnr)+3)
1138       drnucl(iii)=rmax/mxnucl
1139       nrnucl(iii)=nrnucl(iii)+1
1140       do i=1,massnr
1141         r=sqrt(x(i)**2+y(i)**2+z(i)**2)
1142         k=1+int(r/drnucl(iii))
1143         if(k.le.mxnucl)rnucl(k,iii)=rnucl(k,iii)+1
1144       enddo
1145 
1146 c...lorentz trafo
1147 
1148       do i=1,massnr
1149       z(i)=z(i)/cosh(ynuc)
1150       enddo
1151 
1152       return
1153       end
1154 
1155 c-----------------------------------------------------------------------
1156       subroutine conini
1157 c-----------------------------------------------------------------------
1158       include 'epos.inc'
1159 
1160       imax=max(maproj,matarg)
1161       if(idtargin.eq.0)imax=max(imax,40)
1162       do massnr=1,mamxx
1163         dif=0.54
1164         rad=0.
1165         if(massnr.gt.imax.or.massnr.eq.1)then
1166           dif=0
1167           rad=0
1168         elseif(massnr.eq.197)then
1169           dif=0.562
1170           rad=6.5
1171         elseif(massnr.ge.10)then
1172           rad=1.12*massnr**0.33333-0.86*massnr**(-0.33333)
1173         elseif(massnr.ge.3)then
1174           rad=.9*float(massnr)**.3333
1175         elseif(massnr.eq.2)then
1176           dif=4.316
1177           rad=4.68
1178         endif
1179         difnuc(massnr)=dif
1180         radnuc(massnr)=rad
1181       enddo
1182 
1183       end
1184 
1185 c-----------------------------------------------------------------------
1186       subroutine xConNuclDens(iii)
1187 c-----------------------------------------------------------------------
1188 c plots distribution of nucleons in nuclei
1189 c  iii = 1 (proj) or 2 (targ)
1190 c-----------------------------------------------------------------------
1191       include 'epos.inc'
1192       if(model.ne.1)return
1193       massnr=1
1194       if(iii.eq.1)then
1195         massnr=maproj
1196       elseif(iii.eq.2)then
1197         massnr=matarg
1198       endif
1199       if(massnr.eq.1)return
1200       a=1./4.316
1201       b=4.68
1202       write(ifhi,'(a)') '!-----------------------------------------'
1203       write(ifhi,'(a)') '!          nuclear density          '
1204       write(ifhi,'(a)') '!-----------------------------------------'
1205       write(ifhi,'(a)')       'openhisto'
1206       if(massnr.ge.10)write(ifhi,'(a)')'htyp lin xmod lin ymod lin'
1207       if(massnr.lt.10)write(ifhi,'(a)')'htyp lin xmod lin ymod log'
1208       write(ifhi,'(a)')       'text 0 0 "title nuclear density"'
1209       write(ifhi,'(a)')       'text 0.99 -0.15 " r (fm) "'
1210       write(ifhi,'(a)')       'text 0 0 "yaxis rho(r)"'
1211       write(ifhi,'(a,2e11.3)')'xrange',0.,mxnucl*drnucl(iii)
1212       write(ifhi,'(3a)')'yrange',' 0 ',' auto'
1213       write(ifhi,'(a)')       'array 2'
1214       do j=1,mxnucl
1215       r=(j-0.5)*drnucl(iii)
1216       d=0.5*drnucl(iii)
1217       write(ifhi,'(2e12.4)')  r,rnucl(j,iii)/nrnucl(iii)/
1218      *                     (4./3.*pi*((r+d)**3-(r-d)**3))
1219       enddo
1220       write(ifhi,'(a)')       '  endarray'
1221       write(ifhi,'(a)')       'closehisto plot 0-'
1222       write(ifhi,'(a)')       'openhisto'
1223       write(ifhi,'(a)')       'htyp lbo '
1224       write(ifhi,'(a)')       'array 2'
1225       do j=1,mxnucl
1226       r=(j-0.5)*drnucl(iii)
1227       rr=2*r
1228       rho=1
1229       if(massnr.eq.2)then
1230         rho=1.00*((1-exp(-b*a*rr))*exp(-a*rr)/rr)**2
1231       elseif(massnr.eq.197)then
1232         rho=0.16/(1+exp((r-6.5)/0.562))
1233       elseif(massnr.ge.10)then
1234         rho=0.16/(1+exp((r-radnuc(massnr))/difnuc(massnr)))
1235       endif
1236       write(ifhi,'(2e12.4)')  r,rho
1237       enddo
1238       write(ifhi,'(a)')       '  endarray'
1239       write(ifhi,'(a)')       'closehisto plot 0'
1240       end
1241 
1242 c-----------------------------------------------------------------------
1243       subroutine xConThick(iii)
1244 c-----------------------------------------------------------------------
1245       ! plots sigma_pp *T_A (b)  (=average number of collisions)
1246       ! T_A = thickness function
1247       !  iii = 1 (proj) or 2 (targ)
1248       !----------------------------------------------------------------
1249       include 'epos.inc'
1250       parameter(iconimax=20,iconkmax=100)
1251       real thick(2,0:iconimax)
1252       imax=iconimax
1253       kmax=iconkmax
1254       if(model.ne.1)return
1255       ramx=mxnucl*drnucl(iii)
1256       do i=0,imax
1257         x=i/float(imax)*ramx
1258         sum=0
1259         rho0=conrho(iii,0.)
1260         h=ramx/kmax
1261         do k=1,kmax
1262           z=k*h
1263           r=sqrt(x*x+z*z)
1264           rho2=conrho(iii,r)
1265           z=(k-0.5)*h
1266           r=sqrt(x*x+z*z)
1267           rho1=conrho(iii,r)
1268           sum=sum+h/6.*(rho0+4*rho1+rho2)
1269           rho0=rho2
1270         enddo
1271         sum=sum*2 ! integral fro -infty to +infty
1272         thick(1,i)=x
1273         thick(2,i)=sum
1274       enddo
1275       write(ifhi,'(a)')       'openhisto'
1276       write(ifhi,'(a)')       'htyp lbo '
1277       write(ifhi,'(a)')       'txt "xaxis b (fm)" '
1278       write(ifhi,'(a)')       'txt "yaxis [s]?pp! T?A! (b) " '
1279       write(ifhi,'(a)')       'array 2'
1280       do i=0,imax
1281         write(ifhi,'(2e12.4)') thick(1,i),sigine/10.*thick(2,i)
1282       enddo
1283       write(ifhi,'(a)')       '  endarray'
1284       write(ifhi,'(a)')       'closehisto plot 0'
1285 
1286       end
1287 
1288 c-----------------------------------------------------------------------
1289       function conrho(iii,r)
1290 c-----------------------------------------------------------------------
1291       ! nuclear density
1292       !  iii = 1 (proj) or 2 (targ)
1293       !----------------------------------------------------------------
1294       include 'epos.inc'
1295       conrho=1.
1296       if(model.ne.1)return
1297       massnr=1
1298       if(iii.eq.1)then
1299         massnr=maproj
1300       elseif(iii.eq.2)then
1301         massnr=matarg
1302       endif
1303       if(massnr.eq.1)return
1304       a=1./4.316
1305       b=4.68
1306       rr=2*r
1307       rho=1
1308       if(massnr.eq.2.and.rr.gt.0.)then
1309         rho=1.00*((1-exp(-b*a*rr))*exp(-a*rr)/rr)**2
1310       elseif(massnr.eq.197)then
1311         rho=0.16/(1+exp((r-6.5)/0.562))
1312       elseif(massnr.ge.10)then
1313         rho=0.16/(1+exp((r-radnuc(massnr))/difnuc(massnr)))
1314       endif
1315       conrho=rho
1316       end
1317 
1318 
1319 
1320 
1321 
1322 
1323