R=6.68;
a=0.546;
A=208;
σ=6.4 ;(* fm^2 *)
ρ[r_?NumericQ]=((ρ0)/(1+Exp[((r-R)/a)]));
(Subscript[ρ, c])[x_?NumericQ,y_?NumericQ,z_?NumericQ] := ρ0/(1+ Exp[(Sqrt[x^2+y^2+z^2]-R)/a]) ;
stepsize = 1; 
(* rho_0 computation *)
  ρ0 = A ρ0/ (Integrate[  ρ[r] r^2, {r,0,Infinity }]*4Pi)

(* TA thickness function *)

TA [x_?NumericQ,y_?NumericQ] :=(2.0/208.0)*NIntegrate[(Subscript[ρ, c])[x,y,z],{z,0.0,15.0},Method->{Automatic,"SymbolicProcessing"-> 0},PrecisionGoal->5]

(Subscript[ϵ, 2])[b_?NumericQ]:= NIntegrate[(y^2-x^2)TA[x-b/2,y]TA[x+b/2,y],{x,0.0,15.0},{y,0.0,15.0},Method->Automatic,PrecisionGoal->4,AccuracyGoal->20]/NIntegrate[(y^2+x^2)TA[x-b/2,y]TA[x+b/2,y],{x,0.0,15.0},{y,0.0,15.0},Method->Automatic,PrecisionGoal->4,AccuracyGoal->20]

DiscretePlot[(Subscript[ϵ, 2])[b],{b,0,15,stepsize}]

Ellipticity=ParallelTable[{b,(Subscript[ϵ, 2])[b]},{b,0.0,15.0,0.5}]

ListLinePlot[Ellipticity,Frame->True,FrameLabel->{"b [fm]","(Subscript[ϵ, 2])"}]

(Subscript[ϵ, 4])[b_?NumericQ]:= NIntegrate[(6x^2 y^2 - x^4 -y^4)TA[x-b/2,y]TA[x+b/2,y],{x,0.0,15.0},{y,0.0,15.0},Method->Automatic,PrecisionGoal->5,AccuracyGoal->6]/NIntegrate[(x^4+y^4+2x^2 y^2)TA[x-b/2,y]TA[x+b/2,y],{x,0.0,15.0},{y,0.0,15.0},Method->Automatic,PrecisionGoal->5,AccuracyGoal->6]

DiscretePlot[(Subscript[ϵ, 4])[b],{b,0,15,stepsize}]

Quadrangularity=ParallelTable[{b,(Subscript[ϵ, 4])[b]},{b,0.0,15.0,0.5}]

ListLinePlot[Quadrangularity,Frame->True,FrameLabel->{"b [fm]","(Subscript[ϵ, 4])"}]

(*  Two particle distribution Exercise 11 *)

v2 = 0.041;
v3 =0.053;
v4 = 0.034;
v5 =0.012;
Plot[1+2 v2^2 Cos[2Δφ]+2 v3^2 Cos[3Δφ]+2 v4^2 Cos[4Δφ]+2 v5^2 Cos[5Δφ],{Δφ,-Pi/2,3Pi/2},AxesLabel-> {Δφ,p},AxesOrigin->{0,1},Background->None]