New PID class for EMCAL, bayesian analysis done with ESD data, PID information filled...
[u/mrichter/AliRoot.git] / PYTHIA6 / pyquen.F
1 *----------------------------------------------------------------------
2 *
3 *  Filename             : PYQUEN.F
4 *
5 *  First version created: 20-AUG-1997   Author : Igor Lokhtin  
6 *  Last revision        : 23-SEP-2004 
7 *
8 *======================================================================
9 *
10 *  Description : Event generator for simulation of parton rescattering 
11 *                and energy loss in quark-gluon plasma created in heavy 
12 *                ion AA collisons at LHC 
13 *               (implemented as modification of standard pythia jet event) 
14 *
15 *  Method     : I.P.Lokhtin, A.M.Snigirev, Eur.Phys.J. C16 (2000) 527-536;   
16 *               I.P.Lokhtin, A.M.Snigirev, e-print hep-ph/0406038.  
17 *                   
18 *
19 *======================================================================
20
21       SUBROUTINE PYQUEN(A,ifb,bfix)
22       IMPLICIT DOUBLE PRECISION(A-H, O-Z)
23       IMPLICIT INTEGER(I-N)
24       INTEGER PYK,PYCHGE,PYCOMP
25       external pydata  
26       external pyp,pyr,pyk,pyjoin,pyshow
27       external ftaa,funbip 
28       common /pyjets/ n,npad,k(4000,5),p(4000,5),v(4000,5)
29       common /pydat1/ mstu(200),paru(200),mstj(200),parj(200)       
30       common /pysubs/ msel,mselpd,msub(500),kfin(2,-40:40),ckin(200)      
31       common /plpar1/ tau0,T0,TC,sigqq,AW,RA,mvisc,nf
32       common /plglur/ glur(1000,4),kglu(1000,6),nrg,nrgm 
33       common /plquar/ pqua(1000,5),kqua(1000,5),nrq 
34       common /parimp/ b1,psib1,rb1,rb2 
35       common /bintaa/ br 
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) 
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)
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            
162
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 
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  
240       do i=1,500  
241        nis(i)=0 
242        nas(i)=0 
243        nss(i)=0 
244        nus(i)=0 
245       end do                      
246       do i=9,n 
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 
256        elseif(ks.ne.2.and.ksp.ne.2.and.ns.gt.0) then 
257         i1=i1+1               ! last i1 lines not included in strings 
258        end if 
259       end do 
260       i0=n-nes-i1             ! first i0 lines not included in strings 
261       do i=1,ns 
262        nss(i)=nss(i)+i0 
263       end do  
264       
265 * move fragmented particles in bottom of event list  
266       i=i0+1      
267  2    ks=k(i,1)
268       ksp=k(i-1,1) 
269       if(ks.ne.2.and.ksp.ne.2) then 
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
288       else  
289        i=i+1  
290       end if 
291       if(i.le.n-i1) goto 2  
292
293 * define number of additional entries in strings, nas(ns)                     
294       do i=1,nrg 
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           
304       end do 
305       do j=1,ns   
306        do i=1,j   
307         nus(j)=nus(j)+nas(i) 
308        end do 
309       end do 
310             
311 * add emitted gluons in event list  
312       nu=n 
313       n=n+nrg 
314       do i=nu,nu-i1,-1 
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 
321       end do 
322       do ia=ns-1,1,-1  
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 
331       end do 
332
333       do i=1,nrg 
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   
356       end do  
357       
358 * rearrange partons to form strings in event list 
359       do ij=1,1000 
360        ijoin(ij)=0 
361       end do 
362       do i=1,ns 
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)
374       end do  
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
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
487       external plthik, pln, plt, pls, pygauss, gluang 
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)
493       common /plglur/ glur(1000,4),kglu(1000,6),nrg,nrgm  
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  
530        do j=1,4
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 
649           glur(nrg,3)=datan(dpgl)                                    ! phi
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 
749  221   ep0=-1.*tauu*(dlog(max(1.d-10,pyr(0)))+dlog(max(1.d-10,
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
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 **************************************************************************