  We compute the covariance of the galaxy power spectrum multipoles in
perturbation theory, including the effects of nonlinear evolution, nonlinear
and nonlocal bias, radial redshift-space distortions, arbitrary survey window
and shot noise. We rewrite the power spectrum FKP estimator in terms of the
usual windowed galaxy fluctuations and the fluctuations in the number of
galaxies inside the survey volume. We show that this leads to a stronger
super-sample covariance than assumed in the literature and causes a substantial
leakage of Gaussian information. We decompose the covariance matrix into
several contributions that provide an insight into its behavior for different
biased tracers. We show that for realistic surveys, the covariance of power
spectrum multipoles is already dominated by shot noise and super survey mode
coupling in the weakly non-linear regime. Both these effects can be accurately
modeled analytically, making a perturbative treatment of the covariance very
compelling. Our method allows for the covariance to be varied as a function of
cosmology and bias parameters very efficiently, with survey geometry entering
as fixed kernels that can be computed separately using fast fourier transforms
(FFTs). We find excellent agreement between our analytic covariance and that
estimated from BOSS DR12 Patchy mock catalogs in the whole range we tested, up
to $k=0.6$ h/Mpc. This bodes well for application to future surveys such as
DESI and Euclid. The CovaPT code that accompanies this paper is available at
https://github.com/JayWadekar/CovaPT
