ReggeGribovPartonMCInterface/src/epos-jps.f

0001 c-----------------------------------------------------------------------
0002       subroutine jpsifo(npjpsi)
0003 c-----------------------------------------------------------------------
0004 c   forms a jpsi
0005 c-----------------------------------------------------------------------
0006       include 'epos.inc'
0007       include 'epos.incems'
0008       common/geom/rmproj,rmtarg,bmax,bkmx
0009       common/nucl3/phi,bimp
0010       parameter (ndep=129,ndet=129)
0011       common/cdep/xdep(ndep),qdep(ndep),wdep(ndep)
0012       common/cdet/xdet(ndet),qdet(ndet),wdet(ndet)
0013       parameter (nptj=129)
0014       common /cptj/xptj(nptj),qptj(nptj),wptj(nptj)
0015       parameter (mxbim=12)
0016       common/jpsi1/bimmax,kolran,delt,taumi,jpsinu,jpsidr,taudmx
0017       parameter (mxmass=20)
0018
0019       parameter (nxmdk=20)
0020       parameter (ntjpsi=150)
0021       common/jpsi7/xydens(ntjpsi,mxbim,nxmdk,nxmdk),a4min,a4max
0022       common/jpsi8/xys(mxbim,nxmdk,nxmdk),a5min,a5max
0023       common/jpsi9/ami(ntjpsi,mxmass),a6min,a6max
0024
0025       call utpri('jpsifo',ish,ishini,4)
0026       if(ish.ge.6)write(ifch,'(a)')' jpsi formation'
0027
0028 c     trigger
0029 c     -------
0030
0031       ymax=0.5
0032       ymin=-0.5
0033
0034 c     jpsi momenta
0035 c     ------------
0036       id=441
0037       call idmass(id,amass)
0038       s=amass**2
0039  2    rqptj=rangen()*qptj(nptj)
0040       pt=utinvt(nptj,xptj,qptj,rqptj)
0041       phi=2.*pi*rangen()
0042       px=pt*cos(phi)
0043       py=pt*sin(phi)
0044       lo=0
0045     1 lo=lo+1
0046       if(lo.gt.10)call utstop('jpsifo: lo > 10 &')
0047       z=0.19*sqrt(-2*alog(rangen()))*cos(2*pi*rangen()) !1-dim gauss
0048
0049
0050       if(z.gt.1.)goto 1
0051       pz=z*engy/2*ransig()
0052       e=sqrt(s+px**2+py**2+pz**2)
0053       amt=sqrt(amass**2+pt**2)
0054       y=sign(1.,pz)*alog( (e+abs(pz))/amt )
0055       if(y.lt.ymin.or.y.gt.ymax)goto 2
0056
0057 c     fill /cptl/
0058 c     -----------
0059       if(npjpsi.eq.0)then
0060         nptl=nptl+1
0061         npjpsi=nptl
0062       endif
0063       if(npjpsi.gt.mxptl)then
0064         print *,npjpsi,mxptl
0065         call utstop('jpsifo: npjpsi>mxptl&')
0066       endif
0067       istptl(npjpsi)=1
0068       idptl(npjpsi)=id
0069       pptl(1,npjpsi)=px
0070       pptl(2,npjpsi)=py
0071       pptl(3,npjpsi)=pz
0072       pptl(4,npjpsi)=e
0073       pptl(5,npjpsi)=amass
0074       kolran=1+rangen()*kolevt
0075       xorptl(1,npjpsi)=coord(1,kolran)
0076       xorptl(2,npjpsi)=coord(2,kolran)
0077       xorptl(3,npjpsi)=coord(3,kolran)
0078       xorptl(4,npjpsi)=coord(4,kolran)
0079       iorptl(npjpsi)=0
0080       jorptl(npjpsi)=0
0081       tivptl(1,npjpsi)=xorptl(4,npjpsi)
0082       tivptl(2,npjpsi)=ainfin
0083       ifrptl(1,npjpsi)=0
0084       ifrptl(2,npjpsi)=0
0085       if(ish.ge.6) then
0086         call alist("&",npjpsi,npjpsi)
0087         write (ifch,*) xorptl(1,npjpsi)
0088      $       ,xorptl(2,npjpsi),xorptl(3,npjpsi),xorptl(4,npjpsi)
0089      $       ,tivptl(1,npjpsi),tivptl(2,npjpsi)
0090         ii=iproj(kolran)
0091         jj=maproj+itarg(kolran)
0092         call alist("collision&",ii,ii)
0093         call alist("&",jj,jj)
0094       endif
0095       a4min=-15.
0096       a4max= 15.
0097       a5min=-15.
0098       a5max= 15.
0099       a6min= 2.
0100       a6max= 10.
0101       call utprix('jpsifo',ish,ishini,4)
0102       return
0103       end
0104
0105 c-----------------------------------------------------------------------
0106       function sptj(x)
0107 c-----------------------------------------------------------------------
0108 c     jpsi pt-distribution in 200 gev pp
0109 c-----------------------------------------------------------------------
0110       a=0.95
0111       c=1/0.363
0112       z=x/a
0113       sptj=1/a*c**c/utgam1(c)*z**(c-1)*exp(-c*z)
0114       return
0115       end
0116
0117 c-----------------------------------------------------------------------
0118       subroutine jpsian(ifirst)
0119 c-----------------------------------------------------------------------
0120 c     jpsi analysis.
0121 c-----------------------------------------------------------------------
0122       include 'epos.inc'
0123       include 'epos.incems'
0124       parameter (mxbim=12,ntjpsi=150,mxtauc=16)
0125       common/jpsi1/bimmax,kolran,delt,taumi,jpsinu,jpsidr,taudmx
0126       common/jpsi2/jjtot(mxbim),jjnuc(mxbim),jjjtau(mxbim,mxtauc)
0127       common/jpsi3/jjjtot(mxbim,ntjpsi),jjjdro(mxbim,ntjpsi)
0128       common/jpsi4/nnucl(mxbim,ntjpsi),nclose(mxbim,ntjpsi,3)
0129       common/jpsi5/ndrop(mxbim,ntjpsi),jjjnt(mxbim,mxtauc)
0130       parameter (mxmass=20,mxassy=20)
0131       common/jpsi6/ndrp2(mxbim,ntjpsi,mxmass,mxassy)
0132      &     ,ndrop3(mxbim,ntjpsi,mxmass,mxassy)
0133       parameter (nxmdk=20)
0134       common/jpsi7/xydens(ntjpsi,mxbim,nxmdk,nxmdk),a4min,a4max
0135       common/jpsi8/xys(mxbim,nxmdk,nxmdk),a5min,a5max
0136       common/jpsi9/ami(ntjpsi,mxmass),a6min,a6max
0137       common/jpsi10/ndrop0(mxbim,ntjpsi)
0138
0139       double precision tpro,zpro,ttar,ztar,ttaus,detap,detat
0140       common/cttaus/tpro,zpro,ttar,ztar,ttaus,detap,detat /ctain/mtain
0141       common/geom/rmproj,rmtarg,bmax,bkmx
0142       common/nucl3/phi,bimp
0143       parameter (ndep=129,ndet=129)
0144       common/cdep/xdep(ndep),qdep(ndep),wdep(ndep)
0145       common/cdet/xdet(ndet),qdet(ndet),wdet(ndet)
0146
0147       common/c9ptl/tauptl(mxptl),ss0ptl(mxptl)
0148
0149       call utpri('jpsian',ish,ishini,5)
0150
0151       detap=(ypjtl-yhaha)*etafac
0152       detat=-yhaha*etafac
0153       tpro=dcosh(detap)
0154       zpro=dsinh(detap)
0155       ttar=dcosh(detat)
0156       ztar=dsinh(detat)
0157
0158       jpsinu=1
0159       jpsidr=1
0160
0161 c      fac=1                     !    <-------- should be one finally
0162       taudmx=4
0163       rad=sqrt(0.62 / pi)
0164       taud=0
0165       nucia=0
0166       taumi=-2
0167       delt=0.1
0168       bimmax=amin1(rmproj+rmtarg,bmaxim)
0169       delbim=bimmax/mxbim
0170       ii=iproj(kolran)
0171       jj=maproj+itarg(kolran)
0172
0173
0174 c.....first event: delete commom blocks...............................
0175       if(ifirst.eq.1)then
0176         ifirst=0
0177         do nbim=1,mxbim
0178           jjtot(nbim)=0
0179           jjnuc(nbim)=0
0180           do ix=1,nxmdk
0181             do iy=1,nxmdk
0182               xys(nbim,ix,iy)=0.
0183             enddo
0184           enddo
0185           do nt=1,ntjpsi
0186             jjjtot(nbim,nt)=0
0187             nnucl(nbim,nt)=0
0188             jjjdro(nbim,nt)=0
0189             ndrop(nbim,nt)=0
0190             ndrop0(nbim,nt)=0
0191             do kk=1,3
0192               nclose(nbim,nt,kk)=0
0193             enddo
0194             do ix=1,nxmdk
0195               do iy=1,nxmdk
0196                 xydens(nt,nbim,ix,iy)=0.
0197               enddo
0198             enddo
0199           enddo
0200           do mm=1,mxtauc
0201             jjjtau(nbim,mm)=0
0202           enddo
0203         enddo
0204       endif
0205
0206       nbim=1+int(bimevt/delbim)
0207       if(nbim.lt.0.or.nbim.gt.mxbim) goto 5
0208       jjtot(nbim)=jjtot(nbim)+1 !events pro bin
0209
0210       do 1 i=1,nptl
0211         if(idptl(i).eq.441)j=i
0212  1    continue
0213
0214 c     if(jpsidr.eq.1)then
0215 c     write(6,'(a,i5,a,f6.2,a,f6.2,a,f6.2)')'ip=',ip
0216 c     *,'   rin= '
0217 c     *,'   mass= ',pptl(5,ip)
0218 c     enddo
0219 c     endif
0220 c     endif
0221
0222       taumax=0.
0223       ttaus=taumi-delt
0224       do 2 nt=1,ntjpsi
0225         idrin=0
0226         ttaus=ttaus+delt        !increment of time
0227         if(ish.ge.6)write(ifch,*) 'ttaus:-->',ttaus,ii,jj
0228         jpsiex=1
0229         call jtain(j,xj,yj,zj,tj,n,1)
0230         if(n.eq.1.or.n.eq.2.or.n.eq.9)jpsiex=0 !goto 2
0231         if(jpsiex.eq.1)jjjtot(nbim,nt)=jjjtot(nbim,nt)+1
0232
0233 c       nucleons
0234 c       --------
0235
0236         if(jpsinu.eq.1.and.jpsiex.eq.1)then !test jpsi-nucleon collision
0237           do 6 i=1,maproj+matarg
0238             if(i.eq.ii.or.i.eq.jj)goto 6
0239             nnucl(nbim,nt)=nnucl(nbim,nt)+1
0240             t=sngl(ttaus)
0241             x=xorptl(1,i)+(t-xorptl(4,i))*pptl(1,i)/pptl(4,i)
0242             y=xorptl(2,i)+(t-xorptl(4,i))*pptl(2,i)/pptl(4,i)
0243             z=xorptl(3,i)+(t-xorptl(4,i))*pptl(3,i)/pptl(4,i)
0244             pde=(pptl(3,i)+pptl(3,j))/(pptl(4,i)+pptl(4,j))
0245             gam2i=1-pde**2
0246             if(gam2i.eq.0.)goto 6
0247             dist=sqrt((x-xj)**2+(y-yj)**2
0248      &           +1/gam2i*(z-zj-(t-tj)*pde)**2)
0249             if(dist.le.rad)then
0250               nclose(nbim,nt,1)=nclose(nbim,nt,1)+1
0251               nucia=1
0252               if(ish.ge.6)then
0253                 write (ifch,*) "nucl dist:",dist,' dist(sig)='
0254      $               ,rad
0255                 call alist("&",i,i)
0256                 call alist("&",j,j)
0257               endif
0258             elseif(dist.le.rad+1)then
0259               nclose(nbim,nt,2)=nclose(nbim,nt,2)+1
0260             elseif(dist.le.rad+3)then
0261               nclose(nbim,nt,3)=nclose(nbim,nt,3)+1
0262             endif
0263  6        continue
0264         endif
0265
0266 c       particles
0267 c       ---------
0268         do 8 i=1,nptl
0269 c          if ( i.eq.ii.or.i.eq.jj ) goto 8
0270           call jtain(i,x,y,z,t,n,1)
0271           call jtaus(z,tz,sz)
0272           if(n.eq.1.or.n.eq.2.or.n.eq.9)goto 8
0273 c         calculate s
0274 c         -----------
0275           iad=iabs(idptl(i))
0276 c          s=(pptl(4,i)+pptl(4,j))**2-(pptl(3,i)+pptl(3,j))**2
0277 c     $         -(pptl(2,i)+pptl(2,j))**2-(pptl(1,i)+pptl(1,j))**2
0278           if ( iad.eq.120 .or. iad.eq.110 ) then !pion
0279             sig=1.              !
0280           elseif ( iad.eq.121 .or. iad.eq.111 ) then ! rho
0281             sig=1.              !
0282           elseif ( iad.eq.1120 .or. iad.eq.1220 ) then
0283             sig=3.0             ! ???? or 6 ????
0284           else
0285             goto 8
0286           endif
0287           call jtaus(zj,tzj,szj)         !????????????????? OK ?
0288           dist=sqrt((x-xj)**2+(y-yj)**2+(sz-szj)**2)
0289           if ( dist .lt. sqrt(0.1*sig/pi) ) then
0290             istptl(j)=2
0291             if(ish.ge.6)then
0292               write (ifch,*) "dist:",dist,' dist(sig)='
0293      $             ,sqrt(0.1*sig/pi),' sig=',sig
0294               call alist("&",i,i)
0295               call alist("&",j,j)
0296             endif
0297           endif
0298  8      continue
0299
0300 c     droplets
0301 c     --------
0302
0303         if(jpsidr.eq.1)then
0304           call jtaus(zj,tzj,szj)
0305           do 3 i=maproj+matarg+1,nptl
0306 c...........x-y distribution of strings..............................
0307             if(istptl(i).eq.29.and.nt.eq.1)then
0308               call jtain(i,x,y,z,t,n,1)
0309               if(x.gt.a5min.and.x.lt.a5max.and.
0310      &             y.gt.a5min.and.y.lt.a5max)then
0311                 ix=(x-a5min)/(a5max-a5min)*nxmdk + 1
0312                 iy=(y-a5min)/(a5max-a5min)*nxmdk + 1
0313                 xys(nbim,ix,iy)=xys(nbim,ix,iy)+pptl(5,i)
0314               endif
0315             endif
0316             if(istptl(i).gt.10)goto 3
0317             if(i.eq.j)goto 3
0318
0319 c...................................................................
0320
0321             call jtain(i,x,y,z,t,n,1)
0322             if(n.eq.1.or.n.eq.2.or.n.eq.9)goto 3
0323             stop'jpsian: change!!!!        ' !call jintep(i,x,y,z,t,sz,eps,rho)
0324             if(eps.lt.aouni)goto 3 !min-dichte
0325             ndrop(nbim,nt)=ndrop(nbim,nt)+1 !droplets at time nt
0326             ndrop0(nbim,nt)=ndrop0(nbim,nt)+pptl(5,i) !mass
0327             des=0 !?????????????????????????????????
0328             o=sz+des
0329             u=sz-des
0330             r=0 !( xxxx(i) +sngl(ttaus) ) *fac
0331 c..............assym-mass-distribution...............................
0332             assym=log(des/r)
0333             amass=pptl(5,i)
0334             a1min=-5.
0335             a1max=5.
0336             a2min=0.
0337             a2max=40.
0338             if(assym.ge.a1min.and.assym.lt.a1max
0339      &           .and.amass.ge.a2min.and.amass.lt.a2max
0340      &           ) then
0341               nassym=(assym-a1min)/(a1max-a1min)*mxassy+1
0342               namass=(amass-a2min)/(a2max-a2min)*mxmass+1
0343               ndrp2(nbim,nt,namass,nassym)=
0344      &             ndrp2(nbim,nt,namass,nassym)+1
0345             endif
0346
0347 c..............vol-mass-distribution...............................
0348             a3min=-2.
0349             a3max=3.
0350             v=log(pi*r**2.*2.*des)/log(10.)
0351             if(v.ge.a3min.and.v.lt.a3max
0352      &           .and.amass.ge.a2min.and.amass.lt.a2max
0353      &           ) then
0354               nv=(v-a3min)/(a3max-a3min)*mxassy+1
0355               namass=(amass-a2min)/(a2max-a2min)*mxmass+1
0356               ndrop3(nbim,nt,namass,nv)=
0357      &             ndrop3(nbim,nt,namass,nv)+1
0358             endif
0359 c..............x-y distribution of droplet..............................
0360             ix=(x-a4min)/(a4max-a4min)*nxmdk + 1
0361             iy=(y-a4min)/(a4max-a4min)*nxmdk + 1
0362             xydens(nt,nbim,ix,iy)=xydens(nt,nbim,ix,iy)+eps
0363
0364             if(jpsiex.eq.0)goto 3
0365
0366 c           if(mod(nt,10).eq.1)
0367 c           write(6,'(f5.2,i6,5x,3f7.2,5x,4f6.2)')sngl(ttaus),i,
0368 c           *              u,szj,o,sqrt((x-xj)**2+(y-yj)**2),r,v,eps
0369
0370             if(szj.lt.u.or.szj.gt.o)goto 3
0371             if((x-xj)**2+(y-yj)**2.gt.r**2)goto 3
0372
0373 c           write(6,'(a,f5.2,a,i5,a)')'***** t=',sngl(ttaus)
0374 c           *,' -- jpsi in droplet ',i,' *****'
0375
0376             taud=taud+delt
0377             taumax=max(taud,taumax)
0378
0379 c           write (*,*) taud,taumax
0380
0381             idrin=1
0382             jjjdro(nbim,nt)=jjjdro(nbim,nt)+1
0383  3        continue
0384
0385 c     if (idrin.ne.1)taud=max(taud-delt,0.)
0386
0387           if (idrin.ne.1)taud=0.
0388         endif
0389
0390  2    continue                  !end nt-loop
0391
0392       ijmod=2
0393       if(ijmod.eq.2)taud=taumax
0394       if(nucia.eq.1)jjnuc(nbim)=jjnuc(nbim)+1
0395       if(taud.gt.0.)then
0396         do ntaud=1,mxtauc
0397           tauc=ntaud*taudmx/mxtauc
0398           if(taud.gt.tauc)jjjtau(nbim,ntaud)=jjjtau(nbim,ntaud)+1
0399           if(nucia.eq.1.or.taud.gt.tauc)
0400      &         jjjnt(nbim,ntaud)=jjjnt(nbim,ntaud)+1
0401         enddo
0402       else
0403         do ntaud=1,mxtauc
0404           if(nucia.eq.1)
0405      &         jjjnt(nbim,ntaud)=jjjnt(nbim,ntaud)+1
0406         enddo
0407       endif
0408
0409     5 continue
0410       call utprix('jpsian',ish,ishini,5)
0411       return
0412       end
0413
0414 c-----------------------------------------------------------------------
0415       subroutine jtauan(is,im)
0416 c-----------------------------------------------------------------------
0417 c     display collision
0418 c     im = ijk
0419 c          k > 0   -->  post-script
0420 c          j > 0   -->  povray
0421 c            j = 1  -->  time and z in n-n cms
0422 c            j = 2  -->  time and z on hyberbola
0423 c          i > 0   -->  text  ( changes in alist per time step )
0424 c     cut in zeta-x (or y) plane for tau
0425 c-----------------------------------------------------------------------
0426       include 'epos.inc'
0427       double precision tpro,zpro,ttar,ztar,ttaus,detap,detat
0428       common/cttaus/tpro,zpro,ttar,ztar,ttaus,detap,detat /ctain/mtain
0429       parameter (mxbim=12,ntjpsi=150,mxtauc=16)
0430       common/jpsi1/bimmax,kolran,delt,taumi,jpsinu,jpsidr,taudmx
0431       common/jpsi2/jjtot(mxbim),jjnuc(mxbim),jjjtau(mxbim,mxtauc)
0432       common/jpsi3/jjjtot(mxbim,ntjpsi),jjjdro(mxbim,ntjpsi)
0433       common/jpsi4/nnucl(mxbim,ntjpsi),nclose(mxbim,ntjpsi,3)
0434       common/jpsi5/ndrop(mxbim,ntjpsi),jjjnt(mxbim,mxtauc)
0435       parameter (mxmass=20,mxassy=20)
0436       common/jpsi6/ndrp2(mxbim,ntjpsi,mxmass,mxassy)
0437      &     ,ndrop3(mxbim,ntjpsi,mxmass,mxassy)
0438       parameter (nxmdk=20)
0439       common/jpsi7/xydens(ntjpsi,mxbim,nxmdk,nxmdk),a4min,a4max
0440       common/jpsi8/xys(mxbim,nxmdk,nxmdk),a5min,a5max
0441       common/jpsi9/ami(ntjpsi,mxmass),a6min,a6max
0442       common/jpsi10/ndrop0(mxbim,ntjpsi)
0443       character*20 name,nnrr
0444       character*28 filename
0445       character*12 color(20)
0446       character*12 colpo(20)
0447       logical lcalc!,zet
0448       dimension isch(mxptl)
0449 c      zet=.true.
0450 c      zet=.false.
0451       xmin=-10.
0452       xmax=10.
0453       zmin=-10.
0454       zmax=10.
0455
0456 c      zevent=float(nevent*jpsi)
0457       if(mod(im,10).ne.0)then
0458         name='tau-'
0459         n=1
0460  5      l=4
0461         ll=int(log(real(n))/log(10.))+1
0462         do ii=ll,1,-1
0463           l=l+1
0464           name(l:l)=char(48+mod(int(n/10**(ii-1)),10))
0465         enddo
0466         name(l+1:l+3)='.ps'
0467         l=l+3
0468         inquire(file=name(1:l),exist=lcalc)
0469         if (lcalc)then
0470           n=n+1
0471           goto 5
0472         endif
0473         write(*,*) 'jtauan name ',name
0474         ifps=92
0475         open(unit=ifps,file=name(1:l),status='unknown')
0476         WRITE(ifps,'(a)') '%!PS-Adobe-2.0'
0477         WRITE(ifps,'(a)') '%%Title: tt2.fig'
0478         WRITE(ifps,'(a)') '%%Orientation: Portrait'
0479         WRITE(ifps,'(a)') '%%BeginSetup'
0480         WRITE(ifps,'(a)') '%%IncludeFeature: *PageSize A4'
0481         WRITE(ifps,'(a)') '%%EndSetup'
0482         WRITE(ifps,'(a)') '%%EndComments'
0483         WRITE(ifps,*) '/l {lineto} bind def'
0484         WRITE(ifps,*) '/rl {rlineto} bind def'
0485         WRITE(ifps,*) '/m {moveto} bind def'
0486         WRITE(ifps,*) '/rm {rmoveto} bind def'
0487         WRITE(ifps,*) '/s {stroke} bind def'
0488         WRITE(ifps,*) '/gr {grestore} bind def'
0489         WRITE(ifps,*) '/gs {gsave} bind def'
0490         WRITE(ifps,*) '/cp {closepath} bind def'
0491         WRITE(ifps,*) '/tr {translate} bind def'
0492         WRITE(ifps,*) '/sc {scale} bind def'
0493         WRITE(ifps,*) '/sd {setdash} bind def'
0494         WRITE(ifps,*) '/sdo {[.01 .05] 0 sd} bind def'
0495         WRITE(ifps,*) '/sdf {[1 .0] 0 sd} bind def'
0496         WRITE(ifps,*) '/n {newpath} bind def'
0497         WRITE(ifps,*) '/slw {setlinewidth } bind def'
0498         write(ifps,*) '/srgb {setrgbcolor} bind def'
0499         write(ifps,*) '/lgrey      { 0 0.95 0.95 srgb} bind def'
0500         write(ifps,*) '/black      { 0 0 0 srgb} bind def'
0501         write(ifps,*) '/red        { 1 0 0 srgb} bind def  '
0502         write(ifps,*) '/green      { 0 1 0  srgb} bind def  '
0503         write(ifps,*) '/blue       { 0 0 1  srgb} bind def  '
0504         write(ifps,*) '/yellow     { 1 0.5 0  srgb} bind def  '
0505         write(ifps,*) '/turquoise  { 0 1 1  srgb} bind def  '
0506         write(ifps,*) '/purple     { 1 0 1  srgb} bind def  '
0507 c.......write(ifps,*) '/  {   srgb} bind def  '
0508 c.......write(ifps,*) '/  {   srgb} bind def  '
0509         write(ifps,*) '/ef {eofill} bind def'
0510         WRITE(ifps,'(a)') '%%EndProlog'
0511         WRITE(ifps,*) 'gsave'
0512         WRITE(ifps,*) '/Helvetica findfont 10 scalefont setfont'
0513       endif
0514       color(9)='lgrey     '
0515       color(1)='black     '
0516       color(2)='red       '
0517       color(3)='green     '
0518       color(4)='blue      '
0519       color(7)='yellow    '
0520       color(5)='turquoise '
0521       color(6)='purple    '
0522       colpo(1)='Red  '
0523       colpo(2)='Green  '
0524       colpo(3)='Blue  '
0525       colpo(4)='Yellow  '
0526       colpo(5)='Cyan '
0527       colpo(6)='Magenta  '
0528       colpo(7)='Black  '
0529       colpo(8)='Aquamarine '
0530
0531
0532       do i=1,mxptl
0533         isch(i)=0
0534       enddo
0535
0536 c      gray0=.95
0537 c      gray1=0.
0538       iyb=0                     !????????????????????
0539       np=0
0540       nt=-10
0541       deltau=0.1
0542       taumin=-1
0543       ttaus=0
0544       do while (ttaus.lt.20.)
0545         nt=nt+1
0546        ! ttaus=dble(taumin+deltau*(factau**(1.*nt-1.)-1)/(factau-1.))
0547         ttaus=taumin+deltau*nt
0548         tau=ttaus
0549         np=np+1
0550         if(mod(im,10).ne.0)then
0551           write(ifps,'(a,i4)') '%%Page: number ',np
0552           write(ifps,'(a)') 'gsave'
0553           WRITE(ifps,*) '100 700 tr'
0554           scale=0.125
0555           WRITE(ifps,*) 1./scale,1./scale,' sc'
0556           WRITE(ifps,*) scale/2.,' slw'
0557           WRITE(ifps,*) '/Helvetica findfont ',15.*scale
0558      &         ,' scalefont setfont'
0559           write(ifps,*) color(1),' n ',zmin,xmin,' m ( tau:'
0560      $         ,tau,') show '
0561           WRITE(ifps,*) '/Helvetica findfont ',2.*scale
0562      &         ,' scalefont setfont'
0563         endif
0564
0565 *--------------------------------------------------------------------*
0566 *------ povray ------------------------------------------------------*
0567 *--------------------------------------------------------------------*
0568
0569         if (mod(im/100,10).gt.0) then
0570           write (ifch,*)   "-----",np,", tau:",ttaus,"------"
0571         endif
0572         if (mod(im/10,10).gt.0) then
0573           write (nnrr,'(i5)') np
0574           li=6-log(1.*np+0.1)/log(10.)
0575           write (*,*)   "--->"//nnrr(li:5)//"<-----",li,ttaus
0576           ifpo=55
0577           filename="tau."//nnrr(li:5)//".pov"
0578           open(unit=ifpo,file=filename,status='unknown')
0579           write (ifpo,'(a)') '#include "colors.inc";'
0580 c         write (ifpo,'(a)') '#include "shapes.inc" '
0581 c         write (ifpo,'(a)') '#include "textures.inc" '
0582           write (ifpo,'(a)') 'background {color White} '
0583           write (ifpo,'(a)') 'camera {location <0,0,-120> '
0584           write (ifpo,'(a)') '     direction <0,0,2> look_at <0,0,0>} '
0585           write (ifpo,'(a)') 'light_source{<0,300,0> color White} '
0586           write (ifpo,'(a)') 'light_source{<0,5,-90> color White} '
0587           write (ifpo,'(a)') ' '
0588           write (ifpo,'(a)') ' '
0589         endif
0590         do i=1,nptl
0591           if (istptl(i).gt.1) goto 123
0592             if((tivptl(2,i)-tivptl(1,i)).lt.1e-3
0593      $           .and.idptl(i).gt.1000000.and.iyb.eq.0)
0594      $           then
0595               write (*,*) 'tiv1=tiv2 !!!!!!!!',i
0596               tivptl(2,i)=tivptl(1,i)+100.
0597             endif
0598 c...........calculate coordinates....................................
0599             if(mod(im/10,10).eq.1) then !n-n cms frame
0600               if (istptl(i).gt.1
0601      $             .or.ttaus.lt.tivptl(1,i)
0602      $             .or.ttaus.gt.tivptl(2,i)) goto 123
0603               x=xorptl(1,i)+(ttaus-xorptl(4,i))*pptl(1,i)/pptl(4,i)
0604               y=xorptl(2,i)+(ttaus-xorptl(4,i))*pptl(2,i)/pptl(4,i)
0605               z=xorptl(3,i)+(ttaus-xorptl(4,i))*pptl(3,i)/pptl(4,i)
0606             else                !  hyperbola frame
0607               call jtain(i,x,y,z,t,n,0)
0608               call jtaus(z,tz,sz)
0609               z=sz
0610               if(n.ne.0) goto 123
0611             endif
0612 c...........plot sphere or cylinder ................................
0613             if(idptl(i).gt.700000000)
0614      $           then
0615               if(mod(im/10,10).eq.1)then
0616
0617               else
0618                 des=0 !?????????????????????????????
0619                 r=0   !(xxxx(i)+vrad*sngl(ttaus))
0620                 o=sz+des
0621                 u=sz-des
0622                 print *,ttaus,o,u,r,x,y
0623                 ic=4
0624                 if (mod(im/10,10).gt.0) then
0625                   write (ifpo,111) o,x,y,u,x,y,r,colpo(ic)
0626                 endif
0627               endif
0628 c.............text output of changes in time step ...................
0629               if (mod(im/100,10).gt.0) then
0630                 if(isch(i).eq.0)then
0631                   write (ifch,'("> ",$)')
0632                   call alist("&",i,i)
0633                   isch(i)=1
0634                 endif
0635               endif
0636 c$$$              o=sz+des
0637 c$$$              u=sz-des
0638 c$$$              r=(xxxx(i)+vrad*sngl(ttaus))
0639 c$$$              rr2=r**2-(y-yb)**2
0640 c$$$              if(rr2.gt.0.)then
0641 c$$$                write(ifps,*)
0642 c$$$     &               ,' n ',u,x-r,' m ',o,x-r,' l '
0643 c$$$     &               ,o,x+r,' l ',u,x+r,' l cp s '
0644 c$$$                write(ifps,*) ' n ',u,x-r,' m (',i,iorptl(i),') show '
0645 c$$$              endif
0646             else
0647 c.............cylinder................................................
0648               r=0.8
0649               ic=1
0650               if(abs(idptl(i)).lt.999) r=0.5
0651               if(iabs(idptl(i)).eq.1120) ic=2
0652               if(iabs(idptl(i)).eq.1220) ic=3
0653               if(iabs(idptl(i)).eq.441) ic=5
0654               if(mod(im/10,10).gt.0)then
0655                 write (ifpo,110) z,x,y,r,colpo(ic) ! sphere
0656               endif
0657 c.............text...................................................
0658               if(mod(im/100,10).gt.0)then
0659                 if(isch(i).eq.0)then
0660                   write (ifch,'("> ",$)')
0661                   call alist("&",i,i)
0662                   isch(i)=1
0663                 endif
0664               endif
0665             endif
0666             goto 124
0667  123      continue
0668 c.........text................................
0669           if(mod(im/100,10).gt.0)then
0670             if(isch(i).eq.1)then
0671               write (ifch,'("< ",$)')
0672               call alist("&",i,i)
0673               isch(i)=0
0674             endif
0675           endif
0676  124      continue
0677         enddo
0678 c.......
0679  110    format('sphere {<',G12.6,',',g12.6,',',g12.6,'>,',g12.6
0680      $       ,'pigment {color ',a,'}}')
0681  111    format('cylinder {<',
0682      $       G12.6,',',g12.6,',',g12.6,
0683      $       '>,<',
0684      $       G12.6,',',g12.6,',',g12.6,
0685      $       '>,',
0686      $       g12.6,
0687      $       'pigment {color ',a,'}}')
0688         if(mod(im/10,10).gt.0)then
0689           close(unit=ifpo)
0690         endif
0691 *-------------------------------------------------------------------*
0692 *-------   end   povray   ------------------------------------------*
0693 *-------   begin post-script ---------------------------------------*
0694 *-------------------------------------------------------------------*
0695
0696         if(mod(im,10).eq.0) goto 159
0697         yb=-6.
0698         dy=12./12.
0699         yb=yb-dy/2
0700         do iyb=0,11
0701           yb=yb+dy
0702           WRITE(ifps,*) 'gsave'
0703           WRITE(ifps,*) (xmax-xmin)*1.1*float(int(iyb/4))
0704      &         ,-(xmax-xmin)*1.1*mod(iyb,4),' tr'
0705           write(ifps,*) ' n ',zmin,xmin,' m ',zmax,xmin,' l '
0706      &         ,zmax,xmax,' l ',zmin,xmax,' l cp s '
0707 c          ttaus=dble(tau)
0708 c.........particles in layer iyb.............
0709           do 10 i=1,nptl
0710             if (istptl(i).gt.1) goto 10
0711             if((tivptl(2,i)-tivptl(1,i)).lt.1e-3
0712      $           .and.idptl(i).gt.1000000.and.iyb.eq.0)
0713      $           then
0714               write (*,*) 'tiv1=tiv2 !!!!!!!!',i
0715               tivptl(2,i)=tivptl(1,i)+100.
0716             endif
0717             call jtain(i,x,y,z,t,n,0)
0718             call jtaus(z,tz,sz)
0719             if(n.ne.0) goto 10
0720             if(
0721      $           (is.eq.0.or.is.eq.i.or.is.eq.iorptl(i)))then
0722 *.............
0723 *.............   is  is the particle number to observe
0724 *.............   if is=0 then all particles
0725 *.............
0726 c              .and.abs(y-yb).lt.dy/2)then
0727               des=0 !??????????????????????????????????
0728               if(iyb.eq.11.and
0729      $             .abs(tivptl(2,i)-tivptl(1,i)-100.).le.1e-4 ) then
0730                 tivptl(2,i)=tivptl(1,i)+0.01
0731               endif
0732               o=sz+des
0733               u=sz-des
0734               r=0   !(xxxx(i)+vrad*sngl(ttaus))
0735               rr2=r**2-(y-yb)**2
0736               if(rr2.gt.0.)then
0737                 r=sqrt(rr2)
0738 c                write (*,*) i,des,o,u,r,y
0739                 write(ifps,*)  color(mod(i,5)+2)
0740      &               ,' n ',u,x-r,' m ',o,x-r,' l '
0741      &               ,o,x+r,' l ',u,x+r,' l cp s '
0742                 write(ifps,*) ' n ',u,x-r,' m (',i,iorptl(i),') show '
0743               endif
0744             elseif(abs(y-yb).lt.dy/2.and.zmin.lt.sz.and.sz.lt.zmax
0745      &             .and.xmin.lt.x.and.x.lt.xmax)then
0746               r=0.8
0747               ic=1
0748               if(abs(idptl(i)).lt.999)r=0.5
0749               if(abs(idptl(i)).lt.999)ic=2
0750               if(abs(idptl(i)).eq.1120)ic=3
0751               if(abs(idptl(i)).eq.1220)ic=4
0752               if(idptl(i).eq.441) ic=7
0753
0754               io=iorptl(i)
0755               if(is.eq.0.or.io.eq.is)then
0756                 write(ifps,*) ' n ',sz,x,r,0,360,' arc ',color(ic),' s '
0757                 write(ifps,*) ' n ',sz-r,x,' m (',i,io,') show '
0758               endif
0759             endif
0760  10       continue
0761           write(ifps,*) color(1),' n ',zmin,xmin,' m (',yb,') show '
0762           WRITE(ifps,*) 'grestore'
0763         enddo                   !yb bin
0764         write(ifps,'(a)') 'grestore'
0765         write(ifps,*) 'showpage'
0766  159    continue
0767       enddo
0768
0769
0770 c        write(ifps,*) ' n ',y0,x0,' m ',y1,x0,' l ',y1,x1,' l '
0771 c     &    ,y0,x1,' l cp s '
0772 c        write(ifps,*) ' n ',(y0+y1)/2-10.*scale,(x0+x1)/2-5.*scale
0773 c     &    ,' m (',ii,jj,') show '
0774 c        ii=nob-1
0775
0776
0777       if(mod(im,10).ne.0)then
0778         write(ifps,*) 'gr'
0779         write(ifps,'(a)') '%%Trailer'
0780         write(ifps,'(a,i4)') '%%Pages: ',np
0781         write(ifps,'(a)') '%%EOF'
0782         close(unit=ifps)
0783       endif
0784       return
0785       end
0786
0787
0788 c-----------------------------------------------------------------------
0789       subroutine jpsihi
0790 c-----------------------------------------------------------------------
0791 c     histogram
0792 c-----------------------------------------------------------------------
0793       include 'epos.inc'
0794       parameter (mxbim=12,ntjpsi=150,mxtauc=16)
0795       common/jpsi1/bimmax,kolran,delt,taumi,jpsinu,jpsidr,taudmx
0796       common/jpsi2/jjtot(mxbim),jjnuc(mxbim),jjjtau(mxbim,mxtauc)
0797       common/jpsi3/jjjtot(mxbim,ntjpsi),jjjdro(mxbim,ntjpsi)
0798       common/jpsi4/nnucl(mxbim,ntjpsi),nclose(mxbim,ntjpsi,3)
0799       common/jpsi5/ndrop(mxbim,ntjpsi),jjjnt(mxbim,mxtauc)
0800       parameter (mxmass=20,mxassy=20)
0801       common/jpsi6/ndrp2(mxbim,ntjpsi,mxmass,mxassy)
0802      &     ,ndrop3(mxbim,ntjpsi,mxmass,mxassy)
0803       parameter (nxmdk=20)
0804       common/jpsi7/xydens(ntjpsi,mxbim,nxmdk,nxmdk),a4min,a4max
0805       common/jpsi8/xys(mxbim,nxmdk,nxmdk),a5min,a5max
0806       common/jpsi9/ami(ntjpsi,mxmass),a6min,a6max
0807       common/jpsi10/ndrop0(mxbim,ntjpsi)
0808
0809       zevent=float(nevent*jpsi)
0810
0811       write(ifhi,'(a)')    'cd /users/theoric/werner/histo/newdata'
0812       write(ifhi,'(a)')       'newpage'
0813
0814 c     suppression as a function of b
0815 c     ------------------------------
0816
0817       write(ifhi,'(a)')       'zone 1 2 1 openhisto'
0818       write(ifhi,'(a)')       'htyp lin xmod lin ymod lin'
0819       write(ifhi,'(a)')       'text 0 0 "xaxis bmax-b (fm)"'
0820       write(ifhi,'(a)')       'text 0 0 "yaxis J(b) et Jnuc(b) / J"'
0821       write(ifhi,'(a,2e11.3)')'xrange',0.,bimmax
0822       write(ifhi,'(3a)')'yrange',' 0 ',' auto'
0823       write(ifhi,'(a)')       'columnweight 4 column c4 = ( 0 ) '
0824       write(ifhi,'(a)')       'array 4'
0825       do j=mxbim,1,-1
0826       bim=(j-0.5)*bimmax/mxbim
0827       write(ifhi,'(4e12.4)')bimmax-bim,jjtot(j)/zevent,0.,zevent
0828       enddo
0829       write(ifhi,'(a)')       '  endarray'
0830       write(ifhi,'(a)')       'closehisto plot 0-'
0831
0832       write(ifhi,'(a)')       'openhisto'
0833       write(ifhi,'(a)')       'htyp lin xmod lin ymod lin'
0834       write(ifhi,'(a)')       'text 0 0 "xaxis bmax-b (fm)"'
0835       write(ifhi,'(a)') 'text 0 0 " "'
0836       write(ifhi,'(a,2e11.3)')'xrange',0.,bimmax
0837       write(ifhi,'(3a)')'yrange',' 0 ',' auto'
0838       write(ifhi,'(a)')       'columnweight 4 column c4 = ( 0 ) '
0839       write(ifhi,'(a)')       'array 4'
0840       do j=mxbim,1,-1
0841       bim=(j-0.5)*bimmax/mxbim
0842       write(ifhi,'(4e12.4)')bimmax-bim,jjnuc(j)/zevent,0.,zevent
0843       enddo
0844       write(ifhi,'(a)')       '  endarray'
0845       write(ifhi,'(a)')       'closehisto plot 0'
0846
0847
0848       write(ifhi,'(a)')       'openhisto'
0849       write(ifhi,'(a)')       'htyp lin xmod lin ymod lin'
0850       write(ifhi,'(a)')       'text 0 0 "xaxis bmax-b (fm)"'
0851       write(ifhi,'(a)')       'text 0 0 "yaxis survival ratio"'
0852       write(ifhi,'(a,3e11.3)')'xrange',0.,bimmax
0853       write(ifhi,'(3a)')'yrange',' 0.2 ',' auto '
0854       write(ifhi,'(a)')       'columnweight 4 column c4 = ( 0 ) '
0855       write(ifhi,'(a)')       'array 4'
0856       do j=mxbim,1,-1
0857         bim=(j-0.5)*bimmax/mxbim
0858         rat=0
0859         if(jjtot(j).gt.0.)rat=float(jjtot(j)-jjnuc(j))/jjtot(j)
0860         write(ifhi,'(4e12.4)')bimmax-bim,rat,0.,float(jjtot(j))
0861       enddo
0862       write(ifhi,'(a)')       '  endarray'
0863       if(maproj.eq.208.and.matarg.eq.208)then
0864         write(ifhi,'(a)')       'closehisto plot 0-'
0865         write(ifhi,'(a)')       'openhisto'
0866         write(ifhi,'(a)')       'set fmsc 1.0'
0867         write(ifhi,'(a,f4.1,a)')'column c1 = ( ',bimmax,' - c1  )'
0868         write(ifhi,'(a)')       'column c2 = ( c2 * 0.02 )'
0869         write(ifhi,'(a)')       'input na50 ratio-b plot 0'
0870       else
0871         write(ifhi,'(a)')       'closehisto plot 0'
0872       endif
0873
0874 c     nr of jpsi vs t
0875 c     ---------------
0876
0877       write(ifhi,'(a)')       'zone 3 4 1'
0878       do nb=mxbim,1,-1
0879         bim=(nb-0.5)*bimmax/mxbim
0880         write(ifhi,'(a)')       'openhisto'
0881         write(ifhi,'(a)')       'htyp lin xmod lin ymod lin'
0882         write(ifhi,'(a,f5.2,a)')'text .1 .90 "b= ',bim,' fm"'
0883         write(ifhi,'(a)')       'text 0 0 "xaxis time t (fm)"'
0884         write(ifhi,'(a)')       'text 0 0 "yaxis J(b,t) / J"'
0885         write(ifhi,'(a,2e11.3)')'xrange',taumi,taumi+ntjpsi*delt
0886         write(ifhi,'(3a)')      'yrange',' 0 ',' auto'
0887         write(ifhi,'(a)')       'columnweight 4 column c4 = ( 0 ) '
0888         write(ifhi,'(a)')       'array 4'
0889         do j=1,ntjpsi
0890           tau=taumi+(j-0.5)*delt
0891           write(ifhi,'(4e12.4)')tau,float(jjjtot(nb,j))/nevent,0.,nevent
0892         enddo
0893         write(ifhi,'(a)')       '  endarray'
0894         write(ifhi,'(a)')       'closehisto plot 0'
0895       enddo
0896
0897 c     nr of nucleons vs t
0898 c     -------------------
0899
0900       if(jpsinu.eq.1)then
0901         write(ifhi,'(a)')       'zone 3 4 1'
0902         do nb=mxbim,1,-1
0903           bim=(nb-0.5)*bimmax/mxbim
0904           write(ifhi,'(a)')       'openhisto'
0905           write(ifhi,'(a)')       'htyp lin xmod lin ymod lin'
0906           write(ifhi,'(a,f5.2,a)')'text .1 .90 "b= ',bim,' fm"'
0907           write(ifhi,'(a)')       'text 0 0 "xaxis time t (fm)"'
0908           write(ifhi,'(a)')       'text 0 0 "yaxis N(b,t) / J"'
0909           write(ifhi,'(a,2e11.3)')'xrange',taumi,taumi+ntjpsi*delt
0910           write(ifhi,'(3a)')'yrange',' 0 ',' auto'
0911           write(ifhi,'(a)')       'columnweight 4 column c4 = ( 0 ) '
0912           write(ifhi,'(a)')       'array 4'
0913           do j=1,ntjpsi
0914             tau=taumi+(j-0.5)*delt
0915             rat=0
0916             if(jjjtot(nb,j).gt.0)rat=nnucl(nb,j)/float(jjjtot(nb,j))
0917             write(ifhi,'(4e12.4)')tau,rat,0.,float(jjjtot(nb,j))
0918           enddo
0919           write(ifhi,'(a)')       '  endarray'
0920           write(ifhi,'(a)')       'closehisto plot 0'
0921         enddo
0922       endif
0923
0924 c     nr of close nucleons vs t
0925 c     -------------------------
0926
0927       if(jpsinu.eq.1)then
0928         write(ifhi,'(a)')       'zone 3 4 1'
0929         do nb=mxbim,1,-1
0930           bim=(nb-0.5)*bimmax/mxbim
0931           write(ifhi,'(a)')       'openhisto'
0932           write(ifhi,'(a)')       'htyp lin xmod lin ymod lin'
0933           write(ifhi,'(a,f5.2,a)')'text .1 .90 "b= ',bim,' fm"'
0934           write(ifhi,'(a)')       'text 0 0 "xaxis time t (fm)"'
0935           write(ifhi,'(a)')       'text 0 0 "yaxis Nclose(b,t) / J"'
0936           write(ifhi,'(a,2e11.3)')'xrange',taumi,taumi+ntjpsi*delt
0937           write(ifhi,'(3a)')'yrange',' 0 ',' auto '
0938           write(ifhi,'(a)')       'columnweight 4 column c4 = ( 0 ) '
0939           write(ifhi,'(a)')       'array 4'
0940           do j=1,ntjpsi
0941             tau=taumi+(j-0.5)*delt
0942             rat=0
0943             if(jjjtot(nb,j).ne.0)rat=nclose(nb,j,1)/float(jjjtot(nb,j))
0944             write(ifhi,'(4e12.4)')tau,rat,0.,float(jjjtot(nb,j))
0945           enddo
0946           write(ifhi,'(a)')       '  endarray'
0947           write(ifhi,'(a)')       'closehisto plot 0'
0948         enddo
0949       endif
0950
0951 c     number of droplets
0952 c     ------------------
0953
0954       if(jpsidr.eq.1)then
0955       write(ifhi,'(a)')       'zone 3 4 1'
0956       do nb=mxbim,1,-1
0957       bim=(nb-0.5)*bimmax/mxbim
0958       write(ifhi,'(a)')       'openhisto'
0959       write(ifhi,'(a)')       'htyp lfu xmod lin ymod lin'
0960       write(ifhi,'(a,f5.2,a)')'text .1 .90 "b= ',bim,' fm"'
0961       write(ifhi,'(a)')       'text 0 0 "xaxis time t (fm)"'
0962       write(ifhi,'(a)') 'text 0 0 "yaxis D(b,t) / J"'
0963       write(ifhi,'(a,2e11.3)')'xrange',taumi,taumi+ntjpsi*delt
0964       write(ifhi,'(3a)')'yrange',' 0 ',' auto '
0965       write(ifhi,'(a)')       'columnweight 4 column c4 = ( 0 ) '
0966       write(ifhi,'(a)')       'array 4'
0967       do j=1,ntjpsi
0968       tau=taumi+(j-0.5)*delt
0969       rat=0
0970       if(jjjtot(nb,j).ne.0)rat=ndrop(nb,j)/float(jjjtot(nb,j))
0971       write(ifhi,'(4e12.4)')tau,rat,0.,float(jjjtot(nb,j))
0972       enddo
0973       write(ifhi,'(a)')       '  endarray'
0974       write(ifhi,'(a)')       'closehisto plot 0'
0975       enddo
0976       endif
0977
0978 c     number of droplets
0979 c     ------------------
0980
0981       if(jpsidr.eq.1)then
0982         write(ifhi,'(a)')       'zone 3 4 1'
0983         do nb=mxbim,1,-1
0984           bim=(nb-0.5)*bimmax/mxbim
0985           write(ifhi,'(a)')       'openhisto'
0986           write(ifhi,'(a)')       'htyp lfu xmod lin ymod lin'
0987           write(ifhi,'(a,f5.2,a)')'text .1 .90 "b= ',bim,' fm"'
0988           write(ifhi,'(a)')       'text 0 0 "xaxis time t (fm)"'
0989           write(ifhi,'(a)') 'text 0 0 "yaxis mass*D(b,t) / J"'
0990           write(ifhi,'(a,2e11.3)')'xrange',taumi,taumi+ntjpsi*delt
0991           write(ifhi,'(3a)')'yrange',' 0 ',' auto '
0992           write(ifhi,'(a)')       'columnweight 4 column c4 = ( 0 ) '
0993           write(ifhi,'(a)')       'array 4'
0994           do j=1,ntjpsi
0995             tau=taumi+(j-0.5)*delt
0996             rat=0
0997             if(jjjtot(nb,j).ne.0)rat=ndrop0(nb,j)/float(jjjtot(nb,j))
0998             write(ifhi,'(4e12.4)')tau,rat,0.,float(jjjtot(nb,j))
0999           enddo
1000           write(ifhi,'(a)')       '  endarray'
1001           write(ifhi,'(a)')       'closehisto plot 0'
1002         enddo
1003       endif
1004
1005 c$$$c     assymetry and mass of droplets
1006 c$$$c     ------------------
1007 c$$$
1008 c$$$      if(jpsidr.eq.1)then
1009 c$$$         do nb=11,1,-5
1010 c$$$            bim=(nb-0.5)*bimmax/mxbim
1011 c$$$            write(ifhi,'(a)')       'zone 3 4 1'
1012 c$$$            do nt=40,150,10
1013 c$$$               tau=taumi+(nt-0.5)*delt
1014 c$$$
1015 c$$$      write(ifhi,'(a)')       'openhisto'
1016 c$$$      write(ifhi,'(a)')       'htyp lfu xmod lin ymod lin'
1017 c$$$      write(ifhi,*) 'text .1 .90 "b= ',bim,' fm ','t=',tau,'"'
1018 c$$$      write(ifhi,'(a)')       'text 0 0 "xaxis mass "'
1019 c$$$      write(ifhi,'(a)') 'text 0 0 "yaxis assym"'
1020 c$$$      write(ifhi,'(a,2e11.3)')'xrange ',0.,40.
1021 c$$$      write(ifhi,'(a,2e11.3)')'yrange ',-5.,5.
1022 c$$$      write(ifhi,'(a)')       'columnweight 4 column c4 = ( 0 ) '
1023 c$$$      write(ifhi,'(a,i)')       'set ityp2d 5'
1024 c$$$      write(ifhi,'(a,i)')       'array2d ',mxmass
1025 c$$$
1026 c$$$      do i=1,mxassy
1027 c$$$         do j=1,mxmass
1028 c$$$            rat=0.
1029 c$$$            rat=float(ndrp2(nb,nt,j,i))
1030 c$$$     &           /zevent        ! nevent
1031 c$$$     &           /(40./20)      ! mass-bin
1032 c$$$     &           /(10./20)      ! log(assy)-bin
1033 c$$$     &           /1.            ! b-bin
1034 c$$$            write (ifhi,*) rat
1035 c$$$         enddo
1036 c$$$      enddo
1037 c$$$
1038 c$$$      write(ifhi,'(a)')       '  endarray'
1039 c$$$      write(ifhi,'(a)')       'closehisto plot2d'
1040 c$$$      enddo
1041 c$$$      enddo
1042 c$$$      endif
1043 c$$$
1044 c$$$c     volume and mass of droplets
1045 c$$$c     ------------------
1046 c$$$
1047 c$$$      if(jpsidr.eq.1)then
1048 c$$$         do nb=11,1,-5
1049 c$$$            bim=(nb-0.5)*bimmax/mxbim
1050 c$$$            write(ifhi,'(a)')       'zone 3 4 1'
1051 c$$$            do nt=40,150,10
1052 c$$$               tau=taumi+(nt-0.5)*delt
1053 c$$$
1054 c$$$      write(ifhi,'(a)')       'openhisto'
1055 c$$$      write(ifhi,'(a)')       'htyp lfu xmod lin ymod lin'
1056 c$$$      write(ifhi,*) 'text .1 .90 "b= ',bim,' fm ','t=',tau,'"'
1057 c$$$      write(ifhi,'(a)')       'text 0 0 "xaxis mass "'
1058 c$$$      write(ifhi,'(a)') 'text 0 0 "yaxis volume"'
1059 c$$$      write(ifhi,'(a,2e11.3)')'xrange ',0.,40.
1060 c$$$      write(ifhi,'(a,2e11.3)')'yrange ',-2.,3.
1061 c$$$      write(ifhi,'(a)')       'columnweight 4 column c4 = ( 0 ) '
1062 c$$$      write(ifhi,'(a,i)')       'array2d ',mxmass
1063 c$$$
1064 c$$$      do i=1,mxassy
1065 c$$$         do j=1,mxmass
1066 c$$$            rat=0.
1067 c$$$            rat=float(ndrop3(nb,nt,j,i))
1068 c$$$     &           /zevent        ! nevent
1069 c$$$     &           /(40./20)      ! mass-bin
1070 c$$$     &           /(5./20)       ! log(v)-bin
1071 c$$$     &           /1.            ! b-bin
1072 c$$$            write (ifhi,*) rat
1073 c$$$         enddo
1074 c$$$      enddo
1075 c$$$
1076 c$$$      write(ifhi,'(a)')       '  endarray'
1077 c$$$      write(ifhi,'(a)')       'closehisto plot2d'
1078 c$$$      enddo
1079 c$$$      enddo
1080 c$$$      endif
1081 c$$$
1082 c$$$c     xy of droplets
1083 c$$$c     ------------------
1084 c$$$
1085 c$$$      if(jpsidr.eq.1)then
1086 c$$$        write(ifhi,'(a)')       'zone 3 4 1'
1087 c$$$        do nb=11,1,-5
1088 c$$$          bim=(nb-0.5)*bimmax/mxbim
1089 c$$$          do nt=40,150,10
1090 c$$$            tau=taumi+(nt-0.5)*delt
1091 c$$$            write(ifhi,'(a)')       'openhisto'
1092 c$$$            write(ifhi,'(a)')       'htyp lfu xmod lin ymod lin'
1093 c$$$            write(ifhi,*) 'text .1 .90 "xy b= ',bim,' fm ','t=',tau,'"'
1094 c$$$            write(ifhi,'(a)')       'text 0 0 "xaxis x "'
1095 c$$$            write(ifhi,'(a)') 'text 0 0 "yaxis y "'
1096 c$$$            write(ifhi,'(a,2e11.3)')'xrange ',a4min,a4max
1097 c$$$            write(ifhi,'(a,2e11.3)')'yrange ',a4min,a4max
1098 c$$$            write(ifhi,'(a)')       'columnweight 4 column c4 = ( 0 ) '
1099 c$$$            write(ifhi,'(a,i)')       'array2d ',nxmdk
1100 c$$$            do i=1,nxmdk
1101 c$$$              do j=1,nxmdk
1102 c$$$                rat=0.
1103 c$$$                rat=xydens(nt,nb,i,j)
1104 c$$$     &            /zevent       ! nevent
1105 c$$$     &            /((a4max-a4min)/float(nxmdk)) ! x-bin
1106 c$$$     &            /((a4max-a4min)/float(nxmdk)) ! y-bin
1107 c$$$c...............&                 /1.      ! b-bin
1108 c$$$                write (ifhi,*) rat
1109 c$$$              enddo
1110 c$$$            enddo
1111 c$$$            write(ifhi,'(a)')       '  endarray'
1112 c$$$            write(ifhi,'(a)')       'closehisto plot2d'
1113 c$$$          enddo
1114 c$$$        enddo
1115 c$$$      endif
1116 c$$$
1117 c$$$
1118 c$$$c$$$c.....michael-verteilung.........................................
1119 c$$$c$$$      write(ifhi,'(a)')       'zone 3 4 1'
1120 c$$$c$$$      do j=1,ntjpsi
1121 c$$$c$$$        tau=taumi+(j-0.5)*delt
1122 c$$$c$$$
1123 c$$$c$$$        write(ifhi,'(a)')       'openhisto'
1124 c$$$c$$$        write(ifhi,'(a)')       'htyp lin xmod lin ymod log'
1125 c$$$c$$$        write(ifhi,'(a,f5.2,a)')'text .1 .90 "tau= ',tau,' fm"'
1126 c$$$c$$$        write(ifhi,'(a)')       'text 0 0 "xaxis mass t (fm)"'
1127 c$$$c$$$        write(ifhi,'(a)')       'text 0 0 "yaxis Nclose(b,t) / J"'
1128 c$$$c$$$        write(ifhi,'(a,2e11.3)')'xrange',a6min,a6max
1129 c$$$c$$$        write(ifhi,'(3a)')'yrange',' 0 ',' auto '
1130 c$$$c$$$        write(ifhi,'(a)')       'array 2'
1131 c$$$c$$$        do k=1,mxmass
1132 c$$$c$$$          rat=0.
1133 c$$$c$$$          amass=(k-0.5)/(mxmass)*(a6max-a6min)+a6min
1134 c$$$c$$$          rat=ami(j,k)/zevent/(a6max-a6min)*mxmass
1135 c$$$c$$$          write(ifhi,'(2e12.4)')amass,rat
1136 c$$$c$$$        enddo
1137 c$$$c$$$        write(ifhi,'(a)')       '  endarray'
1138 c$$$c$$$        write(ifhi,'(a)')       'closehisto plot 0'
1139 c$$$c$$$      enddo
1140 c$$$
1141 c$$$c     xy of strings
1142 c$$$c     -------------
1143 c$$$
1144 c$$$      if(jpsidr.eq.1)then
1145 c$$$        write(ifhi,'(a)')       'zone 3 4 1'
1146 c$$$        do nb=11,1,-1
1147 c$$$          bim=(nb-0.5)*bimmax/mxbim
1148 c$$$          write(ifhi,'(a)')       'openhisto'
1149 c$$$          write(ifhi,'(a)')       'htyp lfu xmod lin ymod lin'
1150 c$$$          write(ifhi,*) 'text .1 .90 "xy b= ',bim,' fm ','t=',tau,'"'
1151 c$$$          write(ifhi,'(a)')       'text 0 0 "xaxis x "'
1152 c$$$          write(ifhi,'(a)') 'text 0 0 "yaxis y "'
1153 c$$$          write(ifhi,'(a,2e11.3)')'xrange ',a4min,a4max
1154 c$$$          write(ifhi,'(a,2e11.3)')'yrange ',a4min,a4max
1155 c$$$          write(ifhi,'(a)')       'columnweight 4 column c4 = ( 0 ) '
1156 c$$$          write(ifhi,'(a,i)')       'array2d ',nxmdk
1157 c$$$          do i=1,nxmdk
1158 c$$$             do j=1,nxmdk
1159 c$$$                rat=0.
1160 c$$$                rat=xys(nb,i,j)
1161 c$$$     &               /zevent    ! nevent
1162 c$$$     &               /((a5max-a5min)/float(nxmdk)) ! x-bin
1163 c$$$     &               /((a5max-a5min)/float(nxmdk)) ! y-bin
1164 c$$$c...............&                 /1.      ! b-bin
1165 c$$$                write (ifhi,*) rat
1166 c$$$             enddo
1167 c$$$          enddo
1168 c$$$          write(ifhi,'(a)')       '  endarray'
1169 c$$$          write(ifhi,'(a)')       'closehisto plot2d'
1170 c$$$       enddo
1171 c$$$      endif
1172
1173
1174 c     fraction of jpsis in a droplet
1175 c     ------------------------------
1176
1177       if(jpsidr.eq.1)then
1178       write(ifhi,'(a)')       'zone 3 4 1'
1179       do nb=mxbim,1,-1
1180       bim=(nb-0.5)*bimmax/mxbim
1181       write(ifhi,'(a)')       'openhisto'
1182       write(ifhi,'(a)')       'htyp lfu xmod lin ymod lin'
1183       write(ifhi,'(a,f5.2,a)')'text .1 .90 "b= ',bim,' fm"'
1184       write(ifhi,'(a)')       'text 0 0 "xaxis time t (fm)"'
1185       write(ifhi,'(a)') 'text 0 0 "yaxis Jdrop(b,t) / J"'
1186       write(ifhi,'(a,2e11.3)')'xrange',taumi,taumi+ntjpsi*delt
1187       write(ifhi,'(3a)')'yrange',' 0 ',' auto '
1188       write(ifhi,'(a)')       'columnweight 4 column c4 = ( 0 ) '
1189       write(ifhi,'(a)')       'array 4'
1190       do j=1,ntjpsi
1191       tau=taumi+(j-0.5)*delt
1192       rat=0
1193       if(jjjtot(nb,j).ne.0)rat=jjjdro(nb,j)/float(jjjtot(nb,j))
1194       write(ifhi,'(4e12.4)')tau,rat,0.,float(jjjtot(nb,j))
1195       enddo
1196       write(ifhi,'(a)')       '  endarray'
1197       write(ifhi,'(a)')       'closehisto plot 0'
1198       enddo
1199       endif
1200
1201 c     fraction of jpsis with taud gt tauc
1202 c     -----------------------------------
1203
1204       if(jpsidr.eq.1)then
1205       write(ifhi,'(a)')       'zone 2 4 1'
1206       do ntauc=1,mxtauc
1207       tauc=ntauc*(taudmx/mxtauc)
1208       write(ifhi,'(a)')       'openhisto'
1209       write(ifhi,'(a)')       'htyp lfu xmod lin ymod lin'
1210       write(ifhi,'(a,f5.2,a)')'text .1 .90 "tauc= ',tauc,' fm/c"'
1211       write(ifhi,'(a)')       'text 0 0 "xaxis bmax-b (fm)"'
1212       write(ifhi,'(a)')
1213      *'text 0 0 "yaxis J(b, taud) / J(b)"'
1214       write(ifhi,'(a,2e11.3)')'xrange',0.,bimmax
1215       write(ifhi,'(3a)')'yrange',' 0 ',' auto '
1216       write(ifhi,'(a)')       'columnweight 4 column c4 = ( 0 ) '
1217       write(ifhi,'(a)')       'array 4'
1218       do j=mxbim,1,-1
1219       bim=(j-0.5)*bimmax/mxbim
1220       rat=0
1221       if(jjtot(j).gt.0.)rat=float(jjjtau(j,ntauc))/jjtot(j)
1222       write(ifhi,'(4e12.4)')bimmax-bim,rat,0.,float(jjtot(j))
1223       enddo
1224       write(ifhi,'(a)')       '  endarray'
1225       write(ifhi,'(a)')       'closehisto plot 0'
1226       enddo
1227       endif
1228
1229 c     droplet survival ratio
1230 c     --------------
1231
1232       if(jpsidr.eq.1)then
1233       write(ifhi,'(a)')       'zone 2 4 1'
1234       do ntauc=1,mxtauc
1235       tauc=ntauc*(taudmx/mxtauc)
1236       write(ifhi,'(a)')       'openhisto'
1237       write(ifhi,'(a)')       'htyp lfu xmod lin ymod lin'
1238       write(ifhi,'(a,f5.2,a)')'text .1 .90 "tauc= ',tauc,' fm/c"'
1239       write(ifhi,'(a)')       'text 0 0 "xaxis bmax-b (fm)"'
1240       write(ifhi,'(a)')      'text 0 0 "yaxis droplet survival ratio"'
1241       write(ifhi,'(a,2e11.3)')'xrange',0.,bimmax
1242       write(ifhi,'(3a)')'yrange',' 0 ',' auto '
1243       write(ifhi,'(a)')       'columnweight 4 column c4 = ( 0 ) '
1244       write(ifhi,'(a)')       'array 4'
1245       do j=mxbim,1,-1
1246       bim=(j-0.5)*bimmax/mxbim
1247       rat=0
1248       if(jjtot(j).gt.0.)rat=float(jjtot(j)-jjjtau(j,ntauc))/jjtot(j)
1249       write(ifhi,'(4e12.4)')bimmax-bim,rat,0.,float(jjtot(j))
1250       enddo
1251       write(ifhi,'(a)')       '  endarray'
1252       write(ifhi,'(a)')       'closehisto plot 0'
1253       enddo
1254       endif
1255
1256 c     total approx. survival ratio
1257 c     --------------
1258
1259       if(jpsidr.eq.1)then
1260       write(ifhi,'(a)')       'zone 2 4 1'
1261       do ntauc=1,mxtauc
1262       tauc=ntauc*(taudmx/mxtauc)
1263       write(ifhi,'(a)')       'openhisto'
1264       write(ifhi,'(a)')       'htyp lfu xmod lin ymod lin'
1265       write(ifhi,'(a,f5.2,a)')'text .1 .90 "tauc= ',tauc,' fm/c"'
1266       write(ifhi,'(a)')       'text 0 0 "xaxis bmax-b (fm)"'
1267       write(ifhi,'(a)')      'text 0 0 "yaxis tot. ap. survival ratio"'
1268       write(ifhi,'(a,2e11.3)')'xrange',0.,bimmax
1269       write(ifhi,'(3a)')'yrange',' 0 ',' auto '
1270       write(ifhi,'(a)')       'columnweight 4 column c4 = ( 0 ) '
1271       write(ifhi,'(a)')       'array 4'
1272       do j=mxbim,1,-1
1273       bim=(j-0.5)*bimmax/mxbim
1274       rat=0
1275       if(jjtot(j).gt.0.)rat=float(jjtot(j)-jjjtau(j,ntauc))/jjtot(j)
1276       write(ifhi,'(4e12.4)')bimmax-bim,rat,0.,float(jjtot(j))
1277       enddo
1278       write(ifhi,'(a)')       '  endarray'
1279       write(ifhi,'(a)')       'closehisto '
1280       if(maproj.eq.208.and.matarg.eq.208)then
1281       write(ifhi,'(a,a)')   'openhisto htyp lfu input jpbpb',
1282      &        ' j-nucl mult plot 0- '
1283       write(ifhi,'(a)')'openhisto htyp ldo input jpbpb j-nucl plot 0- '
1284       write(ifhi,'(a)')       'openhisto'
1285       write(ifhi,'(a)')       'set fmsc 1.0'
1286       write(ifhi,'(a)')       'column c1 = ( 11.9 - c1  )'
1287       write(ifhi,'(a)')       'column c2 = ( c2 * 0.019 )'
1288       write(ifhi,'(a)')       'input na50 ratio-b plot 0'
1289       elseif(maproj.eq.32.and.matarg.eq.32)then
1290          write(ifhi,'(a,a)')   'openhisto htyp lfu input jss',
1291      &        ' j-nucl mult plot 0- '
1292       write(ifhi,'(a)')'openhisto htyp ldo input jss j-nucl plot 0 '
1293       endif
1294       enddo
1295       endif
1296
1297 c     total survival ratio
1298 c     --------------
1299
1300       if(jpsidr.eq.1)then
1301       write(ifhi,'(a)')       'zone 2 4 1'
1302       do ntauc=1,mxtauc
1303       tauc=ntauc*(taudmx/mxtauc)
1304       write(ifhi,'(a)')       'openhisto'
1305       write(ifhi,'(a)')       'htyp lfu xmod lin ymod lin'
1306       write(ifhi,'(a,f5.2,a)')'text .1 .90 "tauc= ',tauc,' fm/c"'
1307       write(ifhi,'(a)')       'text 0 0 "xaxis bmax-b (fm)"'
1308       write(ifhi,'(a)')      'text 0 0 "yaxis total survival ratio"'
1309       write(ifhi,'(a,2e11.3)')'xrange',0.,bimmax
1310       write(ifhi,'(3a)')'yrange',' 0 ',' auto '
1311       write(ifhi,'(a)')       'columnweight 4 column c4 = ( 0 ) '
1312       write(ifhi,'(a)')       'array 4'
1313       do j=mxbim,1,-1
1314       bim=(j-0.5)*bimmax/mxbim
1315       rat=0
1316       if(jjtot(j).gt.0.)rat=float(jjtot(j)-jjjnt(j,ntauc))/jjtot(j)
1317       write(ifhi,'(4e12.4)')bimmax-bim,rat,0.,float(jjtot(j))
1318       enddo
1319       write(ifhi,'(a)')       '  endarray'
1320       write(ifhi,'(a)')       'closehisto plot 0-'
1321
1322       write(ifhi,'(a)')       'openhisto'
1323       write(ifhi,'(a)')       'htyp ldo xmod lin ymod lin'
1324       write(ifhi,'(a)')       'columnweight 4 column c4 = ( 0 ) '
1325       write(ifhi,'(a)')       'array 4'
1326       do j=mxbim,1,-1
1327       bim=(j-0.5)*bimmax/mxbim
1328       rat=0
1329       if(jjtot(j).gt.0.)rat=float(jjtot(j)-jjnuc(j))/jjtot(j)
1330       write(ifhi,'(4e12.4)')bimmax-bim,rat,0.,float(jjtot(j))
1331       enddo
1332       write(ifhi,'(a)')       '  endarray'
1333       if(maproj.eq.208.and.matarg.eq.208)then
1334       write(ifhi,'(a)')       'closehisto plot 0-'
1335       write(ifhi,'(a)')       'openhisto'
1336       write(ifhi,'(a)')       'set fmsc 1.0'
1337       write(ifhi,'(a)')       'column c1 = ( 11.9 - c1  )'
1338       write(ifhi,'(a)')       'column c2 = ( c2 * 0.019 )'
1339       write(ifhi,'(a)')       'input na50 ratio-b plot 0'
1340       else
1341       write(ifhi,'(a)')       ' closehisto plot 0'
1342       endif
1343       enddo
1344       endif
1345
1346       end
1347