THydjet/hydjet1_1/pyquen1_1.f

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