Convert the 3D luminosity distribution to mass using a bunch of mass particles.
Well, you probably could do that, but I suspect it would be inaccurate and extremely time-consuming. Most Schwarzschild modeling approaches parameterize the 3D luminosity-density distribution using functions from which the gravitational potential can be relatively easily calculated.
For example, one popular approach is to approximate the observed 2D light distribution using the sum of multiple Gaussians (as you mention), which can be deprojected under certain assumptions into the sum of multiple 3D Gaussians to get the 3D luminosity-density distribution. Then you multiply this distribution by a mass-to-light ratio (usually one of the parameters of the model you're ultimately trying to fit to the data) and derive the gravitational potential from that. No need to approximate things with a bunch of particles.
This also allows you to easily add extra hypothetical "dark" components, like dark-matter halos or central supermassive black holes, which have simple analytical descriptions and don't need to modeled with particles.
To my understanding, the orbits of the particles in a Schwarzschild model can be counter-rotating.
Correct. In practice, if the system is (assumed to be) axisymmetric, you can calculate a set of orbits, and then easily generate their counter-rotating counterparts by just reversing the velocity vectors.
Yet, all papers on dynamical modelling using the Schwarzschild method fit models to velocity and velocity dispersion maps instead of a spectral cube.
You're clearly not reading the right papers ;-).
Let me divide this into two parts:
1. Converting spectroscopic data to stellar-kinematic data
Before any modeling takes place, the observed spectroscopic data (long-slit and/or integral-field-unit "spectral cubes") must be analyzed to extract the stellar-kinematic data. This is done by dividing up the data into different spatial bins on the sky (bins along the slit, or 2D bins for a spectral cube) and summing the spectra within each bin, to improve the S/N.
Then you fit each spectrum with a weighted combination of stellar spectra (observed or from models), broadened by a model1 line-of-sight velocity distribution (LOSVD). In very simple cases, the LOSVD is assumed to be just a Gaussian with a specific central velocity $V$ and velocity dispersion $\sigma$. The more common present-day approach is to model the LOSVD with a Gaussian which is multiplied by third- and fourth-order Gauss-Hermite polynomials. The strength of the 3rd-order polynomial is parameterized by $h_{3}$ and measures skew; the 4th-order polynomial is parameterized by $h_{4}$ and measures kurtosis. For galaxies with counter-rotating components, $h_{3}$ and $h_{4}$ provide some information about the non-single-Gaussian nature of the LOSVD, although it's rather fuzzy and poorly constrained.
A more sophisticated approach (requiring high-S/N data) is to "non-parametrically" model the LOSVD with a histogram having twenty, thirty, or more velocity bins. To show how the non-parametric approach works, consider the following stellar-kinematic data for the S0 galaxy NGC 4191:
The first plot, from Fig. B1 of Krajnovic et al. 2011 (MNRAS, 414, 2923), shows $V$ and $\sigma$ maps (from the SAURON IFU). In the velocity field, you can see that the velocity reverses between $r \sim 5$ and 10 arc seconds; at roughly the same radius, $\sigma$ peaks. The paired off-center peaks in $\sigma$ are suggestive of a region where the two counter-rotating components are about equal (as argued by Krajnovic et al. in their paper).
The second plot, based on Fig. 2 of Coccato et al. 2015 (A&A, 581, A65), shows the $\sigma$ map (colored panel, rotated 90 degrees with respect to the figure above) and non-parametric LOSVDs (31 velocity bins) for three of the bins (histograms) from VIRUS-W IFU observations of NGC 4191. You can see how, in the second histogram, the double peaks due to the overlapping, counter-rotating stellar components are clearly represented.
2. Comparing the Schwarzschild model with the kinematic data
After computing the library of stellar orbits for a given galaxy potential and assigning a set of weights to the individual orbits, you "observe" the model by orienting the galaxy to match the observed orientation, dividing it up into projected 2D spatial bins that match the 2D binning used for determining real galaxy's stellar kinematics. Then (for each bin) you compute the weighted LOSVD by integrating along the line of sight through that bin, adding up the radial-velocity contributions from each orbit that intersects the line of sight. If your data consists of $V$, $\sigma$, $h_3$, and $h_4$, then you fit the resulting model LOSVD for that bin with the Gauss-Hermite function to get the model bin's $V$, $\sigma$, $h_3$, and $h_4$. If your data are non-parametric LOSVDs, then you can simply bin the model LOSVD using the same velocity-binning scheme you used for the data.
This set of predicted LOSVD values (either $V$, $\sigma$, $h_3$, and $h_4$ or the non-parametric binned velocity histogram) is what you compare with data's LOSVD values for that spatial bin, repeating for all the spatial bins with kinematic data.
1 "Model" in this case having nothing to do with the Schwarzschild modeling.
