Silicon planar realistic digitization for varying sensor thicknesses with multilple scattering and realistic timing smearing implementation

[angirar]

This PR includes all the up-to-date changes from the sensor thickness studies for multiple scattering implementation in long clusters and realistic timing resolution as studied by Juliette in these slides in Marlin-version.

[madbaron]

I think _currentLayer is 0 here: LoadGeometry() runs from init() before any hit is processed and _currentLayer is initialised to 0 in the constructor. So the thickness of layer 0 selects one global _DigitizedBins for all layers, which defeats the stated purpose (per-sensor-thickness bins)?
If I didn't misread this, I think that since _DigitizedBins is a single member vector filled once at init, per-layer bins probably need a vector<vector> indexed by layer, or the selection moved into ChargeDigitizer() where _currentLayer is valid.

[madbaron]

_layerThickness is std::vector, compared for exact equality against double literals. float(0.075) promotes to 0.07500000298..., so == 0.075 is false.
As a consequence, every layer falls through to the generic 4-bit set.
I think this needs a tolerance (std::fabs(t - 0.075) < eps) or a lookup keyed on rounded microns?

[jules6077]

This is a good catch and should be implemented.

[madbaron]

The CMake option ZSegmented replaces defaulted to ON, is this the intended behaviour?

[madbaron]

Is this used anywhere? Sorry if I missed it!

[madbaron]

The units are inconsistent with the change below

[madbaron]

can this be taken from mcp->getCharge()?

[jules6077]

I think this gives the coulomb charge. We want the absolute value of the particle charge number for this quantity.

[madbaron]

yes!

[madbaron]

