I am following the conventions of http://isites.harvard.edu/fs/docs/icb.topic1146665.files/III-9-RenormalizationGroup.pdf. Consider the QED Lagrangian

$$\mathcal{L}=-\frac{1}{4}Z_3F_{\mu\nu}^2+Z_2\bar{\psi}i\gamma^{\mu}\partial_{\mu}\psi-Z_2Z_mm\bar{\psi}\psi+Z_eZ_2\sqrt{Z_3}e\bar{\psi}\gamma^{\mu}A_{\mu}\psi+\sum^jC_j\mathcal{O}_j$$

where $\mathcal{O}_j=Z_j\partial^n\gamma^mA_{\mu}\ldots{}A_{\nu}\bar{\psi}\ldots\psi$ are operators with all fields evaluated at the same point that have any number of photons, fermions, gamma matrices, factors of the metric... and analytic dependence on derivatives. Everything is written using renormalized fields. Consider in particular the operator

$$\mathcal{O}=Z\bar{\psi}\gamma^{\mu}\partial_{\mu}\psi\bar{\psi}\gamma^{\nu}\partial_{\nu}\psi$$

I want to get the anomalous dimension of this operator at one loop. I know that in order to do that I need to get $Z$ but I am clueless of how to proceed.* *I have the feeling that I have to consider a correlation function but I don't know which. Any indication wouldbe greatly appreciated.