The reaction pi- d -> n n gamma is calculated in chiral perturbation theory so as to facilitate an extraction of the neutron-neutron scattering length (a_nn). We include all diagrams up to O(Q^3). This includes loop effects in the elementary pi- p -> gamma n amplitude and two-body diagrams, both of which were ignored in previous calculations. We find that the chiral expansion for the ratio of the quasi-free (QF) to final-state-interaction (FSI) peaks in the final-state neutron spectrum converges well. Our third-order calculation of the full spectrum is already accurate to better than 5%. Extracting a_nn from the shape of the entire pi- d -> n n gamma spectrum using our calculation in its present stage would thus be possible at the +-0.8 fm level. A fit to the FSI peak only would allow an extraction of a_nn with a theoretical uncertainty of +-0.2 fm. The effects that contribute to these error bars are investigated. The uncertainty in the $nn$ rescattering wave function dominates. This suggests that the quoted theoretical error of +-0.3 fm for the most recent pi- d -> n n gamma measurement may be optimistic. The possibility of constraining the nn rescattering wave function used in our calculation more tightly--and thus reducing the error--is briefly discussed.