You are familiar with $\int_a^b f(x)\,dx$, which is an integral over an interval. One can also define "double" integrals, that is, integrals over a region of the plane. And triple integrals, and so on.
Here they are finding the double integral of a certain product over the unit square. The second term $P(x,y)$ is explained in your post. The first term is $I((x^2+y^2)\lt 1))$. This is an indicator function, which by definition is $1$ where $x^2+y^2\lt 1$, and $0$ elsewhere.
So the product that we are integrating is $1$ inside the quarter-circle, and $0$ outside. The integral over the square of the product is therefore the area of the quarter-circle.
I expect you understand the rest. But for completeness, effectively the program throws a dart repeatedly at the $1\times 1$ square, and records a $1$ for every time the dart lands in the quarter circle. If the number $N$ of dart throws is very large, and $H$ is the number of hits, one would expect $\dfrac{H}{N}$ to be close to the ratio of the area of the quarter-circle to the area of the square, that is, to $\dfrac{\pi}{4}$.
Alternately, one could use older notation, get rid of the indicator function, and integrate over the quarter-circle. But the indicator function version is closer in spirit to the dart throwing model, in the sense that, like in the dart throwing, we record a $1$ whenever we land inside the quarter-circle.
Remark: Although this is an inefficient way of approximating $\pi$, the general idea is extremely useful. It can be used to approximate the volume of extremely complicated multidimensional regions. Tou might look into Monte Carlo Methods.