noOfTrials = 500;
points = RandomPoint[Sphere[], {noOfTrials, 2}];
(*gives the length of a line from pt1 to pt2 that lies inside a
 sphere at origin of radius r*)
dFunc[pt1_, pt2_, r_] := Module[
  {midPt, h},
  midPt = (pt1 + pt2)*0.5;
  h = Norm[midPt];
  If[h >= r, 0, 2 N@Sqrt[r^2 - h^2]]
]
(*simulation*)
sim = Module[
{dr, l, dlTotal},
dr = 10^-2;
 Table[
 dlTotal =
  Plus @@ Table[(dFunc[#1, #2, r + dr] - dFunc[#1, #2, r]) & @@
     pair, {pair, points}];
 {r, dlTotal/(4 Pi r^2 dr)}
 , {r, 0.01, 1, dr}
 ]
];
(*plot*)
ListPlot[sim]