You are not using the full power of the assumption that the radial limits are zero for all points except one
(not merely a.e.). Here are some facts which you might already know:
Lemma 1. There is a universal constant $c>0$ such that for all $0<\epsilon<1$ and all $t\in [-\epsilon,\epsilon]$ we have $P_{1-\epsilon}(t)\ge c\epsilon^{-1}$.
Sketch of proof. Since $t\mapsto P_r(t)$ is symmetric and decreasing in both directions, it suffices to consider $t=\epsilon$. Look at the denominator of the kernel:
$$1-2(1-\epsilon)\cos \epsilon+(1-\epsilon)^2 = 2(1-\epsilon)(1-\cos \epsilon) +\epsilon^2
= O(\epsilon^2)
$$
Since the numerator is $1-r^2\ge \epsilon$, we conclude that the function
$\epsilon\,P_{1-\epsilon}(\epsilon)$ is bounded by a positive constant from below on $(0,1]$. $\Box$
Lemma 2. Suppose $u$ is a positive harmonic function on $U$ and $\lim_{r\to 1} u(re^{i\theta})=0$
for some $\theta\in [0,2\pi)$. Then $\lim_{\epsilon\to 0}\epsilon^{-1}\mu(N_\epsilon(e^{i\theta}))=0$, where
$N_\epsilon$ stands for the $\epsilon$-neighborhood (really, a circular arc with midpoint at $e^{i\theta}$).
Proof. Let $r=1-\epsilon$. Using Lemma 1, obtain
$$
u(re^{i\theta})\ge \int_{-\epsilon}^{\epsilon} u(re^{i(\theta-t)})P_r(t)\,dt
\ge c\epsilon^{-1}\mu(N_\epsilon(e^{i\theta}))
$$
and the conclusion follows. $\Box$
The assumptions in your question imply that $\lim_{\epsilon\to 0}\epsilon^{-1}\mu(N_\epsilon(\xi))=0$ for all $\xi\in T\setminus \{1\}$.
This implies (by a covering argument) that $\mu(T\setminus \{1\})=0$, and therefore $\mu$ is a positive
point mass at $1$.
Let $\gamma_R$ denote the positively oriented contour $[-R,R]$ followed by $Re^{it},0\le t\le \pi.$ If $z$ is in the upper half plane, then for $R>|z|$ we can apply Cauchy's formula to get
$$f(z) = \frac{1}{2\pi i}\int_{\gamma _R} \frac{f(w)}{w-z}\,dw$$
Now let $R\to \infty$ and use the assumed decay rate of $f$ to see
$$f(z) = \frac{1}{2\pi i}\int_\infty^\infty \frac{f(t)}{t-z}\,dt.$$
Claim: If $f$ is replaced by $\bar f$ in the last integral, the integral equals $0.$ Assuming the claim, we have
$$f(z) = \frac{1}{2\pi i}\int_\infty^\infty \frac{f(t)+\bar f (t)}{t-z}\,dt=\frac{1}{2\pi i}\int_\infty^\infty \frac{2\text {Re}f(t)}{t-z}\,dt=\frac{1}{\pi i}\int_\infty^\infty \frac{\text {Re}f(t)}{t-z}\,dt.$$
That's the desired result. I'll leave the proof of the claim for you to ponder for now.
Best Answer
Boundedness is enough. To justify the integral representation of $u$, one usually considers the functions $u_r(z)=u(rz)$ with $0<r<1$. These are smooth on the boundary, so the formula is not an issue: $$u_r(z) = \frac{1}{2\pi} \int_0^{2\pi} u_r(e^{i\theta})\frac{1-|z|^2}{|e^{i\theta}-z|^2}\,d\theta \tag1$$ As $r\to 1$, we have $u_r(z)\to u(z)$ for every $z$ in the open unit disk $\mathbb D$, because $u$ is continuous there. So the issue is how to pass to the limit in the right hand side of (1). Since we do it with $z$ fixed, the kernel $\frac{1}{2\pi}\frac{1-|z|^2}{|e^{i\theta}-z|^2}$ is a nice bounded function of $\theta$; it creates no problems. (One could even wrap it together with $d\theta$ into a probability measure $\omega_z$ on the boundary, called the harmonic measure with respect to $z$.)
Suppose that as $r\to 1$, $u_r(e^{i\theta})$ converges for almost every $\theta$ to some function, which we can denote $u^*(e^{i\theta})$. (The asterisk is optional but it helps to keep things in order, emphasizing the distinction between a function and its limit values on the boundary.) In the presence of a uniform bound $|u_r|\le M$ the dominated convergence theorem applies; hence, $$\frac{1}{2\pi} \int_0^{2\pi} u_r(e^{i\theta})\frac{1-|z|^2}{|e^{i\theta}-z|^2}\,d\theta \to \frac{1}{2\pi} \int_0^{2\pi} u^*(e^{i\theta})\frac{1-|z|^2}{|e^{i\theta}-z|^2}\,d\theta \tag2$$ which answers your question.
Sometimes we don't have a dominating function handy. Instead, it suffices to assume that there exists $p>1$ and $M$ such that $\int |u_r(e^{i\theta})|^p\,d\theta\le M$ for all $r$. Indeed, a bounded sequence in $L^p$ with $1<p<\infty$ has a weakly convergent subsequence. Let $u^*$ be this weak limit. Since the Poisson kernel is a bounded function of $\theta $ (for a fixed $z$), integration against it a continuous linear functional on $L^p$. Thus, (2) holds for a subsequence, which is enough to get the integral representation of $u$. (Using the representation, one can show that the integrals indeed converge as $r\to 1$.)
With $p=1$ the above does not work since $L^1$ is not reflexive. But if $\int |u_r(e^{i\theta})|\,d\theta$ is bounded, then we can find a subsequence that converges in the sense of weak* convergence of measures. Then the integral representation of $u$ involves not the radial limits $u^*$, but some finite signed measure on the boundary, which need not be absolutely continuous.