From a0e56495f7608a554d584e9d761d255983726e36 Mon Sep 17 00:00:00 2001 From: morsch Date: Fri, 1 Oct 2004 11:14:04 +0000 Subject: [PATCH] This version foresees the option to treat vacuum parton showering (final state radiation in Pythia) after, as well as before, medium-induced partonic energy loss. Default option is "vacuum showering before in-medium loss"; the option "vacuum showering after in-medium loss" can be selected just by switching off Pythia final state radiation in main routine, mstj(41)=0 (see updated Pyquen description). (I. Lokhtin) --- PYTHIA6/pyquen.F | 428 +++++++++++++++++++++++++++++++++-------------- 1 file changed, 299 insertions(+), 129 deletions(-) diff --git a/PYTHIA6/pyquen.F b/PYTHIA6/pyquen.F index 8f97a49a6e9..74922ea943d 100644 --- a/PYTHIA6/pyquen.F +++ b/PYTHIA6/pyquen.F @@ -3,7 +3,7 @@ * Filename : PYQUEN.F * * First version created: 20-AUG-1997 Author : Igor Lokhtin -* Last revision : 15-JUN-2004 +* Last revision : 23-SEP-2004 * *====================================================================== * @@ -23,18 +23,19 @@ 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 @@ -91,9 +92,143 @@ 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 @@ -102,144 +237,179 @@ 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 @@ -314,13 +484,13 @@ 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 /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 @@ -357,7 +527,7 @@ * 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 @@ -476,7 +646,7 @@ c tgl=((0.5d0*pi*epartc)**pyr(0))/epartc ! broad-angular spectrum 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 @@ -576,7 +746,7 @@ c tgl=((0.5d0*pi*epartc)**pyr(0))/epartc ! broad-angular spectrum 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 @@ -1019,18 +1189,18 @@ c dc=1.d0 !massless parton 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 + 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) -- 2.43.5