Add optional process support for biased sampling (for focus, point detectors, etc).

[tkittel]

This is related to the issue of "differential cross sections", but loosely speaking (ignoring non-isotropic materials for now) one of the two variables (E,mu) is sampled and the other is used as a differential cross section.

Such features could be used for biasing in our "mini stepping mc" (#101), adding focusing to McStas components, implementing "point detectors", etc.

The features would need new virtual functions on the basic process interface classes, and thus break ABI. It might be an idea to add them when we add other ABI breaking changes (e.g. adding vectorized versions of cross section evaluations, etc).

[tkittel]

Based on my own musings + discussions with @marquezj, @dddijulio, @willend, here are a few ideas for how such interfaces would look. The first (towards a particular mu=cos(scatangle)) would work best with processes where XS varies smoothly, while the latter (towards a mu range) would be needed to deal with other processes (e.g. bragg diffraction) -- and would need either the neutron flight path or the detector extend to be something more than a point. I am not sure yet if the weight fields should contain a probability (or probability density), or also include cross section (using N=0 in the new methods would actually not do any sampling, so the method name could be a bit misleading here). I guess this is why it is best if we do not introduce such functions without actually also implementing something with them at the same time.

      //If processType is Scatter and materialType is Isotropic, it *might* be
      //possible to query the sampling (this can be used to implement biasing,
      //focusing, or efficient tallying). Not all processes might support this,
      //as indicated by the canSampleTowardsXXX() methods - and using the
      //sampleTowardsXXX(..) methods might consequently result in exceptions
      //being raised.

      //Get the scattering contribution towards a particular direction (mu), and
      //sample up to N outgoing energy values in that direction. While the
      //actual number of resulting energy values will typically be the requested
      //N, it might be lower. In particular it should be 0 if there is no
      //probability to scatter in that direction (so also the weight will 0),
      //and if there is only one possible outgoing energy value (e.g. an elastic
      //process), only a single entry might be returned. It is OK to specify
      //N=0, if only the weight is desired.
      struct ScatOutcomesIsotropic_TowardsMu {
        double weight; // p(mu)?, or also include cross section? Should the unit be barn? (i.e. CrossSect?)
        SmallVector<NeutronEnergy,16> ekin;
      };
      virtual bool canSampleTowardsMu() const = 0;
      virtual ScatOutcomesIsotropic_TowardsMu sampleTowardsMu( CachePtr&, RNG&,
                                                               NeutronEnergy,
                                                               CosineScatAngle target_mu,
                                                               std::size_t N = 1 ) const = 0;

      //Similar, but considering a range of mu values. This form is also
      //suitable for sharply peaked differential cross sections (e.g. Bragg
      //diffraction has delta-functions in p(mu).
      struct ScatOutcomesIsotropic_TowardsMuRange {
        double weight;//p(mu_min<=mu<=mu_max)??
        SmallVector<ScatterOutcomeIsotropic,16> outcomes;//list of sampled (mu,ekin) values
      };
      virtual bool canSampleTowardsMuRange() const = 0;
      virtual ScatOutcomesIsotropic_TowardsMu sampleTowardsMuRange( CachePtr&, RNG&,
                                                                    NeutronEnergy,
                                                                    CosineScatAngle target_mu_min,
                                                                    CosineScatAngle target_mu_max,
                                                                    std::size_t N = 1 ) const = 0;

[tkittel]

I am also not super happy with the potential memory allocation in case N>16, so perhaps it could be worth considering to go old-school and instead take a non-const reference to a vector of results instead (and provide the weight in the return value).