Robust Calibration of Non-Perturbative Models with History Matching

Abstract: We apply, for the first time, Bayes Linear Emulation and History Matching to the calibration of non-perturbative models in Monte Carlo event generators. In contrast to the usual approach of "Monte Carlo tuning", History Matching does not result in best-fit plus ellipsoidal parameter uncertainty estimates but instead identifies all parameter space regions that are consistent with data. This approach leads to a systematic and robust quantification of parametric uncertainties in the models, especially in those challenging cases where different, possibly disjoint, regions of parameter space deliver similar results, which are usually not properly treated with current methodology. We highlight the power of this method with the hadronisation models available through Sherpa: the built-in cluster fragmentation Ahadic and string fragmentation through an interface to Pythia.

• The paper introduces Bayes Linear emulation with History Matching to calibrate SHERPA’s AHADIC and PYTHIA 8 hadronisation models across 432 LEP observable bins while explicitly accounting for experimental, model, and Monte Carlo uncertainties.
• History Matching reduces the viable parameter volume by 5.2×10^-4 for AHADIC and 1.5×10^-13 for PYTHIA 8, revealing strongly correlated and multimodal regions that conventional single-point tuning could miss.
• The resulting ensembles quantify parametric and model-choice uncertainty, show comparable global predictions between models, and identify remaining tensions in heavy-flavour fragmentation, identified-hadron yields, and baryon production.

Motivation and context

Calibrating the non-perturbative components of Monte Carlo event generators (MCEGs) remains a methodological weak point in high-energy physics phenomenology. Standard "tuning" workflows — typically polynomial surrogates à la Professor combined with optimisation around a $\chi^2$ minimum — locate a single best-fit point and attach an ellipsoidal uncertainty region around it. This approach has two well-known deficiencies: it does not reliably discover multiple, potentially disjoint, regions of parameter space that describe the data comparably well, and it therefore risks underestimating parametric uncertainty, particularly when calibration results are extrapolated to new kinematic domains.

The paper under review applies Bayes Linear (BL) emulation and History Matching (HM) — well established in epidemiology, climate science, galaxy formation, and reservoir modelling — for the first time to non-perturbative model calibration in particle physics. The target application is hadronisation in the SHERPA event generator, covering both the built-in cluster fragmentation model AHADIC and the Lund string-fragmentation model of PYTHIA 8, accessed through SHERPA's interface. Because both models run on identical perturbative input (NLO matrix elements for $e^+e^-\to q\bar q + \{0,1\}$ jets matched to SHERPA's dipole shower), the authors can quantify not only the parametric uncertainty of each model but also the uncertainty attributable to the choice of hadronisation model itself.

The History Matching framework

HM inverts the logic of optimisation: rather than seeking "good" parameter points, it systematically removes "implausible" ones. An implausibility measure

$$I(x)^2 = \frac{\left(E[f(x)-z]\right)^2}{Var[f(x)-z]}$$

