Repository: igmhub/picca (language: python)

Issue:
Port Pk1d changes to current picca versions
We'll need to port changes in the Pk1d analysis for resolution matrix treatment to current versions of picca.


Discussion / hints:

Fix patch:
diff --git a/bin/picca_Pk1D.py b/bin/picca_Pk1D.py
index e4048b49c..40c602878 100755
--- a/bin/picca_Pk1D.py
+++ b/bin/picca_Pk1D.py
@@ -7,152 +7,371 @@
 from array import array
 import numpy as np
 import fitsio
+import os
 
 from picca import constants
 from picca.data import Delta
-from picca.pk1d import (compute_correction_reso, compute_pk_noise,
+from picca.pk1d import (compute_correction_reso,
+                        compute_correction_reso_matrix, compute_pk_noise,
                         compute_pk_raw, fill_masked_pixels, rebin_diff_noise,
                         split_forest)
 from picca.utils import userprint
+from multiprocessing import Pool
 
 
-def make_tree(tree, max_num_bins):
-    """Makes the ROOT tree to save the data
+def check_linear_binning(delta):
+    """checks if the wavelength binning is linear or log this is stable against masking
 
     Args:
-        tree: TTree
-            The ROOT tree to fill
-        max_num_bins: int
-            Maximum number of bins allowed
+        delta (Delta): delta class as read with Delta.from_...
+
+    Raises:
+        ValueError: Raised if binning is neither linear nor log, or if delta.log_lambda was actually wavelength
 
     Returns:
-        The following variables:
-            z_qso: Quasar redshifts
-            mean_z: Mean redshift of the forest
-            mean_reso: Mean resolution of the forest
-            mean_snr: Mean signal-to-noise ratio in the forest
-            lambda_min_tree: Minimum wavelength (in Angs)
-            lambda_max_tree: Maximum wavelength (in Angs)
-            plate: Plate number of the observation
-            mjd: Modified Julian Date of the observation
-            fiber: Fiberid of the object
-            nb_mask_pix: Number of masked pixels
-            num_bins_tree: Number of bins (storage tree version)
-            k_tree: Fourier modes (tree version)
-            pk_tree: Power spectrum for the different fourier modes (storage
-                tree version)
-            pk_raw_tree: Raw power spectrum for the different fourier modes
-                (storage tree version)
-            pk_noise_tree: Noise power spectrum for the different fourier modes
-                (storage tree version)
-            correction_reso_tree: Resolution for the correlation function
-                (storage tree version)
-            pk_diff_tree: Power spectrum of exposures_diff for the different
-                fourier modes (storage tree version)
+        linear_binning (bool): boolean telling the binning_type
+        pixel_step (float): size of a wavelength bin in the right unit
     """
-    z_qso = array('f', [0.])
-    mean_z = array('f', [0.])
-    mean_reso = array('f', [0.])
-    mean_snr = array('f', [0.])
-    num_masked_pixels_tree = array('f', [0.])
-
-    lambda_min_tree = array('f', [0.])
-    lambda_max_tree = array('f', [0.])
-
-    plate = array('i', [0])
-    mjd = array('i', [0])
-    fiber = array('i', [0])
-
-    num_bins_tree = array('i', [0])
-    k_tree = array('f', max_num_bins * [0.])
-    pk_tree = array('f', max_num_bins * [0.])
-    pk_raw_tree = array('f', max_num_bins * [0.])
-    pk_noise_tree = array('f', max_num_bins * [0.])
-    pk_diff_tree = array('f', max_num_bins * [0.])
-    correction_reso_tree = array('f', max_num_bins * [0.])
-
-    tree.Branch("z_qso", z_qso, "z_qso/F")
-    tree.Branch("mean_z", mean_z, "mean_z/F")
-    tree.Branch("mean_reso", mean_reso, "mean_reso/F")
-    tree.Branch("mean_snr", mean_snr, "mean_snr/F")
-    tree.Branch("lambda_min", lambda_min_tree, "lambda_min/F")
-    tree.Branch("lambda_max", lambda_max_tree, "lambda_max/F")
-    tree.Branch("nb_masked_pixel", num_masked_pixels_tree, "nb_mask_pixel/F")
-
-    tree.Branch("plate", plate, "plate/I")
-    tree.Branch("mjd", mjd, "mjd/I")
-    tree.Branch("fiber", fiber, "fiber/I")
-
-    tree.Branch("NbBin", num_bins_tree, "NbBin/I")
-    tree.Branch("k", k_tree, "k[NbBin]/F")
-    tree.Branch("Pk_raw", pk_raw_tree, "Pk_raw[NbBin]/F")
-    tree.Branch("Pk_noise", pk_noise_tree, "Pk_noise[NbBin]/F")
-    tree.Branch("Pk_diff", pk_diff_tree, "Pk_diff[NbBin]/F")
-    tree.Branch("cor_reso", correction_reso_tree, "cor_reso[NbBin]/F")
-    tree.Branch("Pk", pk_tree, "Pk[NbBin]/F")
-
-    return (z_qso, mean_z, mean_reso, mean_snr, lambda_min_tree,
-            lambda_max_tree, plate, mjd, fiber, num_masked_pixels_tree,
-            num_bins_tree, k_tree, pk_tree, pk_raw_tree, pk_noise_tree,
-            correction_reso_tree, pk_diff_tree)
-
-
-def compute_mean_delta(log_lambda, delta, ivar, z_qso, hist_delta,
-                       hist_delta_rest_frame, hist_delta_obs_frame, hist_ivar,
-                       hist_snr, hist_weighted_delta_rest_frame,
-                       hist_weighted_delta_obs_frame):
-    """Computes the mean delta and stores it in the control histogram variables
 
-    Args:
-        log_lambda: array of floats
-            Logarithm of the wavelength (in Angs)
-        delta: array of floats
-            Mean transmission fluctuation (delta field)
-        ivar: array of floats
-            Inverse variance
-        z_qso: float
-            Redshift of the quasar
-        hist_delta: TProfile2D
-            Root Profile2D histogram to store the mean delta as a function of
-            lambda-lambda_rest_frame
-        hist_delta_rest_frame: TProfile
-            Root Profile histogram to store the mean delta as a function of
-            lambda_rest_frame
-        hist_delta_obs_frame: TProfile
-            Root TProfile histogram to store the mean delta as a function of
-            observed lambda
-        hist_ivar: TH1D
-            Root TH1D histogram to store the inverse variance
-        hist_snr: TH1D
-            Root TH1D histogram to store the signal-to-noise ratio
-        hist_weighted_delta_rest_frame: TProfile
-            Root TProfile histogram to store the mean weighted delta as a
-            function of lambda_rest_frame
-        hist_weighted_delta_obs_frame: TProfile
-            Root TProfile histogram to store the mean weighted delta as a
-            function of observed lambda
-    """
+    diff_lambda = np.diff(10**delta.log_lambda)
+    diff_log_lambda = np.diff(delta.log_lambda)
+    q5_lambda, q25_lambda = np.percentile(diff_lambda, [5, 25])
+    q5_log_lambda, q25_log_lambda = np.percentile(diff_log_lambda, [5, 25])
+    if (q25_lambda - q5_lambda) < 1e-6:
+        #we can assume linear binning for this case
+        linear_binning = True
+        pixel_step = np.min(diff_lambda)
+    elif (q25_log_lambda - q5_log_lambda) < 1e-6 and q5_log_lambda < 0.01:
+        #we can assume log_linear binning for this case
+        linear_binning = False
+        pixel_step = np.min(diff_log_lambda)
+    elif (q5_log_lambda >= 0.01):
+        raise ValueError(
+            "Could not figure out if linear or log wavelength binning was used, probably submitted lambda as log_lambda"
+        )
+    else:
+        raise ValueError(
+            "Could not figure out if linear or log wavelength binning was used"
+        )
+
+    return linear_binning, pixel_step
+
+
+# loop over input files
+num_data=0
+def process_all_files(index_file_args):
+    global num_data
+    file_index, file, args = index_file_args
+    if file_index % 5 == 0:
+        userprint("\rread {} of {} {}".format(file_index, args.len_files,
+                                              num_data),
+                  end="")
+
+    # read fits or ascii file
+    if args.in_format == 'fits':
+        hdul = fitsio.FITS(file)
+        deltas = [Delta.from_fitsio(hdu, pk1d_type=True) for hdu in hdul[1:]]
+    elif args.in_format == 'ascii':
+        ascii_file = open(file, 'r')
+        deltas = [Delta.from_ascii(line) for line in ascii_file]
+
+    #add the check for linear binning on first spectrum only (assuming homogeneity within the file)
+    delta = deltas[0]
+    linear_binning, pixel_step = check_linear_binning(delta)
+    if linear_binning:
+        userprint("\n\nUsing linear binning, results will have units of AA")
+        if (args.disable_reso_matrix or not hasattr(delta, 'resolution_matrix')
+                or delta.resolution_matrix is None):
+            userprint(
+                "Resolution matrix not found or disabled, using Gaussian resolution correction\n"
+            )
+            reso_correction = "Gaussian"
+        else:
+            userprint("Using Resolution matrix for resolution correction\n")
+            reso_correction = "matrix"
+    else:
+        userprint("\n\nUsing log binning, results will have units of km/s")
+        reso_correction = "Gaussian"
+        userprint("Using Gaussian resolution correction\n")
+
+    use_exp_diff_cut=False
+    #add check if diff has ever been set
+    for delta in deltas:
+        #use this cut if some spectra have exposures_differences calculated
+        if not np.sum(delta.exposures_diff)==0:
+            use_exp_diff_cut = True
+            break
+
+    num_data += len(deltas)
+    userprint("\n ndata =  ", num_data)
+    results = None
+
+    for delta in deltas:
+        if (delta.mean_snr <= args.SNR_min
+                or delta.mean_reso >= args.reso_max):
+            continue
+        if use_exp_diff_cut:
+            if np.sum(delta.exposures_diff)==0:
+                continue
+
+        # first pixel in forest
+        selected_pixels = 10**delta.log_lambda > args.lambda_obs_min
+        #this works as selected_pixels returns a bool and argmax points
+        #towards the first occurance for equal values
+        first_pixel_index = (np.argmax(selected_pixels) if
+                             np.any(selected_pixels) else len(selected_pixels))
+
+        # minimum number of pixel in forest
+        min_num_pixels = args.nb_pixel_min
+        if (len(delta.log_lambda) - first_pixel_index) < min_num_pixels:
+            continue
+
+        # Split the forest in n parts
+        max_num_parts = (len(delta.log_lambda) -
+                         first_pixel_index) // min_num_pixels
+        num_parts = min(args.nb_part, max_num_parts)
+
+        #the split_forest function works with either binning, but needs to be uniform
+        if linear_binning:
+            split_array = split_forest(
+                num_parts,
+                pixel_step,
+                10**delta.log_lambda,
+                delta.delta,
+                delta.exposures_diff,
+                delta.ivar,
+                first_pixel_index,
+                reso_matrix=(delta.resolution_matrix
+                             if reso_correction == 'matrix' else None),
+                linear_binning=True)
+            if reso_correction == 'matrix':
+                (mean_z_array, lambda_array, delta_array, exposures_diff_array,
+                 ivar_array, reso_matrix_array) = split_array
+            else:
+                (mean_z_array, lambda_array, delta_array, exposures_diff_array,
+                 ivar_array) = split_array
+        else:
+            (mean_z_array, log_lambda_array, delta_array, exposures_diff_array,
+             ivar_array) = split_forest(num_parts, pixel_step,
+                                        delta.log_lambda, delta.delta,
+                                        delta.exposures_diff, delta.ivar,
+                                        first_pixel_index)
+
+        #the rebin_diff_noise function works with either binning, but needs to be uniform
+        for part_index in range(num_parts):
+            # rebin exposures_diff spectrum
+            if (args.noise_estimate == 'rebin_diff'
+                    or args.noise_estimate == 'mean_rebin_diff'):
+                if linear_binning:
+                    exposures_diff_array[part_index] = rebin_diff_noise(
+                        pixel_step, lambda_array[part_index],
+                        exposures_diff_array[part_index])
+                else:
+                    exposures_diff_array[part_index] = rebin_diff_noise(
+                        pixel_step, log_lambda_array[part_index],
+                        exposures_diff_array[part_index])
+
+            # Fill masked pixels with 0.
+            #the fill_masked_pixels function works with either binning, but needs to be uniform
+            if linear_binning:
+                #the resolution matrix does not need to have pixels filled in any way...
+                (lambda_new, delta_new, exposures_diff_new, ivar_new,
+                 num_masked_pixels) = fill_masked_pixels(
+                     pixel_step, lambda_array[part_index],
+                     delta_array[part_index], exposures_diff_array[part_index],
+                     ivar_array[part_index], args.no_apply_filling)
+            else:
+                (log_lambda_new, delta_new, exposures_diff_new, ivar_new,
+                 num_masked_pixels) = fill_masked_pixels(
+                     pixel_step, log_lambda_array[part_index],
+                     delta_array[part_index], exposures_diff_array[part_index],
+                     ivar_array[part_index], args.no_apply_filling)
+            if num_masked_pixels > args.nb_pixel_masked_max:
+                continue
 
-    for index, _ in enumerate(log_lambda):
-        lambda_ = np.power(10., log_lambda[index])
-        lambda_rf = lambda_ / (1. + z_qso)
-        hist_delta.Fill(lambda_, lambda_rf, delta[index])
-        hist_delta_rest_frame.Fill(lambda_rf, delta[index])
-        hist_delta_obs_frame.Fill(lambda_, delta[index])
-        hist_ivar.Fill(ivar[index])
-        snr_pixel = (delta[index] + 1) * np.sqrt(ivar[index])
-        hist_snr.Fill(snr_pixel)
-        hist_ivar.Fill(ivar[index])
-        if ivar[index] < 1000:
-            hist_weighted_delta_rest_frame.Fill(lambda_rf, delta[index],
-                                                ivar[index])
-            hist_weighted_delta_obs_frame.Fill(lambda_, delta[index],
-                                               ivar[index])
+            # Compute pk_raw, needs uniform binning
+            if linear_binning:
+                k, pk_raw = compute_pk_raw(pixel_step,
+                                           delta_new,
+                                           linear_binning=True)
+            else:
+                k, pk_raw = compute_pk_raw(pixel_step,
+                                           delta_new,
+                                           linear_binning=False)
+
+            # Compute pk_noise
+            run_noise = False
+            if args.noise_estimate == 'pipeline':
+                run_noise = True
+            if linear_binning:
+                pk_noise, pk_diff = compute_pk_noise(pixel_step,
+                                                     ivar_new,
+                                                     exposures_diff_new,
+                                                     run_noise,
+                                                     linear_binning=True,
+                                                     num_noise_exposures=args.num_noise_exp)
+            else:
+                pk_noise, pk_diff = compute_pk_noise(pixel_step,
+                                                     ivar_new,
+                                                     exposures_diff_new,
+                                                     run_noise,
+                                                     linear_binning=False,
+                                                     num_noise_exposures=args.num_noise_exp)
+
+            # Compute resolution correction, needs uniform binning
+            if linear_binning:
+                #in this case all is in AA space
+                if reso_correction == 'matrix':
+                    correction_reso = compute_correction_reso_matrix(
+                        reso_matrix=np.mean(reso_matrix_array[part_index],
+                                            axis=1),
+                        k=k,
+                        delta_pixel=pixel_step,
+                        num_pixel=len(lambda_new))
+                elif reso_correction == 'Gaussian':
+                    #this is roughly converting the mean resolution estimate back to pixels
+                    #and then multiplying with pixel size
+                    mean_reso_AA = pixel_step * delta.mean_reso_pix
+                    correction_reso = compute_correction_reso(
+                        delta_pixel=pixel_step, mean_reso=mean_reso_AA, k=k)
+            else:
+                #in this case all is in velocity space
+                delta_pixel = (pixel_step * np.log(10.) *
+                               constants.speed_light / 1000.)
+                correction_reso = compute_correction_reso(
+                    delta_pixel=delta_pixel, mean_reso=delta.mean_reso, k=k)
+
+            # Compute 1D Pk
+            if args.noise_estimate == 'pipeline':
+                pk = (pk_raw - pk_noise) / correction_reso
+            elif args.noise_estimate == 'mean_pipeline':
+                if args.kmin_noise_avg is None and linear_binning:
+                    #this is roughly the same range as eBOSS analyses for z=2.2
+                    selection = (k > 0) & (k < 1.5)
+                elif args.kmin_noise_avg is None:
+                    selection = (k > 0) & (k < 0.02)
+                else:
+                    selection = (((k > args.kmin_noise_avg) if args.kmax_noise_avg is not None else 1) & 
+                                 ((k < args.kmax_noise_avg) if args.kmax_noise_avg is not None else 1))
+                mean_pk_noise = np.mean(pk_noise[selection])
+                pk = (pk_raw - pk_noise) / correction_reso
+
+
+                
+            elif (args.noise_estimate == 'diff' or args.noise_estimate == 'rebin_diff'):
+                pk = (pk_raw - pk_diff) / correction_reso
+            elif (args.noise_estimate == 'mean_diff' or 'mean_rebin_diff'):
+                if args.kmin_noise_avg is None and linear_binning:
+                    #this is roughly the same range as eBOSS analyses for z=2.2
+                    selection = (k > 0) & (k < 1.5)
+                elif args.kmin_noise_avg is None:
+                    selection = (k > 0) & (k < 0.02)
+                else:
+                    selection = (((k > args.kmin_noise_avg) if args.kmax_noise_avg is not None else 1) & 
+                                 ((k < args.kmax_noise_avg) if args.kmax_noise_avg is not None else 1))
+                mean_pk_diff = np.mean(pk_diff[selection])
+                pk = (pk_raw - mean_pk_diff) / correction_reso
+
+            if args.force_output_in_velocity and linear_binning:
+                #division by 1000 to convert speed_light from m/s to km/s
+                c_kms=constants.speed_light / 1000
+                lambda_mean=np.mean(lambda_new)
+                pk *= c_kms / lambda_mean
+                pk_raw *= c_kms / lambda_mean
+                pk_noise *= c_kms / lambda_mean
+                pk_diff *= c_kms / 1000 / lambda_mean
+                k /= c_kms / lambda_mean
+            
+            # save in fits format
+            if args.out_format == 'fits':
+                header = [{
+                    'name': 'RA',
+                    'value': delta.ra,
+                    'comment': "QSO's Right Ascension [degrees]"
+                }, {
+                    'name': 'DEC',
+                    'value': delta.dec,
+                    'comment': "QSO's Declination [degrees]"
+                }, {
+                    'name': 'Z',
+                    'value': delta.z_qso,
+                    'comment': "QSO's redshift"
+                }, {
+                    'name': 'MEANZ',
+                    'value': mean_z_array[part_index],
+                    'comment': "Absorbers mean redshift"
+                }, {
+                    'name': 'MEANRESO',
+                    'value': delta.mean_reso,
+                    'comment': 'Mean resolution [km/s]'
+                }, {
+                    'name': 'MEANSNR',
+                    'value': delta.mean_snr,
+                    'comment': 'Mean signal to noise ratio'
+                }, {
+                    'name': 'NBMASKPIX',
+                    'value': num_masked_pixels,
+                    'comment': 'Number of masked pixels in the section'
+                }, {
+                    'name': 'LIN_BIN',
+                    'value': linear_binning,
+                    'comment': "analysis was performed on delta with linear binned lambda"
+                }, {
+                    'name': 'LOS_ID',
+                    'value': delta.los_id,
+                    'comment': "line of sight identifier, e.g. THING_ID or TARGETID"
+                },
+                ]
+
+                cols = [k, pk_raw, pk_noise, pk_diff, correction_reso, pk]
+                names = [
+                    'K', 'PK_RAW', 'PK_NOISE', 'PK_DIFF', 'COR_RESO', 'PK'
+                ]
+                comments = [
+                    'Wavenumber', 'Raw power spectrum',
+                    "Noise's power spectrum",
+                    'Noise coadd difference power spectrum',
+                    'Correction resolution function',
+                    'Corrected power spectrum (resolution and noise)'
+                ]
+                if linear_binning and not args.force_output_in_velocity:
+                    baseunit = "AA"
+                else:
+                    baseunit = "km/s"
+                units = [
+                    f'({baseunit})^-1', f'{baseunit}', f'{baseunit}',
+                    f'{baseunit}', f'{baseunit}', f'{baseunit}'
+                ]
+
+                try:
+                    results.write(cols,
+                                  names=names,
+                                  header=header,
+                                  comments=comments,
+                                  units=units)
+                except AttributeError:
+                    userprint("writing to " + args.out_dir + '/Pk1D-' +
+                                           str(file_index) + '.fits.gz')
+                    results = fitsio.FITS((args.out_dir + '/Pk1D-' +
+                                           str(file_index) + '.fits.gz'),
+                                          'rw',
+                                          clobber=True)
+                    results.write(cols,
+                                  names=names,
+                                  header=header,
+                                  comment=comments,
+                                  units=units)
+
+    if (args.out_format == 'fits' and results is not None):
+        results.close()
+
+    return 0
 
 
 def main(cmdargs):
     # pylint: disable-msg=too-many-locals,too-many-branches,too-many-statements
-    """Compute the 1D power spectrum"""
+    """Compute the 1D power spectrum
+    Uses the resolution matrix correction for DESI data"""
     parser = argparse.ArgumentParser(
         formatter_class=argparse.ArgumentDefaultsHelpFormatter,
         description='Compute the 1D power spectrum')
@@ -163,12 +382,11 @@ def main(cmdargs):
                         required=True,
                         help='Output directory')
 
-    parser.add_argument(
-        '--out-format',
-        type=str,
-        default='fits',
-        required=False,
-        help='Output format: root or fits (if root call PyRoot)')
+    parser.add_argument('--out-format',
+                        type=str,
+                        default='fits',
+                        required=False,
+                        help='Output format: ascii or fits')
 
     parser.add_argument('--in-dir',
                         type=str,
@@ -232,7 +450,7 @@ def main(cmdargs):
         default='mean_diff',
         required=False,
         help=('Estimate of Pk_noise '
-              'pipeline/diff/mean_diff/rebin_diff/mean_rebin_diff'))
+              'pipeline/mean_pipeline/diff/mean_diff/rebin_diff/mean_rebin_diff'))
 
     parser.add_argument('--forest-type',
                         type=str,
@@ -240,12 +458,6 @@ def main(cmdargs):
                         required=False,
                         help='Forest used: Lya, SiIV, CIV')
 
-    parser.add_argument('--debug',
-                        action='store_true',
-                        default=False,
-                        required=False,
-                        help='Fill root histograms for debugging')
-
     parser.add_argument(
         '--abs-igm',
         type=str,
@@ -254,283 +466,92 @@ def main(cmdargs):
         help=('Name of the absorption line in picca.constants defining the '
               'redshift of the forest pixels'))
 
-    args = parser.parse_args(cmdargs)
-
-    # Create root file
-    if args.out_format == 'root':
-        # pylint: disable-msg=import-error,import-outside-toplevel
-        # import is done here as ROOT is not a required package for the code
-        # to run, except if args.out_format is set to 'root'
-        from ROOT import TH1D, TFile, TTree, TProfile2D, TProfile
-        store_file = TFile(args.out_dir + "/Testpicca.root", "RECREATE",
-                           "PK 1D studies studies")
-        max_num_bins = 700
-        tree = TTree("Pk1D", "SDSS 1D Power spectrum Ly-a")
-        (z_qso, mean_z, mean_reso, mean_snr, lambda_min_tree, lambda_max_tree,
-         plate, mjd, fiber, num_masked_pixels_tree, num_bins_tree, k_tree,
-         pk_tree, pk_raw_tree, pk_noise_tree, correction_reso_tree,
-         pk_diff_tree) = make_tree(tree, max_num_bins)
-
-        # control histograms
-        if args.forest_type == 'Lya':
-            lambda_min = 1040.
-            lambda_max = 1200.
-        elif args.forest_type == 'SiIV':
-            lambda_min = 1270.
-            lambda_max = 1380.
-        elif args.forest_type == 'CIV':
-            lambda_min = 1410.
-            lambda_max = 1520.
-        hist_delta = TProfile2D('hdelta',
-                                'delta mean as a function of lambda-lambdaRF',
-                                36, 3600., 7200., 16, lambda_min, lambda_max,
-                                -5.0, 5.0)
-        hist_delta_rest_frame = TProfile(
-            'hdelta_RF', 'delta mean as a function of lambdaRF', 320,
-            lambda_min, lambda_max, -5.0, 5.0)
-        hist_delta_obs_frame = TProfile(
-            'hdelta_OBS', 'delta mean as a function of lambdaOBS', 1800, 3600.,
-            7200., -5.0, 5.0)
-        hist_weighted_delta_rest_frame = TProfile(
-            'hdelta_RF_we', 'delta mean weighted as a function of lambdaRF',
-            320, lambda_min, lambda_max, -5.0, 5.0)
-        hist_weighted_delta_obs_frame = TProfile(
-            'hdelta_OBS_we', 'delta mean weighted as a function of lambdaOBS',
-            1800, 3600., 7200., -5.0, 5.0)
-        hist_ivar = TH1D('hivar', '  ivar ', 10000, 0.0, 10000.)
-        hist_snr = TH1D('hsnr', '  snr per pixel ', 100, 0.0, 100.)
-        hist_weighted_delta_rest_frame.Sumw2()
-        hist_weighted_delta_obs_frame.Sumw2()
-
-    # Read deltas
-    if args.in_format == 'fits':
-        files = sorted(glob.glob(args.in_dir + "/*.fits.gz"))
-    elif args.in_format == 'ascii':
-        files = sorted(glob.glob(args.in_dir + "/*.txt"))
+    #additional options
+    parser.add_argument('--num-processors',
+                        type=int,
+                        default=1,
+                        required=False,
+                        help='Number of processors to use for computation')
 
-    num_data = 0
+    parser.add_argument(
+        '--num-noise-exp',
+        default=100,
+        type=int,
+        required=False,
+        help='number of pipeline noise realizations to generate per spectrum')
 
-    # initialize randoms
-    np.random.seed(4)
-
-    # loop over input files
-    for index, file in enumerate(files):
-        if index % 1 == 0:
-            userprint("\rread {} of {} {}".format(index, len(files), num_data),
-                      end="")
-
-        # read fits or ascii file
-        if args.in_format == 'fits':
-            hdul = fitsio.FITS(file)
-            deltas = [
-                Delta.from_fitsio(hdu, pk1d_type=True) for hdu in hdul[1:]
-            ]
-        elif args.in_format == 'ascii':
-            ascii_file = open(file, 'r')
-            deltas = [Delta.from_ascii(line) for line in ascii_file]
-
-        num_data += len(deltas)
-        userprint("\n ndata =  ", num_data)
-        results = None
-
-        # loop over deltas
-        for delta in deltas:
-
-            # Selection over the SNR and the resolution
-            if (delta.mean_snr <= args.SNR_min or
-                    delta.mean_reso >= args.reso_max):
-                continue
+    parser.add_argument('--disable-reso-matrix',
+                        default=False,
+                        action='store_true',
+                        required=False,
+                        help=('do not use the resolution matrix even '
+                              'if it exists and we are on linear binning'))
 
-            # first pixel in forest
-            selected_pixels = 10**delta.log_lambda > args.lambda_obs_min
-            first_pixel_index = (np.argmax(selected_pixels)
-                                 if np.any(selected_pixels) else len(selected_pixels))
+    parser.add_argument(
+        '--force-output-in-velocity',
+        default=False,
+        action='store_true',
+        required=False,
+        help=
+        ('store outputs in units of velocity even for linear binning computations'
+         ))
 
-            # minimum number of pixel in forest
-            min_num_pixels = args.nb_pixel_min
-            if (len(delta.log_lambda) - first_pixel_index) < min_num_pixels:
-                continue
+    parser.add_argument(
+        '--seed',
+        default=4,
+        required=False,
+        type=int,
+        help=
+        ('seed for random number generator, default 4, let system determine seed if set to 0'
+         ))
 
-            # Split in n parts the forest
-            max_num_parts = (len(delta.log_lambda) -
-                             first_pixel_index) // min_num_pixels
-            num_parts = min(args.nb_part, max_num_parts)
-            (mean_z_array, log_lambda_array, delta_array, exposures_diff_array,
-             ivar_array) = split_forest(num_parts, delta.delta_log_lambda,
-                                        delta.log_lambda, delta.delta,
-                                        delta.exposures_diff, delta.ivar,
-                                        first_pixel_index)
-            for index2 in range(num_parts):
 
-                # rebin exposures_diff spectrum
-                if (args.noise_estimate == 'rebin_diff' or
-                        args.noise_estimate == 'mean_rebin_diff'):
-                    exposures_diff_array[index2] = rebin_diff_noise(
-                        delta.delta_log_lambda, log_lambda_array[index2],
-                        exposures_diff_array[index2])
+    parser.add_argument(
+        '--kmin_noise_avg',
+        default=None,
+        required=False,
+        type=float,
+        help=
+        ('minimal mode to take into account when computing noise/diff power average'
+         ))
+    parser.add_argument(
+        '--kmax_noise_avg',
+        default=None,
+        required=False,
+        type=float,
+        help=
+        ('maximal mode to take into account when computing noise/diff power average'
+         ))
 
-                # Fill masked pixels with 0.
-                (log_lambda_new, delta_new, exposures_diff_new, ivar_new,
-                 num_masked_pixels) = fill_masked_pixels(
-                     delta.delta_log_lambda, log_lambda_array[index2],
-                     delta_array[index2], exposures_diff_array[index2],
-                     ivar_array[index2], args.no_apply_filling)
-                if num_masked_pixels > args.nb_pixel_masked_max:
-                    continue
-                if args.out_format == 'root' and args.debug:
-                    compute_mean_delta(log_lambda_new, delta_new, ivar_new,
-                                       delta.z_qso, hist_delta,
-                                       hist_delta_rest_frame,
-                                       hist_delta_obs_frame, hist_ivar,
-                                       hist_snr, hist_weighted_delta_rest_frame,
-                                       hist_weighted_delta_obs_frame)
-
-                # Compute pk_raw
-                k, pk_raw = compute_pk_raw(delta.delta_log_lambda, delta_new)
-
-                # Compute pk_noise
-                run_noise = False
-                if args.noise_estimate == 'pipeline':
-                    run_noise = True
-                pk_noise, pk_diff = compute_pk_noise(delta.delta_log_lambda,
-                                                     ivar_new,
-                                                     exposures_diff_new,
-                                                     run_noise)
+    args = parser.parse_args(cmdargs)
 
-                # Compute resolution correction
-                delta_pixel = (delta.delta_log_lambda * np.log(10.) *
-                               constants.speed_light / 1000.)
-                correction_reso = compute_correction_reso(
-                    delta_pixel, delta.mean_reso, k)
-
-                # Compute 1D Pk
-                if args.noise_estimate == 'pipeline':
-                    pk = (pk_raw - pk_noise) / correction_reso
-                elif (args.noise_estimate == 'diff' or
-                      args.noise_estimate == 'rebin_diff'):
-                    pk = (pk_raw - pk_diff) / correction_reso
-                elif (args.noise_estimate == 'mean_diff' or
-                      args.noise_estimate == 'mean_rebin_diff'):
-                    selection = (k > 0) & (k < 0.02)
-                    if args.noise_estimate == 'mean_rebin_diff':
-                        selection = (k > 0.003) & (k < 0.02)
-                    mean_pk_diff = (sum(pk_diff[selection]) /
-                                    float(len(pk_diff[selection])))
-                    pk = (pk_raw - mean_pk_diff) / correction_reso
-
-                # save in root format
-                if args.out_format == 'root':
-                    z_qso[0] = delta.z_qso
-                    mean_z[0] = mean_z_array[index2]
-                    mean_reso[0] = delta.mean_reso
-                    mean_snr[0] = delta.mean_snr
-                    lambda_min_tree[0] = np.power(10., log_lambda_new[0])
-                    lambda_max_tree[0] = np.power(10., log_lambda_new[-1])
-                    num_masked_pixels_tree[0] = num_masked_pixels
-
-                    plate[0] = delta.plate
-                    mjd[0] = delta.mjd
-                    fiber[0] = delta.fiberid
-
-                    num_bins_tree[0] = min(len(k), max_num_bins)
-                    for index3 in range(num_bins_tree[0]):
-                        k_tree[index3] = k[index3]
-                        pk_raw_tree[index3] = pk_raw[index3]
-                        pk_noise_tree[index3] = pk_noise[index3]
-                        pk_diff_tree[index3] = pk_diff[index3]
-                        pk_tree[index3] = pk[index3]
-                        correction_reso_tree[index3] = correction_reso[index3]
-
-                    tree.Fill()
-
-                # save in fits format
-                if args.out_format == 'fits':
-                    header = [{
-                        'name': 'RA',
-                        'value': delta.ra,
-                        'comment': "QSO's Right Ascension [degrees]"
-                    }, {
-                        'name': 'DEC',
-                        'value': delta.dec,
-                        'comment': "QSO's Declination [degrees]"
-                    }, {
-                        'name': 'Z',
-                        'value': delta.z_qso,
-                        'comment': "QSO's redshift"
-                    }, {
-                        'name': 'MEANZ',
-                        'value': mean_z_array[index2],
-                        'comment': "Absorbers mean redshift"
-                    }, {
-                        'name': 'MEANRESO',
-                        'value': delta.mean_reso,
-                        'comment': 'Mean resolution [km/s]'
-                    }, {
-                        'name': 'MEANSNR',
-                        'value': delta.mean_snr,
-                        'comment': 'Mean signal to noise ratio'
-                    }, {
-                        'name': 'NBMASKPIX',
-                        'value': num_masked_pixels,
-                        'comment': 'Number of masked pixels in the section'
-                    }, {
-                        'name': 'PLATE',
-                        'value': delta.plate,
-                        'comment': "Spectrum's plate id"
-                    }, {
-                        'name':
-                            'MJD',
-                        'value':
-                            delta.mjd,
-                        'comment': ('Modified Julian Date,date the spectrum '
-                                    'was taken')
-                    }, {
-                        'name': 'FIBER',
-                        'value': delta.fiberid,
-                        'comment': "Spectrum's fiber number"
-                    }]
-
-                    cols = [k, pk_raw, pk_noise, pk_diff, correction_reso, pk]
-                    names = [
-                        'k', 'Pk_raw', 'Pk_noise', 'Pk_diff', 'cor_reso', 'Pk'
-                    ]
-                    comments = [
-                        'Wavenumber', 'Raw power spectrum',
-                        "Noise's power spectrum",
-                        'Noise coadd difference power spectrum',
-                        'Correction resolution function',
-                        'Corrected power spectrum (resolution and noise)'
-                    ]
-                    units = [
-                        '(km/s)^-1', 'km/s', 'km/s', 'km/s', 'km/s', 'km/s'
-                    ]
-
-                    try:
-                        results.write(cols,
-                                      names=names,
-                                      header=header,
-                                      comments=comments,
-                                      units=units)
-                    except AttributeError:
-                        results = fitsio.FITS(
-                            (args.out_dir + '/Pk1D-' + str(index) + '.fits.gz'),
-                            'rw',
-                            clobber=True)
-                        results.write(cols,
-                                      names=names,
-                                      header=header,
-                                      comment=comments,
-                                      units=units)
-        if (args.out_format == 'fits' and results is not None):
-            results.close()
-
-    # Store root file results
-    if args.out_format == 'root':
-        store_file.Write()
+    if args.seed==0:
+        seed=None
+    else:
+        seed=args.seed
+    # Read deltas
+    if args.in_format == 'fits':
+        files = sorted(glob.glob(args.in_dir + "/*.fits.gz"))
+    elif args.in_format == 'ascii':
+        files = sorted(glob.glob(args.in_dir + "/*.txt"))
 
+    # initialize randoms
+    np.random.seed(seed)
+    userprint(f"Computing Pk1d for {args.in_dir}")
+    args.len_files = len(files)
+    #create output dir if it does not exist
+    os.makedirs(args.out_dir, exist_ok=True)
+
+    print([[i, f] for i, f in enumerate(files)])
+    if args.num_processors > 1:
+        pool = Pool(args.num_processors)
+        index_file_args = [(i, f, args) for i, f in enumerate(files)]
+        pool.map(process_all_files, index_file_args)
+    else:
+        [process_all_files((i, f, args)) for i, f in enumerate(files)]
     userprint("all done ")
 
 
 if __name__ == '__main__':
-    cmdargs=sys.argv[1:]
+    cmdargs = sys.argv[1:]
     main(cmdargs)
diff --git a/examples/delta_extraction/desi_pk1d_fuji_produce_delta.ini b/examples/delta_extraction/desi_pk1d_fuji_produce_delta.ini
new file mode 100644
index 000000000..d956c056c
--- /dev/null
+++ b/examples/delta_extraction/desi_pk1d_fuji_produce_delta.ini
@@ -0,0 +1,33 @@
+[general]
+overwrite = False
+out dir = $out_fuji_pk1d/Delta_pk1d/
+logging level console = DEBUG
+logging level file = DEBUG
+
+[data]
+type = DesiHealpix
+analysis type = PK 1D
+input directory = /global/cfs/cdirs/desi/spectro/redux/fuji/healpix/
+catalogue = /global/cfs/cdirs/desi/users/edmondc/QSO_catalog/fuji/QSO_cat_fuji_healpix_for_lyalpha.fits
+lambda min = 3500.0
+lambda max = 7500.0
+lambda min rest frame = 1040.0
+lambda max rest frame = 1200.0
+delta lambda = 0.8
+lambda abs IGM = LYA
+wave solution = lin
+#use non-coadded spectra = True     #will use the spectra files, takes way longer for I/O, coadding reasons
+
+[corrections]
+num corrections = 0
+
+[masks]
+num masks = 0
+
+[expected flux]
+type = Dr16ExpectedFlux
+iter out prefix = delta_attributes
+num iterations = 7
+use constant weight = True
+num processors = 4
+order = 0
diff --git a/py/picca/data.py b/py/picca/data.py
index 8e069029d..4e8e028a2 100644
--- a/py/picca/data.py
+++ b/py/picca/data.py
@@ -7,6 +7,7 @@
 import numpy as np
 import iminuit
 import fitsio
+import warnings
 
 from . import constants
 from .utils import userprint, unred
@@ -58,7 +59,7 @@ class QSO(object):
         get_angle_between: Computes the angular separation between two quasars.
     """
 
-    def __init__(self, thingid, ra, dec, z_qso, plate, mjd, fiberid):
+    def __init__(self, los_id, ra, dec, z_qso, plate, mjd, fiberid):
         """Initializes class instance.
 
         Args:
@@ -91,7 +92,10 @@ def __init__(self, thingid, ra, dec, z_qso, plate, mjd, fiberid):
         self.cos_dec = np.cos(dec)
 
         self.z_qso = z_qso
-        self.thingid = thingid
+        self.los_id = los_id
+        #this is for legacy purposes only
+        self.thingid = los_id
+        warnings.warn("currently a thingid entry is created in QSO.__init__, this feature will be removed", DeprecationWarning)
 
         # variables computed in function io.read_objects
         self.weight = None
@@ -931,9 +935,15 @@ class Delta(QSO):
         mean_snr: float
             Mean signal-to-noise ratio in the forest
         mean_reso: float
-            Mean resolution of the forest
+            Mean resolution of the forest in units of velocity (FWHM)
         mean_z: float
             Mean redshift of the forest
+        mean_reso_pix: float
+            Mean resolution of the forest in units of pixels (FWHM)
+        mean_resolution_matrix: array of floats or None
+            Mean (over wavelength) resolution matrix for that forest
+        resolution_matrix: 2d array of floats or None
+            Wavelength dependent resolution matrix for that forest
         delta_log_lambda: float
             Variation of the logarithm of the wavelength between two pixels
         z: array of floats or None
@@ -957,14 +967,15 @@ class Delta(QSO):
 
     """
 
-    def __init__(self, thingid, ra, dec, z_qso, plate, mjd, fiberid, log_lambda,
+    def __init__(self, los_id, ra, dec, z_qso, plate, mjd, fiberid, log_lambda,
                  weights, cont, delta, order, ivar, exposures_diff, mean_snr,
-                 mean_reso, mean_z, delta_log_lambda):
+                 mean_reso, mean_z, resolution_matrix=None,
+                 mean_resolution_matrix=None, mean_reso_pix=None):
         """Initializes class instances.
 
         Args:
-            thingid: integer
-                Thingid of the observation.
+            los_id: integer
+                Thingid or Targetid of the observation.
             ra: float
                 Right-ascension of the quasar (in radians).
             dec: float
@@ -997,10 +1008,16 @@ def __init__(self, thingid, ra, dec, z_qso, plate, mjd, fiberid, log_lambda,
                 Mean resolution of the forest
             mean_z: float
                 Mean redshift of the forest
+            mean_reso_pix: float
+                Mean resolution of the forest in units of pixels (FWHM)
+            mean_resolution_matrix: array of floats or None
+                Mean (over wavelength) resolution matrix for that forest
+            resolution_matrix: 2d array of floats or None
+                Wavelength dependent resolution matrix for that forest
             delta_log_lambda: float
                 Variation of the logarithm of the wavelength between two pixels
         """
-        QSO.__init__(self, thingid, ra, dec, z_qso, plate, mjd, fiberid)
+        QSO.__init__(self, los_id, ra, dec, z_qso, plate, mjd, fiberid)
         self.log_lambda = log_lambda
         self.weights = weights
         self.cont = cont
@@ -1011,7 +1028,9 @@ def __init__(self, thingid, ra, dec, z_qso, plate, mjd, fiberid, log_lambda,
         self.mean_snr = mean_snr
         self.mean_reso = mean_reso
         self.mean_z = mean_z
-        self.delta_log_lambda = delta_log_lambda
+        self.resolution_matrix = resolution_matrix
+        self.mean_resolution_matrix = mean_resolution_matrix
+        self.mean_reso_pix = mean_reso_pix
 
         # variables computed in function io.read_deltas
         self.z = None
@@ -1057,7 +1076,6 @@ def from_fitsio(cls, hdu, pk1d_type=False):
         if 'LOGLAM' in hdu.get_colnames():
             log_lambda = hdu['LOGLAM'][:].astype(float)
         elif 'LAMBDA' in hdu.get_colnames():
-            userprint("no LOGLAM found, trying to read linear_binned_lambda")
             log_lambda = np.log10(hdu['LAMBDA'][:].astype(float))
         else:
             raise KeyError("Did not find LOGLAM or LAMBDA in delta file")
@@ -1067,8 +1085,22 @@ def from_fitsio(cls, hdu, pk1d_type=False):
             exposures_diff = hdu['DIFF'][:].astype(float)
             mean_snr = header['MEANSNR']
             mean_reso = header['MEANRESO']
+            try:
+                mean_reso_pix = header['MEANRESO_PIX']
+            except (KeyError, ValueError):
+                mean_reso_pix = None
+
             mean_z = header['MEANZ']
-            delta_log_lambda = header['DLL']
+            try:
+                #transposing here gives back the actual reso matrix which has been stored transposed
+                resolution_matrix = hdu['RESOMAT'][:].T.astype(float)
+                if resolution_matrix is not None:
+                    mean_resolution_matrix = np.mean(resolution_matrix, axis=1)
+                else:
+                    mean_resolution_matrix = None
+            except (KeyError, ValueError):
+                resolution_matrix = None
+                mean_resolution_matrix = None
             weights = None
             cont = None
         else:
@@ -1076,8 +1108,10 @@ def from_fitsio(cls, hdu, pk1d_type=False):
             exposures_diff = None
             mean_snr = None
             mean_reso = None
-            delta_log_lambda = None
             mean_z = None
+            resolution_matrix = None
+            mean_resolution_matrix = None
+            mean_reso_pix = None
             weights = hdu['WEIGHT'][:].astype(float)
             cont = hdu['CONT'][:].astype(float)
 
@@ -1104,7 +1138,8 @@ def from_fitsio(cls, hdu, pk1d_type=False):
 
         return cls(los_id, ra, dec, z_qso, plate, mjd, fiberid, log_lambda,
                    weights, cont, delta, order, ivar, exposures_diff, mean_snr,
-                   mean_reso, mean_z, delta_log_lambda)
+                   mean_reso, mean_z, resolution_matrix,
+                   mean_resolution_matrix, mean_reso_pix)
 
     @classmethod
     def from_ascii(cls, line):
diff --git a/py/picca/delta_extraction/astronomical_objects/desi_pk1d_forest.py b/py/picca/delta_extraction/astronomical_objects/desi_pk1d_forest.py
index ec7396387..4494dd620 100644
--- a/py/picca/delta_extraction/astronomical_objects/desi_pk1d_forest.py
+++ b/py/picca/delta_extraction/astronomical_objects/desi_pk1d_forest.py
@@ -1,8 +1,12 @@
 """This module defines the class DesiPk1dForest to represent DESI forests
 in the Pk1D analysis
 """
+import numpy as np
+
 from picca.delta_extraction.astronomical_objects.desi_forest import DesiForest
 from picca.delta_extraction.astronomical_objects.pk1d_forest import Pk1dForest
+from picca.delta_extraction.astronomical_objects.forest import Forest
+from picca.delta_extraction.errors import AstronomicalObjectError
 
 
 class DesiPk1dForest(DesiForest, Pk1dForest):
@@ -119,6 +123,9 @@ class DesiPk1dForest(DesiForest, Pk1dForest):
 
     reso: array of floats or None (from Pk1dForest)
     Resolution of the forest
+    
+    resolution_matrix: array of floats or None (from Pk1dForest)
+    Resolution matrix of the forests
     """
     def __init__(self, **kwargs):
         """Initialize instance
@@ -126,10 +133,148 @@ def __init__(self, **kwargs):
         Arguments
         ---------
         **kwargs: dict
-        Dictionary contiaing the information
+        Dictionary containing the information
 
         Raise
         -----
         AstronomicalObjectError if there are missing variables
         """
+
+        self.resolution_matrix = kwargs.get("resolution_matrix")
+        #potentially change this in case we ever want log-binning with DESI Pk1d data
+        #then would need a check of self.wave_solution
+        if self.resolution_matrix is None:
+            raise AstronomicalObjectError(
+                "Error constructing DesiPk1dForest. "
+                "Missing variable 'resolution_matrix'")
+        if "resolution_matrix" in kwargs:
+            del kwargs["resolution_matrix"]
+
+        # call parent constructors
         super().__init__(**kwargs)
+        self.consistency_check()
+
+
+    def consistency_check(self):
+        """Consistency checks after __init__"""
+        super().consistency_check()
+        if self.resolution_matrix.shape[1] != self.flux.shape[0]:
+            raise AstronomicalObjectError(
+                "Error constructing DesiPk1dForest. 'resolution_matrix', "
+                "and 'flux' don't have the "
+                "same size")
+        if "resolution_matrix" not in Forest.mask_fields:
+            Forest.mask_fields += ["resolution_matrix"]
+
+    def coadd(self, other):
+        """Coadd the information of another forest.
+
+        Extends the coadd method of Forest to also include information
+        about the exposures_diff and reso arrays
+
+        Arguments
+        ---------
+        other: Pk1dForest
+        The forest instance to be coadded.
+
+        Raise
+        -----
+        AstronomicalObjectError if other is not a DesiPk1dForest instance
+        """
+        if not isinstance(other, DesiPk1dForest):
+            raise AstronomicalObjectError(
+                "Error coadding DesiPk1dForest. Expected "
+                "DesiPk1dForest instance in other. Found: "
+                f"{type(other)}")
+
+        if other.resolution_matrix.size > 0 and self.resolution_matrix.size > 0:
+            self.resolution_matrix = np.append(self.resolution_matrix,
+                                               other.resolution_matrix,
+                                               axis=1)
+        elif self.resolution_matrix.size == 0:
+            self.resolution_matrix = other.resolution_matrix
+
+        # coadd the deltas by rebinning
+        super().coadd(other)
+
+    def get_data(self):
+        """Get the data to be saved in a fits file.
+
+        Extends the get_data method of Forest to also include data for
+        ivar and exposures_diff.
+
+        Return
+        ------
+        cols: list of arrays
+        Data of the different variables
+
+        names: list of str
+        Names of the different variables
+
+        units: list of str
+        Units of the different variables
+
+        comments: list of str
+        Comments attached to the different variables
+        """
+        cols, names, units, comments = super().get_data()
+
+        #transposing here is necessary to store in fits file
+        cols += [self.resolution_matrix.T]
+        names += ["RESOMAT"]
+        comments += ["Transposed Masked resolution matrix"]
+        units += [""]
+        return cols, names, units, comments
+
+    def rebin(self):
+        """Rebin the arrays and update control variables
+
+        Extends the rebon method of Forest to also rebin exposures_diff and compute
+        the control variable mean_reso.
+
+        Rebinned arrays are flux, ivar, lambda_ or log_lambda,
+        transmission_correctionm, exposures_diff, and reso. Control variables
+        are mean_snr and mean_reso.
+
+        Return
+        ------
+        bins: array of float
+        Binning solution to be used for the rebinning
+
+        rebin_ivar: array of float
+        Rebinned version of ivar
+
+        orig_ivar: array of float
+        Original version of ivar (before applying the function)
+
+        w1: array of bool
+        Masking array for the bins solution
+
+        w2: array of bool
+        Masking array for the rebinned ivar solution
+
+        Raise
+        -----
+        AstronomicalObjectError if Forest.wave_solution is not 'lin' or 'log'
+        """
+        bins, rebin_ivar, orig_ivar, w1, w2 = super().rebin()
+        if len(rebin_ivar) == 0:
+            self.resolution_matrix = np.array([[]])
+            return [], [], [], [], []
+
+        # apply mask due to cuts in bin
+        self.resolution_matrix = self.resolution_matrix[:, w1]
+
+        # rebin resolution_matrix
+        rebin_reso_matrix_aux = np.zeros(
+            (self.resolution_matrix.shape[0], bins.max() + 1))
+        for index, reso_matrix_col in enumerate(self.resolution_matrix):
+            rebin_reso_matrix_aux[index, :] = np.bincount(
+                bins, weights=orig_ivar[w1] * reso_matrix_col)
+        # apply mask due to rebinned inverse vairane
+        self.resolution_matrix = rebin_reso_matrix_aux[:, w2] / rebin_ivar[
+            np.newaxis, w2]
+
+        # return weights and binning solution to be used by child classes if
+        # required
+        return bins, rebin_ivar, orig_ivar, w1, w2
diff --git a/py/picca/delta_extraction/astronomical_objects/forest.py b/py/picca/delta_extraction/astronomical_objects/forest.py
index 1d3556e1e..9e732f7a9 100644
--- a/py/picca/delta_extraction/astronomical_objects/forest.py
+++ b/py/picca/delta_extraction/astronomical_objects/forest.py
@@ -1,6 +1,7 @@
 """This module defines the abstract class Forest from which all
 objects representing a forest must inherit from
 """
+import logging
 import numpy as np
 
 from picca.delta_extraction.astronomical_object import AstronomicalObject
@@ -126,6 +127,8 @@ def __init__(self, **kwargs):
         -----
         AstronomicalObjectError if there are missing variables
         """
+        self.logger = logging.getLogger(__name__)
+
         Forest.class_variable_check()
 
         if Forest.wave_solution == "log":
@@ -548,7 +551,7 @@ def rebin(self):
     def set_class_variables(cls, lambda_min, lambda_max,
                             lambda_min_rest_frame,
                             lambda_max_rest_frame,
-                            pixel_step, wave_solution):
+                            pixel_step, pixel_step_rest_frame, wave_solution):
         """Set class variables
         Arguments
         ---------
@@ -574,9 +577,9 @@ def set_class_variables(cls, lambda_min, lambda_max,
                 np.log10(lambda_max) + pixel_step/2,
                 pixel_step)
             cls.log_lambda_rest_frame_grid = np.arange(
-                np.log10(lambda_min_rest_frame) + pixel_step/2,
+                np.log10(lambda_min_rest_frame) + pixel_step_rest_frame/2,
                 np.log10(lambda_max_rest_frame),
-                pixel_step)
+                pixel_step_rest_frame)
             cls.mask_fields = defaults.get("mask fields log").copy()
         elif wave_solution == "lin":
             cls.lambda_grid = np.arange(
@@ -584,9 +587,9 @@ def set_class_variables(cls, lambda_min, lambda_max,
                 lambda_max + pixel_step/2,
                 pixel_step)
             cls.lambda_rest_frame_grid = np.arange(
-                lambda_min_rest_frame + pixel_step/2,
+                lambda_min_rest_frame + pixel_step_rest_frame/2,
                 lambda_max_rest_frame,
-                pixel_step)
+                pixel_step_rest_frame)
             cls.mask_fields = defaults.get("mask fields lin").copy()
         else:
             raise AstronomicalObjectError("Error in setting Forest class "
diff --git a/py/picca/delta_extraction/astronomical_objects/pk1d_forest.py b/py/picca/delta_extraction/astronomical_objects/pk1d_forest.py
index adc5d5074..f74f992ba 100644
--- a/py/picca/delta_extraction/astronomical_objects/pk1d_forest.py
+++ b/py/picca/delta_extraction/astronomical_objects/pk1d_forest.py
@@ -137,11 +137,18 @@ def __init__(self, **kwargs):
                                           "Missing variable 'reso'")
         del kwargs["reso"]
 
+        self.reso_pix = kwargs.get("reso_pix")
+        if self.reso_pix is None:
+            raise AstronomicalObjectError("Error constructing Pk1dForest. "
+                                          "Missing variable 'reso_pix'")
+        del kwargs["reso_pix"]
+
         # call parent constructor
         super().__init__(**kwargs)
 
         # compute mean quality variables
         self.mean_reso = self.reso.mean()
+        self.mean_reso_pix = self.reso_pix.mean()
         if Forest.wave_solution == "log":
             self.mean_z = (
                 (np.power(10., self.log_lambda[len(self.log_lambda) - 1]) +
@@ -180,6 +187,9 @@ def consistency_check(self):
             Forest.mask_fields += ["exposures_diff"]
         if "reso" not in Forest.mask_fields:
             Forest.mask_fields += ["reso"]
+        if "reso_pix" not in Forest.mask_fields:
+            Forest.mask_fields += ["reso_pix"]
+        
 
     def coadd(self, other):
         """Coadd the information of another forest.
