madspace: evaluate the matrix element on the momenta the event is written with

UpdateNotes.txt

@@ -5,6 +5,23 @@ ANNOUNCEMENT:
    A new LTS (based on the current 3.5.X) is starting now and will act as the stable release for the coming years.
 
 0.3.0 (XX/XX/XX):
+   OM: BUG FIX (mg7): the initial-state mirror paired each event with a matrix
+       element evaluated on the momenta *before* the mirror was applied. madspace
+       covers the beam-swapped copy of an asymmetric initial state by rotating
+       the event by pi about x (py, pz -> -py, -pz), but only the written event
+       was rotated; the weight kept the momenta as generated. The two agree only
+       because |M|^2 is invariant under that rotation -- which a polarised matrix
+       element is not once it is evaluated in a frame that holds the polarised
+       particle at rest, because HELAS then quantises the spin along the frame z
+       axis, the mirror flips it, and the + and - states swap. p p > z{+} j and
+       p p > z{-} j evaluated in the Z rest frame came out 4485 +- 39 pb against
+       6060 +- 54 pb, 23 sigma apart, where the symmetry of a pp collider
+       requires them to be equal; they now agree within 1 sigma. The sum over the
+       three polarisations was already right, which is why no unpolarised result
+       is affected (checked: p p > z j unchanged at 12420 +- 69 pb). The matrix
+       element is now evaluated on the momenta the event is written with, which
+       is both what madevent does and what the systematics reweighting here
+       already assumed, since it reads them back out of the event buffer.
    OM: 'customize_model' now works for every UFO model: generic options (flavour
        scheme, massive leptons, diagonal CKM), commands to fix a parameter or a whole
        block ('set decay all 0'), to change a formula, to choose the values written in

madspace/src/phasespace/integrand.cpp

@@ -637,6 +637,27 @@ NamedVector<Value> Integrand::build_common_part(
     auto flavor_id = args.at("flavor_id");
     auto batch_size_val = fb.batch_size({args.at("weight_before_cuts")});
 
+    // The matrix element is evaluated on the momenta the event is written
+    // with, which for the mirrored half of a beam-swapped subprocess is the
+    // mirrored set (mirror_momenta: py, pz -> -py, -pz, the rotation by pi
+    // about x that moves each leg onto the other beam). This is what madevent
+    // does too -- it flips the momentum array that then goes to both the matrix
+    // element and the event record (super_auto_dsig_group_v4.inc, "Flip momenta
+    // (rotate around x axis)") -- and it is what the systematics reweighting
+    // here already assumes, since that reads the momenta back out of the event
+    // buffer.
+    //
+    // Feeding the unmirrored momenta instead only ever worked because |M|^2 is
+    // invariant under that rotation. A polarised matrix element is not, once it
+    // is evaluated in a frame that holds the polarised particle at rest: HELAS
+    // quantises such a particle along the *frame* z axis (the pp == 0 branch of
+    // vxxxxx) rather than along its own momentum, the mirror flips that axis,
+    // and the + and - states swap. Measured at 48% on g q > z{+} q evaluated in
+    // the Z rest frame.
+    auto momenta_me = args.index_map().contains("momenta_mirror_acc")
+        ? args.at("momenta_mirror_acc")
+        : momenta_acc;
+
     auto scatter_or_drop = [&](Value default_value, Value value) -> Value {
         if (_drop_cuts_and_rescale) {
             return value;
@@ -680,7 +701,7 @@ NamedVector<Value> Integrand::build_common_part(
 
     // Evaluate differential cross section
     ValueVec xs_args{
-        momenta_acc,
+        momenta_me,
         _flavor_remap.size() > 0 ? fb.gather_int(flavor_id, _flavor_remap) : flavor_id,
     };
     xs_args.push_back(x1_acc);