Abstract We study signatures of primordial non-Gaussianity (PNG) in the redshift-space halo field on nonlinear scales using a combination of three summary statistics, namely, the halo mass function (HMF), power spectrum, and bispectrum. The choice of adding the HMF to our previous joint analysis of the power spectrum and bispectrum is driven by a preliminary field-level analysis, in which we train graph neural networks on halo catalogs to infer the PNG f NL parameter. The covariance matrix and the responses of our summaries to changes in model parameters are extracted from a suite of halo catalogs constructed from the Quijote-png N-body simulations. We consider the three main types of PNG: local, equilateral, and orthogonal. Adding the HMF to our previous joint analysis of the power spectrum and bispectrum produces two main effects. First, it reduces the equilateral f NL predicted errors by roughly a factor of 2 while also producing notable, although smaller, improvements for orthogonal PNG. Second, it helps break the degeneracy between the local PNG amplitude, f NL local , and assembly bias, b ϕ , without relying on any external prior assumption. Our final forecasts for the PNG parameters are Δ f NL local = 40 , Δ f NL equil = 200 , Δ f NL ortho = 85 , on a cubic volume of 1 Gpc / h 3 , with a halo number density of n ¯ ∼ 5.1 × 10 − 5 h 3 Mpc − 3 , at z = 1, and considering scales up to k max = 0.5 h Mpc − 1 .