  We explore data analysis techniques for signatures from heavy particle
production during inflation. Heavy particules can be produced by time dependent
masses and couplings, which are ubiquitous in string theory. These localized
excitations induce curvature perturbations with non-zero correlation functions
at all orders. In particular, Flauger et. al. 2016 has shown that the
signal-to-noise as a function of the order $N$ of the correlation function can
peak for $N$ of order $\mathcal{O}(1)$ to $\mathcal{O}(100)$ for an interesting
space of models. As previous non-Gaussianity analyses have focused on
$N=\{3,4\}$, in principle this provides an unexplored data analysis window with
new discovery potential. We derive estimators for arbitrary $N$-point functions
in this model and discuss their properties and covariances. To lowest order,
the heavy particle production phenomenology reduces to a classical Poisson
process, which can be implemented as a search for spherically symmetric
profiles in the curvature perturbations. We explicitly show how to recover this
result from the $N$-point functions and their estimators. Our focus in this
paper is on method development, but we provide an initial data analysis using
WMAP data, which illustrates the particularities of higher $N$-point function
searches.
