I am interested in calculating the p.d.f of the random variable
$$ Y = a + b^\intercal z + c \|z\|_2^2, $$ where $a,c$ are constants, $b \in\mathbb{R}^n$, and $z \sim \mathcal{N}(0,I)$ is a standard normal distribution, where $z \in \mathbb{R}^n$. It is known that the norm squared of a standard normal distribution is a $\chi^2$ distribution, which has p.d.f
$$ f(z) = \frac{1}{2^{n/2}\Gamma(n/2)}z^{n/2 - 1} e^{-z/2}, $$ and $b^\intercal z$ is a normal distribution such that $b^\intercal z \sim \mathcal{N}(0,b^\intercal b)$. I can always rewrite this in terms of a standard normal distribution as $b^\intercal z = \sqrt{b^\intercal b} \ \tilde{z} = \|b\|_2 \ \tilde{z}$, where $\tilde{z} \sim \mathcal{N}(0,1)$ is a standard normal random variable (scalar). This then gives
$$ Y = a + \|b\|_2 \tilde{z} + c \|z\|_2^2, $$ which is the sum of a standard normal random variable and $\chi^2$ random variable with $n$ degrees of freedom. My question is what would be the p.d.f of $Y$, i.e., what is the p.d.f of the sum of a $\chi^2$ distribution and a standard normal distribution?