]>
Commit | Line | Data |
---|---|---|
8b41a278 | 1 | |
2 | *---------------------------------------------------------------------- | |
3 | * | |
4 | * Filename : PYQUEN.F | |
5 | * | |
6 | * Author : Igor Lokhtin (Igor.Lokhtin@cern.ch) | |
7 | * Version : PYQUEN1_5.f | |
8 | * Last revision : 19-DEC-2007 | |
9 | * | |
10 | *====================================================================== | |
11 | * | |
12 | * Description : Event generator for simulation of parton rescattering | |
13 | * and energy loss in expanding quark-gluon plasma created | |
14 | * in ultrarelativistic heavy ion AA collisons | |
15 | * (implemented as modification of standard Pythia jet event) | |
16 | * | |
17 | * Reference: I.P. Lokhtin, A.M. Snigirev, Eur. Phys. J. C 46 (2006) 211 | |
18 | * | |
19 | *====================================================================== | |
20 | ||
21 | SUBROUTINE PYQUEN(A,ifb,bfix) | |
22 | IMPLICIT DOUBLE PRECISION(A-H, O-Z) | |
23 | IMPLICIT INTEGER(I-N) | |
24 | INTEGER PYK,PYCHGE,PYCOMP | |
25 | external pydata | |
26 | external pyp,pyr,pyk,pyjoin,pyshow | |
27 | external funbip,prhoaa,pfunc1 | |
28 | common /pyjets/ n,npad,k(4000,5),p(4000,5),v(4000,5) | |
29 | common /pydat1/ mstu(200),paru(200),mstj(200),parj(200) | |
30 | common /pysubs/ msel,mselpd,msub(500),kfin(2,-40:40),ckin(200) | |
31 | common /plpar1/ tau0,T0,TC,sigqq,AW,RA,mvisc,nf | |
32 | common /plglur/ glur(1000,4),kglu(1000,6),nrg,nrgm | |
33 | common /plquar/ pqua(1000,5),kqua(1000,5),nrq | |
34 | common /parimp/ b1,psib1,rb1,rb2,noquen | |
35 | common /pyqpar/ T0u,tau0u,nfu,ienglu,ianglu | |
36 | common /plfpar/ bgen | |
37 | common /pygeom/ BC | |
38 | common /pythic/ PBAB(110),PTAB(110),PTAAB(110) | |
39 | common /pynup1/ bp,x | |
40 | save/pyjets/,/pydat1/,/pysubs/,/plglur/,/plquar/,/pygeom/, | |
41 | > /pythic/,/plpar1/,/parimp/,/pyqpar/,/plfpar/ | |
42 | dimension ijoik(2),ijoin(1000),ijoin0(1000),nis(500),nss(500), | |
43 | > nas(500),nus(500) | |
44 | ||
45 | * set initial event paramters | |
46 | AW=A ! atomic weight | |
47 | RA=1.15d0*AW**0.333333d0 ! nucleus radius in fm | |
48 | RA3=3.d0*RA | |
49 | mvisc=0 ! flag of QGP viscosity (off here) | |
50 | TC=0.2d0 ! crutical temperature | |
51 | ||
52 | if(nfu.ne.0.and.nfu.ne.1.and.nfu.ne.2.and.nfu.ne.3) nfu=0 | |
53 | nf=nfu ! number of active flavours in QGP | |
54 | if(tau0u.lt.0.01d0.or.tau0u.gt.10.d0) tau0u=0.1d0 | |
55 | tau0=tau0u ! proper time of QGP formation | |
56 | if(T0u.lt.0.2d0.or.T0u.gt.2.d0) T0u=1.d0 | |
57 | T0=T0u*(AW/207.d0)**0.166667d0 ! initial QGP temperatute at b=0 | |
58 | if(ienglu.ne.0.and.ienglu.ne.1.and.ienglu.ne.2) ienglu=0 ! e-loss type | |
59 | if(ianglu.ne.0.and.ianglu.ne.1.and.ianglu.ne.2) ianglu=0 ! angular spec. | |
60 | * | |
61 | pi=3.14159d0 | |
62 | ||
63 | * avoid stopping run if pythia does not conserve energy due to collisional loss | |
64 | mstu(21)=1 | |
65 | ||
66 | * creation of arrays for tabulation of beam/target nuclear thickness function | |
67 | Z2=4.d0*RA | |
68 | Z1=-1.d0*Z2 | |
69 | H=0.01d0*(Z2-Z1) | |
70 | do ib=1,110 | |
71 | BC=RA3*(ib-1)/109.d0 | |
72 | CALL SIMPA(Z1,Z2,H,0.005d0,1.d-8,prhoaa,Z,RES,AIH,AIABS) | |
73 | PBAB(ib)=BC | |
74 | PTAB(ib)=AW*RES | |
75 | end do | |
76 | ||
77 | * calculation of beam/target nuclear overlap function at b=0 | |
78 | * if ifb=1: creation of arrays for tabulation of nuclear overlap function | |
79 | npb=1 | |
80 | if (ifb.eq.1) npb=110 | |
81 | Z1=0.d0 | |
82 | Z2=6.28318d0 | |
83 | H=0.01d0*(Z2-Z1) | |
84 | do ib=1,npb | |
85 | bp=PBAB(ib) | |
86 | CALL SIMPA(Z1,Z2,H,0.05d0,1.d-8,PFUNC1,X,RES,AIH,AIABS) | |
87 | PTAAB(ib)=RES | |
88 | end do | |
89 | ||
90 | * generate impact parameter of A-A collision with jet production | |
91 | if(ifb.eq.0) then | |
92 | if(bfix.lt.0.d0) then | |
93 | write(6,*) 'Impact parameter less than zero!' | |
94 | bfix=0.d0 | |
95 | end if | |
96 | if (bfix.gt.RA3) then | |
97 | write(6,*) 'Impact parameter larger than three nuclear radius!' | |
98 | bfix=RA3 | |
99 | end if | |
100 | b1=bfix | |
101 | else | |
102 | call bipsear(fmax1,xmin1) | |
103 | fmax=fmax1 | |
104 | xmin=xmin1 | |
105 | 11 bb1=xmin*pyr(0) | |
106 | ff1=fmax*pyr(0) | |
107 | fb=funbip(bb1) | |
108 | if(ff1.gt.fb) goto 11 | |
109 | b1=bb1 | |
110 | end if | |
111 | bgen=b1 | |
112 | ||
113 | * generate single event with partonic energy loss | |
114 | nrg=0 | |
115 | ehard=ckin(3) | |
116 | call plinit(ehard) | |
117 | call plevnt(ehard) | |
118 | ||
119 | * reset all in-vacuum radiated guark 4-momenta and codes to zero | |
120 | do i=1,1000 | |
121 | do j=1,5 | |
122 | pqua(i,j)=0.d0 | |
123 | kqua(i,j)=0 | |
124 | end do | |
125 | end do | |
126 | nrq=0 | |
127 | ||
128 | * generate final state shower in vacuum if it was excluded before | |
129 | nrgm=nrg ! fix number of in-medium emitted gluons | |
130 | ip1=0 | |
131 | ip2=0 | |
132 | ip01=0 | |
133 | ip02=0 | |
134 | ip001=0 | |
135 | ip002=0 | |
136 | if(mstj(41).ne.0) goto 5 | |
137 | mstj(41)=1 | |
138 | nn=n | |
139 | do i=9,nn | |
140 | if(k(i,3).eq.7) then | |
141 | ip1=i ! first hard parton (line ip1) | |
142 | kfh1=k(i,1) ! status code of first hard parton | |
143 | qmax1=pyp(i,10) ! transverse momentum of first hard parton | |
144 | end if | |
145 | if(k(i,3).eq.8) then | |
146 | ip2=i ! second hard parton (line ip2) | |
147 | kfh2=k(i,1) ! status code of second hard parton | |
148 | qmax2=pyp(i,10) ! transverse momentum of second hard parton | |
149 | end if | |
150 | end do | |
151 | ||
152 | n1=n | |
153 | call pyshow(ip1,0,qmax1) ! vacuum showering for first hard parton | |
154 | if(n.eq.n1) ip1=0 | |
155 | n2=n | |
156 | call pyshow(ip2,0,qmax2) ! vacuum showering for second hard parton | |
157 | if(n.eq.n2) ip2=0 | |
158 | mstj(41)=0 | |
159 | if(n.eq.nn) goto 5 | |
160 | ||
161 | * find two leading partons after showering | |
162 | do i=nn+1,n | |
163 | if(k(i,3).eq.ip1) ip001=i ! first daughter of first hard parton | |
164 | if(k(i,3).eq.ip2) ip002=i ! first daughter of second hard parton | |
165 | end do | |
166 | ptle1=0.d0 | |
167 | ptle2=0.d0 | |
168 | do i=nn+1,n | |
169 | if (k(i,1).eq.14) goto 3 | |
170 | if(i.ge.ip002.and.ip002.gt.0) then | |
171 | ptl02=pyp(i,10) | |
172 | if(ptl02.gt.ptle2.and.k(i,2).eq.k(ip2,2)) then | |
173 | ip02=i ! leading parton in second shower (line ip02) | |
174 | ptle2=ptl02 ! pt of the leading parton | |
175 | end if | |
176 | elseif(ip001.gt.0) then | |
177 | ptl01=pyp(i,10) | |
178 | if(ptl01.gt.ptle1.and.k(i,2).eq.k(ip1,2)) then | |
179 | ip01=i ! leading parton in first shower (line ip01) | |
180 | ptle1=ptl01 ! pt of the leading parton | |
181 | end if | |
182 | end if | |
183 | 3 continue | |
184 | end do | |
185 | ||
186 | * replace two hard partons by two leading partons in original event record | |
187 | if(ip1.gt.0) then | |
188 | do j=1,5 | |
189 | v(ip1,j)=v(ip01,j) | |
190 | p(ip1,j)=p(ip01,j) | |
191 | end do | |
192 | k(ip1,1)=kfh1 | |
193 | * fix first/last daughter for moving entry | |
194 | do jgl=1,n | |
195 | if(k(jgl,4).eq.ip01) k(jgl,4)=ip1 | |
196 | if(k(jgl,5).eq.ip01) k(jgl,5)=ip1 | |
197 | end do | |
198 | * | |
199 | end if | |
200 | if(ip2.gt.0) then | |
201 | do j=1,5 | |
202 | v(ip2,j)=v(ip02,j) | |
203 | p(ip2,j)=p(ip02,j) | |
204 | end do | |
205 | k(ip2,1)=kfh2 | |
206 | * fix first/last daughter for moving entry | |
207 | do jgl=1,n | |
208 | if(k(jgl,4).eq.ip02) k(jgl,4)=ip2 | |
209 | if(k(jgl,5).eq.ip02) k(jgl,5)=ip2 | |
210 | end do | |
211 | * | |
212 | end if | |
213 | ||
214 | * add final showering gluons to the list of in-medium emitted gluons, | |
215 | * fill the list of emitted quarks by final showering quark pairs, | |
216 | * and remove showering gluons and quarks from the event record | |
217 | do i=nn+1,n | |
218 | if(k(i,1).eq.14.or.i.eq.ip01.or.i.eq.ip02) goto 12 | |
219 | if(k(i,2).ne.21) then ! filling 'plquar' arrays for quarks | |
220 | nrq=nrq+1 | |
221 | do j=1,5 | |
222 | kqua(nrq,j)=k(i,j) | |
223 | pqua(nrq,j)=p(i,j) | |
224 | end do | |
225 | kqua(nrq,1)=2 | |
226 | goto 12 | |
227 | end if | |
228 | if(i.ge.ip002.and.ip002.gt.0) then | |
229 | ish=ip2 | |
230 | else | |
231 | ish=ip1 | |
232 | end if | |
233 | nrg=nrg+1 | |
234 | nur=nrg | |
235 | 7 ishm=kglu(nur-1,6) | |
236 | if(ish.ge.ishm.or.nur.le.2) goto 6 ! adding gluons in 'plglur' arrays | |
237 | do j=1,6 | |
238 | kglu(nur,j)=kglu(nur-1,j) | |
239 | end do | |
240 | do j=1,4 | |
241 | glur(nur,j)=glur(nur-1,j) | |
242 | end do | |
243 | nur=nur-1 | |
244 | goto 7 | |
245 | 6 kglu(nur,1)=2 ! status code | |
246 | kglu(nur,2)=k(i,2) ! particle identificator | |
247 | kglu(nur,3)=k(ish,3) ! parent line number | |
248 | kglu(nur,4)=0 ! special colour info | |
249 | kglu(nur,5)=0 ! special colour info | |
250 | kglu(nur,6)=ish ! associated parton number | |
251 | glur(nur,1)=p(i,4) ! energy | |
252 | glur(nur,2)=pyp(i,10) ! pt | |
253 | glur(nur,3)=pyp(i,15) ! phi | |
254 | glur(nur,4)=pyp(i,19) ! eta | |
255 | 12 continue | |
256 | * remove partons from event list | |
257 | do j=1,5 | |
258 | v(i,j)=0.d0 | |
259 | k(i,j)=0 | |
260 | p(i,j)=0.d0 | |
261 | end do | |
262 | end do | |
263 | n=nn | |
264 | ||
265 | 5 continue | |
266 | ||
267 | * stop generate event if there are no additional gluons | |
268 | if(nrg.lt.1) goto 1 | |
269 | ||
270 | * define number of stirngs (ns) and number of entries in strings before | |
271 | * in-medium radiation (nis(ns)) | |
272 | ns=0 | |
273 | nes=0 | |
274 | i0=0 | |
275 | i1=0 | |
276 | do i=1,500 | |
277 | nis(i)=0 | |
278 | nas(i)=0 | |
279 | nss(i)=0 | |
280 | nus(i)=0 | |
281 | end do | |
282 | do i=9,n | |
283 | ks=k(i,1) | |
284 | ksp=k(i-1,1) | |
285 | if(ks.eq.2) then | |
286 | nis(ns+1)=nis(ns+1)+1 | |
287 | elseif(ks.eq.1.and.nis(ns+1).gt.0) then | |
288 | nis(ns+1)=nis(ns+1)+1 | |
289 | nes=nes+nis(ns+1) ! nes - total number of entries | |
290 | nss(ns+1)=nes | |
291 | ns=ns+1 | |
292 | elseif(ks.ne.2.and.ksp.ne.2.and.ns.gt.0) then | |
293 | i1=i1+1 ! last i1 lines not included in strings | |
294 | end if | |
295 | end do | |
296 | i0=n-nes-i1 ! first i0 lines not included in strings | |
297 | do i=1,ns | |
298 | nss(i)=nss(i)+i0 | |
299 | end do | |
300 | ||
301 | * move fragmented particles in bottom of event list | |
302 | i=i0+1 | |
303 | 2 ks=k(i,1) | |
304 | ksp=k(i-1,1) | |
305 | if(ks.ne.2.and.ksp.ne.2) then | |
306 | n=n+1 | |
307 | do j=1,5 | |
308 | v(n,j)=v(i,j) | |
309 | k(n,j)=k(i,j) | |
310 | p(n,j)=p(i,j) | |
311 | end do | |
312 | * fix first/last daughter for moving entry | |
313 | do jgl=1,n | |
314 | if(k(jgl,4).eq.i) k(jgl,4)=n | |
315 | if(k(jgl,5).eq.i) k(jgl,5)=n | |
316 | end do | |
317 | * | |
318 | do in=i+1,n | |
319 | do j=1,5 | |
320 | v(in-1,j)=v(in,j) | |
321 | k(in-1,j)=k(in,j) | |
322 | p(in-1,j)=p(in,j) | |
323 | end do | |
324 | * fix first/last daughter for moving entry | |
325 | do jgl=1,n | |
326 | if(k(jgl,4).eq.in) k(jgl,4)=in-1 | |
327 | if(k(jgl,5).eq.in) k(jgl,5)=in-1 | |
328 | end do | |
329 | * | |
330 | end do | |
331 | do ip=1,nrg | |
332 | ku=kglu(ip,6) | |
333 | if(ku.gt.i) kglu(ip,6)=ku-1 | |
334 | end do | |
335 | n=n-1 | |
336 | else | |
337 | i=i+1 | |
338 | end if | |
339 | if(i.le.n-i1) goto 2 | |
340 | ||
341 | * define number of additional entries in strings, nas(ns) | |
342 | do i=1,nrg | |
343 | kas=kglu(i,6) | |
344 | if(kas.le.nss(1)) then | |
345 | nas(1)=nas(1)+1 | |
346 | else | |
347 | do j=2,ns | |
348 | if(kas.le.nss(j).and.kas.gt.nss(j-1)) | |
349 | > nas(j)=nas(j)+1 | |
350 | end do | |
351 | end if | |
352 | end do | |
353 | do j=1,ns | |
354 | do i=1,j | |
355 | nus(j)=nus(j)+nas(i) | |
356 | end do | |
357 | end do | |
358 | ||
359 | * add emitted gluons in event list | |
360 | nu=n | |
361 | n=n+nrg | |
362 | do i=nu,nu-i1,-1 | |
363 | is=i+nrg | |
364 | do j=1,5 | |
365 | v(is,j)=v(i,j) | |
366 | k(is,j)=k(i,j) | |
367 | p(is,j)=p(i,j) | |
368 | end do | |
369 | * fix first/last daughter for moving entries | |
370 | do jgl=1,n | |
371 | if(k(jgl,4).eq.i) k(jgl,4)=is | |
372 | if(k(jgl,5).eq.i) k(jgl,5)=is | |
373 | end do | |
374 | * | |
375 | end do | |
376 | ||
377 | do ia=ns-1,1,-1 | |
378 | do i=nss(ia+1)-1,nss(ia),-1 | |
379 | is=i+nus(ia) | |
380 | do j=1,5 | |
381 | v(is,j)=v(i,j) | |
382 | k(is,j)=k(i,j) | |
383 | p(is,j)=p(i,j) | |
384 | end do | |
385 | * fix first/last daughter for moving entries | |
386 | do jgl=1,n | |
387 | if(k(jgl,4).eq.i) k(jgl,4)=is | |
388 | if(k(jgl,5).eq.i) k(jgl,5)=is | |
389 | end do | |
390 | * | |
391 | end do | |
392 | end do | |
393 | ||
394 | do i=1,nrg | |
395 | if(i.le.nus(1)) then | |
396 | ia=nss(1)-1+i | |
397 | else | |
398 | do in=2,ns | |
399 | if(i.le.nus(in).and.i.gt.nus(in-1)) | |
400 | > ia=nss(in)-1+i | |
401 | end do | |
402 | end if | |
403 | eg=glur(i,1) | |
404 | ptg=glur(i,2) | |
405 | phig=glur(i,3) | |
406 | etag=glur(i,4) | |
407 | do j=1,5 | |
408 | v(ia,j)=0.d0 | |
409 | k(ia,j)=kglu(i,j) | |
410 | end do | |
411 | p(ia,1)=ptg*dcos(phig) | |
412 | p(ia,2)=ptg*dsin(phig) | |
413 | p(ia,3)=dsqrt(abs(eg*eg-ptg*ptg)) | |
414 | if(etag.lt.0.d0) p(ia,3)=-1.d0*p(ia,3) | |
415 | p(ia,4)=dsqrt(ptg*ptg+p(ia,3)**2) | |
416 | p(ia,5)=0.d0 | |
417 | end do | |
418 | ||
419 | * rearrange partons to form strings in event list | |
420 | do ij=1,1000 | |
421 | ijoin(ij)=0 | |
422 | ijoin0(ij)=0 | |
423 | end do | |
424 | do i=1,ns | |
425 | njoin=nis(i)+nas(i) | |
426 | if(i.eq.1) then | |
427 | do j=1,njoin | |
428 | ijoin(j)=i0+j | |
429 | end do | |
430 | else | |
431 | do j=1,njoin | |
432 | ijoin(j)=nss(i-1)+nus(i-1)+j | |
433 | end do | |
434 | end if | |
435 | ||
436 | * re-oder additional gluons by z-coordinate along the string | |
437 | if(nas(i).gt.0) then | |
438 | ja=njoin-nas(i) | |
439 | jo1=ijoin(1) | |
440 | jon=ijoin(njoin) | |
441 | etasum=0.d0 | |
442 | detaj=pyp(jo1,19)-pyp(jon,19) | |
443 | do j=ja,njoin-1 | |
444 | jnum=0 | |
445 | etaj=pyp(jo1+j-1,19) | |
446 | etasum=etasum+etaj | |
447 | do jj=ja,njoin-1 | |
448 | etajj=pyp(jo1+jj-1,19) | |
449 | if(detaj.lt.0) then | |
450 | if(etajj.lt.etaj.and.j.ne.jj) jnum=jnum+1 | |
451 | else | |
452 | if(etajj.gt.etaj.and.j.ne.jj) jnum=jnum+1 | |
453 | end if | |
454 | if(etajj.eq.etaj.and.j.lt.jj) jnum=jnum+1 | |
455 | end do | |
456 | ijoin(ja+jnum)=jo1+j-1 | |
457 | end do | |
458 | detas1=abs(pyp(jo1,19)-etasum) | |
459 | detasn=abs(pyp(jon,19)-etasum) | |
460 | if(detasn.gt.detas1) then | |
461 | do j=1,njoin | |
462 | ijoin0(j)=ijoin(j) | |
463 | end do | |
464 | do j=2,nas(i)+1 | |
465 | ijoin(j)=ijoin0(ja+j-2) | |
466 | end do | |
467 | do j=nas(i)+2,njoin-1 | |
468 | ijoin(j)=ijoin0(j-nas(i)) | |
469 | end do | |
470 | end if | |
471 | end if | |
472 | ||
473 | * form strings | |
474 | call pyjoin(njoin,ijoin) | |
475 | ||
476 | end do | |
477 | ||
478 | * add in-vacuum emitted quark pairs | |
479 | if(nrq.lt.2) goto 1 | |
480 | do i=1,nrq,2 | |
481 | n=n+2 | |
482 | do j=1,5 | |
483 | v(n-1,j)=0.d0 | |
484 | k(n-1,j)=kqua(i,j) | |
485 | p(n-1,j)=pqua(i,j) | |
486 | end do | |
487 | in=i+1 | |
488 | 4 ktest=k(n-1,2)+kqua(in,2) | |
489 | if(ktest.eq.0.or.in.eq.nrq) goto 8 | |
490 | in=in+1 | |
491 | goto 4 | |
492 | 8 do j=1,5 | |
493 | v(n,j)=0.d0 | |
494 | k(n,j)=kqua(in,j) | |
495 | p(n,j)=pqua(in,j) | |
496 | end do | |
497 | if(in.gt.i+1) then | |
498 | do j=1,5 | |
499 | kqua(in,j)=kqua(i+1,j) | |
500 | pqua(in,j)=pqua(i+1,j) | |
501 | end do | |
502 | end if | |
503 | end do | |
504 | ||
505 | do ij=1,2 | |
506 | ijoik(ij)=0 | |
507 | end do | |
508 | do i=1,nrq-1,2 | |
509 | k(n+1-i,1)=1 | |
510 | ijoik(1)=n-i | |
511 | ijoik(2)=n+1-i | |
512 | call pyjoin(2,ijoik) | |
513 | end do | |
514 | ||
515 | 1 continue | |
516 | ||
517 | return | |
518 | end | |
519 | ||
520 | ********************************* PLINIT *************************** | |
521 | SUBROUTINE PLINIT(ET) | |
522 | * set time-dependence of plasma parameters | |
523 | IMPLICIT DOUBLE PRECISION(A-H, O-Z) | |
524 | external plvisc | |
525 | common /plpar1/ tau0,T0,TC,sigqq,AW,RA,mvisc,nf | |
526 | common /plpar2/ pln0,taupl,tauh,sigpl,sigh,sigplh,sigqqh,rg,rgn | |
527 | common /plevol/ taup(5000),temp(5000),denp(5000),enep(5000) | |
528 | save /plevol/,/plpar1/,/plpar2/ | |
529 | * | |
530 | pi=3.14159d0 | |
531 | pi2=pi*pi | |
532 | ||
533 | * set number degrees of freedom in QGP | |
534 | hgd=3.d0 | |
535 | rg=(16.d0+10.5d0*nf)/hgd | |
536 | rgn=(16.d0+9.d0*nf)/hgd | |
537 | ||
538 | * set 'fiction' sigma for parton rescattering in QGP | |
539 | sigqq=4.2d0 | |
540 | sigpl=2.25d0*2.25d0*sigqq*(16.d0+4.d0*nf)/(16.d0+9.d0*nf) | |
541 | ||
542 | * set intial plasma temperature, density and energy density in perfect | |
543 | * (if mvisc=0) or viscous (mvisc=1,2) QGP with PLVISC subroitine | |
544 | hst=0.15d0 | |
545 | if(T0.gt.1.5d0.or.mvisc.eq.2) hst=0.25d0 | |
546 | if(T0.gt.1.5d0.and.mvisc.ne.0) hst=0.9d0 | |
547 | T01=T0*5.06d0 | |
548 | TC1=TC*5.06d0 | |
549 | pln0=(16.d0+9.d0*nf)*1.2d0*(T01**3)/pi2 | |
550 | ened0=pi2*(16.d0+10.5d0*nf)*(T01**4)/30.d0 | |
551 | hh=hst*tau0 | |
552 | tau=tau0 ! proper time | |
553 | T=T01 ! temperature | |
554 | den=pln0 ! number density | |
555 | ened=ened0 ! energy density | |
556 | ||
557 | * create array of parameters to configurate QGP time evolution | |
558 | DO I=1,5000 | |
559 | taup(i)=tau ! proper time | |
560 | temp(i)=T/5.06d0 ! temperature | |
561 | denp(i)=den ! number density | |
562 | enep(i)=ened/5.06d0 ! energy density | |
563 | ened1=0.5d0*hh*(1.3333d0*plvisc(T)/(tau*tau)-1.3333d0 | |
564 | > *ened/tau)+ened | |
565 | T1=(30.d0*ened1/((16.d0+10.5d0*nf)*pi2))**0.25d0 | |
566 | tau1=tau+0.5d0*hh | |
567 | ened=hh*(1.3333d0*plvisc(T1)/(tau1*tau1)-1.3333d0 | |
568 | > *ened1/tau1)+ened | |
569 | TPR=T | |
570 | T=(30.d0*ened/((16.d0+10.5d0*nf)*pi2))**0.25d0 | |
571 | den=(16.d0+9.d0*nf)*1.2d0*(T**3)/pi2 | |
572 | tau=tau+hh | |
573 | if(TPR.gt.TC1.and.T.le.TC1) taupl=tau-0.5d0*hh ! QGP lifetime taupl | |
574 | END DO | |
575 | tauh=taupl*rg ! mixed phase lifetime | |
576 | ||
577 | return | |
578 | end | |
579 | ******************************** END PLINIT ************************** | |
580 | ||
581 | ******************************* PLEVNT ****************************** | |
582 | SUBROUTINE PLEVNT(ET) | |
583 | * generate hard parton production vertex and passing with rescattering, | |
584 | * collisional and radiative energy loss of each parton through plasma | |
585 | IMPLICIT DOUBLE PRECISION(A-H, O-Z) | |
586 | IMPLICIT INTEGER(I-N) | |
587 | INTEGER PYK,PYCHGE,PYCOMP | |
588 | external plthik, pln, plt, pls, gauss, gluang | |
589 | external pyp,pyr,pyk | |
590 | common /plpar1/ tau0,T0,TC,sigqq,AW,RA,mvisc,nf | |
591 | common /plpar2/ pln0,taupl,tauh,sigpl,sigh,sigplh,sigqqh,rg,rgn | |
592 | common /thikpa/ fmax,xmin | |
593 | common /pyjets/ n,npad,k(4000,5),p(4000,5),v(4000,5) | |
594 | common /pyqpar/ T0u,tau0u,nfu,ienglu,ianglu | |
595 | common /plglur/ glur(1000,4),kglu(1000,6),nrg,nrgm | |
596 | common /factor/ cfac, kf | |
597 | common /pleave/ taul, temlev | |
598 | common /parimp/ b1, psib1, rb1, rb2, noquen | |
599 | common /plen/ epartc, um | |
600 | common /plos/ elr,rsk | |
601 | common /numje1/ nuj1, nuj2 | |
602 | save/pyjets/,/plglur/,/plpar1/,/plpar2/,/thikpa/,/factor/, | |
603 | < /pleave/,/parimp/,/plen/,/plos/,/numje1/ | |
604 | * | |
605 | pi=3.14159d0 | |
606 | ||
607 | * find minimum of nuclear thikness function with subroutine plsear | |
608 | psib1=pi*(2.d0*pyr(0)-1.d0) | |
609 | call plsear (fmax1,xmin1) | |
610 | fmax=fmax1 | |
611 | xmin=xmin1 | |
612 | ||
613 | * generate vertex of jet production | |
614 | iv=0 | |
615 | 1 rr1=xmin*pyr(0) | |
616 | ff1=fmax*pyr(0) | |
617 | f=plthik(rr1) | |
618 | iv=iv+1 | |
619 | if(ff1.gt.f.and.iv.le.100000) goto 1 | |
620 | r0=rr1 | |
621 | rb1=dsqrt(abs(r0*r0+b1*b1/4.d0+r0*b1*dcos(psib1))) | |
622 | rb2=dsqrt(abs(r0*r0+b1*b1/4.d0-r0*b1*dcos(psib1))) | |
623 | rb1=max(rb1,1.d-4) | |
624 | rb2=max(rb2,1.d-4) | |
625 | * no quenching if noquen=1 or jet production vertex is out of effective dense zone | |
626 | if(noquen.eq.1.or.rb1.gt.RA.or.rb2.gt.RA) goto 7 | |
627 | ||
628 | * find maximum of angular spectrum of radiated gluons with subroutine gluang | |
629 | temin=0.5d0*pi | |
630 | temax=0.5d0*(1.d0+dsqrt(5.d0))*0.0863d0 | |
631 | ftemax=gluang(temax) | |
632 | ||
633 | * reset all radiated gluon 4-momenta and codes to zero ------------------- | |
634 | do i=1,1000 | |
635 | do j=1,4 | |
636 | glur(i,j)=0.d0 | |
637 | kglu(i,j)=0 | |
638 | end do | |
639 | kglu(i,5)=0 | |
640 | kglu(i,6)=0 | |
641 | end do | |
642 | nrg=0 | |
643 | ||
644 | * generate changing 4-momentum of partons due to rescattering and energy loss | |
645 | * (for partons with |eta|<3.5 and pt>3 GeV/c) | |
646 | nuj1=9 ! minimum line number of rescattered parton | |
647 | nuj2=n ! maximum line number of rescattered parton | |
648 | do 2 ip=nuj1,nuj2 ! cycle on travelling partons | |
649 | irasf=0 | |
650 | iraz=0 | |
651 | ks=k(ip,1) ! parton status code | |
652 | kf=k(ip,2) ! parton identificator | |
653 | ka=abs(kf) | |
654 | ko=k(ip,3) ! origin (parent line number) | |
655 | epart=abs(pyp(ip,10)) ! parton transverse momentum | |
656 | etar=pyp(ip,19) ! parton pseudorapidity | |
657 | if(ko.gt.6.and.epart.ge.3.d0.and.abs(etar). | |
658 | > le.7.d0) then | |
659 | if(ka.eq.21.or.ka.eq.1.or.ka.eq.2.or.ka.eq.3. | |
660 | > or.ka.eq.4.or.ka.eq.5.or.ka.eq.6.or.ka.eq.7. | |
661 | > or.ka.eq.8) then | |
662 | if(ks.eq.2.or.ks.eq.1.or.ks.eq.21) then | |
663 | phir=pyp(ip,15) ! parton azimuthal angle | |
664 | tetr=pyp(ip,13) ! parton polar angle | |
665 | yrr=pyp(ip,17) ! parton rapidity | |
666 | stetr=dsin(tetr) ! parton sin(theta) | |
667 | if(abs(stetr).le.1.d-05) then | |
668 | if(stetr.ge.0.d0) then | |
669 | stetr=1.d-05 | |
670 | else | |
671 | stetr=-1.d-05 | |
672 | end if | |
673 | end if | |
674 | phir1=-1.d0*phir | |
675 | tetr1=-1.d0*tetr | |
676 | ||
677 | * set colour factor | |
678 | if(kf.eq.21) then | |
679 | cfac=1.d0 ! for gluon | |
680 | else | |
681 | cfac=0.44444444d0 ! for quark | |
682 | end if | |
683 | ||
684 | * boost from laboratory system to system of hard parton | |
685 | ipar=ip | |
686 | bet0=(r0*dcos(psib1)+0.5d0*b1)/rb1 | |
687 | if(bet0.le.-1.d0) bet0=-0.99999d0 | |
688 | if(bet0.ge.1.d0) bet0=0.99999d0 | |
689 | bet=dacos(bet0) | |
690 | if(psib1.lt.0.d0) bet=-1.d0*bet | |
691 | phip=phir-bet | |
692 | if(phip.gt.pi) phip=phip-2.d0*pi | |
693 | if(phip.lt.-1.d0*pi) phip=phip+2.d0*pi | |
694 | call pyrobo(ip,ip,0.d0,phir1,0.d0,0.d0,0.d0) | |
695 | call pyrobo(ip,ip,tetr1,0.d0,0.d0,0.d0,0.d0) | |
696 | ||
697 | * calculate proper time of parton leaving the effective dense zone | |
698 | aphin=(r0*r0-b1*b1/4.d0)/(rb1*rb2) | |
699 | if(aphin.le.-1.d0) aphin=-0.99999d0 | |
700 | if(aphin.ge.1.d0) aphin=0.99999d0 | |
701 | phin=dacos(aphin) | |
702 | if(psib1.le.0.d0) phin=-1.d0*phin | |
703 | phid=phip-phin | |
704 | if(phid.gt.pi) phid=phid-2.d0*pi | |
705 | if(phid.lt.-1.d0*pi) phid=phid+2.d0*pi | |
706 | taul1=abs(dsqrt(abs(RA*RA-(rb1*dsin(phip))**2))-rb1*dcos(phip)) | |
707 | taul2=abs(dsqrt(abs(RA*RA-(rb2*dsin(phid))**2))-rb2*dcos(phid)) | |
708 | taul=min(taul1,taul2) ! escape time taul | |
709 | temlev=plt(taul,rb1,rb2,yrr) ! QGP temperature at taul | |
710 | if(taul.le.tau0) goto 100 ! escape from QGP if taul<tau0 | |
711 | ||
712 | * start parton rescattering in QGP with proper time iterations | |
713 | tau=tau0 | |
714 | xj=r0*dcos(psib1) | |
715 | yj=r0*dsin(psib1) | |
716 | rj1=rb1 | |
717 | rj2=rb2 | |
718 | 3 tfs=plt(tau,rj1,rj2,yrr) | |
719 | xi=-10.d0*dlog(max(pyr(0),1.d-10))/(sigpl*pln(tau,rj1,rj2,yrr)) | |
720 | vel=abs(p(ip,3))/dsqrt(p(ip,3)**2+p(ip,5)**2) ! parton velocity | |
721 | if(vel.lt.0.3d0) goto 4 | |
722 | tau=tau+xi*vel | |
723 | xj=xj+xi*vel*dcos(phir) | |
724 | yj=yj+xi*vel*dsin(phir) | |
725 | rj1=sqrt(abs(yj**2+(xj+0.5d0*b1)**2)) | |
726 | rj2=sqrt(abs(yj**2+(xj-0.5d0*b1)**2)) | |
727 | if(tfs.le.TC) goto 100 ! escape if temperature drops below TC | |
728 | ||
729 | * transform parton 4-momentum due to next scattering with subroutine pljetr | |
730 | epartc=p(ip,4) ! parton energy | |
731 | um=p(ip,5) ! parton mass | |
732 | sigtr=pls(tfs)*cfac*((p(ip,4)/pyp(ip,8))**2) | |
733 | prob=sigpl/(sigtr/stetr+sigpl) | |
734 | ran=pyr(0) | |
735 | irasf=irasf+1 | |
736 | if(irasf.gt.100000) goto 100 | |
737 | if(ran.lt.prob) goto 3 | |
738 | pltp=plt(tau,rj1,rj2,yrr) | |
739 | pltp3=3.d0*pltp | |
740 | pass=50.6d0/(pln(tau,rj1,rj2,yrr)*sigtr) | |
741 | elr=0.d0 | |
742 | rsk=0.d0 | |
743 | call pljetr(tau,pass,pltp,ipar,epart) | |
744 | irasf=0 | |
745 | ||
746 | * set 4-momentum (in lab system) of next radiated gluon for parton number >8 | |
747 | * and fill arrays of radiated gluons in common block plglur | |
748 | if(nrg.le.1000) then | |
749 | if(abs(elr).gt.0.1d0.and.ip.gt.8) then | |
750 | * generate the angle of emitted gluon | |
751 | if(ianglu.eq.0) then | |
752 | 6 te1=temin*pyr(0) | |
753 | fte1=ftemax*pyr(0) | |
754 | fte2=gluang(te1) | |
755 | if(fte1.gt.fte2) goto 6 | |
756 | tgl=te1 | |
757 | elseif (ianglu.eq.1) then | |
758 | tgl=((0.5d0*pi*epartc)**pyr(0))/epartc | |
759 | else | |
760 | tgl=0.d0 | |
761 | end if | |
762 | pgl=pi*(2.d0*pyr(0)-1.d0) | |
763 | * in comoving system | |
764 | pxgl=abs(elr)*stetr*(dcos(phir)*dcos(tgl)- | |
765 | > dsin(phir)*dsin(tgl)*dsin(pgl)) | |
766 | pygl=abs(elr)*stetr*(dsin(phir)*dcos(tgl)+ | |
767 | > dcos(phir)*dsin(tgl)*dsin(pgl)) | |
768 | pzgl=-1.d0*abs(elr)*stetr*dsin(tgl)*dcos(pgl) | |
769 | ptgl=dsqrt(abs(pxgl*pxgl+pygl*pygl)) | |
770 | psgl=dsqrt(abs(ptgl*ptgl+pzgl*pzgl)) | |
771 | * recalculate in lab system | |
772 | dyg=0.5d0*dlog(max(1.d-9,(psgl+pzgl)/(psgl-pzgl))) | |
773 | pzgl=ptgl*dsinh(yrr+dyg) | |
774 | psgl=dsqrt(abs(ptgl*ptgl+pzgl*pzgl)) | |
775 | * | |
776 | dpgl=pygl/pxgl | |
777 | glur1=abs(elr) ! energy | |
778 | glur3=datan(dpgl) ! phi | |
779 | if(pxgl.lt.0.d0) then | |
780 | if(pygl.ge.0.d0) then | |
781 | glur3=glur3+pi | |
782 | else | |
783 | glur3=glur3-pi | |
784 | end if | |
785 | end if | |
786 | glur4=0.5d0*dlog(max(1.d-9,(psgl+pzgl)/(psgl-pzgl))) ! eta | |
787 | glur2=glur1/dcosh(glur4) ! pt | |
788 | ||
789 | * put in event list radiated gluons with pt > 0.2 GeV only | |
790 | if(glur2.ge.0.2d0) then | |
791 | nrg=nrg+1 | |
792 | * set gluon 4-momentum | |
793 | glur(nrg,1)=glur1 ! energy | |
794 | glur(nrg,2)=glur2 ! pt | |
795 | glur(nrg,3)=glur3 ! phi | |
796 | glur(nrg,4)=glur4 ! eta | |
797 | * set gluon codes | |
798 | kglu(nrg,1)=2 ! status code | |
799 | kglu(nrg,2)=21 ! particle identificator | |
800 | kglu(nrg,3)=k(ipar,3) ! parent line number | |
801 | kglu(nrg,4)=0 ! special colour info | |
802 | kglu(nrg,5)=0 ! special colour info | |
803 | kglu(nrg,6)=ipar ! associated parton number | |
804 | end if | |
805 | end if | |
806 | else | |
807 | write(6,*) 'Warning! Number of emitted gluons is too large!' | |
808 | end if | |
809 | ||
810 | * set parton "thermalization" if pt<T | |
811 | if(abs(p(ip,3)).gt.pltp3) goto 3 | |
812 | 4 continue | |
813 | if(p(ip,3).ge.0.d0) then | |
814 | sigp=1.d0 | |
815 | else | |
816 | sigp=-1.d0 | |
817 | end if | |
818 | 5 iraz=iraz+1 | |
819 | if(iraz.gt.100000) goto 100 | |
820 | ep0=-0.15d0*(dlog(max(1.d-10,pyr(0)))+dlog(max(1.d-10,pyr(0)))+ | |
821 | > dlog(max(1.d-10,pyr(0)))) | |
822 | if(ep0.le.p(ip,5).or.ep0.ge.100.d0) goto 5 | |
823 | pp0=dsqrt(abs(ep0**2-p(ip,5)**2)) | |
824 | probt=pp0/ep0 | |
825 | if(pyr(0).gt.probt) goto 5 | |
826 | ctp0=2.d0*pyr(0)-1.d0 | |
827 | stp0=dsqrt(abs(1.d0-ctp0**2)) | |
828 | php0=pi*(2.d0*pyr(0)-1.d0) | |
829 | p(ip,1)=pp0*stp0*dcos(php0) | |
830 | p(ip,2)=pp0*stp0*dsin(php0) | |
831 | p(ip,3)=sigp*pp0*ctp0 | |
832 | p(ip,4)=dsqrt(p(ip,1)**2+p(ip,2)**2+p(ip,3)**2+p(ip,5)**2) | |
833 | ||
834 | * boost to laboratory system | |
835 | 100 call pyrobo(ip,ip,tetr,phir,0.d0,0.d0,0.d0) | |
836 | end if | |
837 | end if | |
838 | end if | |
839 | 2 continue | |
840 | 7 continue | |
841 | ||
842 | return | |
843 | end | |
844 | ******************************* END PLEVNT ************************* | |
845 | ||
846 | ******************************* PLJETR ***************************** | |
847 | SUBROUTINE PLJETR(tau,y,x,ip,epart) | |
848 | * transform parton 4-momentum due to scattering in plasma at time = tau | |
849 | IMPLICIT DOUBLE PRECISION(A-H, O-Z) | |
850 | IMPLICIT INTEGER(I-N) | |
851 | INTEGER PYK,PYCHGE,PYCOMP | |
852 | external plfun1, pls | |
853 | external pyp,pyr | |
854 | common /plpar1/ tau0,T0,TC,sigqq,AW,RA,mvisc,nf | |
855 | common /plpar2/ pln0,taupl,tauh,sigpl,sigh,sigplh,sigqqh,rg,rgn | |
856 | common /pyjets/ n,npad,k(4000,5),p(4000,5),v(4000,5) | |
857 | common /pyqpar/ T0u,tau0u,nfu,ienglu,ianglu | |
858 | common /pljdat/ ej, z, ygl, alfs, um, epa | |
859 | common /pleave/ taul, temlev | |
860 | common /radcal/ aa, bb | |
861 | common /factor/ cfac, kf | |
862 | common /plos/ elr,rsk | |
863 | save /pyjets/,/plpar1/,/plpar2/,/pyqpar/,/pljdat/,/pleave/,/plos/, | |
864 | < /factor/,/radcal/ | |
865 | * | |
866 | pi=3.14159d0 | |
867 | spi=dsqrt(pi) | |
868 | tauu=x ! redenote temerature tauu=x | |
869 | i=ip ! redenote parton number i=ip | |
870 | iter=0 | |
871 | iraz=0 | |
872 | ||
873 | * boost to system of comoving plasma constituent | |
874 | phir=pyp(i,15) ! parton phi | |
875 | tetr=pyp(i,13) ! parton theta | |
876 | phir1=-1.d0*phir | |
877 | tetr1=-1.d0*tetr | |
878 | call pyrobo(i,i,0.d0,phir1,0.d0,0.d0,0.d0) | |
879 | call pyrobo(i,i,tetr1,0.d0,0.d0,0.d0,0.d0) | |
880 | pp=pyp(i,8) ! parton total momentum | |
881 | ppl=abs(p(i,3)) ! parton pz | |
882 | um=p(i,5) ! parton mass | |
883 | epa=p(i,4) ! parton energy | |
884 | ppt=pyp(i,10) ! parton pt | |
885 | pphi=pyp(i,15) ! parton phi | |
886 | ||
887 | if(ppl.lt.3.d0) goto 222 ! no energy loss if pz<3 GeV/c | |
888 | ||
889 | * generation hard parton-plasma scattering with momentum transfer rsk | |
890 | 221 ep0=-1.*tauu*(dlog(max(1.d-10,pyr(0)))+dlog(max(1.d-10, | |
891 | > pyr(0)))+dlog(max(1.d-10,pyr(0)))) ! energy of 'thermal' parton | |
892 | iter=iter+1 | |
893 | if(ep0.lt.1.d-10.and.iter.le.100000) goto 221 | |
894 | scm=2.*ep0*epa+um*um+ep0*ep0 | |
895 | qm2=(scm-((um+ep0)**2))*(scm-((um-ep0)**2))/scm | |
896 | bub=4.d0*tauu/TC | |
897 | alf=6.d0*pi/((33.d0-2.d0*nf)*dlog(max(bub,1.d-10))) | |
898 | z=pi*4.d0*tauu*tauu*alf*(1.+nf/6.d0) | |
899 | bubs=dsqrt(abs(z))/TC | |
900 | alfs=6.d0*pi/((33.d0-2.d0*nf)*dlog(max(bubs,1.d-10))) | |
901 | phmin2=z | |
902 | phmax2=max(phmin2,qm2) | |
903 | fqmax2=1.d0/(dlog(max(phmin2/(TC*TC),1.d-10)))**2 | |
904 | 12 rn1=pyr(0) | |
905 | tp=1.d0/(rn1/phmax2+(1.d0-rn1)/phmin2) | |
906 | ftp=1.d0/(dlog(max(tp/(TC*TC),1.d-10)))**2 | |
907 | fprob=ftp/fqmax2 | |
908 | rn2=pyr(0) | |
909 | if(fprob.lt.rn2) goto 12 | |
910 | rsk=dsqrt(abs(tp)) | |
911 | if(rsk.gt.ppl) rsk=ppl | |
912 | ||
913 | * calculate radiative energy loss per given scattering with subroutine plfun1 | |
914 | ygl=y*cfac ! mean gluon free path in GeV^{-1} | |
915 | elp=ygl*z ! mimimum radiated energy in LPM regime | |
916 | ej=ppl | |
917 | bb=ej ! maximum radiated energy | |
918 | bbi=max(dsqrt(abs(z)),1.000001d0*elp) | |
919 | aa=min(bb,bbi) ! minimum radiated energy | |
920 | hh=0.00001d0*(bb-aa) | |
921 | REPS=0.01d0 | |
922 | AEPS=1.d-8 | |
923 | CALL SIMPA(aa,bb,hh,REPS,AEPS,plfun1,om,resun,AIH,AIABS) | |
924 | * ! integral over omega for radiative loss | |
925 | call radsear(ermax1,eomin1) | |
926 | ermax=ermax1 | |
927 | eomin=eomin1 | |
928 | 11 resu=eomin*pyr(0)+aa | |
929 | fres=ermax*pyr(0) | |
930 | fres1=plfun1(resu) | |
931 | iraz=iraz+1 | |
932 | if(fres.gt.fres1.and.iraz.lt.100000) goto 11 | |
933 | elr=resu*resun ! energy of radiated gluon | |
934 | ||
935 | * to chancel radiative energy loss (optional case) | |
936 | if(ienglu.eq.2) elr=0.d0 | |
937 | * to chancel collisional energy loss (optional case) | |
938 | if(ienglu.eq.1) rsk=0.d0 | |
939 | ||
940 | * determine the direction of parton moving | |
941 | if(p(i,3).ge.0.d0) then | |
942 | sigp=1.d0 | |
943 | else | |
944 | sigp=-1.d0 | |
945 | end if | |
946 | ||
947 | * calculate new 4-momentum of hard parton | |
948 | phirs=2.d0*pi*pyr(0) | |
949 | epan=epa-rsk*rsk/(2.d0*ep0)-abs(elr) | |
950 | if(epan.lt.0.1d0) then | |
951 | epan=epan+abs(elr) | |
952 | elr=0.d0 | |
953 | if(epan.lt.0.1d0) then | |
954 | rsk=0.d0 | |
955 | epan=epa | |
956 | end if | |
957 | end if | |
958 | pptn=dsqrt(abs(rsk*rsk+(rsk**4)*(1.d0-epa*epa/(ppl*ppl))/ | |
959 | > (4.d0*ep0*ep0)-(rsk**4)*epa/(2.d0*ep0*ppl*ppl)-(rsk**4)/ | |
960 | > (4.d0*ppl*ppl))) | |
961 | ppln=dsqrt(abs(epan*epan-pptn*pptn-p(i,5)**2)) | |
962 | p(i,1)=pptn*dcos(phirs) ! px | |
963 | p(i,2)=pptn*dsin(phirs) ! py | |
964 | p(i,3)=sigp*ppln ! pz | |
965 | p(i,4)=dsqrt(p(i,1)**2+p(i,2)**2+p(i,3)**2+p(i,5)**2) ! E | |
966 | * boost to system of hard parton | |
967 | 222 call pyrobo(i,i,tetr,phir,0.d0,0.d0,0.d0) | |
968 | ||
969 | return | |
970 | end | |
971 | ******************************* END PLJETR ************************** | |
972 | ||
973 | ******************************** PLSEAR *************************** | |
974 | SUBROUTINE PLSEAR (fmax,xmin) | |
975 | * find maximum and 'sufficient minimum' of jet production vertex distribution | |
976 | * xm, fm are outputs. | |
977 | IMPLICIT DOUBLE PRECISION(A-H, O-Z) | |
978 | external plthik | |
979 | common /plpar1/ tau0,T0,TC,sigqq,AW,RA,mvisc,nf | |
980 | save /plpar1/ | |
981 | xmin=3.d0*RA | |
982 | fmax=0.d0 | |
983 | do 10 j=1,1000 | |
984 | x=xmin*(j-1)/999.d0 | |
985 | f=plthik(x) | |
986 | if(f.gt.fmax) then | |
987 | fmax=f | |
988 | end if | |
989 | 10 continue | |
990 | return | |
991 | end | |
992 | ****************************** END PLSEAR ************************** | |
993 | ||
994 | ******************************** RADSEAR *************************** | |
995 | SUBROUTINE RADSEAR (fmax,xmin) | |
996 | * find maximum and 'sufficient minimum' of radiative energy loss distribution | |
997 | * xm, fm are outputs. | |
998 | IMPLICIT DOUBLE PRECISION(A-H, O-Z) | |
999 | external plfun1 | |
1000 | common /radcal/ aa, bb | |
1001 | save /radcal/ | |
1002 | xmin=bb-aa | |
1003 | fmax=0.d0 | |
1004 | do j=1,1000 | |
1005 | x=aa+xmin*(j-1)/999.d0 | |
1006 | f=plfun1(x) | |
1007 | if(f.gt.fmax) then | |
1008 | fmax=f | |
1009 | end if | |
1010 | end do | |
1011 | return | |
1012 | end | |
1013 | ****************************** END RADSEAR ************************** | |
1014 | ||
1015 | ********************************* BIPSEAR *************************** | |
1016 | SUBROUTINE BIPSEAR (fmax,xmin) | |
1017 | * find maximum and 'sufficient minimum' of jet production cross section | |
1018 | * as a function of impact paramater (xm, fm are outputs) | |
1019 | IMPLICIT DOUBLE PRECISION(A-H, O-Z) | |
1020 | external funbip | |
1021 | common /plpar1/ tau0,T0,TC,sigqq,AW,RA,mvisc,nf | |
1022 | save /plpar1/ | |
1023 | xmin=3.d0*RA | |
1024 | fmax=0.d0 | |
1025 | do j=1,1000 | |
1026 | x=xmin*(j-1)/999.d0 | |
1027 | f=funbip(x) | |
1028 | if(f.gt.fmax) then | |
1029 | fmax=f | |
1030 | end if | |
1031 | end do | |
1032 | return | |
1033 | end | |
1034 | ****************************** END RADSEAR ************************** | |
1035 | ||
1036 | **************************** SIMPA ********************************** | |
1037 | SUBROUTINE SIMPA (A1,B1,H1,REPS1,AEPS1,FUNCT,X, | |
1038 | 1 AI,AIH,AIABS) | |
1039 | * calculate intergal of function FUNCT(X) on the interval from A1 to B1 | |
1040 | IMPLICIT DOUBLE PRECISION(A-H, O-Z) | |
1041 | IMPLICIT INTEGER(I-N) | |
1042 | DIMENSION F(7), P(5) | |
1043 | H=dSIGN ( H1, B1-A1 ) | |
1044 | S=dSIGN (1.d0, H ) | |
1045 | A=A1 | |
1046 | B=B1 | |
1047 | AI=0.0d0 | |
1048 | AIH=0.0d0 | |
1049 | AIABS=0.0d0 | |
1050 | P(2)=4.d0 | |
1051 | P(4)=4.d0 | |
1052 | P(3)=2.d0 | |
1053 | P(5)=1.d0 | |
1054 | IF(B-A)1,2,1 | |
1055 | 1 REPS=ABS(REPS1) | |
1056 | AEPS=ABS(AEPS1) | |
1057 | DO 3 K=1,7 | |
1058 | 3 F(K)=10.d16 | |
1059 | X=A | |
1060 | C=0.d0 | |
1061 | F(1)=FUNCT(X)/3.d0 | |
1062 | 4 X0=X | |
1063 | IF( (X0+4.d0*H-B)*S)5,5,6 | |
1064 | 6 H=(B-X0)/4.d0 | |
1065 | IF ( H ) 7,2,7 | |
1066 | 7 DO 8 K=2,7 | |
1067 | 8 F(K)=10.d16 | |
1068 | C=1.d0 | |
1069 | 5 DI2=F (1) | |
1070 | DI3=ABS( F(1) ) | |
1071 | DO 9 K=2,5 | |
1072 | X=X+H | |
1073 | IF((X-B)*S)23,24,24 | |
1074 | 24 X=B | |
1075 | 23 IF(F(K)-10.d16)10,11,10 | |
1076 | 11 F(K)=FUNCT(X)/3.d0 | |
1077 | 10 DI2=DI2+P(K)*F(K) | |
1078 | 9 DI3=DI3+P(K)*ABS(F(K)) | |
1079 | DI1=(F(1)+4.*F(3)+F(5))*2.d0*H | |
1080 | DI2=DI2*H | |
1081 | DI3=DI3*H | |
1082 | IF (REPS) 12,13,12 | |
1083 | 13 IF (AEPS) 12,14,12 | |
1084 | 12 EPS=ABS((AIABS+DI3)*REPS) | |
1085 | IF(EPS-AEPS)15,16,16 | |
1086 | 15 EPS=AEPS | |
1087 | 16 DELTA=ABS(DI2-DI1) | |
1088 | IF(DELTA-EPS)20,21,21 | |
1089 | 20 IF(DELTA-EPS/8.d0)17,14,14 | |
1090 | 17 H=2.d0*H | |
1091 | F(1)=F(5) | |
1092 | F(2)=F(6) | |
1093 | F(3)=F(7) | |
1094 | DO 19 K=4,7 | |
1095 | 19 F(K)=10.d16 | |
1096 | GO TO 18 | |
1097 | 14 F(1)=F(5) | |
1098 | F(3)=F(6) | |
1099 | F(5)=F(7) | |
1100 | F(2)=10.d16 | |
1101 | F(4)=10.d16 | |
1102 | F(6)=10.d16 | |
1103 | F(7)=10.d16 | |
1104 | 18 DI1=DI2+(DI2-DI1)/15.d0 | |
1105 | AI=AI+DI1 | |
1106 | AIH=AIH+DI2 | |
1107 | AIABS=AIABS+DI3 | |
1108 | GO TO 22 | |
1109 | 21 H=H/2.d0 | |
1110 | F(7)=F(5) | |
1111 | F(6)=F(4) | |
1112 | F(5)=F(3) | |
1113 | F(3)=F(2) | |
1114 | F(2)=10.d16 | |
1115 | F(4)=10.d16 | |
1116 | X=X0 | |
1117 | C=0. | |
1118 | GO TO 5 | |
1119 | 22 IF(C)2,4,2 | |
1120 | 2 RETURN | |
1121 | END | |
1122 | ************************* END SIMPA ******************************* | |
1123 | ||
1124 | **************************** SIMPB ********************************** | |
1125 | SUBROUTINE SIMPB (A1,B1,H1,REPS1,AEPS1,FUNCT,X, | |
1126 | 1 AI,AIH,AIABS) | |
1127 | * calculate intergal of function FUNCT(X) on the interval from A1 to B1 | |
1128 | IMPLICIT DOUBLE PRECISION(A-H, O-Z) | |
1129 | IMPLICIT INTEGER(I-N) | |
1130 | DIMENSION F(7), P(5) | |
1131 | H=dSIGN ( H1, B1-A1 ) | |
1132 | S=dSIGN (1.d0, H ) | |
1133 | A=A1 | |
1134 | B=B1 | |
1135 | AI=0.0d0 | |
1136 | AIH=0.0d0 | |
1137 | AIABS=0.0d0 | |
1138 | P(2)=4.d0 | |
1139 | P(4)=4.d0 | |
1140 | P(3)=2.d0 | |
1141 | P(5)=1.d0 | |
1142 | IF(B-A)1,2,1 | |
1143 | 1 REPS=ABS(REPS1) | |
1144 | AEPS=ABS(AEPS1) | |
1145 | DO 3 K=1,7 | |
1146 | 3 F(K)=10.d16 | |
1147 | X=A | |
1148 | C=0.d0 | |
1149 | F(1)=FUNCT(X)/3.d0 | |
1150 | 4 X0=X | |
1151 | IF( (X0+4.d0*H-B)*S)5,5,6 | |
1152 | 6 H=(B-X0)/4.d0 | |
1153 | IF ( H ) 7,2,7 | |
1154 | 7 DO 8 K=2,7 | |
1155 | 8 F(K)=10.d16 | |
1156 | C=1.d0 | |
1157 | 5 DI2=F (1) | |
1158 | DI3=ABS( F(1) ) | |
1159 | DO 9 K=2,5 | |
1160 | X=X+H | |
1161 | IF((X-B)*S)23,24,24 | |
1162 | 24 X=B | |
1163 | 23 IF(F(K)-10.d16)10,11,10 | |
1164 | 11 F(K)=FUNCT(X)/3.d0 | |
1165 | 10 DI2=DI2+P(K)*F(K) | |
1166 | 9 DI3=DI3+P(K)*ABS(F(K)) | |
1167 | DI1=(F(1)+4.*F(3)+F(5))*2.d0*H | |
1168 | DI2=DI2*H | |
1169 | DI3=DI3*H | |
1170 | IF (REPS) 12,13,12 | |
1171 | 13 IF (AEPS) 12,14,12 | |
1172 | 12 EPS=ABS((AIABS+DI3)*REPS) | |
1173 | IF(EPS-AEPS)15,16,16 | |
1174 | 15 EPS=AEPS | |
1175 | 16 DELTA=ABS(DI2-DI1) | |
1176 | IF(DELTA-EPS)20,21,21 | |
1177 | 20 IF(DELTA-EPS/8.d0)17,14,14 | |
1178 | 17 H=2.d0*H | |
1179 | F(1)=F(5) | |
1180 | F(2)=F(6) | |
1181 | F(3)=F(7) | |
1182 | DO 19 K=4,7 | |
1183 | 19 F(K)=10.d16 | |
1184 | GO TO 18 | |
1185 | 14 F(1)=F(5) | |
1186 | F(3)=F(6) | |
1187 | F(5)=F(7) | |
1188 | F(2)=10.d16 | |
1189 | F(4)=10.d16 | |
1190 | F(6)=10.d16 | |
1191 | F(7)=10.d16 | |
1192 | 18 DI1=DI2+(DI2-DI1)/15.d0 | |
1193 | AI=AI+DI1 | |
1194 | AIH=AIH+DI2 | |
1195 | AIABS=AIABS+DI3 | |
1196 | GO TO 22 | |
1197 | 21 H=H/2.d0 | |
1198 | F(7)=F(5) | |
1199 | F(6)=F(4) | |
1200 | F(5)=F(3) | |
1201 | F(3)=F(2) | |
1202 | F(2)=10.d16 | |
1203 | F(4)=10.d16 | |
1204 | X=X0 | |
1205 | C=0. | |
1206 | GO TO 5 | |
1207 | 22 IF(C)2,4,2 | |
1208 | 2 RETURN | |
1209 | END | |
1210 | ************************* END SIMPB ******************************* | |
1211 | ||
1212 | ************************* PARINV ********************************** | |
1213 | SUBROUTINE PARINV(X,A,F,N,R) | |
1214 | * gives interpolation of function F(X) with arrays A(N) and F(N) | |
1215 | IMPLICIT DOUBLE PRECISION(A-H, O-Z) | |
1216 | IMPLICIT INTEGER(I-N) | |
1217 | DIMENSION A(N),F(N) | |
1218 | IF(X.LT.A(1))GO TO 11 | |
1219 | IF(X.GT.A(N))GO TO 4 | |
1220 | K1=1 | |
1221 | K2=N | |
1222 | 2 K3=K2-K1 | |
1223 | IF(K3.LE.1)GO TO 6 | |
1224 | K3=K1+K3/2 | |
1225 | IF(A(K3)-X) 7,8,9 | |
1226 | 7 K1=K3 | |
1227 | GOTO2 | |
1228 | 9 K2=K3 | |
1229 | GOTO2 | |
1230 | 8 P=F(K3) | |
1231 | RETURN | |
1232 | 3 B1=A(K1) | |
1233 | B2=A(K1+1) | |
1234 | B3=A(K1+2) | |
1235 | B4=F(K1) | |
1236 | B5=F(K1+1) | |
1237 | B6=F(K1+2) | |
1238 | R=B4*((X-B2)*(X-B3))/((B1-B2)*(B1-B3))+B5*((X-B1)*(X-B3))/ | |
1239 | 1 ((B2-B1)*(B2-B3))+B6*((X-B1)*(X-B2))/((B3-B1)*(B3-B2)) | |
1240 | RETURN | |
1241 | 6 IF(K2.NE.N)GO TO 3 | |
1242 | K1=N-2 | |
1243 | GOTO3 | |
1244 | 4 C=ABS(X-A(N)) | |
1245 | IF(C.LT.0.1d-7) GO TO 5 | |
1246 | K1=N-2 | |
1247 | 13 CONTINUE | |
1248 | C13 PRINT 41,X | |
1249 | C41 FORMAT(25H X IS OUT OF THE INTERVAL,3H X=,F15.9) | |
1250 | GO TO 3 | |
1251 | 5 R=F(N) | |
1252 | RETURN | |
1253 | 11 C=ABS(X-A(1)) | |
1254 | IF(C.LT.0.1d-7) GO TO 12 | |
1255 | K1=1 | |
1256 | GOTO 13 | |
1257 | 12 R=F(1) | |
1258 | RETURN | |
1259 | END | |
1260 | C************************** END PARINV ************************************* | |
1261 | ||
1262 | * quark-quark scattering differential cross section | |
1263 | double precision FUNCTION PLSIGH(Z) | |
1264 | IMPLICIT DOUBLE PRECISION(A-H, O-Z) | |
1265 | common /plpar1/ tau0,T0,TC,sigqq,AW,RA,mvisc,nf | |
1266 | save /plpar1/ | |
1267 | pi=3.14159d0 | |
1268 | beta=(33.d0-2.d0*nf)/(12.d0*pi) | |
1269 | alfs=1.d0/(beta*dlog(max(1.d-10,z/(TC*TC)))) | |
1270 | PLSIGH=8.d0*pi*alfs*alfs/(9.d0*z*z) | |
1271 | return | |
1272 | end | |
1273 | ||
1274 | * differential radiated gluon spectrum in BDMS model | |
1275 | double precision FUNCTION PLFUN1(or) | |
1276 | IMPLICIT DOUBLE PRECISION(A-H, O-Z) | |
1277 | common /plpar1/ tau0,T0,TC,sigqq,AW,RA,mvisc,nf | |
1278 | common /pljdat/ ej, z, ygl, alfs, um, epa | |
1279 | common /pleave/ taul, temlev | |
1280 | common /factor/ cfac, kf | |
1281 | save /plpar1/,/pljdat/,/pleave/,/factor/ | |
1282 | pi=3.14159d0 | |
1283 | x=min((1.d0-ygl*z/or),or/ej) | |
1284 | if(x.le.0.d0) x=0.d0 | |
1285 | if(x.ge.1.d0) x=0.9999d0 | |
1286 | if(kf.eq.21) then | |
1287 | if(x.ge.0.5d0) x=1.d0-x | |
1288 | spinf=0.5d0*(1.+(1.d0-x)**4+x**4)/(1.d0-x) | |
1289 | else | |
1290 | spinf=1.d0-x+0.5d0*x*x | |
1291 | end if | |
1292 | ak=ygl*z/(or*(1.d0-x)) | |
1293 | al=taul*5.06d0 | |
1294 | uu=0.5d0*al*dsqrt(abs(0.5d0*(1.d0-x+cfac*x*x)*ak* | |
1295 | > dlog(max(16.d0/ak,1.d-10))))/ygl | |
1296 | * if quark production outside the QGP then | |
1297 | * arg=(((dsin(uu)*cosh(uu))**2)+((dcos(uu)*sinh(uu))**2))/(2.d0*uu*uu); | |
1298 | * here quark production inside the QGP | |
1299 | arg=((dcos(uu)*cosh(uu))**2)+((dsin(uu)*sinh(uu))**2) | |
1300 | gl1=(ygl/(cfac*z))**0.3333333d0 | |
1301 | gl2=(um/epa)**1.333333d0 | |
1302 | dc=1.d0/((1.d0+((gl1*gl2*or)**1.5d0))**2) ! massive parton | |
1303 | c dc=1.d0 !massless parton | |
1304 | plfun1=dc*3.d0*alfs*ygl*dlog(max(arg,1.d-20))*spinf/(pi*al*or) | |
1305 | return | |
1306 | end | |
1307 | ||
1308 | * angular distribution of emitted gluons | |
1309 | double precision function gluang(x) | |
1310 | IMPLICIT DOUBLE PRECISION(A-H, O-Z) | |
1311 | s=0.0863d0 | |
1312 | gluang=x*dexp(-1.d0*(x-s)*(x-s)/(2.d0*s*s)) | |
1313 | return | |
1314 | end | |
1315 | ||
1316 | * temperature-dependence of parton-plasma integral cross section | |
1317 | double precision FUNCTION PLS(X) | |
1318 | IMPLICIT DOUBLE PRECISION(A-H, O-Z) | |
1319 | external plsigh | |
1320 | common /plpar1/ tau0,T0,TC,sigqq,AW,RA,mvisc,nf | |
1321 | common /plpar2/ pln0,taupl,tauh,sigpl,sigh,sigplh,sigqqh,rg,rgn | |
1322 | common /plen/ epartc, um | |
1323 | save /plpar1/,/plpar2/,/plen/ | |
1324 | t=X | |
1325 | pi=3.14159d0 | |
1326 | bub=4.d0*t/TC | |
1327 | alf=6.d0*pi/((33.d0-2.d0*nf)*dlog(max(bub,1.d-10))) | |
1328 | ZZ0=4.d0*t*t*pi*alf*(1.d0+nf/6.d0) | |
1329 | scm=4.d0*t*epartc+um*um+4.d0*t*t | |
1330 | ZZ1=max((scm-((um+2.d0*t)**2))*(scm-((um-2.d0*t)**2))/scm,ZZ0) | |
1331 | HH1=0.01d0*ZZ1 | |
1332 | REPS=0.01d0 | |
1333 | AEPS=1.d-8 | |
1334 | CALL SIMPA(ZZ0,ZZ1,HH1,REPS,AEPS,plsigh,ZZ,RESS,AIH,AIABS) | |
1335 | PLS=0.39d0*2.25d0*2.25d0*RESS*(16.d0+4.d0*nf)/(16.d0+9.d0*nf) | |
1336 | return | |
1337 | end | |
1338 | ||
1339 | * temperature-dependence of QGP viscosity (if mvisc=1,2) | |
1340 | double precision FUNCTION PLVISC(X) | |
1341 | IMPLICIT DOUBLE PRECISION(A-H, O-Z) | |
1342 | common /plpar1/ tau0,T0,TC,sigqq,AW,RA,mvisc,nf | |
1343 | save /plpar1/ | |
1344 | pi=3.14159d0 | |
1345 | T=X | |
1346 | TC1=5.06d0*TC | |
1347 | if(X.le.TC1) T=TC1 | |
1348 | if(mvisc.eq.0) then | |
1349 | c=0.d0 | |
1350 | elseif(mvisc.eq.1) then | |
1351 | a=3.4d0*(1.d0+0.12d0*(2.d0*nf+1.d0)) | |
1352 | b=15.d0*(1.d0+0.06d0*nf) | |
1353 | c=4.d0*pi*pi*(10.5d0*nf/a+16.d0/b)/675.d0 | |
1354 | else | |
1355 | c=(1.7d0*nf+1.d0)*0.342d0/(1.d0+nf/6.d0) | |
1356 | end if | |
1357 | bub=4.d0*T/TC1 | |
1358 | alf=6.d0*pi/((33.d0-2.d0*nf)*dlog(max(bub,1.d-10))) | |
1359 | alf1=1.d0/alf | |
1360 | PLVISC=c*(T**3)/(alf*alf*dlog(max(1.d-10,alf1))) | |
1361 | return | |
1362 | end | |
1363 | ||
1364 | * space-time dependence of QGP number density | |
1365 | double precision FUNCTION PLN(X,r1,r2,y) | |
1366 | IMPLICIT DOUBLE PRECISION(A-H, O-Z) | |
1367 | external pythik | |
1368 | common /plpar1/ tau0,T0,TC,sigqq,AW,RA,mvisc,nf | |
1369 | common /plpar2/ pln0,taupl,tauh,sigpl,sigh,sigplh,sigqqh,rg,rgn | |
1370 | common /plevol/ taup(5000),temp(5000),denp(5000),enep(5000) | |
1371 | common /pythic/ PBAB(110),PTAB(110),PTAAB(110) | |
1372 | save /plpar1/,/plpar2/,/plevol/,/pythic/ | |
1373 | pi=3.14159d0 | |
1374 | t=X | |
1375 | if(t.lt.taupl) then | |
1376 | call parinv(t,taup,denp,5000,res) | |
1377 | else | |
1378 | res=1.2d0*(16.d0+9.d0*nf)*((5.06d0*TC)**3)/(pi*pi) | |
1379 | end if | |
1380 | res=res*(pythik(r1)*pythik(r2)*pi*RA*RA/PTAAB(1))**0.75d0 | |
1381 | res=res*dexp(-1.d0*y*y/24.5d0) | |
1382 | PLN=max(1.d-8,res) | |
1383 | return | |
1384 | end | |
1385 | ||
1386 | * space-time dependence of QGP temperature | |
1387 | double precision FUNCTION PLT(X,r1,r2,y) | |
1388 | IMPLICIT DOUBLE PRECISION(A-H, O-Z) | |
1389 | common /plpar1/ tau0,T0,TC,sigqq,AW,RA,mvisc,nf | |
1390 | common /plpar2/ pln0,taupl,tauh,sigpl,sigh,sigplh,sigqqh,rg,rgn | |
1391 | common /plevol/ taup(5000),temp(5000),denp(5000),enep(5000) | |
1392 | common /pythic/ PBAB(110),PTAB(110),PTAAB(110) | |
1393 | save /plpar1/,/plpar2/,/plevol/,/pythic/ | |
1394 | pi=3.14159d0 | |
1395 | t=X | |
1396 | if(t.lt.taupl) then | |
1397 | call parinv(t,taup,temp,5000,res) | |
1398 | else | |
1399 | res=TC | |
1400 | end if | |
1401 | res=res*(pythik(r1)*pythik(r2)*pi*RA*RA/PTAAB(1))**0.25d0 | |
1402 | res=res*(dexp(-1.d0*y*y/24.5d0))**0.333333d0 | |
1403 | PLT=max(1.d-8,res) | |
1404 | return | |
1405 | end | |
1406 | ||
1407 | * impact parameter dependence of jet production cross section | |
1408 | double precision function funbip(x) | |
1409 | IMPLICIT DOUBLE PRECISION(A-H, O-Z) | |
1410 | external ftaa | |
1411 | common /pyint7/ sigt(0:6,0:6,0:5) | |
1412 | save /pyint7/ | |
1413 | br=x | |
1414 | sigin=sigt(0,0,0)-sigt(0,0,1) | |
1415 | taa=ftaa(br) | |
1416 | funbip=taa*br*(1.d0-dexp(-0.1d0*taa*sigin)) | |
1417 | return | |
1418 | end | |
1419 | ||
1420 | * distribution over jet production vertex position | |
1421 | double precision FUNCTION plthik(X) | |
1422 | IMPLICIT DOUBLE PRECISION(A-H, O-Z) | |
1423 | external pythik | |
1424 | common /parimp/ b1, psib1, rb1, rb2, noquen | |
1425 | save /parimp/ | |
1426 | bu=X | |
1427 | r12=dsqrt(abs(bu*bu+b1*b1/4.d0+bu*b1*dcos(psib1))) | |
1428 | r22=dsqrt(abs(bu*bu+b1*b1/4.d0-bu*b1*dcos(psib1))) | |
1429 | PLTHIK=bu*pythik(r12)*pythik(r22) | |
1430 | return | |
1431 | end | |
1432 | ||
1433 | * nuclear overlap function at impact parameter b | |
1434 | double precision function ftaa(r) | |
1435 | IMPLICIT DOUBLE PRECISION(A-H, O-Z) | |
1436 | common /pythic/ PBAB(110),PTAB(110),PTAAB(110) | |
1437 | save /pythic/ | |
1438 | call parinv(r,PBAB,PTAAB,110,RES) | |
1439 | ftaa=RES | |
1440 | return | |
1441 | end | |
1442 | * | |
1443 | double precision function PFUNC1(x) | |
1444 | IMPLICIT DOUBLE PRECISION(A-H, O-Z) | |
1445 | external PFUNC2 | |
1446 | common /pynup1/ bp,xx | |
1447 | common /plpar1/ tau0,T0,TC,sigqq,AW,RA,mvisc,nf | |
1448 | save /plpar1/ | |
1449 | xx=x | |
1450 | EPS=0.05d0 | |
1451 | A=0.d0 | |
1452 | B=3.d0*RA | |
1453 | H=0.01d0*(B-A) | |
1454 | CALL SIMPB(A,B,H,EPS,1.d-8,PFUNC2,Y,RES,AIH,AIABS) | |
1455 | PFUNC1=RES | |
1456 | return | |
1457 | end | |
1458 | * | |
1459 | double precision function PFUNC2(y) | |
1460 | IMPLICIT DOUBLE PRECISION(A-H, O-Z) | |
1461 | external pythik | |
1462 | common /pynup1/ bp,x | |
1463 | r1=sqrt(abs(y*y+bp*bp/4.+y*bp*cos(x))) | |
1464 | r2=sqrt(abs(y*y+bp*bp/4.-y*bp*cos(x))) | |
1465 | PFUNC2=y*pythik(r1)*pythik(r2) | |
1466 | return | |
1467 | end | |
1468 | ||
1469 | * nuclear thickness function | |
1470 | double precision function pythik(r) | |
1471 | IMPLICIT DOUBLE PRECISION(A-H, O-Z) | |
1472 | common /pythic/ PBAB(110),PTAB(110),PTAAB(110) | |
1473 | save /pythic/ | |
1474 | call parinv(r,PBAB,PTAB,110,RES) | |
1475 | pythik=RES | |
1476 | return | |
1477 | end | |
1478 | ||
1479 | * Wood-Saxon nucleon distrubution | |
1480 | double precision function prhoaa(z) | |
1481 | IMPLICIT DOUBLE PRECISION(A-H, O-Z) | |
1482 | common /plpar1/ tau0,T0,TC,sigqq,AW,RA,mvisc,nf | |
1483 | common /pygeom/ BC | |
1484 | save /plpar1/,/pygeom/ | |
1485 | pi=3.14159d0 | |
1486 | df=0.54d0 | |
1487 | r=sqrt(bc*bc+z*z) | |
1488 | rho0=3.d0/(4.d0*pi*RA**3)/(1.d0+(pi*df/RA)**2) | |
1489 | prhoaa=rho0/(1.d0+exp((r-RA)/df)) | |
1490 | return | |
1491 | end | |
1492 | ||
1493 | * function to generate gauss distribution | |
1494 | double precision function gauss(x0,sig) | |
1495 | IMPLICIT DOUBLE PRECISION(A-H, O-Z) | |
1496 | 41 u1=pyr(0) | |
1497 | u2=pyr(0) | |
1498 | v1=2.d0*u1-1.d0 | |
1499 | v2=2.d0*u2-1.d0 | |
1500 | s=v1**2+v2**2 | |
1501 | if(s.gt.1) go to 41 | |
1502 | gauss=v1*dsqrt(-2.d0*dlog(s)/s)*sig+x0 | |
1503 | return | |
1504 | end | |
1505 | ************************************************************************** |