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