There exist a rectangular box and a cylinder in a 3D space. Their sizes, locations and orientations are arbitrary. They may or may not intersect with each other. If they do intersect with each other, theoretically, the volume of their intersection part can be expressed as a 3D integration over a complicated 3D region. Generally speaking, it is very difficult to evaluate this integration analytically. Therefore, I am require to use Monte Carlo method to estimate this integration and obtain an approximate result of the volume of the intersection part.
A rectangular box and a cylinder in 3D space

A rectangular box is described by its eight vertices P1, P2, …, P8. A cylinder is described by the two center points C1 and C2 of its top and bottom faces and the circle radius R of its bottom face.

Parameters given below:
P1: (0.000000e+00, 0.000000e+00, 0.000000e+00)
P2: (1.326481e-02, 8.673519e-01, 4.975186e-01)
P3: (1.011938e+00, 8.806167e-01, 4.477667e-01)
P4: (9.986735e-01, 1.326481e-02, -4.975186e-02)
P5: (4.975186e-02, -4.975186e-01, 8.660254e-01)
P6: (6.301667e-02, 3.698333e-01, 1.363544e+00)
P7: (1.061690e+00, 3.830981e-01, 1.313792e+00)
P8: (1.048425e+00, -4.842538e-01, 8.162735e-01)
C1: (5.557210e-01, -5.721025e-02, 1.089909e+00)
C2: (5.059692e-01, 4.403083e-01, 2.238834e-01)
R: 1.000000e+00
I would love to get some hints on how to come up mathematical expressions and matlab code to solve the question.
I assume that you know some vector notations. Given your points $\vec{P_i}$, the volume of the rectangular box is (https://en.wikipedia.org/wiki/Triple_product) $$V=[(\vec{P_4}-\vec{P_1})\times(\vec{P_2}-\vec{P_1})]\cdot (\vec{P_5}-\vec{P_1})$$ You generate points inside this volume in the following way: for every point you generate three random numbers between $0$ and $1$, say $p,q,r$. The $x,y,z$ coordinates of your point $\vec R_i$ are given by $$\vec R_i=(x,y,z)=(\vec{P_4}-\vec{P_1})p+(\vec{P_2}-\vec{P_1})q +(\vec{P_5}-\vec{P_1})r$$
You know that all $\vec R_i$ are inside the rectangular box, now you need to check if they are in the cylinder as well. For that, the first step is to calculate the angle between the axis of the cylinder and the position of the point with respect to one of the center points. I will use $\vec{C_1}$. The cosine of the angle is $$\cos\theta_1=\frac{(\vec{R_i}-\vec{C_1})\cdot(\vec{C_2}-\vec{C_1})}{\sqrt{(\vec{R_i}-\vec{C_1})\cdot(\vec{R_i}-\vec{C_1})}{\sqrt{(\vec{C_2}-\vec{C_1})\cdot(\vec{C_2}-\vec{C_1})}}}$$ Similarly you can calculate $\cos\theta_2$. From either of these, the distance from the $\vec R_i$ to the cylinder axis is for example $$d=|\vec{R_i}-\vec{C_1}||\sin\theta_1|$$ The conditions that the point is inside the cylinder are $d<R$, and $\theta_{1,2}<\pi/2$. The angle is in radians. The last conditions mean that the point $R_i$ is on the same side of the base $C_1$ as $C_2$, and on the same side of the base $C_2$ as $C_1$ (between bases).
The last step is to repeat the procedure for a large number of points, say $N$, and count how many points are in the cylinder, say $M$. The volume of the intersection is given by $$V_{int}=V\frac{M}{N}$$