  We perform global fit to the quark Sivers function within the transverse
momentum dependent (TMD) factorization formalism in QCD. We simultaneously fit
Sivers asymmetry data from Semi-Inclusive Deep Inelastic Scattering (SIDIS) at
COMPASS, HERMES, and JLab, from Drell-Yan lepton pair production at COMPASS,
and from $W/Z$ boson at RHIC. This extraction is performed at next-to-leading
order (NLO) and next-to-next-to leading logarithmic (NNLL) accuracy. We find
excellent agreement between our extracted asymmetry and the experimental data
for SIDIS and Drell-Yan lepton pair production, while tension arises when
trying to describe the spin asymmetries of $W/Z$ bosons at RHIC. We carefully
assess the situation, and we study in details the impact of the RHIC data and
their implications through different ways of performing the fit. In addition,
we find that the quality of the description of $W/Z$ vector boson asymmetry
data could be strongly sensitive to the DGLAP evolution of Qiu-Sterman
function, besides the usual TMD evolution. We present discussion on this and
the implications for measurements of the transverse-spin asymmetries at the
future Electron Ion Collider.
