TUHKMgen/PYQUEN/pyquen1_5.f

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