  Our goal is to obtain a complete set of angular observables arising in a
generic multi-body process. We show how this can be achieved without the need
to carry out a likelihood fit of the angular distribution to the measured
events. Instead, we apply the method of moments that relies both on the
orthogonality of angular functions and the estimation of integrals by Monte
Carlo techniques. The big advantage of this method is that the joint
distribution of all observables can be easily extracted, even for very few
events. The method of moments is shown to be robust against mismodeling of the
angular distribution. Our main result is an explicit algorithm that accounts
for systematic uncertainties from detector-resolution and acceptance effects.
Finally, we present the necessary process-dependent formulae needed for direct
application of the method to several rare decays of interest.
