For the distribution to be valid, you need to ask for a distribution of a quantity that is a number or vector for a fixed value of r, whereas you are asking for the trajectory. Perhaps you need to sample the function at multiple instances of times, and then ask for the distribution.
It needs to be set, that DSolve, provided all variable are symbolic can find an analytic solution:
(* try this *)
Zt = DSolve[{-z'[t] == ((
r \[Eta]) (\[Rho] \[Alpha] g r z[t] - 2 cos\[Theta] \[Gamma]))/(
8 z[t]), z[0] == L}, z, t]
The solution is real for numeric values of his and for some r and t, but it may only represent a branch of the global solution. Then something like
SmoothHistogram[
RandomVariate[
TransformedDistribution[Re[z[1.]] /. First[Zt] /. subst,
r \[Distributed] NormalDistribution[\[Mu], \[Sigma]] /. subst],
10^5], PlotRange -> All]
works, but shows that the solution is not likely to be a global one, since it has imaginary part.
Hence, you can use ParametricNDSolve too.