@@ -204,7 +214,7 @@ def coadd(self, other):
         self.exposures_diff = np.append(self.exposures_diff,
                                         other.exposures_diff)
         self.reso = np.append(self.reso, other.reso)
-
+        self.reso_pix = np.append(self.reso_pix, other.reso_pix)
         # coadd the deltas by rebinning
         super().coadd(other)
 
@@ -230,13 +240,15 @@ def get_data(self):
         """
         cols, names, units, comments = super().get_data()
 
-        cols += [self.ivar, self.exposures_diff]
-        names += ["IVAR", "DIFF"]
+        cols += [self.ivar, self.exposures_diff, self.reso, self.reso_pix]
+        names += ["IVAR", "DIFF", "RESO", "RESO_PIX"]
         comments += [
             "Inverse variance. Check input spectra for units",
-            "Difference. Check input spectra for units"
+            "Difference. Check input spectra for units",
+            "Resolution estimate (FWHM) for each pixel in units of km/s"
+            "Resolution estimate (sigma) for each pixel in units of pixel size"
         ]
-        units += ["Flux units", "Flux units"]
+        units += ["Flux units", "Flux units", "", ""]
 
         return cols, names, units, comments
 
@@ -261,7 +273,12 @@ def get_header(self):
             {
                 'name': 'MEANRESO',
                 'value': self.mean_reso,
-                'comment': 'Mean resolution'
+                'comment': 'Mean resolution (km/s)'
+            },
+            {
+                'name': 'MEANRESO_PIX',
+                'value': self.mean_reso_pix,
+                'comment': 'Mean resolution (pixels)'
             },
         ]
 
@@ -302,29 +319,37 @@ def rebin(self):
         if len(rebin_ivar) == 0:
             self.exposures_diff = np.array([])
             self.reso = np.array([])
+            self.reso_pix = np.array([])
             return [], [], [], [], []
 
         # apply mask due to cuts in bin
         self.exposures_diff = self.exposures_diff[w1]
         self.reso = self.reso[w1]
+        self.reso_pix = self.reso_pix[w1]
 
         # rebin exposures_diff and reso
         rebin_exposures_diff = np.zeros(bins.max() + 1)
         rebin_reso = np.zeros(bins.max() + 1)
+        rebin_reso_pix = np.zeros(bins.max() + 1)
         rebin_exposures_diff_aux = np.bincount(bins,
                                                weights=orig_ivar[w1] *
                                                self.exposures_diff)
         rebin_reso_aux = np.bincount(bins, weights=orig_ivar[w1] * self.reso)
+        rebin_reso_pix_aux = np.bincount(bins, weights=orig_ivar[w1] * self.reso_pix)
         rebin_exposures_diff[:len(rebin_exposures_diff_aux
                                  )] += rebin_exposures_diff_aux
         rebin_reso[:len(rebin_reso_aux)] += rebin_reso_aux
+        rebin_reso_pix[:len(rebin_reso_pix_aux)] += rebin_reso_pix_aux
+
 
         # apply mask due to rebinned inverse vairane
         self.exposures_diff = rebin_exposures_diff[w2] / rebin_ivar[w2]
         self.reso = rebin_reso[w2] / rebin_ivar[w2]
+        self.reso_pix = rebin_reso_pix[w2] / rebin_ivar[w2]
 
         # finally update control variables
         self.mean_reso = self.reso.mean()
+        self.mean_reso_pix = self.reso_pix.mean()
         if Forest.wave_solution == "log":
             self.mean_z = (
                 (np.power(10., self.log_lambda[len(self.log_lambda) - 1]) +
diff --git a/py/picca/delta_extraction/data.py b/py/picca/delta_extraction/data.py
index 455d889c3..d448aad47 100644
--- a/py/picca/delta_extraction/data.py
+++ b/py/picca/delta_extraction/data.py
@@ -13,6 +13,7 @@
 from picca.delta_extraction.utils import ABSORBER_IGM
 
 accepted_options = ["analysis type", "delta lambda", "delta log lambda",
+                    "delta lambda rest frame",
                     "input directory",
                     "lambda abs IGM",
                     "lambda max", "lambda max rest frame",
@@ -126,11 +127,22 @@ def __parse_config(self, config):
             if pixel_step is None:
                 raise DataError("Missing argument 'delta log lambda' required by "
                                 "Data when 'wave solution' is set to 'log'")
+            pixel_step_rest_frame = config.getfloat("delta log lambda rest frame")
+            if pixel_step_rest_frame is None:
+                pixel_step_rest_frame = pixel_step
+                self.logger.info("'delta log lambda rest frame' not set, using "
+                                 f"the same value as for 'delta log lambda' ({pixel_step_rest_frame})")
+            
         elif wave_solution == "lin":
             pixel_step = config.getfloat("delta lambda")
             if pixel_step is None:
                 raise DataError("Missing argument 'delta lambda' required by "
                                 "Data when 'wave solution' is set to 'lin'")
+            pixel_step_rest_frame = config.getfloat("delta lambda rest frame")
+            if pixel_step_rest_frame is None:
+                pixel_step_rest_frame = pixel_step
+                self.logger.info("'delta lambda rest frame' not set, using"
+                                 f"the same value as for 'delta lambda' ({pixel_step_rest_frame})")
         else:
             raise DataError("Forest.wave_solution must be either "
                             "'log' or 'lin'")
@@ -151,7 +163,8 @@ def __parse_config(self, config):
         Forest.set_class_variables(lambda_min, lambda_max,
                                    lambda_min_rest_frame,
                                    lambda_max_rest_frame,
-                                   pixel_step, wave_solution)
+                                   pixel_step, pixel_step_rest_frame,
+                                   wave_solution)
 
         # instance variables
         self.analysis_type = config.get("analysis type")
@@ -287,6 +300,13 @@ def filter_forests(self):
                 self.logger.progress(
                     f"Rejected forest with los_id {forest.los_id} "
                     "due to finding nan")
+            elif self.analysis_type=='PK 1D' and forest.mean_snr<1:
+                #TODO: add variable for the SNR cut, actually only required for constant weights...
+                self.add_to_rejection_log(forest.get_header(), forest.flux.size,
+                                          f"low SNR ({forest.mean_snr})")
+                self.logger.progress(
+                    f"Rejected forest with los_id {forest.los_id} "
+                    "due to low SNR")
             else:
                 continue
 
diff --git a/py/picca/delta_extraction/data_catalogues/desi_healpix.py b/py/picca/delta_extraction/data_catalogues/desi_healpix.py
index cb9ac29f7..a4671a12c 100644
--- a/py/picca/delta_extraction/data_catalogues/desi_healpix.py
+++ b/py/picca/delta_extraction/data_catalogues/desi_healpix.py
@@ -1,6 +1,7 @@
 """This module defines the class DesiData to load DESI data
 """
 import logging
+import os
 
 import fitsio
 import healpy
@@ -11,7 +12,15 @@
 from picca.delta_extraction.astronomical_objects.forest import Forest
 from picca.delta_extraction.data_catalogues.desi_data import DesiData, defaults, accepted_options
 from picca.delta_extraction.errors import DataError
-from picca.delta_extraction.utils_pk1d import spectral_resolution_desi
+from picca.delta_extraction.utils_pk1d import spectral_resolution_desi, exp_diff_desi
+
+accepted_options = sorted(
+    list(set(accepted_options + ["use non-coadded spectra"])))
+
+defaults.update({
+    "use non-coadded spectra": False,
+})
+
 
 class DesiHealpix(DesiData):
     """Reads the spectra from DESI using healpix mode and formats its data as a
