Introduction: The positron range (PR) effect is a significant factor limiting spatial resolution in positron emission tomography (PET), particularly for high-resolution systems and non-standard isotopes. Objective: This study introduces a novel analytical model to accurately and rapidly describe positron range distributions (PRd) for various PET radioisotopes. Approach: The proposed model explicitly incorporates the Coulomb repulsion effect, the multi-branch nature of certain β⁺ emitters, and the scaling of PR with electronic density. To minimise bias, we used a histogram-free statistical method to derive the cumulative PRd from Monte Carlo (MC) simulated annihilation datasets, avoiding arbitrary histogram binning. A comparative analysis of PR estimates was conducted across three major MC radiation transport algorithm packages: PENELOPE (via PenEasy/PeneloPET), GEANT4 (via GATE), and EGS5 (via PHITS), revealing notable discrepancies between codes, versions, and input configurations, especially at short distances from the source. Main results: The new analytical model demonstrated an excellent reproduction of the simulated data for isotopes including C-11, N-13, O-15, F-18, Cu-64, Ga-68, Rb-82, and I-124, achieving in general coefficients of determination (R²) greater than 0.995 and mean absolute percentage errors (MAPE) < 20%. Compared to previous methods, our model provides a more accurate description of PRd at short and long distances and offers improved R² values. Significance: This work provides a robust framework for generating accurate annihilation point spread function (aPSF) kernels, facilitating improved PR correction in quantitative Nuclear Medical imaging and supporting research with diverse radioisotopes.

Positron range (PR) degrades spatial resolution in positron emission tomography (PET), with its effect being radionuclide- and medium-dependent. Current accurate PR simulators are computationally intensive and cannot be used in clinical time frames. This work presents the Hybrid Monte Carlo (HMC) PR simulator, a GPU-accelerated condensed-history code in CUDA Fortran for 3D annihilation imaging in heterogeneous voxelised geometries derived from CT images and $\mu$-maps. Pre-computed look-up interaction tables (LUTs) covering 24 reference-tissue materials and support for arbitrary ${\beta }^{+}$ emitters up to 3.55 MeV end-point energies (⁶⁸Ga, ⁸²Rb, ¹²⁴I, etc.) allow both GPU and CPU implementations. Maximum speeds reach $4.77·{10}^5$ (Intel i9-14900KF) and $3.32·{10}^8$ (NVIDIA RTX 5080) histories/s, surpassing prior CPU-only PenEasy 2024 by up to 4 orders of magnitude. A fixed step of 0.5 mm was chosen for LUT generation with a uniform $\pm 10\%$ variation each HMC step that allowed reducing aliasing effects. These values are consistent with our previous cross-code benchmark for the isotopes studied. Close agreement with PenEasy 2024 is demonstrated: point-source FWHM and FWTM residuals using ¹⁸F, ⁶⁸Ga, ¹²⁴I and ⁸²Rb are below 0.3 mm at ${10}^7$ histories, and voxelwise $R^2=0.995$, $\mathrm{SSIM}=0.995$ and PSNR $\approx 36.4$ dB are achieved simulating ⁸²Rb uptake from a human-scale [¹⁸F]FDG myocardium acquisition, while decreasing the simulation time from almost 9 h to roughly 1 s using HMC GPU. The resulting annihilation maps enable PET resolution modelling and patient-specific PR correction within clinically practical time frames.