Possible conflicts between Heavy Flavor and Atlas Tuning removed.
[u/mrichter/AliRoot.git] / PYTHIA6 / pyquen.F
index 06c2d4d..9eaebe1 100644 (file)
@@ -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 
 *
 *======================================================================
 *
       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 > 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 
       IMPLICIT DOUBLE PRECISION(A-H, O-Z)
       IMPLICIT INTEGER(I-N)
       INTEGER PYK,PYCHGE,PYCOMP
-      external plthik, pln, plt, pls, gauss, gluang 
+      external plthik, pln, plt, pls, pygauss, 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 
 
 * 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   
@@ -1018,20 +1188,6 @@ c      dc=1.d0                                         !massless parton
        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)