Coding violation fixes
[u/mrichter/AliRoot.git] / THydjet / hydjet1_1 / hydjet1_1.f
1 *----------------------------------------------------------------------
2 *
3 *  Filename             : HYDJET1_1.F
4 *
5 *  Author               : Igor Lokhtin  
6 *  Version              : HYDJET1_1.f
7 *  Last revision        : 27-MAR-2006 
8 *
9 *======================================================================
10 *
11 *  Description : Fast event generator for AA collisons at the LHC            
12 *
13 *  Method I.P. Lokhtin and A.M. Snigirev, Eur. Phys. J C 45 (2006) 211        
14 *
15 *======================================================================
16
17       SUBROUTINE HYDRO(A,ifb,bmin,bmax,bfix,nh)  
18       real hsin,hgauss,hftaa 
19       real AW 
20       real A,bmin,bmax,bfix 
21       integer numjet,numpar 
22       integer ifb,nh,np 
23       external hsin,hgauss,hftaa,numjet,numpar,hyhard,hipsear 
24       external ludata 
25       common /lujets/ n,k(150000,5),p(150000,5),v(150000,5)
26       common /hyjets/ nl,kl(150000,5),pl(150000,5),vl(150000,5) 
27       common /hyipar/ bminh,bmaxh,AW,RA,sigin,np  
28       common /hyfpar/ bgen,nbcol,npart,npyt,nhyd
29       common /hyflow/ ytfl,ylfl,fpart
30       common /hyjpar/ nhsel,ptmin,njet 
31       save/lujets/,/hyjets/,/hyipar/,/hyfpar/,/hyflow/,/hyjpar/ 
32
33       integer nnhyd, khyd
34       real phyd, vhyd
35       common /hyd/ nnhyd, khyd(150000,5),phyd(150000,5),vhyd(150000,5)
36       save /hyd/
37
38 * reset lujets and hyd arrays before event generation 
39       
40       n=0 
41       nnhyd=0 
42       do ncl=1,150000
43        do j=1,5
44         p(ncl,j)=0.
45         phyd(ncl,j)=0.
46         v(ncl,j)=0. 
47         vhyd(ncl,j)=0.  
48         k(ncl,j)=0
49         khyd(ncl,j)=0 
50        enddo
51       end do 
52
53 * set initial beam paramters: atomic weigth and radius in fm 
54       AW=A         
55       RA=1.15*AW**0.333333  
56
57       pi=3.14159
58       
59 * generate impact parameter of A-A collision 
60
61       if(ifb.eq.0) then 
62        if(bfix.lt.0.) then    
63         write(6,*) 'Impact parameter less than zero!'  
64         bfix=0. 
65        end if  
66        if (bfix.gt.2.) then 
67         write(6,*) 'Impact parameter larger than two nuclear radius!'  
68         bfix=2.        
69        end if 
70        b1=bfix*RA
71        bgen=bfix 
72       else        
73        if(bmin.lt.0.) then    
74         write(6,*) 'Impact parameter less than zero!'  
75         bmin=0. 
76        end if 
77        if(bmax.gt.2.) then    
78         write(6,*) 'Impact parameter larger than two nuclear radius!'
79         bmax=2.  
80        end if             
81        bminh=bmin 
82        bmaxh=bmax   
83        call hipsear(fmax1,xmin1) 
84        fmax=fmax1 
85        xmin=xmin1 
86  3     bb1=xmin*rlu(0)+bminh  
87        ff1=fmax*rlu(0) 
88        fb=hsin(bb1) 
89        if(ff1.gt.fb) goto 3    
90        b1=bb1*RA       
91        bgen=bb1                                                                        
92       end if 
93
94 * set flow parameters 
95       Tf=0.14                             ! freeze-out temperature     
96       if (ylfl.lt.0.01.or.ylfl.gt.7.) ylfl=5.
97       etmax=ylfl                          ! longitudinal flow rapidity 
98       if (ytfl.lt.0.01.or.ytfl.gt.3.) ytfl=1.
99       ytmax=ytfl                          ! transverse flow rapidity
100
101 * set inelastic NN cross section, mb 
102       sigin=58.  
103             
104 * calculate number of nucelons-participants 
105       bb=bgen*RA                              ! impact parameter 
106       npart=numpar(bb)                        ! Npart(b) 
107       npar0=411                               ! Npart(Pb,b=0) 
108
109 * calculate number of binary NN sub-collisions 
110       br=max(1.e-10,0.25*bgen*bgen)
111       factor=9.*sigin*AW*AW/(80.*pi*RA*RA)      
112       nbcol=int(factor*hftaa(br))              ! Nsub(b)
113       nbco0=1923                               ! Nsub(Pb,b=0) 
114   
115 * generate total multiplicity in event, np,  
116 * fpart - fraction of soft multiplicity proportional to # of participants, 
117 * fbcol=1-fpart - fraction of multiplicity proprtional to # of NN subcollisions 
118       if(fpart.le.0.or.fpart.gt.1.) fpart=1. 
119       fbcol=1.-fpart 
120       rnp=nh*(fpart*npart+fbcol*nbcol)/(fpart*npar0+fbcol*nbco0)
121       np=int(rnp) 
122       sign=sqrt(rnp)
123  1    if(nhsel.lt.4.and.np.gt.0) np=max(0,int(hgauss(rnp,sign))) 
124       if(np.gt.150000) then 
125        write(6,*) 'Warning, soft multiplicity too large!'   
126        goto 1 
127       end if 
128       
129 * generate hard parton-parton scattering (Q>ptmin) 'njet' times with PYTHIA       
130       if(nhsel.ne.1.and.nhsel.ne.2.and.nhsel.ne.3.and.nhsel.ne.4) 
131      > nhsel=0 
132       njet=0 
133       if(nhsel.ne.0) then 
134        if(ptmin.lt.5.or.ptmin.gt.500.) ptmin=10. 
135        q=ptmin 
136        njet=numjet(q) 
137        call hyhard 
138       end if 
139
140       npyt=n 
141       nhyd=np 
142 c      if(nhsel.lt.3) then 
143 c       nhyd=max(0,np-npyt)
144 c      else 
145 c       nhyd=0 
146 c       np=n
147 c      end if 
148       if(nhyd.eq.0) goto 4 
149
150 * generate sort of hadrons (pions, kaons, nucleons) in jetset7* format 
151       do ip=npyt+1,npyt+np                       !cycle on particles 
152        yy=49.*rlu(0)  
153        if(yy.lt.11.83333333) then
154         kf=211
155        elseif(yy.lt.23.66666667) then
156         kf=-211
157        elseif(yy.lt.35.5) then
158         kf=111
159        elseif(yy.lt.38.375) then
160         kf=321
161        elseif(yy.lt.41.25) then
162         kf=-321 
163        elseif(yy.lt.44.125) then
164         kf=310
165        elseif(yy.lt.47.) then
166         kf=130  
167        elseif(yy.lt.47.5) then
168         kf=2212  
169        elseif(yy.lt.48.) then
170         kf=-2212
171        elseif(yy.lt.48.5) then
172         kf=2112 
173        else  
174         kf=-2112
175        end if
176        n=n+1
177        k(n,1)=1
178        k(n,2)=kf
179        do j=3,5
180         k(n,j)=0
181        end do
182
183        do j=1,5
184         v(n,j)=0.
185        enddo
186         kfa=iabs(kf)
187         p(n,5)=ulmass(kfa) 
188
189 * generate uniform distribution in system of a fluid element
190  2     ep0=-1.*Tf*(log(max(1.e-10,rlu(0)))+log(max(1.e-10,rlu(0)))
191      >   +log(max(1.e-10,rlu(0)))) 
192        if(ep0.le.p(n,5)) go to 2 
193        pp0=sqrt(abs(ep0**2-abs(p(n,5)**2))) 
194        probt=pp0/ep0 
195        if(rlu(0).gt.probt) go to 2 
196        ctp0=2.*rlu(0)-1. 
197        stp0=sqrt(abs(1.-ctp0**2)) 
198        php0=2.*pi*rlu(0) 
199
200 * generate coordinates of a fluid element
201 c       etaf=etmax*(2.*rlu(0)-1.)         ! flat initial eta-spectrum
202        etaf=hgauss(0.,etmax)              ! gaussian initial eta-spectrum
203        phif=2.*pi*rlu(0) 
204        rm1=sqrt(abs(RA*RA-b1*b1/4.*(sin(phif)**2)))+b1*cos(phif)/2. 
205        rm2=sqrt(abs(RA*RA-b1*b1/4.*(sin(phif)**2)))-b1*cos(phif)/2. 
206        RF=min(rm1,rm2) 
207        rrf=RF*RF*rlu(0) 
208
209 * generate four-velocity of a fluid element
210        vradf=sinh(ytmax)
211        sb=RA*RA*(pi-2.*asin(b1/RA/2.))-b1*sqrt(abs(RA*RA-b1*b1/4.)) 
212        reff=sqrt(sqrt(sb/pi)*RA) 
213        urf=vradf*sqrt(abs(rrf))/reff          ! linear transverse profile
214 c       urf=vradf*rrf/reff**2                  !square transverse profile 
215        utf=cosh(etaf)*sqrt(abs(1.+urf*urf)) 
216        uzf=sinh(etaf)*sqrt(abs(1.+urf*urf)) 
217
218 * boost in the laboratory system
219        uipi=pp0*(urf*stp0*cos(phif-php0)+uzf*ctp0) 
220        bfac=uipi/(utf+1.)+ep0 
221        p(n,1)=pp0*stp0*sin(php0)+urf*sin(phif)*bfac 
222        p(n,2)=pp0*stp0*cos(php0)+urf*cos(phif)*bfac                    
223        p(n,3)=pp0*ctp0+uzf*bfac
224        p(n,4)=sqrt(p(n,1)**2+p(n,2)**2+p(n,3)**2+p(n,5)**2) 
225             
226       end do 
227  
228  4    continue
229
230 *      write(*,*) 'NHYD, NPYT, NTOT',nhyd,npyt,nhyd+npyt
231
232 * fill array 'hyd' 
233  
234       nnhyd = nhyd+npyt 
235       
236       do ih=1,n
237        do jh=1,5
238         phyd(ih,jh)=p(ih,jh)
239         khyd(ih,jh)=k(ih,jh)                  
240         vhyd(ih,jh)=v(ih,jh) 
241        end do
242       end do
243
244       return
245       end
246
247 ********************************* HYHARD ***************************
248       SUBROUTINE HYHARD 
249 *     generate 'njet' number of hard parton-parton scatterings with PYTHIA 
250       IMPLICIT DOUBLE PRECISION(A-H, O-Z)
251       IMPLICIT INTEGER(I-N)
252       REAL ptmin,pj,vj,pl,vl,bminh,bmaxh,AW,RA,sigin,bgen 
253       INTEGER PYK,PYCHGE,PYCOMP
254       external pydata 
255       external pyp,pyr,pyk,pyquen  
256       common /pyjets/ n,npad,k(4000,5),p(4000,5),v(4000,5)
257       common /lujets/ nj,kj(150000,5),pj(150000,5),vj(150000,5) 
258       common /hyjets/ nl,kl(150000,5),pl(150000,5),vl(150000,5)            
259       COMMON /PYDAT1/MSTU(200),PARU(200),MSTJ(200),PARJ(200)
260       COMMON /PYDAT2/KCHG(500,4),PMAS(500,4),PARF(2000),VCKM(4,4)
261       COMMON /PYDAT3/MDCY(500,3),MDME(8000,2),BRAT(8000),KFDP(8000,5)
262       COMMON /PYSUBS/MSEL,MSELPD,MSUB(500),KFIN(2,-40:40),CKIN(200)
263       common /pypars/ mstp(200),parp(200),msti(200),pari(200)
264       common /hyjpar/ nhsel,ptmin,njet      
265       common /hyipar/ bminh,bmaxh,AW,RA,sigin,np 
266       common /hyfpar/ bgen,nbcol,npart,npyt,nhyd 
267       save /pyjets/,/lujets/,/hyjets/,/pydat1/,/pydat2/,/pydat3/,
268      +    /pysubs/,/hyjpar/,/hyipar/,/hyfpar/ 
269       
270 * reset array of partons in 'hyjets'
271       nl=0 
272       do i=1,150000
273        do j=1,5
274         pl(i,j)=0.
275        end do
276        do j=1,5
277         vl(i,j)=0.
278        end do
279        do j=1,5
280         kl(i,j)=0
281        end do
282       end do 
283
284 * generate 'njet' PYTHIA events and fill arrays for partons and hadrons 
285       if(njet.ge.1) then 
286        mdcy(pycomp(111),1)=0               ! no pi0 decay 
287        mdcy(pycomp(310),1)=0               ! no K_S0 decay 
288        do ihard=1,njet       
289         mstp(111)=0 
290 c        mstj(41)=0                         ! vacuum showering off 
291         call pyevnt                        ! generate hard scattering
292
293 * generate quenched jets with PYQUEN if nhcel=2 
294         if(nhsel.eq.2.or.nhsel.eq.4) then 
295          ifbp=0 
296          Ap=dble(AW)
297          bfixp=dble(RA*bgen) 
298          call pyquen(Ap,ifbp,bfixp) 
299         end if 
300
301 * fill array of partons
302          nl=nl+n  
303          if(nl.gt.150000-np) goto 51 
304          do i=nl-n+1,nl  
305           ip=i+n-nl                             
306           do j=1,5
307            pl(i,j)=p(ip,j) 
308           end do
309           do j=1,5
310            vl(i,j)=v(ip,j) 
311           end do
312           do j=1,5
313            kl(i,j)=k(ip,j)
314           end do
315           do j=3,5
316            kk=kl(i,j) 
317            if(kk.gt.0) kl(i,j)=kk+nl-n 
318           end do
319          end do 
320  51      continue  
321           
322          call pyexec                         ! hadronization done 
323 c         call pyedit(2)                      ! remove partons & leave hadrons 
324
325 * fill array of final particles
326          nu=nj+n 
327          if(nu.gt.150000-np) then 
328          write(6,*) 'Warning, multiplicity too large! Cut hard part.'   
329           goto 52
330          end if 
331          nj=nu  
332          do i=nj-n+1,nj
333           ip=i+n-nj                            
334           do j=1,5
335            pj(i,j)=p(ip,j) 
336           end do
337           do j=1,5
338            vj(i,j)=v(ip,j) 
339           end do
340           do j=1,5
341            kj(i,j)=k(ip,j)
342           end do
343           do j=3,5
344            kk=kj(i,j) 
345            if(kk.gt.0) then 
346            kj(i,j)=kk+nj-n 
347            end if
348           end do
349          end do               
350       
351        end do 
352  52    njet=ihard-1 
353       end if  
354     
355  
356       return 
357       end 
358 ****************************** END HYHARD **************************      
359       
360 ********************************* HIPSEAR ***************************
361       SUBROUTINE HIPSEAR (fmax,xmin) 
362 * find maximum and 'sufficient minimum' of differential inelasic AA cross 
363 * section as a function of impact paramater (xm, fm are outputs) 
364       real hsin  
365       external hsin  
366       common /hyipar/ bminh,bmaxh,AW,RA,sigin,np  
367       save /hyipar/ 
368      
369       xmin=bmaxh-bminh 
370       fmax=0.
371       do j=1,1000
372       x=bminh+xmin*(j-1)/999.
373       f=hsin(x) 
374        if(f.gt.fmax) then
375         fmax=f
376        endif
377       end do   
378       
379       return
380       end
381 ****************************** END HIPSEAR **************************
382
383 ************************* HARINV **********************************
384       SUBROUTINE HARINV(X,A,F,N,R)                                      
385 * gives interpolation of function F(X) with  arrays A(N) and F(N) 
386       DIMENSION A(N),F(N)                                              
387       IF(X.LT.A(1))GO TO 11                                            
388       IF(X.GT.A(N))GO TO 4                                              
389       K1=1                                                              
390       K2=N                                                              
391  2    K3=K2-K1                                                          
392       IF(K3.LE.1)GO TO 6                                               
393       K3=K1+K3/2                                                        
394       IF(A(K3)-X) 7,8,9                                                 
395  7    K1=K3                                                             
396       GOTO2                                                            
397  9    K2=K3                                                            
398       GOTO2                                                             
399  8    P=F(K3)                                                          
400       RETURN                                                          
401  3    B1=A(K1)                                                          
402       B2=A(K1+1)                                                      
403       B3=A(K1+2)                                                        
404       B4=F(K1)                                                        
405       B5=F(K1+1)                                                        
406       B6=F(K1+2)                                                       
407       R=B4*((X-B2)*(X-B3))/((B1-B2)*(B1-B3))+B5*((X-B1)*(X-B3))/       
408      1 ((B2-B1)*(B2-B3))+B6*((X-B1)*(X-B2))/((B3-B1)*(B3-B2))           
409       RETURN                                                          
410  6    IF(K2.NE.N)GO TO 3                                               
411       K1=N-2                                                            
412       GOTO3                                                            
413  4    C=ABS(X-A(N))                                                     
414       IF(C.LT.0.1) GO TO 5                                           
415       K1=N-2                                                           
416  13   CONTINUE                                                          
417 C13   PRINT 41,X                                                        
418 C41   FORMAT(25H X IS OUT OF THE INTERVAL,3H X=,F15.9)                  
419       GO TO 3                                                           
420  5    R=F(N)                                                           
421       RETURN                                                            
422  11   C=ABS(X-A(1))                                                     
423       IF(C.LT.0.1) GO TO 12                                         
424       K1=1                                                             
425       GOTO 13                                                           
426  12   R=F(1)                                                            
427       RETURN                                                            
428       END                                                              
429 C************************** END HARINV *************************************
430
431 **************************** HSIMPA **********************************
432       SUBROUTINE HSIMPA (A1,B1,H1,REPS1,AEPS1,FUNCT,X,                   
433      1                     AI,AIH,AIABS)                                
434 * calculate intergal of function FUNCT(X) on the interval from A1 to B1 
435       IMPLICIT REAL * 4 ( A-H,O-Z )
436       IMPLICIT INTEGER(I-N)
437       DIMENSION F(7), P(5)                                             
438       H=SIGN ( H1, B1-A1 )                                             
439       S=SIGN (1., H )                                                   
440       A=A1                                                              
441       B=B1                                                              
442       AI=0.0                                                           
443       AIH=0.0                                                           
444       AIABS=0.0                                                        
445       P(2)=4.                                                           
446       P(4)=4.                                                           
447       P(3)=2.                                                           
448       P(5)=1.                                                           
449       IF(B-A)1,2,1                                                      
450  1    REPS=ABS(REPS1)                                                   
451       AEPS=ABS(AEPS1)                                                  
452       DO 3 K=1,7                                                        
453  3    F(K)=10.E16                                                       
454       X=A                                                              
455       C=0.                                                              
456       F(1)=FUNCT(X)/3.                                                  
457  4    X0=X                                                              
458       IF( (X0+4.*H-B)*S)5,5,6                                           
459  6    H=(B-X0)/4.                                                       
460       IF ( H ) 7,2,7                                                   
461  7    DO 8 K=2,7                                                      
462  8    F(K)=10.E16                                                       
463       C=1.                                                           
464  5    DI2=F (1)                                                       
465       DI3=ABS( F(1) )                                                   
466       DO 9 K=2,5                                                       
467       X=X+H                                                           
468       IF((X-B)*S)23,24,24                                              
469  24   X=B                                                              
470  23   IF(F(K)-10.E16)10,11,10                                          
471  11   F(K)=FUNCT(X)/3.                                               
472  10   DI2=DI2+P(K)*F(K)                                                 
473  9    DI3=DI3+P(K)*ABS(F(K))                                            
474       DI1=(F(1)+4.*F(3)+F(5))*2.*H                                      
475       DI2=DI2*H                                                         
476       DI3=DI3*H                                                        
477       IF (REPS) 12,13,12                                               
478  13   IF (AEPS) 12,14,12                                                
479  12   EPS=ABS((AIABS+DI3)*REPS)                                         
480       IF(EPS-AEPS)15,16,16                                              
481  15   EPS=AEPS                                                          
482  16   DELTA=ABS(DI2-DI1)                                               
483       IF(DELTA-EPS)20,21,21                                             
484  20   IF(DELTA-EPS/8.)17,14,14                                          
485  17   H=2.*H                                                            
486       F(1)=F(5)                                                         
487       F(2)=F(6)                                                         
488       F(3)=F(7)                                                         
489       DO 19 K=4,7                                                       
490  19   F(K)=10.E16                                                      
491       GO TO 18                                                         
492  14   F(1)=F(5)                                                         
493       F(3)=F(6)                                                         
494       F(5)=F(7)                                                         
495       F(2)=10.E16                                                       
496       F(4)=10.E16                                                      
497       F(6)=10.E16                                                      
498       F(7)=10.E16                                                      
499  18   DI1=DI2+(DI2-DI1)/15.                                            
500       AI=AI+DI1                                                         
501       AIH=AIH+DI2                                                      
502       AIABS=AIABS+DI3                                                   
503       GO TO 22                                                          
504  21   H=H/2.                                                            
505       F(7)=F(5)                                                        
506       F(6)=F(4)                                                        
507       F(5)=F(3)                                                        
508       F(3)=F(2)                                                         
509       F(2)=10.E16                                                      
510       F(4)=10.E16                                                      
511       X=X0                                                            
512       C=0.                                                             
513       GO TO 5                                                          
514  22   IF(C)2,4,2                                                      
515  2    RETURN                                                        
516       END                                                              
517 ************************* END HSIMPA *******************************
518
519 **************************** HSIMPB **********************************
520       SUBROUTINE HSIMPB (A1,B1,H1,REPS1,AEPS1,FUNCT,X,                   
521      1                     AI,AIH,AIABS)                                
522 * calculate intergal of function FUNCT(X) on the interval from A1 to B1 
523       IMPLICIT REAL * 4 ( A-H,O-Z )
524       IMPLICIT INTEGER(I-N)
525       DIMENSION F(7), P(5)                                             
526       H=SIGN ( H1, B1-A1 )                                             
527       S=SIGN (1., H )                                                   
528       A=A1                                                              
529       B=B1                                                              
530       AI=0.0                                                           
531       AIH=0.0                                                           
532       AIABS=0.0                                                        
533       P(2)=4.                                                           
534       P(4)=4.                                                           
535       P(3)=2.                                                           
536       P(5)=1.                                                           
537       IF(B-A)1,2,1                                                      
538  1    REPS=ABS(REPS1)                                                   
539       AEPS=ABS(AEPS1)                                                  
540       DO 3 K=1,7                                                        
541  3    F(K)=10.E16                                                       
542       X=A                                                              
543       C=0.                                                              
544       F(1)=FUNCT(X)/3.                                                  
545  4    X0=X                                                              
546       IF( (X0+4.*H-B)*S)5,5,6                                           
547  6    H=(B-X0)/4.                                                       
548       IF ( H ) 7,2,7                                                   
549  7    DO 8 K=2,7                                                      
550  8    F(K)=10.E16                                                       
551       C=1.                                                           
552  5    DI2=F (1)                                                       
553       DI3=ABS( F(1) )                                                   
554       DO 9 K=2,5                                                       
555       X=X+H                                                           
556       IF((X-B)*S)23,24,24                                              
557  24   X=B                                                              
558  23   IF(F(K)-10.E16)10,11,10                                          
559  11   F(K)=FUNCT(X)/3.                                               
560  10   DI2=DI2+P(K)*F(K)                                                 
561  9    DI3=DI3+P(K)*ABS(F(K))                                            
562       DI1=(F(1)+4.*F(3)+F(5))*2.*H                                      
563       DI2=DI2*H                                                         
564       DI3=DI3*H                                                        
565       IF (REPS) 12,13,12                                               
566  13   IF (AEPS) 12,14,12                                                
567  12   EPS=ABS((AIABS+DI3)*REPS)                                         
568       IF(EPS-AEPS)15,16,16                                              
569  15   EPS=AEPS                                                          
570  16   DELTA=ABS(DI2-DI1)                                               
571       IF(DELTA-EPS)20,21,21                                             
572  20   IF(DELTA-EPS/8.)17,14,14                                          
573  17   H=2.*H                                                            
574       F(1)=F(5)                                                         
575       F(2)=F(6)                                                         
576       F(3)=F(7)                                                         
577       DO 19 K=4,7                                                       
578  19   F(K)=10.E16                                                      
579       GO TO 18                                                         
580  14   F(1)=F(5)                                                         
581       F(3)=F(6)                                                         
582       F(5)=F(7)                                                         
583       F(2)=10.E16                                                       
584       F(4)=10.E16                                                      
585       F(6)=10.E16                                                      
586       F(7)=10.E16                                                      
587  18   DI1=DI2+(DI2-DI1)/15.                                            
588       AI=AI+DI1                                                         
589       AIH=AIH+DI2                                                      
590       AIABS=AIABS+DI3                                                   
591       GO TO 22                                                          
592  21   H=H/2.                                                            
593       F(7)=F(5)                                                        
594       F(6)=F(4)                                                        
595       F(5)=F(3)                                                        
596       F(3)=F(2)                                                         
597       F(2)=10.E16                                                      
598       F(4)=10.E16                                                      
599       X=X0                                                            
600       C=0.                                                             
601       GO TO 5                                                          
602  22   IF(C)2,4,2                                                      
603  2    RETURN                                                        
604       END                                                              
605 ************************* END HSIMPB *******************************
606
607 * function to calculate differential inelastic AA cross section 
608       real function hsin(x) 
609       external hftaa 
610       real hftaa 
611       common /hyipar/ bminh,bmaxh,AW,RA,sigin,np  
612       save /hyipar/ 
613       sigp=sigin 
614       taapb0=33.16 
615       br=max(1.e-10,0.25*x*x) 
616       hsin=x*(1.-exp(-1.*hftaa(br)*sigp*taapb0)) 
617       return 
618       end 
619
620 * set of functions to calculate number of nucleons-participants at im.par.b 
621       integer function numpar(c) 
622       IMPLICIT REAL * 4 (A-H,O-Z)   
623       real HFUNC1
624       external HFUNC1 
625       common /hynup1/ bp,x    
626       EPS=0.01  
627       A=0. 
628       B=6.28318 
629       H=0.01*(B-A)    
630       bp=c    
631       CALL HSIMPA(A,B,H,EPS,1.E-8,HFUNC1,X,RES,AIH,AIABS)
632       numpar=int(2.*RES) 
633       return 
634       end   
635 *
636       real function HFUNC1(x) 
637       IMPLICIT REAL * 4 (A-H,O-Z) 
638       real HFUNC2 
639       external HFUNC2 
640       common /hyipar/ bminh,bmaxh,AW,RA,sigin,np
641       common /hynup1/ bp,xx
642       save /hyipar/  
643       xx=x 
644       b2=0.5*bp 
645       r1=abs(sqrt(abs(RA*RA-(b2*sin(xx))**2))+b2*cos(xx))
646       r2=abs(sqrt(abs(RA*RA-(b2*sin(xx))**2))-b2*cos(xx))
647       EPS=0.01 
648       A=0. 
649       B=min(r1,r2)
650       H=0.01*(B-A)    
651       CALL HSIMPB(A,B,H,EPS,1.E-8,HFUNC2,Y,RES,AIH,AIABS)
652       HFUNC1=RES 
653       return 
654       end   
655 *      
656       real function HFUNC2(y) 
657       IMPLICIT REAL * 4 (A-H,O-Z) 
658       real hythik 
659       external hythik 
660       common /hyipar/ bminh,bmaxh,AW,RA,sigin,np 
661       common /hynup1/ bp,x 
662       save /hyipar/ 
663       r1=sqrt(abs(y*y+bp*bp/4.+y*bp*cos(x))) 
664       r2=sqrt(abs(y*y+bp*bp/4.-y*bp*cos(x)))
665       s=1.-exp(-0.1*sigin*hythik(r2))
666       HFUNC2=y*hythik(r2)*s 
667       return 
668       end   
669  
670 * to calculate nuclear thikness function 
671        real function hythik(r)   
672        IMPLICIT REAL * 4 (A-H,O-Z) 
673        common /hyipar/ bminh,bmaxh,AW,RA,sigin,np
674        save /hyipar/ 
675        hythik=0.477465*AW*sqrt(abs(RA*RA-r*r))/(RA**3)
676        return
677        end
678
679 * to calculate nuclear overlap function  
680       real function hftaa(br)  
681       common /hyipar/ bminh,bmaxh,AW,RA,sigin,np  
682       save /hyipar/  
683       sbr=sqrt(abs(1.-br)) 
684       taa=1.-br*(1.+(1.-0.25*br)*log(1./br)+2.*(1.-0.25*br)*
685      + (log(1.+sbr)-sbr/(1.+sbr))-0.5*br*(1.-br)/((1.+sbr)**2))    
686       hftaa=taa*((AW/207)**1.33333333) 
687       return 
688       end   
689
690 * function to calculate number of hard parton-parton scatterings with 
691 * Q>x at sqrt{s}=5.5 TeV   
692       integer function numjet(x)      
693       common /hyipar/ bminh,bmaxh,AW,RA,sigin,np
694       common /hyfpar/ bgen,nbcol,npart,npyt,nhyd 
695       save /hyipar/,/hyfpar/ 
696       dimension ptj(30),sgj(30) 
697       data ptj/5.,6.,7.,8.,9.,10.,12.,14.,16.,18.,20.,23.,26.,
698      + 30.,35.,40.,50.,60.,70.,80.,90.,100.,120.,150.,200.,
699      + 250.,300.,400.,450.,500./  
700       data sgj/31.,16.8,9.5,6.1,4.1,2.7,1.4,0.8,0.46,0.3,0.2,
701      + 0.11,6.8e-2,3.9e-2,2.06e-2,1.16e-2,4.46e-3,2.e-3,1.e-3,
702      + 5.3e-4,3.1e-4,1.9e-4,7.7e-5,2.6e-5,5.9e-6,1.6e-6,5.9e-7,
703      + 1.e-7,4.8e-8,2.3e-8/
704       pt=x 
705       i=0 
706  31   i=i+1        
707       if(i.le.30) then 
708        dx=abs(ptj(i)-pt)  
709        if(dx.le.0.1) then 
710         sjet=sgj(i) 
711         goto 32 
712        else
713         goto 31 
714        end if 
715       else 
716        call harinv(pt,ptj,sgj,30,sjet) 
717       end if 
718  32   pjet=sjet/sigin  
719       ijet=0 
720       do i=1,nbcol 
721        if(rlu(0).lt.pjet) ijet=ijet+1 
722       end do 
723       numjet=ijet        
724       return 
725       end   
726
727 * function to generate gauss distribution 
728       real function hgauss(x0,sig)
729  41   u1=rlu(0)
730       u2=rlu(0)
731       v1=2.*u1-1.
732       v2=2.*u2-1.
733       s=v1**2+v2**2
734       if(s.gt.1) go to 41
735       hgauss=v1*sqrt(-2.*log(s)/s)*sig+x0
736       return
737       end
738
739
740 **************************************************************************