- Update to PYQUEN v1.3 (http://lokhtin.web.cern.ch/lokhtin/pyquen/pyquen.txt)
[u/mrichter/AliRoot.git] / PYTHIA6 / pyquen1_3.F
1 *----------------------------------------------------------------------
2 *
3 *  Filename             : PYQUEN.F
4 *
5 *  Author               : Igor Lokhtin  (Igor.Lokhtin@cern.ch)
6 *  Version              : PYQUEN1_3.f 
7 *  Last revision        : 03-OCT-2007 
8 *
9 *======================================================================
10 *
11 *  Description : Event generator for simulation of parton rescattering 
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) 
15 *
16 *  Reference: I.P. Lokhtin, A.M. Snigirev, Eur. Phys. J. C 46 (2006) 211   
17 *                   
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  
25       external pyp,pyr,pyk,pyjoin,pyshow
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
31       common /plglur/ glur(1000,4),kglu(1000,6),nrg,nrgm 
32       common /plquar/ pqua(1000,5),kqua(1000,5),nrq 
33       common /parimp/ b1,psib1,rb1,rb2,noquen 
34       common /pyqpar/ T0u,tau0u,nfu,ienglu,ianglu 
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       mvisc=0                              ! flag of QGP viscosity (off here) 
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 *    
56       pi=3.14159d0 
57
58 * avoid stopping run if pythia does not conserve energy due to collisional loss 
59       mstu(21)=1 
60
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))   
94       T00=T000*((rtaa*sb0)/(rtaa0*sb))**0.25d0 
95       T0=T00*(AW/207.d0)**0.166667d0 
96       
97 * generate single event with partonic energy loss 
98       nrg=0 
99       ehard=ckin(3) 
100       if(b1.le.1.85d0*RA) then
101        call plinit(ehard)  
102        call plevnt(ehard)
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  
122       if(mstj(41).ne.0) goto 5
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 
129         qmax1=pyp(i,10)          ! transverse momentum of first hard parton
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 
134         qmax2=pyp(i,10)          ! transverse momentum of second hard parton 
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   
144       mstj(41)=0 
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            
171
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 
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 *
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  
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 *
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 
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  
261       do i=1,500  
262        nis(i)=0 
263        nas(i)=0 
264        nss(i)=0 
265        nus(i)=0 
266       end do                      
267       do i=9,n 
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
274         nes=nes+nis(ns+1)     ! nes - total number of entries  
275         nss(ns+1)=nes 
276         ns=ns+1 
277        elseif(ks.ne.2.and.ksp.ne.2.and.ns.gt.0) then 
278         i1=i1+1               ! last i1 lines not included in strings 
279        end if 
280       end do 
281       i0=n-nes-i1             ! first i0 lines not included in strings 
282       do i=1,ns 
283        nss(i)=nss(i)+i0 
284       end do  
285       
286 * move fragmented particles in bottom of event list  
287       i=i0+1      
288  2    ks=k(i,1)
289       ksp=k(i-1,1) 
290       if(ks.ne.2.and.ksp.ne.2) then 
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 
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 *
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 
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 *
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
321       else  
322        i=i+1  
323       end if 
324       if(i.le.n-i1) goto 2  
325
326 * define number of additional entries in strings, nas(ns)                     
327       do i=1,nrg 
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           
337       end do 
338       do j=1,ns   
339        do i=1,j   
340         nus(j)=nus(j)+nas(i) 
341        end do 
342       end do 
343             
344 * add emitted gluons in event list  
345       nu=n 
346       n=n+nrg 
347       do i=nu,nu-i1,-1 
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 
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 *
360       end do 
361       do ia=ns-1,1,-1  
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
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 *       
375        end do 
376       end do 
377
378       do i=1,nrg 
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))
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)      
400        p(ia,5)=0.d0   
401       end do  
402       
403 * rearrange partons to form strings in event list 
404       do ij=1,1000 
405        ijoin(ij)=0 
406       end do 
407       do i=1,ns 
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)
419       end do  
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
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 
472       common /plevol/ taup(5000),temp(5000),denp(5000),enep(5000)  
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   
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  
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
533       external plthik, pln, plt, pls, gauss, gluang 
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)
539       common /pyqpar/ T0u,tau0u,nfu,ienglu,ianglu 
540       common /plglur/ glur(1000,4),kglu(1000,6),nrg,nrgm  
541       common /factor/ cfac, kf 
542       common /pleave/ taul, temlev   
543       common /parimp/ b1, psib1, rb1, rb2, noquen 
544       common /plen/ epartc, um 
545       common /plos/ elr,rsk 
546       common /numje1/ nuj1, nuj2  
547       save /pyjets/, /plglur/
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) 
569       if(noquen.eq.1) goto 7
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  
578        do j=1,4
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 
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   
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) 
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 
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 
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 
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 
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         
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 
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                                                  
697           pgl=pi*(2.d0*pyr(0)-1.d0) 
698 * in comoving system 
699           pxgl=abs(elr)*stetr*(dcos(phir)*dcos(tgl)-
700      >      dsin(phir)*dsin(tgl)*dsin(pgl)) 
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))           
705           psgl=dsqrt(abs(ptgl*ptgl+pzgl*pzgl)) 
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 *
711           dpgl=pygl/pxgl        
712           glur1=abs(elr)                                       ! energy 
713           glur3=datan(dpgl)                                    ! phi
714           if(pxgl.lt.0.d0) then 
715            if(pygl.ge.0.d0) then 
716             glur3=glur3+pi 
717            else 
718             glur3=glur3-pi  
719            end if 
720           end if   
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
732 * set gluon codes 
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  
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 
770  100    call pyrobo(ip,ip,tetr,phir,0.d0,0.d0,0.d0)
771        end if 
772        end if 
773        end if 
774  2    continue 
775  7    continue
776  
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)
792       common /pyqpar/ T0u,tau0u,nfu,ienglu,ianglu 
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  
809       tetr=pyp(i,13)                    ! parton theta   
810       phir1=-1.d0*phir 
811       tetr1=-1.d0*tetr 
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)  
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 
824  221   ep0=-1.*tauu*(dlog(max(1.d-10,pyr(0)))+dlog(max(1.d-10,
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
847 * calculate radiative energy loss per given scattering with subroutine plfun1 
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) 
870        if(ienglu.eq.2) elr=0.d0
871 * to chancel collisional energy loss (optional case) 
872        if(ienglu.eq.1) rsk=0.d0 
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)
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  
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))   
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 
900 * boost to system of hard parton 
901  222   call pyrobo(i,i,tetr,phir,0.d0,0.d0,0.d0)
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
915        common /parimp/ b1, psib1, rb1, rb2, noquen 
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
1264        common /parimp/ b1, psib1, rb1, rb2, noquen 
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
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       
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 **************************************************************************