PYTHIA6/pyquen1_3.F

   1 *----------------------------------------------------------------------
   2 *
   3 *  Filename             : PYQUEN.F
   4 *
   5 *  Author               : Igor Lokhtin  (Igor.Lokhtin@cern.ch)
   6 *  Version              : PYQUEN1_3.f
   7 *  Last revision        : 03-OCT-2007
   8 *
   9 *======================================================================
  10 *
  11 *  Description : Event generator for simulation of parton rescattering
  12 *                and energy loss in expanding quark-gluon plasma created
  13 *                in ultrarelativistic heavy ion AA collisons
  14 *               (implemented as modification of standard Pythia jet event)
  15 *
  16 *  Reference: I.P. Lokhtin, A.M. Snigirev, Eur. Phys. J. C 46 (2006) 211
  17 *
  18 *======================================================================
  19
  20       SUBROUTINE PYQUEN(A,ifb,bfix)
  21       IMPLICIT DOUBLE PRECISION(A-H, O-Z)
  22       IMPLICIT INTEGER(I-N)
  23       INTEGER PYK,PYCHGE,PYCOMP
  24       external pydata
  25       external pyp,pyr,pyk,pyjoin,pyshow
  26       external ftaa,funbip
  27       common /pyjets/ n,npad,k(4000,5),p(4000,5),v(4000,5)
  28       common /pydat1/ mstu(200),paru(200),mstj(200),parj(200)
  29       common /pysubs/ msel,mselpd,msub(500),kfin(2,-40:40),ckin(200)
  30       common /plpar1/ tau0,T0,TC,sigqq,AW,RA,mvisc,nf
  31       common /plglur/ glur(1000,4),kglu(1000,6),nrg,nrgm
  32       common /plquar/ pqua(1000,5),kqua(1000,5),nrq
  33       common /parimp/ b1,psib1,rb1,rb2,noquen
  34       common /pyqpar/ T0u,tau0u,nfu,ienglu,ianglu
  35       common /bintaa/ br
  36       common /plfpar/ bgen
  37       save /pyjets/, /pydat1/, /pysubs/, /plglur/, /plquar/
  38       dimension ijoik(2),ijoin(1000),nis(500),nss(500),nas(500),nus(500)
  39
  40 * set initial event paramters
  41       AW=A                                 ! atomic weight
  42       RA=1.15d0*AW**0.333333d0             ! nucleus radius in fm
  43       RA2=2.d0*RA
  44       mvisc=0                              ! flag of QGP viscosity (off here)
  45       TC=0.2d0                             ! crutical temperature
  46
  47       if(nfu.ne.0.and.nfu.ne.1.and.nfu.ne.2.and.nfu.ne.3) nfu=0
  48       nf=nfu                               ! number of active flavours in QGP
  49       if(tau0u.lt.0.01d0.or.tau0u.gt.10.d0) tau0u=0.1d0
  50       tau0=tau0u                          ! proper time of QGP formation
  51       if(T0u.lt.0.2d0.or.T0u.gt.2.d0) T0u=1.d0
  52       T000=T0u                            ! initial QGP temperatute
  53       if(ienglu.ne.0.and.ienglu.ne.1.and.ienglu.ne.2) ienglu=0 ! e-loss type
  54       if(ianglu.ne.0.and.ianglu.ne.1.and.ianglu.ne.2) ianglu=0 ! angular spec.
  55 *
  56       pi=3.14159d0
  57
  58 * avoid stopping run if pythia does not conserve energy due to collisional loss
  59       mstu(21)=1
  60
  61 * generate impact parameter of A-A collision with jet production
  62       if(ifb.eq.0) then
  63        if(bfix.lt.0.d0) then
  64         write(6,*) 'Impact parameter less than zero!'
  65         bfix=0.d0
  66        end if
  67        if (bfix.gt.RA2) then
  68         write(6,*) 'Impact parameter larger than two nuclear radius!'
  69         bfix=RA2
  70        end if
  71        b1=bfix
  72       else
  73        call bipsear(fmax1,xmin1)
  74        fmax=fmax1
  75        xmin=xmin1
  76  11    bb1=xmin*pyr(0)
  77        ff1=fmax*pyr(0)
  78        fb=funbip(bb1)
  79        if(ff1.gt.fb) goto 11
  80        b1=bb1
  81       end if
  82       bgen=b1
  83
  84 * calculate initial QGP temperature as function of centrality
  85       sb0=pi*RA*RA
  86       sb=RA*RA*(pi-2.d0*dasin(0.5d0*b1/RA))-b1*dsqrt(abs(RA*RA-
  87      >   b1*b1/4.d0))
  88       rtaa0=9.d0*AW*AW/(8.d0*sb0)
  89       br=max(1.d-10,b1*b1/(4.d0*RA*RA))
  90       call simpa(0.d0,20.d0,0.001d0,0.001d0,1.d-08,ftaa,xx,rest,
  91      >           aih,aiabs)
  92       rtaa=rtaa0*(1.d0-br*(1.d0+(1.d0-0.25d0*br)*dlog(1.d0/br)+
  93      >     4.d0*rest/pi))
  94       T00=T000*((rtaa*sb0)/(rtaa0*sb))**0.25d0
  95       T0=T00*(AW/207.d0)**0.166667d0
  96
  97 * generate single event with partonic energy loss
  98       nrg=0
  99       ehard=ckin(3)
 100       if(b1.le.1.85d0*RA) then
 101        call plinit(ehard)
 102        call plevnt(ehard)
 103       end if
 104
 105 * reset all in-vacuum radiated guark 4-momenta and codes to zero
 106       do i=1,1000
 107        do j=1,5
 108         pqua(i,j)=0.d0
 109         kqua(i,j)=0
 110        end do
 111       end do
 112       nrq=0
 113
 114 * generate final state shower in vacuum if it was excluded before
 115       nrgm=nrg                        ! fix number of in-medium emitted gluons
 116       ip1=0
 117       ip2=0
 118       ip01=0
 119       ip02=0
 120       ip001=0
 121       ip002=0
 122       if(mstj(41).ne.0) goto 5
 123       mstj(41)=1
 124       nn=n
 125       do i=9,nn
 126        if(k(i,3).eq.7) then
 127         ip1=i                    ! first hard parton (line ip1)
 128         kfh1=k(i,1)              ! status code of first hard parton
 129         qmax1=pyp(i,10)          ! transverse momentum of first hard parton
 130        end if
 131        if(k(i,3).eq.8) then
 132         ip2=i                    ! second hard parton (line ip2)
 133         kfh2=k(i,1)              ! status code of second hard parton
 134         qmax2=pyp(i,10)          ! transverse momentum of second hard parton
 135        end if
 136       end do
 137
 138       n1=n
 139       call pyshow(ip1,0,qmax1)    ! vacuum showering for first hard parton
 140       if(n.eq.n1) ip1=0
 141       n2=n
 142       call pyshow(ip2,0,qmax2)    ! vacuum showering for second hard parton
 143       if(n.eq.n2) ip2=0
 144       mstj(41)=0
 145       if(n.eq.nn) goto 5
 146
 147 * find two leading partons after showering
 148       do i=nn+1,n
 149        if(k(i,3).eq.ip1) ip001=i   ! first daughter of first hard parton
 150        if(k(i,3).eq.ip2) ip002=i   ! first daughter of second hard parton
 151       end do
 152       ptle1=0.d0
 153       ptle2=0.d0
 154       do i=nn+1,n
 155        if (k(i,1).eq.14) goto 3
 156        if(i.ge.ip002.and.ip002.gt.0) then
 157         ptl02=pyp(i,10)
 158         if(ptl02.gt.ptle2.and.k(i,2).eq.k(ip2,2)) then
 159          ip02=i                   ! leading parton in second shower (line ip02)
 160          ptle2=ptl02              ! pt of the leading parton
 161         end if
 162        elseif(ip001.gt.0) then
 163         ptl01=pyp(i,10)
 164         if(ptl01.gt.ptle1.and.k(i,2).eq.k(ip1,2)) then
 165          ip01=i                   ! leading parton in first shower (line ip01)
 166          ptle1=ptl01              ! pt of the leading parton
 167         end if
 168        end if
 169  3     continue
 170       end do
 171
 172 * replace two hard partons by two leading partons in original event record
 173       if(ip1.gt.0) then
 174        do j=1,5
 175         v(ip1,j)=v(ip01,j)
 176         p(ip1,j)=p(ip01,j)
 177        end do
 178        k(ip1,1)=kfh1
 179 * fix first/last daughter for moving entry
 180         do jgl=1,n
 181          if(k(jgl,4).eq.ip01) k(jgl,4)=ip1
 182          if(k(jgl,5).eq.ip01) k(jgl,5)=ip1
 183         end do
 184 *
 185       end if
 186       if(ip2.gt.0) then
 187        do j=1,5
 188         v(ip2,j)=v(ip02,j)
 189         p(ip2,j)=p(ip02,j)
 190        end do
 191        k(ip2,1)=kfh2
 192 * fix first/last daughter for moving entry
 193         do jgl=1,n
 194          if(k(jgl,4).eq.ip02) k(jgl,4)=ip2
 195          if(k(jgl,5).eq.ip02) k(jgl,5)=ip2
 196         end do
 197 *
 198       end if
 199
 200 * add final showering gluons to the list of in-medium emitted gluons,
 201 * fill the list of emitted quarks by final showering quark pairs,
 202 * and remove showering gluons and quarks from the event record
 203       do i=nn+1,n
 204        if(k(i,1).eq.14.or.i.eq.ip01.or.i.eq.ip02) goto 12
 205        if(k(i,2).ne.21) then               ! filling 'plquar' arrays for quarks
 206         nrq=nrq+1
 207         do j=1,5
 208          kqua(nrq,j)=k(i,j)
 209          pqua(nrq,j)=p(i,j)
 210         end do
 211          kqua(nrq,1)=2
 212         goto 12
 213        end if
 214        if(i.ge.ip002.and.ip002.gt.0) then
 215         ish=ip2
 216        else
 217         ish=ip1
 218        end if
 219        nrg=nrg+1
 220        nur=nrg
 221  7     ishm=kglu(nur-1,6)
 222        if(ish.ge.ishm.or.nur.le.2) goto 6   ! adding gluons in 'plglur' arrays
 223        do j=1,6
 224         kglu(nur,j)=kglu(nur-1,j)
 225        end do
 226        do j=1,4
 227         glur(nur,j)=glur(nur-1,j)
 228        end do
 229        nur=nur-1
 230        goto 7
 231  6     kglu(nur,1)=2                              ! status code
 232        kglu(nur,2)=k(i,2)                         ! particle identificator
 233        kglu(nur,3)=k(ish,3)                       ! parent line number
 234        kglu(nur,4)=0                              ! special colour info
 235        kglu(nur,5)=0                              ! special colour info
 236        kglu(nur,6)=ish                            ! associated parton number
 237        glur(nur,1)=p(i,4)                         ! energy
 238        glur(nur,2)=pyp(i,10)                      ! pt
 239        glur(nur,3)=pyp(i,15)                      ! phi
 240        glur(nur,4)=pyp(i,19)                      ! eta
 241  12    continue
 242        do j=1,5                              ! remove partons from event list
 243         v(i,j)=0.d0
 244         k(i,j)=0
 245         p(i,j)=0.d0
 246        end do
 247       end do
 248       n=nn
 249
 250  5    continue
 251
 252 * stop generate event if there are no additional gluons
 253       if(nrg.lt.1) goto 1
 254
 255 * define number of stirngs (ns) and number of entries in strings before
 256 * in-medium radiation (nis(ns))
 257       ns=0
 258       nes=0
 259       i0=0
 260       i1=0
 261       do i=1,500
 262        nis(i)=0
 263        nas(i)=0
 264        nss(i)=0
 265        nus(i)=0
 266       end do
 267       do i=9,n
 268        ks=k(i,1)
 269        ksp=k(i-1,1)
 270        if(ks.eq.2) then
 271         nis(ns+1)=nis(ns+1)+1
 272        elseif(ks.eq.1.and.nis(ns+1).gt.0) then
 273         nis(ns+1)=nis(ns+1)+1
 274         nes=nes+nis(ns+1)     ! nes - total number of entries
 275         nss(ns+1)=nes
 276         ns=ns+1
 277        elseif(ks.ne.2.and.ksp.ne.2.and.ns.gt.0) then
 278         i1=i1+1               ! last i1 lines not included in strings
 279        end if
 280       end do
 281       i0=n-nes-i1             ! first i0 lines not included in strings
 282       do i=1,ns
 283        nss(i)=nss(i)+i0
 284       end do
 285
 286 * move fragmented particles in bottom of event list
 287       i=i0+1
 288  2    ks=k(i,1)
 289       ksp=k(i-1,1)
 290       if(ks.ne.2.and.ksp.ne.2) then
 291        n=n+1
 292        do j=1,5
 293         v(n,j)=v(i,j)
 294         k(n,j)=k(i,j)
 295         p(n,j)=p(i,j)
 296        end do
 297 * fix first/last daughter for moving entry
 298        do jgl=1,n
 299         if(k(jgl,4).eq.i) k(jgl,4)=n
 300         if(k(jgl,5).eq.i) k(jgl,5)=n
 301        end do
 302 *
 303        do in=i+1,n
 304         do j=1,5
 305          v(in-1,j)=v(in,j)
 306          k(in-1,j)=k(in,j)
 307          p(in-1,j)=p(in,j)
 308         end do
 309 * fix first/last daughter for moving entry
 310         do jgl=1,n
 311          if(k(jgl,4).eq.in) k(jgl,4)=in-1
 312          if(k(jgl,5).eq.in) k(jgl,5)=in-1
 313         end do
 314 *
 315        end do
 316        do ip=1,nrg
 317         ku=kglu(ip,6)
 318         if(ku.gt.i) kglu(ip,6)=ku-1
 319        end do
 320        n=n-1
 321       else
 322        i=i+1
 323       end if
 324       if(i.le.n-i1) goto 2
 325
 326 * define number of additional entries in strings, nas(ns)
 327       do i=1,nrg
 328        kas=kglu(i,6)
 329        if(kas.le.nss(1)) then
 330         nas(1)=nas(1)+1
 331        else
 332         do j=2,ns
 333          if(kas.le.nss(j).and.kas.gt.nss(j-1))
 334      >   nas(j)=nas(j)+1
 335         end do
 336        end if
 337       end do
 338       do j=1,ns
 339        do i=1,j
 340         nus(j)=nus(j)+nas(i)
 341        end do
 342       end do
 343
 344 * add emitted gluons in event list
 345       nu=n
 346       n=n+nrg
 347       do i=nu,nu-i1,-1
 348        is=i+nrg
 349        do j=1,5
 350         v(is,j)=v(i,j)
 351         k(is,j)=k(i,j)
 352         p(is,j)=p(i,j)
 353        end do
 354 * fix first/last daughter for moving entries
 355        do jgl=1,n
 356         if(k(jgl,4).eq.i) k(jgl,4)=is
 357         if(k(jgl,5).eq.i) k(jgl,5)=is
 358        end do
 359 *
 360       end do
 361       do ia=ns-1,1,-1
 362        do i=nss(ia+1)-1,nss(ia),-1
 363         is=i+nus(ia)
 364         do j=1,5
 365          v(is,j)=v(i,j)
 366          k(is,j)=k(i,j)
 367          p(is,j)=p(i,j)
 368         end do
 369 * fix first/last daughter for moving entries
 370         do jgl=1,n
 371          if(k(jgl,4).eq.i) k(jgl,4)=is
 372          if(k(jgl,5).eq.i) k(jgl,5)=is
 373         end do
 374 *
 375        end do
 376       end do
 377
 378       do i=1,nrg
 379        if(i.le.nus(1)) then
 380        ia=nss(1)-1+i
 381        else
 382         do in=2,ns
 383          if(i.le.nus(in).and.i.gt.nus(in-1))
 384      >   ia=nss(in)-1+i
 385         end do
 386        end if
 387        eg=glur(i,1)
 388        ptg=glur(i,2)
 389        phig=glur(i,3)
 390        etag=glur(i,4)
 391        do j=1,5
 392         v(ia,j)=0.d0
 393         k(ia,j)=kglu(i,j)
 394        end do
 395        p(ia,1)=ptg*dcos(phig)
 396        p(ia,2)=ptg*dsin(phig)
 397        p(ia,3)=dsqrt(abs(eg*eg-ptg*ptg))
 398        if(etag.lt.0.d0) p(ia,3)=-1.d0*p(ia,3)
 399        p(ia,4)=dsqrt(ptg*ptg+p(ia,3)**2)
 400        p(ia,5)=0.d0
 401       end do
 402
 403 * rearrange partons to form strings in event list
 404       do ij=1,1000
 405        ijoin(ij)=0
 406       end do
 407       do i=1,ns
 408        njoin=nis(i)+nas(i)
 409        if(i.eq.1) then
 410         do j=1,njoin
 411          ijoin(j)=i0+j
 412         end do
 413        else
 414         do j=1,njoin
 415          ijoin(j)=nss(i-1)+nus(i-1)+j
 416         end do
 417        end if
 418        call pyjoin(njoin,ijoin)
 419       end do
 420
 421 * add in-vacuum emitted quark pairs
 422       if(nrq.lt.2) goto 1
 423       do i=1,nrq,2
 424        n=n+2
 425        do j=1,5
 426         v(n-1,j)=0.d0
 427         k(n-1,j)=kqua(i,j)
 428         p(n-1,j)=pqua(i,j)
 429        end do
 430        in=i+1
 431  4     ktest=k(n-1,2)+kqua(in,2)
 432        if(ktest.eq.0.or.in.eq.nrq) goto 8
 433        in=in+1
 434        goto 4
 435  8     do j=1,5
 436         v(n,j)=0.d0
 437         k(n,j)=kqua(in,j)
 438         p(n,j)=pqua(in,j)
 439        end do
 440        if(in.gt.i+1) then
 441         do j=1,5
 442          kqua(in,j)=kqua(i+1,j)
 443          pqua(in,j)=pqua(i+1,j)
 444         end do
 445        end if
 446       end do
 447
 448       do ij=1,2
 449        ijoik(ij)=0
 450       end do
 451       do i=1,nrq-1,2
 452        k(n+1-i,1)=1
 453        ijoik(1)=n-i
 454        ijoik(2)=n+1-i
 455        call pyjoin(2,ijoik)
 456       end do
 457
 458  1    continue
 459
 460       return
 461       end
 462
 463 ********************************* PLINIT ***************************
 464       SUBROUTINE PLINIT(ET)
 465 * set nucleus thikness and plasma parameters
 466       IMPLICIT DOUBLE PRECISION(A-H, O-Z)
 467       IMPLICIT INTEGER(I-N)
 468       INTEGER PYK,PYCHGE,PYCOMP
 469       external plvisc
 470       common /plpar1/ tau0,T0,TC,sigqq,AW,RA,mvisc,nf
 471       common /plpar2/ pln0,taupl,tauh,sigpl,sigh,sigplh,sigqqh,rg,rgn
 472       common /plevol/ taup(5000),temp(5000),denp(5000),enep(5000)
 473       save /plevol/
 474 *
 475       pi=3.14159d0
 476       pi2=pi*pi
 477
 478 * set number degrees of freedom in QGP
 479       hgd=3.d0
 480       rg=(16.d0+10.5d0*nf)/hgd
 481       rgn=(16.d0+9.d0*nf)/hgd
 482
 483 * set 'fiction' sigma for parton rescattering in QGP
 484       sigqq=4.2d0
 485       sigpl=2.25d0*2.25d0*sigqq*(16.d0+4.d0*nf)/(16.d0+9.d0*nf)
 486
 487 * set intial plasma temperature, density and energy density in perfect
 488 * (if mvisc=0) or viscous (mvisc=1,2) QGP with PLVISC subroitine
 489       hst=0.15d0
 490       if(T0.gt.1.5d0.or.mvisc.eq.2) hst=0.25d0
 491       if(T0.gt.1.5d0.and.mvisc.ne.0) hst=0.9d0
 492       T01=T0*5.06d0
 493       TC1=TC*5.06d0
 494       pln0=(16.d0+9.d0*nf)*1.2d0*(T01**3)/pi2
 495       ened0=pi2*(16.d0+10.5d0*nf)*(T01**4)/30.d0
 496       hh=hst*tau0
 497       tau=tau0                          ! proper time
 498       T=T01                             ! temperature
 499       den=pln0                          ! number density
 500       ened=ened0                        ! energy density
 501
 502 * create array of parameters to configurate QGP time evolution
 503       DO I=1,5000
 504        taup(i)=tau                      ! proper time
 505        temp(i)=T/5.06d0                 ! temperature
 506        denp(i)=den                      ! number density
 507        enep(i)=ened/5.06d0              ! energy density
 508        ened1=0.5d0*hh*(1.3333d0*plvisc(T)/(tau*tau)-1.3333d0
 509      >       *ened/tau)+ened
 510        T1=(30.d0*ened1/((16.d0+10.5d0*nf)*pi2))**0.25d0
 511        tau1=tau+0.5d0*hh
 512        ened=hh*(1.3333d0*plvisc(T1)/(tau1*tau1)-1.3333d0
 513      >      *ened1/tau1)+ened
 514        TPR=T
 515        T=(30.d0*ened/((16.d0+10.5d0*nf)*pi2))**0.25d0
 516        den=(16.d0+9.d0*nf)*1.2d0*(T**3)/pi2
 517        tau=tau+hh
 518        if(TPR.gt.TC1.and.T.le.TC1) taupl=tau-0.5d0*hh  ! QGP lifetime taupl
 519       END DO
 520       tauh=taupl*rg                                    ! mixed phase lifetime
 521
 522       return
 523       end
 524 ******************************** END PLINIT **************************
 525
 526 ******************************* PLEVNT ******************************
 527       SUBROUTINE PLEVNT(ET)
 528 * generate hard parton production vertex and passing with rescattering,
 529 * collisional and radiative energy loss of each parton through plasma
 530       IMPLICIT DOUBLE PRECISION(A-H, O-Z)
 531       IMPLICIT INTEGER(I-N)
 532       INTEGER PYK,PYCHGE,PYCOMP
 533       external plthik, pln, plt, pls, gauss, gluang
 534       external pyp,pyr,pyk
 535       common /plpar1/ tau0,T0,TC,sigqq,AW,RA,mvisc,nf
 536       common /plpar2/ pln0,taupl,tauh,sigpl,sigh,sigplh,sigqqh,rg,rgn
 537       common /thikpa/ fmax,xmin
 538       common /pyjets/ n,npad,k(4000,5),p(4000,5),v(4000,5)
 539       common /pyqpar/ T0u,tau0u,nfu,ienglu,ianglu
 540       common /plglur/ glur(1000,4),kglu(1000,6),nrg,nrgm
 541       common /factor/ cfac, kf
 542       common /pleave/ taul, temlev
 543       common /parimp/ b1, psib1, rb1, rb2, noquen
 544       common /plen/ epartc, um
 545       common /plos/ elr,rsk
 546       common /numje1/ nuj1, nuj2
 547       save /pyjets/, /plglur/
 548 *
 549       pi=3.14159d0
 550
 551 * find minimum of nuclear thikness function with subroutine plsear
 552       psib1=pi*(2.d0*pyr(0)-1.d0)
 553       call plsear (fmax1,xmin1)
 554       fmax=fmax1
 555       xmin=xmin1
 556
 557 * generate vertex of jet production
 558       iv=0
 559  1    rr1=xmin*pyr(0)
 560       ff1=fmax*pyr(0)
 561       f=plthik(rr1)
 562       iv=iv+1
 563       if(ff1.gt.f.and.iv.le.100000) goto 1
 564       r0=rr1
 565       rb1=dsqrt(abs(r0*r0+b1*b1/4.d0+r0*b1*dcos(psib1)))
 566       rb2=dsqrt(abs(r0*r0+b1*b1/4.d0-r0*b1*dcos(psib1)))
 567       rb1=max(rb1,1.d-4)
 568       rb2=max(rb2,1.d-4)
 569       if(noquen.eq.1) goto 7
 570
 571 * find maximum of angular spectrum of radiated gluons with subroutine gluang
 572       temin=0.5d0*pi
 573       temax=0.5d0*(1.d0+dsqrt(5.d0))*0.0863d0
 574       ftemax=gluang(temax)
 575
 576 * reset all radiated gluon 4-momenta and codes to zero -------------------
 577       do i=1,1000
 578        do j=1,4
 579         glur(i,j)=0.d0
 580         kglu(i,j)=0
 581        end do
 582        kglu(i,5)=0
 583        kglu(i,6)=0
 584       end do
 585       nrg=0
 586
 587 * generate changing 4-momentum of partons due to rescattering and energy loss
 588 * (for partons with |eta|<3.5 and pt>3 GeV/c)
 589       nuj1=9                            ! minimum line number of rescattered parton
 590       nuj2=n                            ! maximum line number of rescattered parton
 591       do 2 ip=nuj1,nuj2                 ! cycle on travelling partons
 592        irasf=0
 593        iraz=0
 594        ks=k(ip,1)                       ! parton status code
 595        kf=k(ip,2)                       ! parton identificator
 596        ka=abs(kf)
 597        ko=k(ip,3)                       ! origin (parent line number)
 598        epart=abs(pyp(ip,10))            ! parton transverse momentum
 599        etar=pyp(ip,19)                  ! parton pseudorapidity
 600        if(ko.gt.6.and.epart.ge.3.d0.and.abs(etar).
 601      >  le.3.5d0) then
 602        if(ka.eq.21.or.ka.eq.1.or.ka.eq.2.or.ka.eq.3.
 603      >  or.ka.eq.4.or.ka.eq.5.or.ka.eq.6.or.ka.eq.7.
 604      >  or.ka.eq.8) then
 605        if(ks.eq.2.or.ks.eq.1.or.ks.eq.21) then
 606         phir=pyp(ip,15)                 ! parton azimuthal angle
 607         tetr=pyp(ip,13)                 ! parton polar angle
 608         yrr=pyp(ip,17)                  ! parton rapidity
 609         stetr=dsin(tetr)                ! parton sin(theta)
 610         if(abs(stetr).le.1.d-05) then
 611          if(stetr.ge.0.d0) then
 612           stetr=1.d-05
 613          else
 614           stetr=-1.d-05
 615          end if
 616         end if
 617         phir1=-1.d0*phir
 618         tetr1=-1.d0*tetr
 619
 620 * set colour factor
 621         if(kf.eq.21) then
 622          cfac=1.d0                      ! for gluon
 623         else
 624          cfac=0.44444444d0              ! for quark
 625         end if
 626
 627 * boost from laboratory system to system of hard parton
 628         ipar=ip
 629         bet0=(r0*dcos(psib1)+0.5d0*b1)/rb1
 630         if(bet0.le.-1.d0) bet0=-0.99999d0
 631         if(bet0.ge.1.d0) bet0=0.99999d0
 632         bet=dacos(bet0)
 633         if(psib1.lt.0.d0) bet=-1.d0*bet
 634         phip=phir-bet
 635         if(phip.gt.pi) phip=phip-2.d0*pi
 636         if(phip.lt.-1.d0*pi) phip=phip+2.d0*pi
 637         call pyrobo(ip,ip,0.d0,phir1,0.d0,0.d0,0.d0)
 638         call pyrobo(ip,ip,tetr1,0.d0,0.d0,0.d0,0.d0)
 639
 640 * calculate proper time of parton leaving QGP
 641         aphin=(r0*r0-b1*b1/4.d0)/(rb1*rb2)
 642         if(aphin.le.-1.d0) aphin=-0.99999d0
 643         if(aphin.ge.1.d0) aphin=0.99999d0
 644         phin=dacos(aphin)
 645         if(psib1.le.0.d0) phin=-1.d0*phin
 646         phid=phip-phin
 647         if(phid.gt.pi) phid=phid-2.d0*pi
 648         if(phid.lt.-1.d0*pi) phid=phid+2.d0*pi
 649         taul1=abs(dsqrt(abs(RA*RA-(rb1*dsin(phip))**2))-rb1*dcos(phip))
 650         taul2=abs(dsqrt(abs(RA*RA-(rb2*dsin(phid))**2))-rb2*dcos(phid))
 651         taul=min(taul1,taul2)             ! escape time taul
 652         temlev=plt(taul)                  ! QGP temperature at taul
 653         if(taul.le.tau0) goto 100         ! escape from QGP if taul<tau0
 654
 655 * start parton rescattering in QGP with proper time iterations
 656         tau=tau0
 657  3      tfs=plt(tau)
 658         xi=-10.d0*dlog(max(pyr(0),1.d-10))/(sigpl*pln(tau))
 659         vel=abs(p(ip,3))/dsqrt(p(ip,3)**2+p(ip,5)**2) ! parton velocity
 660         if(vel.lt.0.3d0) goto 4
 661         tau=tau+xi*vel
 662         if(tau.ge.taul.or.tfs.le.TC) goto 100  ! escape if tau>taul or >taupl
 663
 664 * transform parton 4-momentum due to next scattering with subroutine pljetr
 665         epartc=p(ip,4)                         ! parton energy
 666         um=p(ip,5)                             ! parton mass
 667         sigtr=pls(tfs)*cfac*((p(ip,4)/pyp(ip,8))**2)
 668         prob=sigpl/(sigtr/stetr+sigpl)
 669         ran=pyr(0)
 670         irasf=irasf+1
 671         if(irasf.gt.100000) goto 100
 672         if(ran.lt.prob) goto 3
 673         pltp=plt(tau)
 674         pltp3=3.d0*pltp
 675         pass=50.6d0/(pln(tau)*sigtr)
 676         elr=0.d0
 677         rsk=0.d0
 678         call pljetr(tau,pass,pltp,ipar,epart)
 679         irasf=0
 680
 681 * set 4-momentum (in lab system) of next radiated gluon for parton number >8
 682 * and fill arrays of radiated gluons in common block plglur
 683         if(nrg.le.1000) then
 684          if(abs(elr).gt.0.1d0.and.ip.gt.8) then
 685 * generate the angle of emitted gluon
 686           if(ianglu.eq.0) then
 687  6         te1=temin*pyr(0)
 688            fte1=ftemax*pyr(0)
 689            fte2=gluang(te1)
 690            if(fte1.gt.fte2) goto 6
 691            tgl=te1
 692           elseif (ianglu.eq.1) then
 693            tgl=((0.5d0*pi*epartc)**pyr(0))/epartc
 694           else
 695            tgl=0.d0
 696           end if
 697           pgl=pi*(2.d0*pyr(0)-1.d0)
 698 * in comoving system
 699           pxgl=abs(elr)*stetr*(dcos(phir)*dcos(tgl)-
 700      >      dsin(phir)*dsin(tgl)*dsin(pgl))
 701           pygl=abs(elr)*stetr*(dsin(phir)*dcos(tgl)+
 702      >      dcos(phir)*dsin(tgl)*dsin(pgl))
 703           pzgl=-1.d0*abs(elr)*stetr*dsin(tgl)*dcos(pgl)
 704           ptgl=dsqrt(abs(pxgl*pxgl+pygl*pygl))
 705           psgl=dsqrt(abs(ptgl*ptgl+pzgl*pzgl))
 706 * recalculate in lab system
 707           dyg=0.5d0*dlog(max(1.d-9,(psgl+pzgl)/(psgl-pzgl)))
 708           pzgl=ptgl*dsinh(yrr+dyg)
 709           psgl=dsqrt(abs(ptgl*ptgl+pzgl*pzgl))
 710 *
 711           dpgl=pygl/pxgl
 712           glur1=abs(elr)                                       ! energy
 713           glur3=datan(dpgl)                                    ! phi
 714           if(pxgl.lt.0.d0) then
 715            if(pygl.ge.0.d0) then
 716             glur3=glur3+pi
 717            else
 718             glur3=glur3-pi
 719            end if
 720           end if
 721           glur4=0.5d0*dlog(max(1.d-9,(psgl+pzgl)/(psgl-pzgl))) ! eta
 722           glur2=glur1/dcosh(glur4)                             ! pt
 723
 724 * put in event list radiated gluons with pt > 0.2 GeV only
 725           if(glur2.ge.0.2d0) then
 726            nrg=nrg+1
 727 * set gluon 4-momentum
 728            glur(nrg,1)=glur1                                   ! energy
 729            glur(nrg,2)=glur2                                   ! pt
 730            glur(nrg,3)=glur3                                   ! phi
 731            glur(nrg,4)=glur4                                   ! eta
 732 * set gluon codes
 733            kglu(nrg,1)=2                           ! status code
 734            kglu(nrg,2)=21                          ! particle identificator
 735            kglu(nrg,3)=k(ipar,3)                   ! parent line number
 736            kglu(nrg,4)=0                           ! special colour info
 737            kglu(nrg,5)=0                           ! special colour info
 738            kglu(nrg,6)=ipar                        ! associated parton number
 739           end if
 740          end if
 741         else
 742          write(6,*) 'Warning! Number of emitted gluons is too large!'
 743         end if
 744
 745 * set parton "thermalization" if pt<T
 746         if(abs(p(ip,3)).gt.pltp3) goto 3
 747  4      continue
 748         if(p(ip,3).ge.0.d0) then
 749          sigp=1.d0
 750         else
 751          sigp=-1.d0
 752         end if
 753  5      iraz=iraz+1
 754         if(iraz.gt.100000) goto 100
 755         ep0=-0.15d0*(dlog(max(1.d-10,pyr(0)))+dlog(max(1.d-10,pyr(0)))+
 756      >  dlog(max(1.d-10,pyr(0))))
 757         if(ep0.le.p(ip,5).or.ep0.ge.100.d0) goto 5
 758         pp0=dsqrt(abs(ep0**2-p(ip,5)**2))
 759         probt=pp0/ep0
 760         if(pyr(0).gt.probt) goto 5
 761         ctp0=2.d0*pyr(0)-1.d0
 762         stp0=dsqrt(abs(1.d0-ctp0**2))
 763         php0=pi*(2.d0*pyr(0)-1.d0)
 764         p(ip,1)=pp0*stp0*dcos(php0)
 765         p(ip,2)=pp0*stp0*dsin(php0)
 766         p(ip,3)=sigp*pp0*ctp0
 767         p(ip,4)=dsqrt(p(ip,1)**2+p(ip,2)**2+p(ip,3)**2+p(ip,5)**2)
 768
 769 * boost to laboratory system
 770  100    call pyrobo(ip,ip,tetr,phir,0.d0,0.d0,0.d0)
 771        end if
 772        end if
 773        end if
 774  2    continue
 775  7    continue
 776
 777       return
 778       end
 779 ******************************* END PLEVNT *************************
 780
 781 ******************************* PLJETR *****************************
 782       SUBROUTINE PLJETR(tau,y,x,ip,epart)
 783 * transform parton 4-momentum due to scattering in plasma at time = tau
 784       IMPLICIT DOUBLE PRECISION(A-H, O-Z)
 785       IMPLICIT INTEGER(I-N)
 786       INTEGER PYK,PYCHGE,PYCOMP
 787       external plfun1, pls
 788       external pyp,pyr
 789       common /plpar1/ tau0,T0,TC,sigqq,AW,RA,mvisc,nf
 790       common /plpar2/ pln0,taupl,tauh,sigpl,sigh,sigplh,sigqqh,rg,rgn
 791       common /pyjets/ n,npad,k(4000,5),p(4000,5),v(4000,5)
 792       common /pyqpar/ T0u,tau0u,nfu,ienglu,ianglu
 793       common /pljdat/ ej, z, ygl, alfs, um, epa
 794       common /pleave/ taul, temlev
 795       common /radcal/ aa, bb
 796       common /factor/ cfac, kf
 797       common /plos/ elr,rsk
 798       save /pyjets/
 799 *
 800       pi=3.14159d0
 801       spi=dsqrt(pi)
 802       tauu=x                            ! redenote temerature tauu=x
 803       i=ip                              ! redenote parton number i=ip
 804       iter=0
 805       iraz=0
 806
 807 * boost to system of comoving plasma constituent
 808       phir=pyp(i,15)                    ! parton phi
 809       tetr=pyp(i,13)                    ! parton theta
 810       phir1=-1.d0*phir
 811       tetr1=-1.d0*tetr
 812       call pyrobo(i,i,0.d0,phir1,0.d0,0.d0,0.d0)
 813       call pyrobo(i,i,tetr1,0.d0,0.d0,0.d0,0.d0)
 814       pp=pyp(i,8)                       ! parton total momentum
 815       ppl=abs(p(i,3))                   ! parton pz
 816       um=p(i,5)                         ! parton mass
 817       epa=p(i,4)                        ! parton energy
 818       ppt=pyp(i,10)                     ! parton pt
 819       pphi=pyp(i,15)                    ! parton phi
 820
 821       if(ppl.lt.3.d0) goto 222          ! no energy loss if pz<3 GeV/c
 822
 823 * generation hard parton-plasma scattering with momentum transfer rsk
 824  221   ep0=-1.*tauu*(dlog(max(1.d-10,pyr(0)))+dlog(max(1.d-10,
 825      >   pyr(0)))+dlog(max(1.d-10,pyr(0))))     ! energy of 'thermal' parton
 826        iter=iter+1
 827        if(ep0.lt.1.d-10.and.iter.le.100000) goto 221
 828        scm=2.*ep0*epa+um*um+ep0*ep0
 829        qm2=(scm-((um+ep0)**2))*(scm-((um-ep0)**2))/scm
 830        bub=4.d0*tauu/TC
 831        alf=6.d0*pi/((33.d0-2.d0*nf)*dlog(max(bub,1.d-10)))
 832        z=pi*4.d0*tauu*tauu*alf*(1.+nf/6.d0)
 833        bubs=dsqrt(abs(z))/TC
 834        alfs=6.d0*pi/((33.d0-2.d0*nf)*dlog(max(bubs,1.d-10)))
 835        phmin2=z
 836        phmax2=max(phmin2,qm2)
 837        fqmax2=1.d0/(dlog(max(phmin2/(TC*TC),1.d-10)))**2
 838  12    rn1=pyr(0)
 839        tp=1.d0/(rn1/phmax2+(1.d0-rn1)/phmin2)
 840        ftp=1.d0/(dlog(max(tp/(TC*TC),1.d-10)))**2
 841        fprob=ftp/fqmax2
 842        rn2=pyr(0)
 843        if(fprob.lt.rn2) goto 12
 844        rsk=dsqrt(abs(tp))
 845        if(rsk.gt.ppl) rsk=ppl
 846
 847 * calculate radiative energy loss per given scattering with subroutine plfun1
 848        ygl=y*cfac                      ! mean gluon free path in GeV^{-1}
 849        elp=ygl*z                       ! mimimum radiated energy in LPM regime
 850        ej=ppl
 851        bb=ej                           ! maximum radiated energy
 852        bbi=max(dsqrt(abs(z)),1.000001d0*elp)
 853        aa=min(bb,bbi)                  ! minimum radiated energy
 854        hh=0.00001d0*(bb-aa)
 855        REPS=0.01d0
 856        AEPS=1.d-8
 857        CALL SIMPA(aa,bb,hh,REPS,AEPS,plfun1,om,resun,AIH,AIABS)
 858 *                                      ! integral over omega for radiative loss
 859        call radsear(ermax1,eomin1)
 860        ermax=ermax1
 861        eomin=eomin1
 862  11    resu=eomin*pyr(0)+aa
 863        fres=ermax*pyr(0)
 864        fres1=plfun1(resu)
 865        iraz=iraz+1
 866        if(fres.gt.fres1.and.iraz.lt.100000) goto 11
 867        elr=resu*resun                   ! energy of radiated gluon
 868
 869 * to chancel radiative energy loss (optional case)
 870        if(ienglu.eq.2) elr=0.d0
 871 * to chancel collisional energy loss (optional case)
 872        if(ienglu.eq.1) rsk=0.d0
 873
 874 * determine the direction of parton moving
 875        if(p(i,3).ge.0.d0) then
 876         sigp=1.d0
 877        else
 878         sigp=-1.d0
 879        end if
 880
 881 * calculate new 4-momentum of hard parton
 882        phirs=2.d0*pi*pyr(0)
 883        epan=epa-rsk*rsk/(2.d0*ep0)-abs(elr)
 884        if(epan.lt.0.1d0) then
 885         epan=epan+abs(elr)
 886         elr=0.d0
 887         if(epan.lt.0.1d0) then
 888          rsk=0.d0
 889          epan=epa
 890         end if
 891        end if
 892        pptn=dsqrt(abs(rsk*rsk+(rsk**4)*(1.d0-epa*epa/(ppl*ppl))/
 893      >      (4.d0*ep0*ep0)-(rsk**4)*epa/(2.d0*ep0*ppl*ppl)-(rsk**4)/
 894      >      (4.d0*ppl*ppl)))
 895        ppln=dsqrt(abs(epan*epan-pptn*pptn-p(i,5)**2))
 896        p(i,1)=pptn*dcos(phirs)                                 ! px
 897        p(i,2)=pptn*dsin(phirs)                                 ! py
 898        p(i,3)=sigp*ppln                                        ! pz
 899        p(i,4)=dsqrt(p(i,1)**2+p(i,2)**2+p(i,3)**2+p(i,5)**2)   ! E
 900 * boost to system of hard parton
 901  222   call pyrobo(i,i,tetr,phir,0.d0,0.d0,0.d0)
 902
 903       return
 904       end
 905 ******************************* END PLJETR **************************
 906
 907 ******************************** PLSEAR ***************************
 908        SUBROUTINE PLSEAR (fmax,xmin)
 909 * finding maximum and 'sufficient minimum of nucleus thikness function.
 910 * xm, fm are outputs.
 911        IMPLICIT DOUBLE PRECISION(A-H, O-Z)
 912        IMPLICIT INTEGER(I-N)
 913        INTEGER PYK,PYCHGE,PYCOMP
 914        external plthik
 915        common /parimp/ b1, psib1, rb1, rb2, noquen
 916        common /plpar1/ tau0,T0,TC,sigqq,AW,RA,mvisc,nf
 917
 918        rm1=dsqrt(abs(RA*RA-b1*b1/4.d0*(dsin(psib1)**2)))+
 919      >     b1*dcos(psib1)/2.d0
 920        rm2=dsqrt(abs(RA*RA-b1*b1/4.d0*(dsin(psib1)**2)))-
 921      >     b1*dcos(psib1)/2.d0
 922        xmin=min(rm1,rm2)
 923        fmax=0.d0
 924        do 10 j=1,1000
 925         x=xmin*(j-1)/999.d0
 926         f=plthik(x)
 927         if(f.gt.fmax) then
 928          fmax=f
 929         end if
 930   10   continue
 931
 932        return
 933        end
 934 ****************************** END PLSEAR **************************
 935
 936 ******************************** RADSEAR ***************************
 937        SUBROUTINE RADSEAR (fmax,xmin)
 938 * find maximum and 'sufficient minimum of radiative energy loss distribution
 939 * xm, fm are outputs.
 940       IMPLICIT DOUBLE PRECISION(A-H, O-Z)
 941       IMPLICIT INTEGER(I-N)
 942       INTEGER PYK,PYCHGE,PYCOMP
 943       external plfun1
 944       common /radcal/ aa, bb
 945
 946       xmin=bb-aa
 947       fmax=0.d0
 948       do j=1,1000
 949        x=aa+xmin*(j-1)/999.d0
 950        f=plfun1(x)
 951        if(f.gt.fmax) then
 952         fmax=f
 953        end if
 954       end do
 955
 956       return
 957       end
 958 ****************************** END RADSEAR **************************
 959
 960 ********************************* BIPSEAR ***************************
 961       SUBROUTINE BIPSEAR (fmax,xmin)
 962 * find maximum and 'sufficient minimum' of jet production cross section
 963 * as a function of impact paramater (xm, fm are outputs)
 964       IMPLICIT DOUBLE PRECISION(A-H, O-Z)
 965       IMPLICIT INTEGER(I-N)
 966       INTEGER PYK,PYCHGE,PYCOMP
 967       external funbip
 968       common /plpar1/ tau0,T0,TC,sigqq,AW,RA,mvisc,nf
 969
 970       xmin=2.*RA
 971       fmax=0.d0
 972       do j=1,1000
 973        x=xmin*(j-1)/999.d0
 974        f=funbip(x)
 975        if(f.gt.fmax) then
 976         fmax=f
 977        end if
 978       end do
 979
 980       return
 981       end
 982 ****************************** END RADSEAR **************************
 983
 984 **************************** SIMPA **********************************
 985       SUBROUTINE SIMPA (A1,B1,H1,REPS1,AEPS1,FUNCT,X,
 986      1                     AI,AIH,AIABS)
 987 * calculate intergal of function FUNCT(X) on the interval from A1 to B1
 988       IMPLICIT DOUBLE PRECISION(A-H, O-Z)
 989       IMPLICIT INTEGER(I-N)
 990       DIMENSION F(7), P(5)
 991       H=dSIGN ( H1, B1-A1 )
 992       S=dSIGN (1.d0, H )
 993       A=A1
 994       B=B1
 995       AI=0.0d0
 996       AIH=0.0d0
 997       AIABS=0.0d0
 998       P(2)=4.d0
 999       P(4)=4.d0
