Complicated contour integral to be solved. Anyone can help to solve the following integral? 
$$I=\int_{0}^{\infty} dp p^{-1}e^{-2p^{2}M^{-2}}\sin(pr)\frac{M^2}{M^2+p^2}$$
at this stage I am able to write the integral as 
$$I=\int_{0}^{\infty} dp p^{-1}e^{-2p^{2}M^{-2}}\frac{M^2(\exp (i p r)-\exp (-i p r))}{2i(p+iM)(p-iM)}$$
if the integral was without $p^{-1}e^{-2p^{2}M^{-2}}$ part, one could use the general contour integral formula. But now things are not as simple. So the question is how to solve the integral?

Note: $M > 0$ and $r > 0$ (i.e. +ve constants) 
 A: Adopting @Winther's suggestion, substitute $q=p/M$ and write $R=Mr$, so the integral we seek is
$$ I(R) = \int_0^{\infty} e^{-2q^2} \frac{\sin{Rq}}{q} \frac{dq}{1+q^2}. $$
Then
$$ -I'''(R) + I'(R) = \int_0^{\infty} e^{-2q^2} \cos{Rq} \, dq \tag{1} $$
The integral on the right can be done by the usual methods (more differentiation under the integral sign, for example, or completing the square &c.) to give $ \sqrt{\frac{\pi}{8}} e^{-R^2/8} $. Integrating from $0$ to $R$, $I''(0)=I(0)=0$, so
$$ -I''(R)+I(R) = \frac{\pi}{2} \operatorname{erf}{\left( R/\sqrt{8} \right)} $$
A basis of solutions is given by $ u_1(R) = e^R $ and $u_2(R) = e^{-R}$, and so the Wronskian is -2, and variation of parameters gives
$$ I(R) = \int_0^R \frac{1}{2}(-e^{R-r}+e^{r-R}) \frac{\pi}{2} \operatorname{erf}{\left( r/\sqrt{8} \right)} \, dr + A\sinh{R} \\
= -\frac{\pi}{2}\int_0^R \sinh{(R-r)} \operatorname{erf}{\left( r/\sqrt{8} \right)} \, dr + A\sinh{R}
$$
$$ A=I'(0) = \int_0^{\infty} \frac{e^{-2q^2}}{1+q^2} \, dq = \frac{e^2\pi}{2}(1-\operatorname{erf}(\sqrt{2})) $$
(this can be computed using the Schwinger parametrisation, for example.)
The final integral can be done using integration by parts repeatedly (or interchanging the order of integration), and we end up with the final answer
$$ I(R) = \frac{e^2\pi}{4} \sinh{R} -\frac{e^2\pi}{4} e^{-R} \operatorname{erf}\left(\frac{R-4}{\sqrt{8}}\right)+\frac{\pi}{2}  \operatorname{erf}\left(\frac{R}{\sqrt{8}}\right)-\frac{e^2\pi}{4}  e^{R} \operatorname{erf}\left(\frac{R+4}{\sqrt{8}}\right) $$

Edited to add: One can find the differential equation as follows: notice that
$$ \frac{d}{dR} \sin{Rq} = q\cos{Rq}, $$
so one differentiation under the integral sign gets rid of the $q^{-1}$. I still didn't much fancy evaluating
$$ I'(R) = \int_0^{\infty} e^{-2q^2} \cos{Rq} \frac{dq}{1+q^2}, $$
so I then used that
$$ \left( -\frac{d^2}{dR^2} + 1 \right) \cos{Rq} = (q^2+1)\cos{Rq} $$
to get rid of the $\frac{1}{1+q^2}$, which gives the equation (1).
Now, to evaluate
$$ J = \int_0^{\infty} e^{-2q^2} \cos{Rq} \, dq = \frac{1}{2} \int_{-\infty}^{\infty} e^{-2q^2} e^{iRq} \, dq, $$
we could look in a table of Fourier transforms, or complete the square:
$$ -2q^2 + iRq = -2\left( q - \frac{iR}{4} \right)^2 - \frac{R^2}{8}, $$
and then we have a standard Gaussian integral times $e^{-R^2/8}$. (We do need to shift the integration contour back to the real line, of course.)
A: Thanks @cheppers ok I re-derived most of the parts, can you please elaborate on this part please: 
A basis of solutions is given by $ u_1(R) = e^R $ and $u_2(R) = e^{-R}$, and so the Wronskian is -2, and variation of parameters gives
$$ I(R) = \int_0^R \frac{1}{2}(-e^{R-r}+e^{r-R}) \frac{\pi}{2} \operatorname{erf}{\left( r/\sqrt{8} \right)} \, dr + A\sinh{R} \\
= -\frac{\pi}{2}\int_0^R \sinh{(R-r)} \operatorname{erf}{\left( r/\sqrt{8} \right)} \, dr + A\sinh{R}
$$
$$ A=I'(0) = \int_0^{\infty} \frac{e^{-2q^2}}{1+q^2} \, dq = \frac{e^2\pi}{2}(1-\operatorname{erf}(\sqrt{2})) $$
(this can be computed using the Schwinger parametrisation, for example.)
The final integral can be done using integration by parts repeatedly (or interchanging the order of integration), and we end up with the final answer
$$ I(R) = \frac{e^2\pi}{4} \sinh{R} -\frac{e^2\pi}{4} e^{-R} \operatorname{erf}\left(\frac{R-4}{\sqrt{8}}\right)+\frac{\pi}{2}  \operatorname{erf}\left(\frac{R}{\sqrt{8}}\right)-\frac{e^2\pi}{4}  e^{R} \operatorname{erf}\left(\frac{R+4}{\sqrt{8}}\right) $$
