I'm about to derive the closed form of the integral with parameters $c_i>0, i \in \{1,2,\dots,2N\}$, which is $$ I=\int_0^{\pi/2}\prod_{i=1}^{2N}\left(\frac{\sin^2\theta}{\sin^2\theta+c_i}\right)d\theta. $$ By letting $t=\cos^2\theta$, we can further obtain $$ \int_0^{\pi/2}\prod_{i=1}^{2N}\left(\frac{\sin^2\theta}{\sin^2\theta+c_i}\right)d\theta \\ =\frac{1}{2}\int_0^1\frac{1}{\sqrt{t(1-t)}}\prod_{i=1}^{2N}\left(\frac{t}{t+c_i}\right)dt, $$ which is hard for me to solve.
However, as a special case, when there are only two different items in $c_i ,i \in \{1,2,\dots,2N\}$, namely $c_1$ and $c_2$, the integral reduces to $$ \frac{1}{2}\int_0^1\frac{1}{\sqrt{t(1-t)}}\left(\frac{t}{t+c_1}\right)^N\left(\frac{t}{t+c_2}\right)^Ndt \\ =\frac{\left(c_1c_2\right)^{-N}}{2}\int_0^1t^\left(2N-\frac{1}{2}\right)\left(1-t\right)^{-\frac{1}{2}}\left(1+\frac{t}{c_1}\right)^{-N}\left(1+\frac{t}{c_2}\right)^{-N}dt \\ =\frac{\left(c_1c_2\right)^{-N}}{2}\frac{\Gamma\left(2N+\frac{1}{2}\right)\Gamma\left(\frac{1}{2}\right)}{\Gamma\left(2N+1\right)}F_1\left(2N+\frac{1}{2};N,N;2N+1;-c_1^{-1},-c_2^{-1}\right), $$ where Eq. 9 in Appell Hypergeometric Function is used and $F_1$ is Appell's first kind of hypergeometric function.
I'm curious that in general cases where all $c_i$'s are different, how to derive $I$.
Furthermore, I wonder whether it is possible to derive $\lim_{c_i\rightarrow +\infty,i\in {1,2,\dots,2N}} I$.
$$I_m=\int_0^{\pi/2}d\theta\,\prod_{i=1}^{m}\frac{\sin^2(\theta)}{\sin^2(\theta)+c_i}\,$$ $$\theta=\sin ^{-1}\left(\sqrt{x}\right)\quad\implies\quad I_m=\frac 12\int_0^{1}\frac{dx}{ \sqrt{x(1-x)} }\,\prod _{i=1}^{m} \frac{x}{x+c_i}$$ Using partial fraction decomposition, $$I_m=\frac 12\int_0^{1}\frac{1}{\sqrt{x(1-x)}}\,\left(1+\sum _{i=1}^{m}\frac {a_i }{x+c_i}\right)dx$$ $$I_m=\frac \pi 2\left(1+\sum _{i=1}^{m}\frac{a_i}{\sqrt{c_i (c_i+1)}}\right)$$
If all $c_i=c$
$$I_m=\frac{\sqrt{\pi }}{2}\,\,\frac{\Gamma \left(m+\frac{1}{2}\right)}{\Gamma (m+1)}\frac 1{c^m}\,\,\, _2F_1\left(m,m+\frac{1}{2};m+1;-\frac{1}{c}\right)$$
and asymptotically if $c \to \infty$,
$$I_m= \frac{\sqrt{\pi }}{2}\,\frac{\Gamma \left(m+\frac{1}{2}\right)}{\Gamma (m+1)}\,\frac 1{c^m}\Bigg(1-\frac{m (2 m+1)}{2 (m+1)}\frac 1c+O\left(\frac{1}{c^2}\right)\Bigg)$$