This paper presents fully automated hp-mesh refinement strategies applied to diffusion equations. In hp strategies, both the mesh size and the polynomial order can vary locally. Numerical results show that exponential convergence rates are achieved for a fraction of the number of unknowns needed with uniform refinement and h-adaptive strategies. The treatment of adaptivity in the multigroup case and the derivation of goal-oriented estimators for neutronics calculations are described. The smoothness of the multigroup components can vary greatly as a function of the energy group; this fact called for the development of group-dependent adapted spatial meshes. The goal-oriented process combines the standard hp adaptation technique with a goal-oriented adaptivity based on the simultaneous solution of an adjoint problem in order to compute quantities of interest, such as reaction rates in subdomains and pointwise fluxes or currents. These algorithms are tested for various multigroup one-dimensional and two-dimensional diffusion problems, and the numerical results confirm the exponential convergence rates predicted theoretically.