Deriving a PDF for an annulus distribution

An alternative, and faster, way to generate samples with the desired distribution using TransformedDistribution:

ClearAll[distF]
distF[r_, σ_] := TransformedDistribution[{(r + d ) Sin[θ], (r + d ) Cos[θ]}, 
  {Distributed[d, NormalDistribution[0, σ]], Distributed[θ,
                UniformDistribution[{0, 2 Pi}]]}]
ListPlot[RandomVariate[distF[5, .5] , 10000], PlotStyle -> PointSize[Tiny]]

data = RandomVariate[distF[5, .5] , 1000000];
empdist = SmoothKernelDistribution[data];
ContourPlot[PDF[empdist, {x, y}], {x, -8, 8}, {y, -8, 8}, PlotPoints -> 50]

Table[fTorusRand[5, 0.5], {i, 1, 1000000}]; // AbsoluteTiming // First

13.499060

RandomVariate[distF[5, .5] , 1000000]; // AbsoluteTiming // First

3.012873

Sorry for the trouble - I have found a way:

data = Table[fTorusRand[5, 0.5], {i, 1, 100000}];
empD1 = SmoothKernelDistribution[data];
ContourPlot[PDF[empD1, {x, y}], {x, -8, 8}, {y, -8, 8}]

Gives me what I want:

Edit: So using the Jacobian here, I can get an exact PDF. Still messing around with the algebra though: $x = (r+d) \sin(\theta)$ and $y = (r+d) \cos(\theta)$. So I can implicitly differentiate these to find the Jacobians...

This is a naive approach, so I wonder if I am overlooking some complication here; I would have thought that, given the description of your distribution, you could consider the "solid of revolution" generated by rotating the PDF of a normal distribution with the required parameters around the vertical $z$ axis, so something like the following:

PDF[NormalDistribution[6, 1], x] /. x -> Sqrt[x^2 + y^2]

(* Out: E^(-(1/2) (-6+Sqrt[x^2+y^2])^2)/Sqrt[2 Pi] *)

Plot3D[
  Evaluate[PDF[NormalDistribution[6, 1], x] /. x -> Sqrt[x^2 + y^2]],
  {x, -10, 10}, {y, -10, 10}, PlotPoints -> 75
]

ContourPlot[
  Evaluate[PDF[NormalDistribution[6, 1], x] /. x -> Sqrt[x^2 + y^2]],
  {x, -10, 10}, {y, -10, 10}, PlotPoints -> 75
]

Of course, as it is the expression is not properly normalized, so an appropriate factor should be added to make sure that its integral is equal to 1.