c c 13may2004: slight modification to allow quantum (as opposed to Boltzmann) statistics c c change since 16may2003: added v4 calculation c c this code written late December 2002. I want to get this as neat as c possible, and calculate as much stuff "in one go" as I can. c It does the integrals in the paper draft of Fabrice and me, calculating c o spectra[pi/K/p] c o v2(pT)[pi/K/p] c o spatial-correlation-tensor elements (e.g. <\tilde{x}^2>) and HBT(pT,phi_p)[pi only] c o spacetime separations b/t pi/K/p (to be implemented) c c It makes two output files per particle type c o PID-vs-phi.BW, containing tensor elements, average x,y,t, HBT radii, and N vs phi c o PID-vs-pT.BW containing N, v2, and FourierCoefficients of HBT radii vs pT c c program one_case c c point of this program is to plot R(phi) or R(mt) for the transverse c HBT radii, using the routines in the file chi2_ellipse.for c implicit none double precision pT,phi_p,mT integer npar,ipar,ipid parameter (npar=8) double precision xval dimension xval(npar) character*10 name(npar) data (name(ipar),ipar=1,npar)/'Ry_max','Rx_max','T','rho0_max', + 'rhoa_max','tau','a_skin','del-tau'/ double precision pi/3.1415927/ double precision mass_pion,mass_proton,mass_kaon,pT_Randy,mass parameter(mass_pion=0.138d0,mass_kaon=0.494d0,mass_proton=0.938d0) double precision pTmin,pTmax,delpT double precision phi_p_value,phi_p_max,phi_p_min,phi_p_step integer just_one_phi_p c these next variables are just the FourierCoefficients of the oscillating radii double precision Ro2_FC(3),Rs2_FC(3),Ros2_FC(3),Rl2_FC(3) integer index,iorder,count double precision fsum_sum,v2sum,v2 double precision v4sum,v4,v6,v6sum c these variables are the calculated quantities coming from BWcalc double precision fsum,xtilde2_ave,ytilde2_ave,xtilde_ytilde_ave, + xtilde_ttilde_ave,ytilde_ttilde_ave,ttilde2_ave,ztilde2_ave, + xave,yave,tave,Rout2,Rside2,Ros2,Rlong2 common/BW_answers/fsum,xtilde2_ave,ytilde2_ave, + xtilde_ytilde_ave,xtilde_ttilde_ave,ytilde_ttilde_ave, + ttilde2_ave,ztilde2_ave,xave,yave,tave, + Rout2,Rside2,Ros2,Rlong2 double precision AveCos2phi_b,Norm_phi_b ! added 19 dec 03 - Voloshin wanted it calculated common/CalculateAveCos2phib/AveCos2phi_b,Norm_phi_b integer*4 maxterm,iqstat,boson,fermion parameter(boson=1,fermion=-1) write(6,*)'--------------------------------------------------' write(6,*)'This is the all-new truly-boost-invariant BlastWave' write(6,*)'--------------------------------------------------' do ipar=1,npar write(6,*)'Enter value for ',name(ipar) read(5,*)xval(ipar) enddo c 13may2004 mal c write(6,*)'I can now treat protons as fermions and', + 'pions/kaons as bosons' write(6,*)'To what order in the quantum-statistics expansion' write(6,*)' shall we go? (n=1 is classical/Boltzman stats)' read(5,*)maxterm open(unit=19,file='define_title.kumac',status='unknown') write(19,177)char(39),(xval(ipar),ipar=1,8),char(39) close(unit=19) 177 format(1x,'title ',a1,'Ry=',f5.2,' Rx=',f5.2, + ' T=',f5.2,' rho0=',f4.2,' rhoa=',f5.3, + ' tau0=',f4.2,' askin=',f5.3,' Dtau=',f5.3,a1) c--- write(6,*)'Do calculation for all phi_p (0)' write(6,*)' or only one angle (1)?' write(6,*)' OR skip the BW numbers and make an x-y plot (2)?' write(6,*)' OR just output versus mass (only one phi!) (3)?' read(5,*)just_one_phi_p if (just_one_phi_p.eq.3) then call calculate_avept_versus_mass(name,xval,npar) stop elseif (just_one_phi_p.eq.2) then write(6,*)'Input pT and phi_p' read(5,*)pT,phi_p write(6,*)'Which particle? 1=pion 2=kaon 3=proton' read(5,*)ipid if (ipid.eq.1) then mass = mass_pion iqstat = boson elseif (ipid.eq.2) then mass = mass_kaon iqstat = boson elseif (ipid.eq.3) then mass = mass_proton iqstat = fermion else stop 'Not a valid choice - quitting' endif call make_xy_plot(xval,pT,phi_p,mass,iqstat,maxterm) stop elseif (just_one_phi_p.eq.1) then write(6,*)'Enter value of phi_p' read(5,*)phi_p_value phi_p_min = phi_p_value phi_p_max = phi_p_value phi_p_step = 10.0d0 else phi_p_min = 0.0d0 phi_p_max = 2.0d0*pi-pi/8.0d0 phi_p_step = pi/8.0d0 endif write(6,*)'Enter min/max pT, and pT step' read(5,*)pTmin,pTmax,delpT c I don't use the kaons in the BW paper anyway, so I just take them totally out c as I do a huge batch of BW jobs - 27may2004 c do ipid=1,3 write(6,*)'NOTE: Kaon calculations disabled (to save time)' do ipid=1,3,2 write(6,*)'Doing ipid = ',ipid if (ipid.eq.1) then mass = mass_pion iqstat = boson elseif (ipid.eq.2) then mass = mass_kaon iqstat = boson elseif (ipid.eq.3) then mass = mass_proton iqstat = fermion endif call open_files(ipid,name,xval,npar,just_one_phi_p,phi_p_value, + maxterm) do pT=pTmin,pTmax,delpT mT = dsqrt(mass*mass + pT*pT) do index=1,3 Ro2_FC(index) = 0.0d0 Rs2_FC(index) = 0.0d0 Ros2_FC(index) = 0.0d0 Rl2_FC(index) = 0.0d0 enddo count=0 fsum_sum = 0.0d0 v2sum = 0.0d0 v4sum = 0.0d0 v6sum = 0.0d0 * * IMPT: MUST zero these quantities at THIS point (so that we get integrated over \phi_p AveCos2phi_b = 0.0d0 ! 19 Dec 03 - Sergei V. asked me to check this... Norm_phi_b = 0.0d0 ! this quantity is kept so that we calculate the *average* * do phi_p=phi_p_min,phi_p_max,phi_p_step count = count+1 call BWcalc(xval,pT,phi_p,mass,iqstat,maxterm) ! calculated values passed back in commonblock do index=1,3 iorder = (index-1)*2 Ro2_FC(index) = Ro2_FC(index) + + Rout2*dcos(dble(iorder)*phi_p) Rs2_FC(index) = Rs2_FC(index) + + Rside2*dcos(dble(iorder)*phi_p) Ros2_FC(index) = Ros2_FC(index) + + Ros2*dsin(dble(iorder)*phi_p) Rl2_FC(index) = Rl2_FC(index) + + Rlong2*dcos(dble(iorder)*phi_p) enddo fsum_sum = fsum_sum + fsum ! fsum is the {1}_{0,0} integral in the paper v2sum = v2sum + fsum*dcos(2.0d0*phi_p) v4sum = v4sum + fsum*dcos(4.0d0*phi_p) v6sum = v6sum + fsum*dcos(6.0d0*phi_p) write(22,922) + pT,phi_p,mT*fsum, + xtilde2_ave,ytilde2_ave,xtilde_ytilde_ave, + xtilde_ttilde_ave,ytilde_ttilde_ave,ttilde2_ave, + ttilde2_ave,xave,yave,tave, + Rout2,Rside2,Ros2,Rlong2 enddo ! phi_p do index=1,3 Ro2_FC(index) = Ro2_FC(index)/dble(count) Rs2_FC(index) = Rs2_FC(index)/dble(count) Ros2_FC(index) = Ros2_FC(index)/dble(count) Rl2_FC(index) = Rl2_FC(index)/dble(count) enddo v2 = v2sum/fsum_sum v4 = v4sum/fsum_sum v6 = v6sum/fsum_sum write(6,*)'pT and :',pT, + AveCos2phi_b/Norm_phi_b write(21,921) + pT,mT*fsum_sum,v2, + (Ro2_FC(index),index=1,3), + (Rs2_FC(index),index=1,3), + (Ros2_FC(index),index=1,3), + (Rl2_FC(index),index=1,3),v4,v6, ! v4 added 19oct2003 + AveCos2phi_b/Norm_phi_b enddo ! pT enddo ! ipid 921 format(18(e9.3,1x)) 922 format(17(e9.3,1x)) end c----------------------------------------------------------------------------------- c----------------------------------------------------------------------------------- subroutine BWcalc(xval,pT,phi_p,mass,iqstat,maxterm) implicit none double precision xval(*) * * index parameter * ----- --------- * 1 Ry * 2 Rx * 3 T * 4 rho_0 * 5 rho_a * 6 tau (peak freezeout time in propertime) * 7 a_skin (*fractional* skin depth) * 8 del-tau (duration of freezeout in proper time) integer*4 iqstat ! this will be +1 if particle is boson and -1 if it is fermion integer*4 maxterm ! this is the maximum term we should calculate before truncating ! the series expansion of the quantum statistics ! if maxterm=1, then we only use the first term, which is the ! same as Boltzmann (classical) statistics integer*4 nterm double precision beta_N,signfactor double precision Ry,Rx,T,rho_0,rho_a,tau,a_skin,del_tau double precision Rout2,Rside2,Ros2,Rlong2 double precision phi_p,rho,phi_b,dphi,spatial_density integer iphi double precision fsum,xsum,ysum,x2sum,y2sum,xysum,phi_s,r,x,y double precision tsum,xtsum,ytsum,t2sum,z2sum double precision fvalue,fvalue_r,xave,yave,x2ave,y2ave,xyave,f double precision xtilde2_ave,ytilde2_ave,xtilde_ytilde_ave double precision tave,t2ave,z2ave,xtave,ytave double precision ttilde2_ave,ztilde2_ave,xtilde_ttilde_ave,ytilde_ttilde_ave double precision r_step,phi_step,pi,Rmax parameter(pi=3.1415927d0) double precision cphi,sphi,c2phi,s2phi double precision pT,mT2,pT2,mT double precision mass,mass2,beta_t2,beta_t double precision G00,G01,G02,G20,K0,K1,K2 double precision bessk0,bessk1,bessk double precision factor double precision alpha,beta,r_ellipse common/BW_answers/fsum,xtilde2_ave,ytilde2_ave, + xtilde_ytilde_ave,xtilde_ttilde_ave,ytilde_ttilde_ave, + ttilde2_ave,ztilde2_ave,xave,yave,tave, + Rout2,Rside2,Ros2,Rlong2 double precision eta2 ! (Ry/Rx)**2 double precision AveCos2phi_b,Norm_phi_b ! added 19 dec 03 - Voloshin wanted it calculated common/CalculateAveCos2phib/AveCos2phi_b,Norm_phi_b mass2=mass**2 Ry = dabs(xval(1)) Rx = dabs(xval(2)) T = dabs(xval(3)) rho_0 = dabs(xval(4)) rho_a = xval(5) tau = xval(6) a_skin = xval(7) del_tau = xval(8) pT2 = pT*pT mT2 = mass2+pT2 mT = dsqrt(mT2) beta_t2 = pT2/mT2 beta_t = dsqrt(beta_t2) if (Ry.gt.Rx) then Rmax = Ry else Rmax = Rx endif Rmax = Rmax*(1.0+8.0d0*a_skin) r_step = Rmax/200.0 phi_step = pi/dfloat(256) ! dividing by power of 2 ensures that ! 90 degrees is hit as well as zero deg. fsum = 0.0 xsum = 0.0 ysum = 0.0 x2sum = 0.0 y2sum = 0.0 xysum = 0.0 tsum = 0.0 xtsum = 0.0 ytsum = 0.0 t2sum = 0.0 z2sum = 0.0 ***no more s2 !!! eta = ((1.0d0+2.0d0*s2)/(1.0d0-2.0d0*s2))**(1.0d0/3.0d0) eta2 = (Ry/Rx)**2 c this outer loop over "nterm" is for the quantum statistics - mal 13may2004 do nterm=1,maxterm signfactor = dfloat(iqstat**(nterm+1)) do phi_s = 0.0,2.0d0*pi-phi_step/2.0,phi_step ! integral over phi_s do r = 0.0,Rmax,r_step ! note x = r*dcos(phi_s) y = r*dsin(phi_s) c===== the "f" calculated here in this sub-block is what we see in the definition of c===== the curly brackets in the paper: {B'}_kj. It is everything except for the B' and the G_jk c===== we also figure out here what is rho, and thus alpha and beta c application of the Theta function: r_ellipse = dsqrt(y*y+eta2*x*x)/Ry spatial_density = + 1.0d0/(1.0d0+dexp((r_ellipse-1.0d0)/a_skin)) ! Woods-Saxon phi_b = datan(dtan(phi_s)/eta2) c c of course the problem with atan() is that it is multiple-valued... c how to get the "right" answer (and it DOES matter) ? c Well, we know that phi_b and phi_s cannot be VERY much different, c so if they begin to differ by pi or more, then we are at the c wrong answer 1 dphi = phi_b - phi_s if (dabs(dphi).gt.3.1415927d0/2.0d0) then phi_b = phi_b - 3.1415927d0 * dphi/dabs(dphi) goto 1 endif rho = (rho_0 + rho_a*dcos(2.0d0*phi_b))*r_ellipse alpha = pT*sinh(rho)/T beta = mT*cosh(rho)/T fvalue = dexp(dfloat(nterm)*alpha*dcos(phi_b-phi_p))* ! note the "nterm" - mal 13may2004 + spatial_density c===== c===== end of sub-block c===== beta_N = beta * dfloat(nterm) ! this is now the argument used below - mal 13may2004 K0 = bessk0(beta_N) K1 = bessk1(beta_N) * K2 = bessk(2,beta_N) K2 = (2.0d0/beta_N)*K1 + K0 ! Gradshteyn&Ryzhik 8.486(17) - checked numerically 26dec2002 mal G00 = 2.0d0*K1 G01 = 2.0d0*(K1/beta_N + K0) G02 = 2.0d0*(K2/beta_N + K1) G20 = 2.0d0*K2/beta_N * fvalue_r = fvalue*r ! because integral is blah*r*dr*dphi AveCos2phi_b = AveCos2phi_b + + dcos(2.0d0*phi_b)*fvalue_r*G00 ! 19 Dec 03 Norm_phi_b = Norm_phi_b + fvalue_r*G00 * ABOVE IS RIGHT / BELOW IS WRONG - I'M JUST PLAYING... * AveCos2phi_b = AveCos2phi_b + * + dcos(2.0d0*phi_b)*fvalue_r ! 19 Dec 03 * Norm_phi_b = Norm_phi_b + fvalue_r c after all the loops, these will be the {B'}_{i,j} factors, which are just c integrals if we use classical (Boltzmann) statistics, and are an infinite c (in principle) sum of integrals in quantum statisitcs c e.g. fsum = {1}_{0,0} fsum = fsum + fvalue_r * G00 * signfactor xsum = xsum + fvalue_r * x * G00 * signfactor ysum = ysum + fvalue_r * y * G00 * signfactor x2sum = x2sum + fvalue_r * x*x * G00 * signfactor y2sum = y2sum + fvalue_r * y*y * G00 * signfactor xysum = xysum + fvalue_r * x*y * G00 * signfactor * tsum = tsum + fvalue_r*G01 * signfactor xtsum = xtsum + fvalue_r*G01*x * signfactor ytsum = ytsum + fvalue_r*G01*y * signfactor * t2sum = t2sum + fvalue_r*G02 * signfactor * z2sum = z2sum + fvalue_r*G20 * signfactor enddo ! end integral over r enddo ! end integral over phi_s enddo ! end integral over nterm * factor = (del_tau**2 + tau**2)/tau tsum = tsum * factor xtsum = xtsum * factor ytsum = ytsum * factor * factor = 3.0d0*del_tau**2 + tau**2 t2sum = t2sum * factor z2sum = z2sum * factor * xave = xsum/fsum yave = ysum/fsum tave = tsum/fsum x2ave = x2sum/fsum y2ave = y2sum/fsum t2ave = t2sum/fsum z2ave = z2sum/fsum xyave = xysum/fsum xtave = xtsum/fsum ytave = ytsum/fsum xtilde2_ave = x2ave-xave*xave ytilde2_ave = y2ave-yave*yave ttilde2_ave = t2ave-tave*tave ztilde2_ave = z2ave - 0.0d0 xtilde_ytilde_ave = xyave - xave*yave xtilde_ttilde_ave = xtave - xave*tave ytilde_ttilde_ave = ytave - yave*tave cphi = dcos(phi_p) sphi = dsin(phi_p) c2phi = dcos(2.0d0*phi_p) s2phi = dsin(2.0d0*phi_p) Rout2 = 0.5d0*(xtilde2_ave+ytilde2_ave) 1 + 0.5d0*(xtilde2_ave-ytilde2_ave)*c2phi 2 + xtilde_ytilde_ave*s2phi 3 - 2.0d0*beta_t*(xtilde_ttilde_ave*cphi + ytilde_ttilde_ave*sphi) 4 + beta_t2*ttilde2_ave Rside2 = 0.5d0*(xtilde2_ave+ytilde2_ave) 1 - 0.5d0*(xtilde2_ave-ytilde2_ave)*c2phi 2 - xtilde_ytilde_ave*s2phi Ros2 = xtilde_ytilde_ave*c2phi 1 + 0.5d0*(ytilde2_ave-xtilde2_ave)*s2phi 2 + beta_t*(xtilde_ttilde_ave*sphi - ytilde_ttilde_ave*cphi) Rlong2 = ztilde2_ave return end c----------------------------------------------------------------------------------- c----------------------------------------------------------------------------------- subroutine make_xy_plot(xval,pT,phi_p,mass,iqstat,maxterm) * * This routine is a bastardization of the routine above (BWcalc). * It makes an x-versus-y plot of the emission zone for a given pT and phi_p. * Because it makes a Cartesian plot, it steps in x and y, instead of r and phi_s * as BWcalc. There, it was appropriate to step in r and phi_s, to avoid biasing * effects as FourierCoefficients are calculated. Here, we step in x and y to * avoid binning effects in the x-y histogram * * What is plotted on the 'z' axis is the inteGRAND of {1}_00, i.e. G00*exp(alpha*cos(phi_b-phi_p))*Omega * (I think it's appropriate to include the G00=2*K1(beta), since it DOES matter, as beta depends * on x and y (so it's not "just" an overall normalization), and G00 is the coefficient in the * calculation of and and ). I do leave out any H_i, since that IS "just" * an overall normalization implicit none double precision xval(*) * * index parameter * ----- --------- * 1 Ry * 2 Rx * 3 T * 4 rho_0 * 5 rho_a * 6 tau (peak freezeout time in propertime) * 7 a_skin (*fractional* skin depth) * 8 del-tau (duration of freezeout in proper time) double precision Ry,Rx,T,rho_0,rho_a,tau,a_skin,del_tau integer*4 iqstat ! will be -1 if it's a fermion and +1 if it's a boson integer*4 maxterm ! this is the maximum term we should calculate before truncating ! the series expansion of the quantum statistics ! if maxterm=1, then we only use the first term, which is the ! same as Boltzmann (classical) statistics integer*4 nterm double precision beta_N,signfactor double precision phi_p,rho,phi_b,dphi,spatial_density integer iphi double precision phi_s,r,x,y,xystep double precision fvalue,fvalue_r,f double precision pi,Rmax parameter(pi=3.1415927d0) double precision pT,mT2,pT2,mT double precision mass,mass2,beta_t2,beta_t double precision G00,G01,G02,G20,K0,K1,K2 double precision bessk0,bessk1,bessk * double precision factor double precision alpha,beta,r_ellipse double precision eta2 ! (Ry/Rx)**2 real*8 h common /pawc/h(100000) character*35 htitle integer*4 nbins,icycle,lrec,istat real*4 plotthis,x4,y4,Rmax4 c---- the following is for superimposing many small objects (pp sources) in a larger source (AuAu source) integer*4 isup,ncenters,icenter,ncenters_max parameter(ncenters_max=30000) real*8 xcenter(ncenters_max),ycenter(ncenters_max) ! these are where the protons are placed inside Au real*8 Rcomp,deltax,Rmaxsmall real*4 x4comp,y4comp c--- call hlimit(200000) mass2=mass**2 Ry = dabs(xval(1)) Rx = dabs(xval(2)) T = dabs(xval(3)) rho_0 = dabs(xval(4)) rho_a = xval(5) tau = xval(6) a_skin = xval(7) del_tau = xval(8) pT2 = pT*pT mT2 = mass2+pT2 mT = dsqrt(mT2) beta_t2 = pT2/mT2 beta_t = dsqrt(beta_t2) if (Ry.gt.Rx) then Rmax = Ry else Rmax = Rx endif Rmax = Rmax*(1.0+8.0d0*a_skin) Rmax = Rmax*1.1 htitle = 'phi_p = *.*** pT = *.*** pions' write(htitle(9:13),'(f5.3)')phi_p write(htitle(21:25),'(f5.3)')pT nbins = 100 Rmax4 = sngl(Rmax) call hbook2(10,htitle,nbins,-Rmax4,Rmax4,nbins,-Rmax4,Rmax4,0.) xystep = 2.0d0*Rmax / dfloat(nbins) c--- added April 2004 (in trento) --- c c allow option to superimpose many pp collisions to make a AA collision which is just N*(p+p) c write(6,*)'Hey wanna superimpose many of these objects in ', + 'one large composit one? (1/0)' write(6,*)'[I only deal with round composite objects BTW]' read(5,*)isup if (isup.ne.0) then write(6,*)'What is radius of composite object?' read(5,*)Rcomp write(htitle(28:34),'(a)')'overlay' c now, calculate the centers of all the protons which you will overlay c they need to be inside a radius which is R_composite - R_proton c gotta do this again to get the smaller object radius since Rmax now is something else :(( if (Ry.gt.Rx) then Rmaxsmall = Ry else Rmaxsmall = Rx endif nbins = 400 Rmax4 = sngl(1.1*Rcomp) call + hbook2(20,htitle,nbins,-Rmax4,Rmax4,nbins,-Rmax4,Rmax4,0.) write(6,*)'I will overlay smaller objects inside the', + ' composite in a grid - what should be delta-x between', + ' their centers?' write(6,*)'[something like 1/4 the radius of small object]' read(5,*)deltax ncenters = 0 do x=-Rcomp+deltax/2.0d0,Rcomp-deltax/2.0d0,deltax do y=-Rcomp+deltax/2.0d0,Rcomp-deltax/2.0d0,deltax r=dsqrt(x*x+y*y) if (r.lt.(Rcomp-Rmaxsmall)) then ! put a proton here... ncenters=ncenters+1 if (ncenters.gt.ncenters_max) then stop 'increase ncenters_max' endif xcenter(ncenters) = x ycenter(ncenters) = y endif enddo enddo write(6,*)'OK, I will place ',ncenters,' small objects' endif c--------- end of addition april 2004 ----------- ***no more s2 !!! eta = ((1.0d0+2.0d0*s2)/(1.0d0-2.0d0*s2))**(1.0d0/3.0d0) eta2 = (Ry/Rx)**2 do nterm=1,maxterm signfactor = dfloat(iqstat**(nterm+1)) do x = -Rmax+xystep/2.0d0,Rmax-xystep/2.0d0,xystep x4 = sngl(x) do y = -Rmax+xystep/2.0d0,Rmax-xystep/2.0d0,xystep y4 = sngl(y) r = dsqrt(x*x+y*y) phi_s = datan2(y,x) c===== the "f" calculated here in this sub-block is what we see in the definition of c===== the curly brackets in the paper: {B'}_kj. It is everything except for the B' and the G_jk c===== we also figure out here what is rho, and thus alpha and beta c application of the Theta function: r_ellipse = dsqrt(y*y+eta2*x*x)/Ry spatial_density = + 1.0d0/(1.0d0+dexp((r_ellipse-1.0d0)/a_skin)) ! Woods-Saxon phi_b = datan(dtan(phi_s)/eta2) c c of course the problem with atan() is that it is multiple-valued... c how to get the "right" answer (and it DOES matter) ? c Well, we know that phi_b and phi_s cannot be VERY much different, c so if they begin to differ by pi or more, then we are at the c wrong answer 1 dphi = phi_b - phi_s if (dabs(dphi).gt.3.1415927d0/2.0d0) then phi_b = phi_b - 3.1415927d0 * dphi/dabs(dphi) goto 1 endif rho = (rho_0 + rho_a*dcos(2.0d0*phi_b))*r_ellipse alpha = pT*sinh(rho)/T beta = mT*cosh(rho)/T fvalue = dexp(dfloat(nterm)*alpha*dcos(phi_b-phi_p))* ! note "nterm" in there - mal 13may2004 + spatial_density c===== c===== end of sub-block c===== beta_N = beta*dfloat(nterm) K1 = bessk1(beta_N) G00 = 2.0d0*K1 plotthis = sngl(fvalue*G00) call hfill(10,x4,y4,plotthis) c--- april 2004 c--- OK, now if we are supposed to superimpose pp in AuAu, let's do it here. c--- loop over centers if (isup.ne.0) then do icenter=1,ncenters x4comp = sngl(xcenter(icenter)+x) y4comp = sngl(ycenter(icenter)+y) call hfill(20,x4comp,y4comp,plotthis) enddo endif c--- that was easy!! enddo ! y enddo ! x enddo ! nterm * icycle = 0 lrec = 1024 call hropen(31,'temp2','xyplot.hst',' N',lrec,istat) if (istat.ne.0) write(6,*)'**! error opening hbook file!' if (istat.ne.0) stop icycle = 0 call hrout(10,icycle,' ') if (isup.ne.0) call hrout(20,icycle,' ') call hrend('temp2') close(unit=31) return end c----------------------------------------------------------------------------------- c----------------------------------------------------------------------------------- subroutine open_files(ipid,name,xval,npar,one_phi_only,phival, + maxterm) implicit none integer*4 maxterm integer ipid,npar,ipar,one_phi_only character*10 name(*) double precision xval(*),phival character*180 line if (ipid.eq.1) then open(unit=21,file='pions-vs-pT.BW',status='unknown') open(unit=22,file='pions-vs-phi.BW',status='unknown') write(21,'(a)')'% Boost-invariant results for PIONS' write(22,'(a)')'% Boost-invariant results for PIONS' elseif (ipid.eq.2) then open(unit=21,file='kaons-vs-pT.BW',status='unknown') open(unit=22,file='kaons-vs-phi.BW',status='unknown') write(21,'(a)')'% Boost-invariant results for KAONS' write(22,'(a)')'% Boost-invariant results for KAONS' elseif (ipid.eq.3) then open(unit=21,file='protons-vs-pT.BW',status='unknown') open(unit=22,file='protons-vs-phi.BW',status='unknown') write(21,'(a)')'% Boost-invariant results for PROTONS' write(22,'(a)')'% Boost-invariant results for PROTONS' else stop 'Huh - what pid?' endif write(21,'(a,i2)')'% Nterms into quantum expansion was ',maxterm write(21,'(a)')'% (1 would mean Boltzmann statistics)' write(22,'(a,i2)')'% Nterms into quantum expansion was ',maxterm write(22,'(a)')'% (1 would mean Boltzmann statistics)' if (one_phi_only.eq.1) then write(21,174)phival write(22,174)phival endif write(21,172) write(22,172) do ipar=1,npar write(21,171)name(ipar),xval(ipar) write(22,171)name(ipar),xval(ipar) enddo write(21,172) write(22,172) 171 format('% ',a10,' = ',f6.3) 172 format('%') 174 format('% NOTE !!!! -- Calculation is for phi_p=',f6.3,' ONLY!') * unit 21 (pT dep stuff) has 18 quantities line = '% pT N v2' // + ' R2o0 R2o2 R2o4' // + ' R2s0 R2s2 R2s4' // + ' R2os0 R2os2 R2os4' // + ' R2l0 R2l2 R2l4' // + ' v4 v6 s2' write(21,'(a)')line * unit 22 (pT and phi dep stuff) has 17 quantities line = + '% pT phi_p N ' // + ' ' // + ' ' // + ' Ro2 Rs2 Ros2 Rl2' write(22,'(a)')line return end subroutine calculate_avept_versus_mass(name,xval,npar) implicit none integer npar,ipar character*10 name(*) double precision xval(*) double precision mass,massmin,massmax,massstep double precision pT,pTmin,pTmax,pTstep character*170 line double precision sum,wsum integer nsteps double precision ratio,ratio_old,frac_increase double precision mT,dndpT c these variables are the calculated quantities coming from BWcalc double precision fsum,xtilde2_ave,ytilde2_ave,xtilde_ytilde_ave, + xtilde_ttilde_ave,ytilde_ttilde_ave,ttilde2_ave,ztilde2_ave, + xave,yave,tave,Rout2,Rside2,Ros2,Rlong2 common/BW_answers/fsum,xtilde2_ave,ytilde2_ave, + xtilde_ytilde_ave,xtilde_ttilde_ave,ytilde_ttilde_ave, + ttilde2_ave,ztilde2_ave,xave,yave,tave, + Rout2,Rside2,Ros2,Rlong2 write(6,*)'Input mass min and max and stepsize (GeV/c^2)' read(5,*)massmin,massmax,massstep * write(6,*)'Input pT min and max and stepsize (GeV/c)' * read(5,*)ptmin,ptmax,ptstep open(unit=23,file='pTave-vs-mass.BW',status='unknown') write(23,172) write(23,173) write(23,172) do ipar=1,npar write(23,171)name(ipar),xval(ipar) enddo write(23,172) 171 format('% ',a10,' = ',f6.3) 172 format('%') 173 format('% Average pT versus mass from Boost-Invariant BlastWave') line = '% mass ' write(23,'(a)')line pTstep = 0.04 do mass=massmin,massmax,massstep write(6,*)'Working on mass=',mass sum = 0.0d0 wsum = 0.0d0 ratio_old = 0.0d0 pT = pTstep nsteps = 0 10 continue ! pT loop - loop until wsum/sum changes by less than 1% nsteps = nsteps+1 mT = dsqrt(mass*mass+pT*pT) *** do pT = pTmin,pTmax,pTstep call BWcalc(xval,pT,0.0d0,mass,1,1) ! note classical statistics used here! dndpT = fsum*mT*pT sum = sum + dndpT ! fsum is calculated by BWcalc and passed as common wsum = wsum + dndpT*pT ratio = wsum/sum frac_increase = (ratio-ratio_old)/ratio ratio_old = ratio pT = pT + pTstep if (frac_increase.gt.0.001) goto 10 *** enddo write(23,181)mass,ratio write(6,*)nsteps,' steps for mass=',mass,' - stopped at pT=',pT enddo 181 format(1x,f8.4,2x,f8.4) return end