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