We present a framework for $q_T$ resummation at N$^3$LL+NNLO accuracy for arbitrary color-singlet processes based on a factorization theorem in SCET. Our implementation CuTe-MCFM is fully differential in the Born kinematics and matches to large-$q_T$ fixed-order predictions at relative order $alpha_s^2$. It provides an efficient way to estimate uncertainties from fixed-order truncation, resummation, and parton distribution functions. In addition to $W^pm$, $Z$ and $H$ production, also the diboson processes $gammagamma,Zgamma,ZH$ and $W^pm H$ are available, including decays. We discuss and exemplify the framework with several direct comparisons to experimental measurements as well as inclusive benchmark results. In particular, we present novel results for $gammagamma$ and $Zgamma$ at N$^3$LL+NNLO and discuss in detail the power corrections induced by photon isolation requirements.