Here a lot gets confused. The reason is, as quite often, the abuse of the idea of photons. Let's take the case of black-body radiation in a cavity which is among the very few examples, where the quantum theory of free fields is sufficient for an accurat description of an observable phenomenon, namely black-body radiation.
It starts with the free electromagnetic field, which reads in radiation gauge, i.e., the gauge, where A^0=0 and \vec{\nabla} \cdot \vec{A}=0 on the classical level, and with this gauge constraint you quantize the radiation field. This is most conveniently done in terms of Fock space of the momentum-single-photon space for photons in a cavity, which we choose as a cubic box of length L and impose periodic boundary conditions (you could as well work with rigid boundary conditions, but that's less convenient). The field operator reads
\hat{A}^{\mu}(x)=\frac{1}{\sqrt{L^3}} \sum_{\lambda=\pm 1} \sum_{\vec{p}} [\hat{a}(\vec{p},\lambda) \vec{\epsilon}(\vec{p},\lambda) \exp(-[\mathrm{i} \omega t+\mathrm{i} \vec{x} \cdot \vec{k}]+\text{h.c.}].
Here p^0=|\vec{p}| (massless (!) photons). The \lambda labels the helicity, i.e., sums over the two circular polarized modes of the em. fields, and \vec{p} \in \frac{2 \pi}{L}.
The energy (Hamilton) and momentum operators are calculated from the normal ordered Belinfante energy-momentum tensor (or equivalently from the canonical energy momentum tensor, because total energy and momentum do not depend on this):
\hat{P}^{\mu} =\sum_{\vec{p},\lambda} \hat{N}(\vec{p},\lambda) p^{\mu}, \quad \text{with} \quad \hat{N}(\vec{p},\lambda)=\hat{a}^{\dagger}(\vec{p},\lambda) \hat{a}(\vec{p},\lambda).
The Hilbert space is spanned by the Fock-space basis vektors
|\{N(\vec{p},\lambda) \} \rangle = \prod_{\vec{p},\lambda} \frac{1}{\sqrt{N(\vec{p},\lambda)!}} (\hat{a}^{\dagger}(\vec{p},\lambda))^{N(\vec{p},\lambda)} |\Omega \rangle,
where |\Omega \rangle is the vacuum state, which obeys \hat{a}(\vec{p},\lambda)|\Omega \rangle =0 for all (\vec{p},\lambda).
Now we can look at the equilibrium state. Since photons carry no conserved charges, we work with the canonical ensemble. The covariant statistical operator for a non-rotating state (0 total angular momentum) reads
\hat{R}=\frac{1}{Z} \exp(-u \cdot \hat{P}).
Here, u is a time-like vector, appearing as Lagrange parameters for the maximum-entropy principle to keep the average energy and momentum as constraints.
The partition sum is calculated with the above Fock-space (occupation number) basis
Z=\mathrm{Tr} \exp(-u \cdot \hat{P})=\prod_{\vec{p},\lambda} \frac{1}{1-\exp(-u \cdot p)}.
In the limit of a very big container this partition sum can be expressed in the form
\ln Z=-2 \frac{V}{(2 \pi)^3} \int_{\mathbb{R}^3} \ln[1-\exp(-u \cdot p)].
The factor 2 comes from the two helicity states for each mode of the em. field.
The partition sum can be evaluated exactly. Without loss of generality we can assume that \vec{u}=u_3 \vec{e}_z and introducing spherical coordinates gives
\ln Z=-2 \frac{V}{(2 \pi)^2} \int_0^{\infty} \mathrm{d} p \int_{-1}^1 \int_{-1}^1\mathrm{d} \eta \ln[1-\exp[-\beta p(u^0-u_3 \eta)]] = \frac{V u_0 \pi^2}{45 (u \cdot u)^2}.
On the first glance this looks not manifestly covariant, but the reason is that we evaluated the partition sum in the lab frame, where the cavity with the radiation inside is moving.
We can make the partition sum manifestly covariant, by expressing everything with respect to the restframe of the cavity (which is at the same the of course the frame where the radiation inside the cavity has total 0 momentu either). To set end we set u=(\beta,0,0,0). In the general case u=\beta/\sqrt{1-v^2}(1,\vec{v}). Then we have to take into account the fact that the length contraction gives V=V_0 \sqrt{1-v^2}=V_0 \frac{\beta}{u^0}.
This gives the canonical potential
\Omega=-\ln Z=\frac{\pi^2 V_0}{45 (u \cdot u)^{3/2}},
where we have written u \cdot u=\beta^2. This has the advantage that we can take derivatives with independent variables u_{\mu} to get the average energy and momentum. At the end of the calculation we use u \cdot u=\beta^2 again to express everything in the scalar quantities V_0 and \beta=1/T (T: temperature). This gives finally
\langle P^{\mu}_{\text{rad}} \rangle=-\frac{\partial}{\partial u_{\mu}}=\frac{\pi^2 V_0}{45 \beta^4} \frac{1}{\sqrt{1-v^2}} (1,\vec{v}).
This shows that the total mass of the radiation is given by its energy in the rest frame, which is
m_{\text{rad}}=E_{0,\text{rad}}=\frac{\pi^2 V_0}{45 \beta^4},
as is well known from textbooks.
This reflects the true meaning of the famous formula "E=m", which of course refers to the rest frame of the moving massive body.
The em. field itself is, of course, a massless field according to the Standard Model of particles, and this assumption is one among the best experimentally tested ones in whole physics.
