We consider a parabolic system in divergence form with measurable coefficients in a non-smooth bounded domain when the associated nonhomogeneous term belongs to a weighted Orlicz space. We generalize the Calder\'{o}n-Zygmund theorem for the weak solution of such a system as an optimal estimate in weighted Orlicz spaces, by essentially proving that the spatial gradient is as integrable as the nonhomogeneous term under a possibly optimal assumption on the coefficients and a minimal geometric assumption on the boundary of the domain