This mirrors the classical microcanonical phase volume.

We have

$Z(\beta):=\frac{1}{h^{3N}N!}\int_{\Gamma_{\mathcal M}} \mathrm{e}^{-\beta\ H({\bf q},{\bf p})}\ \mathrm d\Gamma $

$U\equiv\langle H\rangle=-\dfrac{\partial}{\partial\beta}\log\,Z(\beta)$

The usage of $Z(\beta)$ is very similar to that of a characteristic function/moment generating function for a probability density $f(x)$:

$E[e^{-ikX}]:=\int f(x)\ \mathrm e^{-ikx}\ \mathrm dx$

in probability theory, where

$E[X^n]=\left(i\frac{\partial}{\partial k}\right)_{k=0}^nE[e^{-ikX}].$

In fact, to compute expectation values with the partition function is even a little more straight forward, since in that case the probability function itself has exponential form and the temperature acts like the otherwise auxiliary parameter $k$. Notice a second difference: The parameter $\beta$ is multiplied by the energy $H(q,p)$, not directly by the variables of integration $q,p$.

See density of states

$ Z(\beta)=\int{\mathrm e}^{-\beta\,E}{\mathrm d}\varphi(E)=\int D(E)\ \mathrm{e}^{-\beta\ E}\ \mathrm d E $

$U\equiv\langle H\rangle=-\dfrac{\partial}{\partial\beta}\log\,Z(\beta)=\int E\,\dfrac{{\mathrm e}^{-\beta\,E}}{Z(\beta)}\,D(E)\,{\mathrm d}E$

Note that in QM, there is no phase spaces to smoothly integrate over.

Experimental observation says that certain objects (the ones involved in experiments measuring the energy spectrum from a box, the sun, etc.) the itegrand $E\,\dfrac{{\mathrm e}^{-\beta\,E}}{Z(\beta)}\,D(E)$ ought to be replaces by $\propto E\,\dfrac{1}{\mathrm e^{\beta\hbar\omega}-1}\,E^2$.

For a single particle and $E\propto \omega\propto |k|$, it can indeed by argued that $D(E)\propto E^2$, see density of states. The temperature dependency of the energy distribution however isn't observed to follow a Boltzmann Distribution but instead a statistics known from the Grand canonical partition function.

Planck's ad hoc step is to say that any radiation ray of energy $E$ is actually partitioned into $n$ parts of $(\hbar\omega)$. Here $n$ isn't fixed but partioning is instead also thermailzed (that's basically the Definition of a black body) and thus follows a Boltzmann distribution in $(n\hbar\omega)$. This amounts to raplacing certain functions of energy (they are energy densities, if we introduce spatial variables into the model)

Roughly

$\int {\mathrm d}E\cdot p_E \, \mapsto \, \sum_{n=0}^\infty\int{\mathrm d}\,(n\hbar\omega)\,p_\omega(n)\,\mapsto\,\sum_{n=0}^\infty\int{\mathrm d}\,(n\hbar|k|)\,|k|^2\,p_k(n)$

If we compare the classical expression with the empirical result, the replacement must be

$E\,\dfrac{{\mathrm e}^{-\beta\,E}}{Z(\beta)}\mapsto\sum_{n=0}^\infty\dfrac{(n(\hbar\omega))\,{\mathrm e}^{-\beta\,n\,(\hbar\omega)}}{\sum_{n=0}^\infty{\mathrm e}^{-\beta\,n\,(\hbar\omega)}}=(\hbar\omega)\dfrac{1}{\mathrm e^{\beta\hbar\omega}-1}$

The common derivation of Plancks law makes the above step and then neglects spatial inhomogenies and just introduces a characteristic length $L$. Then we can use the characteristic speed (of light) $c$ to get a characteristic frequency $c/L$. Thus we can intorduce powers of $\omega$ via the unitless expression $\omega\left/\right.\dfrac{c}{L}$. If $D$ is quadratic (the case of photon gas), we get

Classical canonical ensemble

Classical density of states