Fix madspace issue 292

madspace/include/madspace/phasespace/color_ordered_mapping.hpp

@@ -1,5 +1,7 @@
 #pragma once
 
+#include <array>
+#include <optional>
 #include <vector>
 
 #include "madspace/phasespace/base.hpp"
@@ -53,25 +55,41 @@ namespace madspace {
 class ColorOrderedMapping : public Mapping {
 public:
     /**
-     * @param color_order       0-based permutation of `{0, ..., n-1}`
-     *                          (`n = n_out + 2`) giving the chain order;
-     *                          particles 0 and 1 are the incoming beams. See
-     *                          @ref Topology.
-     * @param t_invariant_power Exponent of the `1/|t|^p` sampling of every
-     *                          momentum transfer; see @ref Invariant.
-     * @param s_invariant_power Exponent of the `1/s^p` sampling of every
-     *                          time-like invariant; see @ref Invariant.
+     * Settings of the chain. Besides these, the chain follows AmpliCol's
+     * gen23 [2]: a side of exactly two particles (with the other side not
+     * empty) is one 2->3 step with the opposite beam as the previous emission,
+     * the composite and rest masses are sampled as `1/sqrt(s)` (flat in the
+     * mass), and the s23 random numbers of the 2->3 steps are remapped with the
+     * arcsine map of @ref TwoToThreeParticleScattering.
+     */
+    struct Options {
+        /// 0-based permutation of `{0, ..., n-1}` (`n = n_out + 2`) giving the
+        /// chain order; particles 0 and 1 are the incoming beams. See
+        /// @ref Topology.
+        std::vector<std::size_t> color_order;
+        /// Exponent of the `1/|x|^p` sampling of the momentum transfers and
+        /// the s23 invariants of the chain; see @ref Invariant. Unset: 1 with
+        /// cuts (AmpliCol) and 0.8 without, where these ranges reach zero and
+        /// `1/|x|` would pile the points up at the regulator of
+        /// @ref Invariant.
+        std::optional<double> chain_invariant_power = std::nullopt;
+        /// In every 2->3 step, keep the azimuthal solution that passes the pt
+        /// cuts, and pick with the discrete choice only when both pass, as
+        /// AmpliCol does. The weight assumes uniformly sampled discrete
+        /// choices: turn it off when they are sampled with a learned
+        /// distribution.
+        bool cut_aware_branch = true;
+    };
+
+    /**
+     * @param options           Chain order and settings; see @ref Options.
      * @param pt_min            Per-outgoing-particle minimum transverse
      *                          momentum; empty disables the cut. See @ref Cuts.
      * @param m_inv_min         Symmetric matrix of minimum pair invariant
      *                          masses; empty disables the cut. See @ref Cuts.
      * @param dr_min            Symmetric matrix of minimum pair `delta_r`
      *                          separations; empty disables the cut. See
      *                          @ref Cuts.
-     * @param arcsine_s23       Remap the random number of every 2->3
-     *                          invariant with the arcsine map in the azimuth,
-     *                          removing the `1/|sin(phi)|` edge peak of the
-     *                          weight; see @ref TwoToThreeParticleScattering.
      * @param y_max             Per-outgoing-particle bound on the absolute lab
      *                          rapidity, negative (or empty) for none. Used
      *                          where a block scatters the two beams in the
@@ -87,13 +105,10 @@ class ColorOrderedMapping : public Mapping {
      *                          conditions.
      */
     ColorOrderedMapping(
-        const std::vector<std::size_t>& color_order,
-        double t_invariant_power = 0.8,
-        double s_invariant_power = 0.8,
+        const Options& options,
         const std::vector<double>& pt_min = {},
         const std::vector<std::vector<double>>& m_inv_min = {},
         const std::vector<std::vector<double>>& dr_min = {},
-        bool arcsine_s23 = true,
         const std::vector<double>& y_max = {}
     );
 
@@ -125,13 +140,10 @@ class ColorOrderedMapping : public Mapping {
     // The public constructor with the rapidity bounds worked out once.
     ColorOrderedMapping(
         const RapidityBounds& bounds,
-        const std::vector<std::size_t>& color_order,
-        double t_invariant_power,
-        double s_invariant_power,
+        const Options& options,
         const std::vector<double>& pt_min,
         const std::vector<std::vector<double>>& m_inv_min,
-        const std::vector<std::vector<double>>& dr_min,
-        bool arcsine_s23
+        const std::vector<std::vector<double>>& dr_min
     );
 
     Result build_forward_impl(
@@ -151,6 +163,22 @@ class ColorOrderedMapping : public Mapping {
     // 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;
+    // Floor on the squared mass of a set of outgoing particles (masses and
+    // cuts), and floor on its transverse mass: the sum of those of the
+    // particles with their pt cuts, at least the square root of the former.
+    std::array<Value, 2> set_floors(
+        FunctionBuilder& fb, const std::vector<std::size_t>& set, const ValueVec& m_out
+    ) const;
+    // Range of the squared mass of the rest system of a walk over the set s
+    // after its peel j, from a system of mass prev_mass.
+    std::array<Value, 2> rest_mass_range(
+        FunctionBuilder& fb,
+        const std::vector<std::size_t>& s,
+        std::size_t j,
+        Value prev_mass,
+        Value e_cm,
+        const ValueVec& m_out
+    ) const;
 
     // 0-indexed outgoing-particle indices (values in {0,...,n_out-1}).
     // _set1 contains the outgoing particles attached to beam 0's side,
@@ -179,12 +207,14 @@ class ColorOrderedMapping : public Mapping {
     bool _has_cut;
     bool _rapidity_window;
 
-    Invariant _uniform_invariant;
+    Invariant _mass_invariant;
     TwoToTwoParticleScattering _com_scattering;
     TwoToTwoParticleScattering _lab_scattering;
     // first peel of a single chain, whose incoming momenta are the beams
     TwoToTwoParticleScattering _chain_scattering;
     TwoToThreeParticleScattering _two_to_three;
+    // the 2->3 step of a side of two particles, with an incoming recoil
+    TwoToThreeParticleScattering _two_to_three_beam;
     DoubleT _double_t;
 };
 

madspace/include/madspace/phasespace/phasespace.hpp

@@ -80,14 +80,21 @@ class PhaseSpaceMapping : public Mapping {
      * @param leptonic         If true, skip the PDF convolution (no `x1`/`x2`
      *                         sampling).
      * @param invariant_power  Exponent of the `1/s^p` sampling of every
-     *                         time-like invariant. See @ref Invariant.
+     *                         time-like invariant. See @ref Invariant. The
+     *                         chain of @ref TChannelMode::color_ordered takes
+     *                         its powers from @p color_ordered.
      * @param t_channel_mode   The @ref TChannelMode strategy.
      * @param cuts             Fiducial cuts to fold into the sampling. See
      *                         @ref Cuts.
      * @param permutations     One outgoing-particle permutation per channel.
      *                         More than one enables the `permutation_index`
      *                         condition.
-     * @param color_order      Color ordering for @ref TChannelMode::color_ordered.
+     * @param color_ordered    Chain order and settings of
+     *                         @ref TChannelMode::color_ordered, required there
+     *                         and not allowed for the other modes. An unset
+     *                         `chain_invariant_power` is 1 with cuts and
+     *                         @p invariant_power without. See
+     *                         @ref ColorOrderedMapping::Options.
      * @param beam_rapidity    Rapidity of the beams' centre-of-mass frame in the
      *                         lab frame the momenta are returned in:
      *                         0.5 ln(E1 / E2) for beam energies E1 (along +z)
@@ -100,6 +107,14 @@ class PhaseSpaceMapping : public Mapping {
      *                         beams' frame (leg 1 from beam 2, leg 2 from beam 1)
      *                         before the lab boost and the cuts. With identical
      *                         beams that is a rotation of the lab event too.
+     * @param s_hat_power      Exponent of the `1/s_hat^p` sampling of the
+     *                         partonic invariant s_hat in the PDF convolution
+     *                         (see @ref Luminosity), in every mode. Unset, it is
+     *                         1.2 for @ref TChannelMode::color_ordered when
+     *                         masses or cuts bound s_hat from below (as
+     *                         AmpliCol), and @p invariant_power otherwise. A
+     *                         resonant s_hat is sampled with its Breit-Wigner
+     *                         either way.
      */
     PhaseSpaceMapping(
         const Topology& topology,
@@ -109,9 +124,10 @@ class PhaseSpaceMapping : public Mapping {
         TChannelMode t_channel_mode = propagator,
         const std::optional<Cuts>& cuts = std::nullopt,
         const std::vector<std::vector<std::size_t>>& permutations = {},
-        const std::optional<std::vector<std::size_t>>& color_order = std::nullopt,
+        const std::optional<ColorOrderedMapping::Options>& color_ordered = std::nullopt,
         double beam_rapidity = 0.,
-        bool mirror_beams = false
+        bool mirror_beams = false,
+        std::optional<double> s_hat_power = std::nullopt
     );
 
     /**
@@ -121,11 +137,18 @@ class PhaseSpaceMapping : public Mapping {
      * @param cm_energy        Total collision energy.
      * @param leptonic         If true, skip the PDF convolution.
      * @param invariant_power  Exponent of the `1/s^p` invariant sampling. See
-     *                         @ref Invariant.
+     *                         @ref Invariant. The chain of
+     *                         @ref TChannelMode::color_ordered takes its
+     *                         powers from @p color_ordered.
      * @param mode             The @ref TChannelMode strategy (default
      *                         @ref TChannelMode::rambo).
      * @param cuts             Fiducial cuts. See @ref Cuts.
-     * @param color_order      Color ordering for @ref TChannelMode::color_ordered.
+     * @param color_ordered    Chain order and settings of
+     *                         @ref TChannelMode::color_ordered, required there
+     *                         and not allowed for the other modes. An unset
+     *                         `chain_invariant_power` is 1 with cuts and
+     *                         @p invariant_power without. See
+     *                         @ref ColorOrderedMapping::Options.
      * @param beam_rapidity    Rapidity of the beams' centre-of-mass frame in the
      *                         lab frame the momenta are returned in:
      *                         0.5 ln(E1 / E2) for beam energies E1 (along +z)
@@ -138,6 +161,14 @@ class PhaseSpaceMapping : public Mapping {
      *                         beams' frame (leg 1 from beam 2, leg 2 from beam 1)
      *                         before the lab boost and the cuts. With identical
      *                         beams that is a rotation of the lab event too.
+     * @param s_hat_power      Exponent of the `1/s_hat^p` sampling of the
+     *                         partonic invariant s_hat in the PDF convolution
+     *                         (see @ref Luminosity), in every mode. Unset, it is
+     *                         1.2 for @ref TChannelMode::color_ordered when
+     *                         masses or cuts bound s_hat from below (as
+     *                         AmpliCol), and @p invariant_power otherwise. A
+     *                         resonant s_hat is sampled with its Breit-Wigner
+     *                         either way.
      */
     PhaseSpaceMapping(
         const std::vector<double>& external_masses,
@@ -146,9 +177,10 @@ class PhaseSpaceMapping : public Mapping {
         double invariant_power = 0.8,
         TChannelMode mode = rambo,
         const std::optional<Cuts>& cuts = std::nullopt,
-        const std::optional<std::vector<std::size_t>>& color_order = std::nullopt,
+        const std::optional<ColorOrderedMapping::Options>& color_ordered = std::nullopt,
         double beam_rapidity = 0.,
-        bool mirror_beams = false
+        bool mirror_beams = false,
+        std::optional<double> s_hat_power = std::nullopt
     );
 
     /**

madspace/include/madspace/phasespace/three_particle.hpp

@@ -174,6 +174,23 @@ class TwoToThreeParticleScattering : public Mapping {
      *                          holds it more precisely than pa + pb - p3 (a
      *                          soft system next to the beams) keeps that
      *                          precision.
+     * @param incoming_recoil   If true, `momentum3` is an incoming momentum
+     *                          k (the other beam), and the second invariant is
+     *                          the momentum transfer @f$(p_2 - k)^2 \le 0@f$
+     *                          instead of @f$(p_2 + p_3)^2@f$; its absolute
+     *                          value is sampled. With @p has_cut it is bounded
+     *                          by @f$|(p_2 - k)^2| \ge \texttt{etmin\_2}^2@f$,
+     *                          the `drcut` and `s23_min_cut` conditions are
+     *                          not used. Requires @p p12_condition. AmpliCol
+     *                          builds a side of two particles this way, with
+     *                          the opposite beam as the previous emission.
+     * @param cut_aware_branch  With @p has_cut, form both azimuthal solutions
+     *                          and keep the one passing the transverse-energy
+     *                          floors `etmin_1`, `etmin_2` (AmpliCol's
+     *                          gen23_one_step): `discrete_choice` only picks
+     *                          when both pass. The weight assumes
+     *                          `discrete_choice` is sampled uniformly; with a
+     *                          learned distribution of it, leave this off.
      */
     TwoToThreeParticleScattering(
         double t_invariant_power = 0,
@@ -184,7 +201,9 @@ class TwoToThreeParticleScattering : public Mapping {
         double s_width = 0,
         bool has_cut = false,
         bool arcsine_s23 = true,
-        bool p12_condition = false
+        bool p12_condition = false,
+        bool incoming_recoil = false,
+        bool cut_aware_branch = false
     );
 
     /// Number of discrete inputs (1): which of the two two-body solutions.
@@ -204,12 +223,23 @@ class TwoToThreeParticleScattering : public Mapping {
 
     std::array<Value, 3>
     split_conditions(FunctionBuilder& fb, const NamedVector<Value>& conditions) const;
+    // Sampled and kinematic range of the second invariant: s23 itself, or
+    // |(p2 - k)^2| for an incoming recoil. Arguments: the kinematic range of
+    // s23 and the conditions.
+    std::array<Value, 4> sampling_ranges(
+        FunctionBuilder& fb,
+        Value s23_min,
+        Value s23_max,
+        const NamedVector<Value>& conditions
+    ) const;
     Invariant _t_invariant;
     Invariant _s_invariant;
     double _s_power, _s_mass;
     bool _has_cut;
     bool _arcsine_s23;
     bool _p12_condition;
+    bool _incoming_recoil;
+    bool _cut_aware_branch;
 };
 
 } // namespace madspace

madspace/instruction_set.yaml

@@ -1212,6 +1212,53 @@ two_to_three_particle_scattering_inverse:
       desc: "Jacobian of the mapping"
   desc: "Inverse of `two_to_three_particle_scattering`. Building block of @ref TwoToThreeParticleScattering."
 
+two_to_three_branch_choice:
+  inputs:
+    - name: phi_choice_in
+      type: [int]
+      desc: "sampled branch selector, used when both solutions pass"
+    - name: pa
+      type: [float, 4]
+      desc: "first incoming momentum"
+    - name: p12
+      type: [float, 4]
+      desc: "momentum of the system of the first and second outgoing particles"
+    - name: p3
+      type: [float, 4]
+      desc: "momentum of the already-generated third particle"
+    - name: u
+      type: [float]
+      desc: "position of s23 in its kinematic range"
+    - name: u_c
+      type: [float]
+      desc: "1 - u"
+    - name: t1_abs
+      type: [float]
+      desc: "absolute value of the first momentum transfer"
+    - name: m1
+      type: [float]
+      desc: "mass of the first outgoing particle"
+    - name: m2
+      type: [float]
+      desc: "mass of the second outgoing particle"
+    - name: etmin_1
+      type: [float]
+      desc: "transverse-energy floor of the first outgoing momentum"
+    - name: etmin_2
+      type: [float]
+      desc: "transverse-energy floor of the second outgoing momentum"
+  outputs:
+    - name: phi_choice
+      type: [int]
+      desc: "branch selector to use: a passing solution"
+    - name: factor
+      type: [float]
+      desc: "weight factor for a uniformly sampled selector: 1 if both solutions pass, 1/2 if one does, 0 if none"
+    - name: inv_factor
+      type: [float]
+      desc: "inverse of the factor, 0 if no solution passes"
+  desc: "Choice of the azimuthal solution of the 2 to 3 step that passes the transverse-energy cuts, as in AmpliCol. Building block of @ref TwoToThreeParticleScattering."
+
 double_t_scattering:
   inputs:
     - name: r_phi

madspace/src/compgraphs/function.cpp

@@ -348,11 +348,16 @@ FunctionBuilder::instruction(InstructionPtr instruction, const ValueVec& args) {
         auto& const_arg = args.at(1 - var_index);
         switch (opcode) {
         case opcodes::add:
-        case opcodes::sub:
             if (std::get<double>(const_arg.literal_value) == 0.) {
                 return {var_arg};
             }
             break;
+        case opcodes::sub:
+            // x - 0 = x, but 0 - x is not x
+            if (var_index == 0 && std::get<double>(const_arg.literal_value) == 0.) {
+                return {var_arg};
+            }
+            break;
         case opcodes::mul:
             if (std::get<double>(const_arg.literal_value) == 1.) {
                 return {var_arg};