GollumFit: Neutrino Likelihood Analysis

• GollumFit is an open-source framework for binned-likelihood analyses in neutrino astrophysics, integrating both physics parameters and detailed detector systematics.
• Its architecture employs analytic, gradient, and spline reweighting modes to efficiently update Monte Carlo weights and evaluate likelihoods in high-dimensional parameter spaces.
• The framework significantly speeds up analyses by eliminating on-the-fly Monte Carlo resampling, enabling fast, scalable fits for experiments like IceCube and beyond.

GollumFit is an open-source C++/Python framework specifically designed for performing binned-likelihood analyses on high-energy ($E \gtrsim 100\,\mathrm{GeV}$) neutrino telescope data. Developed with computational efficiency and experimental flexibility in mind, GollumFit provides robust treatment of both universal and experiment-specific systematic uncertainties. Its architecture supports fitting over $\mathcal{O}(10\!-\!40)$ model parameters—physics and nuisance systematics—without requiring regeneration of Monte Carlo samples at each likelihood evaluation, addressing key analysis bottlenecks in large-scale neutrino telescope experiments such as IceCube and offering extensibility to other observatories.

1. Scope, Design Motivation, and Parameterization

GollumFit targets binned-likelihood maximization problems in next-generation neutrino telescopes. Its core justification arises from the need to maximize a complex likelihood over dozens of highly correlated parameters—physics parameters of interest and experiment- or detector-specific nuisance parameters—while avoiding the prohibitive computational cost of on-the-fly Monte Carlo resampling. The framework integrates two parameter classes:

• Common parameters: applicable to any neutrino telescope, including atmospheric flux uncertainties, cross-section attenuation factors, hadronic yield uncertainties with full covariance, cosmic-ray spectrum parameters, and broken-power-law astrophysical flux models.
• IceCube-specific parameters: encompassing digital optical module (DOM) efficiency, hole-ice properties, and “SnowStorm” ice model parameters that encode amplitude and phase with covariance, and additional experiment operational systematics.

Typical use-cases include diffuse astrophysical vs. atmospheric neutrino flux separation and searches for sterile neutrinos, but the codebase is explicitly modular, permitting straightforward extension to other experiments such as KM3NeT or Baikal-GVD by supplying new MC and systematic parameterizations (Collaboration, 4 Jun 2025).

2. Framework Architecture and Workflow

The central object is gollumfit::GollumFit, constructed from three primary inputs: SteeringParams (binning configurations, parallelization flags, systematics file paths), DataPaths (data, MC, and systematics derivatives file locations), and gollumfit::Priors (parameter prior definitions). Each parameter—physics ($\theta_k$) or systematics ($\eta_j$)—is mapped to a parameter-specific “weighter” object, a templated C++ class responsible for updating MC event weights ($w_e$) in response to parameter changes.

Three modes for this reweighting are supported:

1. Analytic: closed-form formulas, e.g., global scaling,
2. Gradient: linear expansion using precomputed MC–vs–parameter derivatives,
3. Spline: multidimensional interpolation for non-analytic systematics.

At run time, GollumFit applies weight corrections per MC event according to the current parameter vector $(\vec{\theta}, \vec{\eta})$, bins the reweighted events to compute $\mu_i(\vec{\theta},\vec{\eta})$, and evaluates the likelihood. Minimization leverages the L-BFGS-B algorithm, with automatic differentiation support via the PhysTools package for efficient analytic gradient computation.

A code-level workflow is illustrated below:

