We present a contour integral method (CIM) for the efficient computation of outputs from parametric linear systems in control form over specified time intervals and under a user-prescribed accuracy. The CIM framework employs a quadrature rule to approximate the inverse Laplace transform, executed along a modified integration contour. For the parametric case, we demonstrate how this method integrates effectively with projection-based model order reduction (MOR) where the projection spaces are constructed via a greedy-type algorithm which is guided by an appropriately derived error estimate and adheres to Hermite interpolation conditions. This methodology significantly reduces computational expenses when evaluating the input-output relations for varied parameters across the parametric domain, as well as for a wide class of input functions and initial conditions sufficiently captured by a low-dimensional subspace. We validate the precision and efficacy of the approach using numerical experiments on a range of benchmark test problems for both non-parametric and parametric linear control systems.