]> git.uio.no Git - u/mrichter/AliRoot.git/blob - EPOS/epos167/epos-ico-158.f
following coverity, fix for case of invalid DET string passed to GetNalignable
[u/mrichter/AliRoot.git] / EPOS / epos167 / epos-ico-158.f
1 c-----------------------------------------------------------------------
2 c  activate making inicon table by the follwing commands in optns file:
3 c-----------------------------------------------------------------------
4 c        set nevent 1
5 c        set bminim 3      (for example)
6 c        set bmaxim 5      (for example)
7 c        set ninicon 1000  (for example)
8 c        set icocore 1  !or 2
9 c        make icotable
10 c        switch  splitting on
11 c        switch  fusion on
12 c        switch clusterdecay off
13 c-----------------------------------------------------------------------
14 c  inicon table means storage of IcoE, IcoV, IcoF
15 c     IcoE ... energy density in the comoving frame
16 c     IcoV ... flow velocity in hyperbolic coordinates 
17 c               v=(vx,vy,veta=vz/tau) 
18 c                with vz being the velocity in the frame moving with rapidity eta
19 c     IcoV ... net flavor density in the comoving frame
20 c  the indices of these fields ix,iy,iz also refer to hyperbolic coordinates
21 c     (x,y,eta) at given tau
22 c     the corresponding grid is defined in epos.incico
23 c-----------------------------------------------------------------------
24
25 c-----------------------------------------------------------------------
26       subroutine IniCon(nin)
27 c-----------------------------------------------------------------------
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
118 c-----------------------------------------------------------------------
119       subroutine IcoStr(iflag)
120 c-----------------------------------------------------------------------
121 c     energy momentum tensor and flavor current 
122 c     and corresponding contractions from string method
123 c
124 c     iflag = 1 initialization   
125 c           = 2 sum up density
126 c           = 3 normalization     
127 c     
128 c     output: common blocks /Ico3/  /Ico4/  /Ico5/ 
129 c-----------------------------------------------------------------------
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
143 c......................................................................
144 c                             initialization      
145 c......................................................................
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
173 c......................................................................
174 c                        fill arrays      
175 c......................................................................
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
243 c............................................................................
244 c                 normalization and diagonalization
245 c............................................................................
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       
307 c------------------------------------------------------------------------------------      
308       subroutine DiagTmunu(Tmunu,u,v,gamma,eps,nonz,ix,iy,iz)
309 c------------------------------------------------------------------------------------      
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