Repository navigation
Crossing: integrate the helicity PRs (#278 #279 #263 #255 #261 #217), MadSpin with crossing - #294
Open
oliviermattelaer wants to merge 139 commits into
Open
oliviermattelaer wants to merge 139 commits into
oliviermattelaer wants to merge 139 commits into
Conversation
Port of the Fortran option to the C++ backend. helas.h gets vxxxxxr(momenta, vmass, nhel, nsv, flv, vc, ipar, iref): vxxxxx followed by eps -> eps - (r.eps)/(r.p) p with r the momentum of particle iref, skipped for a massive vector, for |nhel| != 1, for iref == ipar and in the lanes where r.p = 0 (fpternary/cxternary on SIMD, as vxxxxx does). The writer emits it for every massless external vector, with the reference of helas_call_writers.get_axial_gauge_refs, when the output line has --axial_gauge=True (standalone and mg7). Validated on g g > g g g, u u~ > g g g, e+ e- > u u~ g g and g g > t t~ g against the default gauge, on common momenta, scalar and simd_128: FPTYPE=d median 7e-16 (max 1.5e-11 on a near-collinear RAMBO point), f and m equally consistent. g g > t t~ g differs at 1.6e-7 with the physical top width (the standalone output keeps the t-channel width, which breaks gauge invariance) and agrees to 2e-16 at WT = 0. No gain yet: madmatrix computes every amplitude of every helicity, so the vanishing amplitudes are still evaluated and the shift costs 3-4% (g g > g g g 84.3k vs 88.1k events/s). The gain needs helicity recycling with per-helicity zero-amplitude removal, as in madevent. Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
Port of the dialect half of a1c1120b4 / b62f23c05 (claude/helicity-recycling- madmatrix-db046c) onto main's hel_recycle, without the PR #49 changes those commits sat on (chunked -O0 subroutines, the AMP gather, C-parity, the DAG compatibility mask). What is fortran-specific in the recycler is only its text layer: how a call is found and taken apart, which argument is the helicity, the leg, the result, how a slot is named, how a line continues. That is now Dialect; everything between parse and write (the DAG, the unfolding, the good-helicity and dead-amplitude filters, the amplitude split plan) is untouched and shared. CppDialect reads what madmatrix writes: a template argument list before the real one (and m_pars->MT), a result that is not the last argument, amplitudes that all write one scratch and so carry their diagram in a //HRAMP tag, and the per-row colour fold (//HRFOLD blocks with a ${row} hole) stamped out once per helicity row. The P1N split is planned by plan_amp_split, shared by the fortran split_amps and the C++ renderer. New here: VXXXXXR is in the external routines of both dialects, and the C++ leg is the argument after the output slot, since vxxxxxr ends with its reference leg rather than its own. Fortran output is unchanged: replaying gen_ximprove's recycler on the 24 matrix1_orig.f of 12 madevent outputs (gluon, quark, e+e- and top processes, default and axial gauge) gives the same matrix1_optim.f byte for byte as main's module, including the two that differ from the files on disk (a pre-existing NWAVEFUNCS leak between recyclers, fixed separately). Tests: TestCppDialect, TestCppAmplitudeSplit (from db046c) plus the vxxxxxr leg. Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
…state It is read back as NWAVEFUNCS but was never reset, so a matrix element recycled after a bigger one in the same process -- gen_ximprove loops over every P directory in one -- declared the bigger one's wavefunction array: e+ e- > u u~ after e+ e- > u u~ g g g got NWAVEFUNCS=14 instead of 12. Harmless (an oversized array) but it made the output depend on what was recycled before it. Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
…ty zero amplitudes Port of the madmatrix side of a1c1120b4 / b62f23c05 (db046c) onto main's backend layout, through the C++ dialect of hel_recycle (77a83b1). backend/{cpu,simd}/SigmaKin.cc gain, CPU only: - calculate_all_jamps: every wavefunction once per COPY, the amplitudes the warm-up found vanishing at a row left out, the others folded into per-row tables (jamps, multichannel numerators); - load_jamps: one row handed back to the helicity loop, with the bookkeeping calculate_jamps does for its own row (numerators, jamp2 for the color choice). Both the BLAS and the per-helicity loop call it. They are compiled only when the P directory has HelRecycling.inc (__has_include), so an output without --hel_recycling, or whose warm-up failed, builds exactly as before. The page setup both functions share moved from calculate_jamps to JampsPreamble.h (a fragment, not a header). Zero amplitudes are found the way madevent finds them, not by a second implementation: a probe build (make BACKEND=scalar FPTYPE=d HRPROBE=1) has every amplitude record |amp|^2, and an amplitude is alive at a helicity once it is above 1e-12 of that helicity's largest at the same point -- a gauge cancellation is a rounding residual, never 0. sigmaKin_setGoodHel prints the good helicities and the dead (helicity, amplitude) pairs under MG_DUMP_GOODHEL in madevent's protocol; ProcessExporterMadMatrixStandalone. finalize builds and runs it, and write_recycled hands them to the recycler. The rows are those of the probe's own sigmaKin_getGoodHel, so the table can never be short of one the runtime asks for (at these parameters; hr_unbaked_row warns otherwise). The recycler input (hr_orig: one call per line, amplitudes tagged //HRAMP, their per-row fold in //HRFOLD blocks) is written BEFORE the ordinary calls, while the amplitude numbers are still intact, and the numbers are restored after -- the ordinary pass renumbers every amplitude to the one scratch, which is what forced db046c to key a second pass on positions. Flavor-mask guards, shared color sub-expressions and the divergence-test variant of an external are left out of hr_orig (multi-line constructs the reader cannot unfold). Refused, with a warning, for the FD gauge, split amplitude orders and directories with more than one matrix element. Also: do_output asks ALOHA for the P1N routines for this output too, and the madmatrix ALOHA writer recognises their identity denominator (the integer 1, which it compared to the string '1', then wrote '/ ( 1 )', which a cxsmpl<double> cannot be divided by). g g > g g g, FPTYPE=d, 2048 points vs plain: scalar and simd_128 median 2.2e-16, max 9.7e-14; with --axial_gauge max 1.5e-11 (the near-collinear point plain axial also has). simd_128 throughput 86.5k -> 237.8k evt/s (2.75x), 249.7k with --axial_gauge (2.89x). Flag off: every generated process file identical to before. Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
One function holding the whole unrolled block is what makes a big process unbuildable: the compiler's per-function cost is superlinear. g g > g g g g recycles into 136510 statements, which build in 808 s inline and in 29 s cut into noinline hr_chunk_N functions of 200 statements that calculate_all_jamps calls in order (52 s at 1000, 76 s at 3000, for the same run time; inline runs 13% faster, 14.8k vs 13.1k events/s, not worth 13 minutes). The locals of calculate_all_jamps they read are passed in (HR_CHUNK_STATE); the model tables are file-static in the same translation unit. A cut only falls at brace depth zero after a finished statement -- `if( storeChannelWeights )` has its body on the next line (db046c's rule and its tests). Below 15000 statements the block stays inline: chunks cost ~9% at run time (g g > g g g: 239k -> 218k events/s, the same with 8 chunks or 38) and its 8148 statements build in 6 s inline. --hel_recycling_chunk=<n> overrides, 0 keeps any block inline. Measured (simd_128, FPTYPE=d, |M|^2 vs plain on 2048 points): g g > g g g g, axial gauge 2.20x, max 1.3e-12 g g > t t~ g g 1.70x (1.76x axial), max 7e-15 Tests: TestRecycledChunkBoundaries (from db046c), TestRecycledChunkPolicy, TestHelRecyclingWarmup. Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
…as madevent sigmaKin_getGoodHel kept every helicity whose |M|^2 was not exactly 0 at one of the sampled events. A helicity whose amplitudes only cancel to rounding is never exactly 0, so it was kept and evaluated at every event for nothing: the all-plus and one-minus helicities of a gluonic process, 12 of the 32 of g g > g g g. madevent's filter (matrix_madevent_group_v4.inc) keeps a helicity when, at one event, |M|^2 > limhel/ncomb * sum over helicities, limhel being its hidden run_card parameter (1e-8 by default); this is that rule, with the per-event sum over the helicities of the same flavor (the union over flavors and events is unchanged). CPU, SIMD and GPU backends; the GPU one is inspected, not compiled here (no CUDA/HIP on this machine). g g > g g g now keeps 20 rows, and the --hel_recycling warm-up reports exactly madevent's counts (108 vanishing (helicity, amplitude) pairs, 378 in the axial gauge; u u~ > g g g 12 rows and 84). |M|^2 is bit-identical to before on 2048 points for g g > g g g, g g > t t~ (g), p p > e+ e-, u u~ > d d~ g, e+ e- > u u~ g g and u u~ > g g g: what was dropped was below one ulp. simd_128 FPTYPE=d, events/s relative to the previous plain build: g g > g g g plain 1.61x, --hel_recycling 3.13x, + --axial_gauge 3.84x u u~ > g g g plain 1.31x, --hel_recycling 2.29x, + --axial_gauge 2.96x e+ e- > u u~ g g unchanged (its dropped helicities were exactly 0 already) Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
…/madmatrix-axial-hel-recycling # Conflicts: # tests/unit_tests/iolibs/test_helas_call_writers.py
…reference The C++ side of be2519f (#216): a massive vector, which get_axial_gauge_refs gives no reference since the shift needs p^2 = 0, was written as vxxxxxr with itself as reference -- a no-op, but p p > w+ w+ j j QCD=0, which has nothing to gauge, came out different from the plain output. Now its EvaluateDiagrams.inc is identical to the plain one in all three P directories, and g g > g g g is unchanged. Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
Two review findings on #217. The warm-up ran `check_sa.exe matrix`, which feeds nevt copies of ONE phase-space point, so the baked helicity rows and the dead (helicity, amplitude) pairs rested on a single point. The run-time scan samples 16 distinct points of the real sample, and a row it finds good that was not baked contributes nothing (hr_unbaked_row only warns). The warm-up now runs `check_sa.exe perf 1 64 1`, and the probe build (-DMG_HR_PROBE) samples 64 events instead of 16 and keeps every row above 1e-14 of the sum instead of 1e-8: a row baked for nothing only costs time, a missing one loses physics, and a rounding residual (~1e-30) is still far below. Counts unchanged on g g > g g g (20 rows, 108 / 378 axial), u u~ > g g g, g g > t t~ g and p p > e+ e-; |M|^2 vs plain as before and no unbaked row in any run. The generated diagram calls of a --hel_recycling output carry MG_HR_PROBE_AMP( n ) after every amplitude, which only the cpu and simd SigmaKin.cc defined: a cuda/hip build of such an output did not compile. The GPU SigmaKin.cc now defines it empty (recycling and the probe are CPU only). Not compiled here (no CUDA/HIP on this machine). Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
Review finding on #217. HelicityRecycler kept bad_amps and bad_amps_perhel as the lists its callers pass in (gen_ximprove and madmatrix's write_recycled both pass sorted lists), and good_helicity asks `(hel, amp) in bad_hel_amp` once per unfolded copy of every amplitude: a linear scan each time. Both are only ever tested for membership, so they are now sets, built once in __init__, which speeds up the madevent and the madmatrix recycling alike. madmatrix g g > g g g g --hel_recycling --axial_gauge (9956 dead pairs): the whole output, warm-up build included, 15.4 s -> 7.2 s, with HelRecycling.inc, HelRecyclingCalls.inc and EvaluateDiagrams.inc byte-identical. Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
…e per session Review finding on #217. axial_gauge_requested kept a module-level `warned` flag that was never reset, so in one MG5 session only the first output with --axial_gauge in the FD gauge said the option was ignored; every later one dropped it silently. axial_gauge_requested is now a pure check, asked per matrix element as before, and warns only when called with warn=True, which each exporter does once when it is built: ProcessExporterFortran (standalone_fortran and, through super(), madevent) and ProcessExporterMadMatrix. In the FD gauge madevent no longer prints the sde_strategy advice for an option it is not using. Checked in one session: set gauge FD, then output standalone and output madevent with --axial_gauge (one warning each) and output madevent without it (none); without FD, madevent still writes VXXXXXR with the sde_strategy note. Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
…l-hel-recycling # Conflicts: # madgraph/madevent/hel_recycle.py # tests/unit_tests/madevent/test_hel_recycle.py
…l-hel-recycling # Conflicts: # madmatrix/model_handling.py
From claude/helrecycle-dag-memo (Aug 2026), which never reached main; ported onto the Dialect version of the recycler. All output-byte-identical. - DAG.find_path was a recursive DFS asked once per (dependency, external) pair -- 114 million calls on g g > g g g g g. The graph is two levels deep by construction (store_wav's edges go straight to the externals, its caller's ext_deps being a transitive closure already), so reachability is a set membership: ext_closure per node, unioned; self.graph goes. - dependencies / old_names / kill_old scanned all_wavs (every wavefunction ever built, dead ones included): bucketed by old_name. kill_old empties its bucket. - good_helicity's "covered by a good combination" rebuilt set(comb) per comb per call: the combinations are bitmasks, ANDed per dependency and cached (set_good_wav_combs, rebuilt in get_gwc where they change). - get_arguments (the fortran dialect's parser) is asked again for every object unfolded out of a line: a 32-slot cache, handing out copies. - all_hel.index -> External.all_hel_index, built in set_helicity_table. The split_amps bitmask selection (a86fb4e5d) came with plan_amp_split already, and the bad-amplitude sets (2c2e1a9) with 432eeca. Recycler alone, replayed on madevent inputs, old -> new, matrix1_optim.f byte-identical: g g > g g g g g 53.0 s -> 6.2 s (22 MB, 47700 dead pairs), g g > g g g g axial 0.85 -> 0.21 s, g g > t t~ g g 0.30 -> 0.08 s. madmatrix --hel_recycling outputs (g g > g g g, u u~ > g g g axial, g g > t t~ g, p p > e+ e-, g g > g g g g axial): HelRecycling*.inc byte-identical. Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
… page
umami_matrix_element (cpu and simd backends) rounded every batch up to two
SIMD pages ("for some reason"), which only mixed precision on SIMD needs
(FPTYPE=m/e, two neppV pages per kernel step). A one-event call to the
scalar backend therefore computed two events, and simd_128 two pages.
The padding events were left at zero (momenta and g_s), which evaluate to
nan. The good-helicity filter runs at the first call and keeps a helicity
if its |M|^2 is not 0 for one of the first events: a padding lane among
them kept every helicity for the rest of the run, the vanishing ones of a
q q~ initial state included. The padding lanes are now copies of the first
event of their page (as the GPU backend already does), so the filter only
sees real kinematics; the two fill loops (flavor-sorted or not) become one.
One |M|^2 per call (Apple M5 Pro, 1 thread, FPTYPE=d, us per point):
before after (f2py -O3)
scalar g g > t t~ 4.09 2.29 2.62
scalar u u~ > t t~ 2.61 0.97 0.80
scalar g g > t t~ g g 324 174 223
scalar u u~ > t t~ g g 110 28.9 36.1
simd128 u u~ > t t~ g g 123 32.9
A batched run whose first call had a single event: u u~ > t t~ g g
simd_128 31.2 -> 16.2 us. |M|^2 unchanged (same values to 1e-15).
madmatrix_good_helicity_count() (extern "C", cpu/simd, not part of UMAMI)
reports how many helicities the filter kept. The acceptance test
test_standalone_umami_padding uses it on p p > e+ e- (scalar/d,
simd_128/d, simd_128/m): a one-event first call must keep the helicities a
16-event first call keeps (4, not 16; it fails with 16 != 4 without the
fix), and a batch with alternating flavors and partial pages must give the
|M|^2 of single-event calls.
The density spin modes (madspin/full, PA, onshell) evaluated their
production and decay density matrices through the f2py modules of an
`output standalone_fortran --density=1`. They now use the libraries of a
madmatrix standalone output.
madmatrix:
- backend/common/Density.{h,cc} (C++ backends): madmatrix_density, the
port of the Fortran GET_DENSITY (same POS/ALLOW_HEL convention, divided by
the same iden, beam polarisation folded in, split orders paired as the
color sum does), madmatrix_matrix_element (the SMATRIX normalisation) and
madmatrix_set_beampol. Helicity rows that vanish on the first 20 events of
a (flavor, open legs) are skipped, the Fortran SMATRIX filter.
- calculate_jamps takes the helicity row instead of its index (a local cHel
shadows the table, EvaluateDiagrams.inc is unchanged): GET_DENSITY builds
combinations a polarised process leaves out of its helicity table.
check_sa perf of g g > t t~ g g is unchanged, scalar and simd_128.
- An external leg marked off-shell ('*', MadSpin's decaying particles) is
evaluated at the invariant mass of its momentum, as the Fortran output
does (offshell_mass). It used the pole mass: off shell the top density was
2-4% off for one helicity and MadSpin's maximum weights ~8% low. Exact in
the scalar backend; SIMD/GPU keep the pole mass (one mass per call), the
same for the on-shell momenta mg7 generates.
- make LIBTAG=...: goes into every library name. Two outputs in one process
shared libmadmatrix_common_<backend>.so by name, and the second silently
ran with the first one's parameters (decay |M|^2 = 0). make PGOAL=lib
builds the libraries without check_sa.exe.
- The standalone exporter writes SubProcesses/madmatrix_processes.json
(library, legs, processes, off-shell legs, every pdg combination of each
flavor index).
- madmatrix/density.py: ctypes binding with the f2py module's method names,
pdg -> (library, flavor index) dispatch with the FLAVOR_ORDER_REPAIR leg
reordering, refuses a library name already loaded from another path and a
non-scalar build of a process with off-shell legs.
MadSpin: option density_backend (madmatrix by default, fortran). A
loop-induced production takes fortran (MadLoop), and a madmatrix output
that fails falls back to it.
Validation: densities and |M|^2 agree with the f2py ones to 1e-15 on shell,
1e-13 off shell (p p > t* t~*, t* > b w+, w+ > l+ vl, merged flavors); with
decay_generator madevent the decayed events of fortran and madmatrix are
byte-identical. p p > t t~, both tops leptonic, 10k events: wall time
66.0 -> 58.8 s (1 core), 31.3 -> 23.5 s (6 cores), from the matrix-element
generation and compilation (~13 -> 6.3 s); per point the density kernel is
2.3x faster.
Tests: tests/unit_tests/iolibs/test_madmatrix_density.py,
test_madmatrix_density_matches_fortran (tests/acceptance_tests/test_cmd.py).
…ecoil floor Measured on p p > t t~, both tops leptonic, 10^4 events, spinmode madspin, reproducible decay pools (decay_generator madevent), 1 core; gain = time of max-weight scan + decay loop before / after, outputs compared event by event. - Event.copy() (lhe_parser): what Event(str(event)) gives -- the LHE fields of the particles and of the event -- without rendering and re-parsing the text, which MadSpin did per trial (production copy, decay copy, reshuffling probes). Momenta keep their doubles instead of 11 digits. 13.6 us instead of 45 us per copy; gain 1.06 joint, 1.13 sequential, 1.11 with 145 weights. - Event.__str__: blank lines are collapsed in the free-text parts only (the same text, byte for byte), not by a regular expression over the particles and every weight: 21 us instead of 38 us; gain 1.08 with 145 weights. - Joint offshell: the density basis (get_iden, get_pdir, leg positions) and the on-shell |M_prod|^2 depend on the production event only, and are worked out once per event instead of per trial (_hoisted_production). Gain 1.06. - project_massless_partons is skipped where there is nothing to project: only aMC@NLO puts massless partons on a Monte-Carlo mass shell, so not for an LO production (decided from the banner's run_card) nor for a decay, which MadSpin generated itself. It triggered on any m^2 != 0 and rescaled rounding errors for ~40 us per density matrix. Gain 1.16 joint, 1.18 sequential. - Decay densities are evaluated on the decay at rest, the madmatrix library boosting it along its parent (madmatrix_density_boosted); add_decays boosts the accepted decays, so rejected trials no longer boost in python. Gain 1.01-1.05 (the boost itself is 2 us). Old libraries keep the python boost. - mass_shuffle: Newton iteration inlined (newtonmethod evaluated f twice per step), E^2 - m^2 computed once, three unused sums dropped; bit-identical, 8.6 us instead of 14 us (a failed solve, a vanishing derivative or an overflow still fails the reshuffling, as around newtonmethod). calculate_matrix_element takes the single get_momenta assignment when no merged pdg repeats in the final state (get_all_momenta redid the mapping to find nothing to permute) and reads its option once. - Correctness, joint offshell max-weight scan: draws 2..400 of a production event reused the production density of the first draw, while every draw reshuffles the production to its own mass set and the decay loop evaluates it on every trial -- the bound was measured on another weight. Now evaluated per draw: overweight events 3 -> 1 (largest factor 3.05 -> 1.06) for 1.6% more trials. - offshell_min_recoil (madspin/full, default 0.5; the study of 2026-09-30, left uncommitted then): evaluated per draw, the bound of the previous item is the true one, and with a massless particle recoiling against the resonances the true one is unbounded -- on u u~ > z g a Z drawn up to sqrt(shat) pushes the gluon soft, the weight grows like 1/chi (RAMBO's momentum scale factor) and its mean over the Breit-Wigner diverges. The scan fix alone took test_madspin_mixed_flavor_decay_log_summary from 16 to 161-505 trials/event. Below chi = 0.5 the production is not reshuffled: a zero weight in the joint accept/reject, a mass-set redraw in the sequential one. Only when a massless final-state particle recoils: t t~ and the other processes without one are unchanged. mass_shuffle(return_chi=True), reshuffle_production(min_recoil=). Together (all of the above), 1f9d40b against this commit, interleaved, medians of 2-3 runs: - madevent pools, 1 core: scan + loop 38.5 -> 27.4 s sequential (whole run 69.8 -> 59.0 s); joint 32.8 -> 28.3 s, of which one more pool refill (+4.7 s: the scan fix's 1.6% more trials cross a refill boundary), so 27.4 -> 18.1 s without the refills (whole run 63.4 -> 59.6 s). - mg7 pools (the default), 6 cores: scan + loop 7.0 -> 5.0 s joint (auto), 7.3 -> 5.5 s sequential; wall 16.5 -> 15.6 s and 17.2 -> 14.5 s (these pools are not reproducible: trial counts differ by up to 4%). The decayed events are those of the previous code (momenta within 3e-12), except where the scan fix changes the bound. Tests: TestEventCopyAndText (test_lhe_parser), TestDecayLoopShortcuts, TestJointScanProductionDensity, TestOffshellRecoilFloor (test_madspin), the boosted density in test_madmatrix_density_matches_fortran; two MadSpin tests asserted the text round trip itself and now assert the copy. Unit suites madspin (678), various (754), iolibs (192) and the 15 MadSpin acceptance tests (aMC@NLO input, gridpacks and use_old_dir included) pass; test_madspin_mixed_flavor_decay_log_summary runs at 11.1-13.9 trials/event.
…BO points Review suggestion (PR #261): the helicity filter ran at the first umami_matrix_element call, so the helicities kept for the whole run depended on what that call held. Padding was one way to get it wrong (the previous commit); special kinematics is another, and worse: a first call with one u u~ > e+ e- event in the exact forward direction, where the (1 - cos theta) helicities vanish, dropped them, and every later |M|^2 was wrong by factors 5 to 100 (5.71e-4 -> 1.18e-4, 2.09e-3 -> 1.75e-5). umami_initialize (cpu and simd backends) now evaluates the matrix element once, per library, on 16 RAMBO points (sqrt(s) = max(1.5 TeV, 2 x the sum of the masses); a decay at rest), which runs the filter there; it already scans every flavor. RAMBO draws from a private generator so that the points of check_sa (rambo::rn) do not move: rambo() and get_momenta() take the random-number source, the old signatures forward rn. rambo exits on a closed phase space, so 2 > 1 and decays below threshold are checked first and keep the first-call choice, as does a point that is not on shell. The GPU backend still chooses at the first call (now on real events). test_standalone_umami_padding becomes test_standalone_umami_first_call: a first call of one generic event, of 16 events, or of one forward event must keep the same helicities and give the same |M|^2 afterwards (scalar/d, simd_128/d, simd_128/m). Without this commit: 2 != 4 helicities for the forward first call. Also: rambo.h's rn() silenced its unused argument with a self-assignment (-Wself-assign, now also printed by the library build); (void) instead.
The cpu and simd backends choose the good helicities in umami_initialize, on 16 RAMBO points (previous commit); the GPU backend still chose them at the first umami_matrix_element call, so a first call of special points (one u u~ > e+ e- event in the exact forward direction) dropped the helicities vanishing there for the rest of the run. backend/gpu/umami.cc now does the same as the cpu and simd backends, with device buffers: the RAMBO points (same private generator, same 2 > 1, threshold and on-shell guards) are made on the host, copied to the device, and umami_matrix_element runs once on them on the default stream (no UMAMI_OUT_GPU_STREAM), which the filter's own copies use too; the stream is synchronized before the buffers are freed. It runs once per library, on the device current in umami_initialize, whose parameters initProc has just set (madspace activates the device first). The first-call path stays as the fallback where RAMBO has no phase space. madmatrix_good_helicity_count() (extern "C", not part of UMAMI) and sigmaKin_getNGoodHel() are added to the GPU backend for parity (cNGoodHel, the host copy of dcNGoodHel). Its comment, in all three backends, now says where the filter runs. test_standalone_umami_first_call_gpu runs the checks of test_standalone_umami_first_call (now in _check_umami_first_call) for BACKEND=cuda and/or hip, FPTYPE=d and m, with device buffers allocated through the CUDA/HIP runtime from the helper. Skipped without nvcc and hipcc on the PATH, or without a device.
The GPU CI ran mg7 (pp_ttx, pp_jj) and the interference checks, but no acceptance test: test_standalone_umami_first_call_gpu (previous commit), which skips without nvcc/hipcc or a device, ran nowhere. New job acceptance in gpu_runner_ci.yml (after build_madspace, for its venv; no madspace or PDF set used), with .github/gpu_ci/acceptance_gpu.sh: the environment of the other GPU jobs (modules, venv, GPU_ARCH, default mg7 configuration), MAKEFLAGS=-j<cpus> for the tests' own make calls, then tests/test_manager.py -p A -t0 for each test of DEFAULT_TESTS (or GPU_TESTS), now test_standalone_umami_first_call_gpu. test_manager exits 1 on a skip, so on a GPU node a skipped test fails the job; summary.txt gives, per test, ok/FAILED, unittest's "Ran N tests" and verdict, and the skip reason (the runner's verdict says OK for a skip). stop_runner waits for it; README row added. Checked locally with a fake nvidia-smi: the skip on a Mac gives exit 1 and "test_standalone_umami_first_call_gpu=FAILED (Ran 1 test in 0.177s, OK, skipped 'neither nvcc nor hipcc on the PATH')". _check_umami_first_call puts the end of the make output in its failure message: on the GPU CI the build log is otherwise lost.
…the Fortran Two failures of the CI of PR #263, both from the madmatrix density backend: - madmatrix/density.py imported numpy at module level, and MadSpin imports it whatever its spinmode. test_manager imports every test module to find the ones it runs, so every job without numpy (all but the MadSpin unit job) died importing tests/unit_tests/iolibs/test_madmatrix_density.py: IOtest_bundle_1, IOtest_loop_sqso, unittest_bundle_1/2/5, unittest_new_coverage and unittest_write_model. And MadSpin itself needed numpy just to be imported, where main imports it inside the functions that use it. numpy is now imported when present: MadDensityLibrary refuses to load without it, and the density unit tests skip. - test_complex_mass_scheme (acceptancetest_madevent_bundle_6): the parameters madmatrix writes in the complex-mass scheme do not compile, and the failure only comes at the build, which the fallback of output_density_me cannot catch. density_backend now takes the Fortran for a complex-mass sample, as for a loop-induced production. Tests: TestDensityBackendChoice. With numpy blocked, the iolibs density tests skip and test_manager collects every unit test module (the unittest_write_model command). test_complex_mass_scheme and test_madmatrix_density_matches_fortran pass.
…l-hel-recycling # Conflicts: # madgraph/iolibs/template_files/madmatrix/backend/cpu/SigmaKin.cc # madgraph/iolibs/template_files/madmatrix/backend/simd/SigmaKin.cc
…xes, faster trials Follow-up of the review of PR #263. - madmatrix (model_handling): computeDependentCouplings_fromG recomputed, event by event, every aS-dependent parameter of the model, not only those the dependent couplings of the process need. With loop_sm(-no_b_mass) that is every R2/UV counterterm parameter, written with COND, reglog and mdl_complexi, which the generated code does not define: the output did not compile, and MadSpin, whose density backend is madmatrix, aborted on every aMC@NLO sample whose banner imports a loop_sm model (e.g. tests/input_files/ttbar_nlo.lhe.gz). Only the parameters the couplings need, directly or through one another, are written now (for g g > t t~ in loop_sm: none, GC_4 = -G and GC_5 = cI * G). That sample decays with madmatrix to the same events as with the Fortran density matrices. - MadSpin output_density_me builds the madmatrix libraries right after writing them, while the processes are still the current ones, and falls back to the Fortran output when the build fails (it only caught a failing output, and the build in compile() had no fallback). - offshell_min_recoil: a recoiling particle is massless when the model says so, not only when its LHE mass is zero. aMC@NLO writes its partons on the shower's Monte-Carlo mass shell (0.32 GeV quarks, a 0.75 GeV gluon with HERWIG6): those recoils were never floored, and their weight grew ~1/chi. - density_debug: a trial the recoil floor refuses (weight zero) is not compared with the full matrix element; the comparison divided by its zero. - ms_dir/use_old_dir: the directory is stamped with the MadSpin cache format and the MadGraph version that built it (ms_version.dat), and one holding cached material with another stamp, or none, is refused (MadSpinStaleDirectory): its maximum weights were bounded by another code (this PR changed the offshell scan and added the recoil floor) and its libraries may not be the ones this version loads. - CI: tests/unit_tests/iolibs/test_madmatrix_density.py (unittest_bundle_3) and test_madmatrix_density_matches_fortran (acceptancetest_32). - Density.o is compiled only in the builds madmatrix/density.py asks for (DENSITY=1, the scalar libraries MadSpin loads), not in every library of every madmatrix output (mg7, check_sa). - Speed (t t~ leptonic, madevent pools, 1 core, scan + decay loop): joint 19.1 -> 16.2 s (x1.18), sequential 20.7 -> 19.0 s (x1.09), medians of 3 interleaved runs, the same trials and byte-identical decayed events: * get_density: the flavour structure of an event (matrix-element order, tag, where each leg sits in the event, the pdgs to pass) is cached per (pid, status) pattern, shared by the trial copies of a production event and the decays of a channel; get_pdir, get_mapping and get_pdg ran on every call. * get_onshell_evt_and_wgt checks the production tag once per production event, and takes sqrt(shat) from the source of a trial copy. * mass_shuffle skips its two Lorentz maps when the momenta already sum to rest (every decay reshuffle): bit-identical, 8.4 -> 5.7 us per call. Tests: TestDensityOutputFallback, TestReusedDirectoryVersion, the shower-mass case in TestOffshellRecoilFloor; test_make_variables expects DENSITY=1.
…o on a redraw - reshuffle_production applies min_recoil only when one of the particles it rescales -- the top-level final state of split_event_by_onshell_propagator, without a virtuality of its own -- is massless by recoil_massless (MadSpin: a zero LHE mass, or massless in the model). _offshell_min_recoil looked at every status-1 line, so the massless decay products of a resonance already decayed in the event (t > b w+, w+ > u d~ written by madevent), which move rigidly with it and never go soft, switched the floor on: with the top decayed in the event and the t~ by MadSpin, 0.54% of the feasible t~ masses -- the upper Breit-Wigner tail near threshold, with bounded weights -- were refused. - The retry of reshuffle_production after an infeasible mass set dropped min_recoil, so the redrawn set skipped the floor (reachable with a massive particle outside the mass budget, e.g. an undecayed one): it now keeps the floor and the masslessness test. - _offshell_min_recoil() returns the option and _recoil_massless(particle) the test. Tests: test_floor_ignores_the_daughters_of_a_decayed_resonance and test_a_redrawn_mass_set_keeps_the_floor (both fail on the previous code), the floor decisions now through reshuffle_production.
…erformance-eval-122592 # Conflicts: # .github/gpu_ci/README.md # .github/workflows/gpu_runner_ci.yml
…erformance-eval-122592 # Conflicts: # .github/gpu_ci/README.md # .github/workflows/gpu_runner_ci.yml
…ity integration
A pdgs array of the wrong length matched the uncrossed branches on its leading codes (LPDGS copies only MIN(NPDG, maxpart) of them) while the folded-crossing fallback refused it, so folded and expanded outputs disagreed on such a call. Every path now answers NaN with the ERROR line.
…acters without an apostrophe FortranWriter re-split the overlong lines and read the apostrophe in "MadSpin's" as an open string, printing a '//' fragment in every generated f2py_wrapper.f. Comment only; pre-existing on origin/main.
… row, closed over the frame-dependent legs The good-helicity scan kept a row only above limhel/ncomb of the helicity sum on its RAMBO points, all at sqrt(s) = max(1500, 2*mtotal). A row suppressed by powers of m/E looks negligible there: for heft g g > b b~ two b helicity-flip rows are 2e-11 of the sum at 1500 GeV but about 4% each near threshold, where the cross section lives. They were dropped and the mg7 cross section came out 11% low (1.856e8 against 2.090e8 pb). The relative test of #217 only met #261's fixed-energy scan in the merge. With a massive external leg (ProcessTables::hel_massive) only an exact 0 is now dropped, and in every backend the list is closed over the legs whose helicity states mix under a boost (hel_frame_leg, from frame_dependent_helicity_legs, the C++ twin of the fortran warm-up's HMV), so that one list holds in every frame and for every crossing. Massless processes keep the relative test and its gluonic dedup. heft g g > b b~: 16 good rows instead of 14, 2.0905e8 +- 1.7e5 pb. test_generation_heft_mg7 now asserts 3% instead of 10%, and test_massive_leg_keeps_mass_suppressed_rows checks the 16 rows.
…rs.cc compile Under the complex-mass scheme Parameters.cc takes the sqrt of complex parameters (CMASS_mdl_MT = sqrt(MT^2 - i*MT*WT), mdl_cw, ...), and no sqrt accepted cxsmpl, so every mg7/madmatrix build failed (8 errors); MadSpin's mg7 decay generator broke test_complex_mass_scheme wherever madspace is built. The principal root goes through std::complex, next to conj in all four mgOnGpuCxtypes.h copies. HRDCOD=1 (constexpr_sqrt on cxsmpl) stays unsupported. Regression test: test_complex_mass_scheme_mg7_compile (no madspace).
debugging = True made every TestCmdLoop setUp rmtree the shared <MG5DIR>/TEST_AMC, so two concurrent test_manager runs in one checkout deleted each other's process directories (missing PS.input, missing Source/make_opts). Follow TestCmdShell2/MECmdShell: read unittest.debug in setUp, use tempfile.mkdtemp otherwise, and remove it in tearDown.
The good helicities are chosen once, on scan points at sqrt(s) = max(1500, 2*mtotal). Before 611c68d the relative limhel cut dropped the b helicity-flip rows there (4e-15 of the sum), and sm g g > b b~ came out 4% low at sqrt(s) = 1.3*2*MB and 29% low at 1.02*2*MB. The new test compares the C++ |M|^2 with the Combridge formula at 8 points from 1.02*2*MB to 600*MB to 1e-8: it fails on 9007da0 and passes on HEAD.
…se space With --crossing_table=all, t > b e+ ve also scans its crossings at their own masses, e.g. b~ > t~ e+ ve with sqrt(s) = m_b. RAMBO then STOPs the whole probe. The STOP exits with status 0, so the exporter found no good helicity and silently kept the compute-all matrix.f: recycling was lost for the base rows too. A 1->N configuration lighter than its products has no physical point and needs no row, so hel_warmup now skips it, as the C++ scan does. An empty warm-up measurement is now reported with a warning that shows the end of the probe output. Already present on 3a87b59 (PR #49 warm-up driver), not introduced by the integration.
…t, and the gpu init-scan note is its own paragraph
…hat the madmatrix and acceptance jobs share it
… only on C-parity-violating processes, the CPU C-parity dedup draws other rows on a C-symmetric one
…xtra crossings, the PDG-keyed smatrixhel/density keep answering the generated and folded subprocesses only
…ty of the beam, so a small parton x no longer loses digits to the beam mass get_betaz built beta from E*pz products of the massive beam (about 1e9 for x ~ 1e-4) whose difference is the m**2 term, losing up to 8 digits; with beta -> 1 the lab E and pz of an event then moved by up to 1e-5 (6e-3 at x = 1e-6) under rounding-level changes of the CM momenta. get_rapz takes log((E+pz)_in/(E+pz)_out) and zboost_with_rap scales E+pz and E-pz by exp(-+y), the same boost evaluated at machine precision. Folded and expanded crossing output of p p > w+ j now agree to 2e-11 in lab rapidity (was 2e-5).
…entity-only table, no router or cross-group dependent) is written on the plain path, without the CROSSUSE decode and per-crossing counters it would never use
…ged, so SMATRIXHEL answers 0 for it as the expanded output does, instead of wrapping it onto code-NCOMB
…leg-order repair instead of being frozen, so an order that needs a POS leg moved answers the right rho in the expanded and folded outputs; an order still unresolved falls through to the crossed fallback and then NaN, never a silent 0
…lone output zeroes its AMP again, so a diagram the flavor mask skips (u d > u d in q q > q q) reads 0 instead of leftover stack data in the density, GET_ALL_INTER and GET_VALUE
…nstead of 240 s, since the wall time follows the machine load (about 290 s at a load average of 100 for about 100 s of CPU), and checks that 'set nevents 500' was applied on the '500 written' count of the MadSpin summary line rather than on the timeout
… no longer raises AttributeError on every file it did not have to write plain (MadSpin run_bridge, test_hepmc_decay), and the plain-written .gz fallback can gzip its output
…ython gets its own name, so a second reweighting (another sample, or the next acceptance test) no longer binds to the first one's liball_rw_me_2me / liballrw_mevirt_2me
The dynamic loader resolves a library by its install name (SONAME on
Linux) to the one already loaded, whatever directory the new one comes
from, and nothing loaded can be unloaded. test_nlo_reweighting_new_parameter
run after test_nlo_reweighting therefore had its f2py module linked to the
ttbar rw_me library ("Symbol not found: _f77_m12_get_nhel_entry_, Expected
in <previous tmpdir>/rw_me/.../liball_rw_me_2me.dylib"); with matching
symbols it would silently have used the previous matrix elements.
- module-level loaded_f2py_names records each imported directory name;
get_f2py_suffix picks '_r<N>' for a name already loaded, and
compile_SubProcess_dir then passes PROCNAME (tree), ROOTNAME (loop) and
F2PYMOD with that suffix; the suffix map is saved in rwgt.pkl.
- makefile_sa_f2py and the MadLoop standalone makefile take the f2py
module name from F2PYMOD (default unchanged: all_matrix$(MENUM)py).
- load_module imports the suffixed name and forgets the stale parent
packages, which may point to another rwgt_dir.
Ported from the uncommitted follow-up on claude/beautiful-williamson-542354,
extended to rw_mevirt.
Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
… logical relpath would not resolve (macOS /var -> /private/var), so the MC@NLO Pythia8 run directory under the temp dir finds its xmldoc/tunes and test_gen_evt_onlygen writes events.hepmc.gz again
…tjet-config is not configured instead of a TypeError in Popen(None), and test_rivet_from_file / test_contur_from_file (and the mg7 twin) skip when pythia8, rivet, yoda, fastjet (and contur) are not installed, as the other mg7 post-processing tests do
…second copy of the param_card as tmp.dat in the MG5 root (or the check output path); compute_widths reads and updates the copy in its temporary directory, so the root file was never read and was left behind by every 'check cms' with recompute_width always/first_time
…e, so --hel_recycling no longer refuses the directory With --crossing_table=all, t > b e+ ve also records crossed decays of a light parent (b~ > t~ e+ ve, ...). scan_points found no RAMBO point at their masses and gave them the base points under the base id, so sigmaKin_getGoodHel scanned them off shell, the probe printed HEL/OFFSHELL and output.py wrote the directory without recycling, although the Fortran warm-up (b7f7b89) skips such a decay and recycles. A crossing lighter than its products (or a 2 > 1 lighter than its beams) is never evaluated on shell and needs no row: scan_points now tests rambo::has_phase_space and marks it with sigmaKin_setNoPhaseSpace, and getGoodHel (cpu, simd, gpu) leaves it out of the scan, in the probe and at run time alike, so the recycled table and the run-time list stay the same. The refusal now only covers a crossing RAMBO does not converge on. t > b e+ ve: 4 of 16 rows recycled in C++, as in Fortran, |M|^2 identical to the plain build.
…ds it now (scan_points in umami.cc and the --hel_recycling probe's off-shell list), no longer the limhel threshold, which stopped depending on it when the threshold became hel_massive ? 0 : limhel
… as the C++ good-helicity scan does since 611c68d Both --hel_recycling warm-ups sample fixed energies (C++ sqrt(s) = max(1500, 2*mtotal), Fortran 750-4500 GeV), where a contribution suppressed by powers of m/E looks negligible though it is not near threshold. 611c68d made the C++ row test exact for hel_massive, but two relative cuts were left: - the C++ probe (hr_probe_close, cpu and simd) dropped an amplitude below 1e-12 of the largest one of its row; with hel_massive it now drops only an exact 0, like MG_HR_EXACT_ZERO under --axial_gauge; - the Fortran standalone warm-up kept a row above LIMHEL/NCOMB of the sum; with a massive leg of the scanned crossing (RELHEL) it now keeps every non-zero row before the HMV closure. C-parity pairs are still taken only among rows above LIMHEL (GOODREL), so a row kept only as non-zero stays computed instead of being copied from a partner never compared. No wrong matrix element was found with the old cuts; this is hardening. Measured cost: C++ dropped pairs are unchanged (e+ e- > w+ w- z 135, u u~ > w+ w- g 54, others 0: all exact zeros). Fortran rows are unchanged for g g > t t~ g(g), u u~ > t t~ g, e+ e- > t t~ g, u u~ > w+ w- g, e+ e- > w+ w- z, u u~ > z g g, g g > b b~ b b~ and heft g g > h g, b b~; heft g g > h h keeps 4 rows instead of 2 (the Jz=2 rows are 1e-32 round-off at the boosted points), as the C++ scan already does.
…t parameters and propagator denominators) Under the complex-mass scheme the alphaS-dependent parameters are complex (mdl_sqrt__aS = sqrt(aS), G, mdl_G__exp__2). The writer of computeDependentCouplings_fromG only rewrote 'constexpr double' lines, so the constexpr cxsmpl<double> UFO assignments were kept verbatim: they read the undeclared aS, redefined G and were not constant expressions. They are now commented out (before G) or made per-event complex (const cxtype_sv, and cxtype in the SIMD-float loop), as the real ones are made fptype_sv. The CMS branch of the C++ ALOHA writer also cast the (complex) propagator mass to a real fptype_denom, which does not compile and would have dropped the -i*M*W of the complex mass; it is now cast to cxtype_denom_sv like the coupling. p p > t t~ j (the reported reproducer), g g > t t~ g, u u~ > t t~ g and g u > t t~ u now build (scalar/simd_128, FPTYPE d/f/m) and agree with the CMS Fortran standalone at the same point to 1e-8. HRDCOD=1 stays unsupported under CMS (constexpr_sqrt of cxsmpl). Regression test: test_complex_mass_scheme_qcd_mg7_compile.
… on crossed interference subprocesses (p p > j j and p p > t t~ j QED^2==2), as the helicity sum and with signed helicity draws (--rnd-hel) No GPU CI check reached the signed (nampso > 1) helicity selection on crossed extended ids: interference_checks.py has only single-subprocess split orders, and the crossed checks of madmatrix_checks.py have no squared-order constraint. crossed_gpu_vs_cpu gains signed=True: |M|^2, helicity and colour compared only above 1e-10 of the largest |M|^2 (vanishing interferences leave rounding residues, which must stay residues on the GPU), every crossing row with |M|^2 != 0, negative draws, and some draw that differs from the helicity sum (helicities of both signs: P1_gQ_ttxQ; P1_QQ_QQ brings the 11 rows). Checked with the host HIP emulation; reporting the base code in the signed GPU branch fails it (313 helicities differ).
…py runs (crossed MSSM, massive, axial-gauge and signed interference crossings, |M|^2 = 0 with good helicities), no longer only those of #255
…helicity-integration # Conflicts: # MadSpin/decay.py # madmatrix/output.py
This branch has not been deployed
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment
1 participant
Add this suggestion to a batch that can be applied as a single commit.This suggestion is invalid because no changes were made to the code.Suggestions cannot be applied while the pull request is closed.Suggestions cannot be applied while viewing a subset of changes.Only one suggestion per line can be applied in a batch.Add this suggestion to a batch that can be applied as a single commit.Applying suggestions on deleted lines is not supported.You must change the existing code in this line in order to create a valid suggestion.Outdated suggestions cannot be applied.This suggestion has been applied or marked resolved.Suggestions cannot be applied from pending reviews.Suggestions cannot be applied on multi-line comments.Suggestions cannot be applied while the pull request is queued to merge.Suggestion cannot be applied right now. Please check back later.
Follow-up to #49. Merges the open PRs that touch the helicity code #49 changes, and reconciles them with #49's design. Conflicts were settled in this order: correct physics, then one good-helicity list for all crossings (so helicity recycling works with crossing), then polarised computations, then correct helicity/colour selection.
Merged
pdgsarray.py_get_all_intercallable. Follow-ups: split orders summed (noINTERoverflow); recycledSMATRIXHEL(HEL=-1)now sums the helicities (it returned 0).check_sa --flavorstakes crossed ids).--hel_recycling, ported onto Crossing symmetry #49's crossing: one recycled table for every crossing; a zero amplitude is dropped only where it vanishes for every crossing; a row missing from the table is computed, not counted as zero; C++ axial gauge works on crossed legs.Decisions taken with the user
GET_ALL_INTER_IDXtakesNHELin the caller's leg order, likePandPOS.py_get_density/py_get_all_interfall back to the crossing, assmatrixheldoes. A pdg list no process serves now gives NaN and an error line.*legs keep crossing: the crossed external call uses the off-shell mass of the momentum in its slot.onshell_v1,onshell,madspin/fullandPAmodes. It also keeps crossing with polarised beams when the density backend is madmatrix.madspin_v1(msP/msF) stays without crossing.About 75 further fixes came out of the verification, in helicity selection, the good-helicity scans, polarised helicity codes, and the GPU/CPU behaviour for a zero |M|².
Tests
p p > t* t~ j, t~ > b~ w-: within 6e-16 for C++, recycled C++ and Fortran.p p > w+ j: 1 directory with crossing, 6 without. The decayed events are byte-identical with crossing on and off, for both density backends.Not verified
mainsince madspace: always draw the beam orientation before the mapping #245 (missing PDF set;run_card.tomlkey; DLPack read-only array), independently of this PR.beam_rapidity/mirror_beams.