from gollumfit import SteeringParams, DataPaths, Priors, GollumFit steering = SteeringParams( ebins=40, cbins=40, elims=(1e2,1e8), clims=(-1,1), nthreads=4 ) data_paths = DataPaths( mc_file="/path/to/icecube_mc.root", data_file="/path/to/icecube_data.root", sys_dir="/path/to/spline_and_gradients/" ) priors = Priors() priors.add_gaussian("convNorm", 1.0, 0.2) priors.add_gaussian("astroNorm", 1e-18, 0.5e-18) priors.add_gaussian("astroDeltaGamma", 0.0, 0.3) priors.add_uniform("astroPivot", 4.0, 8.0) gf = GollumFit(steering, data_paths, priors) init_guess = priors.get_central_values_dict() result = gf.MinLLH(init_guess) print("Best-fit parameters:", result.best_fit) print("Min NLL:", result.minimum) print("Eval count:", result.n_evaluations) hist, ebins, cbins = gf.GetExpectation(result.best_fit) import matplotlib.pyplot as plt plt.pcolormesh(ebins, cbins, hist.T) plt.xlabel("log10(E/GeV)"); plt.ylabel("cos(theta)") plt.title("Best-fit MC expectation") plt.show()

(Collaboration, 4 Jun 2025)

3. Binned-Likelihood Formalism and Optimization

GollumFit implements a likelihood formalism designed for binned event data. Letting $i$ denote analysis bins (typically two-dimensional in $\log_{10}E$ vs. $\cos\theta$, with optional topological categorization), $x_i$ the observed event count in bin $i$, and $$\mu_i(\vec{\theta},\vec{\eta}) = \sum_{e \in i} w_e(\vec{\theta}, \vec{\eta})$$ the expected count from the reweighted MC, the framework maximizes the posterior

$$\mathcal{L}_{\rm post}(\vec{\theta}, \vec{\eta}) = \mathcal{L}(\vec{\theta}, \vec{\eta}) \times \Pi(\vec{\theta}, \vec{\eta}),$$

where

$$\mathcal{L}(\vec{\theta},\vec{\eta}) = \prod_{i} \mathcal{L}_{\rm eff}(\mu_i(\vec{\theta},\vec{\eta}), x_i)$$

is the product of effective Poisson likelihoods $\mathcal{L}_{\rm eff}$—as proposed by Argüelles, Schneider, and Yuan (JHEP 06 (2019) 030)—that explicitly accounts for MC sample statistical uncertainty. Priors are incorporated multiplicatively:

$$\Pi(\vec{\theta},\vec{\eta}) = \prod_k \pi_{\rm phys}(\theta_k) \times \prod_j \pi_{\rm nuisance}(\eta_j),$$

with typically Gaussian or uniform forms (e.g., $\pi(\eta)\propto\exp[-(\eta-\eta_0)^2/(2\sigma^2)]$).

Minimization proceeds by targeting the negative log-posterior:

$$-\ln\mathcal{L}_{\rm post} = -\sum_i\ln\mathcal{L}_{\rm eff}(\mu_i, x_i) - \sum_k\ln\pi(\theta_k) - \sum_j\ln\pi(\eta_j)$$

Hypothesis testing (e.g., for sterile-neutrino models) uses a log-likelihood ratio test statistic:

$$\mathrm{TS} = -2 \ln \frac{\mathcal{L}_{\rm post}(\mathrm{null})} {\mathcal{L}_{\rm post}(\mathrm{best\!-\!fit})}$$

(Collaboration, 4 Jun 2025).

4. Systematic Uncertainty Handling and Parameter Classes

Parameter coverage in GollumFit explicitly distinguishes “common” telescopic and IceCube-specific systematics, as summarized below.

Parameter Class	Quantity/Type	Reweighting Mode
convNorm	Flux normalization (universal)	Analytic
Atmospheric	$\rho_{\rm atm}$, $\sigma_{K\!-\!{\rm Air}}$	Analytic / Spline
Hadronic Yield	DAEMONFLUX ($K$, $\pi$, $p$, $n$), 10 params	Analytic / Covariance
CR Spectrum	GSF$_{1..6}$, 6 params, full covariance	Analytic
Astro Flux	$\Phi_{\rm astro}$, $\Delta\gamma$, pivot	Analytic / Spline
Prompt Flux	Normalization, $\nu/\bar{\nu}$ ratio	Analytic
Cross-section	$\nu$ Att., $\bar{\nu}$ Att.	Analytic
IceCube: DOM eff.	Spline	Spline
IceCube: Hole ice	Spline	Spline
SnowStorm ice	5 amplitude + 4 phase (covariant grad.)	Gradient / Spline

