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.

What you've posted isn't exactly in Mathematica format (and I suspect that you probably want lowercase Greek letters rather than uppercase). I think you'll get more help if you use the Code sample icon (left-most unlabeled icon - a blank page with < >) and/or attach a Mathematica notebook.