Energy-dependent self-adaptive mesh refinement algorithms are developed for a continuous Bubnov-Galerkin spatial discretization of the multi-group neutron diffusion equation using NURBS-based isogeometric analysis (IGA). The spatially self-adaptive algorithms employ both mesh (h) and polynomial degree (p) refinement. Constraint-based equations are established across irregular interfaces with hanging-nodes; they are based upon master-slave relationships and the conservative interpolation between surface meshes. A similar Galerkin projection is employed in the conservative interpolation between volume meshes to evaluate group-to-group source terms over energy-dependent meshes; and to evaluate interpolation-based error measures. Enforcing continuity over an irregular mesh does introduce discretization errors. However, local mesh refinement allows for a better allocation of computational resources; and thus, more accuracy per degree of freedom. Two a posteriori interpolation-based error measures are proposed. The first heuristically minimizes local contributions to the discretization error, which becomes competitive for global quantities of interest (QoIs). However, for localized QoIs, over energy-dependent meshes, certain multi-group components may become under-resolved. The second employs duality arguments to minimize important error contributions, which consistently and reliably reduces the error in the QoI.