All parameters are treated as free (with priors), coupled directly to the weight recalculation chain. This unified treatment allows for high-dimensional correlated systematics, and onboard analytic/gradient/spline mapping as detailed in the code’s metadata. (Collaboration, 4 Jun 2025)

5. Performance Characteristics and Scaling

GollumFit achieves computational efficiency through several mechanisms. Null-injection performance (Figure 1, (Collaboration, 4 Jun 2025)) demonstrates recovery of all prior-centered values to within $<1\,\sigma$ for 38 simultaneously floating nuisance parameters, even from randomly offset initializations. Likelihood evaluation time scales linearly with the MC event count.

• With FastMC compression at $k=0.25$, approximately $0.8 \times 10^6$ “meta-events” can be fitted with $\mathcal{O}(10~\mathrm{ms})$ per likelihood evaluation on an Intel Xeon Platinum core.
• Uncompressed full samples ($1.3 \times 10^7$ events) require roughly $150\,\mathrm{ms}$ per evaluation.
• Complete fits with $814,847$ FastMC-compressed events and 38 parameters converge in tens of seconds per run on one core (Figure 2).
• Removal of certain systematics classes produces fit-time speedups of 10–50%.

The framework supports multicore parallelization for further speedup of event loop and reweight operations. (Collaboration, 4 Jun 2025)

6. Additional Features and Utilities

GollumFit provides a comprehensive ancillary feature set:

• FastMC: Combines MC events into statistically equivalent “meta-events” to accelerate reweight loops. User-tunable compression through $k\sim0.25$ (API: ConstructFastMC, WriteCompact).
• GetExpectation / GetExpectationEvents: Retrieve binned or unbinned weighted MC under arbitrary parameters, enabling pseudo-data, pull, and diagnostic studies.
• EvalLLHGradient: Provides analytic or automatic-differentiation gradients for external samplers, supporting methods such as nested sampling or Hamiltonian Monte Carlo.
• Containerized release: Full Docker/Singularity images for reproducible deployment and HPC portability.
• Extensible dimensionality: Support for additional bin axes (e.g., right ascension, event topology) directly in SteeringParams.

(Collaboration, 4 Jun 2025)

7. Relationship to Established Neutrino Analysis Frameworks

Comparison with contemporaneous frameworks reveals distinct design choices in GollumFit:

• ROOT/RooFit and BAT offer general statistical models but lack native support for event-by-event reweighting in high-dimensional systematics spaces characteristic of modern neutrino observatories.
• IceTray is central to IceCube’s event simulation/reconstruction pipeline, but does not provide gradient-enabled, high-level binned-likelihood optimization or FastMC compression.
• Golem and the JUNO analysis pack are tightly coupled to individual experiments; by contrast, GollumFit offers modular, open-source treatment of both common (flux/model) and detector-specific uncertainties.
• For published diffuse-flux neutrino analyses in IceCube (e.g., Phys. Rev. Lett. 125 (2020) 141801, IceCube:2024uzv), GollumFit has accelerated workflow pipelines and remains, as of publication, the only publicly available package of comparable scope.

In summary, GollumFit occupies a unique position among statistical analysis frameworks for neutrino telescope data, serving as a fast, flexible, binned-likelihood fitter capable of handling $\mathcal{O}(30\!-\!40)$ highly correlated nuisance parameters, supporting comprehensive systematic uncertainty propagation, and providing scalable computational tools for the next generation of neutrino astrophysics research (Collaboration, 4 Jun 2025).