Blast-wave Model

Klaus Reygers, reygers@physi.uni-heidelberg,de, 17 June 2015

Clear["Global`*"]

Parameter

Rp:=6 (* fm *)

Calculate surface velocity from average transverse velocity

βs=((2+n)/2)β;

βs/.{β->0.651, n->0.712}

Transverse rapidity

ρ[r_?NumericQ, β_?NumericQ, n_?NumericQ]:=ArcTanh[(((r/R))^n) ((2+n)/2)β]/.R->Rp

Plot[ρ[r, 0.6, 1], {r, 0., Rp},FrameLabel->{Style["radius (fm)",16],Style["transverse rapidity ",16]},Frame->True]

Transverse velocoty profile as a function of r

vProfile=Plot[(((r/R))^n) ((2+n)/2)β /.{β->0.651, n->0.712, R->Rp},{r,0,6},Frame->True,FrameLabel->{Style["radius (fm)",16], Style["transverse velocity (Subscript[β, T])",16]}]

Export[NotebookDirectory[]<>"blastwave_velocity_profile.png",vProfile]

Blastwave formula for particle pT spectra f = 1/(2 PI pT) dN/dpTdy

f[pt_?NumericQ, A_?NumericQ,β_?NumericQ, T_?NumericQ, n_?NumericQ, m_?NumericQ]:=A NIntegrate[r (Sqrt[(m^2)+(pt^2)])BesselI[0,((pt Sinh[ArcTanh[(((r/Rp))^n) ((2+n)/2)β]])/T)]BesselK[1,(((Sqrt[(m^2)+(pt^2)])Cosh[ArcTanh[(((r/Rp))^n) ((2+n)/2)β]])/T)],{r,0,Rp}, AccuracyGoal->4, Method->"LocalAdaptive"]

Parameters for protons and pions (for Pb-Pb, 2.76 TeV, 0-5% most central, from arXiv:1303.0737)

parsProton={m->0.938,dndy->34,β->0.651, T->0.095, n->0.712};

parsProton2={m->0.938,dndy->34,β->0., T->0.16, n->0.712};

parsPion={m->0.140,dndy->733,β->0.651, T->0.095, n->0.712};

parsPion2={m->0.240,dndy->733,β->0, T->0.095, n->0.712};

AProton=NIntegrate[2 π pt f[pt,1, β, T, n, m]/.parsProton,{pt, 0., 10}];

AProton2=NIntegrate[2 π pt f[pt,1, β, T, n, m]/.parsProton2,{pt, 0., 10}];

APion=NIntegrate[2 π pt f[pt,1, β, T, n, m]/.parsPion,{pt, 0., 10}];

APion2=NIntegrate[2 π pt f[pt,1, β, T, n, m]/.parsPion2,{pt, 0., 10}];

Spectra

col1=ColorData[97,"ColorList"][[1]];

col2=ColorData[97,"ColorList"][[2]];

pionplot=LogPlot[{f[pt,dndy/APion, β, T, n, m]/.parsPion,f[pt,dndy/APion2, β, T, n, m]/.parsPion2}, {pt, 0.,2},Frame->True, FrameLabel->{Style["(Subscript[p, T]) (GeV/c)",16], Style["\!\(\*FractionBox[\(1\), \(2\\ π\\ \*SubscriptBox[\(p\), \(T\)]\)]\)(((d^2) N)/((Subscript[dp, T]) dy))",16]},Epilog->{Inset[Style["⟨(Subscript[β, T])⟩ = 0.65",12,col1],{1.5,Log[50]}],Inset[Style["⟨(Subscript[β, T])⟩ = 0",12,col2],{1,Log[0.1]}],Inset[Style["pions",12],{0.2, Log[1/20000]}]}]

Export[NotebookDirectory[]<>"blastwave_pions.png",pionplot]

protonplot=LogPlot[{f[pt,dndy/AProton, β, T, n, m]/.parsProton,f[pt,dndy/AProton2, β, T, n, m]/.parsProton2}, {pt, 0.,2},Frame->True, FrameLabel->{Style["(Subscript[p, T]) (GeV/c)",16], Style["\!\(\*FractionBox[\(1\), \(2\\ π\\ \*SubscriptBox[\(p\), \(T\)]\)]\)(((d^2) N)/((Subscript[dp, T]) dy))",16]},Epilog->{Inset[Style["⟨(Subscript[β, T])⟩ = 0.65",12,col1],{1.5,Log[5]}],Inset[Style["⟨(Subscript[β, T])⟩ = 0",12,col2],{1,Log[0.7]}],Inset[Style["protons",12],{0.3, Log[1.5/100]}]}]

Export[NotebookDirectory[]<>"blastwave_protons.png",protonplot]