@@ -62,8 +71,32 @@ def __init__(self, config):
         """
         self.logger = logging.getLogger(__name__)
 
+        self.use_non_coadded_spectra = None
+        self.__parse_config(config)
+        #init of DesiData needs to come last, as it contains the actual data reading and thus needs all config
         super().__init__(config)
 
+
+
+    def __parse_config(self, config):
+        """Parse the configuration options
+
+        Arguments
+        ---------
+        config: configparser.SectionProxy
+        Parsed options to initialize class
+
+        Raise
+        -----
+        DataError upon missing required variables
+        """
+        self.use_non_coadded_spectra = config.getboolean(
+            "use non-coadded spectra")
+        if self.use_non_coadded_spectra is None:
+            raise DataError(
+                "Missing argument 'use non-coadded spectra' required by DesiHealpix"
+            )
+
     def read_data(self):
         """Read the data.
 
@@ -83,9 +116,11 @@ def read_data(self):
         """
         in_nside = 64
 
-
         healpix = [
-            healpy.ang2pix(in_nside, np.pi / 2 - row["DEC"], row["RA"], nest=True)
+            healpy.ang2pix(in_nside,
+                           np.pi / 2 - row["DEC"],
+                           row["RA"],
+                           nest=True)
             for row in self.catalogue
         ]
         self.catalogue["HEALPIX"] = healpix
@@ -95,17 +130,21 @@ def read_data(self):
         forests_by_targetid = {}
         is_sv = True
         for (index,
-             (healpix, survey)), group in zip(enumerate(grouped_catalogue.groups.keys),
-                                    grouped_catalogue.groups):
+             (healpix,
+              survey)), group in zip(enumerate(grouped_catalogue.groups.keys),
+                                     grouped_catalogue.groups):
 
             if survey not in ["sv", "sv1", "sv2", "sv3"]:
                 is_sv = False
 
             input_directory = f'{self.input_directory}/{survey}/dark'
+            coadd_name = "spectra" if self.use_non_coadded_spectra else "coadd"
             filename = (
-                f"{input_directory}/{healpix//100}/{healpix}/coadd-{survey}-"
+                f"{input_directory}/{healpix//100}/{healpix}/{coadd_name}-{survey}-"
                 f"dark-{healpix}.fits")
 
+            #TODO: not sure if we want the dark survey to be hard coded in here, probably won't run on anything else, but still
+
             self.logger.progress(
                 f"Read {index} of {len(grouped_catalogue.groups.keys)}. "
                 f"num_data: {len(forests_by_targetid)}")
@@ -152,6 +191,9 @@ def read_file(self, filename, catalogue, forests_by_targetid):
         colors = ["B", "R"]
         if "Z_FLUX" in hdul:
             colors.append("Z")
+
+        reso_from_truth = False
+        no_scores_available = False
         for color in colors:
             spec = {}
             try:
@@ -163,14 +205,33 @@ def read_file(self, filename, catalogue, forests_by_targetid):
                 for key in ["FLUX", "IVAR"]:
                     spec[key][w] = 0.
                 if self.analysis_type == "PK 1D":
+                    if self.use_non_coadded_spectra and "SCORES" in hdul:
+                        # Calibration factor given in https://desi.lbl.gov/trac/browser/code/desimodel/trunk/data/tsnr/
+                        spec['TEFF_LYA'] = 11.80090901380597 * hdul['SCORES'][f'TSNR2_LYA_{color}'].read()
+                    else:
+                        spec['TEFF_LYA'] = np.ones(spec["FLUX"].shape[0])
+                        if self.use_non_coadded_spectra and not no_scores_available:
+                            self.logger.warning("SCORES are missing, Teff information (and thus DIFF) will be garbage")
+                        no_scores_available=True
+
                     if f"{color}_RESOLUTION" in hdul:
                         spec["RESO"] = hdul[f"{color}_RESOLUTION"].read()
                     else:
-                        raise DataError(
-                            "Error while reading {color} band from "
-                            "{filename}. Analysis type is  'PK 1D', "
-                            "but file does not contain HDU "
-                            f"'{color}_RESOLUTION' ")
+                        basename_truth=os.path.basename(filename).replace('spectra-','truth-')
+                        pathname_truth=os.path.dirname(filename)
+                        filename_truth=f"{pathname_truth}/{basename_truth}"
+                        if os.path.exists(filename_truth):
+                            with fitsio.FITS(filename_truth) as hdul_truth:
+                                spec["RESO"] = hdul_truth[f"{color}_RESOLUTION"].read()
+                            if not reso_from_truth:
+                                self.logger.debug("no resolution in files, reading from truth files")
+                            reso_from_truth=True
+                        else:
+                            raise DataError(
+                                "Error while reading {color} band from "
+                                "{filename}. Analysis type is  'PK 1D', "
+                                "but file does not contain HDU "
+                                f"'{color}_RESOLUTION' ")
                 spectrographs_data[color] = spec
             except OSError:
                 self.logger.warning(
@@ -196,10 +257,16 @@ def read_file(self, filename, catalogue, forests_by_targetid):
                 w_t = w_t[0]
             # Construct DesiForest instance
             # Fluxes from the different spectrographs will be coadded
-            for spec_name, spec in spectrographs_data.items():
-                ivar = spec['IVAR'][w_t].copy()
-                flux = spec['FLUX'][w_t].copy()
-
+            for spec in spectrographs_data.values():
+                if self.use_non_coadded_spectra:
+                    ivar = np.atleast_2d(spec['IVAR'][w_t])
+                    ivar_coadded_flux = np.atleast_2d(ivar*spec['FLUX'][w_t]).sum(axis=0)
+                    ivar = ivar.sum(axis=0)
+                    flux = (ivar_coadded_flux / ivar)
+                else:
+                    flux = spec['FLUX'][w_t].copy()
+                    ivar = spec['IVAR'][w_t].copy()
+                
                 args = {
                     "flux": flux,
                     "ivar": ivar,
@@ -219,17 +286,29 @@ def read_file(self, filename, catalogue, forests_by_targetid):
                 if self.analysis_type == "BAO 3D":
                     forest = DesiForest(**args)
                 elif self.analysis_type == "PK 1D":
-                    reso_sum = spec['RESO'][w_t].copy()
-                    reso_in_km_per_s = spectral_resolution_desi(
+                    if self.use_non_coadded_spectra and not no_scores_available:
+                        exposures_diff = exp_diff_desi(spec, w_t)
+                    else:
+                        exposures_diff = np.zeros(spec['WAVELENGTH'].shape)
+                    if not reso_from_truth:
+                        if len(spec['RESO'][w_t].shape)<3:
+                            reso_sum = spec['RESO'][w_t].copy()
+                        else:
+                            reso_sum = spec['RESO'][w_t].sum(axis=0)
+                    else:
+                        reso_sum = spec['RESO'][:, :]
+                    reso_in_pix, reso_in_km_per_s = spectral_resolution_desi(
                         reso_sum, spec['WAVELENGTH'])
-                    exposures_diff = np.zeros(spec['WAVELENGTH'].shape)
-
                     args["exposures_diff"] = exposures_diff
                     args["reso"] = reso_in_km_per_s
+                    args["resolution_matrix"] = reso_sum
+                    args["reso_pix"] = reso_in_pix
+
                     forest = DesiPk1dForest(**args)
                 else:
-                    raise DataError("Unkown analysis type. Expected 'BAO 3D'"
-                                    f"or 'PK 1D'. Found '{self.analysis_type}'")
+                    raise DataError(
+                        "Unkown analysis type. Expected 'BAO 3D'"
+                        f"or 'PK 1D'. Found '{self.analysis_type}'")
 
                 # rebin arrays
                 # this needs to happen after all arrays are initialized by
diff --git a/py/picca/delta_extraction/data_catalogues/desi_tile.py b/py/picca/delta_extraction/data_catalogues/desi_tile.py
index c7ba3e3f6..ceaedac26 100644
--- a/py/picca/delta_extraction/data_catalogues/desi_tile.py
+++ b/py/picca/delta_extraction/data_catalogues/desi_tile.py
@@ -13,16 +13,20 @@
 from picca.delta_extraction.astronomical_objects.forest import Forest
 from picca.delta_extraction.data_catalogues.desi_data import DesiData, defaults, accepted_options
 from picca.delta_extraction.errors import DataError
-from picca.delta_extraction.utils_pk1d import spectral_resolution_desi
+from picca.delta_extraction.utils_pk1d import spectral_resolution_desi, exp_diff_desi
 
-accepted_options = sorted(list(set(accepted_options+[
-    "use all", "use single nights"])))
+accepted_options = sorted(
+    list(
+        set(accepted_options +
+            ["use all", "use single nights", "use non-coadded spectra"])))
 
 defaults.update({
     "use all": False,
     "use single nights": False,
+    "use non-coadded spectra": False,
 })
 
+
 class DesiTile(DesiData):
     """Reads the spectra from DESI using tile mode and formats its data as a
     list of Forest instances.
