Below is my confirmation that the distribution should be proportional to $E^{n-1} e^{-E/kT}$. I originally did these calculations for massive particles, which lead to a different result than what the OP was expecting, and so I posted it as an answer. However, since my results now agree with the OP's, this no longer constitutes a possible answer, and I don't know what's going on.
The Boltzmann distribution in velocitymomentum space for massless particles in $n$ dimensions is $$ f(\vec{v}) \, d^n\vec{v} \propto e^{-\beta mv^2/2} \, d^n\vec{v}, $$$$ f(\vec{p}) \, d^n\vec{p} \propto e^{-\beta c p} \, d^n\vec{p}, $$ with $\beta = 1/kT$. To transform this into a distribution with respect to energy, we change to polar coordinates in velocitymomentum space, in which $$ d^n \vec{v} = v^{n-1} \, dv \, d\Omega_v $$$$ d^n \vec{p} = p^{n-1} \, dp \, d\Omega_p $$ with $d \Omega_v$$d \Omega_p$ representing the solid angle in velocity space. Integrating over solid angle just yields a constant that can be folded into the normalization, so the distribution with respect to speeds is $$ f(v) \, dv \propto v^{n-1} e^{-\beta m v^2/2} \, dv $$$$ f(p) \, dp \propto p^{n-1} e^{-\beta cp} \, dp $$ Finally, to transform this into a distribution with respect to energy, we note that $v = \sqrt{2E/m}$$p = E/c$ and so $dv = 1/\sqrt{2mE} \, dE$$dp = dE/c$; and thus $$ f(E) \, dE \propto E^{(n-1)/2} e^{-\beta E} \, \frac{dE}{\sqrt{E}} = E^{n/2 - 1} e^{-\beta E} \, dE. $$$$ f(E) \, dE \propto E^{n-1} e^{-\beta E} \, dE. $$ Thus, the distribution should be $f(E) \propto \sqrt{E} e^{-E/kT}$$f(E) \propto E^2 e^{-E/kT}$ in 3D and $f(E) \propto e^{-E/kT}$$f(E) \propto E e^{-E/kT}$ in 2D. So your 2D code is working as expected, but your 3D code does not seem to be yielding the expected result.
Note that this pre-factor of $E^{n/2 - 1}$ in the distribution is just (proportional to) the density of states with respect to energy.