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