The maximal characterization of real Hardy space.

Question

Denote by $H_r^1(\mathbb{R}^n)$ the real Hardy space (integrable functions which have atomic decompositions). Let $\Phi(x)$ be a Schwarz function and define for each $f\in L^1$ the maximal function $$ Mf(x)=\sup_{\epsilon>0}|\Phi_{\epsilon}*f(x)|\quad\textrm{where}\quad \Phi_{\epsilon}(x)=\frac{1}{\epsilon^n}\Phi(\frac{1}{\epsilon}) $$ Then, as claimed in problem 6 of chapter 2, Functional Analysis by Stein, $f\in H_r^1$ iff $Mf\in L^1$. I have difficulty in proving the "if" direction. Recall the technique called Calderon-Zygmund decomposition which appears in the context of Stein's book in the construction of the atomic decomposition of $L^p$ function with compact support. There we make use of a different maximal function, the truncated Hardy-Littlewood function $f^{\dagger}$. We decompose its upper level set into dyadic cubes $\{f^{\dagger}>2^k\}=\bigcup_j Q_j^k$ and estimate meanvalue $$ m_j^k=\frac{1}{m(Q_j^k)}\left|\int_{Q_j^k}f(x)dx\right| $$ However, this argument cannot apply directly to the problem. In the definition of $Mf$, the absolute value is taken outside the integral $\Phi_{\epsilon}*f$, thus the meanvalue $m_j^k$ cannot be easily estimated as before. Some modifications must be made, though I haven't figure them out. Could anyone provide an outline of the proof?