@@ -78,10 +82,13 @@ def __init__(self, config):
         # load variables from config
         self.use_all = None
         self.use_single_nights = None
-        self.__parse_config(config)
+        self.use_non_coadded_spectra = None
 
+        self.__parse_config(config)
+        #init of DesiData needs to come last, as it contains the actual data reading and thus needs all config
         super().__init__(config)
 
+
     def __parse_config(self, config):
         """Parse the configuration options
 
@@ -100,7 +107,15 @@ def __parse_config(self, config):
 
         self.use_single_nights = config.getboolean("use single nights")
         if self.use_single_nights is None:
-            raise DataError("Missing argument 'use single nights' required by DesiTile")
+            raise DataError(
+                "Missing argument 'use single nights' required by DesiTile")
+
+        self.use_non_coadded_spectra = config.getboolean(
+            "use non-coadded spectra")
+        if self.use_non_coadded_spectra is None:
+            raise DataError(
+                "Missing argument 'use non-coadded spectra' required by DesiTile"
+            )
 
     def read_data(self):
         """Read the spectra and formats its data as Forest instances.
@@ -118,8 +133,8 @@ def read_data(self):
         DataError if the analysis type is PK 1D and resolution data is not present
         DataError if no quasars were found
         """
-        if np.any((self.catalogue['TILEID'] < 60000) &
-                  (self.catalogue['TILEID'] >= 1000)):
+        if np.any((self.catalogue['TILEID'] < 60000)
+                  & (self.catalogue['TILEID'] >= 1000)):
             is_sv = False
         else:
             is_sv = True
@@ -127,9 +142,13 @@ def read_data(self):
         forests_by_targetid = {}
         num_data = 0
 
+        coadd_name = "spectra" if self.use_non_coadded_spectra else "coadd"
+
         if self.use_single_nights or "cumulative" in self.input_directory:
-            files_in = sorted(glob.glob(os.path.join(self.input_directory, "**/coadd-*.fits"),
-                              recursive=True))
+            files_in = sorted(
+                glob.glob(os.path.join(self.input_directory,
+                                       f"**/{coadd_name}-*.fits"),
+                          recursive=True))
 
             if "cumulative" in self.input_directory:
                 petal_tile_night = [
@@ -143,15 +162,21 @@ def read_data(self):
                 ]
         else:
             if self.use_all:
-                files_in = sorted(glob.glob(os.path.join(self.input_directory, "**/all/**/coadd-*.fits"),
-                             recursive=True))
+                files_in = sorted(
+                    glob.glob(os.path.join(self.input_directory,
+                                           f"**/all/**/{coadd_name}-*.fits"),
+                              recursive=True))
             else:
-                files_in = sorted(glob.glob(os.path.join(self.input_directory, "**/deep/**/coadd-*.fits"),
-                             recursive=True))
+                files_in = sorted(
+                    glob.glob(os.path.join(self.input_directory,
+                                           f"**/deep/**/{coadd_name}-*.fits"),
+                              recursive=True))
+
             petal_tile = [
                 f"{entry['PETAL_LOC']}-{entry['TILEID']}"
                 for entry in self.catalogue
             ]
+
         # this uniqueness check is to ensure each petal/tile/night combination
         # only appears once in the filelist
         petal_tile_night_unique = np.unique(petal_tile_night)
@@ -169,7 +194,8 @@ def read_data(self):
             try:
                 hdul = fitsio.FITS(filename)
             except IOError:
-                self.logger.warning(f"Error reading file {filename}. Ignoring file")
+                self.logger.warning(
+                    f"Error reading file {filename}. Ignoring file")
                 continue
 
             fibermap = hdul['FIBERMAP'].read()
@@ -202,6 +228,7 @@ def read_data(self):
             petal_spec = fibermap['PETAL_LOC'][0]
 
             targetid_spec = fibermap['TARGETID']
+            no_scores_available = False
 
             spectrographs_data = {}
             for color in colors:
@@ -212,6 +239,17 @@ def read_data(self):
                     spec['IVAR'] = (hdul[f'{color}_IVAR'].read() *
                                     (hdul[f'{color}_MASK'].read() == 0))
                     if self.analysis_type == "PK 1D":
+                        if self.use_non_coadded_spectra and "SCORES" in hdul:
+                            spec['TEFF_LYA'] = 11.80090901380597 * hdul[
+                                'SCORES'][f'TSNR2_LYA_{color}'].read()
+                        else:
+                            spec['TEFF_LYA'] = np.ones(spec["FLUX"].shape[0])
+                            if self.use_non_coadded_spectra:
+                                self.logger.warning(
+                                    "SCORES are missing, Teff information (and thus DIFF) will be garbage"
+                                )
+                            no_scores_available = True
+
                         if f"{color}_RESOLUTION" in hdul:
                             spec["RESO"] = hdul[f"{color}_RESOLUTION"].read()
                         else:
@@ -256,9 +294,15 @@ def read_data(self):
                     w_t = w_t[0]
 
                 for spec in spectrographs_data.values():
-                    ivar = spec['IVAR'][w_t].copy()
-                    flux = spec['FLUX'][w_t].copy()
-
+                    if self.use_non_coadded_spectra:
+                        ivar = np.atleast_2d(spec['IVAR'][w_t])
+                        ivar_coadded_flux = np.atleast_2d(
+                            ivar * spec['FLUX'][w_t]).sum(axis=0)
+                        ivar = ivar.sum(axis=0)
+                        flux = (ivar_coadded_flux / ivar)
+                    else:
+                        flux = spec['FLUX'][w_t].copy()
+                        ivar = spec['IVAR'][w_t].copy()
                     args = {
                         "flux": flux,
                         "ivar": ivar,
@@ -281,18 +325,28 @@ def read_data(self):
                     if self.analysis_type == "BAO 3D":
                         forest = DesiForest(**args)
                     elif self.analysis_type == "PK 1D":
-                        reso_sum = spec['RESO'][w_t].copy()
-                        reso_in_km_per_s = np.real(
+                        if self.use_non_coadded_spectra and not no_scores_available:
+                            exposures_diff = exp_diff_desi(spec, w_t)
+                        else:
+                            exposures_diff = np.zeros(spec['WAVELENGTH'].shape)
+                        if len(spec['RESO'][w_t].shape) < 3:
+                            reso_sum = spec['RESO'][w_t].copy()
+                        else:
+                            reso_sum = spec['RESO'][w_t].sum(axis=0)
+                        reso_in_pix, reso_in_km_per_s = np.real(
                             spectral_resolution_desi(reso_sum,
                                                      spec['WAVELENGTH']))
-                        exposures_diff = np.zeros(spec['log_lambda'].shape)
 
                         args["exposures_diff"] = exposures_diff
                         args["reso"] = reso_in_km_per_s
+                        args["resolution_matrix"] = reso_sum
+                        args["reso_pix"] = reso_in_pix
+
                         forest = DesiPk1dForest(**args)
                     else:
-                        raise DataError("Unkown analysis type. Expected 'BAO 3D'"
-                                        f"or 'PK 1D'. Found '{self.analysis_type}'")
+                        raise DataError(
+                            "Unkown analysis type. Expected 'BAO 3D'"
+                            f"or 'PK 1D'. Found '{self.analysis_type}'")
 
                     # rebin arrays
                     # this needs to happen after all arrays are initialized by
@@ -306,7 +360,8 @@ def read_data(self):
                         forests_by_targetid[targetid] = forest
 
                 num_data += 1
-        self.logger.progress("Found {} quasars in input files".format(num_data))
+        self.logger.progress(
+            "Found {} quasars in input files".format(num_data))
 
         if num_data == 0:
             raise DataError("No Quasars found, stopping here")
diff --git a/py/picca/delta_extraction/data_catalogues/sdss_data.py b/py/picca/delta_extraction/data_catalogues/sdss_data.py
index d187e3981..d12625a20 100644
--- a/py/picca/delta_extraction/data_catalogues/sdss_data.py
+++ b/py/picca/delta_extraction/data_catalogues/sdss_data.py
@@ -181,7 +181,8 @@ def read_from_spec(self, catalogue):
                         "mjd": mjd,
                         "fiberid": fiberid,
                         "exposures_diff": exposures_diff,
-                        "reso": reso
+                        "reso": reso,
+                        "reso_pix": wdisp
                     })
             else:
                 raise DataError(f"analysis_type = {self.analysis_type}")
diff --git a/py/picca/delta_extraction/expected_fluxes/dr16_expected_flux.py b/py/picca/delta_extraction/expected_fluxes/dr16_expected_flux.py
index 3255f613b..c48e755ec 100644
--- a/py/picca/delta_extraction/expected_fluxes/dr16_expected_flux.py
+++ b/py/picca/delta_extraction/expected_fluxes/dr16_expected_flux.py
@@ -428,6 +428,7 @@ def __parse_config(self, config):
         if self.use_ivar_as_weight is None:
             raise ExpectedFluxError(
                 "Missing argument 'use ivar as weight' required by Dr16ExpectedFlux")
+        
 
     def compute_continuum(self, forest):
         """Compute the forest continuum.
