First Fourier coefficient of weight $k$ holomorphic cusp form. Let $f$ a weight $k$ holomorphic Hecke cusp form with $\|f\|^2=\langle f,f\rangle=1$ with fourier expansion
$$f(z)=\sum_{r\geq 1}a_f(r)e(rz)$$
Let $\displaystyle\lambda_f(r)=\frac{a_f(r)r^{(-k+1)/2}}{a_f(1)}$, then $\lambda_f(r)$ are the eigenvalues of the normalized Hecke operators $n^{(1-k)/2}T_n$.
We define  $$L(s,\operatorname{sym}^2(f))=\zeta(2s)\sum_{n=1}^\infty \frac{\lambda_f(n^2)}{n^s}$$
I am trying to prove the following (this appears in different papers like this one in  Luo-Sarnak)

$$|a_f(1)|^2=\frac{(4\pi)^{k-1}}{\Gamma(k)}\frac{2\pi^2}{L(1,\operatorname{sym}^2(f))}.$$

What I was able to prove so far is that 
$$a_f(1)=\frac{\zeta(2s)}{L(s,\operatorname{sym}^2 (f))}\frac{L(k+s-1,\operatorname{sym}^2 (f))}{\zeta(2(k+s-1))}$$
If I put $s=1$, the problem reduces to find $L(1,\operatorname{sym}^2(f))$ and $L(k,\operatorname{sym}^2(f))$, but I don't think this is a good approach. Any ideas or references would be appreciated.

Thanks to @Peter Humphries in telling me the approach in his comment.
Denote $\Gamma=SL_2(\mathbb{Z})$ and consider $$F(s)=\int_{\Gamma\setminus\mathbb{H}}|f(z)|^2y^kE(z,s)d\mu(z)$$
where $E(z,s)=\sum_{\gamma\in\Gamma_\infty\setminus\Gamma}(\operatorname{Im}(\gamma z))^s$ is the usual Eisenstein series.
First, since $\operatorname{Res}_{s=1}(E(z,s))=\frac{3}{\pi}$, we have 
$$\operatorname{Res}_{s=1}F(s)=\frac{3}{\pi}\int_{\Gamma\setminus\mathbb{H}}|f(z)|^2y^kd\mu(z)=\frac{3}{\pi}\langle f,f\rangle=\frac{3}{\pi}$$
Now, let's compute this residue in a different manner: 
By modularity, it's easy to check  that $|f(z)|^2y^k=|f(z)|^2(\operatorname{Im}(z))^k$ is $\Gamma$-invariant, i.e.  $|f(z)|^2(\operatorname{Im}(z))^k=|f(\gamma z)|^2(\operatorname{Im}(\gamma z))^k$. So, we can unfold the integral obtaining 
$$F(s)=\int_{\Gamma_\infty\setminus\mathbb{H}}|f(z)|^2y^ky^sd\mu(z)=\int_{0}^\infty\int_{0}^1|f(z)|^2y^ky^sy^{-2}dxdy$$
Using the Fourier expansion of $f(z)$ and integrating  with respect to the $x$ variable, we have $$\int_{0}^1|f(z)|^2 dx=\sum_{r}|a_{f}(r)|^2e^{-4\pi ry}$$ 
So, $$F(s)=\sum_{r}|a_f(r)|^2\int_{0}^\infty y^{k+s-2}e^{-4\pi ry}dy$$
After a change of variable $t=4\pi ry$, we can see that $\int_{0}^\infty y^{k+s-2}e^{-4\pi ry}dy=\frac{\Gamma(k+s-1)}{(4\pi r)^{k+s-1}}$.
Then, \begin{align*}F(s)&=\sum_{r}|a_f(r)|^2\frac{\Gamma(k+s-1)}{(4\pi r)^{k+s-1}}\\
&=|a_f(1)|^2\frac{\Gamma(k+s-1)}{(4\pi )^{k+s-1}}\sum_{r}\frac{\lambda_f(r)^2}{r^s}\\
&=|a_f(1)|^2\frac{\Gamma(k+s-1)}{(4\pi )^{k+s-1}}\frac{\zeta(s)}{\zeta(2s)}L(s,\operatorname{sym}^2f)\end{align*}
Then $$\operatorname{Res}_{s=1}F(s)=|a_f(1)|^2\frac{\Gamma(k)}{(4\pi )^{k}}\frac{1}{\zeta(2)}L(1,\operatorname{sym}^2f)$$
Equating $$\frac{3}{\pi}=|a_f(1)|^2\frac{\Gamma(k)}{(4\pi )^{k}}\frac{1}{\zeta(2)}L(1,\operatorname{sym}^2f),$$
we get our result.
