How to extract subscripts in mathematica and apply a function on them?

For a bit of context, I am trying to evaluate an integration over Haar randm unitaries. I couldn't find any inbuilt function, so I am trying to implement a simple code for the Wingarten functtion in mathematica. I have to take average over following type of matrix elements. $$e^{-\beta (1)} u_{2,1} u_{3,1} \left(u_{1,1}\right){}^*\left(u_{2,1}\right){}^*$$ or $$e^{-2 \beta (1)} u_{2,1}^2 \left(\left(u_{1,1}\right){}^*\right){}^2$$

where the exponential pre-factor can be ignored. I am aware of the formula for such average as below,

So, my goal is to extract the values of a,b,c,d and their primed versions for each matrix elements and do the averaging using this formula. I tried doing pattern matching, but I couldn't figure out a proper way to identify and distinguish those subscripts. Besides, i am trying to have a general function e.g I may have terms that look like

$$e^{-\beta (1)} u_{2,1} u_{3,1} \left(u_{1,1}\right){}^* \left(u_{2,1}\right){}^* u_{1,2} u_{4,1}$$

or something with higher power of terms. There are formulas similar to B2 (product of delta) for those cases too, but I can't figure out a systematic way to extract these subscripts and evaluate the integral using such formula in Mathematica. I thought of defining $U$ with $U[[1,2]]$ elements instead of subscript variables, but there also I suppose pattern matching is the only way to extract the indices when encountering terms like below, $$e^{-\beta (1)} u[[2,1]] u[[3,1]] \left(u[[1,1]]\right){}^* \left(u[[2,1]]\right){}^*$$

Any help will be much appreciated. I am not very well-versed with Mathematica, please excuse me if this is a straightforward technique to implement. I couldn't generalise the existing answers on extracting subscript variables as my case involves power of variables and I want to distinguish between subscript of conjugates. Thank you in advance.