Back to home page

Project CMSSW displayed by LXR

 
 

    


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

0001 c----------------------------------------------------------------------
0002       subroutine ahydro
0003 c----------------------------------------------------------------------
0004       include 'epos.inc'
0005       common/geom/rmproj,rmtarg,bmax,bkmx
0006       common/cranphi/ranphi
0007       ntry=0
0008  1    ntry=ntry+1
0009       if(ntry.eq.100)stop'in ahydro: infine loop (080719)   '
0010       nevt=0
0011       nptl=0
0012       rmproj=1.19*maproj**(1./3.)-1.61*maproj**(-1./3.)+fctrmx*.54
0013       rmtarg=1.19*matarg**(1./3.)-1.61*matarg**(-1./3.)+fctrmx*.54
0014       b1=bminim
0015       b2=min(rmproj+rmtarg,bmaxim)
0016       bimevt=sqrt(b1**2+(b2**2-b1**2)*rangen())
0017       phievt=0
0018       ranphi=0
0019       call InitializeHyperbola !does not affect results but necessary
0020       call HydroFO(ier)
0021       if(ier.eq.1)goto1
0022       do n=1,nptl
0023         iorptl(n)=0
0024         jorptl(n)=0
0025         istptl(n)=0
0026         ifrptl(1,n)=0
0027         ifrptl(2,n)=0
0028         tivptl(1,n)=xorptl(4,n)
0029         call idtau(idptl(n),pptl(4,n),pptl(5,n),taugm)
0030         r=rangen()
0031         tivptl(2,n)=tivptl(1,n)+taugm*(-alog(r))
0032         radptl(n)=0.
0033         dezptl(n)=0.
0034         itsptl(n)=0
0035         rinptl(n)=0
0036       enddo
0037       end
0038 
0039 c----------------------------------------------------------------------
0040       subroutine amicro(iret)
0041 c----------------------------------------------------------------------
0042 c  microcanonical decay of cluster specified via keu...ket, tecm, volu
0043 c----------------------------------------------------------------------
0044       include 'epos.inc'
0045       parameter(maxp=500)
0046       common/confg/np,amass(maxp),ident(maxp),pcm(5,maxp),wtxlog,wtlog
0047       double precision seedp
0048       data ncntmi/0/
0049       save ncntmi
0050       call utpri('amicro',ish,ishini,4)
0051       ncntmi=ncntmi+1
0052 
0053       if(ncntmi.eq.1)then
0054         call ranfgt(seedp)      !not to change the seed ...
0055         if(hydt.ne.'---')then
0056           call HydroTable2(0)
0057           call DefineParticles
0058         else
0059           call ManiParticles(93,0)
0060         endif
0061         call ranfst(seedp)      ! ... after this initialization
0062       endif
0063 
0064       iret=0
0065       nevt=0
0066       nptl=0
0067 c      dez=0.5e-4
0068       call InitializeHyperbola !does not affect results but necessary
0069 
0070 c  50  continue
0071       call GraCan
0072       energ=0
0073       do i=1,np
0074         energ=energ+pcm(4,i)
0075       enddo
0076       !print*,'+++++',energ,tecm,tecm/volu
0077 ccc      if(abs(energ-tecm).gt.0.1) goto50   !uncomment for energy conservation
0078 
0079       do n=1,np
0080         nptl=nptl+1
0081         if(nptl.gt.mxptl)call utstop('StaHadShort: mxptl too small&')
0082         idptl(nptl)=ident(n)
0083         do j=1,4
0084           pptl(j,nptl)=pcm(j,n)
0085           xorptl(j,nptl)=0
0086         enddo
0087         pptl(5,nptl)=amass(n)
0088         ityptl(nptl)=19
0089       enddo
0090 
0091       do n=1,nptl
0092         iorptl(n)=0
0093         jorptl(n)=0
0094         istptl(n)=0
0095         ifrptl(1,n)=0
0096         ifrptl(2,n)=0
0097         tivptl(1,n)=0
0098         call idtau(idptl(n),pptl(4,n),pptl(5,n),taugm)
0099         r=rangen()
0100         tivptl(2,n)=tivptl(1,n)+taugm*(-alog(r))
0101         radptl(n)=0.
0102         itsptl(n)=0
0103         rinptl(n)=-9999
0104       enddo
0105 
0106       call utprix('amicro',ish,ishini,4)
0107       return
0108       end
0109 
0110 c-----------------------------------------------------------------------
0111       subroutine GraCan
0112 c-----------------------------------------------------------------------
0113       ! decays cluster specified via tfo and volu using
0114       !    hadron gas (ioclude=4) or recombination (ioclude=5)
0115       ! pcm(5,n) is the usual energy (ioclude=4)
0116       ! or the effective energy, based on quark mass sum (ioclude=5)
0117       !        so for the flow boost use this energy
0118       !----------------------------------------------------------------
0119       include 'epos.inc'
0120       include 'epos.inchy'
0121       common/copt/istat
0122       parameter(maxp=500)
0123       common/confg/np,amass(maxp),ident(maxp),pcm(5,maxp),wtxlog,wtlog
0124       parameter (mspez=54)
0125       common/cspez1/nspez,ispez(mspez),aspez(2,mspez),gspez(mspez)
0126       common/cspez2/kspez,jspez(mspez)
0127       common/cspez3/fuga(mspez)
0128       common/cspez4/ffstat(2,0:mspez+2) /ctfo/tfo
0129       real u(3)
0130       io3=ioclude-3
0131       if(io3.lt.1.or.io3.gt.2)stop'in GraCan: wrong ioclude (140808) '
0132       yie= volu * ffstat(io3,nspez) * ffstat(io3,mspez+2)
0133       np=yie
0134       if(rangen().le.yie-np)np=np+1
0135       do n=1,np
0136         r=rangen()*yie
0137         i=0
0138         do while(yie*ffstat(io3,i)/ffstat(io3,nspez).lt.r)
0139          i=i+1
0140         enddo
0141         ident(n)=ispez(i)
0142         amass(n)=aspez(1,i)
0143         fug=fuga(i)
0144         js=jspez(i)
0145         x=RanBoseFermi(aspez(io3,i),fug,js,tfo,istat)
0146         p=x*tfo
0147         e=sqrt(amass(n)**2+p**2)
0148         ex=sqrt(aspez(io3,i)**2+p**2)
0149         u(3)=2.*rangen()-1.
0150         phi=2.*pi*rangen()
0151         u(1)=sqrt(1.-u(3)**2)*cos(phi)
0152         u(2)=sqrt(1.-u(3)**2)*sin(phi)
0153         pcm(1,n)=p*u(1)
0154         pcm(2,n)=p*u(2)
0155         pcm(3,n)=p*u(3)
0156         pcm(4,n)=e
0157         pcm(5,n)=ex
0158       enddo
0159       end
0160 
0161 c----------------------------------------------------------------------
0162       subroutine HydroTable2(ichk)
0163 c----------------------------------------------------------------------
0164       include 'epos.inc'
0165       include 'epos.inchy'
0166       parameter (mspez=54)
0167       common/cspez1/nspez,ispez(mspez),aspez(2,mspez),gspez(mspez)
0168       common/cspez2/kspez,jspez(mspez) /cstep/netastep /crapi/delrapi
0169       common/ctempcrit/tempcrit,epscrit
0170       character tabname*550
0171       
0172       if(hydt.eq.'---')stop'in HydroTable2. no table found (150808) '
0173       write(tabname,'(a)')fnnx(1:nfnnx)//'epos.ini'//hydt//' '
0174       call HydroTable(ichk,tabname)
0175       epscri(3)=epscrit
0176       kspez=1
0177       nspez=54
0178       end
0179 
0180 c------------------------------------------------------------------------------
0181       subroutine HydroTable(ichk,tabname) ! 13 Aug 08 (same as gyx.f except inc)
0182 c------------------------------------------------------------------------------
0183       include 'epos.inc'
0184       include 'epos.inchy'
0185       character txtdum*40, tabname*550
0186       common/ctempcrit/tempcrit,epscrit/ctfo/tfo
0187       nchar=index(tabname,' ')-1
0188       write(*,'(a,$)')
0189      * ' reading table '//tabname(1:nchar)//' ...'
0190       open(unit=3,file=tabname(1:nchar),status='old'
0191      * ,err=99)
0192       read(3,'(a)')txtdum
0193       read(3,*)maprojx,matargx,engyx,epscrit,tempcrit
0194       tfo=tempcrit
0195       if(ichk.eq.1)then
0196        if(maprojx.ne.maproj)stop'HydroTable: maprojx.ne.maproj.       '
0197        if(matargx.ne.matarg)stop'HydroTable: matargx.ne.matarg.       '
0198        if(engyx.ne.engy)stop'HydroTable: engyx.ne.engy.       '
0199       endif
0200       read(3,*)ncenthy,netahy,ntauhy,nphihy,nradhy
0201       if(ncenthx.lt.ncenthy)stop'HydroTable: ncenthx too small.   '
0202       if(netahx.lt.netahy)  stop'HydroTable: netahx too small.   '
0203       if(ntauhx.lt.ntauhy)  stop'HydroTable: ntauhx too small.   '
0204       if(nphihx.lt.nphihy)  stop'HydroTable: nphihx too small.   '
0205       if(nradhx.lt.nradhy)  stop'HydroTable: nradhx too small.   '
0206       read(3,*)(ntauhoc(ncent),ncent=1,ncenthy)
0207       read(3,*)(centhy(ncent),ncent=1,ncenthy)
0208      *         ,(etahy(neta),  neta=  1,netahy)
0209      *         ,(phihy(nphi),  nphi=  1,nphihy)
0210      *         ,(radhy(nrad),  nrad=  1,nradhy)
0211       read(3,*)((tauhoc(ncent,ntau),ncent=1,ncenthy), ntau=1,ntauhy)
0212       read(3,*)((((epsii(ncent,neta,nphi,nrad),ncent=1,ncenthy)
0213      *          ,neta=1,netahy),nphi=1,nphihy),nrad=1,nradhy)
0214       read(3,*)meos,meosmu
0215       read(3,*)((eos(i,m),i=1,3),m=1,meos)
0216      *          ,((eosmu(i,m),i=1,17),m=1,meosmu)
0217       read(3,*)mcenty
0218       read(3,*)(bcenty(m),m=1,mcenty),(pcenty(m),m=1,mcenty)
0219       read(3,*)nbimp
0220       read(3,*)(bimpar(1,m),m=1,nbimp),(bimpar(2,m),m=1,nbimp)
0221       read(3,*)((((rar(ncent,neta,ntau,nphi),ncent=1,ncenthy)
0222      *          ,neta=1,netahy),ntau=1,ntauhy),nphi=1,nphihy)
0223       read(3,*)(((((var(i,ncent,neta,ntau,nphi),i=1,3),ncent=1,ncenthy)
0224      *           ,neta=1,netahy),ntau=1,ntauhy),nphi=1,nphihy)
0225       close(3)
0226       do ntau=1,ntauhx
0227         tauhy(ntau)=tauhoc(1,ntau)
0228       enddo
0229       do ncent=2,ncenthy
0230         do ntau=1,ntauhoc(ncent)
0231           if(abs(tauhoc(ncent,ntau)-tauhy(ntau)).gt.1e-4)
0232      *     stop'in HydroTable: different tau grids.       '
0233         enddo
0234       enddo
0235       print*,' done'
0236       return
0237    99 print*,'HydroTable: error opening hydro table'
0238       print*,'      file=',tabname(1:nchar)
0239       stop'070817'
0240       end
0241 
0242 c----------------------------------------------------------------------------------
0243       subroutine SymmTab(ncent1,ncent2)   ! 13 Aug 08 (same as gyx.f except inc)
0244 c----------------------------------------------------------------------------------
0245       include 'epos.inchy'
0246       common/cderivs/
0247      .  ddrdtau(ncenthx,-netahx+1:netahx-1,ntauhx,nphihx)
0248      . ,ddrdphi(ncenthx,-netahx+1:netahx-1,ntauhx,nphihx)
0249      . ,ddrdeta(ncenthx,-netahx+1:netahx-1,ntauhx,nphihx)
0250       common/cgmma/gaa(ncenthx,-netahx+1:netahx-1,ntauhx,nphihx)
0251 
0252                        !..............................................
0253       write(*,'(a,$)')' making symmetric tables ...'
0254                        !''''''''''''''''''''''''''''''''''''''''''''''
0255 
0256       do ncent=ncent1,ncent2
0257        do meta=1-netahy,netahy-1
0258         do ntau=1,ntauhx
0259          do nphi=1,nphihy
0260            waa(ncent,meta,ntau,nphi)=0
0261            raa(ncent,meta,ntau,nphi)=0
0262            vaa(1,ncent,meta,ntau,nphi)=0
0263            vaa(2,ncent,meta,ntau,nphi)=0
0264            vaa(3,ncent,meta,ntau,nphi)=0
0265            ddrdphi(ncent,meta,ntau,nphi)=0
0266            ddrdtau(ncent,meta,ntau,nphi)=0
0267            ddrdeta(ncent,meta,ntau,nphi)=0
0268            gaa(ncent,meta,ntau,nphi)=0
0269          enddo
0270         enddo
0271        enddo
0272       enddo
0273 
0274       do ncent=ncent1,ncent2
0275       ntauho=ntauhoc(ncent)
0276 
0277       do neta=1,netahy
0278        meta=neta-1
0279        zetahy(meta)=etahy(neta)
0280        zetahy(-meta)=-etahy(neta)
0281       enddo
0282       if(mod(nphihy,2).ne.1)stop'in Symm (2103200820)   '
0283       do neta=1,netahy
0284         meta=neta-1
0285         do ntau=1,ntauho
0286           do nphi=1,nphihy
0287            !-------------------------
0288            !  symmetry: eta -> -eta ;  x -> -x ;  y -> y
0289            !    or   eta -> -eta ;  r -> r ;  phi -> pi-phi
0290            !    vtg -> -vtg
0291            !    v3 -> -v3
0292            !    drdphi -> -drdphi
0293            !    drdeta -> -drdeta
0294            !-------------------------
0295            nphimed=nphihy/2+1
0296            if(nphi.le.nphimed)then
0297             jphi=nphimed+1-nphi
0298            else
0299             jphi=nphihy-(nphi-nphimed)
0300            endif
0301            if(meta.ne.0)then
0302            raa(  ncent, meta,ntau,nphi)= rar(  ncent,neta,ntau,nphi)
0303            vaa(1,ncent, meta,ntau,nphi)= var(1,ncent,neta,ntau,nphi)
0304            vaa(2,ncent, meta,ntau,nphi)= var(2,ncent,neta,ntau,nphi)
0305            vaa(3,ncent, meta,ntau,nphi)= var(3,ncent,neta,ntau,nphi)
0306            raa(  ncent,-meta,ntau,nphi)= rar(  ncent,neta,ntau,jphi)
0307            vaa(1,ncent,-meta,ntau,nphi)= var(1,ncent,neta,ntau,jphi)
0308            vaa(2,ncent,-meta,ntau,nphi)=-var(2,ncent,neta,ntau,jphi)
0309            vaa(3,ncent,-meta,ntau,nphi)=-var(3,ncent,neta,ntau,jphi)
0310            else
0311            raa(ncent, meta,ntau,nphi)=
0312      .     0.5*(rar(ncent,neta,ntau,nphi)+rar(ncent,neta,ntau,jphi))
0313            vaa(1,ncent,meta,ntau,nphi)=
0314      .     0.5*(var(1,ncent,neta,ntau,nphi)+var(1,ncent,neta,ntau,jphi))
0315            vaa(2,ncent, meta,ntau,nphi)=
0316      .     0.5*(var(2,ncent,neta,ntau,nphi)-var(2,ncent,neta,ntau,jphi))
0317            vaa(3,ncent, meta,ntau,nphi)=
0318      .     0.5*(var(3,ncent,neta,ntau,nphi)-var(3,ncent,neta,ntau,jphi))
0319            endif
0320           enddo
0321         enddo
0322       enddo
0323 
0324       do meta=1-netahy,netahy-1
0325         do nphi=1,nphihy
0326             n=ntauho
0327             do while(raa(ncent,meta,n,nphi).eq.0.0.and.n.gt.2)
0328               n=n-1
0329             enddo
0330             n=n+1
0331             n=min(n,ntauho)
0332             ntauhec(ncent,meta,nphi)=n
0333         enddo
0334       enddo
0335 
0336       dphi=phihy(2)-phihy(1)
0337       dtau=tauhy(2)-tauhy(1)
0338       deta=etahy(2)-etahy(1)
0339       do meta=1-netahy,netahy-1
0340         mem=meta-1
0341         mem=max(mem,1-netahy)
0342         mep=meta+1
0343         mep=min(mep,netahy-1)
0344         do nphi=1,nphihy
0345           npp=nphi+1
0346           npp=min(npp,nphihy)
0347           npm=nphi-1
0348           npm=max(npm,1)
0349           do ntau=1,ntauhec(ncent,meta,nphi)
0350             ntm=ntau-1
0351             ntm=max(ntm,1)
0352             ntp=ntau+1
0353             ntp=min(ntp,ntauhec(ncent,meta,nphi))
0354             ddrdphi(ncent,meta,ntau,nphi)=(raa(ncent,meta,ntau,npp )
0355      .              -raa(ncent,meta,ntau,npm )) / ((npp-npm)*dphi)
0356             ddrdtau(ncent,meta,ntau,nphi)=(raa(ncent,meta,ntp ,nphi)
0357      .              -raa(ncent,meta,ntm ,nphi)) / ((ntp-ntm)*dtau)
0358             ddrdeta(ncent,meta,ntau,nphi)=(raa(ncent,mep ,ntau,nphi)
0359      .              -raa(ncent,mem ,ntau,nphi)) / ((mep-mem)*deta)
0360           enddo
0361         enddo
0362       enddo
0363 
0364       do meta=1-netahy,netahy-1
0365         do nphi=1,nphihy
0366           do ntau=1,ntauhec(ncent,meta,nphi)
0367             vv=sqrt(vaa(1,ncent,meta,ntau,nphi)**2
0368      .             +vaa(2,ncent, meta,ntau,nphi)**2
0369      .             +vaa(3,ncent, meta,ntau,nphi)**2)
0370             gaa(ncent,meta,ntau,nphi)=1./sqrt((1-vv)*(1+vv))
0371           enddo
0372         enddo
0373       enddo
0374 
0375       do meta=1-netahy,netahy-1
0376        do nphi=1,nphihy
0377 c        phi=phihy(nphi)
0378         do ntau=1,ntauhec(ncent,meta,nphi)
0379           tau=tauhy(ntau)
0380           rad=raa(ncent,meta,ntau,nphi)
0381           vrd=vaa(1,ncent,meta,ntau,nphi)
0382           vtg=vaa(2,ncent,meta,ntau,nphi)
0383           v3=vaa(3,ncent,meta,ntau,nphi)
0384           vv=sqrt(vrd**2+vtg**2+v3**2)
0385           gm=1./sqrt((1-vv)*(1+vv))
0386           dVs=      - ddrdtau(ncent,meta,ntau,nphi) * tau * rad *gm
0387      .              +                    rad * tau       *gm*vrd
0388      .              + ddrdphi(ncent,meta,ntau,nphi) * tau       *gm*vtg
0389      .              - ddrdeta(ncent,meta,ntau,nphi) * rad       *gm*v3
0390           waa(ncent,meta,ntau,nphi)=abs(dVs)
0391         enddo !tau
0392        enddo !phi
0393       enddo !neta
0394 
0395       enddo !ncent
0396 
0397       print*,' done'
0398 
0399       end
0400 
0401 c--------------------------------------------------------------------------------------
0402       subroutine DefineParticles   ! 14 Aug 08 (same as gyx.f except inc)
0403 c--------------------------------------------------------------------------------------
0404       include 'epos.inc'
0405       include 'epos.inchy'
0406       common/copt/istat
0407       parameter (mspez=54)
0408       common/cspez1/nspez,ispez(mspez),aspez(2,mspez),gspez(mspez)
0409       common/cspez2/kspez,jspez(mspez)
0410       common/cspez3/fuga(mspez)
0411       common/cspez4/ffstat(2,0:mspez+2) /ctfo/tfo
0412       parameter (nlag=15)
0413       real xlag(nlag),wlag(nlag)
0414 c      parameter (klax=5)
0415       common/cspez7/klass(mspez)
0416       data klass/
0417      *  0, 1, 1, 2, 2, 4*0
0418      *, 9*0
0419      *, 8*0
0420      *, 8*0
0421      *,10*0
0422      *,10*0/
0423       data ispez/
0424      *   110,  120, -120,  130, -130,  230, -230,  220,  330
0425      *,  111,  121, -121,  131, -131,  231, -231,  221,  331
0426      *, 1120,-1120, 1220,-1220, 1130,-1130, 2130,-2130
0427      *, 1230,-1230, 2230,-2230, 1330,-1330, 2330,-2330
0428      *, 1111,-1111, 1121,-1121, 1221,-1221, 2221,-2221, 1131,-1131
0429      *, 1231,-1231, 2231,-2231, 1331,-1331, 2331,-2331, 3331,-3331 /
0430       data jspez/
0431      *  -1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1 , 36*1/
0432       data gspez/
0433      *  9*1.
0434      *, 9*3.
0435      *, 8*2.
0436      *, 8*2.
0437      *,10*4.
0438      *,10*4./
0439       data (aspez(1,m),m=1,mspez)/
0440      *    1.349766E-01   !pi               0
0441      *,2* 1.3957018E-01  !pi               +
0442      *,2* 4.93677E-01    !K                +
0443      *,2* 4.97648E-01    !K                0
0444      *,   5.4751E-01     !eta              0
0445      *,   9.5778E-01     !eta'(958)        0
0446      *,3* 7.755E-01      !rho(770)         0
0447      *,2* 8.9166E-01     !K*(892)          +
0448      *,2* 8.9600E-01     !K*(892)          0
0449      *,   7.8265E-01     !omega(782)       0
0450      *,   1.019460E+00   !phi(1020)        0
0451      *,2* 9.3827203E-01  !p                +
0452      *,2* 9.3956536E-01  !n                0
0453      *,2* 1.18937E+00    !Sigma            +
0454      *,2* 1.115683E+00   !Lambda           0
0455      *,2* 1.192642E+00   !Sigma            0
0456      *,2* 1.197449E+00   !Sigma            -
0457      *,2* 1.31483E+00    !Xi               0
0458      *,2* 1.32131E+00    !Xi               -
0459      *,8* 1.2320E+00     !Delta(1232)      -
0460      *,2* 1.3828E+00     !Sigma(1385)      +
0461      *,2* 1.3837E+00     !Sigma(1385)      0
0462      *,2* 1.3872E+00     !Sigma(1385)      -
0463      *,2* 1.53180E+00    !Xi(1530)         0
0464      *,2* 1.5350E+00     !Xi(1530)         -
0465      *,2* 1.67245E+00    !Omega            -
0466      */
0467 c      data (aspez(1,m),m=1,mspez)/
0468 c     *     0.13496, 2* 0.13957, 2* 0.49367, 2* 0.49767, 0.54880,0.95760
0469 c     *, 3* 0.77000, 2* 0.88810, 2* 0.89220,    0.78260, 1.01960
0470 c     *, 2* 0.93828, 2* 0.93957, 2* 1.18940, 2* 1.11560
0471 c     *, 2* 1.19250, 2* 1.19740, 2* 1.31490, 2* 1.32130
0472 c     *, 8* 1.23200, 2* 1.38230, 2* 1.38200, 2* 1.38750
0473 c     *, 2* 1.53180, 2* 1.53500, 2* 1.67220   /
0474       parameter (nquark=3)
0475       real amaq(0:nquark)
0476       data amaq /0., 0.337, 0.337, 0.486/
0477       do m=1,mspez
0478         id=ispez(m)
0479         call idflav(id,i1,i2,i3,jspin,index)
0480         amx=amaq(abs(i1))+amaq(abs(i2))+amaq(abs(i3))
0481         aspez(2,m)=amx
0482       enddo
0483 
0484       write(*,'(a,$)')' DefineParticles  ...'
0485 
0486       istat=1  !0=Boltzmann, 1=Bose/Fermi
0487       pi=3.1415927
0488       hbar=0.197327
0489 
0490       call gaulag(xlag,wlag,nlag,0.)
0491       do n=1,nlag
0492       wlag(n)=wlag(n)*exp(xlag(n))
0493       enddo
0494 
0495       !write(*,'(4x,a,$)')'fuga:'
0496       do m=kspez,nspez
0497         id=ispez(m)
0498         chem=0
0499         if(meosmu.gt.0)then
0500          ihi=idxHiranoTable(id)
0501          if(ihi.gt.0)then
0502           k=1
0503           do while(eosmu(1,k).gt.tfo)
0504            k=k+1
0505           enddo
0506           f=(tfo-eosmu(1,k))/(eosmu(1,k-1)-eosmu(1,k))
0507           chem=eosmu(ihi,k)*(1-f) + eosmu(ihi,k-1)*f
0508          endif
0509         endif
0510         fuga(m)=exp(chem/tfo)
0511         !if(m.lt.6)write(*,'(i5,a,f6.2,$)')id,':',fuga(m)
0512       enddo
0513 
0514       !write(*,'(/25x,a,$)')'yie:'
0515       facphase=1./(2*pi*hbar)**3
0516       do m=0,kspez-1
0517         ffstat(1,m)=0
0518         ffstat(2,m)=0
0519       enddo
0520        eesum=0
0521        hhsum=0
0522        do m=kspez,nspez
0523         id=ispez(m)
0524         am=aspez(1,m)
0525         amx=aspez(2,m)
0526         esum=0
0527         fsum=0
0528         gsum=0
0529         hsum=0
0530         do n=1,nlag
0531           x=xlag(n)  ! p/tfo
0532           e=sqrt(am**2+x**2*tfo**2)
0533           w=exp(-sqrt(am**2/tfo**2+x**2))  * fuga(m)
0534           fsum=fsum+wlag(n)*x**2*w /(1+istat*jspez(m)*w)
0535           esum=esum+wlag(n)*x**2*w /(1+istat*jspez(m)*w) *e
0536           wx=exp(-sqrt(amx**2/tfo**2+x**2)) !???? * fuga(m)
0537           gsum=gsum+wlag(n)*x**2*wx /(1+istat*jspez(m)*wx)
0538           hsum=hsum+wlag(n)*x**2*wx /(1+istat*jspez(m)*wx) *e
0539         enddo
0540         esum=esum * facphase * gspez(m) * 4 * pi * tfo**2  * tfo
0541         fsum=fsum * facphase * gspez(m) * 4 * pi * tfo**2  * tfo
0542         gsum=gsum * facphase * gspez(m) * 4 * pi * tfo**2  * tfo
0543         hsum=hsum * facphase * gspez(m) * 4 * pi * tfo**2  * tfo
0544         ffstat(1,m)=ffstat(1,m-1)+fsum
0545         ffstat(2,m)=ffstat(2,m-1)+gsum
0546         eesum=eesum+esum
0547         hhsum=hhsum+hsum
0548         !if(m.lt.6)write(*,'(i5,a,f8.5,$)')id,':',fsum
0549       enddo
0550       ffstat(1,mspez+1)=eesum
0551       ffstat(1,mspez+2)=1
0552       ffstat(2,mspez+1)=hhsum
0553       ffstat(2,mspez+2)=eesum/hhsum
0554 
0555       print*,'  E/V=',eesum,'  Ereco/V=',hhsum,'    done'
0556       !do n=1,nlag
0557       !print*,n,xlag(n),wlag(n)
0558       !enddo
0559       end
0560 
0561 c--------------------------------------------------------------------------------------
0562       subroutine RestFrameFO(nsimu,bimp
0563      . ,jcent,ier,iprint,ianl,isto) ! 14 Aug 08 (same as gyx.f except inc)
0564 c--------------------------------------------------------------------------------------
0565       include 'epos.inc'
0566       include 'epos.inchy'
0567       common/copt/istat
0568       parameter (mspez=54,klax=5)
0569       common/cspez1/nspez,ispez(mspez),aspez(2,mspez),gspez(mspez)
0570       common/cspez2/kspez,jspez(mspez)  /cho/netaho
0571       common/cpro/probdi(20,100),dprob,mxprob
0572       common/crap/nraphi,jfac  /cstep/netastep
0573       common/cspez3/fuga(mspez)
0574       common/cspez4/ffstat(2,0:mspez+2)  /ctfo/tfo
0575       real u(3),q(4),waii(2),wbii(2),wcii(2)
0576       real wrii(2),wwii(2),wxii(2),wyii(2),wzii(2)
0577       parameter(numiv=100,rapmax=5,rapmin=-rapmax,ptvmax=2)
0578       common/cana1/rapeta(klax,-netahx+1:netahx-1,numiv)
0579       common/cana1b/rapar(klax,numiv),v2rapar(klax,numiv)
0580       common/cana1c/phaar(klax,numiv,numiv)
0581       common/cana1d/sapar(klax,numiv),v2sapar(klax,numiv)
0582       common/cana1e/ptvar(klax,numiv),ptwar(klax,numiv)
0583       common/cspez7/klass(mspez)
0584       common/cderivs/
0585      .  ddrdtau(ncenthx,-netahx+1:netahx-1,ntauhx,nphihx)
0586      . ,ddrdphi(ncenthx,-netahx+1:netahx-1,ntauhx,nphihx)
0587      . ,ddrdeta(ncenthx,-netahx+1:netahx-1,ntauhx,nphihx)
0588       ier=1
0589       if(bimp.gt.centhy(ncenthy)+0.5)return
0590       ier=0
0591       pi=3.1415927
0592 c      hbar=0.197327
0593       delrax= 2*rapmax   / numiv
0594       delptv= ptvmax   / numiv
0595       do i=1,100
0596        probdi(jcent,i)=0
0597       enddo
0598       mxprob=jcent
0599       dprob=0.1
0600       ncent=2
0601       do while(ncent.lt.ncenthy.and.centhy(ncent).lt.bimp)
0602         ncent=ncent+1
0603       enddo
0604       n1=ncent-1
0605       n2=ncent
0606       frac=(bimp-centhy(n1))/(centhy(n2)-centhy(n1))
0607       g1=1-frac
0608       g2=frac
0609       ntauho=max(ntauhoc(n1),ntauhoc(n2))
0610 
0611       if(iprint.eq.1)then
0612        print*,'++++++++',ncent, bimp, centhy(n1),centhy(n2),g1,g2
0613        print*,'++++++++',ntauho,netaho,jfac*netastep,nsimu
0614        write(*,'(a,$)')' Rest Frame freeze out ...    '
0615       endif
0616 
0617       dphi=phihy(2)-phihy(1)
0618       dtau=tauhy(2)-tauhy(1)
0619 c      deta=etahy(2)-etahy(1)
0620       dleta=(etahy(2)-etahy(1))/jfac
0621       dall=dphi*dtau*dleta*0.125
0622       do nphi=1,nphihy
0623       if(iprint.eq.1.and.
0624      .mod(nphihy+1-nphi,10).eq.0)write(*,'(i3,$)')(nphihy+1-nphi)/10
0625       phi=phihy(nphi)
0626       fphi=2
0627       if(nphi.eq.1.or.nphi.eq.nphihy)fphi=1
0628       do ntau=1,ntauho
0629       tau=tauhy(ntau)
0630       ftau=2
0631       if(ntau.eq.1.or.ntau.eq.ntauho)ftau=1
0632       do meta=-netaho+netastep,netaho,netastep
0633         do ii=1,2
0634         if(ii.eq.1)then
0635           mt=meta-netastep
0636         else                 !if(ii.eq.2)
0637           mt=meta
0638         endif
0639         wwii(ii)=g1*waa(n1,mt,ntau,nphi)    +g2*waa(n2,mt,ntau,nphi)
0640         wwii(ii)=max(0.,wwii(ii))
0641         wxii(ii)=g1*vaa(1,n1,mt,ntau,nphi)  +g2*vaa(1,n2,mt,ntau,nphi)
0642         wyii(ii)=g1*vaa(2,n1,mt,ntau,nphi)  +g2*vaa(2,n2,mt,ntau,nphi)
0643         wzii(ii)=g1*vaa(3,n1,mt,ntau,nphi)  +g2*vaa(3,n2,mt,ntau,nphi)
0644         wrii(ii)=g1*raa(n1,mt,ntau,nphi)    +g2*raa(n2,mt,ntau,nphi)
0645         waii(ii)=g1*ddrdtau(n1,mt,ntau,nphi)+g2*ddrdtau(n2,mt,ntau,nphi)
0646         wbii(ii)=g1*ddrdphi(n1,mt,ntau,nphi)+g2*ddrdphi(n2,mt,ntau,nphi)
0647         wcii(ii)=g1*ddrdeta(n1,mt,ntau,nphi)+g2*ddrdeta(n2,mt,ntau,nphi)
0648         enddo
0649         jmax=jfac*netastep
0650         do j=0,jmax
0651           f=2
0652           if(j.eq.0.or.j.eq.jmax)f=1
0653           fall=fphi*ftau*f
0654           dVs=wwii(1)+j/float(jmax)*(wwii(2)-wwii(1))
0655           dVs=dVs*dall*fall
0656           vr=wxii(1)+j/float(jmax)*(wxii(2)-wxii(1))
0657           vt=wyii(1)+j/float(jmax)*(wyii(2)-wyii(1))
0658           v3=wzii(1)+j/float(jmax)*(wzii(2)-wzii(1))
0659           v1=vr*cos(phi)+vt*sin(phi)
0660           v2=vr*sin(phi)-vt*cos(phi)
0661           gmx=max(1e-8, 1.-vr*vr-vt*vt-v3*v3)
0662           gm=1./sqrt(gmx)
0663           rad=wrii(1)+j/float(jmax)*(wrii(2)-wrii(1))
0664           ieta=(meta-netastep)*jfac+j
0665           eta=ieta*dleta
0666           finc=10
0667           volu=abs(dVs)*finc
0668           io3=ioclude-3
0669           if(io3.lt.1.or.io3.gt.2)
0670      .     stop'in RestFrameFO: wrong ioclude (150808) '
0671           yie= volu * ffstat(io3,nspez) * ffstat(io3,mspez+2)
0672           do nsim=1,nsimu
0673             np=yie
0674             if(rangen().le.yie-np)np=np+1
0675             if(np.gt.0)then
0676               do n=1,np
0677                 r=rangen()*yie
0678                 i=0
0679                 do while(yie*ffstat(io3,i)/ffstat(io3,nspez).lt.r)
0680                  i=i+1
0681                 enddo
0682                 kss=klass(i)
0683                 id=ispez(i)
0684                 am=aspez(1,i)
0685                 fug=fuga(i)
0686                 js=jspez(i)
0687                 x=RanBoseFermi(aspez(io3,i),fug,js,tfo,istat)
0688                 p=x*tfo
0689                 e=sqrt(p**2+aspez(1,i)**2)
0690                 ex=sqrt(p**2+aspez(io3,i)**2)
0691                 !print*,id,e
0692                 u(3)=2.*rangen()-1.
0693                 angle=2.*pi*rangen()
0694                 u(1)=sqrt(1.-u(3)**2)*cos(angle)
0695                 u(2)=sqrt(1.-u(3)**2)*sin(angle)
0696                 q(1)=p*u(1)
0697                 q(2)=p*u(2)
0698                 q(3)=p*u(3)
0699                 q(4)=ex
0700                 call utlob3(-1, v1*gm , v2*gm , v3*gm , gm ,1e0
0701      .                   , q(1), q(2), q(3), q(4))
0702                 w1=FOFactor(wrii(1),waii(1),wbii(1),wcii(1),tau,phi,q)
0703                 w2=FOFactor(wrii(2),waii(2),wbii(2),wcii(2),tau,phi,q)
0704                 fof=w1+j/float(jmax)*(w2-w1)
0705                 fof=fof*dall*fall
0706                 probab=fof/volu/e
0707                 ij=1+probab/dprob
0708                 if(ij.ge.1.and.ij.le.100)
0709      .           probdi(jcent,ij)=probdi(jcent,ij)+1
0710                 if(rangen().le.probab)then !accept
0711                   if(io3.eq.2)q(4)=e
0712                   rap=eta +  0.5*log((q(4)+q(3))/(q(4)-q(3)))
0713                   pt2=q(1)**2+q(2)**2
0714                   ptr=sqrt(pt2)
0715                   pha=sign(1.,q(2))*acos(q(1)/ptr)
0716                   if(pha.lt.0.)pha=pha+2*pi
0717                   if(ianl.eq.1)then !~~~~~~~~~~~~~~~~~~~~~~~~~~
0718                     pz=sqrt(am**2+ptr**2)*sinh(rap)
0719                     pz2p=pz/sqrt(pz**2+ptr**2)
0720                     nrap=1+(rap-rapmin)/delrax
0721                     npha=1+pha/(2*pi/numiv)
0722                     sap=0.501
0723                     if(pz2p.le.-1.)then
0724                      nsap=0
0725                     elseif(pz2p.ge.1.)then
0726                      nsap=numiv+1
0727                     else
0728                      sap=0.5*log((1+pz2p)/(1-pz2p))
0729                      nsap=1+(sap-rapmin)/delrax
0730                     endif
0731                     nptv=1+ptr/delptv
0732                     if(nrap.ge.1.and.nrap.le.numiv.and.kss.gt.0)then
0733                       rapeta(kss,meta,nrap)=rapeta(kss,meta,nrap)+1
0734                       rapar(kss,nrap)=rapar(kss,nrap)+1
0735                       if(pt2.gt.0.)
0736      .                v2rapar(kss,nrap)=v2rapar(kss,nrap)
0737      .                 +(q(1)**2-q(2)**2)/pt2
0738                       if(npha.ge.1.and.npha.le.numiv)then
0739                         phaar(kss,nrap,npha)=phaar(kss,nrap,npha)+1
0740                       endif
0741                     endif
0742                     if(nsap.ge.1.and.nsap.le.numiv.and.kss.gt.0)then
0743                       sapar(kss,nsap)=sapar(kss,nsap)+1
0744                       if(pt2.gt.0.)
0745      .                v2sapar(kss,nsap)=v2sapar(kss,nsap)
0746      .                 +(q(1)**2-q(2)**2)/pt2
0747                     endif
0748                     if(nptv.ge.1.and.nptv.le.numiv.and.kss.gt.0
0749      .              .and.abs(rap).le.0.5)then
0750                       ptvar(kss,nptv)=ptvar(kss,nptv)+1
0751                     endif
0752                     if(nptv.ge.1.and.nptv.le.numiv.and.kss.gt.0
0753      .              .and.abs(sap).le.0.5)then
0754                       ptwar(kss,nptv)=ptwar(kss,nptv)+1
0755                     endif
0756                   endif !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
0757                   if(isto.eq.1)
0758      .            call FOStore(id,ptr,pha,am,rap,rad,phi,tau,eta)
0759                 endif
0760               enddo !np
0761             endif !np.gt.0
0762           enddo !nsim
0763         enddo !j
0764 
0765       enddo !meta
0766       enddo !ntau
0767       enddo !nphi
0768 
0769       if(iprint.eq.1)print*,' done'
0770       end
0771 
0772 c--------------------------------------------------------------------------------------
0773       subroutine FOStore(id,ptr,pha,am,rap,rad,phi,tau,eta)
0774 c--------------------------------------------------------------------------------------
0775       include 'epos.inc'
0776       include 'epos.inchy'
0777       common/cranphi/ranphi
0778       data ncntfoa/0/
0779       save ncntfoa
0780       ncntfoa=ncntfoa+1
0781 c      if(ncntfoa.eq.1)nptlb=nptl
0782       phinull=phievt+ranphi
0783       nptl=nptl+1
0784       if(nptl.gt.mxptl)
0785      . call utstop('FOStore: mxptl too small&')
0786       idptl(nptl)=id
0787       pptl(1,nptl)=ptr*cos(pha+phinull)
0788       pptl(2,nptl)=ptr*sin(pha+phinull)
0789       pptl(3,nptl)=sqrt(am**2+ptr**2)*sinh(rap)
0790       pptl(4,nptl)=sqrt(am**2+ptr**2)*cosh(rap)
0791       pptl(5,nptl)=am
0792       ityptl(nptl)=60
0793       xorptl(1,nptl)=rad*cos(phi+phinull)
0794       xorptl(2,nptl)=rad*sin(phi+phinull)
0795       xorptl(3,nptl)=tau*sinh(eta)
0796       xorptl(4,nptl)=tau*cosh(eta)
0797       end
0798 
0799 c--------------------------------------------------------------------------------------
0800       function FOFactor(wr,wa,wb,wc,tau,phi,q) ! 13 Aug 08 (same as gyx.f)
0801 c--------------------------------------------------------------------------------------
0802       real q(4)
0803       rad=wr
0804       qrd= cos(phi)*q(1)+sin(phi)*q(2)
0805       qtg= sin(phi)*q(1)-cos(phi)*q(2)
0806       fof=      - wa * tau * rad *q(4)
0807      .          +      rad * tau *qrd
0808      .          + wb * tau       *qtg
0809      .          - wc * rad       *q(3)
0810       FOFactor=fof
0811       end
0812 
0813 c---------------------------------------------------------------------------------
0814       function RanBoseFermi(am,fug,js,tfo,istat) ! 13 Aug 08 (same as gyx.f)
0815 c---------------------------------------------------------------------------------
0816 c   generates randomly x=p/tfo according to Bose/Fermi
0817 c------------------------------------------------------------------------------
0818   1   x=RanTherm(am/tfo)
0819       if(istat.eq.0)goto2  !Boltznann
0820       p=x*tfo
0821       e=sqrt(p**2+am**2)
0822       w=exp(-e/tfo)*fug
0823       if(js.eq.-1)then !bosons
0824         w0=exp(-am/tfo)*fug
0825         pacc=(1.-w0)/(1.-w)
0826       elseif(js.eq.1)then !fermions
0827         pacc=1./(1.+w)
0828       else
0829         stop'in RanBoseFermi: unknown statistics (080726) '
0830       endif
0831       !print*,'+++++',am,js,pacc
0832       if(rangen().gt.pacc)goto1
0833   2   RanBoseFermi=x
0834       end
0835 
0836 c------------------------------------------------------------------------------
0837       function RanTherm(a) ! 13 Aug 08 (same as gyx.f)
0838 c------------------------------------------------------------------------------
0839 c   generates a random number according to f(x) ~ x**2*exp(-sqrt(x*2+a**2))
0840 c   in the interval from zero to infinity
0841 c------------------------------------------------------------------------------
0842       !ntry=0
0843   1   i=2
0844       if(rangen().le.a**3/(a**3+3*a**2+6*a+6))i=1
0845       !ntry=ntry+1
0846       if(i.eq.1)then
0847         x=a*rangen()**(1./3.)
0848         if(rangen().gt.exp(a-sqrt(x**2+a**2)))goto1
0849       elseif(i.eq.2)then
0850         !f(z)~a**2+2*a*z+z**2)*exp(-z)  from zero to infty
0851         r=rangen()
0852         if(r.lt.a**2/(a**2+2*a+2))then
0853            z=-log(rangen())
0854         elseif(r.lt.(a**2+2*a)/(a**2+2*a+2))then
0855            z=-log(rangen())-log(rangen())
0856         else
0857            z=-log(rangen())-log(rangen())-log(rangen())
0858         endif
0859         x=a+z
0860         if(rangen().gt.exp(x-sqrt(x**2+a**2)))goto1
0861       endif
0862       RanTherm=x
0863       !ia=a*100
0864       !print*,ia,'   ',ntry,'    ',x
0865       end
0866 
0867 c---------------------------------------------------------------------------------
0868       subroutine DefineRapScale ! 13 Aug 08 (same as gyx.f except inc)
0869 c---------------------------------------------------------------------------------
0870       include 'epos.inchy'
0871       common/crap/nraphi,jfac /cho/netaho /cstep/netastep /crapi/delrapi
0872       deta=etahy(2)-etahy(1)
0873       write(*,'(a,$)')' DefineRapScale  ...'
0874       jfac=1
0875       delrap=deta
0876       do while(delrap.gt.delrapi)
0877         jfac=jfac+1
0878         delrap=deta/jfac
0879       enddo
0880       write(*,'(f7.2,i3,f7.2,5x,$)')deta,jfac,delrap
0881       nraphi=5
0882       rapmn=-delrap*nraphi
0883       do while(rapmn.gt.-2.)
0884         nraphi=nraphi+1
0885         if(nraphi.gt.nraphx)stop'(1539070608)   '
0886         rapmn=-delrap*nraphi
0887       enddo
0888       netaho=netahy-1
0889       do while(mod(netaho,netastep).ne.0)
0890       netaho=netaho-1
0891       enddo
0892       !print*,'$$$$$$ nraphi,rapmn,netaho:',nraphi,rapmn,netaho
0893       print*,' done'
0894       end
0895 
0896 c----------------------------------------------------------------------
0897       subroutine HydroFO(ier)
0898 c----------------------------------------------------------------------
0899       include 'epos.inc'
0900       include 'epos.inchy'
0901       common/cen/ncentr
0902       parameter (klax=5)
0903       parameter(numiv=100)!,rapmax=5,rapmin=-rapmax)
0904       common/cana1b/rapar(klax,numiv),v2rapar(klax,numiv)
0905       common/cana1d/sapar(klax,numiv),v2sapar(klax,numiv)
0906       common/cstep/netastep /crapi/delrapi
0907       double precision seedp
0908 
0909       data icntcf /0/
0910       save icntcf
0911       icntcf=icntcf+1
0912       if(icntcf.eq.1)then
0913         call ranfgt(seedp)      !not to change the seed ...
0914         do n=1,klax
0915          do nu=1,numiv
0916           rapar(n,nu)=0
0917           v2rapar(n,nu)=0
0918           sapar(n,nu)=0
0919           v2sapar(n,nu)=0
0920          enddo
0921         enddo
0922         netastep=1
0923         delrapi=0.2
0924         call HydroTable2(0)
0925         call SymmTab(1,ncenthy)
0926         call DefineParticles
0927         call DefineRapScale
0928         call ranfst(seedp)      ! ... after this initialization
0929       endif
0930       if(iappl.eq.9)then
0931         call CheckBimp(ier)
0932         if(ier.eq.1)return
0933         call GetNpart
0934       endif
0935       call RestFrameFO(1,bimevt,1,ier,0,1,1)
0936       end
0937 
0938 c----------------------------------------------------------------------
0939       subroutine xEnergy
0940 c----------------------------------------------------------------------
0941       include 'epos.inc'
0942       parameter (mspez=54)
0943       common/cspez4/ffstat(2,0:mspez+2)
0944       common/ccsum/eesum,hhsum
0945       eesum=ffstat(1,mspez+1)
0946       write(ifhi,'(a)')       '!##################################'
0947       write(ifhi,'(a,i3)')    '!   energy     '
0948       write(ifhi,'(a)')       '!##################################'
0949       write(ifhi,'(a)') ' openhisto htyp pgs xmod lin ymod lin '
0950       write(ifhi,'(a)') ' array 2'
0951       write(ifhi,'(2e11.3)')eesum*volu,0.005
0952       write(ifhi,'(a)') ' endarray closehisto plot 0'
0953       end
0954 
0955 c----------------------------------------------------------------------
0956       subroutine GetNpart
0957 c----------------------------------------------------------------------
0958       include 'epos.inc'
0959       include 'epos.inchy'
0960       mcent=2
0961       do while(mcent.lt.mcenty.and.bcenty(mcent).lt.bimevt)
0962         mcent=mcent+1
0963       enddo
0964       m1=mcent-1
0965       m2=mcent
0966       hrac=(bimevt-bcenty(m1))/(bcenty(m2)-bcenty(m1))
0967       h1=1-hrac
0968       h2=hrac
0969       npartic=h1*pcenty(m1)+h2*pcenty(m2)
0970       npjevt=npartic/2
0971       ntgevt=npartic-npjevt
0972       !print*,'++++++++'
0973       !.    ,mcent, bimevt, bcenty(m1),bcenty(m2),h1,h2,npartic
0974       end
0975 
0976 c----------------------------------------------------------------------
0977       subroutine CheckBimp(ier)
0978 c----------------------------------------------------------------------
0979       include 'epos.inc'
0980       include 'epos.inchy'
0981       data icntcb /0/
0982       save icntcb, wref,bref
0983       icntcb=icntcb+1
0984       if(icntcb.eq.1)then
0985         nmax=0
0986         wmax=0
0987         do n=1,nbimp
0988          if(bimpar(2,n).gt.wmax)then
0989            wmax=bimpar(2,n)
0990            nmax=n
0991          endif
0992         enddo
0993         nref=nmax*0.75
0994         bref=bimpar(1,nref)
0995         wref=bimpar(2,nref)
0996       endif
0997       ier=0
0998       if(bimevt.lt.bref)then
0999         q=1
1000       else
1001         w=wref/bref*bimevt
1002         n=2
1003         do while(n.lt.nbimp.and.bimpar(1,n).lt.bimevt)
1004           n=n+1
1005         enddo
1006         n1=n-1
1007         n2=n
1008         frac=(bimevt-bimpar(1,n1))/(bimpar(1,n2)-bimpar(1,n1))
1009         g1=1-frac
1010         g2=frac
1011         wx=g1*bimpar(2,n1)+g2*bimpar(2,n2)
1012         wx=max(0.,wx)
1013         q=wx/w
1014       endif
1015       if(rangen().gt.q)ier=1
1016       !if(ier.eq.0)print*,'+++++',bimevt,q
1017 
1018       end
1019 
1020 c----------------------------------------------------------------------
1021       subroutine xCoopFryPt(kss)
1022 c----------------------------------------------------------------------
1023       include 'epos.inc'
1024       parameter(numiv=100,ptvmax=2,klax=5)
1025       common/cana1e/ptvar(klax,numiv),ptwar(klax,numiv)
1026       delptv= ptvmax   / numiv
1027       pi=3.1415927
1028       write(ifhi,'(a)')       '!##################################'
1029       write(ifhi,'(a,i3)')    '!   pt    '
1030       write(ifhi,'(a)')       '!##################################'
1031       write(ifhi,'(a)') 'openhisto htyp lin xmod lin ymod log '
1032       write(ifhi,'(a,f7.2)') 'xrange 0 2 '
1033       write(ifhi,'(a)') 'txt  "xaxis p?t!"'
1034       write(ifhi,'(a)') 'txt  "yaxis dn / 2[p] p?t! dp?t!"'
1035       write(ifhi,'(a)') 'array 2'
1036       do n=1,numiv
1037         x=(n-0.5)*delptv
1038         y=ptvar(kss,n)/float(nevent)/delptv/2./pi/x
1039         write(ifhi,'(2e14.6)')x,y
1040       enddo
1041       write(ifhi,'(a)') 'endarray closehisto plot 0'
1042 c      write(ifhi,'(a)') 'openhisto htyp lin xmod lin ymod lin '
1043 c      write(ifhi,'(a,f7.2)') 'xrange 0 2 '
1044 c      write(ifhi,'(a)') 'array 2'
1045 c      do n=1,numiv
1046 c        x=(n-0.5)*delptv
1047 c        y=0
1048 c        if(ptwar(kss,n).gt.0.)y=ptwar(kss,n)/float(nevent)/delptv/2./pi/x
1049 c        write(ifhi,'(2e14.6)')x,y
1050 c      enddo
1051 c      write(ifhi,'(a)') 'endarray closehisto plot 0'
1052       end
1053 
1054 c----------------------------------------------------------------------
1055       subroutine xCoopFryRap(kss)
1056 c----------------------------------------------------------------------
1057       include 'epos.inc'
1058       parameter(numiv=100,rapmax=5,rapmin=-rapmax,klax=5)
1059       common/cana1b/rapar(klax,numiv),v2rapar(klax,numiv)
1060       common/cana1d/sapar(klax,numiv),v2sapar(klax,numiv)
1061       delrax= 2*rapmax   / numiv
1062       write(ifhi,'(a)')       '!##################################'
1063       write(ifhi,'(a,i3)')    '!   rapidity distribution     '
1064       write(ifhi,'(a)')       '!##################################'
1065       write(ifhi,'(a)') ' openhisto htyp lin xmod lin ymod lin '
1066       write(ifhi,'(a)') ' xrange -5 5'
1067       write(ifhi,'(a)') ' txt  "xaxis y "'
1068       write(ifhi,'(a)') ' txt  "yaxis dn/dy "'
1069       write(ifhi,'(a)') ' array 2'
1070       do n=1,numiv
1071         x=rapmin+(n-0.5)*delrax
1072         y=rapar(kss,n)/float(nevent)/delrax
1073         write(ifhi,'(2e11.3)')x,y
1074       enddo
1075       write(ifhi,'(a)') ' endarray closehisto plot 0'
1076 c      write(ifhi,'(a)') ' openhisto htyp lin xmod lin ymod lin '
1077 c      write(ifhi,'(a)') ' xrange -5 5'
1078 c      write(ifhi,'(a)') ' txt  "xaxis  [c] "'
1079 c      write(ifhi,'(a)') ' txt  "yaxis dn/d[c]"'
1080 c      write(ifhi,'(a)') ' array 2'
1081 c      do n=1,numiv
1082 c        x=rapmin+(n-0.5)*delrax
1083 c        y=sapar(kss,n)/float(nevent)/delrax
1084 c        write(ifhi,'(2e11.3)')x,y
1085 c      enddo
1086 c      write(ifhi,'(a)') ' endarray closehisto plot 0'
1087       end
1088 
1089 c----------------------------------------------------------------------
1090       subroutine xCoopFryV2(kss)
1091 c----------------------------------------------------------------------
1092       include 'epos.inc'
1093       parameter(numiv=100,rapmax=5,rapmin=-rapmax,klax=5)
1094       common/cana1b/rapar(klax,numiv),v2rapar(klax,numiv)
1095       common/cana1d/sapar(klax,numiv),v2sapar(klax,numiv)
1096       delrax= 2*rapmax   / numiv
1097       write(ifhi,'(a)')       '!##################################'
1098       write(ifhi,'(a,i3)')    '!   v2     '
1099       write(ifhi,'(a)')       '!##################################'
1100       write(ifhi,'(a)') ' openhisto htyp lin xmod lin ymod lin '
1101       write(ifhi,'(a,f7.2)') ' xrange -5 5 yrange 0 auto'
1102       write(ifhi,'(a)') ' txt  "xaxis y "'
1103       write(ifhi,'(a)') ' txt  "yaxis v?2!"'
1104       write(ifhi,'(a)') ' array 2'
1105       do n=1,numiv
1106         x=rapmin+(n-0.5)*delrax
1107         y=0
1108         if(rapar(kss,n).gt.0.)y=v2rapar(kss,n)/rapar(kss,n)
1109         write(ifhi,'(2e11.3)')x,y
1110       enddo
1111       write(ifhi,'(a)') ' endarray closehisto plot 0'
1112 c      write(ifhi,'(a)') ' openhisto htyp lin xmod lin ymod lin '
1113 c      write(ifhi,'(a,f7.2)') ' xrange -5 5 yrange 0 auto'
1114 c      write(ifhi,'(a)') ' txt  "xaxis  [c]"'
1115 c      write(ifhi,'(a)') ' txt  "yaxis v?2!"'
1116 c      write(ifhi,'(a)') ' array 2'
1117 c      do n=1,numiv
1118 c        x=rapmin+(n-0.5)*delrax
1119 c        y=0
1120 c        if(sapar(kss,n).gt.0.)y=v2sapar(kss,n)/sapar(kss,n)
1121 c        write(ifhi,'(2e11.3)')x,y
1122 c      enddo
1123 c      write(ifhi,'(a)') ' endarray closehisto plot 0'
1124       end
1125 
1126 c----------------------------------------------------------------------
1127       subroutine DropletDecay(ip,iret)
1128 c----------------------------------------------------------------------
1129       ! statistical hadronization with imposed flow
1130       !   (to be distiguished from real hydro flow as in StaHadDF)
1131       !---------------------------------------------------------------
1132       include 'epos.inc'
1133       include 'epos.inchy'
1134       common/cxyzt/xptl(mxptl),yptl(mxptl),zptl(mxptl),tptl(mxptl)
1135      *,optl(mxptl),uptl(mxptl),sptl(mxptl),rptl(mxptl,3)
1136       parameter(maxp=500)
1137       common/confg/np,amass(maxp),ident(maxp),pcm(5,maxp),wtxlog,wtlog
1138       common/citer/iter,itermx
1139       integer jc(nflav,2)
1140       double precision p(5),c(5)!,yyrmax
1141       parameter(maxit=50000)
1142       common/count/nacc,nrej,naccit(maxit),nptot,npit(maxit)
1143       dimension uu(4),pe(5),pa(5)
1144       common/ylong/ylong(maxp)
1145       common/vradi/vradi(maxp),vtang(maxp),phifop(maxp),radfop(maxp)
1146      . ,taufop(maxp)
1147       common/xxxspecsy/ndrop(-4:4,-4:4,-4:4)
1148       common/cdelzet/delzet,delsgr /cvocell/vocell
1149       common/cen/ncentr /ctauhu/ntauhu(ncenthx,1-netahx:netahx-1)
1150       common/cranphi/ranphi /ctfo/tfo
1151       double precision seedp
1152       !data vv2 /0./ nvv2 /0/  vv3 /0./
1153       !save vv2,nvv2,vv3
1154       data icntsta /0/
1155       save icntsta
1156       icntsta=icntsta+1
1157       if(ioclude.eq.1)stop'ioclude.eq.1 no longer supported.    '
1158       if(ioclude.eq.2)stop'ioclude.eq.2 no longer supported.    '
1159 
1160       call utpri('DropletDecay',ish,ishini,4)
1161 
1162       if(ish.ge.3)then
1163       write(ifch,140)
1164   140 format(/' ----------------------------------'/
1165      *        '    droplet decay'/
1166      *        ' ----------------------------------')
1167       write(ifch,*)'droplet:'
1168       call alist('&',ip,ip)
1169       endif
1170       call ManiParticles(1,1) !store parameters, to be restored at the end
1171 
1172       if(icntsta.eq.1)then
1173         call ranfgt(seedp)      !not to change the seed ...
1174         call ManiParticles(92,0)
1175         !redefine parameters, for this subroutine; for the moment same parameters
1176         ! for remnants and droplets, this could easily be changed
1177         call ManiParticles(1,2) !store local parameters
1178         call ranfst(seedp)      ! ... after this initialization
1179       else
1180         call ManiParticles(-1,2) !restore local parameters
1181       endif
1182 
1183       iret=0
1184       do j=1,5
1185       c(j)=pptl(j,ip)
1186       enddo
1187 
1188       call idquac(ip,nqi,nsi,nai,jc)
1189       keu=jc(1,1)-jc(1,2)
1190       ked=jc(2,1)-jc(2,2)
1191       kes=jc(3,1)-jc(3,2)
1192       kec=jc(4,1)-jc(4,2)
1193       keb=jc(5,1)-jc(5,2)
1194       ket=jc(6,1)-jc(6,2)
1195       !print*,'droplet uds=',keu,ked,kes,'   E=',pptl(5,ip)
1196 
1197     !~~~~~define masses~~~~~~~~~~~~~~~~
1198       amin=utamnu(keu,ked,kes,kec,keb,ket,5)
1199       aumin=amuseg+yrmaxi
1200       ipo=ip
1201       if(ityptl(ip).eq.60)ipo=iorptl(ip)
1202       tecmor=pptl(5,ipo)
1203       tecm=pptl(5,ip)
1204       tecmxx=tecm
1205 
1206       if(iappl.eq.4.or.iorsdf.ne.3
1207      &.or.ityptl(ip).eq.40.or.ityptl(ip).eq.50)then !not for droplets from remnants
1208         yrmax=0
1209       else
1210         yrmax=yrmaxi             !define in ainit
1211       endif
1212 
1213 
1214     !~~~~~redefine energy in case of imposed radial flow~~~~~~~~~~~~~~~~
1215       fradfli=1.
1216       if(yrmax.gt.1e-5)fradfli=fradflii                !define in ainit
1217       if(tecm*fradfli.lt.amin.and.tecm.gt.0.)fradfli=1.1*amin/tecm  !if flow too large, do something anyway (saturation of flow)
1218       if(yrmax.gt.1e-5.and.tecmor.gt.aumin
1219      &.and.tecm*fradfli.gt.amin) then
1220          tecm=tecm*fradfli
1221       else
1222         yrmax=0.
1223         fradfli=1.
1224       endif
1225 
1226     !~~~~~redefine energy in case of long coll flow
1227       if(iappl.eq.4.or.iorsdf.ne.3
1228      &.or.ityptl(ip).eq.40.or.ityptl(ip).eq.50)then !not for droplets from remnants
1229         yco=0
1230       else
1231        if(ylongmx.lt.0.)then
1232         yco=delzet * 1.75
1233        else
1234         yco=ylongmx
1235        endif
1236       endif
1237       tecmx=tecm
1238       if(yco.gt.0..and.tecmor.gt.aumin) then
1239         tecm=tecm/sinh(yco)*yco
1240       else
1241         yco=0.
1242       endif
1243       !print*,'========= cluster energy: ',pptl(5,ip),tecmx,tecm
1244 
1245     !~~~~~define volu~~~~~~~~~~~~~~~~
1246       volu=tecm/epscri(3)
1247 
1248     !~~~~~~~~~decay~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1249         ibarini=(keu+ked+kes+kec+keb+ket)/3 !for baryon number conservation
1250         niter=0
1251     4   niter=niter+1
1252         call GraCan
1253         energ=0
1254         ibaryon=ibarini
1255         do i=1,np
1256           energ=energ+pcm(4,i)
1257           if(abs(ident(i)).gt.1000)ibaryon=ibaryon-sign(1,ident(i))
1258         enddo
1259         !print*,'++4+++',energ,tecm,tecm/volu,ibaryon
1260         if((abs(energ-tecm).gt.0.1.or.ibaryon.ne.0)
1261      &       .and.niter.le.maxit) goto 4 !comment for energy non-conservation
1262         if(niter.gt.maxit)then
1263          iret=1
1264          goto1000
1265         endif
1266 
1267     !~~~~~~~~~~long coll flow -> particles~~~~~~~~~~~~~~~~
1268       tecmxxx=tecm
1269       if(yco.gt.0.) then
1270         errlim=0.0001
1271         tecm=tecmx
1272         niter=0
1273  611    energ=0.
1274         niter=niter+1
1275         do i=1,np
1276           ylong(i)=(2*rangen()-1)*yco
1277           uu(3)=sinh(ylong(i))
1278           uu(4)=cosh(ylong(i))
1279           energ=energ+uu(4)*pcm(4,i)+uu(3)*pcm(3,i)
1280         enddo
1281         if(abs(energ-tecm).gt.0.1)then
1282           goto 611
1283         elseif(niter.ge.1000)then
1284           if(ish.ge.1)write(ifch,*)'Long Flow failed:'
1285      &                             ,energ,tecm
1286           yco=0
1287           tecm=tecmxxx
1288           goto 400
1289         endif
1290         !print*,'===== energy after flow boosts',energ,'   soll: ',tecm
1291         do j=1,4
1292           pe(j)=0.
1293         enddo
1294         do i=1,np
1295           uu(1)= 0
1296           uu(2)= 0
1297           uu(3)= sinh(ylong(i))
1298           uu(4)= cosh(ylong(i))
1299           call utlob3(-1,uu(1),uu(2),uu(3),uu(4),1e0
1300      *         , pcm(1,i), pcm(2,i), pcm(3,i), pcm(4,i))
1301           do j=1,4
1302           pe(j)=pe(j)+pcm(j,i)
1303           enddo
1304         enddo
1305         pe(5)=sqrt(pe(4)**2-pe(3)**2-pe(2)**2-pe(1)**2)
1306         !write(6,'(a,5e11.3)')'flow boosts',pe
1307         do j=1,4
1308           pa(j)=0.
1309         enddo
1310         do i=1,np
1311           call utlob3(1,pe(1),pe(2),pe(3),pe(4),pe(5)
1312      *         , pcm(1,i), pcm(2,i), pcm(3,i), pcm(4,i))
1313           do j=1,4
1314             pa(j)=pa(j)+pcm(j,i)
1315           enddo
1316         enddo
1317         pa(5)=sqrt(pa(4)**2-pa(3)**2-pa(2)**2-pa(1)**2)
1318         !write(6,'(a,5e11.3)')' cms boost ',pa
1319         esoll=tecm
1320         scal=1.
1321         do ipass=1,200
1322           sum=0.
1323           do  j=1,np
1324             do k=1,3
1325               pcm(k,j)=scal*pcm(k,j)
1326             enddo
1327             pcm(4,j)=sqrt(pcm(1,j)**2+pcm(2,j)**2+pcm(3,j)**2
1328      *           +amass(j)**2)
1329             sum=sum+pcm(4,j)
1330           enddo
1331           scal=esoll/sum
1332           !write(6,*)'ipass,scal,e,esoll:'
1333           !    $         ,ipass,scal,sum,esoll
1334           if(abs(scal-1.).le.errlim) goto301
1335         enddo
1336  301    continue
1337         do j=1,4
1338           pa(j)=0.
1339         enddo
1340         do i=1,np
1341           do j=1,4
1342             pa(j)=pa(j)+pcm(j,i)
1343           enddo
1344         enddo
1345         pa(5)=sqrt(pa(4)**2-pa(3)**2-pa(2)**2-pa(1)**2)
1346         !write(6,'(a,5e11.3)')' rescaling ',pa
1347       endif
1348  400  continue
1349 
1350     !~~~~~~~~~~radial flow -> particles~~~~~~~~~~~~~~~~~~
1351       if(fradfli.lt.1.) then
1352         errlim=0.0001
1353         tecm=tecmxx
1354         phinull=phievt+ranphi
1355         do n=1,np
1356           vradi(n)=tanh(sqrt(rangen())*yrmax)
1357           vtang(n)=0.
1358           phifop(n)=rangen()*2*pi
1359           radfop(n)=0.
1360           taufop(n)=0.
1361         enddo
1362         energ=0.
1363         do n=1,np
1364           co=cos(phifop(n)+phinull)
1365           si=sin(phifop(n)+phinull)
1366           v1=vradi(n)*co+vtang(n)*si
1367           v2=vradi(n)*si-vtang(n)*co
1368           v3=0
1369           a=1.-v1*v1-v2*v2-v3*v3
1370           if(a.le.0.)stop'in DropletDecay (20032008)        '
1371           gm=1./sqrt(a)
1372 
1373           uu(1)=v1*gm
1374           uu(2)=v2*gm
1375           uu(3)=v3*gm
1376           uu(4)=gm
1377           call utlob3(-1,uu(1),uu(2),uu(3),uu(4),1e0
1378      *         , pcm(1,n), pcm(2,n), pcm(3,n), pcm(5,n))
1379           energ=energ+sqrt(amass(n)**2
1380      *                     +pcm(1,n)**2+pcm(2,n)**2+pcm(3,n)**2)
1381         enddo
1382         esoll=tecm
1383         !print*,tecm,energ
1384         iiscale=1
1385         if(iiscale.eq.1)then
1386           scal=1.
1387           do ipass=1,200
1388             sum=0.
1389             do  j=1,np
1390               do k=1,3
1391                 pcm(k,j)=scal*pcm(k,j)
1392               enddo
1393               pcm(4,j)=sqrt(pcm(1,j)**2+pcm(2,j)**2+pcm(3,j)**2
1394      *           +amass(j)**2)
1395               sum=sum+pcm(4,j)
1396             enddo
1397             scal=esoll/sum
1398             if(abs(scal-1.).le.errlim) goto300
1399           enddo
1400   300     continue
1401         endif
1402       else
1403         phinull=0.
1404         do n=1,np
1405           vradi(n)=0.
1406           vtang(n)=0.
1407           phifop(n)=0.
1408           radfop(n)=0.
1409           taufop(n)=0.
1410         enddo
1411       endif
1412     !~~~~~~~~~~~~~~~
1413 
1414       nptlb=nptl
1415       do n=1,np
1416         nptl=nptl+1
1417         if(nptl.gt.mxptl)
1418      .   call utstop('DropletDecay: mxptl too small&')
1419         idptl(nptl)=ident(n)
1420         do j=1,4
1421           p(j)=pcm(j,n)
1422         enddo
1423         p(5)=amass(n)
1424         call utlob2(-1,c(1),c(2),c(3),c(4),c(5),p(1),p(2),p(3),p(4),10)
1425         do j=1,5
1426           pptl(j,nptl)=p(j)
1427         enddo
1428         if(fradfli.lt.1.) then
1429           ityptl(nptl)=60
1430         elseif(ityptl(ip).eq.40.or.ityptl(ip).eq.50)then
1431           ityptl(nptl)=ityptl(ip)+1
1432         else
1433           ityptl(nptl)=19
1434         endif
1435         if(ityptl(ip).eq.60)then
1436          if(ityptl(nptl).eq.60)then
1437           ipo=iorptl(ip)
1438           phi=phifop(n)
1439           tau=taufop(n)
1440           r=radfop(n)
1441           !---add r-randomness
1442           !dr=5
1443           !do while(dr.lt.-2.or.dr.gt.2.)
1444           ! dr=sqrt(3.)*(rangen()+rangen()+rangen()+rangen()-2)
1445           !enddo
1446           !r=r+dr
1447           !---------
1448           zeta=0.5*log((p(4)+p(3))/(p(4)-p(3)))
1449           !deleta=etahy(2)-etahy(1)
1450           !zeta=zetaor-0.5*delzet+delzet*rangen()
1451           z=tau*sinh(zeta)
1452           t=tau*cosh(zeta)
1453           xorptl(1,nptl)=r*cos(phi+phinull)
1454           xorptl(2,nptl)=r*sin(phi+phinull)
1455           xorptl(3,nptl)=z
1456           xorptl(4,nptl)=t
1457          else
1458           xorptl(1,nptl)=xorptl(1,ip)
1459           xorptl(2,nptl)=xorptl(2,ip)
1460           xorptl(3,nptl)=xorptl(3,ip)
1461           xorptl(4,nptl)=xorptl(4,ip)
1462          endif
1463         endif
1464       enddo
1465       if(ish.ge.3)then
1466         write(ifch,*)'decay products:'
1467         call alist('&',nptlb+1,nptl)
1468         if(ish.ge.5)then
1469           write(ifch,*)'momentum sum:'
1470           do kk=1,5
1471           pptl(kk,nptl+1)=0
1472           do ii=nptlb+1,nptl
1473           pptl(kk,nptl+1)=pptl(kk,nptl+1)+pptl(kk,ii)
1474           enddo
1475           pptl(kk,nptl+2)=c(kk)
1476           enddo
1477           call alist('&',nptl+1,nptl+2)
1478         endif
1479       endif
1480 
1481  1000 continue
1482       call ManiParticles(-1,1)
1483       call utprix('DropletDecay',ish,ishini,4)
1484       end
1485 
1486 c--------------------------------------------------------------------
1487       subroutine ManiParticles(isens,isto)
1488 c--------------------------------------------------------------------
1489       ! isens=1 -> store,  isens=-1 -> restore,  isens=9... -> redefine
1490       !-----------------------------------------------------------
1491       include 'epos.inc'
1492       include 'epos.inchy'
1493       parameter (mspez=54,msto=2)
1494       common/cspez1/nspez,ispez(mspez),aspez(2,mspez),gspez(mspez)
1495       common/cspez2/kspez,jspez(mspez)
1496       common/cspez3/fuga(mspez)
1497       common/cspez4/ffstat(2,0:mspez+2)
1498       common/ctfo/tfo
1499       common/sp1/ffstatSave(2,0:mspez+2,msto),fugaSave(mspez,msto)
1500       common/sp2/iocludeSave(msto),kspezSave(msto),nspezSave(msto)
1501       common/sp3/tfoSave(msto),meosmuSave(msto),epscri3Save(msto)
1502       if(isens.eq.1)then !~~~~~~~~~~~~~~store~~~~~~~~~~~~~~
1503         iocludeSave(isto)=ioclude
1504         kspezSave(isto)=kspez
1505         nspezSave(isto)=nspez
1506         if(kspez.gt.0)then
1507          do m=kspez,nspez
1508           ffstatSave(1,m,isto)=ffstat(1,m)
1509           ffstatSave(2,m,isto)=ffstat(2,m)
1510           fugaSave(m,isto)=fuga(m)
1511          enddo
1512         endif
1513         do i=1,2
1514          do j=1,2
1515           ffstatSave(i,mspez+j,isto)=ffstat(i,mspez+j)
1516          enddo
1517         enddo
1518         tfoSave(isto)=tfo
1519         meosmuSave(isto)=meosmu
1520         epscri3Save(isto)=epscri(3)
1521       elseif(isens.eq.-1)then !~~~~~~~~~~~~~~restore~~~~~~~~~~~~~~
1522         ioclude=iocludeSave(isto)
1523         kspez=kspezSave(isto)
1524         nspez=nspezSave(isto)
1525         if(kspez.gt.0)then
1526           do m=kspez,nspez
1527            ffstat(1,m)=ffstatSave(1,m,isto)
1528            ffstat(2,m)=ffstatSave(2,m,isto)
1529            fuga(m)=fugaSave(m,isto)
1530           enddo
1531         endif
1532         do i=1,2
1533          do j=1,2
1534           ffstat(i,mspez+j)=ffstatSave(i,mspez+j,isto)
1535          enddo
1536         enddo
1537         tfo=tfoSave(isto)
1538         meosmu=meosmuSave(isto)
1539         epscri(3)=epscri3Save(isto)
1540       elseif(isens.eq.92)then !~~~~~~~~~~~~~~redefine set2~~~~~~~~~~~~~~
1541         !used for remnant and droplet decay --- may be modified!!!
1542         ioclude=5             !choice of hadronization (4 or 5)
1543         kspez=1               !choice of
1544         nspez=54              !   particles
1545         tfo=0.130             !freeze out temperature used in this routine
1546         meosmu=0              !no chemical potentials used in this routine
1547         epscri(3)=0.0765      !should be consitent with tfo see EoS table below
1548         call DefineParticles
1549       elseif(isens.eq.93)then !~~~~~~~~~~~~~~redefine set3~~~~~~~~~~~~~
1550         !used for application micro
1551         kspez=1               !choice of
1552         nspez=54              !   particles
1553         tfo=0.130             !freeze out temperature used in this routine
1554         meosmu=0              !no chemical potentials used in this routine
1555         epscri(3)=0.0765      !should be consitent with tfo see EoS table below
1556         call DefineParticles
1557       endif
1558       !--------------------------!
1559       ! epscri(3)  !    tfo      !
1560       !--------------------------!
1561       !   0.045    !   0.11954   !
1562       !   0.050    !   0.12157   !
1563       !   0.055    !   0.12327   !
1564       !   0.060    !   0.12496   !
1565       !   0.065    !   0.12665   !
1566       !   0.070    !   0.12801   !
1567       !   0.075    !   0.1297    !
1568       !   0.080    !   0.13072   !
1569       !   0.085    !   0.13207   !
1570       !   0.090    !   0.13343   !
1571       !   0.095    !   0.13444   !
1572       !   0.100    !   0.13546   !
1573       !   0.105    !   0.13648   !
1574       !   0.110    !   0.13749   !
1575       !   0.115    !   0.13851   !
1576       !   0.120    !   0.13952   !
1577       !   0.125    !   0.1402    !
1578       !--------------------------!
1579       end
1580 
1581 c------------------------------------------------------------------------------
1582       subroutine xSpaceTime
1583 c------------------------------------------------------------------------------
1584       include 'epos.inc'
1585       if(iSpaceTime.eq.1.and.ioclude.gt.1)then
1586          call xCoreCorona(0,0)
1587          do m=0,4,2
1588           do meta=-m,m,max(1,2*m)
1589            call xFoMass(meta)
1590            call xFoRadius(meta)
1591            call xFoRadVelocity(meta)
1592            call xFoTanVelocity(meta)
1593           enddo
1594          enddo
1595          do m=0,4,2
1596           do meta=-m,m,max(1,2*m)
1597            call xFreezeOutTauX(meta)
1598           enddo
1599          enddo
1600          call xFeff
1601          !call xFreezeOutTauEta
1602          !call xFreezeOutTZ
1603          call xEos
1604       elseif(iSpaceTime.eq.1)then
1605          call xCoreCorona(0,0)
1606          !stop'bjinta: space-time plots require ioclude>1.          '
1607       endif
1608       end
1609 
1610 c------------------------------------------------------------------------------
1611       subroutine xEos
1612 c------------------------------------------------------------------------------
1613       include'epos.inc'
1614       include'epos.inchy'
1615       common/latt/temp(14),press(14),epsi(14)
1616       call Lattice
1617 
1618       write(ifhi,'(a)')       '!##################################'
1619       write(ifhi,'(a,i3)')    '!   eos 1     '
1620       write(ifhi,'(a)')       '!##################################'
1621       write(ifhi,'(a)') ' openhisto htyp lin xmod lin ymod lin '
1622       write(ifhi,'(a)') ' xrange 0.   3.  yrange 0 auto'
1623       write(ifhi,'(a)') ' txt  "xaxis [e] (GeV/fm^3!)"'
1624       write(ifhi,'(a)') ' txt  "yaxis p([e]) (GeV/fm^3!)"'
1625       write(ifhi,'(a)') ' array 2'
1626       do m=1,meos
1627        if(eos(1,m).le.3.)write(ifhi,'(2e11.3)')eos(1,m),eos(2,m)
1628       enddo
1629       write(ifhi,'(a)') ' endarray closehisto plot 0-'
1630       write(ifhi,'(a)') ' openhisto htyp poc  '
1631       write(ifhi,'(a)') ' array 2'
1632       do m=2,13
1633        if(epsi(m).le.3.)write(ifhi,'(2e11.3)')epsi(m),press(m)
1634       enddo
1635       write(ifhi,'(a)') 'endarray closehisto plot 0'
1636       write(ifhi,'(a)')       '!##################################'
1637       write(ifhi,'(a,i3)')    '!   eos 2     '
1638       write(ifhi,'(a)')       '!##################################'
1639       write(ifhi,'(a)') ' openhisto htyp lin xmod lin ymod lin '
1640       write(ifhi,'(a)') ' xrange 0.   3.  yrange 0 auto'
1641       write(ifhi,'(a)') ' txt  "xaxis [e] (GeV/fm^3!)"'
1642       write(ifhi,'(a)') ' txt  "yaxis T([e]) (GeV)"'
1643       write(ifhi,'(a)') ' array 2'
1644       do m=1,meos
1645        if(eos(1,m).le.3.)write(ifhi,'(2e11.3)')eos(1,m),eos(3,m)
1646       enddo
1647       write(ifhi,'(a)') ' endarray closehisto plot 0-'
1648       write(ifhi,'(a)') ' openhisto htyp poc  '
1649       write(ifhi,'(a)') ' array 2'
1650       do m=2,13
1651        if(epsi(m).le.3.)write(ifhi,'(2e11.3)')epsi(m),temp(m)
1652       enddo
1653       write(ifhi,'(a)') 'endarray closehisto plot 0'
1654       !write(ifhi,'(a)')       '!##################################'
1655       !write(ifhi,'(a,i3)')    '!   eosmu     '
1656       !write(ifhi,'(a)')       '!##################################'
1657       !write(ifhi,'(a)') ' openhisto htyp lin xmod lin ymod lin '
1658       !write(ifhi,'(a)') ' xrange 0.   0.2  yrange 0 auto'
1659       !write(ifhi,'(a)') ' txt  "xaxis T (GeV)"'
1660       !write(ifhi,'(a)') ' txt  "yaxis [m]?[p]!(T) (GeV)"'
1661       !write(ifhi,'(a)') ' array 2'
1662       !do m=1,meosmu
1663       ! if(eos(1,m).le.3.)write(ifhi,'(2e11.3)')eosmu(1,m),eosmu(2,m)
1664       !enddo
1665       !write(ifhi,'(a)') 'endarray closehisto plot 0'
1666       end
1667 
1668 c------------------------------------------------------------------------------
1669       subroutine Lattice
1670 c------------------------------------------------------------------------------
1671       parameter (hc=0.1973,Tc=0.170/hc)
1672       common/latt/temp(14),press(14),epsi(14)
1673       real t2tc(14),p2T4(14),e2T4(14)
1674       ! T/Tc  no units
1675       data t2tc /0.80,0.87,0.96,1.02,1.07,1.14,1.20,1.28
1676      *          ,1.35,1.52,1.70,1.90,2.24,2.55/
1677       ! p/T^4,  no units
1678       data p2T4 /0.05,0.15,0.39,0.60,0.86,1.12,1.40,1.66
1679      *          ,1.91,2.32,2.65,2.89,3.19,3.41/
1680       do i=1,14
1681        temp(i)=t2tc(i)*Tc*hc                ! in GeV
1682        press(i)=p2T4(i)*(t2tc(i)*Tc)**4*hc  ! in GeV/fm^3
1683       enddo
1684       do i=2,13
1685         f1=p2T4(i-1)*(t2tc(i-1)*Tc)**4
1686         f2=p2T4(i  )*(t2tc(i  )*Tc)**4
1687         f3=p2T4(i+1)*(t2tc(i+1)*Tc)**4
1688         a=(t2tc(i  )*Tc - t2tc(i-1)*Tc)
1689         b=(t2tc(i+1)*Tc - t2tc(i  )*Tc)
1690         s=(f2-f1)/a*b/(a+b) + (f3-f2)/b*a/(a+b)
1691         s2T3=s / (t2tc(i)*Tc)**3
1692         e2T4(i)=  s2T3 - p2T4(i)
1693         epsi(i)=e2T4(i)*(t2tc(i)*Tc)**4*hc  ! in GeV/fm^3
1694       enddo
1695       end
1696 
1697 c------------------------------------------------------------------------------
1698       subroutine xFreezeOutTauX(meta)
1699 c------------------------------------------------------------------------------
1700       include 'epos.inc'
1701       include 'epos.inchy'
1702       common/cxyzt/xptl(mxptl),yptl(mxptl),zptl(mxptl),tptl(mxptl)
1703      *,optl(mxptl),uptl(mxptl),sptl(mxptl),rptl(mxptl,3)
1704       common/cen/ncentr /ctauhu/ntauhu(ncenthx,1-netahx:netahx-1)
1705       !..........................................................................
1706       nhis=1
1707       npl=0
1708       nplx=0
1709       deleta=etahy(2)-etahy(1)
1710       eta1=zetahy(meta)-deleta/2
1711       eta2=zetahy(meta)+deleta/2
1712       etaav=zetahy(meta)
1713       taumax=tauhoc(ncentr,ntauhu(ncentr,meta))+4
1714       do n=1,nptl
1715         if(ityptl(n).eq.60
1716      * .and.istptl(n).ne.12.and.istptl(n).ne.11)then
1717          if(istptl(iorptl(n)).eq.11)then
1718           tau=0
1719           tau2=xorptl(4,n)**2-xorptl(3,n)**2
1720           if(tau2.gt.0.)tau=sqrt(tau2)
1721           if(tau.lt.taumax)then
1722            rap=
1723      *     .5*alog((xorptl(4,n)+xorptl(3,n))/(xorptl(4,n)-xorptl(3,n)))
1724            if(rap.ge.eta1.and.rap.le.eta2)then
1725            if(abs(xorptl(2,n)).le.2)then
1726             npl=npl+1
1727             nplx=nplx+1
1728             if(npl.eq.1)then
1729               if(nplx.gt.1)
1730      *        write(ifhi,'(a)')       '  endarray closehisto plot 0-'
1731               write(ifhi,'(a)')       '!-------------------------------'
1732               write(ifhi,'(a)')       '!   tau-x      '
1733               write(ifhi,'(a)')       '!-------------------------------'
1734               write(ifhi,'(a)')       '!newpage'
1735               write(ifhi,'(a,i1)')    'openhisto name t-r-0-',nhis
1736               write(ifhi,'(a)')       'htyp prl xmod lin ymod lin'
1737               write(ifhi,'(a)')       'xrange -10 10'
1738               write(ifhi,'(a,f5.1)')  'yrange 0 ',taumax
1739               write(ifhi,'(a)')       'txt  "xaxis x (fm)"'
1740               write(ifhi,'(a)')       'txt  "yaxis [t] (fm/c)"'
1741               write(ifhi,'(a,f4.1,a)')'text 0.35 0.90 "[c]=',etaav,'"'
1742               write(ifhi,'(a)')    'text 0.02 0.9 """B#y"B#"L#2fm""'
1743           write(ifhi,'(a,f4.1,a)')'text 0.65 0.9 "',centhy(ncentr),'fm"'
1744               write(ifhi,'(a)')       'array 2'
1745             endif
1746             write(ifhi,'(2e11.3)')xorptl(1,n),tau
1747             if(npl.eq.1000)then
1748               nhis=nhis+1
1749               npl=0
1750             endif
1751            endif
1752            endif
1753           endif
1754          endif
1755         endif
1756       enddo
1757       if(nplx.gt.0)write(ifhi,'(a)')  '  endarray closehisto plot 0-'
1758     !..........................................................................
1759        nhis=1
1760       npl=0
1761       nplx=0
1762       npli=20
1763       do n=1,nptl
1764         if(dezptl(n).lt.1e3.and.n.le.maxfra
1765      *  .and.(istptl(n).eq.0.or.istptl(n).eq.1))then
1766            rap=1e10
1767            xp=tptl(n)+zptl(n)
1768            xm=tptl(n)-zptl(n)
1769            if(xp.gt.0.0.and.xm.gt.0.0)rap=.5*alog(xp/xm)
1770            if(rap.ge.-0.5.and.rap.le.0.5)then
1771            if(abs(xorptl(2,n)).le.2)then
1772             npl=npl+npli
1773             nplx=nplx+1
1774             if(npl.eq.npli)then
1775               if(nplx.gt.1)
1776      *        write(ifhi,'(a)')       '  endarray closehisto plot 0-'
1777               write(ifhi,'(a)')       '!-------------------------------'
1778               write(ifhi,'(a)')       '!   tau-x    corona  '
1779               write(ifhi,'(a)')       '!-------------------------------'
1780               write(ifhi,'(a,i1)')    'openhisto name t-r-0-cor-',nhis
1781               write(ifhi,'(a)')       'htyp pgl '
1782               write(ifhi,'(a)')       'array 2'
1783             endif
1784             x=xorptl(1,n)
1785             y=xorptl(2,n)
1786             z=xorptl(3,n)
1787             t=xorptl(4,n)
1788             tau=0
1789             tau2=t**2-z**2
1790             if(tau2.gt.0.)tau=sqrt(tau2)
1791             if(abs(y).le.2)
1792      .      write(ifhi,'(2e11.3)')x,tau
1793             dt=0.1
1794             do k=1,npli-1
1795             x=x+pptl(1,n)/pptl(4,n)*dt
1796             y=y+pptl(2,n)/pptl(4,n)*dt
1797             z=z+pptl(3,n)/pptl(4,n)*dt
1798             t=t+dt
1799             tau=0
1800             tau2=t**2-z**2
1801             if(tau2.gt.0.)tau=sqrt(tau2)
1802             if(abs(y).le.2)
1803      .            write(ifhi,'(2e11.3)')x,tau
1804             enddo
1805             if(npl.eq.1000)then
1806               nhis=nhis+1
1807               npl=0
1808             endif
1809            endif
1810            endif
1811         endif
1812       enddo
1813       write(ifhi,'(a)')  '  endarray closehisto plot 0'
1814     !..........................................................................
1815        end
1816 
1817 c------------------------------------------------------------------------------
1818       subroutine xFeff
1819 c------------------------------------------------------------------------------
1820       include 'epos.inc'
1821       include 'epos.inchy'
1822       common/cen/ncentr
1823       write(ifhi,'(a)') '!-------------------------------'
1824       write(ifhi,'(a)') '!   feff      '
1825       write(ifhi,'(a)') '!-------------------------------'
1826       write(ifhi,'(a)') '!newpage'
1827       write(ifhi,'(a)') 'openhisto name feff'
1828       write(ifhi,'(a)') 'htyp lin xmod lin ymod lin'
1829       write(ifhi,'(a)') 'xrange -6 6'
1830       write(ifhi,'(a)') 'yrange 0 1.2'
1831       write(ifhi,'(a)') 'txt  "xaxis [c] "'
1832       write(ifhi,'(a)') 'txt  "yaxis f?eff!"'
1833       write(ifhi,'(a,f4.1,a)') 'text 0.65 0.9 "',centhy(ncentr),'fm"'
1834       write(ifhi,'(a)') 'array 2'
1835       do meta=-netahy+1,netahy-1
1836         neta=abs(meta)+1
1837         eta=etahy(neta)
1838         if(meta.lt.0)eta=-eta
1839         write(ifhi,'(2e11.3)') eta,feff(ncentr,meta)
1840       enddo
1841       write(ifhi,'(a)')    '  endarray closehisto plot 0'
1842       end
1843 
1844 c------------------------------------------------------------------------------
1845       subroutine xFreezeOutTauEta
1846 c------------------------------------------------------------------------------
1847       include 'epos.inc'
1848       include 'epos.inchy'
1849       common/cxyzt/xptl(mxptl),yptl(mxptl),zptl(mxptl),tptl(mxptl)
1850      *,optl(mxptl),uptl(mxptl),sptl(mxptl),rptl(mxptl,3)
1851       common/cen/ncentr /ctauhu/ntauhu(ncenthx,1-netahx:netahx-1)
1852       taumax=tauhoc(ncentr,ntauhu(ncentr,0))+4
1853     !..........................................................................
1854       nhis=1
1855       npl=0
1856       do n=1,nptl
1857         if(ityptl(n).eq.60
1858      * .and.istptl(n).ne.12.and.istptl(n).ne.11)then
1859          if(istptl(iorptl(n)).eq.11)then
1860           tau=0
1861           tau2=xorptl(4,n)**2-xorptl(3,n)**2
1862           if(tau2.gt.0.)tau=sqrt(tau2)
1863           if(tau.lt.taumax)then
1864             npl=npl+1
1865             if(npl.eq.1)then
1866               write(ifhi,'(a)')      '!-------------------------------'
1867               write(ifhi,'(a)')      '!   tau-eta      '
1868               write(ifhi,'(a)')      '!-------------------------------'
1869               write(ifhi,'(a)')      '!newpage'
1870               write(ifhi,'(a,i1)')   'openhisto name t-eta-',nhis
1871               write(ifhi,'(a)')      'htyp prl xmod lin ymod lin'
1872               write(ifhi,'(a)')      'xrange -4 4'
1873               write(ifhi,'(a,f5.1)')  'yrange 0 ',taumax
1874               write(ifhi,'(a)')    'txt  "xaxis [c] "'
1875               write(ifhi,'(a)')    'txt  "yaxis [t] (fm/c)"'
1876       write(ifhi,'(a,f4.1,a)')'text 0.65 0.9 "',centhy(ncentr),'fm"'
1877               write(ifhi,'(a)')       'array 2'
1878             endif
1879             eta=
1880      *     .5*alog((xorptl(4,n)+xorptl(3,n))/(xorptl(4,n)-xorptl(3,n)))
1881             write(ifhi,'(2e11.3)') eta,tau
1882             if(npl.eq.1000)then
1883               write(ifhi,'(a)')    '  endarray closehisto plot 0-'
1884               nhis=nhis+1
1885               npl=0
1886             endif
1887           endif
1888          endif
1889         endif
1890       enddo
1891       if(npl.ne.0)write(ifhi,'(a)')  '  endarray closehisto plot 0'
1892       if(npl.eq.0)stop'xFreezeOutTZ: no particles!!!!!            '
1893       end
1894 
1895 c------------------------------------------------------------------------------
1896       subroutine xFreezeOutTZ
1897 c------------------------------------------------------------------------------
1898       include 'epos.inc'
1899       include 'epos.inchy'
1900       common/cxyzt/xptl(mxptl),yptl(mxptl),zptl(mxptl),tptl(mxptl)
1901      *,optl(mxptl),uptl(mxptl),sptl(mxptl),rptl(mxptl,3)
1902       common/cen/ncentr
1903     !..........................................................................
1904       nhis=1
1905       npl=0
1906       do n=1,nptl
1907         if(ityptl(n).eq.60
1908      * .and.istptl(n).ne.12.and.istptl(n).ne.11)then
1909          if(istptl(iorptl(n)).eq.11)then
1910             npl=npl+1
1911             if(npl.eq.1)then
1912               write(ifhi,'(a)')      '!-------------------------------'
1913               write(ifhi,'(a)')      '!   t-z      '
1914               write(ifhi,'(a)')      '!-------------------------------'
1915               write(ifhi,'(a)')      '!newpage'
1916               write(ifhi,'(a,i1)')   'openhisto name t-z-',nhis
1917               write(ifhi,'(a)')      'htyp prl xmod lin ymod lin'
1918               write(ifhi,'(a)')      'xrange -25 25'
1919               write(ifhi,'(a)')      'yrange 0 25 '
1920               write(ifhi,'(a)')    'txt  "xaxis z (fm)"'
1921               write(ifhi,'(a)')    'txt  "yaxis t (fm/c)"'
1922       write(ifhi,'(a,f4.1,a)')'text 0.70 0.22 "',centhy(ncentr),'fm"'
1923               write(ifhi,'(a)')       'array 2'
1924             endif
1925             write(ifhi,'(2e11.3)') xorptl(3,n),xorptl(4,n)
1926             if(npl.eq.1000)then
1927               write(ifhi,'(a)')    '  endarray closehisto plot 0-'
1928               nhis=nhis+1
1929               npl=0
1930             endif
1931          endif
1932         endif
1933       enddo
1934       if(npl.ne.0)write(ifhi,'(a)')  '  endarray closehisto plot 0'
1935       if(npl.eq.0)stop'xFreezeOutTZ: no particles!!!!!            '
1936       end
1937 
1938 c------------------------------------------------------------------------------
1939       subroutine xFoMass(meta)
1940 c------------------------------------------------------------------------------
1941       include 'epos.inc'
1942       include 'epos.inchy'
1943       common/cen/ncentr /ctauhu/ntauhu(ncenthx,1-netahx:netahx-1)
1944       taumax=tauhoc(ncentr,ntauhu(ncentr,meta))+2
1945 c      netahyxx=netahy
1946       write(ifhi,'(a)')    '!----------------------------------------'
1947       write(ifhi,'(a,i3)') '!   hydro freeze out rate     '
1948       write(ifhi,'(a)')    '!----------------------------------------'
1949       write(ifhi,'(a)')       '!newpage'
1950       do ii=1,3
1951         if(ii.eq.1)nphi=2
1952         if(ii.eq.2)nphi=1+nphihy/4
1953         if(ii.eq.3)nphi=1+nphihy/2
1954         write(ifhi,'(a,3i1)')'openhisto htyp lin name w-',meta,nphi,ii
1955         if(ii.eq.1)then !----------------------
1956          write(ifhi,'(a,f4.1)')'xmod lin xrange 0. ',taumax
1957          write(ifhi,'(a)')    'txt  "xaxis [t] (fm/c)"'
1958          write(ifhi,'(a)') 'ymod lin yrange auto auto '
1959          write(ifhi,'(a,f4.1,a)')'text 0.1 0.9 "  [c]=',zetahy(meta),'"'
1960          write(ifhi,'(a)')'txt "yaxis w "'
1961          write(ifhi,'(a,f4.1,a)')'text 0.65 0.9 "',centhy(ncentr),'fm"'
1962         endif       !-------------------------------
1963         write(ifhi,'(a)')'array 2'
1964         deltau=tauhoc(ncentr,2)-tauhoc(ncentr,1)
1965         do ntau=2,ntauhu(ncentr,meta)
1966          dy=waa(ncentr,meta,ntau,nphi)
1967          write(ifhi,'(2e13.5)')tauhoc(ncentr,ntau)-deltau/2,dy
1968         enddo
1969         write(ifhi,'(a)') 'endarray closehisto '
1970         if(ii.ne.3)write(ifhi,'(a)') 'plot 0-'
1971         if(ii.eq.3)write(ifhi,'(a)') 'plot 0'
1972       enddo
1973       end
1974 
1975 c------------------------------------------------------------------------------
1976       subroutine xFoRadius(meta)
1977 c------------------------------------------------------------------------------
1978       include 'epos.inc'
1979       include 'epos.inchy'
1980       common/cen/ncentr /ctauhu/ntauhu(ncenthx,1-netahx:netahx-1)
1981       taumax=tauhoc(ncentr,ntauhu(ncentr,meta))+2
1982 c      netahyxx=netahy
1983       write(ifhi,'(a)')    '!----------------------------------------'
1984       write(ifhi,'(a,i3)') '!        hydro freeze out radius     '
1985       write(ifhi,'(a)')    '!----------------------------------------'
1986       write(ifhi,'(a)')    '!newpage'
1987       do ii=1,3
1988         if(ii.eq.1)nphi=2
1989         if(ii.eq.2)nphi=1+nphihy/4
1990         if(ii.eq.3)nphi=1+nphihy/2
1991         write(ifhi,'(a,3i1)')'openhisto htyp lin name r-',meta,nphi,ii
1992         if(ii.eq.1)then !----------------------
1993          write(ifhi,'(a,f4.1)')'xmod lin xrange 0. ',taumax
1994          write(ifhi,'(a)')'txt  "xaxis [t] (fm/c)"'
1995          write(ifhi,'(a)') 'ymod lin yrange auto auto '
1996          write(ifhi,'(a,f4.1,a)')'text 0.1 0.9 "  [c]=',zetahy(meta),'"'
1997          write(ifhi,'(a)')'txt "yaxis r  (fm) "'
1998          write(ifhi,'(a,f4.1,a)')'text 0.65 0.9 "',centhy(ncentr),'fm"'
1999         endif !-------------------------------
2000         write(ifhi,'(a)')'array 2'
2001         do ntau=2,ntauhu(ncentr,meta)
2002          write(ifhi,'(2e13.5)')tauhoc(ncentr,ntau)
2003      .     ,raa(ncentr,meta,ntau,nphi)
2004          enddo
2005         write(ifhi,'(a)') 'endarray closehisto '
2006         if(ii.ne.3)write(ifhi,'(a)') 'plot 0-'
2007         if(ii.eq.3)write(ifhi,'(a)') 'plot 0'
2008       enddo
2009       end
2010 
2011 c------------------------------------------------------------------------------
2012       subroutine xFoRadVelocity(meta)
2013 c------------------------------------------------------------------------------
2014       include 'epos.inc'
2015       include 'epos.inchy'
2016       common/cen/ncentr /ctauhu/ntauhu(ncenthx,1-netahx:netahx-1)
2017       taumax=tauhoc(ncentr,ntauhu(ncentr,meta))+2
2018 c      netahyxx=netahy
2019       write(ifhi,'(a)')    '!----------------------------------------'
2020       write(ifhi,'(a,i3)') '!     hydro freeze out rad velocity      '
2021       write(ifhi,'(a)')    '!----------------------------------------'
2022       write(ifhi,'(a)')    '!newpage'
2023       do ii=1,3
2024         if(ii.eq.1)nphi=2
2025         if(ii.eq.2)nphi=1+nphihy/4
2026         if(ii.eq.3)nphi=1+nphihy/2
2027         write(ifhi,'(a,3i1)')'openhisto htyp lin name y-',meta,nphi,ii
2028         if(ii.eq.1)then !----------------------
2029          write(ifhi,'(a,f4.1)')'xmod lin xrange 0. ',taumax
2030          write(ifhi,'(a)')'txt "xaxis [t] (fm/c)"'
2031          write(ifhi,'(a)') 'ymod lin yrange auto auto '
2032          write(ifhi,'(a,f4.1,a)')'text 0.1 0.9 "  [c]=',zetahy(meta),'"'
2033          write(ifhi,'(a)')'txt "yaxis v?rad!  "'
2034          write(ifhi,'(a,f4.1,a)')'text 0.65 0.9 "',centhy(ncentr),'fm"'
2035         endif !-------------------------------
2036         write(ifhi,'(a)')'array 2'
2037         do ntau=2,ntauhu(ncentr,meta)
2038          write(ifhi,'(2e13.5)')tauhoc(ncentr,ntau)
2039      .     ,vaa(1,ncentr,meta,ntau,nphi)
2040         enddo
2041         write(ifhi,'(a)') 'endarray closehisto '
2042         if(ii.ne.3)write(ifhi,'(a)') 'plot 0-'
2043         if(ii.eq.3)write(ifhi,'(a)') 'plot 0'
2044       enddo
2045       end
2046 c------------------------------------------------------------------------------
2047       subroutine xFoTanVelocity(meta)
2048 c------------------------------------------------------------------------------
2049       include 'epos.inc'
2050       include 'epos.inchy'
2051       common/cen/ncentr /ctauhu/ntauhu(ncenthx,1-netahx:netahx-1)
2052       taumax=tauhoc(ncentr,ntauhu(ncentr,meta))+2
2053 c      netahyxx=netahy
2054       write(ifhi,'(a)')    '!----------------------------------------'
2055       write(ifhi,'(a,i3)') '!    hydro freeze out tang velocity      '
2056       write(ifhi,'(a)')    '!----------------------------------------'
2057       write(ifhi,'(a)')    '!newpage'
2058       do ii=1,4
2059         if(ii.eq.1)nphi=2
2060         if(ii.eq.2)nphi=1+nphihy/8
2061         if(ii.eq.3)nphi=1+nphihy/8*3
2062         if(ii.eq.4)nphi=1+nphihy/2
2063         write(ifhi,'(a,3i1)')'openhisto htyp lin name y-',meta,nphi,ii
2064         if(ii.eq.1)then !----------------------
2065          write(ifhi,'(a,f4.1)')'xmod lin xrange 0. ',taumax
2066          write(ifhi,'(a)')'txt "xaxis [t] (fm/c)"'
2067          write(ifhi,'(a)') 'ymod lin yrange auto auto '
2068          write(ifhi,'(a,f4.1,a)')'text 0.1 0.9 "  [c]=',zetahy(meta),'"'
2069          write(ifhi,'(a)')'txt "yaxis v?tan!  "'
2070           write(ifhi,'(a,f4.1,a)')'text 0.65 0.9 "',centhy(ncentr),'fm"'
2071         endif !-------------------------------
2072         write(ifhi,'(a)')'array 2'
2073         do ntau=2,ntauhu(ncentr,meta)
2074          write(ifhi,'(2e13.5)')tauhoc(ncentr,ntau)
2075      .     ,vaa(2,ncentr,meta,ntau,nphi)
2076         enddo
2077         write(ifhi,'(a)') 'endarray closehisto '
2078         if(ii.ne.4)write(ifhi,'(a)') 'plot 0-'
2079         if(ii.eq.4)write(ifhi,'(a)') 'plot 0'
2080       enddo
2081       end
2082 
2083 c-----------------------------------------------------------------------
2084       subroutine xCoreCorona(iii,jjj)
2085 c-----------------------------------------------------------------------
2086 c     space-time evolution of core and corona
2087 c
2088 c     cluster ............   ist=11  ity=60
2089 c     core particles .....   ist=0   ity=60
2090 c     corona particles ...   ist=0   ity/=60
2091 c
2092 c    iii=1: plot also binary collisions
2093 c    jjj>0: multiplicity trigger (useful for pp)
2094 c-----------------------------------------------------------------------
2095       include 'epos.inc'
2096       include 'epos.incems'
2097       include 'epos.inchy'
2098       common/cen/ncentr
2099       common/cxyzt/xptl(mxptl),yptl(mxptl),zptl(mxptl),tptl(mxptl)
2100      *,optl(mxptl),uptl(mxptl),sptl(mxptl),rptl(mxptl,3)
2101       common/cdelzet/delzet,delsgr
2102       parameter (myy=48,mrr=100)
2103       real yy(myy),rr(mrr)
2104       common/cranphi/ranphi
2105 
2106       phinll=phievt+ranphi
2107       phi=   phievt+ranphi
2108       !print*,'   EventPhi=',phievt,'   RandomPhi=',ranphi
2109       rapmax=6
2110       radmax=10
2111       r1=0.0
2112       if(maproj.gt.1)r1=radnuc(maproj)
2113       r2=0.0
2114       if(matarg.gt.1)r2=radnuc(matarg)
2115       if(maproj.eq.1.and.matarg.gt.1)r1=r2
2116       if(maproj.gt.1.and.matarg.eq.1)r2=r1
2117       a=7.8
2118       b=bimevt/2
2119 c      n1=koievt
2120       n2=0
2121       do k=1,koll
2122        if(itpr(k).gt.0)n2=n2+1
2123       enddo
2124 c      n3=nglevt
2125 
2126       if(jjj.gt.0)then
2127       multy1=0
2128        do i=maproj+matarg+1,nptl
2129         if(istptl(i).eq.0)then
2130           amt=pptl(5,i)**2+pptl(1,i)**2+pptl(2,i)**2
2131           rap=1000
2132           if(amt.gt.0..and.pptl(4,i).gt.0.)then
2133             amt=sqrt(amt)
2134             rap=sign(1.,pptl(3,i))*alog((pptl(4,i)+abs(pptl(3,i)))/amt)
2135           endif
2136           ch=0
2137           if(abs(idptl(i)).ge.100.and.abs(idptl(i)).lt.10000)then
2138             call idchrg(idptl(i),ch)
2139             if(abs(ch).gt.0.1.and.abs(rap).le.1.)multy1=multy1+1
2140           endif
2141          endif
2142         enddo
2143         ih1=jjj/100
2144         ih2=mod(jjj,100)
2145         if(0.5*multy1.lt.ih1.or.0.5*multy1.gt.ih2)return
2146       endif
2147       xmax=1+int(a*1.5)
2148 
2149       write(ifhi,'(a)')       '!---------------------------------'
2150       write(ifhi,'(a)')       '!   core particles                '
2151       write(ifhi,'(a)')       '!---------------------------------'
2152       write(ifhi,'(a)')       '!newpage'
2153       write(ifhi,'(a)')      'openhisto name st1'
2154       write(ifhi,'(a)')       'htyp prv xmod lin ymod lin'
2155       write(ifhi,'(a,2e11.3)')'xrange',-xmax,xmax
2156       write(ifhi,'(a,2e11.3)')'yrange',-a,a
2157       write(ifhi,'(a)')    'text 0.85 -0.25 " x (fm)"'
2158       write(ifhi,'(a)')    'txt "yaxis y (fm)"'
2159       write(ifhi,'(a)')       'text 0.05 0.90 "core"'
2160       write(ifhi,'(a)')       'text 0.05 0.80 "corona"'
2161       write(ifhi,'(a,f4.1,a)')'text 0.82 0.07 "[c]=0"'
2162       write(ifhi,'(a,f4.1,a)')'text 0.65 0.9 "',centhy(ncentr),'fm"'
2163       write(ifhi,'(a)')       'array 2'
2164       ncore=0
2165       do i=1,nptl
2166        if(abs(sptl(i)).lt.0.5*delsgr.and.istptl(i).eq.2)then
2167          write(ifhi,'(2e11.3)')xptl(i)*cos(phi)+yptl(i)*sin(phi)
2168      .      ,                      -xptl(i)*sin(phi)+yptl(i)*cos(phi)
2169          ncore=ncore+1
2170        endif
2171       enddo
2172       write(ifhi,'(a)')    '  endarray'
2173       write(ifhi,'(a)')    'closehisto plot 0-'
2174       write(ifhi,'(a)')       '!----------------------------------'
2175       write(ifhi,'(a)')       '!   corona particles               '
2176       write(ifhi,'(a)')       '!----------------------------------'
2177       write(ifhi,'(a)')      'openhisto name st2'
2178       write(ifhi,'(a)')       'htyp pgv xmod lin ymod lin'
2179       write(ifhi,'(a)')       'array 2'
2180       ncorona=0
2181       do i=1,nptl
2182        if(abs(sptl(i)).lt.0.5*delsgr.and.istptl(i).eq.0
2183      .  .and.ityptl(i).ne.60.and.ityptl(i).ne.19)then
2184          write(ifhi,'(2e11.3)')xptl(i)*cos(phi)+yptl(i)*sin(phi)
2185      .      ,                      -xptl(i)*sin(phi)+yptl(i)*cos(phi)
2186         ncorona=ncorona+1
2187        endif
2188       enddo
2189       write(ifhi,'(a)')    '  endarray'
2190       write(ifhi,'(a)')    'closehisto'
2191       !print*,'b=',bimevt,'   ncorona:ncore =  ',ncorona,':',ncore
2192          if(iii.eq.1)then
2193       write(ifhi,'(a)')       '!----------------------------------'
2194       write(ifhi,'(a)')       '!   binary collisions             '
2195       write(ifhi,'(a)')       '!----------------------------------'
2196       write(ifhi,'(a)')    'plot 0-'
2197       write(ifhi,'(a)')      'openhisto name coo'
2198       write(ifhi,'(a)')       'htyp pbl xmod lin ymod lin'
2199       write(ifhi,'(a,2e11.3)')'xrange',-xmax,xmax
2200       write(ifhi,'(a,2e11.3)')'yrange',-a,a
2201       write(ifhi,'(a)')       'array 2'
2202       do k=1,koll
2203        if(itpr(k).gt.0)then
2204          write(ifhi,'(2e11.3)')coord(1,k)*cos(phi)+coord(2,k)*sin(phi)
2205      * ,                      -coord(1,k)*sin(phi)+coord(2,k)*cos(phi)
2206        endif
2207       enddo
2208       write(ifhi,'(a)')    '  endarray'
2209       write(ifhi,'(a)')    'closehisto'
2210          endif
2211       if(r1.ne.0.0)then
2212         write(ifhi,'(a)')    'plot 0-'
2213         write(ifhi,'(a)')       '!----------------------------------'
2214         write(ifhi,'(a)')       '!   hard spheres             '
2215         write(ifhi,'(a)')       '!----------------------------------'
2216         write(ifhi,'(a)')   'openhisto name stc1 htyp lyu'
2217         write(ifhi,'(a)')   'array 2'
2218         do j=-50,50
2219          phi=j/50.*pi*0.55
2220          write(ifhi,'(2e11.3)')r1*cos(phi)-b,r1*sin(phi)
2221         enddo
2222         write(ifhi,'(a)')    '  endarray'
2223         write(ifhi,'(a)')    'closehisto'
2224       endif
2225       if(r2.ne.0.0)then
2226         write(ifhi,'(a)')    'plot 0-'
2227         write(ifhi,'(a)')   'openhisto name stc1 htyp lyu'
2228         write(ifhi,'(a)')   'array 2'
2229         do j=-50,50
2230          phi=j/50.*pi*0.55
2231          write(ifhi,'(2e11.3)')-r1*cos(phi)+b,r1*sin(phi)
2232         enddo
2233         write(ifhi,'(a)')    '  endarray'
2234         write(ifhi,'(a)')    'closehisto'
2235       endif
2236 
2237       write(ifhi,'(a)')    'plot 0'
2238 
2239    !........................................................................................
2240       if(ioclude.le.1)return
2241    !........................................................................................
2242       delrap=2*rapmax/float(myy)
2243       do m=1,myy
2244         yy(m)=0
2245       enddo
2246       do n=1,nptl
2247         if(dezptl(n).lt.1e3.and.n.le.maxfra.and.istptl(n).eq.2)then
2248         routp=-sin(phinll)*xptl(n)+cos(phinll)*yptl(n)
2249         rinp = cos(phinll)*xptl(n)+sin(phinll)*yptl(n)
2250         if(abs(rinp).le.1.and.abs(routp).le.1.)then
2251           rapx=dezptl(n)
2252           eco=0
2253           amt=pptl(5,n)**2+pptl(1,n)**2+pptl(2,n)**2
2254           if(amt.gt.0..and.pptl(4,n)+abs(pptl(3,n)).gt.0.d0)then
2255             amt=sqrt(amt)
2256             rap=sign(1.,pptl(3,n))*alog((pptl(4,n)+abs(pptl(3,n)))/amt)
2257             eco=amt*cosh(rap-rapx)
2258           endif
2259           m=(rapx+rapmax)/delrap+1
2260           if(m.gt.myy)m=myy
2261           yy(m)=yy(m)+eco
2262         endif
2263         endif
2264       enddo
2265       write(ifhi,'(a)')'!---------------------------------------------'
2266       write(ifhi,'(a)')'!    core segment energy per d[c]dxdy         '
2267       write(ifhi,'(a)')'!           vs space-time rapidity rapx       '
2268       write(ifhi,'(a)')'!   (same as histogram rapx eco... in optns)  '
2269       write(ifhi,'(a)')'!---------------------------------------------'
2270       write(ifhi,'(a)')       '!newpage'
2271       write(ifhi,'(a)')    'openhisto name rapx'
2272       write(ifhi,'(a)')       'htyp lin xmod lin ymod lin'
2273       write(ifhi,'(a,2f7.3)') 'xrange ',-rapmax,rapmax
2274       write(ifhi,'(a)')        'yrange 0 auto '
2275       write(ifhi,'(a)')  'txt "title initial energy          "'
2276       write(ifhi,'(a,f4.1,a)')'text 0.05 0.70 "x=0"'
2277       write(ifhi,'(a,f4.1,a)')'text 0.05 0.60 "y=0"'
2278       write(ifhi,'(a,f4.1,a)')'text 0.65 0.9 "',centhy(ncentr),'fm"'
2279       write(ifhi,'(a)')  'txt  "xaxis space-time rapidity [c] "'
2280       write(ifhi,'(a)')  'txt  "yaxis dE/d[c]dxdy "'
2281       write(ifhi,'(a)')       'array 2'
2282       do m=1,myy
2283         write(ifhi,'(2e11.3)')-rapmax+(m-0.5)*delrap, yy(m)/4./delrap
2284       enddo
2285       write(ifhi,'(a)')  '  endarray closehisto plot 0-'
2286    !........................................................................................
2287       delrap=2*rapmax/float(myy)
2288       do m=1,myy
2289         yy(m)=0
2290       enddo
2291       do n=1,nptl
2292         if(dezptl(n).lt.1e3.and.n.le.maxfra.and.istptl(n)/2.eq.0)then
2293         routp=-sin(phinll)*xptl(n)+cos(phinll)*yptl(n)
2294         rinp = cos(phinll)*xptl(n)+sin(phinll)*yptl(n)
2295         if(abs(rinp).le.1.and.abs(routp).le.1.)then
2296           rapx=dezptl(n)
2297           eco=0
2298           amt=pptl(5,n)**2+pptl(1,n)**2+pptl(2,n)**2
2299           if(amt.gt.0..and.pptl(4,n)+abs(pptl(3,n)).gt.0.d0)then
2300             amt=sqrt(amt)
2301             rap=sign(1.,pptl(3,n))*alog((pptl(4,n)+abs(pptl(3,n)))/amt)
2302             eco=amt*cosh(rap-rapx)
2303           endif
2304           m=(rapx+rapmax)/delrap+1
2305           if(m.ge.1.and.m.le.myy)yy(m)=yy(m)+eco
2306         endif
2307         endif
2308       enddo
2309       write(ifhi,'(a)')'!---------------------------------------------'
2310       write(ifhi,'(a)')'!    corona segment energy per d[c]dxdy       '
2311       write(ifhi,'(a)')'!           vs space-time rapidity rapx       '
2312       write(ifhi,'(a)')'!   (same as histogram rapx eco... in optns)  '
2313       write(ifhi,'(a)')'!---------------------------------------------'
2314       write(ifhi,'(a)')    'openhisto name rapx'
2315       write(ifhi,'(a)')       'htyp lin xmod lin ymod lin'
2316       write(ifhi,'(a,2f7.3)') 'xrange ',-rapmax,rapmax
2317       write(ifhi,'(a)')        'yrange 0 auto '
2318       write(ifhi,'(a)')  'txt "title initial energy          "'
2319       write(ifhi,'(a)')  'txt  "xaxis space-time rapidity [c] "'
2320       write(ifhi,'(a)')  'txt  "yaxis dE/d[c]dxdy "'
2321       write(ifhi,'(a)')       'array 2'
2322       do m=1,myy
2323         write(ifhi,'(2e11.3)')-rapmax+(m-0.5)*delrap, yy(m)/4./delrap
2324       enddo
2325       write(ifhi,'(a)')  '  endarray closehisto plot 0-'
2326       call xEiniEta(1)
2327       write(ifhi,'(a)')  'plot 0'
2328    !........................................................................................
2329       delrad=2*radmax/float(mrr)
2330       do m=1,mrr
2331         rr(m)=0
2332       enddo
2333       do n=1,nptl
2334         if(dezptl(n).lt.1e3.and.n.le.maxfra.and.istptl(n).eq.2)then
2335         routp=-sin(phinll)*xptl(n)+cos(phinll)*yptl(n)
2336         rinp = cos(phinll)*xptl(n)+sin(phinll)*yptl(n)
2337         rapx=dezptl(n)
2338         if(abs(rapx).le.1.and.abs(routp).le.1.)then
2339           eco=0
2340           amt=pptl(5,n)**2+pptl(1,n)**2+pptl(2,n)**2
2341           if(amt.gt.0..and.pptl(4,n)+abs(pptl(3,n)).gt.0.d0)then
2342             amt=sqrt(amt)
2343             rap=sign(1.,pptl(3,n))*alog((pptl(4,n)+abs(pptl(3,n)))/amt)
2344             eco=amt*cosh(rap-rapx)
2345           endif
2346           m=(rinp+radmax)/delrad+1
2347           if(m.gt.mrr)m=mrr
2348           rr(m)=rr(m)+eco
2349         endif
2350         endif
2351       enddo
2352       write(ifhi,'(a)')'!---------------------------------------------'
2353       write(ifhi,'(a)')'!    core segment energy per d[c]dxdy vs x    '
2354       write(ifhi,'(a)')'!   (same as histogram rinp eco... in optns)  '
2355       write(ifhi,'(a)')'!---------------------------------------------'
2356       write(ifhi,'(a)')       '!newpage'
2357       write(ifhi,'(a)')    'openhisto name rapx'
2358       write(ifhi,'(a)')       'htyp lin xmod lin ymod lin'
2359       write(ifhi,'(a,2f7.3)') 'xrange ',-radmax,radmax
2360       write(ifhi,'(a)')        'yrange 0 auto '
2361       write(ifhi,'(a)')  'txt "title initial energy          "'
2362       write(ifhi,'(a,f4.1,a)')'text 0.05 0.70 "[c]=0"'
2363       write(ifhi,'(a,f4.1,a)')'text 0.05 0.60 "y=0"'
2364       write(ifhi,'(a,f4.1,a)')'text 0.65 0.9 "',centhy(ncentr),'fm"'
2365       write(ifhi,'(a)')  'txt  "xaxis x (fm)"'
2366       write(ifhi,'(a)')  'txt  "yaxis dE/d[c]dxdy "'
2367       write(ifhi,'(a)')       'array 2'
2368       do m=1,mrr
2369         write(ifhi,'(2e11.3)')-radmax+(m-0.5)*delrad, rr(m)/4./delrad
2370       enddo
2371       write(ifhi,'(a)')  '  endarray closehisto plot 0-'
2372    !........................................................................................
2373       delrad=2*radmax/float(mrr)
2374       do m=1,mrr
2375         rr(m)=0
2376       enddo
2377       do n=1,nptl
2378         if(dezptl(n).lt.1e3.and.n.le.maxfra.and.istptl(n)/2.eq.0)then
2379         routp=-sin(phinll)*xptl(n)+cos(phinll)*yptl(n)
2380         rinp = cos(phinll)*xptl(n)+sin(phinll)*yptl(n)
2381         rapx=dezptl(n)
2382         if(abs(rapx).le.1.and.abs(routp).le.1.)then
2383           eco=0
2384           amt=pptl(5,n)**2+pptl(1,n)**2+pptl(2,n)**2
2385           if(amt.gt.0..and.pptl(4,n)+abs(pptl(3,n)).gt.0.d0)then
2386             amt=sqrt(amt)
2387             rap=sign(1.,pptl(3,n))*alog((pptl(4,n)+abs(pptl(3,n)))/amt)
2388             eco=amt*cosh(rap-rapx)
2389           endif
2390           m=(rinp+radmax)/delrad+1
2391           if(m.gt.mrr)m=mrr
2392           rr(m)=rr(m)+eco
2393         endif
2394         endif
2395       enddo
2396       write(ifhi,'(a)')'!---------------------------------------------'
2397       write(ifhi,'(a)')'!    corona segment energy per d[c]dxdy vs x  '
2398       write(ifhi,'(a)')'!   (same as histogram rinp eco... in optns)  '
2399       write(ifhi,'(a)')'!---------------------------------------------'
2400       write(ifhi,'(a)')    'openhisto name rapx'
2401       write(ifhi,'(a)')       'htyp lin xmod lin ymod lin'
2402       write(ifhi,'(a)')       'array 2'
2403       do m=1,mrr
2404         write(ifhi,'(2e11.3)')-radmax+(m-0.5)*delrad, rr(m)/4./delrad
2405       enddo
2406       write(ifhi,'(a)')  '  endarray closehisto plot 0-'
2407       call xEiniX(1)
2408       write(ifhi,'(a)')  'plot 0'
2409    !........................................................................................
2410       delrad=2*radmax/float(mrr)
2411       do m=1,mrr
2412         rr(m)=0
2413       enddo
2414       do n=1,nptl
2415         if(dezptl(n).lt.1e3.and.n.le.maxfra.and.istptl(n).eq.2)then
2416         routp=-sin(phinll)*xptl(n)+cos(phinll)*yptl(n)
2417         rinp = cos(phinll)*xptl(n)+sin(phinll)*yptl(n)
2418         rapx=dezptl(n)
2419         if(abs(rapx).le.1.and.abs(rinp).le.1.)then
2420           eco=0
2421           amt=pptl(5,n)**2+pptl(1,n)**2+pptl(2,n)**2
2422           if(amt.gt.0..and.pptl(4,n)+abs(pptl(3,n)).gt.0.d0)then
2423             amt=sqrt(amt)
2424             rap=sign(1.,pptl(3,n))*alog((pptl(4,n)+abs(pptl(3,n)))/amt)
2425             eco=amt*cosh(rap-rapx)
2426           endif
2427           m=(routp+radmax)/delrad+1
2428           if(m.gt.mrr)m=mrr
2429           rr(m)=rr(m)+eco
2430         endif
2431         endif
2432       enddo
2433       write(ifhi,'(a)')'!---------------------------------------------'
2434       write(ifhi,'(a)')'!    core segment energy per d[c]dxdy vs y    '
2435       write(ifhi,'(a)')'!   (same as histogram routp eco... in optns)  '
2436       write(ifhi,'(a)')'!---------------------------------------------'
2437       write(ifhi,'(a)')       '!newpage'
2438       write(ifhi,'(a)')    'openhisto name rout'
2439       write(ifhi,'(a)')       'htyp lin xmod lin ymod lin'
2440       write(ifhi,'(a,2f7.3)') 'xrange ',-radmax,radmax
2441       write(ifhi,'(a)')        'yrange 0 auto '
2442       write(ifhi,'(a)')  'txt "title initial energy          "'
2443       write(ifhi,'(a,f4.1,a)')'text 0.05 0.70 "[c]=0"'
2444       write(ifhi,'(a,f4.1,a)')'text 0.05 0.60 "x=0"'
2445       write(ifhi,'(a,f4.1,a)')'text 0.65 0.9 "',centhy(ncentr),'fm"'
2446       write(ifhi,'(a)')  'txt  "xaxis y (fm)"'
2447       write(ifhi,'(a)')  'txt  "yaxis dE/d[c]dxdy "'
2448       write(ifhi,'(a)')       'array 2'
2449       do m=1,mrr
2450         write(ifhi,'(2e11.3)')-radmax+(m-0.5)*delrad, rr(m)/4./delrad
2451       enddo
2452       write(ifhi,'(a)')  '  endarray closehisto plot 0-'
2453    !........................................................................................
2454       delrad=2*radmax/float(mrr)
2455       do m=1,mrr
2456         rr(m)=0
2457       enddo
2458       do n=1,nptl
2459         if(dezptl(n).lt.1e3.and.n.le.maxfra.and.istptl(n)/2.eq.0)then
2460         routp=-sin(phinll)*xptl(n)+cos(phinll)*yptl(n)
2461         rinp = cos(phinll)*xptl(n)+sin(phinll)*yptl(n)
2462         rapx=dezptl(n)
2463         if(abs(rapx).le.1.and.abs(rinp).le.1.)then
2464           eco=0
2465           amt=pptl(5,n)**2+pptl(1,n)**2+pptl(2,n)**2
2466           if(amt.gt.0..and.pptl(4,n)+abs(pptl(3,n)).gt.0.d0)then
2467             amt=sqrt(amt)
2468             rap=sign(1.,pptl(3,n))*alog((pptl(4,n)+abs(pptl(3,n)))/amt)
2469             eco=amt*cosh(rap-rapx)
2470           endif
2471           m=(routp+radmax)/delrad+1
2472           if(m.gt.mrr)m=mrr
2473           rr(m)=rr(m)+eco
2474         endif
2475         endif
2476       enddo
2477       write(ifhi,'(a)')'!---------------------------------------------'
2478       write(ifhi,'(a)')'!    corona segment energy per d[c]dxdy vs y  '
2479       write(ifhi,'(a)')'!   (same as histogram routp eco... in optns)  '
2480       write(ifhi,'(a)')'!---------------------------------------------'
2481       write(ifhi,'(a)')    'openhisto name rout'
2482       write(ifhi,'(a)')       'htyp lin xmod lin ymod lin'
2483       write(ifhi,'(a)')       'array 2'
2484       do m=1,mrr
2485         write(ifhi,'(2e11.3)')-radmax+(m-0.5)*delrad, rr(m)/4./delrad
2486       enddo
2487       write(ifhi,'(a)')  '  endarray closehisto plot 0-'
2488       call xEiniY(1)
2489       write(ifhi,'(a)')  'plot 0'
2490 
2491       end
2492 
2493 c------------------------------------------------------------------------------
2494       subroutine xEini(ii)
2495 c------------------------------------------------------------------------------
2496       include 'epos.inc'
2497       include 'epos.inchy'
2498       common/cen/ncentr
2499 
2500       entry xEiniX(ii)
2501       write(ifhi,'(a)')    '!----------------------------------------'
2502       write(ifhi,'(a,i3)') '!        hydro initial energy vs x  '
2503       write(ifhi,'(a)')    '!----------------------------------------'
2504       write(ifhi,'(a)')'openhisto  array 2'
2505       do nr=nradhy,2,-1
2506        write(ifhi,'(2e13.5)')
2507      .     -radhy(nr),epsii(ncentr,ii,1,nr)*tauhoc(ncentr,1)
2508       enddo
2509       do nr=1,nradhy
2510        write(ifhi,'(2e13.5)')
2511      .     radhy(nr),epsii(ncentr,ii,1,nr)*tauhoc(ncentr,1)
2512       enddo
2513       write(ifhi,'(a)') 'endarray closehisto '
2514       return
2515 
2516       entry xEiniY(ii)
2517       write(ifhi,'(a)')    '!----------------------------------------'
2518       write(ifhi,'(a,i3)') '!        hydro initial energy vs y     '
2519       write(ifhi,'(a)')    '!----------------------------------------'
2520       write(ifhi,'(a)')'openhisto  array 2'
2521       do nr=nradhy,2,-1
2522        write(ifhi,'(2e13.5)')
2523      .    -radhy(nr),epsii(ncentr,ii,nphihy,nr)*tauhoc(ncentr,1)
2524       enddo
2525       do nr=1,nradhy
2526        write(ifhi,'(2e13.5)')
2527      .    radhy(nr),epsii(ncentr,ii,nphihy,nr)*tauhoc(ncentr,1)
2528       enddo
2529       write(ifhi,'(a)') 'endarray closehisto '
2530       return
2531 
2532       entry xEiniEta(ii)
2533       write(ifhi,'(a)')    '!----------------------------------------'
2534       write(ifhi,'(a,i3)') '!        hydro initial energy vs y     '
2535       write(ifhi,'(a)')    '!----------------------------------------'
2536       write(ifhi,'(a)')'openhisto  array 2'
2537       do neta=netahy,2,-1
2538        write(ifhi,'(2e13.5)')
2539      .   -etahy(neta),epsii(ncentr,neta,1,1)*tauhoc(ncentr,1)
2540       enddo
2541       do neta=1,netahy
2542        write(ifhi,'(2e13.5)')
2543      .   etahy(neta),epsii(ncentr,neta,1,1)*tauhoc(ncentr,1)
2544       enddo
2545       write(ifhi,'(a)') 'endarray closehisto '
2546       return
2547 
2548       end
2549 
2550 c-----------------------------------------------------------------------------
2551       integer function idxHiranoTable(id)
2552 c-----------------------------------------------------------------------------
2553       ida=abs(id)
2554       ihi=0
2555       if(ida.eq.120.or.id.eq.110)   ihi=2       !pion
2556       if(id.eq.220)                 ihi=3       !eta
2557       if(ida.eq.121.or.id.eq.111)   ihi=4       !rho
2558       if(id.eq.221)                 ihi=5       !omega
2559       !if(id.eq.)                        ihi=6         !sigma
2560       if(id.eq.330)                 ihi=7       !eta prime
2561       !if(id.eq.)                        ihi=8         !f_0
2562       !if(id.eq.)                        ihi=9         !a_0
2563       if(id.eq.331)                 ihi=10      !phi
2564       !if(id.eq.)                        ihi=11          !h_1
2565       if(ida.eq.130.or.ida.eq.230)  ihi=12      !K
2566       if(ida.eq.131.or.ida.eq.231)  ihi=13      !K star
2567       if(ida.eq.1120.or.ida.eq.1220)ihi=14      !N
2568       if(ida.eq.1111.or.ida.eq.1121)ihi=15      !Delta
2569       if(ida.eq.1221.or.ida.eq.2221)ihi=15      !Delta
2570       if(ida.eq.2130)               ihi=16      !Lambda
2571       if(ida.eq.1130.or.ida.eq.1230)ihi=17      !Sigma
2572       if(ida.eq.2230)               ihi=17      !Sigma
2573       idxHiranoTable=ihi
2574       end
2575 
2576 c----------------------------------------------------------------------
2577       subroutine InitializeHyperbola
2578 c----------------------------------------------------------------------
2579       include 'epos.inc'
2580       double precision tpro,zpro,ttar,ztar,ttaus,detap,detat
2581       common /cttaus/  tpro,zpro,ttar,ztar,ttaus,detap,detat
2582       ttaus=1
2583       ypjtl=6
2584       yhaha=3
2585       etapro=(ypjtl-yhaha)*etafac
2586       etatar=-yhaha*etafac
2587       detap=dble(etapro)
2588       detat=dble(etatar)
2589       tpro=dcosh(detap)
2590       zpro=dsinh(detap)
2591       ttar=dcosh(detat)
2592       ztar=dsinh(detat)
2593       end
2594 
2595 c----------------------------------------------------------------------
2596       integer function idGet()
2597 c----------------------------------------------------------------------
2598       include 'epos.inc'
2599       integer jc(nflav,2),ic(2)
2600       if(ish.ge.5)write(ifch,'(a/6i4)')
2601      *' keu ked kes kec keb ket:',keu,ked,kes,kec,keb,ket
2602       do i=1,6
2603       jc(i,1)=0
2604       jc(i,2)=0
2605       enddo
2606       if(keu.ge.0)jc(1,1)=keu
2607       if(ked.ge.0)jc(2,1)=ked
2608       if(kes.ge.0)jc(3,1)=kes
2609       if(kec.ge.0)jc(4,1)=kec
2610       if(keb.ge.0)jc(5,1)=keb
2611       if(ket.ge.0)jc(6,1)=ket
2612       if(keu.lt.0)jc(1,2)=-keu
2613       if(ked.lt.0)jc(2,2)=-ked
2614       if(kes.lt.0)jc(3,2)=-kes
2615       if(kec.lt.0)jc(4,2)=-kec
2616       if(keb.lt.0)jc(5,2)=-keb
2617       if(ket.lt.0)jc(6,2)=-ket
2618       if(ish.ge.5)write(ifch,'(a/6i4/6i4)')' jc:',jc
2619       idr=0
2620       do  nf=1,nflav
2621         do  ij=1,2
2622           if(jc(nf,ij).ge.10)idr=7*10**8
2623         enddo
2624       enddo
2625       if(idr/10**8.ne.7)then
2626         call idenco(jc,ic,ireten)
2627         if(ic(1).eq.0.and.ic(2).eq.0)then
2628           ic(1)=100000
2629           ic(2)=100000
2630         endif
2631 
2632         idr=8*10**8+ic(1)*100+ic(2)/100
2633         if(ish.ge.5)write(ifch,'(a,i9)')' id:',idr
2634       else
2635         idr=idr
2636      *       +mod(jc(1,1)+jc(2,1)+jc(3,1)+jc(4,1),10**4)*10**4
2637      *       +mod(jc(1,2)+jc(2,2)+jc(3,2)+jc(4,2),10**4)
2638         call idtrbi(jc,ibptl(1,1),ibptl(2,1),ibptl(3,1),ibptl(4,1))
2639       endif
2640       idGet=idr
2641       end
2642