Extracting constant terms: is there a direct way? $\DeclareMathOperator\CT{CT}$
Let $\CT_t(f(t))$ denote the constant term of the Laurent polynomial of $f(t)$.
Define the two functions $F(x_1,\dots,x_n)$ and $G(y)$ by
$$F:=\prod_{i=1}^nx_i^{-1}(1-x_i)^{-2}\prod_{1\leq i<j\leq n}
(1-x_i-x_j)^{-1} \qquad \text{and} \qquad 
G:=n!\cdot y^{-n}e^{(n+1)y+y^2/2}.$$
I like to ask:

QUESTION. Is the following true? It would be great if there is a direct way to compare these two.
$$\CT_{x_1}\CT_{x_2}\cdots\CT_{x_n}\left(F(x_1,\dots,x_n)\right)=\CT_y\left(G(y)\right).$$

NOTE 1. The sequence on the right-hand side is available at the OEIS as A301741 with an explicit evaluation.
NOTE 2. Incidentally, we also have (a consequence of Han's formula proved here)
$$\CT_{x_1}\CT_{x_2}\cdots\CT_{x_n}\left(F(x_1,\dots,x_n)\right)=
\sum_{\lambda\vdash n}f^{\lambda}\prod_{u\in\lambda}
\frac{(n+2)^{h_u}+n^{h_u}}{(n+2)^{h_u}-n^{h_u}};$$
where $h_u$ is the hook-length of cell $u$ (in the Young diagram of $\lambda$) and $f^{\lambda}$ is the number of Standard Young Tableau of shape $\lambda$ (given by the hook-length formula).
NOTE 3. A cute analogue: let $f:=\prod_{i=1}^nx_i^{-1}(1-x_i)^{-1}\prod_{1\leq i<j\leq n}(1+x_i+x_j)$ and $g:=n!\cdot y^{-n}e^{ny-y^2/2}$. Then,
$$\CT_{x_1}\CT_{x_2}\cdots\CT_{x_n}\left(f(x_1,\dots,x_n)\right)
=\CT_y(g(y)).$$
Proof. Fedor's reasoning applies (it'd be nice to employ Richard's too)
\begin{align*}
{\rm CT}\, f&=
[x_1\ldots x_n] \prod_i (1-x_i)^{-1}\prod_{i<j}(1+x_i+x_j) \\
&=[x_1\ldots x_n]\prod_i\exp(x_i)\prod_{i<j}\exp(x_i+x_j-x_ix_j)\\
&=[x_1\ldots x_n] \exp\left(\sum x_i+\sum_{i<j}(x_i+x_j-x_ix_j)\right)\\
&=[x_1\ldots x_n]\exp\left(n\cdot S-S^2/2\right).
\end{align*}
 A: By the exponential formula, the constant term of $G(y)$ equals $\sum_w
(n+1)^{\mathrm{fix}(w)}$, where $w$ ranges over all involutions in the
symmetric group $S_n$ and $\mathrm{fix}(w)$ is the number of fixed
points of $w$. In other words, this is the number of involutions in
$S_n$ whose fixed points are colored by one of the colors
$1,2,\dots,n+1$.
The constant term of $F$ is the coefficient of $x_1\cdots x_n$ in
$\hat{F} = x_1\cdots x_n F$. Hence we can ignore all exponents greater
than 1 in the expansion of $\hat{F}$, so we want the coefficient of
$x_1\cdots x_n$ in
$$ \prod_{i=1}^n (1+x_i)^2 \prod_{1\leq i<j\leq n}
    (1+x_i+x_j+x_ix_j+x_ix_j).\ \ \ (*) $$
Given an involution $w\in S_n$ with each fixed point colored by one of
the colors $1,2,\dots,n+1$, we can obtain a monomial $x_1\cdots x_n$
as follows: if the fixed point $a$ is colored 1, choose $x_a$ from the
first factor of $(1+x_a)^2$ in (*). If $a$ is colored 2, choose $x_a$
from the second factor of $(1+x_a)^2$. If $a$ is colored $k>2$ then
consider the $(k-2)$nd factor of $\prod_{1\leq i<j\leq
n}(1+x_i+x_j+x_ix_j+x_ix_j)$ (where we order the ${n\choose 2}$
factors in some way) containing a term $x_a$. Let $x_b$ be the other
variable appearing in this factor. Say that this is the $(m-2)$nd
factor containing $x_b$. If $b$ is not colored $m$, then choose the
term $x_a$ from the factor $1+x_a+x_b+x_ax_b+x_ax_b$. If $b$ is
colored $m$, then choose the first of the two terms equal to $x_ax_b$.
It remains to account for the 2-cycles. If $(a,b)$ is a 2-cycle, then
choose the second term equal to $x_ax_b$ from the factor
$1+x_a+x_b+x_ax_b+x_ax_b$. From all remaining factors we choose the
term 1. This sets up a bijection between the
colored involutions and the monomials $x_1\cdots x_n$ appearing in the
expansion of $\hat{F}$.
A: For power series $u(x_1,\ldots,x_n),v(x_1,\ldots,x_n)$ call $u,v$ similar and write $u\sim v$ if all monomials $\prod x_i^{c_i}$ with $c_i\in \{0,1\}$ have equal coefficients in $u,v$. In other words, if $u$ is congruent to $v$ modulo the ideal generated by $x_i^2$'s. Note that this similarity respects addition and multiplucation, and that $(1-x_i)^{-1}\sim \exp(x_i)$ and $(1-x_i-x_j)^{-1}\sim 1+x_i+x_j+2x_ix_j\sim\exp(x_i+x_j+x_ix_j)$. Thus
\begin{align*}
{\rm CT}\, F&=
[x_1\ldots x_n] \prod_i (1-x_i)^{-2}\prod_{i<j}(1-x_i-x_j)^{-1}\\&=
[x_1\ldots x_n]\prod_i\exp(2x_i)\prod_{i<j}\exp(x_i+x_j+x_ix_j)\\&=[x_1\ldots x_n] \exp\left(
\sum 2x_i+\sum_{i<j}(x_i+x_j+x_ix_j)
\right)\\
&=[x_1\ldots x_n]\exp\left((n+1)S+S^2/2\right),
\end{align*}
where $S=x_1+\ldots+x_n$ (since $S^2/2\sim \sum_{i<j} x_ix_j$). Now if we expand $\exp\left((n+1)S+S^2/2\right)$ as a power series in $S$, we get $[x_1\ldots x_n]S^n=n!$ and $[x_1\ldots x_n]S^k=0$ for $k\ne n$, thus your identity.
