EPOS/epos167/epos-jps-164.f

   1 c-----------------------------------------------------------------------
   2       subroutine jpsifo(npjpsi)
   3 c-----------------------------------------------------------------------
   4 c   forms a jpsi
   5 c-----------------------------------------------------------------------
   6       include 'epos.inc'
   7       include 'epos.incems'
   8       common/geom/rmproj,rmtarg,bmax,bkmx
   9       common/nucl3/phi,bimp
  10       parameter (ndep=129,ndet=129)
  11       common/cdep/xdep(ndep),qdep(ndep),wdep(ndep)
  12       common/cdet/xdet(ndet),qdet(ndet),wdet(ndet)
  13       parameter (nptj=129)
  14       common /cptj/xptj(nptj),qptj(nptj),wptj(nptj)
  15       parameter (mxbim=12)
  16       common/jpsi1/bimmax,kolran,delt,taumi,jpsinu,jpsidr,taudmx
  17       parameter (mxmass=20)
  18
  19       parameter (nxmdk=20)
  20       parameter (ntjpsi=150)
  21       common/jpsi7/xydens(ntjpsi,mxbim,nxmdk,nxmdk),a4min,a4max
  22       common/jpsi8/xys(mxbim,nxmdk,nxmdk),a5min,a5max
  23       common/jpsi9/ami(ntjpsi,mxmass),a6min,a6max
  24
  25       call utpri('jpsifo',ish,ishini,4)
  26       if(ish.ge.6)write(ifch,'(a)')' jpsi formation'
  27
  28 c     trigger
  29 c     -------
  30
  31       ymax=0.5
  32       ymin=-0.5
  33
  34 c     jpsi momenta
  35 c     ------------
  36       id=441
  37       call idmass(id,amass)
  38       s=amass**2
  39  2    rqptj=rangen()*qptj(nptj)
  40       pt=utinvt(nptj,xptj,qptj,rqptj)
  41       phi=2.*pi*rangen()
  42       px=pt*cos(phi)
  43       py=pt*sin(phi)
  44       lo=0
  45     1 lo=lo+1
  46       if(lo.gt.10)call utstop('jpsifo: lo > 10 &')
  47       z=0.19*sqrt(-2*alog(rangen()))*cos(2*pi*rangen()) !1-dim gauss
  48
  49
  50       if(z.gt.1.)goto 1
  51       pz=z*engy/2*ransig()
  52       e=sqrt(s+px**2+py**2+pz**2)
  53       amt=sqrt(amass**2+pt**2)
  54       y=sign(1.,pz)*alog( (e+abs(pz))/amt )
  55       if(y.lt.ymin.or.y.gt.ymax)goto 2
  56
  57 c     fill /cptl/
  58 c     -----------
  59       if(npjpsi.eq.0)then
  60         nptl=nptl+1
  61         npjpsi=nptl
  62       endif
  63       if(npjpsi.gt.mxptl)then
  64         print *,npjpsi,mxptl
  65         call utstop('jpsifo: npjpsi>mxptl&')
  66       endif
  67       istptl(npjpsi)=1
  68       idptl(npjpsi)=id
  69       pptl(1,npjpsi)=px
  70       pptl(2,npjpsi)=py
  71       pptl(3,npjpsi)=pz
  72       pptl(4,npjpsi)=e
  73       pptl(5,npjpsi)=amass
  74       kolran=1+rangen()*kolevt
  75       xorptl(1,npjpsi)=coord(1,kolran)
  76       xorptl(2,npjpsi)=coord(2,kolran)
  77       xorptl(3,npjpsi)=coord(3,kolran)
  78       xorptl(4,npjpsi)=coord(4,kolran)
  79       iorptl(npjpsi)=0
  80       jorptl(npjpsi)=0
  81       tivptl(1,npjpsi)=xorptl(4,npjpsi)
  82       tivptl(2,npjpsi)=ainfin
  83       ifrptl(1,npjpsi)=0
  84       ifrptl(2,npjpsi)=0
  85       if(ish.ge.6) then
  86         call alist("&",npjpsi,npjpsi)
  87         write (ifch,*) xorptl(1,npjpsi)
  88      $       ,xorptl(2,npjpsi),xorptl(3,npjpsi),xorptl(4,npjpsi)
  89      $       ,tivptl(1,npjpsi),tivptl(2,npjpsi)
  90         ii=iproj(kolran)
  91         jj=maproj+itarg(kolran)
  92         call alist("collision&",ii,ii)
  93         call alist("&",jj,jj)
  94       endif
  95       a4min=-15.
  96       a4max= 15.
  97       a5min=-15.
  98       a5max= 15.
  99       a6min= 2.
 100       a6max= 10.
 101       call utprix('jpsifo',ish,ishini,4)
 102       return
 103       end
 104
 105 c-----------------------------------------------------------------------
 106       function sptj(x)
 107 c-----------------------------------------------------------------------
 108 c     jpsi pt-distribution in 200 gev pp
 109 c-----------------------------------------------------------------------
 110       a=0.95
 111       c=1/0.363
 112       z=x/a
 113       sptj=1/a*c**c/utgam1(c)*z**(c-1)*exp(-c*z)
 114       return
 115       end
 116
 117 c-----------------------------------------------------------------------
 118       subroutine jpsian(ifirst)
 119 c-----------------------------------------------------------------------
 120 c     jpsi analysis.
 121 c-----------------------------------------------------------------------
 122       include 'epos.inc'
 123       include 'epos.incems'
 124       parameter (mxbim=12,ntjpsi=150,mxtauc=16)
 125       common/jpsi1/bimmax,kolran,delt,taumi,jpsinu,jpsidr,taudmx
 126       common/jpsi2/jjtot(mxbim),jjnuc(mxbim),jjjtau(mxbim,mxtauc)
 127       common/jpsi3/jjjtot(mxbim,ntjpsi),jjjdro(mxbim,ntjpsi)
 128       common/jpsi4/nnucl(mxbim,ntjpsi),nclose(mxbim,ntjpsi,3)
 129       common/jpsi5/ndrop(mxbim,ntjpsi),jjjnt(mxbim,mxtauc)
 130       parameter (mxmass=20,mxassy=20)
 131       common/jpsi6/ndrp2(mxbim,ntjpsi,mxmass,mxassy)
 132      &     ,ndrop3(mxbim,ntjpsi,mxmass,mxassy)
 133       parameter (nxmdk=20)
 134       common/jpsi7/xydens(ntjpsi,mxbim,nxmdk,nxmdk),a4min,a4max
 135       common/jpsi8/xys(mxbim,nxmdk,nxmdk),a5min,a5max
 136       common/jpsi9/ami(ntjpsi,mxmass),a6min,a6max
 137       common/jpsi10/ndrop0(mxbim,ntjpsi)
 138
 139       double precision tpro,zpro,ttar,ztar,ttaus,detap,detat
 140       common/cttaus/tpro,zpro,ttar,ztar,ttaus,detap,detat /ctain/mtain
 141       common/geom/rmproj,rmtarg,bmax,bkmx
 142       common/nucl3/phi,bimp
 143       parameter (ndep=129,ndet=129)
 144       common/cdep/xdep(ndep),qdep(ndep),wdep(ndep)
 145       common/cdet/xdet(ndet),qdet(ndet),wdet(ndet)
 146
 147       common/c9ptl/tauptl(mxptl),ss0ptl(mxptl)
 148
 149       call utpri('jpsian',ish,ishini,5)
 150
 151       detap=(ypjtl-yhaha)*etafac
 152       detat=-yhaha*etafac
 153       tpro=dcosh(detap)
 154       zpro=dsinh(detap)
 155       ttar=dcosh(detat)
 156       ztar=dsinh(detat)
 157
 158       jpsinu=1
 159       jpsidr=1
 160
 161       fac=1                     !    <-------- should be one finally
 162       taudmx=4
 163       rad=sqrt(0.62 / pi)
 164       taud=0
 165       nucia=0
 166       taumi=-2
 167       delt=0.1
 168       bimmax=amin1(rmproj+rmtarg,bmaxim)
 169       delbim=bimmax/mxbim
 170       ii=iproj(kolran)
 171       jj=maproj+itarg(kolran)
 172
 173
 174 c.....first event: delete commom blocks...............................
 175       if(ifirst.eq.1)then
 176         ifirst=0
 177         do nbim=1,mxbim
 178           jjtot(nbim)=0
 179           jjnuc(nbim)=0
 180           do ix=1,nxmdk
 181             do iy=1,nxmdk
 182               xys(nbim,ix,iy)=0.
 183             enddo
 184           enddo
 185           do nt=1,ntjpsi
 186             jjjtot(nbim,nt)=0
 187             nnucl(nbim,nt)=0
 188             jjjdro(nbim,nt)=0
 189             ndrop(nbim,nt)=0
 190             ndrop0(nbim,nt)=0
 191             do kk=1,3
 192               nclose(nbim,nt,kk)=0
 193             enddo
 194             do ix=1,nxmdk
 195               do iy=1,nxmdk
 196                 xydens(nt,nbim,ix,iy)=0.
 197               enddo
 198             enddo
 199           enddo
 200           do mm=1,mxtauc
 201             jjjtau(nbim,mm)=0
 202           enddo
 203         enddo
 204       endif
 205
 206       nbim=1+int(bimevt/delbim)
 207       if(nbim.lt.0.or.nbim.gt.mxbim) goto 5
 208       jjtot(nbim)=jjtot(nbim)+1 !events pro bin
 209
 210       do 1 i=1,nptl
 211         if(idptl(i).eq.441)j=i
 212  1    continue
 213
 214 c     if(jpsidr.eq.1)then
 215 c     write(6,'(a,i5,a,f6.2,a,f6.2,a,f6.2)')'ip=',ip
 216 c     *,'   rin= '
 217 c     *,'   mass= ',pptl(5,ip)
 218 c     enddo
 219 c     endif
 220 c     endif
 221
 222       taumax=0.
 223       ttaus=taumi-delt
 224       do 2 nt=1,ntjpsi
 225         idrin=0
 226         ttaus=ttaus+delt        !increment of time
 227         if(ish.ge.6)write(ifch,*) 'ttaus:-->',ttaus,ii,jj
 228         jpsiex=1
 229         call jtain(j,xj,yj,zj,tj,n,1)
 230         if(n.eq.1.or.n.eq.2.or.n.eq.9)jpsiex=0 !goto 2
 231         if(jpsiex.eq.1)jjjtot(nbim,nt)=jjjtot(nbim,nt)+1
 232
 233 c       nucleons
 234 c       --------
 235
 236         if(jpsinu.eq.1.and.jpsiex.eq.1)then !test jpsi-nucleon collision
 237           do 6 i=1,maproj+matarg
 238             if(i.eq.ii.or.i.eq.jj)goto 6
 239             nnucl(nbim,nt)=nnucl(nbim,nt)+1
 240             t=sngl(ttaus)
 241             x=xorptl(1,i)+(t-xorptl(4,i))*pptl(1,i)/pptl(4,i)
 242             y=xorptl(2,i)+(t-xorptl(4,i))*pptl(2,i)/pptl(4,i)
 243             z=xorptl(3,i)+(t-xorptl(4,i))*pptl(3,i)/pptl(4,i)
 244             pde=(pptl(3,i)+pptl(3,j))/(pptl(4,i)+pptl(4,j))
 245             gam2i=1-pde**2
 246             if(gam2i.eq.0.)goto 6
 247             dist=sqrt((x-xj)**2+(y-yj)**2
 248      &           +1/gam2i*(z-zj-(t-tj)*pde)**2)
 249             if(dist.le.rad)then
 250               nclose(nbim,nt,1)=nclose(nbim,nt,1)+1
 251               nucia=1
 252               if(ish.ge.6)then
 253                 write (ifch,*) "nucl dist:",dist,' dist(sig)='
 254      $               ,rad
 255                 call alist("&",i,i)
 256                 call alist("&",j,j)
 257               endif
 258             elseif(dist.le.rad+1)then
 259               nclose(nbim,nt,2)=nclose(nbim,nt,2)+1
 260             elseif(dist.le.rad+3)then
 261               nclose(nbim,nt,3)=nclose(nbim,nt,3)+1
 262             endif
 263  6        continue
 264         endif
 265
 266 c       particles
 267 c       ---------
 268         do 8 i=1,nptl
 269 c          if ( i.eq.ii.or.i.eq.jj ) goto 8
 270           call jtain(i,x,y,z,t,n,1)
 271           call jtaus(z,tz,sz)
 272           if(n.eq.1.or.n.eq.2.or.n.eq.9)goto 8
 273 c         calculate s
 274 c         -----------
 275           iad=iabs(idptl(i))
 276 c          s=(pptl(4,i)+pptl(4,j))**2-(pptl(3,i)+pptl(3,j))**2
 277 c     $         -(pptl(2,i)+pptl(2,j))**2-(pptl(1,i)+pptl(1,j))**2
 278           if ( iad.eq.120 .or. iad.eq.110 ) then !pion
 279             sig=1.              !
 280           elseif ( iad.eq.121 .or. iad.eq.111 ) then ! rho
 281             sig=1.              !
 282           elseif ( iad.eq.1120 .or. iad.eq.1220 ) then
 283             sig=3.0             ! ???? or 6 ????
 284           else
 285             goto 8
 286           endif
 287           call jtaus(zj,tzj,szj)         !????????????????? OK ?
 288           dist=sqrt((x-xj)**2+(y-yj)**2+(sz-szj)**2)
 289           if ( dist .lt. sqrt(0.1*sig/pi) ) then
 290             istptl(j)=2
 291             if(ish.ge.6)then
 292               write (ifch,*) "dist:",dist,' dist(sig)='
 293      $             ,sqrt(0.1*sig/pi),' sig=',sig
 294               call alist("&",i,i)
 295               call alist("&",j,j)
 296             endif
 297           endif
 298  8      continue
 299
 300 c     droplets
 301 c     --------
 302
 303         if(jpsidr.eq.1)then
 304           call jtaus(zj,tzj,szj)
 305           do 3 i=maproj+matarg+1,nptl
 306 c...........x-y distribution of strings..............................
 307             if(istptl(i).eq.29.and.nt.eq.1)then
 308               call jtain(i,x,y,z,t,n,1)
 309               if(x.gt.a5min.and.x.lt.a5max.and.
 310      &             y.gt.a5min.and.y.lt.a5max)then
 311                 ix=(x-a5min)/(a5max-a5min)*nxmdk + 1
 312                 iy=(y-a5min)/(a5max-a5min)*nxmdk + 1
 313                 xys(nbim,ix,iy)=xys(nbim,ix,iy)+pptl(5,i)
 314               endif
 315             endif
 316             if(istptl(i).gt.10)goto 3
 317             if(i.eq.j)goto 3
 318
 319 c...................................................................
 320  384        continue
 321             call jtain(i,x,y,z,t,n,1)
 322             if(n.eq.1.or.n.eq.2.or.n.eq.9)goto 3
 323             stop'jpsian: change!!!!        ' !call jintep(i,x,y,z,t,sz,eps,rho)
 324             if(eps.lt.aouni)goto 3 !min-dichte
 325             ndrop(nbim,nt)=ndrop(nbim,nt)+1 !droplets at time nt
 326             ndrop0(nbim,nt)=ndrop0(nbim,nt)+pptl(5,i) !mass
 327             des=0 !?????????????????????????????????
 328             o=sz+des
 329             u=sz-des
 330             r=( xxxx(i) +sngl(ttaus) ) *fac
 331 c..............assym-mass-distribution...............................
 332             assym=log(des/r)
 333             amass=pptl(5,i)
 334             a1min=-5.
 335             a1max=5.
 336             a2min=0.
 337             a2max=40.
 338             if(assym.ge.a1min.and.assym.lt.a1max
 339      &           .and.amass.ge.a2min.and.amass.lt.a2max
 340      &           ) then
 341               nassym=(assym-a1min)/(a1max-a1min)*mxassy+1
 342               namass=(amass-a2min)/(a2max-a2min)*mxmass+1
 343               ndrp2(nbim,nt,namass,nassym)=
 344      &             ndrp2(nbim,nt,namass,nassym)+1
 345             endif
 346
 347 c..............vol-mass-distribution...............................
 348             a3min=-2.
 349             a3max=3.
 350             v=log(pi*r**2.*2.*des)/log(10.)
 351             if(v.ge.a3min.and.v.lt.a3max
 352      &           .and.amass.ge.a2min.and.amass.lt.a2max
 353      &           ) then
 354               nv=(v-a3min)/(a3max-a3min)*mxassy+1
 355               namass=(amass-a2min)/(a2max-a2min)*mxmass+1
 356               ndrop3(nbim,nt,namass,nv)=
 357      &             ndrop3(nbim,nt,namass,nv)+1
 358             endif
 359 c..............x-y distribution of droplet..............................
 360             ix=(x-a4min)/(a4max-a4min)*nxmdk + 1
 361             iy=(y-a4min)/(a4max-a4min)*nxmdk + 1
 362             xydens(nt,nbim,ix,iy)=xydens(nt,nbim,ix,iy)+eps
 363
 364             if(jpsiex.eq.0)goto 3
 365
 366 c           if(mod(nt,10).eq.1)
 367 c           write(6,'(f5.2,i6,5x,3f7.2,5x,4f6.2)')sngl(ttaus),i,
 368 c           *              u,szj,o,sqrt((x-xj)**2+(y-yj)**2),r,v,eps
 369
 370             if(szj.lt.u.or.szj.gt.o)goto 3
 371             if((x-xj)**2+(y-yj)**2.gt.r**2)goto 3
 372
 373 c           write(6,'(a,f5.2,a,i5,a)')'***** t=',sngl(ttaus)
 374 c           *,' -- jpsi in droplet ',i,' *****'
 375
 376             taud=taud+delt
 377             taumax=max(taud,taumax)
 378
 379 c           write (*,*) taud,taumax
 380
 381             idrin=1
 382             jjjdro(nbim,nt)=jjjdro(nbim,nt)+1
 383  3        continue
 384
 385 c     if (idrin.ne.1)taud=max(taud-delt,0.)
 386
 387           if (idrin.ne.1)taud=0.
 388         endif
 389
 390  2    continue                  !end nt-loop
 391
 392       ijmod=2
 393       if(ijmod.eq.2)taud=taumax
 394       if(nucia.eq.1)jjnuc(nbim)=jjnuc(nbim)+1
 395       if(taud.gt.0.)then
 396         do ntaud=1,mxtauc
 397           tauc=ntaud*taudmx/mxtauc
 398           if(taud.gt.tauc)jjjtau(nbim,ntaud)=jjjtau(nbim,ntaud)+1
 399           if(nucia.eq.1.or.taud.gt.tauc)
 400      &         jjjnt(nbim,ntaud)=jjjnt(nbim,ntaud)+1
 401         enddo
 402       else
 403         do ntaud=1,mxtauc
 404           if(nucia.eq.1)
 405      &         jjjnt(nbim,ntaud)=jjjnt(nbim,ntaud)+1
 406         enddo
 407       endif
 408
 409     5 continue
 410       call utprix('jpsian',ish,ishini,5)
 411       return
 412       end
 413
 414 c-----------------------------------------------------------------------
 415       subroutine jtauan(is,im)
 416 c-----------------------------------------------------------------------
 417 c     display collision
 418 c     im = ijk
 419 c          k > 0   -->  post-script
 420 c          j > 0   -->  povray
 421 c            j = 1  -->  time and z in n-n cms
 422 c            j = 2  -->  time and z on hyberbola
 423 c          i > 0   -->  text  ( changes in alist per time step )
 424 c     cut in zeta-x (or y) plane for tau
 425 c-----------------------------------------------------------------------
 426       include 'epos.inc'
 427       double precision tpro,zpro,ttar,ztar,ttaus,detap,detat
 428       common/cttaus/tpro,zpro,ttar,ztar,ttaus,detap,detat /ctain/mtain
 429       parameter (mxbim=12,ntjpsi=150,mxtauc=16)
 430       common/jpsi1/bimmax,kolran,delt,taumi,jpsinu,jpsidr,taudmx
 431       common/jpsi2/jjtot(mxbim),jjnuc(mxbim),jjjtau(mxbim,mxtauc)
 432       common/jpsi3/jjjtot(mxbim,ntjpsi),jjjdro(mxbim,ntjpsi)
 433       common/jpsi4/nnucl(mxbim,ntjpsi),nclose(mxbim,ntjpsi,3)
 434       common/jpsi5/ndrop(mxbim,ntjpsi),jjjnt(mxbim,mxtauc)
 435       parameter (mxmass=20,mxassy=20)
 436       common/jpsi6/ndrp2(mxbim,ntjpsi,mxmass,mxassy)
 437      &     ,ndrop3(mxbim,ntjpsi,mxmass,mxassy)
 438       parameter (nxmdk=20)
 439       common/jpsi7/xydens(ntjpsi,mxbim,nxmdk,nxmdk),a4min,a4max
 440       common/jpsi8/xys(mxbim,nxmdk,nxmdk),a5min,a5max
 441       common/jpsi9/ami(ntjpsi,mxmass),a6min,a6max
 442       common/jpsi10/ndrop0(mxbim,ntjpsi)
 443       character*20 name,nnrr
 444       character*28 filename
 445       character*12 color(20)
 446       character*12 colpo(20)
 447       logical lcalc!,zet
 448       dimension isch(mxptl)
 449 c      zet=.true.
 450 c      zet=.false.
 451       xmin=-10.
 452       xmax=10.
 453       zmin=-10.
 454       zmax=10.
 455
 456 c      zevent=float(nevent*jpsi)
 457       if(mod(im,10).ne.0)then
 458         name='tau-'
 459         n=1
 460  5      l=4
 461         ll=int(log(real(n))/log(10.))+1
 462         do ii=ll,1,-1
 463           l=l+1
 464           name(l:l)=char(48+mod(int(n/10**(ii-1)),10))
 465         enddo
 466         name(l+1:l+3)='.ps'
 467         l=l+3
 468         inquire(file=name(1:l),exist=lcalc)
 469         if (lcalc)then
 470           n=n+1
 471           goto 5
 472         endif
 473         write(*,*) 'jtauan name ',name
 474         ifps=92
 475         open(unit=ifps,file=name(1:l),status='unknown')
 476         WRITE(ifps,'(a)') '%!PS-Adobe-2.0'
 477         WRITE(ifps,'(a)') '%%Title: tt2.fig'
 478         WRITE(ifps,'(a)') '%%Orientation: Portrait'
 479         WRITE(ifps,'(a)') '%%BeginSetup'
 480         WRITE(ifps,'(a)') '%%IncludeFeature: *PageSize A4'
 481         WRITE(ifps,'(a)') '%%EndSetup'
 482         WRITE(ifps,'(a)') '%%EndComments'
 483         WRITE(ifps,*) '/l {lineto} bind def'
 484         WRITE(ifps,*) '/rl {rlineto} bind def'
 485         WRITE(ifps,*) '/m {moveto} bind def'
 486         WRITE(ifps,*) '/rm {rmoveto} bind def'
 487         WRITE(ifps,*) '/s {stroke} bind def'
 488         WRITE(ifps,*) '/gr {grestore} bind def'
 489         WRITE(ifps,*) '/gs {gsave} bind def'
 490         WRITE(ifps,*) '/cp {closepath} bind def'
 491         WRITE(ifps,*) '/tr {translate} bind def'
 492         WRITE(ifps,*) '/sc {scale} bind def'
 493         WRITE(ifps,*) '/sd {setdash} bind def'
 494         WRITE(ifps,*) '/sdo {[.01 .05] 0 sd} bind def'
 495         WRITE(ifps,*) '/sdf {[1 .0] 0 sd} bind def'
 496         WRITE(ifps,*) '/n {newpath} bind def'
 497         WRITE(ifps,*) '/slw {setlinewidth } bind def'
 498         write(ifps,*) '/srgb {setrgbcolor} bind def'
 499         write(ifps,*) '/lgrey      { 0 0.95 0.95 srgb} bind def'
 500         write(ifps,*) '/black      { 0 0 0 srgb} bind def'
 501         write(ifps,*) '/red        { 1 0 0 srgb} bind def  '
 502         write(ifps,*) '/green      { 0 1 0  srgb} bind def  '
 503         write(ifps,*) '/blue       { 0 0 1  srgb} bind def  '
 504         write(ifps,*) '/yellow     { 1 0.5 0  srgb} bind def  '
 505         write(ifps,*) '/turquoise  { 0 1 1  srgb} bind def  '
 506         write(ifps,*) '/purple     { 1 0 1  srgb} bind def  '
 507 c.......write(ifps,*) '/  {   srgb} bind def  '
 508 c.......write(ifps,*) '/  {   srgb} bind def  '
 509         write(ifps,*) '/ef {eofill} bind def'
 510         WRITE(ifps,'(a)') '%%EndProlog'
 511         WRITE(ifps,*) 'gsave'
 512         WRITE(ifps,*) '/Helvetica findfont 10 scalefont setfont'
 513       endif
 514       color(9)='lgrey     '
 515       color(1)='black     '
 516       color(2)='red       '
 517       color(3)='green     '
 518       color(4)='blue      '
 519       color(7)='yellow    '
 520       color(5)='turquoise '
 521       color(6)='purple    '
 522       colpo(1)='Red  '
 523       colpo(2)='Green  '
 524       colpo(3)='Blue  '
 525       colpo(4)='Yellow  '
 526       colpo(5)='Cyan '
 527       colpo(6)='Magenta  '
 528       colpo(7)='Black  '
 529       colpo(8)='Aquamarine '
 530
 531
 532       do i=1,mxptl
 533         isch(i)=0
 534       enddo
 535
 536 c      gray0=.95
 537 c      gray1=0.
 538       iyb=0                     !????????????????????
 539       np=0
 540       nt=-10
 541       deltau=0.1
 542       taumin=-1
 543       ttaus=0
 544       do while (ttaus.lt.20.)
 545         nt=nt+1
 546        ! ttaus=dble(taumin+deltau*(factau**(1.*nt-1.)-1)/(factau-1.))
 547         ttaus=taumin+deltau*nt
 548         tau=ttaus
 549         np=np+1
 550         if(mod(im,10).ne.0)then
 551           write(ifps,'(a,i4)') '%%Page: number ',np
 552           write(ifps,'(a)') 'gsave'
 553           WRITE(ifps,*) '100 700 tr'
 554           scale=0.125
 555           WRITE(ifps,*) 1./scale,1./scale,' sc'
 556           WRITE(ifps,*) scale/2.,' slw'
 557           WRITE(ifps,*) '/Helvetica findfont ',15.*scale
 558      &         ,' scalefont setfont'
 559           write(ifps,*) color(1),' n ',zmin,xmin,' m ( tau:'
 560      $         ,tau,') show '
 561           WRITE(ifps,*) '/Helvetica findfont ',2.*scale
 562      &         ,' scalefont setfont'
 563         endif
 564
 565 *--------------------------------------------------------------------*
 566 *------ povray ------------------------------------------------------*
 567 *--------------------------------------------------------------------*
 568
 569         if (mod(im/100,10).gt.0) then
 570           write (ifch,*)   "-----",np,", tau:",ttaus,"------"
 571         endif
 572         if (mod(im/10,10).gt.0) then
 573           write (nnrr,'(i5)') np
 574           li=6-log(1.*np+0.1)/log(10.)
 575           write (*,*)   "--->"//nnrr(li:5)//"<-----",li,ttaus
 576           ifpo=55
 577           filename="tau."//nnrr(li:5)//".pov"
 578           open(unit=ifpo,file=filename,status='unknown')
 579           write (ifpo,'(a)') '#include "colors.inc";'
 580 c         write (ifpo,'(a)') '#include "shapes.inc" '
 581 c         write (ifpo,'(a)') '#include "textures.inc" '
 582           write (ifpo,'(a)') 'background {color White} '
 583           write (ifpo,'(a)') 'camera {location <0,0,-120> '
 584           write (ifpo,'(a)') '     direction <0,0,2> look_at <0,0,0>} '
 585           write (ifpo,'(a)') 'light_source{<0,300,0> color White} '
 586           write (ifpo,'(a)') 'light_source{<0,5,-90> color White} '
 587           write (ifpo,'(a)') ' '
 588           write (ifpo,'(a)') ' '
 589         endif
 590         do i=1,nptl
 591           if (istptl(i).gt.1) goto 123
 592             if((tivptl(2,i)-tivptl(1,i)).lt.1e-3
 593      $           .and.idptl(i).gt.1000000.and.iyb.eq.0)
 594      $           then
 595               write (*,*) 'tiv1=tiv2 !!!!!!!!',i
 596               tivptl(2,i)=tivptl(1,i)+100.
 597             endif
 598 c...........calculate coordinates....................................
 599             if(mod(im/10,10).eq.1) then !n-n cms frame
 600               if (istptl(i).gt.1
 601      $             .or.ttaus.lt.tivptl(1,i)
 602      $             .or.ttaus.gt.tivptl(2,i)) goto 123
 603               x=xorptl(1,i)+(ttaus-xorptl(4,i))*pptl(1,i)/pptl(4,i)
 604               y=xorptl(2,i)+(ttaus-xorptl(4,i))*pptl(2,i)/pptl(4,i)
 605               z=xorptl(3,i)+(ttaus-xorptl(4,i))*pptl(3,i)/pptl(4,i)
 606             else                !  hyperbola frame
 607               call jtain(i,x,y,z,t,n,0)
 608               call jtaus(z,tz,sz)
 609               z=sz
 610               if(n.ne.0) goto 123
 611             endif
 612 c...........plot sphere or cylinder ................................
 613             if(idptl(i).gt.700000000)
 614      $           then
 615               if(mod(im/10,10).eq.1)then
 616
 617               else
 618                 des=0 !?????????????????????????????
 619                 r=0   !(xxxx(i)+vrad*sngl(ttaus))
 620                 o=sz+des
 621                 u=sz-des
 622                 print *,ttaus,o,u,r,x,y
 623                 ic=4
 624                 if (mod(im/10,10).gt.0) then
 625                   write (ifpo,111) o,x,y,u,x,y,r,colpo(ic)
 626                 endif
 627               endif
 628 c.............text output of changes in time step ...................
 629               if (mod(im/100,10).gt.0) then
 630                 if(isch(i).eq.0)then
 631                   write (ifch,'("> ",$)')
 632                   call alist("&",i,i)
 633                   isch(i)=1
 634                 endif
 635               endif
 636 c$$$              o=sz+des
 637 c$$$              u=sz-des
 638 c$$$              r=(xxxx(i)+vrad*sngl(ttaus))
 639 c$$$              rr2=r**2-(y-yb)**2
 640 c$$$              if(rr2.gt.0.)then
 641 c$$$                write(ifps,*)
 642 c$$$     &               ,' n ',u,x-r,' m ',o,x-r,' l '
 643 c$$$     &               ,o,x+r,' l ',u,x+r,' l cp s '
 644 c$$$                write(ifps,*) ' n ',u,x-r,' m (',i,iorptl(i),') show '
 645 c$$$              endif
 646             else
 647 c.............cylinder................................................
 648               r=0.8
 649               ic=1
 650               if(abs(idptl(i)).lt.999) r=0.5
 651               if(iabs(idptl(i)).eq.1120) ic=2
 652               if(iabs(idptl(i)).eq.1220) ic=3
 653               if(iabs(idptl(i)).eq.441) ic=5
 654               if(mod(im/10,10).gt.0)then
 655                 write (ifpo,110) z,x,y,r,colpo(ic) ! sphere
 656               endif
 657 c.............text...................................................
 658               if(mod(im/100,10).gt.0)then
 659                 if(isch(i).eq.0)then
 660                   write (ifch,'("> ",$)')
 661                   call alist("&",i,i)
 662                   isch(i)=1
 663                 endif
 664               endif
 665             endif
 666             goto 124
 667  123      continue
 668 c.........text................................
 669           if(mod(im/100,10).gt.0)then
 670             if(isch(i).eq.1)then
 671               write (ifch,'("< ",$)')
 672               call alist("&",i,i)
 673               isch(i)=0
 674             endif
 675           endif
 676  124      continue
 677         enddo
 678 c.......
 679  110    format('sphere {<',G12.6,',',g12.6,',',g12.6,'>,',g12.6
 680      $       ,'pigment {color ',a,'}}')
 681  111    format('cylinder {<',
 682      $       G12.6,',',g12.6,',',g12.6,
 683      $       '>,<',
 684      $       G12.6,',',g12.6,',',g12.6,
 685      $       '>,',
 686      $       g12.6,
 687      $       'pigment {color ',a,'}}')
 688         if(mod(im/10,10).gt.0)then
 689           close(unit=ifpo)
 690         endif
 691 *-------------------------------------------------------------------*
 692 *-------   end   povray   ------------------------------------------*
 693 *-------   begin post-script ---------------------------------------*
 694 *-------------------------------------------------------------------*
 695  145    continue
 696         if(mod(im,10).eq.0) goto 159
 697         yb=-6.
 698         dy=12./12.
 699         yb=yb-dy/2
 700         do iyb=0,11
 701           yb=yb+dy
 702           WRITE(ifps,*) 'gsave'
 703           WRITE(ifps,*) (xmax-xmin)*1.1*float(int(iyb/4))
 704      &         ,-(xmax-xmin)*1.1*mod(iyb,4),' tr'
 705           write(ifps,*) ' n ',zmin,xmin,' m ',zmax,xmin,' l '
 706      &         ,zmax,xmax,' l ',zmin,xmax,' l cp s '
 707 c          ttaus=dble(tau)
 708 c.........particles in layer iyb.............
 709           do i=1,nptl
 710             if (istptl(i).gt.1) goto 10
 711             if((tivptl(2,i)-tivptl(1,i)).lt.1e-3
 712      $           .and.idptl(i).gt.1000000.and.iyb.eq.0)
 713      $           then
 714               write (*,*) 'tiv1=tiv2 !!!!!!!!',i
 715               tivptl(2,i)=tivptl(1,i)+100.
 716             endif
 717             call jtain(i,x,y,z,t,n,0)
 718             call jtaus(z,tz,sz)
 719             if(n.ne.0) goto 10
 720             if(
 721      $           (is.eq.0.or.is.eq.i.or.is.eq.iorptl(i)))then
 722 *.............
 723 *.............   is  is the particle number to observe
 724 *.............   if is=0 then all particles
 725 *.............
 726 c              .and.abs(y-yb).lt.dy/2)then
 727               des=0 !??????????????????????????????????
 728               if(iyb.eq.11.and
 729      $             .abs(tivptl(2,i)-tivptl(1,i)-100.).le.1e-4 ) then
 730                 tivptl(2,i)=tivptl(1,i)+0.01
 731               endif
 732               o=sz+des
 733               u=sz-des
 734               r=0   !(xxxx(i)+vrad*sngl(ttaus))
 735               rr2=r**2-(y-yb)**2
 736               if(rr2.gt.0.)then
 737                 r=sqrt(rr2)
 738 c                write (*,*) i,des,o,u,r,y
 739                 write(ifps,*)  color(mod(i,5)+2)
 740      &               ,' n ',u,x-r,' m ',o,x-r,' l '
 741      &               ,o,x+r,' l ',u,x+r,' l cp s '
 742                 write(ifps,*) ' n ',u,x-r,' m (',i,iorptl(i),') show '
 743               endif
 744             elseif(abs(y-yb).lt.dy/2.and.zmin.lt.sz.and.sz.lt.zmax
 745      &             .and.xmin.lt.x.and.x.lt.xmax)then
 746               r=0.8
 747               ic=1
 748               if(abs(idptl(i)).lt.999)r=0.5
 749               if(abs(idptl(i)).lt.999)ic=2
 750               if(abs(idptl(i)).eq.1120)ic=3
 751               if(abs(idptl(i)).eq.1220)ic=4
 752               if(idptl(i).eq.441) ic=7
 753
 754               io=iorptl(i)
 755               if(is.eq.0.or.io.eq.is)then
 756                 write(ifps,*) ' n ',sz,x,r,0,360,' arc ',color(ic),' s '
 757                 write(ifps,*) ' n ',sz-r,x,' m (',i,io,') show '
 758               endif
 759             endif
 760  10       enddo
 761           write(ifps,*) color(1),' n ',zmin,xmin,' m (',yb,') show '
 762           WRITE(ifps,*) 'grestore'
 763         enddo                   !yb bin
 764         write(ifps,'(a)') 'grestore'
 765         write(ifps,*) 'showpage'
 766  159    continue
 767       enddo
 768
 769
 770 c        write(ifps,*) ' n ',y0,x0,' m ',y1,x0,' l ',y1,x1,' l '
 771 c     &    ,y0,x1,' l cp s '
 772 c        write(ifps,*) ' n ',(y0+y1)/2-10.*scale,(x0+x1)/2-5.*scale
 773 c     &    ,' m (',ii,jj,') show '
 774 c        ii=nob-1
 775
 776
 777       if(mod(im,10).ne.0)then
 778         write(ifps,*) 'gr'
 779         write(ifps,'(a)') '%%Trailer'
 780         write(ifps,'(a,i4)') '%%Pages: ',np
 781         write(ifps,'(a)') '%%EOF'
 782         close(unit=ifps)
 783       endif
 784       return
 785       end
 786
 787
 788 c-----------------------------------------------------------------------
 789       subroutine jpsihi
 790 c-----------------------------------------------------------------------
 791 c     histogram
 792 c-----------------------------------------------------------------------
 793       include 'epos.inc'
 794       parameter (mxbim=12,ntjpsi=150,mxtauc=16)
 795       common/jpsi1/bimmax,kolran,delt,taumi,jpsinu,jpsidr,taudmx
 796       common/jpsi2/jjtot(mxbim),jjnuc(mxbim),jjjtau(mxbim,mxtauc)
 797       common/jpsi3/jjjtot(mxbim,ntjpsi),jjjdro(mxbim,ntjpsi)
 798       common/jpsi4/nnucl(mxbim,ntjpsi),nclose(mxbim,ntjpsi,3)
 799       common/jpsi5/ndrop(mxbim,ntjpsi),jjjnt(mxbim,mxtauc)
 800       parameter (mxmass=20,mxassy=20)
 801       common/jpsi6/ndrp2(mxbim,ntjpsi,mxmass,mxassy)
 802      &     ,ndrop3(mxbim,ntjpsi,mxmass,mxassy)
 803       parameter (nxmdk=20)
 804       common/jpsi7/xydens(ntjpsi,mxbim,nxmdk,nxmdk),a4min,a4max
 805       common/jpsi8/xys(mxbim,nxmdk,nxmdk),a5min,a5max
 806       common/jpsi9/ami(ntjpsi,mxmass),a6min,a6max
 807       common/jpsi10/ndrop0(mxbim,ntjpsi)
 808
 809       zevent=float(nevent*jpsi)
 810
 811       write(ifhi,'(a)')    'cd /users/theoric/werner/histo/newdata'
 812       write(ifhi,'(a)')       'newpage'
 813
 814 c     suppression as a function of b
 815 c     ------------------------------
 816
 817       write(ifhi,'(a)')       'zone 1 2 1 openhisto'
 818       write(ifhi,'(a)')       'htyp lin xmod lin ymod lin'
 819       write(ifhi,'(a)')       'text 0 0 "xaxis bmax-b (fm)"'
 820       write(ifhi,'(a)')       'text 0 0 "yaxis J(b) et Jnuc(b) / J"'
 821       write(ifhi,'(a,2e11.3)')'xrange',0.,bimmax
 822       write(ifhi,'(3a)')'yrange',' 0 ',' auto'
 823       write(ifhi,'(a)')       'columnweight 4 column c4 = ( 0 ) '
 824       write(ifhi,'(a)')       'array 4'
 825       do j=mxbim,1,-1
 826       bim=(j-0.5)*bimmax/mxbim
 827       write(ifhi,'(4e12.4)')bimmax-bim,jjtot(j)/zevent,0.,zevent
 828       enddo
 829       write(ifhi,'(a)')       '  endarray'
 830       write(ifhi,'(a)')       'closehisto plot 0-'
 831
 832       write(ifhi,'(a)')       'openhisto'
 833       write(ifhi,'(a)')       'htyp lin xmod lin ymod lin'
 834       write(ifhi,'(a)')       'text 0 0 "xaxis bmax-b (fm)"'
 835       write(ifhi,'(a)') 'text 0 0 " "'
 836       write(ifhi,'(a,2e11.3)')'xrange',0.,bimmax
 837       write(ifhi,'(3a)')'yrange',' 0 ',' auto'
 838       write(ifhi,'(a)')       'columnweight 4 column c4 = ( 0 ) '
 839       write(ifhi,'(a)')       'array 4'
 840       do j=mxbim,1,-1
 841       bim=(j-0.5)*bimmax/mxbim
 842       write(ifhi,'(4e12.4)')bimmax-bim,jjnuc(j)/zevent,0.,zevent
 843       enddo
 844       write(ifhi,'(a)')       '  endarray'
 845       write(ifhi,'(a)')       'closehisto plot 0'
 846
 847
 848       write(ifhi,'(a)')       'openhisto'
 849       write(ifhi,'(a)')       'htyp lin xmod lin ymod lin'
 850       write(ifhi,'(a)')       'text 0 0 "xaxis bmax-b (fm)"'
 851       write(ifhi,'(a)')       'text 0 0 "yaxis survival ratio"'
 852       write(ifhi,'(a,3e11.3)')'xrange',0.,bimmax
 853       write(ifhi,'(3a)')'yrange',' 0.2 ',' auto '
 854       write(ifhi,'(a)')       'columnweight 4 column c4 = ( 0 ) '
 855       write(ifhi,'(a)')       'array 4'
 856       do j=mxbim,1,-1
 857         bim=(j-0.5)*bimmax/mxbim
 858         rat=0
 859         if(jjtot(j).gt.0.)rat=float(jjtot(j)-jjnuc(j))/jjtot(j)
 860         write(ifhi,'(4e12.4)')bimmax-bim,rat,0.,float(jjtot(j))
 861       enddo
 862       write(ifhi,'(a)')       '  endarray'
 863       if(maproj.eq.208.and.matarg.eq.208)then
 864         write(ifhi,'(a)')       'closehisto plot 0-'
 865         write(ifhi,'(a)')       'openhisto'
 866         write(ifhi,'(a)')       'set fmsc 1.0'
 867         write(ifhi,'(a,f4.1,a)')'column c1 = ( ',bimmax,' - c1  )'
 868         write(ifhi,'(a)')       'column c2 = ( c2 * 0.02 )'
 869         write(ifhi,'(a)')       'input na50 ratio-b plot 0'
 870       else
 871         write(ifhi,'(a)')       'closehisto plot 0'
 872       endif
 873
 874 c     nr of jpsi vs t
 875 c     ---------------
 876
 877       write(ifhi,'(a)')       'zone 3 4 1'
 878       do nb=mxbim,1,-1
 879         bim=(nb-0.5)*bimmax/mxbim
 880         write(ifhi,'(a)')       'openhisto'
 881         write(ifhi,'(a)')       'htyp lin xmod lin ymod lin'
 882         write(ifhi,'(a,f5.2,a)')'text .1 .90 "b= ',bim,' fm"'
 883         write(ifhi,'(a)')       'text 0 0 "xaxis time t (fm)"'
 884         write(ifhi,'(a)')       'text 0 0 "yaxis J(b,t) / J"'
 885         write(ifhi,'(a,2e11.3)')'xrange',taumi,taumi+ntjpsi*delt
 886         write(ifhi,'(3a)')      'yrange',' 0 ',' auto'
 887         write(ifhi,'(a)')       'columnweight 4 column c4 = ( 0 ) '
 888         write(ifhi,'(a)')       'array 4'
 889         do j=1,ntjpsi
 890           tau=taumi+(j-0.5)*delt
 891           write(ifhi,'(4e12.4)')tau,float(jjjtot(nb,j))/nevent,0.,nevent
 892         enddo
 893         write(ifhi,'(a)')       '  endarray'
 894         write(ifhi,'(a)')       'closehisto plot 0'
 895       enddo
 896
 897 c     nr of nucleons vs t
 898 c     -------------------
 899
 900       if(jpsinu.eq.1)then
 901         write(ifhi,'(a)')       'zone 3 4 1'
 902         do nb=mxbim,1,-1
 903           bim=(nb-0.5)*bimmax/mxbim
 904           write(ifhi,'(a)')       'openhisto'
 905           write(ifhi,'(a)')       'htyp lin xmod lin ymod lin'
 906           write(ifhi,'(a,f5.2,a)')'text .1 .90 "b= ',bim,' fm"'
 907           write(ifhi,'(a)')       'text 0 0 "xaxis time t (fm)"'
 908           write(ifhi,'(a)')       'text 0 0 "yaxis N(b,t) / J"'
 909           write(ifhi,'(a,2e11.3)')'xrange',taumi,taumi+ntjpsi*delt
 910           write(ifhi,'(3a)')'yrange',' 0 ',' auto'
 911           write(ifhi,'(a)')       'columnweight 4 column c4 = ( 0 ) '
 912           write(ifhi,'(a)')       'array 4'
 913           do j=1,ntjpsi
 914             tau=taumi+(j-0.5)*delt
 915             rat=0
 916             if(jjjtot(nb,j).gt.0)rat=nnucl(nb,j)/float(jjjtot(nb,j))
 917             write(ifhi,'(4e12.4)')tau,rat,0.,float(jjjtot(nb,j))
 918           enddo
 919           write(ifhi,'(a)')       '  endarray'
 920           write(ifhi,'(a)')       'closehisto plot 0'
 921         enddo
 922       endif
 923
 924 c     nr of close nucleons vs t
 925 c     -------------------------
 926
 927       if(jpsinu.eq.1)then
 928         write(ifhi,'(a)')       'zone 3 4 1'
 929         do nb=mxbim,1,-1
 930           bim=(nb-0.5)*bimmax/mxbim
 931           write(ifhi,'(a)')       'openhisto'
 932           write(ifhi,'(a)')       'htyp lin xmod lin ymod lin'
 933           write(ifhi,'(a,f5.2,a)')'text .1 .90 "b= ',bim,' fm"'
 934           write(ifhi,'(a)')       'text 0 0 "xaxis time t (fm)"'
 935           write(ifhi,'(a)')       'text 0 0 "yaxis Nclose(b,t) / J"'
 936           write(ifhi,'(a,2e11.3)')'xrange',taumi,taumi+ntjpsi*delt
 937           write(ifhi,'(3a)')'yrange',' 0 ',' auto '
 938           write(ifhi,'(a)')       'columnweight 4 column c4 = ( 0 ) '
 939           write(ifhi,'(a)')       'array 4'
 940           do j=1,ntjpsi
 941             tau=taumi+(j-0.5)*delt
 942             rat=0
 943             if(jjjtot(nb,j).ne.0)rat=nclose(nb,j,1)/float(jjjtot(nb,j))
 944             write(ifhi,'(4e12.4)')tau,rat,0.,float(jjjtot(nb,j))
 945           enddo
 946           write(ifhi,'(a)')       '  endarray'
 947           write(ifhi,'(a)')       'closehisto plot 0'
 948         enddo
 949       endif
 950
 951 c     number of droplets
 952 c     ------------------
 953
 954       if(jpsidr.eq.1)then
 955       write(ifhi,'(a)')       'zone 3 4 1'
 956       do nb=mxbim,1,-1
 957       bim=(nb-0.5)*bimmax/mxbim
 958       write(ifhi,'(a)')       'openhisto'
 959       write(ifhi,'(a)')       'htyp lfu xmod lin ymod lin'
 960       write(ifhi,'(a,f5.2,a)')'text .1 .90 "b= ',bim,' fm"'
 961       write(ifhi,'(a)')       'text 0 0 "xaxis time t (fm)"'
 962       write(ifhi,'(a)') 'text 0 0 "yaxis D(b,t) / J"'
 963       write(ifhi,'(a,2e11.3)')'xrange',taumi,taumi+ntjpsi*delt
 964       write(ifhi,'(3a)')'yrange',' 0 ',' auto '
 965       write(ifhi,'(a)')       'columnweight 4 column c4 = ( 0 ) '
 966       write(ifhi,'(a)')       'array 4'
 967       do j=1,ntjpsi
 968       tau=taumi+(j-0.5)*delt
 969       rat=0
 970       if(jjjtot(nb,j).ne.0)rat=ndrop(nb,j)/float(jjjtot(nb,j))
 971       write(ifhi,'(4e12.4)')tau,rat,0.,float(jjjtot(nb,j))
 972       enddo
 973       write(ifhi,'(a)')       '  endarray'
 974       write(ifhi,'(a)')       'closehisto plot 0'
 975       enddo
 976       endif
 977
 978 c     number of droplets
 979 c     ------------------
 980
 981       if(jpsidr.eq.1)then
 982         write(ifhi,'(a)')       'zone 3 4 1'
 983         do nb=mxbim,1,-1
 984           bim=(nb-0.5)*bimmax/mxbim
 985           write(ifhi,'(a)')       'openhisto'
 986           write(ifhi,'(a)')       'htyp lfu xmod lin ymod lin'
 987           write(ifhi,'(a,f5.2,a)')'text .1 .90 "b= ',bim,' fm"'
 988           write(ifhi,'(a)')       'text 0 0 "xaxis time t (fm)"'
 989           write(ifhi,'(a)') 'text 0 0 "yaxis mass*D(b,t) / J"'
 990           write(ifhi,'(a,2e11.3)')'xrange',taumi,taumi+ntjpsi*delt
 991           write(ifhi,'(3a)')'yrange',' 0 ',' auto '
 992           write(ifhi,'(a)')       'columnweight 4 column c4 = ( 0 ) '
 993           write(ifhi,'(a)')       'array 4'
 994           do j=1,ntjpsi
 995             tau=taumi+(j-0.5)*delt
 996             rat=0
 997             if(jjjtot(nb,j).ne.0)rat=ndrop0(nb,j)/float(jjjtot(nb,j))
 998             write(ifhi,'(4e12.4)')tau,rat,0.,float(jjjtot(nb,j))
 999           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