Muon collider MAIA detector postprocessing

[jmduarte]

• Update postprocessing to work with muon collider MAIA detector data
• One file for testing: https://uaf-3.t2.ucsd.edu/~atuna/muoncollider/data/mlpf/ttbar/v04/ttbar_reco_10000.slcio.edm4hep.root
• Check equivalent to script before rebase
• Investigate bias of target definition

@farakiko @alexandertuna

[jmduarte]

Bug caught by Claude

Root cause: the merge clobbers repeat hosts

In postprocessing.py:923, vec1 is rebuilt from the original gpdata.gen_features on every iteration, but the result is written back into the running `pt_arr/eta_arr/phi_arr/energy_arr:

vec1 = vector.obj(pt=gpdata.gen_features["pt"][idx_gp_bestcluster], ...)   # original, every time
vec = vec0 + vec1
energy_arr[idx_gp_bestcluster] = vec.energy                                 # accumulated target

So the second unmatched particle merging into a given host computes E_host_original + E_new and throws away everything earlier merges added. In your 5 events: 145 particles merged into 61 hosts, 14 hosts received more than one, and one host received 30 merges — 29 of them discarded. That's 1981 GeV, 5.2% of all visible gen energy, silently gone.

Changing vec1 to read from the running arrays recovers exactly the predicted 1981.4 GeV, and the ceiling goes back to 1.000

[jmduarte]

Following up on the "Investigate bias of target definition" item. Summary of what I found on the MAIA ttbar sample (10 events, 42 matched jets) — three separate issues, two of them pre-existing on main and affecting CLIC/CLD as well.

The starting symptom: target jet $p_T$ / truth jet $p_T$ had max 0.990 and not a single jet above 1.0. A hard ceiling rather than a smear means the target is a strict subset of truth by construction, so it was a question of finding where energy leaks rather than a resolution effect.

1. Merge overwrite (629d67f) — pre-existing, from #463

In assign_genparticles_to_obj_and_merge, unmatched genparticles are merged into the genparticle owning their best cluster. vec1 was rebuilt from gpdata.gen_features on every iteration while the result was written back to pt_arr/energy_arr, so a second merge into the same host computed E_host_original + E_new and discarded whatever earlier merges had accumulated.

On this sample: 145 particles merged into 61 hosts, 14 hosts received more than one, and one host received 30 merges — 29 discarded. That silently dropped 1981 GeV, 5.2% of the visible gen energy.

This is much worse at a 10 TeV muon collider than at CLIC/CLD simply because the jets are denser (60–120 constituents at 1–5 TeV), so clusters are shared far more and multi-merges are routine.

Note the assert below the merge loop was meant to catch exactly this, but it is one-sided:

assert (np.sum(gen_features_new["energy"]) - np.sum(gpdata.gen_features["energy"])) < 1e-2

The merge can only ever lose energy, so the expression is always negative and the assert can never fire (it evaluates to −444.7 here). I left it alone, because making it two-sided would immediately trip on the remaining continue path a few lines above, which silently drops genparticles that cannot be associated with a unique cluster (237 particles, 445 GeV, highest 10.95 GeV $p_T$). Worth fixing together.

2. gp_to_hit weights were not energies (65a9078) — pre-existing

The weights mixed three incompatible scales: tracker hits contributed a literal 1.0, while calorimeter hits contributed raw CalohitMCTruthLink weights, which are proportional to deposited energy but with a constant that differs per subdetector (it absorbs the sampling fraction). Summing the weights on a hit and dividing by the reconstructed hit energy gives ~151 in the ECAL and ~603 in the HCAL.

So anything reading these as an energy was wrong. mask_visible_hit claimed to be "at least 10% of their energy", but the quantity it cut on had a median of 20–230 — it was effectively "has at least one linked hit".

Each link now takes the genparticle's share of the total weight on its hit, scaled by that hit's calibrated energy. That is correct whatever convention the weights follow, so it should hold for CLIC/CLD too. E_deposited/E_gen now has median 0.66 with 99.6% of values in [0, 1.5], i.e. a real energy fraction. gp_in_calo becomes meaningful for the same reason.

3. The 10% cut alone misses MIPs (51b34ae)

Once the units are right, the fractional cut is a shower-containment criterion and drops 6 of 7 muons: a muon deposits a roughly constant few GeV whatever its momentum (3.4–23 GeV here, median 4.2), so its fraction falls with energy (0.006–0.04).

gp_in_tracker cannot serve as the escape hatch — none of these muons have a track linked above its 20% threshold, so OR-ing with it would have recovered none of them. Instead the mask is now OR-ed with an absolute deposit threshold, which recovers all 7. Both thresholds are module-level constants:

visible_energy_fraction = 0.10
visible_energy_deposit = 1.0  # GeV

The absolute one is a detector-dependent scale that has to sit below the MIP peak and above the noise. 1 GeV is in the middle of a broad plateau — anything from 0.1 to 2 GeV keeps 7/7 muons and 84–85% of the status-1 energy — but it will want revisiting for CLIC/CLD, which have thinner calorimeters and a lower MIP peak.

Net effect

	mean	median	max
original	0.848	0.894	0.990
+ merge fix	0.877	0.933	1.000
+ corrected weights, 10% only	0.841	0.896	1.000
+ OR MIP deposit > 1 GeV	0.876	0.929	1.000

The merge fix is strictly monotonic per jet — 21 improved, 21 unchanged, none degraded — and recovers 4904 GeV of target jet energy. The response ends up close to where the merge fix alone put it, but now via a mask that actually means what it says rather than one that was accidentally passing nearly everything.

Open questions

1. The absolute term admits leaked neutrals. The 20 particles it lets in deposit only ~2.6% of their generated energy. That is right for MIPs, whose momentum comes from the track, but 9 of the 20 are photons that leaked (e.g. a 237 GeV photon depositing 3.2 GeV) and have no track — so MLPF would be asked to predict energy nothing measured. Restricting the absolute term to charged particles excludes them and gives 0.834. I did not apply that; it seems like a call for the group.

2. There is a ~15% floor that no mask can reach. 186 particles carrying 13455 GeV deposit exactly zero, mostly photons down the beampipe. Since the truth jets are clustered from all status-1 particles including those, the target can never reach 1.0. Probably wants to be addressed on the truth side (defining the reference within acceptance) rather than in the mask.

3. CLIC/CLD validation. Item 2 above touches every Key4hep detector, not just MAIA. Someone should run a before/after on a CLD sample before this merges — I have only exercised the MAIA path.

Caveat throughout: 42 matched jets from 10 events is thin, so treat the individual numbers as indicative.

cc @jpata @farakiko @alexandertuna

[jmduarte]

CLD validation, as promised in the previous comment. Test file from scripts/fetch_test_data_cld.sh (reco_p8_ee_ttbar_ecm365_300000.root), 50 events, --detector cld. Runs clean, no assertion failures.

CLD does not have the pathology

Worth stating up front, because it changes how to read the MAIA numbers. On main, CLD target response is already 0.9923, median exactly 1.000, with 98 of 304 jets above 1.0. There is no hard ceiling.

So the bias I was chasing is specific to dense TeV-scale muon-collider jets, where clusters are shared heavily and the merge overwrite compounds. It is not a general defect in the target definition.

Effect of this PR on CLD, per commit

289 jets matched in all four variants:

variant	mean	median	max	>1.0	Δ
main	0.9923	1.0000	1.022	98
+ merge fix (629d67f)	0.9926	1.0000	1.037	99	+0.0003
+ corrected weights (65a9078)	0.9815	0.9987	1.047	78	−0.0112
+ MIP OR (51b34ae)	0.9879	0.9991	1.047	81	+0.0064

Net −0.0045 (−0.45%): 6 jets improved, 211 unchanged, 85 degraded.

The merge fix is a no-op on CLD (+0.0003), exactly as expected — sparse 365 GeV jets almost never produce a host that takes more than one merge, so the bug never fires. Compare MAIA, where one host absorbed 30 merges.

The regression is entirely the mask now actually cutting at 10%, about half of it recovered by the MIP OR. I want to be straight about what this means rather than dress it up: on CLD, main's broken mask was accidentally producing a better-looking response by passing essentially anything with a linked hit, and correcting the units costs 0.45%. Whether genparticles that deposit less than 10% of their energy belong in the target is a physics decision, not a bug — flagging it here rather than deciding it.

The MIP threshold is detector-dependent, and 1.0 GeV was too close to the edge

Median muon deposit is 4.22 GeV on MAIA but 1.98 GeV on CLD — the shallower calorimeter roughly halves it. Usable range for visible_energy_deposit:

• MAIA: 0.1 – 2.0 GeV
• CLD: 0.05 – 1.0 GeV (at 2.0 GeV it cuts into CLD's MIP peak, keeping 15/20 muons instead of 19/20)

1.0 GeV sat mid-plateau for MAIA but right at CLD's upper edge, so I lowered the default to 0.5 GeV (a9ff5e3), which is inside both with room either side. It costs nothing: CLD output is bit-identical (0 of 302 jets change) and MAIA moves 2 of 42 jets by +0.0002.

This parameter will want revisiting for any detector with a different calorimeter depth. It is a module-level constant now, so that is at least a one-line change.

One thing the jet numbers hide

On CLD, without the MIP OR only 10 of 20 muons survive the mask, versus 19 with it. That barely moves jet $p_T$, so it would not show up in any of the tables above — but it would show up in muon reconstruction. Same story as MAIA, where the fractional cut alone dropped 6 of 7.

Two notes on the fetch script

• It calls wget, which is not present on macOS by default; the script has set -e but the failure lands inside a pushd subshell, so it exits 0 having downloaded nothing and silently produces no parquet. curl -L with the same URLs works.
• The closing message points at tests/validate_parquet_gt.py, which does not exist — tests/ has validate_cms_gt.py and validate_inclusive_hits.py.

Still open

The two items from the previous comment stand, and CLD does not change either: whether the absolute term should be restricted to charged particles (it admits leaked neutrals whose energy nothing measured), and the ~15% of MAIA status-1 energy carried by particles that deposit exactly zero, which caps the achievable response there and is probably a truth-side acceptance question.

Numbers are from 50 CLD events / 10 MAIA events, so treat them as indicative.

[jmduarte]

Correction to something I claimed in both comments above, plus a question about the MAIA sample.

The "~15% floor" was the wrong number

I said the ~15% of status-1 energy deposited by nothing at all "sets a floor on the achievable response". That is wrong in the way that matters, because jets are clustered with min_pt=5 and so what counts is $p_T$, not energy:

zero-deposit status-1 particles	energy	scalar $p_T$
MAIA	13455 GeV (14.8%)	1814 GeV (3.0%)
CLD	64 GeV (1.0%)	35 GeV (0.7%)

The 14.8% is dominated by four photons at $\eta$ = 8.8–9.2 ($\theta \approx 0.01$–$0.03^\circ$) carrying 783 / 2514 / 2541 / 4924 GeV. Their $p_T$ is 0.39 / 0.78 / 0.49 / 1.26 GeV — all below the 5 GeV jet cut, so they never enter a gen jet and contribute nothing to the response plots. The honest figure for the jet-level floor is ~3%, not 15%.

It is also not a forward-acceptance effect

I had assumed the forward region was the culprit. It is not:

• 99.5% of MAIA's zero-deposit $p_T$ is at $|\eta| < 2.5$, i.e. in the barrel and endcap.
• The nozzle region proper ($2.44 < |\eta| < 4$, $\theta$ = 2–10°) accounts for 43 GeV, 0.04% of the energy.
• CLD is equally blind further forward — 44 of its 48 particles at $|\eta| > 4$ also deposit nothing. The difference is not coverage but what is there: multi-TeV beam-collinear radiation at 10 TeV $\mu\mu$ versus 2 GeV total at CLD. At $\theta \approx 0.17$ mrad nothing instruments that anyway, well inside even a BeamCal.

I also checked the file: there is no forward calorimeter collection being missed, only ECal/HCal barrel+endcap and the Yoke.

The actual in-acceptance effect: status-1 particles with no simulator endpoint flags

This is what I would like a second opinion on, @alexandertuna.

In $|\eta| < 2.5$, 49 generator status-1 particles carry simulatorStatus = 0x80000000 — no endpoint bits set at all — and deposit exactly zero energy. They account for 1792 GeV of $p_T$, 2.9% of the status-1 scalar $p_T$, which is essentially the whole in-acceptance zero-deposit effect (1806 GeV). Some are hard: one carries 256 GeV of $p_T$, and there are 329 and 307 GeV pions in the barrel depositing nothing.

Using the EDM4hep bit definitions (30=CreatedInSimulation, 29=Backscatter, 28=VertexIsNotEndpointOfParent, 27=DecayedInTracker, 26=DecayedInCalorimeter, 25=LeftDetector, 24=Stopped):

• Every one of the 2452 particles that does deposit energy has at least one endpoint bit set (bits 24–27). Zero exceptions, on both detectors.
• Bit 30 (CreatedInSimulation) is set on none of them, so these are generator particles, not simulation products — consistent with generatorStatus == 1.
• The other zero-deposit particles were simulated: 87 have 0x85000000 (Stopped + DecayedInCalorimeter) and 17 have 0x89000000 (Stopped + DecayedInTracker). Together they carry only ~14 GeV of $p_T$, so they are negligible.

The reading that fits is that those 49 particles were never propagated through Geant4. A 329 GeV pion in the barrel depositing exactly zero is not physical, and it is hard to see how it acquires no endpoint flag if it was tracked.

One more oddity in the same direction: bit 31 is set on every MAIA particle and on no CLD particle. EDM4hep only defines bits up to 30, so bit 31 has no meaning in the schema. CLD's values are the clean ones (e.g. 0x05000000 = Stopped + DecayedInCalorimeter, versus MAIA's 0x85000000 for the same physical case). That smells like an artifact of the slcio → edm4hep conversion rather than anything physical, but it would be worth confirming it is harmless.

So: is a subset of status-1 particles expected to go unsimulated in this production — a generator-level filter before Geant4, or something in the BIB-overlay chain? If it is expected, then the target definition is behaving correctly and ~3% is simply not reconstructable. If it is not expected, the sample may want regenerating, and the MAIA response numbers in my earlier comments would improve accordingly.

For reference the affected particles' nearest neighbours in angle are the parent partons (generatorStatus 51/52/44 — gluons, a strange quark, an anti-top), so they look like ordinary shower products rather than anything exotic.

Numbers from 10 MAIA events and 20 CLD events.