I’m a bit confused about some of your questions–for instance, PyRing already performs the primary statistical calculations for constraining QNMs from the data, and question 2 seems purely to be about plotting, not GW data analysis, but perhaps I am misunderstanding it. However, if you haven’t already read the ringdown analysis section (VIII A) in the latest LVK testing GR paper, then this might provide some of the information you’re looking for. You might also have a look at the corresponding section (VII A) in the previous LVC testing GR paper, particularly since it analyzes GW150914, which it appears that you may be considering, given the paper you link. Both of these papers provide some code for plotting, etc. as part of their data releases (here for the latest paper), which could potentially be useful for you.
Additionally, to reply to 1. and 3. more explicitly (though most of this information is in those papers), for 1., it appears that you are interested in combining together the constraints from multiple events. Here there are two methods used, either simply multiplying likelihoods, assuming that the values for the deviation parameters are the same for all events, or inferring the distribution of deviations in the detected population hierarchically, as described in this paper. Unfortunately, I don’t think that the code to do the hierarchical analysis is public, but you could ask Max Isi if he’d be willing to share it if you’re interested.
For 3., if you’re analyzing a real signal, there won’t be a single GR prediction for the frequency and damping time (of the least-damped QNM), just a posterior distribution of predictions obtained from the full IMR analysis of the signal, which I presume is what you’re referring to. You can see the median and 90% credible intervals for GW150914 and many other events in Table IX of the previous testing GR paper mentioned above and for the more recent events in Table XII of the latest testing GR paper. The posterior samples for these quantities are released in the data releases, e.g., here for the previous testing GR paper including the GW150914 results.
If you want to compute these yourself, then you can start from the posterior samples for the final (redshifted) mass and spin provided by the LVK (linked from GWOSC, e.g., for the latest results on GW150914, the GWTC-2.1 samples). You then can apply either the classic fits for the frequency and damping time for the 220 mode (or another mode) from Appendix E of Berti, Cardoso, and Will, or the more recent ones from London and Fauchon-Jones (implemented in the positive package), or compute the QNM frequencies directly (solving the eigenvalue problem) using, e.g., Stein’s qnm package.