diffxsec — working with differential cross sections

git clone <this repo> && cd diffxsec
pip install -r requirements.txt
python3 scripts/run_miniboone_ccqe.py    # ~40 s, writes plots/ and results.json
pytest -q                                # closure tests

A small, self-contained project that reproduces the method of the MiniBooNE ν_μ CCQE double-differential measurement (PRD 81, 092005 (2010), arXiv:1002.2680). It also shows how to build 1D, 2D and 3D differential cross sections from weighted events, the same way you would from a CAF.

1. Why the CAF has xsec but no differential xsec

A per-event xsec is a number the generator attaches to an event: the cross section of that process at that event's E_ν. A differential cross section is not an event property. It is a binned density you build by:

1. choosing observables x = (x₁, x₂, …),
2. histogramming events in x,
3. dividing by the bin volume and by the "exposure":

            N_i                       N_i  = signal events in bin i (weighted, eff.-corrected, bkg-subtracted)
dσ/dx_i = -------------------------   Φ    = energy-integrated flux [ν/cm²]
          Φ · N_targets · Δx_i         N_t  = target nucleons;  Δx_i = bin volume

• 1D: Δx = width, giving dσ/dQ² [cm²/GeV²].
• 2D: Δx = ΔT·Δcos, giving d²σ/dT dcosθ [cm²/GeV].
• 3D: Δx = ΔT·Δcos·ΔT_p, giving d³σ/… [cm²/GeV²].

Because Φ is integrated over the whole spectrum, these are flux-integrated (flux-averaged) cross sections: model ⊗ flux. Comparing to a model means folding the model with the same flux, which is what generator.py does.

The closure check: summing dσ/dx·Δx over all bins must give the flux-averaged total cross section in that phase space.

2. Layout

file	what
diffxsec/physics.py	Llewellyn Smith free-nucleon CCQE dσ/dQ², dipole form factors
diffxsec/flux.py	binned flux; from_txt() for a real table; miniboone_like() stand-in
diffxsec/generator.py	weighted-event CCQE on carbon: Fermi gas (k_F=220 MeV), E_b=34 MeV, Pauli blocking; MiniBooNE's E_ν^QE, Q²_QE
diffxsec/xsec.py	diff_xsec() for N-dimensional binning, closure, χ², MiniBooNE binning
scripts/run_miniboone_ccqe.py	the reproduction: models, (pseudo-)data, M_A fit, plots
scripts/caf_diff_xsec_template.py	the same recipe applied to a DUNE CAF (edit branch names)
tests/test_closure.py	free-nucleon sanity, MC vs analytic closure, ∫dσ = σ (run by GitHub Actions)
data/	where the real MiniBooNE flux / d²σ tables go (see data/README.md)

3. Weighted-event recipe (generator side)

Each event is drawn with E_ν ~ flux pdf and carries weight w = dσ (cm²) for its sampled kinematics. Then:

<σ>_flux    = mean(w)
dσ/dx_i     = Σ_{e∈i} w_e / N_gen / Δx_i
stat err_i  = sqrt(Σ_{e∈i} w_e²) / N_gen / Δx_i

xsec.diff_xsec(ev, ["Tmu","cosmu","Tp"], [e1,e2,e3]) handles any dimension.

4. Results

check	value
closure, no nuclear effects: MC vs analytic ⟨σ_free⟩	8.668e-39 vs 8.680e-39 (0.15%)
RFG, M_A = 1.03 GeV, flux-averaged σ per neutron	7.58e-39 cm²
RFG, M_A = 1.35 GeV	9.49e-39 cm²
MiniBooNE quoted flux-integrated total (data)	9.429e-39 cm²
fraction of the 1.35 total inside MiniBooNE's (T_μ, cosθ) bins	8.20e-39 (86%)

This reproduces MiniBooNE's headline: an RFG with the then-standard M_A≈1.03 GeV falls ~20% below the measured normalisation, while M_A≈1.35 GeV matches. Today that excess is understood as mostly multinucleon (2p2h) contributions, not a larger axial mass.

Plots in plots/:

1. 1_d2sigma_map.png: the 2D d²σ/dT_μ dcosθ_μ map in MiniBooNE binning.
2. 2_cos_slices.png: T_μ distributions in angle slices, the way the paper shows the 2D result.
3. 3_dsigma_dQ2QE.png: single-differential dσ/dQ²_QE from the same events, with the 1.35/1.03 ratio showing M_A's effect grows with Q².
4. 4_triple_diff.png: triple differential with T_p added. The proton spectrum shifts with muon kinematics, which is the information a 2D muon-only measurement integrates away (Fermi motion and Pauli blocking show up here).
5. 5_MA_scan.png: Δχ² vs M_A, full covariance (stat + 10.7% correlated normalisation) and shape-only.

5. Caveats

• Flux and data are stand-ins. The flux table (Table V) and the d²σ table (Table VI) weren't fetchable from this environment, so:
  • the flux is a Gamma-shaped spectrum tuned to MiniBooNE's quoted mean energy (788 MeV, 0–3 GeV); absolute numbers depend on this at the ~10% level;
  • the "data" are pseudo-data generated at M_A=1.35 with MiniBooNE-like statistics (146,070 events) plus one coherent 10.7% normalisation shift. The M_A fit returning 1.35 is therefore a closure test of the fitting machinery, not a measurement. In 2_cos_slices.png the points sit low relative to the 1.35 model because that random normalisation draw came out low; this is what a correlated normalisation error looks like.
• To use the real paper numbers, drop in:
  • data/miniboone_flux.txt: E_lo E_hi flux per 50 MeV bin, then run with --flux data/miniboone_flux.txt;
  • data/miniboone_d2sigma.csv: Tmu_lo,Tmu_hi,cos_lo,cos_hi,value,error with value and error in 1e-41 cm²/GeV (the paper's units). It's picked up automatically. These come from the paper appendix or the MiniBooNE data-release page.
• The real MiniBooNE error includes shape systematics (flux shape, background model, detector), so its χ² surface is much broader than the stat+norm one here. Model MC statistics are not in the covariance (χ²_min ≈ 170/126).
• The nuclear model is a simple global Fermi gas, not MiniBooNE's exact RFG with the κ Pauli-blocking scaling, and not GENIE/NuWro/NEUT.

6. Applying this to a CAF

See scripts/caf_diff_xsec_template.py. Key points:

• Choose observables you can reconstruct (lepton p and cosθ, hadronic energy, E_ν^reco, …). For atmospherics the angle is relative to the true ν direction, not the detector z axis.
• Normalise either with flux × targets (method A) or with the flux-averaged generator cross section (method B). Check what the CAF xsec branch actually means in your production before using method B.
• For systematics, compute dσ/dx for each reweight universe and build the covariance as Cov_ij = ⟨(x_i − x̄_i)(x_j − x̄_j)⟩ over universes. This is the same thing MiniBooNE did with their multisims.