As answered by RaphaelB4 it seems best first to simulate a random point on $\mathcal{S}$ and then to apply the rejection method, see f.i. L. Devroye, Non Uniform Random Variate Generation (1986). For the second step you need the maximum $m$ of $x \to \Delta(X,p)$ and then you have to generate random $U \in [0,m]$. (I have not investigated if $m \leq 2$.) For high dimension $d$ the simulation of a random point $y \in \mathcal{S}$ is not trivial. One possibility is to generate $d$ $N(0,1)$ independent values $(y_1,\ldots,y_d)$. Then $\frac{1}{\sqrt{\sum_{i=1}^d y_i^2}}(y_1,\ldots,y_d)$ is a random point on $\mathcal{S}$.

As answered by RaphaelB4 it seems best first to simulate a random point on $\mathcal{S}$ and then to apply the rejection method, see f.i. L. Devroye, Non Uniform Random Variate Generation (1986). For the second step you need the maximum $m$ of $x \to \Delta(X,p)$ and then you have to generate random $U \in [0,m]$. (I have not investigated if $m \leq 2$.) For high dimension $d$ the simulation of a random point $y \in \mathcal{S}$ is not trivial. One possibility is to generate $d$ $N(0,1)$-distributed independent values $(y_1,\ldots,y_d)$. Then $\frac{1}{\sqrt{\sum_{i=1}^d y_i^2}}(y_1,\ldots,y_d)$ is a random point on $\mathcal{S}$.