I have prepared a branch ("package_cleanup", https://github.com/MuonColliderSoft/MuonCVXDDigitiser/tree/package_cleanup, based on this PR) that takes care of the comments above and a few other things.

If you agree to the changes, you should be able to fast forward to the new tip and get everything in this PR (I can also open another PR and obsolete this one if you prefer - the contribution history will be preserved no matter what).

Here's an AI-assisted summary of the changes, with a couple of questions.

---

1. Per-thickness charge digi is never picked up

As submitted, the per-thickness charge digitisation never activates. Three independent defects stack up in src/MuonCVXDDigitiser.cc:

1. The selection runs before any hit exists. The _DigitizedBins assignment sits at the end of LoadGeometry(), which is called from processRunHeader(), always before the first hit is decoded. _currentLayer is therefore still 0 from the constructor, so _layerThickness[_currentLayer] is layer 0's thickness and one global table is chosen for every layer.
2. The thickness comparison can never be true. _layerThickness is std::vector<float>; float(0.075) widens to 0.07500000298…, which is != the double literal 0.075. Same for 0.1 / 0.2 / 0.4. Every layer falls through to the generic 4-bit set even for layer 0.
3. isVertex is read uninitialised. init() assigns exactly one of isVertex / isInnerTracker / isOuterTracker and leaves the other two indeterminate.

Individually each is small; together they mean the feature the PR exists to add is inert.

Change: _DigitizedBins becomes vector<vector<double>> — one table per layer, built at geometry load (where _currentLayer is meaningless) and indexed by _currentLayer inside ChargeDigitizer() (where it is valid). Matching is keyed on lround(thickness_mm * 1000) microns via a std::map instead of float equality. Per your answer, the guard stays isVertex, so barrel and endcap are covered.

Note for the PR description: production MAIA_v0 is 50 µm, which was not one of the four tabulated thicknesses. I have copied over the generic 4 bit bins to 50 µm, but please let me know if you'd prefer I undid this change.

---

2. A second bug, found while making the slice thickness configurable

The multiple-scattering stepping loop ran one slice too many at every tabulated thickness. while (z_traveled < sensorT) accumulated z_traveled += z_segment in floating point, and because _layerThickness is a float, 50 µm / 5 µm evaluates to 10.000000149 — so after 10 steps z_traveled was still fractionally short and a whole 11th slice ran.

sensor	slices intended	actually run
50 µm	10	11
75 µm	15	16
100 µm	20	21
200 µm	40	41
400 µm	80	81

Each extra slice also advanced pos by another pathL_segment, so the track was propagated past the sensor face — about 10% too far at 50 µm — before the exit point was derived from it.

This is worth knowing before comparing against the Marlin-version numbers in the slides: if those were produced with this code, they carry the extra slice.

Change: step an integer number of slices, n = ceil(sensorT / sliceT × (1 − 1e-6)), with the slice thickness exposed as a new MSSliceThickness parameter (default 0.005 mm). The small relative tolerance absorbs the float representation noise while still rounding a genuine 10.4 up to 11.

---

3. Silent behaviour changes that would have caught users out

Both of these change results for existing steering files with no error and no warning.

ZSegmented default flipped from on to off

Moving the compile-time switch to a runtime parameter is right, but the default changed:

	master	PR as submitted
MuonCVXDDigitiser ladder length	isVertex && isBarrel (geometry-derived)	ZSegmented, default false
DetElemSlidingWindow sensor offset	CMake ZSEGMENTED, default ON	zSegmented, default false

For the vertex barrel this switches _layerLadderLength from lengthSensor × sensorsPerLadder to lengthSensor, which changes GetPixelsInaColumn() and the local-y mapping — i.e. reconstructed hit positions move.

Change: ZSegmented becomes a tri-state int — -1 auto, 0 off, 1 on — defaulting to -1, where auto reproduces master's geometry-derived logic exactly. You keep the ability to force it on for a tracker barrel or off for a test, without the silent flip.

TimeSmearingSigma silently changed meaning

It stopped setting the smearing width but survived as an on/off gate (if (_timeSmearingSigma > 0)), so any steering file with a tuned value silently switched to the new realistic model.

Change: new TimeSmearingModel selector — 0 none, 1 constant sigma from TimeSmearingSigma, 2 realistic (default). I picked realistic-by-default, so this is a deliberate break that needs a release note: disabling smearing is now TimeSmearingModel=0, and the old constant behaviour is TimeSmearingModel=1.

IsBarrel as a user parameter

MuonCVXDRealDigitiser asked the user to declare whether the sub-detector is a barrel, while MuonCVXDDigitiser::init() already derives it from SubDetectorName. Two processors that can disagree about the same detector is a footgun.

Change: derive it the same way in MuonCVXDRealDigitiser and drop the parameter. DetElemSlidingWindow correspondingly goes back to a single resolved flag rather than the two the PR added.

---

4. The larger modelling change

This one goes beyond the PR's own scope and is the main thing worth discussing.

The problem

FindLocalPosition() takes the track direction from mcp->getMomentum() — the MC particle's momentum at its production vertex — rather than hit->getMomentum(), the momentum recorded at the hit. Measured on muonGun_pT_0_50_sim_0.slcio:

angle between vertex and hit momentum	median	p75	p90	max
VertexBarrel hits	3.5°	65.9°	134.7°	174.5°

A quarter of hits had the assumed incidence direction wrong by more than 65°. That direction sets tanx/tany, which set the trail tilt, the track length and hence the deposited charge. It also sets _currentParticleMomentum, which enters the fluctuation model through β².

Separately, hit->getPathLength() — the true path in the sensitive volume — was printed at DEBUG6 and never used; the trail always spanned the full thickness. In this sample:

pathLength / straight-line crossing	p05	median	p95
	0.962	1.012	1.609

4.7% of hits travel <90% of the assumed length (stopped, produced inside, clipped a corner); 14.1% travel >110%. The 1.2% median excess is the scattering and curvature Geant4 already simulated.

The change

Four commits:

1. Local direction and momentum from the hit, falling back to the MC particle only when the hit carries no momentum. This one is unconditional — the vertex momentum is simply the wrong quantity, in either mode.
2. Trail geometry from getPathLength() — a straight segment centred on the hit position along the local direction, capped at the full-thickness crossing and clipped to the slab (ProduceSignalPoints() computes drift as halfThickness − z and would go negative otherwise). Property-tested over 200k random configurations: no point escapes the sensor.
3. dEmean anchored to hit->getEDep() instead of the _energyLoss × trackLength parametrisation, with the parametrisation kept as a fallback. SampleFluctuations() still supplies the segment-to-segment Landau structure; only the mean is pinned to truth.
4. Both of those gated behind a new ResimulateIonisation flag, so the re-simulation chain they replaced is still reachable for optimisation studies rather than being deleted.

Introduced ResimulateIonisation (int, default 0) selects where the ionisation trail comes from:

	0 — Geant4 (default)	1 — re-simulation
trail length	hit->getPathLength(), falling back to the crossing	full sensor crossing along the local direction
trail extent	centred on the hit, capped at the crossing, clipped to the slab	entry face → exit face
dE per segment	hit->getEDep() / N	EnergyLoss × trackLength / N
1/n² top-up to the G4 deposit	active	skipped

Mode 1 is not simply "master again": commit 1 above still applies, so the trail follows the direction at the hit rather than at the production vertex, and the top-up loop is skipped so the total stays independent of the Geant4 deposit. What it does restore exactly is master's trail geometry and its parametrised energy — segment count, track length, segment depth and every ionisation-point coordinate reproduce the pre-gating formulas bit-for-bit (checked algebraically and numerically over a range of incidence angles, including the MaxTrackLength cap).

DoMultipleScattering is orthogonal and applies in both modes — the deflected direction still orients the trail in either — but in mode 0 it no longer sets the trail's length, which now comes from Geant4.

The effect

Charge closure against the true Geant4 deposit:

	sum(reco EDep)/sum(sim EDep)	per-event median
master	0.824	0.883
+ direction & path length	0.922	1.147
+ anchored dEmean	0.936	1.037

Master was systematically losing ~12% of the deposited energy. The residual 6% is threshold losses and charge digitisation, which is expected. The bottom row is now the default (ResimulateIonisation=0); the flag exists so this progression can be walked back towards the parametrisation, on the same build, when tuning it.

Direct measurement of _eSum / true EDep after the existing 1/n² top-up loop: median 0.9985, 68% of hits below 1.0, but mean 1.134 with a tail to 2.86 — a small number of hits where a large Landau excursion overshoots the true deposit and the one-sided loop cannot correct it. Whether that tail should be capped is a physics question worth your view; it lives exactly where delta rays live and may be legitimate.

Consequence for multiple scattering

Geant4TrackerWeightedAction averages the momentum over the step, so the local direction already contains all scattering up to mid-sensor. What re-simulating MS still adds is only the residual half-thickness wiggle:

p	sensor	θ₀	true lateral displacement	fraction of a 25 µm pixel
0.1 GeV	400 µm	10.6 mrad	2.45 µm	10%
0.1 GeV	50 µm	3.4 mrad	0.10 µm	0.4%
1 GeV	400 µm	0.71 mrad	0.16 µm	0.7%
1 GeV	50 µm	0.23 mrad	0.007 µm	0.03%

So intra-sensor MS is a sub-pixel effect except in the thick-sensor, low-momentum corner — which is precisely the study this PR is for, so the feature earns its place; it just shouldn't be expected to move anything in production. DoMultipleScattering still defaults to off, and with ResimulateIonisation=0 the re-simulated deflection no longer feeds the trail length at all — set both to 1 to exercise the chain the PR wrote.

---

5. Open questions

1. Should the MS slice loop stay at all? ResimulateIonisation=1 is where it does the most work — but the gating does not touch how it is implemented, and the implementation has a problem: it rotates the whole ionisation trail by the final accumulated angle about the hit position, rather than following the kinked path. That overestimates the MS-induced transverse extent by exactly √3 (x·θ₀ instead of x·θ₀/√3) and destroys the physical displacement–angle correlation (ρ = √3/2). A closed-form PDG sampling would be exact for a slab, need no slices, and would make the slice-thickness and log-term issues disappear:

   θ_plane = z₂·θ₀
   y_plane = z₁·x·θ₀/√12 + z₂·x·θ₀/2      (z₁, z₂ independent standard normals)
   

2. The PDG log term is evaluated per slice, not for the total thickness. Accumulating n slices this way underestimates θ₀ by 12% (50 µm, 10 slices) to 21% (400 µm, 80 slices). Was that the intent, and were the numbers in the slides produced this way?
3. Should re-simulation mode still top up to the Geant4 deposit? The 1/n² loop at the end of ProduceIonisationPoints() tops the sampled charge up to hit->getEDep(). It is active on master, and it is currently gated off in mode 1 so that the parametrisation can be judged on its own. Leaving it on would make mode 1 a closer replica of master; turning it off, as now, makes the mode genuinely independent of Geant4's energy. This is a one-line choice (if (anchorToG4)), and your view on which is more useful for the studies would settle it.
4. The pos / entry / exit bookkeeping in the MS block is currently dead — _currentLocalPosition is reset to the original position and _currentEntryPoint/_currentExitPoint are never read anywhere. Was the kinked path meant to displace the ionisation trail, or only to rotate it? The answer decides whether that code becomes live or gets deleted.

---

6. Smaller fixes included

• Reproducibility: MS used std::mt19937_64 seeded from std::random_device, outside Marlin's control — verified that two runs of the same job with the same RandomSeed produced different physics. Now uses CLHEP RandGauss::shoot, like every other smearing stage.
• Charge number: q_charge was hard-coded to 1; now |q| from the MC particle. Neutrals are skipped rather than scattered — with z = 0 the PDG log term is log(0) and θ₀ evaluates to 0 × −inf = NaN, which would have propagated silently into hit positions.
• Unit expressions: SampleFluctuations arguments and hcharge used * dd4hep::keV / dd4hep::GeV and / dd4hep::GeV where the correct conversions are / dd4hep::MeV and * dd4hep::GeV. Harmless under dd4hep's default TGeo units (GeV = 1) but wrong by 10³ under DD4HEP_USE_GEANT4_UNITS. Worth stressing that the hcharge form is a no-op in the default build: the 1/n² top-up loop behaved identically before and after, so none of the charge-closure numbers above are attributable to it.
• Per-hit debug block (deltaR, incidentTheta, beta) was computed for every hit regardless of log level; now guarded by streamlog_level(DEBUG6). Also abs() → std::fabs() on a double, and an unreachable branch removed.
• Timing constants: float → double throughout TimeSmearer (the override members are double), the override sentinel > -1.0 → >= 0. (a -0.5 override would previously have been accepted as a negative sigma), and the sigma_TDC comment corrected — 0.025/√12 is a 25 ps LSB, not the 25 ns the comment claimed.
• Naming: zSegmented/isBarrel → ZSegmented/IsBarrel to match every other parameter. Both are new in this PR, so no shipped steering file breaks.
• Whitespace, indentation, an unused #include <array>, a dead commented block containing a $$ syntax typo, and a no-op double c = 1; in the MS formula.

---

7. Validation

All runs: muonGun_pT_0_50_sim_0.slcio (MAIA_v0, 10 events), VertexBarrelCollection, RandomSeed=1234567, LayerIDs = 0 1 2 4 6.

test	result
master vs branch, TimeSmearingModel=1, before the modelling change	bit-for-bit identical, all 59 hits, every field
per-thickness bins, MAIA 50 µm, 4 bits	identical — no table for 50 µm, generic set used
per-thickness bins, 100 µm geometry, 4 bits, barrel	differs in EDep only, 4/4 hits
per-thickness bins, 100 µm geometry, 4 bits, endcap	differs in EDep only, 8/8 hits
per-thickness bins, 100 µm geometry, 5 bits	identical — tables are 4-bit only
MS reproducibility, same seed, two runs	identical after the fix; differed before it
trail clipping, 200k random configurations	no ionisation point outside the sensor
ResimulateIonisation=1 trail vs. the pre-gating formulas	segment count, track length, segment depth and every point coordinate agree to ~1e-17 mm, over a range of incidence angles and both signs of u_z
full build in the container image	both libraries link; no compiler diagnostic from MuonCVXDDigitiser.cc
Marlin -x with the built library	ResimulateIonisation registered, default 0

The first row is the important one: every cleanup and default-restoration commit was proven to change nothing observable before the modelling work began. After the modelling change that no longer holds by design — digitised charge and position move for roughly a quarter of hits, all in the low-momentum and partial-traversal population.

Caveat on the gating: ResimulateIonisation=1 has been verified by construction — the formulas were shown to reproduce the pre-gating trail exactly, the code compiles and the parameter registers — but the digitiser has not yet been run end-to-end in that mode on the validation sample. The numbers in section 4 are all ResimulateIonisation=0. A repeat of the first-row comparison with the flag set is the obvious next check.

Steering-file gotcha found during validation: MAIA_v0's vertex barrel has 5 layers whose cellID layer values are 0, 1, 2, 4, 6. Supplying LayerIDs = 0…7 is silently accepted, and layerMapping() then returns 4 and 6 as indices into 5-element vectors. No error, no crash, plausible-looking output. See follow-up E2 below.

---

8. Deliberately not included

Seven pre-existing issues on master, unrelated to this PR, proposed as a separate PR against master so they don't inflate this review. (There is no E5 — a numbering gap, nothing missing.)

	issue
E1	_layerLadderLength = 2 * z_layout.zHalfSensitive is missing * dd4hep::cm / dd4hep::mm — a cm value stored in an mm field. Inner/outer tracker only.
E2	layerMapping() returns −1 for an unknown ID, and FindLocalPosition() then indexes _layerHalfPhi[-1] before the existing _currentLayer == -1 guard is reached.
E3	_cutOnDeltaRays is passed by non-const reference into SampleFluctuations(), which writes to it — the member drifts across calls. Fixable entirely inside MyG4UniversalFluctuationForSi.
E4	phi = std::atan(y/x) should be atan2 (debug output only).
E6	TimeDigitizer() caches numBins and binWidth in function-local statics, so a job running a barrel and an endcap instance with different TimeMaximum would digitise the second with the first's bin width.
E7	MuonCVXDRealDigitiser.h doc block is stale: four registered parameters undocumented, three documented under names that no longer exist.

---

Appendix — commits

31bb9e2 gating re-sim path
44772ab anchor mean edep to G4
48a225c build ionisation trail from G4 path length
fd0c504 take local direction instead of MC vertex
45ec86c make slice thickness configurable
a53470b use mcpart q instead of 1
d9f3ee4 use randgauss to ensure reproducibility
d9f886b select digi bins by sensor thickness
4abe0e1 derive barrel/vertex flags from SubDetectorName
8ad7347 add selector for timesmearingmodel
896ba25 make zsegmented tri-state defaulting to geo behaviour
43678ea match init order to decl order
e6d5eac default init dettypes
1c3cf6d get rid of c
243ee62 indentation and cleanup
32a4b3f Unify naming to MuonCVXDRealDigitiser
3cb63dd avoid double>float roundtrips, fix sentinel
c08a0be guard hit debug and fix abs/atan
48c8ac0 fix unit conversion

New or changed processor parameters:

parameter	processor	change
ZSegmented	both	bool → tri-state int, default -1 (auto)
IsBarrel	RealDigitiser	removed — derived from SubDetectorName
TimeSmearingModel	Digitiser	new, default 2 (realistic)
TimeSmearingSigma	Digitiser	now only used when TimeSmearingModel=1
MSSliceThickness	Digitiser	new, default 0.005 mm
ResimulateIonisation	Digitiser	new, default 0 — Geant4 path length and EDep; 1 restores the re-simulated trail and the EnergyLoss parametrisation