I have computed the moments $m_n=\int p(x)x^ndx$ of a distribution to be $m_0=1$, $m_1=1$, $m_2=2$ and, in general, $$m_n=\frac{(2n-2)!}{(n-1)!^2}.$$
The moment generating function is $$f(t)=1+\sum_{n=0}^\infty \frac{(2n)!}{n!^2(n+1)!}t^{n+1}=1+te^{2t}\left(I_0(2t)-I_1(2t)\right),$$ in terms of Bessel functions. An inverse Fourier transform $\int e^{-itx}f(it)dt$ then gave me the function $$ \rho(x)=\delta(x)+\frac{1}{\pi x^{3/2}\sqrt{4-x}}.$$
The function $\rho(x)$ indeed has all the correct moments, except for the norm, $m_0$, because it is in fact not normalizable.
This is weird. If the moment generating function converges everywhere, the probability distribution in uniquely determined (right?), but this $\rho$ is not a probability distribution.
Thr general problem of determining a probability distribution from its moments is usually called the "moment problem". See here for some references: https://mathoverflow.net/questions/3525/when-are-probability-distributions-completely-determined-by-their-moments
In your case the result of the inversion suggests that the distribution is a mixture of a point mass and a continuous distribution (if the continuous part has support on ($\epsilon, 4$)). This might explain the spike in the histogram.