*---------------------------------------------------------------------- * * Filename : HYDRO.F * * First version created: 03-AUG-1996 Author : Igor Lokhtin * Last revision : 25-FEB-2005 * *====================================================================== * * Description : Fast event generator for simulation of flow effects in * central and semi-central heavy ion AA collisons at LHC * * Method : N.A.Kruglov, I.P.Lokhtin, L.I.Sarycheva, A.M.Snigirev, * Z. Phys. C 76 (1997) 99-105; * I.P.Lokhtin, L.I.Sarycheva, A.M.Snigirev, * Phys. Lett. B 537 (2002) 261-267; * I.P.Lokhtin, A.M.Snigirev, e-print hep-ph/0312204. * *====================================================================== SUBROUTINE HYDRO(A,ifb,bmin,bmax,bfix,ny) real crosin, gauss, ftaa real AW real A, bmin, bmax, bfix integer ifb, ny external crosin, gauss, ftaa external ludata common /lujets/ n,k(150000,5),p(150000,5),v(150000,5) common /hyipar/ bminh, bmaxh, AW, RA common /hyfpar/ bgen common /hyflow/ ytfl save /lujets/ * reset lujets massives before event generation n=0 do ncl=1,150000 do j=1,5 p(ncl,j)=0. enddo do j=1,5 k(ncl,j)=0 enddo end do * set initial beam paramters: atomic weigth and radius in fm AW=A RA=1.15*AW**0.333333 * pi=3.14159 * generate impact parameter of A-A collision if(ifb.eq.0) then if(bfix.lt.0.) then write(6,*) 'Impact parameter less than zero!' bfix=0. end if if (bfix.gt.2.) then write(6,*) 'Impact parameter larger than two nuclear radius!' bfix=2. end if b1=bfix*RA bgen=bfix else if(bmin.lt.0.) then write(6,*) 'Impact parameter less than zero!' bmin=0. end if if(bmax.gt.2.) then write(6,*) 'Impact parameter larger than two nuclear radius!' bmax=2. end if bminh=bmin bmaxh=bmax call bipsear(fmax1,xmin1) fmax=fmax1 xmin=xmin1 3 bb1=xmin*rlu(0)+bminh ff1=fmax*rlu(0) fb=crosin(bb1) if(ff1.gt.fb) goto 3 b1=bb1*RA bgen=bb1 end if * set flow parameters Tf=0.14 ! freeze-out temperature etmax=5. ! longitudinal flow rapidity if (ytfl.lt.0.01.or.ytfl.gt.3.) ytfl=1. ytmax=ytfl ! transverse flow rapidity * calculate nuclear overlap function for uniform nuceon distribution br=max(1.e-10,0.25*bgen*bgen) rtaa=ftaa(br) * generate total multiplicity in event, np, * with input from mean charged particle density at y=0 for central Pb-Pb events rnp=ny*rtaa*9.625 np=int(rnp) sign=sqrt(rnp) 1 np=int(gauss(rnp,sign)) if(np.gt.150000) goto 1 * generate sort of hadrons (pions, kaons, nucleons) in jetset7* format do ip=1,np !cycle on particles yy=49.*rlu(0) if(yy.lt.11.83333333) then kf=211 elseif(yy.lt.23.66666667) then kf=-211 elseif(yy.lt.35.5) then kf=111 elseif(yy.lt.38.375) then kf=321 elseif(yy.lt.41.25) then kf=-321 elseif(yy.lt.44.125) then kf=310 elseif(yy.lt.47.) then kf=130 elseif(yy.lt.47.5) then kf=2212 elseif(yy.lt.48.) then kf=-2212 elseif(yy.lt.48.5) then kf=2112 else kf=-2112 endif n=n+1 k(n,1)=1 k(n,2)=kf do j=3,5 k(n,j)=0 enddo do j=1,5 v(n,j)=0. enddo kfa=iabs(kf) p(n,5)=ulmass(kfa) * generate uniform distribution in system of a fluid element 2 ep0=-1.*Tf*(log(max(1.e-10,rlu(0)))+log(max(1.e-10,rlu(0)))+ > log(max(1.e-10,rlu(0)))) if(ep0.le.p(n,5)) go to 2 pp0=sqrt(abs(ep0**2-abs(p(n,5)**2))) probt=pp0/ep0 if(rlu(0).gt.probt) go to 2 ctp0=2.*rlu(0)-1. stp0=sqrt(abs(1.-ctp0**2)) php0=2.*pi*rlu(0) * generate coordinates of a fluid element c etaf=etmax*(2.*rlu(0)-1.) ! flat initial eta-spectrum etaf=gauss(0.,etmax) ! gaussian initial eta-spectrum phif=2.*pi*rlu(0) rm1=sqrt(abs(RA*RA-b1*b1/4.*(sin(phif)**2)))+b1*cos(phif)/2. rm2=sqrt(abs(RA*RA-b1*b1/4.*(sin(phif)**2)))-b1*cos(phif)/2. RF=min(rm1,rm2) rrf=RF*RF*rlu(0) * generate four-velocity of a fluid element vradf=sinh(ytmax) sb=RA*RA*(pi-2.*asin(b1/RA/2.))-b1*sqrt(abs(RA*RA-b1*b1/4.)) reff=sqrt(sqrt(sb/pi)*RA) urf=vradf*sqrt(abs(rrf))/reff ! linear transverse profile c urf=vradf*rrf/RA**2 !square transverse profile utf=cosh(etaf)*sqrt(abs(1.+urf*urf)) uzf=sinh(etaf)*sqrt(abs(1.+urf*urf)) * boost in the laboratory system uipi=pp0*(urf*stp0*cos(phif-php0)+uzf*ctp0) bfac=uipi/(utf+1.)+ep0 p(n,1)=pp0*stp0*sin(php0)+urf*sin(phif)*bfac p(n,2)=pp0*stp0*cos(php0)+urf*cos(phif)*bfac p(n,3)=pp0*ctp0+uzf*bfac p(n,4)=sqrt(p(n,1)**2+p(n,2)**2+p(n,3)**2+p(n,5)**2) end do return end ********************************* BIPSEAR *************************** SUBROUTINE BIPSEAR (fmax,xmin) * find maximum and 'sufficient minimum' of differential inelasic AA cross * section as a function of impact paramater (xm, fm are outputs) real crosin external crosin common /hyipar/ bminh, bmaxh, AW, RA xmin=bmaxh-bminh fmax=0. do j=1,1000 x=bminh+xmin*(j-1)/999. f=crosin(x) if(f.gt.fmax) then fmax=f endif end do return end ****************************** END BIPSEAR ************************** * function to calculate differential inelastic AA cross section real function crosin(x) external ftaa real ftaa common /hyipar/ bminh, bmaxh, AW, RA sigp=58. taapb0=33.16 br=max(1.e-10,0.25*x*x) crosin=x*(1.-exp(-1.*ftaa(br)*sigp*taapb0)) return end * function to calculate nuclear overlap function real function ftaa(br) common /hyipar/ bminh, bmaxh, AW, RA sbr=sqrt(abs(1.-br)) taa=1.-br*(1.+(1.-0.25*br)*log(1./br)+2.*(1.-0.25*br)* + (log(1.+sbr)-sbr/(1.+sbr))-0.5*br*(1.-br)/((1.+sbr)**2)) ftaa=taa*((AW/207)**1.33333333) return end * function to generate gauss distribution real function gauss(x0,sig) 41 u1=rlu(0) u2=rlu(0) v1=2.*u1-1. v2=2.*u2-1. s=v1**2+v2**2 if(s.gt.1) go to 41 gauss=v1*sqrt(-2.*log(s)/s)*sig+x0 return end **************************************************************************