residual background in the same region. The mass peak regions were defined as $m_{\rm inv}=0.10$-$0.17$ and $0.48$-$0.62$ GeV/$c^2$ for $\pi^0$ and $\eta$, respectively.

Acceptance and reconstruction efficiency (efficiency hereafter) are estimated using a [geant3]{}-based [@Brun:1978fy] Monte-Carlo simulation of the PHENIX detector. The simulation was tuned to reproduce the observed mass peaks and widths of $\pi^0$ and $\eta$ in the real data. To account for the effect of underlying events the simulated mesons were embedded in real data in each centrality, then analyzed with the same methods as the real data. Final efficiencies also account for branching ratios of the analyzed decay modes and for the ERT-A trigger efficiency in the corresponding data sample.

Invariant yields of $\pi^0$ and $\eta$ are obtained as follows: $$\frac{1}{N_{\rm event}}\frac{d^2N}{2\pi p_T dp_T dy} = 
\frac{N_{\rm raw}}{2\pi p_T N_{\rm event} \epsilon_{\rm rec} 
\Delta p_T \Delta y},$$

where $N_{\rm raw}$ is the particle raw yield and $\epsilon_{\rm rec}$ is the efficiency (including acceptance and all other corrections), $N_{\rm event}$ is the number of analyzed events.

Systematic uncertainties are classified into three types. Type A represents uncertainties, that are entirely -uncorrelated; these are added in quadrature to the statistical uncertainty. Type B uncertainties are -correlated, but different from 