PYTHIA6/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,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, 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        if (z < 0.) write(6,*) "Warning in PLJETR: z < 0.", z
 900        bubs=dsqrt(abs(z))/TC
 901        alfs=6.d0*pi/((33.d0-2.d0*nf)*dlog(max(bubs,1.d-10)))
 902        phmin2=z
 903        phmax2=max(phmin2,qm2)
 904        fqmax2=1.d0/(dlog(max(phmin2/(TC*TC),1.d-10)))**2
 905  12    rn1=pyr(0)
 906        tp=1.d0/(rn1/phmax2+(1.d0-rn1)/phmin2)
 907        ftp=1.d0/(dlog(max(tp/(TC*TC),1.d-10)))**2
 908        fprob=ftp/fqmax2
 909        rn2=pyr(0)
 910        if(fprob.lt.rn2) goto 12
 911        rsk=dsqrt(abs(tp))
 912        if(rsk.gt.ppl) rsk=ppl
 913
 914 * calculate radiative energy loss per given scattering with subroutine plfun1
 915        ygl=y*cfac                      ! mean gluon free path in GeV^{-1}
 916        elp=ygl*z                       ! mimimum radiated energy in LPM regime
 917        ej=ppl
 918        bb=ej                           ! maximum radiated energy
 919        bbi=max(dsqrt(abs(z)),1.000001d0*elp)
 920        aa=min(bb,bbi)                  ! minimum radiated energy
 921        hh=0.00001d0*(bb-aa)
 922        REPS=0.01d0
 923        AEPS=1.d-8
 924        CALL SIMPA(aa,bb,hh,REPS,AEPS,plfun1,om,resun,AIH,AIABS)
 925 *                                      ! integral over omega for radiative loss
 926        call radsear(ermax1,eomin1)
 927        ermax=ermax1
 928        eomin=eomin1
 929  11    resu=eomin*pyr(0)+aa
 930        fres=ermax*pyr(0)
 931        fres1=plfun1(resu)
 932        iraz=iraz+1
 933        if(fres.gt.fres1.and.iraz.lt.100000) goto 11
 934        elr=resu*resun                   ! energy of radiated gluon
 935
 936 * to chancel radiative energy loss (optional case)
 937        if(ienglu.eq.2) elr=0.d0
 938 * to chancel collisional energy loss (optional case)
 939        if(ienglu.eq.1) rsk=0.d0
 940
 941 * determine the direction of parton moving
 942        if(p(i,3).ge.0.d0) then
 943         sigp=1.d0
 944        else
 945         sigp=-1.d0
 946        end if
 947
 948 * calculate new 4-momentum of hard parton
 949        phirs=2.d0*pi*pyr(0)
 950        epan=epa-rsk*rsk/(2.d0*ep0)-abs(elr)
 951        if(epan.lt.0.1d0) then
 952         epan=epan+abs(elr)
 953         elr=0.d0
 954         if(epan.lt.0.1d0) then
 955          rsk=0.d0
 956          epan=epa
 957         end if
 958        end if
 959        pptn=dsqrt(abs(rsk*rsk+(rsk**4)*(1.d0-epa*epa/(ppl*ppl))/
 960      >      (4.d0*ep0*ep0)-(rsk**4)*epa/(2.d0*ep0*ppl*ppl)-(rsk**4)/
 961      >      (4.d0*ppl*ppl)))
 962        ppln=dsqrt(abs(epan*epan-pptn*pptn-p(i,5)**2))
 963        p(i,1)=pptn*dcos(phirs)                                 ! px
 964        p(i,2)=pptn*dsin(phirs)                                 ! py
 965        p(i,3)=sigp*ppln                                        ! pz
 966        p(i,4)=dsqrt(p(i,1)**2+p(i,2)**2+p(i,3)**2+p(i,5)**2)   ! E
 967 * boost to system of hard parton
 968  222   call pyrobo(i,i,tetr,phir,0.d0,0.d0,0.d0)
 969
 970       return
 971       end
 972 ******************************* END PLJETR **************************
 973
 974 ******************************** PLSEAR ***************************
 975        SUBROUTINE PLSEAR (fmax,xmin)
 976 * find maximum and 'sufficient minimum' of jet production vertex distribution
 977 * xm, fm are outputs.
 978        IMPLICIT DOUBLE PRECISION(A-H, O-Z)
 979        external plthik
 980        common /plpar1/ tau0,T0,TC,sigqq,AW,RA,mvisc,nf
 981        save /plpar1/
 982        xmin=3.d0*RA
 983        fmax=0.d0
 984        do 10 j=1,1000
 985         x=xmin*(j-1)/999.d0
 986         f=plthik(x)
 987         if(f.gt.fmax) then
 988          fmax=f
 989         end if
 990   10   continue
 991        return
 992        end
 993 ****************************** END PLSEAR **************************
 994
 995 ******************************** RADSEAR ***************************
 996        SUBROUTINE RADSEAR (fmax,xmin)
 997 * find maximum and 'sufficient minimum' of radiative energy loss distribution
 998 * xm, fm are outputs.
 999       IMPLICIT DOUBLE PRECISION(A-H, O-Z)
