Feat madspace ramon

madspace/include/madspace/compgraphs/function_builder_mixin.inc

@@ -300,6 +300,16 @@ std::array<Value, 3> t_inv_value_and_min_max(Value pa, Value pb, Value p1, Value
     return {output_vector[0], output_vector[1], output_vector[2]};
 }
 
+std::array<Value, 2> t_inv_min_max_cut(Value pa, Value pb, Value m1, Value m2, Value t_min_cut) {
+    auto output_vector = instruction("t_inv_min_max_cut", {pa, pb, m1, m2, t_min_cut});
+    return {output_vector[0], output_vector[1]};
+}
+
+std::array<Value, 3> t_inv_value_and_min_max_cut(Value pa, Value pb, Value p1, Value p2, Value t_min_cut) {
+    auto output_vector = instruction("t_inv_value_and_min_max_cut", {pa, pb, p1, p2, t_min_cut});
+    return {output_vector[0], output_vector[1], output_vector[2]};
+}
+
 std::array<Value, 2> t1_inv_min_max_doublet(Value pa, Value pb, Value m1, Value mir_min) {
     auto output_vector = instruction("t1_inv_min_max_doublet", {pa, pb, m1, mir_min});
     return {output_vector[0], output_vector[1]};
@@ -330,6 +340,16 @@ std::array<Value, 3> s23_value_and_min_max(Value pa, Value pb, Value p3, Value t
     return {output_vector[0], output_vector[1], output_vector[2]};
 }
 
+std::array<Value, 2> s23_min_max_cut(Value pa, Value pb, Value p3, Value t1_abs, Value m1, Value m2, Value s23_min_cut) {
+    auto output_vector = instruction("s23_min_max_cut", {pa, pb, p3, t1_abs, m1, m2, s23_min_cut});
+    return {output_vector[0], output_vector[1]};
+}
+
+std::array<Value, 3> s23_value_and_min_max_cut(Value pa, Value pb, Value p3, Value t1_abs, Value p1, Value p2, Value s23_min_cut) {
+    auto output_vector = instruction("s23_value_and_min_max_cut", {pa, pb, p3, t1_abs, p1, p2, s23_min_cut});
+    return {output_vector[0], output_vector[1], output_vector[2]};
+}
+
 Value invariants_from_momenta(Value p_ext, Value factors) {
     return instruction("invariants_from_momenta", {p_ext, factors})[0];
 }

madspace/include/madspace/compgraphs/opcode_mixin.inc

@@ -67,86 +67,90 @@ three_body_decay = 65,
 three_body_decay_inverse = 66,
 t_inv_min_max = 67,
 t_inv_value_and_min_max = 68,
-t1_inv_min_max_doublet = 69,
-t1_inv_value_and_min_max_doublet = 70,
-t2_inv_min_max_doublet = 71,
-t2_inv_value_and_min_max_doublet = 72,
-s23_min_max = 73,
-s23_value_and_min_max = 74,
-invariants_from_momenta = 75,
-sde2_channel_weights = 76,
-subchannel_weights = 77,
-apply_subchannel_weights = 78,
-pt_eta_phi_x = 79,
-mirror_momenta = 80,
-momenta_to_x1x2 = 81,
-uniform_invariant = 82,
-uniform_invariant_inverse = 83,
-breit_wigner_invariant = 84,
-breit_wigner_invariant_inverse = 85,
-stable_invariant = 86,
-stable_invariant_inverse = 87,
-stable_invariant_nu = 88,
-stable_invariant_nu_inverse = 89,
-fast_rambo_massless = 90,
-fast_rambo_massless_inverse = 91,
-fast_rambo_massless_com = 92,
-fast_rambo_massive = 93,
-fast_rambo_massive_inverse = 94,
-fast_rambo_massive_com = 95,
-cut_unphysical = 96,
-cut_one = 97,
-cut_all = 98,
-cut_any = 99,
-scale_transverse_energy = 100,
-scale_transverse_mass = 101,
-scale_half_transverse_mass = 102,
-scale_partonic_energy = 103,
-chili_forward = 104,
-chili_inverse = 105,
-matrix_element = 106,
-collect_channel_weights = 107,
-interpolate_pdf = 108,
-interpolate_alpha_s = 109,
-matmul = 110,
-relu = 111,
-leaky_relu = 112,
-elu = 113,
-gelu = 114,
-sigmoid = 115,
-softplus = 116,
-rqs_reshape = 117,
-rqs_find_bin = 118,
-rqs_forward = 119,
-rqs_inverse = 120,
-softmax = 121,
-softmax_prior = 122,
-sample_discrete = 123,
-sample_discrete_inverse = 124,
-sample_discrete_probs = 125,
-sample_discrete_probs_inverse = 126,
-discrete_histogram = 127,
-permute_momenta = 128,
-gather = 129,
-gather_int = 130,
-select_int = 131,
-select = 132,
-select_vector = 133,
-argsort = 134,
-quantile = 135,
-one_hot = 136,
-madnis_abs_weight = 137,
-madnis_softclip = 138,
-madnis_variance = 139,
-madnis_single_channel_variance = 140,
-madnis_multi_channel_variance = 141,
-nonzero = 142,
-batch_gather = 143,
-batch_scatter = 144,
-random = 145,
-random_int = 146,
-unweight = 147,
-vegas_forward = 148,
-vegas_inverse = 149,
-vegas_histogram = 150,
-histogram = 151
+t_inv_min_max_cut = 69,
+t_inv_value_and_min_max_cut = 70,
+t1_inv_min_max_doublet = 71,
+t1_inv_value_and_min_max_doublet = 72,
+t2_inv_min_max_doublet = 73,
+t2_inv_value_and_min_max_doublet = 74,
+s23_min_max = 75,
+s23_value_and_min_max = 76,
+s23_min_max_cut = 77,
+s23_value_and_min_max_cut = 78,
+invariants_from_momenta = 79,
+sde2_channel_weights = 80,
+subchannel_weights = 81,
+apply_subchannel_weights = 82,
+pt_eta_phi_x = 83,
+mirror_momenta = 84,
+momenta_to_x1x2 = 85,
+uniform_invariant = 86,
+uniform_invariant_inverse = 87,
+breit_wigner_invariant = 88,
+breit_wigner_invariant_inverse = 89,
+stable_invariant = 90,
+stable_invariant_inverse = 91,
+stable_invariant_nu = 92,
+stable_invariant_nu_inverse = 93,
+fast_rambo_massless = 94,
+fast_rambo_massless_inverse = 95,
+fast_rambo_massless_com = 96,
+fast_rambo_massive = 97,
+fast_rambo_massive_inverse = 98,
+fast_rambo_massive_com = 99,
+cut_unphysical = 100,
+cut_one = 101,
+cut_all = 102,
+cut_any = 103,
+scale_transverse_energy = 104,
+scale_transverse_mass = 105,
+scale_half_transverse_mass = 106,
+scale_partonic_energy = 107,
+chili_forward = 108,
+chili_inverse = 109,
+matrix_element = 110,
+collect_channel_weights = 111,
+interpolate_pdf = 112,
+interpolate_alpha_s = 113,
+matmul = 114,
+relu = 115,
+leaky_relu = 116,
+elu = 117,
+gelu = 118,
+sigmoid = 119,
+softplus = 120,
+rqs_reshape = 121,
+rqs_find_bin = 122,
+rqs_forward = 123,
+rqs_inverse = 124,
+softmax = 125,
+softmax_prior = 126,
+sample_discrete = 127,
+sample_discrete_inverse = 128,
+sample_discrete_probs = 129,
+sample_discrete_probs_inverse = 130,
+discrete_histogram = 131,
+permute_momenta = 132,
+gather = 133,
+gather_int = 134,
+select_int = 135,
+select = 136,
+select_vector = 137,
+argsort = 138,
+quantile = 139,
+one_hot = 140,
+madnis_abs_weight = 141,
+madnis_softclip = 142,
+madnis_variance = 143,
+madnis_single_channel_variance = 144,
+madnis_multi_channel_variance = 145,
+nonzero = 146,
+batch_gather = 147,
+batch_scatter = 148,
+random = 149,
+random_int = 150,
+unweight = 151,
+vegas_forward = 152,
+vegas_inverse = 153,
+vegas_histogram = 154,
+histogram = 155

madspace/include/madspace/phasespace/base.hpp

@@ -50,6 +50,7 @@ class Mapping {
     const NamedVector<Type>& output_types() const { return _output_types; }
     const NamedVector<Type>& condition_types() const { return _condition_types; }
     const std::string& name() const { return _name; }
+    virtual std::size_t discrete_dim() const { return 0; }
 
 protected:
     // TODO: make parameters const ref

madspace/include/madspace/phasespace/color_ordered_mapping.hpp

@@ -14,13 +14,31 @@ class ColorOrderedMapping : public Mapping {
 public:
     // color_order: 0-indexed permutation of {0, ..., n-1} (n = n_out + 2).
     // Particles 0 and 1 are the two incoming beams.
+    //
+    // Optional cuts (all indexed by 0-based outgoing-particle index, i.e. the
+    // same indexing as the masses m_out and the entries of the color order
+    // minus 2):
+    //   * pt_min[i]        : minimum transverse momentum of outgoing particle i.
+    //   * m_inv_min[i][j] : minimum invariant mass of the pair (i, j).
+    //   * dr_min[i][j]     : minimum delta-R separation of the pair (i, j).
+    // Passing any non-empty cut container enables cut-aware sampling. The cuts
+    // are translated into invariant-space bounds exactly as in the Fortran
+    // reference (phase_space_gen23): per-subset invariant-mass floors
+    // (invm_min) on the sampled s-channel masses, the adjacent-pair floor on
+    // the 2->3 s23 invariant, and a |t| floor (pt^2) on each peeled particle.
+    // Empty containers (the default) reproduce the previous cut-free behaviour
+    // exactly.
     ColorOrderedMapping(
         const std::vector<std::size_t>& color_order,
         double t_invariant_power = 0.8,
-        double s_invariant_power = 0.8
+        double s_invariant_power = 0.8,
+        const std::vector<double>& pt_min = {},
+        const std::vector<std::vector<double>>& m_inv_min = {},
+        const std::vector<std::vector<double>>& dr_min = {}
     );
 
     std::size_t random_dim() const { return _random_dim; }
+    std::size_t discrete_dim() const override { return _discrete_dim; }
 
 private:
     Result build_forward_impl(
@@ -34,13 +52,24 @@ class ColorOrderedMapping : public Mapping {
         const NamedVector<Value>& conditions
     ) const override;
 
+    // pt^2 of outgoing particle i (0 if no pt cut on it).
+    double pt2(std::size_t i) const;
+    // Cut-derived invariant-mass^2 floor (gen23 invm_min, without the mass^2
+    // term, which is applied separately) for a subset of outgoing particles.
+    // Returns 0 when cuts are disabled or the subset has fewer than 2 members.
+    double cut_floor(const std::vector<std::size_t>& subset) const;
+
     // 0-indexed outgoing-particle indices (values in {0,...,n_out-1}).
     // _set1 contains the outgoing particles attached to beam 0's side,
     // _set2 those attached to beam 1's side, in peel order.
     std::vector<std::size_t> _set1;
     std::vector<std::size_t> _set2;
     std::size_t _n_out;
     std::size_t _random_dim;
+    // Number of discrete two-solution choices (one per 2->3 peel). These are
+    // supplied/recovered as a separate batch_int channel, not through the
+    // continuous random_dim() block (opt-in r_disc).
+    std::size_t _discrete_dim;
     // True iff exactly one of (set1, set2) has size 1 (and the other >= 2).
     // In that case the central block is DoubleT instead of 2->2.
     bool _use_double_t;
@@ -50,6 +79,11 @@ class ColorOrderedMapping : public Mapping {
     // produced as a single t-channel chain seeded directly off the beams.
     bool _use_single_chain;
 
+    // Cut configuration (empty => all bounds resolve to 0 = no cut).
+    std::vector<double> _pt_min;
+    std::vector<std::vector<double>> _m_inv_min;
+    std::vector<std::vector<double>> _dr_min;
+
     Invariant _uniform_invariant;
     TwoToTwoParticleScattering _com_scattering;
     TwoToTwoParticleScattering _lab_scattering;

madspace/include/madspace/phasespace/cuts.hpp

@@ -23,6 +23,8 @@ class Cuts : public FunctionGenerator {
     double sqrt_s_min() const;
     std::vector<double> eta_max() const;
     std::vector<double> pt_min() const;
+    std::vector<std::vector<double>> m_inv_min() const;
+    std::vector<std::vector<double>> dr_min() const;
 
 private:
     NamedVector<Value> build_function_impl(

madspace/include/madspace/phasespace/observable.hpp

@@ -48,6 +48,8 @@ class Observable : public FunctionGenerator {
         const std::string& name = ""
     );
     ObservableOption observable() const { return _observable; }
+    const nested_vector2<me_int_t>& indices() const { return _indices; }
+    bool sum_momenta() const { return _sum_momenta; }
     std::vector<std::size_t> simple_observable_indices() const {
         if (_sum_momenta || _sum_observable || _indices.size() != 1) {
             return {};

madspace/include/madspace/phasespace/phasespace.hpp

@@ -2,6 +2,7 @@
 
 #include "madspace/phasespace/base.hpp"
 #include "madspace/phasespace/chili.hpp"
+#include "madspace/phasespace/color_ordered_mapping.hpp"
 #include "madspace/phasespace/cuts.hpp"
 #include "madspace/phasespace/invariants.hpp"
 #include "madspace/phasespace/luminosity.hpp"
@@ -14,7 +15,7 @@ namespace madspace {
 
 class PhaseSpaceMapping : public Mapping {
 public:
-    enum TChannelMode { propagator, rambo, chili };
+    enum TChannelMode { propagator, rambo, chili, color_ordered };
 
     PhaseSpaceMapping(
         const Topology& topology,
@@ -23,7 +24,8 @@ class PhaseSpaceMapping : public Mapping {
         double invariant_power = 0.8,
         TChannelMode t_channel_mode = propagator,
         const std::optional<Cuts>& cuts = std::nullopt,
-        const std::vector<std::vector<std::size_t>>& permutations = {}
+        const std::vector<std::vector<std::size_t>>& permutations = {},
+        const std::optional<std::vector<std::size_t>>& color_order = std::nullopt
     );
 
     PhaseSpaceMapping(
@@ -32,12 +34,14 @@ class PhaseSpaceMapping : public Mapping {
         bool leptonic = false,
         double invariant_power = 0.8,
         TChannelMode mode = rambo,
-        const std::optional<Cuts>& cuts = std::nullopt
+        const std::optional<Cuts>& cuts = std::nullopt,
+        const std::optional<std::vector<std::size_t>>& color_order = std::nullopt
     );
 
     std::size_t random_dim() const {
         return 3 * _topology.outgoing_masses().size() - (_leptonic ? 4 : 2);
     }
+    std::size_t discrete_dim() const override { return _n_discrete; }
     std::size_t particle_count() const {
         return _topology.outgoing_masses().size() + 2;
     }
@@ -61,8 +65,14 @@ class PhaseSpaceMapping : public Mapping {
     double _sqrt_s_lab;
     bool _leptonic;
     bool _map_luminosity;
+    std::size_t _n_discrete;
     std::vector<Invariant> _s_invariants;
-    std::variant<TPropagatorMapping, FastRamboMapping, ChiliMapping, std::monostate>
+    std::variant<
+        TPropagatorMapping,
+        FastRamboMapping,
+        ChiliMapping,
+        ColorOrderedMapping,
+        std::monostate>
         _t_mapping;
     std::vector<std::variant<TwoBodyDecay, ThreeBodyDecay, FastRamboMapping>> _s_decays;
     nested_vector2<me_int_t> _permutations;

madspace/include/madspace/phasespace/t_propagator_mapping.hpp

@@ -12,7 +12,9 @@ namespace madspace {
 class TPropagatorMapping : public Mapping {
 public:
     TPropagatorMapping(
-        const std::vector<std::size_t>& integration_order, double invariant_power = 0.8
+        const std::vector<std::size_t>& integration_order,
+        double invariant_power = 0.8,
+        const std::vector<double>& pt_min = {}
     );
     std::size_t random_dim() const { return 3 * _integration_order.size() - 1; }
 
@@ -28,8 +30,12 @@ class TPropagatorMapping : public Mapping {
         const NamedVector<Value>& conditions
     ) const override;
 
+    // pt^2 of the outgoing particle at position i (0 if no pt cut).
+    double pt2(std::size_t i) const;
+
     std::vector<std::size_t> _integration_order;
     std::vector<bool> _sample_sides;
+    std::vector<double> _pt_min;
     Invariant _uniform_invariant;
     TwoToTwoParticleScattering _com_scattering;
     TwoToTwoParticleScattering _lab_scattering;

madspace/include/madspace/phasespace/three_particle.hpp

@@ -36,6 +36,9 @@ class TwoToThreeParticleScattering : public Mapping {
         double s_width = 0
     );
 
+    // one discrete input: which of the 2 two-body solutions to take
+    std::size_t discrete_dim() const override { return 1; }
+
 private:
     Result build_forward_impl(
         FunctionBuilder& fb,