This contribution presents a hierarchical multigrid approach for the solution of large-scale finite cell problems on both uniform grids and multi-level hp-discretizations. The proposed scheme takes advantage of the hierarchical basis functions utilized in the finite cell method and the multi-level hp-method, which is attributed to the use of high-order integrated Legendre basis functions and overlay meshes, to yield a simple and elegant multigrid scheme. This simplicity is reflected in the fact that transitioning between multigrid levels only involves the inclusion or exclusion of specific basis functions. All restriction and prolongation operators, therefore, reduce to binary matrices that do not need to be explicitly assembled or applied, saving computational time and memory. Elementwise and patchwise additive Schwarz smoothing techniques are used to mitigate the influence of the cut cells on the conditioning of the linear systems, while maintaining the parallelizability of the solver. The effectiveness of the scheme is numerically verified in various examples and convergence rates that are independent of the cut configuration, mesh size, refinement level and, in certain scenarios, even the polynomial order are shown. A series of numerical examples demonstrate the applicability of the scheme for solving large immersed systems with multiple millions and even billions of unknowns on massively parallel machines.