diff --git a/bin/picca_xcf.py b/bin/picca_xcf.py
index 63c693356..c96b2d266 100755
--- a/bin/picca_xcf.py
+++ b/bin/picca_xcf.py
@@ -15,7 +15,7 @@
 
 from picca import constants, xcf, io, prep_del, utils
 from picca.data import Forest
-from picca.utils import userprint
+from picca.utils import userprint, get_qso_weights
 
 
 def corr_func(healpixs):
@@ -389,6 +389,10 @@ def main(cmdargs):
     z[w] /= weights_list.sum(axis=0)[w]
     num_pairs = num_pairs_list.sum(axis=0)
 
+    # Get stacked QSO weights for metal matrix computations
+    z_qso_cat = np.array([qso.z for qso in healpix for healpix in xcf.objs.values()])
+    z_qso, weights_qso = get_qso_weights(z_qso_cat, args.z_ref, args.z_evol_obj)
+
     results = fitsio.FITS(args.out, 'rw', clobber=True)
     header = [{
         'name': 'RPMIN',
@@ -466,6 +470,13 @@ def main(cmdargs):
                   header=header2,
                   extname='COR')
 
+    header3 = [{}]
+    results.write([z_qso, weights_qso],
+                  names=['Z_BIN', 'QSO_WEIGHTS'],
+                  comment=['Redshift', 'QSO weights'],
+                  header=header3,
+                  extname='QSO_WEIGHTS')
+
     results.close()
 
     t3 = time.time()
diff --git a/py/picca/utils.py b/py/picca/utils.py
index 95ce0579d..7ad5fdf3f 100644
--- a/py/picca/utils.py
+++ b/py/picca/utils.py
@@ -491,3 +491,33 @@ def unred(wave, ebv, R_V=3.1, LMC2=False, AVGLMC=False):
     corr = 1. / (10.**(0.4 * curve))
 
     return corr
+
+
+def get_qso_weights(z_qso_cat, z_ref, z_evol, zbins=1000):
+    """Computes the stacked QSO weights for the fast metal matrix computation
+
+    Parameters
+    ----------
+    z_qso_cat : array
+        Array of quasar redshifts read from the QSO catalog
+    z_ref : float
+        Reference redshift
+    z_evol : float
+        Redshift evolution parameter
+    zbins : int, optional
+        Number of bins in the histogram, by default 1000
+
+    Returns
+    -------
+    (array, array)
+        QSO redshift bins and their corresponding weights
+    """
+    weights_qso_cat = ((1. + z_qso_cat) / (1. + z_ref))**(z_evol - 1.)
+
+    histo_w, zbins = np.histogram(z_qso_cat, bins=zbins, weights=weights_qso_cat)
+    histo_wz, _ = np.histogram(z_qso_cat, bins=zbins, weights=weights_qso_cat*z_qso_cat)
+    selection = histo_w > 0
+    z_qso = histo_wz[selection] / histo_w[selection]  # weighted mean in bins
+    weights_qso = histo_w[selection]
+
+    return z_qso, weights_qso