SUBROUTINE dst(P,dP,y,pT,mT,sign, * T,vT,et0,tau,R,aT,mu,dmu,nmu,mult,n1) c c Calculates the direct invariant one-particle distribution c and derivatives with respect to T, vT, et0, muB, tau, R, aT, ycm, fs c for the model described in LA-UR-96-0782, also in http://xxx.lanl.gov c with designation nucl-th/9603007. c c Integration over eta is done by n1(3)-point Gauss-Jacobi quadratures c with variable xe = eta/et0 and weight function We = 1 - xe^2. c Integration over xp=x/R is done by n1(2)-point Gauss-Jacobi quads c with variable xx = xp/sqrt(1-xe^2), weight Wx = sqrt(1-xx^2). c Here x is the transverse coordinate pointing in the direction of c the pT of the particle under consideration. c Integration over yp=y/R is done by n1(1)-point Gauss-Legendre quads c with variable xy = yp/ypmax, where c ypmax = sqrt(1-xe^2-xp^2) = sqrt(1-xe^2)*sqrt(1-xx^2). c Here y is the other transverse coordinate. Azimuthal symmetry means c the model is unchanged under y -> -y, so the xy integral is from c 0 to 1 and then the result is doubled. c c Notice that the variable y used in the program is not a coordinate, c rather it is the particle's source-frame rapidity y = y_p - y_s. c IMPLICIT NONE INTEGER i,j,k,ne,nx,ny,ii,nmu,imu,ivT,iaT,Ifirst,n1(6) PARAMETER (ne=30, nx=20, ny=20) REAL P,dP(10),y,pT,mT,T,vT,et0,tau,R,aT,mu(4),dmu(4,9),eta,xp,yp REAL xpmax(ne),ypmax,Weix(nx),s(10),s4,fact,sign REAL xe(ne),we(ne),xx(nx),wx(nx),xy(ny),wy(ny),mult(4),vT20,a,b REAL numer,denom,shep,chep,vT2,gamT,expa,gamT0(ne,nx,ny),rhop2 REAL aT0,sqr,sqr0(ne,nx,ny),num1,num2,pi2,fac PARAMETER (pi2 = 6.283185307, fac = 1./pi2**3) SAVE xe,we,xx,wx,xy,wy,xpmax,Weix,gamT0,vT20,Ifirst,sqr0,aT0 DATA vT20,aT0 /2.,-2./ DATA Ifirst /1/ c c Find abcissas and weights for integration (first call only) c if (Ifirst.eq.1) then Ifirst = 0 if (n1(1).gt.ny) then write(*,*) * 'too many points for y-integration, n1(1) > ny in dst' stop endif if (n1(2).gt.nx) then write(*,*) * 'too many points for x-integration, n1(2) > nx in dst' stop endif if (n1(3).gt.ne) then write(*,*) * 'too many points for eta-integration, n1(3) > ne in dst' stop endif call gaujac(xe,we,n1(3),1.,1.) call gaujac(xx,wx,n1(2),.5,.5) call gauleg(0.,1.,xy,wy,n1(1)) do i=1,n1(3) xpmax(i) = sqrt(1.-xe(i)**2) enddo do i=1,n1(2) Weix(i) = sqrt(1.-xx(i)**2) enddo endif c c If vT and aT have not changed, no need to calculate gamT and sqr c vT2 = vT*vT if (vT2.eq.vT20) then ivT = 1 else ivT = 0 vT20 = vT2 endif if (aT.eq.aT0) then iaT = 1 else iaT = 0 aT0 = aT endif c c Integrate over eta (spacetime rapidity) c do ii=1,10 s(ii) = 0. enddo do i=1,n1(3) eta = et0*xe(i) chep = exp(eta-y) shep = .5*(chep - 1./chep) chep = .5*(chep + 1./chep) c c Integrate over xp = x/R c do j=1,n1(2) xp = xpmax(i)*xx(j) ypmax = xpmax(i)*Weix(j) c c Integrate over yp = y/R c do k=1,n1(1) yp = ypmax*xy(k) rhop2 = xp*xp + yp*yp c c transverse flow gamma if (ivT.eq.1) then gamT = gamT0(i,j,k) else gamT = 1./sqrt(1.-vT2*rhop2) gamT0(i,j,k) = gamT endif c c transverse freezeout curvature if (iaT.eq.1) then sqr = sqr0(i,j,k) else sqr = sqrt(1. + aT*rhop2) sqr0(i,j,k) = sqr endif c c p.u(x) and p.n(x) a = gamT*mT*chep/T b = gamT*vT*pT*xp/T num1 = we(i)*wx(j)*wy(k)*sqr*mT num2 = we(i)*wx(j)*wy(k)*pT*xp*tau/R numer = num1*chep - aT*num2 c c Sum over particles in each isospin multiplet c do imu=1,nmu expa = exp(a - b - mu(imu)) denom = (expa + sign)/mult(imu) fact = numer/(denom*(1.+sign/expa)) c P s(1) = s(1) + numer/denom c T*dP/dT s(2) = s(2) + (a - b + dmu(imu,1))*fact c vT*dP/dvT s(3) = s(3) + (gamT**2*(b-vT2*rhop2*a)+dmu(imu,2))*fact c et0*dP/det0 - P s4 = shep*(num1/denom - mT*gamT*fact/T) s(4) = s(4) + eta*s4 c dP/dmuB s(5) = s(5) + dmu(imu,4)*fact c tau*dP/dtau s(6) = s(6) + dmu(imu,5)*fact + + (num1*chep - 2.*aT*num2)/denom c R*dP/dR s(7) = s(7) + dmu(imu,6)*fact + + (2.*num1*chep - aT*num2)/denom c dP/daT s(8) = s(8) + dmu(imu,7)*fact + + (.5*rhop2*num1*chep/sqr**2 - num2)/denom c dP/dycm s(9) = s(9) + dmu(imu,8)*fact + s4 c fs*dP/dfs - |S|*fs^|S|P s(10) = s(10) + dmu(imu,9)*fact enddo enddo enddo enddo c c Overall normalization c fact = 2.*tau*R*R*et0*fac P = fact*s(1) do ii=1,9 dP(ii) = fact*s(ii+1) enddo return END SUBROUTINE distpim(P,dP,y,pT,T,vT,et0,muB,tau,R,aT,muS,muI, * dmuS,dmuI,rw,fs,n1) c c Calculates the negative pion distribution and its derivatives c with respect to T, vT, et0, muB (=muB/T), tau, R, aT, ycm, fs. c c For two-body resonance decays, Mran=0 and c a resonance mres decays into a pion mpi and another particle mX. c c All three-body decays considered have at least two pions in the c final state. So for these, mres decays into two pions and mX. c The other pion and mX can have an invariant mass anywhere in c mpi + mX < M < mres - mpi c So Mran = (mres - 2mpi - mX)/2 c Mave = (mres + mX)/2 c c sgn = -1. for bosons and +1. for fermions c nmu is the number of different baryon (nB), strangeness (nS) and c isospin (nI) states of mres. For each state, mult is the spin c multiplicity of mres multiplied by branching ratio of the decay. c c Integration over left-over invariant mass M in 3-particle decays c is done by n1(6)-point Gauss-Legendre quadratures c Integration over rapidity of the resonance yres is done by c n1(5)-point Gauss-Legendre quadratures c Integration over mTres is done by n1(4)-point Gauss-Chebyshev quads c with variable xk = (mTres-mT_ave)/del_mT and weight function c Wk = 1/sqrt(1-xk^2) c = sqrt(E2)*del_mT/sqrt((pT*pTres)^2-(E0*mres-mTres*E)^2) c Note: \int dmTres f(mTres)/sqrt((pT*pTres)^2-(E0*mres-mTres*E)^2) c = \int_-1^1 dxk/sqrt(1-xk^2) f(mTres(xk))/sqrt(E2) c When pT=0, del_mT=0, so mTres(xk) = mT_ave, and E2(mTres) = E2(mT_ave). c Consequently the xk integral can be done analytically: c \int_-1^1 dxk/sqrt(1-xk^2) = pi c IMPLICIT NONE INTEGER Nd,Ndecay,i,j,k,ii,iM,nj,nk,nka,imu,nM PARAMETER (Nd=18, nM=10,nj=20,nk=20) INTEGER nmu(Nd),nS(4,Nd),nB(4,Nd),niM(Nd),n1(6),Ifirst REAL mres(Nd),E0(Nd),p0(Nd),sgn(Nd),mult(4,Nd),mlt(4) REAL mX(Nd),Mave(Nd),Mran(Nd),sa(10),sM(10),s(10) REAL xiM(nM),wiM(nM),xj(nj),wj(nj),xk(nk),wk,rw,mpi,mpi2 REAL P,dP(10),Pd,dPd(10),y,pT,mT,T,vT,et0,muB,tau,R,muS,muI,aT REAL nI(4,Nd),dmuS(9),dmuI(9),mu(4,Nd),dmu(4,Nd,9),mud(4),ctau(3) REAL muB0,muS0,muI0,dmud(4,9),mres2,C,M,p0M,M2,E02,del_y,fs REAL yres,del_mT,E2,pTres,mTres,pi2,pT2,E,mT_ave,fact,sqrtE2 PARAMETER (pi2 = 6.283185307, mpi = 139.6, mpi2 = mpi**2) c SAVE Ifirst,muB0,muS0,muI0,mres,E0,p0,sgn,niM,nmu,mult,nI,nS,nB SAVE mX,Mave,Mran,xiM,wiM,xj,wj,xk,wk,Ndecay,mu,dmu,ctau DATA Ifirst /1/ DATA muB0,muS0,muI0 /3*0./ DATA mres /547.5, 547.5, 781.9, 781.9, 770., 892., 1232., * 1385., 1385., 1407., * 1116.,1194.,1321., * 958.,1440.,1440.,1525.,1525./ DATA mX /135., 0., 135., mpi, mpi, 494., 939., * 1116., 1194., 1194., * 939.,939.,1116., * 549.,939.,939.,939.,939./ DATA Mave /341.25, 273.75, 458.45, 0., 0., 0., 0., 0., 0., 0., * 0.,0.,0., * 753.5,0.,1189.5,0.,1232./ DATA Mran /66.65, 134.15, 183.85, 0., 0., 0., 0., 0., 0., 0., * 0.,0.,0., * 65.5,0.,111.5,0.,154./ DATA sgn /-1.,-1.,-1.,-1.,-1.,-1.,1.,1.,1.,1., * 1.,1.,1., * -1.,1.,1.,1.,1./ DATA nmu /1,1,1,1,2,2,4,2,4,2,1,2,1,1,2,4,2,4/ DATA mult /.236,0.,0.,0., .049,0.,0.,0., 2.66,0.,0.,0., * .07,0.,0.,0., 3.,3.,0.,0., 2.,2.,0.,0., 3.98,1.33,3.98,1.33, * 3.52,3.52,0.,0., .24,.24,.24,.24, .67,.67,0.,0., * 1.28,0.,0.,0., 2.,.97,0.,0., 2.,0.,0.,0., * .652,0.,0.,0., .8,.8,0.,0., .416,.65,.416,.65, 2.,2.,0.,0., * .922,1.74,.922,1.74/ DATA nI /0.,0.,0.,0., 0.,0.,0.,0., 0.,0.,0.,0., 0.,0.,0.,0., * -1.,0.,0.,0., -.5,-.5,0.,0., -1.5,-.5,-1.5,-.5, * -1.,-1.,0.,0., -1.,0.,-1.,0., 0.,0.,0.,0., * 0.,0.,0.,0., -1.,-1.,0.,0., -.5,0.,0.,0., * 0.,0.,0.,0., -.5,-.5,0.,0., .5,-.5,.5,-.5, -.5,-.5,0.,0., * .5,-.5,.5,-.5/ DATA nS /0,0,0,0, 0,0,0,0, 0,0,0,0, 0,0,0,0, * 0,0,0,0, 1,-1,0,0, 0,0,0,0, * -1,1,0,0, -1,-1,1,1, -1,1,0,0, * -1,0,0,0, -1,1,0,0, -1,0,0,0, * 0,0,0,0, 0,0,0,0, 0,0,0,0, 0,0,0,0, 0,0,0,0/ DATA nB /0,0,0,0, 0,0,0,0, 0,0,0,0, 0,0,0,0, * 0,0,0,0, 0,0,0,0, 1,1,-1,-1, * 1,-1,0,0, 1,1,-1,-1, 1,-1,0,0, * 1,0,0,0, 1,-1,0,0, 1,0,0,0, * 0,0,0,0, 1,-1,0,0, 1,1,-1,-1, 1,-1,0,0, 1,1,-1,-1/ DATA ctau /7.89,3.4,4.91/ c pT2 = pT*pT mT = sqrt(mpi2+pT2) c c Only call this part when the chemical potentials change c if (muB.ne.muB0.or.muS.ne.muS0.or.muI.ne.muI0.or.Ifirst.eq.1) then muB0 = muB muS0 = muS muI0 = muI c c Find abscissas and weights for integration (first time only) c if (Ifirst.eq.1) then Ifirst = 0 if (n1(4).gt.nk) then write(*,*) * 'too many points for mT-integration, n1(4) > nk in distpim' stop endif if (n1(5).gt.nj) then write(*,*) * 'too many points for y-integration, n1(5) > nj in distpim' stop endif if (n1(6).gt.nM) then write(*,*) * 'too many points for M-integration, n1(6) > nM in distpim' stop endif do i=1,n1(4) xk(i) = cos(.5*pi2*(float(i)-.5)/float(n1(4))) enddo wk = .5*pi2/float(n1(4)) call gauleg(-1.,1.,xj,wj,n1(5)) call gauleg(-1.,1.,xiM,wiM,n1(6)) if (rw.lt.0.) then Ndecay = 0 elseif (rw.eq.0.) then Ndecay = 10 else Ndecay = 13 endif c c Calculate E0 and p0 for two-body decays: c do i=1,Ndecay if (Mran(i).eq.0.) then niM(i) = 1 mres2 = mres(i)**2 E0(i) = .5*(mres2 + mpi2 - mX(i)**2)/mres(i) p0(i) = .5*sqrt( (mres2-(mpi+mX(i))**2) * *(mres2-(mpi-mX(i))**2)/mres2 ) else niM(i) = n1(6) endif enddo endif c c Calculate the chemical potential and its derivatives c There is an implicit assumption of positive baryon number c do i=1,Ndecay do imu=1,nmu(i) mu(imu,i) = nI(imu,i)*muI + nB(imu,i)*muB + nS(imu,i)*muS do ii=1,9 dmu(imu,i,ii) = nI(imu,i)*dmuI(ii) + nS(imu,i)*dmuS(ii) enddo dmu(imu,i,4) = dmu(imu,i,4) + float(nB(imu,i))*sign(1.,muB) enddo enddo endif c c c Direct pi- c mud(1) = -muI mlt(1) = 1. do ii=1,9 dmud(1,ii) = -dmuI(ii) enddo CALL dst(P,dP,y,pT,mT,-1.,T,vT,et0,tau,R,aT,mud,dmud,1,mlt,n1) c c c Pi- from resonance decays c do ii=1,10 sa(ii) = 0. enddo do i=1,Ndecay mres2 = mres(i)*mres(i) do imu=1,nmu(i) mud(imu) = mu(imu,i) mlt(imu) = mult(imu,i) do ii=1,9 dmud(imu,ii) = dmu(imu,i,ii) enddo enddo c c In three-body decays, integrate over mass^2 of the other 2 particles c C = 0. do ii=1,10 sM(ii) = 0. enddo do iM=1,niM(i) if (Mran(i).ne.0.) then M = Mave(i) + Mran(i)*xiM(iM) M2 = M*M p0(i) =.5*sqrt((mres2-(mpi+M)**2)*(mres2-(mpi-M)**2)/mres2) p0M = sqrt((M2-(mpi+mX(i))**2)*(M2-(mpi-mX(i))**2))/M E0(i) = .5*(mres2 + mpi2 - M2)/mres(i) endif E02 = E0(i)*E0(i) del_y = log((sqrt(E02+pT2) + p0(i))/mT) c c Integrate over yres (rapidity of the resonance) c do ii=1,10 s(ii) = 0. enddo do j=1,n1(5) yres = y + del_y*xj(j) E = mT*cosh(y-yres) E2 = E*E - pT2 mT_ave = mres(i)*E0(i)*E/E2 if (pT.eq.0.) then del_mT = 0. sqrtE2 = E nka = 1 else del_mT = mres(i)*pT*sqrt(E02 - E2)/E2 sqrtE2 = sqrt(E2) nka = n1(4) endif c c Integrate over mTres c do k=1,nka mTres = mT_ave + del_mT*xk(k) pTres = sqrt(mTres*mTres - mres2) c c Pi- from single-generation resonance decays c CALL dst(Pd,dPd,yres,pTres,mTres,sgn(i), * T,vT,et0,tau,R,aT,mud,dmud,nmu(i),mlt,n1) c fact = wj(j)*mTres/sqrtE2 if ((i.gt.10).and.(i.le.13).and.(pTres.gt.0.)) * fact = fact*(1.-exp(-mres(i)*rw/(pTres*ctau(i-10)))) s(1) = s(1) + fact*Pd do ii=2,10 s(ii) = s(ii) + fact*dPd(ii-1) enddo enddo enddo c c C is the normalization of P(M) for 3-body decays c if (Mran(i).ne.0.) then fact = wiM(iM)*p0M*del_y C = C + wiM(iM)*p0(i)*p0M else fact = del_y/p0(i) endif do ii=1,10 sM(ii) = sM(ii) + fact*s(ii) enddo enddo if (Mran(i).eq.0.) C=1. fact = mres(i)/C c c sa(10) = fs*dsa(1)/dfs (incomplete chemical equilibrium means fs<1) c if (nS(1,i).ne.0) then fact = fact*fs**abs(nS(1,i)) sa(10) = sa(10) + abs(nS(1,i))*fact*sM(1) endif do ii=1,10 sa(ii) = sa(ii) + fact*sM(ii) enddo enddo fact = wk/pi2 if (pT.eq.0.) fact = .5 P = P + fact*sa(1) do ii=1,9 dP(ii) = dP(ii) + fact*sa(ii+1) enddo c c Give correct Jacobians for minimization variables c dP(2) = dP(2)*2.*(1.-vT) dP(3) = dP(3) + P dP(7) = dP(7)*(1.-aT*aT) if (rw.gt.0.) then call distpims(Pd,dPd,y,pT,T,vT,et0,muB,tau,R,aT,muS,muI, * dmuS,dmuI,rw,fs,n1) P = P + Pd do ii =1,9 dP(ii) = dP(ii) + dPd(ii) enddo endif dP(10) = 0. return END SUBROUTINE distpims(P,dP,y,pT,T,vT,et0,muB,tau,R,aT,muS,muI, * dmuS,dmuI,rw,fs,n1) c c Calculates the negative pion distribution and its derivatives c coming from two-generation decays of the following form: c mresa -> mXa + mres -> p + mX c c Explicity, the decays considered are: c sigma0 -> gamma + lambda -> p + pi- c hyper -> pi + lambda -> p + pi- c sigma(1385) -> pi + lambda -> p + pi- c lam/sig(1398) -> pi + sigma -> p + pi c anti_lam/sig(1398) -> pi + anti_sigma- -> anti_n + pi- c c The following approximations have been made to cut down on c integration times: c In hyperon decay, both hyper- and hyper0 are taken to have I=0. c Thus the two isospin states are considered a single resonance. c The same is done with the 3 sigma(1385) isospin states (set I=0). c For decay into sigmas, lambda(1405) and all the isospin states of c sigma(1385) are combined into a resonance lam/sig(1398). c Some of these would decay into sigma0 which decays to lambda and c then to p. We pretend sigma0 goes directly to p + pi-. c c Integration over rapidity of the resonance yres is done by c n1(5)-point Gauss-Legendre quadratures c Integration over mTres is done by n1(4)-point Gauss-Chebyshev quads c IMPLICIT NONE INTEGER Nd,Ndecay,i,j,k,j2,k2,ii,nj,nk,nka,imu PARAMETER (Nd=4,nj=20,nk=20) INTEGER nmu(Nd),nS(2,Nd),nB(2,Nd),n1(6),Ifirst REAL mres(Nd),mX(Nd),E0(Nd),p0(Nd),sgn(Nd),mult(2,Nd),mlt(4) REAL mresa(Nd),mXa(Nd),E0a(Nd),p0a(Nd),pTres2,ctau(Nd) REAL E2a,mTa_av,del_mTa,mTresa,pTresa,fact,E02a,del_ya,yresa,Ea REAL sa(10),sb(10),sc(10),fs,xj(nj),wj(nj),xk(nk),wk,rw REAL P,dP(10),Pd,dPd(10),y,pT,mT,T,vT,et0,muB,tau,R,muS,muI,aT REAL nI(2,Nd),dmuS(9),dmuI(9),mu(2,Nd),dmu(2,Nd,9),mud(4) REAL muB0,muS0,muI0,dmud(4,9),mres2,mres2a,E02,del_y,mpi,mpi2 REAL yres,del_mT,E2,pTres,mTres,pi2,pT2,E,mT_ave,sqrtE2,sqrtE2a PARAMETER (pi2 = 6.283185307, mpi=139.6, mpi2 = mpi**2) c SAVE Ifirst,muB0,muS0,muI0,mresa,mXa,sgn,nmu,mult,nI,nS,nB,ctau SAVE mres,mX,xk,wk,xj,wj,E0,p0,E0a,p0a,mu,dmu DATA Ifirst /1/ DATA muB0,muS0,muI0 /3*0./ DATA mresa /1197.,1318.,1385.,1398./ DATA mXa /0.,137.3,138.1,138.1/ DATA sgn /1., 1., 1., 1./ DATA nmu /1,1,1,2/ DATA mult /2.,0., 2.,0., 6.75,0., 1.88,.554/ DATA nI /0.,0., 0.,0., 0.,0., 0.,0./ DATA nS /-1,0, -2,0, -1,0, -1,1/ DATA nB /1,0, 1,0, 1,0, 1,-1/ DATA ctau /7.89,14.7,7.89,5.3/ DATA mres /1116.,1116.,1116.,1193.5/ DATA mX /939.,939.,939.,939./ c pT2 = pT*pT mT = sqrt(mpi2+pT2) c c Only call this part when the chemical potentials change c if (muB.ne.muB0.or.muS.ne.muS0.or.muI.ne.muI0.or.Ifirst.eq.1) then muB0 = muB muS0 = muS muI0 = muI c c Find abscissas and weights for ingetration (first time only) c if (Ifirst.eq.1) then Ifirst = 0 if (n1(4).gt.nk) then write(*,*) * 'too many points for mT-integration, n1(4) > nk in distpros' stop endif if (n1(5).gt.nj) then write(*,*) * 'too many points for y-integration, n1(5) > nj in distpros' stop endif do i=1,n1(4) xk(i) = cos(.5*pi2*(float(i)-.5)/float(n1(4))) enddo wk = .5*pi2/float(n1(4)) call gauleg(-1.,1.,xj,wj,n1(5)) Ndecay = Nd c c Calculate E0 and p0 for two-body decays: c do i=1,Ndecay mres2 = mres(i)**2 E0(i) = .5*(mres2 + mpi2 - mX(i)**2)/mres(i) p0(i) = .5*sqrt( (mres2-(mpi+mX(i))**2) * *(mres2-(mpi-mX(i))**2)/mres2 ) mres2a = mresa(i)**2 E0a(i) = .5*(mres2a + mres(i)**2 - mXa(i)**2)/mresa(i) p0a(i) = .5*sqrt( (mres2a-(mres(i)+mXa(i))**2) * *(mres2a-(mres(i)-mXa(i))**2)/mres2a ) enddo endif c c Calculate the chemical potential and its derivatives c There is an implicit assumption of positive baryon number c do i=1,Ndecay do imu=1,nmu(i) mu(imu,i) = nI(imu,i)*muI + nB(imu,i)*muB + nS(imu,i)*muS do ii=1,9 dmu(imu,i,ii) = nI(imu,i)*dmuI(ii) + nS(imu,i)*dmuS(ii) enddo dmu(imu,i,4) = dmu(imu,i,4) + nB(imu,i)*sign(1.,muB) enddo enddo endif c c c Sum over resonance decays c do ii=1,10 sa(ii) = 0. enddo do i=1,Ndecay mres2 = mres(i)*mres(i) do imu=1,nmu(i) mud(imu) = mu(imu,i) mlt(imu) = mult(imu,i) do ii=1,9 dmud(imu,ii) = dmu(imu,i,ii) enddo enddo E02 = E0(i)*E0(i) E02a = E0a(i)*E0a(i) del_y = log((sqrt(E02+pT2) + p0(i))/mT) c c Integrate over yres (rapidity of the daughter resonance) c do ii=1,10 sb(ii) = 0. enddo do j=1,n1(5) yres = y + del_y*xj(j) E = mT*cosh(y-yres) E2 = E*E - pT2 mT_ave = mres(i)*E0(i)*E/E2 if (pT.eq.0.) then del_mT = 0. sqrtE2 = E nka = 1 else del_mT = mres(i)*pT*sqrt(E02 - E2)/E2 sqrtE2 = sqrt(E2) nka = n1(4) endif c c Integrate over mTres (daughter resonance) c do k=1,nka mTres = mT_ave + del_mT*xk(k) pTres2 = mTres*mTres - mres2 pTres = sqrt(pTres2) del_ya = log((sqrt(E02a+pTres2)+p0a(i))/mTres) c c Integrate over yresa (rapidity of the parent resonance) c do ii=1,10 sc(ii) = 0. enddo do j2=1,n1(5) yresa = yres + del_ya*xj(j2) Ea = mTres*cosh(yres-yresa) E2a = Ea*Ea - pTres2 mTa_av = mresa(i)*E0a(i)*Ea/E2a del_mTa = mresa(i)*pTres*sqrt(E02a-E2a)/E2a sqrtE2a = sqrt(E2a) c c Integrate over mTresa (parent resonance) c do k2=1,n1(4) mTresa = mTa_av + del_mTa*xk(k2) pTresa = sqrt(mTresa*mTresa - mresa(i)**2) c c Freezeout distribution of the parent resonance c CALL dst(Pd,dPd,yresa,pTresa,mTresa,sgn(i), * T,vT,et0,tau,R,aT,mud,dmud,nmu(i),mlt,n1) fact = wj(j2)*mTresa/sqrtE2a sc(1) = sc(1) + fact*Pd do ii=2,10 sc(ii) = sc(ii) + fact*dPd(ii-1) enddo enddo enddo fact = del_ya*wj(j)*mTres/sqrtE2 * *(1.-exp(-mres(i)*rw/(pTres*ctau(i)))) do ii=1,10 sb(ii) = sb(ii) + fact*sc(ii) enddo enddo enddo fact = mres(i)*mresa(i)*del_y/(p0(i)*p0a(i))*fs**abs(nS(1,i)) sa(10) = sa(10) + abs(nS(1,i))*fact*sb(1) do ii=1,10 sa(ii) = sa(ii) + fact*sb(ii) enddo enddo fact = (wk/pi2)**2 if (pT.eq.0.) fact = .5*wk/pi2 P = fact*sa(1) do ii=1,9 dP(ii) = fact*sa(ii+1) enddo c c Jacobians for minimization variables c dP(2) = dP(2)*2.*(1.-vT) dP(3) = dP(3) + P dP(7) = dP(7)*(1.-aT*aT) dP(10) = 0. return END SUBROUTINE distkp(P,dP,y,pT,T,vT,et0,muB,tau,R,aT,muS,muI, * dmuS,dmuI,rw,fs,n1) c c Calculates the positive kaon distribution and its derivatives c with respect to T, vT, et0, muB (=muB/T), tau, R, aT, ycm, fs. c c A resonance mres decays into a kaon mK and another particle mX. c sgn = -1. for bosons and +1. for fermions c nmu is the number of different baryon (nB), strangeness (nS) and c isospin (nI) states of mres. For each state, mult is the spin c multiplicity of mres multiplied by branching ratio of the decay. c c Integration over rapidity of the resonance yres is done by c n1(5)-point Gauss-Legendre quadratures c Integration over mTres is done by n1(4)-point Gauss-Chebyshev quads c with variable xk = (mTres-mT_ave)/del_mT and weight function c Wk = 1/sqrt(1-xk^2) c = sqrt(E2)*del_mT/sqrt((pT*pTres)^2-(E0*mres-mTres*E)^2) c Note: \int dmTres f(mTres)/sqrt((pT*pTres)^2-(E0*mres-mTres*E)^2) c = \int_-1^1 dxk/sqrt(1-xk^2) f(mTres(xk))/sqrt(E2) c When pT=0, del_mT=0, so mTres(xk) = mT_ave, and E2(mTres) = E2(mT_ave). c Consequently the xk integral can be done analytically: c \int_-1^1 dxk/sqrt(1-xk^2) = pi c IMPLICIT NONE INTEGER Nd,Ndecay,i,j,k,ii,nj,nk,nka,imu,Ifirst,n1(6) PARAMETER (Nd=1, nj=20, nk=20) INTEGER nmu(Nd),nS(2,Nd) REAL mres(Nd),E0(Nd),p0(Nd),sgn(Nd),mult(2,Nd),mlt(4),mX(Nd) REAL sa(10),s(10),xj(nj),wj(nj),xk(nk),wk,rw,fs REAL P,dP(10),Pd,dPd(10),y,pT,mT,T,vT,et0,tau,R,aT,muS,muI,muB REAL nI(2,Nd),dmuS(9),dmuI(9),mu(2,Nd),dmu(2,Nd,9),mud(4) REAL muB0,muS0,muI0,dmud(4,9),mres2,E02,del_y,mK,mK2 REAL yres,del_mT,E2,pTres,mTres,pi2,pT2,E,mT_ave,fact,sqrtE2 PARAMETER (pi2 = 6.283185307, mK=494., mK2 = mK**2) c SAVE Ifirst,muB0,muS0,muI0,mres,E0,p0,sgn,nmu,mult SAVE mX,nI,nS,xj,wj,xk,wk,Ndecay,mu,dmu DATA Ifirst /1/ DATA muB0,muS0,muI0 /3*0./ DATA mres /892./ DATA mX /137.3/ DATA sgn /-1./ DATA nmu /2/ DATA mult /1.,2./ DATA nI /.5,-.5/ DATA nS /1,1/ c pT2 = pT*pT mT = sqrt(mK2+pT2) c c Only call this part when the chemical potentials change c if (muB.ne.muB0.or.muS.ne.muS0.or.muI.ne.muI0.or.Ifirst.eq.1) then muB0 = muB muS0 = muS muI0 = muI c c Calculate weights and abscissas for integration (first time only) c if (Ifirst.eq.1) then Ifirst = 0 if (n1(4).gt.nk) then write(*,*) * 'too many points for mT-integration, n1(4) > nk in distkp' stop endif if (n1(5).gt.nj) then write(*,*) * 'too many points for y-integration, n1(5) > nj in distkp' stop endif do i=1,n1(4) xk(i) = cos(.5*pi2*(float(i)-.5)/float(n1(4))) enddo wk = .5*pi2/float(n1(4)) call gauleg(-1.,1.,xj,wj,n1(5)) if (rw.lt.0.) then Ndecay = 0 else Ndecay = Nd endif c c Calculate E0 and p0 for two-body decays: c do i=1,Ndecay mres2 = mres(i)**2 E0(i) = .5*(mres2 + mK2 - mX(i)**2)/mres(i) p0(i) = .5*sqrt( (mres2-(mK+mX(i))**2) * *(mres2-(mK-mX(i))**2)/mres2 ) enddo endif c c Calculate the chemical potential and its derivatives c do i=1,Ndecay do imu=1,nmu(i) mu(imu,i) = nI(imu,i)*muI + nS(imu,i)*muS do ii=1,9 dmu(imu,i,ii) = nI(imu,i)*dmuI(ii) + nS(imu,i)*dmuS(ii) enddo enddo enddo endif c c c Direct K+ c mud(1) = muS + .5*muI mlt(1) = 1. do ii=1,9 dmud(1,ii) = dmuS(ii) + .5*dmuI(ii) enddo CALL dst(P,dP,y,pT,mT,-1.,T,vT,et0,tau,R,aT,mud,dmud,1,mlt,n1) P = fs*P dP(9) = dP(9) + P c c c K+ from resonance decays c do ii=1,10 sa(ii) = 0. enddo do i=1,Ndecay mres2 = mres(i)*mres(i) do imu=1,nmu(i) mud(imu) = mu(imu,i) mlt(imu) = mult(imu,i) do ii=1,9 dmud(imu,ii) = dmu(imu,i,ii) enddo enddo E02 = E0(i)*E0(i) del_y = log((sqrt(E02+pT2) + p0(i))/mT) c c Integrate over yres (rapidity of the resonance) c do ii=1,10 s(ii) = 0. enddo do j=1,n1(5) yres = y + del_y*xj(j) E = mT*cosh(y-yres) E2 = E*E - pT2 mT_ave = mres(i)*E0(i)*E/E2 if (pT.eq.0.) then del_mT = 0. sqrtE2 = E nka = 1 else del_mT = mres(i)*pT*sqrt(E02 - E2)/E2 sqrtE2 = sqrt(E2) nka = n1(4) endif c c Integrate over mTres c do k=1,nka mTres = mT_ave + del_mT*xk(k) pTres = sqrt(mTres*mTres - mres2) c c K+ from two-body resonance decays c CALL dst(Pd,dPd,yres,pTres,mTres,sgn(i), * T,vT,et0,tau,R,aT,mud,dmud,nmu(i),mlt,n1) c fact = wj(j)*mTres/sqrtE2 s(1) = s(1) + fact*Pd do ii=2,10 s(ii) = s(ii) + fact*dPd(ii-1) enddo enddo enddo fact = mres(i)*del_y/p0(i) c c sa(10) = fs*dsa(1)/dfs (incomplete chemical equilibrium means fs<1) c if (nS(1,i).ne.0) then fact = fact*fs**abs(nS(1,i)) sa(10) = sa(10) + abs(nS(1,i))*fact*s(1) endif do ii=1,10 sa(ii) = sa(ii) + fact*s(ii) enddo enddo fact = wk/pi2 if (pT.eq.0.) fact = .5 P = P + fact*sa(1) do ii=1,9 dP(ii) = dP(ii) + fact*sa(ii+1) enddo c c Jacobians for minimization variables c dP(2) = dP(2)*2.*(1.-vT) dP(3) = dP(3) + P dP(7) = dP(7)*(1.-aT*aT) dP(10) = 0. return END SUBROUTINE distpro(P,dP,y,pT,T,vT,et0,muB,tau,R,aT,muS,muI, * dmuS,dmuI,rw,fs,n1) c c Calculates the proton distribution and its derivatives c with respect to T, vT, et0, muB (=muB/T), tau, R, aT, ycm, fs. c c A resonance mres decays into a proton mp and another particle mX. c sgn = -1. for bosons and +1. for fermions c nmu is the number of different baryon (nB), strangeness (nS) and c isospin (nI) states of mres. For each state, mult is the spin c multiplicity of mres multiplied by branching ratio of the decay. c c Integration over rapidity of the resonance yres is done by c n1(5)-point Gauss-Legendre quadratures c Integration over mTres is done by n1(4)-point Gauss-Chebyshev quads c with variable xk = (mTres-mT_ave)/del_mT and weight function c Wk = 1/sqrt(1-xk^2) c = sqrt(E2)*del_mT/sqrt((pT*pTres)^2-(E0*mres-mTres*E)^2) c Note: \int dmTres f(mTres)/sqrt((pT*pTres)^2-(E0*mres-mTres*E)^2) c = \int_-1^1 dxk/sqrt(1-xk^2) f(mTres(xk))/sqrt(E2) c When pT=0, del_mT=0, so mTres(xk) = mT_ave, and E2(mTres) = E2(mT_ave). c Consequently the xk integral can be done analytically: c \int_-1^1 dxk/sqrt(1-xk^2) = pi c IMPLICIT NONE INTEGER Nd,Ndecay,i,j,k,ii,nj,nk,nka,imu,Ifirst,n1(6) PARAMETER (Nd=3, nj=20, nk=20) INTEGER nmu(Nd),nS(3,Nd),nB(3,Nd) REAL mres(Nd),E0(Nd),p0(Nd),sgn(Nd),mult(3,Nd),mlt(4),mX(Nd) REAL sa(10),s(10),xj(nj),wj(nj),xk(nk),wk,rw,fs REAL P,dP(10),Pd,dPd(10),y,pT,mT,T,vT,et0,tau,R,aT,muS,muI,muB REAL nI(3,Nd),dmuS(9),dmuI(9),mu(3,Nd),dmu(3,Nd,9),mud(4) REAL muB0,muS0,muI0,dmud(4,9),mres2,E02,del_y,mp,mp2,ctau(2) REAL yres,del_mT,E2,pTres,mTres,pi2,pT2,E,mT_ave,fact,sqrtE2 PARAMETER (pi2 = 6.283185307, mp=939., mp2 = mp**2) c SAVE Ifirst,muB0,muS0,muI0,mres,mX,E0,p0,sgn,nmu,mult SAVE nI,nS,nB,xj,wj,xk,wk,Ndecay,mu,dmu,ctau DATA Ifirst /1/ DATA muB0,muS0,muI0 /3*0./ DATA mres /1232.,1116.,1189./ DATA mX /137.3,139.6,135./ DATA sgn /1.,1.,1./ DATA nmu /3,1,1/ DATA mult /3.98,2.65,1.33, 1.28,0.,0., 1.03,0.,0./ DATA nI /1.5,.5,-.5, 0.,0.,0., 1.,0.,0./ DATA nS /0,0,0, -1,0,0, -1,0,0/ DATA nB /1,1,1, 1,0,0, 1,0,0/ DATA ctau /7.89,2.396/ c pT2 = pT*pT mT = sqrt(mp2+pT2) c c Only call this part when the chemical potentials change c if (muB.ne.muB0.or.muS.ne.muS0.or.muI.ne.muI0.or.Ifirst.eq.1) then muB0 = muB muS0 = muS muI0 = muI c c Calcualte abscissas and weights for integration (first time only) c if (Ifirst.eq.1) then Ifirst = 0 if (n1(4).gt.nk) then write(*,*) * 'too many points for mT-integration, n1(4) > nk in distpro' stop endif if (n1(5).gt.nj) then write(*,*) * 'too many points for y-integration, n1(5) > nj in distpro' stop endif do i=1,n1(4) xk(i) = cos(.5*pi2*(float(i)-.5)/float(n1(4))) enddo wk = .5*pi2/float(n1(4)) call gauleg(-1.,1.,xj,wj,n1(5)) if (rw.lt.0.) then Ndecay = 0 elseif (rw.eq.0.) then Ndecay = 1 else Ndecay = Nd endif c c Calculate E0 and p0 for two-body decays: c do i=1,Ndecay mres2 = mres(i)**2 E0(i) = .5*(mres2 + mp2 - mX(i)**2)/mres(i) p0(i) = .5*sqrt( (mres2-(mp+mX(i))**2) * *(mres2-(mp-mX(i))**2)/mres2 ) enddo endif c c Calculate the chemical potential and its derivatives c There is an implicit assumption of positive baryon number c do i=1,Ndecay do imu=1,nmu(i) mu(imu,i) = nI(imu,i)*muI + nB(imu,i)*muB + nS(imu,i)*muS do ii=1,9 dmu(imu,i,ii) = nI(imu,i)*dmuI(ii) + nS(imu,i)*dmuS(ii) enddo dmu(imu,i,4) = dmu(imu,i,4) + float(nB(imu,i))*sign(1.,muB) enddo enddo endif c c c Direct protons c mud(1) = muB + .5*muI mlt(1) = 2. do ii=1,9 dmud(1,ii) = .5*dmuI(ii) enddo dmud(1,4) = dmud(1,4) + 1. CALL dst(P,dP,y,pT,mT,1.,T,vT,et0,tau,R,aT,mud,dmud,1,mlt,n1) c c c Sum over resonance decays c do ii=1,10 sa(ii) = 0. enddo do i=1,Ndecay mres2 = mres(i)*mres(i) do imu=1,nmu(i) mud(imu) = mu(imu,i) mlt(imu) = mult(imu,i) do ii=1,9 dmud(imu,ii) = dmu(imu,i,ii) enddo enddo E02 = E0(i)*E0(i) del_y = log((sqrt(E02+pT2) + p0(i))/mT) c c Integrate over yres (rapidity of the resonance) c do ii=1,10 s(ii) = 0. enddo do j=1,n1(5) yres = y + del_y*xj(j) E = mT*cosh(y-yres) E2 = E*E - pT2 mT_ave = mres(i)*E0(i)*E/E2 if (pT.eq.0.) then del_mT = 0. sqrtE2 = E nka = 1 else del_mT = mres(i)*pT*sqrt(E02 - E2)/E2 sqrtE2 = sqrt(E2) nka = n1(4) endif c c Integrate over mTres c do k=1,nka mTres = mT_ave + del_mT*xk(k) pTres = sqrt(mTres*mTres - mres2) c c Freezeout distribution of the resonance c CALL dst(Pd,dPd,yres,pTres,mTres,sgn(i), * T,vT,et0,tau,R,aT,mud,dmud,nmu(i),mlt,n1) fact = wj(j)*mTres/sqrtE2 if ((i.gt.1).and.(pTres.gt.0.)) * fact = fact*(1.-exp(-mres(i)*rw/(pTres*ctau(i-1)))) s(1) = s(1) + fact*Pd do ii=2,10 s(ii) = s(ii) + fact*dPd(ii-1) enddo enddo enddo fact = mres(i)*del_y/p0(i) c c sa(10) = fs*dsa(1)/dfs (incomplete chemical equilibrium means fs<1) c if (nS(1,i).ne.0) then fact = fact*fs**abs(nS(1,i)) sa(10) = sa(10) + abs(nS(1,i))*fact*s(1) endif do ii=1,10 sa(ii) = sa(ii) + fact*s(ii) enddo enddo fact = wk/pi2 if (pT.eq.0.) fact = .5 P = P + fact*sa(1) do ii=1,9 dP(ii) = dP(ii) + fact*sa(ii+1) enddo c c Jacobians for minimization variables c dP(2) = dP(2)*2.*(1.-vT) dP(3) = dP(3) + P dP(7) = dP(7)*(1.-aT*aT) if (rw.gt.0.) then call distpros(Pd,dPd,y,pT,T,vT,et0,muB,tau,R,aT,muS,muI, * dmuS,dmuI,rw,fs,n1) P = P + Pd do ii =1,9 dP(ii) = dP(ii) + dPd(ii) enddo endif dP(10) = 0. return END SUBROUTINE distpros(P,dP,y,pT,T,vT,et0,muB,tau,R,aT,muS,muI, * dmuS,dmuI,rw,fs,n1) c c Calculates the proton distribution and its derivatives c coming from two-generation decays of the following form: c mresa -> mXa + mres -> p + mX c c Explicity, the decays considered are: c sigma0 -> gamma + lambda -> p + pi- c hyper -> pi + lambda -> p + pi- c sigma(1385) -> pi + lambda -> p + pi- c lam/sig(1398) -> pi + sigma -> p + pi c c The following approximations have been made to cut down on c integration times: c In hyperon decay, both hyper- and hyper0 are taken to have I=0. c Thus the two isospin states are considered a single resonance. c The same is done with the 3 sigma(1385) isospin states (set I=0). c For decay into sigmas, lambda(1405) and all the isospin states of c sigma(1385) are combined into a resonance lam/sig(1398). c Some of these would decay into sigma0 which decays to lambda and c then to p. We pretend sigma0 goes directly to p + pi-. c c Integration over rapidity of the resonance yres is done by c n1(5)-point Gauss-Legendre quadratures c Integration over mTres is done by n1(4)-point Gauss-Chebyshev quads c IMPLICIT NONE INTEGER Nd,Ndecay,i,j,k,j2,k2,ii,nj,nk,nka,imu PARAMETER (Nd=4,nj=20,nk=20) INTEGER nmu(Nd),nS(1,Nd),nB(1,Nd),n1(6),Ifirst REAL mres(Nd),mX(Nd),E0(Nd),p0(Nd),sgn(Nd),mult(1,Nd),mlt(4) REAL mresa(Nd),mXa(Nd),E0a(Nd),p0a(Nd),pTres2,ctau(Nd) REAL E2a,mTa_av,del_mTa,mTresa,pTresa,fact,E02a,del_ya,yresa,Ea REAL sa(10),sb(10),sc(10),fs,xj(nj),wj(nj),xk(nk),wk,rw REAL P,dP(10),Pd,dPd(10),y,pT,mT,T,vT,et0,muB,tau,R,muS,muI,aT REAL nI(1,Nd),dmuS(9),dmuI(9),mu(1,Nd),dmu(1,Nd,9),mud(4) REAL muB0,muS0,muI0,dmud(4,9),mres2,mres2a,E02,del_y,mp,mp2 REAL yres,del_mT,E2,pTres,mTres,pi2,pT2,E,mT_ave,sqrtE2,sqrtE2a PARAMETER (pi2 = 6.283185307, mp=939., mp2 = mp**2) c SAVE Ifirst,muB0,muS0,muI0,mresa,mXa,sgn,nmu,mult,nI,nS,nB,ctau SAVE mres,mX,xk,wk,xj,wj,E0,p0,E0a,p0a,mu,dmu DATA Ifirst /1/ DATA muB0,muS0,muI0 /3*0./ DATA mresa /1197.,1318.,1385.,1398./ DATA mXa /0.,137.3,138.1,138.1/ DATA sgn /1., 1., 1., 1./ DATA nmu /1,1,1,1/ DATA mult /2.,2.,6.75,1.32/ DATA nI /0.,0.,0.,0./ DATA nS /-1,-2,-1,-1/ DATA nB /1,1,1,1/ DATA ctau /7.89,14.7,7.89,5.14/ DATA mres /1116.,1116.,1116.,1191./ DATA mX /139.6,139.6,139.6,137./ c pT2 = pT*pT mT = sqrt(mp2+pT2) c c Only call this part when the chemical potentials change c if (muB.ne.muB0.or.muS.ne.muS0.or.muI.ne.muI0.or.Ifirst.eq.1) then muB0 = muB muS0 = muS muI0 = muI c c Find abscissas and weights for ingetration (first time only) c if (Ifirst.eq.1) then Ifirst = 0 if (n1(4).gt.nk) then write(*,*) * 'too many points for mT-integration, n1(4) > nk in distpros' stop endif if (n1(5).gt.nj) then write(*,*) * 'too many points for y-integration, n1(5) > nj in distpros' stop endif do i=1,n1(4) xk(i) = cos(.5*pi2*(float(i)-.5)/float(n1(4))) enddo wk = .5*pi2/float(n1(4)) call gauleg(-1.,1.,xj,wj,n1(5)) Ndecay = Nd c c Calculate E0 and p0 for two-body decays: c do i=1,Ndecay mres2 = mres(i)**2 E0(i) = .5*(mres2 + mp2 - mX(i)**2)/mres(i) p0(i) = .5*sqrt( (mres2-(mp+mX(i))**2) * *(mres2-(mp-mX(i))**2)/mres2 ) mres2a = mresa(i)**2 E0a(i) = .5*(mres2a + mres(i)**2 - mXa(i)**2)/mresa(i) p0a(i) = .5*sqrt( (mres2a-(mres(i)+mXa(i))**2) * *(mres2a-(mres(i)-mXa(i))**2)/mres2a ) enddo endif c c Calculate the chemical potential and its derivatives c There is an implicit assumption of positive baryon number c do i=1,Ndecay do imu=1,nmu(i) mu(imu,i) = nI(imu,i)*muI + nB(imu,i)*muB + nS(imu,i)*muS do ii=1,9 dmu(imu,i,ii) = nI(imu,i)*dmuI(ii) + nS(imu,i)*dmuS(ii) enddo dmu(imu,i,4) = dmu(imu,i,4) + nB(imu,i)*sign(1.,muB) enddo enddo endif c c c Sum over resonance decays c do ii=1,10 sa(ii) = 0. enddo do i=1,Ndecay mres2 = mres(i)*mres(i) do imu=1,nmu(i) mud(imu) = mu(imu,i) mlt(imu) = mult(imu,i) do ii=1,9 dmud(imu,ii) = dmu(imu,i,ii) enddo enddo E02 = E0(i)*E0(i) E02a = E0a(i)*E0a(i) del_y = log((sqrt(E02+pT2) + p0(i))/mT) c c Integrate over yres (rapidity of the daughter resonance) c do ii=1,10 sb(ii) = 0. enddo do j=1,n1(5) yres = y + del_y*xj(j) E = mT*cosh(y-yres) E2 = E*E - pT2 mT_ave = mres(i)*E0(i)*E/E2 if (pT.eq.0.) then del_mT = 0. sqrtE2 = E nka = 1 else del_mT = mres(i)*pT*sqrt(E02 - E2)/E2 sqrtE2 = sqrt(E2) nka = n1(4) endif c c Integrate over mTres (daughter resonance) c do k=1,nka mTres = mT_ave + del_mT*xk(k) pTres2 = mTres*mTres - mres2 pTres = sqrt(pTres2) del_ya = log((sqrt(E02a+pTres2)+p0a(i))/mTres) c c Integrate over yresa (rapidity of the parent resonance) c do ii=1,10 sc(ii) = 0. enddo do j2=1,n1(5) yresa = yres + del_ya*xj(j2) Ea = mTres*cosh(yres-yresa) E2a = Ea*Ea - pTres2 mTa_av = mresa(i)*E0a(i)*Ea/E2a del_mTa = mresa(i)*pTres*sqrt(E02a-E2a)/E2a sqrtE2a = sqrt(E2a) c c Integrate over mTresa (parent resonance) c do k2=1,n1(4) mTresa = mTa_av + del_mTa*xk(k2) pTresa = sqrt(mTresa*mTresa - mresa(i)**2) c c Freezeout distribution of the parent resonance c CALL dst(Pd,dPd,yresa,pTresa,mTresa,sgn(i), * T,vT,et0,tau,R,aT,mud,dmud,nmu(i),mlt,n1) fact = wj(j2)*mTresa/sqrtE2a sc(1) = sc(1) + fact*Pd do ii=2,10 sc(ii) = sc(ii) + fact*dPd(ii-1) enddo enddo enddo fact = del_ya*wj(j)*mTres/sqrtE2 * *(1.-exp(-mres(i)*rw/(pTres*ctau(i)))) do ii=1,10 sb(ii) = sb(ii) + fact*sc(ii) enddo enddo enddo fact = mres(i)*mresa(i)*del_y/(p0(i)*p0a(i))*fs**abs(nS(1,i)) sa(10) = sa(10) + abs(nS(1,i))*fact*sb(1) do ii=1,10 sa(ii) = sa(ii) + fact*sb(ii) enddo enddo fact = (wk/pi2)**2 if (pT.eq.0.) fact = .5*wk/pi2 P = fact*sa(1) do ii=1,9 dP(ii) = fact*sa(ii+1) enddo c c Jacobians for minimization variables c dP(2) = dP(2)*2.*(1.-vT) dP(3) = dP(3) + P dP(7) = dP(7)*(1.-aT*aT) dP(10) = 0. return END SUBROUTINE dst2(P,dP,Y,pTx,pTy,mT,q,sgn, * T,vT,et0,tau,R,aT,mu,dmu,nmu,mult,n2,pqf) c c Calculates the fourier transform term of the 2-particle distribution c and derivatives with respect to T, vT, et0, muB, tau, R, aT, ycm, fs. c P(1) is the real part. P(2) is the imaginary part. c c Integration over eta is done by n2(3)-point Gauss-Jacobi quadratures c with variable xe = eta/et0 and weight function We = 1 - xe^2 c Integration over xp=x/R is done by n2(2)-point Gauss-Jacobi quads c with variable xx = xp/sqrt(1-xe^2), weight Wx = sqrt(1-xx^2) c Integration over yp=y/R is done by n2(1)-point Gauss-Legendre quads c with variable xy = yp/ypmax c ypmax = sqrt(1-xe^2-xp^2) = sqrt(1-xe^2)*sqrt(1-xx^2) c IMPLICIT NONE INTEGER i,j,k,ne,nx,ny,ii,nmu,imu,ivT,iy,it,n2(6),Ifirst,iaT,iq PARAMETER (ne=30, nx=20, ny=20) REAL P(2),dP(9,2),y,mT,T,vT,et0,tau,R,mu(4),dmu(4,9),eta,xp,yp REAL xpmax(ne),ypmax,Weix(nx),s(10,2),pi2,fact REAL xe(ne),we(ne),xx(nx),wx(nx),xy(ny),wy(ny),mult(4),vT20,a,b REAL numer,denom,shep,chep,vT2,gamT,expa,gamT0(ne,nx,ny),rhop2,sgn REAL chep0(ne),shep0(ne),y0,et00,q(4),tau0,time,z,time0(ne),z0(ne) REAL cosqx,sinqx,s1,s2,pTx,pTy,aT,aT0,sqr,sqr0(ne,nx,ny),qx0,qx,qy REAL num1,num2,q0(4),R0,cosqx0(ne,nx,ny),cosqy0(ne,nx,ny) REAL sinqx0(ne,nx,ny),sinqy0(ne,nx,ny),cosqy,sinqy,pqf(4) REAL c1,c2,dc1,dc2,reqx,imqx,s3,dreqx,dimqx,fac PARAMETER (pi2 = 6.283185307, fac = 1./pi2**3) SAVE xe,we,xx,wx,xy,wy,xpmax,Weix,gamT0,vT20,y0,chep0,shep0,et00 SAVE tau0,time0,z0,Ifirst,aT0,sqr0,R0,q0,cosqx0,sinqx0 SAVE cosqy0,sinqy0 DATA vT20,y0,et00,tau0,aT0,R0 /2.,100000.,-1.,-1.,-2.,-1./ DATA q0 /4*31415927./ DATA Ifirst /1/ c c Find abcissas and weights for integration (first call only) c if (Ifirst.eq.1) then Ifirst = 0 if (n2(1).gt.ny) then write(*,*) * 'too many points for y-integration, n2(1) > ny in dst2' stop endif if (n2(2).gt.nx) then write(*,*) * 'too many points for x-integration, n2(2) > nx in dst2' stop endif if (n2(3).gt.ne) then write(*,*) * 'too many points for eta-integration, n2(3) > ne in dst2' stop endif call gaujac(xe,we,n2(3),1.,1.) call gaujac(xx,wx,n2(2),.5,.5) call gauleg(-1.,1.,xy,wy,n2(1)) do i=1,n2(3) xpmax(i) = sqrt(1.-xe(i)**2) enddo do i=1,n2(2) Weix(i) = sqrt(1.-xx(i)**2) enddo endif c c If parameters have not changed, do not calculate gamT, sqr, z, time c vT2 = vT*vT if (vT2.eq.vT20) then ivT = 1 else ivT = 0 vT20 = vT2 endif if (aT.eq.aT0) then iaT = 1 else iaT = 0 aT0 = aT endif if (et0.eq.et00) then if (y.eq.y0) then iy = 1 else iy = 0 y0 = y endif if (tau.eq.tau0) then it = 1 else it = 0 tau0 = tau endif else iy = 0 it = 0 et00 = et0 y0 = y tau0 = tau endif if ((q(4).eq.q0(4)).and.(q(3).eq.q0(3)).and.(q(2).eq.q0(2)) * .and.(q(1).eq.q0(1)).and.(R.eq.R0)) then iq = 1 else iq = 0 do i=1,4 q0(i) = q(i) enddo R0 = R endif c c Integrate over eta (spacetime rapidity) c do ii=1,10 s(ii,1) = 0. s(ii,2) = 0. enddo do i=1,n2(3) eta = et0*xe(i) if (iy.eq.1) then chep = chep0(i) shep = shep0(i) else chep = exp(eta-y) shep = .5*(chep - 1./chep) chep = .5*(chep + 1./chep) chep0(i) = chep shep0(i) = shep endif if (it.eq.1) then time = time0(i) z = z0(i) else time = tau*cosh(eta) z = tau*sinh(eta) time0(i) = time z0(i) = z endif c c Integrate over xp = x/R c do j=1,n2(2) xp = xpmax(i)*xx(j) ypmax = xpmax(i)*Weix(j) c c Integrate over yp = y/R c do k=1,n2(1) yp = ypmax*xy(k) rhop2 = xp*xp + yp*yp c c transverse flow gamma if (ivT.eq.1) then gamT = gamT0(i,j,k) else gamT = 1./sqrt(1.-vT2*rhop2) gamT0(i,j,k) = gamT endif c c transverse freezeout curvature if (iaT.eq.1) then sqr = sqr0(i,j,k) else sqr = sqrt(1. + aT*rhop2) sqr0(i,j,k) = sqr endif c c p.u(x), p.n(x), and q.x a = gamT*mT*chep/T b = gamT*vT*(xp*pTx+yp*pTy)/T num1 = we(i)*wx(j)*wy(k)*sqr*mT num2 = we(i)*wx(j)*wy(k)*(pTx*xp + pTy*yp)*tau/R numer = num1*chep - aT*num2 qx0 = sqr*(q(4)*time - q(3)*z) qx = q(1)*xp*R qy = q(2)*yp*R c c Calculation of fourier factors if ((iq.eq.1).and.(it.eq.1).and.(iaT.eq.1)) then cosqx = cosqx0(i,j,k) sinqx = sinqx0(i,j,k) cosqy = cosqy0(i,j,k) sinqy = sinqy0(i,j,k) else cosqx = cos(qx0 - qx) sinqx = sin(qx0 - qx) cosqy = cos(qy) sinqy = sin(qy) cosqx0(i,j,k) = cosqx sinqx0(i,j,k) = sinqx cosqy0(i,j,k) = cosqy sinqy0(i,j,k) = sinqy endif c1 = pqf(1)*cosqy + pqf(2)*sinqy c2 = pqf(3)*cosqy + pqf(4)*sinqy dc1 = -pqf(1)*sinqy + pqf(2)*cosqy dc2 = -pqf(3)*sinqy + pqf(4)*cosqy reqx = c1*cosqx - c2*sinqx imqx = c1*sinqx + c2*cosqx dreqx = -c1*sinqx - c2*cosqx dimqx = c1*cosqx -c2*sinqx c c Sum over different members of isospin multiplets c do imu=1,nmu expa = exp(a - b - mu(imu)) denom = (expa + sgn)/mult(imu) fact = numer/(denom*(1.+sgn/expa)) c P s(1,1) = s(1,1) + reqx*numer/denom s(1,2) = s(1,2) + imqx*numer/denom c T*dP/dT s1 = (a - b + dmu(imu,1))*fact s(2,1) = s(2,1) + reqx*s1 s(2,2) = s(2,2) + imqx*s1 c vT*dP/dvT s1 = (gamT**2*(b-vT2*rhop2*a)+dmu(imu,2))*fact s(3,1) = s(3,1) + reqx*s1 s(3,2) = s(3,2) + imqx*s1 c dP/dycm s1 = shep*(num1/denom - mT*gamT*fact/T) s(9,1) = s(9,1) + reqx*(s1 + dmu(imu,8)*fact) s(9,2) = s(9,2) + imqx*(s1 + dmu(imu,8)*fact) c et0*dP/det0 s2 = sqr*(q(4)*z-q(3)*time)*eta*numer/denom s(4,1) = s(4,1) + reqx*eta*s1 + dreqx*s2 s(4,2) = s(4,2) + imqx*eta*s1 + dimqx*s2 c dP/dmuB s(5,1) = s(5,1) + reqx*dmu(imu,4)*fact s(5,2) = s(5,2) + imqx*dmu(imu,4)*fact c tau*dP/dtau s1 = dmu(imu,5)*fact + (num1*chep-2.*aT*num2)/denom s2 = qx0*numer/denom s(6,1) = s(6,1) + reqx*s1 + dreqx*s2 s(6,2) = s(6,2) + imqx*s1 + dimqx*s2 c R*dP/dR s1 = dmu(imu,6)*fact + (2.*num1*chep - aT*num2)/denom s2 = -qx*numer/denom s3 = qy*numer/denom s(7,1) = s(7,1) + reqx*s1 + dreqx*s2 + + (dc1*cosqx - dc2*sinqx)*s3 s(7,2) = s(7,2) + imqx*s1 + dimqx*s2 + + (dc1*sinqx + dc2*cosqx)*s3 c dP/daT s1 = dmu(imu,7)*fact + + (.5*rhop2*num1*chep/sqr**2 - num2)/denom s2 = .5*rhop2*qx0*numer/(denom*sqr**2) s(8,1) = s(8,1) + reqx*s1 + dreqx*s2 s(8,2) = s(8,2) + imqx*s1 + dimqx*s2 c fs*dP/dfs s(10,1) = s(10,1) + reqx*dmu(imu,9)*fact s(10,2) = s(10,2) + imqx*dmu(imu,9)*fact enddo enddo enddo enddo fact = tau*R*R*et0*fac P(1) = fact*s(1,1) P(2) = fact*s(1,2) do ii=1,9 dP(ii,1) = fact*s(ii+1,1) dP(ii,2) = fact*s(ii+1,2) enddo return END SUBROUTINE corpim(P,dP,y,KT,dy,qout,qside,T,vT,et0,muB,tau,R,aT * ,lam,muS,muI,dmuS,dmuI,rw,fs,n1,n2) c c Calculates the negative pion correlation func and its derivatives c with respect to T, vT, et0, muB (=muB/T), tau, R, aT, ycm, fs, lam. c c For two-body resonance decays, Mran=0 and c a resonance mres decays into a pion mpi and another particle mX. c c All three-body decays considered have at least two pions in the c final state. So for these, mres decays into two pions and mX. c The other pion and mX can have an invariant mass anywhere in c mpi + mX < M < mres - mpi c So Mran = (mres - 2mpi - mX)/2 c Mave = (mres + mX)/2 c c sgn = -1. for bosons and +1. for fermions c nmu is the number of different baryon (nB), strangeness (nS) and c isospin (nI) states of mres. For each state, mult is the spin c multiplicity of mres multiplied by branching ratio of the decay. c c Integration over left-over invariant mass M in 3-particle decays c is done by n2(6)-point Gauss-Legendre quadratures c Integration over rapidity of the resonance yres is done by c n2(5)-point Gauss-Legendre quadratures c Integration over mTres is done by n2(4)-point Gauss-Chebyshev quads c with variable xk = (mTres-mT_ave)/del_mT and weight function c Wk = 1/sqrt(1-xk^2) c = sqrt(E2)*del_mT/sqrt((pT*pTres)^2-(E0*mres-mTres*E)^2) c Note: \int dmTres f(mTres)/sqrt((pT*pTres)^2-(E0*mres-mTres*E)^2) c = \int_-1^1 dxk/sqrt(1-xk^2) f(mTres(xk))/sqrt(E2) c When pT=0, del_mT=0, so mTres(xk) = mT_ave, and E2(mTres) = E2(mT_ave). c Consequently the xk integral can be done analytically: c \int_-1^1 dxk/sqrt(1-xk^2) = pi c IMPLICIT NONE INTEGER Nd,Ndecay,i,j,k,ii,iM,nj,nk,nka,imu,Ifirst,nM PARAMETER (Nd=15, nM=10, nj=20, nk=20) INTEGER nmu(Nd),nS(4,Nd),nB(4,Nd),niM(Nd),n1(6),n2(6) REAL mres(Nd),E0(Nd),p0(Nd),sgn(Nd),mult(4,Nd),mlt(4),gam(Nd) REAL mX(Nd),Mave(Nd),Mran(Nd),sa(10,2),sM(10,2),s(10,2) REAL xiM(nM),wiM(nM),xj(nj),wj(nj),xk(nk),wk,rw,Pr(2),dPr(9,2),P REAL dP(10),Pd(2),dPd(9,2),y,KT,mT,T,vT,et0,muB,tau,R,aT,muS,muI REAL nI(4,Nd),dmuS(9),dmuI(9),mu(4,Nd),dmu(4,Nd,9),mud(4),mpi,mpi2 REAL muB0,muS0,muI0,dmud(4,9),mres2,C,M,p0M,M2,E02,del_y,moff REAL yres,del_mT,E2p,pTres,mTres,pi2,KT2,E,mT_ave,fact,sqrtE2 REAL dy,qout,qside,lam,y1,y2,pT1,pT2,mT1,mT2,K0,K3,E1,E2,pz1,pz2 REAL YK,q(4),chyr,shyr,pTx,pTy,den,pq,pqy,P1,P2,dP1(10),dP2(10) REAL pqf(4),fs PARAMETER (pi2 = 6.283185307, mpi = 139.6, mpi2 = mpi**2) c SAVE Ifirst,muB0,muS0,muI0,mres,E0,p0,sgn,niM,nmu,mult,nI,nS,nB SAVE mX,Mave,Mran,xiM,wiM,xj,wj,xk,wk,Ndecay,mu,dmu,gam DATA Ifirst /1/ DATA muB0,muS0,muI0 /3*0./ DATA mres /547.5, 547.5, 781.9, 781.9, 770., 892., 1232., * 1385., 1385., 1407., * 958.,1440.,1440.,1525.,1525./ DATA mX /135., 0., 135., mpi, mpi, 494., 939., * 1116., 1194., 1194., * 549.,939.,939.,939.,939./ DATA Mave /341.25, 273.75, 458.45, 0., 0., 0., 0., 0., 0., 0., * 753.5,0.,1189.5,0.,1232./ DATA Mran /66.65, 134.15, 183.85, 0., 0., 0., 0., 0., 0., 0., * 65.5,0.,111.5,0.,154./ DATA sgn /-1.,-1.,-1.,-1.,-1.,-1.,1.,1.,1.,1., * -1.,1.,1.,1.,1./ DATA nmu /1,1,1,1,2,2,4,2,4,2,1,2,4,2,4/ DATA mult /.236,0.,0.,0., .049,0.,0.,0., 2.66,0.,0.,0., * .07,0.,0.,0., 3.,3.,0.,0., 2.,2.,0.,0., 3.98,1.33,3.98,1.33, * 3.52,3.52,0.,0., .24,.24,.24,.24, .67,.67,0.,0., * .652,0.,0.,0., .8,.8,0.,0., .416,.65,.416,.65, 2.,2.,0.,0., * .922,1.74,.922,1.74/ DATA nI /0.,0.,0.,0., 0.,0.,0.,0., 0.,0.,0.,0., 0.,0.,0.,0., * -1.,0.,0.,0., -.5,-.5,0.,0., -1.5,-.5,-1.5,-.5, * -1.,-1.,0.,0., -1.,0.,-1.,0., 0.,0.,0.,0., * 0.,0.,0.,0., -.5,-.5,0.,0., .5,-.5,.5,-.5, -.5,-.5,0.,0., * .5,-.5,.5,-.5/ DATA nS /0,0,0,0, 0,0,0,0, 0,0,0,0, 0,0,0,0, * 0,0,0,0, 1,-1,0,0, 0,0,0,0, * -1,1,0,0, -1,-1,1,1, -1,1,0,0, * 0,0,0,0, 0,0,0,0, 0,0,0,0, 0,0,0,0, 0,0,0,0/ DATA nB /0,0,0,0, 0,0,0,0, 0,0,0,0, 0,0,0,0, * 0,0,0,0, 0,0,0,0, 1,1,-1,-1, * 1,-1,0,0, 1,1,-1,-1, 1,-1,0,0, * 0,0,0,0, 1,-1,0,0, 1,1,-1,-1, 1,-1,0,0, 1,1,-1,-1/ DATA gam /.0012,.0012,8.4,8.4,151.,50.,120.,38.,38.,50., * .24,240.,240.,140.,140./ c c Find off-shell mass and momentum components for K c KT2 = KT*KT y1 = y + .5*dy y2 = y - .5*dy pT1 = sqrt((KT+.5*qout)**2 + .25*qside**2) pT2 = sqrt((KT-.5*qout)**2 + .25*qside**2) mT1 = sqrt(mpi2 + pT1*pT1) mT2 = sqrt(mpi2 + pT2*pT2) E1 = mT1*cosh(y1) E2 = mT2*cosh(y2) pz1 = mT1*sinh(y1) pz2 = mT2*sinh(y2) K0 = .5*(E1 + E2) K3 = .5*(pz1 + pz2) mT = sqrt(K0*K0 - K3*K3) YK = .5*log((K0+K3)/(K0-K3)) q(1) = qout q(2) = qside q(3) = pz1 - pz2 q(4) = E1 - E2 moff = K0*K0-K3*K3-KT*KT if (moff.le.0.) then write(*,*) 'WARNING: moff = ',moff write(100,*) 'WARNING: moff = ',moff P = 0. do ii=1,10 dP(ii) = 0. enddo return endif moff = sqrt(moff) c c Find decay energy and momenta with offshell mass c do i=1,Nd if (Mran(i).eq.0.) then E0(i) = .5*(mres(i)**2 + moff**2 - mX(i)**2)/mres(i) p0(i) = (mres(i)**2-(moff+mX(i))**2)* * (mres(i)**2-(moff-mX(i))**2) if (p0(i).lt.0.) then write(*,*) 'mres,moff,mX',mres(i),moff,mX(i) write(100,*) 'mres,moff,mX',mres(i),moff,mX(i) else p0(i) =.5*sqrt(p0(i))/mres(i) endif else Mran(i) = .5*(mres(i)-moff-mX(i)-mpi) Mave(i) = .5*(mres(i)-moff+mX(i)+mpi) if (Mran(i).lt.0.) then write(*,*) 'mres,moff,mX',mres(i),moff,mX(i) write(100,*) 'mres,moff,mX',mres(i),moff,mX(i) endif endif enddo c c Only call this part when the chemical potentials change c if (muB.ne.muB0.or.muS.ne.muS0.or.muI.ne.muI0.or.Ifirst.eq.1) then muB0 = muB muS0 = muS muI0 = muI c c Find abscissas and weights for integration (first time only) c if (Ifirst.eq.1) then Ifirst = 0 if (n2(4).gt.nk) then write(*,*) * 'too many points for mT-integration, n2(4) > nk in corpim' stop endif if (n2(5).gt.nj) then write(*,*) * 'too many points for y-integration, n2(5) > nj in corpim' stop endif if (abs(n2(6)).gt.nM) then write(*,*) * 'too many points for M-integration, n2(6) > nM in corpim' stop endif do i=1,n2(4) xk(i) = cos(.5*pi2*(float(i)-.5)/float(n2(4))) enddo wk = .5*pi2/float(n2(4)) call gauleg(-1.,1.,xj,wj,n2(5)) call gauleg(-1.,1.,xiM,wiM,abs(n2(6))) if (rw.lt.0.) then Ndecay = 0 else Ndecay = 10 endif c c Specify which resonances have 3-body decays c do i=1,Ndecay if (Mran(i).eq.0.) then niM(i) = 1 else niM(i) = abs(n2(6)) endif enddo endif c c Calculate the chemical potential and its derivatives c There is an implicit assumption of positive baryon number c do i=1,Ndecay do imu=1,nmu(i) mu(imu,i) = nI(imu,i)*muI + nB(imu,i)*muB + nS(imu,i)*muS do ii=1,9 dmu(imu,i,ii) = nI(imu,i)*dmuI(ii) + nS(imu,i)*dmuS(ii) enddo dmu(imu,i,4) = dmu(imu,i,4) +float(nB(imu,i))*sign(1.,muB) enddo enddo endif c c c Direct pi- c mud(1) = -muI mlt(1) = 1. do ii=1,9 dmud(1,ii) = -dmuI(ii) enddo pqf(1) = 1. pqf(2) = 0. pqf(3) = 0. pqf(4) = 0. CALL dst2(Pd,dPd,YK,KT,0.,mT,q,-1., * T,vT,et0,tau,R,aT,mud,dmud,1,mlt,n2,pqf) c c c Pi- from resonance decays c do ii=1,10 sa(ii,1) = 0. sa(ii,2) = 0. enddo do i=1,Ndecay if ((p0(i).ge.0.).and.(Mran(i).ge.0.)) then mres2 = mres(i)*mres(i) do imu=1,nmu(i) mud(imu) = mu(imu,i) mlt(imu) = mult(imu,i) do ii=1,9 dmud(imu,ii) = dmu(imu,i,ii) enddo enddo c c In three-body decays, integrate over mass^2 of the other 2 particles c C = 0. do ii=1,10 sM(ii,1) = 0. sM(ii,2) = 0. enddo do iM=1,niM(i) if (Mran(i).ne.0.) then M = Mave(i) + Mran(i)*xiM(iM) M2 = M*M p0(i) = * .5*sqrt((mres2-(moff+M)**2)*(mres2-(moff-M)**2)/mres2) p0M = sqrt((M2-(mpi+mX(i))**2)*(M2-(mpi-mX(i))**2))/M E0(i) = .5*(mres2 + moff*moff - M2)/mres(i) endif E02 = E0(i)*E0(i) del_y = log((sqrt(E02+KT2) + p0(i))/mT) c c Integrate over yres (rapidity of the resonance) c do ii=1,10 s(ii,1) = 0. s(ii,2) = 0. enddo do j=1,n2(5) yres = YK + del_y*xj(j) chyr = cosh(yres) shyr = sinh(yres) E = mT*cosh(YK-yres) E2p = E*E - KT2 mT_ave = mres(i)*E0(i)*E/E2p if (KT.eq.0.) then del_mT = 0. sqrtE2 = E nka = 1 else del_mT = mres(i)*KT*sqrt(E02 - E2p)/E2p sqrtE2 = sqrt(E2p) nka = n2(4) endif c c Integrate over mTres c do k=1,nka mTres = mT_ave + del_mT*xk(k) pTres = sqrt(mTres*mTres - mres2) if (KT.eq.0.) then pTx = 0. pTy = pTres else pTx = (mTres*E - E0(i)*mres(i))/KT pTy = sqrt(pTres**2-pTx**2) endif pq = ( mTres*(chyr*q(4)-shyr*q(3)) - pTx*q(1) ) * /(mres(i)*gam(i)) pqy = pTy*q(2)/(mres(i)*gam(i)) den = (1.-pq*pq+pqy*pqy)**2 + 4.*pq*pq pqf(1) = (1. + pq*pq + pqy*pqy)/den pqf(2) = -pqy*(1.-pq*pq+pqy*pqy)/den pqf(3) = pq*(1.+pq*pq-pqy*pqy)/den pqf(4) = -2.*pq*pqy/den CALL dst2(Pr,dPr,yres,pTx,pTy,mTres,q,sgn(i), * T,vT,et0,tau,R,aT,mud,dmud,nmu(i),mlt,n2,pqf) c P fact = wj(j)*mTres/sqrtE2 s(1,1) = s(1,1) + fact*Pr(1) s(1,2) = s(1,2) + fact*Pr(2) c dP do ii=2,10 s(ii,1) = s(ii,1) + fact*dPr(ii-1,1) s(ii,2) = s(ii,2) + fact*dPr(ii-1,2) enddo enddo enddo c c C is the normalization of P(M) for 3-body decays c if (Mran(i).ne.0.) then fact = wiM(iM)*p0M*del_y C = C + wiM(iM)*p0(i)*p0M else fact = del_y/p0(i) endif do ii=1,10 sM(ii,1) = sM(ii,1) + fact*s(ii,1) sM(ii,2) = sM(ii,2) + fact*s(ii,2) enddo enddo if (Mran(i).eq.0.) C=1. fact = mres(i)/C c c sa(10,i) = fs*dsa(1,i)/dfs (incomplete chemical equilibrium means fs<1) c if (nS(1,i).ne.0) then fact = fact*fs**abs(nS(1,i)) sa(10,1) = sa(10,1) + abs(nS(1,i))*fact*sM(1,1) sa(10,2) = sa(10,2) + abs(nS(1,i))*fact*sM(1,2) endif do ii=1,10 sa(ii,1) = sa(ii,1) + fact*sM(ii,1) sa(ii,2) = sa(ii,2) + fact*sM(ii,2) enddo endif enddo fact = wk/pi2 if (KT.eq.0.) fact = .5 Pd(1) = Pd(1) + fact*sa(1,1) Pd(2) = Pd(2) + fact*sa(1,2) do ii=1,9 dPd(ii,1) = dPd(ii,1) + fact*sa(ii+1,1) dPd(ii,2) = dPd(ii,2) + fact*sa(ii+1,2) enddo c c Give correct Jacobians for minimization variables c dPd(2,1) = dPd(2,1)*2.*(1.-vT) dPd(2,2) = dPd(2,2)*2.*(1.-vT) dPd(3,1) = dPd(3,1) + Pd(1) dPd(3,2) = dPd(3,2) + Pd(2) dPd(7,1) = dPd(7,1)*(1.-aT*aT) dPd(7,2) = dPd(7,2)*(1.-aT*aT) c c c Calculate one-pion distributions c CALL distpim(P1,dP1,y1,pT1,T,vT,et0,muB,tau,R,aT, * muS,muI,dmuS,dmuI,rw,fs,n1) CALL distpim(P2,dP2,y2,pT2,T,vT,et0,muB,tau,R,aT, * muS,muI,dmuS,dmuI,rw,fs,n1) c c Calculate correlation function and its derivatives c fact = lam/(P1*P2) P = 1. + fact*(Pd(1)**2 + Pd(2)**2) do ii=1,9 dP(ii) = 2.*fact*(Pd(1)*dPd(ii,1) + Pd(2)*dPd(ii,2)) + - (P-1.)*(dP1(ii)/P1 + dP2(ii)/P2) enddo c lam*dP/dlam dP(10) = P - 1. return END SUBROUTINE corkp(P,dP,y,KT,dy,qout,qside,T,vT,et0,muB,tau,R,aT * ,lam,muS,muI,dmuS,dmuI,rw,fs,n1,n2) c c Calculates the positive kaon correlation func and its derivatives c with respect to T, vT, et0, muB (=muB/T), tau, R, aT, ycm, fs, lamK. c c A resonance mres decays into a kaon mK and another particle mX. c sgn = -1. for bosons and +1. for fermions c nmu is the number of different baryon (nB), strangeness (nS) and c isospin (nI) states of mres. For each state, mult is the spin c multiplicity of mres multiplied by branching ratio of the decay. c c Integration over rapidity of the resonance yres is done by c n2(5)-point Gauss-Legendre quadratures c Integration over mTres is done by n2(4)-point Gauss-Chebyshev quads c with variable xk = (mTres-mT_ave)/del_mT and weight function c Wk = 1/sqrt(1-xk^2) c = sqrt(E2)*del_mT/sqrt((pT*pTres)^2-(E0*mres-mTres*E)^2) c Note: \int dmTres f(mTres)/sqrt((pT*pTres)^2-(E0*mres-mTres*E)^2) c = \int_-1^1 dxk/sqrt(1-xk^2) f(mTres(xk))/sqrt(E2) c When pT=0, del_mT=0, so mTres(xk) = mT_ave, and E2(mTres) = E2(mT_ave). c Consequently the xk integral can be done analytically: c \int_-1^1 dxk/sqrt(1-xk^2) = pi c IMPLICIT NONE INTEGER Nd,Ndecay,i,j,k,ii,nj,nk,nka,imu,Ifirst PARAMETER (Nd=1, nj=20, nk=20) INTEGER nmu(Nd),nS(2,Nd),n1(6),n2(6) REAL mres(Nd),E0(Nd),p0(Nd),sgn(Nd),mult(2,Nd),mlt(4),gam(Nd) REAL mX(Nd),sa(10,2),s(10,2) REAL xj(nj),wj(nj),xk(nk),wk,rw,mK,Pr(2),dPr(9,2),P,dP(10) REAL Pd(2),dPd(9,2),y,KT,mT,T,vT,et0,muB,tau,R,aT,muS,muI REAL nI(2,Nd),dmuS(9),dmuI(9),mu(2,Nd),dmu(2,Nd,9),mud(4) REAL muB0,muS0,muI0,dmud(4,9),mres2,E02,del_y,pqf(4),fs REAL yres,del_mT,E2p,pTres,mTres,pi2,KT2,E,mT_ave,fact,sqrtE2 REAL dy,qout,qside,lam,y1,y2,pT1,pT2,mT1,mT2,K0,K3,E1,E2,pz1,pz2 REAL YK,q(4),chyr,shyr,pTx,pTy,den,pq,pqy,P1,P2,dP1(10),dP2(10) PARAMETER (pi2 = 6.283185307) c SAVE Ifirst,muB0,muS0,muI0,mres,E0,p0,sgn,nmu,mult SAVE mX,nI,nS,xj,wj,xk,wk,Ndecay,mu,dmu,gam DATA Ifirst /1/ DATA muB0,muS0,muI0 /3*0./ DATA mres /892./ DATA sgn /-1./ DATA nmu /2/ DATA mult /1.,2./ DATA gam /51./ DATA mX /139.6/ DATA nI /.5,-.5/ DATA nS /1,1/ c c Calculate off-shell mass and momentum components of K c KT2 = KT*KT y1 = y + .5*dy y2 = y - .5*dy pT1 = sqrt((KT+.5*qout)**2 + .25*qside**2) pT2 = sqrt((KT-.5*qout)**2 + .25*qside**2) mT1 = sqrt(494.*494. + pT1*pT1) mT2 = sqrt(494.*494. + pT2*pT2) E1 = mT1*cosh(y1) E2 = mT2*cosh(y2) pz1 = mT1*sinh(y1) pz2 = mT2*sinh(y2) K0 = .5*(E1 + E2) K3 = .5*(pz1 + pz2) mT = sqrt(K0*K0 - K3*K3) YK = .5*log((K0+K3)/(K0-K3)) q(1) = qout q(2) = qside q(3) = pz1 - pz2 q(4) = E1 - E2 mK = K0*K0-K3*K3-KT*KT if (mK.le.0.) then write(*,*) 'WARNING: mK = ',mK write(100,*) 'WARNING: mK = ',mK P = 0. P = 0. do ii=1,10 dP(ii) = 0. enddo return endif mK = sqrt(mK) c c Find decay energy and momenta with offshell mass c do i=1,Nd E0(i) = .5*(mres(i)**2 + mK**2 - mX(i)**2)/mres(i) p0(i) = (mres(i)**2-(mK+mX(i))**2)* * (mres(i)**2-(mK-mX(i))**2) if (p0(i).lt.0.) then write(*,*) 'mres,mK,mX',mres(i),mK,mX(i) write(100,*) 'mres,mK,mX',mres(i),mK,mX(i) else p0(i) =.5*sqrt(p0(i))/mres(i) endif enddo c c Only call this part when the chemical potentials change c if (muB.ne.muB0.or.muS.ne.muS0.or.muI.ne.muI0.or.Ifirst.eq.1) then muB0 = muB muS0 = muS muI0 = muI if (Ifirst.eq.1) then Ifirst = 0 if (n2(4).gt.nk) then write(*,*) * 'too many points for mT-integration, n2(4) > nk in corpim' stop endif if (n2(5).gt.nj) then write(*,*) * 'too many points for y-integration, n2(5) > nj in corpim' stop endif do i=1,n2(4) xk(i) = cos(.5*pi2*(float(i)-.5)/float(n2(4))) enddo wk = .5*pi2/float(n2(4)) call gauleg(-1.,1.,xj,wj,n2(5)) if (rw.lt.0.) then Ndecay = 0 else Ndecay = Nd endif endif c c Calculate the chemical potential and its derivatives c do i=1,Ndecay do imu=1,nmu(i) mu(imu,i) = nI(imu,i)*muI + nS(imu,i)*muS do ii=1,9 dmu(imu,i,ii) = nI(imu,i)*dmuI(ii) + nS(imu,i)*dmuS(ii) enddo enddo enddo endif c c c Direct K+ c mud(1) = muS + .5*muI mlt(1) = 1. do ii=1,9 dmud(1,ii) = dmuS(ii) + .5*dmuI(ii) enddo pqf(1) = 1. pqf(2) = 0. pqf(3) = 0. pqf(4) = 0. CALL dst2(Pd,dPd,YK,KT,0.,mT,q,-1., * T,vT,et0,tau,R,aT,mud,dmud,1,mlt,n2,pqf) Pd(1) = fs*Pd(1) Pd(2) = fs*Pd(2) dPd(9,1) = dPd(9,1) + Pd(1) dPd(9,2) = dPd(9,2) + Pd(2) c c c K+ from resonance decays c do ii=1,10 sa(ii,1) = 0. sa(ii,2) = 0. enddo do i=1,Ndecay if (p0(i).ge.0.) then mres2 = mres(i)*mres(i) do imu=1,nmu(i) mud(imu) = mu(imu,i) mlt(imu) = mult(imu,i) do ii=1,9 dmud(imu,ii) = dmu(imu,i,ii) enddo enddo E02 = E0(i)*E0(i) del_y = log((sqrt(E02+KT2) + p0(i))/mT) c c Integrate over yres (rapidity of the resonance) c do ii=1,10 s(ii,1) = 0. s(ii,2) = 0. enddo do j=1,n2(5) yres = YK + del_y*xj(j) chyr = cosh(yres) shyr = sinh(yres) E = mT*cosh(YK-yres) E2p = E*E - KT2 mT_ave = mres(i)*E0(i)*E/E2p if (KT.eq.0.) then del_mT = 0. sqrtE2 = E nka = 1 else del_mT = mres(i)*KT*sqrt(E02 - E2p)/E2p sqrtE2 = sqrt(E2p) nka = n2(4) endif c c Integrate over mTres c do k=1,nka mTres = mT_ave + del_mT*xk(k) pTres = sqrt(mTres*mTres - mres2) if (KT.eq.0.) then pTx = 0. pTy = pTres else pTx = (mTres*E - E0(i)*mres(i))/KT pTy = sqrt(pTres**2-pTx**2) endif pq = ( mTres*(chyr*q(4)-shyr*q(3)) - pTx*q(1) ) * /(mres(i)*gam(i)) pqy = pTy*q(2)/(mres(i)*gam(i)) den = (1.-pq*pq+pqy*pqy)**2 + 4.*pq*pq pqf(1) = (1. + pq*pq + pqy*pqy)/den pqf(2) = -pqy*(1.-pq*pq+pqy*pqy)/den pqf(3) = pq*(1.+pq*pq-pqy*pqy)/den pqf(4) = -2.*pq*pqy/den CALL dst2(Pr,dPr,yres,pTx,pTy,mTres,q,sgn(i), * T,vT,et0,tau,R,aT,mud,dmud,nmu(i),mlt,n2,pqf) c P fact = wj(j)*mTres/sqrtE2 s(1,1) = s(1,1) + fact*Pr(1) s(1,2) = s(1,2) + fact*Pr(2) c dP do ii=2,10 s(ii,1) = s(ii,1) + fact*dPr(ii-1,1) s(ii,2) = s(ii,2) + fact*dPr(ii-1,2) enddo enddo enddo fact = mres(i)*del_y/p0(i) c c sa(10,i) = fs*dsa(1,i)/dfs (incomplete chemical equilibrium means fs<1) c if (nS(1,i).ne.0) then fact = fact*fs**abs(nS(1,i)) sa(10,1) = sa(10,1) + abs(nS(1,i))*fact*s(1,1) sa(10,2) = sa(10,2) + abs(nS(1,i))*fact*s(1,2) endif do ii=1,10 sa(ii,1) = sa(ii,1) + fact*s(ii,1) sa(ii,2) = sa(ii,2) + fact*s(ii,2) enddo endif enddo fact = wk/pi2 if (KT.eq.0.) fact = .5 Pd(1) = Pd(1) + fact*sa(1,1) Pd(2) = Pd(2) + fact*sa(1,2) do ii=1,9 dPd(ii,1) = dPd(ii,1) + fact*sa(ii+1,1) dPd(ii,2) = dPd(ii,2) + fact*sa(ii+1,2) enddo c c Give correct Jacobians for minimization variables c dPd(2,1) = dPd(2,1)*2.*(1.-vT) dPd(2,2) = dPd(2,2)*2.*(1.-vT) dPd(3,1) = dPd(3,1) + Pd(1) dPd(3,2) = dPd(3,2) + Pd(2) dPd(7,1) = dPd(7,1)*(1.-aT*aT) dPd(7,2) = dPd(7,2)*(1.-aT*aT) c dPd(9,1) = dPd(9,1) + Pd(1) c dPd(9,2) = dPd(9,2) + Pd(2) c c c Calculate one-kaon distributions c CALL distkp(P1,dP1,y1,pT1,T,vT,et0,muB,tau,R,aT, * muS,muI,dmuS,dmuI,rw,fs,n1) CALL distkp(P2,dP2,y2,pT2,T,vT,et0,muB,tau,R,aT, * muS,muI,dmuS,dmuI,rw,fs,n1) c c Calculate correlation function and its derivatives c fact = lam/(P1*P2) P = 1. + fact*(Pd(1)**2 + Pd(2)**2) do ii=1,9 dP(ii) = 2.*fact*(Pd(1)*dPd(ii,1) + Pd(2)*dPd(ii,2)) + - (P-1.)*(dP1(ii)/P1 + dP2(ii)/P2) enddo c dP(2) = dP(2)*2.*(1.-vT) c dP(7) = dP(7)*(1.-aT*aT) c lam*dP/dlam dP(10) = P - 1. return END SUBROUTINE convrt(mass,ycor,YK,KT,qo,qs,ql,dqy,opt) c c Convert qlong into dqy = y1-y2 c If opt = 3, ql = qlong in the YCMS c If opt = 2, ql = qlong in the frame defined by YK c IMPLICIT NONE INTEGER opt REAL mass,ycor,YK,KT,qo,qs,dqy,mT1,mT2,sh5dy,CHY,SHY,ql c mT1 = sqrt((KT+.5*qo)**2 + .25*qs**2 + mass**2) mT2 = sqrt((KT-.5*qo)**2 + .25*qs**2 + mass**2) if (abs(opt).eq.3) then sh5dy = ql/(mT1+mT2) else CHY = cosh(YK-ycor) SHY = sinh(YK-ycor) sh5dy = ( ql*(mT1+mT2)*CHY - SHY*(mT1-mT2)* * sqrt(ql**2 +(mT1+mT2)**2 +4.*mT1*mT2*SHY**2) ) * /((mT1+mT2)**2 +4.*mT1*mT2*SHY**2) endif dqy = 2.*log(sh5dy + sqrt(sh5dy**2 + 1.)) return END SUBROUTINE gauleg(x1,x2,x,w,n) INTEGER n REAL x1,x2,x(n),w(n) REAL EPS PARAMETER (EPS=3.e-14) INTEGER i,j,m REAL p1,p2,p3,pp,xl,xm,z,z1 m=(n+1)/2 xm=0.5*(x2+x1) xl=0.5*(x2-x1) do 12 i=1,m z=cos(3.141592654*(i-.25)/(n+.5)) 1 continue p1=1. p2=0. do 11 j=1,n p3=p2 p2=p1 p1=((2.*j-1.)*z*p2-(j-1.)*p3)/j 11 continue pp=n*(z*p1-p2)/(z*z-1.) z1=z z=z1-p1/pp if(abs(z-z1).gt.EPS)goto 1 x(i)=xm-xl*z x(n+1-i)=xm+xl*z w(i)=2.*xl/((1.-z*z)*pp*pp) w(n+1-i)=w(i) 12 continue return END SUBROUTINE gaujac(x,w,n,alf,bet) INTEGER n,MAXIT REAL alf,bet,w(n),x(n) REAL EPS PARAMETER (EPS=3.e-14,MAXIT=10) CU USES gammln INTEGER i,its,j REAL alfbet,an,bn,r1,r2,r3,gammln REAL a,b,c,p1,p2,p3,pp,temp,z,z1 do 13 i=1,n if(i.eq.1)then an=alf/n bn=bet/n r1=(1.+alf)*(2.78/(4.+n*n)+.768*an/n) r2=1.+1.48*an+.96*bn+.452*an*an+.83*an*bn z=1.-r1/r2 else if(i.eq.2)then r1=(4.1+alf)/((1.+alf)*(1.+.156*alf)) r2=1.+.06*(n-8.)*(1.+.12*alf)/n r3=1.+.012*bet*(1.+.25*abs(alf))/n z=z-(1.-z)*r1*r2*r3 else if(i.eq.3)then r1=(1.67+.28*alf)/(1.+.37*alf) r2=1.+.22*(n-8.)/n r3=1.+8.*bet/((6.28+bet)*n*n) z=z-(x(1)-z)*r1*r2*r3 else if(i.eq.n-1)then r1=(1.+.235*bet)/(.766+.119*bet) r2=1./(1.+.639*(n-4.)/(1.+.71*(n-4.))) r3=1./(1.+20.*alf/((7.5+alf)*n*n)) z=z+(z-x(n-3))*r1*r2*r3 else if(i.eq.n)then r1=(1.+.37*bet)/(1.67+.28*bet) r2=1./(1.+.22*(n-8.)/n) r3=1./(1.+8.*alf/((6.28+alf)*n*n)) z=z+(z-x(n-2))*r1*r2*r3 else z=3.*x(i-1)-3.*x(i-2)+x(i-3) endif alfbet=alf+bet do 12 its=1,MAXIT temp=2.+alfbet p1=(alf-bet+temp*z)/2. p2=1. do 11 j=2,n p3=p2 p2=p1 temp=2*j+alfbet a=2*j*(j+alfbet)*(temp-2.) b=(temp-1.)*(alf*alf-bet*bet+temp*(temp-2.)*z) c=2.*(j-1+alf)*(j-1+bet)*temp p1=(b*p2-c*p3)/a 11 continue pp=(n*(alf-bet-temp*z)*p1+2.*(n+alf)*(n+bet)*p2)/(temp* *(1.-z*z)) z1=z z=z1-p1/pp if(abs(z-z1).le.EPS)goto 1 12 continue write(*,*) 'WARNING: too many iterations in gaujac' write(100,*) 'WARNING: too many iterations in gaujac' 1 x(i)=z w(i)=exp(gammln(alf+n)+gammln(bet+n)-gammln(n+1.)-gammln(n+ *alfbet+1.))*temp*2.**alfbet/(pp*p2) 13 continue return END SUBROUTINE gaulag(x,w,n,alf) INTEGER n,MAXIT REAL alf,w(n),x(n) REAL EPS PARAMETER (EPS=3.e-14,MAXIT=10) CU USES gammln INTEGER i,its,j REAL ai,gammln REAL p1,p2,p3,pp,z,z1 do 13 i=1,n if(i.eq.1)then z=(1.+alf)*(3.+.92*alf)/(1.+2.4*n+1.8*alf) else if(i.eq.2)then z=z+(15.+6.25*alf)/(1.+.9*alf+2.5*n) else ai=i-2 z=z+((1.+2.55*ai)/(1.9*ai)+1.26*ai*alf/(1.+3.5*ai))* *(z-x(i-2))/(1.+.3*alf) endif do 12 its=1,MAXIT p1=1. p2=0. do 11 j=1,n p3=p2 p2=p1 p1=((2*j-1+alf-z)*p2-(j-1+alf)*p3)/j 11 continue pp=(n*p1-(n+alf)*p2)/z z1=z z=z1-p1/pp if(abs(z-z1).le.EPS)goto 1 12 continue write(*,*) 'WARNING: too many iterations in gaulag' write(100,*) 'WARNING: too many iterations in gaulag' 1 x(i)=z w(i)=-exp(gammln(alf+n)-gammln(float(n)))/(pp*n*p2) 13 continue return END FUNCTION gammln(xx) REAL gammln,xx INTEGER j REAL ser,stp,tmp,x,y,cof(6) SAVE cof,stp DATA cof,stp/76.18009172947146,-86.50532032941677, *24.01409824083091,-1.231739572450155,.1208650973866179e-2, *-.5395239384953e-5,2.5066282746310005/ x=xx y=x tmp=x+5.5 tmp=(x+0.5)*log(tmp)-tmp ser=1.000000000190015 do 11 j=1,6 y=y+1. ser=ser+cof(j)/y 11 continue gammln=tmp+log(stp*ser/x) return END