compares simulator output $f(x)$ to observations $z$ relative to the combined uncertainty from observational error $e$ and model discrepancy $\epsilon(x)$. Points with $I(x)$ exceeding a threshold (typically $I=3$, following Pukelsheim's rule) are ruled out. Crucially, the framework makes explicit the structural relation $f(x) = z + e + \epsilon(x)$, forcing the practitioner to specify both observational and model-discrepancy variances rather than absorbing them into an effective $e^+e^-\to q\bar q + \{0,1\}$0.

To make the search computationally feasible in a roughly 20-dimensional parameter space — where direct exploration would require over half a million simulator runs — the simulator is replaced by BL emulators of the form $e^+e^-\to q\bar q + \{0,1\}$1, requiring only first- and second-moment specifications and updated via the Bayes Linear update equations rather than a full Gaussian-process posterior. Emulator uncertainty enters the implausibility measure, and successive "waves" of the algorithm — design, simulate, emulate, rule out — progressively shrink the non-implausible region until emulator uncertainty is subdominant. Two properties of HM are emphasised as methodologically decisive: (i) outputs need only be emulated selectively at each wave, since a point ruled out by one output stays ruled out; and (ii) if the non-implausible space empties, this signals a genuine model–data conflict rather than silently returning the "least bad" posterior, as MCMC- or ABC-based methods would.

Application setup

The calibration uses $e^+e^-\to q\bar q + \{0,1\}$2 observable bins from eleven LEP-era analyses (ALEPH, DELPHI, OPAL, SLD) plus PDG hadron multiplicities, accessed via Rivet. AHADIC contributes $e^+e^-\to q\bar q + \{0,1\}$3 parameters and the PYTHIA 8 string model $e^+e^-\to q\bar q + \{0,1\}$4. Each wave runs $e^+e^-\to q\bar q + \{0,1\}$5 parameter sets, with $e^+e^-\to q\bar q + \{0,1\}$6 (AHADIC) or $e^+e^-\to q\bar q + \{0,1\}$7 (PYTHIA 8) used for emulator training and the remainder held out for diagnostics. Only a subset of outputs — $e^+e^-\to q\bar q + \{0,1\}$8–$e^+e^-\to q\bar q + \{0,1\}$9 per wave, selected via principal variable analysis so that at least $$I(x)^2 = \frac{\left(E[f(x)-z]\right)^2}{Var[f(x)-z]}$$0 of each observable's variability is represented — is emulated at any wave; over the full match, $$I(x)^2 = \frac{\left(E[f(x)-z]\right)^2}{Var[f(x)-z]}$$1 (AHADIC) and $$I(x)^2 = \frac{\left(E[f(x)-z]\right)^2}{Var[f(x)-z]}$$2 (PYTHIA 8) distinct outputs were emulated. Emulator quality is verified through standardised diagnostics checking predictive agreement, implausibility classification consistency, and absence of systematic bias.

Two uncertainty contributions deserve emphasis. First, the finite Monte Carlo statistics of the generator are treated as an $$I(x)^2 = \frac{\left(E[f(x)-z]\right)^2}{Var[f(x)-z]}$$3-independent contribution to the model discrepancy; the authors note that in just over $$I(x)^2 = \frac{\left(E[f(x)-z]\right)^2}{Var[f(x)-z]}$$4 of observations this stochastic variability exceeds the observational error, and in extreme cases is five times larger — ignoring it would spuriously rule out viable parameter space. Second, a model discrepancy of $$I(x)^2 = \frac{\left(E[f(x)-z]\right)^2}{Var[f(x)-z]}$$5–$$I(x)^2 = \frac{\left(E[f(x)-z]\right)^2}{Var[f(x)-z]}$$6 of the observation magnitude is included, reflecting the inherent imperfection of the hadronisation models.

Structure of the non-implausible parameter space

The compression achieved by HM is substantial. For AHADIC, the non-implausible volume falls from $$I(x)^2 = \frac{\left(E[f(x)-z]\right)^2}{Var[f(x)-z]}$$7 to $$I(x)^2 = \frac{\left(E[f(x)-z]\right)^2}{Var[f(x)-z]}$$8 over three waves; for PYTHIA 8, five waves reduce the volume from $$I(x)^2 = \frac{\left(E[f(x)-z]\right)^2}{Var[f(x)-z]}$$9 to $f(x)$0 — a reduction of more than thirteen orders of magnitude. Additional verification waves reduced the volume by only $f(x)$1 and $f(x)$2, respectively, justifying termination. At the late waves, more than $f(x)$3 of emulators had low uncertainty, and wave members matched at least $f(x)$4 of the $f(x)$5 observations.

The most consequential finding concerns the topology of the surviving space. Two-dimensional projections reveal strongly correlated and frequently multi-modal structures: for AHADIC, BARYON_FRACTION versus P_QQ1_by_P_QQ0 forms a "banana-shaped" distribution with two distinct maxima at opposite ends, and similar bimodality appears for BARYON_FRACTION versus P_QS_by_P_QQ_norm. For PYTHIA 8, the well-known aLund–bLund correlation is confirmed, while aLund–Sigma is tightly constrained and unimodal. The authors state the implication directly: a tune that settles in one mode would miss the equally viable region around the other — precisely the failure mode of single-point calibration that HM is designed to expose.

Observable predictions and model comparison

Propagating the final-wave parameter sets through the generator yields uncertainty envelopes that constitute the paper's estimate of non-perturbative parametric uncertainty. For inclusive observables and global event shapes (thrust, thrust major/minor, charged-particle multiplicity), the two models give nearly indistinguishable predictions, with envelopes comparable to or moderately larger than experimental uncertainties. Differences emerge in heavy-flavour fragmentation — AHADIC shows larger variation and better coverage of the $f(x)$6-quark fragmentation function, while PYTHIA 8 slightly undershoots the data centrally — and in identified-hadron yields: AHADIC overproduces $f(x)$7 and $f(x)$8 by roughly $f(x)$9 and underestimates $z$0, $z$1, and $z$2 yields by approximately $z$3, while PYTHIA 8 overshoots $z$4 by $z$5 and undershoots $z$6 comparably. Notably, the final AHADIC wave exhibits large spread in $z$7, $z$8, and $z$9 baryon yields despite being matched to the corresponding data — a residual tension the authors leave unresolved.

The reduced $e$0 distributions of all $e$1 final-wave members for both models peak near $e$2 with no member exceeding $e$3 (AHADIC) or $e$4 (PYTHIA 8), confirming that the entire non-implausible set consists of tunes of comparable quality. In the distribution tails — particularly the differential 2-, 3-, and 4-jet rates — parameter uncertainties exceed the combined measurement and simulation statistical error, and the authors state plainly that tighter constraints there require both more accurate data and higher simulation statistics.

Uncertainty propagation and limitations

The non-implausible region is explicitly not a posterior: all surviving points carry equal weight, and no posterior mode exists. The paper therefore devotes attention to downstream use — constructing pseudo-posteriors via importance weighting against an inflated Gaussian proposal, subselecting representative ensembles for expensive secondary analyses, and, if a single "best" point is demanded, using internal discrepancy to retain the parametric uncertainty the match affords.

The authors identify several open issues. The per-bin emulation of histogram observables ignores their underlying density structure; a parametrisation through hyperparameters, or a fully multivariate emulator with explicit between-bin covariances, could reduce the number of waves but requires careful prior specification and uncertainty propagation. The stochastic variability of the generator, while accounted for, remains inescapable and limits constraining power in distribution tails; systematic enhancement of rare configurations is suggested but not implemented. Finally, some LEP observations were excluded on grounds of expert judgement regarding reliability, and validating those exclusions against the final non-implausible space remains undone.

Conclusion

This work demonstrates that Bayes Linear emulation with History Matching is computationally tractable for full-scale hadronisation calibration, compressing the parameter spaces of two distinct non-perturbative models by up to thirteen orders of magnitude against 432 LEP observable bins while exposing multi-modal, strongly correlated viable regions that single-point tuning would miss. The resulting non-implausible ensembles provide parametric and model-choice uncertainty estimates for hadronisation in SHERPA. The method's practical adoption in phenomenological analyses, however, depends on on-the-fly reweighting algorithms for non-perturbative parameters, which the authors identify as an area of active development, and on extending the approach to the underlying-event and remnant-fragmentation components relevant for hadronic collisions.