The classical Selberg zeta function is a mediator between spectral entities and dynamical entities of hyperbolic surfaces, as it is defined by means of the geodesic length spectrum and encodes in its zeros the spectral parameters of the Laplacian of the considered hyperbolic surface. We will consider the Selberg zeta function of infinite-area, geometrically finite hyperbolic orbisurfaces with twists by finite-dimensional unitary representations and hence for vector-valued situations. We will present a factorization formula in terms of the Weierstrass product of the Laplace resonances, Barnes G-functions, gamma functions and the singularity degrees of the representation. Similar to the classical, untwisted case, this provides a spectral and geometric interpretation of the zeros and poles of the Selberg zeta function, but this time by spectral and geometric entities of the orbisurface and by the representation. We will see that this factorization formula generalizes the factorization result by Borthwick, Judge and Perry to hyperbolic orbisurfaces with orbifold singularities as well as to unitary twists. We will further see that the presence of orbifold singularities yields a separate, previously unobserved contribution to the factorization formula, even in the untwisted case. This is joint work with Moritz Doll.