Interpolating four-dimensional data

Before doing anything, you have to remove any duplicate $(x,y)$ pairs from your data. Interpolating functions don't like double-valued functions. 10 of your points are like these two:

data[[{3197, 3250}]]
(* {{10.20913`, -0.18115`, -0.15105`, 50.42686`}, 
    {10.20913`, -0.18115`, -0.15105`, 50.42685`}} *)

So just throw out the duplicates via

data = DeleteDuplicatesBy[data, #[[;; 2]] &];
xy = data[[All, ;; 2]];

How can I generate, let's say 100000, random (x,y) initial conditions inside the oval region of the above plot?

You can define the oval region using ConvexHullMesh, and then use RandomPoint to generate your points.

xy = data[[All, ;; 2]];
pts = RandomPoint[ConvexHullMesh[xy], 100000];
Graphics[{Point[pts], Red, Point[xy]}]

Assuming that we have generated 100000 (x,y) initial conditions, inside the oval region, how can we exploit the original data in order to interpolate the momenta $p_x$ and $p_y$ of the additional initial conditions?

Now just make your interpolation functions. You need to set the InterpolationOrder to 1 because the $(x,y)$ points don't lie on a rectangular grid. Check the documentation for Interpolation to see why you have to rearrange the list structure.

px = Interpolation[{{#1, #2}, #3} & @@@ data, 
  InterpolationOrder -> 1]
py = Interpolation[{{#1, #2}, #4} & @@@ data, 
   InterpolationOrder -> 1];

and generate the new data via

newData = {#1, #2, px[##], py[##]} & @@@ pts;

Here is the new data (in blue) along with the old data (in red) for the px variable

Graphics3D[{
  Red, PointSize@Medium, Point[data[[All, ;; 3]]],
  Blue, PointSize@Small, Point[
   newData[[All, ;; 3]]]}]

Edit If you are working in an older version of Mathematica, before there was a ConvexHullMesh or RandomPoint, you can still make this work, but the only way I could think was to use rejection sampling to get the new data points.

positionDuplicates[list_] := 
    Flatten[Rest /@ GatherBy[Range@Length[list], list[[#]] &]]; 
data = ReplacePart[data, 
    Thread[positionDuplicates[data[[All, ;; 2]]] -> Sequence[]]]; 
xy = data[[All, ;; 2]]; 
px = Interpolation[{{#1, #2}, #3} & @@@ data, InterpolationOrder -> 1]; 
py = Interpolation[{{#1, #2}, #4} & @@@ data, InterpolationOrder -> 1]; 
range = {Min@#, Max@#} & /@ Transpose[xy]; 
pointsBeta = Transpose[{ 
    RandomReal[First@range, 100000], 
    RandomReal[Last@range, 100000] 
    } 
]; 
pgon = Polygon[ 
    data[[Graphics`Mesh`ConvexHull[xy], ;; 2]] 
]; 
pts = Select[pointsBeta, Graphics`Mesh`InPolygonQ[pgon, #] &];
newData = {#1, #2, px[##], py[##]} & @@@ pts; 

Graphics3D[{ 
    Red, PointSize@Medium, Point[data[[All, ;; 3]]], 
    Blue, PointSize@Small, Point[ newData[[All, ;; 3]]]
    }, Method -> {"ShrinkWrap" -> True}]

should give the same result, but only giving 80,000 points instead of 100,000. Code borrowed and adapted from this answer and this answer.

For the random generation part, once you have an interpolation function you can actually just let the interpolation function tell you when points are out of the region:

points = Select[ RandomReal[{-2, 2}, {200, 3}], 
                Norm[#[[1 ;; 2]]] < 1 &];
int = Interpolation[points, InterpolationOrder -> 1]
random = DeleteCases[
   Quiet@Check[Append[#, int @@ #], "Fail"] & /@
    RandomReal[{Min@Flatten[points], Max@Flatten[points]}, {2000, 2}],
     "Fail"];
Show[{
  ConvexHullMesh[points[[All, 1 ;; 2]]],
  Graphics@Point[random[[All, 1 ;; 2]]]}]

note I used ConvexHullMesh here only for the illustration graphic.