I have been trying to work on this for a bit but it has been difficult since the code is just pure spaghetti. Really not trying to roast you for writing bad code and not documenting anything or removing things you're not using, but I think the combination of all of this is what is making it difficult for anyone to debug. Cleaning it up a lot and documenting it may make it significantly easier for yourself and others to fix. Anyways there's two things that stood out as I have been debugging it.

1: You're deploying nested integrations calls across multidimensional integrations at all levels. Every single time you evaluate any single part of the outer integration you are calling the inner integration to sweep across the entity of every single part of every corner of its dimensions. Mathematically that just does not make sense, but if it is being done intentionally (and I just do not understand why) then I would say that you can just calculate the integrals by hand since their boundaries appear to always be fixed for the whole file. Not sure if this is a helpful metaphor but its like asking why does using your fighter jet take more petrol to get to the market than your motorcycle. Both will get you there, but one of them makes more sense. Since the boundaries are always fixed it would be worth it to rewrite your code to a much simpler logic than abstractly calling fully nested multidimensional integrations if this calculation is in fact correct.

2: A bunch of the epsilons are never used. This may be causing some numerical issue if something is 0 but is not supposed to be. Though I do not know if this is the case since there's no documentation.