In this work, we present numerical analysis for a distributed optimal control problem, with box constraint on the control, governed by a subdiffusion equation which involves a fractional derivative of order $alphain(0,1)$ in time. The fully discrete scheme is obtained by applying the conforming linear Galerkin finite element method in space, L1 scheme/backward Euler convolution quadrature in time, and the control variable by a variational type discretization. With a space mesh size $h$ and time stepsize $tau$, we establish the following order of convergence for the numerical solutions of the optimal control problem: $O(tau^{min({1}/{2}+alpha-epsilon,1)}+h^2)$ in the discrete $L^2(0,T;L^2(Omega))$ norm and $O(tau^{alpha-epsilon}+ell_h^2h^2)$ in the discrete $L^infty(0,T;L^2(Omega))$ norm, with any small $epsilon>0$ and $ell_h=ln(2+1/h)$. The analysis relies essentially on the maximal $L^p$-regularity and its discrete analogue for the subdiffusion problem. Numerical experiments are provided to support the theoretical results.