- Update to PYQUEN v1.3 (http://lokhtin.web.cern.ch/lokhtin/pyquen/pyquen.txt)

- Possibility to set the pyquen parameters.

diff --git a/PYTHIA6/AliPythia.cxx b/PYTHIA6/AliPythia.cxx
index 7e6f0702ada76beef433f7ab69f4d8bb4572ea7b..9173c056a66f04a4b56d8566ac17ec5194b70794 100644 (file)
--- a/PYTHIA6/AliPythia.cxx
+++ b/PYTHIA6/AliPythia.cxx
@@ -18,9 +18,10 @@
 
 #include "AliPythia.h"
 #include "AliPythiaRndm.h"
-#include "../FASTSIM/AliFastGlauber.h"
-#include "../FASTSIM/AliQuenchingWeights.h"
+#include "AliFastGlauber.h"
+#include "AliQuenchingWeights.h"
 #include "TVector3.h"
+#include "PyquenCommon.h"
 
 ClassImp(AliPythia)
 
@@ -1209,9 +1210,25 @@ void  AliPythia::Quench()
 void AliPythia::Pyquen(Double_t a, Int_t ibf, Double_t b)
 {
     // Igor Lokthine's quenching routine
+    // http://lokhtin.web.cern.ch/lokhtin/pyquen/pyquen.txt
+
     pyquen(a, ibf, b);
 }
 
+void AliPythia::SetPyquenParameters(Double_t t0, Double_t tau0, Int_t nf, Int_t iengl, Int_t iangl)
+{
+    // Set the parameters for the PYQUEN package.
+    // See comments in PyquenCommon.h
+    
+    
+    PYQPAR.t0    = t0;
+    PYQPAR.tau0  = tau0;
+    PYQPAR.nf    = nf;
+    PYQPAR.iengl = iengl;
+    PYQPAR.iangl = iangl;
+}
+
+
 void AliPythia::Pyevnw()
 {
     // New multiple interaction scenario

diff --git a/PYTHIA6/AliPythia.h b/PYTHIA6/AliPythia.h
index 3b4645de1e1c57d62c1f50ec26afd1449a51f202..ddb05087fa155983d3e7e48d570f9d4e04f2aaa5 100644 (file)
--- a/PYTHIA6/AliPythia.h
+++ b/PYTHIA6/AliPythia.h
@@ -45,6 +45,7 @@ class AliPythia : public TPythia6, public AliRndm
     virtual void Pyshow(Int_t ip1, Int_t ip2, Double_t qmax);
     virtual void Pyrobo(Int_t imi, Int_t ima, Double_t the, Double_t phi, Double_t bex, Double_t bey, Double_t bez);
     virtual void InitQuenching(Float_t bmin, Float_t bmax, Float_t k, Int_t iECMethod, Float_t zmax = 0.97, Int_t ngmax = 30);
+    virtual void SetPyquenParameters(Double_t t0, Double_t tau0, Int_t nf, Int_t iengl, Int_t iangl);
     virtual void Pyquen(Double_t a, Int_t ibf, Double_t b);
     virtual void GetQuenchingParameters(Double_t& xp, Double_t& yp, Double_t z[4]);
     // return instance of the singleton

diff --git a/PYTHIA6/PyquenCommon.h b/PYTHIA6/PyquenCommon.h
new file mode 100644 (file)
index 0000000..1f4a0a2
--- /dev/null
+++ b/PYTHIA6/PyquenCommon.h
@@ -0,0 +1,58 @@
+#ifndef ROOT_PyqCommon
+#define ROOT_PyqCommon
+
+#ifndef __CFORTRAN_LOADED
+//*KEEP,cfortran.
+#include "cfortran.h"
+//*KEND.
+#endif
+
+extern "C" {
+
+//  common /pyqpar/ T0u,tau0u,nfu,ienglu,ianglu
+    
+typedef struct {
+    Double_t t0;
+    Double_t tau0;
+    Int_t    nf;
+    Int_t    iengl;
+    Int_t    iangl;
+} PyqparCommon;
+
+#define PYQPAR COMMON_BLOCK(PYQPAR,pyqpar)
+ COMMON_BLOCK_DEF(PyqparCommon,PYQPAR);
+}
+
+/*
+Parameters in COMMON BLOCK PYQPAR can be varied by user: 
+
+
+COMMON /pyqpar/ T0,tau0,nf,ienglu,ianglu
+
+T0 - initial temparature of quark-gluon plasma  
+(allowed range is 0.2 GeV < T0 < 2 GeV, default value is T0=1 GeV);
+
+tau0 - proper time of quark-gluon plasma formation
+(allowed range is 0.01 < tau0 < 10 fm/c, default value is tau0=0.1 fm/c)
+
+nf - number of active quark flavours in quark-gluon plasma
+(nf=0, 1, 2 or 3, default value is nf=0);
+
+ienglu - flag to fix type of medium-induced partonic energy loss
+(ienglu=0 - radiative and collisional loss,
+ienglu=1 - radiative loss only, ienglu=2 - collisional loss only,
+default value is ienglu=0);
+
+ianglu - flag to fix type of angular distribution of emitted gluons
+(ianglu=0 - small-angular, ianglu=1 - wide-angular, ianglu=2 - collinear,
+default value is ianglu-0).  
+
+NOTE! If specified by user value of such parameter is out of allowed range, 
+the default value is used in PYQUEN run. 
+
+NOTE! Default parameters of quark-gluon plasma (T0, tau0, nf) were selected as 
+an estimation for LHC heavy ion beam energies. 
+*/
+
+
+#endif

diff --git a/PYTHIA6/libpythia6.pkg b/PYTHIA6/libpythia6.pkg
index 5d07cf3b887b333163ac1e22356b6512b70de2ec..4e9e5861e506a95cfb651696cecf6843868d5423 100644 (file)
--- a/PYTHIA6/libpythia6.pkg
+++ b/PYTHIA6/libpythia6.pkg
@@ -6,17 +6,22 @@ DHDR:=pythia6LinkDef.h
 
 EXPORT:=
 
+PACKFFLAGS:=$(FFLAGS) -DPYTHIA
+
 CSRCS:= \
 main.c pythia6_common_address.c 
 
+EINCLUDE:= PYTHIA6/charybdis
+
 FSRCS:= \
 pythia6214.f \
 pdf_alice.F\
 pythia6_common_block_address.F \
 tpythia6_called_from_cc.F\
 openDecayTable.F\
-pyquen.F\
-lhaglue.F
+pyquen1_3.F\
+lhaglue.F\
+charybdis/charybdis1003.F
 
 ifeq ($(ALICE_TARGET),linux)
    GCC_VERSION    := $(shell $(CXX) -v 2>&1 | \

diff --git a/PYTHIA6/pyquen.F b/PYTHIA6/pyquen1_3.F
similarity index 85%
rename from PYTHIA6/pyquen.F
rename to PYTHIA6/pyquen1_3.F
index 9eaebe179baee05138a4f057997b6433033e61a0..24a5f62f7eac3f72e2195a15c414c98247a62540 100644 (file)
--- a/PYTHIA6/pyquen.F
+++ b/PYTHIA6/pyquen1_3.F
@@ -2,20 +2,19 @@
 *
 *  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)
@@ -31,7 +30,8 @@
       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/
@@ -41,16 +41,23 @@
       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    
@@ -84,10 +91,11 @@
      >           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)  
@@ -111,19 +119,19 @@
       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
       
@@ -133,6 +141,7 @@
       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  
@@ -166,7 +175,13 @@
         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  
@@ -174,6 +189,12 @@
         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,
@@ -250,7 +271,7 @@
         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 
@@ -273,12 +294,24 @@
         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) 
