*----------------------------------------------------------------------
*
-* Filename : PYQUEN.F
+* Filename : PYQUEN.F
*
-* Author : Igor Lokhtin (Igor.Lokhtin@cern.ch)
-* Version : PYQUEN1_5.f
-* Last revision : 19-DEC-2007
+* Author : Igor Lokhtin (Igor.Lokhtin@cern.ch)
+* Version : PYQUEN1_5.f, v.1.5.1
+* Last revision v.5.1 : 19-DEC-2007
+* Last revision v.5.1.1 : 06-MAY-2010
*
*======================================================================
*
RA3=3.d0*RA
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
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
- if(ks.eq.2.or.ks.eq.1.or.ks.eq.21) then
+ if(ks.eq.2.or.ks.eq.1) then
phir=pyp(ip,15) ! parton azimuthal angle
tetr=pyp(ip,13) ! parton polar angle
yrr=pyp(ip,17) ! parton rapidity
if(irasf.gt.100000) goto 100
if(ran.lt.prob) goto 3
pltp=plt(tau,rj1,rj2,yrr)
+ if(pltp.le.TC) goto 100 ! escape if temperature drops below TC
pltp3=3.d0*pltp
pass=50.6d0/(pln(tau,rj1,rj2,yrr)*sigtr)
elr=0.d0
qm2=(scm-((um+ep0)**2))*(scm-((um-ep0)**2))/scm
bub=4.d0*tauu/TC
alf=6.d0*pi/((33.d0-2.d0*nf)*dlog(max(bub,1.d-10)))
- z=pi*4.d0*tauu*tauu*alf*(1.+nf/6.d0)
- if (z < 0.) write(6,*) "Warning in PLJETR: z < 0.", z
- bubs=dsqrt(abs(z))/TC
+ z=abs(pi*4.d0*tauu*tauu*alf*(1.+nf/6.d0))
+ bubs=dsqrt(z)/TC
alfs=6.d0*pi/((33.d0-2.d0*nf)*dlog(max(bubs,1.d-10)))
phmin2=z
phmax2=max(phmin2,qm2)
elp=ygl*z ! mimimum radiated energy in LPM regime
ej=ppl
bb=ej ! maximum radiated energy
- bbi=max(dsqrt(abs(z)),1.000001d0*elp)
+ bbi=max(dsqrt(z),1.000001d0*elp)
aa=min(bb,bbi) ! minimum radiated energy
hh=0.00001d0*(bb-aa)
REPS=0.01d0
* if quark production outside the QGP then
* arg=(((dsin(uu)*cosh(uu))**2)+((dcos(uu)*sinh(uu))**2))/(2.d0*uu*uu);
* here quark production inside the QGP
-c AM: better numerical stability
-
- if (uu .lt. 5.) then
- arg = ((dcos(uu)*cosh(uu))**2)+((dsin(uu)*sinh(uu))**2)
- xlogarg = dlog(max(arg,1.d-20))
- else
- xlogarg = log(0.25) + 2. * uu
- endif
-
-
- gl1=(ygl/(cfac*abs(z)))**0.3333333d0
+ arg=((dcos(uu)*cosh(uu))**2)+((dsin(uu)*sinh(uu))**2)
+ gl1=(ygl/(cfac*z))**0.3333333d0
gl2=(um/epa)**1.333333d0
- dc=1.d0/((1.d0+((gl1*gl2*or)**1.5d0))**2) ! massive parton
+ dc=1.d0/((1.d0+((gl1*gl2*or)**1.5d0))**2) ! massive parton
c dc=1.d0 !massless parton
-c plfun1 = dc*3.d0*alfs*ygl*dlog(max(arg,1.d-20))*spinf/(pi*al*or)
- plfun1 = dc*3.d0*alfs*ygl*xlogarg*spinf/(pi*al*or)
+ plfun1=dc*3.d0*alfs*ygl*dlog(max(arg,1.d-20))*spinf/(pi*al*or)
return
end