1000       external plfun1
1001       common /radcal/ aa, bb
1002       save /radcal/
1003       xmin=bb-aa
1004       fmax=0.d0
1005       do j=1,1000
1006        x=aa+xmin*(j-1)/999.d0
1007        f=plfun1(x)
1008        if(f.gt.fmax) then
1009         fmax=f
1010        end if
1011       end do
1012       return
1013       end
1014 ****************************** END RADSEAR **************************
1015
1016 ********************************* BIPSEAR ***************************
1017       SUBROUTINE BIPSEAR (fmax,xmin)
1018 * find maximum and 'sufficient minimum' of jet production cross section
1019 * as a function of impact paramater (xm, fm are outputs)
1020       IMPLICIT DOUBLE PRECISION(A-H, O-Z)
1021       external funbip
1022       common /plpar1/ tau0,T0,TC,sigqq,AW,RA,mvisc,nf
1023       save /plpar1/
1024       xmin=3.d0*RA
1025       fmax=0.d0
1026       do j=1,1000
1027        x=xmin*(j-1)/999.d0
1028        f=funbip(x)
1029        if(f.gt.fmax) then
1030         fmax=f
1031        end if
1032       end do
1033       return
1034       end
1035 ****************************** END RADSEAR **************************
1036
1037 **************************** SIMPA **********************************
1038       SUBROUTINE SIMPA (A1,B1,H1,REPS1,AEPS1,FUNCT,X,
1039      1                     AI,AIH,AIABS)
1040 * calculate intergal of function FUNCT(X) on the interval from A1 to B1
1041       IMPLICIT DOUBLE PRECISION(A-H, O-Z)
1042       IMPLICIT INTEGER(I-N)
1043       DIMENSION F(7), P(5)
1044       H=dSIGN ( H1, B1-A1 )
1045       S=dSIGN (1.d0, H )
1046       A=A1
1047       B=B1
1048       AI=0.0d0
1049       AIH=0.0d0
1050       AIABS=0.0d0
1051       P(2)=4.d0
1052       P(4)=4.d0
1053       P(3)=2.d0
1054       P(5)=1.d0
1055       IF(B-A)1,2,1
1056  1    REPS=ABS(REPS1)
1057       AEPS=ABS(AEPS1)
1058       DO 3 K=1,7
1059  3    F(K)=10.d16
1060       X=A
1061       C=0.d0
1062       F(1)=FUNCT(X)/3.d0
1063  4    X0=X
1064       IF( (X0+4.d0*H-B)*S)5,5,6
1065  6    H=(B-X0)/4.d0
1066       IF ( H ) 7,2,7
1067  7    DO 8 K=2,7
1068  8    F(K)=10.d16
1069       C=1.d0
1070  5    DI2=F (1)
1071       DI3=ABS( F(1) )
1072       DO 9 K=2,5
1073       X=X+H
1074       IF((X-B)*S)23,24,24
1075  24   X=B
1076  23   IF(F(K)-10.d16)10,11,10
1077  11   F(K)=FUNCT(X)/3.d0
1078  10   DI2=DI2+P(K)*F(K)
1079  9    DI3=DI3+P(K)*ABS(F(K))
1080       DI1=(F(1)+4.*F(3)+F(5))*2.d0*H
1081       DI2=DI2*H
1082       DI3=DI3*H
1083       IF (REPS) 12,13,12
1084  13   IF (AEPS) 12,14,12
1085  12   EPS=ABS((AIABS+DI3)*REPS)
1086       IF(EPS-AEPS)15,16,16
1087  15   EPS=AEPS
1088  16   DELTA=ABS(DI2-DI1)
1089       IF(DELTA-EPS)20,21,21
1090  20   IF(DELTA-EPS/8.d0)17,14,14
1091  17   H=2.d0*H
1092       F(1)=F(5)
1093       F(2)=F(6)
1094       F(3)=F(7)
1095       DO 19 K=4,7
1096  19   F(K)=10.d16
1097       GO TO 18
1098  14   F(1)=F(5)
1099       F(3)=F(6)
1100       F(5)=F(7)
1101       F(2)=10.d16
1102       F(4)=10.d16
1103       F(6)=10.d16
1104       F(7)=10.d16
1105  18   DI1=DI2+(DI2-DI1)/15.d0
1106       AI=AI+DI1
1107       AIH=AIH+DI2
1108       AIABS=AIABS+DI3
1109       GO TO 22
1110  21   H=H/2.d0
1111       F(7)=F(5)
1112       F(6)=F(4)
1113       F(5)=F(3)
1114       F(3)=F(2)
1115       F(2)=10.d16
1116       F(4)=10.d16
1117       X=X0
1118       C=0.
1119       GO TO 5
1120  22   IF(C)2,4,2
1121  2    RETURN
1122       END
1123 ************************* END SIMPA *******************************
1124
1125 **************************** SIMPB **********************************
1126       SUBROUTINE SIMPB (A1,B1,H1,REPS1,AEPS1,FUNCT,X,
1127      1                     AI,AIH,AIABS)
1128 * calculate intergal of function FUNCT(X) on the interval from A1 to B1
1129       IMPLICIT DOUBLE PRECISION(A-H, O-Z)
1130       IMPLICIT INTEGER(I-N)
1131       DIMENSION F(7), P(5)
1132       H=dSIGN ( H1, B1-A1 )
1133       S=dSIGN (1.d0, H )
1134       A=A1
1135       B=B1
1136       AI=0.0d0
1137       AIH=0.0d0
1138       AIABS=0.0d0
1139       P(2)=4.d0
1140       P(4)=4.d0
1141       P(3)=2.d0
1142       P(5)=1.d0
1143       IF(B-A)1,2,1
1144  1    REPS=ABS(REPS1)
1145       AEPS=ABS(AEPS1)
1146       DO 3 K=1,7
1147  3    F(K)=10.d16
1148       X=A
1149       C=0.d0
1150       F(1)=FUNCT(X)/3.d0
1151  4    X0=X
1152       IF( (X0+4.d0*H-B)*S)5,5,6
1153  6    H=(B-X0)/4.d0
1154       IF ( H ) 7,2,7
1155  7    DO 8 K=2,7
1156  8    F(K)=10.d16
1157       C=1.d0
1158  5    DI2=F (1)
1159       DI3=ABS( F(1) )
1160       DO 9 K=2,5
1161       X=X+H
1162       IF((X-B)*S)23,24,24
1163  24   X=B
1164  23   IF(F(K)-10.d16)10,11,10
1165  11   F(K)=FUNCT(X)/3.d0
1166  10   DI2=DI2+P(K)*F(K)
1167  9    DI3=DI3+P(K)*ABS(F(K))
1168       DI1=(F(1)+4.*F(3)+F(5))*2.d0*H
1169       DI2=DI2*H
1170       DI3=DI3*H
1171       IF (REPS) 12,13,12
1172  13   IF (AEPS) 12,14,12
1173  12   EPS=ABS((AIABS+DI3)*REPS)
1174       IF(EPS-AEPS)15,16,16
1175  15   EPS=AEPS
1176  16   DELTA=ABS(DI2-DI1)
1177       IF(DELTA-EPS)20,21,21
1178  20   IF(DELTA-EPS/8.d0)17,14,14
1179  17   H=2.d0*H
1180       F(1)=F(5)
1181       F(2)=F(6)
1182       F(3)=F(7)
1183       DO 19 K=4,7
1184  19   F(K)=10.d16
1185       GO TO 18
1186  14   F(1)=F(5)
1187       F(3)=F(6)
1188       F(5)=F(7)
1189       F(2)=10.d16
1190       F(4)=10.d16
1191       F(6)=10.d16
1192       F(7)=10.d16
1193  18   DI1=DI2+(DI2-DI1)/15.d0
1194       AI=AI+DI1
1195       AIH=AIH+DI2
1196       AIABS=AIABS+DI3
1197       GO TO 22
1198  21   H=H/2.d0
1199       F(7)=F(5)
1200       F(6)=F(4)
1201       F(5)=F(3)
1202       F(3)=F(2)
1203       F(2)=10.d16
1204       F(4)=10.d16
1205       X=X0
1206       C=0.
1207       GO TO 5
1208  22   IF(C)2,4,2
1209  2    RETURN
1210       END
1211 ************************* END SIMPB *******************************
1212
1213 ************************* PARINV **********************************
1214       SUBROUTINE PARINV(X,A,F,N,R)
1215 * gives interpolation of function F(X) with  arrays A(N) and F(N)
1216       IMPLICIT DOUBLE PRECISION(A-H, O-Z)
1217       IMPLICIT INTEGER(I-N)
1218       DIMENSION A(N),F(N)
1219       IF(X.LT.A(1))GO TO 11
1220       IF(X.GT.A(N))GO TO 4
1221       K1=1
1222       K2=N
1223  2    K3=K2-K1
1224       IF(K3.LE.1)GO TO 6
1225       K3=K1+K3/2
1226       IF(A(K3)-X) 7,8,9
1227  7    K1=K3
1228       GOTO2
1229  9    K2=K3
1230       GOTO2
1231  8    P=F(K3)
1232       RETURN
1233  3    B1=A(K1)
1234       B2=A(K1+1)
1235       B3=A(K1+2)
1236       B4=F(K1)
1237       B5=F(K1+1)
1238       B6=F(K1+2)
1239       R=B4*((X-B2)*(X-B3))/((B1-B2)*(B1-B3))+B5*((X-B1)*(X-B3))/
1240      1 ((B2-B1)*(B2-B3))+B6*((X-B1)*(X-B2))/((B3-B1)*(B3-B2))
1241       RETURN
1242  6    IF(K2.NE.N)GO TO 3
1243       K1=N-2
1244       GOTO3
1245  4    C=ABS(X-A(N))
1246       IF(C.LT.0.1d-7) GO TO 5
1247       K1=N-2
1248  13   CONTINUE
1249 C13   PRINT 41,X
1250 C41   FORMAT(25H X IS OUT OF THE INTERVAL,3H X=,F15.9)
1251       GO TO 3
1252  5    R=F(N)
1253       RETURN
1254  11   C=ABS(X-A(1))
1255       IF(C.LT.0.1d-7) GO TO 12
1256       K1=1
1257       GOTO 13
1258  12   R=F(1)
1259       RETURN
1260       END
1261 C************************** END PARINV *************************************
1262
1263 * quark-quark scattering differential cross section
1264        double precision FUNCTION PLSIGH(Z)
1265        IMPLICIT DOUBLE PRECISION(A-H, O-Z)
1266        common /plpar1/ tau0,T0,TC,sigqq,AW,RA,mvisc,nf
1267        save /plpar1/
1268        pi=3.14159d0
1269        beta=(33.d0-2.d0*nf)/(12.d0*pi)
1270        alfs=1.d0/(beta*dlog(max(1.d-10,z/(TC*TC))))
1271        PLSIGH=8.d0*pi*alfs*alfs/(9.d0*z*z)
1272        return
1273        end
1274
1275 * differential radiated gluon spectrum in BDMS model
1276       double precision FUNCTION PLFUN1(or)
1277       IMPLICIT DOUBLE PRECISION(A-H, O-Z)
1278       common /plpar1/ tau0,T0,TC,sigqq,AW,RA,mvisc,nf
1279       common /pljdat/ ej, z, ygl, alfs, um, epa
1280       common /pleave/ taul, temlev
1281       common /factor/ cfac, kf
1282       save /plpar1/,/pljdat/,/pleave/,/factor/
1283       pi=3.14159d0
1284       x=min((1.d0-ygl*z/or),or/ej)
1285       if(x.le.0.d0) x=0.d0
1286       if(x.ge.1.d0) x=0.9999d0
1287       if(kf.eq.21) then
1288        if(x.ge.0.5d0) x=1.d0-x
1289        spinf=0.5d0*(1.+(1.d0-x)**4+x**4)/(1.d0-x)
1290       else
1291        spinf=1.d0-x+0.5d0*x*x
1292       end if
1293       ak=ygl*z/(or*(1.d0-x))
1294       al=taul*5.06d0
1295       uu=0.5d0*al*dsqrt(abs(0.5d0*(1.d0-x+cfac*x*x)*ak*
1296      >   dlog(max(16.d0/ak,1.d-10))))/ygl
1297 * if  quark production outside the QGP then
1298 * arg=(((dsin(uu)*cosh(uu))**2)+((dcos(uu)*sinh(uu))**2))/(2.d0*uu*uu);
1299 * here quark production inside the QGP
1300 c AM: better numerical stability
1301
1302       if (uu .lt. 5.) then
1303          arg     = ((dcos(uu)*cosh(uu))**2)+((dsin(uu)*sinh(uu))**2)
1304          xlogarg = dlog(max(arg,1.d-20))
1305       else
1306          xlogarg = log(0.25) + 2. * uu
1307       endif
1308
1309
1310       gl1=(ygl/(cfac*abs(z)))**0.3333333d0
1311       gl2=(um/epa)**1.333333d0
1312       dc=1.d0/((1.d0+((gl1*gl2*or)**1.5d0))**2)        ! massive parton
1313 c      dc=1.d0                                         !massless parton
1314 c      plfun1 = dc*3.d0*alfs*ygl*dlog(max(arg,1.d-20))*spinf/(pi*al*or)
1315       plfun1   = dc*3.d0*alfs*ygl*xlogarg*spinf/(pi*al*or)
1316       return
1317       end
1318
1319 * angular distribution of emitted gluons
1320       double precision function gluang(x)
1321       IMPLICIT DOUBLE PRECISION(A-H, O-Z)
1322       s=0.0863d0
1323       gluang=x*dexp(-1.d0*(x-s)*(x-s)/(2.d0*s*s))
1324       return
1325       end
1326
1327 * temperature-dependence of parton-plasma integral cross section
1328        double precision FUNCTION PLS(X)
1329        IMPLICIT DOUBLE PRECISION(A-H, O-Z)
1330        external plsigh
1331        common /plpar1/ tau0,T0,TC,sigqq,AW,RA,mvisc,nf
1332        common /plpar2/ pln0,taupl,tauh,sigpl,sigh,sigplh,sigqqh,rg,rgn
1333        common /plen/ epartc, um
1334        save /plpar1/,/plpar2/,/plen/
1335        t=X
1336        pi=3.14159d0
1337        bub=4.d0*t/TC
1338        alf=6.d0*pi/((33.d0-2.d0*nf)*dlog(max(bub,1.d-10)))
1339        ZZ0=4.d0*t*t*pi*alf*(1.d0+nf/6.d0)
1340        scm=4.d0*t*epartc+um*um+4.d0*t*t
1341        ZZ1=max((scm-((um+2.d0*t)**2))*(scm-((um-2.d0*t)**2))/scm,ZZ0)
1342        HH1=0.01d0*ZZ1
1343        REPS=0.01d0
1344        AEPS=1.d-8
1345        CALL SIMPA(ZZ0,ZZ1,HH1,REPS,AEPS,plsigh,ZZ,RESS,AIH,AIABS)
1346        PLS=0.39d0*2.25d0*2.25d0*RESS*(16.d0+4.d0*nf)/(16.d0+9.d0*nf)
1347        return
1348        end
1349
1350 * temperature-dependence of QGP viscosity (if mvisc=1,2)
1351       double precision FUNCTION PLVISC(X)
1352       IMPLICIT DOUBLE PRECISION(A-H, O-Z)
1353       common /plpar1/ tau0,T0,TC,sigqq,AW,RA,mvisc,nf
1354       save /plpar1/
1355       pi=3.14159d0
1356       T=X
1357       TC1=5.06d0*TC
1358       if(X.le.TC1) T=TC1
1359       if(mvisc.eq.0) then
1360        c=0.d0
1361       elseif(mvisc.eq.1) then
1362        a=3.4d0*(1.d0+0.12d0*(2.d0*nf+1.d0))
1363        b=15.d0*(1.d0+0.06d0*nf)
1364        c=4.d0*pi*pi*(10.5d0*nf/a+16.d0/b)/675.d0
1365       else
1366        c=(1.7d0*nf+1.d0)*0.342d0/(1.d0+nf/6.d0)
1367       end if
1368       bub=4.d0*T/TC1
1369       alf=6.d0*pi/((33.d0-2.d0*nf)*dlog(max(bub,1.d-10)))
1370       alf1=1.d0/alf
1371       PLVISC=c*(T**3)/(alf*alf*dlog(max(1.d-10,alf1)))
1372       return
1373       end
1374
1375 * space-time dependence of QGP number density
1376        double precision FUNCTION PLN(X,r1,r2,y)
1377        IMPLICIT DOUBLE PRECISION(A-H, O-Z)
1378        external pythik
1379        common /plpar1/ tau0,T0,TC,sigqq,AW,RA,mvisc,nf
1380        common /plpar2/ pln0,taupl,tauh,sigpl,sigh,sigplh,sigqqh,rg,rgn
1381        common /plevol/ taup(5000),temp(5000),denp(5000),enep(5000)
1382        common /pythic/ PBAB(110),PTAB(110),PTAAB(110)
1383        save /plpar1/,/plpar2/,/plevol/,/pythic/
1384        pi=3.14159d0
1385        t=X
1386        if(t.lt.taupl) then
1387         call parinv(t,taup,denp,5000,res)
1388        else
1389         res=1.2d0*(16.d0+9.d0*nf)*((5.06d0*TC)**3)/(pi*pi)
1390        end if
1391        res=res*(pythik(r1)*pythik(r2)*pi*RA*RA/PTAAB(1))**0.75d0
1392        res=res*dexp(-1.d0*y*y/24.5d0)
1393        PLN=max(1.d-8,res)
1394        return
1395        end
1396
1397 * space-time dependence of QGP temperature
1398        double precision FUNCTION PLT(X,r1,r2,y)
1399        IMPLICIT DOUBLE PRECISION(A-H, O-Z)
1400        common /plpar1/ tau0,T0,TC,sigqq,AW,RA,mvisc,nf
1401        common /plpar2/ pln0,taupl,tauh,sigpl,sigh,sigplh,sigqqh,rg,rgn
1402        common /plevol/ taup(5000),temp(5000),denp(5000),enep(5000)
1403        common /pythic/ PBAB(110),PTAB(110),PTAAB(110)
1404        save /plpar1/,/plpar2/,/plevol/,/pythic/
1405        pi=3.14159d0
1406        t=X
1407        if(t.lt.taupl) then
1408         call parinv(t,taup,temp,5000,res)
1409        else
1410         res=TC
1411        end if
1412         res=res*(pythik(r1)*pythik(r2)*pi*RA*RA/PTAAB(1))**0.25d0
1413         res=res*(dexp(-1.d0*y*y/24.5d0))**0.333333d0
1414         PLT=max(1.d-8,res)
1415        return
1416        end
1417
1418 * impact parameter dependence of jet production cross section
1419       double precision function funbip(x)
1420       IMPLICIT DOUBLE PRECISION(A-H, O-Z)
1421       external ftaa
1422       common /pyint7/ sigt(0:6,0:6,0:5)
1423       save /pyint7/
1424       br=x
1425       sigin=sigt(0,0,0)-sigt(0,0,1)
1426       taa=ftaa(br)
1427       funbip=taa*br*(1.d0-dexp(-0.1d0*taa*sigin))
1428       return
1429       end
1430
1431 * distribution over jet production vertex position
1432       double precision FUNCTION plthik(X)
1433       IMPLICIT DOUBLE PRECISION(A-H, O-Z)
1434       external pythik
1435       common /parimp/ b1, psib1, rb1, rb2, noquen
1436       save /parimp/
1437       bu=X
1438       r12=dsqrt(abs(bu*bu+b1*b1/4.d0+bu*b1*dcos(psib1)))
1439       r22=dsqrt(abs(bu*bu+b1*b1/4.d0-bu*b1*dcos(psib1)))
1440       PLTHIK=bu*pythik(r12)*pythik(r22)
1441       return
1442       end
1443
1444 * nuclear overlap function at impact parameter b
1445       double precision function ftaa(r)
1446       IMPLICIT DOUBLE PRECISION(A-H, O-Z)
1447       common /pythic/ PBAB(110),PTAB(110),PTAAB(110)
1448       save /pythic/
1449       call parinv(r,PBAB,PTAAB,110,RES)
1450       ftaa=RES
1451       return
1452       end
1453 *
1454       double precision function PFUNC1(x)
1455       IMPLICIT DOUBLE PRECISION(A-H, O-Z)
1456       external PFUNC2
1457       common /pynup1/ bp,xx
1458       common /plpar1/ tau0,T0,TC,sigqq,AW,RA,mvisc,nf
1459       save /plpar1/
1460       xx=x
1461       EPS=0.05d0
1462       A=0.d0
1463       B=3.d0*RA
1464       H=0.01d0*(B-A)
1465       CALL SIMPB(A,B,H,EPS,1.d-8,PFUNC2,Y,RES,AIH,AIABS)
1466       PFUNC1=RES
1467       return
1468       end
1469 *
1470       double precision function PFUNC2(y)
1471       IMPLICIT DOUBLE PRECISION(A-H, O-Z)
1472       external pythik
1473       common /pynup1/ bp,x
1474       r1=sqrt(abs(y*y+bp*bp/4.+y*bp*cos(x)))
1475       r2=sqrt(abs(y*y+bp*bp/4.-y*bp*cos(x)))
1476       PFUNC2=y*pythik(r1)*pythik(r2)
1477       return
1478       end
1479
1480 * nuclear thickness function
1481       double precision function pythik(r)
1482       IMPLICIT DOUBLE PRECISION(A-H, O-Z)
1483       common /pythic/ PBAB(110),PTAB(110),PTAAB(110)
1484       save /pythic/
1485       call parinv(r,PBAB,PTAB,110,RES)
1486       pythik=RES
1487       return
1488       end
1489
1490 * Wood-Saxon nucleon distrubution
1491       double precision function prhoaa(z)
1492       IMPLICIT DOUBLE PRECISION(A-H, O-Z)
1493       common /plpar1/ tau0,T0,TC,sigqq,AW,RA,mvisc,nf
1494       common /pygeom/ BC
1495       save /plpar1/,/pygeom/
1496       pi=3.14159d0
1497       df=0.54d0
1498       r=sqrt(bc*bc+z*z)
1499       rho0=3.d0/(4.d0*pi*RA**3)/(1.d0+(pi*df/RA)**2)
1500       prhoaa=rho0/(1.d0+exp((r-RA)/df))
1501       return
1502       end
1503
1504 * function to generate gauss distribution
1505       double precision function gauss(x0,sig)
1506       IMPLICIT DOUBLE PRECISION(A-H, O-Z)
1507  41   u1=pyr(0)
1508       u2=pyr(0)
1509       v1=2.d0*u1-1.d0
1510       v2=2.d0*u2-1.d0
1511       s=v1**2+v2**2
1512       if(s.gt.1) go to 41
1513       gauss=v1*dsqrt(-2.d0*dlog(s)/s)*sig+x0
1514       return
1515       end
1516 **************************************************************************