This is my solution. Might still contain mistakes. Hope this helps others too!
The bessel function
Using the following properties
\begin{equation}
\int_0^\infty dq q J_v(q \rho_1) J_v(q \rho_2)= \frac{\delta(r_1-r_2)}{r_1}
\end{equation}
\begin{equation}
J'_v(z)=- J_{v+1}(z)+\frac{v}{z}J_{v}(z)
\end{equation}
\begin{equation}
J_0'(z)=-J_1(z)
\end{equation}
\begin{equation}
J_{-n}=(-1)^n J_{n}.
\end{equation}
Therefore;
\begin{equation}
q^3 J'_{-1}(q \rho_1) J'_{-1}(q \rho_2)= q^3 J_0(q \rho_1) J_0(q \rho_2)+ q \frac{1}{\rho_1 \rho_2} J_{-1}(q\rho_1)J_{-1}(q\rho_2)+q^2\left( \frac{1}{\rho_2}J_0(q \rho_1) J'_0(q \rho_2)+\frac{1}{\rho_1}J'_0(q \rho_1) J_0(q \rho_2) \right)
\end{equation}
This gives
\begin{equation}
\int_0^\infty dq q^3 J_0(q \rho_1) J_0(q \rho_2)=\left(\partial_{\rho_1}\partial_{\rho_2}-\frac{1}{\rho_1 \rho_2}-\frac{\partial_{\rho_1}}{\rho_1}-\frac{\partial_{\rho_2}}{\rho_2}\right)\frac{\delta(\rho_1-\rho_2)}{\rho_1}
\end{equation}
Solving the integral
\begin{equation}
I=\int_0^\infty d \rho f(\rho) \left(-\frac{\partial^2_{\rho}}{\rho}+\frac{3\partial_{\rho}}{\rho^2}-\frac{3}{\rho^3} \right)f(\rho).
\end{equation}
define $x\equiv \rho^2$.
\begin{equation}
I=-\int_0^\infty dx f(x)\left[ 2\partial^2_x-2\frac{\partial_x}{x}+\frac{3}{2x^2}\right] f(x).
\end{equation}
Using the identities
\begin{equation}
\partial_x(x L^1_p(x))=(p+1) L_{p}^0
\end{equation}
\begin{equation}
\partial_x(L_p^0)=-L_{p-1} (\text{for}~ p\geq1) ~(=0~ \text{otherwise})
\end{equation}
the derivatives to $f(x)=x L_p^1(x) e^{-x/2}$ are
\begin{equation}
\partial_x f(x)=(- \frac{x}{2} L_p^1(x)+(p+1)L_p^0(x) ) e^{x/2}
\end{equation}
\begin{equation}
\partial^2_x f(x)=(\frac{x}{4} L_p^1(x)-(p+1)L_p^0(x)-L^1_{p-1}(!)) e^{x/2}
\end{equation}
where the $(!)$ points out that this term is only there is $p\geq1$.
This gives the following integral
\begin{equation}
I=-\int_0^\infty dx e^{-x} \left[L_p^1(x)L_p^1(x)(\frac{3}{2}+x+\frac{x^2}{2}) -x (p+1) L_p^0(x) L_p^1(x)- 2 x L_{p-1}^1(x) L_{p}^1(x)\right]
\end{equation}
Orthogonality conditions
Making us of the orthogonality conditions
\begin{equation}
\int_0^\infty dx x^a e^{-x} L_n^a L_m^a =\frac{(n+a)!}{n!} \delta_{n,m}
\end{equation}
\begin{equation}
\int_0^\infty dx x^{a+1} e^{-x} (L_n^a(x))^2 =\frac{(n+a)!}{n!} (2n+a+1)
\end{equation}
and
\begin{equation}
L_p^1=\sum_{i=0}^p L_{i}^0(x)
\end{equation}
\begin{equation}
L_n^a=L_n^{a+1}-L_{n-1}^{a+1}.
\end{equation}
then
\begin{equation}
I=-(3p/2+(p+1)+(p+1)(2p+1+1)/2-(p+1)(p+1)-0)=-\frac{5p}{2}-1
\end{equation}