Optical Glauber Calculation for A+A Collisions

Klaus Reygers, reygers@physi.uni-heidelberg,de, 21 May 2015

Clear["Global`*"]

Parameters

User input

SigmaNN:=6.4 (* Nucleon-nucleon cross section in fm^2: Pb+Pb at LHC: 6.4, Au+Au at RHIC 4.2 *)

A:=208 (* Pb:208, Au: 197 *)

R:=6.62 (* Radius in fm: Pb: 6.62, Au: 6.38 *)

a:= 0.546 (* diffuseness parameter in fm *)

bmax:=20 (* maximum impact parameter in fm *)

bstep:=0.2(* stepsize in fm in calculation of Npart(b), etc *)

Define volume that fully subtends the nucleus (and thus also the overlap region)

xmax= R+5 a (* transverse area that fully includes the larger of the two nuclei *)

ymax =  xmax; (* transverse area that fully includes the larger of the two nuclei *)

zmax=xmax;

Calculation of thickness function T

Define nuclear geometry

rho:= rho0/(1+Exp[(r-Rp)/ap]) (* nuclear density distribution as a function of the radius r *)

norm=Integrate[4 Pi r^2 rho/.rho0->1, {r, 0, Infinity},Assumptions->{ap>0,Rp>0}]/.{Rp->R,ap->a}; (* normalization *)

rhonorm[s_?NumericQ,z_?NumericQ]=rho /. {r-> Sqrt[s^2+z^2] ,Rp->R,ap->a, rho0->A/norm} ;(*3D nucleon density distribution with correct normalization *)

t[s_?NumericQ]:=NIntegrate[rhonorm[s,z],{z,-zmax,zmax},AccuracyGoal->4, Method->"LocalAdaptive"] (* nuclear thickness function*)

Monitor[tint=Interpolation[Table[{s,t[s]},{s,0.,20,bstep}]] ,s];

Plot[tint[s],{s,0,10.},AxesLabel->{"transverse radius s (fm)", "nucleons/area T(s)"}]

Calculate nuclear overlap function

tab[b_?NumericQ]:=NIntegrate[tint[(Sqrt[(((x+b/2))^2)+(y^2)])] tint[(Sqrt[(((x-b/2))^2)+(y^2)])],{x,-xmax,xmax},{y,-ymax,ymax},AccuracyGoal->4, Method->"LocalAdaptive"] (* define nuclear overlap function TAB *)

Monitor[tabint=Interpolation[Table[{b,tab[b]},{b,0.,bmax,bstep}]] ,b];

Plot[tabint[b], {b,0,bmax}, PlotPoints->20, AxesOrigin->{0,0}, AxesLabel->{"b (fm)", "(Subscript[T, AB])(b)"}]

Calculate interaction probability and dσ/db

pint[b_?NumericQ]:=1-Exp[-SigmaNN tabint[b]]

Plot[2 Pi b pint[b], {b,0,bmax}, PlotPoints->20, AxesOrigin->{0,0}, AxesLabel->{"b (fm)", "dσ/db (fm)"}]

Calculate Fraction of Total Cross Section as a Function of bmax

sigtot=NIntegrate[2 Pi b pint[b],{b,0,bmax}] (* total cross section in fm^2 *)

sigfrac[bmax_?NumericQ]:=NIntegrate[1/sigtot 2 Pi b pint[b],{b,0,bmax}]

Plot[sigfrac[bmax],{bmax,0,18},PlotPoints->20,AxesOrigin->{0,0}, AxesLabel->{"(Subscript[b, max]) (fm)", "fraction of (Subscript[σ, tot])"}]

Calculate number of participants

npart[b_?NumericQ]:= NIntegrate[tint[(Sqrt[(((x-b/2))^2)+(y^2)])] (1-Exp[-tint[(Sqrt[(((x+b/2))^2)+(y^2)])] SigmaNN])+ tint[(Sqrt[(((x+b/2))^2)+(y^2)])] (1-Exp[-tint[(Sqrt[(((x-b/2))^2)+(y^2)])] SigmaNN]),{x,-xmax,xmax},{y,-ymax,ymax}, AccuracyGoal->6, Method->"LocalAdaptive"] (* total number of participants *)

Monitor[npartint=Interpolation[Table[{b,npart[b]},{b,0.,bmax,bstep}]] ,b];(* interpolate Npart(b) *)

npartinelint[b_?NumericQ]:=npartint[b]/pint[b]

Plot[npartint[b], {b,0,bmax}, PlotPoints->20, AxesOrigin->{0,0}, AxesLabel->{"b (fm)", "(Subscript[N, part])(b)"}]

Table: Impact parameter b, Ncoll(b), Npart(b), fraction of total cross section

table=Table[{b,NumberForm[SigmaNN tabint[b]/pint[b],{5,1}],NumberForm[npartinelint[b],{4,1}],NumberForm[100 sigfrac[b],3]},{b,0,bmax,1}] ;

TableForm[table,TableHeadings-> {None,{"b (fm)","Ncoll","Npart","σ/(Subscript[σ, tot]) (%)"}},TableAlignments->"."]