diff --git a/py/picca/delta_extraction/masks/absorber_mask.py b/py/picca/delta_extraction/masks/absorber_mask.py
index 46d6d4c06..21fe7535a 100644
--- a/py/picca/delta_extraction/masks/absorber_mask.py
+++ b/py/picca/delta_extraction/masks/absorber_mask.py
@@ -15,6 +15,7 @@
 
 accepted_options = ["absorber mask width", "filename"]
 
+
 class AbsorberMask(Mask):
     """Class to mask Absorbers
 
@@ -58,7 +59,8 @@ def __init__(self, config):
 
         los_id_name = config.get("los_id name")
         if los_id_name is None:
-            raise MaskError("Missing argument 'los_id name' required by AbsorberMask")
+            raise MaskError(
+                "Missing argument 'los_id name' required by AbsorberMask")
 
         columns_list = [los_id_name, "LAMBDA_ABS"]
         try:
@@ -81,17 +83,20 @@ def __init__(self, config):
         for los_id in np.unique(cat[los_id_name]):
             w = (los_id == cat[los_id_name])
             self.los_ids[los_id] = list(cat[los_id_name][w])
-        num_absorbers = np.sum([len(los_id) for los_id in self.los_ids.values()])
+        num_absorbers = np.sum(
+            [len(los_id) for los_id in self.los_ids.values()])
 
         self.logger.progress(" In catalog: {} absorbers".format(num_absorbers))
-        self.logger.progress(" In catalog: {} forests have absorbers\n".format(len(self.los_ids)))
+        self.logger.progress(" In catalog: {} forests have absorbers\n".format(
+            len(self.los_ids)))
 
         # setup transmission limit
         # transmissions below this number are masked
         self.absorber_mask_width = config.getfloat("absorber mask width")
         if self.absorber_mask_width is None:
-            raise MaskError("Missing argument 'absorber mask width' required by "
-                            "AbsorbersMask")
+            raise MaskError(
+                "Missing argument 'absorber mask width' required by "
+                "AbsorbersMask")
 
     def apply_mask(self, forest):
         """Applies the mask. The mask is done by removing the affected
@@ -111,16 +116,20 @@ def apply_mask(self, forest):
         elif Forest.wave_solution == "lin":
             log_labmda = np.log10(forest.lambda_)
         else:
-            raise MaskError("Forest.wave_solution must be either 'log' or 'lin'")
+            raise MaskError(
+                "Forest.wave_solution must be either 'log' or 'lin'")
 
         # load DLAs
         if self.los_ids.get(forest.los_id) is not None:
             # find out which pixels to mask
             w = np.ones(log_lambda.size, dtype=bool)
             for lambda_absorber in self.los_ids.get(forest.los_id):
-                w &= (np.fabs(1.e4 * (log_lambda - np.log10(lambda_absorber))) >
-                      self.absorber_mask_width)
+                w &= (np.fabs(1.e4 * (log_lambda - np.log10(lambda_absorber)))
+                      > self.absorber_mask_width)
 
             # do the actual masking
             for param in Forest.mask_fields:
-                setattr(forest, param, getattr(forest, param)[w])
+                if param in ['resolution_matrix']:
+                    setattr(forest, param, getattr(forest, param)[:, w])
+                else:
+                    setattr(forest, param, getattr(forest, param)[w])
\ No newline at end of file
diff --git a/py/picca/delta_extraction/masks/bal_mask.py b/py/picca/delta_extraction/masks/bal_mask.py
index 088fa0653..901d1900d 100644
--- a/py/picca/delta_extraction/masks/bal_mask.py
+++ b/py/picca/delta_extraction/masks/bal_mask.py
@@ -34,6 +34,7 @@
     "lOI": 1039
 }
 
+
 class BalMask(Mask):
     """Class to mask BALs
 
@@ -86,37 +87,36 @@ def __init__(self, config):
 
         los_id_name = config.get("los_id name")
         if los_id_name is None:
-            raise MaskError("Missing argument 'los_id name' required by BalMask")
+            raise MaskError(
+                "Missing argument 'los_id name' required by BalMask")
         elif los_id_name == "THING_ID":
             ext_name = 'BALCAT'
         elif los_id_name == "TARGETID":
             ext_name = 'ZCATALOG'
         else:
-            raise MaskError("Unrecognized los_id name. Expected one of 'THING_ID' "
-                            f" or 'TARGETID'. Found {los_id_name}")
+            raise MaskError(
+                "Unrecognized los_id name. Expected one of 'THING_ID' "
+                f" or 'TARGETID'. Found {los_id_name}")
 
         # setup bal index limit
         self.bal_index_type = config.getfloat("bal index type")
         if self.bal_index_type is None:
-            self.bal_index_type = MaskError("Missing argument 'bal index type' "
-                                            "required by BalMask")
+            self.bal_index_type = MaskError(
+                "Missing argument 'bal index type' "
+                "required by BalMask")
         if self.bal_index_type == "ai":
-            columns_list = [
-                los_id_name, 'VMIN_CIV_450', 'VMAX_CIV_450'
-            ]
+            columns_list = [los_id_name, 'VMIN_CIV_450', 'VMAX_CIV_450']
         elif self.bal_index_type == "bi":
-            columns_list = [
-                los_id_name, 'VMIN_CIV_2000', 'VMAX_CIV_2000'
-            ]
+            columns_list = [los_id_name, 'VMIN_CIV_2000', 'VMAX_CIV_2000']
         else:
-            self.bal_index_type = MaskError("In BalMask, unrecognized value "
-                                            "for 'bal_index_type'. Expected one "
-                                            "of 'ai' or 'bi'. Found "
-                                            f"{self.bal_index_type}")
+            self.bal_index_type = MaskError(
+                "In BalMask, unrecognized value "
+                "for 'bal_index_type'. Expected one "
+                "of 'ai' or 'bi'. Found "
+                f"{self.bal_index_type}")
 
         self.logger.progress(f"Reading BAL catalog from: {filename}")
 
-
         try:
             hdul = fitsio.FITS(filename)
             self.cat = {col: hdul[ext_name][col][:] for col in columns_list}
@@ -154,8 +154,9 @@ def add_bal_rest_frame(self, los_id, los_id_name):
         else:  # AI, the default
             velocity_list = ['VMIN_CIV_450', 'VMAX_CIV_450']
 
-        mask_rest_frame_bal = Table(names=['log_lambda_min', 'log_lambda_max',
-                                           'lambda_min', 'lambda_max'],
+        mask_rest_frame_bal = Table(names=[
+            'log_lambda_min', 'log_lambda_max', 'lambda_min', 'lambda_max'
+        ],
                                     dtype=['f4', 'f4', 'f4', 'f4'])
         min_velocities = []  # list of minimum velocities
         max_velocities = []  # list of maximum velocities
@@ -181,9 +182,12 @@ def add_bal_rest_frame(self, los_id, los_id_name):
             for line in lines.values():
                 log_lambda_min = np.log10(line * (1 - min_vel / SPEED_LIGHT))
                 log_lambda_max = np.log10(line * (1 - max_vel / SPEED_LIGHT))
-                mask_rest_frame_bal.add_row([log_lambda_min, log_lambda_max,
-                                             10**log_lambda_min, 10**log_lambda_max,
-                                             ])
+                mask_rest_frame_bal.add_row([
+                    log_lambda_min,
+                    log_lambda_max,
+                    10**log_lambda_min,
+                    10**log_lambda_max,
+                ])
 
         return mask_rest_frame_bal
 
@@ -207,18 +211,25 @@ def apply_mask(self, forest):
             if Forest.wave_solution == "log":
                 w = np.ones(forest.log_lambda.size, dtype=bool)
                 for mask_range in mask_table:
-                    rest_frame_log_lambda = forest.log_lambda - np.log10(1. + forest.z)
-                    w &= ((rest_frame_log_lambda < mask_range['log_lambda_min']) |
-                          (rest_frame_log_lambda > mask_range['log_lambda_max']))
+                    rest_frame_log_lambda = forest.log_lambda - np.log10(
+                        1. + forest.z)
+                    w &= (
+                        (rest_frame_log_lambda < mask_range['log_lambda_min'])
+                        |
+                        (rest_frame_log_lambda > mask_range['log_lambda_max']))
             elif Forest.wave_solution == "lin":
                 w = np.ones(forest.lambda_.size, dtype=bool)
                 for mask_range in mask_table:
-                    rest_frame_lambda = forest.lambda_/(1. + forest.z)
+                    rest_frame_lambda = forest.lambda_ / (1. + forest.z)
                     w &= ((rest_frame_lambda < mask_range['lambda_min']) |
                           (rest_frame_lambda > mask_range['lambda_max']))
             else:
-                raise MaskError("Forest.wave_solution must be either 'log' or 'lin'")
+                raise MaskError(
+                    "Forest.wave_solution must be either 'log' or 'lin'")
 
             # do the actual masking
             for param in Forest.mask_fields:
-                setattr(forest, param, getattr(forest, param)[w])
+                if param in ['resolution_matrix']:
+                    setattr(forest, param, getattr(forest, param)[:, w])
+                else:
+                    setattr(forest, param, getattr(forest, param)[w])
\ No newline at end of file
diff --git a/py/picca/delta_extraction/masks/dla_mask.py b/py/picca/delta_extraction/masks/dla_mask.py
index 3dca3c942..a7024dc9c 100644
--- a/py/picca/delta_extraction/masks/dla_mask.py
+++ b/py/picca/delta_extraction/masks/dla_mask.py
@@ -22,6 +22,7 @@
 num_points = 10000
 gaussian_dist = np.random.normal(size=num_points) * np.sqrt(2)
 
+
 class DlaMask(Mask):
     """Class to mask DLAs
 
