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