]> git.uio.no Git - u/mrichter/AliRoot.git/blame - EPOS/epos167/epos-jps-164.f
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
Fixes needed by gfortran, it doesn't perform lazy evaluation (Piotr)
[u/mrichter/AliRoot.git] / EPOS / epos167 / epos-jps-164.f
CommitLineData
9ef1c2d9 1c-----------------------------------------------------------------------
2 subroutine jpsifo(npjpsi)
3c-----------------------------------------------------------------------
4c forms a jpsi
5c-----------------------------------------------------------------------
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
28c trigger
29c -------
30
31 ymax=0.5
32 ymin=-0.5
33
34c jpsi momenta
35c ------------
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
57c fill /cptl/
58c -----------
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
105c-----------------------------------------------------------------------
106 function sptj(x)
107c-----------------------------------------------------------------------
108c jpsi pt-distribution in 200 gev pp
109c-----------------------------------------------------------------------
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
117c-----------------------------------------------------------------------
118 subroutine jpsian(ifirst)
119c-----------------------------------------------------------------------
120c jpsi analysis.
121c-----------------------------------------------------------------------
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
174c.....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
214c if(jpsidr.eq.1)then
215c write(6,'(a,i5,a,f6.2,a,f6.2,a,f6.2)')'ip=',ip
216c *,' rin= '
217c *,' mass= ',pptl(5,ip)
218c enddo
219c endif
220c 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
233c nucleons
234c --------
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
266c particles
267c ---------
268 do 8 i=1,nptl
269c 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
273c calculate s
274c -----------
275 iad=iabs(idptl(i))
276c s=(pptl(4,i)+pptl(4,j))**2-(pptl(3,i)+pptl(3,j))**2
277c $ -(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
300c droplets
301c --------
302
303 if(jpsidr.eq.1)then
304 call jtaus(zj,tzj,szj)
305 do 3 i=maproj+matarg+1,nptl
306c...........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
319c...................................................................
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
331c..............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
347c..............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
359c..............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
366c if(mod(nt,10).eq.1)
367c write(6,'(f5.2,i6,5x,3f7.2,5x,4f6.2)')sngl(ttaus),i,
368c * 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
373c write(6,'(a,f5.2,a,i5,a)')'***** t=',sngl(ttaus)
374c *,' -- jpsi in droplet ',i,' *****'
375
376 taud=taud+delt
377 taumax=max(taud,taumax)
378
379c write (*,*) taud,taumax
380
381 idrin=1
382 jjjdro(nbim,nt)=jjjdro(nbim,nt)+1
383 3 continue
384
385c 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
414c-----------------------------------------------------------------------
415 subroutine jtauan(is,im)
416c-----------------------------------------------------------------------
417c display collision
418c im = ijk
419c k > 0 --> post-script
420c j > 0 --> povray
421c j = 1 --> time and z in n-n cms
422c j = 2 --> time and z on hyberbola
423c i > 0 --> text ( changes in alist per time step )
424c cut in zeta-x (or y) plane for tau
425c-----------------------------------------------------------------------
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)
449c zet=.true.
450c zet=.false.
451 xmin=-10.
452 xmax=10.
453 zmin=-10.
454 zmax=10.
455
456c 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 '
507c.......write(ifps,*) '/ { srgb} bind def '
508c.......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
536c gray0=.95
537c 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";'
580c write (ifpo,'(a)') '#include "shapes.inc" '
581c 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
598c...........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
612c...........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
628c.............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
636c$$$ o=sz+des
637c$$$ u=sz-des
638c$$$ r=(xxxx(i)+vrad*sngl(ttaus))
639c$$$ rr2=r**2-(y-yb)**2
640c$$$ if(rr2.gt.0.)then
641c$$$ write(ifps,*)
642c$$$ & ,' n ',u,x-r,' m ',o,x-r,' l '
643c$$$ & ,o,x+r,' l ',u,x+r,' l cp s '
644c$$$ write(ifps,*) ' n ',u,x-r,' m (',i,iorptl(i),') show '
645c$$$ endif
646 else
647c.............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
657c.............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
668c.........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
678c.......
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 '
707c ttaus=dble(tau)
708c.........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*.............
726c .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)
738c 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
770c write(ifps,*) ' n ',y0,x0,' m ',y1,x0,' l ',y1,x1,' l '
771c & ,y0,x1,' l cp s '
772c write(ifps,*) ' n ',(y0+y1)/2-10.*scale,(x0+x1)/2-5.*scale
773c & ,' m (',ii,jj,') show '
774c 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
788c-----------------------------------------------------------------------
789 subroutine jpsihi
790c-----------------------------------------------------------------------
791c histogram
792c-----------------------------------------------------------------------
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
814c suppression as a function of b
815c ------------------------------
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
874c nr of jpsi vs t
875c ---------------
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
897c nr of nucleons vs t
898c -------------------
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
924c nr of close nucleons vs t
925c -------------------------
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
951c number of droplets
952c ------------------
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
978c number of droplets
979c ------------------
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
1005c$$$c assymetry and mass of droplets
1006c$$$c ------------------
1007c$$$
1008c$$$ if(jpsidr.eq.1)then
1009c$$$ do nb=11,1,-5
1010c$$$ bim=(nb-0.5)*bimmax/mxbim
1011c$$$ write(ifhi,'(a)') 'zone 3 4 1'
1012c$$$ do nt=40,150,10
1013c$$$ tau=taumi+(nt-0.5)*delt
1014c$$$
1015c$$$ write(ifhi,'(a)') 'openhisto'
1016c$$$ write(ifhi,'(a)') 'htyp lfu xmod lin ymod lin'
1017c$$$ write(ifhi,*) 'text .1 .90 "b= ',bim,' fm ','t=',tau,'"'
1018c$$$ write(ifhi,'(a)') 'text 0 0 "xaxis mass "'
1019c$$$ write(ifhi,'(a)') 'text 0 0 "yaxis assym"'
1020c$$$ write(ifhi,'(a,2e11.3)')'xrange ',0.,40.
1021c$$$ write(ifhi,'(a,2e11.3)')'yrange ',-5.,5.
1022c$$$ write(ifhi,'(a)') 'columnweight 4 column c4 = ( 0 ) '
1023c$$$ write(ifhi,'(a,i)') 'set ityp2d 5'
1024c$$$ write(ifhi,'(a,i)') 'array2d ',mxmass
1025c$$$
1026c$$$ do i=1,mxassy
1027c$$$ do j=1,mxmass
1028c$$$ rat=0.
1029c$$$ rat=float(ndrp2(nb,nt,j,i))
1030c$$$ & /zevent ! nevent
1031c$$$ & /(40./20) ! mass-bin
1032c$$$ & /(10./20) ! log(assy)-bin
1033c$$$ & /1. ! b-bin
1034c$$$ write (ifhi,*) rat
1035c$$$ enddo
1036c$$$ enddo
1037c$$$
1038c$$$ write(ifhi,'(a)') ' endarray'
1039c$$$ write(ifhi,'(a)') 'closehisto plot2d'
1040c$$$ enddo
1041c$$$ enddo
1042c$$$ endif
1043c$$$
1044c$$$c volume and mass of droplets
1045c$$$c ------------------
1046c$$$
1047c$$$ if(jpsidr.eq.1)then
1048c$$$ do nb=11,1,-5
1049c$$$ bim=(nb-0.5)*bimmax/mxbim
1050c$$$ write(ifhi,'(a)') 'zone 3 4 1'
1051c$$$ do nt=40,150,10
1052c$$$ tau=taumi+(nt-0.5)*delt
1053c$$$
1054c$$$ write(ifhi,'(a)') 'openhisto'
1055c$$$ write(ifhi,'(a)') 'htyp lfu xmod lin ymod lin'
1056c$$$ write(ifhi,*) 'text .1 .90 "b= ',bim,' fm ','t=',tau,'"'
1057c$$$ write(ifhi,'(a)') 'text 0 0 "xaxis mass "'
1058c$$$ write(ifhi,'(a)') 'text 0 0 "yaxis volume"'
1059c$$$ write(ifhi,'(a,2e11.3)')'xrange ',0.,40.
1060c$$$ write(ifhi,'(a,2e11.3)')'yrange ',-2.,3.
1061c$$$ write(ifhi,'(a)') 'columnweight 4 column c4 = ( 0 ) '
1062c$$$ write(ifhi,'(a,i)') 'array2d ',mxmass
1063c$$$
1064c$$$ do i=1,mxassy
1065c$$$ do j=1,mxmass
1066c$$$ rat=0.
1067c$$$ rat=float(ndrop3(nb,nt,j,i))
1068c$$$ & /zevent ! nevent
1069c$$$ & /(40./20) ! mass-bin
1070c$$$ & /(5./20) ! log(v)-bin
1071c$$$ & /1. ! b-bin
1072c$$$ write (ifhi,*) rat
1073c$$$ enddo
1074c$$$ enddo
1075c$$$
1076c$$$ write(ifhi,'(a)') ' endarray'
1077c$$$ write(ifhi,'(a)') 'closehisto plot2d'
1078c$$$ enddo
1079c$$$ enddo
1080c$$$ endif
1081c$$$
1082c$$$c xy of droplets
1083c$$$c ------------------
1084c$$$
1085c$$$ if(jpsidr.eq.1)then
1086c$$$ write(ifhi,'(a)') 'zone 3 4 1'
1087c$$$ do nb=11,1,-5
1088c$$$ bim=(nb-0.5)*bimmax/mxbim
1089c$$$ do nt=40,150,10
1090c$$$ tau=taumi+(nt-0.5)*delt
1091c$$$ write(ifhi,'(a)') 'openhisto'
1092c$$$ write(ifhi,'(a)') 'htyp lfu xmod lin ymod lin'
1093c$$$ write(ifhi,*) 'text .1 .90 "xy b= ',bim,' fm ','t=',tau,'"'
1094c$$$ write(ifhi,'(a)') 'text 0 0 "xaxis x "'
1095c$$$ write(ifhi,'(a)') 'text 0 0 "yaxis y "'
1096c$$$ write(ifhi,'(a,2e11.3)')'xrange ',a4min,a4max
1097c$$$ write(ifhi,'(a,2e11.3)')'yrange ',a4min,a4max
1098c$$$ write(ifhi,'(a)') 'columnweight 4 column c4 = ( 0 ) '
1099c$$$ write(ifhi,'(a,i)') 'array2d ',nxmdk
1100c$$$ do i=1,nxmdk
1101c$$$ do j=1,nxmdk
1102c$$$ rat=0.
1103c$$$ rat=xydens(nt,nb,i,j)
1104c$$$ & /zevent ! nevent
1105c$$$ & /((a4max-a4min)/float(nxmdk)) ! x-bin
1106c$$$ & /((a4max-a4min)/float(nxmdk)) ! y-bin
1107c$$$c...............& /1. ! b-bin
1108c$$$ write (ifhi,*) rat
1109c$$$ enddo
1110c$$$ enddo
1111c$$$ write(ifhi,'(a)') ' endarray'
1112c$$$ write(ifhi,'(a)') 'closehisto plot2d'
1113c$$$ enddo
1114c$$$ enddo
1115c$$$ endif
1116c$$$
1117c$$$
1118c$$$c$$$c.....michael-verteilung.........................................
1119c$$$c$$$ write(ifhi,'(a)') 'zone 3 4 1'
1120c$$$c$$$ do j=1,ntjpsi
1121c$$$c$$$ tau=taumi+(j-0.5)*delt
1122c$$$c$$$
1123c$$$c$$$ write(ifhi,'(a)') 'openhisto'
1124c$$$c$$$ write(ifhi,'(a)') 'htyp lin xmod lin ymod log'
1125c$$$c$$$ write(ifhi,'(a,f5.2,a)')'text .1 .90 "tau= ',tau,' fm"'
1126c$$$c$$$ write(ifhi,'(a)') 'text 0 0 "xaxis mass t (fm)"'
1127c$$$c$$$ write(ifhi,'(a)') 'text 0 0 "yaxis Nclose(b,t) / J"'
1128c$$$c$$$ write(ifhi,'(a,2e11.3)')'xrange',a6min,a6max
1129c$$$c$$$ write(ifhi,'(3a)')'yrange',' 0 ',' auto '
1130c$$$c$$$ write(ifhi,'(a)') 'array 2'
1131c$$$c$$$ do k=1,mxmass
1132c$$$c$$$ rat=0.
1133c$$$c$$$ amass=(k-0.5)/(mxmass)*(a6max-a6min)+a6min
1134c$$$c$$$ rat=ami(j,k)/zevent/(a6max-a6min)*mxmass
1135c$$$c$$$ write(ifhi,'(2e12.4)')amass,rat
1136c$$$c$$$ enddo
1137c$$$c$$$ write(ifhi,'(a)') ' endarray'
1138c$$$c$$$ write(ifhi,'(a)') 'closehisto plot 0'
1139c$$$c$$$ enddo
1140c$$$
1141c$$$c xy of strings
1142c$$$c -------------
1143c$$$
1144c$$$ if(jpsidr.eq.1)then
1145c$$$ write(ifhi,'(a)') 'zone 3 4 1'
1146c$$$ do nb=11,1,-1
1147c$$$ bim=(nb-0.5)*bimmax/mxbim
1148c$$$ write(ifhi,'(a)') 'openhisto'
1149c$$$ write(ifhi,'(a)') 'htyp lfu xmod lin ymod lin'
1150c$$$ write(ifhi,*) 'text .1 .90 "xy b= ',bim,' fm ','t=',tau,'"'
1151c$$$ write(ifhi,'(a)') 'text 0 0 "xaxis x "'
1152c$$$ write(ifhi,'(a)') 'text 0 0 "yaxis y "'
1153c$$$ write(ifhi,'(a,2e11.3)')'xrange ',a4min,a4max
1154c$$$ write(ifhi,'(a,2e11.3)')'yrange ',a4min,a4max
1155c$$$ write(ifhi,'(a)') 'columnweight 4 column c4 = ( 0 ) '
1156c$$$ write(ifhi,'(a,i)') 'array2d ',nxmdk
1157c$$$ do i=1,nxmdk
1158c$$$ do j=1,nxmdk
1159c$$$ rat=0.
1160c$$$ rat=xys(nb,i,j)
1161c$$$ & /zevent ! nevent
1162c$$$ & /((a5max-a5min)/float(nxmdk)) ! x-bin
1163c$$$ & /((a5max-a5min)/float(nxmdk)) ! y-bin
1164c$$$c...............& /1. ! b-bin
1165c$$$ write (ifhi,*) rat
1166c$$$ enddo
1167c$$$ enddo
1168c$$$ write(ifhi,'(a)') ' endarray'
1169c$$$ write(ifhi,'(a)') 'closehisto plot2d'
1170c$$$ enddo
1171c$$$ endif
1172
1173
1174c fraction of jpsis in a droplet
1175c ------------------------------
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
1201c fraction of jpsis with taud gt tauc
1202c -----------------------------------
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
1229c droplet survival ratio
1230c --------------
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
1256c total approx. survival ratio
1257c --------------
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
1297c total survival ratio
1298c --------------
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