@@ -73,7 +74,8 @@ def __init__(self, config):
 
         los_id_name = config.get("los_id name")
         if los_id_name is None:
-            raise MaskError("Missing argument 'los_id name' required by DlaMask")
+            raise MaskError(
+                "Missing argument 'los_id name' required by DlaMask")
 
         self.logger.progress(f"Reading DLA catalog from: {filename}")
         columns_list = [los_id_name, "Z", "NHI"]
@@ -98,7 +100,8 @@ def __init__(self, config):
         num_dlas = np.sum([len(los_id) for los_id in self.los_ids.values()])
 
         self.logger.progress('In catalog: {} DLAs'.format(num_dlas))
-        self.logger.progress('In catalog: {} forests have a DLA\n'.format(len(self.los_ids)))
+        self.logger.progress('In catalog: {} forests have a DLA\n'.format(
+            len(self.los_ids)))
 
         # setup transmission limit
         # transmissions below this number are masked
@@ -112,7 +115,8 @@ def __init__(self, config):
         if mask_file is not None:
             try:
                 self.mask = Table.read(mask_file,
-                                       names=('type', 'wave_min', 'wave_max', 'frame'),
+                                       names=('type', 'wave_min', 'wave_max',
+                                              'frame'),
                                        format='ascii')
                 self.mask = self.mask['frame'] == 'RF_DLA'
             except (OSError, ValueError):
@@ -139,28 +143,32 @@ def apply_mask(self, forest):
         elif Forest.wave_solution == "lin":
             lambda_ = forest.lambda_
         else:
-            raise MaskError("Forest.wave_solution must be either 'log' or 'lin'")
+            raise MaskError(
+                "Forest.wave_solution must be either 'log' or 'lin'")
 
         # load DLAs
         if self.los_ids.get(forest.los_id) is not None:
             dla_transmission = np.ones(len(lambda_))
             for (z_abs, nhi) in self.los_ids.get(forest.los_id):
-                dla_transmission *= DlaProfile(lambda_, z_abs, nhi).transmission
+                dla_transmission *= DlaProfile(lambda_, z_abs,
+                                               nhi).transmission
 
             # find out which pixels to mask
             w = dla_transmission > self.dla_mask_limit
             if len(self.mask) > 0:
                 for mask_range in self.mask:
                     for (z_abs, nhi) in self.los_ids.get(forest.los_id):
-                        w &= ((lambda_ / (1. + z_abs) <
-                               mask_range['wave_min']) |
-                              (lambda_ / (1. + z_abs) >
-                               mask_range['wave_max']))
+                        w &= ((lambda_ / (1. + z_abs) < mask_range['wave_min'])
+                              | (lambda_ /
+                                 (1. + z_abs) > mask_range['wave_max']))
 
             # do the actual masking
             forest.transmission_correction *= dla_transmission
             for param in Forest.mask_fields:
-                setattr(forest, param, getattr(forest, param)[w])
+                if param in ['resolution_matrix']:
+                    setattr(forest, param, getattr(forest, param)[:, w])
+                else:
+                    setattr(forest, param, getattr(forest, param)[w])
 
 
 class DlaProfile:
@@ -206,10 +214,8 @@ def __init__(self, lambda_, z_abs, nhi):
         self.z_abs = z_abs
         self.nhi = nhi
 
-        self.transmission = self.profile_lya_absorption(lambda_,
-                                                        z_abs, nhi)
-        self.transmission *= self.profile_lyb_absorption(
-            lambda_, z_abs, nhi)
+        self.transmission = self.profile_lya_absorption(lambda_, z_abs, nhi)
+        self.transmission *= self.profile_lyb_absorption(lambda_, z_abs, nhi)
 
     @staticmethod
     def profile_lya_absorption(lambda_, z_abs, nhi):
diff --git a/py/picca/delta_extraction/masks/lines_mask.py b/py/picca/delta_extraction/masks/lines_mask.py
index e0eb8dd2e..db8d8f09b 100644
--- a/py/picca/delta_extraction/masks/lines_mask.py
+++ b/py/picca/delta_extraction/masks/lines_mask.py
@@ -12,6 +12,7 @@
 
 accepted_options = ["filename"]
 
+
 class LinesMask(Mask):
     """Class to mask (sky) lines
 
@@ -44,7 +45,8 @@ def __init__(self, config):
 
         mask_file = config.get("filename")
         if mask_file is None:
-            raise MaskError("Missing argument 'filename' required by LinesMask")
+            raise MaskError(
+                "Missing argument 'filename' required by LinesMask")
         try:
             mask = Table.read(mask_file,
                               names=('type', 'wave_min', 'wave_max', 'frame'),
@@ -80,7 +82,8 @@ def apply_mask(self, forest):
             for mask_range in self.mask_obs_frame:
                 w &= ((forest.log_lambda < mask_range['log_wave_min']) |
                       (forest.log_lambda > mask_range['log_wave_max']))
-            log_lambda_rest_frame = forest.log_lambda - np.log10(1.0 + forest.z)
+            log_lambda_rest_frame = forest.log_lambda - np.log10(1.0 +
+                                                                 forest.z)
             for mask_range in self.mask_rest_frame:
                 w &= ((log_lambda_rest_frame < mask_range['log_wave_min']) |
                       (log_lambda_rest_frame > mask_range['log_wave_max']))
@@ -89,7 +92,7 @@ def apply_mask(self, forest):
             for mask_range in self.mask_obs_frame:
                 w &= ((forest.lambda_ < mask_range['wave_min']) |
                       (forest.lambda_ > mask_range['wave_max']))
-            lambda_rest_frame = forest.lambda_/(1.0 + forest.z)
+            lambda_rest_frame = forest.lambda_ / (1.0 + forest.z)
             for mask_range in self.mask_rest_frame:
                 w &= ((lambda_rest_frame < mask_range['wave_min']) |
                       (lambda_rest_frame > mask_range['wave_max']))
@@ -100,4 +103,7 @@ def apply_mask(self, forest):
 
         # do the actual masking
         for param in Forest.mask_fields:
-            setattr(forest, param, getattr(forest, param)[w])
+            if param in ['resolution_matrix']:
+                setattr(forest, param, getattr(forest, param)[:, w])
+            else:
+                setattr(forest, param, getattr(forest, param)[w])
\ No newline at end of file
diff --git a/py/picca/delta_extraction/utils_pk1d.py b/py/picca/delta_extraction/utils_pk1d.py
index d777639f0..1eb27bb0f 100644
--- a/py/picca/delta_extraction/utils_pk1d.py
+++ b/py/picca/delta_extraction/utils_pk1d.py
@@ -16,6 +16,7 @@
 # create logger
 module_logger = logging.getLogger(__name__)
 
+
 def exp_diff(hdul, log_lambda):
     """Compute the difference between exposures.
 
@@ -59,7 +60,8 @@ def exp_diff(hdul, log_lambda):
 
             # exclude masks 25 (COMBINEREJ), 23 (BRIGHTSKY)?
             rebin_ivar_exp = np.bincount(log_lambda_bins,
-                                         weights=ivar_exp * (mask & 2**25 == 0))
+                                         weights=ivar_exp *
+                                         (mask & 2**25 == 0))
             rebin_flux_exp = np.bincount(log_lambda_bins,
                                          weights=(ivar_exp * flux_exp *
                                                   (mask & 2**25 == 0)))
@@ -68,8 +70,10 @@ def exp_diff(hdul, log_lambda):
                 flux_total_odd[:len(rebin_ivar_exp) - 1] += rebin_flux_exp[:-1]
                 ivar_total_odd[:len(rebin_ivar_exp) - 1] += rebin_ivar_exp[:-1]
             else:
-                flux_total_even[:len(rebin_ivar_exp) - 1] += rebin_flux_exp[:-1]
-                ivar_total_even[:len(rebin_ivar_exp) - 1] += rebin_ivar_exp[:-1]
+                flux_total_even[:len(rebin_ivar_exp) -
+                                1] += rebin_flux_exp[:-1]
+                ivar_total_even[:len(rebin_ivar_exp) -
+                                1] += rebin_ivar_exp[:-1]
 
     w = ivar_total_odd > 0
     flux_total_odd[w] /= ivar_total_odd[w]
@@ -87,6 +91,77 @@ def exp_diff(hdul, log_lambda):
     return exposures_diff
 
 