@@ -318,6 +351,12 @@
         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 
@@ -327,6 +366,12 @@
          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 
 
@@ -350,8 +395,8 @@
        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  
       
@@ -424,7 +469,7 @@
       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 
@@ -442,7 +487,8 @@
 * 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
@@ -484,20 +530,21 @@
       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              
 
@@ -519,6 +566,7 @@
       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 
@@ -537,18 +585,20 @@
       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    
@@ -556,7 +606,14 @@
         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 
 
@@ -576,10 +633,10 @@
         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
@@ -624,46 +681,63 @@
 * 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 
@@ -693,12 +767,13 @@ c          tgl=((0.5d0*pi*epartc)**pyr(0))/epartc ! broad-angular spectrum
         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 ************************* 
@@ -714,6 +789,7 @@ c          tgl=((0.5d0*pi*epartc)**pyr(0))/epartc ! broad-angular spectrum
       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 
@@ -730,12 +806,11 @@ c          tgl=((0.5d0*pi*epartc)**pyr(0))/epartc ! broad-angular spectrum
 
 * 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 
@@ -769,7 +844,7 @@ c          tgl=((0.5d0*pi*epartc)**pyr(0))/epartc ! broad-angular spectrum
        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                           
@@ -792,9 +867,9 @@ c          tgl=((0.5d0*pi*epartc)**pyr(0))/epartc ! broad-angular spectrum
        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 
@@ -805,18 +880,25 @@ c       rsk=0.d0
 
 * 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
@@ -830,7 +912,7 @@ c       rsk=0.d0
        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)))+
@@ -1179,7 +1261,7 @@ c      dc=1.d0                                         !massless parton
        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) 
@@ -1188,6 +1270,20 @@ 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)