I presume that you are referring to the different waveforms available in LALSimulation (and thus also through, e.g., PyCBC). If so, I’m not sure exactly where the specific number of 72 comes from, since there are currently 114 implemented (see lalsimulation.NumApproximants). Nevertheless, I can answer how one picks a particular waveform model to analyze a given event. Unfortunately, there is not a publicly available guide to the various waveform models. However, you can usually find the details about the physics included in a given model by searching for its name. See also this post of mine for details about some of the waveform models used in LVK analyses, giving some details not discussed below and up to date as of the beginning of this year.
Many of the waveforms in LALSimulation are old and just kept for backwards compatibility. Some are even only applicable to the inspiral phase (e.g., the post-Newtonian waveforms TaylorF2, SpinTaylorT4, etc.), and thus should not be used to analyze a BBH, except possibly for very low-mass ones, where the merger-ringdown phase is at high frequencies where the current detectors are not very sensitive. However, as I mention in the link above, there are two standard waveforms used in LVK analyses of BBHs in both the GWTC-2.1 and GWTC-3 papers: IMRPhenomXPHM and SEOBNRv4PHM. Both of these waveforms are for quasicircular binaries and include the effects of precession as well as a selection of higher harmonics. IMRPhenomXPHM is a phenomenological model in the frequency domain, so it is relatively fast and often used for data analysis studies (which are usually carried out in the frequency domain). SEOBNRv4PHM is a time-domain effective-one-body (EOB) model which is significantly slower than IMRPhenomXPHM. However, since the two models are constructed in a significantly different way, it was good to use both of them for data analysis to assess the effects of systematic uncertainties in the modeling.
There are also some improved versions of these models now available, notably IMRPhenomXO4a, which has precession tuned to numerical relativity (see this paper) and the asymmetry between the \ell = 2, m = \pm 2 modes for precessing systems (see this paper). There is also an option to use a numerical solution of the post-Newtonian SpinTaylor equations in the IMRPhenomXPHM precession (accessed using a LALSimulation parameter—I can let you know the details if you are interested; there is a paper in preparation). There is also SEOBNRv5PHM (available through the new waveforms interface, so through ExternalPython in, e.g., lalsim-inspiral), which contains additional higher modes, compared to SEOBNRv4PHM, and also has significant improvements in speed and accuracy.
There are also separate aligned-spin and/or only dominant mode versions of these models (IMRPhenomXAS, IMRPhenomXHM, IMRPhenomXP; SEOBNRv4, SEOBNRv4HM [available in both
reduced order model (_ROM) and post-adiabatic (_PA) versions], SEOBNRv4P). Specifically, in each set these are aligned-spin only dominant mode, aligned spin with higher modes, and precessing only dominant mode (in the coprecessing frame).
There are also other models that are worth considering, notably the time-domain numerical relativity surrogate model NRSur7dq4, which directly interpolates quasicircular precessing numerical relativity simulations. It is thus limited in length by the length of the numerical simulations, so it is only applicable to sufficiently high-mass signals, where the signal in the sensitive band of current GW detectors is not that long. It is additionally limited in the mass ratios for which it provides accurate results, at most 6:1, extrapolating from the training range of waveforms with mass ratios up to 4:1. (It also is only trained on waveforms with spin magnitudes up to 0.8, but extrapolates well up to extremal spins.) This is the most accurate waveform model in its region of validity, including all the spherical harmonic modes up to \ell = 4, though it is not as fast as IMRPhenomXPHM (though not as slow as SEOBNRv4PHM; see the IMRPhenomXPHM paper for timing comparisons). There are also some aligned-spin surrogate models not included in LALSimulation—see here
There is also another EOB model, TEOBResumS (see here for references), where the latest version is not in LALSimulation. It also includes precession and higher modes in the quasicircular case, and also models eccentricity for aligned-spin systems, including higher modes.
If you want to analyze a higher-mass binary black hole, then IMRPhenomXPHM is likely the first choice, due to its speed. You might then want to consider additional models to assess waveform systematics, particularly for higher-SNR signals, especially ones with nonnegligible spins, precession, and/or higher mass ratios. For low-mass binaries, where the higher-order modes are less prominent and the analysis time is longer (due to the long signals), it might be sufficient to use IMRPhenomXP, but for a high-SNR signal, IMRPhenomXPHM will probably be necessary to obtain sufficiently accurate results. Older, fast but less accurate waveform models like IMRPhenomPv2 (single-spin precession and no higher modes) are sometimes still used when fast results are needed, particularly for low-mass systems.
For BNS and NSBH systems, there are specialized waveform models including the effects of the neutron star, notably its tidal deformability. There are versions of TEOBResumS applicable to both BNS and NSBH systems, as well as tidal versions of both the Phenom and SEOBNR models, e.g., for BNSs, IMRPhenomXP_NRTidalv2 (paper will appear on arXiv soon; the NRTidalv3 version has been developed, but is not yet in the public version of LALSimulation) and SEOBNRv4_ROM_NRTidalv2 (based on the frequency domain reduced order model of the aligned-spin SEOBNRv4 model); see here for a description of the NRTidalv2 model. There is also the SEOBNRv4T aligned-spin model and its frequency-domain surrogate SEOBNRv4T_surrogate. There are also the dedicated aligned-spin NSBH models IMRPhenomNSBH and SEOBNRv4_ROM_NRTidalv2_NSBH which include the effect of tidal disruption that is important for certain parts of the NSBH parameter space. However, for more unequal-mass NSBH systems, where there is no tidal disruption of the neutron star, the matter effects are negligible and one can analyze them with (precessing) BBH waveforms with good accuracy, as was done for GW200105_162426 and GW200115_042309 by the LVK.
As far as deciding whether an event should be analyzed with BBH, BNS, or NSBH waveforms, one considers the masses inferred from the detection pipelines and then from preliminary parameter estimation. If both masses are well above 3 solar masses, the BBH waveforms are appropriate. If either or both have significant support below 3 solar masses, then analyses with BNS and/or NSBH waveforms are called for.
I hope that this has answered your questions, but this is a large topic, so feel free to ask for more information if necessary. I cannot promise to reply quickly, but other commenters might also have useful insights here.