Chi square fit, fitting data to a straight line: Incomplete gamma function as goodness-of-fit estimation

On page 661 of the book "Numerical Recipes in C:The art of scientific computing", we fit a set of N data points $$(x_{i}, y_i)$$ to a straight line model $$y(x)=a+bx.$$ We do this by minimizing the chi-square merit function with fitting parameters a,b, as in eq. (15.2.2) of the book $$\chi²(a,b)=\sum_{i=1}^N\big(\frac{y_i-a-bx_i}{\sigma_i}\big)^2.$$ Here $$\sigma$$ is the uncertainity for each $$y_i$$ measurement. The goodness-of-fit of the data is then calculated using the incomplete gamma function as in equation eq. (6.2.3) $$Q=\text{gammaq}\left(\frac{(N-2)}{2}, \chi^2/2\right)= \\ \frac{1}{\Gamma((N-2)/2)}\int_{\chi^2/2}^{\infty}e^{-t} t^{(N-2)/2-1}dt,$$ where $$\Gamma$$ is the complete gamma function. Since this is way more abstract to me than other fitting routines, I want to ask if someone has a deeper understanding of why we use the incomplete gamma function gammaq as the goodness-of-fit estimation. What connects the $$\chi^2$$ fitting with the gammaq function, e.g. how do we know that the incomplete gamma function returns us a gof estimation for chi-square fitting?

$$\chi^2$$ distribution is a special case of the Gamma distribution : a $$\chi^2$$ distribution with $$n$$ degrees of freedom is the same as a Gamma distribution with shape parameter $$n/2$$ and scale parameter $$2$$ : $$\chi^2_n \sim \text{Gamma}(n/2,2)$$, or equivalently $$\chi^2_n/2 \sim \text{Gamma}(n/2,1)$$ by the scaling property of the Gamma distribution.
So the quantity $$Q$$ is just the tail probability of a $$\chi^2$$ distribution with $$N-2$$ degrees of freedom: $$Q = P(X > \chi^2_{obs})$$ where $$\chi^2_{obs}$$ is the observed value and $$X \sim \chi^2_{N-2}$$. It tells you what is the probability to observe a similar or larger value than $$\chi^2_{obs}$$ under the hypothesis that the model is correct. (you have $$N-2$$ degrees of freedom because of the two free parameters in the fit)