  We extend the formalism to calculate non-Gaussianity of primordial curvature
perturbations produced by preheating in the presence of a light scalar field.
The calculation is carried out in the separate universe approximation using the
non-perturbative delta N formalism and lattice field theory simulations.
Initial conditions for simulations are drawn from a statistical ensemble
determined by modes that left the horizon during inflation, with the
time-dependence of Hubble rate during inflation taken into account. Our results
show that cosmic variance, i.e., the contribution from modes with wavelength
longer than the size of the observable universe today, plays a key role in
determining the dominant contribution. We illustrate our formalism by applying
it to an observationally-viable preheating model motivated by non-minimal
coupling to gravity, and study its full parameter dependence.
