Gamma-Distribution-like integral

Let $\gamma$ be Hankle's contour. It is well-known that: $$\dfrac{1}{\Gamma(z)}=\dfrac{i}{2\pi}\int_{\gamma}(-t)^{-z}e^{-t}{\rm d}t$$ STEP 1: Considering the integral: \begin{align*} J(a)& =\int_{0}^{\infty}\dfrac{x^xe^{-ax}}{\Gamma(x+1)}{\rm d}x\\ & =\dfrac{i}{2\pi}\int_0^{\infty}e^{-ax}{\rm d}x\int_{\gamma}(-t)^{-x-1}e^{-tx}{\rm d}t\\ & =-\dfrac{i}{2\pi}\int_{\gamma}\dfrac{{\rm d}t}{t(a+t+\log(-t))}\\ & =-\dfrac{1}{1+W_{-1}(-e^{-a})}, \end{align*} where $W_{-1}(z)$ is Lambert W function. The last equation is right due to residue theorem.

STEP 2: Hence \begin{align*} I & = e\int_{0}^{\infty}{\rm d}x\int_{1}^{\infty}\dfrac{x^xe^{-a(1+x)}}{\Gamma(1+x)}{\rm d}a\\ & =e\int_1^{\infty}e^{-a}J(a){\rm d}a\\ & =-e\int_1^{\infty}\dfrac{e^{-a}}{1+W_{-1}(-e^{-a})}{\rm d}a\\ & =-e\int_0^{-1/e}\dfrac{{\rm d}t}{1+W_{-1}(t)}\\ & =e\int_{-\infty}^{-1}e^W_{-1}{\rm d}W=1 \end{align*} The last equation is right due to $W'_{-1}(z)(1+W_{-1}(z))e^{W_{-1}(z)}=1.$

Corollary：$$\int_0^{\infty}\frac{x^{x+p-1}e^{-x}}{\Gamma(x+p+1)}dx=\frac{1}{p}.$$ Here Proposition 5.5 at Page 27

The same question has been discussed on MathOverflow, see here