| cb220f83 |
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 | ************************************************************************** |
git://git.uio.no / u / mrichter / AliRoot.git / blame