R=6.68;
a=0.546;
A=208;
σ=6.4 ;(* fm^2 *)
ρ[r_]=((ρ0)/(1+Exp[((r-R)/a)]));
stepsize = 1;

(* rho_0 computation *)
  ρ0 = A ρ0/ (Integrate[  ρ[r] r^2, {r,0,Infinity }]*4Pi)

(* TA thickness function *)

TA [x,y] =Integrate[ρ [Sqrt[x^2+y^2+z^2]] / A  ,{z,-Infinity,Infinity}];
Plot3D[NIntegrate[ρ[Sqrt[x^2+y^2+z^2] ]/ A  ,{z,-Infinity,Infinity}],{x,-10,10},{y,-10,10},PlotRange->All,AxesLabel->{x,y,TA}]

(* TAA  integrand *)
b=0;
Plot3D[(NIntegrate[ρ [Sqrt[(x+b/2)^2+y^2+z^2]] /A ,{z,-Infinity,Infinity}])*(NIntegrate[ρ [Sqrt[(x-b/2)^2+y^2+z^2]] /A ,{z,-Infinity,Infinity}]),{x,-10,10},{y,-10,10},PlotRange->All,AxesLabel->{x,y,TAA Integrand}]

integrandNucleiOverlapFunction=Interpolation[Flatten[ Table[ {{{x,y,b},Nof=NIntegrate[ρ[Sqrt[(x+b/2)^2+y^2+z^2]]/A,{z,-∞,∞},Method-> Automatic]NIntegrate[ρ[Sqrt[(x-b/2)^2+y^2+z^2]]/A,{z,-∞,∞},Method-> Automatic]},{{-x,y,b},Nof},{{x,-y,b},Nof},{{-x,-y,b},Nof}},{x,0.0001,10.0001,stepsize},{y,0.0001,10.0001,stepsize},{b,0,20,stepsize}],3]];

Animate[Plot3D[integrandNucleiOverlapFunction[x,y,b],{x,-10,10},{y,-10,10},AxesLabel->Automatic] ,{b,0,12,4}]

(* TAA  nuclei overlap function *)

TAA =Table[ NIntegrate[integrandNucleiOverlapFunction[x,y,b],{x,-10.0001,10.0001},{y,-10.0001,10.0001}],{b,0,20,stepsize}];

Ncoll = A^2*σ *TAA;
ClearAll[b];
yData=Ncoll;
xData=Range[0,20,stepsize];
ListPlot[Transpose[{xData,yData}],AxesLabel->{b,Ncol}]

(* Collison probability *)
CollProb = 1- Exp[- Ncoll];

ImpactParamter = Table[b,{b,0,20,stepsize}];
dCrossSection =   CollProb Take[ImpactParamter * 2 Pi, Length[CollProb]];
yData=dCrossSection;
xData=Range[0,20,stepsize];
ClearAll[b,σ]
ListPlot[Transpose[{xData,yData}],AxesLabel->{b,dσ}]

TotalCrossSection= stepsize * With[{l=Length[dCrossSection]},Sum[dCrossSection[[i]],{i,1,l}]] fm^2