* Filename : PYQUEN.F
*
* First version created: 20-AUG-1997 Author : Igor Lokhtin
-* Last revision : 15-JUN-2004
+* Last revision : 23-SEP-2004
*
*======================================================================
*
IMPLICIT INTEGER(I-N)
INTEGER PYK,PYCHGE,PYCOMP
external pydata
- external pyp,pyr,pyk,pyjoin
+ external pyp,pyr,pyk,pyjoin,pyshow
external ftaa,funbip
common /pyjets/ n,npad,k(4000,5),p(4000,5),v(4000,5)
common /pydat1/ mstu(200),paru(200),mstj(200),parj(200)
common /pysubs/ msel,mselpd,msub(500),kfin(2,-40:40),ckin(200)
common /plpar1/ tau0,T0,TC,sigqq,AW,RA,mvisc,nf
- common /plglur/ glur(1000,4),kglu(1000,6),nrg
+ 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 /bintaa/ br
- common /plfpar/ bgen
- save /pyjets/, /pydat1/, /pysubs/, /plglur/
- dimension ijoin(1000),nis(100),nss(100),nas(100),nus(100)
+ common /plfpar/ bgen
+ save /pyjets/, /pydat1/, /pysubs/, /plglur/, /plquar/
+ dimension ijoik(2),ijoin(1000),nis(500),nss(500),nas(500),nus(500)
* set initial event paramters
AW=A ! atomic weight
if(b1.le.1.85d0*RA) then
call plinit(ehard)
call plevnt(ehard)
- end if
+ end if
+
+* reset all in-vacuum radiated guark 4-momenta and codes to zero
+ do i=1,1000
+ do j=1,5
+ pqua(i,j)=0.d0
+ kqua(i,j)=0
+ end do
+ end do
+ nrq=0
+
+* generate final state shower in vacuum if it was excluded before
+ nrgm=nrg ! fix number of in-medium emitted gluons
+ ip1=0
+ ip2=0
+ ip01=0
+ ip02=0
+ ip001=0
+ ip002=0
+ 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
+ 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
+ end if
+ end do
+
+ n1=n
+ call pyshow(ip1,0,qmax1) ! vacuum showering for first hard parton
+ if(n.eq.n1) ip1=0
+ n2=n
+ call pyshow(ip2,0,qmax2) ! vacuum showering for second hard parton
+ if(n.eq.n2) ip2=0
+ if(n.eq.nn) goto 5
+
+* find two leading partons after showering
+ do i=nn+1,n
+ if(k(i,3).eq.ip1) ip001=i ! first daughter of first hard parton
+ if(k(i,3).eq.ip2) ip002=i ! first daughter of second hard parton
+ end do
+ ptle1=0.d0
+ ptle2=0.d0
+ do i=nn+1,n
+ if (k(i,1).eq.14) goto 3
+ if(i.ge.ip002.and.ip002.gt.0) then
+ ptl02=pyp(i,10)
+ if(ptl02.gt.ptle2.and.k(i,2).eq.k(ip2,2)) then
+ ip02=i ! leading parton in second shower (line ip02)
+ ptle2=ptl02 ! pt of the leading parton
+ end if
+ elseif(ip001.gt.0) then
+ ptl01=pyp(i,10)
+ if(ptl01.gt.ptle1.and.k(i,2).eq.k(ip1,2)) then
+ ip01=i ! leading parton in first shower (line ip01)
+ ptle1=ptl01 ! pt of the leading parton
+ end if
+ end if
+ 3 continue
+ end do
-* stop generate event if there are no in-medium emitted gluons
+* replace two hard partons by two leading partons in original event record
+ if(ip1.gt.0) then
+ do j=1,5
+ v(ip1,j)=v(ip01,j)
+ p(ip1,j)=p(ip01,j)
+ end do
+ k(ip1,1)=kfh1
+ end if
+ if(ip2.gt.0) then
+ do j=1,5
+ v(ip2,j)=v(ip02,j)
+ p(ip2,j)=p(ip02,j)
+ end do
+ k(ip2,1)=kfh2
+ end if
+
+* add final showering gluons to the list of in-medium emitted gluons,
+* fill the list of emitted quarks by final showering quark pairs,
+* and remove showering gluons and quarks from the event record
+ do i=nn+1,n
+ if(k(i,1).eq.14.or.i.eq.ip01.or.i.eq.ip02) goto 12
+ if(k(i,2).ne.21) then ! filling 'plquar' arrays for quarks
+ nrq=nrq+1
+ do j=1,5
+ kqua(nrq,j)=k(i,j)
+ pqua(nrq,j)=p(i,j)
+ end do
+ kqua(nrq,1)=2
+ goto 12
+ end if
+ if(i.ge.ip002.and.ip002.gt.0) then
+ ish=ip2
+ else
+ ish=ip1
+ end if
+ nrg=nrg+1
+ nur=nrg
+ 7 ishm=kglu(nur-1,6)
+ if(ish.ge.ishm.or.nur.le.2) goto 6 ! adding gluons in 'plglur' arrays
+ do j=1,6
+ kglu(nur,j)=kglu(nur-1,j)
+ end do
+ do j=1,4
+ glur(nur,j)=glur(nur-1,j)
+ end do
+ nur=nur-1
+ goto 7
+ 6 kglu(nur,1)=2 ! status code
+ kglu(nur,2)=k(i,2) ! particle identificator
+ kglu(nur,3)=k(ish,3) ! parent line number
+ kglu(nur,4)=0 ! special colour info
+ kglu(nur,5)=0 ! special colour info
+ kglu(nur,6)=ish ! associated parton number
+ glur(nur,1)=p(i,4) ! energy
+ glur(nur,2)=pyp(i,10) ! pt
+ glur(nur,3)=pyp(i,15) ! phi
+ glur(nur,4)=pyp(i,19) ! eta
+ 12 continue
+ do j=1,5 ! remove partons from event list
+ v(i,j)=0.d0
+ k(i,j)=0
+ p(i,j)=0.d0
+ end do
+ end do
+ n=nn
+
+ 5 continue
+
+* stop generate event if there are no additional gluons
if(nrg.lt.1) goto 1
* define number of stirngs (ns) and number of entries in strings before
nes=0
i0=0
i1=0
- do i=1,100
- nis(i)=0
- nas(i)=0
- nss(i)=0
- nus(i)=0
+ do i=1,500
+ nis(i)=0
+ nas(i)=0
+ nss(i)=0
+ nus(i)=0
end do
do i=9,n
- ks=k(i,1)
- ksp=k(i-1,1)
- if(ks.eq.2) then
- 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
- nss(ns+1)=nes
- ns=ns+1
+ ks=k(i,1)
+ ksp=k(i-1,1)
+ if(ks.eq.2) then
+ 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
+ nss(ns+1)=nes
+ ns=ns+1
elseif(ks.ne.2.and.ksp.ne.2.and.ns.gt.0) then
- i1=i1+1 ! last i1 lines not included in strings
+ i1=i1+1 ! last i1 lines not included in strings
end if
end do
- i0=n-nes-i1 ! first i0 lines not included in strings
+ i0=n-nes-i1 ! first i0 lines not included in strings
do i=1,ns
- nss(i)=nss(i)+i0
+ nss(i)=nss(i)+i0
end do
* move fragmented particles in bottom of event list
- icount = 0
i=i0+1
2 ks=k(i,1)
- icount = icount + 1
- if (icount.gt.200) stop
ksp=k(i-1,1)
if(ks.ne.2.and.ksp.ne.2) then
- n=n+1
- do j=1,5
- v(n,j)=v(i,j)
- k(n,j)=k(i,j)
- p(n,j)=p(i,j)
- 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
- end do
- do ip=1,nrg
- ku=kglu(ip,6)
- if(ku.gt.i) kglu(ip,6)=ku-1
- end do
- n=n-1
+ n=n+1
+ do j=1,5
+ v(n,j)=v(i,j)
+ k(n,j)=k(i,j)
+ p(n,j)=p(i,j)
+ 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
+ end do
+ do ip=1,nrg
+ ku=kglu(ip,6)
+ if(ku.gt.i) kglu(ip,6)=ku-1
+ end do
+ n=n-1
else
- i=i+1
+ i=i+1
end if
if(i.le.n-i1) goto 2
+
* define number of additional entries in strings, nas(ns)
do i=1,nrg
- kas=kglu(i,6)
- if(kas.le.nss(1)) then
- nas(1)=nas(1)+1
- else
- do j=2,ns
- if(kas.le.nss(j).and.kas.gt.nss(j-1))
- > nas(j)=nas(j)+1
- end do
- end if
+ kas=kglu(i,6)
+ if(kas.le.nss(1)) then
+ nas(1)=nas(1)+1
+ else
+ do j=2,ns
+ if(kas.le.nss(j).and.kas.gt.nss(j-1))
+ > nas(j)=nas(j)+1
+ end do
+ end if
end do
do j=1,ns
- do i=1,j
- nus(j)=nus(j)+nas(i)
- end do
+ do i=1,j
+ nus(j)=nus(j)+nas(i)
+ end do
end do
-
-* add emitted gluons in event list
+
+* add emitted gluons in event list
nu=n
n=n+nrg
do i=nu,nu-i1,-1
- is=i+nrg
- do j=1,5
- v(is,j)=v(i,j)
- k(is,j)=k(i,j)
- p(is,j)=p(i,j)
- end do
+ is=i+nrg
+ do j=1,5
+ v(is,j)=v(i,j)
+ k(is,j)=k(i,j)
+ p(is,j)=p(i,j)
+ end do
end do
do ia=ns-1,1,-1
- do i=nss(ia+1)-1,nss(ia),-1
- is=i+nus(ia)
- do j=1,5
- v(is,j)=v(i,j)
- k(is,j)=k(i,j)
- p(is,j)=p(i,j)
- end do
- end do
+ do i=nss(ia+1)-1,nss(ia),-1
+ is=i+nus(ia)
+ do j=1,5
+ v(is,j)=v(i,j)
+ k(is,j)=k(i,j)
+ p(is,j)=p(i,j)
+ end do
+ end do
end do
-
+
do i=1,nrg
- if(i.le.nus(1)) then
- ia=nss(1)-1+i
- else
- do in=2,ns
- if(i.le.nus(in).and.i.gt.nus(in-1))
- > ia=nss(in)-1+i
- end do
- end if
- eg=glur(i,1)
- ptg=glur(i,2)
- phig=glur(i,3)
- etag=glur(i,4)
- do j=1,5
- v(ia,j)=0.d0
- k(ia,j)=kglu(i,j)
- 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
- p(ia,5)=0.d0
+ if(i.le.nus(1)) then
+ ia=nss(1)-1+i
+ else
+ do in=2,ns
+ if(i.le.nus(in).and.i.gt.nus(in-1))
+ > ia=nss(in)-1+i
+ end do
+ end if
+ eg=glur(i,1)
+ ptg=glur(i,2)
+ phig=glur(i,3)
+ etag=glur(i,4)
+ do j=1,5
+ v(ia,j)=0.d0
+ k(ia,j)=kglu(i,j)
+ 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
+ p(ia,5)=0.d0
end do
-* rearrange partons to form strings
+* rearrange partons to form strings in event list
do ij=1,1000
- ijoin(ij)=0
+ ijoin(ij)=0
end do
do i=1,ns
- njoin=nis(i)+nas(i)
- if(i.eq.1) then
- do j=1,njoin
- ijoin(j)=i0+j
- end do
- else
- do j=1,njoin
- ijoin(j)=nss(i-1)+nus(i-1)+j
- end do
- end if
- call pyjoin(njoin,ijoin)
+ njoin=nis(i)+nas(i)
+ if(i.eq.1) then
+ do j=1,njoin
+ ijoin(j)=i0+j
+ end do
+ else
+ do j=1,njoin
+ ijoin(j)=nss(i-1)+nus(i-1)+j
+ end do
+ end if
+ call pyjoin(njoin,ijoin)
end do
-
+
+* add in-vacuum emitted quark pairs
+ if(nrq.lt.2) goto 1
+ do i=1,nrq,2
+ n=n+2
+ do j=1,5
+ v(n-1,j)=0.d0
+ k(n-1,j)=kqua(i,j)
+ p(n-1,j)=pqua(i,j)
+ end do
+ in=i+1
+ 4 ktest=k(n-1,2)+kqua(in,2)
+ if(ktest.eq.0.or.in.eq.nrq) goto 8
+ in=in+1
+ goto 4
+ 8 do j=1,5
+ v(n,j)=0.d0
+ k(n,j)=kqua(in,j)
+ p(n,j)=pqua(in,j)
+ end do
+ if(in.gt.i+1) then
+ do j=1,5
+ kqua(in,j)=kqua(i+1,j)
+ pqua(in,j)=pqua(i+1,j)
+ end do
+ end if
+ end do
+
+ do ij=1,2
+ ijoik(ij)=0
+ end do
+ do i=1,nrq-1,2
+ k(n+1-i,1)=1
+ ijoik(1)=n-i
+ ijoik(2)=n+1-i
+ call pyjoin(2,ijoik)
+ end do
+
1 continue
return
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 /plglur/ glur(1000,4),kglu(1000,6),nrg
+ common /plglur/ glur(1000,4),kglu(1000,6),nrg,nrgm
common /factor/ cfac, kf
common /pleave/ taul, temlev
common /parimp/ b1, psib1, rb1, rb2
* reset all radiated gluon 4-momenta and codes to zero -------------------
do i=1,1000
- do j=1,4
+ do j=1,4
glur(i,j)=0.d0
kglu(i,j)=0
end do
psgl=dsqrt(abs(ptgl*ptgl+pzgl*pzgl))
dpgl=pygl/pxgl
glur(nrg,1)=abs(elr) ! energy
- glur(nrg,3)=datan(dpgl) ! phi
+ glur(nrg,3)=datan(dpgl) ! phi
if(pxgl.lt.0.d0) then
if(pygl.ge.0.d0) then
glur(nrg,3)=glur(nrg,3)+pi
if(ppl.lt.3.d0) goto 222 ! no energy loss if pz<3 GeV/c
* generation hard parton-plasma scattering with momentum transfer rsk
- 221 ep0=-1.*tauu*(dlog(max(1.d-10,pyr(0)))+dlog(max(1.d-10,
+ 221 ep0=-1.*tauu*(dlog(max(1.d-10,pyr(0)))+dlog(max(1.d-10,
> pyr(0)))+dlog(max(1.d-10,pyr(0)))) ! energy of 'thermal' parton
iter=iter+1
if(ep0.lt.1.d-10.and.iter.le.100000) goto 221
return
end
-* function to generate gauss distribution
-c double precision function gauss(x0,sig)
-c IMPLICIT DOUBLE PRECISION(A-H, O-Z)
-c IMPLICIT INTEGER(I-N)
-c 41 u1=pyr(0)
-c u2=pyr(0)
-c v1=2.d0*u1-1.d0
-c v2=2.d0*u2-1.d0
-c s=v1**2+v2**2
-c if(s.gt.1) go to 41
-c gauss=v1*dsqrt(-2.d0*dlog(s)/s)*sig+x0
-c return
-c end
-
* function to calculate angular distribution of emitted gluons
double precision function gluang(x)
IMPLICIT DOUBLE PRECISION(A-H, O-Z)
git://git.uio.no / u / mrichter / AliRoot.git / blobdiff