Updates needed for pyquen production.
authormorsch <morsch@f7af4fe6-9843-0410-8265-dc069ae4e863>
Fri, 16 Dec 2011 08:51:58 +0000 (08:51 +0000)
committermorsch <morsch@f7af4fe6-9843-0410-8265-dc069ae4e863>
Fri, 16 Dec 2011 08:51:58 +0000 (08:51 +0000)
Leticia Cunqueiro

PYTHIA6/CMakelibpythia6.4.21.pkg
PYTHIA6/QPYTHIA/pyquen1_5.F

index f37bf21..0b2cba5 100644 (file)
@@ -35,7 +35,9 @@ set ( EXPORT )
 
 set ( CSRCS  pythia6.4.21/main.c pythia6.4.21/pythia6_common_address.c)
 
-set ( FSRCS  pythia6.4.21/pythia6_common_block_address.F pythia6.4.21/tpythia6_called_from_cc.F pythia6.4.21/pythia-6.4.21.f pythia6.4.21/pydummy_6.4.21.f)
+set ( EINCLUDE PYTHIA6/QPYTHIA)
+
+set ( FSRCS  pythia6.4.21/pythia6_common_block_address.F pythia6.4.21/tpythia6_called_from_cc.F pythia6.4.21/pythia-6.4.21.f QPYTHIA/pyquen1_5.F)
 
 if( ALICE_TARGET STREQUAL "win32gcc")
        
index 06aea33..bed1692 100644 (file)
@@ -1,11 +1,12 @@
 
 *----------------------------------------------------------------------
 *
-*  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
 *
 *======================================================================
 *
@@ -48,7 +49,7 @@
       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 
@@ -1297,22 +1298,12 @@ C************************** END PARINV *************************************
 * 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