+def exp_diff_desi(spec_dict, mask_targetid):
+    """Computes the difference between exposures.
+
+    More precisely computes de semidifference between two customized coadded
+    spectra obtained from weighted averages of the even-number exposures, for
+    the first spectrum, and of the odd-number exposures, for the second one
+    (see section 3.2 of Chabanier et al. 2019).
+
+    Args:
+        spec_dict: dict
+            spec dictionary from desi_healpix/tile
+        mask targetid: array of ints
+            Targetids to select for calculating the exp differences
+
+    Returns:
+        The difference between exposures
+    """
+    teff_lya = np.atleast_1d(spec_dict["TEFF_LYA"][mask_targetid])
+    argsort = np.flip(np.argsort(teff_lya))
+    flux = np.atleast_2d(spec_dict["FLUX"][mask_targetid])[argsort, :]
+    ivar = np.atleast_2d(spec_dict["IVAR"][mask_targetid])[argsort, :]
+    teff_lya = teff_lya[argsort]
+
+    num_exp = len(flux)
+    if (num_exp < 2):
+        module_logger.debug("Not enough exposures for diff, Spectra rejected")
+        return None
+    elif (num_exp > 100):
+        module_logger.debug("More than 100 exposures, potentially wrong file type and using wavelength axis here, skipping?")
+        return None
+    flux_total_odd = np.zeros(flux.shape[1])
+    ivar_total_odd = np.zeros(flux.shape[1])
+    flux_total_even = np.zeros(flux.shape[1])
+    ivar_total_even = np.zeros(flux.shape[1])
+    teff_even = 0
+    teff_odd = 0
+    teff_total = np.sum(teff_lya)
+    teff_last = teff_lya[-1]
+    for index_exp in range(2 * (num_exp // 2)):
+        flexp = flux[index_exp]
+        ivexp = ivar[index_exp]
+        teff_lya_exp = teff_lya[index_exp]
+        if index_exp % 2 == 1:
+            flux_total_odd += flexp * ivexp
+            ivar_total_odd += ivexp
+            teff_odd += teff_lya_exp
+        else:
+            flux_total_even += flexp * ivexp
+            ivar_total_even += ivexp
+            teff_even += teff_lya_exp
+
+    w = ivar_total_odd > 0
+    flux_total_odd[w] /= ivar_total_odd[w]
+    w = ivar_total_even > 0
+    flux_total_even[w] /= ivar_total_even[w]
+
+    alpha = 1
+    if (num_exp % 2 == 1):
+        n_even = (num_exp - 1) // 2
+        alpha_N_old = np.sqrt(4. * n_even * (num_exp - n_even)) / num_exp
+        # alpha_N = np.sqrt(4.*t_even*(t_exp-t_even))/t_exp
+        alpha_C_new = np.sqrt((teff_total - teff_last) / teff_total)
+        alpha_N_new = np.sqrt(
+            (teff_total - teff_last) * (teff_total + teff_last)) / teff_total
+        alpha = alpha_N_new
+    diff = 0.5 * (flux_total_even - flux_total_odd) * alpha
+
+    return diff
+
+
+
 def spectral_resolution(wdisp,
                         with_correction=False,
                         fiberid=None,
@@ -126,11 +201,11 @@ def spectral_resolution(wdisp,
         # fiberids greater than 500 corresponds to the second spectrograph
         fiberid = fiberid % 500
         if fiberid < 100:
-            correction = (1. + (correction - 1) * .25 + (correction - 1) * .75 *
-                          (fiberid) / 100.)
+            correction = (1. + (correction - 1) * .25 +
+                          (correction - 1) * .75 * (fiberid) / 100.)
         elif fiberid > 400:
-            correction = (1. + (correction - 1) * .25 + (correction - 1) * .75 *
-                          (500 - fiberid) / 100.)
+            correction = (1. + (correction - 1) * .25 +
+                          (correction - 1) * .75 * (500 - fiberid) / 100.)
 
         # apply the correction
         reso *= correction
@@ -154,14 +229,29 @@ def spectral_resolution_desi(reso_matrix, lambda_):
     reso_in_km_per_s: array
     The spectral resolution
     """
-    delta_lambda = ((lambda_[-1] - lambda_[0]) /
-                    float(len(lambda_) - 1))
+    delta_log_lambda = np.diff(np.log10(lambda_))
+    #note that this would be the same result as before (except for the missing bug) in
+    #case of log-uniform binning, but for linear binning pixel size chenges wrt lambda
+    delta_log_lambda = np.append(
+        delta_log_lambda,
+        [delta_log_lambda[-1] + (delta_log_lambda[-1] - delta_log_lambda[-2])])
     reso = np.clip(reso_matrix, 1.0e-6, 1.0e6)
-    rms_in_pixel = (np.sqrt(1.0 / 2.0 / np.log(
-        reso[len(reso) // 2][:] / reso[len(reso) // 2 - 1][:])) + np.sqrt(
-            4.0 / 2.0 / np.log(
-                reso[len(reso) // 2][:] / reso[len(reso) // 2 - 2][:]))) / 2.0
-
-    reso_in_km_per_s = (rms_in_pixel * SPEED_LIGHT * delta_lambda)
-
-    return reso_in_km_per_s
+    #assume reso = A*exp(-(x-central_pixel_pos)**2 / 2 / sigma**2)
+    #=> sigma = sqrt((x-central_pixel_pos)/2)**2 / log(A/reso)
+    #   A = reso(central_pixel_pos)  
+    # the following averages over estimates for four symmetric values of x
+    rms_in_pixel = (
+        (np.sqrt(1.0 / 2.0 / np.log(
+            reso[reso.shape[0] // 2, :] / reso[reso.shape[0] // 2 - 1, :])) +
+         np.sqrt(4.0 / 2.0 / np.log(
+             reso[reso.shape[0] // 2, :] / reso[reso.shape[0] // 2 - 2, :])) +
+         np.sqrt(1.0 / 2.0 / np.log(
+             reso[reso.shape[0] // 2, :] / reso[reso.shape[0] // 2 + 1, :])) +
+         np.sqrt(4.0 / 2.0 / np.log(
+             reso[reso.shape[0] // 2, :] / reso[reso.shape[0] // 2 + 2, :])))
+        / 4.0) #this is rms
+
+    reso_in_km_per_s = (rms_in_pixel * SPEED_LIGHT * delta_log_lambda *
+                        np.log(10.0))   #this is FWHM
+
+    return rms_in_pixel, reso_in_km_per_s
\ No newline at end of file
diff --git a/py/picca/io.py b/py/picca/io.py
index d6fb04ce6..5a25859a0 100644
--- a/py/picca/io.py
+++ b/py/picca/io.py
@@ -1328,7 +1328,7 @@ def read_from_minisv_desi(in_dir, catalog, pk1d=None, usesinglenights=False, use
 
                 if pk1d is not None:
                     reso_sum = spec['RESO'][w_t].copy()
-                    reso_in_km_per_s = np.real(
+                    reso_in_pixel, reso_in_km_per_s = np.real(
                         spectral_resolution_desi(reso_sum, spec['log_lambda']))
                     exposures_diff = np.zeros(spec['log_lambda'].shape)
                 else:
diff --git a/py/picca/pk1d.py b/py/picca/pk1d.py
index c72a9bfde..9c0c9fc28 100644
--- a/py/picca/pk1d.py
+++ b/py/picca/pk1d.py
@@ -11,29 +11,31 @@
 See the respective docstrings for more details
 """
 import numpy as np
-from scipy.fftpack import fft
+from numpy.fft import rfft, rfftfreq
 
 from . import constants
 from .utils import userprint
 
 
 def split_forest(num_parts,
-                 delta_log_lambda,
-                 log_lambda,
+                 pixel_step,
+                 lambda_or_log_lambda,
                  delta,
                  exposures_diff,
                  ivar,
                  first_pixel_index,
-                 abs_igm="LYA"):
+                 abs_igm="LYA",
+                 reso_matrix=None,
+                 linear_binning=False):
     """Splits the forest in n parts
 
     Args:
         num_parts: int
             Number of parts
-        delta_log_lambda: float
-            Variation of the logarithm of the wavelength between two pixels
-        log_lambda: array of float
-            Logarith of the wavelength (in Angs)
+        pixel_step: float
+            Variation of the wavelength (or log(wavelength)) between two pixels
+        lambda_or_log_lambda: array of float
+            Wavelength (in Angs) (or its log, but needs to be consistent with delta_lambda_or_log_lambda)
         delta: array of float
             Mean transmission fluctuation (delta field)
         exposures_diff: array of float
@@ -47,11 +49,14 @@ def split_forest(num_parts,
         abs_igm: string - default: "LYA"
             Name of the absorption in picca.constants defining the
             redshift of the forest pixels
+        reso_matrix: 2d-array of floats
+            The resolution matrix used for corrections
+        linear_binning: assume linear wavelength binning, lambda_or_log_lambda vectors will be actual lambda (else they will be log_lambda)
 
     Returns:
         The following variables:
             mean_z_array: Array with the mean redshift the parts of the forest
-            log_lambda_array: Array with logarith of the wavelength for the
+            lambda_or_log_lambda_array: Array with logarith of the wavelength for the
                 parts of the forest
             delta_array: Array with the deltas for the parts of the forest
             exposures_diff_array: Array with the exposures_diff for the parts of
@@ -59,53 +64,77 @@ def split_forest(num_parts,
             ivar_array: Array with the ivar for the parts of the forest
 
     """
-    log_lambda_limit = [log_lambda[first_pixel_index]]
-    num_bins = (len(log_lambda) - first_pixel_index) // num_parts
+    lambda_or_log_lambda_limit = [lambda_or_log_lambda[first_pixel_index]]
+    num_bins = (len(lambda_or_log_lambda) - first_pixel_index) // num_parts
 
     mean_z_array = []
-    log_lambda_array = []
+    lambda_or_log_lambda_array = []
     delta_array = []
     exposures_diff_array = []
     ivar_array = []
+    if reso_matrix is not None:
+        reso_matrix_array = []
 
     for index in range(1, num_parts):
-        log_lambda_limit.append(log_lambda[num_bins * index +
-                                           first_pixel_index])
+        lambda_or_log_lambda_limit.append(
+            lambda_or_log_lambda[num_bins * index + first_pixel_index])
 
-    log_lambda_limit.append(log_lambda[len(log_lambda) - 1] +
-                            0.1 * delta_log_lambda)
+    lambda_or_log_lambda_limit.append(
+        lambda_or_log_lambda[len(lambda_or_log_lambda) - 1] +
+        0.1 * pixel_step)
 
     for index in range(num_parts):
-        selection = ((log_lambda >= log_lambda_limit[index]) &
-                     (log_lambda < log_lambda_limit[index + 1]))
+        selection = (
+            (lambda_or_log_lambda >= lambda_or_log_lambda_limit[index]) &
+            (lambda_or_log_lambda < lambda_or_log_lambda_limit[index + 1]))
 
-        log_lambda_part = log_lambda[selection].copy()
+        lambda_or_log_lambda_part = lambda_or_log_lambda[selection].copy()
         lambda_abs_igm = constants.ABSORBER_IGM[abs_igm]
-        mean_z = (np.power(10., log_lambda_part[len(log_lambda_part) - 1]) +
-                  np.power(10., log_lambda_part[0])) / 2. / lambda_abs_igm - 1.0
+
+        if linear_binning:
+            mean_z = np.mean(lambda_or_log_lambda_part) / lambda_abs_igm - 1.0
+        else:
+            mean_z = np.mean(10**
+                             lambda_or_log_lambda_part) / lambda_abs_igm - 1.0
+
+        if reso_matrix is not None:
+            reso_matrix_part = reso_matrix[:, selection].copy()
 
         mean_z_array.append(mean_z)
-        log_lambda_array.append(log_lambda_part)
+        lambda_or_log_lambda_array.append(lambda_or_log_lambda_part)
         delta_array.append(delta[selection].copy())
         exposures_diff_array.append(exposures_diff[selection].copy())
         ivar_array.append(ivar[selection].copy())
+        if reso_matrix is not None:
+            reso_matrix_array.append(reso_matrix_part)
 
-    return (mean_z_array, log_lambda_array, delta_array, exposures_diff_array,
-            ivar_array)
+    out = [
+        mean_z_array, lambda_or_log_lambda_array, delta_array,
+        exposures_diff_array, ivar_array
+    ]
+    if reso_matrix is not None:
+        out.append(reso_matrix_array)
+    return out
 
 
-def rebin_diff_noise(delta_log_lambda, log_lambda, exposures_diff):
+def rebin_diff_noise(pixel_step, lambda_or_log_lambda,
+                     exposures_diff):
     """Rebin the semidifference between two customized coadded spectra to
     construct the noise array
 
+    Note that inputs can be either linear or log-lambda spaced units (but 
+    pixel_step and lambda_or_log_lambda need the same unit)
+
     The rebinning is done by combining 3 of the original pixels into analysis
     pixels.
 
     Args:
-        delta_log_lambda: float
-            Variation of the logarithm of the wavelength between two pixels
-        log_lambda: array of floats
+        pixel_step: float
+            Variation of the logarithm of the wavelength between two pixels 
+            for linear binnings this would need to be the wavelength difference
+        lambda_or_log_lambda: array of floats
             Array containing the logarithm of the wavelengths (in Angs)
+            for linear binnings this would need to be just wavelength
         exposures_diff: array of floats
             Semidifference between two customized coadded spectra obtained from
             weighted averages of the even-number exposures, for the first
@@ -118,11 +147,11 @@ def rebin_diff_noise(delta_log_lambda, log_lambda, exposures_diff):
     if exposures_diff.size < rebin:
         userprint("Warning: exposures_diff.size too small for rebin")
         return exposures_diff
-    rebin_delta_log_lambda = rebin * delta_log_lambda
+    rebin_delta_lambda_or_log_lambda = rebin * pixel_step
 
     # rebin not mixing pixels separated by masks
-    bins = np.floor((log_lambda - log_lambda.min()) / rebin_delta_log_lambda +
-                    0.5).astype(int)
+    bins = np.floor((lambda_or_log_lambda - lambda_or_log_lambda.min()) /
+                    rebin_delta_lambda_or_log_lambda + 0.5).astype(int)
 
     rebin_exposure_diff = np.bincount(bins.astype(int), weights=exposures_diff)
     rebin_counts = np.bincount(bins.astype(int))
@@ -145,14 +174,18 @@ def rebin_diff_noise(delta_log_lambda, log_lambda, exposures_diff):
     return noise
 
 
-def fill_masked_pixels(delta_log_lambda, log_lambda, delta, exposures_diff,
-                       ivar, no_apply_filling):
+def fill_masked_pixels(delta_lambda_or_log_lambda, lambda_or_log_lambda, delta,
+                       exposures_diff, ivar, no_apply_filling):
     """Fills the masked pixels with zeros
 
+    Note that inputs can be either linear or log-lambda spaced units (but 
+    delta_lambda_or_log_lambda and lambda_or_log_lambda need the same unit);
+    spectrum needs to be uniformly binned in that unit
+
     Args:
-        delta_log_lambda: float
+        delta_lambda_or_log_lambda: float
             Variation of the logarithm of the wavelength between two pixels
-        log_lambda: array of floats
+        lambda_or_log_lambda: array of floats
             Array containing the logarithm of the wavelengths (in Angs)
         delta: array of floats
             Mean transmission fluctuation (delta field)
@@ -167,7 +200,7 @@ def fill_masked_pixels(delta_log_lambda, log_lambda, delta, exposures_diff,
 
     Returns:
         The following variables:
-            log_lambda_new: Array containing the logarithm of the wavelengths
+            lambda_or_log_lambda_new: Array containing the logarithm of the wavelengths
                 (in Angs)
             delta_new: Mean transmission fluctuation (delta field)
             exposures_diff_new: Semidifference between two customized coadded
@@ -178,22 +211,23 @@ def fill_masked_pixels(delta_log_lambda, log_lambda, delta, exposures_diff,
             num_masked_pixels: Number of masked pixels
     """
     if no_apply_filling:
-        return log_lambda, delta, exposures_diff, ivar, 0
+        return lambda_or_log_lambda, delta, exposures_diff, ivar, 0
 
-    log_lambda_index = log_lambda.copy()
-    log_lambda_index -= log_lambda[0]
-    log_lambda_index /= delta_log_lambda
-    log_lambda_index += 0.5
-    log_lambda_index = np.array(log_lambda_index, dtype=int)
-    index_all = range(log_lambda_index[-1] + 1)
-    index_ok = np.in1d(index_all, log_lambda_index)
+    lambda_or_log_lambda_index = lambda_or_log_lambda.copy()
+    lambda_or_log_lambda_index -= lambda_or_log_lambda[0]
+    lambda_or_log_lambda_index /= delta_lambda_or_log_lambda
+    lambda_or_log_lambda_index += 0.5
+    lambda_or_log_lambda_index = np.array(lambda_or_log_lambda_index,
+                                          dtype=int)
+    index_all = range(lambda_or_log_lambda_index[-1] + 1)
+    index_ok = np.in1d(index_all, lambda_or_log_lambda_index)
 
     delta_new = np.zeros(len(index_all))
     delta_new[index_ok] = delta
 
-    log_lambda_new = np.array(index_all, dtype=float)
-    log_lambda_new *= delta_log_lambda
-    log_lambda_new += log_lambda[0]
+    lambda_or_log_lambda_new = np.array(index_all, dtype=float)
+    lambda_or_log_lambda_new *= delta_lambda_or_log_lambda
+    lambda_or_log_lambda_new += lambda_or_log_lambda[0]
 
     exposures_diff_new = np.zeros(len(index_all))
     exposures_diff_new[index_ok] = exposures_diff
@@ -201,51 +235,59 @@ def fill_masked_pixels(delta_log_lambda, log_lambda, delta, exposures_diff,
     ivar_new = np.zeros(len(index_all), dtype=float)
     ivar_new[index_ok] = ivar
 
-    num_masked_pixels = len(index_all) - len(log_lambda_index)
+    num_masked_pixels = len(index_all) - len(lambda_or_log_lambda_index)
 
-    return (log_lambda_new, delta_new, exposures_diff_new, ivar_new,
+    return (lambda_or_log_lambda_new, delta_new, exposures_diff_new, ivar_new,
             num_masked_pixels)
 
 
-def compute_pk_raw(delta_log_lambda, delta):
+def compute_pk_raw(delta_lambda_or_log_lambda, delta, linear_binning=False):
     """Computes the raw power spectrum
 
     Args:
-        delta_log_lambda: float
-            Variation of the logarithm of the wavelength between two pixels
+        delta_lambda_or_log_lambda: float
+            Variation of (the logarithm of) the wavelength between two pixels
         delta: array of floats
             Mean transmission fluctuation (delta field)
+        linear_binning: if set then inputs need to be in AA, outputs will be 1/AA
+                        else inputs will be in log(AA) and outputs in s/km
 
     Returns:
         The following variables
             k: the Fourier modes the Power Spectrum is measured on
             pk: the Power Spectrum
     """
-    # spectral length in km/s
-    length_lambda = (delta_log_lambda * constants.SPEED_LIGHT * np.log(10.) *
-                     len(delta))
+    if linear_binning:  # spectral length in AA
+        length_lambda = (delta_lambda_or_log_lambda * len(delta))
+    else:  # spectral length in km/s
+        length_lambda = (delta_lambda_or_log_lambda * constants.SPEED_LIGHT *
+                         np.log(10.) * len(delta))
 
     # make 1D FFT
     num_pixels = len(delta)
-    num_bins_fft = num_pixels // 2 + 1
-    fft_delta = fft(delta)
+    fft_delta = rfft(delta)
 
     # compute power spectrum
-    fft_delta = fft_delta[:num_bins_fft]
-    pk = (fft_delta.real**2 + fft_delta.imag**2) * length_lambda / num_pixels**2
-    k = np.arange(num_bins_fft, dtype=float) * 2 * np.pi / length_lambda
+    pk = (fft_delta.real**2 +
+          fft_delta.imag**2) * length_lambda / num_pixels**2
+    k = 2 * np.pi * rfftfreq(num_pixels, length_lambda / num_pixels)
 
     return k, pk
 
 
-def compute_pk_noise(delta_log_lambda, ivar, exposures_diff, run_noise):
+def compute_pk_noise(delta_lambda_or_log_lambda,
+                     ivar,
+                     exposures_diff,
+                     run_noise,
+                     num_noise_exposures=10,
+                     linear_binning=False):
     """Computes the noise power spectrum
 
     Two noise power spectrum are computed: one using the pipeline noise and
     another one using the noise derived from exposures_diff
 
     Args:
-        delta_log_lambda: float
+        delta_lambda_or_log_lambda: float
             Variation of the logarithm of the wavelength between two pixels
         ivar: array of floats
             Array containing the inverse variance
@@ -256,6 +298,8 @@ def compute_pk_noise(delta_log_lambda, ivar, exposures_diff, run_noise):
         run_noise: boolean
             If False the noise power spectrum using the pipeline noise is not
             computed and an array filled with zeros is returned instead
+        num_noise_exposures: int
+            Number of exposures to average for noise power estimate
 
     Returns:
         The following variables
@@ -266,7 +310,6 @@ def compute_pk_noise(delta_log_lambda, ivar, exposures_diff, run_noise):
     num_pixels = len(ivar)
     num_bins_fft = num_pixels // 2 + 1
 
-    num_noise_exposures = 10
     pk_noise = np.zeros(num_bins_fft)
     error = np.zeros(num_pixels)
     w = ivar > 0
@@ -276,12 +319,14 @@ def compute_pk_noise(delta_log_lambda, ivar, exposures_diff, run_noise):
         for _ in range(num_noise_exposures):
             delta_exp = np.zeros(num_pixels)
             delta_exp[w] = np.random.normal(0., error[w])
-            _, pk_exp = compute_pk_raw(delta_log_lambda, delta_exp)
+            _, pk_exp = compute_pk_raw(delta_lambda_or_log_lambda,
+                                       delta_exp,
+                                       linear_binning=linear_binning)
             pk_noise += pk_exp
 
         pk_noise /= float(num_noise_exposures)
 
-    _, pk_diff = compute_pk_raw(delta_log_lambda, exposures_diff)
+    _, pk_diff = compute_pk_raw(delta_lambda_or_log_lambda, exposures_diff, linear_binning=linear_binning)
 
     return pk_noise, pk_diff
 
@@ -292,7 +337,7 @@ def compute_correction_reso(delta_pixel, mean_reso, k):
     Args:
         delta_pixel: float
             Variation of the logarithm of the wavelength between two pixels
-            (in km/s)
+            (in km/s or Ang depending on the units of k submitted)
         mean_reso: float
             Mean resolution of the forest
         k: array of floats
@@ -304,12 +349,54 @@ def compute_correction_reso(delta_pixel, mean_reso, k):
     num_bins_fft = len(k)
     correction = np.ones(num_bins_fft)
 
-    sinc = np.ones(num_bins_fft)
-    sinc[k > 0.] = (np.sin(k[k > 0.] * delta_pixel / 2.0) /
-                    (k[k > 0.] * delta_pixel / 2.0))**2
+    pixelization_factor = np.sinc(k * delta_pixel / (2 * np.pi))**2
 
     correction *= np.exp(-(k * mean_reso)**2)
-    correction *= sinc
+    correction *= pixelization_factor
+    return correction
+
+
+def compute_correction_reso_matrix(reso_matrix, k, delta_pixel, num_pixel):
+    """Computes the resolution correction based on the resolution matrix using linear binning
+
+    Args:
+        delta_pixel: float
+            Variation of the logarithm of the wavelength between two pixels
+            (in km/s or Ang depending on the units of k submitted)
+        num_pixel: int
+            Length  of the spectrum in pixels
+        mean_reso: float
+            Mean resolution of the forest
+        k: array of floats
+            Fourier modes
+
+    Returns:
+        The resolution correction
+    """
+    #this allows either computing the power for each pixel seperately or for the mean
+    reso_matrix = np.atleast_2d(reso_matrix)
+
+    W2arr = []
+    #first compute the power in the resmat for each pixel, then average
+    for resmat in reso_matrix:
+        r = np.append(resmat, np.zeros(num_pixel - resmat.size))
+        k_resmat, W2 = compute_pk_raw(delta_pixel, r, linear_binning=True)
+        try:
+            assert np.all(k_resmat == k)
+        except AssertionError:
+            raise ("for some reason the resolution matrix correction has "
+                   "different k scaling than the pk")
+        W2 /= W2[0]
+        W2arr.append(W2)
+
+    Wres2 = np.mean(W2arr, axis=0)
+    pixelization_factor = np.sinc(k * delta_pixel / (2 * np.pi))**2
+
+
+    #the following assumes that the resolution matrix is storing the actual resolution convolved with the pixelization kernel along each matrix axis
+    correction = Wres2
+    correction /= pixelization_factor
+
     return correction
 
 
@@ -354,7 +441,6 @@ class Pk1D:
         __init__: Initialize class instance.
         from_fitsio: Initialize instance from a fits file.
     """
-
     def __init__(self,
                  ra,
                  dec,
diff --git a/py/picca/prep_pk1d.py b/py/picca/prep_pk1d.py
index ed40f65f2..f56667aba 100644
--- a/py/picca/prep_pk1d.py
+++ b/py/picca/prep_pk1d.py
@@ -52,7 +52,8 @@ def exp_diff(hdul, log_lambda):
 
             # exclude masks 25 (COMBINEREJ), 23 (BRIGHTSKY)?
             rebin_ivar_exp = np.bincount(log_lambda_bins,
-                                         weights=ivar_exp * (mask & 2**25 == 0))
+                                         weights=ivar_exp *
+                                         (mask & 2**25 == 0))
             rebin_flux_exp = np.bincount(log_lambda_bins,
                                          weights=(ivar_exp * flux_exp *
                                                   (mask & 2**25 == 0)))
@@ -61,8 +62,10 @@ def exp_diff(hdul, log_lambda):
                 flux_total_odd[:len(rebin_ivar_exp) - 1] += rebin_flux_exp[:-1]
                 ivar_total_odd[:len(rebin_ivar_exp) - 1] += rebin_ivar_exp[:-1]
             else:
-                flux_total_even[:len(rebin_ivar_exp) - 1] += rebin_flux_exp[:-1]
-                ivar_total_even[:len(rebin_ivar_exp) - 1] += rebin_ivar_exp[:-1]
+                flux_total_even[:len(rebin_ivar_exp) -
+                                1] += rebin_flux_exp[:-1]
+                ivar_total_even[:len(rebin_ivar_exp) -
+                                1] += rebin_ivar_exp[:-1]
 
     w = ivar_total_odd > 0
     flux_total_odd[w] /= ivar_total_odd[w]
@@ -80,6 +83,72 @@ def exp_diff(hdul, log_lambda):
     return exposures_diff
 
 
+def exp_diff_desi(hdul, mask_targetid):
+    """Computes the difference between exposures.
+
+    More precisely computes de semidifference between two customized coadded
+    spectra obtained from weighted averages of the even-number exposures, for
+    the first spectrum, and of the odd-number exposures, for the second one
+    (see section 3.2 of Chabanier et al. 2019).
+
+    Args:
+        hdul: fitsio.fitslib.FITS
+            Header Data Unit List opened by fitsio
+        mask targetid: array of ints
+            Targetids to select for calculating the exp differences
+
+    Returns:
+        The difference between exposures
+    """
+    argsort = np.flip(np.argsort(hdul["TEFF_LYA"][mask_targetid][:]))
+    flux = hdul["FL"][mask_targetid][argsort, :]
+    ivar = hdul["IV"][mask_targetid][argsort, :]
+    teff_lya = hdul["TEFF_LYA"][mask_targetid][argsort]
+
+    num_exp = len(flux)
+    if (num_exp < 2):
+        userprint("DBG : not enough exposures for diff, spectra rejected")
+        return None
+    flux_total_odd = np.zeros(flux.shape[1])
+    ivar_total_odd = np.zeros(flux.shape[1])
+    flux_total_even = np.zeros(flux.shape[1])
+    ivar_total_even = np.zeros(flux.shape[1])
+    teff_even = 0
+    teff_odd = 0
+    teff_total = np.sum(teff_lya)
+    teff_last = teff_lya[-1]
+    for index_exp in range(2 * (num_exp // 2)):
+        flexp = flux[index_exp]
+        ivexp = ivar[index_exp]
+        teff_lya_exp = teff_lya[index_exp]
+        if index_exp % 2 == 1:
+            flux_total_odd += flexp * ivexp
+            ivar_total_odd += ivexp
+            teff_odd += teff_lya_exp
+        else:
+            flux_total_even += flexp * ivexp
+            ivar_total_even += ivexp
+            teff_even += teff_lya_exp
+
+    w = ivar_total_odd > 0
+    flux_total_odd[w] /= ivar_total_odd[w]
+    w = ivar_total_even > 0
+    flux_total_even[w] /= ivar_total_even[w]
+
+    alpha = 1
+    if (num_exp % 2 == 1):
+        n_even = (num_exp - 1) // 2
+        alpha_N_old = np.sqrt(4. * n_even * (num_exp - n_even)) / num_exp
+        # alpha_N = np.sqrt(4.*t_even*(t_exp-t_even))/t_exp
+        alpha_C_new = np.sqrt((teff_total - teff_last) / teff_total)
+        alpha_N_new = np.sqrt(
+            (teff_total - teff_last) * (teff_total + teff_last)) / teff_total
+        alpha = alpha_N_new
+    diff = 0.5 * (flux_total_even - flux_total_odd) * alpha
+
+    return diff
+
+
 def spectral_resolution(wdisp,
                         with_correction=False,
                         fiberid=None,
@@ -112,11 +181,11 @@ def spectral_resolution(wdisp,
         # fiberids greater than 500 corresponds to the second spectrograph
         fiberid = fiberid % 500
         if fiberid < 100:
-            correction = (1. + (correction - 1) * .25 + (correction - 1) * .75 *
-                          (fiberid) / 100.)
+            correction = (1. + (correction - 1) * .25 +
+                          (correction - 1) * .75 * (fiberid) / 100.)
         elif fiberid > 400:
-            correction = (1. + (correction - 1) * .25 + (correction - 1) * .75 *
-                          (500 - fiberid) / 100.)
+            correction = (1. + (correction - 1) * .25 +
+                          (correction - 1) * .75 * (500 - fiberid) / 100.)
 
         # apply the correction
         reso *= correction
@@ -126,6 +195,7 @@ def spectral_resolution(wdisp,
 
 def spectral_resolution_desi(reso_matrix, log_lambda):
     """Computes the spectral resolution for DESI spectra
+    Note that this is only giving rough estimates, it relies on a Gaussian resolution matrix
 
     Args:
         reso_matrix: array
@@ -139,12 +209,19 @@ def spectral_resolution_desi(reso_matrix, log_lambda):
     delta_log_lambda = ((log_lambda[-1] - log_lambda[0]) /
                         float(len(log_lambda) - 1))
     reso = np.clip(reso_matrix, 1.0e-6, 1.0e6)
-    rms_in_pixel = (np.sqrt(1.0 / 2.0 / np.log(
-        reso[len(reso) // 2][:] / reso[len(reso) // 2 - 1][:])) + np.sqrt(
-            4.0 / 2.0 / np.log(
-                reso[len(reso) // 2][:] / reso[len(reso) // 2 - 2][:]))) / 2.0
-
-    reso_in_km_per_s = (rms_in_pixel * SPEED_LIGHT * delta_log_lambda *
-                        np.log(10.0))
 
-    return reso_in_km_per_s
+    #
+    rms_in_pixel = (
+        np.sqrt(1.0 / 2.0 / np.log(
+            reso[len(reso) // 2][:] / reso[len(reso) // 2 - 1][:])) +
+        np.sqrt(4.0 / 2.0 /
+                np.log(reso[len(reso) // 2][:] / reso[len(reso) // 2 - 2][:]))
+        + np.sqrt(1.0 / 2.0 / np.log(
+            reso[len(reso) // 2][:] / reso[len(reso) // 2 + 1][:])) +
+        np.sqrt(4.0 / 2.0 / np.log(
+            reso[len(reso) // 2][:] / reso[len(reso) // 2 + 2][:]))) / 4.0
+
+    avg_reso_in_km_per_s = (rms_in_pixel * SPEED_LIGHT * delta_log_lambda *
+                            np.log(10.0))
+
+    return rms_in_pixel, avg_reso_in_km_per_s


Test patch:
diff --git a/py/picca/tests/delta_extraction/astronomical_object_tests.py b/py/picca/tests/delta_extraction/astronomical_object_tests.py
index 7123861db..bd960e909 100644
--- a/py/picca/tests/delta_extraction/astronomical_object_tests.py
+++ b/py/picca/tests/delta_extraction/astronomical_object_tests.py
@@ -28,10 +28,18 @@
 EXPOSURES_DIFF2 = np.ones(10) * 3
 EXPOSURES_DIFF_REBIN = np.ones(5)
 EXPOSURES_DIFF_COADD = np.ones(5) * 2
+
 RESO = np.ones(10)
 RESO2 = np.ones(10) * 3
 RESO_REBIN = np.ones(5)
 RESO_COADD = np.ones(5) * 2
+
+
+RESO_PIX = np.ones(10)
+RESO_PIX2 = np.ones(10) * 3
+RESO_PIX_REBIN = np.ones(5)
+RESO_PIX_COADD = np.ones(5) * 2
+
 LOG_LAMBDA = np.array([
     3.5562825, 3.5563225, 3.5565825, 3.5566225, 3.5568825, 3.5569225,
     3.5571825, 3.5572225, 3.5574825, 3.5575225
@@ -42,6 +50,12 @@
     [3610.0, 3610.4, 3650.0, 3650.4, 3670.0, 3670.4, 3680.0, 3680.4, 3700.0, 3700.4])
 LAMBDA_REBIN = np.array([3610, 3650, 3670, 3680, 3700])
 
+RESOMAT = np.ones([7, 10])
+RESOMAT2 = np.ones([7, 10])*3
+
+RESOMAT_REBIN = np.ones([7, 5])
+RESOMAT_COADD = np.ones([7, 5])*2
+
 THINGID = 100000000
 TARGETID = 100000000
 
@@ -151,12 +165,15 @@
 kwargs_pk1d_forest.update({
     "exposures_diff": EXPOSURES_DIFF,
     "reso": RESO,
+    "reso_pix": RESO_PIX
 })
 
 kwargs_pk1d_forest2 = kwargs_forest2.copy()
 kwargs_pk1d_forest2.update({
     "exposures_diff": EXPOSURES_DIFF2,
     "reso": RESO2,
+    "reso_pix": RESO_PIX2
+
 })
 
 kwargs_pk1d_forest_log = kwargs_pk1d_forest.copy()
@@ -183,6 +200,7 @@
 kwargs_pk1d_forest_rebin.update({
     "exposures_diff": EXPOSURES_DIFF_REBIN,
     "reso": RESO_REBIN,
+    "reso_pix": RESO_PIX_REBIN
 })
 
 kwargs_pk1d_forest_log_rebin = kwargs_pk1d_forest_rebin.copy()
@@ -199,6 +217,7 @@
 kwargs_pk1d_forest_coadd.update({
     "exposures_diff": EXPOSURES_DIFF_COADD,
     "reso": RESO_COADD,
+    "reso_pix": RESO_PIX_COADD
 })
 
 kwargs_pk1d_forest_log_coadd = kwargs_pk1d_forest_coadd.copy()
@@ -251,17 +270,29 @@
 # define contructors for DesiPk1dForest
 kwargs_desi_pk1d_forest = kwargs_desi_forest.copy()
 kwargs_desi_pk1d_forest.update(kwargs_pk1d_forest_lin)
+kwargs_desi_pk1d_forest.update({
+    "resolution_matrix": RESOMAT,
+})
 del kwargs_desi_pk1d_forest["los_id"]
 
 kwargs_desi_pk1d_forest2 = kwargs_desi_forest2.copy()
 kwargs_desi_pk1d_forest2.update(kwargs_pk1d_forest_lin2)
+kwargs_desi_pk1d_forest2.update({
+    "resolution_matrix": RESOMAT2,
+})
 del kwargs_desi_pk1d_forest2["los_id"]
 
 kwargs_desi_pk1d_forest_rebin = kwargs_pk1d_forest_rebin.copy()
 kwargs_desi_pk1d_forest_rebin.update(kwargs_desi_forest_rebin)
+kwargs_desi_pk1d_forest_rebin.update({
+    "resolution_matrix": RESOMAT_REBIN,
+})
 
 kwargs_desi_pk1d_forest_coadd = kwargs_pk1d_forest_coadd.copy()
 kwargs_desi_pk1d_forest_coadd.update(kwargs_desi_forest_coadd)
+kwargs_desi_pk1d_forest_coadd.update({
+    "resolution_matrix": RESOMAT_COADD,
+})
 
 # define contructors for SdssForest
 kwargs_sdss_forest = kwargs_forest_log.copy()
@@ -414,8 +445,10 @@ def assert_forest_object(self, test_obj, kwargs):
         self.assertTrue(np.allclose(test_obj.flux, flux))
         self.assertTrue(np.allclose(test_obj.ivar, ivar))
 
-        if isinstance(test_obj, Pk1dForest):
-            self.assertTrue(len(Forest.mask_fields) == 6)
+        if isinstance(test_obj, DesiPk1dForest):
+            self.assertTrue(len(Forest.mask_fields) == 8)
+        elif isinstance(test_obj, Pk1dForest):
+            self.assertTrue(len(Forest.mask_fields) == 7)
         else:
             self.assertTrue(len(Forest.mask_fields) == 4)
         self.assertTrue(Forest.mask_fields[0] == "flux")
@@ -454,6 +487,8 @@ def assert_forest_object(self, test_obj, kwargs):
                 np.allclose(test_obj.exposures_diff,
                             kwargs.get("exposures_diff")))
             self.assertTrue(np.allclose(test_obj.reso, kwargs.get("reso")))
+            self.assertTrue(np.allclose(test_obj.reso_pix, kwargs.get("reso_pix")))
+
             if Forest.wave_solution == "log":
                 log_lambda = kwargs.get("log_lambda")
                 mean_z = ((np.power(10., log_lambda[len(log_lambda) - 1]) +
@@ -569,7 +604,11 @@ def assert_get_header(self, test_obj):
             self.assertTrue(header[index + 2].get("name") == "MEANRESO")
             self.assertTrue(header[index +
                                    2].get("value") == test_obj.mean_reso)
-            index += 2
+            self.assertTrue(header[index + 3].get("name") == "MEANRESO_PIX")
+            self.assertTrue(header[index +
+                                   3].get("value") == test_obj.mean_reso_pix)
+            
+            index += 3
         if isinstance(test_obj, SdssForest):
             self.assertTrue(header[index + 1].get("name") == "THING_ID")
             self.assertTrue(header[index + 1].get("value") == test_obj.thingid)
@@ -840,12 +879,14 @@ def test_desi_pk1d_forest(self):
                     3610, 3610.4, 3650, 3650.4, 3670, 3670.4, 3680, 3680.4,
                     3700, 3700.4
                 ]),
-            "targetid":
-                100000000,
             "exposures_diff":
                 np.ones(10),
             "reso":
                 np.ones(10),
+            "reso_pix":
+                np.ones(10),
+            "resolution_matrix":
+                np.ones([7, 10])
         }
         with self.assertRaises(AstronomicalObjectError):
             DesiPk1dForest(**kwargs)
@@ -869,10 +910,42 @@ def test_desi_pk1d_forest(self):
                 ]),
             "targetid":
                 100000000,
+            "resolution_matrix":
+                np.ones([7, 10])
         }
         with self.assertRaises(AstronomicalObjectError):
             DesiPk1dForest(**kwargs)
 
+        # create a DesiForest with missing DesiPk1dForest variables
+        kwargs = {
+            "ra":
+                0.15,
+            "dec":
+                0.0,
+            "z":
+                2.1,
+            "flux":
+                np.ones(15),
+            "ivar":
+                np.ones(15) * 4,
+            "lambda":
+                np.array([
+                    3610, 3610.4, 3650, 3650.4, 3670, 3670.4, 3680, 3680.4,
+                    3700, 3700.4
+                ]),
+            "exposures_diff":
+                np.ones(10),
+            "targetid":
+                100000000,
+            "reso":
+                np.ones(10),
+            "reso_pix":
+                np.ones(10),
+        }
+        with self.assertRaises(AstronomicalObjectError):
+            DesiPk1dForest(**kwargs)
+
+
         # create forest with missing Forest variables
         kwargs = {
             "los_id": 9999,
@@ -886,6 +959,7 @@ def test_desi_pk1d_forest(self):
             "fiber": 0,
             "exposures_diff": np.ones(10),
             "reso": np.ones(10),
+            "reso_pix": np.ones(15)
         }
         with self.assertRaises(AstronomicalObjectError):
             DesiPk1dForest(**kwargs)
@@ -1433,6 +1507,8 @@ def test_sdss_pk1d_forest(self):
                 np.ones(15),
             "reso":
                 np.ones(15),
+            "reso_pix":
+                np.ones(15),
         }
         with self.assertRaises(AstronomicalObjectError):
             SdssPk1dForest(**kwargs)
@@ -1465,6 +1541,7 @@ def test_sdss_pk1d_forest(self):
             "mjd": 0,
             "exposures_diff": np.ones(15),
             "reso": np.ones(15),
+            "reso_pix": np.ones(15),
         }
         with self.assertRaises(AstronomicalObjectError):
             SdssPk1dForest(**kwargs)
diff --git a/py/picca/tests/delta_extraction/data_tests.py b/py/picca/tests/delta_extraction/data_tests.py
index 11c62ed4d..654f674e8 100644
--- a/py/picca/tests/delta_extraction/data_tests.py
+++ b/py/picca/tests/delta_extraction/data_tests.py
@@ -88,6 +88,7 @@ def test_data_filter_forests(self):
                             "rejection log file": "rejection_log.fits.gz",
                             "wave solution": "log",
                             "delta log lambda": 3e-4,
+                            "delta log lambda rest frame": 3e-4,
                             "input directory": f"{THIS_DIR}/data/",
                         }})
         for key, value in defaults_data.items():
@@ -110,6 +111,7 @@ def test_data_filter_forests(self):
                                    "rejection log file": "rejection_log.fits.gz",
                                    "wave solution": "log",
                                    "delta log lambda": 3e-4,
+                                   "delta log lambda rest frame": 3e-4,
                                    "input directory": f"{THIS_DIR}/data/",
                          }})
         for key, value in defaults_data.items():
diff --git a/py/picca/tests/delta_extraction/test_utils.py b/py/picca/tests/delta_extraction/test_utils.py
index bf8a8b105..33564b6f4 100644
--- a/py/picca/tests/delta_extraction/test_utils.py
+++ b/py/picca/tests/delta_extraction/test_utils.py
@@ -40,8 +40,9 @@ def setup_forest(wave_solution, rebin=1):
         pixel_step = 1e-4 * rebin
     elif wave_solution == "lin":
         pixel_step = 1.
+    pixel_step_rf = pixel_step
 
-    Forest.set_class_variables(3600.0, 5500.0, 1040.0, 1200.0, pixel_step,
+    Forest.set_class_variables(3600.0, 5500.0, 1040.0, 1200.0, pixel_step, pixel_step_rf,
                                wave_solution)
 
 setup_forest("log", rebin=3)
diff --git a/py/picca/tests/test_helpers.py b/py/picca/tests/test_helpers.py
index 00ce6786b..59c38ad12 100644
--- a/py/picca/tests/test_helpers.py
+++ b/py/picca/tests/test_helpers.py
@@ -66,14 +66,14 @@ def compare_fits(self, path1, path2, nameRun="", rel_tolerance=1e-05, abs_tolera
             for el in r_m:
                 name = el['name']
                 if len(name) > 5 and name[:5] == "TTYPE":
-                    ld_m += [el['value'].replace(" ", "")]
+                    ld_m += [el['value'].replace(" ", "").upper()]
             ###
             r_b = b[i].read_header().records()
             ld_b = []
             for el in r_b:
                 name = el['name']
                 if len(name) > 5 and name[:5] == "TTYPE":
-                    ld_b += [el['value'].replace(" ", "")]
+                    ld_b += [el['value'].replace(" ", "").upper()]
 
             self.assertListEqual(ld_m, ld_b, "{}".format(nameRun))