1000       P(3)=2.d0
1001       P(5)=1.d0
1002       IF(B-A)1,2,1
1003  1    REPS=ABS(REPS1)
1004       AEPS=ABS(AEPS1)
1005       DO 3 K=1,7
1006  3    F(K)=10.d16
1007       X=A
1008       C=0.d0
1009       F(1)=FUNCT(X)/3.d0
1010  4    X0=X
1011       IF( (X0+4.d0*H-B)*S)5,5,6
1012  6    H=(B-X0)/4.d0
1013       IF ( H ) 7,2,7
1014  7    DO 8 K=2,7
1015  8    F(K)=10.d16
1016       C=1.d0
1017  5    DI2=F (1)
1018       DI3=ABS( F(1) )
1019       DO 9 K=2,5
1020       X=X+H
1021       IF((X-B)*S)23,24,24
1022  24   X=B
1023  23   IF(F(K)-10.d16)10,11,10
1024  11   F(K)=FUNCT(X)/3.d0
1025  10   DI2=DI2+P(K)*F(K)
1026  9    DI3=DI3+P(K)*ABS(F(K))
1027       DI1=(F(1)+4.*F(3)+F(5))*2.d0*H
1028       DI2=DI2*H
1029       DI3=DI3*H
1030       IF (REPS) 12,13,12
1031  13   IF (AEPS) 12,14,12
1032  12   EPS=ABS((AIABS+DI3)*REPS)
1033       IF(EPS-AEPS)15,16,16
1034  15   EPS=AEPS
1035  16   DELTA=ABS(DI2-DI1)
1036       IF(DELTA-EPS)20,21,21
1037  20   IF(DELTA-EPS/8.d0)17,14,14
1038  17   H=2.d0*H
1039       F(1)=F(5)
1040       F(2)=F(6)
1041       F(3)=F(7)
1042       DO 19 K=4,7
1043  19   F(K)=10.d16
1044       GO TO 18
1045  14   F(1)=F(5)
1046       F(3)=F(6)
1047       F(5)=F(7)
1048       F(2)=10.d16
1049       F(4)=10.d16
1050       F(6)=10.d16
1051       F(7)=10.d16
1052  18   DI1=DI2+(DI2-DI1)/15.d0
1053       AI=AI+DI1
1054       AIH=AIH+DI2
1055       AIABS=AIABS+DI3
1056       GO TO 22
1057  21   H=H/2.d0
1058       F(7)=F(5)
1059       F(6)=F(4)
1060       F(5)=F(3)
1061       F(3)=F(2)
1062       F(2)=10.d16
1063       F(4)=10.d16
1064       X=X0
1065       C=0.
1066       GO TO 5
1067  22   IF(C)2,4,2
1068  2    RETURN
1069       END
1070 ************************* END SIMPA *******************************
1071
1072 ************************* PARINV **********************************
1073       SUBROUTINE PARINV(X,A,F,N,R)
1074 * gives interpolation of function F(X) with  arrays A(N) and F(N)
1075       IMPLICIT DOUBLE PRECISION(A-H, O-Z)
1076       IMPLICIT INTEGER(I-N)
1077       DIMENSION A(N),F(N)
1078       IF(X.LT.A(1))GO TO 11
1079       IF(X.GT.A(N))GO TO 4
1080       K1=1
1081       K2=N
1082  2    K3=K2-K1
1083       IF(K3.LE.1)GO TO 6
1084       K3=K1+K3/2
1085       IF(A(K3)-X) 7,8,9
1086  7    K1=K3
1087       GOTO2
1088  9    K2=K3
1089       GOTO2
1090  8    P=F(K3)
1091       RETURN
1092  3    B1=A(K1)
1093       B2=A(K1+1)
1094       B3=A(K1+2)
1095       B4=F(K1)
1096       B5=F(K1+1)
1097       B6=F(K1+2)
1098       R=B4*((X-B2)*(X-B3))/((B1-B2)*(B1-B3))+B5*((X-B1)*(X-B3))/
1099      1 ((B2-B1)*(B2-B3))+B6*((X-B1)*(X-B2))/((B3-B1)*(B3-B2))
1100       RETURN
1101  6    IF(K2.NE.N)GO TO 3
1102       K1=N-2
1103       GOTO3
1104  4    C=ABS(X-A(N))
1105       IF(C.LT.0.1d-7) GO TO 5
1106       K1=N-2
1107  13   CONTINUE
1108 C13   PRINT 41,X
1109 C41   FORMAT(25H X IS OUT OF THE INTERVAL,3H X=,F15.9)
1110       GO TO 3
1111  5    R=F(N)
1112       RETURN
1113  11   C=ABS(X-A(1))
1114       IF(C.LT.0.1d-7) GO TO 12
1115       K1=1
1116       GOTO 13
1117  12   R=F(1)
1118       RETURN
1119       END
1120 C************************** END PARINV *************************************
1121
1122 * function to calculate quark-quark scattering differential cross section
1123        double precision FUNCTION PLSIGH(Z)
1124        IMPLICIT DOUBLE PRECISION(A-H, O-Z)
1125        IMPLICIT INTEGER(I-N)
1126        INTEGER PYK,PYCHGE,PYCOMP
1127        common /plpar1/ tau0,T0,TC,sigqq,AW,RA,mvisc,nf
1128        pi=3.14159d0
1129        beta=(33.d0-2.d0*nf)/(12.d0*pi)
1130        alfs=1.d0/(beta*dlog(max(1.d-10,z/(TC*TC))))
1131        PLSIGH=8.d0*pi*alfs*alfs/(9.d0*z*z)
1132        return
1133        end
1134
1135 * function to calculate differential radiated gluon spectrum in BDMS model
1136       double precision FUNCTION PLFUN1(or)
1137       IMPLICIT DOUBLE PRECISION(A-H, O-Z)
1138       IMPLICIT INTEGER(I-N)
1139       INTEGER PYK,PYCHGE,PYCOMP
1140       common /plpar1/ tau0,T0,TC,sigqq,AW,RA,mvisc,nf
1141       common /pljdat/ ej, z, ygl, alfs, um, epa
1142       common /pleave/ taul, temlev
1143       common /factor/ cfac, kf
1144       pi=3.14159d0
1145       x=min((1.d0-ygl*z/or),or/ej)
1146       if(x.le.0.d0) x=0.d0
1147       if(x.ge.1.d0) x=0.9999d0
1148       if(kf.eq.21) then
1149        if(x.ge.0.5d0) x=1.-x
1150        spinf=0.5*(1.+(1.d0-x)**4+x**4)/(1.d0-x)
1151       else
1152        spinf=1.d0-x+0.5d0*x*x
1153       end if
1154       ak=ygl*z/(or*(1.d0-x))
1155       al=taul*5.06d0
1156       uu=0.5d0*al*dsqrt(abs(0.5d0*(1.d0-x+cfac*x*x)*ak*
1157      >   dlog(max(16.d0/ak,1.d-10))))/ygl
1158 * if  quark production outside the QGP then
1159 * arg=(((dsin(uu)*cosh(uu))**2)+((dcos(uu)*sinh(uu))**2))/(2.d0*uu*uu);
1160 * here quark production inside the QGP
1161       arg=((dcos(uu)*cosh(uu))**2)+((dsin(uu)*sinh(uu))**2)
1162       gl1=(ygl/(cfac*z))**0.3333333d0
1163       gl2=(um/epa)**1.333333d0
1164       dc=1.d0/((1.d0+((gl1*gl2*or)**1.5d0))**2)       ! massive parton
1165 c      dc=1.d0                                         !massless parton
1166       plfun1=dc*3.d0*alfs*ygl*dlog(max(arg,1.d-20))*spinf/(pi*al*or)
1167       return
1168       end
1169
1170 * function to calculate time-dependence of QGP viscosity (if mvisc=1,2)
1171       double precision FUNCTION PLVISC(X)
1172       IMPLICIT DOUBLE PRECISION(A-H, O-Z)
1173       IMPLICIT INTEGER(I-N)
1174       INTEGER PYK,PYCHGE,PYCOMP
1175       common /plpar1/ tau0,T0,TC,sigqq,AW,RA,mvisc,nf
1176       pi=3.14159d0
1177       T=X
1178       TC1=5.06d0*TC
1179       if(X.le.TC1) T=TC1
1180       if(mvisc.eq.0) then
1181        c=0.d0
1182       elseif(mvisc.eq.1) then
1183        a=3.4d0*(1.d0+0.12d0*(2.d0*nf+1.d0))
1184        b=15.d0*(1.d0+0.06d0*nf)
1185        c=4.d0*pi*pi*(10.5d0*nf/a+16.d0/b)/675.d0
1186       else
1187        c=(1.7d0*nf+1.d0)*0.342d0/(1.d0+nf/6.d0)
1188       end if
1189       bub=4.d0*T/TC1
1190       alf=6.d0*pi/((33.d0-2.d0*nf)*dlog(max(bub,1.d-10)))
1191       alf1=1.d0/alf
1192       PLVISC=c*(T**3)/(alf*alf*dlog(max(1.d-10,alf1)))
1193       return
1194       end
1195
1196 * function to calculate time-dependence of QGP number density
1197        double precision FUNCTION PLN(X)
1198        IMPLICIT DOUBLE PRECISION(A-H, O-Z)
1199        IMPLICIT INTEGER(I-N)
1200        INTEGER PYK,PYCHGE,PYCOMP
1201        common /plpar1/ tau0,T0,TC,sigqq,AW,RA,mvisc,nf
1202        common /plpar2/ pln0,taupl,tauh,sigpl,sigh,sigplh,sigqqh,rg,rgn
1203        common /plevol/ taup(5000),temp(5000),denp(5000),enep(5000)
1204        save /plevol/
1205        pi2=3.14159d0*3.14159d0
1206        t=X
1207        if(t.lt.taupl) then
1208         call parinv(t,taup,denp,5000,res)
1209        else
1210         res=1.2d0*(16.d0+9.d0*nf)*((5.06d0*TC)**3)/pi2
1211        end if
1212        PLN=res
1213        return
1214        end
1215
1216 * function to calculate time-dependence of QGP temperature
1217        double precision FUNCTION PLT(X)
1218        IMPLICIT DOUBLE PRECISION(A-H, O-Z)
1219        IMPLICIT INTEGER(I-N)
1220        INTEGER PYK,PYCHGE,PYCOMP
1221        common /plpar1/ tau0,T0,TC,sigqq,AW,RA,mvisc,nf
1222        common /plpar2/ pln0,taupl,tauh,sigpl,sigh,sigplh,sigqqh,rg,rgn
1223        common /plevol/ taup(5000),temp(5000),denp(5000),enep(5000)
1224        save /plevol/
1225        t=X
1226        if(t.lt.taupl) then
1227         call parinv(t,taup,temp,5000,res)
1228        else
1229         res=TC
1230        end if
1231         PLT=res
1232        return
1233        end
1234
1235 * function to caculate time-dependence of parton-plasma integral cross section
1236        double precision FUNCTION PLS(X)
1237        IMPLICIT DOUBLE PRECISION(A-H, O-Z)
1238        IMPLICIT INTEGER(I-N)
1239        INTEGER PYK,PYCHGE,PYCOMP
1240        external plsigh
1241        common /plpar1/ tau0,T0,TC,sigqq,AW,RA,mvisc,nf
1242        common /plpar2/ pln0,taupl,tauh,sigpl,sigh,sigplh,sigqqh,rg,rgn
1243        common /plen/ epartc, um
1244        t=X
1245        pi=3.14159d0
1246        bub=4.d0*t/TC
1247        alf=6.d0*pi/((33.d0-2.d0*nf)*dlog(max(bub,1.d-10)))
1248        ZZ0=4.d0*t*t*pi*alf*(1.d0+nf/6.d0)
1249        scm=4.d0*t*epartc+um*um+4.d0*t*t
1250        ZZ1=max((scm-((um+2.d0*t)**2))*(scm-((um-2.d0*t)**2))/scm,ZZ0)
1251        HH1=0.01d0*ZZ1
1252        REPS=0.01d0
1253        AEPS=1.d-8
1254        CALL SIMPA(ZZ0,ZZ1,HH1,REPS,AEPS,plsigh,ZZ,RESS,AIH,AIABS)
1255        PLS=0.39d0*2.25d0*2.25d0*RESS*(16.d0+4.d0*nf)/(16.d0+9.d0*nf)
1256        return
1257        end
1258
1259 * function to calculate nuclear thikness function
1260        double precision FUNCTION PLTHIK(X)
1261        IMPLICIT DOUBLE PRECISION(A-H, O-Z)
1262        IMPLICIT INTEGER(I-N)
1263        INTEGER PYK,PYCHGE,PYCOMP
1264        common /parimp/ b1, psib1, rb1, rb2, noquen
1265        common /plpar1/ tau0,T0,TC,sigqq,AW,RA,mvisc,nf
1266        bu=X
1267        r12=bu*bu+b1*b1/4.d0+bu*b1*dcos(psib1)
1268        r22=bu*bu+b1*b1/4.d0-bu*b1*dcos(psib1)
1269        PLTHIK=dsqrt(abs((RA*RA-r12)*(RA*RA-r22)))*bu
1270        return
1271        end
1272
1273 * function to generate gauss distribution
1274       double precision function gauss(x0,sig)
1275       IMPLICIT DOUBLE PRECISION(A-H, O-Z)
1276       IMPLICIT INTEGER(I-N)
1277  41   u1=pyr(0)
1278       u2=pyr(0)
1279       v1=2.d0*u1-1.d0
1280       v2=2.d0*u2-1.d0
1281       s=v1**2+v2**2
1282       if(s.gt.1) go to 41
1283       gauss=v1*dsqrt(-2.d0*dlog(s)/s)*sig+x0
1284       return
1285       end
1286
1287 * function to calculate angular distribution of emitted gluons
1288       double precision function gluang(x)
1289       IMPLICIT DOUBLE PRECISION(A-H, O-Z)
1290       IMPLICIT INTEGER(I-N)
1291       s=0.0863d0
1292       gluang=x*dexp(-1.d0*(x-s)*(x-s)/(2.d0*s*s))
1293       return
1294       end
1295
1296 * function to calculate jet production vs. centrality
1297       double precision function funbip(x)
1298       IMPLICIT DOUBLE PRECISION(A-H, O-Z)
1299       IMPLICIT INTEGER(I-N)
1300       INTEGER PYK,PYCHGE,PYCOMP
1301       common /plpar1/ tau0,T0,TC,sigqq,AW,RA,mvisc,nf
1302       dimension bip(15), bipr(15), pjet(15)
1303       data bip/0.d0,0.5d0,1.5d0,2.5d0,3.5d0,4.5d0,5.5d0,6.5d0,7.5d0,
1304      >         8.5d0,9.5d0,10.5d0,11.5d0,12.5d0,13.5d0/
1305       data pjet/200000.d0,217558.d0,625570.d0,949850.d0,1.17128d+06,
1306      >   1.30123d+06,1.32297d+06,1.18483d+06,1.02584d+06,839982.d0,
1307      >   621238.d0,399300.d0,227456.d0,113982.d0,41043.d0/
1308       bu=x
1309       do i=1,15
1310        bipr(i)=bip(i)*RA/6.8d0
1311       end do
1312       call parinv (bu,bipr,pjet,15,res)
1313       funbip=res
1314       return
1315       end
1316
1317 * function integrated at calculation of initial QGP temperature vs. centrality
1318       double precision function ftaa(x)
1319       IMPLICIT DOUBLE PRECISION(A-H, O-Z)
1320       IMPLICIT INTEGER(I-N)
1321       INTEGER PYK,PYCHGE,PYCOMP
1322       common/bintaa/ br
1323       a=1.d0+x*x
1324       ftaa=(1.d0-br*x*x/a)*dlog(1.d0+x*x*(1.d0-br))/(a*a)
1325       return
1326       end
1327 **************************************************************************