In the heavy quark limit inclusive production rate of a heavy meson can be factorized, in which the nonperturbative effect related to the heavy meson can be characterized by matrix elements defined in the heavy quark effective theory. Using this factorization, predictions for the full spin density matrix of a spin-1 and spin-2 meson can be obtained and they are characterized only by one coefficient representing the nonperturbative effect. Predictions for spin-1 heavy meson are compared with experiment performed at $e^+e^-$ colliders in the energy range from $sqrt{s}=10.5$GeV to $sqrt{s}=91$GeV, a complete agreement is found for $D^*$- and $B^*$-meson. For $D^{**}$ meson, our prediction suffers a large correction, as indicated by experimental data. There exists another approach by taking heavy mesons as bound systems, in which the total angular momentum of the light degrees of freedom is 1/2 and 3/2 for spin-1 and spin-2 meson respectively, then the diagonal parts of spin density matrices can be obtained. However, there are distinct differences in the predictions from the two approaches and they are discussed in detail.