Merge mg5amcnlo 3.7.3 and madspin_density into main - #82
Merged
Conversation
The scan's cost is dominated by reading decay events and evaluating their matrix elements -- the same work the joint scan does, and irreducible: to bound each slot's weight you need every slot's decay and density. Profiling a 4-top case put ~62% there, so no contraction trick brings the scan near the joint one. But the probe events are independent, exactly like the unweighting, so the scan parallelises the same way. get_sequential_maxwgt now reads the probe events into memory and, for nb_core > 1, forks one worker per contiguous slice (reusing the unweighting machinery: fork so self is inherited rather than pickled, per-worker RNG streams, _reopen_decay_pool for an independent pool view, results marshalled as JSON). Each worker returns its per-event maximum-weight vectors; the parent concatenates them -- order independent, since _combine_maxwgt takes the max and spread over all events -- and combines per slot as before. Probe weights are cast to python float so they survive the JSON round-trip. p p > t t~ [PA], nb_core=4: scan wall-clock 3.65s -> 1.86s, same maximum weights (1.042, 2.69) and efficiency as serial, cross section 483.81 vs joint 483.78, dilepton Delta-phi within 0.01 sigma. A unit test checks the event-range split concatenates to the whole scan. Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com>
The progress line printed "Efficiency: <1/self.efficiency>", which in
sequential mode is the decay events drawn per accepted event summed over
every decaying particle -- one number whose meaning was not obvious.
In sequential mode print it per particle instead, e.g.
decaying event number 5000. Decay events per accepted event,
per particle: p0=1.04 p1=3.20 [17.1 s]
so each particle's cost is visible (their sum is the old number). It is
emitted once per event -- gated on the slot-0 pass, since every accepted
event goes through slot 0 exactly once -- not once per slot. The joint path
keeps a single number, now labelled "Trials per event".
Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com>
Under the parallel unweighting (and now the parallel max-weight scan) every forked worker printed the "decaying event number ..." and "Event x/y" progress lines, so an 18-core run interleaved 18 copies. Gate both on the shard tag: print only from the worker whose _shard_tag is 0, and from the serial path (where it is unset). The other workers stay silent. Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com>
onshell and PA are implemented and validated; madspin still falls back to the joint accept/reject. Capture how to lift it. The obstacle: madspin evaluates the production density at reshuffled (offshell) momenta that couple all decay masses, so rho is not fixed while the chain is built. Olivier's fix: draw every decaying particle's invariant mass before the loop, reshuffle the production once up front, and reuse the fixed offshell rho -- which also moves the production-validity check earlier. Because the masses are then fixed, a decay that cannot be reshuffled to its mass forces a restart of the whole set rather than a local redraw. Also record the subtlety the fixed rho does not remove: joint madspin uses onshell denominators with an offshell numerator, whereas the per-particle method normalises by the offshell traces. To match, each slot must be normalised by the onshell decay ME (one extra evaluation per decay) and the telescoping re-derived; do not enable madspin sequential before that is done and A/B'd against joint madspin. Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com>
…ping
The onshell-denominator / offshell-numerator mismatch is not an open
problem: the per-particle weight
W_k = <rho_off, D1_off (x)...(x) Dk_off (x) I...> *
(jac_prod jac_1...jac_k) / (|M_prod|^2_on |M_1,dec|^2_on...|M_k,dec|^2_on)
accepted as W_k/W_{k-1}, telescopes to the joint madspin weight. In the
ratio the production ME, the production jacobian and the slots below k
cancel, leaving [P_k/P_{k-1}] * jac_dec_k / |M_k,dec|^2_on -- the existing
contraction with offshell rho and Dk, normalised by the ONSHELL decay ME.
This unifies with the built code: Dhat_k = D_k^off / Tr(D_k^on) for madspin
vs D_k^on / Tr(D_k^on) for onshell -- same denominator, offshell numerator.
Cost: one extra ME evaluation per decay (density before and after the
reshuffle).
Also record why the up-front mass draw is exact: the BW prior cancels the
BW jacobian carried in the weight, so the accepted mass distribution is the
physical marginal even though the mass is fixed per chain.
Still gate behind an A/B against joint madspin on a genuinely offshell
(large-width) process, since ttbar's near-onshell tops could hide a slip.
Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com>
Implements the full-offshell per-particle accept/reject per the design in
MADSPIN_SEQUENTIAL_PLAN.md section 10: _offshell_production draws every
virtuality up front, reshuffles a copy of the production to that mass set,
and evaluates the fixed offshell rho there; the offshell branch of
sequential_accept_reject weights each slot by
(N_k/N_{k-1}) * jac_bw_k * Tr(D_k^off)/|M_k|^2_on, with the per-chain
production reshuffling jacobian on slot 0. The offshell decay density is
taken on a copy so the drawn decay stays onshell and the final
add_decays + single reshuffle_production rebuild consistent kinematics.
PA and onshell are unchanged (71 unit tests pass, PA/onshell A/B still
0.008% on cross section).
Left gated off in _sequential_active: on ttbar it needs ~340 decay-ME
evaluations per event (slot 0 ~313, slot 1 ~27) versus joint madspin's
~122 (madspin is inherently peaked -- joint itself is 61 trials/event).
The per-mass-set reshuffling jacobian and the offshell tail land in slot
0's per-angle accept/reject, and with the mass fixed per chain an unlucky
draw cannot be escaped by redrawing angles. Making it worthwhile needs a
mass-set-level accept/reject before the per-angle loop; physics
correctness is also still unverified (the run is too slow to finish an
A/B). Recorded in the plan.
Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com>
…gated off)
Restructures the offshell (madspin/full) path so the mass set is
accept/rejected before the per-angle loop: w_mass = Tr(rho_off) *
jac_reshuffle * prod jac_bw_k carries every factor that depends on the
masses but not the decay angles, and the per-angle weights reduce to
(N_k/N_{k-1}) * Tr(D_k^off)/|M_k|^2_on. This correctly isolates the
production reshuffling jacobian -- its bound is modest (C_mass ~ 14 on
ttbar), where before it inflated slot 0 to ~322.
But it does not make sequential madspin competitive, for a structural
reason: the intrinsic tail is the per-angle offshell decay reweighting
Tr(D^off)/|M|^2_on (reweighting a pool decay to the offshell mass), whose
per-slot bound on ttbar is ~124 and ~161 -- each *worse* than joint
madspin's full-weight bound of ~61. The offshell tails are anti-correlated
across decays, so the joint test captures a cancellation the per-particle
factorisation loses. Sequential madspin needs ~280 decay-ME evaluations
per event vs the joint ~122; slower, structurally.
So madspin/full stay gated off in _sequential_active; PA and onshell are
unchanged and validated (71 unit tests pass). The implementation is kept as
the correct foundation and the finding is recorded in the plan.
Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com>
…ssue, not the factorisation madspin_v1 (weighted decays) does the same physics in ~7.8 trials/event, while joint density madspin is ~61. So density madspin (joint and sequential) is ~8x less efficient than it should be: a real inefficiency in the density-madspin weight, not a property of the per-particle decomposition (my earlier "structural/anti-correlation" conclusion had a wrong baseline). The sequential per-slot scan localises it: C_mass ~ 14 (production fine) but C_0 ~ 124, C_1 ~ 161 -- the tail is the per-decay offshell reweighting Tr(D^off)/|M|^2_on. The decay pool is generated in the onshell frame and the density mode reshuffles each decay to an offshell mass over BW_cut=15 widths (the code warns >25 breaks NWA) and unweights that; madspin_v1 keeps the decay weighted and avoids the cost. Open item is therefore the density-madspin offshell decay reweighting (shared joint/sequential), not the sequential factorisation. Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com>
The earlier conclusion that sequential madspin is slower than the joint test was an artifact of the inclusive `w+ > all all` decay, whose top non-resonant contributions blow up density madspin as a whole (61 trials/event joint) via the offshell reweighting tail -- independent of which accept/reject is used. For a physical resonant decay (t > w+ b, w+ > l+ vl) the mass-set-level restructure makes sequential madspin both correct and faster than the joint test, validated end to end on p p > t t~, same production events: efficiency 5.6 decay-ME evaluations/event (slot 0 = 2.1, slot 1 = 3.5) vs joint 8.9 (4.46 trials x 2 decays); cross section 23.7742 vs joint 23.7750 (0.003%); dilepton Delta-phi within 0.32 sigma on the full 10000-event sample. Re-enable madspin/full in _sequential_active. PA/onshell unchanged, 71 unit tests pass. The `w+ > all all` non-resonant blow-up is a separate density-madspin weight issue (shared with the joint test) recorded in the plan. Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com>
…tion Two follow-ups to the sequential accept/reject. 1. Max-weight bound margin 1.05 -> 1.10 in _combine_maxwgt. The sequential accept/reject cannot carry a per-slot overweight forward (it redraws a slot until it accepts, rather than doing staged unweighting), so a weight above the bound biases the sample directly. A slightly wider margin lowers that overflow probability. Shared with the joint bound, which only becomes a touch more conservative. 2. Reserve 50% more decays for the maximum-weight scan. The scan draws nevents_for_max (Nevents_for_max_weight * max_weight_ps_point) decays *per slot*, so a pdg with several slots (e.g. two tops) consumes several times the bare reservation, and the sequential offshell scan draws a few more on every restart. Without margin the scan drains the pool and forces a slow mid-scan refill (observed: >2*75*400 top decays). A 1.5x reservation keeps the scan inside its pool. Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com>
…-weight scan A forked max-weight-scan worker that ran its decay pool dry refilled a huge number of events (900k+ in a 4-top run, ~130s), because _scan_maxwgt_range passed the *global* remaining scan count (nevents - i) as nb_remain. That number sizes the decay-pool refill, and the refill already multiplies by nb_core to build one shared pool for every worker -- so the per-worker refill came out nb_core times too large. Pass the worker's own remaining range (stop - i) instead. Measured on p p > t t~ t t~: the mid-scan refill drops from ~927k to ~67k events (from ~130s to a few seconds), same maximum weights. Also correct the reservation comment from the previous commit: the measurement showed the scan does *no* restarts for tops (it draws exactly nevents_for_max per slot). The 50% reservation margin is still useful, but for a different reason -- the parallel scan stripes the pool across workers and an uneven split otherwise exhausts one worker's slice. Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com>
… file The parallel max-weight scan exhausted its decay pool and refilled even when the pool was several times larger than the scan needs. Cause: the decay pool is split into nb_core files (nb_unweight_output = the madspin nb_core, e.g. 18), but _scan_maxwgt_parallel reduced nb_core to the number of non-empty event ranges (e.g. 15 for 75 events over 18 cores). That made len(paths) != nb_core in _reopen_decay_pool, so every worker fell onto the striding branch -- and that branch opened evtfile.name, which for a split pool is only its *first* file. Each worker then strided 1/nb_core of a single file (~1/270 of the pool) and exhausted almost immediately. Two fixes: - _scan_maxwgt_parallel keeps the original nb_core as the pool-addressing count, so each worker opens its own file (paths[shard_id]); trailing empty shards are just not launched. - _reopen_decay_pool's fallback, for a split pool whose file count does not match nb_core, strides the whole *chained* pool instead of only its first file (correctness for any residual mismatch, e.g. nb_core > nevents). Measured on p p > t t~ and p p > t t~ t t~, 18 cores: the mid-scan pool refills go from present (and huge) to zero, same maximum weights. Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com>
…iple of nb_core) Round the number of max-weight probe events up to a multiple of nb_core, so the parallel scan splits evenly -- every worker gets the same number of events instead of one worker carrying an extra one whose pool slice can run short. Reduce max_weight_ps_point in proportion so the sampling budget (nevents * nb_ps_point) stays roughly the same. Default 75 events x 400 ps on 18 cores becomes 90 x 333 (90 = 5*18, total 29970 ~= 30000). Serial (nb_core=1) and already-even cases are untouched. Measured on p p > t t~, 18 cores: probe splits 5 events per worker, no mid-scan pool refill, same maximum weights. Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com>
… too The joint accept/reject's max-weight estimate (get_maxwgt_for_onshell) was serial while the sequential one already forked across cores. Give it the same treatment. - Generalise _scan_maxwgt_parallel to take the worker entry and its extra args, so it drives both scans (it forks one worker per contiguous slice of the probe events, each with its own RNG and reopened decay pool, and concatenates the per-event results). - Factor the joint per-event loop into _joint_maxwgt_range and add _joint_maxwgt_shard_entry (mirroring the sequential range/shard-entry). - get_maxwgt_for_onshell now reads the probe events into memory, applies the same even-split rounding (nevents up to a multiple of nb_core, nb_ps_point down to keep the budget), and dispatches serial vs parallel. Refill sizing uses the worker's own remaining range (stop - i). Validated on p p > t t~ [PA, joint], leptonic, 2000 events: serial 75x500 vs parallel 90x417 on 18 cores give the same cross section (23.7142) and the same unweighting efficiency (2.81 vs 2.77 trials/event) -- i.e. the same max-weight bound -- with no mid-scan pool refill. 71 unit tests pass. Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com>
"thread done. [X s]" was indistinguishable across the parallel unweighting workers. Print the worker's shard number and the worker count instead -- "worker 1 of 4 done. [1.5 s]" -- so the interleaved lines can be told apart; the serial path (no shard tag) prints "decay unweighting done". Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com>
MadSpin: bound the mass stage per production event, for any 2 -> N
…cation `sigma` came out identical for `BW_cut = 15`, `BW_cut = 1` or any other value. MadSpin samples each resonance's virtuality only inside +- BW_cut widths of the pole but normalises with the *full* width -- sigma_prod * BR on the density side, the param-card BR of the chain on the v1 side -- so it reported the whole rate while producing only the part of the Breit-Wigner inside the window. Measured against a truth sample (`p p > t t~ j, t > w+ b, t~ > w- b~`, where MG5's `gForceBW = 1` branch sets `cut_bw = .true.` and the truncation is in the integrated cross section), truth/MadSpin moves BW_cut = 15: 0.9628 -> 1.0052 (649.35 +- 2.19 pb against 646.02) BW_cut = 1: 0.4956 -> 0.9977 (334.28 +- 1.13 pb against 335.05) The second row is what settles it: the old code reported the same 674.44 pb for a sample it had cut in half. The factor is the sampler's own normalisation, not an approximation of it. Both generators draw m^2 flat in R = atan((m^2-M^2)/(M.Gamma)), whose full range is pi: `_mass_window` returns this as its gap/pi jacobian and `generate_inv_mass_sch` computes it as `bwdelf`. Integrating the density the code samples from over the window it samples in is closed-form, so `bw_retained_fraction` evaluates it exactly rather than falling back on the linearised 2/pi * atan(2N) (0.97869 against 0.97879 for a top at N = 15). What no self-consistent calculation can supply is the numerator -- it would need the integral over the part of the Breit-Wigner that was never sampled -- which leaves a residual of a few tenths of a percent, +0.4 % to +1.0 % for a t t~ pair at BW_cut = 15. Documented at the function and in doc/madspin_sequential_plan.md section 16. Which resonances differs by path, and the reason is the normalisation: * density (`madspin`, `full`, `PA`): the top-level virtualities only. The nested W of `t > w+ b, w+ > l+ vl` is not redrawn by MadSpin -- it comes from the MG5-generated decay events and is only boosted and rotated -- so its truncation is already inside the measured partial width the BR is built from. * v1 (`madspin_v1`): every resonance of the chain, nested included. `merge_itree` marks every decay-side s-channel invariant free and the driver BW-samples each; the v1 BR is the param card's, which carries no truncation. Measured on `t > w+ b, w+ > all all`: 0.95785 (density, top^2) against 0.91614 (v1, top^2 . W^2). * `onshell`, `onshell_v1`, `none`: no correction. They sample no virtuality, so inventing a loss for them would be the same error with the sign flipped. Verified: identical sigma at BW_cut = 15 and 1. * 2 -> 1 production: no correction. sqrt(shat) fixes the virtuality and nothing is drawn -- verified on `p p > w+`, which is why the acceptance-test cross sections are unchanged. The factor goes into `branching_ratio` before anything reads it, so it reaches both the <init> block and every event weight and sigma = mean(w) stays true under IDWTUP = -4. Composes with the overweight carry (which still shows up as mean(w) above XSECUP by the amount it reports) and with decay_output = weighted (whose mean(w) self-check now compares against the truncated reference: 1.48 sigma on a 100-event run). The parallel comparator learned to read the factor from the log and compares sigma_prod * BR with it divided back out, so the five-mode invariant survives; `assert_bw_truncation_matches_spinmode` asserts the per-mode factors themselves. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
…st it moves Finishes f9ee215 by measuring what it claims, on fresh samples. **The headline.** `p p > t t~ j`, 30 000 production events, `spinmode PA`: `sigma(BW_cut = 15) = 645.741 pb` against `sigma(BW_cut = 1) = 334.910 pb`, ratio 1.928103 against the 1.928102 the closed form predicts. Seven digits, because the factor is the sampler's own normalisation and not a fit to it. The old code reported the same number for both. `onshell` reports 674.1565 pb at both cuts -- unchanged, and equal to the production sample's own. **The physics.** Against truth samples `p p > t t~ j, t > w+ b, t~ > w- b~` at matched `bwcutoff`, truth/MadSpin moves BW_cut = 15: 0.9688 -> 1.0114 +- 0.0024 (truth 653.12 +- 1.14) BW_cut = 1: 0.4971 -> 1.0006 +- 0.0027 (truth 335.11 +- 0.68) The second row is new information the earlier study said it needed and could not have: it ran truth at `bwcutoff = 15` only, and bounded the residual to +0.4 %..+1.0 % by attributing it to the decay numerator `m.Gamma(m)/(m_t.Gamma_t)` running 0.52 to 1.71 across the window while the factor holds it flat. That explanation predicts the residual collapses when the window narrows, and it does: +1.1 % at 15, +0.06 % at 1. It is a window-width effect, not an offset, and it vanishes exactly where the correction is largest. Folded against the 5M-event truth of RESULTS.md (651.8 +- 0.22) rather than this 30k one, the `BW_cut = 15` residual is +0.90 % +- 0.05 %, inside the documented band; the +1.1 % above and the +0.5 % of the commit message are two 30k measurements 1.5 sigma apart about it. **The per-path asymmetry, measured rather than argued.** Same sample, same cut, `t > w+ b, w+ > e+ ve`: density (`PA`) reports 0.95785 = top^2, `madspin_v1` reports 0.91614 = top^2 . W^2, and `madspin_v1` on the un-nested `t > w+ b` reports 0.95785 -- agreeing with the density path when there is nothing nested to disagree about. **Blast radius: one number.** A sweep of every cross-section, BR and event-weight assertion in `tests/` finds exactly one that moves: `test_complex_mass_scheme`'s post-`decay_events` target, `440.779 -> 431.39` (-2.13 %, one top at 0.9786983 -- `t~` is not in the card and MadSpin does not auto-conjugate). Measured: production 442.887 +- 4.815, decayed 433.4528, ratio 0.9786983 to seven digits. It did **not** fail. `4*err1` on a 100-event run is +-4.3 % and swallowed a 2.13 % shift, so the count of *broken* tests is zero and the count of tests carrying a now-wrong expectation is one. Updated for the second reason: a tolerance wide enough to hide a systematic is not a check on it. The two targets were only ever equal because `BR(t -> w+ b)` is 1 to seven digits. The six `p p > w+`/`p p > w-` acceptance cross-sections are untouched, as the 2 -> 1 guard intends -- and that guard is load-bearing, not decoration: at 2 -> 2 the W factor would have moved `100521.5` by 2212 pb against an `error` of 800. `test_wj_production_with_ms_decay` is the one other affected path and already omits its cross references. Suite: 444 unit tests OK; `test_complex_mass_scheme` OK. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
MadSpin: the reported cross-section now carries the Breit-Wigner truncation
collect_events cached one open file + one mmap per input file, in a dict local to each copy worker. Since the event refs are shuffled before being split into chunks, every worker touches every input file and therefore opens all of them. mmap.mmap() dup()s the descriptor, so the cost was 2 * workers * nb_input_files descriptors. This was harmless while the workers were multiprocessing.Process, as each child had its own RLIMIT_NOFILE budget. Once they became threads of a single process the cost became multiplicative, and `p p > w+ j [QCD]` (12 event files) on an 18-core Mac needs 432 descriptors against the macOS soft limit of 256: OSError: [Errno 24] Too many open files Replace the per-worker cache with a single FDPool shared by all workers. Reads go through os.pread, which ignores the file offset and so lets every thread share one descriptor per file without locking around the read; entries are pinned while in use and evicted LRU once the pool reaches its cap. The cap is derived from RLIMIT_NOFILE, which is raised towards the (much larger) hard limit first. The disk-backed path had the same unbounded cache, single-writer but over an unbounded number of files, and now shares the pool too. Peak descriptor use over 40 files with 18 workers drops from 666 to 62, and the output stays byte-identical for a given seed. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
Review follow-up. FDPool treated max_open as advisory: when every pooled entry was pinned by a concurrent read, _evict_locked gave up and _acquire opened anyway, so the pool grew to max_open + workers rather than staying bounded. Wait on a condition variable instead. This cannot deadlock, since a reader holds at most one pin and always releases it, so an entry always becomes evictable. Measured with 32 concurrent readers against a cap of 4: 32 open before, 4 after. The k-way merge of the disk-backed path opens every shuffle run at once, and the run count is unbounded (one per EXTERNAL_RUN_RECORD_CAPACITY events). That is a pre-existing overrun in the path chosen for the largest jobs, and the previous commit's budget ignored those descriptors entirely: 139 runs against a limit of 128 fails on 3.7.3 and still failed after the first fix. Merge in passes of at most a bounded fan-in first, and split one budget between the input pool and that fan-in. A pass is only reached above fan_in * EXTERNAL_RUN_RECORD_CAPACITY events, so real jobs never take one; the output is unchanged either way, as merging sorted runs in passes gives the same key order as one big merge. _pread_exact returned a short read silently, which would write a truncated event; raise instead. The regression test now skips when the inherited hard limit is already below the limit it wants to impose, rather than failing on setrlimit. os.pread and resource are POSIX-only, which is not a new constraint here: MG5_aMC already requires make, gfortran and a POSIX shell. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
An overweight has been CARRIED rather than clipped since section 14, so it is
no longer a silent bias anywhere. What is left is how loudly to say it -- and
one region deserves less noise than the rest, because there is nothing to fix
in it.
MadSpin evaluates the production side with every resonance ON its pole. An
event that has not got the invariant mass to put them all there is asking the
factorisation for something it does not have: its Breit-Wigner windows stop
being set by BW_cut and start being cut off by the energy budget, and the
jacobian of the reshuffle onto a drawn mass set diverges because there is no
recoil left to absorb it. _near_nwa_threshold names that region,
sqrt(shat) < sum_r pole_r + _NWA_THRESHOLD_WIDTHS * sum_r Gamma_r
over the particles the event actually decays, counted with multiplicity.
sqrt(shat) because that is the quantity MadSpin itself spends as the mass-draw
budget, in _upfront_production and on the joint path alike; the margin in
summed widths because the width is the only scale that says how far off its
pole a resonance may go. _NWA_THRESHOLD_WIDTHS = 1.0.
Measured, p p > t t~ at 6.5+6.5 TeV, spinmode madspin, BW_cut 15:
50 000 production events x 400 free mass sets each (2.0e7 draws) put every one
of the 239 over-bound draws, on every one of the 14 events producing one,
inside 0.24 summed widths of 2 m_t -- a factor four inside this margin, in a
region holding 0.31 % of the sample. An independent 500 000-event joint run
says the same of all 17 of its overflows. p p > t t~ j does NOT: its 265
overflows in 300 000 events are nowhere near threshold in sqrt(shat) or in
m(t t~), and they keep the warning, which is the intended behaviour.
The end-of-run lines drop from warning to info only when EVERY carried
overweight is in the region, and quote both halves either way. The total stays
the first number on the line and no count, factor or cross-section shift moves:
this is a report, not an accounting change.
Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
Section 15 measured where the offshell overweights are, got `p p > t t~` right and `p p > t t~ j` wrong: it said "same divergence of J, reached by the mass draw rather than by the energy budget". Measured on the offending trials, J is 1.07 on them (0.98 to 2.47) against a sample-wide J_corner median of 1.22 -- the reshuffle does nothing. The whole factor is the matrix element. Once the production process has a final state besides the resonances, its matrix element carries a propagator of the resonance itself on the line the jet is radiated from, at (p_r + p_j)^2 = m_r^2 + 2 p_r.p_j. On the pole that is >= M_r^2 for any real jet, so it is unreachable; drawing m_r below the pole opens it, and there the production matrix element is a Breit-Wigner peak regulated only by M_r Gamma_r. In 2 -> 2 the only internal resonance line is t-channel and spacelike, so the region does not exist -- which is the whole of the 26x. All 265 carried overweights of the 300 000-event run sit within 5.2 M_t Gamma_t of that pole (51 % within one, the largest -- factor 48.9 -- at 0.12) against 3.30 for the trials that merely came close to the bound. Reachable inside BW_cut = 15 for 3.5 % of the production events, 0.9 % at 10, 0.03 % at 5; rerun there the count falls 265 -> 97 -> 13 and the tail collapses 48.9 -> 46.6 -> 6.9. No bound is changed. Raising it needs 49x and costs 49x. A per-event bound J_corner(e) x max(w/J^P) -- section 15's mass-stage construction -- is free in 2 -> 2 (3.43 trials/event against 3.46, all 17 overflows gone, measured over all 1.7e6 trials) and a 62x slowdown in 2 -> 3, because there J^P is 1 and dividing by it shrinks nothing: it fixes the case that does not matter. That is written up as a recommendation, not taken. What is changed is the report. _near_production_resonance tags the second region the way the previous commit tagged threshold, and the overweight line splits three ways -- threshold, production resonance, neither -- dropping to info only when the third is empty. On p p > t t~ j it tags 265 of 265 and the log goes from one WARNING to none, with the 300 000 written events byte-identical to the base tip's and every number on the line unchanged; p p > t t~ still goes through the threshold branch. Section 17 has the measurements, including what the population costs: the carry restores +0.245 % of the cross-section and 2.4 % of the effective statistics, and before PR #375 clipping was leaving the low tail of the top lineshape 25 % low. tests/test_manager.py test_madspin -t0: 473 tests, green (462 + 11 new). Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
…-- refuted Four readings of the hypothesis, on the run's own instrumented sample, with the 431 near-bound trials (0.3 < w/C < 1) as the control that is also heavy and 40 000 ordinary mass draws as the baseline. AUC between over-bound and near-bound, where 0.500 means the variable says nothing: sqrt(shat) 0.523 (readings 1 and 3) sqrt(shat) - sum m_i' 0.525 (reading 2) sum m_i'/sqrt(shat) 0.474 |chi - 1| 0.490 (reading 4) d = |(p_r+p_j)^2 - M^2|/(M.Gamma) 0.133 (1)/(3): the over-bound trials ARE mildly shifted down in sqrt(shat) -- the median alone hid that, 634 GeV against the sample's 716 -- but their own minimum is 387.8 against the sample's 368.9, only 1.5 % are below 400 GeV, and the lowest decile holds 16.2 % of them. A 1.6x enrichment, which the mechanism predicts: a softer jet reaches the pole more often. (2): refuted backwards. The over-bound never come within 61 GeV of the reshuffle boundary and never fill more than 84 % of sqrt(shat); ordinary draws reach 21 GeV and 94 %. They are FURTHER from their own boundary than a random draw is. (4): |chi - 1| separates heavy from ordinary strongly (AUC 0.885) but over-bound from near-bound not at all (0.490). It is the marker of "a mass was drawn low", which is what makes J^P 1.07 instead of 1.00 -- and 1.07 is not 48.9. Nothing adds anything beyond d: Spearman r(ln w/C, d) = -0.688, and the partial correlations given d are +0.14, +0.15 and -0.08, from raw correlations of +0.04 to +0.05. What residual there is has the wrong sign for the hypothesis -- within d quartiles, HIGHER sqrt(shat) gives a higher weight. And "drawn low" is not the mechanism in disguise: Spearman(d, min m') = 0.007 over ordinary draws, and of the draws that put a resonance more than 8 widths below its pole only 0.87 % land inside d < 5. Necessary, and short of sufficient by a factor 115; the extra condition is on the jet. Doc only, in section 17. No code change; 473 tests still green. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
…-2 partner
Sourcery flagged that `_near_production_resonance` lets any status-2 particle
serve as the partner `k`, so two decayed resonances could be mistaken for a
resonance-plus-jet production pole, and suggested restricting `k` to status 1.
Checked numerically. The test fires only when (p_r + p_k)^2 lands in
M_r^2 +- 10 M_r Gamma_r, and (p_r + p_k)^2 >= (m_r + m_k)^2 for any two
physical momenta, so a status-2 partner needs
m_r + m_k <= sqrt(M_r^2 + 10 M_r Gamma_r)
while MadSpin never samples a mass more than BW_cut widths below its pole
(decay.py: m_min = max(mpole - BW_cut*w, 0.5)). Sweeping t/W/Z/h in pairs with
both masses at that floor: t+t~ needs BW_cut >= 55.3 (303 GeV against a window
stopping at 181.7), t+W >= 20.5, and the tightest pair Z+W misses by 1.6 GeV at
the default BW_cut = 15, opening only from 15.4 -- past the point MadSpin's own
check already calls too large for the narrow-width approximation.
The suggestion is not applied. `status` records whether MadSpin attached a
decay, not whether the particle was radiated off the r line, so the guard would
make the same momenta tag or not depending on the user's decay card; and where
the false positive does open, the same pairing is the genuine W* -> W Z
production resonance, so it would cost a true positive at exactly that setting.
No behaviour change: a comment at the pairing loop and three tests that pin the
boundary instead of assuming it -- the Z+W walk across BW_cut 5 -> 25, the top
pair at any usable BW_cut, and a light partner tagging identically whether or
not it was decayed.
Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
The early return added by 083fb7e divided by sum w unconditionally and sat above the z >= _OVERWEIGHT_MIN_Z branch that exists because sum w is not a usable denominator for a sample whose weights cancel, so an all-zero sample raised instead of reporting. The shift the line quotes is now computed once per branch and named, and the smallness test is made against that same number: sum w when it is a legitimate denominator, sum |w| when it is not. It therefore never divides by a cross-section consistent with zero, and a sample with no scale at all gives nan, which compares false and keeps the full note rather than being declared small on an undefined ratio. The test is on the magnitude: a sample with counter-events can lose several per-cent of its cross-section with every carried factor still above 1, and that is not "very small". The 0.5 is now _OVERWEIGHT_QUIET_PERCENT, next to the other named thresholds here, and the region split is computed after the early return so that it is structurally clear it plays no part in the quiet decision. Also drops the now-unused `joint` local and restores the space the shortened production-resonance sentence lost before "The other N are NOT ...". Tests: the report now has two regimes and the suite describes both. The fixtures of TestNwaThresholdReport and TestProductionResonanceReport are scaled to +3% so they still reach the split they are about, and a new TestOverweightQuietBelowThreshold pins the quiet one -- including the unexplained overweight that used to warn, the boundary, the negative shift, the cancelling sample on either side of the threshold, and that the head of the line is character for character the same across it. Three assertions that had been passing on text the verbosity change deleted ('invalid by construction', 'Lowering BW_cut', a bare '-3') are retargeted at text the messages still contain. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
MadSpin: the 2 -> 3 joint overweights are the production's own resonance, not the jacobian
MadSpin: tag overweights inside the NWA-invalid threshold region
…n it fits Three follow-ups to the pooled reader: - A lock-free fast path. Once the soft limit is raised, the whole input set almost always fits, and then there is nothing to evict and nothing to synchronise: DirectFDReader holds every file open and reads straight from the descriptor. This is now the normal case. - No sharing machinery in the single-writer external path. The disk-backed copy has exactly one writer, so when it does have to bound the descriptors it uses LRUFDReader -- the same LRU without the lock, the pin counting or the condition variable. The single-worker in-memory branch uses it too. FDPool is left for the only case that needs it: several copy workers sharing a reader too small for the input set. - Sizing from the actual counts instead of fixed 320/256/64 caps. fd_budget() raises the soft limit only as far as the files in hand require, so a small job no longer disturbs the limit at all, and plan_external_fds() splits the result between the merge fan-in and the reader rather than capping both. Measured on the workload this was reported against -- 5M events of p p > e+ ve [real=QCD] as 1000 files x 5000 events, 5.5 GB, under the macOS soft limit of 256. Before, the reader held 256 of the 1000 files and roughly three of every four events paid an open()+close(); now it holds all 1000: events.lhe.gz 815 s -> 379 s (2.2x, mean of 2 and 2 runs) events.lhe 716 s -> 263 s (2.7x, mean of 3 and 6 runs) nb_core does not enter into it: above 128 input files the disk-backed path is selected and it is single-writer by construction. Output is byte-identical across both versions and both core counts (one md5 over 5 runs, 5000000 events). Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
…imit Fix "Too many open files" in collect_events on many-core machines
# Conflicts: # .github/workflows/acceptancetest.yml # .github/workflows/warm_cache.yml # UpdateNotes.txt # VERSION # madgraph/interface/madgraph_interface.py # madgraph/various/misc.py # tests/acceptance_tests/test_cmd.py # tests/acceptance_tests/test_cmd_madevent.py # tests/unit_tests/various/test_banner.py
Conflict resolutions of note: * Standalone Fortran templates (matrix_standalone_v4 / _splitOrders_ / _msF_, check_sa.f): main's FLAV_IDX / FLAVOR threading and madspin_density's beam-polarisation reweighting now coexist. In _msF_ the reweighting is kept before the helicity sum, which is madspin_density's bug fix. * MadSpin/decay.py: main's per-flavor max_weight and madspin_density's overweight `carry` are combined on the same per-flavor bound. Every matrix-element call now goes through options.beampol_me(), including two sites outside the conflict that main still had on the old (0.5, 0.5) convention. * MadSpin/interface_madspin.py: madspin_density's sequential-maxweight and parallel-unweighting paths plus main's merged-particle PDG handling. The cached f2py lambda takes pdg_list per call *and* dispatches production / decay to mymod[0] / mymod[1]. * MadSpinInterface._revert_merged gains a class-level default: main set it as an instance attribute on one path only, while madspin_density's new tests drive the density path from a bare stub. * .github/workflows: MG5's gating and ubuntu-24.04-only runners are kept; madspin_density's new jobs are imported onto them. * tests/input_files/IOTestsComparison: resolved to the current reference. The goldens still need a regeneration pass -- the merged templates declare an extra JJ loop variable that the stored references do not have. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
In makefile_sa_f_sp the static-library rule for libmodel ran make in
Source/DHELAS, so libmodel.$(libext) was never produced by its own rule.
The two .$(dylibext) rules just below it are correct (DHELAS -> DHELAS,
MODEL -> MODEL), which is what makes this a copy-paste slip rather than
an intentional dependency.
It is normally masked because lib/libmodel.a already exists by the time
the standalone check binary is linked. From a clean tree it is not:
$ rm -f lib/libmodel.a lib/libdhelas.a Source/*/*.o
$ make -C SubProcesses/P1_epem_mupmum
ld: library 'model' not found
collect2: error: ld returned 1 exit status
With the rule pointing at Source/MODEL, the same clean build produces
lib/libmodel.a and links check.
Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
…ensity acceptancetest_85, _91 and _92 (test_gen_evt_onlygen and the herwig6 / pythia6 stdhep tests) were commented out on main and re-enabled by the madspin_density merge. They are legacy stdhep paths that are not expected to pass here, so remove the enabled copies again. The commented-out stubs main already carried are left untouched, so the acceptancetest.yml job set is now main's plus only the genuinely new jobs the merge brought in. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
calculate_matrix_element kept main's model_init block, which assumes
mymod is a single f2py extension. The density/onshell modes keep the
production and decay extensions in self.f2py_module[0]/[1], so mymod is
a list there and the block died with
AttributeError: 'list' object has no attribute 'initialise'
as soon as the max-weight scan evaluated a production matrix element.
Those modules are already initialised through initialise_f2py_module,
guarded by model_init_prod / model_init_decay, so restrict this block to
the single-module v1 path it was written for.
Merge fallout: main added the block on a path only the v1 modes reached,
and madspin_density added the list-based module storage.
Reproduced and fixed with:
MadSpin/madspin <card> # spinmode madspin, tree-level ttbar LHE
which now runs to completion and writes events_decayed.lhe.
Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
ProcessExporterFortranSA.make() fell back to 'make ../lib/libdhelas.a'
and '../lib/libmodel.a' when the plain 'make' failed. Those are only
real targets when the makefile was configured with the default static
libext; with 'dynamic' set the makefile knows only
'../lib/libdhelas.$(libext)' and make stops with
make: *** No rule to make target '../lib/libdhelas.a'. Stop.
That message was then re-raised as the user-visible failure, hiding why
the plain 'make' failed in the first place.
Build through the libext-agnostic 'libdhelas' / 'libmodel' phony targets
instead -- both Source/makefile templates provide them -- and re-raise
the original error when that does not help either, so the reported
message is the actual compilation failure.
Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
madspin_density moved the end-of-run accounting into _apply_accounting
and logged the headline line with logger.info. MadSpin's do_launch
carries @misc.mute_logger(), and both the madevent and the mg7
post-processing adapters run MadSpin with the decay loggers raised above
INFO, so the line was dropped from the run log entirely.
main relies on it being there: test_madspin_mixed_flavor_decay_log_summary
and its _mg7 variant both assert on it, and the former even documents the
expected form as
CRITICAL: MadSpin unweight efficiency: 0.3697
Restore the level main used. Both tests were failing on the merge with
"density-mode summary line not found in log".
Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
…base name
combine_name has two naming schemes. The same-base one (FFV1_2) already
picks the right slot for the FLV_Coupling flag: '%(propa)s' for a
wavefunction, '%(tags)s' for an amplitude. The cross-base fallback that
FFV2_FFS1 takes appended '%(propa)s' unconditionally.
For an amplitude (outgoing == 0) HelasAmplitude.get_helas_call_dict fills
'propa' with '' and puts the flag into 'tags', so the call site emitted
FFV2_FFS1_0 while ALOHA wrote the routine as FFV2_FFS1M_0:
check gauge e+ e- > ve ve~ w+ w-
NameError: name 'FFV2_FFS1_0' is not defined.
Did you mean: 'FFV2_FFS1M_0'?
Every backend shares this function, so the same mismatch shows up in
Fortran as an undefined _ffv2_ffs3_0_ at link time. Guard the fallback the
same way the first scheme is guarded; non-FLV models have tags == '' and
are unaffected.
Pre-existing: no CI job ran the check-gauge tests before
acceptancetest_check_gauge arrived with the madspin_density merge.
Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
These two tests came from upstream, where apply_flavor_grouping does not
exist and 'check gauge p p > w+ w-' evaluates four light-quark
subprocesses. In MG7 grouping is on by default, so the four are carried by
the single merged matrix element Q Qx > w+ w- and the gauge block reports
Q Qx > w+ w- ... Passed
Summary: 1/1 passed, 0/1 failed
which is the shipped behaviour, not a regression: the per-flavor coverage
lives in the flavor-grouping block, which compares the merged matrix
element against the unmerged one for every flavor and both orderings
(8/8). Assert that too, so the coverage the gauge block no longer provides
is still pinned by the test.
Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
unittest_10 runs 'tests.unit_tests.madspin.test_madspin' but only checks
out the repo -- no pip cache, so no numpy on the runner. That was fine
until the madspin_density merge, whose new test classes
(TestPureInterferenceRestriction, TestScanMaxwgtDecomposition,
TestSequentialAcceptReject, ...) work on arrays and import numpy:
Ran 502 tests in 58.167s
FAILED ( errors=67)
ModuleNotFoundError: No module named 'numpy'
Upstream covers those tests from unittest_madspin_sequential, which does
restore the pip cache, so neither branch saw this on its own.
Restore the pip cache in unittest_10 as well rather than dropping the
module from it: that keeps main's coverage exactly as it was. With numpy
available the same 502 tests pass.
Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
_build_gridpack wrote a bare 'output <dir>'. That means madevent upstream,
where this test comes from, but MG5's default output mode is 'mg7' in
MadGraph7, so it produced a madspace/cudacpp tree carrying a
Cards/run_card.toml and no Cards/run_card.dat:
FileNotFoundError: [Errno 2] No such file or directory:
'/tmp/ro_gridpack_xxxxxxxx/PROC/Cards/run_card.dat'
The two assertions before it (mg5_aMC returned 0, the process directory
exists) both passed, which is why the failure only surfaced on the
run_card read.
Everything the test does afterwards -- flipping run_card.dat to
gridpack=True, bin/generate_events, run_01_gridpack.tar.gz, the extracted
madevent/ tree and bin/internal/restore_data -- is madevent-only, so ask
for that exporter explicitly.
Runs green locally: 3 concurrent workers off the frozen gridpack, 184s.
Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
test_pythia8_delphes_parallel wrote a bare 'output %s -f'. That is the
madevent exporter upstream, where the test comes from, but MG5's default
output mode is 'mg7' in MadGraph7, which writes no HTML/ directory at all.
load_result reads HTML/results.pkl, so the run surfaced as
TypeError: 'NoneType' object is not subscriptable
out of check_parton_output, with nothing pointing at the output mode.
This was the only bare 'output' left in test_cmd_madevent.py; every other
madevent test in the file already spells it 'output madevent'.
Also check mg5_aMC's return code before reading the results: the run had
already failed by then, and without the check that failure could only ever
show up as the TypeError above.
Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
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
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.
Brings
mainup to date with two upstream mg5amcnlo branches, in this order:mg5amcnlo/3.7.3(tipfb8ece13a)mg5amcnlo/madspin_density(tip02ee49a59)madspin_densitydoes not contain3.7.3's tip, so the order matters and both merges are real (no fast-forward).Conflict resolutions worth reviewing
Standalone Fortran templates
matrix_standalone_v4.inc,matrix_standalone_splitOrders_v4.inc,matrix_standalone_msF_v4.inc,check_sa.f.Both branches edited the same helicity loops:
mainthreadsFLAV_IDX/FLAVORthroughMATRIXandGET_ALL_INTER,madspin_densityinserts the beam-polarisation reweighting. Both are kept. In_msF_the reweighting stays before the helicity sum — that ismadspin_density's bug fix, and dropping it would leave the v1 decay distribution unpolarised.MadSpin/decay.pymain's per-flavormax_weightandmadspin_density's overweightcarryare combined on the same per-flavor bound (decay_mw_for_event), so the safety net and the accept/reject threshold cannot drift apart.options.beampol_me(). Two sites outside the conflict regions still hadmain's old(0.5, 0.5)convention and were converted — leaving them would have applied a spurious polarisation factor to unpolarised runs.MadSpin/interface_madspin.pymadspin_density's sequential-maxweight and parallel-unweighting paths, plusmain's merged-particle PDG handling.pdg_listper call (main's merged-ID fix) and dispatches production/decay tomymod[0]/mymod[1](madspin_density's split ME dirs)._revert_mergedbecomes a class-level default.mainset it as an instance attribute on a single path;madspin_density's new tests drive the density path from a bare stub. Without the default, 3 tests error withAttributeError.CI workflows
MG7's gating and
ubuntu-24.04-only runners are kept throughout; upstream's new jobs are imported onto them (readonly_gridpack,acceptancetest_check_gauge,acceptancetest_PA_decay,unittest_madspin_sequential, …). The scheduledwarm_cacherun now only warms the heptools cache (keep-warm-heptools) instead of rebuilding the chain, and no longer deletes caches.acceptancetest_85/_91/_92(test_gen_evt_onlygen, hw6/py6 stdhep) were commented out onmainand re-enabled bymadspin_density; they are imported as enabled. Easy to re-comment if they are not expected to pass here.Test status
./tests/test_manager.py -t0— 1489 tests, 5 failures, 1 error.The 1 error is environmental:
test_DensityMatrixObservables22needsscipy, which is not installed on the machine that ran this.The 5 failures are all IOTest golden drift, and all the same single line:
The merged templates declare the
JJloop variable the beam-polarisation loop needs. Per instruction the conflicted goldens were resolved to the current reference and are awaiting a regeneration pass.Note that the goldens are currently in a mixed state: the 18 conflicted ones kept
main's content, while 39 others (improve_ps.f,loop_matrix.f) auto-merged and tookmadspin_density's. A regeneration pass overtests/input_files/IOTestsComparison/will settle both.All MadSpin and lhe_parser unit tests pass (531/531).