]> git.uio.no Git - u/mrichter/AliRoot.git/blame - EPOS/epos167/epos-yyy-168.f
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
Update master to aliroot
[u/mrichter/AliRoot.git] / EPOS / epos167 / epos-yyy-168.f
CommitLineData
9ef1c2d9 1c----------------------------------------------------------------------
2 subroutine hnbaaa156(ip,iret)
3c----------------------------------------------------------------------
4c microcanonical decay of cluster ip via loop over hnbmet
5c----------------------------------------------------------------------
6 include 'epos.inc'
7 common/cxyzt/xptl(mxptl),yptl(mxptl),zptl(mxptl),tptl(mxptl)
8 *,optl(mxptl),uptl(mxptl),sptl(mxptl),rptl(mxptl,3)
9 parameter(maxp=500)
10 common/confg/np,amass(maxp),ident(maxp),pcm(5,maxp),wtxlog,wtlog
11 common/citer/iter,itermx
12 double precision tpro,zpro,ttar,ztar,ttaus,detap,detat
13 common /cttaus/ tpro,zpro,ttar,ztar,ttaus,detap,detat
14 integer jc(nflav,2)
15 double precision p(5),c(5)
16 parameter(maxit=50000)
17 common/count/nacc,nrej,naccit(maxit),nptot,npit(maxit)
18 dimension be(4),pe(5),pa(5)
19 common/yradx/yrad(maxp),phirad(maxp)
20 common/xxxspecsy/ndrop(-4:4,-4:4,-4:4)
21 common/cdelzet/delzet,delsgr /cvocell/vocell
22 call utpri('hnbaaa',ish,ishini,4)
23
24 if(ish.ge.3)then
25 write(ifch,140)sngl(ttaus)
26 140 format(/' ----------------------------------'/
27 *' droplet decay at tau =',f6.2/
28 *' ----------------------------------')
29 write(ifch,*)'droplet:'
30 call alist('&',ip,ip)
31 endif
32
33 iret=0
34 do j=1,5
35 c(j)=pptl(j,ip)
36 enddo
37
38 call idquac(ip,nqi,nsi,nai,jc)
39 keu=jc(1,1)-jc(1,2)
40 ked=jc(2,1)-jc(2,2)
41 kes=jc(3,1)-jc(3,2)
42 kec=jc(4,1)-jc(4,2)
43 keb=jc(5,1)-jc(5,2)
44 ket=jc(6,1)-jc(6,2)
45 !print*,'droplet uds=',keu,ked,kes,' E=',pptl(5,ip)
46
47 volu=4./3.*pi*radptl(ip)**3
48
49 if(volu.le.0)call utstop('hnbaaa: volume = 0&')
50 if(volu.lt.0.01)then
51 call utmsg('hnbaaa')
52 write(ifch,*)'***** very small volume:',volu
53 write(ifch,*)
54 * 'id:',idptl(ip),' r:',radptl(ip),' m:',pptl(5,ip)
55 call utmsgf
56 endif
57
58 !~~~~~redefine energy in case of radial flow
59 amin=utamnu(keu,ked,kes,kec,keb,ket,4) !utamnu(...,4) and not utamnu(...,5)
60 aumin=amuseg !otherwise droplet from remnant decay
61 ipo=ip !could be too light after flow
62 if(ityptl(ip).eq.60)ipo=iorptl(ip)
63 tecmor=pptl(5,ipo)
64 if(iappl.eq.4.or.iorsdf.ne.3
65 &.or.ityptl(ip).eq.40.or.ityptl(ip).eq.50)then !not for droplets from remnants
66 yrmax=0
67 else
68 yrmax=max(0.0,yradmx+yradmi*alog10(engy/200.))
69 if(maproj.eq.1.and.matarg.eq.1)then
70 yrmax=max(0.0,yradpp+yradpi*alog10(engy/1800.))
71 !aumin=amin
72 !if(yrmax.gt.0.2)print*,'===',tecmor,aamin,yrmax
73 endif
74 endif
75 fradflo=1.
76 if(yrmax.ne.0.)
77 &fradflo=1./((sinh(yrmax)*yrmax-cosh(yrmax)+1.)/(yrmax**2/2.))
78 tecm=pptl(5,ip)
79 tecmxx=tecm
80 if(yrmax.gt.0..and.tecmor.gt.aumin
81 & .and.tecm*fradflo.gt.amin) then
82 ! redefine energy to account for collective flow
83 ! \int_0^yrmax f(y) d2y = E_new (effective mass)
84 ! \int_0^yrmax cosh(y) f(y) d2y = E_old
85 tecm=tecm*fradflo
86 if(tecm.lt.amin)stop'aaahnb: small mass. should not happen. '
87 else
88 yrmax=0.
89 endif
90
91 !~~~~~redefine energy in case of long coll flow
92 if(iappl.eq.4.or.iorsdf.ne.3
93 &.or.ityptl(ip).eq.40.or.ityptl(ip).eq.50)then !not for droplets from remnants
94 yco=0
95 else
96 if(ylongmx.lt.0.)then
97 yco=delzet * 1.75
98 else
99 yco=ylongmx
100 endif
101 endif
102 tecmx=tecm
103 if(yco.gt.0..and.tecmor.gt.aumin) then
104 tecm=tecm/sinh(yco)*yco
105 else
106 yco=0.
107 endif
108 !print*,'========= cluster energy: ',pptl(5,ip),tecmx,tecm
109
110 !~~~~~~~~~redefine volume~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
111
112ccc vocri=tecm/epscri(ioclude) !??????????????????????????????
113ccc volu=max(vocri,vocell) !????????????????????????????
114
115 !~~~~~~~~~decay~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
116 call hnbini(iret)
117 !if(iret.ne.0)write(ifch,*)'***** unsucessfull hnbini *****'
118 if(iret.ne.0)goto1000
119 if(ioinct.ge.1)goto1
120
121 do iter=1,itermx
122 naccit(iter)=0
123 call hnbmet
124 enddo
125
1261 continue
127
128 if(ioceau.eq.1.and.iappl.eq.1)call xhnbte(ip)
129
130 !~~~~~~~~~~long coll flow -> particles~~~~~~~~~~~~~~~~
131 if(yco.gt.0.) then
132 errlim=0.0001
133 tecm=tecmx
134 611 energ=0.
135 do i=1,np
136 yrad(i)=(2*rangen()-1)*yco
137 be(3)=sinh(yrad(i))
138 be(4)=cosh(yrad(i))
139 energ=energ+be(4)*pcm(4,i)-be(3)*pcm(3,i)
140 enddo
141 if(abs(energ-tecm).gt.0.1) goto 611
142 !print*,'===== energy after flow boosts',energ,' soll: ',tecm
143 do j=1,4
144 pe(j)=0.
145 enddo
146 do i=1,np
147 be(1)= 0
148 be(2)= 0
149 be(3)= sinh(yrad(i))
150 be(4)= cosh(yrad(i))
151 call utlob3(1,be(1),be(2),be(3),be(4),1e0
152 * , pcm(1,i), pcm(2,i), pcm(3,i), pcm(4,i))
153 do j=1,4
154 pe(j)=pe(j)+pcm(j,i)
155 enddo
156 enddo
157 pe(5)=sqrt(pe(4)**2-pe(3)**2-pe(2)**2-pe(1)**2)
158 !write(6,'(a,5e11.3)')'flow boosts',pe
159 do j=1,4
160 pa(j)=0.
161 enddo
162 do i=1,np
163 call utlob3(1,pe(1),pe(2),pe(3),pe(4),pe(5)
164 * , pcm(1,i), pcm(2,i), pcm(3,i), pcm(4,i))
165 do j=1,4
166 pa(j)=pa(j)+pcm(j,i)
167 enddo
168 enddo
169 pa(5)=sqrt(pa(4)**2-pa(3)**2-pa(2)**2-pa(1)**2)
170 !write(6,'(a,5e11.3)')' cms boost ',pa
171 esoll=tecm
172 scal=1.
173 do ipass=1,200
174 sum=0.
175 do j=1,np
176 do k=1,3
177 pcm(k,j)=scal*pcm(k,j)
178 enddo
179 pcm(4,j)=sqrt(pcm(1,j)**2+pcm(2,j)**2+pcm(3,j)**2
180 * +amass(j)**2)
181 sum=sum+pcm(4,j)
182 enddo
183 scal=esoll/sum
184 !write(6,*)'ipass,scal,e,esoll:'
185 ! $ ,ipass,scal,sum,esoll
186 if(abs(scal-1.).le.errlim) goto301
187 enddo
188 301 continue
189 do j=1,4
190 pa(j)=0.
191 enddo
192 do i=1,np
193 do j=1,4
194 pa(j)=pa(j)+pcm(j,i)
195 enddo
196 enddo
197 pa(5)=sqrt(pa(4)**2-pa(3)**2-pa(2)**2-pa(1)**2)
198 !write(6,'(a,5e11.3)')' rescaling ',pa
199 endif
200
201 !~~~~~~~~~~radial flow -> particles~~~~~~~~~~~~~~~~~~
202 if(yrmax.gt.0.) then
203 fecc=0
204 aa=1
205 bb=0
206 cc=0
207 dd=1
208 if(ityptl(ip).eq.60)then
209 ipo=iorptl(ip)
210 xx=uptl(ipo) ! <x**2>
211 yy=optl(ipo) ! <y**2>
212 xy=desptl(ipo) ! <x*y>
213 dta=0.5*abs(xx-yy)
214 ev1=(xx+yy)/2+sqrt(dta**2+xy**2)
215 ev2=(xx+yy)/2-sqrt(dta**2+xy**2)
216 if(xy.lt.0..and.dta.ne.0.)then
217 theta=0.5*atan(-xy/dta)
218 elseif(xy.gt.0..and.dta.ne.0.)then
219 theta=-0.5*atan(xy/dta)
220 else
221 theta=0
222 endif
223 !eccx=(yy-xx)/(yy+xx)
224 yy=ev1
225 xx=ev2
226 ecc=(yy-xx)/(yy+xx)
227 !print*,eccx,ecc,theta
228 fecc=facecc*ecc
229 fecc=min(0.3,fecc)
230 phiclu=phievt+theta
231 aa=cos(phiclu)
232 bb=sin(phiclu)
233 cc=-sin(phiclu)
234 dd=cos(phiclu)
235 endif
236 errlim=0.0001
237 tecm=tecmxx
238 610 energ=0.
239 do i=1,np
240 yrad(i)=sqrt(rangen())
241 phirad(i)=2.*pi*rangen()
242 bex=-dsinh(dble(yrad(i)*yrmax))*cos(phirad(i))*(1+fecc)
243 bey=-dsinh(dble(yrad(i)*yrmax))*sin(phirad(i))*(1-fecc)
244 be(1)=aa*bex+cc*bey
245 be(2)=bb*bex+dd*bey
246 be(3)=-0d0
247 be(4)=sqrt(1+be(1)**2+be(2)**2)
248 bp=0d0
249 do k=1,3
250 bp=bp+pcm(k,i)*be(k)
251 enddo
252 en=be(4)*pcm(4,i)+bp
253 energ=energ+en
254 enddo
255 if(abs(energ-tecm).gt.0.1) goto 610
256 energ=0.
257 do i=1,np
258 bex=dsinh(dble(yrad(i)*yrmax))*cos(phirad(i))*(1+fecc)
259 bey=dsinh(dble(yrad(i)*yrmax))*sin(phirad(i))*(1-fecc)
260 be(1)=aa*bex+cc*bey
261 be(2)=bb*bex+dd*bey
262 be(3)=0d0
263 be(4)=sqrt(1+be(1)**2+be(2)**2)
264 call utlob3(1,be(1),be(2),be(3),be(4),1e0
265 * , pcm(1,i), pcm(2,i), pcm(3,i), pcm(4,i))
266 energ=energ+pcm(4,i)
267 enddo
268 esoll=tecm
269 scal=1.
270 do ipass=1,200
271 sum=0.
272 do j=1,np
273 do k=1,3
274 pcm(k,j)=scal*pcm(k,j)
275 enddo
276 pcm(4,j)=sqrt(pcm(1,j)**2+pcm(2,j)**2+pcm(3,j)**2
277 * +amass(j)**2)
278 sum=sum+pcm(4,j)
279 enddo
280 scal=esoll/sum
281 !write(6,*)'ipass,scal,e,esoll:'
282 ! $ ,ipass,scal,sum,esoll
283 if(abs(scal-1.).le.errlim) goto300
284 enddo
285 300 continue
286 else
287 do n=1,np
288 yrad(n)=0.
289 phirad(n)=0.
290 enddo
291 endif
292 !~~~~~~~~~~~~~~~
293
294 nptlb=nptl
295 do n=1,np
296 nptl=nptl+1
297 if(nptl.gt.mxptl)call utstop('hnbptl: mxptl too small&')
298 idptl(nptl)=ident(n)
299 do j=1,4
300 p(j)=pcm(j,n)
301 enddo
302 p(5)=amass(n)
303 call utlob2(-1,c(1),c(2),c(3),c(4),c(5),p(1),p(2),p(3),p(4),10)
304 do j=1,5
305 pptl(j,nptl)=p(j)
306 enddo
307 if(tecmor.gt.aumin)then
308 ityptl(nptl)=60
309 else
310 ityptl(nptl)=19
311 endif
312 if(ityptl(ip).eq.60)then
313 if(ityptl(nptl).eq.60)then
314 ipo=iorptl(ip)
315 xx=uptl(ipo) ! <x**2>
316 yy=optl(ipo) ! <y**2>
317 rini=sqrt(5./3.*(xx+yy)) !<r**2>=3/5*R**2 for sphere of radius R
318 r=1.15*rini*yrad(n) !yrad=y/ymax
319 tau=2.25/sqrt(yrad(n)**2+0.04)-0.75
320 z=xorptl(3,ipo)
321 t=xorptl(4,ipo)
322 !zeta=0.5*log((t+z)/(t-z))-0.5*delzet+2*0.5*delzet*rangen()
323 zeta=0.5*log((p(4)+p(3))/(p(4)-p(3)))
324 z=tau*sinh(zeta)
325 t=tau*cosh(zeta)
326 xorptl(1,nptl)=xorptl(1,ipo)+r*cos(phirad(n))
327 xorptl(2,nptl)=xorptl(2,ipo)+r*sin(phirad(n))
328 xorptl(3,nptl)=z
329 xorptl(4,nptl)=t
330 else
331 xorptl(1,nptl)=xorptl(1,ip)
332 xorptl(2,nptl)=xorptl(2,ip)
333 xorptl(3,nptl)=xorptl(3,ip)
334 xorptl(4,nptl)=xorptl(4,ip)
335 endif
336 endif
337 enddo
338
339 if(ish.ge.3)then
340 write(ifch,*)'decay products:'
341 call alist('&',nptlb+1,nptl)
342 if(ish.ge.5)then
343 write(ifch,*)'momentum sum:'
344 do kk=1,5
345 pptl(kk,nptl+1)=0
346 do ii=nptlb+1,nptl
347 pptl(kk,nptl+1)=pptl(kk,nptl+1)+pptl(kk,ii)
348 enddo
349 pptl(kk,nptl+2)=c(kk)
350 enddo
351 call alist('&',nptl+1,nptl+2)
352 endif
353 endif
354
355 1000 continue
356 call utprix('hnbaaa',ish,ishini,4)
357 return
358 end