ReggeGribovPartonMCInterface/src/epos-ico.f

0001 c-----------------------------------------------------------------------
0002 c  activate making inicon table by the follwing commands in optns file:
0003 c-----------------------------------------------------------------------
0004 c        set nevent 1
0005 c        set bminim 3      (for example)
0006 c        set bmaxim 5      (for example)
0007 c        set ninicon 1000  (for example)
0008 c        set icocore 1  !or 2
0009 c        make icotable
0010 c        switch  splitting on
0011 c        switch  fusion on
0012 c        switch clusterdecay off
0013 c-----------------------------------------------------------------------
0014 c  inicon table means storage of IcoE, IcoV, IcoF
0015 c     IcoE ... energy density in the comoving frame
0016 c     IcoV ... flow velocity in hyperbolic coordinates
0017 c               v=(vx,vy,veta=vz/tau)
0018 c                with the velocity in the frame moving with rapidity eta
0019 c     IcoV ... net flavor density in the comoving frame
0020 c  the indices of these fields ix,iy,iz also refer to hyperbolic coordinates
0021 c     (x,y,eta) at given tau
0022 c     the corresponding grid is defined in epos.incico
0023 c-----------------------------------------------------------------------
0024
0025 c-----------------------------------------------------------------------
0026       subroutine IniCon(nin)
0027 c-----------------------------------------------------------------------
0028       include 'epos.inc'
0029       include 'epos.incico'
0030       double precision seedf
0031       character*80 fn
0032       logical table
0033
0034       ii=index(fnhi(1:nfnhi),".histo")-1
0035       fn=fnhi(1:ii)//".ico"
0036       inquire(file=fn(1:ii+4),exist=table)
0037
0038       if(icotabr.eq.1)then  !read from file
0039
0040         if(.not.table)then
0041           write(ifmt,'(//10x,3a//)')'STOP: file '
0042      *             ,fn(1:ii+4),' not found !!!'
0043           stop
0044         endif
0045         write(ifmt,'(3a)')'read from ',fn(1:ii+4),' ...'
0046         open(97,file=fn,status='old')
0047         read(97,*) iversn
0048         read(97,*) laprojx,maprojx,latargx,matargx
0049         read(97,*) engyx
0050         read(97,*) bminimx,bmaximx
0051         read(97,*) tauicox
0052         read(97,*) iabs_ninicon_dum
0053         read(97,*) nxicox,nyicox,nzicox
0054         read(97,*)xminicox,xmaxicox,yminicox,ymaxicox,zminicox,zmaxicox
0055         if(laprojx.ne.laproj.or.maprojx.ne.maproj
0056      * .or.latargx.ne.latarg.or.matargx.ne.matarg)stop'(1)'
0057         if(engyx.ne.engy)    stop'variables dont match (2)    '
0058         if(bminimx.ne.bminim.or.bmaximx.ne.bmaxim)
0059      *                       stop'variables dont match (3)    '
0060         if(tauicox.ne.tauico)stop'variables dont match (4)    '
0061         if(nxicox.ne.nxico.or.nyicox.ne.nyico
0062      * .or.nzicox.ne.nzico)stop'(5)'
0063         if( xminicox.ne.xminico.or.xmaxicox.ne.xmaxico
0064      * .or.yminicox.ne.yminico.or.ymaxicox.ne.ymaxico
0065      * .or.zminicox.ne.zminico.or.zmaxicox.ne.zmaxico)stop'(6)'
0066         read(97,*) IcoE,IcoV,IcoF
0067         close(97)
0068
0069       elseif(icotabr.eq.0)then  ! calculate
0070
0071         !if(icotabm.eq.1.and.table)then
0072           !write(ifmt,'(//10x,3a//)')'STOP: file '
0073           !*             ,fn(1:ii+4),' already exists !!!'
0074           !stop
0075         !endif
0076
0077         if(nin.eq.1)then
0078           write(ifmt,'(2a,i7,a)')' calculate inicon table ',
0079      &    'based on',iabs(ninicon),' ico-events ...'
0080           call IcoStr(1)
0081         endif
0082
0083         call ranfgt(seedf)
0084         if(nin.le.iabs(ninicon).and.mod(nin,modsho).eq.0)
0085      &        write(ifmt,'(a,i7,a,f6.2,a,d27.16)')' ico-event ',nin
0086      &                   ,'   bimevt:',bimevt,'   seedf:',seedf
0087         if(nin.le.iabs(ninicon).and.jwseed.eq.1)then
0088           open(unit=1,file=fnch(1:nfnch-5)//'see',status='unknown')
0089           write(1,'(a,i7,a,f6.2,a,d27.16)')' ico-event ',nin
0090      &                   ,'   bimevt:',bimevt,'   seedf:',seedf
0091          close(1)
0092         endif
0093         !if(nin.le.iabs(ninicon).and.mod(nin,modsho).eq.0)
0094         !&         print*,'+++++ time: ',timefin-timeini
0095         seedc=seedf
0096         call IcoStr(2)
0097
0098         if(nin.eq.iabs(ninicon))then
0099           write(ifmt,*)'normalize and diagonalize engy-mom tensor'
0100            call IcoStr(3)
0101
0102            if(icotabm.eq.1)then  ! make table
0103             write(ifmt,'(2a)')' write inicon table into ',fn(1:ii+4)
0104             open(97,file=fn,status='unknown')
0105             write(97,*) iversn
0106             write(97,*) laproj,maproj,latarg,matarg
0107             write(97,*) engy
0108             write(97,*) bminim,bmaxim
0109             write(97,*) tauico
0110             write(97,*) iabs(ninicon)
0111             write(97,*) nxico,nyico,nzico
0112             write(97,*) xminico,xmaxico,yminico,ymaxico,zminico,zmaxico
0113             write(97,*) IcoE,IcoV,IcoF
0114             close(97)
0115           endif
0116          endif
0117
0118       else
0119
0120         stop'IniCon: wrong choice'
0121
0122       endif
0123
0124       end
0125
0126 c-----------------------------------------------------------------------
0127       subroutine IcoStr(iflag)
0128 c-----------------------------------------------------------------------
0129 c     energy momentum tensor and flavor current
0130 c     and corresponding contractions from string method
0131 c
0132 c     iflag = 1 initialization
0133 c           = 2 sum up density
0134 c           = 3 normalization
0135 c
0136 c     output: common blocks /Ico3/  /Ico4/  /Ico5/
0137 c-----------------------------------------------------------------------
0138       include 'epos.inc'
0139       include 'epos.incico'
0140       double precision IcoT(4,4,nxico,nyico,nzico)
0141      *                ,IcoC(3,4,nxico,nyico,nzico)
0142       common/cxyzt/xptl(mxptl),yptl(mxptl),zptl(mxptl),tptl(mxptl)
0143      *,optl(mxptl),uptl(mxptl),sptl(mxptl),rptl(mxptl,3)
0144       common /Ico10/ntot,ish0,irs
0145       double precision p(4),v(4),u(4),eps,T(4,4),gamma,rap,rapx
0146       integer ic(2),jc(nflav,2)
0147       common/cranphi/ranphi
0148       logical ok
0149       goto (1,10,100),iflag
0150       return
0151
0152 c......................................................................
0153 c                             initialization
0154 c......................................................................
0155
0156  1    do ix=1,nxico
0157         do iy=1,nyico
0158           do iz=1,nzico
0159             IcoE(ix,iy,iz)=0.
0160             do i=1,3
0161               IcoV(i,ix,iy,iz)=0.
0162             enddo
0163             do i=1,3
0164               IcoF(i,ix,iy,iz)=0.
0165             enddo
0166             do i=1,4
0167             do j=1,4
0168               IcoT(i,j,ix,iy,iz)=0.
0169             enddo
0170             enddo
0171             do i=1,3
0172             do j=1,4
0173               IcoC(i,j,ix,iy,iz)=0.
0174             enddo
0175             enddo
0176           enddo
0177         enddo
0178       enddo
0179       ntot=0
0180       return
0181
0182 c......................................................................
0183 c                        fill arrays
0184 c......................................................................
0185
0186  10   ntot=ntot+1
0187       dxico=(xmaxico-xminico)/nxico
0188       dyico=(ymaxico-yminico)/nyico
0189       dzico=(zmaxico-zminico )/nzico
0190       phinll=phievt+ranphi
0191
0192       do j=1,nptl
0193        ok=.true.
0194        if(icocore.eq.2
0195      . .and.(ityptl(j).lt.20.or.ityptl(j).ge.40))ok=.false.
0196        if(ok
0197      . .and.dezptl(j).lt.1e3.and.j.le.maxfra.and.istptl(j).eq.2)then
0198         aa=cos(phinll)
0199         bb=sin(phinll)
0200         x=xptl(j)*aa+yptl(j)*bb
0201         cc=-sin(phinll)
0202         dd=cos(phinll)
0203         y=xptl(j)*cc+yptl(j)*dd
0204         z=dezptl(j)
0205         ix=1+(x-xminico)/dxico
0206         iy=1+(y-yminico)/dyico
0207         iz=1+(z-zminico)/dzico
0208         if(ix.ge.1.and.ix.le.nxico.and.
0209      .  iy.ge.1.and.iy.le.nyico.and.
0210      .  iz.ge.1.and.iz.le.nzico)then
0211           !~~~~~~~~~~~~
0212           ! T^i^k = \int d3p p^i p^k /E f(\vec x,\vec p)
0213           ! N^k =   \int d3p p^k /E f_net(\vec x,\vec p)
0214           !~~~~~~~~~~~~
0215           rapx=zminico+(iz-0.5)*dzico
0216           amt2=pptl(5,j)**2+pptl(1,j)**2+pptl(2,j)**2
0217           pp=pptl(4,j)+abs(pptl(3,j))
0218           if(amt2.gt.0..and.pp.gt.0.)then
0219            amt=sqrt(amt2)
0220            rap=sign(1.,pptl(3,j))*alog(pp/amt)
0221            p(1)=pptl(1,j)
0222            p(2)=pptl(2,j)
0223            p(3)=amt*sinh(rap-rapx)
0224            p(4)=amt*cosh(rap-rapx)
0225            do i=1,4
0226             do k=1,4
0227               IcoT(i,k,ix,iy,iz)=IcoT(i,k,ix,iy,iz)+p(i)*p(k)/p(4)
0228             enddo
0229            enddo
0230            id=idptl(j)
0231            ida=iabs(id/10)
0232            ids=id/iabs(id)
0233            if(ida.ne.111.and.ida.ne.222.and.ida.ne.333)id=id/10*10
0234            if(ida.eq.111.or. ida.eq.222.or. ida.eq.333)id=id/10*10+ids
0235            if(ida.eq.213)id=1230*ids
0236            ic(1)=idtrai(1,id,1)
0237            ic(2)=idtrai(2,id,1)
0238            call iddeco(ic,jc)
0239            do i=1,3
0240             fi=jc(i,1)-jc(i,2)
0241             do k=1,4
0242              IcoC(i,k,ix,iy,iz)=IcoC(i,k,ix,iy,iz)+p(k)/p(4)*fi
0243             enddo
0244            enddo
0245           endif
0246           !~~~~~~~~~~~~
0247         endif
0248        endif
0249       enddo
0250
0251       return
0252
0253 c............................................................................
0254 c                 normalization and diagonalization
0255 c............................................................................
0256
0257  100  continue
0258
0259       !~~~normalize
0260
0261       vol= (xmaxico-xminico)/nxico
0262      .    *(ymaxico-yminico)/nyico
0263      .    *(zmaxico-zminico)/nzico  *tauico
0264       fac=1./ntot/vol
0265
0266       do ix=1,nxico
0267         do iy=1,nyico
0268           do iz=1,nzico
0269             do i=1,4
0270              do k=1,4
0271               IcoT(i,k,ix,iy,iz)=IcoT(i,k,ix,iy,iz)*fac
0272              enddo
0273             enddo
0274             do i=1,3
0275              do k=1,4
0276               IcoC(i,k,ix,iy,iz)=IcoC(i,k,ix,iy,iz)*fac
0277              enddo
0278             enddo
0279           enddo
0280         enddo
0281       enddo
0282
0283       !~~~diagonalize T
0284
0285       do ix=1,nxico
0286         do iy=1,nyico
0287           do iz=1,nzico
0288             do i=1,4
0289             do k=1,4
0290              T(i,k)=IcoT(i,k,ix,iy,iz)
0291             enddo
0292             enddo
0293             call DiagTmunu(T,u,v,gamma,eps,nonz,ix,iy,iz)
0294             if(nonz.eq.1)then
0295               IcoE(ix,iy,iz)=eps
0296               IcoV(1,ix,iy,iz)=v(1)
0297               IcoV(2,ix,iy,iz)=v(2)
0298                !rapx=zminico+(iz-0.5)*dzico
0299                !rap=rapx+log(gamma*(1+v(3)))
0300                !IcoV(3,ix,iy,iz)=tanh(rap)
0301               IcoV(3,ix,iy,iz)=v(3) / tauico
0302               do i=1,3
0303                IcoF(i,ix,iy,iz)=IcoC(i,4,ix,iy,iz)*u(4)
0304                do  k=1,3
0305                 IcoF(i,ix,iy,iz)=IcoF(i,ix,iy,iz)
0306      .                  -IcoC(i,k,ix,iy,iz)*u(k)
0307                enddo
0308               enddo
0309             endif
0310           enddo
0311         enddo
0312       enddo
0313
0314       return
0315       end
0316
0317 c------------------------------------------------------------------------------------
0318       subroutine DiagTmunu(Tmunu,u,v,gamma,eps,nonz,ix,iy,iz)
0319 c------------------------------------------------------------------------------------
0320       ! solve lambda * T * lambda = diag(eps,P,P,P), for symmetric T
0321       !
0322       !  lambda =  g      g*vx        g*vy        g*vz          =symmetric!
0323       !            g*vx   a*vx*vx+1   a*vy*vx     a*vz*vx
0324       !            g*vy   a*vx*vy     a*vy*vy+1   a*vz*vy
0325       !            g*vz   a*vx*vz     a*vy*vz     a*vz*vz+1
0326       ! with g=gamma, a=g*g/(g+1)
0327       !
0328       ! so T*lambda(v)=lambda(-v)*diag(eps,P,P,P)
0329       ! first column: four equations
0330       !                    eps=T00+T0k*vk
0331       !                -eps*vi=Ti0+Tik*vk
0332       ! solved iteratively to get eps, vi
0333       ! returns u,v,eps if nonz=1 (otherwise zero tensor)
0334       !  (by T.Kodama)
0335       !  (K. Werner: modifs)
0336
0337       include 'epos.incico'
0338       double precision Tmunu(4,4), u(4),eps,gamma,g,a, Lor(4,4),w(4),sum
0339       double precision v(4),vx(3),tt(4,4),err,sg(4)
0340
0341       sg(4)=1d0
0342       v(4)=1d0
0343       do i=1,3
0344        sg(i)=-1d0
0345       enddo
0346       sum=0d0
0347       do i=1,4
0348        do k=1,4
0349        sum=sum+dabs(Tmunu(i,k))
0350        end do
0351       end do
0352       nonz=0
0353       if(sum.eq.0.0d0)return
0354       nonz=1
0355
0356       do k=1,3
0357        v(k)=0.
0358       end do
0359       eps=0
0360
0361       do lrep=1,100
0362        epsx=eps
0363        do i=1,3
0364         vx(i)=v(i)
0365        end do
0366        eps=Tmunu(4,4)
0367        do k=1,3
0368         eps=eps-Tmunu(4,k)*vx(k)
0369        end do
0370        if(eps.le.0d0)then
0371          print*,'DiagTmunu: sum(abs(Tmunu))=',sum,'   eps=',eps
0372          print*,'Tmunu(4,nu):',(Tmunu(4,nu),nu=1,4)
0373          if(eps.gt.-1e-5)then
0374            nonz=0
0375            return
0376          else
0377            print*,'STOP in DiagTmunu: negative energy density.  '
0378            stop'(3003200808)'
0379          endif
0380        endif
0381        do i=1,3
0382         Tv=0d0
0383         do k=1,3
0384          Tv=Tv+Tmunu(i,k)*vx(k)
0385         end do
0386         v(i)=(Tmunu(i,4)-Tv)/eps
0387        end do
0388        if(lrep.gt.60)then
0389         do i=1,3
0390          v(i)=0.5d0*(vx(i)+v(i))
0391         enddo
0392        endif
0393        !print*,'Tmunu: ',lrep,abs(eps-epsx),(abs(v(i)-vx(i)),i=1,3)
0394        err=1d-6
0395        if(lrep.gt.50)err=1d-5
0396        if(lrep.gt.89)err=1d-4
0397        if(abs(eps-epsx).lt.err.and.abs(v(1)-vx(1)).lt.err
0398      . .and.abs(v(2)-vx(2)).lt.err.and.abs(v(3)-vx(3)).lt.err)goto1
0399         do i=1,4
0400           w(i)=0
0401           do k=1,4
0402           w(i)=w(i)+Tmunu(i,k)*v(k)*sg(k)
0403           enddo
0404           w(i)=w(i)-eps*v(i)
0405         enddo
0406         if(lrep.gt.95
0407      ..and.w(1)*w(1)+w(2)*w(2)+w(3)*w(3)+w(4)*w(4).lt.err)goto1
0408       end do
0409
0410   1   v2=0.d0
0411       do i=1,3
0412        v2=v2+v(i)**2
0413       end do
0414       if(v2.ge.1.)stop'DiagTmunu: v2 ge 1.          '
0415       gamma=1./sqrt(abs(1.-v2))
0416       u(4)=gamma
0417       u(1)=v(1)*gamma
0418       u(2)=v(2)*gamma
0419       u(3)=v(3)*gamma
0420
0421       !~~~check
0422       g=gamma
0423       a=g*g/(g+1d0)
0424       Lor(4,4)=g
0425       do k=1,3
0426        Lor(4,k)=-g*v(k)
0427        Lor(k,4)=Lor(4,k)
0428       enddo
0429       do i=1,3
0430       do k=1,3
0431        Lor(i,k)=a*v(i)*v(k)
0432       enddo
0433       enddo
0434       do k=1,3
0435        Lor(k,k)=Lor(k,k)+1
0436       enddo
0437       do i=1,4
0438       do k=1,4
0439         tt(i,k)=0d0
0440         do m=1,4
0441         do n=1,4
0442           tt(i,k)=tt(i,k)+Lor(i,m)*Tmunu(m,n)*Lor(n,k)
0443         enddo
0444         enddo
0445       enddo
0446       enddo
0447       err=err*1d2
0448       if(tt(1,4).gt.err.or.tt(2,4).gt.err.or.tt(3,4).gt.err)then
0449         print*,'************** DiagTmunu: nonzero T14 or T24 or T34'
0450      .  ,'  after ',lrep,' iterations'
0451         print*,'d_eps or d_v(i): ',abs(eps-epsx),(abs(v(i)-vx(i)),i=1,3)
0452         print*,'Tmunu ( ix=',ix-nxico/2-1,'  iy=',iy-nyico/2-1,' ) :'
0453      .  ,'  iz=',iz-nzico/2-1
0454         do i=1,4
0455         write(*,'(4f10.5,2x,4f10.5)')
0456      .  (sngl(tmunu(i,k)),k=1,4),(sngl(tt(i,k)),k=1,4)
0457         enddo
0458       endif
0459
0460       end
0461