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