File indexing completed on 2024-04-06 12:14:01
0001
0002 subroutine ahydro
0003
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
0040 subroutine amicro(iret)
0041
0042
0043
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
0068 call InitializeHyperbola !does not affect results but necessary
0069
0070
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
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
0111 subroutine GraCan
0112
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
0162 subroutine HydroTable2(ichk)
0163
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
0181 subroutine HydroTable(ichk,tabname) ! 13 Aug 08 (same as gyx.f except inc)
0182
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
0243 subroutine SymmTab(ncent1,ncent2) ! 13 Aug 08 (same as gyx.f except inc)
0244
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
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
0402 subroutine DefineParticles ! 14 Aug 08 (same as gyx.f except inc)
0403
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
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
