PYTHIA6/QPYTHIA/pyquen1_5.F

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