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