  We present a high-precision QED calculation, with 0.1% theoretical accuracy,
of two photon production in $e^+ e^-$ annihilation, as required by more and
more accurate luminosity monitoring at flavour factories. The accuracy of the
approach, which is based on the matching of exact next-to-leading order
corrections with a QED Parton Shower algorithm, is demonstrated through a
detailed analysis of the impact of the various sources of radiative corrections
to the experimentally relevant observables. The calculation is implemented in
the latest version of the event generator BabaYaga, available for precision
simulations of photon pair production at $e^+ e^-$ colliders of moderately high
energies.
