*
* Filename : PYQUEN.F
*
-* First version created: 20-AUG-1997 Author : Igor Lokhtin
-* Last revision : 23-SEP-2004
+* Author : Igor Lokhtin (Igor.Lokhtin@cern.ch)
+* Version : PYQUEN1_3.f
+* Last revision : 03-OCT-2007
*
*======================================================================
*
* Description : Event generator for simulation of parton rescattering
-* and energy loss in quark-gluon plasma created in heavy
-* ion AA collisons at LHC
-* (implemented as modification of standard pythia jet event)
+* and energy loss in expanding quark-gluon plasma created
+* in ultrarelativistic heavy ion AA collisons
+* (implemented as modification of standard Pythia jet event)
*
-* Method : I.P.Lokhtin, A.M.Snigirev, Eur.Phys.J. C16 (2000) 527-536;
-* I.P.Lokhtin, A.M.Snigirev, e-print hep-ph/0406038.
+* Reference: I.P. Lokhtin, A.M. Snigirev, Eur. Phys. J. C 46 (2006) 211
*
-*
*======================================================================
SUBROUTINE PYQUEN(A,ifb,bfix)
common /plpar1/ tau0,T0,TC,sigqq,AW,RA,mvisc,nf
common /plglur/ glur(1000,4),kglu(1000,6),nrg,nrgm
common /plquar/ pqua(1000,5),kqua(1000,5),nrq
- common /parimp/ b1,psib1,rb1,rb2
+ common /parimp/ b1,psib1,rb1,rb2,noquen
+ common /pyqpar/ T0u,tau0u,nfu,ienglu,ianglu
common /bintaa/ br
common /plfpar/ bgen
save /pyjets/, /pydat1/, /pysubs/, /plglur/, /plquar/
AW=A ! atomic weight
RA=1.15d0*AW**0.333333d0 ! nucleus radius in fm
RA2=2.d0*RA
- nf=0 ! number of active flavours in QGP
- TC=0.2d0 ! crutical temperature
- tau0=0.1d0 ! proper time of QGP formation
mvisc=0 ! flag of QGP viscosity (off here)
-*
+ TC=0.2d0 ! crutical temperature
+
+ if(nfu.ne.0.and.nfu.ne.1.and.nfu.ne.2.and.nfu.ne.3) nfu=0
+ nf=nfu ! number of active flavours in QGP
+ if(tau0u.lt.0.01d0.or.tau0u.gt.10.d0) tau0u=0.1d0
+ tau0=tau0u ! proper time of QGP formation
+ if(T0u.lt.0.2d0.or.T0u.gt.2.d0) T0u=1.d0
+ T000=T0u ! initial QGP temperatute
+ if(ienglu.ne.0.and.ienglu.ne.1.and.ienglu.ne.2) ienglu=0 ! e-loss type
+ if(ianglu.ne.0.and.ianglu.ne.1.and.ianglu.ne.2) ianglu=0 ! angular spec.
+*
pi=3.14159d0
* avoid stopping run if pythia does not conserve energy due to collisional loss
mstu(21)=1
-
+
* generate impact parameter of A-A collision with jet production
if(ifb.eq.0) then
if(bfix.lt.0.d0) then
> aih,aiabs)
rtaa=rtaa0*(1.d0-br*(1.d0+(1.d0-0.25d0*br)*dlog(1.d0/br)+
> 4.d0*rest/pi))
- T00=((rtaa*sb0)/(rtaa0*sb))**0.25d0
+ T00=T000*((rtaa*sb0)/(rtaa0*sb))**0.25d0
T0=T00*(AW/207.d0)**0.166667d0
* generate single event with partonic energy loss
+ nrg=0
ehard=ckin(3)
if(b1.le.1.85d0*RA) then
call plinit(ehard)
ip02=0
ip001=0
ip002=0
- if(mstj(41).ne.0) goto 5
+ if(mstj(41).ne.0) goto 5
mstj(41)=1
nn=n
do i=9,nn
if(k(i,3).eq.7) then
ip1=i ! first hard parton (line ip1)
kfh1=k(i,1) ! status code of first hard parton
- qmax1=p(i,4) ! energy of first hard parton
+ qmax1=pyp(i,10) ! transverse momentum of first hard parton
end if
if(k(i,3).eq.8) then
ip2=i ! second hard parton (line ip2)
kfh2=k(i,1) ! status code of second hard parton
- qmax2=p(i,4) ! energy of second hard parton
+ qmax2=pyp(i,10) ! transverse momentum of second hard parton
end if
end do
n2=n
call pyshow(ip2,0,qmax2) ! vacuum showering for second hard parton
if(n.eq.n2) ip2=0
+ mstj(41)=0
if(n.eq.nn) goto 5
* find two leading partons after showering
v(ip1,j)=v(ip01,j)
p(ip1,j)=p(ip01,j)
end do
- k(ip1,1)=kfh1
+ k(ip1,1)=kfh1
+* fix first/last daughter for moving entry
+ do jgl=1,n
+ if(k(jgl,4).eq.ip01) k(jgl,4)=ip1
+ if(k(jgl,5).eq.ip01) k(jgl,5)=ip1
+ end do
+*
end if
if(ip2.gt.0) then
do j=1,5
p(ip2,j)=p(ip02,j)
end do
k(ip2,1)=kfh2
+* fix first/last daughter for moving entry
+ do jgl=1,n
+ if(k(jgl,4).eq.ip02) k(jgl,4)=ip2
+ if(k(jgl,5).eq.ip02) k(jgl,5)=ip2
+ end do
+*
end if
* add final showering gluons to the list of in-medium emitted gluons,
nis(ns+1)=nis(ns+1)+1
elseif(ks.eq.1.and.nis(ns+1).gt.0) then
nis(ns+1)=nis(ns+1)+1
- nes=nes+nis(ns+1) ! nes - total number if entries
+ nes=nes+nis(ns+1) ! nes - total number of entries
nss(ns+1)=nes
ns=ns+1
elseif(ks.ne.2.and.ksp.ne.2.and.ns.gt.0) then
k(n,j)=k(i,j)
p(n,j)=p(i,j)
end do
+* fix first/last daughter for moving entry
+ do jgl=1,n
+ if(k(jgl,4).eq.i) k(jgl,4)=n
+ if(k(jgl,5).eq.i) k(jgl,5)=n
+ end do
+*
do in=i+1,n
do j=1,5
v(in-1,j)=v(in,j)
k(in-1,j)=k(in,j)
p(in-1,j)=p(in,j)
end do
+* fix first/last daughter for moving entry
+ do jgl=1,n
+ if(k(jgl,4).eq.in) k(jgl,4)=in-1
+ if(k(jgl,5).eq.in) k(jgl,5)=in-1
+ end do
+*
end do
do ip=1,nrg
ku=kglu(ip,6)
k(is,j)=k(i,j)
p(is,j)=p(i,j)
end do
+* fix first/last daughter for moving entries
+ do jgl=1,n
+ if(k(jgl,4).eq.i) k(jgl,4)=is
+ if(k(jgl,5).eq.i) k(jgl,5)=is
+ end do
+*
end do
do ia=ns-1,1,-1
do i=nss(ia+1)-1,nss(ia),-1
k(is,j)=k(i,j)
p(is,j)=p(i,j)
end do
+* fix first/last daughter for moving entries
+ do jgl=1,n
+ if(k(jgl,4).eq.i) k(jgl,4)=is
+ if(k(jgl,5).eq.i) k(jgl,5)=is
+ end do
+*
end do
end do
p(ia,1)=ptg*dcos(phig)
p(ia,2)=ptg*dsin(phig)
p(ia,3)=dsqrt(abs(eg*eg-ptg*ptg))
- if(etag.lt.0.d0) p(ia,3)=-1.*p(ia,3)
- p(ia,4)=eg
+ if(etag.lt.0.d0) p(ia,3)=-1.d0*p(ia,3)
+ p(ia,4)=dsqrt(ptg*ptg+p(ia,3)**2)
p(ia,5)=0.d0
end do
external plvisc
common /plpar1/ tau0,T0,TC,sigqq,AW,RA,mvisc,nf
common /plpar2/ pln0,taupl,tauh,sigpl,sigh,sigplh,sigqqh,rg,rgn
- common /plevol/ taup(5000),temp(5000),denp(5000),enep(5000)
+ common /plevol/ taup(5000),temp(5000),denp(5000),enep(5000)
save /plevol/
*
pi=3.14159d0
* set intial plasma temperature, density and energy density in perfect
* (if mvisc=0) or viscous (mvisc=1,2) QGP with PLVISC subroitine
hst=0.15d0
- if(mvisc.eq.2.and.T0.gt.0.6d0) hst=0.25d0
+ if(T0.gt.1.5d0.or.mvisc.eq.2) hst=0.25d0
+ if(T0.gt.1.5d0.and.mvisc.ne.0) hst=0.9d0
T01=T0*5.06d0
TC1=TC*5.06d0
pln0=(16.d0+9.d0*nf)*1.2d0*(T01**3)/pi2
IMPLICIT DOUBLE PRECISION(A-H, O-Z)
IMPLICIT INTEGER(I-N)
INTEGER PYK,PYCHGE,PYCOMP
- external plthik, pln, plt, pls, pygauss, gluang
+ external plthik, pln, plt, pls, gauss, gluang
external pyp,pyr,pyk
common /plpar1/ tau0,T0,TC,sigqq,AW,RA,mvisc,nf
common /plpar2/ pln0,taupl,tauh,sigpl,sigh,sigplh,sigqqh,rg,rgn
common /thikpa/ fmax,xmin
common /pyjets/ n,npad,k(4000,5),p(4000,5),v(4000,5)
+ common /pyqpar/ T0u,tau0u,nfu,ienglu,ianglu
common /plglur/ glur(1000,4),kglu(1000,6),nrg,nrgm
common /factor/ cfac, kf
common /pleave/ taul, temlev
- common /parimp/ b1, psib1, rb1, rb2
+ common /parimp/ b1, psib1, rb1, rb2, noquen
common /plen/ epartc, um
common /plos/ elr,rsk
common /numje1/ nuj1, nuj2
- save /pyjets/, /plglur/
+ save /pyjets/, /plglur/
*
pi=3.14159d0
rb2=dsqrt(abs(r0*r0+b1*b1/4.d0-r0*b1*dcos(psib1)))
rb1=max(rb1,1.d-4)
rb2=max(rb2,1.d-4)
+ if(noquen.eq.1) goto 7
* find maximum of angular spectrum of radiated gluons with subroutine gluang
temin=0.5d0*pi
nrg=0
* generate changing 4-momentum of partons due to rescattering and energy loss
-* (for partons with |eta|<3 and p>5 GeV/c)
- nuj1=7 ! minimum number of rescattered parton
- nuj2=n ! maximum number of rescattered parton
+* (for partons with |eta|<3.5 and pt>3 GeV/c)
+ nuj1=9 ! minimum line number of rescattered parton
+ nuj2=n ! maximum line number of rescattered parton
do 2 ip=nuj1,nuj2 ! cycle on travelling partons
irasf=0
iraz=0
ks=k(ip,1) ! parton status code
kf=k(ip,2) ! parton identificator
ka=abs(kf)
- epart=abs(pyp(ip,10)) ! parton total momentum
- etar=pyp(ip,19) ! parton pseudorapidity
- if(epart.ge.5.d0.and.abs(etar).le.3.d0) then
+ ko=k(ip,3) ! origin (parent line number)
+ epart=abs(pyp(ip,10)) ! parton transverse momentum
+ etar=pyp(ip,19) ! parton pseudorapidity
+ if(ko.gt.6.and.epart.ge.3.d0.and.abs(etar).
+ > le.3.5d0) then
if(ka.eq.21.or.ka.eq.1.or.ka.eq.2.or.ka.eq.3.
> or.ka.eq.4.or.ka.eq.5.or.ka.eq.6.or.ka.eq.7.
> or.ka.eq.8) then
phir=pyp(ip,15) ! parton azimuthal angle
tetr=pyp(ip,13) ! parton polar angle
yrr=pyp(ip,17) ! parton rapidity
- stetr=max(dsin(tetr),1.d-04) ! parton sin(theta)
+ stetr=dsin(tetr) ! parton sin(theta)
+ if(abs(stetr).le.1.d-05) then
+ if(stetr.ge.0.d0) then
+ stetr=1.d-05
+ else
+ stetr=-1.d-05
+ end if
+ end if
phir1=-1.d0*phir
tetr1=-1.d0*tetr
if(psib1.lt.0.d0) bet=-1.d0*bet
phip=phir-bet
if(phip.gt.pi) phip=phip-2.d0*pi
- if(phip.lt.-1.d0*pi) phip=phip+2.d0*pi
- call pyrobo(0,0,0.d0,phir1,0.d0,0.d0,0.d0)
- call pyrobo(0,0,tetr1,0.d0,0.d0,0.d0,0.d0)
-
+ if(phip.lt.-1.d0*pi) phip=phip+2.d0*pi
+ call pyrobo(ip,ip,0.d0,phir1,0.d0,0.d0,0.d0)
+ call pyrobo(ip,ip,tetr1,0.d0,0.d0,0.d0,0.d0)
+
* calculate proper time of parton leaving QGP
aphin=(r0*r0-b1*b1/4.d0)/(rb1*rb2)
if(aphin.le.-1.d0) aphin=-0.99999d0
* set 4-momentum (in lab system) of next radiated gluon for parton number >8
* and fill arrays of radiated gluons in common block plglur
if(nrg.le.1000) then
- if(abs(elr).gt.0.1d0.and.ip.gt.8) then
- nrg=nrg+1
- 6 te1=temin*pyr(0)
- fte1=ftemax*pyr(0)
- fte2=gluang(te1)
- if(fte1.gt.fte2) goto 6
- tgl=te1 ! gaussian angular spectrum
-c tgl=0.d0 ! collinear angular spectrum
-c tgl=((0.5d0*pi*epartc)**pyr(0))/epartc ! broad-angular spectrum
+ if(abs(elr).gt.0.1d0.and.ip.gt.8) then
+* generate the angle of emitted gluon
+ if(ianglu.eq.0) then
+ 6 te1=temin*pyr(0)
+ fte1=ftemax*pyr(0)
+ fte2=gluang(te1)
+ if(fte1.gt.fte2) goto 6
+ tgl=te1
+ elseif (ianglu.eq.1) then
+ tgl=((0.5d0*pi*epartc)**pyr(0))/epartc
+ else
+ tgl=0.d0
+ end if
pgl=pi*(2.d0*pyr(0)-1.d0)
- pxgl=abs(elr)*(dcos(phir)*dcos(tgl)/dcosh(yrr)+
- > dcos(phir)*dsin(tgl)*dcos(pgl)*dsinh(yrr)/dcosh(yrr)-
+* in comoving system
+ pxgl=abs(elr)*stetr*(dcos(phir)*dcos(tgl)-
> dsin(phir)*dsin(tgl)*dsin(pgl))
- pygl=abs(elr)*(dsin(phir)*dcos(tgl)/dcosh(yrr)+
- > dsin(phir)*dsin(tgl)*dcos(pgl)*dsinh(yrr)/dcosh(yrr)+
- > dcos(phir)*dsin(tgl)*dsin(pgl))
- pzgl=abs(elr)*(dsinh(yrr)*dcos(tgl)-dsin(tgl)*dcos(pgl))
- > /dcosh(yrr)
- ptgl=dsqrt(abs(pxgl*pxgl+pygl*pygl))
+ pygl=abs(elr)*stetr*(dsin(phir)*dcos(tgl)+
+ > dcos(phir)*dsin(tgl)*dsin(pgl))
+ pzgl=-1.d0*abs(elr)*stetr*dsin(tgl)*dcos(pgl)
+ ptgl=dsqrt(abs(pxgl*pxgl+pygl*pygl))
psgl=dsqrt(abs(ptgl*ptgl+pzgl*pzgl))
+* recalculate in lab system
+ dyg=0.5d0*dlog(max(1.d-9,(psgl+pzgl)/(psgl-pzgl)))
+ pzgl=ptgl*dsinh(yrr+dyg)
+ psgl=dsqrt(abs(ptgl*ptgl+pzgl*pzgl))
+*
dpgl=pygl/pxgl
- glur(nrg,1)=abs(elr) ! energy
- glur(nrg,3)=datan(dpgl) ! phi
+ glur1=abs(elr) ! energy
+ glur3=datan(dpgl) ! phi
if(pxgl.lt.0.d0) then
if(pygl.ge.0.d0) then
- glur(nrg,3)=glur(nrg,3)+pi
+ glur3=glur3+pi
else
- glur(nrg,3)=glur(nrg,3)-pi
+ glur3=glur3-pi
end if
end if
- glur(nrg,4)=0.5d0*dlog(max(1.d-9,(psgl+pzgl)/(psgl-pzgl))) ! eta
- glur(nrg,2)=glur(nrg,1)/dcosh(glur(nrg,4)) ! pt
+ glur4=0.5d0*dlog(max(1.d-9,(psgl+pzgl)/(psgl-pzgl))) ! eta
+ glur2=glur1/dcosh(glur4) ! pt
+
+* put in event list radiated gluons with pt > 0.2 GeV only
+ if(glur2.ge.0.2d0) then
+ nrg=nrg+1
+* set gluon 4-momentum
+ glur(nrg,1)=glur1 ! energy
+ glur(nrg,2)=glur2 ! pt
+ glur(nrg,3)=glur3 ! phi
+ glur(nrg,4)=glur4 ! eta
* set gluon codes
- kglu(nrg,1)=2 ! status code
- kglu(nrg,2)=21 ! particle identificator
- kglu(nrg,3)=k(ipar,3) ! parent line number
- kglu(nrg,4)=0 ! special colour info
- kglu(nrg,5)=0 ! special colour info
- kglu(nrg,6)=ipar ! associated parton number
- end if
+ kglu(nrg,1)=2 ! status code
+ kglu(nrg,2)=21 ! particle identificator
+ kglu(nrg,3)=k(ipar,3) ! parent line number
+ kglu(nrg,4)=0 ! special colour info
+ kglu(nrg,5)=0 ! special colour info
+ kglu(nrg,6)=ipar ! associated parton number
+ end if
+ end if
else
write(6,*) 'Warning! Number of emitted gluons is too large!'
end if
p(ip,4)=dsqrt(p(ip,1)**2+p(ip,2)**2+p(ip,3)**2+p(ip,5)**2)
* boost to laboratory system
- 100 call pyrobo(0,0,tetr,phir,0.d0,0.d0,0.d0)
+ 100 call pyrobo(ip,ip,tetr,phir,0.d0,0.d0,0.d0)
end if
end if
end if
2 continue
-
+ 7 continue
+
return
end
******************************* END PLEVNT *************************
common /plpar1/ tau0,T0,TC,sigqq,AW,RA,mvisc,nf
common /plpar2/ pln0,taupl,tauh,sigpl,sigh,sigplh,sigqqh,rg,rgn
common /pyjets/ n,npad,k(4000,5),p(4000,5),v(4000,5)
+ common /pyqpar/ T0u,tau0u,nfu,ienglu,ianglu
common /pljdat/ ej, z, ygl, alfs, um, epa
common /pleave/ taul, temlev
common /radcal/ aa, bb
* boost to system of comoving plasma constituent
phir=pyp(i,15) ! parton phi
- tetr=pyp(i,13) ! parton theta
- stetr=max(dsin(tetr),1.d-08) ! parton sin(theta)
+ tetr=pyp(i,13) ! parton theta
phir1=-1.d0*phir
tetr1=-1.d0*tetr
- call pyrobo(0,0,0.d0,phir1,0.d0,0.d0,0.d0)
- call pyrobo(0,0,tetr1,0.d0,0.d0,0.d0,0.d0)
+ call pyrobo(i,i,0.d0,phir1,0.d0,0.d0,0.d0)
+ call pyrobo(i,i,tetr1,0.d0,0.d0,0.d0,0.d0)
pp=pyp(i,8) ! parton total momentum
ppl=abs(p(i,3)) ! parton pz
um=p(i,5) ! parton mass
rsk=dsqrt(abs(tp))
if(rsk.gt.ppl) rsk=ppl
-* calculate radiative energy loss per given scattering with subroutin plfun1
+* calculate radiative energy loss per given scattering with subroutine plfun1
ygl=y*cfac ! mean gluon free path in GeV^{-1}
elp=ygl*z ! mimimum radiated energy in LPM regime
ej=ppl
elr=resu*resun ! energy of radiated gluon
* to chancel radiative energy loss (optional case)
-c elr=0.d0
+ if(ienglu.eq.2) elr=0.d0
* to chancel collisional energy loss (optional case)
-c rsk=0.d0
+ if(ienglu.eq.1) rsk=0.d0
* determine the direction of parton moving
if(p(i,3).ge.0.d0) then
* calculate new 4-momentum of hard parton
phirs=2.d0*pi*pyr(0)
- epan=max(epa-rsk*rsk/(2.d0*ep0)-abs(elr),1.d0)
+ epan=epa-rsk*rsk/(2.d0*ep0)-abs(elr)
+ if(epan.lt.0.1d0) then
+ epan=epan+abs(elr)
+ elr=0.d0
+ if(epan.lt.0.1d0) then
+ rsk=0.d0
+ epan=epa
+ end if
+ end if
pptn=dsqrt(abs(rsk*rsk+(rsk**4)*(1.d0-epa*epa/(ppl*ppl))/
> (4.d0*ep0*ep0)-(rsk**4)*epa/(2.d0*ep0*ppl*ppl)-(rsk**4)/
> (4.d0*ppl*ppl)))
ppln=dsqrt(abs(epan*epan-pptn*pptn-p(i,5)**2))
- p(i,1)=pptn*dcos(phirs) ! px
- p(i,2)=pptn*dsin(phirs) ! py
- p(i,3)=sigp*ppln ! pz
- p(i,4)=epan ! E
-
+ p(i,1)=pptn*dcos(phirs) ! px
+ p(i,2)=pptn*dsin(phirs) ! py
+ p(i,3)=sigp*ppln ! pz
+ p(i,4)=dsqrt(p(i,1)**2+p(i,2)**2+p(i,3)**2+p(i,5)**2) ! E
* boost to system of hard parton
- 222 call pyrobo(0,0,tetr,phir,0.d0,0.d0,0.d0)
+ 222 call pyrobo(i,i,tetr,phir,0.d0,0.d0,0.d0)
return
end
IMPLICIT INTEGER(I-N)
INTEGER PYK,PYCHGE,PYCOMP
external plthik
- common /parimp/ b1, psib1, rb1, rb2
+ common /parimp/ b1, psib1, rb1, rb2, noquen
common /plpar1/ tau0,T0,TC,sigqq,AW,RA,mvisc,nf
rm1=dsqrt(abs(RA*RA-b1*b1/4.d0*(dsin(psib1)**2)))+
IMPLICIT DOUBLE PRECISION(A-H, O-Z)
IMPLICIT INTEGER(I-N)
INTEGER PYK,PYCHGE,PYCOMP
- common /parimp/ b1, psib1, rb1, rb2
+ common /parimp/ b1, psib1, rb1, rb2, noquen
common /plpar1/ tau0,T0,TC,sigqq,AW,RA,mvisc,nf
bu=X
r12=bu*bu+b1*b1/4.d0+bu*b1*dcos(psib1)
return
end
+* function to generate gauss distribution
+ double precision function gauss(x0,sig)
+ IMPLICIT DOUBLE PRECISION(A-H, O-Z)
+ IMPLICIT INTEGER(I-N)
+ 41 u1=pyr(0)
+ u2=pyr(0)
+ v1=2.d0*u1-1.d0
+ v2=2.d0*u2-1.d0
+ s=v1**2+v2**2
+ if(s.gt.1) go to 41
+ gauss=v1*dsqrt(-2.d0*dlog(s)/s)*sig+x0
+ return
+ end
+
* function to calculate angular distribution of emitted gluons
double precision function gluang(x)
IMPLICIT DOUBLE PRECISION(A-H, O-Z)