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