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