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