We consider an elliptic self-adjoint first order pseudodifferential operator acting on columns of m complex-valued half-densities over a connected compact n-dimensional manifold without boundary. The eigenvalues of the principal symbol are assumed to be simple but no assumptions are made on their sign, so the operator is not necessarily semi-bounded. We study the spectral function, i.e. the sum of squares of Euclidean norms of eigenfunctions evaluated at a given point of the manifold, with summation carried out over all eigenvalues between zero and a positive λ. We derive a two-term asymptotic formula for the spectral function as λ tends to plus infinity. In doing this we establish that all previous publications on the subject give incorrect or incomplete formulae for the second asymptotic coefficient. We then restrict our study to the case when m=2, n=3, the operator is differential and has trace-free principal symbol, and address the question: is our operator a massless Dirac operator? We prove that it is a massless Dirac operator if and only if the following two conditions are satisfied at every point of the manifold: a) the subprincipal symbol is proportional to the identity matrix and b) the second asymptotic coefficient of the spectral function is zero.