Integrating over Euler Angles I have a $6\times6$ matrix having its elements being functions of Euler's angles (ZXZ rotation scheme) representing a tensor physical property. To find the average of the tensor property, I need to integrate this matrix over space in a spherical coordinate system. I found at many places that two of the Euler's angles are same as the spherical coordinates $(\theta, \phi)$, but I am not able to imagine how to deal with the third Euler's angle in this integration. Would appreciate any help.
 A: I believe that the orientational average of a function $f(\phi, \theta, \chi)$ of the three Euler angles $\phi$, $\theta$, $\chi$ may be written
$$
\langle f\rangle = 
\frac{1}{8\pi^2}
\int_0^{2\pi} d\phi \, 
\int_0^{\pi} \sin\theta d\theta \, \int_0^{2\pi} d\chi \, 
f(\phi, \theta, \chi)
$$
Here $\theta$ is the angle of rotation about the $x$ axis,
which is applied in between the two rotations about the $z$ axes. You can see an example of this being used
in this paper
(in that paper the intermediate rotation is about the $y$ axis,
but I don't believe that affects the volume elements appearing
in the integral).
[EDIT following OP comments]
I'm just going to make a few more remarks here,
following the comments made on my original answer.
(Generally, extended discussion in the comments is discouraged, and I can see the danger of going far beyond the original question).
From your comments, and your other question on Math SE, it seems that you wish to orientationally average a $6\times 6$ stiffness tensor for a crystal
with some symmetry. 
Certainly, you can use symmetry to restrict the range of integration. You need to identify symmetry elements (rotations) of the crystal, let's say $\mathcal{R}$, which satisfy $f(\Omega)=f(\mathcal{R}\Omega)$ where $\Omega$ is short for $(\phi,\theta,\chi)$. Then you need to identify the corresponding angle ranges, 
i.e. the sub-regions of angular integration that are mapped into each other by this same operation $\mathcal{R}$. It may be easiest to represent $\mathcal{R}$ as a rotation about an axis by a specified angle, but this is easily converted into other forms. If you are doing this numerically, it would make sense to check (numerically) once or twice that you get the same answer from the full integration as from the symmetry-reduced one.
However, for a symmetric crystal, surely the angular dependence of a fourth-rank tensor is known exactly?
For example, for cubic crystals, this open source paper by KM Knowles and PR Howie, J Elast, 120, 87 (2015) describes in eqn (5) how the stiffness tensor $C_{ijk\ell}$ (written in its full form, not in the reduced Voigt $6\times6$ notation) depends on orientation, which is represented as a $3\times3$ rotation matrix with components $a_{ij}$. (These can be written, if you like, in terms of the Euler angles).
Once you know the elements of $C_{ijk\ell}$ in a convenient set of axes based on the crystal unit cell,
computing the orientational average of these rotation matrix elements is relatively straightforward, and for an isotropic orientational average (which seems to be what you want, i.e. you haven't mentioned an orientational distribution function) it can surely be evaluated analytically.
An example is given in section 7 of that paper.
Of course, your crystal may not be cubic. But for any crystal symmetry, it should be possible to apply the same principles, and it may have already been done for your crystal, if you search in the literature.
