Update master to aliroot
[u/mrichter/AliRoot.git] / EPOS / epos167 / epos-dro-168.f
CommitLineData
9ef1c2d9 1c----------------------------------------------------------------------
2 subroutine amicro
3c----------------------------------------------------------------------
4c microcanonical decay of cluster specified via keu...ket, tecm, volu
5c----------------------------------------------------------------------
6 include 'epos.inc'
7 common/drop9/mkeu,mked,mkes,mkec,mkeb,mket
8 $ ,iwkeu,iwked,iwkes,iwkec,iwkeb,iwket
9 common/cxyzt/xptl(mxptl),yptl(mxptl),zptl(mxptl),tptl(mxptl)
10 *,optl(mxptl),uptl(mxptl),sptl(mxptl),rptl(mxptl,3)
11 double precision tpro,zpro,ttar,ztar,ttaus,detap,detat
60ed90b3 12 common /cttaus/ tpro,zpro,ttar,ztar,ttaus,detap,detat
9ef1c2d9 13 integer jc(nflav,2),ic(2)
14
15 call utpri('amicro',ish,ishini,4)
16
17 nevt=0
18 nptl=0
19
20 ttaus=1
21ctp060829 taus=sngl(ttaus)
22 ypjtl=6
23 yhaha=3
24 etapro=(ypjtl-yhaha)*etafac
25 etatar=-yhaha*etafac
26 detap=dble(etapro)
27 detat=dble(etatar)
28 tpro=dcosh(detap)
29 zpro=dsinh(detap)
30 ttar=dcosh(detat)
31 ztar=dsinh(detat)
32
33 if(ish.ge.5)write(ifch,'(a/6i4)')
34 *' keu ked kes kec keb ket:',keu,ked,kes,kec,keb,ket
35 do i=1,6
36 jc(i,1)=0
37 jc(i,2)=0
38 enddo
39 10 if(mkeu.ne.0.or.iwkeu.ne.0)then
40 12 x=-log(rangen())*iwkeu+mkeu
41 if(rangen().lt.0.5) x=real(mkeu)-(x-real(mkeu))
42 if( exp(-((x-mkeu)/iwkeu)**2) / exp(-abs(x-mkeu)/iwkeu)/2
43 $ .lt. rangen() ) goto 12
44 keu=nint(x)
45 endif
46 if(mked.ne.0.or.iwked.ne.0)then
47 13 x=-log(rangen())*iwked+mked
48 if(rangen().lt.0.5) x=real(mked)-(x-real(mked))
49 if( exp(-((x-mked)/iwked)**2) / exp(-abs(x-mked)/iwked)/2
50 $ .lt. rangen() ) goto 13
51 ked=nint(x)
52 endif
53 if(mkes.ne.0.or.iwkes.ne.0)then
54 15 x=-log(rangen())*iwkes+mkes
55 if(rangen().lt.0.5) x=real(mkes)-(x-real(mkes))
56 if( exp(-((x-mkes)/iwkes)**2) / exp(-abs(x-mkes)/iwkes)/2
57 $ .lt. rangen() ) goto 15
58 kes=nint(x)
59 endif
60ed90b3 60 if(real((keu+ked+kes)/3).ne.real(keu+kes+ked)/3.) goto 10
9ef1c2d9 61 if(keu.ge.0)jc(1,1)=keu
62 if(ked.ge.0)jc(2,1)=ked
63 if(kes.ge.0)jc(3,1)=kes
64 if(kec.ge.0)jc(4,1)=kec
65 if(keb.ge.0)jc(5,1)=keb
66 if(ket.ge.0)jc(6,1)=ket
67 if(keu.lt.0)jc(1,2)=-keu
68 if(ked.lt.0)jc(2,2)=-ked
69 if(kes.lt.0)jc(3,2)=-kes
70 if(kec.lt.0)jc(4,2)=-kec
71 if(keb.lt.0)jc(5,2)=-keb
72 if(ket.lt.0)jc(6,2)=-ket
73 if(ish.ge.5)write(ifch,'(a/6i4/6i4)')' jc:',jc
74 idr=0
75 do nf=1,nflav
76 do ij=1,2
77 if(jc(nf,ij).ge.10)idr=7*10**8
78 enddo
79 enddo
80 if(idr/10**8.ne.7)then
81 call idenco(jc,ic,ireten)
82c if(ireten.eq.1)call utstop('ametro: idenco ret code = 1&')
83 if(ic(1).eq.0.and.ic(2).eq.0)then
84 ic(1)=100000
85 ic(2)=100000
86 endif
87
88 idr=8*10**8+ic(1)*100+ic(2)/100
89 if(ish.ge.5)write(ifch,'(a,i9)')' id:',idr
90 dez=0.5e-4
60ed90b3 91 else
9ef1c2d9 92 idr=idr
93 * +mod(jc(1,1)+jc(2,1)+jc(3,1)+jc(4,1),10**4)*10**4
94 * +mod(jc(1,2)+jc(2,2)+jc(3,2)+jc(4,2),10**4)
95 call idtrbi(jc,ibptl(1,1),ibptl(2,1),ibptl(3,1),ibptl(4,1))
96 dez=0.5e-4
97 endif
98
99 pptl(1,1)=0
100 pptl(2,1)=0
101 pptl(3,1)=0
102 pptl(4,1)=tecm
103 pptl(5,1)=tecm
104
105 g=1
106
107 nptl=nptl+1
108 idptl(nptl)=idr
109 pptl(1,nptl)=0
110 pptl(2,nptl)=0
111 pptl(3,nptl)=0
112 pptl(4,nptl)=tecm*g
113 pptl(5,nptl)=tecm*g
114 radptl(nptl)=(3*volu/pi/4.)**0.33333
115 istptl(nptl)=10
116 tivptl(2,1)=1.
60ed90b3 117
9ef1c2d9 118 nptlb=nptl
119 call hnbaaa(nptl,iret)
120 if(iret.eq.1)stop'STOP: sr amicro: hnbaaa return code = 1. '
121 istptl(nptl)=istptl(nptl)+1
122 ifrptl(1,nptl)=nptlb+1
123 ifrptl(2,nptl)=nptl
124 x=xorptl(1,nptl)
125 y=xorptl(2,nptl)
126 z=xorptl(3,nptl)
127 t=xorptl(4,nptl)
128 do n=nptlb+1,nptl
129 iorptl(n)=nptlb
130 jorptl(n)=0
131 istptl(n)=0
132 ifrptl(1,n)=0
133 ifrptl(2,n)=0
60ed90b3 134 xorptl(1,n)=x
135 xorptl(2,n)=y
9ef1c2d9 136 xorptl(3,n)=z
137 xorptl(4,n)=t
138 tivptl(1,n)=t
139 call idtau(idptl(n),pptl(4,n),pptl(5,n),taugm)
140 r=rangen()
141 tivptl(2,n)=t+taugm*(-alog(r))
142 ityptl(n)=60
143 radptl(n)=0.
144 itsptl(n)=0
145 rinptl(n)=-9999
146 enddo
147
148 if(ioceau.eq.1)call xhnbte(0)
149 if(iociau.eq.1)call xhnbti(0)
150
151 call utprix('amicro',ish,ishini,4)
152 return
153 end
154
155c-----------------------------------------------------------------------
156 subroutine hgcaaa
157c-----------------------------------------------------------------------
158c hadronic resonance gas in grand canonical treatment
60ed90b3 159c returns T, chemical potentials and hadronic yield
160c (hadron chemical potentials as combinations of quark chemical potentials)
9ef1c2d9 161c
162c input:
163c iostat: 1: Boltzmann approximation, 0: quantum statistics /metr3/
164c tecm: droplet energy /confg/
165c volu: droplet volume /confg/
166c keu ked kes kec keb ket: net flavor number /drop5/
167c
168c output:
169c tem : temperature [GeV] /cgchg/
170c chem(1:nflav): quark chem. pot. [GeV] /cflav/
171c chemgc(1:nspecs): hadron chem. pot. [GeV] /cgchg/
172c ptlngc(1:nspecs): hadron number /cgchg/
173c rmsngc(1:nspecs): standard deviation of hadron number /cgchg/
174c
175c exact treatment (iostat=0):
176c for massive hadrons : first in Boltzmann approximation with analytical
177c expressions for particle and energy densities,
178c then by using quantum statistics in integral form,
179c extracting mu and T using numerical integration
180c and an iterative procedure to solve for mu, T
181c for massless hadrons : using analytic expressions for massles particles
182c and employing the same algorithm as for massive
183c-----------------------------------------------------------------------
184 include 'epos.inc'
60ed90b3 185 parameter (mspecs=56)
9ef1c2d9 186 common/cspecs/nspecs,ispecs(mspecs),aspecs(mspecs),gspecs(mspecs)
187 common/cgchg/rmsngc(mspecs),ptlngc(mspecs),chemgc(mspecs),tem
188 common/cbol/rmsbol(mspecs),ptlbol(mspecs),chebol(mspecs),tembol
189 common/cflavs/nflavs,kef(nflav),chem(nflav)
190 common/ciakt/gen,iafs,ians,genm
191 common/cnrit/nrit
192 gen=10.0**(-epsgc)
193 genm=gen/10.
194
195 isho=ish
196 if(ishsub/100.eq.51)ish=mod(ishsub,100)
197
198 iug=(1+iospec)/2*2-1
199
200
201c initialization
202c --------------
203 kef(1)=keu
204 kef(2)=ked
205 kef(3)=kes
206 kef(4)=kec
207 kef(5)=keb
208 kef(6)=ket
209
210 if(iug.eq.1)nflavs=1
211 if(iug.eq.3)nflavs=2
212 if(iug.eq.5)nflavs=2
213 if(iug.eq.7)nflavs=3
214 if(iug.eq.9)nflavs=3
215 if(iug.eq.11)nflavs=3
216 tem=0.0
217 do i=1,nflavs
218 chem(i)=0.0
219 enddo
220 call hgchac(0)
221 do i=1,nspecs
222 ptlngc(i)=0.0
223 rmsngc(i)=0.0
224 enddo
225 nrit=0
226
227 if(ish.ge.5)then
228 write(ifch,*)('-',l=1,10)
229 *,' entry sr hgcaaa ',('-',l=1,30)
230 write(ifch,'(1x,a,2x,3i3)')
231 *'>>> grand canonical hadron gas for droplet with u d s content:'
232 *,keu,ked,kes
233 write(ifch,'(1x,a,2x,f7.3,2x,a,2x,f7.3)')
234 *'mass [GeV]:',tecm,'volume [fm^3]:',volu
235 endif
236
237 if(iug.eq.1.and.keu.ne.0.and.ish.ge.5)then
238 write(ifch,*)'inversion impossible !!!'
239 write(ifch,*)'keu=0 required for this option'
240 write(ifch,*)'T = mu(i) = 0 returned'
241 if(ish.ge.5)write(ifch,*)('-',i=1,30)
242 *,' exit sr hgcaaa ',('-',i=1,10)
243 return
244 endif
245 if(iug.eq.3.and.(keu+ked).ne.0.and.ish.ge.5)then
246 write(ifch,*)'inversion impossible !!!'
247 write(ifch,*)'keu+ked=0 required for this option'
248 write(ifch,*)'T = mu(i) = 0 returned'
249 if(ish.ge.5)write(ifch,*)('-',i=1,30)
250 *,' exit sr hgcaaa ',('-',i=1,10)
251 return
252 endif
253 kf=keu+ked+kes+kec+keb+ket
254 kf=abs(kf)
255 if(kf.ne.0)then
256 if(mod(kf,3).ne.0.and.ish.ge.5)then
257 write(ifch,*)'inversion impossible !!!'
258 write(ifch,*)'sum must be multiple of three'
259 write(ifch,*)'T = mu(i) = 0 returned'
260 if(ish.ge.5)write(ifch,*)('-',i=1,30)
261 *,' exit sr hgcaaa ',('-',i=1,10)
262 return
263 endif
264 endif
265
266
267c initial T (m=0, baryon free)
268c -------------------------------
269 gfac=0.0
270
271 if(iostat.eq.0.and.iospec.eq.iug)then
272 do i=1,nspecs
273 igsp=int(gspecs(i))
274 if(mod(igsp,2).eq.0)then
275 gfac=gfac+7.*gspecs(i)/8.
276 else
277 gfac=gfac+gspecs(i)
278 endif
279 enddo
60ed90b3 280 if(iabs(ispecs(nspecs)).lt.10)gfac=gfac+16.
9ef1c2d9 281 tem=(tecm/volu*hquer**3*30./pi**2/gfac)**.25
282 else
283 do i=1,nspecs
284 gfac=gfac+gspecs(i)
285 enddo
286 if(iabs(ispecs(nspecs)).lt.10)gfac=gfac+16.
287 tem=(tecm/volu*hquer**3*pi**2/gfac/3.)**.25
288 tem=2.*tem
289 endif
290
291 if(ish.ge.5)write(ifch,1)'initial T :',tem
2921 format(1x,a,3x,f9.6)
293
294 if(ish.ge.5)write(ifch,*)'iospec: ',iospec
295
296 if(ish.ge.5.and.iospec.ne.iug)then
297 write(ifch,*)'inversion in Boltzmann approx. :'
60ed90b3 298 elseif(ish.ge.5.and.iospec.eq.iug)then
9ef1c2d9 299 write(ifch,*)'inversion for massless hadrons :'
300 endif
301
60ed90b3 302 if(ish.ge.5)then
9ef1c2d9 303 if(nflavs.eq.1)write(ifch,'(3x,a,8x,a)')
304 *'T:','chemu:'
305 if(nflavs.eq.2)write(ifch,'(3x,a,8x,a,5x,a)')
306 *'T:','chemu:','chemd:'
307 if(nflavs.eq.3)write(ifch,'(3x,a,8x,a,5x,a,5x,a)')
308 *'T:','chemu:','chemd:','chems:'
309 endif
310
311 k=1
31210 continue
313 if(ish.ge.9.and.mod(k,10).eq.0)
314 *write(ifch,*)'hgc iteration:',k
315 if(ish.ge.9)call hgccch(1)
316
317c search for temperature (chem=const)
318c -----------------------------------
319 idt=0
320 temo=tem
321
322 if(iospec.eq.iug)then
323
324c massless particles
325c ------------------
326 if(iostat.eq.0)then
327 if(ish.ge.9)
328 *write(ifch,*)'iteration (massless):',k
329 call hgctm0
330 elseif(iostat.eq.1)then
331 if(ish.ge.9)
332 *write(ifch,*)'iteration (Boltzmann, massless):',k
333 call hgctbo(ibna)
334 if(ibna.eq.1)then
335 tem=temo
336 goto20
337 endif
338 endif
339
340 else
341
342c Boltzmann approxiamtion (massive particles)
343c -------------------------------------------
344 if(ish.ge.9)
345 *write(ifch,*)'iteration (Boltzmann, massive):',k
346 call hgctbo(ibna)
347 if(ibna.eq.1)then
348 tem=temo
349 goto20
350 endif
351
352 endif
353
354 if(tem.le.1.e-6.and.ish.ge.5)then
355 write(ifch,*)'inversion imposssible'
356 write(ifch,*)'T:',tem
357 if(ioinco.ge.1)call hnbmin(keu,ked,kes,kec)
358 if(ish.ge.5)write(ifch,*)('-',i=1,30)
359 *,' exit sr hgcaaa ',('-',i=1,10)
360 ish=isho
361 return
362 endif
363
364 dt=abs(temo-tem)
365 if(dt.le.gen*temo.or.dt.le.genm)idt=1
366
367c search for chemical potentials (tem=const)
368c ------------------------------------------
369 idch=0
370 ibna=0
371
372 do iafs=1,nflavs
373 chemo=chem(iafs)
374
375 if(iospec.eq.iug)then
376
377c massless particles
378c ------------------
379 if(iostat.eq.0)then
380 call hgccm0
381 elseif(iostat.eq.1)then
382 call hgccbo(ibna)
383 endif
384
385 else
386
387c Boltzmann approxiamtion (massive particles)
388c -------------------------------------------
389 call hgccbo(ibna)
390
391 endif
392
393 dch=abs(chemo-chem(iafs))
394 if(ish.ge.9)write(ifch,*)'dch:',dch
395 if(dch.le.abs(gen*chemo).or.dch.le.genm)idch=idch+1
396 if(ibna.eq.1)then
397 chem(iafs)=chemo
398 call hgchac(0)
399 goto20
400 endif
401
402 enddo
403
404
405c new hadron chem. potentials
406c ---------------------------
407 call hgchac(0)
408
409
410 if(ish.ge.5.and.nflavs.eq.1)
411 *write(ifch,'(1x,f8.6,2x,f9.6)')
412 *tem,chem(1)
413 if(ish.ge.5.and.nflavs.eq.2)
414 *write(ifch,'(1x,f8.6,2x,f9.6,2x,f9.6)')
415 *tem,chem(1),chem(2)
416 if(ish.ge.5.and.nflavs.eq.3)
417 *write(ifch,'(1x,f8.6,2x,f9.6,2x,f9.6,2x,f9.6)')
418 *tem,chem(1),chem(2),chem(3)
419 if(idch.eq.nflavs.and.idt.eq.1)goto20
420
421
422 k=k+1
423
424 if(k.gt.300)then
425 if(ish.ge.5)
426 *write(ifch,*)'failure in approximate solution'
427 goto20
428 endif
429
430 goto10
431
43220 continue
433 if(ish.ge.9)call hgccch(0)
434 if(ish.ge.5)write(ifch,'(1x,a,1x,f9.6)')' T :',tem
435 do i=1,nflavs
436 if(i.eq.1.and.ish.ge.5)
437 *write(ifch,'(1x,a,1x,f9.6)')'chemu:',chem(1)
438 if(i.eq.2.and.ish.ge.5)
439 *write(ifch,'(1x,a,1x,f9.6)')'chemd:',chem(2)
440 if(i.eq.3.and.ish.ge.5)
441 *write(ifch,'(1x,a,1x,f9.6)')'chems:',chem(3)
442 enddo
60ed90b3 443
9ef1c2d9 444
445c checking results
446c ----------------
447 if(ish.ge.5)call hgcchb
448
449c particle yield
450c --------------
451 call hgcpyi(1)
452
453c checking flavor conservation
454c ----------------------------
455 if(ish.ge.5)call hgccfc
456
457 if(iug.eq.iospec.and.iostat.eq.0)then
458 if(ish.ge.5)write(ifch,*)
459 *'approximation and exact treatment equal'
460 if(ish.ge.5)write(ifch,*)('-',i=1,30)
461 *,' exit sr hgcaaa ',('-',i=1,10)
462 ish=isho
463 return
464 endif
465
466c continue or return approximate values
467c -------------------------------------
468 do i=1,nspecs
469 rmsbol(i)=rmsngc(i)
470 ptlbol(i)=ptlngc(i)
471 chebol(i)=chemgc(i)
472 enddo
473 tembol=tem
474 if(iostat.eq.1)then
475 if(ish.ge.5)write(ifch,*)('-',i=1,30)
476 *,' exit sr hgcaaa ',('-',i=1,10)
477 ish=isho
478 return
479 endif
60ed90b3 480
9ef1c2d9 481
482c quantum statistics
483c ------------------
484 if(ish.ge.5)write(ifch,*)'quantum statistics:'
60ed90b3 485 if(ish.ge.5.and.nflavs.eq.1)write(ifch,'(3x,a,8x,a)')
9ef1c2d9 486 *'T:','chemu:'
60ed90b3 487 if(ish.ge.5.and.nflavs.eq.2)write(ifch,'(3x,a,8x,a,6x,a)')
9ef1c2d9 488 *'T:','chemu:','chemd:'
60ed90b3 489 if(ish.ge.5.and.nflavs.eq.3)write(ifch,'(3x,a,8x,a,6x,a,6x,a)')
9ef1c2d9 490 *'T:','chemu:','chemd:','chems:'
491 k=1
492
49330 continue
494 if(ish.ge.9.and.mod(k,10).eq.0)
495 *write(ifch,*)'hgc iteration:',k
496
497c new temperature
498c ---------------
499 idt=0
500 temo=tem
501 call hgctex
502 if(ish.ge.5.and.nflavs.eq.1)
503 *write(ifch,'(1x,f10.8,2x,f10.7)')
504 *tem,chem(1)
505 if(ish.ge.5.and.nflavs.eq.2)
506 *write(ifch,'(1x,f10.8,2x,f10.7,2x,f10.7)')
507 *tem,chem(1),chem(2)
508 if(ish.ge.5.and.nflavs.eq.3)
509 *write(ifch,'(1x,f10.8,2x,f10.7,2x,f10.7,2x,f10.7)')
510 *tem,chem(1),chem(2),chem(3)
511
512 if(tem.le.1.e-6.and.ish.ge.5)then
513 write(ifch,*)'inversion imposssible'
514 write(ifch,*)'T:',tem
515 call hnbmin(keu,ked,kes,kec)
516 if(ish.ge.5)write(ifch,*)('-',i=1,30)
517 *,' exit sr hgcaaa ',('-',i=1,10)
518 ish=isho
519 return
520 endif
521
522 dt=abs(temo-tem)
523 if(dt.le.gen*temo.or.dt.le.genm)idt=1
524 if(ish.ge.9)write(ifch,*)'dtem:',dt
525
526c new quark chem. potentials
527c --------------------------
528 idch=0
529 do iafs=1,nflavs
530 chemo=chem(iafs)
531 call hgccex
532 dch=abs(chemo-chem(iafs))
533 if(ish.ge.9)write(ifch,*)'dche:',dch
534 if(dch.le.abs(gen*chemo).or.dch.le.genm)idch=idch+1
535 enddo
536
537c new hadron chem. potentials
538c ---------------------------
539 call hgchac(0)
540
541 if(idch.eq.nflavs.and.idt.eq.1)then
54250 continue
543 if(ish.ge.5)write(ifch,*)'results:'
544 if(ish.ge.5)write(ifch,51)' T :',tem
545 if(nflavs.ge.1.and.ish.ge.5)write(ifch,51)'chemu:',chem(1)
546 if(nflavs.ge.2.and.ish.ge.5)write(ifch,51)'chemd:',chem(2)
547 if(nflavs.ge.3.and.ish.ge.5)write(ifch,51)'chems:',chem(3)
54851 format(1x,a,3x,f9.6)
549
550c checking results
551c ----------------
60ed90b3 552 if(ish.ge.5)call hgcchh(i)
9ef1c2d9 553
554c particle yield
555c --------------
556 call hgcpyi(0)
557
558c checking flavor conservation
559c ----------------------------
560 call hgccfc
561
562 if(ish.ge.5)write(ifch,*)('-',i=1,30)
563 *,' exit sr hgcaaa ',('-',i=1,10)
564 ish=isho
60ed90b3 565 return
9ef1c2d9 566 endif
567
568 if(k.gt.300)then
569 if(ish.ge.5)
570 *write(ifch,*)'failure in exact solution'
571 if(ish.ge.5)write(ifch,*)'results:'
572 if(ish.ge.5)write(ifch,51)' T :',tem
573 if(nflavs.ge.1.and.ish.ge.5)write(ifch,51)'chemu:',chem(1)
574 if(nflavs.ge.2.and.ish.ge.5)write(ifch,51)'chemd:',chem(2)
575 if(nflavs.ge.3.and.ish.ge.5)write(ifch,51)'chems:',chem(3)
576
577c particle yield
578c --------------
579 call hgcpyi(0)
580
581 if(ish.ge.5)write(ifch,*)('-',i=1,30)
582 *,' exit sr hgcaaa ',('-',i=1,10)
583 ish=isho
584 return
585
586 endif
587
588 k=k+1
589 goto30
590
591 end
592
593
594c---------------------------------------------------------------------
595 function hgcbi0(x)
596c---------------------------------------------------------------------
597 DOUBLE PRECISION p1,p2,p3,p4,p5,p6,p7,q1,q2,q3,q4,q5,q6,q7,q8,q9,y
598 SAVE p1,p2,p3,p4,p5,p6,p7,q1,q2,q3,q4,q5,q6,q7,q8,q9
599 DATA p1,p2,p3,p4,p5,p6,p7/1.0d0,3.5156229d0,3.0899424d0,
600 *1.2067492d0,0.2659732d0,0.360768d-1,0.45813d-2/
601 DATA q1,q2,q3,q4,q5,q6,q7,q8,q9/0.39894228d0,0.1328592d-1,
602 *0.225319d-2,-0.157565d-2,0.916281d-2,-0.2057706d-1,0.2635537d-1,
603 *-0.1647633d-1,0.392377d-2/
604 if (abs(x).lt.3.75) then
605 y=dble((x/3.75)**2)
606 hgcbi0=sngl(p1+y*(p2+y*(p3+y*(p4+y*(p5+y*(p6+y*p7))))))
607 else
608 ax=abs(x)
609 y=dble(3.75/ax)
610 hgcbi0=(exp(ax)/sqrt(ax))*(q1+y*(q2+y*(q3+y*(q4+y*(q5+y*(q6+y*
611 *(q7+y*(q8+y*q9))))))))
612 endif
613 return
614 end
615
616
617c------------------------------------------------------------------------
618 function hgcbi1(x)
619c------------------------------------------------------------------------
620 DOUBLE PRECISION p1,p2,p3,p4,p5,p6,p7,q1,q2,q3,q4,q5,q6,q7,q8,q9,y
621 SAVE p1,p2,p3,p4,p5,p6,p7,q1,q2,q3,q4,q5,q6,q7,q8,q9
622 DATA p1,p2,p3,p4,p5,p6,p7/0.5d0,0.87890594d0,0.51498869d0,
623 *0.15084934d0,0.2658733d-1,0.301532d-2,0.32411d-3/
624 DATA q1,q2,q3,q4,q5,q6,q7,q8,q9/0.39894228d0,-0.3988024d-1,
625 *-0.362018d-2,0.163801d-2,-0.1031555d-1,0.2282967d-1,-0.2895312d-1,
626 *0.1787654d-1,-0.420059d-2/
627 if (abs(x).lt.3.75) then
628 y=dble((x/3.75)**2)
629 hgcbi1=x*(p1+y*(p2+y*(p3+y*(p4+y*(p5+y*(p6+y*p7))))))
630 else
631 ax=abs(x)
632 y=dble(3.75/ax)
633 hgcbi1=(exp(ax)/sqrt(ax))*(q1+y*(q2+y*(q3+y*(q4+y*(q5+y*(q6+y*
634 *(q7+y*(q8+y*q9))))))))
635 if(x.lt.0.)hgcbi1=-hgcbi1
636 endif
637 return
638 END
639
640
641c---------------------------------------------------------------------
642 function hgcbk0(x)
643c------------------------------------------------------------------------
644 DOUBLE PRECISION p1,p2,p3,p4,p5,p6,p7,q1,q2,q3,q4,q5,q6,q7,y
645 SAVE p1,p2,p3,p4,p5,p6,p7,q1,q2,q3,q4,q5,q6,q7
646 DATA p1,p2,p3,p4,p5,p6,p7/-0.57721566d0,0.42278420d0,0.23069756d0,
647 *0.3488590d-1,0.262698d-2,0.10750d-3,0.74d-5/
648 DATA q1,q2,q3,q4,q5,q6,q7/1.25331414d0,-0.7832358d-1,0.2189568d-1,
649 *-0.1062446d-1,0.587872d-2,-0.251540d-2,0.53208d-3/
650 if (x.le.2.0) then
651 y=dble(x*x/4.0)
652 hgcbk0=(-log(x/2.0)*hgcbi0(x))+(p1+y*(p2+y*(p3+y*(p4+y*(p5+y*
653 *(p6+y*p7))))))
654 else
655 y=dble(2.0/x)
656 hgcbk0=(exp(-x)/sqrt(x))*(q1+y*(q2+y*(q3+y*(q4+y*(q5+y*(q6+y*
657 *q7))))))
658 endif
659 return
660 END
661
662
663c---------------------------------------------------------------
664 function hgcbk1(x)
665c--------------------------------------------------------------------
666 DOUBLE PRECISION p1,p2,p3,p4,p5,p6,p7,q1,q2,q3,q4,q5,q6,q7,y
667 SAVE p1,p2,p3,p4,p5,p6,p7,q1,q2,q3,q4,q5,q6,q7
668 DATA p1,p2,p3,p4,p5,p6,p7/1.0d0,0.15443144d0,-0.67278579d0,
669 *-0.18156897d0,-0.1919402d-1,-0.110404d-2,-0.4686d-4/
670 DATA q1,q2,q3,q4,q5,q6,q7/1.25331414d0,0.23498619d0,-0.3655620d-1,
671 *0.1504268d-1,-0.780353d-2,0.325614d-2,-0.68245d-3/
672 if (x.le.2.0) then
673 y=dble(x*x/4.0)
674 hgcbk1=(log(x/2.0)*hgcbi1(x))+(1.0/x)*(p1+y*(p2+y*(p3+y*(p4+y*
675 *(p5+y*(p6+y*p7))))))
676 else
677 y=dble(2.0/x)
678 hgcbk1=(exp(-x)/sqrt(x))*(q1+y*(q2+y*(q3+y*(q4+y*(q5+y*(q6+y*
679 *q7))))))
680 endif
681 return
682 END
683
684
685c-------------------------------------------------------------------
686 function hgcbk(n,x)
60ed90b3 687c------------------------------------------------------------------
9ef1c2d9 688 tox=2.0/x
689 bkm=hgcbk0(x)
690 bk=hgcbk1(x)
691 do 11 j=1,n-1
692 bkp=bkm+j*tox*bk
693 bkm=bk
694 bk=bkp
69511 continue
696 hgcbk=bk
697 return
698 END
699
700
60ed90b3 701c----------------------------------------------------------------
9ef1c2d9 702 subroutine hgccbo(iba)
703c----------------------------------------------------------------
704c returns new chem(iafs) for boltzmann statistics
705c input:
706c tem
707c kef/volu
708c output:
709c chem(iafs)
710c-----------------------------------------------------------------------
711 common/cnsta/pi,pii,hquer,prom,piom,ainfin
712 common/drop6/tecm,volu
60ed90b3 713 parameter (mspecs=56)
9ef1c2d9 714 common/cspecs/nspecs,ispecs(mspecs),aspecs(mspecs),gspecs(mspecs)
715 common/cgchg/rmsngc(mspecs),ptlngc(mspecs),chemgc(mspecs),tem
716 parameter(nflav=6)
717 common/cflavs/nflavs,kef(nflav),chem(nflav)
718 common/cflac/ifok(nflav,mspecs),ifoa(nflav)
719 common/ciakt/gen,iafs,ians,genm
720 external hgcbk
721 k=1
722 iba=0
723 c1=-0.5
724 c2=0.5
725 goto11
726
727c new chemical potential
728c ----------------------
72910 chem(iafs)=c1+0.5*(c2-c1)
73011 continue
731 fd=0.0
732 call hgchac(0)
60ed90b3 733
9ef1c2d9 734 do i=1,nspecs
735
736 if(ifok(iafs,i).ne.0)then
737 if((chemgc(i)/tem).gt.70.)then
738 hpd=1.e30
739 else
740 hpd=exp(chemgc(i)/tem)
741 endif
742 if(aspecs(i).ne.0.)then
743 fk2=hgcbk(2,aspecs(i)/tem)
744 hpd=hpd*gspecs(i)*aspecs(i)**2*tem*fk2
745 */2./pi**2/hquer**3
746 else
747 hpd=hpd*gspecs(i)*tem**3/pi**2/hquer**3
748 endif
749 hfd=ifok(iafs,i)*hpd
750 fd=fd+hfd
751 endif
752
753 enddo
754
755 dfd=abs(fd-(kef(iafs)/volu))
756 if(dfd.le.abs(gen*(kef(iafs)/volu)).or.dfd.le.genm)return
757c if(abs(fd).ge.100.)then
758c iba=1
759c return
760c endif
761
762
763 if(fd.gt.(kef(iafs)/volu))then
764 c2=chem(iafs)
765 else
766 c1=chem(iafs)
767 endif
768
769 k=k+1
770 if(k.gt.300)return
771
772 goto10
60ed90b3 773
9ef1c2d9 774 end
775
776
777c----------------------------------------------------------------------
778 subroutine hgccch(iii)
779c----------------------------------------------------------------------
780c checks convergence of iterative algorithm
781c plots iteration values for T and mu_i
782c----------------------------------------------------------------------
783 include 'epos.inc'
60ed90b3 784 parameter (mspecs=56)
9ef1c2d9 785 common/cgchg/rmsngc(mspecs),ptlngc(mspecs),chemgc(mspecs),tem
786 common/cflavs/nflavs,kef(nflav),chem(nflav)
787 parameter (nbin=500)
788 common/cdatc/data(nbin),datb(nbin),datc(nbin),datd(nbin)
789 *,date(nbin),datf(nbin),datg(nbin),dath(nbin),dati(nbin)
790 common/cnrit/nrit
791 character cen*4,cvol*4,cu*3,cd*3,cs*3
792
793 if(iii.gt.0)then
794
795 nrit=nrit+1
796 data(nrit)=nrit
797 datb(nrit)=tem
798 datc(nrit)=chem(1)
799 datd(nrit)=chem(2)
800 date(nrit)=chem(3)
801
802 elseif(iii.eq.0)then
803
804 nrit=nrit+1
805 data(nrit)=nrit
806 datb(nrit)=tem
807 datc(nrit)=chem(1)
808 datd(nrit)=chem(2)
809 date(nrit)=chem(3)
810 do i=1,nrit
811 datf(i)=datb(nrit)
812 datg(i)=datc(nrit)
813 dath(i)=datd(nrit)
814 dati(i)=date(nrit)
815 enddo
816
817 x1=data(1)
818 x2=data(nrit)
819 write(cen,'(f4.1)')tecm
820 write(cvol,'(f4.1)')volu
821 write(cu,'(i3)')keu
822 write(cd,'(i3)')ked
823 write(cs,'(i3)')kes
60ed90b3 824
9ef1c2d9 825
826 write(ifhi,'(a)') 'newpage zone 1 4 1 openhisto'
827 write(ifhi,'(a)') 'htyp lin xmod lin ymod lin'
828 write(ifhi,'(a)') 'text 0 0 "xaxis Iteration"'
829 write(ifhi,'(a)') 'text 0 0 "yaxis T (GeV)"'
830 write(ifhi,'(a)') 'text 0.15 0.9 "E= '//cen//'"'
831 write(ifhi,'(a)') 'text 0.4 0.9 "V= '//cvol//'"'
832 write(ifhi,'(a,2e11.3)')'xrange',x1,x2
833 write(ifhi,'(3a)')'yrange',' auto',' auto'
834 write(ifhi,'(a)') 'array 2'
835 do j=1,nrit
836 write(ifhi,'(2e12.4)')data(j),datb(j)
837 enddo
838 write(ifhi,'(a)') ' endarray'
839 write(ifhi,'(a)') 'closehisto plot 0-'
840
841 write(ifhi,'(a)') 'openhisto'
842 write(ifhi,'(a)') 'htyp lin xmod lin ymod lin'
843 write(ifhi,'(a,2e11.3)')'xrange',x1,x2
844 write(ifhi,'(3a)')'yrange',' auto',' auto'
845 write(ifhi,'(a)') 'array 2'
846 do j=1,nrit
847 write(ifhi,'(2e12.4)')data(j),datf(j)
848 enddo
849 write(ifhi,'(a)') ' endarray'
850 write(ifhi,'(a)') 'closehisto plot 0'
60ed90b3 851
9ef1c2d9 852 write(ifhi,'(a)') 'openhisto'
853 write(ifhi,'(a)') 'htyp lin xmod lin ymod lin'
854 write(ifhi,'(a)') 'text 0 0 "xaxis Iteration"'
855 write(ifhi,'(a)') 'text 0 0 "yaxis [m]^1! (GeV)"'
856 write(ifhi,'(a)') 'text 0.15 0.9 "Q^1!= '//cu//'"'
857 write(ifhi,'(a,2e11.3)')'xrange',x1,x2
858 write(ifhi,'(3a)')'yrange',' auto',' auto'
859 write(ifhi,'(a)') 'array 2'
860 do j=1,nrit
861 write(ifhi,'(2e12.4)')data(j),datc(j)
862 enddo
863 write(ifhi,'(a)') ' endarray'
864 write(ifhi,'(a)') 'closehisto plot 0-'
865
866 write(ifhi,'(a)') 'openhisto'
867 write(ifhi,'(a)') 'htyp lin xmod lin ymod lin'
868 write(ifhi,'(a,2e11.3)')'xrange',x1,x2
869 write(ifhi,'(3a)')'yrange',' auto',' auto'
870 write(ifhi,'(a)') 'array 2'
871 do j=1,nrit
872 write(ifhi,'(2e12.4)')data(j),datg(j)
873 enddo
874 write(ifhi,'(a)') ' endarray'
875 write(ifhi,'(a)') 'closehisto plot 0'
876
877 write(ifhi,'(a)') 'openhisto'
878 write(ifhi,'(a)') 'htyp lin xmod lin ymod lin'
879 write(ifhi,'(a)') 'text 0 0 "xaxis Iteration"'
880 write(ifhi,'(a)') 'text 0 0 "yaxis [m]^2! (GeV)"'
881 write(ifhi,'(a)') 'text 0.15 0.9 "Q^2!= '//cd//'"'
882 write(ifhi,'(a,2e11.3)')'xrange',x1,x2
883 write(ifhi,'(3a)')'yrange',' auto',' auto'
884 write(ifhi,'(a)') 'array 2'
885 do j=1,nrit
886 write(ifhi,'(2e12.4)')data(j),datd(j)
887 enddo
888 write(ifhi,'(a)') ' endarray'
889 write(ifhi,'(a)') 'closehisto plot 0-'
890
891 write(ifhi,'(a)') 'openhisto'
892 write(ifhi,'(a)') 'htyp lin xmod lin ymod lin'
893 write(ifhi,'(a,2e11.3)')'xrange',x1,x2
894 write(ifhi,'(3a)')'yrange',' auto',' auto'
895 write(ifhi,'(a)') 'array 2'
896 do j=1,nrit
897 write(ifhi,'(2e12.4)')data(j),dath(j)
898 enddo
899 write(ifhi,'(a)') ' endarray'
900 write(ifhi,'(a)') 'closehisto plot 0'
901
902 write(ifhi,'(a)') 'openhisto'
903 write(ifhi,'(a)') 'htyp lin xmod lin ymod lin'
904 write(ifhi,'(a)') 'text 0 0 "xaxis Iteration"'
905 write(ifhi,'(a)') 'text 0 0 "yaxis [m]^3! (GeV)"'
906 write(ifhi,'(a)') 'text 0.15 0.9 "Q^3!= '//cs//'"'
907 write(ifhi,'(a,2e11.3)')'xrange',x1,x2
908 write(ifhi,'(3a)')'yrange',' auto',' auto'
909 write(ifhi,'(a)') 'array 2'
910 do j=1,nrit
911 write(ifhi,'(2e12.4)')data(j),date(j)
912 enddo
913 write(ifhi,'(a)') ' endarray'
914 write(ifhi,'(a)') 'closehisto plot 0-'
915
916 write(ifhi,'(a)') 'openhisto'
917 write(ifhi,'(a)') 'htyp lin xmod lin ymod lin'
918 write(ifhi,'(a,2e11.3)')'xrange',x1,x2
919 write(ifhi,'(3a)')'yrange',' auto',' auto'
920 write(ifhi,'(a)') 'array 2'
921 do j=1,nrit
922 write(ifhi,'(2e12.4)')data(j),dati(j)
923 enddo
924 write(ifhi,'(a)') ' endarray'
925 write(ifhi,'(a)') 'closehisto plot 0'
926
927 endif
60ed90b3 928
9ef1c2d9 929 return
60ed90b3 930
9ef1c2d9 931 end
932
933c-----------------------------------------------------------------------
934 subroutine hgccex
935c-----------------------------------------------------------------------
936c returns new chem(iafs) for massive quantum statistics
937c input:
938c tem
939c kef/volu
940c output:
941c chem(iafs)
942c-----------------------------------------------------------------------
943 include 'epos.inc'
60ed90b3 944 parameter (mspecs=56)
9ef1c2d9 945 common/cspecs/nspecs,ispecs(mspecs),aspecs(mspecs),gspecs(mspecs)
946 common/cflavs/nflavs,kef(nflav),chem(nflav)
947 common/cflac/ifok(nflav,mspecs),ifoa(nflav)
948 common/ciakt/gen,iafs,ians,genm
949 external hgcfhn
950
951 k=1
952
953 c1=-0.5
954 c2=0.5
955 goto11
956
957c new chemical potential
958c ----------------------
95910 chem(iafs)=c1+0.5*(c2-c1)
96011 continue
961
962 fd=0.0
963 do ians=1,nspecs
964 if(ifok(iafs,ians).ne.0)then
965
966 call hgchac(0)
60ed90b3 967 call hgclim(a,b)
9ef1c2d9 968 if(b.eq.0.0)then
969 hpd=0.0
970 else
971 call uttraq(hgcfhn,a,b,hpd)
972 endif
973 hpd=hpd*gspecs(ians)/2./pi**2/hquer**3
974 fd=fd+hpd*ifok(iafs,ians)
975
976 endif
977 enddo
978
979 dfd=abs(fd-(kef(iafs)/volu))
980 if(dfd.le.abs(gen*(kef(iafs)/volu)).or.dfd.le.genm)return
981
982 if(fd.gt.(kef(iafs)/volu))then
983 c2=chem(iafs)
984 else
985 c1=chem(iafs)
986 endif
987
988 k=k+1
989 if(k.gt.300)then
990 if(ish.ge.5)
991 *write(ifch,*)'failure at cex at iafs:',iafs
992 return
993 endif
994
995 goto10
996
997 end
998
999
1000c------------------------------------------------------------------
1001 subroutine hgccfc
1002c------------------------------------------------------------------
1003c checks flavor conservation in particle yield
1004c------------------------------------------------------------------
1005 include 'epos.inc'
60ed90b3 1006 parameter (mspecs=56)
9ef1c2d9 1007 common/cspecs/nspecs,ispecs(mspecs),aspecs(mspecs),gspecs(mspecs)
1008 common/cgchg/rmsngc(mspecs),ptlngc(mspecs),chemgc(mspecs),tem
1009 common/cflavs/nflavs,kef(nflav),chem(nflav)
1010 common/cflac/ifok(nflav,mspecs),ifoa(nflav)
1011
1012 if(ish.ge.5)write(ifch,*)'checking flavor conservation'
1013 do i=1,nflavs
1014 ckef=0.0
1015 do ii=1,nspecs
1016 ckef=ckef+ifok(i,ii)*ptlngc(ii)
1017 enddo
1018 dkef=abs(ckef-kef(i))
1019 if(dkef.le.1.e-2)then
1020 if(i.eq.1.and.ish.ge.5)write(ifch,*)'u conserved'
1021 if(i.eq.2.and.ish.ge.5)write(ifch,*)'d conserved'
60ed90b3 1022 if(i.eq.3.and.ish.ge.5)write(ifch,*)'s conserved'
9ef1c2d9 1023 else
1024 if(i.eq.1.and.ish.ge.5)write(ifch,*)'u not conserved'
1025 if(i.eq.2.and.ish.ge.5)write(ifch,*)'d not conserved'
1026 if(i.eq.3.and.ish.ge.5)write(ifch,*)'s not conserved'
1027 if(ish.ge.5)write(ifch,*)'df:',dkef
1028 endif
1029 enddo
1030
60ed90b3 1031 return
9ef1c2d9 1032 end
1033
1034c----------------------------------------------------------------
1035 subroutine hgcchb
1036c----------------------------------------------------------------
1037c checks results by numerical integration
1038c----------------------------------------------------------------
1039 include 'epos.inc'
60ed90b3 1040 parameter (mspecs=56)
9ef1c2d9 1041 common/cspecs/nspecs,ispecs(mspecs),aspecs(mspecs),gspecs(mspecs)
1042 common/cgchg/rmsngc(mspecs),ptlngc(mspecs),chemgc(mspecs),tem
1043 common/cflavs/nflavs,kef(nflav),chem(nflav)
1044 common/cflac/ifok(nflav,mspecs),ifoa(nflav)
1045 common/ciakt/gen,iafs,ians,genm
1046 external hgcfbe
1047 external hgcfbn
1048 if(ish.ge.5)write(ifch,*)
1049 *'check by numer. calc. of expect. values:'
1050 iced=0
1051 ceden=0.0
1052 do ians=1,nspecs
1053 call hgclim(a,b)
1054 if(b.eq.0.0)then
1055 cedh=0.0
1056 else
1057 call uttraq(hgcfbe,a,b,cedh)
1058 endif
1059 if(ish.ge.9)write(ifch,*)'cedh:',cedh
1060 ced=cedh*gspecs(ians)/2./pi**2/hquer**3
1061 ceden=ceden+ced
1062 enddo
1063
1064 if(iabs(ispecs(nspecs)).lt.10)
1065 *ceden=ceden+(8.*pi**2*tem**4/15.+bag4rt**4)/hquer**3
1066
1067 if(ish.ge.5)write(ifch,*)'energy density :',ceden
1068 ded=abs((tecm/volu)-ceden)
1069 if((tecm/volu)*gen.ge.ded.or.ded.le.gen)iced=1
1070 icfd=0
1071
1072 do i=1,nflavs
1073 cfd=0.0
1074 do ians=1,nspecs
1075 call hgclim(a,b)
1076 if(b.eq.0.0)then
1077 hpd=0.0
1078 else
1079 call uttraq(hgcfbn,a,b,hpd)
1080 endif
1081 hfd=ifok(i,ians)*hpd*gspecs(ians)/2./pi**2/hquer**3
1082 if(ish.ge.9)write(ifch,*)'hfd:',hfd
1083 cfd=cfd+hfd
1084 enddo
1085 if(i.eq.1.and.ish.ge.5)write(ifch,5)'flavor density u :',cfd
1086 if(i.eq.2.and.ish.ge.5)write(ifch,5)'flavor density d :',cfd
1087 if(i.eq.3.and.ish.ge.5)write(ifch,5)'flavor density s :',cfd
10885 format(1x,a,1x,f12.6)
1089 dfd=abs(cfd-(kef(i)/volu))
1090 if(abs(gen*(kef(i)/volu)).ge.dfd.or.dfd.le.gen)
1091 *icfd=icfd+1
1092 enddo
1093
1094 if(iced.eq.1.and.icfd.eq.nflavs)then
1095 if(ish.ge.5)write(ifch,*)'results agree'
1096 else
1097 if(ish.ge.5)write(ifch,*)'results disagree'
1098 endif
1099
1100 return
1101 end
1102
1103c----------------------------------------------------------------
1104 subroutine hgcchh(icorr)
1105c----------------------------------------------------------------
1106c checks results by numerical integration
1107c----------------------------------------------------------------
1108 include 'epos.inc'
60ed90b3 1109 parameter (mspecs=56)
9ef1c2d9 1110 common/cspecs/nspecs,ispecs(mspecs),aspecs(mspecs),gspecs(mspecs)
1111 common/cgchg/rmsngc(mspecs),ptlngc(mspecs),chemgc(mspecs),tem
1112 common/cflavs/nflavs,kef(nflav),chem(nflav)
1113 common/cflac/ifok(nflav,mspecs),ifoa(nflav)
1114 common/ciakt/gen,iafs,ians,genm
1115 external hgcfhe
1116 external hgcfhn
1117 icorr=0
1118 if(ish.ge.5)write(ifch,*)
1119 *'check by numer. calc. of expect. values:'
1120
1121 iced=0
1122 ceden=0.0
1123 do ians=1,nspecs
1124 call hgclim(a,b)
1125 if(b.eq.0.0)then
1126 cedh=0.0
1127 else
1128 call uttraq(hgcfhe,a,b,cedh)
1129 endif
1130 if(ish.ge.9)write(ifch,*)'cedh:',cedh
1131 ced=cedh*gspecs(ians)/2./pi**2/hquer**3
1132 ceden=ceden+ced
1133 enddo
1134
1135 if(iabs(ispecs(nspecs)).lt.10)
1136 *ceden=ceden+(8.*pi**2*tem**4/15.+bag4rt**4)/hquer**3
1137
1138 if(ish.ge.5)write(ifch,*)'energy density :',ceden
1139 ded=abs((tecm/volu)-ceden)
1140 if((tecm/volu)*gen.ge.ded.or.ded.le.gen)iced=1
1141
1142 icfd=0
1143
1144 do i=1,nflavs
1145 cfd=0.0
1146 do ians=1,nspecs
1147 call hgclim(a,b)
1148 if(b.eq.0.0)then
1149 hpd=0.0
1150 else
1151 call uttraq(hgcfhn,a,b,hpd)
1152 endif
1153 hfd=ifok(i,ians)*hpd*gspecs(ians)/2./pi**2/hquer**3
1154 if(ish.ge.9)write(ifch,*)'hfd:',hfd
60ed90b3 1155 cfd=cfd+hfd
9ef1c2d9 1156 enddo
1157 if(i.eq.1.and.ish.ge.5)write(ifch,5)'flavor density u :',cfd
1158 if(i.eq.2.and.ish.ge.5)write(ifch,5)'flavor density d :',cfd
1159 if(i.eq.3.and.ish.ge.5)write(ifch,5)'flavor density s :',cfd
11605 format(1x,a,1x,f9.6)
1161 dfd=abs(cfd-(kef(i)/volu))
1162 if(abs(gen*(kef(i)/volu)).ge.dfd.or.dfd.le.gen)
1163 *icfd=icfd+1
1164 enddo
1165
1166 if(iced.eq.1.and.icfd.eq.nflavs)then
1167 if(ish.ge.5)write(ifch,*)'results agree'
1168 icorr=1
1169 else
1170 if(ish.ge.5)write(ifch,*)'results disagree'
1171 endif
1172
60ed90b3 1173 return
9ef1c2d9 1174 end
1175
1176
1177c--------------------------------------------------------------------
1178 subroutine hgccm0
1179c--------------------------------------------------------------------
1180c returns new quark chemical potentials for massless quantum statistics
1181c input:
1182c tem
1183c kef/volu
1184c output:
1185c chem
1186c---------------------------------------------------------------------
1187 include 'epos.inc'
60ed90b3 1188 parameter (mspecs=56)
9ef1c2d9 1189 common/cspecs/nspecs,ispecs(mspecs),aspecs(mspecs),gspecs(mspecs)
1190 common/cgchg/rmsngc(mspecs),ptlngc(mspecs),chemgc(mspecs),tem
1191 common/cflavs/nflavs,kef(nflav),chem(nflav)
1192 common/cflac/ifok(nflav,mspecs),ifoa(nflav)
1193 common/ciakt/gen,iafs,ians,genm
1194 external hgcfhn
1195 k=1
1196 z3=1.2020569
1197
1198 c1=-0.5
1199 c2=0.5
1200 goto11
1201
1202c new chemical potential
1203c ----------------------
120410 chem(iafs)=c1+0.5*(c2-c1)
120511 continue
1206
1207 fd=0.0
1208 call hgchac(0)
1209
1210 do i=1,nspecs
1211 if(ifok(iafs,i).ne.0)then
1212
1213 igsp=int(gspecs(i))
1214 if(mod(igsp,2).eq.0)then
1215
1216 if(ispecs(i).gt.0)then
1217 hpd=gspecs(i)*(chemgc(i)*tem**2+chemgc(i)**3/pi**2)/6./hquer**3
1218 else
1219 hpd=0.0
1220 endif
1221
1222c else
1223c if(ispecs(i).gt.0)then
1224c hpd=gspecs(i)*(chemgc(i)*tem**2/3.-chemgc(i)**3/pi**2/6.)/hquer**3
1225c else
1226c hpd=0.0
1227c endif
1228c endif
1229
1230c n=1
1231c0 xx=n*abs(chemgc(i))/tem
1232c if(xx.le.60.)then
1233c hpd=hpd+(-1.)**(n+1)/n**3/exp(xx)
1234c n=n+1
1235c goto20
1236c endif
1237c hpd=hpd*gspecs(i)*tem**3/pi**2/hquer**3
1238c if(chemgc(i).eq.abs(chemgc(i)))then
60ed90b3 1239c hpd=gspecs(i)*(chemgc(i)*tem**2+chemgc(i)**3/pi**2)/6./hquer**3
9ef1c2d9 1240c *-hpd
1241c endif
1242
1243c else
1244c hpd=3.*gspecs(i)*tem**3*z3/4./pi**2/hquer**3
1245c endif
60ed90b3 1246
9ef1c2d9 1247 else
1248
1249 hpd=gspecs(i)*tem**3*z3/pi**2/hquer**3
1250
1251 endif
1252
1253 hfd=hpd*ifok(iafs,i)
1254 fd=fd+hfd
1255
1256 endif
1257 enddo
1258
1259 dfd=abs(fd-(kef(iafs)/volu))
1260 if(dfd.le.abs(gen*(kef(iafs)/volu)).or.dfd.le.genm)return
1261
1262 if(fd.gt.(kef(iafs)/volu))then
1263 c2=chem(iafs)
1264 else
1265 c1=chem(iafs)
1266 endif
1267
1268 k=k+1
1269 if(k.gt.300)then
1270 if(ish.ge.5)
1271 *write(ifch,*)'failure at cm0 at iafs:',iafs
1272 return
1273 endif
1274 goto10
1275 end
1276
1277c-----------------------------------------------------------------------
1278 function hgcfbe(x)
1279c-----------------------------------------------------------------------
1280c integrand of energy density
1281c------------------------------------------------------------------------
60ed90b3 1282 parameter (mspecs=56)
9ef1c2d9 1283 common/cspecs/nspecs,ispecs(mspecs),aspecs(mspecs),gspecs(mspecs)
1284 common/cgchg/rmsngc(mspecs),ptlngc(mspecs),chemgc(mspecs),tem
1285 common/ciakt/gen,iafs,ians,genm
1286 eex=81.
1287 hgcfbe=0.0
1288 sq=sqrt(x**2+aspecs(ians)**2)
1289 if(tem.ne.0.0)eex=(sq-chemgc(ians))/tem
1290 if(eex.gt.60.)return
1291 if(eex.lt.-60)then
1292 hgcfbe=1.e25
1293 return
1294 endif
1295
1296 hgcfbe=sq*x**2*exp(-eex)
1297
1298 return
1299 end
1300
1301c-----------------------------------------------------------------
1302 function hgcfbf(x)
1303c-----------------------------------------------------------------
1304c integrand of mean square variance of energy
1305c----------------------------------------------------------------
60ed90b3 1306 parameter (mspecs=56)
9ef1c2d9 1307 common/cspecs/nspecs,ispecs(mspecs),aspecs(mspecs),gspecs(mspecs)
1308 common/cgchg/rmsngc(mspecs),ptlngc(mspecs),chemgc(mspecs),tem
1309 common/ciakt/gen,iafs,ians,genm
1310 eex=61
1311 hgcfbf=0.0
1312
1313 sq=sqrt(x**2+aspecs(ians)**2)
1314 if(tem.ne.0.0)eex=(sq-chemgc(ians))/tem
1315 if(eex.gt.60.)return
1316 if(eex.lt.-60)then
1317 hgcfbf=1.e25
1318 return
1319 endif
1320
1321 hgcfbf=(aspecs(ians)**2+x**2)*x**2*exp(-eex)
1322
1323 return
1324 end
1325
1326c-----------------------------------------------------------------
1327 function hgcfbn(x)
1328c-----------------------------------------------------------------
1329c integrand of hadron density
1330c-----------------------------------------------------------------
60ed90b3 1331 parameter (mspecs=56)
9ef1c2d9 1332 common/cspecs/nspecs,ispecs(mspecs),aspecs(mspecs),gspecs(mspecs)
1333 common/cgchg/rmsngc(mspecs),ptlngc(mspecs),chemgc(mspecs),tem
1334 common/ciakt/gen,iafs,ians,genm
1335 eex=81.
1336 hgcfbn=0.0
1337
1338 sq=sqrt(x**2+aspecs(ians)**2)
1339 if(tem.ne.0.0)eex=(sq-chemgc(ians))/tem
1340 if(eex.gt.80.)return
1341 if(eex.lt.-60)then
1342 hgcfbn=1.e25
1343 return
1344 endif
1345
1346 hgcfbn=x**2*exp(-eex)
1347
1348 return
1349 end
1350
1351c-----------------------------------------------------------------------
1352 function hgcfhe(x)
1353c-----------------------------------------------------------------------
1354c integrand of energy density
1355c------------------------------------------------------------------------
60ed90b3 1356 parameter (mspecs=56)
9ef1c2d9 1357 common/cspecs/nspecs,ispecs(mspecs),aspecs(mspecs),gspecs(mspecs)
1358 common/cgchg/rmsngc(mspecs),ptlngc(mspecs),chemgc(mspecs),tem
1359 common/ciakt/gen,iafs,ians,genm
1360 eex=81.
1361 hgcfhe=0.0
1362 igsp=int(gspecs(ians))
1363
1364 sq=sqrt(x**2+aspecs(ians)**2)
1365 if(tem.ne.0.0)eex=(sq-chemgc(ians))/tem
1366 if(eex.gt.80.)return
1367
1368 if(mod(igsp,2).ne.0)then
1369 d=-1.0
60ed90b3 1370 if(eex.lt.1.e-10)return
9ef1c2d9 1371 else
1372 d=1.0
60ed90b3 1373 endif
9ef1c2d9 1374
1375 hgcfhe=sq*x**2/(exp(eex)+d)
1376
60ed90b3 1377 return
9ef1c2d9 1378 end
1379
1380c-----------------------------------------------------------------
1381 function hgcfhf(x)
1382c-----------------------------------------------------------------
1383c integrand of mean square variance of energy
1384c----------------------------------------------------------------
60ed90b3 1385 parameter (mspecs=56)
9ef1c2d9 1386 common/cspecs/nspecs,ispecs(mspecs),aspecs(mspecs),gspecs(mspecs)
1387 common/cgchg/rmsngc(mspecs),ptlngc(mspecs),chemgc(mspecs),tem
1388 common/ciakt/gen,iafs,ians,genm
1389 eex=61
1390 hgcfhf=0.0
1391 igsp=int(gspecs(ians))
1392
1393 sq=sqrt(x**2+aspecs(ians)**2)
1394 if(tem.ne.0.0)eex=(sq-chemgc(ians))/tem
1395 if(eex.gt.60.)return
1396 if(eex.lt.(-60.))return
1397
1398 if(mod(igsp,2).ne.0)then
1399 d=-1.0
1400 if(eex.lt.1.0e-10.and.eex.gt.(-1.0e-10))return
1401 else
1402 d=1.0
1403 endif
1404
1405 hgcfhf=(aspecs(ians)**2+x**2)*x**2/(exp(eex)+2.0*d+exp(-eex))
1406
1407 return
1408 end
1409
1410c-----------------------------------------------------------------
1411 function hgcfhn(x)
1412c-----------------------------------------------------------------
1413c integrand of hadron density
1414c-----------------------------------------------------------------
60ed90b3 1415 parameter (mspecs=56)
9ef1c2d9 1416 common/cspecs/nspecs,ispecs(mspecs),aspecs(mspecs),gspecs(mspecs)
1417 common/cgchg/rmsngc(mspecs),ptlngc(mspecs),chemgc(mspecs),tem
1418 common/ciakt/gen,iafs,ians,genm
1419 eex=81.
1420 hgcfhn=0.0
1421 igsp=int(gspecs(ians))
1422
1423 sq=sqrt(x**2+aspecs(ians)**2)
1424 if(tem.ne.0.0)eex=(sq-chemgc(ians))/tem
1425 if(eex.gt.80.)return
1426
1427 if(mod(igsp,2).ne.0)then
1428 d=-1.0
60ed90b3 1429 if(eex.lt.1.e-10)return
9ef1c2d9 1430 else
1431 d=1.0
60ed90b3 1432 endif
9ef1c2d9 1433
1434 hgcfhn=x**2/(exp(eex)+d)
1435
1436 return
1437 end
1438
1439c-----------------------------------------------------------------
1440 function hgcfhw(x)
1441c-----------------------------------------------------------------
1442c integrand of mean square variance of hadron yield
1443c----------------------------------------------------------------
60ed90b3 1444 parameter (mspecs=56)
9ef1c2d9 1445 common/cspecs/nspecs,ispecs(mspecs),aspecs(mspecs),gspecs(mspecs)
1446 common/cgchg/rmsngc(mspecs),ptlngc(mspecs),chemgc(mspecs),tem
1447 common/ciakt/gen,iafs,ians,genm
1448 eex=61
1449 hgcfhw=0.0
1450 igsp=int(gspecs(ians))
1451
1452 sq=sqrt(x**2+aspecs(ians)**2)
1453 if(tem.ne.0.0)eex=(sq-chemgc(ians))/tem
60ed90b3 1454 if(eex.gt.60.)return
1455 if(eex.lt.(-60.))return
9ef1c2d9 1456
1457 if(mod(igsp,2).ne.0)then
1458 d=-1.0
60ed90b3 1459 if(eex.lt.1.0e-10.and.eex.gt.(-1.0e-10))return
9ef1c2d9 1460 else
1461 d=1.0
60ed90b3 1462 endif
9ef1c2d9 1463
1464 hgcfhw=x**2/(exp(eex)+2.0*d+exp(-eex))
1465
1466 return
1467 end
1468
1469
1470c-----------------------------------------------------------------
1471 subroutine hgchac(iboco)
1472c------------------------------------------------------------------
1473c returns hadronic chemical potentials as combinations of quark
1474c chemical potentials
1475c----------------------------------------------------------------------
1476 include 'epos.inc'
60ed90b3 1477 parameter (mspecs=56)
9ef1c2d9 1478 common/cspecs/nspecs,ispecs(mspecs),aspecs(mspecs),gspecs(mspecs)
1479 common/cgchg/rmsngc(mspecs),ptlngc(mspecs),chemgc(mspecs),tem
1480 common/cflavs/nflavs,kef(nflav),chem(nflav)
1481 common/cflac/ifok(nflav,mspecs),ifoa(nflav)
1482
1483 do i=1,nspecs
1484 chemgc(i)=0.0
1485 do ii=1,nflavs
1486 chemgc(i)=chemgc(i)+ifok(ii,i)*chem(ii)
1487 if(ish.ge.9)write(ifch,*)'mu_i:',chem(ii),' k_i:',ifok(ii,i)
1488 enddo
1489 if(ish.ge.9)write(ifch,*)'mu_nu:',chemgc(i)
1490 igsp=int(gspecs(i))
1491 if(mod(igsp,2).ne.0.and.chemgc(i).gt.aspecs(i).and.iboco.eq.0)
1492 *chemgc(i)=aspecs(i)
1493 enddo
1494
1495 return
1496 end
1497
1498
1499c-----------------------------------------------------------------------
1500 subroutine hgclim(a,b)
1501c----------------------------------------------------------------------
1502c returns integration limits for numerical evaluation of particle
1503c and energy densities using quantum statistics
1504c----------------------------------------------------------------------
1505 include 'epos.inc'
60ed90b3 1506 parameter (mspecs=56)
9ef1c2d9 1507 common/cspecs/nspecs,ispecs(mspecs),aspecs(mspecs),gspecs(mspecs)
1508 common/cgchg/rmsngc(mspecs),ptlngc(mspecs),chemgc(mspecs),tem
1509 common/ciakt/gen,iafs,ians,genm
1510
1511 igsp=int(gspecs(ians))
1512
1513 if(mod(igsp,2).ne.0)then
1514 a=0.001
1515 else
1516 a=0.0
1517 endif
1518
1519 b=0.0
1520 bb=(chemgc(ians)+tem*80.)**2-aspecs(ians)**2
1521 if(ish.ge.9)write(ifch,*)'bb:',bb
1522 if(bb.ge.0.0)b=sqrt(bb)
1523 if(bb.lt.0.0)then
1524 if(ish.ge.9)write(ifch,*)'failure at hgclim, bb=',bb
1525 if(ish.ge.9)write(ifch,'(1x,a,i5,a,2x,f12.6,1x,a,2x,f9.6)')
1526 *'mu(',ispecs(ians),'):',chemgc(ians),' T:',tem
1527 endif
1528 if(ish.ge.9)write(ifch,*)'ians:',ians,' a:',a,' b:',b
1529 return
1530 end
1531
1532c------------------------------------------------------------------------
1533 subroutine hgcnbi(iret)
1534c-----------------------------------------------------------------------
1535c uses hgcaaa results to generate initial hadron set, nlattc, iozero
1536c input:
1537c ptlngc(1:nspecs): particle number expectation values /cgchg/
1538c output:
1539c nump: number of hadrons /chnbin/
1540c ihadro(1:nump): hadron ids /chnbin/
1541c nlattc: lattice size /clatt/
1542c iozero: zero weight /metr1/
1543c-----------------------------------------------------------------------
1544 include 'epos.inc'
1545 parameter(maxp=500)
1546 common/chnbin/nump,ihadro(maxp)
60ed90b3 1547 parameter (mspecs=56)
9ef1c2d9 1548 common/cspecs/nspecs,ispecs(mspecs),aspecs(mspecs),gspecs(mspecs)
1549 common/cgchg/rmsngc(mspecs),ptlngc(mspecs),chemgc(mspecs),tem
1550 common/cgctot/rmstot,ptltot
1551 common/camgc/amgc,samgc,amtot
1552 common/cflavs/nflavs,kef(nflav),chem(nflav)
1553 common/cflac/ifok(nflav,mspecs),ifoa(nflav)
1554 common/clatt/nlattc,npmax
1555 common/cgcnb/nptlgc(mspecs)
1556 common/ctaue/taue
1557 common/cgck/k(nflav),kp(nflav),kps(nflav)
1558 *,idp(maxp),ida(mspecs),idb(mspecs)
1559 integer hgcndn
1560
1561 iret=0
1562 isho=ish
1563 if(ishsub/100.eq.50)ish=mod(ishsub,100)
1564
1565 if(ish.ge.7)write(ifch,*)('-',l=1,10)
1566 *,' entry sr hgcnbi ',('-',l=1,30)
1567
1568
1569 nh=nint(ptltot)
1570 iug=(1+iospec)/2*2-1
1571 if(iug.lt.9)call utstop('hgcnbi: iospec < 9&')
1572
1573c determine nlattc
1574c ----------------
1575 if(ionlat.eq.1)then
1576 s1=ptltot+2.*rmstot
1577 s2=1.3*ptltot
1578 s=max(s1,s2,6.)
1579 nlattc=nint(s)
1580 elseif(ionlat.eq.2)then
1581 s1=ptltot+3.*rmstot
1582 s2=1.5*ptltot
1583 s=max(s1,s2,6.)
1584 nlattc=nint(s)
1585 elseif(ionlat.eq.3)then
1586 s1=ptltot+4.*rmstot
1587 s2=2.*ptltot
1588 s=max(s1,s2,6.)
1589 nlattc=nint(s)
1590 elseif(ionlat.eq.0)then
1591 nlattc=8*(tecm/10)*(1/(tecm/volu))**0.2*(nspecs/3.)**0.3
1592 if(aspecs(1).lt.0.010)nlattc=nlattc*3
1593 nlattc=max(nlattc,20)
1594 endif
1595
1596 if(ish.ge.7)write(ifch,*)'nlattc:',nlattc
1597
1598c determine iozero
1599c ----------------
1600 if(iozero.eq.-1)then
1601 iozero=nspecs
1602 elseif(iozero.eq.-2)then
1603 iozero=nspecs*int(sqrt(volu/tecm))
1604 endif
1605
1606c modify iozero for testing
1607c -------------------------
1608 if(iozevt.gt.0)then
1609 iozero=(nrevt/iozevt+1)*iozinc !nrevt=event number - 1 !!
1610 write(ifch,*)'nrevt+1:',nrevt+1,' iozero:',iozero
1611 endif
1612
1613c initial hadron set
1614c ------------------
1615 ammin=2.*aspecs(1)
1616 if(tecm.lt.ammin)then
1617 write(ifch,*)'impossible to generate hadron configuration'
1618 call utstop('hgcnbi: tecm less than two pi0 masses&')
1619 endif
60ed90b3 1620
9ef1c2d9 1621 kk=1
1622100 continue
1623
1624 if(kk.gt.20)then
1625 iret=1
1626 if(ish.ge.7)then
1627 write(ifch,*)'failed to generate hadron set for'
1628 *,' event:',nrevt+1
1629 write(ifch,*)'u d s :',keu,ked,kes,' E:',tecm
1630 write(ifch,*)('-',i=1,30)
1631 *,' exit sr hgcnbi ',('-',i=1,10)
1632 endif
1633 ish=isho
1634 return
1635 endif
1636
1637 amtot=0.0
1638 do i=1,nspecs
1639 nptlgc(i)=0
1640 enddo
1641 do ii=1,nflavs
1642 k(ii)=kef(ii)
1643 enddo
1644
1645 if(ish.ge.7)write(ifch,*)
60ed90b3 1646 *'sample hadron multiplicities and total mass:'
9ef1c2d9 1647
1648 kbar=keu+ked+kes
1649 kpar=iabs(keu)+iabs(ked)+iabs(kes)
1650 nbar=kbar/3.
1651 if(ish.ge.7)write(ifch,*)'baryon number:',nbar,' parton number:'
1652 *,kpar
1653
1654 nn=2
1655 if(ioinco.ne.2)then
1656 nn=hgcndn(0)
1657 else
1658 nn=nh
1659 endif
1660 nb=iabs(nbar)
1661 if(ish.ge.7)write(ifch,*)'<n>:',nh,' n_sample:',nn,' n_bar:',nb
1662 if(nn.gt.nb.and.nb.ne.0.and.nb.ge.nh)nn=nb
1663 if(nn.lt.nb.and.nb.ne.0)nn=nb
1664 km=kpar-iabs(kbar)
1665 nt=km/2+nb
1666 if(nt.gt.nn)nn=nt
1667 nn=max(nn,2)
1668
1669 if(ioinco.eq.2)then
1670 nit=15*taue
1671 else
1672 itpn=100
1673 nit=nn*itpn
1674 endif
1675 nbb=0
1676 n=0
1677
1678c start with nb protons
1679 nptlgc(19)=nptlgc(19)+nb
1680 n=nb
1681 amtot=amtot+nb*aspecs(19)
1682 do ii=1,nflavs
1683 k(ii)=k(ii)-ifok(ii,19)*nb
1684 enddo
1685 nbb=nbb+nb
1686
1687
1688 do it=1,nit
1689
1690 xsp=nspecs
1691 x0=0.5
1692 xib=x0+xsp*rangen()
1693 ib=nint(xib)
1694 if(ib.gt.nspecs)ib=nspecs
1695 if(ib.lt.1)ib=1
1696 kb=ifok(1,ib)+ifok(2,ib)+ifok(3,ib)
1697 if(rangen().lt.0.5.and.nptlgc(ib).ge.1)then
1698 ni=-1
1699 else
1700 ni=1
1701 endif
1702 as=1.0
1703 if(nptlgc(ib).eq.0)as=0.5
1704 if(nptlgc(ib).eq.1.and.ni.eq.(-1))as=2.0
1705 if(ish.ge.9)write(ifch,*)
1706 *'id:',ispecs(ib),' <i>:',ptlngc(ib),' ni:',ni
1707
1708 if(ni.ne.0)then
1709
1710 if(ptlngc(ib).gt.5.0)then
1711
1712 pnla=hgcpnl(ib,0)
1713 pnlb=hgcpnl(ib,ni)
1714 pnlog=-pnla+pnlb
1715 if(ish.ge.9)write(ifch,*)'pnlog:',pnlog
1716 if(pnlog.lt.60)then
1717 pn=exp(pnlog)
60ed90b3 1718 else
9ef1c2d9 1719 pn=1.1
1720 endif
1721
1722 else
1723
1724 if(ni.eq.1)then
1725 pn=ptlngc(ib)/(nptlgc(ib)+1)
1726 elseif(ni.eq.(-1).and.ptlngc(ib).gt.1.e-20)then
1727 pn=nptlgc(ib)/ptlngc(ib)
1728 elseif(nptlgc(ib).gt.0)then
1729 pn=1.1
1730 else
1731 pn=0.0
1732 endif
1733
1734 endif
1735
1736 pm=1.0
1737 if(ioinfl.ge.0)then
1738 pmla=hgcpml(ib,0,ib,0)
1739 pmlb=hgcpml(ib,ni,ib,0)
1740 pmlog=-pmla+pmlb
1741 if(ish.ge.9)write(ifch,*)'pmlog:',pmlog
1742 if(pmlog.lt.60)then
1743 pm=exp(pmlog)
1744 else
1745 pm=1.1
1746 endif
1747 endif
1748
1749 p=pn*pm*as
1750 r=rangen()
1751 if(r.le.p)then
1752 nptlgc(ib)=nptlgc(ib)+ni
1753 n=n+ni
1754 amtot=amtot+ni*aspecs(ib)
1755 do ii=1,nflavs
1756 k(ii)=k(ii)-ifok(ii,ib)*ni
1757 enddo
1758 if(kb.ne.0)nbb=nbb+ni
1759 if(ish.ge.7.and.ni.gt.0)write(ifch,*)'add:'
1760 if(ish.ge.7.and.ni.lt.0)write(ifch,*)'remove:'
1761 if(ish.ge.7)write(ifch,*)'id:',ispecs(ib),' <n_i>:',ptlngc(ib)
1762 *,' n_i:',nptlgc(ib)
1763 if(ish.ge.7)write(ifch,*)'<n>:',nn,' it:',it
1764 if(ish.ge.7)write(ifch,*)'<M>:',amgc,' M:',amtot
1765 if(ish.ge.7)write(ifch,*)'p:',p,' r:',r
1766 if(ish.ge.7)write(ifch,*)'flav defect: u:',k(1),' d:'
1767 *,k(2),' s:',k(3)
1768 if(n.ge.nn.and.ioinco.ne.2)goto102
1769 endif
1770
1771 endif
1772
1773 enddo
1774
1775
1776102 continue
60ed90b3 1777
9ef1c2d9 1778 ndd=0
1779c if(nbb.lt.nb)then
1780c nba=nb-nbb
1781c if(nbar.gt.0)then
60ed90b3 1782c if(ish.ge.7)write(ifch,*)'add protons: nba:',nba
9ef1c2d9 1783c nptlgc(19)=nptlgc(19)+nba
1784c n=n+nba
1785c amtot=amtot+aspecs(19)*nba
1786c elseif(nbar.lt.0)then
60ed90b3 1787c if(ish.ge.7)write(ifch,*)'add aprotons: nba:',nba
9ef1c2d9 1788c nptlgc(20)=nptlgc(20)+nba
1789c n=n+nba
1790c amtot=amtot+aspecs(20)*nba
1791c endif
1792c endif
1793 if(n.lt.nn.and.ioinco.ne.2)then
1794 ndd=nn-n
1795 nd=mod(ndd,4)
1796 xn=n
1797 xnn=nn
1798 xl=(xnn-xn)/4.
1799 l=aint(xl)
1800 if(ish.ge.7)write(ifch,*)'add pions/etas: ndd:',ndd
1801 *,' l:',l,' nd:',nd
1802 if(l.ge.1)then
1803 do j=1,l
1804 nptlgc(1)=nptlgc(1)+1
1805 nptlgc(2)=nptlgc(2)+1
1806 nptlgc(3)=nptlgc(3)+1
1807 nptlgc(8)=nptlgc(8)+1
1808 amtot=amtot+aspecs(1)+aspecs(2)+aspecs(3)+aspecs(8)
1809 enddo
1810 endif
1811 if(nd.eq.1)then
1812 nptlgc(1)=nptlgc(1)+1
1813 amtot=amtot+aspecs(1)
1814 elseif(nd.eq.2)then
1815 nptlgc(2)=nptlgc(2)+1
1816 nptlgc(3)=nptlgc(3)+1
1817 amtot=amtot+aspecs(2)+aspecs(3)
1818 elseif(nd.eq.3)then
1819 nptlgc(2)=nptlgc(2)+1
1820 nptlgc(3)=nptlgc(3)+1
1821 nptlgc(1)=nptlgc(1)+1
1822 amtot=amtot+aspecs(2)+aspecs(3)+aspecs(1)
1823 endif
1824 endif
1825
1826 if(n.eq.0.and.ioinco.eq.2)then
1827 nptlgc(2)=nptlgc(2)+1
1828 nptlgc(3)=nptlgc(3)+1
1829 amtot=amtot+aspecs(2)+aspecs(3)
1830 elseif(n.eq.1.and.ioinco.eq.2)then
1831 nptlgc(1)=nptlgc(1)+1
1832 amtot=amtot+aspecs(1)
1833 endif
1834
1835 if(amtot.ge.tecm.and.ioinfl.ge.0)then
1836 if(ish.ge.7)write(ifch,*)
1837 *'total mass exceeded , redo configuration'
1838 kk=kk+1
1839 goto100
1840 endif
60ed90b3 1841
9ef1c2d9 1842
1843 iii=0
1844 if(ish.ge.7)then
1845 write(ifch,*)'u d s :',keu,ked,kes,' E:',tecm
1846 write(ifch,*)
1847 *'hadron set without flavor conservation:'
1848 endif
1849 do i=1,nspecs
1850 n=nptlgc(i)
1851 if(n.ge.1)then
1852 do j=1,n
1853 iii=iii+1
1854 if(iii.gt.maxp)stop'iii>maxp in hgcnbi'
1855 idp(iii)=ispecs(i)
1856 enddo
1857 endif
1858 enddo
1859 if(ish.ge.7)then
1860 write(ifch,'(1x,10i6)')(idp(i),i=1,iii)
1861 write(ifch,*)'flav defect: u:',k(1),' d:'
1862 *,k(2),' s:',k(3)
1863 write(ifch,*)'M:',amtot,' <M>:',amgc
1864 endif
1865 if(ioinfl.le.0)goto1000
1866
1867 ll=1
1868 llmax=nn*25
1869 ior=1
1870
1871120 if(k(1).ne.0.or.k(2).ne.0.or.k(3).ne.0)then
1872
1873 if(kk.gt.6)ior=0
1874
1875 if(ish.ge.7)write(ifch,*)
1876 *'remaining flavor defect before operation:',ll
1877 if(ish.ge.7)write(ifch,*)'flav defect: u:',k(1),' d:'
1878 *,k(2),' s:',k(3)
1879
1880 nida=0
1881 do i=1,nspecs
1882 if(nptlgc(i).gt.0)then
1883 nida=nida+1
1884 ida(nida)=i
1885 endif
1886 enddo
1887
1888 if(nida.eq.0)then
1889 if(ish.ge.7)write(ifch,*)'no proposals in a , redo'
1890 kk=kk+1
1891 goto100
1892 endif
1893
1894
1895 xna=0.5+nida*rangen()
1896 na=nint(xna)
1897 if(na.gt.nida)na=nida
1898 if(na.lt.1)na=1
1899 ia=ida(na)
1900 if(ish.ge.7)write(ifch,*)'nida:',nida,' ia:',ia
1901
1902 nidb=0
1903 do ii=1,nflavs
1904 kp(ii)=k(ii)+ifok(ii,ia)
1905 kps(ii)=isign(1,kp(ii))
1906 enddo
1907 if(ish.ge.7)write(ifch,*)
1908 *' assemble: u:',kp(1),' d:',kp(2),' s:',kp(3)
1909 do i=1,nspecs
1910 iacc=0
1911 naccsp=0
1912 naccmi=1
1913 do ii=1,nflavs
1914 naccsp=naccsp+iabs(ifok(ii,i))
1915 if(kp(ii).ne.0)then
1916 if(kps(ii)*ifok(ii,i).le.kps(ii)*kp(ii)
1917 *.and.kps(ii)*ifok(ii,i).gt.0)iacc=iacc+iabs(ifok(ii,i))
1918 endif
1919 enddo
1920 if(kp(1).eq.0.and.kp(2).eq.0.and.kp(3).eq.0)naccmi=0
1921 if(iacc.eq.naccsp.and.naccsp.ge.naccmi)then
1922 nidb=nidb+1
1923 idb(nidb)=i
1924 endif
1925 enddo
1926
1927 if(nidb.eq.0)then
1928 if(ish.ge.7)write(ifch,*)'no proposals in b , redo'
1929 kk=kk+1
1930 goto100
1931 endif
1932
1933 xnb=0.5+nidb*rangen()
1934 nb=nint(xnb)
1935 if(nb.gt.nidb)nb=nidb
1936 if(nb.lt.1)nb=1
1937 ib=idb(nb)
1938 if(ish.ge.7)write(ifch,*)'nidb:',nidb,' ib:',ib
1939 if(ish.ge.7)write(ifch,*)
1940 *'proposal:',ispecs(ia),' --> ',ispecs(ib)
1941
1942 asym=1.0
1943
1944c if(asym.gt.0.0)then
1945
1946 if(ptlngc(ia).gt.5.0)then
1947 pnali=hgcpnl(ia,0)
1948 pnalf=hgcpnl(ia,-1)
1949 pnalog=-pnali+pnalf
1950 if(ish.ge.7)write(ifch,*)'pnalog:',pnalog
1951 if(pnalog.lt.60)then
1952 pna=exp(pnalog)
1953 else
1954 pna=1.1
1955 endif
1956 else
1957 if(ptlngc(ia).gt.1.e-20)then
1958 pna=nptlgc(ia)/ptlngc(ia)
1959 elseif(nptlgc(ia).gt.0)then
1960 pna=1.1
1961 else
1962 pna=0.0
1963 endif
1964 endif
60ed90b3 1965
9ef1c2d9 1966 if(ptlngc(ib).gt.5.0)then
1967 pnbli=hgcpnl(ib,0)
1968 pnblf=hgcpnl(ib,1)
1969 pnblog=-pnbli+pnblf
1970 if(ish.ge.7)write(ifch,*)'pnblog:',pnblog
1971 if(pnblog.lt.60)then
1972 pnb=exp(pnblog)
1973 else
1974 pnb=1.1
1975 endif
60ed90b3 1976 else
9ef1c2d9 1977 pnb=ptlngc(ib)/(nptlgc(ib)+1)
1978 endif
1979
1980
1981 pmli=hgcpml(ia,0,ib,0)
1982 pmlf=hgcpml(ia,-1,ib,1)
1983 pmlog=-pmli+pmlf
1984 if(ish.ge.7)write(ifch,*)'pmlog:',pmlog
1985 if(pmlog.lt.60)then
1986 pm=exp(pmlog)
1987 else
1988 pm=1.1
1989 endif
1990
1991 p=pna*pnb*pm*asym
1992 if(ior.eq.0)then
1993 r=0.0
1994 else
1995 r=rangen()
1996 endif
60ed90b3 1997
9ef1c2d9 1998c else
1999
2000c r=1.0
2001c p=0.0
2002
2003c endif
2004
2005 if(r.lt.p)then
2006 if(ish.ge.7)write(ifch,*)'p:',p,' r:',r,' asymmetry:',asym
2007 if(ish.ge.7)write(ifch,*)'remove ',ispecs(ia),' add ',ispecs(ib)
2008 *,' proposal accepted'
2009 nptlgc(ia)=nptlgc(ia)-1
2010 nptlgc(ib)=nptlgc(ib)+1
2011 amtot=amtot-aspecs(ia)+aspecs(ib)
2012 do ii=1,nflavs
2013 k(ii)=k(ii)+ifok(ii,ia)-ifok(ii,ib)
2014 enddo
2015 endif
2016
60ed90b3 2017
9ef1c2d9 2018 if(k(1).ne.0.or.k(2).ne.0.or.k(3).ne.0)then
2019 ll=ll+1
2020 if(ll.le.llmax)then
2021 goto120
2022 else
2023 if(ish.ge.7)write(ifch,*)'failed to remove defect, redo'
2024 kk=kk+1
2025 goto100
2026 endif
2027 endif
2028
2029 endif
2030
20311000 continue
2032
2033 nump=0
2034 kcu=0
2035 kcd=0
2036 kcs=0
2037 do i=1,nspecs
2038 n=nptlgc(i)
2039 if(n.ge.1)then
2040 do j=1,n
2041 nump=nump+1
2042 ihadro(nump)=ispecs(i)
2043 kcu=kcu+ifok(1,i)
2044 kcd=kcd+ifok(2,i)
2045 kcs=kcs+ifok(3,i)
2046 enddo
2047 endif
2048 enddo
2049
2050 if(ioinfl.gt.0)then
2051 if(kcu.ne.keu.or.kcd.ne.ked.or.kcs.ne.kes)then
2052 if(ish.ge.7)write(ifch,*)
2053 *'failed to remove flavor defect, redo configuration'
2054 kk=kk+1
2055 goto100
2056 endif
2057 endif
2058
2059 if(ioinct.ge.1)then
2060 chitot=0.0
2061 nutot=nspecs
2062 do i=1,nspecs
2063 chi=0.0
2064 if(rmsngc(i).gt.1.e-10)chi=(ptlngc(i)-nptlgc(i))/rmsngc(i)
2065 chitot=chitot+chi**2
2066 enddo
2067 call xhgccc(chitot)
60ed90b3 2068
9ef1c2d9 2069 u=0
2070 d=0
2071 s=0
2072 do i=1,nspecs
2073 u=u+ifok(1,i)*nptlgc(i)
2074 d=d+ifok(2,i)*nptlgc(i)
2075 s=s+ifok(3,i)*nptlgc(i)
2076 enddo
2077 call xhgcfl(u,d,s,0)
2078 call xhgcam(amtot,0)
2079 endif
2080
2081 if(ish.ge.7)then
2082 write(ifch,*)
2083 *'initial hadron set for droplet decay:'
2084 write(ifch,'(1x,10i6)')(ihadro(i),i=1,nump)
2085 endif
2086 if(nump.ge.nlattc)then
2087 nlattc=nump+1
2088 if(ish.ge.7)then
2089 write(ifch,*)'initial set > nlattc !'
2090 write(ifch,*)'new nlattc:',nlattc
2091 endif
2092 endif
2093 if(ish.ge.7)then
2094 write(ifch,*)'keu:',kef(1),' kcu:',kcu,' ku:',k(1)
2095 write(ifch,*)'ked:',kef(2),' kcd:',kcd,' kd:',k(2)
2096 write(ifch,*)'kes:',kef(3),' kcs:',kcs,' ks:',k(3)
2097 write(ifch,*)' nh:',nh,' nump:',nump
2098 write(ifch,*)' nu:',nutot,' chi^2:',chitot
2099 write(ifch,*)'iozero:',iozero,' iomom:',iomom
2100 write(ifch,*)
2101 *'total mass:',amtot,' droplet mass:',tecm
2102 write(ifch,*)'trials needed:',kk
2103 *,' operations needed:',ll
2104 write(ifch,*)'iterations:',it,' pions added:',ndd
2105 write(ifch,*)('-',i=1,30)
2106 *,' exit sr hgcnbi ',('-',i=1,10)
2107 endif
2108 ish=isho
2109 return
2110
2111 end
2112
2113c--------------------------------------------------------------------
2114 integer function hgcndn(i)
2115c--------------------------------------------------------------------
2116c returns random multiplicity from gaussian distribution for species i
2117c---------------------------------------------------------------------
2118 include 'epos.inc'
60ed90b3 2119 parameter (mspecs=56)
9ef1c2d9 2120 common/cgchg/rmsngc(mspecs),ptlngc(mspecs),chemgc(mspecs),tem
2121 common/cgctot/rmstot,ptltot
2122 common/clatt/nlattc,npmax
2123 a=iowidn
2124 kk=0
2125
2126 if(i.eq.0)then
2127
21281 continue
2129 kk=kk+1
2130 p=0.0
2131 nmin=2
2132 nh=nint(ptltot)
2133 nmax=nlattc
2134 xn=1.5+(nmax-nmin)*rangen()
2135 n=nint(xn)
2136 x=(n-ptltot)**2/2.0
2137 y=-70.
2138 if(rmstot.gt.1.e-15)y=-x/rmstot**2*a**2
2139 if(y.lt.70.)p=exp(y)
2140 if(rmstot.gt.1.e-15.and.iowidn.lt.0)p=p/sqrt(2.*pi)/rmstot
2141 if(p.ge.rangen())then
2142 hgcndn=n
2143 if(ish.ge.9)write(ifch,*)'hgcndn: k:',kk,' n:',hgcndn
2144 return
2145 else
2146 if(kk.le.25)goto1
2147 hgcndn=max(2,nh)
2148 if(ish.ge.9)write(ifch,*)'hgcndn: k:',kk,' n:',hgcndn
2149 return
2150 endif
2151
2152 else
60ed90b3 2153
9ef1c2d9 21542 continue
2155 kk=kk+1
2156 p=0.0
2157 nmin=0
2158 nh=nint(ptlngc(i))
2159 nmax=2*nh
2160 nmax=max(2,nmax)
2161 xn=-0.5+(nmax-nmin)*rangen()
2162 n=nint(xn)
2163 x=(n-ptlngc(i))**2/2.0
2164 if(x.lt.1.e-30)then
2165 p=1.
2166 else
2167 y=-70.
2168 if(rmsngc(i).gt.1.e-15)y=-x/rmsngc(i)**2
2169 if(y.lt.70.)p=exp(y)
2170 if(rmsngc(i).gt.1.e-15.and.iowidn.lt.0)
2171 *p=p/sqrt(2.*pi)/rmsngc(i)
2172 endif
2173 if(p.ge.rangen())then
2174 hgcndn=n
2175 if(ish.ge.9)write(ifch,*)'hgcndn: k:',kk,' n:',hgcndn
2176 return
2177 else
2178 if(kk.le.25)goto2
2179 hgcndn=nh
2180 if(ish.ge.9)write(ifch,*)'hgcndn: k:',kk,' n:',hgcndn
2181 return
2182 endif
2183
2184 endif
2185
2186 end
2187
2188c--------------------------------------------------------------------
2189 function hgcpml(i1,n1,i2,n2)
2190c--------------------------------------------------------------------
2191 include 'epos.inc'
60ed90b3 2192 parameter (mspecs=56)
9ef1c2d9 2193 common/cspecs/nspecs,ispecs(mspecs),aspecs(mspecs),gspecs(mspecs)
2194 common/camgc/amgc,samgc,amtot
2195 common/cgcnb/nptlgc(mspecs)
2196 if(ish.ge.9)write(ifch,*)'i1:',i1,' i2:',i2
2197 if(ish.ge.9)write(ifch,*)'n1:',n1,' n2:',n2
2198 hgcpml=-1.e30
2199 ampr=n1*aspecs(i1)+n2*aspecs(i2)
2200 if((amtot+ampr).lt.tecm.and.(amtot+ampr).ge.0
2201 *.and.nptlgc(i1).ge.(-n1).and.nptlgc(i2).ge.(-n2))then
2202 hgcpml=0.0
2203 pl=(amtot-amgc+ampr)**2/2.0
2204 if(pl.lt.1.e-30)then
2205 hgcpml=0.0
2206 return
2207 endif
2208 if(samgc.gt.1.e-15)hgcpml=-pl/samgc**2
2209 endif
2210 if(ish.ge.9)write(ifch,*)'hgcpml:',hgcpml
2211 return
2212 end
2213
2214c--------------------------------------------------------------------
2215 function hgcpnl(i,n)
2216c--------------------------------------------------------------------
2217 include 'epos.inc'
60ed90b3 2218 parameter (mspecs=56)
9ef1c2d9 2219 common/cgchg/rmsngc(mspecs),ptlngc(mspecs),chemgc(mspecs),tem
2220 common/cgcnb/nptlgc(mspecs)
2221 if(ish.ge.9)write(ifch,*)'i:',i,' n:',n
2222 hgcpnl=-1.e30
2223 if(nptlgc(i).ge.(-n))then
2224 pl=(nptlgc(i)-ptlngc(i)+n)**2/2.0
2225 if(pl.lt.1.e-30)then
2226 hgcpnl=0.0
2227 return
2228 endif
2229 if(rmsngc(i).gt.1.e-15)hgcpnl=-pl/rmsngc(i)**2
2230 endif
2231 if(ish.ge.9)write(ifch,*)'hgcpnl:',hgcpnl
2232 return
2233 end
2234
2235
2236c--------------------------------------------------------------------
2237 subroutine hgcpen
2238c--------------------------------------------------------------------
2239c returns array for twodimensional plot of energy- and flavor-
2240c density
2241c--------------------------------------------------------------------
2242c xpar1,xpar2 temperature range
2243c xpar3 # of bins for temperature
2244c xpar4,xpar5 chem.pot. range
2245c xpar6 # of bins for chem.pot.
2246c xpar7 max. density
2247c xpar8 strange chem.pot.
2248c--------------------------------------------------------------------
2249 include 'epos.inc'
60ed90b3 2250 parameter (mspecs=56)
9ef1c2d9 2251 common/cspecs/nspecs,ispecs(mspecs),aspecs(mspecs),gspecs(mspecs)
2252 common/cgchg/rmsngc(mspecs),ptlngc(mspecs),chemgc(mspecs),tem
2253 common/cflavs/nflavs,kef(nflav),chem(nflav)
2254 common/cflac/ifok(nflav,mspecs),ifoa(nflav)
2255 common/ciakt/gen,iafs,ians,genm
2256 parameter (nbin=100)
2257 real edensi(nbin,nbin),qdensi(nbin,nbin)
2258 external hgcfhe
2259 external hgcfhn
2260 external hgcfbe
2261 external hgcfbn
2262
2263 iug=(1+iospec)/2*2-1
2264
2265c initialization
2266c --------------
2267
2268 if(iug.eq.1)nflavs=1
2269 if(iug.eq.3)nflavs=2
2270 if(iug.eq.5)nflavs=2
2271 if(iug.eq.7)nflavs=3
2272 if(iug.eq.9)nflavs=3
2273 if(iug.eq.11)nflavs=3
2274 tem=0.0
2275 do i=1,nflavs
2276 chem(i)=0.0
2277 enddo
2278 call hgchac(0)
2279 do i=1,nspecs
2280 ptlngc(i)=0.0
2281 rmsngc(i)=0.0
2282 enddo
2283
2284 nbt=nint(xpar3)
2285 nbc=nint(xpar6)
2286 nbc=min(nbc,100)
2287 nbt=min(nbt,100)
2288 dt=(xpar2-xpar1)/nbt
2289 dc=(xpar5-xpar4)/nbc
2290 ymax=xpar7
2291 cs=xpar8
2292
2293
2294 t0=xpar1+dt/2.
2295 c0=xpar4+dc/2
2296 do i=1,nbc
2297 chem(1)=c0+(i-1)*dc
2298 chem(2)=chem(1)
2299 chem(3)=cs
2300 chem(4)=0.0
2301 chem(5)=0.0
2302 chem(6)=0.0
2303 call hgchac(0)
2304 do ii=1,nbt
2305 tem=t0+(ii-1)*dt
2306 if(ish.ge.5)write(ifch,*)' mu:',chem(1),' T:',tem
2307
2308 qd=0.0
2309 ed=0.0
2310
2311 do ians=1,nspecs
2312
2313 call hgclim(a,b)
2314
2315 if(b.eq.0.0)then
2316 hden=0.0
2317 elseif(iostat.eq.0)then
2318 call uttraq(hgcfhn,a,b,hden)
2319 elseif(iostat.eq.1)then
2320 call uttraq(hgcfbn,a,b,hden)
2321 endif
2322 hd=hden*gspecs(ians)/2./pi**2/hquer**3
60ed90b3 2323
9ef1c2d9 2324 if(ish.ge.7)write(ifch,*)'i:',ians,' n_u:',ifok(1,ians),' hd:',hd
2325
60ed90b3 2326 qd=qd+ifok(1,ians)*hd+ifok(2,ians)*hd
9ef1c2d9 2327 if(qd.gt.ymax)qd=ymax
2328c if(qd.gt.ymax)qd=0.0
2329 if(qd.lt.-ymax)qd=-ymax
2330c if(qd.lt.-ymax)qd=0.0
60ed90b3 2331
9ef1c2d9 2332
2333 if(b.eq.0.0)then
2334 edi=0.0
2335 elseif(iostat.eq.0)then
2336 call uttraq(hgcfhe,a,b,edi)
2337 elseif(iostat.eq.1)then
2338 call uttraq(hgcfbe,a,b,edi)
2339 endif
2340 edi=edi*gspecs(ians)/2./pi**2/hquer**3
2341
2342 if(ish.ge.7)write(ifch,*)'i:',ians,' mu:',chemgc(ians)
2343 * ,' edi:',edi
60ed90b3 2344
9ef1c2d9 2345 ed=ed+edi
2346 if(ed.gt.ymax)ed=ymax
2347c if(ed.gt.ymax)ed=0.0
2348 enddo
2349
2350 if(ish.ge.5)write(ifch,*)' ed:',ed,' qd:',qd
2351 edensi(i,ii)=ed
2352 qdensi(i,ii)=qd
60ed90b3 2353
9ef1c2d9 2354 enddo
2355 enddo
2356
2357 write(ifhi,'(a)') 'openhisto'
2358 write(ifhi,'(a,2e11.3)')'xrange',xpar1,xpar2
2359 write(ifhi,'(a,2e11.3)')'yrange',xpar4,xpar5
2360 write(ifhi,'(a)') 'set ityp2d 5'
2361 write(ifhi,'(a,i4)') 'array2d',nbt
2362 do j=1,nbc
2363 do jj=1,nbt
2364 write(ifhi,'(e11.3)') edensi(j,jj)
2365 enddo
2366 enddo
2367 write(ifhi,'(a)') ' endarray'
2368 write(ifhi,'(a)') 'closehisto plot2d'
2369
2370 write(ifhi,'(a)') 'openhisto'
2371 write(ifhi,'(a,2e11.3)')'xrange',xpar1,xpar2
2372 write(ifhi,'(a,2e11.3)')'yrange',xpar4,xpar5
2373 write(ifhi,'(a)') 'set ityp2d 5'
2374 write(ifhi,'(a,i4)') 'array2d',nbt
2375 do j=1,nbc
2376 do jj=1,nbt
2377 write(ifhi,'(e11.3)') qdensi(j,jj)
2378 enddo
2379 enddo
2380 write(ifhi,'(a)') ' endarray'
2381 write(ifhi,'(a)') 'closehisto plot2d'
2382
2383 return
2384 end
2385
2386c--------------------------------------------------------------------
2387 subroutine hgcpfl
2388c--------------------------------------------------------------------
2389c returns array for twodimensional plot of energy- and flavor-
2390c density fluctuations
2391c--------------------------------------------------------------------
2392c xpar1,xpar2 temperature range
2393c xpar3 # of bins for temperature
2394c xpar4,xpar5 chem.pot. range
2395c xpar6 # of bins for chem.pot.
2396c xpar7 max. density
2397c xpar8 strange chem.pot.
2398c--------------------------------------------------------------------
2399 include 'epos.inc'
60ed90b3 2400 parameter (mspecs=56)
9ef1c2d9 2401 common/cspecs/nspecs,ispecs(mspecs),aspecs(mspecs),gspecs(mspecs)
2402 common/cgchg/rmsngc(mspecs),ptlngc(mspecs),chemgc(mspecs),tem
2403 common/cflavs/nflavs,kef(nflav),chem(nflav)
2404 common/ciakt/gen,iafs,ians,genm
2405 parameter (nbin=100)
2406 real efl(nbin,nbin),qfl(nbin,nbin),v(nbin),wn(nbin),we(nbin)
2407 external hgcfhf
2408 external hgcfhe
2409 external hgcfhn
2410 external hgcfhw
2411 external hgcfbf
2412 external hgcfbe
2413 external hgcfbn
2414
2415 iug=(1+iospec)/2*2-1
2416
2417c initialization
2418c --------------
2419
2420 if(iug.eq.1)nflavs=1
2421 if(iug.eq.3)nflavs=2
2422 if(iug.eq.5)nflavs=2
2423 if(iug.eq.7)nflavs=3
2424 if(iug.eq.9)nflavs=3
2425 if(iug.eq.11)nflavs=3
2426 tem=0.0
2427 do i=1,nflavs
2428 chem(i)=0.0
2429 enddo
2430 call hgchac(0)
2431 do i=1,nspecs
2432 ptlngc(i)=0.0
2433 rmsngc(i)=0.0
2434 enddo
2435
2436 nbt=nint(xpar3)
2437 nbv=nint(xpar6)
2438 nbv=min(nbv,100)
2439 nbt=min(nbt,100)
2440 dt=(xpar2-xpar1)/nbt
2441 dv=(xpar5-xpar4)/nbv
2442 ymax=1.e20
2443 chem(1)=xpar7
2444 chem(2)=xpar7
2445 chem(3)=xpar8
2446 call hgchac(0)
2447
2448
2449 t0=xpar1+dt/2.
2450 v0=xpar4
2451 do i=1,nbv
2452 volu=v0+(i-1)*dv
2453 do ii=1,nbt
2454 tem=t0+(ii-1)*dt
2455 if(ish.ge.5)write(ifch,*)'volu:',volu,' tem:',tem
2456
2457 ev=0.0
2458 ee=0.0
2459 qv=0.0
2460 qe=0.0
2461
2462 do ians=1,nspecs
2463
2464 call hgclim(a,b)
2465
2466 if(b.eq.0.0)then
2467 hn=0.0
2468 hv=0.0
2469 elseif(iostat.eq.0)then
2470 call uttraq(hgcfhn,a,b,hn)
2471 call uttraq(hgcfhw,a,b,hv)
2472 elseif(iostat.eq.1)then
2473 call uttraq(hgcfbn,a,b,hn)
2474 hv=hn
2475 endif
2476 hn=hn*volu*gspecs(ians)/2./pi**2/hquer**3
2477 hv=hv*volu*gspecs(ians)/2./pi**2/hquer**3
2478 if(ish.ge.5)write(ifch,*)'hn:',hn,' hv:',hv
60ed90b3 2479
9ef1c2d9 2480 hn=max(hn,1.e-15)
2481 qv=qv+hv
2482 qe=qe+hn
2483
2484
2485 if(qv.gt.ymax)qv=ymax
2486 if(qe.gt.ymax)qe=ymax
60ed90b3 2487
9ef1c2d9 2488
2489 if(b.eq.0.0)then
2490 eei=0.0
2491 evi=0.0
2492 elseif(iostat.eq.0)then
2493 call uttraq(hgcfhe,a,b,eei)
2494 call uttraq(hgcfhf,a,b,evi)
2495 elseif(iostat.eq.1)then
2496 call uttraq(hgcfbe,a,b,eei)
2497 call uttraq(hgcfbf,a,b,evi)
2498 endif
2499 eei=eei*volu*gspecs(ians)/2./pi**2/hquer**3
2500 evi=evi*volu*gspecs(ians)/2./pi**2/hquer**3
2501 if(ish.ge.5)write(ifch,*)'eei:',eei,' evi:',evi
2502
60ed90b3 2503
9ef1c2d9 2504 eei=max(eei,1.e-15)
2505 ev=ev+evi
2506 ee=ee+eei
2507 if(ev.gt.ymax)ev=ymax
2508 if(ee.gt.ymax)ee=ymax
2509 enddo
2510 if(ish.ge.5)write(ifch,*)'qv:',qv,' ev:',ev
2511
2512 qfl(i,ii)=0.
2513 efl(i,ii)=0.
2514 if(ev.gt.0.0.and.ee.gt.1.e-15)efl(i,ii)=sqrt(ev)/ee
2515 if(qv.gt.0.0.and.ee.gt.1.e-15)qfl(i,ii)=sqrt(qv)/qe
2516 if(tem.eq.0.195)then
2517 we(i)=efl(i,ii)
2518 wn(i)=qfl(i,ii)
2519 v(i)=volu
2520 endif
60ed90b3 2521
9ef1c2d9 2522 enddo
2523 enddo
2524
2525 write(ifhi,'(a)') 'openhisto'
2526 write(ifhi,'(a,2e11.3)')'xrange',xpar1,xpar2
2527 write(ifhi,'(a,2e11.3)')'yrange',xpar4,xpar5
2528 write(ifhi,'(a)') 'set ityp2d 5'
2529 write(ifhi,'(a,i4)') 'array2d',nbt
2530 do j=1,nbv
2531 do jj=1,nbt
2532 write(ifhi,'(e11.3)') efl(j,jj)
2533 enddo
2534 enddo
2535 write(ifhi,'(a)') ' endarray'
2536 write(ifhi,'(a)') 'closehisto plot2d'
2537
2538 write(ifhi,'(a)') 'openhisto'
2539 write(ifhi,'(a,2e11.3)')'xrange',xpar1,xpar2
2540 write(ifhi,'(a,2e11.3)')'yrange',xpar4,xpar5
2541 write(ifhi,'(a)') 'set ityp2d 5'
2542 write(ifhi,'(a,i4)') 'array2d',nbt
2543 do j=1,nbv
2544 do jj=1,nbt
2545 write(ifhi,'(e11.3)') qfl(j,jj)
2546 enddo
2547 enddo
2548 write(ifhi,'(a)') ' endarray'
2549 write(ifhi,'(a)') 'closehisto plot2d'
2550
2551 write(ifhi,'(a)') 'newpage zone 1 2 1'
2552 write(ifhi,'(a)') 'openhisto'
2553 write(ifhi,'(a,2e11.3)')'xrange',xpar4,xpar5
2554 write(ifhi,'(a)') 'htyp lfu xmod lin ymod lin'
2555 write(ifhi,'(a,i4)') 'array 2'
2556 do j=1,nbv
2557 write(ifhi,'(2e13.5)')v(j),we(j)
2558 enddo
2559 write(ifhi,'(a)') ' endarray'
60ed90b3 2560 write(ifhi,'(a)') 'closehisto plot 0'
9ef1c2d9 2561
2562 write(ifhi,'(a)') 'openhisto'
2563 write(ifhi,'(a,2e11.3)')'xrange',xpar4,xpar5
2564 write(ifhi,'(a)') 'htyp lfu xmod lin ymod lin'
2565 write(ifhi,'(a,i4)') 'array 2'
2566 do j=1,nbv
2567 write(ifhi,'(2e13.5)')v(j),wn(j)
2568 enddo
2569 write(ifhi,'(a)') ' endarray'
2570 write(ifhi,'(a)') 'closehisto plot 0'
2571
2572
2573 return
2574 end
2575
2576
2577c------------------------------------------------------------------
2578 subroutine hgcpyi(ist)
2579c------------------------------------------------------------------
60ed90b3 2580c returns particle yield
9ef1c2d9 2581c input:
2582c tem : temperature
2583c chemgc: chemical potentials
2584c output:
2585c ptlngc: expectation value of particle number for each species
2586c rmsngc: standard deviation of ptlngc
2587c ptltot: total particle number
2588c rmstot: standard deviation of ptltot
2589c works for hadrons and partons
2590c ist=1 boltzmann statistics
2591c ist=0 quantum statistics
2592c--------------------------------------------------------------------
2593 include 'epos.inc'
60ed90b3 2594 parameter (mspecs=56)
9ef1c2d9 2595 common/cspecs/nspecs,ispecs(mspecs),aspecs(mspecs),gspecs(mspecs)
2596 common/cgchg/rmsngc(mspecs),ptlngc(mspecs),chemgc(mspecs),tem
2597 common/cgctot/rmstot,ptltot
2598 common/camgc/amgc,samgc,amtot
2599 common/ciakt/gen,iafs,ians,genm
2600 external hgcfhw
2601 external hgcfhn
2602
2603 if(iabs(ispecs(nspecs)).lt.10)then
2604
2605c parton yield
2606c ------------
2607 if(ish.ge.5)write(ifch,*)'parton yield:'
2608 gln=16.*1.20206*tem**3/pi**2*volu/hquer**3
2609 sdg=sqrt(gln) !!???
2610 if(ish.ge.5)write(ifch,'(1x,a,f10.4,2x,a,f9.4,a)')
2611 *'<N( 0)> :',gln,' sigma :',sdg,' (qm-statistics!)'
2612 ptltot=gln
2613 rmstot=0.0
2614 vartot=gln
2615
2616 else
2617
2618 if(ish.ge.5)write(ifch,*)'hadronic yield:'
2619 ptltot=0.0
2620 rmstot=0.0
2621 vartot=0.0
2622
2623 endif
2624
2625 amgc=0.0
2626 samgc=0.0
2627
2628 do ians=1,nspecs
2629
2630c hadronic yield
2631c --------------
2632 if(ist.eq.0)then
2633
2634 call hgclim(a,b)
2635 if(b.eq.0.0)then
2636 hden=0.0
2637 else
2638 call uttraq(hgcfhn,a,b,hden)
2639 endif
2640 ptlngc(ians)=hden*volu*gspecs(ians)/2./pi**2/hquer**3
2641
2642 else
2643
2644 if((chemgc(ians)/tem).gt.70.)then
2645 hpd=1.e30
2646 else
2647 hpd=exp(chemgc(ians)/tem)
2648 endif
2649 if(aspecs(ians).ne.0.)then
2650 fk2=hgcbk(2,aspecs(ians)/tem)
2651 hpd=hpd*gspecs(ians)*aspecs(ians)**2*tem*fk2
2652 */2./pi**2/hquer**3
2653 else
2654 hpd=hpd*gspecs(ians)*tem**3/pi**2/hquer**3
2655 endif
2656 ptlngc(ians)=hpd*volu
2657
2658 endif
2659
2660 ptltot=ptltot+ptlngc(ians)
2661 amgc=amgc+ptlngc(ians)*aspecs(ians)
2662 if(amgc.ge.tecm)amgc=tecm*0.9
2663
2664c standard deviation
2665c ------------------
2666 rmsngc(ians)=0.0
2667
2668 if(ist.eq.0)then
60ed90b3 2669
9ef1c2d9 2670 call uttraq(hgcfhw,a,b,var)
2671 var=var*gspecs(ians)*volu/2./pi**2/hquer**3
2672 vartot=vartot+var
2673 if(var.ge.0.0)rmsngc(ians)=sqrt(var)
2674 samgc=samgc+var*aspecs(ians)
60ed90b3 2675
9ef1c2d9 2676 else
2677
2678 if(ptlngc(ians).ge.0.0)rmsngc(ians)=sqrt(ptlngc(ians))
2679 vartot=vartot+ptlngc(ians)
2680 samgc=samgc+ptlngc(ians)*aspecs(ians)
2681
2682 endif
2683
2684
2685 if(ish.ge.7)write(ifch,'(2x,a,i5,a,2x,f8.4,5x,a,3x,f8.4)')
2686 *'m(',ispecs(ians),') :',aspecs(ians),'mu :',chemgc(ians)
2687 if(ish.ge.5)write(ifch,'(1x,a,i5,a,2x,f8.4,2x,a,2x,f10.4)')
2688 *'<N(',ispecs(ians),')> :',ptlngc(ians),'sigma :',rmsngc(ians)
2689
2690 enddo
60ed90b3 2691
9ef1c2d9 2692 if(vartot.ge.0.0)rmstot=sqrt(vartot)
2693 if(samgc.ge.0.0)samgc=sqrt(samgc)
2694 if(amgc.ge.tecm)samgc=sqrt(amgc)
2695 if(ish.ge.5)write(ifch,'(1x,a,2x,f8.4,2x,a,2x,f10.4)')
2696 *'<N( all)> :',ptltot,'sigma :',rmstot
2697 if(ish.ge.5)write(ifch,'(1x,a,2x,f8.4,2x,a,2x,f10.4)')
2698 *'<M_tot> :',amgc,'sigma :',samgc
2699
60ed90b3 2700 return
9ef1c2d9 2701 end
2702
2703c------------------------------------------------------------------------
2704 subroutine hgctbo(iba)
2705c------------------------------------------------------------------------
2706c returns new tem using boltzmann statistics in analytic form
2707c input:
2708c chemgc
2709c tecm/volu
2710c output:
2711c tem
2712c----------------------------------------------------------------------
2713 include 'epos.inc'
60ed90b3 2714 parameter (mspecs=56)
9ef1c2d9 2715 common/cspecs/nspecs,ispecs(mspecs),aspecs(mspecs),gspecs(mspecs)
2716 common/cgchg/rmsngc(mspecs),ptlngc(mspecs),chemgc(mspecs),tem
2717 common/ciakt/gen,iafs,ians,genm
2718 external hgcbk
2719 external hgcbk1
2720 iba=0
2721 k=1
2722 t1=0.0
2723 t2=1.0
60ed90b3 2724
9ef1c2d9 2725 goto15
2726
272710 tem=t1+.5*(t2-t1)
2728 if(tem.le.1.e-7)return
272915 eden=0.0
2730
2731 do i=1,nspecs
2732
2733 if(aspecs(i).ne.0)then
2734 if(tem.ne.0.)arr=aspecs(i)/tem
2735 cba=(aspecs(i)/tem+12.*tem/aspecs(i)-3.*chemgc(i)/aspecs(i))
2736 **hgcbk(2,arr)+(3.-chemgc(i)/tem)*hgcbk1(arr)
2737 else
2738 cba=4.*tem-chemgc(i)
2739 endif
2740
2741 if(cba.lt.0.0)then
2742 iba=1
2743 return
2744 endif
2745
2746 if(tem.ne.0.)x=chemgc(i)/tem
2747
2748 if(x.le.70.)then
2749 y=exp(x)
2750 else
2751 y=1.e30
2752 endif
2753
2754 if(aspecs(i).ne.0.)then
2755 edi=y*(3./arr*hgcbk(2,arr)+hgcbk1(arr))
2756 **gspecs(i)*aspecs(i)**3*tem/2./pi**2/hquer**3
2757 else
2758 edi=y*3.*gspecs(i)*tem**4/pi**2/hquer**3
2759 endif
2760
2761 eden=eden+edi
2762
2763 enddo
2764
2765 if(iabs(ispecs(nspecs)).lt.10)
2766 *eden=eden+(8.*pi**2*tem**4/15.+bag4rt**4)/hquer**3
2767
2768 de=abs(eden-(tecm/volu))
2769 if(de.le.gen*(tecm/volu).or.de.le.genm)return
2770c if(eden.ge.100.)return
2771
2772 if(eden.gt.(tecm/volu))then
2773 t2=tem
2774 else
2775 t1=tem
2776 endif
2777
2778 if(k.gt.300)return
60ed90b3 2779
9ef1c2d9 2780 k=k+1
2781 goto10
2782 end
2783
2784c----------------------------------------------------------------------
2785 subroutine hgctex
2786c----------------------------------------------------------------------
2787c returns new tem using massive quantum statistics in integral form
2788c input:
2789c chemgc
2790c tecm/volu
2791c output:
2792c tem
2793c----------------------------------------------------------------------
2794 include 'epos.inc'
60ed90b3 2795 parameter (mspecs=56)
9ef1c2d9 2796 common/cspecs/nspecs,ispecs(mspecs),aspecs(mspecs),gspecs(mspecs)
2797 common/cgchg/rmsngc(mspecs),ptlngc(mspecs),chemgc(mspecs),tem
2798 common/ciakt/gen,iafs,ians,genm
2799 external hgcfhe
2800 k=1
2801 t1=0.0
2802 t2=tem+0.1
2803 goto15
2804
2805c new temperature
2806c ---------------
280710 tem=t1+.5*(t2-t1)
280815 continue
2809 if(tem.le.1.e-6)return
2810 eden=0.0
60ed90b3 2811
9ef1c2d9 2812 do ians=1,nspecs
2813 call hgclim(a,b)
2814 if(b.eq.0.0)then
2815 edi=0.0
2816 else
2817 call uttraq(hgcfhe,a,b,edi)
2818 endif
2819 edi=edi*gspecs(ians)/2./pi**2/hquer**3
2820 eden=eden+edi
2821 enddo
2822
2823 if(iabs(ispecs(nspecs)).lt.10)
2824 *eden=eden+(8.*pi**2*tem**4/15.+bag4rt**4)/hquer**3
2825
2826 de=abs(eden-(tecm/volu))
2827 if(de.le.gen*(tecm/volu).or.de.le.genm)return
2828
2829 if(eden.gt.(tecm/volu))then
2830 t2=tem
2831 else
2832 t1=tem
2833 endif
2834
2835 if(k.gt.300)then
2836 if(ish.ge.5)
2837 *write(ifch,*)'failure in tex'
2838 return
2839 endif
60ed90b3 2840
9ef1c2d9 2841 k=k+1
2842 goto10
2843 end
2844
2845c-----------------------------------------------------------------
2846 subroutine hgctm0
2847c-----------------------------------------------------------------
2848c returns new tem using massless quantum statistics in analytic form
2849c input:
2850c chemgc
2851c tecm/volu
2852c output:
2853c tem
2854c----------------------------------------------------------------------
2855
2856 include 'epos.inc'
60ed90b3 2857 parameter (mspecs=56)
9ef1c2d9 2858 common/cspecs/nspecs,ispecs(mspecs),aspecs(mspecs),gspecs(mspecs)
2859 common/cgchg/rmsngc(mspecs),ptlngc(mspecs),chemgc(mspecs),tem
2860 common/ciakt/gen,iafs,ians,genm
2861
2862 k=1
60ed90b3 2863
9ef1c2d9 2864 t1=0.0
2865 t2=1.0
286610 tem=t1+.5*(t2-t1)
2867 if(tem.le.1.e-6)return
2868 eden=0.0
2869
2870 do i=1,nspecs
2871
2872 igsp=int(gspecs(i))
2873 if(mod(igsp,2).eq.0)then
2874 edhm0=7./240.*pi**2*tem**4+chemgc(i)**2*tem**2/8.
2875 *+chemgc(i)**4/pi**2/16.
2876 else
2877 edhm0=pi**2*tem**4/30.+chemgc(i)**2*tem**2/4.
2878 *-chemgc(i)**4/pi**2/16.
2879 endif
2880 edi=edhm0*gspecs(i)/hquer**3
2881
2882
2883 eden=eden+edi
2884 enddo
2885
2886 if(iabs(ispecs(nspecs)).lt.10)
2887 *eden=eden+(8.*pi**2*tem**4/15.+bag4rt**4)/hquer**3
2888
2889 de=abs(eden-(tecm/volu))
2890 if(de.le.gen*(tecm/volu).or.de.le.genm)return
2891
2892 if(eden.gt.(tecm/volu))then
2893 t2=tem
2894 else
2895 t1=tem
2896 endif
2897
2898 if(k.gt.300)then
2899 if(ish.ge.5)
2900 *write(ifch,*)'failure in tm0'
2901 return
2902 endif
60ed90b3 2903
9ef1c2d9 2904 k=k+1
2905 goto10
2906 end
2907
2908c----------------------------------------------------------------------
2909 subroutine hnbxxx(ip,iret)
2910c----------------------------------------------------------------------
2911c decays droplet very fast ... and hopefully not too badly
2912c----------------------------------------------------------------------
2913 include 'epos.inc'
2914 integer jc(nflav,2),jc1(nflav,2)
2915 integer ifl(nflav)
2916 double precision p(5),c(5)
2917 real u(3),pin(4),pout(4),poutx(4)
60ed90b3 2918 parameter (mspecs=56)
9ef1c2d9 2919 common/cspecs/nspecs,ispecs(mspecs),aspecs(mspecs),gspecs(mspecs)
2920 !-------------------------------------------------------------
2921 ! 110, 120, -120, 130, -130, 230, -230, 220, 330
2922 !, 111, 121, -121, 131, -131, 231, -231, 221, 331
2923 !, 1120,-1120, 1220,-1220, 1130,-1130, 2130,-2130
2924 !, 1230,-1230, 2230,-2230, 1330,-1330, 2330,-2330
2925 !, 1111,-1111, 1121,-1121, 1221,-1221, 2221,-2221, 1131,-1131
60ed90b3 2926 !, 1231,-1231, 2231,-2231, 1331,-1331, 2331,-2331, 3331,-3331
9ef1c2d9 2927 !-------------------------------------------------------------
2928 integer ianti(mspecs)
2929 data ianti/ 1, 3, 2, 5, 4, 7, 6, 8, 9
2930 * ,10,12,11,14,13,16,15,17,18
2931 * ,20,19,22,21,24,23,26,25
2932 * ,28,27,30,29,32,31,34,33
2933 * ,36,35,38,37,40,39,42,41,44,43
2934 * ,46,45,48,47,50,49,52,51,54,53,0,0/
2935 common/xxxspecs/wtot,wspecs(mspecs),zspecs(mspecs)
2936 parameter(mxdrop=35,mxe=10)
2937 common/xxxspecsy/ndrop(-4:4,-4:4,-4:4)
2938 common/xxxspecsx/ee(mxe),wwspecs(mxdrop,mxe,mspecs)
2939 common/cspec2/jspecs(2,nflav,mspecs)
2940 common/cspec4/lkfoi(8,-3:3,-3:3,-3:3,-3:3) !-charm
2941ctp060829 parameter (mxidh=3331)
2942 real wbar(-1:1)
2943
2944 call utpri('hnbxxx',ish,ishini,4)
2945 iret=0
2946
2947 call idquac(ip,nqi,nsi,nai,jc)
2948 keu=jc(1,1)-jc(1,2)
2949 ked=jc(2,1)-jc(2,2)
2950 kes=jc(3,1)-jc(3,2)
2951 n=ndrop(keu,ked,kes)
2952 if(n.eq.0)stop'hnbxxx: n=0'
60ed90b3 2953
9ef1c2d9 2954c...fill wspecs
2955
2956 e=pptl(5,ip)
2957 k=1
2958 do while(k.lt.mxe.and.ee(k).lt.e)
2959 k=k+1
2960 enddo
2961 k=max(k,2)
2962 xi=(e-ee(k-1))/(ee(k)-ee(k-1))
2963 do i=1,nspecs
2964 ii=i
2965 if(n.lt.0)ii=ianti(i)
2966 w1=wwspecs(abs(n),k-1,ii)
2967 w2=wwspecs(abs(n),k,ii)
60ed90b3 2968 wspecs(i)=w1+xi*(w2-w1)
9ef1c2d9 2969 enddo
60ed90b3 2970
9ef1c2d9 2971 if(ish.ge.5)then
2972 write(ifch,*)'keu,ked,kes,n:',keu,ked,kes,n
2973 write(ifch,'(9x, 9f6.3)')(wspecs(i),i=1,9)
2974 write(ifch,'(9x, 9f6.3)')(wspecs(i),i=10,18)
2975 write(ifch,'(9x, 8f6.3)')(wspecs(i),i=19,26)
2976 write(ifch,'(9x, 8f6.3)')(wspecs(i),i=27,34)
2977 write(ifch,'(9x,10f6.3)')(wspecs(i),i=35,44)
2978 write(ifch,'(9x,10f6.3)')(wspecs(i),i=45,54)
2979 endif
2980
2981c...initializations
2982
2983 do j=1,5
2984 c(j)=pptl(j,ip)
2985 enddo
2986 do j=1,4
2987 pin(j)=pptl(j,ip)
2988 pout(j)=0
2989 poutx(j)=0
2990 enddo
2991 wtot=0
2992 do i=1,nspecs
2993 wtot=wtot+wspecs(i)
2994 zspecs(i)=0
2995 enddo
2996 wbar(1)=0
2997 wbar(-1)=0
2998 do i=19,53,2
2999 wbar(1)=wbar(1)+wspecs(i)
3000 wbar(-1)=wbar(-1)+wspecs(i+1)
3001 enddo
3002 wbar(1)=wbar(1) /wtot
3003 wbar(-1)=wbar(-1)/wtot
3004 w12=0
3005 do i=19,34
3006 w12=w12+wspecs(i)
3007 enddo
3008 w32=0 !35,36,41,42,53,54 excluded
60ed90b3 3009 do i=37,40
9ef1c2d9 3010 w32=w32+wspecs(i)
3011 enddo
60ed90b3 3012 do i=43,52
9ef1c2d9 3013 w32=w32+wspecs(i)
3014 enddo
3015 sum=w12+w32
3016 w12=w12/sum
3017 w32=w32/sum
3018 w0=0
3019 do i=1,9
3020 w0=w0+wspecs(i)
3021 enddo
3022 w1=0
3023 do i=10,18
3024 w1=w1+wspecs(i)
3025 enddo
3026 sum=w0+w1
3027 w0=w0/sum
3028 w1=w1/sum
3029 nptlb=nptl
3030 call idquac(ip,nq,ns,na,jc)
3031 nbar=nq/3
3032 nbarini=nbar
3033 do nf=1,nflav
3034 ifl(nf)=jc(nf,1)-jc(nf,2)
3035 enddo
60ed90b3 3036
3037c...print
3038
9ef1c2d9 3039 if(ish.ge.5)then
3040 write(ifch,*)'ip=',ip,' id=',idptl(ip),' e=',sngl(c(5))
3041 write(ifch,*)'jc=',jc
3042 write(ifch,*)'nq=',nq,' na=',na,' nbar=',nbar
3043 write(ifch,*)'wtot=',wtot
3044 endif
3045
3046c...generate number of hadrons
60ed90b3 3047
9ef1c2d9 3048 wfac=1.05 !mean increased by factor wfac
3049 aa=wtot*wfac
3050 if(aa.le.70.)then
3051 776 nhad=0
3052 pr=1
3053 sum=1
3054 wfac=1.05 !mean increased by factor wfac
3055 rr=rangen()*exp(aa)
3056 do while (sum.lt.rr)
3057 nhad=nhad+1
3058 if(nhad.gt.10*aa)goto776
3059 pr=pr*aa/nhad
3060 sum=sum+pr
3061 !print*,'r:',rr,' n:',nhad,' sum pr:',sum
3062 enddo
3063 nhad=max(2,nhad)
3064 else
3065 778 nhad=aa-30
3066 sum=0
3067 rr=rangen()
3068 do while (sum.lt.rr)
3069 nhad=nhad+1
3070 if(nhad.gt.10*aa)goto778
3071 sum=sum+exp(nhad-aa)*(aa/nhad)**nhad/sqrt(2*pi*nhad)
3072 !print*,'r:',rr,' n:',nhad,' sum pr:',sum
3073 enddo
3074 nhad=max(2,nhad)
3075 endif
3076 if(ish.ge.5)write(ifch,*)'-----> ',nhad,' hadrons'
60ed90b3 3077
3078c...generate first n-2 hadrons
9ef1c2d9 3079
3080 do n=1,nhad-2
3081 1 sum=0
3082 i=0
3083 rr=rangen()*wtot
3084 do while (sum.lt.rr)
3085 i=i+1
3086 sum=sum+wspecs(i)
3087 enddo
3088 if(ispecs(i).gt.1000)then
3089 nbari=1
3090 elseif(ispecs(i).lt.-1000)then
3091 nbari=-1
3092 else
60ed90b3 3093 nbari=0
9ef1c2d9 3094 endif
3095 if(nbari*nbarini.gt.0
3096 * .and.nbari*(nbar-nbari).lt.0.and.rangen().gt.wbar(-nbari))then
3097 if(ish.ge.5)
3098 * write(ifch,*)'-----',nptl,nbar,ispecs(i),wbar(-nbari),rr
3099 goto1
3100 elseif(nbari*(nbar-nbari).lt.0)then
3101 if(ish.ge.5)write(ifch,*)'+++++',wbar(-nbari)
60ed90b3 3102 endif
3103 nbar=nbar-nbari
9ef1c2d9 3104 nptl=nptl+1
3105 id=ispecs(i)
60ed90b3 3106 idptl(nptl)=id
9ef1c2d9 3107 call idquac(nptl,nq,ns,na,jc1)
3108 do nf=1,nflav
3109 ifl(nf)=ifl(nf)-jc1(nf,1)+jc1(nf,2)
3110 enddo
3111 if(ish.ge.5)
60ed90b3 3112 * write(ifch,*)'nptl=',nptl,' id=',id,' ifl=',ifl
9ef1c2d9 3113 enddo
3114 do nf=1,nflav
3115 if(ifl(nf).ge.0)then
3116 jc(nf,1)=ifl(nf)
3117 jc(nf,2)=0
3118 else
3119 jc(nf,1)=0
3120 jc(nf,2)=-ifl(nf)
60ed90b3 3121 endif
9ef1c2d9 3122 enddo
3123 if(ish.ge.5)then
3124 write(ifch,*)'jc=',jc
3125 write(ifch,*)'hadrons:',(idptl(n),n=nptlb+1,nptl)
3126 * ,' --> nbar=',nbar
3127 endif
60ed90b3 3128
9ef1c2d9 3129c...last two hadrons
60ed90b3 3130
9ef1c2d9 3131 if(nbar.ne.0)then
3132 do n=1,abs(nbar)
60ed90b3 3133 ii=nbar/abs(nbar)
9ef1c2d9 3134 i1=idraflz(jc,(3-ii)/2)
3135 i2=idraflz(jc,(3-ii)/2)
3136 i3=idraflz(jc,(3-ii)/2)
3137 if(i1.eq.i2.and.i2.eq.i3)then
3138 id=ii*(i1*1000+i2*100+i3*10+1)
60ed90b3 3139 else
9ef1c2d9 3140 if(i2.lt.i1)then
60ed90b3 3141 ix=i1
9ef1c2d9 3142 i1=i2
3143 i2=ix
60ed90b3 3144 endif
9ef1c2d9 3145 if(i3.lt.i2)then
60ed90b3 3146 ix=i2
9ef1c2d9 3147 i2=i3
3148 i3=ix
60ed90b3 3149 endif
9ef1c2d9 3150 if(i2.lt.i1)then
60ed90b3 3151 ix=i1
9ef1c2d9 3152 i1=i2
3153 i2=ix
60ed90b3 3154 endif
9ef1c2d9 3155 ispin=0
3156 if(rangen().lt.w32)ispin=1
3157 id=ii*(i1*1000+i2*100+i3*10+ispin)
60ed90b3 3158 endif
9ef1c2d9 3159 nptl=nptl+1
60ed90b3 3160 idptl(nptl)=id
9ef1c2d9 3161 if(ish.ge.5)
60ed90b3 3162 * write(ifch,*)'nptl=',nptl,' baryon=',id,' jc=',jc
9ef1c2d9 3163 enddo
3164 endif
3165
3166 call idquacjc(jc,nqu,naq)
3167 do while (nqu.gt.0.or.nptl.eq.nptlb)
3168 if(nqu.eq.0)then
3169 i1=1.5+rangen()
3170 i2=1.5+rangen()
3171 jc(i1,1)=jc(i1,1)+1
3172 jc(i2,2)=jc(i2,2)+1
3173 else
3174 i1=idraflz(jc,1)
3175 i2=idraflz(jc,2)
3176 endif
3177 ii=1
3178 if(i2.lt.i1)then
60ed90b3 3179 ix=i1
9ef1c2d9 3180 i1=i2
3181 i2=ix
3182 ii=-1
3183 endif
3184 ispin=0
3185 if(rangen().lt.w1)ispin=1
3186 id=ii*(i1*100+i2*10+ispin)
3187 nptl=nptl+1
60ed90b3 3188 idptl(nptl)=id
9ef1c2d9 3189 if(ish.ge.5)write(ifch,*)'nptl=',nptl,' nqu=',nqu
60ed90b3 3190 & ,' naq=',naq,' --> meson',id
9ef1c2d9 3191 call idquacjc(jc,nqu,naq)
3192 enddo
3193
3194 nmiss=nhad-nptl+nptlb
3195 if(ish.ge.5)then
3196 write(ifch,*)nmiss,' hadron(s) missing'
3197 write(ifch,*)'hadrons:',(idptl(n),n=nptlb+1,nptl)
3198 endif
60ed90b3 3199
3200c nsechad=lkfoi(1,ifl(1),ifl(2),ifl(3),ifl(4))
3201c if(nsechad.gt.0)then
9ef1c2d9 3202c i2x=min(nsechad,1+rangen()*nsechad)
3203c i2=lkfoi(1+i2x,ifl(1),ifl(2),ifl(3),ifl(4))
3204c !print*,'secnd chosen hadron:',ispecs(i2),wzspecs(i2)
60ed90b3 3205
3206c... generate momenta
9ef1c2d9 3207
3208 if(ish.ge.5)write(ifch,*)'hadron momenta:'
3209 do n=nptlb+1,nptl
3210 id=idptl(n)
3211 call idmass(id,am)
3212 !prepare proposal function
3213 ! f_prop: f0(x)=0, x<am
60ed90b3 3214 ! f1(x)=const=f2(b), am<x<b,
9ef1c2d9 3215 ! f2(x)=x**2*exp(7-a*x), x>b
3216 a=11
3217 b=0.9
3218 b=max(b,am+0.1)
3219 !relative weights ... consider here f1/exp(7), f2/exp(7)
3220 c1=(b-am)*b**2*exp(-a*b)
3221 c2=(b**2/a+2*b/a**2+2/a**3)*exp(-a*b)
3222 if(ish.ge.5)write(ifch,*)'c1 c2:',c1,c2
3223 rrr=1
3224 fff=0
3225 do while(rrr.gt.fff)
3226 !proposal
3227 i=2
3228 if(rangen().lt.c1/(c1+c2))i=1
3229 if(ish.ge.5)write(ifch,*)'i=',i
3230 if(i.eq.1)then
3231 x=am+rangen()*(b-am)
3232 fx=b**2*exp(-a*b) *exp(7.)
3233 else
3234 r=alog(rangen())
3235 !find root of log((x**2/a+2*x/a**2+2/a**3)*exp(-a*x)/c2)-r
3236 x=20
3237 do k=1,4
3238 f=alog(x**2/a+2*x/a**2+2/a**3)-a*x-alog(c2)-r
3239 fp=(2*x/a+2/a**2)/(x**2/a+2*x/a**2+2/a**3)-a
3240 xb=x
3241 x=x-f/fp
3242 if(ish.ge.5)write(ifch,*)'k,xb,f,fp,x: ',k,xb,f,fp,x
3243 enddo
3244 fx=x**2*exp(-a*x) *exp(7.)
3245 endif
3246 !acceptance
3247 rrr=rangen()
3248 if(x.le.2)then
3249 fff=sqrt(x-am)*x*exp(-4*x-1.75*x**2) /fx
3250 else
3251 fff=sqrt(x-am)*x*exp(7-11*x) /fx
3252 endif
3253 if(ish.ge.5)write(ifch,*)'fff,rrr',fff,rrr
3254 enddo
3255 e=x
3256 pa=sqrt(e**2-am**2)
3257 u(3)=2.*rangen()-1.
3258 phi=2.*pi*rangen()
3259 u(1)=sqrt(1.-u(3)**2)*cos(phi)
3260 u(2)=sqrt(1.-u(3)**2)*sin(phi)
3261 do j=1,3
3262 pptl(j,n)=pa*u(j)
3263 pout(j)=pout(j)+pptl(j,n)
3264 enddo
3265 pptl(4,n)=e
3266 pptl(5,n)=am
3267 pout(4)=pout(4)+e
3268 if(ish.ge.5)
3269 * write(ifch,*)pptl(1,n),pptl(2,n),pptl(3,n),pptl(4,n),pptl(5,n)
3270 enddo
3271
3272c...check total energy ... rescale (maybe)
3273
3274 if(ish.ge.5)then
3275 write(ifch,*)'pout(cms)=',pout
3276 write(ifch,*)'energy_in(cms)= ',sngl(c(5))
3277 endif
3278
3279c...boost
3280
3281 do n=nptlb+1,nptl
3282 do j=1,4
3283 p(j)=pptl(j,n)
3284 enddo
3285 call utlob2(-1,c(1),c(2),c(3),c(4),c(5),p(1),p(2),p(3),p(4),10)
3286 do j=1,4
3287 pptl(j,n)=p(j)
3288 poutx(j)=poutx(j)+pptl(j,n)
3289 enddo
3290 enddo
3291 if(ish.ge.5)then
3292 write(ifch,*)'pout(lab)=',poutx
3293 write(ifch,*)'pin(lab)= ',pin
3294 endif
3295 call utprix('hnbxxx',ish,ishini,4)
3296 return
3297 end
60ed90b3 3298
9ef1c2d9 3299c----------------------------------------------------------------------
3300 subroutine hnbxxxini
3301c----------------------------------------------------------------------
3302 include 'epos.inc'
3303 logical lcalc
60ed90b3 3304 parameter (mspecs=56)
9ef1c2d9 3305 common/cspecs/nspecs,ispecs(mspecs),aspecs(mspecs),gspecs(mspecs)
3306 !-------------------------------------------------------------
3307 ! 110, 120, -120, 130, -130, 230, -230, 220, 330
3308 !, 111, 121, -121, 131, -131, 231, -231, 221, 331
3309 !, 1120,-1120, 1220,-1220, 1130,-1130, 2130,-2130
3310 !, 1230,-1230, 2230,-2230, 1330,-1330, 2330,-2330
3311 !, 1111,-1111, 1121,-1121, 1221,-1221, 2221,-2221, 1131,-1131
60ed90b3 3312 !, 1231,-1231, 2231,-2231, 1331,-1331, 2331,-2331, 3331,-3331
9ef1c2d9 3313 !-------------------------------------------------------------
3314 common/xxxspecs/wtot,wspecs(mspecs),zspecs(mspecs)
60ed90b3 3315 integer ittspecs(mspecs)
9ef1c2d9 3316 parameter(mxdrop=35,mxe=10)
3317 common/xxxspecsy/ndrop(-4:4,-4:4,-4:4)
3318 common/xxxspecsx/ee(mxe),wwspecs(mxdrop,mxe,mspecs)
3319 inquire(file=fndr(1:nfndr),exist=lcalc)
3320 if(lcalc)then
3321 if(inicnt.eq.1)then
3322 do ku=-4,4
3323 do kd=-4,4
3324 do ks=-4,4
3325 ndrop(ku,kd,ks)=0
3326 enddo
3327 enddo
3328 enddo
3329 write(ifmt,'(3a)')'read from ',fndr(1:nfndr),' ...'
3330 open(1,file=fndr(1:nfndr),status='old')
3331 read(1,*)mxxdrop
60ed90b3 3332 if(mxxdrop.ne.mxdrop)stop'hnbxxxini: wrong nr of droplets'
9ef1c2d9 3333 do n=1,mxdrop
3334 read(1,*)ku,kd,ks
3335 if(abs(ku).gt.4)stop'hnbxxxini: ku out of range'
3336 if(abs(kd).gt.4)stop'hnbxxxini: kd out of range'
3337 if(abs(ks).gt.4)stop'hnbxxxini: ks out of range'
3338 ndrop(ku,kd,ks)=n
3339 ndrop(-ku,-kd,-ks)=-n
3340 enddo
3341 read(1,*)(ittspecs(i),i=1,nspecs)
3342 do i=1,nspecs
3343 if(ittspecs(i).ne.ispecs(i))stop'hnbxxxini: wrong id table'
3344 enddo
3345 read(1,*)ee
3346 do n=1,mxdrop
3347 read(1,*)((wwspecs(n,k,i),i=1,nspecs),k=1,mxe)
3348 enddo
3349 close(1)
3350 endif
3351 else
3352 stop'hnbxxxini: file not found. '
3353 endif
60ed90b3 3354 end
9ef1c2d9 3355
3356c----------------------------------------------------------------------
3357 subroutine hnbaaa(ip,iret)
3358c----------------------------------------------------------------------
60ed90b3 3359 include 'epos.inc'
3360 if(ioclude.eq.1)call hnbaaa156(ip,iret)
9ef1c2d9 3361 if(ioclude.eq.2)stop'ioclude.eq.2 no longer supported. '
3362 if(ioclude.eq.3)call hnbaaanew(ip,iret)
3363 end
60ed90b3 3364
9ef1c2d9 3365c----------------------------------------------------------------------
3366 subroutine hnbaaanew(ip,iret)
3367c----------------------------------------------------------------------
3368c microcanonical decay of cluster ip via loop over hnbmet
3369c----------------------------------------------------------------------
3370 include 'epos.inc'
3371 include 'epos.inchy'
3372 common/cxyzt/xptl(mxptl),yptl(mxptl),zptl(mxptl),tptl(mxptl)
3373 *,optl(mxptl),uptl(mxptl),sptl(mxptl),rptl(mxptl,3)
3374 parameter(maxp=500)
3375 common/confg/np,amass(maxp),ident(maxp),pcm(5,maxp),wtxlog,wtlog
3376 common/citer/iter,itermx
3377 integer jc(nflav,2)
3378 double precision p(5),c(5)
3379 parameter(maxit=50000)
3380 common/count/nacc,nrej,naccit(maxit),nptot,npit(maxit)
3381 dimension uu(4),pe(5),pa(5)
3382 common/yrad/yrad(maxp),phifop(maxp),radfop(maxp),taufop(maxp)
3383 common/xxxspecsy/ndrop(-4:4,-4:4,-4:4)
3384 common/cdelzet/delzet,delsgr /cvocell/vocell
3385 common/cen/ncentr /ctauhac/ntauhac(ncenthy,netahy)
3386 common/cranphi/ranphi,ranecc,weiecc
3387 character txt*40
60ed90b3 3388 data icnthnb /0/ !vv2 /0./ nvv2 /0/ vv3 /0./
9ef1c2d9 3389 !save vv2,nvv2,vv3
3390 save icnthnb
60ed90b3 3391
9ef1c2d9 3392 call utpri('hnbaaa',ish,ishini,4)
60ed90b3 3393
9ef1c2d9 3394 if(ish.ge.3)then
3395 write(ifch,140)
3396 140 format(/' ----------------------------------'/
3397 *' droplet decay'/
3398 *' ----------------------------------')
3399 write(ifch,*)'droplet:'
3400 call alist('&',ip,ip)
3401 endif
60ed90b3 3402
9ef1c2d9 3403 iret=0
3404 do j=1,5
3405 c(j)=pptl(j,ip)
3406 enddo
60ed90b3 3407
9ef1c2d9 3408 call idquac(ip,nqi,nsi,nai,jc)
3409 keu=jc(1,1)-jc(1,2)
3410 ked=jc(2,1)-jc(2,2)
3411 kes=jc(3,1)-jc(3,2)
3412 kec=jc(4,1)-jc(4,2)
3413 keb=jc(5,1)-jc(5,2)
3414 ket=jc(6,1)-jc(6,2)
3415 !print*,'droplet uds=',keu,ked,kes,' E=',pptl(5,ip)
3416
3417 volu=4./3.*pi*radptl(ip)**3
3418
3419 if(volu.le.0)call utstop('hnbaaa: volume = 0&')
3420 if(volu.lt.0.01)then
3421 call utmsg('hnbaaa')
3422 write(ifch,*)'***** very small volume:',volu
3423 write(ifch,*)
3424 * 'id:',idptl(ip),' r:',radptl(ip),' m:',pptl(5,ip)
3425 call utmsgf
3426 endif
60ed90b3 3427
3428 !~~~~~~~~~read in freeze out surface properties from hydro~~~~~~~~~~~~
3429 if(iorsdf.eq.3.and.ityptl(ip).eq.60)then
3430 icnthnb=icnthnb+1
9ef1c2d9 3431 if(icnthnb.eq.1)then
60ed90b3 3432 !here we use epos.iniXXX (like epos.ini1fc)
9ef1c2d9 3433 !instead of the generic filename epos.inihy
3434 open(3,file=fnnx(1:nfnnx)//fnhy(1:nfnhy),status='old',err=99)
3435 read(3,*)txt
3436 read(3,*)maprojx,matargx,engyx,epscrix
3437 if(maprojx.ne.maproj)stop'hnbaaa: maprojx.ne.maproj. '
3438 if(matargx.ne.matarg)stop'hnbaaanew: matargx.ne.matarg. '
3439 if(engyx.ne.engy)stop'hnbaaanew: engyx.ne.engy. '
3440 if(epscrix.ne.epscri(ioclude))then
3441 print*,'hnbaaanew: epscrix.ne.epscri(',ioclude,'). '
3442 stop
60ed90b3 3443 endif
9ef1c2d9 3444 read(3,*)ncenthyx,netahyx,ntauhyx,nphihyx,nradhyx
3445 if(ncenthyx.ne.ncenthy)stop'hnbaaanew: ncenthyx.ne.ncenthy. '
3446 if(netahyx.ne.netahy)stop'hnbaaanew: netahyx.ne.netahy. '
3447 if(ntauhyx.ne.ntauhy)stop'hnbaaanew: ntauhyx.ne.ntauhy. '
3448 if(nphihyx.ne.nphihy)stop'hnbaaanew: nphihyx.ne.nphihy. '
3449 if(nradhyx.ne.nradhy)stop'hnbaaanew: nradhyx.ne.nradhy. '
3450 read(3,*)ntauhoc
3451 read(3,*)centhy,etahy,phihy,radhy
3452 read(3,*)tauhoc
3453 read(3,*)epsii
3454 read(3,*)rom,yom,wom
3455 close(3)
3456 endif
3457 endif
3458
60ed90b3 3459 !~~~~~~~~~define womi yomi romi~~~~~~~~~~~~
3460 if(iorsdf.eq.3.and.icnthnb.eq.1)then
9ef1c2d9 3461 if(ioclude.eq.3)then
3462 do ncent=1,ncenthy
3463 do neta=1,netahy
3464 do ntau=1,ntauhoc(ncent)
3465 do i=1,2
3466 womi(ncent,neta,ntau,i)=wom(ncent,neta,ntau,i)
3467 yomi(ncent,neta,ntau,i)=yom(ncent,neta,ntau,i)
3468 romi(ncent,neta,ntau,i)=rom(ncent,neta,ntau,i)
3469 enddo
3470 enddo
3471 enddo
3472 enddo
3473 elseif(ioclude.eq.2)then
3474 stop'in hnbaaanew: ioclude=2 not supported any more.'
3475 else
60ed90b3 3476 stop'in hnbaaa: invalid ioclude. '
9ef1c2d9 3477 endif
3478 endif
3479
3480 !~~tau partition function paut(ncent,neta,ntau)
3481 !~~phi partition function pauf(ncent,neta,ntau,nphi)
60ed90b3 3482 if(iorsdf.eq.3.and.icnthnb.eq.1)then
9ef1c2d9 3483 do ncent=1,ncenthy
3484 do neta=1,netahy
3485 womax=0
3486 ntauhac(ncent,neta)=0
60ed90b3 3487 do ntau=1,ntauhoc(ncent)
3488 if(womi(ncent,neta,ntau,1).gt.womax )then
9ef1c2d9 3489 womax=womi(ncent,neta,ntau,1)
3490 ntauhac(ncent,neta)=ntau
3491 endif
3492 !print*,ncent,neta,ntau,ntauhac(ncent,neta)
3493 !. , womi(ncent,neta,ntau,1)
3494 enddo
3495 paut(ncent,neta,1)=0
3496 do ntau=2,ntauhac(ncent,neta)
3497 paut(ncent,neta,ntau)=womi(ncent,neta,ntau,1)/womax
3498 !print*,ncent,neta,ntau,paut(ncent,neta,ntau)
3499 enddo
3500 !~~phi~~~~
3501 do ntau=2,ntauhac(ncent,neta)
3502 pauf(ncent,neta,ntau,1)=0
3503 do nphi=2,nphihy
3504 dphi=phihy(nphi)-phihy(nphi-1)
3505 c1=cos(2*phihy(nphi-1))
3506 c2=cos(2*phihy(nphi))
3507 w1=womi(ncent,neta,ntau ,1)+c1*womi(ncent,neta,ntau ,2)
60ed90b3 3508 . -womi(ncent,neta,ntau-1,1)-c1*womi(ncent,neta,ntau-1,2)
9ef1c2d9 3509 w2=womi(ncent,neta,ntau ,1)+c2*womi(ncent,neta,ntau ,2)
60ed90b3 3510 . -womi(ncent,neta,ntau-1,1)-c2*womi(ncent,neta,ntau-1,2)
9ef1c2d9 3511 pauf(ncent,neta,ntau,nphi)
60ed90b3 3512 . =pauf(ncent,neta,ntau,nphi-1)+0.5*(w1+w2)*dphi
9ef1c2d9 3513 enddo
3514 w=pauf(ncent,neta,ntau,nphihy)
3515 if(w.eq.0.)stop'hnbaaanew: w.eq.0. '
3516 do nphi=2,nphihy
3517 pauf(ncent,neta,ntau,nphi)=pauf(ncent,neta,ntau,nphi)/w
3518 enddo
3519 enddo
3520 enddo
3521 enddo
3522 endif
60ed90b3 3523
9ef1c2d9 3524 !~~~~~~~~~determine ncentr~~~~~~~~~~~~
3525 if(iorsdf.eq.3.and.ityptl(ip).eq.60)then !!!fusion on!!!
3526 ncentr=0
3527 dbmin=1000.
3528 do ncent=1,ncenthy
3529 db=abs(bimevt-centhy(ncent))
3530 if(db.lt.dbmin)then
3531 dbmin=db
3532 ncentr=ncent
60ed90b3 3533 endif
9ef1c2d9 3534 enddo
3535 !print*,ncentr,bimevt,centhy(ncentr)
3536 endif
3537
60ed90b3 3538 !~~~~~define masses~~~~~~~~~~~~~~~~
9ef1c2d9 3539 amin=utamnu(keu,ked,kes,kec,keb,ket,5)
3540 aumin=amuseg
3541 ipo=ip
3542 if(ityptl(ip).eq.60)ipo=iorptl(ip)
3543 tecmor=pptl(5,ipo)
3544 tecm=pptl(5,ip)
3545 tecmxx=tecm
3546
3547 !~~~~~~~~~determine netar~~~~~~~~~~~~
60ed90b3 3548 if(iorsdf.eq.3.and.ityptl(ip).eq.60)then
9ef1c2d9 3549 z=xorptl(3,ipo)
3550 t=xorptl(4,ipo)
3551 !print*,z,t,ityptl(ip),tecm
3552 if(t+z.le.0..or.t-z.le.0.)then
3553 zetaor=1000.
3554 else
3555 zetaor=abs(0.5*log((t+z)/(t-z)))
3556 endif
3557 netar=0
3558 detamin=1000.
3559 do neta=1,netahy
3560 deta=abs(zetaor-etahy(neta))
3561 if(deta.lt.detamin)then
3562 detamin=deta
3563 netar=neta
60ed90b3 3564 endif
9ef1c2d9 3565 enddo
3566 !print*,netar,zetaor,etahy(netar)
3567 endif
3568
3569 fradflo=1.
60ed90b3 3570
3571 !~~~~~redefine energy in case of radial flow~~~~~~~~~~~~~~~~
9ef1c2d9 3572 if(iappl.ne.4.and.iorsdf.eq.3.and.tecmor.gt.aumin
3573 ..and.ityptl(ip).eq.60)then
3574 am=0.
3575 do ntau=2,ntauhac(ncentr,netar)
3576 d=paut(ncentr,netar,ntau)-paut(ncentr,netar,ntau-1)
3577 am=am+d
3578 enddo
3579 en=0.
3580 do ntau=2,ntauhac(ncentr,netar)
3581 d=paut(ncentr,netar,ntau)-paut(ncentr,netar,ntau-1)
3582 y=( yomi(ncentr,netar,ntau-1,1)
3583 . +yomi(ncentr,netar,ntau,1)) / 2
3584 en=en+d*cosh(y)
3585 enddo
3586 fradflo=1.
3587 if(en.ne.0.)fradflo=am/en
3588 if(tecm*fradflo.lt.amin)fradflo=1.
60ed90b3 3589 endif
9ef1c2d9 3590 tecm=tecm*fradflo
60ed90b3 3591
3592 !~~~~~redefine energy in case of long coll flow
9ef1c2d9 3593 if(iappl.eq.4.or.iorsdf.ne.3
60ed90b3 3594 &.or.ityptl(ip).eq.40.or.ityptl(ip).eq.50)then !not for droplets from remnants
9ef1c2d9 3595 yco=0
3596 else
3597 if(ylongmx.lt.0.)then
3598 yco=delzet * 1.75
3599 else
3600 yco=ylongmx
3601 endif
60ed90b3 3602 endif
9ef1c2d9