Yes, there is a Gauss-Bonnet-Chern theorem for orbifolds, under a certain technical restriction. The proof is by reduction to the classical case, extending all definitions in the only meaningful way. It has been proven by Satake a long time ago:
http://projecteuclid.org/DPubS/Repository/1.0/Disseminate?view=body&id=pdf_1&handle=euclid.jmsj/1261153826
Here is an outline, in my own terms. For me, an orbifold $O$ is a smooth Deligne-Mumford stack (see here: http://arxiv.org/pdf/math/0203100). There is the coarse moduli space $|O|$. If $O$ is a quotient $M//G$ of some proper action of a discrete group, then $|O|$ is the space of orbits. In the definition of orbifolds that is used in $3$-dimensional topology, $|O|$ would be the underlying space which has some atlasses.
On orbifolds, we have the following structures available: there is the tangent bundle $TO \to O$; $TO$ is an orbifold as well; this is a vector bundle in the category of stacks; the quotient $|TO| \to |O|$ is not quite a vector bundle.
Let me make the following assumption:
There exists a compact oriented manifold $N$ and a finite covering $N \to O$ of degree $d$.
Under these circumstances, we can define the Euler number of the orbifold $\chi(O) \in \mathbb{Q}$ as follows:
$$\chi(O) := \frac{1}{d} \chi(N).$$
Using the multiplicativity of the ordinary Euler number, it is easy to see that $\chi(O)$ only depends on $O$ and not on $N$.
Now I define the geometric side of the G-B-C theorem.
Orbifolds have tangent bundles and we can talk about differential forms on orbifolds, de Rham complex, connnections on vector bundles and the Levi-Civitta connection in the same way as for manifolds. The cleanest way to define connections might be the framework of connections on principal bundles. The last piece of differential geometry needed is Chern-Weil theory. Conclusion: to any Riemann metric on an $2n$-dimensional oriented orbifold (i.e., the tangent bundle is oriented), we get an Euler-Form $e$, a closed $2n$-form.
This form represents an element in the cohomology of the orbifold, but we have to say what this means. If $O=M //G$ were a global quotient, then $H^* (O) := H^* (EG \times_G M)$, the Borel-equivariant cohomology. If $O$ is not a global quotient, you have to replace the Borel space by the homotopy type of the orbifold.
How do you integrate the $2n$-form? Well, given a closed oriented $2n$-manifold $N$ and a map $f:N \to O$ as before,
we define
$$\int_O e := \frac{1}{d} \int_N f^* e,$$
which is the only reasonable way to define integration over orbifolds. Since $N$ is a manifold, we can apply the classical Gauss-Bonnet-Chern theorem and find that
$$\chi(O)=\int_O e.$$
When $E\to M$ is an oriented vector bundle of rank $2n$ over a
compact manifold $M$, it has a well-defined de Rham Euler class $e(E)$
in $H^{2n}_{dR}(M)$, and a representative $2n$-form for $e(E)$
can be computed as follows:
Fix a positive definite inner product $\langle,\rangle$ on $E$. (Since any two such inner products are equivalent under automorphisms of $E$, it doesn't matter which one.) Let $\nabla$ be a $\langle,\rangle$-orthogonal connection on $E$, and let $K^\nabla$ denote the curvature of $\nabla$, regarded as a $2$-form with values in ${\frak{so}}(E)$. Let $e(\nabla) = \mathrm{Pf}(K^\nabla/2\pi)$, which is a well-defined $2n$-form on $M$. (N.B.: The definition of $\mathrm{Pf}$ requires both the inner product and the orientation of $E$.) Then $e(E)=\bigl[e(\nabla)\bigr]\in H^{2n}_{dR}(M)$. (That $\bigl[e(\nabla)\bigr]$ is independent of the choice of $\nabla$
is one of the first results proved in Chern-Weil theory.)
If $M$ is a compact, oriented $2n$-manifold, then the value of $e(E)$ on $[M]$, the fundamental class of $M$, can be computed as follows: Let $Y$ be a section of $E$ that has only a finite number of zeroes. (By Whitney transversality, the generic section of $E$ satisfies this condition.) Using the orientations of $E$ and $M$, one defines the index
of an isolated zero $z$ of $Y$, which is an integer $\iota_Y(z)$. Then
$$
e(E)\bigl([M]\bigr) = \int_M e(\nabla) = \sum_{z\in Z}\ \iota_Y(z).
$$
By the Poincaré-Hopf Theorem, when $E = TM$ (as oriented bundles),
this sum is equal to the Euler characteristic of $M$, which explains
why $e(E)$ is called the Euler class of $E$.
To prove this result, one chooses an orthogonal connection $\nabla$ on $E$ that depends on $Y$ and whose Euler form $e(\nabla)$ is easy to evaluate explicitly. This can be done as follows:
Let $Z\subset M$ be the (finite) zero set of $Y$ and let $U\subset M$ be an open neighborhood of $Z$ that consists of disjoint, smoothly embedded $2n$-balls, one around each element of $Z$. Let $\phi$ be a smooth function on $M$ that is identically equal to $1$ on an open neighborhood of $Z$ and whose support $K\subset U$ is a disjoint union of closed, smoothly embedded balls,
one around each element of $Z$.
Now, $E$ is trivial over $U$,
so choose a positively oriented $\langle,\rangle$-orthonormal basis of sections
$s_1,\ldots, s_{2n}$ of $E$ over $U$ and let $\nabla_1$ be the (flat) connection on $E_U$ for which these sections are parallel. Let $\bar Y = Y/|Y|$ be the normalized unit section defined on $M\setminus Z$, and define a second connection on $U$ by
$$
\nabla_2 s = \nabla_1 s
+ (1{-}\phi)\bigl( \bar Y\otimes \langle s,\nabla_1 \bar Y\ \rangle
- \langle s,\bar Y\ \rangle\ \nabla_1\bar Y\ \bigr).
$$
It is easily verified that this formula does define a connection on $U$,
that $\nabla_2$ is $\langle,\rangle$-orthogonal, and that,
outside of $K$ (i.e., where $\phi\equiv0$), the vector field $\bar Y$ is
$\nabla_2$-parallel. (Note that, because $\phi\equiv1$ near $Z$, where $Y$ is
not defined, this formula does extend smoothly across $Z$, agreeing with $\nabla_1$
on a neighborhood of $Z$.)
Meanwhile, on $M\setminus K$, write $E$ as a $\langle,\rangle$-orthogonal direct sum $E = \mathbb{R}\cdot Y \oplus E'$ and choose a $\langle,\rangle$-compatible connection $\nabla_3$ on $E$ over $M\setminus K$ that preserves this splitting. Using a partition of unity subordinate to the open cover of $M$ defined by $U$ and $M\setminus K$, construct a $\langle,\rangle$-orthogonal connection $\nabla$ that agrees with $\nabla_2$ on $K$ and with $\nabla_3$ on $M\setminus U$, and for which the line bundle $\mathbb{R}\cdot Y\subset E$
is parallel on $M\setminus K$.
Now, on $M\setminus K$, the curvature of $\nabla$ takes values in ${\frak{so}}(E')\subset{\frak{so}}(E)$, so $e(\nabla) = \textrm{Pf}\bigl(K^{\nabla}/(2\pi)\bigr)$ vanishes identically outside of $K$. It suffices now to evaluate the
integral of $e(\nabla)$ over a single component $B$ of $K$, which may be
assumed to be the unit ball in $\mathbb{R}^{2n}$, so restrict attention to
this case. Write $\bar Y = s_1 u_1 +\cdots + s_{2n}\ u_{2n}$, and note that,
on $B$, one has, by definition,
$$
\nabla s_i = \nabla_2 s_i
= \sum_{j=1}^{2n} s_j\otimes (1{-}\phi)(u_j\ du_i - u_i\ du_j),
$$
i.e., the connection $1$-forms of $\nabla$ in this basis are
$\omega_{ij} = (1{-}\phi)(u_i\ du_j - u_j\ du_i)$.
Using the identity ${u_1}^2 +\cdots + {u_{2n}}^2 = 1$,
the curvature forms are easily computed to satisfy
$$
\Omega_{ij} = d\omega_{ij} + \omega_{ik}\wedge\omega_{kj}
= (1{-}\phi^2)\ du_i\wedge du_j - d\phi\wedge(u_i\ du_j - u_j\ du_i).
$$
At this point, you have to know the definition of the Pfaffian. (I'll wait while you look it up.) Using the fact that the result has to be invariant under $SO(2n)$-rotations, it is easy to show that, on $B$, one has
$$
e(\nabla) = \textrm{Pf}\left(\frac{\Omega}{2\pi}\right)
= c_n (1{-}\phi^2)^{n-1} d\phi\wedge u^*(\Upsilon)
$$
for some universal constant $c_n$ and where $\Upsilon$ is the $SO(2n)$-invariant $2n$-form
of unit volume on $S^{2n-1}\subset\mathbb{R}^{2n}$ and $u: B\setminus z\to S^{2n-1}$
is $u = (u_1,\ldots, u_{2n})$. (I'll let you evaluate the constant $c_n$. This can be done in a number of ways, but it's essentially a combinatorial exercise.) In particular,
$$
e(\nabla) = d\bigl(P_n(1{-}\phi)\ u^*(\Upsilon)\bigr)
$$
where $P_n(t)$ is the polynomial in $t$ (of degree $2n{-}1)$) that vanishes at $t=0$ and
satisfies $P'_n(1{-}t) = -c_n(1{-}t^2)^{n-1}$. By Stokes' Theorem, one has
$$
\int_B e(\nabla) = P_n(1) \int_{\partial B}u^*(\Upsilon)
= P_n(1)\ \textrm{deg}(u:\partial B\to S^{2n-1}) = P_n(1)\ \iota_Y(z)
$$
Thus, it follows that there is a universal constant $C_n = P_n(1)$ such that
$$
e(E)\bigl([M]\bigr) = \int_M e(\nabla) = C_n \ \sum_{z\in Z}\ \iota_Y(z).
$$
Now, $C_1 = 1$, as one can show by computing with the constant curvature metric on the $2$-sphere and noting that, for $Y$ the gradient of the height function on $S^2$, the sum of the indices is $2$. Finally, by properties of the Pfaffian and the index, one sees that,
if such a formula as above is to hold, then one must have $C_{n+m}=C_nC_m$. Thus, $C_n=1$ for all $n$, and the formula is proved.
(Of course, one can avoid these final tricks for evaluating the constants by carefully doing the combinatorial exercise, but I'll leave that to the curious.)
Best Answer
One can find some applications in Arthur Besse's book "Einstein manifolds", $\S6\ D$.
Proposition. If a closed manifold $M^4$ admits an Einstein metric, then $\chi(M)\ge 0$. Moreover $\chi(M)=0$ iff $M$ admits a flat metric.
Proof. Since $Rc(g)=\lambda g$, we have $|Rc|^2=4\lambda^2=r^2/4$. This implies $\chi(M)=1/32\pi^2 \int|Rm|^2dV\ge 0$ with equality iff $Rm=0$. $\square$
Combining with another formula for the signature of a 4-manifold $$ \sigma(M)=\frac{1}{12\pi^2}\int_M(|W_+|^2-|W_-|^2)dV, $$ where $W_{\pm}$ are compontens of the Weyl tensor, one can prove a stronger inequality: for a closed Einstein manifold $M^4$ $$ \chi(M)\ge 3/2|\sigma(M)|. $$
The equality case was studied by Nigel Hitchin, Compact four-dimensional Einstein manifolds, (1974). He proved that if $(M^4, g)$ is an Einstein manifold with $\chi(M)= 3/2|\sigma(M)|$, then $M$ is Ricci flat and either universal cover of $M$ is K3-surface, or $M$ is flat.
These are important results, since topological obstructions for the existence of Einstein metrics are very scarce. In particular in dimensions $\ge 5$ there are no known such obstructions.