After reading Bateman & Diamond's and Pollack's books, I realize that Hooley's extension of Brun's pure sieve is the way to go. In the pure sieve, Brun derives an upper bound sieve via the following inequality
$$
\sum_{\substack{d|(n,P(z))\\\nu(d)\le\ell}}\mu(d)\ge
\begin{cases}
1 & (n,P(z))=1 \\
0 & \text{otherwise}
\end{cases}
$$
for even $\ell$ and $P(z)=\prod_{p<z}p$, where $\nu(d)$ counts the number of distinct prime factors in $d$. Compared to Brun's original attempt to improve his sieve, Hooley's trick is much simpler. In particular, he defines a sequence of squarefree numbers $P_1,P_2,\cdots,P_r$ such that $P(z)=P_1P_2\cdots P_r$ and a sequence of even numbers $\ell_1,\ell_2,\cdots,\ell_r$. This allows for a more reflexible range of parameters to be chosen:
$$
\sum_{d|(n,P)}\mu(d)=\prod_{1\le j\le r}\sum_{d|(n,P_j)}\mu(d)\le\prod_j\sum_{\substack{d|(n,P_j)\\\nu(d)\le\ell_j}}\mu(d)
$$
Plugging this into $N_z$, we have
\begin{aligned}
N_z
&\le\prod_j\sum_{\substack{d|P_j\\\nu(d)\le\ell_j}}\mu(d)|\mathcal A_d| \\
&=X\underbrace{\prod_j\sum_{\substack{d|P_j\\\nu(d)\le\ell_j}}{\mu(d)b(d)\over d}}_Q+\underbrace{\sum_{\substack{d_1,d_2,\cdots,d_r\\d_j|P_j,\nu(d)\le\ell_j}}|r(d_1d_2\cdots p_r)|}_R
\end{aligned}
To continue, we went on estimating $Q$ and $R$.
An upper bound for $Q$
For convenience, we let $p^-(d)$ denote the least prime divisor of $d$. Then
Lemma: Suppose $m$ is squarefree, $h(d)$ is multiplicative. Then for any $k\ge0$ there is
$$
\sum_{d|m}h(d)=\sum_{\substack{d|m\\\nu(d)\le k}}h(d)+\sum_{\substack{d|m\\\nu(d)=k+1}}h(d)\prod_{p<p^-(d)}(1+h(p))
$$
Proof. See Lemma 12.3 of Bateman & Diamond.
Applying this lemma to $Q$, we have
\begin{aligned}
\sum_{\substack{d|P_j\\\nu(d)\le\ell_j}}{\mu(d)b(d)\over d}
&=\sum_{d|P_j}{\mu(d)b(d)\over d}-\sum_{\substack{d|P_j\\\nu(d)=\ell_j+1}}{\mu(d)b(d)\over d}\prod_{p<p^-(d)}\left(1-{b(p)\over p}\right) \\
&=\underbrace{\prod_{p|P_j}\left(1-{b(p)\over p}\right)}_{V_j}+\sum_{\substack{d|P_j\\\nu(d)=\ell_j+1}}{b(d)\over d}\prod_{p<p^-(d)}\left(1-{b(p)\over p}\right)\le V_j+K_j
\end{aligned}
where the second equal symbol follows from the fact that $\mu(d)=-1$ in the second sum when $\ell_j$ is even, and $K_j$ is defined as follows
$$
K_j=\sum_{\substack{d|P_j\\\nu(d)=\ell_j+1}}{b(d)\over d}
$$
If we were to write $d=p_1p_2\cdots p_{\ell_j+1}$ such that $p_1>p_2>\cdots>p_{\ell_j+1}\ge2$. Then
\begin{aligned}
K_j
&=\sum_{\substack{1\le k\le\ell_j+1\\p_k|P_j\\p_k\text{ decreasing}}}{b(p_1)b(p_2)\cdots b(p_{\ell_j+1})\over p_1p_2\cdots p_{\ell_j+1}} \\
&\le{1\over(\ell_j+1)!}\sum_{\substack{1\le k\le\ell_j+1\\p_k|P_j}}{b(p_1)b(p_2)\cdots b(p_{\ell_j+1})\over p_1p_2\cdots p_{\ell_j+1}} \\
&={1\over(\ell_j+1)!}\left(\sum_{p|P_j}{b(p)\over p}\right)^{\ell_j+1}
\le{(\log V_j^{-1})^{\ell_j+1}\over(\ell_j+1)!}
\end{aligned}
where the last inequality follow from the elementary fact that for each $0\le y<1$
$$
y\le y+{y^2\over2}+{y^3\over3}+\cdots=\log(1-y)^{-1}
$$
Finally, we define
$$
V(z)=V_1V_2\cdots V_r=\prod_{p<z}\left(1-{b(p)\over p}\right)
$$
Then $Q$ can be written as
\begin{aligned}
Q
&\le V(z)\prod_j\left(1+{1\over V_j}{(\log V_j^{-1})^{\ell_j+1}\over(\ell_j+1)!}\right) \\
&\le V(z)\exp\left\{\sum_j{1\over V_j}{(\log V_j^{-1})^{\ell_j+1}\over(\ell_j+1)!}\right\}\triangleq V(z)\exp E
\end{aligned}
Now, we can start choosing our parameters $P_j$ and $\ell_j$.
Hooley's choice of $P_j$ and $\ell_j$ and the Iwaniec condition
For any $z\ge2$, we will reach one if we continue taking square roots. As a result, we can partition $[2,z)$ into following intervals:
$$
[2,z)=[z^{1/2},z)\cup[z^{1/4},z^{1/2})\cup\cdots\cup[z^{2^{-r}},z^{2^{1-r}})\cup[2,z^{2^{-r}})
$$
This indicates that we can let $P_j$ be the product of all primes within $[z^{2^{-j}},2^{2^{1-j}})$ for $1\le j\le r-1$ and $P_r$ be the product of all primes in $[2,z^{2^{1-r}})$. To estimate $E$, we first require $b(p)$ to satsify a weaker Iwaniec condition:
$$
\prod_{a\le p<b}\left(1-{b(p)\over p}\right)^{-1}\le A\left(\log b\over\log a\right)^\kappa\tag{$\Omega$}
$$
Then we have for $1\le j\le r-1$
$$
V_j^{-1}=\prod_{z^{2^{-j}}\le p<2^{2^{1-j}}}\left(1-{b(p)\over p}\right)^{-1}\le A(\log2)^\kappa=L=\text{const.}
$$
and
$$
V_r^{-1}=\prod_{2\le p<z^{2^{1-r}}}\left(1-{b(p)\over p}\right)^{-1}\le A\left(2^{1-r}\log z\over\log2\right)^\kappa
$$
Since ($\Omega$) does not allow us to bound $V_r^{-1}$ above by a constant, we can set $\ell_r=+\infty$ to eliminate this term from $E$, yielding
$$
E\le L\sum_{1\le j\le r-1}{(\log L)^{\ell_j+1}\over(\ell_j+1)!}
$$
Since $\ell_j$ needs to be even, we can plainly set $\ell_j=2j$ for $1\le j\le r-1$ to get
$$
E\le L\sum_{j\ge1}{(\log L)^{2j+1}\over(2j+1)!}\le L\sum_{m\ge0}{(\log L)^m\over m!}=L^2
$$
This indicates that $E$ is bounded above by a constant, so $Q\ll V(z)$. Now, we move our attention to the remainder term $R$.
Estimating the remainder $R$
Under the assumption $|r(d)|\le b(d)$, we have
$$
R\le\sum_{k|P_r}b(k)\prod_{1\le j\le r-1}\sum_{\substack{d|P_j\\\nu(d)\le2j}}b(d)=\prod_{p<z^{2^{1-r}}}(1+b(p))\prod_{1\le j\le r-1}\color{blue}{\sum_{\substack{d|P_j\\\nu(d)\le2j}}b(d)}
$$
For blue part, we may use the fact that $d\le z^{2^{2-j}j}$ when $\nu(d)\le2j$ to deduce
$$
\sum_{\substack{d|P_j\\\nu(d)\le2j}}b(d)\le z^{2^{2-j}j}\prod_{2^{2^{-j}}<p\le z^{2^{1-j}}}\left(1+{b(p)\over p}\right)
$$
Since for any $0\le y<1$, $1+y\le1+y+y^2+y^3+\cdots=(1-y)^{-1}$, so the produce in the above expression is bounded by $V_j$. Thus, we obtain the following bound for $R$
$$
R\le z^\alpha\prod_{p\le z^{2^{1-r}}}(1+b(p))\prod_{1\le j\le r-1}V_j^{-1}
$$
where $\alpha=\sum_{j\ge1}2^{2-j}j$. Since $V_r<1$ and $b(p)<p$, we get
$$
R\le{z^\alpha\over V(z)}\prod_{p\le z^{2^{1-r}}}p={z^\alpha\over V(z)}\exp[\vartheta(z^{2^{2-r}})]
$$
Since $z^{2^{1-r}}\to0$ as $r\to+\infty$, we can set $r$ large enough so that $z^{2^{1-r}}<\log\log X$. Eventually, using Chebyshev's bounds and ($\Omega$), we see that there exists a constant $\beta>0$ such that $R\le z^{\alpha}(\log X)^\beta$, so we can conclude our task with the following result:
Theorem (Brun): Under the aforementioned hypotheses, there exists fixed $0<A<1$ such that
$$
N_z\ll X\prod_{p<X^A}\left(1-{b(p)\over p}\right)
$$
Bibliography
Bateman, P. T., & Diamond, H. G. (2004). Analytic number theory: An introductory course. World Scientific.
Halberstam, H., & Richert, H.-E. (1974). Sieve methods. Academic Press.
Pollack, P. (2009). Not always buried deep: A second course in elementary number theory. American Mathematical Society.
Best Answer
For practical purposes, one only needs to care about the case where $\alpha$ is even, which indicates that
$$ n+\alpha\equiv0\pmod p\Rightarrow \begin{cases} \text{1 solution} & p=2 \\ \text{2 solutions} & p>2 \end{cases}\Rightarrow\omega(p)\text{ solutions}. $$
By the Chinese remainder theorem, we know for sure that $n+\alpha\equiv0\pmod d$ has $\omega(d)$ solutions whenever $d$ is squarefree. If we were to divide the range $[1,x]$ into intervals of size $d$. Then
$$ \left||A_d|-\omega(d)\left\lfloor\frac xd\right\rfloor\right|\le\omega(d). $$
Because $|y-\lfloor y\rfloor|\le1$, we see immediately that
$$ |R_d|=\left||A_d|-{\omega(d)\over d}x\right|\le2\omega(d), $$
which still satisfies the conditions for Brun's sieve (In fact, $|R_d|\le L\omega(d)$ for any fixed $L>0$ will work).