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 **************************************************************************