Instead of treating this as a technical problem, I would recommend that you consider the underlying conceptual issue: QuadGK uses adaptive quadrature, which is tricky to differentiate (it may be differentiable locally, but not continuous).
I would recommend using simple Gaussian quadrature instead, eg
In particular, if your interval is finite, you should use Gauss-Legendre with fixed nodes/weights, and transform according to θ₁, θ₂ etc. Then integration is simply a dot product, which is differentiable.