Update master to aliroot
[u/mrichter/AliRoot.git] / EPOS / epos167 / epos-ico-158.f
CommitLineData
9ef1c2d9 1c-----------------------------------------------------------------------
2c activate making inicon table by the follwing commands in optns file:
3c-----------------------------------------------------------------------
4c set nevent 1
5c set bminim 3 (for example)
6c set bmaxim 5 (for example)
7c set ninicon 1000 (for example)
8c set icocore 1 !or 2
9c make icotable
10c switch splitting on
11c switch fusion on
12c switch clusterdecay off
13c-----------------------------------------------------------------------
14c inicon table means storage of IcoE, IcoV, IcoF
15c IcoE ... energy density in the comoving frame
16c IcoV ... flow velocity in hyperbolic coordinates
17c v=(vx,vy,veta=vz/tau)
18c with vz being the velocity in the frame moving with rapidity eta
19c IcoV ... net flavor density in the comoving frame
20c the indices of these fields ix,iy,iz also refer to hyperbolic coordinates
21c (x,y,eta) at given tau
22c the corresponding grid is defined in epos.incico
23c-----------------------------------------------------------------------
24
25c-----------------------------------------------------------------------
26 subroutine IniCon(nin)
27c-----------------------------------------------------------------------
28 include 'epos.inc'
29 include 'epos.incico'
30 double precision seedf
31 character*80 fn
32 logical table
33
34 ii=index(fnhi(1:nfnhi),".histo")-1
35 fn=fnhi(1:ii)//".ico"
36 inquire(file=fn(1:ii+4),exist=table)
37
38 if(icotabr.eq.1)then !read from file
39
40 if(.not.table)then
41 write(ifmt,'(//10x,3a//)')'STOP: file '
42 * ,fn(1:ii+4),' not found !!!'
43 stop
44 endif
45 write(ifmt,'(3a)')'read from ',fn(1:ii+4),' ...'
46 open(97,file=fn,status='old')
47 read(97,*) laprojx,maprojx,latargx,matargx
48 read(97,*) engyx
49 read(97,*) bminimx,bmaximx
50 read(97,*) tauicox
51 read(97,*) iabs_ninicon
52 read(97,*) nxicox,nyicox,nzicox
53 read(97,*)xminicox,xmaxicox,yminicox,ymaxicox,zminicox,zmaxicox
54 if(laprojx.ne.laproj.or.maprojx.ne.maproj
55 * .or.latargx.ne.latarg.or.matargx.ne.matarg)stop'(1)'
56 if(engyx.ne.engy) stop'variables dont match (2) '
57 if(bminimx.ne.bminim.or.bmaximx.ne.bmaxim)
58 * stop'variables dont match (3) '
59 if(tauicox.ne.tauico)stop'variables dont match (4) '
60 if(nxicox.ne.nxico.or.nyicox.ne.nyico
61 * .or.nzicox.ne.nzico)stop'(5)'
62 if( xminicox.ne.xminico.or.xmaxicox.ne.xmaxico
63 * .or.yminicox.ne.yminico.or.ymaxicox.ne.ymaxico
64 * .or.zminicox.ne.zminico.or.zmaxicox.ne.zmaxico)stop'(6)'
65 read(97,*) IcoE,IcoV,IcoF
66 close(97)
67
68 elseif(icotabr.eq.0)then ! calculate
69
70 !if(icotabm.eq.1.and.table)then
71 !write(ifmt,'(//10x,3a//)')'STOP: file '
72 !* ,fn(1:ii+4),' already exists !!!'
73 !stop
74 !endif
75
76 if(nin.eq.1)then
77 write(ifmt,'(2a,i7,a)')' calculate inicon table ',
78 & 'based on',iabs(ninicon),' ico-events ...'
79 call IcoStr(1)
80 endif
81
82 call ranfgt(seedf)
83 if(nin.le.iabs(ninicon).and.mod(nin,modsho).eq.0)
84 & write(ifmt,'(a,i7,a,f6.2,a,d27.16)')' ico-event ',nin
85 & ,' bimevt:',bimevt,' seedf:',seedf
86 !if(nin.le.iabs(ninicon).and.mod(nin,modsho).eq.0)
87 !& print*,'+++++ time: ',timefin-timeini
88 seedc=seedf
89 call IcoStr(2)
90
91 if(nin.eq.iabs(ninicon))then
92 write(ifmt,*)'normalize and diagonalize engy-mom tensor'
93 call IcoStr(3)
94
95 if(icotabm.eq.1)then ! make table
96 write(ifmt,'(2a)')' write inicon table into ',fn(1:ii+4)
97 open(97,file=fn,status='unknown')
98 write(97,*) laproj,maproj,latarg,matarg
99 write(97,*) engy
100 write(97,*) bminim,bmaxim
101 write(97,*) tauico
102 write(97,*) iabs(ninicon)
103 write(97,*) nxico,nyico,nzico
104 write(97,*) xminico,xmaxico,yminico,ymaxico,zminico,zmaxico
105 write(97,*) IcoE,IcoV,IcoF
106 close(97)
107 endif
108 endif
109
110 else
111
112 stop'IniCon: wrong choice'
113
114 endif
115
116 end
117
118c-----------------------------------------------------------------------
119 subroutine IcoStr(iflag)
120c-----------------------------------------------------------------------
121c energy momentum tensor and flavor current
122c and corresponding contractions from string method
123c
124c iflag = 1 initialization
125c = 2 sum up density
126c = 3 normalization
127c
128c output: common blocks /Ico3/ /Ico4/ /Ico5/
129c-----------------------------------------------------------------------
130 include 'epos.inc'
131 include 'epos.incico'
132 double precision IcoT(4,4,nxico,nyico,nzico)
133 * ,IcoC(3,4,nxico,nyico,nzico)
134 common/cxyzt/xptl(mxptl),yptl(mxptl),zptl(mxptl),tptl(mxptl)
135 *,optl(mxptl),uptl(mxptl),sptl(mxptl),rptl(mxptl,3)
136 common /Ico10/ntot,ish0,irs
137 double precision p(4),v(4),u(4),eps,T(4,4),gamma,rap,rapx
138 integer ic(2),jc(nflav,2)
139 logical ok
140 goto (1,10,100),iflag
141 return
142
143c......................................................................
144c initialization
145c......................................................................
146
147 1 do ix=1,nxico
148 do iy=1,nyico
149 do iz=1,nzico
150 IcoE(ix,iy,iz)=0.
151 do i=1,3
152 IcoV(i,ix,iy,iz)=0.
153 enddo
154 do i=1,3
155 IcoF(i,ix,iy,iz)=0.
156 enddo
157 do i=1,4
158 do j=1,4
159 IcoT(i,j,ix,iy,iz)=0.
160 enddo
161 enddo
162 do i=1,3
163 do j=1,4
164 IcoC(i,j,ix,iy,iz)=0.
165 enddo
166 enddo
167 enddo
168 enddo
169 enddo
170 ntot=0
171 return
172
173c......................................................................
174c fill arrays
175c......................................................................
176
177 10 ntot=ntot+1
178 dxico=(xmaxico-xminico)/nxico
179 dyico=(ymaxico-yminico)/nyico
180 dzico=(zmaxico-zminico )/nzico
181
182 do j=1,nptl
183 ok=.true.
184 if(icocore.eq.2
185 . .and.(ityptl(j).lt.20.or.ityptl(j).ge.40))ok=.false.
186 if(ok
187 . .and.dezptl(j).lt.1e3.and.j.le.maxfra.and.istptl(j).eq.2)then
188 aa=cos(phievt)
189 bb=sin(phievt)
190 x=xptl(j)*aa+yptl(j)*bb
191 cc=-sin(phievt)
192 dd=cos(phievt)
193 y=xptl(j)*cc+yptl(j)*dd
194 z=dezptl(j)
195 ix=1+(x-xminico)/dxico
196 iy=1+(y-yminico)/dyico
197 iz=1+(z-zminico)/dzico
198 if(ix.ge.1.and.ix.le.nxico.and.
199 . iy.ge.1.and.iy.le.nyico.and.
200 . iz.ge.1.and.iz.le.nzico)then
201 !~~~~~~~~~~~~
202 ! T^i^k = \int d3p p^i p^k /E f(\vec x,\vec p)
203 ! N^k = \int d3p p^k /E f_net(\vec x,\vec p)
204 !~~~~~~~~~~~~
205 rapx=zminico+(iz-0.5)*dzico
206 amt2=pptl(5,j)**2+pptl(1,j)**2+pptl(2,j)**2
207 pp=pptl(4,j)+abs(pptl(3,j))
208 if(amt2.gt.0..and.pp.gt.0.)then
209 amt=sqrt(amt2)
210 rap=sign(1.,pptl(3,j))*alog(pp/amt)
211 p(1)=pptl(1,j)
212 p(2)=pptl(2,j)
213 p(3)=amt*sinh(rap-rapx)
214 p(4)=amt*cosh(rap-rapx)
215 do i=1,4
216 do k=1,4
217 IcoT(i,k,ix,iy,iz)=IcoT(i,k,ix,iy,iz)+p(i)*p(k)/p(4)
218 enddo
219 enddo
220 id=idptl(j)
221 ida=iabs(id/10)
222 ids=id/iabs(id)
223 if(ida.ne.111.and.ida.ne.222.and.ida.ne.333)id=id/10*10
224 if(ida.eq.111.or. ida.eq.222.or. ida.eq.333)id=id/10*10+ids
225 if(ida.eq.213)id=1230*ids
226 ic(1)=idtrai(1,id,1)
227 ic(2)=idtrai(2,id,1)
228 call iddeco(ic,jc)
229 do i=1,3
230 fi=jc(i,1)-jc(i,2)
231 do k=1,4
232 IcoC(i,k,ix,iy,iz)=IcoC(i,k,ix,iy,iz)+p(k)/p(4)*fi
233 enddo
234 enddo
235 endif
236 !~~~~~~~~~~~~
237 endif
238 endif
239 enddo
240
241 return
242
243c............................................................................
244c normalization and diagonalization
245c............................................................................
246
247 100 continue
248
249 !~~~normalize
250
251 vol= (xmaxico-xminico)/nxico
252 . *(ymaxico-yminico)/nyico
253 . *(zmaxico-zminico)/nzico *tauico
254 fac=1./ntot/vol
255
256 do ix=1,nxico
257 do iy=1,nyico
258 do iz=1,nzico
259 do i=1,4
260 do k=1,4
261 IcoT(i,k,ix,iy,iz)=IcoT(i,k,ix,iy,iz)*fac
262 enddo
263 enddo
264 do i=1,3
265 do k=1,4
266 IcoC(i,k,ix,iy,iz)=IcoC(i,k,ix,iy,iz)*fac
267 enddo
268 enddo
269 enddo
270 enddo
271 enddo
272
273 !~~~diagonalize T
274
275 do ix=1,nxico
276 do iy=1,nyico
277 do iz=1,nzico
278 do i=1,4
279 do k=1,4
280 T(i,k)=IcoT(i,k,ix,iy,iz)
281 enddo
282 enddo
283 call DiagTmunu(T,u,v,gamma,eps,nonz,ix,iy,iz)
284 if(nonz.eq.1)then
285 IcoE(ix,iy,iz)=eps
286 IcoV(1,ix,iy,iz)=v(1)
287 IcoV(2,ix,iy,iz)=v(2)
288 !rapx=zminico+(iz-0.5)*dzico
289 !rap=rapx+log(gamma*(1+v(3)))
290 !IcoV(3,ix,iy,iz)=tanh(rap)
291 IcoV(3,ix,iy,iz)=v(3) / tauico
292 do i=1,3
293 IcoF(i,ix,iy,iz)=IcoC(i,4,ix,iy,iz)*u(4)
294 do k=1,3
295 IcoF(i,ix,iy,iz)=IcoF(i,ix,iy,iz)
296 . -IcoC(i,k,ix,iy,iz)*u(k)
297 enddo
298 enddo
299 endif
300 enddo
301 enddo
302 enddo
303
304 return
305 end
306
307c------------------------------------------------------------------------------------
308 subroutine DiagTmunu(Tmunu,u,v,gamma,eps,nonz,ix,iy,iz)
309c------------------------------------------------------------------------------------
310 ! solve lambda * T * lambda = diag(eps,P,P,P), for symmetric T
311 !
312 ! lambda = g g*vx g*vy g*vz =symmetric!
313 ! g*vx a*vx*vx+1 a*vy*vx a*vz*vx
314 ! g*vy a*vx*vy a*vy*vy+1 a*vz*vy
315 ! g*vz a*vx*vz a*vy*vz a*vz*vz+1
316 ! with g=gamma, a=g*g/(g+1)
317 !
318 ! so T*lambda(v)=lambda(-v)*diag(eps,P,P,P)
319 ! first column: four equations
320 ! eps=T00+T0k*vk
321 ! -eps*vi=Ti0+Tik*vk
322 ! solved iteratively to get eps, vi
323 ! returns u,v,eps if nonz=1 (otherwise zero tensor)
324 ! (by T.Kodama)
325 ! (K. Werner: modifs)
326
327 include 'epos.incico'
328 double precision Tmunu(4,4), u(4),eps,gamma,g,a, Lor(4,4),w(4),sum
329 double precision v(4),vx(3),tt(4,4),err,sg(4)
330
331 sg(4)=1d0
332 v(4)=1d0
333 do i=1,3
334 sg(i)=-1d0
335 enddo
336 sum=0d0
337 do i=1,4
338 do k=1,4
339 sum=sum+dabs(Tmunu(i,k))
340 end do
341 end do
342 nonz=0
343 if(sum.eq.0.0d0)return
344 nonz=1
345
346 do k=1,3
347 v(k)=0.
348 end do
349 eps=0
350
351 do lrep=1,100
352 epsx=eps
353 do i=1,3
354 vx(i)=v(i)
355 end do
356 eps=Tmunu(4,4)
357 do k=1,3
358 eps=eps-Tmunu(4,k)*vx(k)
359 end do
360 if(eps.le.0d0)then
361 print*,'DiagTmunu: sum(abs(Tmunu))=',sum,' eps=',eps
362 stop'negative energy density. '
363 endif
364 do i=1,3
365 Tv=0d0
366 do k=1,3
367 Tv=Tv+Tmunu(i,k)*vx(k)
368 end do
369 v(i)=(Tmunu(i,4)-Tv)/eps
370 end do
371 if(lrep.gt.60)then
372 do i=1,3
373 v(i)=0.5d0*(vx(i)+v(i))
374 enddo
375 endif
376 !print*,'Tmunu: ',lrep,abs(eps-epsx),(abs(v(i)-vx(i)),i=1,3)
377 err=1d-6
378 if(lrep.gt.50)err=1d-5
379 if(lrep.gt.89)err=1d-4
380 if(abs(eps-epsx).lt.err.and.abs(v(1)-vx(1)).lt.err
381 . .and.abs(v(2)-vx(2)).lt.err.and.abs(v(3)-vx(3)).lt.err)goto1
382 do i=1,4
383 w(i)=0
384 do k=1,4
385 w(i)=w(i)+Tmunu(i,k)*v(k)*sg(k)
386 enddo
387 w(i)=w(i)-eps*v(i)
388 enddo
389 if(lrep.gt.95
390 ..and.w(1)*w(1)+w(2)*w(2)+w(3)*w(3)+w(4)*w(4).lt.err)goto1
391 end do
392
393 1 v2=0.d0
394 do i=1,3
395 v2=v2+v(i)**2
396 end do
397 if(v2.ge.1.)stop'DiagTmunu: v2 ge 1. '
398 gamma=1./sqrt(abs(1.-v2))
399 u(4)=gamma
400 u(1)=v(1)*gamma
401 u(2)=v(2)*gamma
402 u(3)=v(3)*gamma
403
404 !~~~check
405 g=gamma
406 a=g*g/(g+1d0)
407 Lor(4,4)=g
408 do k=1,3
409 Lor(4,k)=-g*v(k)
410 Lor(k,4)=Lor(4,k)
411 enddo
412 do i=1,3
413 do k=1,3
414 Lor(i,k)=a*v(i)*v(k)
415 enddo
416 enddo
417 do k=1,3
418 Lor(k,k)=Lor(k,k)+1
419 enddo
420 do i=1,4
421 do k=1,4
422 tt(i,k)=0d0
423 do m=1,4
424 do n=1,4
425 tt(i,k)=tt(i,k)+Lor(i,m)*Tmunu(m,n)*Lor(n,k)
426 enddo
427 enddo
428 enddo
429 enddo
430 err=err*1d2
431 if(tt(1,4).gt.err.or.tt(2,4).gt.err.or.tt(3,4).gt.err)then
432 print*,'************** DiagTmunu: nonzero T14 or T24 or T34'
433 . ,' after ',lrep,' iterations'
434 print*,'d_eps or d_v(i): ',abs(eps-epsx),(abs(v(i)-vx(i)),i=1,3)
435 print*,'Tmunu ( ix=',ix-nxico/2-1,' iy=',iy-nyico/2-1,' ) :'
436 . ,' iz=',iz-nzico/2-1
437 do i=1,4
438 write(*,'(4f10.5,2x,4f10.5)')
439 . (sngl(tmunu(i,k)),k=1,4),(sngl(tt(i,k)),k=1,4)
440 enddo
441 endif
442
443 end
444