PROGRAM plotfrz c c For the model described in LA-UR-96-0782, also in http://xxx.lanl.gov c with designation nucl-th/9603007. c c Plots of Ed^3N/dp^3 (or d^2N/(dydpT)) for fixed y c fort.11 = pi- c fort.12 = pi+ c fort.13 = K+ c fort.14 = K- c fort.15 = protons c fort.16 = pbar c fort.20 = pi+/pi- c c Plots dN/dy as a function of y c fort.21 = pi- c fort.22 = pi+ c fort.23 = K+ c fort.24 = K- c fort.25 = protons c fort.26 = pbar c c Plots 1D correlations c side out long c pi- fort.31 fort.41 fort.51 c pi+ fort.32 fort.42 fort.52 c K+ fort.33 fort.43 fort.53 c K- fort.34 fort.44 fort.54 c c Plots out-long correlation c pi- fort.71 c pi+ fort.72 c K+ fort.73 c K- fort.74 c IMPLICIT NONE INTEGER npT,i,ny,mdist(6),mdndy(6),mcorr(4),mqplot(3),nq,nq2 INTEGER n1(6),n2(6),option(4),j,nfile,nq1 REAL T,vT,et0,y,pTmax,dpT,pT,f,dndy,dy,ymax,dqy REAL masspi,massp,massK,muB,muS,muI,tau,R,pi2,g,rw,dq,qc REAL YK,KT,qside,qout,qz,qmax,lam,P,dP(10),dmuS(9),dmuI(9) REAL mTm,aT,ycm,gevfac,normpi,normk,Ycor,qmin,lamK,dq2,qc2,fs REAL normp,ymin,mq2d,q1min,q1max,q2min,q2max PARAMETER (masspi = 139.6, massK = 494., massp = 939.) PARAMETER (pi2 =6.283185307) EXTERNAL distpim,distkp,distpro c c Read in parameters c OPEN (unit=10,file='plotfrz.in',status='old') read(10,*) T,vT,et0 read(10,*) muB,tau,R read(10,*) aT,ycm,fs read(10,*) lam,lamK,muS,muI read(10,*) normpi,normK,normp read(10,*) mdist read(10,*) y,pTmax,npT read(10,*) mdndy read(10,*) ymin,ymax,ny read(10,*) mcorr read(10,*) mqplot read(10,*) qmin,qmax,nq read(10,*) mq2D read(10,*) q1min,q1max,q2min,q2max,nq1,nq2 read(10,*) YK,KT,qside,qout,dqy read(10,*) ycor, rw read(10,*) option read(10,*) n1 read(10,*) n2 CLOSE(10) y = y-ycm muS = muS/T muI = muI/T tau = tau/197.33 R = R/197.33 qz = dqy if (option(2).eq.1) then gevfac = 1000. else gevfac = 1. endif c dpT = pTmax/float(npT) mTm = -dpT pT = -dpT c c c Plot dN/dydpT c do i=0,npT mTm = mTm + dpT pT = pT + dpT c pi- if (mdist(1).eq.1) then if (option(3).eq.0) pT = sqrt(mTm*(mTm + 2.*masspi)) CALL distpim(P,dP,y,pT,T,vT,et0,muB,tau,R,aT, * muS,muI,dmuS,dmuI,rw,fs,n1) g = P if (option(3).eq.0) then write(11,'(f10.3,e15.6)') mTm/gevfac,P*normpi*gevfac**2 else write(11,'(f10.3,e15.6)')pT/gevfac,pi2*pT*normpi*P*gevfac endif endif c pi+ if (mdist(2).eq.1) then if (option(3).eq.0) pT = sqrt(mTm*(mTm + 2.*masspi)) muB = -muB muS = -muS muI = -muI CALL distpim(P,dP,y,pT,T,vT,et0,muB,tau,R,aT, * muS,muI,dmuS,dmuI,rw,fs,n1) if (option(3).eq.0) then write(12,'(f10.3,e15.6)') mTm/gevfac,P*normpi*gevfac**2 else write(12,'(f10.3,e15.6)')pT/gevfac,pi2*pT*normpi*P*gevfac endif muB = -muB muS = -muS muI = -muI c pi+/pi- if (mdist(1).eq.1) then write(20,'(f10.3,e15.6)') pT/gevfac,P/g endif endif c K+ if (mdist(3).eq.1) then if (option(3).eq.0) pT = sqrt(mTm*(mTm + 2.*massK)) CALL distkp(P,dP,y,pT,T,vT,et0,muB,tau,R,aT, * muS,muI,dmuS,dmuI,rw,fs,n1) if (option(3).eq.0) then write(13,'(f10.3,e15.6)') mTm/gevfac,P*normk*gevfac**2 else write(13,'(f10.3,e15.6)') pT/gevfac,pi2*pT*normK*P*gevfac endif endif c K- if (mdist(4).eq.1) then if (option(3).eq.0) pT = sqrt(mTm*(mTm + 2.*massK)) muS = -muS muI = -muI CALL distkp(P,dP,y,pT,T,vT,et0,muB,tau,R,aT, * muS,muI,dmuS,dmuI,rw,fs,n1) if (option(3).eq.0) then write(14,'(f10.3,e15.6)') mTm/gevfac,P*normk*gevfac**2 else write(14,'(f10.3,e15.6)') pT/gevfac,pi2*pT*normK*P*gevfac endif muS = -muS muI = -muI endif c protons if (mdist(5).eq.1) then if (option(3).eq.0) pT = sqrt(mTm*(mTm + 2.*massp)) CALL distpro(P,dP,y,pT,T,vT,et0,muB,tau,R,aT, * muS,muI,dmuS,dmuI,rw,fs,n1) if (option(3).eq.0) then write(15,'(f10.3,e15.6)') mTm/gevfac,P*normp*gevfac**2 else write(15,'(f10.3,e15.6)') pT/gevfac,pi2*pT*normp*P*gevfac endif endif c antiprotons if (mdist(5).eq.1) then if (option(3).eq.0) pT = sqrt(mTm*(mTm + 2.*massp)) muB = -muB muS = -muS muI = -muI CALL distpro(P,dP,y,pT,T,vT,et0,muB,tau,R,aT, * muS,muI,dmuS,dmuI,rw,fs,n1) if (option(3).eq.0) then write(16,'(f10.3,e15.6)') mTm/gevfac,P*normp*gevfac**2 else write(16,'(f10.3,e15.6)') pT/gevfac,pi2*pT*normp*P*gevfac endif muB = -muB muS = -muS muI = -muI endif enddo c c c Plot dndy c dy = (ymax-ymin)/float(ny) y = ymin -dy -ycm do i=0,ny y = y + dy c pi- if (mdndy(1).eq.1) then f = pi2*dndy(y,T,vT,et0,muB,tau,R,aT,muS,muI,rw,fs,n1, * masspi,distpim) write(21,'(f10.3,e15.6)') y+ycm,f endif c pi+ if (mdndy(2).eq.1) then muB = -muB muS = -muS muI = -muI f = pi2*dndy(y,T,vT,et0,muB,tau,R,aT,muS,muI,rw,fs,n1, * masspi,distpim) write(22,'(f10.3,e15.6)') y+ycm,f muB = -muB muS = -muS muI = -muI endif c K+ if (mdndy(3).eq.1) then f = pi2*dndy(y,T,vT,et0,muB,tau,R,aT,muS,muI,rw,fs,n1, * massK,distkp) write(23,'(f10.3,e15.6)') y+ycm,f endif c K- if (mdndy(4).eq.1) then muS = -muS muI = -muI f = pi2*dndy(y,T,vT,et0,muB,tau,R,aT,muS,muI,rw,fs,n1, * massK,distkp) write(24,'(f10.3,e15.6)') y+ycm,f muS = -muS muI = -muI endif c protons if (mdndy(5).eq.1) then f = pi2*dndy(y,T,vT,et0,muB,tau,R,aT,muS,muI,rw,fs,n1, * massp,distpro) write(25,'(f10.3,e15.6)') y+ycm,f endif c pbars if (mdndy(6).eq.1) then muB = -muB muS = -muS muI = -muI f = pi2*dndy(y,T,vT,et0,muB,tau,R,aT,muS,muI,rw,fs,n1, * massp,distpro) write(26,'(f10.3,e15.6)') y+ycm,f muB = -muB muS = -muS muI = -muI endif enddo c c c Plot 1-D projection of HBT correlation function c dq = (qmax-qmin)/float(nq) qc = qmin -dq do i=0,nq qc = qc + dq if (qc.eq.0.) qc = .1 c pi- if (mcorr(1).eq.1) then c side if (mqplot(1).eq.1) then if (option(4).ge.2) * call convrt(masspi,Ycor,YK,KT,qout,qc,qz,dqy,option(4)) CALL corpim(P,dP,YK-ycm,KT,dqy,qout,qc,T,vT,et0,muB, * tau,R,aT,lam,muS,muI,dmuS,dmuI,rw,fs,n1,n2) write(*,*) qc,P write(31,'(f10.3,e15.6)') qc,P endif c out if (mqplot(2).eq.1) then if (option(4).ge.2) * call convrt(masspi,Ycor,YK,KT,qc,qside,qz,dqy,option(4)) CALL corpim(P,dP,YK-ycm,KT,dqy,qc,qside,T,vT,et0,muB, * tau,R,aT,lam,muS,muI,dmuS,dmuI,rw,fs,n1,n2) write(*,*) qc,P write(41,'(f10.3,e15.6)') qc,P endif c long if (mqplot(3).eq.1) then if (option(4).ge.2) then call convrt(masspi,Ycor,YK,KT, * qout,qside,qc,dqy,option(4)) else dqy = qc endif CALL corpim(P,dP,YK-ycm,KT,dqy,qout,qside,T,vT,et0,muB, * tau,R,aT,lam,muS,muI,dmuS,dmuI,rw,fs,n1,n2) write(*,*) qc,P write(51,'(f10.3,e15.6)') qc,P endif endif c pi+ if (mcorr(2).eq.1) then muI = -muI muS = -muS muB = -muB c side if (mqplot(1).eq.1) then if (option(4).ge.2) * call convrt(masspi,Ycor,YK,KT,qout,qc,qz,dqy,option(4)) CALL corpim(P,dP,YK-ycm,KT,dqy,qout,qc,T,vT,et0,muB, * tau,R,aT,lam,muS,muI,dmuS,dmuI,rw,fs,n1,n2) write(*,*) qc,P write(32,'(f10.3,e15.6)') qc,P endif c out if (mqplot(2).eq.1) then if (option(4).ge.2) * call convrt(masspi,Ycor,YK,KT,qout,qc,qz,dqy,option(4)) CALL corpim(P,dP,YK-ycm,KT,dqy,qc,qside,T,vT,et0,muB, * tau,R,aT,lam,muS,muI,dmuS,dmuI,rw,fs,n1,n2) write(*,*) qc,P write(42,'(f10.3,e15.6)') qc,P endif c long if (mqplot(3).eq.1) then if (option(4).ge.2) then call convrt(masspi,Ycor,YK,KT, * qout,qside,qc,dqy,option(4)) else dqy = qc endif CALL corpim(P,dP,YK-ycm,KT,dqy,qout,qside,T,vT,et0,muB, * tau,R,aT,lam,muS,muI,dmuS,dmuI,rw,fs,n1,n2) write(*,*) qc,P write(52,'(f10.3,e15.6)') qc,P endif muI = -muI muS = -muS muB = -muB endif c K+ if (mcorr(3).eq.1) then c side if (mqplot(1).eq.1) then if (option(4).ge.2) * call convrt(massK,Ycor,YK,KT,qout,qc,qz,dqy,option(4)) CALL corkp(P,dP,YK-ycm,KT,dqy,qout,qc,T,vT,et0,muB, * tau,R,aT,lamK,muS,muI,dmuS,dmuI,rw,fs,n1,n2) write(*,*) qc,P write(33,'(f10.3,e15.6)') qc,P endif c out if (mqplot(2).eq.1) then if (option(4).ge.2) * call convrt(massK,Ycor,YK,KT,qc,qside,qz,dqy,option(4)) CALL corkp(P,dP,YK-ycm,KT,dqy,qc,qside,T,vT,et0,muB, * tau,R,aT,lamK,muS,muI,dmuS,dmuI,rw,fs,n1,n2) write(*,*) qc,P write(43,'(f10.3,e15.6)') qc,P endif c long if (mqplot(3).eq.1) then if (option(4).ge.2) then call convrt(massK,Ycor,YK,KT, * qout,qside,qc,dqy,option(4)) else dqy = qc endif CALL corkp(P,dP,YK-ycm,KT,dqy,qout,qside,T,vT,et0,muB, * tau,R,aT,lam,muS,muI,dmuS,dmuI,rw,fs,n1,n2) write(*,*) qc,P write(53,'(f10.3,e15.6)') qc,P endif endif c K- if (mcorr(4).eq.1) then muI = -muI muS = -muS muB = -muB c side if (mqplot(1).eq.1) then if (option(4).ge.2) * call convrt(massK,Ycor,YK,KT,qout,qc,qz,dqy,option(4)) CALL corkp(P,dP,YK-ycm,KT,dqy,qout,qc,T,vT,et0,muB, * tau,R,aT,lamK,muS,muI,dmuS,dmuI,rw,fs,n1,n2) write(*,*) qc,P write(34,'(f10.3,e15.6)') qc,P endif c out if (mqplot(2).eq.1) then if (option(4).ge.2) * call convrt(massK,Ycor,YK,KT,qc,qside,qz,dqy,option(4)) CALL corkp(P,dP,YK-ycm,KT,dqy,qc,qside,T,vT,et0,muB, * tau,R,aT,lamK,muS,muI,dmuS,dmuI,rw,fs,n1,n2) write(*,*) qc,P write(44,'(f10.3,e15.6)') qc,P endif c long if (mqplot(3).eq.1) then if (option(4).ge.2) then call convrt(massK,Ycor,YK,KT, * qout,qside,qc,dqy,option(4)) else dqy = qc endif CALL corkp(P,dP,YK-ycm,KT,dqy,qout,qside,T,vT,et0,muB, * tau,R,aT,lam,muS,muI,dmuS,dmuI,rw,fs,n1,n2) write(*,*) qc,P write(54,'(f10.3,e15.6)') qc,P endif muI = -muI muS = -muS muB = -muB endif if ((qc.eq..1).and.(dq.ne..1)) qc = 0. enddo c c c Plot 2-D out-long projection of HBT correlation function c if (mq2D.eq.1) then dq = (q1max-q1min)/float(nq1) dq2 = (q2max-q2min)/float(nq2) qc = q1min -dq do i=0,nq1 nfile = 200+i qc = qc + dq qc2 = -dq2 do j=0,nq2 qc2 = qc2 + dq2 c pi- if (mcorr(1).eq.1) then call convrt(masspi,Ycor,YK,KT,qc2,qside,qc,dqy,option(4)) CALL corpim(P,dP,YK-ycm,KT,dqy,qc2,qside,T,vT,et0,muB, * tau,R,aT,lam,muS,muI,dmuS,dmuI,rw,fs,n1,n2) write(*,*) qc,qc2,P write(71,'(f10.3,e15.6)') qc,qc2,P write(nfile,*) P endif c pi+ if (mcorr(2).eq.1) then muI = -muI muS = -muS muB = -muB call convrt(masspi,Ycor,YK,KT,qc2,qside,qc,dqy,option(4)) CALL corpim(P,dP,YK-ycm,KT,dqy,qc2,qside,T,vT,et0,muB, * tau,R,aT,lam,muS,muI,dmuS,dmuI,rw,fs,n1,n2) write(*,*) qc,qc2,P write(72,'(f10.3,e15.6)') qc,qc2,P write(nfile,*) P muI = -muI muS = -muS muB = -muB endif c K+ if (mcorr(3).eq.1) then call convrt(massK,Ycor,YK,KT,qc2,qside,qc,dqy,option(4)) CALL corkp(P,dP,YK-ycm,KT,dqy,qc2,qside,T,vT,et0,muB, * tau,R,aT,lamK,muS,muI,dmuS,dmuI,rw,fs,n1,n2) write(*,*) qc,qc2,P write(73,'(f10.3,e15.6)') qc,qc2,P write(nfile,*) P endif c K- if (mcorr(4).eq.1) then muI = -muI muS = -muS muB = -muB call convrt(massK,Ycor,YK,KT,qc2,qside,qc,dqy,option(4)) CALL corkp(P,dP,YK-ycm,KT,dqy,qc2,qside,T,vT,et0,muB, * tau,R,aT,lamK,muS,muI,dmuS,dmuI,rw,fs,n1,n2) write(*,*) qc,qc2,P write(74,'(f10.3,e15.6)') qc,qc2,P write(nfile,*) P muI = -muI muS = -muS muB = -muB endif enddo if (mcorr(1).eq.1) write(71,*) if (mcorr(2).eq.1) write(72,*) if (mcorr(3).eq.1) write(73,*) if (mcorr(4).eq.1) write(74,*) enddo endif END FUNCTION dndy(y,T,vT,et0,muB,tau,R,aT,muS,muI,rw,fs,n1,mass,func) c c Calculates dndy by integrating the distribution func over mTdmT. c The numerical mT integration is done using nl-point Gauss-Laguerre c quadrature with the integration variable w = dxj*(mT-mass) c The coefficient dxj is found by determining the slope of the log c of the distribution between mT-mass = 1500. and 2000. c IMPLICIT NONE INTEGER j,n1(6),nl,Ifirst PARAMETER (nl=12) REAL dndy,y,T,vT,et0,muB,tau,R,aT,muS,muI,fs REAL rw,mass,xl(nl),wl(nl),pT,dxj,mT,P,dP(10),dmuS(9),dmuI(9) REAL func DATA dmuS,dmuI /18*0./ DATA Ifirst /1/ SAVE xl,wl,dmuS,dmuI,Ifirst c c Find abcissas and weights for integration (first call only) c if (Ifirst.eq.1) then Ifirst = 0 call gaulag(xl,wl,nl,0.) endif c c Integrate over xl(j) = dxj*(mT-mass) c dndy = 0. do j=-1,nl c c j=-1 and j=0 determine the log of distribution at mT-m=500 and 1000. c From these values, the slope dxj is determined. c if (j.eq.-1) then mT = 500. + mass elseif (j.eq.0) then mT = 1000. + mass else mT = xl(j)/dxj + mass endif pT = sqrt(mT**2-mass**2) CALL func(P,dP,y,pT,T,vT,et0,muB,tau,R,aT, * muS,muI,dmuS,dmuI,rw,fs,n1) if (j.eq.-1) then dxj = log(mT*P) elseif (j.eq.0) then dxj = (dxj-log(mT*P))/500. else dndy = dndy + wl(j)*exp(xl(j))*mT*P endif enddo dndy = dndy/dxj return END