Obtaining the temperature from Bose-Einstein and Fermi-Dirac distribution Lets say you are given a distribution function $f(p)$ and you want to define a temperature, $T_f$, for this distribution. (I assume $\mu = 0$.) 
It is then natural to define a temperature the following way:
\begin{equation}
T_f \equiv \frac{ \int d^3p \ G(p) f(p)}{\int d^3p \ f(p)},
\end{equation}
where $G(p)$ is defined by the following equation
\begin{equation}
T = \frac{ \int d^3p \ G(p) f_{eq}(p,T)}{\int d^3p \ f_{eq}(p,T)},
\end{equation}
where $f_{eq}(p,T)$ is the equilibrium thermal distribution function. 
I know that if $f_{eq}$ is given by the Maxwell-Boltzmann distribution, then $$G_{MB}(p) = \frac{p^2}{3E},$$
where $E = \sqrt{p^2 + m^2}$.
What I need is to find an expression for $G(p)$ if $f_{eq}$ is the Bose-Einstein or Fermi-Dirac distribution
$$ f_{eq} = \frac{1}{e^{E(p)/T} \pm 1}.$$
I do not need an analytic expression for $G(p)$, an integral that I can solve numerically is sufficient. It seems to me that this should be possible to do, but I just can't think of how. 
 A: I think such a function may only exist in the Maxwell-Boltzmann limit. Here's why:
For simplicity let us parametrize everything in terms of $\beta = 1/T$ and denote $Z(\beta) = \int{d^3p\; f_{eq}(p, \beta)}$. Rewrite the latter as 
$$
Z(\beta) =  4\pi \int_0^\infty{dp\;\frac{p^2}{e^{\beta E_p}\pm 1}} = 4\pi \int_m^\infty{dE\;\frac{E\sqrt{E^2-m^2}}{e^{\beta E}\pm 1}} = \\
=\frac{4\pi}{3} \int_m^\infty{dE\;\left[\frac{d}{dE}(E^2-m^2)^{3/2}\right]\frac{1}{e^{\beta E}\pm 1}}
$$
and upon integrating by parts,
$$
Z(\beta) = \beta \frac{4\pi}{3} \int_m^\infty{dE\; (E^2-m^2)^{3/2} \frac{e^{\beta E}}{\left(e^{\beta E}\pm 1\right)^2}}
$$
Now, for given $E$ let $e^{\beta E}/\left(e^{\beta E}\pm 1\right)^2$ be the Laplace transform of $\Lambda_\pm(\epsilon; E)$, such that
$$
\frac{e^{\beta E}}{\left(e^{\beta E}\pm 1\right)^2} = \int_0^\infty{d\epsilon \; \Lambda_\pm(\epsilon; E)e^{-\beta \epsilon}}
$$
(I know the phrasing is awkward, but I'm trying to avoid complex plane integration issues with the inverse Laplace transform.) If $\Lambda_\pm(\epsilon; E)$ exist, substitute in $Z(\beta)$, rearrange, and obtain
$$
\frac{Z(\beta)}{4\pi\beta} = \frac{1}{3}\int_0^\infty{d\epsilon \; \left[ \int_m^\infty{dE\; (E^2-m^2)^{3/2}\Lambda_\pm(\epsilon; E) }\right]e^{-\beta \epsilon}}
$$
Basically, this gives us an expression for the Laplace transform of $Z(\beta)/(4\pi\beta)$. Keeping it in mind, let us look for a function $G(p )$ such that 
$$
\frac{1}{\beta} = \frac{1}{Z(\beta)} \int{d^3p\;G(p )f_{eq}(p,\beta)}
$$
or
$$
\frac{Z(\beta)}{\beta} = \int{d^3p\;G(p )f_{eq}(p,\beta)} = 4\pi\int_0^\infty{dp\; \frac{p^2G(p )}{e^{\beta E_p}\pm 1}  } = 4\pi\int_m^\infty{dE\; \frac{{\bar G}(E )E \sqrt{E^2 - m^2}}{e^{\beta E}\pm 1}  } 
$$
where ${\bar G}(E) = G(p )$. As before, let $1/\left(e^{\beta E}\pm 1\right)$ be the Laplace transform of $\Lambda^0_\pm(\epsilon; E)$, such that
$$
\frac{1}{e^{\beta E}\pm 1} = \int_0^\infty{d\epsilon \; \Lambda^0_\pm(\epsilon; E)e^{-\beta \epsilon}}
$$
and obtain
$$
\frac{Z(\beta)}{4\pi\beta} = \int_0^\infty{d\epsilon \; \left[ \int_m^\infty{dE\; E(E^2-m^2)^{1/2}{\bar G}(E)\Lambda^0_\pm(\epsilon; E) }\right]e^{-\beta \epsilon}}
$$
This is yet another expression for the Laplace transform of $Z(\beta)/(4\pi\beta)$. Identifying with the one obtained previously gives 
$$
\int_m^\infty{dE\; E(E^2-m^2)^{1/2}{\bar G}(E)\Lambda^0_\pm(\epsilon; E) } = \frac{1}{3}  \int_m^\infty{dE\; (E^2-m^2)^{3/2}\Lambda_\pm(\epsilon; E) }
$$
or
$$
\int_m^\infty{dE\; (E^2-m^2)^{1/2} \left[ E \;{\bar G}(E)\Lambda^0_\pm(\epsilon; E)  -  \frac{1}{3} (E^2-m^2) \Lambda_\pm(\epsilon; E) \right]} = 0
$$
But ${\bar G}(E)$ has to satisfy this identity for all $\epsilon \ge 0$. This effectively implies that $\Lambda_\pm(\epsilon; E)  = \chi(E)  \Lambda^0_\pm(\epsilon; E)$ for some suitable $\chi(E)$ and in turn means 
$$
\frac{e^{\beta E}}{\left(e^{\beta E}\pm 1\right)^2} = \chi(E) \frac{1}{e^{\beta E}\pm 1}\;\;\; \Rightarrow \;\;\; \chi(E) = \frac{e^{\beta E}}{e^{\beta E}\pm 1}
$$
Since the lhs above is always temperature-independent while the rhs is only in the low-temperature limit, it seems that a proper function ${\bar G}(E)$ can exist only in the same limit, when 
$\Lambda_\pm(\epsilon; E)  \rightarrow  \Lambda^0_\pm(\epsilon; E)$ and
$$
{\bar G}(E) = \frac{E^2 - m^2}{3E} = \frac{p^2}{3E}
$$
as you already know.
