Proving an entire function that maps real to real, upper to upper, is a linear function

Basically, the conditions of the problem is same as this one. That is given an entire function which maps real line to real line, $\operatorname{Im}(z)>0$ to $\operatorname{Im}(f(z))>0$. I need to prove $f(z)$ is a linear function.

So far I have found two proofs, but I can understand neither of them.

The first one is the second part of zhw's answer to the question linked above. However, the answer involves Poisson integral and measure, which is too advanced for me to understand the details.

The second one is this paper. The Theorem 3 is basically a generalized problem, and the proof looks very short and elementary. However, I failed to get the idea of the proof. It looks like a Fourier series, or a poisson kernel, but not the exact definition. And there is only comparisons and I did not understand how the equality is established.

I am asking for an elementary(using complex analysis course stuff) proof, or explanation. Thanks.

Step I. Schwarz Reflection Principle implies that $f(\overline{z})=\overline{f(z)}$. Hence, if $\,\mathrm{Im}\, z>0,$ then $\mathrm{Im}\, f(z)>0.$

Step II. $f'(x)\ne 0$, for all $x\in\mathbb R$.

Assume that $f'(x_0)=0$, for some $x_0\in\mathbb R$. Then $g(x)=f(x)-f(x_0)$ satisfies whichever properties $f$ is required to satisfy in the OP. Now, $g(x_0)=g'(x_0)=0$, and thus $\frac{1}{2\pi i}\int_{|z-x_0|=r}g'/g\ge 2$, where $r>0$ is chosen so that $g(x_0-r), g(x_0+r)\ne 0$. But the total argument change of $g'/g$ can not exceed $\pi$ in each of the two half planes.

Step III. $f$ vanishes at unique $x_0\in\mathbb R$. Uniqueness follows from Step II. If $f(x)\ne 0$, for all $x\in\mathbb R$, then $f(z)\ne 0$, for all $z\in\mathbb C$, and hence, the exists an entire $g$, such that $f=\mathrm{e}^g$. Since $\,\mathrm{Im}\, z>0,$ implies $\mathrm{Im}\, f(z)>0,$ then $\,\mathrm{Im}\, z>0,$ implies $\mathrm{Im}\, g(z)\in \big(2\pi k^+, 2\pi(k^++1)\big),$ for some integer $k^+$, and similarly $\,\mathrm{Im}\, z<0,$ implies $\mathrm{Im}\, g(z)\in \big(2\pi (k^--1), 2\pi k^-\big)$, for some integer $k^-$. In which case $\mathrm{Im}\,g$ would be bounded, and so would be $g$. Without loss of generality, $x_0=0$.

Step III implies that $f$ is of the form $f(z)=(z-x_0)\mathrm{e}^{g(z)}$.

In the upper half plane, $f(z)=\mathrm{e}^{g(z)+\log (z-x_0)}$. (Here $\log (z-x_0)$ is the branch, defined in $\mathbb C\setminus (-\infty,x_0],$ with $\mathrm{Im}\,\log(z-x_0)\in (0,\pi)$, when $\mathrm{Im}\,(z-x_0)>0$.) Again, $\mathrm{Im}\,f(z)>0$, implies that there exists an integer $k$ such that $\mathrm{Im}\,\big(g(z)+\log (z-x_0)\big)\in \big(2k\pi,(2k+1)\pi\big)$, and since $\mathrm{Im}\,\log (z-x_0)\in (0,\pi),$ then $\mathrm{Im}\,g(z)\in \big((2k-1)\pi,(2k+1)\pi\big)$. Similarly, in the lower half plane, $\mathrm{Im}\,g(z)\in \big((2k'-1)\pi,(2k'+1)\pi\big)$, for some integer $k'$. Thus $g$ is bounded and hence constant, which implies that $f(z)=a(z-x_0)$, for some $a\ne 0$ real.