In this work we develop a Hessian-based sampling method for the construction of goal-oriented reduced order models with high-dimensional parameter inputs. Model reduction is known very challenging for high-dimensional parametric problems whose solutions also live in high-dimensional manifolds. However, the manifold of some quantity of interest (QoI) depending on the parametric solutions may be low-dimensional. We use the Hessian of the QoI with respect to the parameter to detect this low-dimensionality, and draw training samples by projecting the high-dimensional parameter to a low-dimensional subspace spanned by the eigenvectors of the Hessian corresponding to its dominating eigenvalues. Instead of forming the full Hessian, which is computationally intractable for a high-dimensional parameter, we employ a randomized algorithm to efficiently compute the dominating eigenpairs of the Hessian whose cost does not depend on the nominal dimension of the parameter but only on the intrinsic dimension of the QoI. We demonstrate that the Hessian-based sampling leads to much smaller errors of the reduced basis approximation for the QoI compared to a random sampling for a diffusion equation with random input obeying either uniform or Gaussian distributions.