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