PYTHIA6/pyquen.F

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