The question also asks how to show that $\mathcal{H}\mu$ is defined almost everywhere. I don't have a reference for this, but the standard proof can be modified without too much difficulty. The maximal operator $\mathcal{H}^*\mu(x)\equiv\sup_{\epsilon > 0}\lvert\mathcal{H}_\epsilon\mu(x)\rvert$ is weak (1,1) continuous,\left\lvert\left\lbrace x\colon\mathcal{H}^*\mu(x) > \lambda\right\rbrace\right\rvert\le C\lVert\mu\rVert/\lambda,for all $\lambda > 0$ and a fixed constant $C$ (I'm using $\lvert\cdot\rvert$ to denote the Lebesgue measure). I'm working from the notes Interpolation, Maximal Operators, and the Hilbert Transform, by Michael Wong (Theorem 8.7). These notes look at the case where $\mu$ is absolutely continuous, but the proof carries across to the general case with no big changes.
By weak continuity, if there exists measures $\mu_n$ with $\lVert\mu-\mu_n\rVert\to0$ such that $\mathcal{H}\mu_n$ all exist almost everywhere, then $\mathcal{H}\mu$ exists almost everywhere.If $\mu$ is absolutely continuous with differentiable density, then $\mathcal{H}\mu$ exists everywhere in the standard way. As the differentiable functions are dense in $L^1$, this extends to all absolutely continuous measures.Only the case for singular measures $\mu$ remain. So, there exists a measurable $S\subseteq\mathbb{R}$ with zero Lebesgue measure with $\mu(S^c)=0$. Then we can choose compact sets $K_n\subseteq S$ with $\mu(S\setminus K_n)\to0$. Note that the measures $\mu_n\equiv 1_{K_n}\cdot\mu$ are supported on the compact sets $K_n$ and $\lVert\mu-\mu_n\rVert\to0$. It follows that $\mathcal{H}\mu_n(x)$ is defined everywhere outside of $K_n$ (in fact, $\mathcal{H}_\epsilon\mu_n(x)$ is a constant function of $\epsilon$ for $\epsilon$ small enough that $B_\epsilon(x)\cap K_n=\emptyset$). So, $\mathcal{H}\mu_n$ is defined almost everywhere and, hence, so is $\mathcal{H}\mu$.
Showing that the two definitions agree almost everywhere is easy! Using the truncated transform $$\mathcal{H}_\epsilon\mu(x)=\frac1\pi\int_{\lvert y-x\rvert > \epsilon}\frac{d\mu(y)}{x-y}$$ then, by definition, $\mathcal{H}\mu(x)=\lim_{\epsilon\to0}\mathcal{H}_\epsilon\mu(x)$ for all $x$ at which the limit exists. Convolve the identity $$\Re\left(\frac{1}{x+ih}\right)=\frac{x}{x^2+h^2}=\int_0^11_{\left\lbrace\lvert x\rvert > h\sqrt{t/(1-t)}\right\rbrace}\frac1x\,dt$$ with $\frac1\pi d\mu$ to obtain, $$\Re\left(\mathcal{H}\mu\right)(x+ih)=\int_0^1\mathcal{H}_{h\sqrt{t/(1-t)}}\mu(x)\,dt.$$ The integrand on the right hand side tends to $\mathcal{H}\mu(x)$ as $h\to0$, whenever this is defined, so bounded convergence gives $\Re(\mathcal{H}\mu)(x+ih)\to\mathcal{H}\mu(x)$.