  We study the $e^+e^-\to\gamma\omega J/\psi$ process using $11.6 ~\rm fb^{-1}$
$e^+ e^-$ annihilation data taken at center-of-mass energies from
$\sqrt{s}=4.008~\rm GeV$ to $4.600~\rm GeV$ with the BESIII detector at the
BEPCII storage ring. The $X(3872)$ resonance is observed for the first time in
the $\omega J/\psi$ system with a significance of more than $5\sigma$. The
relative decay ratio of $X(3872)\to\omega J/\psi$ and $\pi^+\pi^- J/\psi$ is
measured to be $\mathcal{R}=1.6^{+0.4}_{-0.3}\pm0.2$, where the first error is
statistical and the second systematic (the same hereafter). The
$\sqrt{s}$-dependent cross section of $e^+e^-\to\gamma X(3872)$ is also
measured and investigated, and it can be described by a single Breit-Wigner
resonance, referred to as the $Y(4200)$, with a mass of
$4200.6^{+7.9}_{-13.3}\pm3.0~{\rm MeV}/c^2$ and a width of
$115^{+38}_{-26}\pm12~{\rm MeV}$. In addition, to describe the $\omega J/\psi$
mass distribution above $3.9~\rm GeV/c^2$, we need at least one additional
Breit-Wigner resonance, labeled as $X(3915)$, in the fit. The mass and width of
the $X(3915)$ are measured to be $3926.4\pm2.2\pm1.2~\rm MeV/c^2$ and
$3.8\pm7.5\pm2.6~\rm MeV$, or $3932.6\pm8.7\pm4.7~\rm MeV/c^2$ and
$59.7\pm15.5\pm3.7~\rm MeV$, respectively, depending on the fit models. The
resonant parameters of the $X(3915)$ agree with those of the $Y(3940)$ in $B\to
K\omega J/\psi$ and of the $X(3915)$ in $\gamma\gamma\to\omega J/\psi$ by the
Belle and BABAR experiments within errors.
