1 *----------------------------------------------------------------------
5 * Author : Igor Lokhtin (Igor.Lokhtin@cern.ch)
6 * Version : PYQUEN1_3.f
7 * Last revision : 03-OCT-2007
9 *======================================================================
11 * Description : Event generator for simulation of parton rescattering
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)
16 * Reference: I.P. Lokhtin, A.M. Snigirev, Eur. Phys. J. C 46 (2006) 211
18 *======================================================================
20 SUBROUTINE PYQUEN(A,ifb,bfix)
21 IMPLICIT DOUBLE PRECISION(A-H, O-Z)
23 INTEGER PYK,PYCHGE,PYCOMP
25 external pyp,pyr,pyk,pyjoin,pyshow
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
31 common /plglur/ glur(1000,4),kglu(1000,6),nrg,nrgm
32 common /plquar/ pqua(1000,5),kqua(1000,5),nrq
33 common /parimp/ b1,psib1,rb1,rb2,noquen
34 common /pyqpar/ T0u,tau0u,nfu,ienglu,ianglu
37 save /pyjets/, /pydat1/, /pysubs/, /plglur/, /plquar/
38 dimension ijoik(2),ijoin(1000),nis(500),nss(500),nas(500),nus(500)
40 * set initial event paramters
42 RA=1.15d0*AW**0.333333d0 ! nucleus radius in fm
44 mvisc=0 ! flag of QGP viscosity (off here)
45 TC=0.2d0 ! crutical temperature
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.
58 * avoid stopping run if pythia does not conserve energy due to collisional loss
61 * generate impact parameter of A-A collision with jet production
64 write(6,*) 'Impact parameter less than zero!'
68 write(6,*) 'Impact parameter larger than two nuclear radius!'
73 call bipsear(fmax1,xmin1)
84 * calculate initial QGP temperature as function of centrality
86 sb=RA*RA*(pi-2.d0*dasin(0.5d0*b1/RA))-b1*dsqrt(abs(RA*RA-
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,
92 rtaa=rtaa0*(1.d0-br*(1.d0+(1.d0-0.25d0*br)*dlog(1.d0/br)+
94 T00=T000*((rtaa*sb0)/(rtaa0*sb))**0.25d0
95 T0=T00*(AW/207.d0)**0.166667d0
97 * generate single event with partonic energy loss
100 if(b1.le.1.85d0*RA) then
105 * reset all in-vacuum radiated guark 4-momenta and codes to zero
114 * generate final state shower in vacuum if it was excluded before
115 nrgm=nrg ! fix number of in-medium emitted gluons
122 if(mstj(41).ne.0) goto 5
127 ip1=i ! first hard parton (line ip1)
128 kfh1=k(i,1) ! status code of first hard parton
129 qmax1=pyp(i,10) ! transverse momentum of first hard parton
132 ip2=i ! second hard parton (line ip2)
133 kfh2=k(i,1) ! status code of second hard parton
134 qmax2=pyp(i,10) ! transverse momentum of second hard parton
139 call pyshow(ip1,0,qmax1) ! vacuum showering for first hard parton
142 call pyshow(ip2,0,qmax2) ! vacuum showering for second hard parton
147 * find two leading partons after showering
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
155 if (k(i,1).eq.14) goto 3
156 if(i.ge.ip002.and.ip002.gt.0) then
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
162 elseif(ip001.gt.0) then
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
172 * replace two hard partons by two leading partons in original event record
179 * fix first/last daughter for moving entry
181 if(k(jgl,4).eq.ip01) k(jgl,4)=ip1
182 if(k(jgl,5).eq.ip01) k(jgl,5)=ip1
192 * fix first/last daughter for moving entry
194 if(k(jgl,4).eq.ip02) k(jgl,4)=ip2
195 if(k(jgl,5).eq.ip02) k(jgl,5)=ip2
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
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
214 if(i.ge.ip002.and.ip002.gt.0) then
222 if(ish.ge.ishm.or.nur.le.2) goto 6 ! adding gluons in 'plglur' arrays
224 kglu(nur,j)=kglu(nur-1,j)
227 glur(nur,j)=glur(nur-1,j)
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
242 do j=1,5 ! remove partons from event list
252 * stop generate event if there are no additional gluons
255 * define number of stirngs (ns) and number of entries in strings before
256 * in-medium radiation (nis(ns))
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
274 nes=nes+nis(ns+1) ! nes - total number of entries
277 elseif(ks.ne.2.and.ksp.ne.2.and.ns.gt.0) then
278 i1=i1+1 ! last i1 lines not included in strings
281 i0=n-nes-i1 ! first i0 lines not included in strings
286 * move fragmented particles in bottom of event list
290 if(ks.ne.2.and.ksp.ne.2) then
297 * fix first/last daughter for moving entry
299 if(k(jgl,4).eq.i) k(jgl,4)=n
300 if(k(jgl,5).eq.i) k(jgl,5)=n
309 * fix first/last daughter for moving entry
311 if(k(jgl,4).eq.in) k(jgl,4)=in-1
312 if(k(jgl,5).eq.in) k(jgl,5)=in-1
318 if(ku.gt.i) kglu(ip,6)=ku-1
326 * define number of additional entries in strings, nas(ns)
329 if(kas.le.nss(1)) then
333 if(kas.le.nss(j).and.kas.gt.nss(j-1))
344 * add emitted gluons in event list
354 * fix first/last daughter for moving entries
356 if(k(jgl,4).eq.i) k(jgl,4)=is
357 if(k(jgl,5).eq.i) k(jgl,5)=is
362 do i=nss(ia+1)-1,nss(ia),-1
369 * fix first/last daughter for moving entries
371 if(k(jgl,4).eq.i) k(jgl,4)=is
372 if(k(jgl,5).eq.i) k(jgl,5)=is
383 if(i.le.nus(in).and.i.gt.nus(in-1))
395 p(ia,1)=ptg*dcos(phig)
396 p(ia,2)=ptg*dsin(phig)
397 p(ia,3)=dsqrt(abs(eg*eg-ptg*ptg))
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)
403 * rearrange partons to form strings in event list
415 ijoin(j)=nss(i-1)+nus(i-1)+j
418 call pyjoin(njoin,ijoin)
421 * add in-vacuum emitted quark pairs
431 4 ktest=k(n-1,2)+kqua(in,2)
432 if(ktest.eq.0.or.in.eq.nrq) goto 8
442 kqua(in,j)=kqua(i+1,j)
443 pqua(in,j)=pqua(i+1,j)
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
470 common /plpar1/ tau0,T0,TC,sigqq,AW,RA,mvisc,nf
471 common /plpar2/ pln0,taupl,tauh,sigpl,sigh,sigplh,sigqqh,rg,rgn
472 common /plevol/ taup(5000),temp(5000),denp(5000),enep(5000)
478 * set number degrees of freedom in QGP
480 rg=(16.d0+10.5d0*nf)/hgd
481 rgn=(16.d0+9.d0*nf)/hgd
483 * set 'fiction' sigma for parton rescattering in QGP
485 sigpl=2.25d0*2.25d0*sigqq*(16.d0+4.d0*nf)/(16.d0+9.d0*nf)
487 * set intial plasma temperature, density and energy density in perfect
488 * (if mvisc=0) or viscous (mvisc=1,2) QGP with PLVISC subroitine
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
494 pln0=(16.d0+9.d0*nf)*1.2d0*(T01**3)/pi2
495 ened0=pi2*(16.d0+10.5d0*nf)*(T01**4)/30.d0
497 tau=tau0 ! proper time
499 den=pln0 ! number density
500 ened=ened0 ! energy density
502 * create array of parameters to configurate QGP time evolution
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
510 T1=(30.d0*ened1/((16.d0+10.5d0*nf)*pi2))**0.25d0
512 ened=hh*(1.3333d0*plvisc(T1)/(tau1*tau1)-1.3333d0
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
518 if(TPR.gt.TC1.and.T.le.TC1) taupl=tau-0.5d0*hh ! QGP lifetime taupl
520 tauh=taupl*rg ! mixed phase lifetime
524 ******************************** END PLINIT **************************
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
533 external plthik, pln, plt, pls, gauss, gluang
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)
539 common /pyqpar/ T0u,tau0u,nfu,ienglu,ianglu
540 common /plglur/ glur(1000,4),kglu(1000,6),nrg,nrgm
541 common /factor/ cfac, kf
542 common /pleave/ taul, temlev
543 common /parimp/ b1, psib1, rb1, rb2, noquen
544 common /plen/ epartc, um
545 common /plos/ elr,rsk
546 common /numje1/ nuj1, nuj2
547 save /pyjets/, /plglur/
551 * find minimum of nuclear thikness function with subroutine plsear
552 psib1=pi*(2.d0*pyr(0)-1.d0)
553 call plsear (fmax1,xmin1)
557 * generate vertex of jet production
563 if(ff1.gt.f.and.iv.le.100000) goto 1
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)))
569 if(noquen.eq.1) goto 7
571 * find maximum of angular spectrum of radiated gluons with subroutine gluang
573 temax=0.5d0*(1.d0+dsqrt(5.d0))*0.0863d0
576 * reset all radiated gluon 4-momenta and codes to zero -------------------
587 * generate changing 4-momentum of partons due to rescattering and energy loss
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
591 do 2 ip=nuj1,nuj2 ! cycle on travelling partons
594 ks=k(ip,1) ! parton status code
595 kf=k(ip,2) ! parton identificator
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).
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.
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
609 stetr=dsin(tetr) ! parton sin(theta)
610 if(abs(stetr).le.1.d-05) then
611 if(stetr.ge.0.d0) then
622 cfac=1.d0 ! for gluon
624 cfac=0.44444444d0 ! for quark
627 * boost from laboratory system to system of hard parton
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
633 if(psib1.lt.0.d0) bet=-1.d0*bet
635 if(phip.gt.pi) phip=phip-2.d0*pi
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)
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
645 if(psib1.le.0.d0) phin=-1.d0*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
655 * start parton rescattering in QGP with proper time iterations
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
662 if(tau.ge.taul.or.tfs.le.TC) goto 100 ! escape if tau>taul or >taupl
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)
671 if(irasf.gt.100000) goto 100
672 if(ran.lt.prob) goto 3
675 pass=50.6d0/(pln(tau)*sigtr)
678 call pljetr(tau,pass,pltp,ipar,epart)
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
684 if(abs(elr).gt.0.1d0.and.ip.gt.8) then
685 * generate the angle of emitted gluon
690 if(fte1.gt.fte2) goto 6
692 elseif (ianglu.eq.1) then
693 tgl=((0.5d0*pi*epartc)**pyr(0))/epartc
697 pgl=pi*(2.d0*pyr(0)-1.d0)
699 pxgl=abs(elr)*stetr*(dcos(phir)*dcos(tgl)-
700 > dsin(phir)*dsin(tgl)*dsin(pgl))
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))
705 psgl=dsqrt(abs(ptgl*ptgl+pzgl*pzgl))
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))
712 glur1=abs(elr) ! energy
713 glur3=datan(dpgl) ! phi
714 if(pxgl.lt.0.d0) then
715 if(pygl.ge.0.d0) then
721 glur4=0.5d0*dlog(max(1.d-9,(psgl+pzgl)/(psgl-pzgl))) ! eta
722 glur2=glur1/dcosh(glur4) ! pt
724 * put in event list radiated gluons with pt > 0.2 GeV only
725 if(glur2.ge.0.2d0) then
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
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
742 write(6,*) 'Warning! Number of emitted gluons is too large!'
745 * set parton "thermalization" if pt<T
746 if(abs(p(ip,3)).gt.pltp3) goto 3
748 if(p(ip,3).ge.0.d0) then
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))
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)
769 * boost to laboratory system
770 100 call pyrobo(ip,ip,tetr,phir,0.d0,0.d0,0.d0)
779 ******************************* END PLEVNT *************************
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
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)
792 common /pyqpar/ T0u,tau0u,nfu,ienglu,ianglu
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
802 tauu=x ! redenote temerature tauu=x
803 i=ip ! redenote parton number i=ip
807 * boost to system of comoving plasma constituent
808 phir=pyp(i,15) ! parton phi
809 tetr=pyp(i,13) ! parton theta
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)
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
821 if(ppl.lt.3.d0) goto 222 ! no energy loss if pz<3 GeV/c
823 * generation hard parton-plasma scattering with momentum transfer rsk
824 221 ep0=-1.*tauu*(dlog(max(1.d-10,pyr(0)))+dlog(max(1.d-10,
825 > pyr(0)))+dlog(max(1.d-10,pyr(0)))) ! energy of 'thermal' parton
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
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)))
836 phmax2=max(phmin2,qm2)
837 fqmax2=1.d0/(dlog(max(phmin2/(TC*TC),1.d-10)))**2
839 tp=1.d0/(rn1/phmax2+(1.d0-rn1)/phmin2)
840 ftp=1.d0/(dlog(max(tp/(TC*TC),1.d-10)))**2
843 if(fprob.lt.rn2) goto 12
845 if(rsk.gt.ppl) rsk=ppl
847 * calculate radiative energy loss per given scattering with subroutine plfun1
848 ygl=y*cfac ! mean gluon free path in GeV^{-1}
849 elp=ygl*z ! mimimum radiated energy in LPM regime
851 bb=ej ! maximum radiated energy
852 bbi=max(dsqrt(abs(z)),1.000001d0*elp)
853 aa=min(bb,bbi) ! minimum radiated energy
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)
862 11 resu=eomin*pyr(0)+aa
866 if(fres.gt.fres1.and.iraz.lt.100000) goto 11
867 elr=resu*resun ! energy of radiated gluon
869 * to chancel radiative energy loss (optional case)
870 if(ienglu.eq.2) elr=0.d0
871 * to chancel collisional energy loss (optional case)
872 if(ienglu.eq.1) rsk=0.d0
874 * determine the direction of parton moving
875 if(p(i,3).ge.0.d0) then
881 * calculate new 4-momentum of hard parton
883 epan=epa-rsk*rsk/(2.d0*ep0)-abs(elr)
884 if(epan.lt.0.1d0) then
887 if(epan.lt.0.1d0) then
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)/
895 ppln=dsqrt(abs(epan*epan-pptn*pptn-p(i,5)**2))
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
900 * boost to system of hard parton
901 222 call pyrobo(i,i,tetr,phir,0.d0,0.d0,0.d0)
905 ******************************* END PLJETR **************************
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
915 common /parimp/ b1, psib1, rb1, rb2, noquen
916 common /plpar1/ tau0,T0,TC,sigqq,AW,RA,mvisc,nf
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
934 ****************************** END PLSEAR **************************
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
944 common /radcal/ aa, bb
949 x=aa+xmin*(j-1)/999.d0
958 ****************************** END RADSEAR **************************
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
968 common /plpar1/ tau0,T0,TC,sigqq,AW,RA,mvisc,nf
982 ****************************** END RADSEAR **************************
984 **************************** SIMPA **********************************
985 SUBROUTINE SIMPA (A1,B1,H1,REPS1,AEPS1,FUNCT,X,
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)
991 H=dSIGN ( H1, B1-A1 )
1011 IF( (X0+4.d0*H-B)*S)5,5,6
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
1031 13 IF (AEPS) 12,14,12
1032 12 EPS=ABS((AIABS+DI3)*REPS)
1033 IF(EPS-AEPS)15,16,16
1035 16 DELTA=ABS(DI2-DI1)
1036 IF(DELTA-EPS)20,21,21
1037 20 IF(DELTA-EPS/8.d0)17,14,14
1052 18 DI1=DI2+(DI2-DI1)/15.d0
1070 ************************* END SIMPA *******************************
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)
1078 IF(X.LT.A(1))GO TO 11
1079 IF(X.GT.A(N))GO TO 4
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))
1101 6 IF(K2.NE.N)GO TO 3
1105 IF(C.LT.0.1d-7) GO TO 5
1109 C41 FORMAT(25H X IS OUT OF THE INTERVAL,3H X=,F15.9)
1114 IF(C.LT.0.1d-7) GO TO 12
1120 C************************** END PARINV *************************************
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
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)
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
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
1149 if(x.ge.0.5d0) x=1.-x
1150 spinf=0.5*(1.+(1.d0-x)**4+x**4)/(1.d0-x)
1152 spinf=1.d0-x+0.5d0*x*x
1154 ak=ygl*z/(or*(1.d0-x))
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
1165 c dc=1.d0 !massless parton
1166 plfun1=dc*3.d0*alfs*ygl*dlog(max(arg,1.d-20))*spinf/(pi*al*or)
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
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
1187 c=(1.7d0*nf+1.d0)*0.342d0/(1.d0+nf/6.d0)
1190 alf=6.d0*pi/((33.d0-2.d0*nf)*dlog(max(bub,1.d-10)))
1192 PLVISC=c*(T**3)/(alf*alf*dlog(max(1.d-10,alf1)))
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)
1205 pi2=3.14159d0*3.14159d0
1208 call parinv(t,taup,denp,5000,res)
1210 res=1.2d0*(16.d0+9.d0*nf)*((5.06d0*TC)**3)/pi2
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)
1227 call parinv(t,taup,temp,5000,res)
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
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
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)
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)
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
1264 common /parimp/ b1, psib1, rb1, rb2, noquen
1265 common /plpar1/ tau0,T0,TC,sigqq,AW,RA,mvisc,nf
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
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)
1283 gauss=v1*dsqrt(-2.d0*dlog(s)/s)*sig+x0
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)
1292 gluang=x*dexp(-1.d0*(x-s)*(x-s)/(2.d0*s*s))
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/
1310 bipr(i)=bip(i)*RA/6.8d0
1312 call parinv (bu,bipr,pjet,15,res)
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
1324 ftaa=(1.d0-br*x*x/a)*dlog(1.d0+x*x*(1.d0-br))/(a*a)
1327 **************************************************************************