Repository: igmhub/picca (language: python)

Issue:
Add constant weights for true cont
This adds the option to use constant weights for the true continuum analysis, which works without the raw statistics file. Also fixes a bug in the current true cont calculation.


Discussion / hints:

Fix patch:
diff --git a/docs/delta_extraction/data_model/data_model-Data.png b/docs/delta_extraction/data_model/data_model-Data.png
index 9f3fd4c21..69c87bf91 100644
Binary files a/docs/delta_extraction/data_model/data_model-Data.png and b/docs/delta_extraction/data_model/data_model-Data.png differ
diff --git a/docs/delta_extraction/data_model/data_model-ExpectedFlux.png b/docs/delta_extraction/data_model/data_model-ExpectedFlux.png
index c07f4e276..de4807727 100644
Binary files a/docs/delta_extraction/data_model/data_model-ExpectedFlux.png and b/docs/delta_extraction/data_model/data_model-ExpectedFlux.png differ
diff --git a/docs/delta_extraction/data_model/data_model-FullDataModel.png b/docs/delta_extraction/data_model/data_model-FullDataModel.png
index 11dcca92c..6c37412a0 100644
Binary files a/docs/delta_extraction/data_model/data_model-FullDataModel.png and b/docs/delta_extraction/data_model/data_model-FullDataModel.png differ
diff --git a/docs/delta_extraction/data_model/data_model-Main.png b/docs/delta_extraction/data_model/data_model-Main.png
index 04a763e66..dffa6987a 100644
Binary files a/docs/delta_extraction/data_model/data_model-Main.png and b/docs/delta_extraction/data_model/data_model-Main.png differ
diff --git a/docs/delta_extraction/data_model/data_model.drawio b/docs/delta_extraction/data_model/data_model.drawio
index e0d9ee827..825c43f90 100644
--- a/docs/delta_extraction/data_model/data_model.drawio
+++ b/docs/delta_extraction/data_model/data_model.drawio
@@ -1,6 +1,6 @@
-<mxfile host="app.diagrams.net" modified="2022-03-30T10:00:08.532Z" agent="5.0 (X11)" etag="ZxQMu9UQknNjdmgtxlQK" version="17.2.6" type="github" pages="7">
+<mxfile host="app.diagrams.net" modified="2022-03-30T16:34:56.373Z" agent="5.0 (X11)" etag="s_WS4omrb5KVoFcUA0wD" version="17.2.6" type="github" pages="7">
   <diagram id="nPyiXVHnBq5R0qJwF2HZ" name="Main">
-    <mxGraphModel dx="4543" dy="1163" grid="1" gridSize="10" guides="1" tooltips="1" connect="1" arrows="1" fold="1" page="1" pageScale="1" pageWidth="827" pageHeight="1169" math="0" shadow="0">
+    <mxGraphModel dx="3903" dy="923" grid="1" gridSize="10" guides="1" tooltips="1" connect="1" arrows="1" fold="1" page="1" pageScale="1" pageWidth="827" pageHeight="1169" math="0" shadow="0">
       <root>
         <mxCell id="0" />
         <mxCell id="1" parent="0" />
@@ -91,7 +91,7 @@
         <mxCell id="4YIj42Dbc83F9fyeYHnL-46" value="Forest" style="swimlane;fontStyle=3;childLayout=stackLayout;horizontal=1;startSize=26;fillColor=#ffcd28;horizontalStack=0;resizeParent=1;resizeParentMax=0;resizeLast=0;collapsible=1;marginBottom=0;swimlaneFillColor=#ffffff;rounded=0;shadow=0;glass=0;comic=0;strokeColor=#000000;gradientColor=#ffa500;swimlaneLine=1;gradientDirection=north;" parent="1" vertex="1">
           <mxGeometry x="-660" y="870" width="270" height="501" as="geometry" />
         </mxCell>
-        <mxCell id="4YIj42Dbc83F9fyeYHnL-47" value="+ bad_continuum_reason[0..1]: str&#xa;&gt; blinding: str&#xa;+ continuum[0..1]: array of float&#xa;+ deltas[0..1]: array of float&#xa;+ flux: array of float&#xa;+ ivar: array of float&#xa;+ lambda_[0..1]: array of float&#xa;&gt; lambda_grid[0..1]: array of float&#xa;&gt; lambda_rest_frame_grid[0..1]: array of float&#xa;+ log_lambda[0..1]: array of float&#xa;&gt; log_lambda_grid[0..1]: array of float&#xa;&gt; log_lambda_rest_frame_grid[0..1]: array of float&#xa;&gt; mask_fields: list of str&#xa;+ mean_snr: float&#xa;+ transmission_correction: array of float&#xa;&gt; wave_solution: &quot;lin&quot; or &quot;log&quot;&#xa;+ weights[0..1]: array of float" style="text;strokeColor=#000000;fillColor=#ffcd28;align=left;verticalAlign=middle;overflow=hidden;rotatable=0;points=[[0,0.5],[1,0.5]];portConstraint=eastwest;rounded=0;shadow=0;glass=0;comic=0;fontStyle=0;gradientColor=none;swimlaneFillColor=#ffffff;gradientDirection=north;startSize=26;" parent="4YIj42Dbc83F9fyeYHnL-46" vertex="1">
+        <mxCell id="4YIj42Dbc83F9fyeYHnL-47" value="+ bad_continuum_reason[0..1]: str&#xa;&gt; blinding: str&#xa;+ continuum[0..1]: array of float&#xa;+ deltas[0..1]: array of float&#xa;+ flux: array of float&#xa;+ ivar: array of float&#xa;+ log_lambda: array of float&#xa;&gt; log_lambda_grid: array of float&#xa;&gt; log_lambda_rest_frame_grid: array of float&#xa;&gt; mask_fields: list of str&#xa;+ mean_snr: float&#xa;+ transmission_correction: array of float&#xa;&gt; wave_solution: &quot;lin&quot; or &quot;log&quot;&#xa;+ weights[0..1]: array of float" style="text;strokeColor=#000000;fillColor=#ffcd28;align=left;verticalAlign=middle;overflow=hidden;rotatable=0;points=[[0,0.5],[1,0.5]];portConstraint=eastwest;rounded=0;shadow=0;glass=0;comic=0;fontStyle=0;gradientColor=none;swimlaneFillColor=#ffffff;gradientDirection=north;startSize=26;" parent="4YIj42Dbc83F9fyeYHnL-46" vertex="1">
           <mxGeometry y="26" width="270" height="244" as="geometry" />
         </mxCell>
         <mxCell id="4YIj42Dbc83F9fyeYHnL-48" value="" style="text;strokeColor=none;fillColor=#ffcd28;align=left;verticalAlign=top;spacingLeft=4;spacingRight=4;overflow=hidden;rotatable=0;points=[[0,0.5],[1,0.5]];portConstraint=eastwest;gradientColor=none;" parent="4YIj42Dbc83F9fyeYHnL-46" vertex="1">
@@ -409,7 +409,7 @@
         <mxCell id="4YIj42Dbc83F9fyeYHnL-184" value="Data" style="swimlane;fontStyle=3;childLayout=stackLayout;horizontal=1;startSize=26;fillColor=#ffcd28;horizontalStack=0;resizeParent=1;resizeParentMax=0;resizeLast=0;collapsible=1;marginBottom=0;swimlaneFillColor=#ffffff;rounded=0;shadow=0;glass=0;comic=0;strokeColor=#000000;gradientColor=#ffa500;swimlaneLine=1;gradientDirection=north;" parent="1" vertex="1">
           <mxGeometry x="-570" y="492" width="390" height="340" as="geometry" />
         </mxCell>
-        <mxCell id="4YIj42Dbc83F9fyeYHnL-185" value=" + analysis_type: str&#xa; + forests: list of Forest&#xa; + logger: logging.Logger&#xa; + min_num_pix: int&#xa; + min_snr: float&#xa; + out_dir: str &#xa; + rejection_log_file: str&#xa; + rejection_log_initialized: bool&#xa; + rejection_log_cols: list of list&#xa; + rejection_log_comments: list of list&#xa; + rejection_log_names: list of list" style="text;strokeColor=#000000;fillColor=#ffcd28;align=left;verticalAlign=middle;overflow=hidden;rotatable=0;points=[[0,0.5],[1,0.5]];portConstraint=eastwest;rounded=0;shadow=0;glass=0;comic=0;fontStyle=0;gradientColor=none;swimlaneFillColor=#ffffff;gradientDirection=north;startSize=26;" parent="4YIj42Dbc83F9fyeYHnL-184" vertex="1">
+        <mxCell id="4YIj42Dbc83F9fyeYHnL-185" value=" + analysis_type: str&#xa; + forests: list of Forest&#xa; + logger: logging.Logger&#xa; + min_num_pix: int&#xa; + min_snr: float&#xa; + out_dir: str&#xa; + rejection_log_file: str&#xa; + rejection_log_initialized: bool&#xa; + rejection_log_cols: list of list&#xa; + rejection_log_comments: list of list&#xa; + rejection_log_names: list of list" style="text;strokeColor=#000000;fillColor=#ffcd28;align=left;verticalAlign=middle;overflow=hidden;rotatable=0;points=[[0,0.5],[1,0.5]];portConstraint=eastwest;rounded=0;shadow=0;glass=0;comic=0;fontStyle=0;gradientColor=none;swimlaneFillColor=#ffffff;gradientDirection=north;startSize=26;" parent="4YIj42Dbc83F9fyeYHnL-184" vertex="1">
           <mxGeometry y="26" width="390" height="154" as="geometry" />
         </mxCell>
         <mxCell id="4YIj42Dbc83F9fyeYHnL-186" value="" style="line;strokeWidth=1;fillColor=none;align=left;verticalAlign=middle;spacingTop=-1;spacingLeft=3;spacingRight=3;rotatable=0;labelPosition=right;points=[];portConstraint=eastwest;rounded=0;shadow=0;glass=0;comic=0;" parent="4YIj42Dbc83F9fyeYHnL-184" vertex="1">
@@ -808,7 +808,7 @@
     </mxGraphModel>
   </diagram>
   <diagram name="Simplified Main" id="joHP-pMOU07-pkZl2VE5">
-    <mxGraphModel dx="3716" dy="1163" grid="1" gridSize="10" guides="1" tooltips="1" connect="1" arrows="1" fold="1" page="1" pageScale="1" pageWidth="827" pageHeight="1169" math="0" shadow="0">
+    <mxGraphModel dx="2672" dy="741" grid="1" gridSize="10" guides="1" tooltips="1" connect="1" arrows="1" fold="1" page="1" pageScale="1" pageWidth="827" pageHeight="1169" math="0" shadow="0">
       <root>
         <mxCell id="75n_6tslTMfhbzDpBMPq-0" />
         <mxCell id="75n_6tslTMfhbzDpBMPq-1" parent="75n_6tslTMfhbzDpBMPq-0" />
@@ -1432,23 +1432,23 @@
             <mxPoint x="10" as="offset" />
           </mxGeometry>
         </mxCell>
-        <mxCell id="B4_5931JTfAuX9vA9zEw-0" style="edgeStyle=orthogonalEdgeStyle;rounded=0;orthogonalLoop=1;jettySize=auto;html=1;exitX=0.25;exitY=2;exitDx=0;exitDy=0;endArrow=none;endFill=0;entryX=0.5;entryY=0;entryDx=0;entryDy=0;exitPerimeter=0;" parent="75n_6tslTMfhbzDpBMPq-1" source="B4_5931JTfAuX9vA9zEw-2" target="75n_6tslTMfhbzDpBMPq-23" edge="1">
+        <mxCell id="B4_5931JTfAuX9vA9zEw-0" style="edgeStyle=orthogonalEdgeStyle;rounded=0;orthogonalLoop=1;jettySize=auto;html=1;exitX=0.25;exitY=2;exitDx=0;exitDy=0;endArrow=none;endFill=0;entryX=0.5;entryY=0;entryDx=0;entryDy=0;exitPerimeter=0;" edge="1" parent="75n_6tslTMfhbzDpBMPq-1" source="B4_5931JTfAuX9vA9zEw-2" target="75n_6tslTMfhbzDpBMPq-23">
           <mxGeometry relative="1" as="geometry">
             <mxPoint x="-529.6" y="423.5160000000001" as="targetPoint" />
             <Array as="points" />
           </mxGeometry>
         </mxCell>
-        <mxCell id="B4_5931JTfAuX9vA9zEw-1" value="find_bins" style="swimlane;fontStyle=1;childLayout=stackLayout;horizontal=1;startSize=26;horizontalStack=0;resizeParent=1;resizeParentMax=0;resizeLast=0;collapsible=1;marginBottom=0;swimlaneFillColor=#ffffff;rounded=0;shadow=0;glass=0;comic=0;strokeColor=#9673a6;swimlaneLine=1;gradientDirection=north;fillColor=#e1d5e7;gradientColor=#CC99FF;" parent="75n_6tslTMfhbzDpBMPq-1" vertex="1">
+        <mxCell id="B4_5931JTfAuX9vA9zEw-1" value="find_bins" style="swimlane;fontStyle=1;childLayout=stackLayout;horizontal=1;startSize=26;horizontalStack=0;resizeParent=1;resizeParentMax=0;resizeLast=0;collapsible=1;marginBottom=0;swimlaneFillColor=#ffffff;rounded=0;shadow=0;glass=0;comic=0;strokeColor=#9673a6;swimlaneLine=1;gradientDirection=north;fillColor=#e1d5e7;gradientColor=#CC99FF;" vertex="1" parent="75n_6tslTMfhbzDpBMPq-1">
           <mxGeometry x="-730" y="410" width="120" height="27" as="geometry" />
         </mxCell>
-        <mxCell id="B4_5931JTfAuX9vA9zEw-2" value="" style="line;strokeWidth=1;fillColor=none;align=left;verticalAlign=middle;spacingTop=-1;spacingLeft=3;spacingRight=3;rotatable=0;labelPosition=right;points=[];portConstraint=eastwest;rounded=0;shadow=0;glass=0;comic=0;" parent="B4_5931JTfAuX9vA9zEw-1" vertex="1">
+        <mxCell id="B4_5931JTfAuX9vA9zEw-2" value="" style="line;strokeWidth=1;fillColor=none;align=left;verticalAlign=middle;spacingTop=-1;spacingLeft=3;spacingRight=3;rotatable=0;labelPosition=right;points=[];portConstraint=eastwest;rounded=0;shadow=0;glass=0;comic=0;" vertex="1" parent="B4_5931JTfAuX9vA9zEw-1">
           <mxGeometry y="26" width="120" height="1" as="geometry" />
         </mxCell>
       </root>
     </mxGraphModel>
   </diagram>
   <diagram id="_gUzTnbLqX_FSo8O-JG7" name="Correction">
-    <mxGraphModel dx="2889" dy="2332" grid="1" gridSize="10" guides="1" tooltips="1" connect="1" arrows="1" fold="1" page="1" pageScale="1" pageWidth="827" pageHeight="1169" math="0" shadow="0">
+    <mxGraphModel dx="1845" dy="1910" grid="1" gridSize="10" guides="1" tooltips="1" connect="1" arrows="1" fold="1" page="1" pageScale="1" pageWidth="827" pageHeight="1169" math="0" shadow="0">
       <root>
         <mxCell id="c4hXGpoebVqtJSqK0afD-0" />
         <mxCell id="c4hXGpoebVqtJSqK0afD-1" parent="c4hXGpoebVqtJSqK0afD-0" />
@@ -1860,7 +1860,7 @@
     </mxGraphModel>
   </diagram>
   <diagram id="RSB3thF5JrANvTTKpqBh" name="Mask">
-    <mxGraphModel dx="2889" dy="1163" grid="1" gridSize="10" guides="1" tooltips="1" connect="1" arrows="1" fold="1" page="1" pageScale="1" pageWidth="827" pageHeight="1169" math="0" shadow="0">
+    <mxGraphModel dx="1845" dy="741" grid="1" gridSize="10" guides="1" tooltips="1" connect="1" arrows="1" fold="1" page="1" pageScale="1" pageWidth="827" pageHeight="1169" math="0" shadow="0">
       <root>
         <mxCell id="Cd3ZDb_xAI04ongL0umg-0" />
         <mxCell id="Cd3ZDb_xAI04ongL0umg-1" parent="Cd3ZDb_xAI04ongL0umg-0" />
@@ -2139,7 +2139,7 @@
     </mxGraphModel>
   </diagram>
   <diagram id="2t2efU7mqRpUEueg9p6R" name="Data">
-    <mxGraphModel dx="2889" dy="1163" grid="1" gridSize="10" guides="1" tooltips="1" connect="1" arrows="1" fold="1" page="1" pageScale="1" pageWidth="827" pageHeight="1169" math="0" shadow="0">
+    <mxGraphModel dx="2249" dy="923" grid="1" gridSize="10" guides="1" tooltips="1" connect="1" arrows="1" fold="1" page="1" pageScale="1" pageWidth="827" pageHeight="1169" math="0" shadow="0">
       <root>
         <mxCell id="RnRJP-5kd4kx1agwHSpP-0" />
         <mxCell id="RnRJP-5kd4kx1agwHSpP-1" parent="RnRJP-5kd4kx1agwHSpP-0" />
@@ -2192,7 +2192,7 @@
         <mxCell id="RnRJP-5kd4kx1agwHSpP-102" value="Forest" style="swimlane;fontStyle=3;childLayout=stackLayout;horizontal=1;startSize=26;fillColor=#ffcd28;horizontalStack=0;resizeParent=1;resizeParentMax=0;resizeLast=0;collapsible=1;marginBottom=0;swimlaneFillColor=#ffffff;rounded=0;shadow=0;glass=0;comic=0;strokeColor=#000000;gradientColor=#ffa500;swimlaneLine=1;gradientDirection=north;" parent="RnRJP-5kd4kx1agwHSpP-1" vertex="1">
           <mxGeometry x="1360" y="641" width="250" height="510" as="geometry" />
         </mxCell>
-        <mxCell id="RnRJP-5kd4kx1agwHSpP-103" value="+ bad_continuum_reason[0..1]: str&#xa;&gt; blinding: str&#xa;+ continuum[0..1]: array of float&#xa;+ deltas[0..1]: array of float&#xa;+ flux: array of float&#xa;+ ivar: array of float&#xa;+ lambda_[0..1]: array of float&#xa;&gt; lambda_grid[0..1]: array of float&#xa;&gt; lambda_rest_frame_grid[0..1]: array of float&#xa;+ log_lambda[0..1]: array of float&#xa;&gt; log_lambda_grid[0..1]: array of float&#xa;&gt; log_lambda_rest_frame_grid[0..1]: array of float&#xa;&gt; mask_fields: list of str&#xa;+ mean_snr: float&#xa;+ transmission_correction: array of float&#xa;&gt; wave_solution: &quot;lin&quot; or &quot;log&quot;&#xa;+ weights[0..1]: array of float" style="text;strokeColor=#000000;fillColor=#ffcd28;align=left;verticalAlign=middle;overflow=hidden;rotatable=0;points=[[0,0.5],[1,0.5]];portConstraint=eastwest;rounded=0;shadow=0;glass=0;comic=0;fontStyle=0;gradientColor=none;swimlaneFillColor=#ffffff;gradientDirection=north;startSize=26;" parent="RnRJP-5kd4kx1agwHSpP-102" vertex="1">
+        <mxCell id="RnRJP-5kd4kx1agwHSpP-103" value="+ bad_continuum_reason[0..1]: str&#xa;&gt; blinding: str&#xa;+ continuum[0..1]: array of float&#xa;+ deltas[0..1]: array of float&#xa;+ flux: array of float&#xa;+ ivar: array of float&#xa;+ log_lambda: array of float&#xa;&gt; log_lambda_grid: array of float&#xa;&gt; log_lambda_rest_frame_grid: array of float&#xa;&gt; mask_fields: list of str&#xa;+ mean_snr: float&#xa;+ transmission_correction: array of float&#xa;&gt; wave_solution: &quot;lin&quot; or &quot;log&quot;&#xa;+ weights[0..1]: array of float" style="text;strokeColor=#000000;fillColor=#ffcd28;align=left;verticalAlign=middle;overflow=hidden;rotatable=0;points=[[0,0.5],[1,0.5]];portConstraint=eastwest;rounded=0;shadow=0;glass=0;comic=0;fontStyle=0;gradientColor=none;swimlaneFillColor=#ffffff;gradientDirection=north;startSize=26;" parent="RnRJP-5kd4kx1agwHSpP-102" vertex="1">
           <mxGeometry y="26" width="250" height="244" as="geometry" />
         </mxCell>
         <mxCell id="RnRJP-5kd4kx1agwHSpP-104" value="" style="text;strokeColor=none;fillColor=#ffcd28;align=left;verticalAlign=top;spacingLeft=4;spacingRight=4;overflow=hidden;rotatable=0;points=[[0,0.5],[1,0.5]];portConstraint=eastwest;gradientColor=none;" parent="RnRJP-5kd4kx1agwHSpP-102" vertex="1">
@@ -3173,7 +3173,7 @@
     </mxGraphModel>
   </diagram>
   <diagram id="je2KSbhj7LrZFA1j369E" name="ExpectedFlux">
-    <mxGraphModel dx="4543" dy="1163" grid="1" gridSize="10" guides="1" tooltips="1" connect="1" arrows="1" fold="1" page="1" pageScale="1" pageWidth="827" pageHeight="1169" math="0" shadow="0">
+    <mxGraphModel dx="3499" dy="741" grid="1" gridSize="10" guides="1" tooltips="1" connect="1" arrows="1" fold="1" page="1" pageScale="1" pageWidth="827" pageHeight="1169" math="0" shadow="0">
       <root>
         <mxCell id="tOnHf5iAFq_aPN0IRysd-0" />
         <mxCell id="tOnHf5iAFq_aPN0IRysd-1" parent="tOnHf5iAFq_aPN0IRysd-0" />
@@ -3282,41 +3282,41 @@
           </mxGeometry>
         </mxCell>
         <mxCell id="-0v5vkKTASBL-IAy4xzU-7" value="DR16ExpectedFlux" style="swimlane;fontStyle=1;childLayout=stackLayout;horizontal=1;startSize=26;fillColor=#ffcd28;horizontalStack=0;resizeParent=1;resizeParentMax=0;resizeLast=0;collapsible=1;marginBottom=0;swimlaneFillColor=#ffffff;rounded=0;shadow=0;glass=0;comic=0;strokeColor=#000000;gradientColor=#ffa500;swimlaneLine=1;gradientDirection=north;" parent="tOnHf5iAFq_aPN0IRysd-1" vertex="1">
-          <mxGeometry x="-1160" y="840" width="390" height="580" as="geometry" />
+          <mxGeometry x="-1160" y="840" width="390" height="490" as="geometry" />
         </mxCell>
-        <mxCell id="-0v5vkKTASBL-IAy4xzU-8" value="+ continuum_fit_parameters: dict&#xa;+ delta_lambda[0..1]: float&#xa;+ delta_log_lambda[0..1]: float&#xa;+ get_eta: scipy.interpolate.interp1d&#xa;+ get_fudge: scipy.interpolate.interp1d&#xa;+ get_mean_cont: scipy.interpolate.interp1d&#xa;+ get_stack_delta[0..1]: scipy.interpolate.interp1d &#xa;+ get_var_lss: scipy.interpolate.interp1d&#xa;+ iter_out_prefix: str&#xa;+ lambda_[0..1]: array of float&#xa;+ limit_eta: tuple of floats&#xa;+ limit_var_lss: tuple of floats&#xa;+ log_lambda[0..1]: array of float&#xa;+ logger: logging.Logger&#xa;+ order: int&#xa;+ out_dir: str&#xa;+ num_bins_variance: int&#xa;+ num_iterations: int&#xa;+ num_processors: int&#xa;+ use_constant_weight: boolean&#xa;+ use_ivar_as_weight: boolean&#xa;" style="text;strokeColor=none;fillColor=#ffcd28;align=left;verticalAlign=top;spacingLeft=4;spacingRight=4;overflow=hidden;rotatable=0;points=[[0,0.5],[1,0.5]];portConstraint=eastwest;gradientColor=none;" parent="-0v5vkKTASBL-IAy4xzU-7" vertex="1">
-          <mxGeometry y="26" width="390" height="354" as="geometry" />
+        <mxCell id="-0v5vkKTASBL-IAy4xzU-8" value="+ continuum_fit_parameters: dict&#xa;+ delta_lambda[0..1]: float&#xa;+ delta_log_lambda[0..1]: float&#xa;+ get_eta: scipy.interpolate.interp1d&#xa;+ get_fudge: scipy.interpolate.interp1d&#xa;+ get_mean_cont: scipy.interpolate.interp1d&#xa;+ get_stack_delta[0..1]: scipy.interpolate.interp1d &#xa;+ get_var_lss: scipy.interpolate.interp1d&#xa;+ iter_out_prefix: str&#xa;+ lambda_[0..1]: array of float&#xa;+ limit_eta: tuple of floats&#xa;+ limit_var_lss: tuple of floats&#xa;+ log_lambda_var_func_grid[0..1]: array of float&#xa;+ logger: logging.Logger&#xa;+ order: int&#xa;+ out_dir: str&#xa;+ num_bins_variance: int&#xa;+ num_iterations: int&#xa;+ num_processors: int&#xa;+ use_constant_weight: boolean&#xa;+ use_ivar_as_weight: boolean&#xa;" style="text;strokeColor=none;fillColor=#ffcd28;align=left;verticalAlign=top;spacingLeft=4;spacingRight=4;overflow=hidden;rotatable=0;points=[[0,0.5],[1,0.5]];portConstraint=eastwest;gradientColor=none;" parent="-0v5vkKTASBL-IAy4xzU-7" vertex="1">
+          <mxGeometry y="26" width="390" height="304" as="geometry" />
         </mxCell>
         <mxCell id="-0v5vkKTASBL-IAy4xzU-9" value="" style="text;strokeColor=#000000;fillColor=#ffcd28;align=left;verticalAlign=middle;overflow=hidden;rotatable=0;points=[[0,0.5],[1,0.5]];portConstraint=eastwest;rounded=0;shadow=0;glass=0;comic=0;fontStyle=0;gradientColor=none;swimlaneFillColor=#ffffff;gradientDirection=north;startSize=26;" parent="-0v5vkKTASBL-IAy4xzU-7" vertex="1">
-          <mxGeometry y="380" width="390" height="26" as="geometry" />
+          <mxGeometry y="330" width="390" height="26" as="geometry" />
         </mxCell>
         <mxCell id="-0v5vkKTASBL-IAy4xzU-10" value="" style="line;strokeWidth=1;fillColor=none;align=left;verticalAlign=middle;spacingTop=-1;spacingLeft=3;spacingRight=3;rotatable=0;labelPosition=right;points=[];portConstraint=eastwest;rounded=0;shadow=0;glass=0;comic=0;" parent="-0v5vkKTASBL-IAy4xzU-7" vertex="1">
-          <mxGeometry y="406" width="390" height="1" as="geometry" />
+          <mxGeometry y="356" width="390" height="1" as="geometry" />
         </mxCell>
-        <mxCell id="-0v5vkKTASBL-IAy4xzU-11" value=" + __init__(config: configparser.SectionProxy)&#xa; + _initialize_arrays_lin()&#xa; + _initialize_arrays_log()&#xa; + _parse_config(config: configparser.SectionProxy)&#xa; + compute_continuum(forests: list of Forest)&#xa; + compute_delta_stack(forests: list of Forest, stack_from_deltas: bool)&#xa; + compute_mean_cont_lin(forests: list of Forest)&#xa; + compute_mean_cont_log(forests: list of Forest)&#xa; + compute_expected_flux(forests: list of Forest)&#xa; + compute_var_stats(forests: list of Forest)&#xa;" style="text;strokeColor=#000000;fillColor=#ffcd28;align=left;verticalAlign=middle;overflow=hidden;rotatable=0;points=[[0,0.5],[1,0.5]];portConstraint=eastwest;rounded=0;shadow=0;glass=0;comic=0;fontStyle=0;gradientColor=none;swimlaneFillColor=#ffffff;gradientDirection=north;startSize=26;" parent="-0v5vkKTASBL-IAy4xzU-7" vertex="1">
-          <mxGeometry y="407" width="390" height="173" as="geometry" />
+        <mxCell id="-0v5vkKTASBL-IAy4xzU-11" value=" + __init__(config: configparser.SectionProxy)&#xa; + _initialize_arrays()&#xa; + _parse_config(config: configparser.SectionProxy)&#xa; + compute_continuum(forests: list of Forest)&#xa; + compute_delta_stack(forests: list of Forest, stack_from_deltas: bool)&#xa; + compute_mean_cont(forests: list of Forest)&#xa; + compute_expected_flux(forests: list of Forest)&#xa; + compute_var_stats(forests: list of Forest)&#xa;" style="text;strokeColor=#000000;fillColor=#ffcd28;align=left;verticalAlign=middle;overflow=hidden;rotatable=0;points=[[0,0.5],[1,0.5]];portConstraint=eastwest;rounded=0;shadow=0;glass=0;comic=0;fontStyle=0;gradientColor=none;swimlaneFillColor=#ffffff;gradientDirection=north;startSize=26;" parent="-0v5vkKTASBL-IAy4xzU-7" vertex="1">
+          <mxGeometry y="357" width="390" height="133" as="geometry" />
         </mxCell>
         <mxCell id="bhYFi5XTDl4pGI8U_RH3-0" value="TrueMockContinuum" style="swimlane;fontStyle=1;childLayout=stackLayout;horizontal=1;startSize=26;fillColor=#ffcd28;horizontalStack=0;resizeParent=1;resizeParentMax=0;resizeLast=0;collapsible=1;marginBottom=0;swimlaneFillColor=#ffffff;rounded=0;shadow=0;glass=0;comic=0;strokeColor=#000000;gradientColor=#ffa500;swimlaneLine=1;gradientDirection=north;" parent="tOnHf5iAFq_aPN0IRysd-1" vertex="1">
-          <mxGeometry x="-1480" y="840" width="260" height="258" as="geometry" />
+          <mxGeometry x="-1480" y="840" width="280" height="258" as="geometry" />
         </mxCell>
         <mxCell id="bhYFi5XTDl4pGI8U_RH3-1" value="+ get_mean_flux: scipy.interpolate.interp1d&#xa;+ get_var_lss: scipy.interpolate.interp1d&#xa;+ input_directory: str&#xa;+ iter_out_prefix: str&#xa;+ num_processors: int" style="text;strokeColor=none;fillColor=#ffcd28;align=left;verticalAlign=top;spacingLeft=4;spacingRight=4;overflow=hidden;rotatable=0;points=[[0,0.5],[1,0.5]];portConstraint=eastwest;gradientColor=none;" parent="bhYFi5XTDl4pGI8U_RH3-0" vertex="1">
-          <mxGeometry y="26" width="260" height="84" as="geometry" />
+          <mxGeometry y="26" width="280" height="84" as="geometry" />
         </mxCell>
         <mxCell id="bhYFi5XTDl4pGI8U_RH3-2" value="" style="text;strokeColor=#000000;fillColor=#ffcd28;align=left;verticalAlign=middle;overflow=hidden;rotatable=0;points=[[0,0.5],[1,0.5]];portConstraint=eastwest;rounded=0;shadow=0;glass=0;comic=0;fontStyle=0;gradientColor=none;swimlaneFillColor=#ffffff;gradientDirection=north;startSize=26;" parent="bhYFi5XTDl4pGI8U_RH3-0" vertex="1">
-          <mxGeometry y="110" width="260" height="26" as="geometry" />
+          <mxGeometry y="110" width="280" height="26" as="geometry" />
         </mxCell>
         <mxCell id="bhYFi5XTDl4pGI8U_RH3-3" value="" style="line;strokeWidth=1;fillColor=none;align=left;verticalAlign=middle;spacingTop=-1;spacingLeft=3;spacingRight=3;rotatable=0;labelPosition=right;points=[];portConstraint=eastwest;rounded=0;shadow=0;glass=0;comic=0;" parent="bhYFi5XTDl4pGI8U_RH3-0" vertex="1">
-          <mxGeometry y="136" width="260" height="1" as="geometry" />
+          <mxGeometry y="136" width="280" height="1" as="geometry" />
         </mxCell>
-        <mxCell id="bhYFi5XTDl4pGI8U_RH3-4" value="+ __init__(config: configpaser.SectionProxy)&#xa;+ _parse_config(config: configparser.SectionProxy)&#xa;+ compute_expected_flux(forests: list of Forest)&#xa;+ compute_mean_cont_lin(forests: list of Forest)&#xa;+ compute_mean_cont_log(forests: list of Forest)&#xa;+ read_true_continuum(forest: Forest)&#xa;+ read_var_lss()&#xa;+ save_delta_attributes()" style="text;strokeColor=#000000;fillColor=#ffcd28;align=left;verticalAlign=middle;overflow=hidden;rotatable=0;points=[[0,0.5],[1,0.5]];portConstraint=eastwest;rounded=0;shadow=0;glass=0;comic=0;fontStyle=0;gradientColor=none;swimlaneFillColor=#ffffff;gradientDirection=north;startSize=26;" parent="bhYFi5XTDl4pGI8U_RH3-0" vertex="1">
-          <mxGeometry y="137" width="260" height="121" as="geometry" />
+        <mxCell id="bhYFi5XTDl4pGI8U_RH3-4" value="+ __init__(config: configpaser.SectionProxy)&#xa;+ _parse_config(config: configparser.SectionProxy)&#xa;+ compute_expected_flux(forests: list of Forest)&#xa;+ compute_mean_cont(forests: list of Forest)&#xa;+ read_true_continuum(forest: Forest)&#xa;+ read_var_lss()&#xa;+ save_delta_attributes()" style="text;strokeColor=#000000;fillColor=#ffcd28;align=left;verticalAlign=middle;overflow=hidden;rotatable=0;points=[[0,0.5],[1,0.5]];portConstraint=eastwest;rounded=0;shadow=0;glass=0;comic=0;fontStyle=0;gradientColor=none;swimlaneFillColor=#ffffff;gradientDirection=north;startSize=26;" parent="bhYFi5XTDl4pGI8U_RH3-0" vertex="1">
+          <mxGeometry y="137" width="280" height="121" as="geometry" />
         </mxCell>
         <mxCell id="bhYFi5XTDl4pGI8U_RH3-5" value="" style="endArrow=block;endSize=16;endFill=0;html=1;exitX=0.5;exitY=0;exitDx=0;exitDy=0;entryX=0.483;entryY=1.037;entryDx=0;entryDy=0;entryPerimeter=0;" parent="tOnHf5iAFq_aPN0IRysd-1" source="bhYFi5XTDl4pGI8U_RH3-0" target="u-n-N70yIythKyl8P2P_-12" edge="1">
           <mxGeometry width="160" relative="1" as="geometry">
             <mxPoint x="-990" y="812" as="sourcePoint" />
             <mxPoint x="-1270" y="760" as="targetPoint" />
             <Array as="points">
-              <mxPoint x="-1350" y="780" />
+              <mxPoint x="-1340" y="780" />
               <mxPoint x="-1190" y="780" />
             </Array>
           </mxGeometry>
@@ -3370,12 +3370,13 @@
             <mxPoint x="-1450.24" y="540" as="targetPoint" />
           </mxGeometry>
         </mxCell>
-        <mxCell id="wSkfboFQg1evl7dp9Hlh-0" value="&lt;font color=&quot;#0000ff&quot;&gt;&amp;nbsp;Tested by&amp;nbsp;&lt;/font&gt;" style="endArrow=none;html=1;edgeStyle=orthogonalEdgeStyle;strokeColor=#0000FF;entryX=0.38;entryY=0.998;entryDx=0;entryDy=0;entryPerimeter=0;" parent="tOnHf5iAFq_aPN0IRysd-1" source="-0v5vkKTASBL-IAy4xzU-8" target="u-n-N70yIythKyl8P2P_-17" edge="1">
+        <mxCell id="wSkfboFQg1evl7dp9Hlh-0" value="&lt;font color=&quot;#0000ff&quot;&gt;&amp;nbsp;Tested by&amp;nbsp;&lt;/font&gt;" style="endArrow=none;html=1;edgeStyle=orthogonalEdgeStyle;strokeColor=#0000FF;entryX=0.38;entryY=0.998;entryDx=0;entryDy=0;entryPerimeter=0;exitX=0.003;exitY=-1;exitDx=0;exitDy=0;exitPerimeter=0;" parent="tOnHf5iAFq_aPN0IRysd-1" source="-0v5vkKTASBL-IAy4xzU-10" target="u-n-N70yIythKyl8P2P_-17" edge="1">
           <mxGeometry x="-0.8508" relative="1" as="geometry">
-            <mxPoint x="-1160" y="1193" as="sourcePoint" />
+            <mxPoint x="-1190" y="1193" as="sourcePoint" />
             <mxPoint x="-1400" y="673.7180000000001" as="targetPoint" />
             <mxPoint as="offset" />
             <Array as="points">
+              <mxPoint x="-1159" y="1193" />
               <mxPoint x="-1596" y="1193" />
             </Array>
           </mxGeometry>
@@ -3454,7 +3455,7 @@
     </mxGraphModel>
   </diagram>
   <diagram id="iPTn10C6GmmupZz-RSug" name="FullDataModel">
-    <mxGraphModel dx="3716" dy="2332" grid="1" gridSize="10" guides="1" tooltips="1" connect="1" arrows="1" fold="1" page="1" pageScale="1" pageWidth="827" pageHeight="1169" math="0" shadow="0">
+    <mxGraphModel dx="3076" dy="2092" grid="1" gridSize="10" guides="1" tooltips="1" connect="1" arrows="1" fold="1" page="1" pageScale="1" pageWidth="827" pageHeight="1169" math="0" shadow="0">
       <root>
         <mxCell id="A5sfwNRnfNnv_qPRvbdW-0" />
         <mxCell id="A5sfwNRnfNnv_qPRvbdW-1" parent="A5sfwNRnfNnv_qPRvbdW-0" />
@@ -3734,7 +3735,7 @@
         <mxCell id="A5sfwNRnfNnv_qPRvbdW-44" value="Forest" style="swimlane;fontStyle=3;childLayout=stackLayout;horizontal=1;startSize=26;fillColor=#ffcd28;horizontalStack=0;resizeParent=1;resizeParentMax=0;resizeLast=0;collapsible=1;marginBottom=0;swimlaneFillColor=#ffffff;rounded=0;shadow=0;glass=0;comic=0;strokeColor=#000000;gradientColor=#ffa500;swimlaneLine=1;gradientDirection=north;" parent="A5sfwNRnfNnv_qPRvbdW-1" vertex="1">
           <mxGeometry x="1280" y="560" width="250" height="517" as="geometry" />
         </mxCell>
-        <mxCell id="A5sfwNRnfNnv_qPRvbdW-45" value="+ bad_continuum_reason[0..1]: str&#xa;&gt; blinding: str&#xa;+ continuum[0..1]: array of float&#xa;+ deltas[0..1]: array of float&#xa;+ flux: array of float&#xa;+ ivar: array of float&#xa;+ lambda_[0..1]: array of float&#xa;&gt; lambda_grid[0..1]: array of float&#xa;&gt; lambda_rest_frame_grid[0..1]: array of float&#xa;+ log_lambda[0..1]: array of float&#xa;&gt; log_lambda_grid[0..1]: array of float&#xa;&gt; log_lambda_rest_frame_grid[0..1]: array of float&#xa;&gt; mask_fields: list of str&#xa;+ mean_snr: float&#xa;+ transmission_correction: array of float&#xa;&gt; wave_solution: &quot;lin&quot; or &quot;log&quot;&#xa;+ weights[0..1]: array of float" style="text;strokeColor=#000000;fillColor=#ffcd28;align=left;verticalAlign=middle;overflow=hidden;rotatable=0;points=[[0,0.5],[1,0.5]];portConstraint=eastwest;rounded=0;shadow=0;glass=0;comic=0;fontStyle=0;gradientColor=none;swimlaneFillColor=#ffffff;gradientDirection=north;startSize=26;" parent="A5sfwNRnfNnv_qPRvbdW-44" vertex="1">
+        <mxCell id="A5sfwNRnfNnv_qPRvbdW-45" value="+ bad_continuum_reason[0..1]: str&#xa;&gt; blinding: str&#xa;+ continuum[0..1]: array of float&#xa;+ deltas[0..1]: array of float&#xa;+ flux: array of float&#xa;+ ivar: array of float&#xa;+ log_lambda: array of float&#xa;&gt; log_lambda_grid: array of float&#xa;&gt; log_lambda_rest_frame_grid: array of float&#xa;&gt; mask_fields: list of str&#xa;+ mean_snr: float&#xa;+ transmission_correction: array of float&#xa;&gt; wave_solution: &quot;lin&quot; or &quot;log&quot;&#xa;+ weights[0..1]: array of float" style="text;strokeColor=#000000;fillColor=#ffcd28;align=left;verticalAlign=middle;overflow=hidden;rotatable=0;points=[[0,0.5],[1,0.5]];portConstraint=eastwest;rounded=0;shadow=0;glass=0;comic=0;fontStyle=0;gradientColor=none;swimlaneFillColor=#ffffff;gradientDirection=north;startSize=26;" parent="A5sfwNRnfNnv_qPRvbdW-44" vertex="1">
           <mxGeometry y="26" width="250" height="254" as="geometry" />
         </mxCell>
         <mxCell id="A5sfwNRnfNnv_qPRvbdW-46" value="" style="text;strokeColor=none;fillColor=#ffcd28;align=left;verticalAlign=top;spacingLeft=4;spacingRight=4;overflow=hidden;rotatable=0;points=[[0,0.5],[1,0.5]];portConstraint=eastwest;gradientColor=none;" parent="A5sfwNRnfNnv_qPRvbdW-44" vertex="1">
@@ -3848,7 +3849,7 @@
         <mxCell id="A5sfwNRnfNnv_qPRvbdW-82" value="" style="line;strokeWidth=1;fillColor=none;align=left;verticalAlign=middle;spacingTop=-1;spacingLeft=3;spacingRight=3;rotatable=0;labelPosition=right;points=[];portConstraint=eastwest;rounded=0;shadow=0;glass=0;comic=0;" parent="A5sfwNRnfNnv_qPRvbdW-79" vertex="1">
           <mxGeometry y="376" width="390" height="1" as="geometry" />
         </mxCell>
-        <mxCell id="A5sfwNRnfNnv_qPRvbdW-83" value=" + __init__(config: configparser.SectionProxy)&#xa; + _initialize_arrays_lin()&#xa; + _initialize_arrays_log()&#xa; + _parse_config(config: configparser.SectionProxy)&#xa; + compute_continuum(forests: list of Forest)&#xa; + compute_delta_stack(forests: list of Forest, stack_from_deltas: bool)&#xa; + compute_mean_cont_lin(forests: list of Forest)&#xa; + compute_mean_cont_log(forests: list of Forest)&#xa; + compute_expected_flux(forests: list of Forest)&#xa; + compute_var_stats(forests: list of Forest)&#xa;" style="text;strokeColor=#000000;fillColor=#ffcd28;align=left;verticalAlign=middle;overflow=hidden;rotatable=0;points=[[0,0.5],[1,0.5]];portConstraint=eastwest;rounded=0;shadow=0;glass=0;comic=0;fontStyle=0;gradientColor=none;swimlaneFillColor=#ffffff;gradientDirection=north;startSize=26;" parent="A5sfwNRnfNnv_qPRvbdW-79" vertex="1">
+        <mxCell id="A5sfwNRnfNnv_qPRvbdW-83" value=" + __init__(config: configparser.SectionProxy)&#xa; + _initialize_arrays()&#xa; + _parse_config(config: configparser.SectionProxy)&#xa; + compute_continuum(forests: list of Forest)&#xa; + compute_delta_stack(forests: list of Forest, stack_from_deltas: bool)&#xa; + compute_mean_cont(forests: list of Forest)&#xa; + compute_expected_flux(forests: list of Forest)&#xa; + compute_var_stats(forests: list of Forest)&#xa;" style="text;strokeColor=#000000;fillColor=#ffcd28;align=left;verticalAlign=middle;overflow=hidden;rotatable=0;points=[[0,0.5],[1,0.5]];portConstraint=eastwest;rounded=0;shadow=0;glass=0;comic=0;fontStyle=0;gradientColor=none;swimlaneFillColor=#ffffff;gradientDirection=north;startSize=26;" parent="A5sfwNRnfNnv_qPRvbdW-79" vertex="1">
           <mxGeometry y="377" width="390" height="173" as="geometry" />
         </mxCell>
         <mxCell id="A5sfwNRnfNnv_qPRvbdW-84" value="TrueMockContinuum" style="swimlane;fontStyle=1;childLayout=stackLayout;horizontal=1;startSize=26;fillColor=#ffcd28;horizontalStack=0;resizeParent=1;resizeParentMax=0;resizeLast=0;collapsible=1;marginBottom=0;swimlaneFillColor=#ffffff;rounded=0;shadow=0;glass=0;comic=0;strokeColor=#000000;gradientColor=#ffa500;swimlaneLine=1;gradientDirection=north;" parent="A5sfwNRnfNnv_qPRvbdW-1" vertex="1">
@@ -3863,7 +3864,7 @@
         <mxCell id="A5sfwNRnfNnv_qPRvbdW-87" value="" style="line;strokeWidth=1;fillColor=none;align=left;verticalAlign=middle;spacingTop=-1;spacingLeft=3;spacingRight=3;rotatable=0;labelPosition=right;points=[];portConstraint=eastwest;rounded=0;shadow=0;glass=0;comic=0;" parent="A5sfwNRnfNnv_qPRvbdW-84" vertex="1">
           <mxGeometry y="136" width="280" height="1" as="geometry" />
         </mxCell>
-        <mxCell id="A5sfwNRnfNnv_qPRvbdW-88" value="+ __init__(config: configpaser.SectionProxy) &#xa;+ _parse_config(config: configparser.SectionProxy)&#xa;+ compute_expected_flux(forests: list of Forest)&#xa;+ compute_mean_cont_lin(forests: list of Forest)&#xa;+ compute_mean_cont_log(forests: list of Forest)&#xa;+ read_true_continuum(forest: Forest)&#xa;+ read_var_lss()&#xa;+ save_delta_attributes()" style="text;strokeColor=#000000;fillColor=#ffcd28;align=left;verticalAlign=middle;overflow=hidden;rotatable=0;points=[[0,0.5],[1,0.5]];portConstraint=eastwest;rounded=0;shadow=0;glass=0;comic=0;fontStyle=0;gradientColor=none;swimlaneFillColor=#ffffff;gradientDirection=north;startSize=26;" parent="A5sfwNRnfNnv_qPRvbdW-84" vertex="1">
+        <mxCell id="A5sfwNRnfNnv_qPRvbdW-88" value="+ __init__(config: configpaser.SectionProxy) &#xa;+ _parse_config(config: configparser.SectionProxy)&#xa;+ compute_expected_flux(forests: list of Forest)&#xa;+ compute_mean_cont(forests: list of Forest)&#xa;+ read_true_continuum(forest: Forest)&#xa;+ read_var_lss()&#xa;+ save_delta_attributes()" style="text;strokeColor=#000000;fillColor=#ffcd28;align=left;verticalAlign=middle;overflow=hidden;rotatable=0;points=[[0,0.5],[1,0.5]];portConstraint=eastwest;rounded=0;shadow=0;glass=0;comic=0;fontStyle=0;gradientColor=none;swimlaneFillColor=#ffffff;gradientDirection=north;startSize=26;" parent="A5sfwNRnfNnv_qPRvbdW-84" vertex="1">
           <mxGeometry y="137" width="280" height="133" as="geometry" />
         </mxCell>
         <mxCell id="A5sfwNRnfNnv_qPRvbdW-89" value="" style="endArrow=block;endSize=16;endFill=0;html=1;exitX=0.5;exitY=0;exitDx=0;exitDy=0;entryX=0.5;entryY=1;entryDx=0;entryDy=0;entryPerimeter=0;" parent="A5sfwNRnfNnv_qPRvbdW-1" source="A5sfwNRnfNnv_qPRvbdW-84" target="A5sfwNRnfNnv_qPRvbdW-94" edge="1">
diff --git a/py/picca/delta_extraction/astronomical_objects/desi_forest.py b/py/picca/delta_extraction/astronomical_objects/desi_forest.py
index b2841b112..a820730da 100644
--- a/py/picca/delta_extraction/astronomical_objects/desi_forest.py
+++ b/py/picca/delta_extraction/astronomical_objects/desi_forest.py
@@ -19,78 +19,11 @@ class DesiForest(Forest):
 
     Class Attributes
     ----------------
-    blinding: str (from Forest)
-    Name of the blinding strategy used
-
-    lambda_grid: array of float or None (from Forest)
-    Common grid in lambda_ based on the specified minimum and maximum
-    wavelengths, and delta_lambda.
-
-    lambda_rest_frame_grid: array of float or None (from Forest)
-    Same as lambda_grid but for rest-frame wavelengths
-
-    log_lambda_grid: array of float or None (from Forest)
-    Common grid in log_lambda based on the specified minimum and maximum
-    wavelengths, and delta_log_lambda.
-
-    log_lambda_rest_frame_grid: array of float or None (from Forest)
-    Same as log_lambda_grid but for rest-frame wavelengths.
-
-    mask_fields: list of str (from Forest)
-    Names of the fields that are affected by masking. In general it will
-    be "flux" and "ivar" but some child classes might add more.
-
-    wave_solution: "lin" or "log" (from Forest)
-    Determines whether the wavelength solution has linear spacing ("lin") or
-    logarithmic spacing ("log").
-
+    (see Forest in py/picca/delta_extraction/astronomical_objects/forest.py)
+    
     Attributes
     ----------
-    dec: float (from AstronomicalObject)
-    Declination (in rad)
-
-    healpix: int (from AstronomicalObject)
-    Healpix number associated with (ra, dec)
-
-    los_id: longint (from AstronomicalObject)
-    Line-of-sight id. Same as targetid
-
-    ra: float (from AstronomicalObject)
-    Right ascention (in rad)
-
-    z: float (from AstronomicalObject)
-    Redshift
-
-    bad_continuum_reason: str or None
-    Reason as to why the continuum fit is not acceptable. None for acceptable
-    contiuum.
-
-    continuum: array of float or None (from Forest)
-    Quasar continuum. None for no information
-
-    deltas: array of float or None (from Forest)
-    Flux-transmission field (delta field). None for no information
-
-    flux: array of float (from Forest)
-    Flux
-
-    ivar: array of float (from Forest)
-    Inverse variance
-
-    lambda_: array of float or None (from Forest)
-    Wavelength (in Angstroms)
-
-    log_lambda: array of float or None (from Forest)
-    Logarithm of the wavelength (in Angstroms)
-
-    mean_snr: float (from Forest)
-    Mean signal-to-noise of the forest
-
-    transmission_correction: array of float (from Forest)
-    Transmission correction.
-
-    weights: array of float or None (from Forest)
-    Weights associated to the delta field. None for no information
+    (see Forest in py/picca/delta_extraction/astronomical_objects/forest.py)
 
     night: list of int
     Identifier of the night where the observation was made. None for no info
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 4494dd620..35165f165 100644
--- a/py/picca/delta_extraction/astronomical_objects/desi_pk1d_forest.py
+++ b/py/picca/delta_extraction/astronomical_objects/desi_pk1d_forest.py
@@ -24,107 +24,17 @@ class DesiPk1dForest(DesiForest, Pk1dForest):
     get_header (from DesiForest, Pk1dForest)
     __init__
 
-
     Class Attributes
     ----------------
-    blinding: str (from Forest)
-    Name of the blinding strategy used
-    
-    lambda_grid: array of float or None (from Forest)
-    Common grid in lambda_ based on the specified minimum and maximum
-    wavelengths, and delta_lambda.
-
-    lambda_rest_frame_grid: array of float or None (from Forest)
-    Same as lambda_grid but for rest-frame wavelengths
-
-    log_lambda_grid: array of float or None (from Forest)
-    Common grid in log_lambda based on the specified minimum and maximum
-    wavelengths, and delta_log_lambda.
-
-    log_lambda_rest_frame_grid: array of float or None (from Forest)
-    Same as log_lambda_grid but for rest-frame wavelengths.
-
-    mask_fields: list of str (from Forest)
-    Names of the fields that are affected by masking. In general it will
-    be "flux" and "ivar" but some child classes might add more.
-
-    wave_solution: "lin" or "log" (from Forest)
-    Determines whether the wavelength solution has linear spacing ("lin") or
-    logarithmic spacing ("log").
-
-    lambda_abs_igm: float (from Pk1dForest)
-    Wavelength of the IGM absorber
+    (see DesiForest in py/picca/delta_extraction/astronomical_objects/desi_forest.py)
+    (see Pk1dForest in py/picca/delta_extraction/astronomical_objects/pk1d_forest.py)
 
     Attributes
     ----------
-    dec: float (from AstronomicalObject)
-    Declination (in rad)
-
-    healpix: int (from AstronomicalObject)
-    Healpix number associated with (ra, dec)
-
-    los_id: longint (from AstronomicalObject)
-    Line-of-sight id. Same as targetid
-
-    ra: float (from AstronomicalObject)
-    Right ascention (in rad)
-
-    z: float (from AstronomicalObject)
-    Redshift
-
-    bad_continuum_reason: str or None
-    Reason as to why the continuum fit is not acceptable. None for acceptable
-    contiuum.
-
-    continuum: array of float or None (from Forest)
-    Quasar continuum. None for no information
-
-    deltas: array of float or None (from Forest)
-    Flux-transmission field (delta field). None for no information
-
-    flux: array of float (from Forest)
-    Flux
-
-    ivar: array of float (from Forest)
-    Inverse variance
-
-    lambda_: array of float or None (from Forest)
-    Wavelength (in Angstroms)
-
-    log_lambda: array of float or None (from Forest)
-    Logarithm of the wavelength (in Angstroms)
-
-    mean_snr: float (from Forest)
-    Mean signal-to-noise of the forest
-
-    transmission_correction: array of float (from Forest)
-    Transmission correction.
-
-    weights: array of float or None (from Forest)
-    Weights associated to the delta field. None for no information
-
-    night: list of int (from DesiForest)
-    Identifier of the night where the observation was made. None for no info
-
-    petal: list of int (from DesiForest)
-    Identifier of the spectrograph used in the observation. None for no info
-
-    targetid: int (from DesiForest)
-    Targetid of the object
-
-    tile: list of int (from DesiForest)
-    Identifier of the tile used in the observation. None for no info
-
-    exposures_diff: array of floats (from Pk1dForest)
-    Difference between exposures
-
-    mean_z: float
-    Mean redshift of the forest (from Pk1dForest)
-
-    reso: array of floats or None (from Pk1dForest)
-    Resolution of the forest
+    (see DesiForest in py/picca/delta_extraction/astronomical_objects/desi_forest.py)
+    (see Pk1dForest in py/picca/delta_extraction/astronomical_objects/pk1d_forest.py)
     
-    resolution_matrix: array of floats or None (from Pk1dForest)
+    resolution_matrix: 2d-array of floats or None
     Resolution matrix of the forests
     """
     def __init__(self, **kwargs):
diff --git a/py/picca/delta_extraction/astronomical_objects/forest.py b/py/picca/delta_extraction/astronomical_objects/forest.py
index 9e732f7a9..16a050e5c 100644
--- a/py/picca/delta_extraction/astronomical_objects/forest.py
+++ b/py/picca/delta_extraction/astronomical_objects/forest.py
@@ -9,8 +9,7 @@
 from picca.delta_extraction.utils import find_bins
 
 defaults = {
-    "mask fields log": ["flux", "ivar", "transmission_correction", "log_lambda"],
-    "mask fields lin": ["flux", "ivar", "transmission_correction", "lambda_"]
+    "mask fields": ["flux", "ivar", "transmission_correction", "log_lambda"],
 }
 
 class Forest(AstronomicalObject):
@@ -34,24 +33,16 @@ class Forest(AstronomicalObject):
     blinding: str
     Name of the blinding strategy used
 
-    lambda_grid: array of float or None
-    Common grid in lambda_ based on the specified minimum and maximum
-    wavelengths, and delta_lambda.
-
-    lambda_rest_frame_grid: array of float or None (from Forest)
-    Same as lambda_grid but for rest-frame wavelengths
-
     log_lambda_grid: array of float or None
     Common grid in log_lambda based on the specified minimum and maximum
-    wavelengths, and delta_log_lambda.
+    wavelengths, and pixelisation.
 
     log_lambda_rest_frame_grid: array of float or None
     Same as log_lambda_grid but for rest-frame wavelengths.
 
     mask_fields: list of str
     Names of the fields that are affected by masking. In general it will
-    be "flux", "ivar", "transmission_correction" and either "log_lambda" if
-    Forest.wave_solution is "log" or "lambda_" if Forests.wave_solution is "lin",
+    be "flux", "ivar", "transmission_correction", and "log_lambda",
     but some child classes might add more.
 
     wave_solution: "lin" or "log"
@@ -60,22 +51,6 @@ class Forest(AstronomicalObject):
 
     Attributes
     ----------
-    dec: float (from AstronomicalObject)
-    Declination (in rad)
-
-    healpix: int (from AstronomicalObject)
-    Healpix number associated with (ra, dec)
-
-    los_id: longint (from AstronomicalObject)
-    Line-of-sight id. Same as thingid for SDSS data, and same as targetid for
-    DESI data
-
-    ra: float (from AstronomicalObject)
-    Right ascention (in rad)
-
-    z: float (from AstronomicalObject)
-    Redshift
-
     bad_continuum_reason: str or None
     Reason as to why the continuum fit is not acceptable. None for acceptable
     contiuum.
@@ -92,11 +67,10 @@ class Forest(AstronomicalObject):
     ivar: array of float
     Inverse variance
 
-    lambda_: array of float
-    Wavelength (in Angstroms)
-
     log_lambda: array of float or None
-    Logarithm of the wavelength (in Angstroms)
+    Logarithm of the wavelength (in Angstroms). Differs from log_lambda_grid
+    as the particular instance might not have full wavelength coverage or
+    might have some missing pixels (because they are masked)
 
     mean_snr: float
     Mean signal-to-noise of the forest
@@ -108,11 +82,9 @@ class Forest(AstronomicalObject):
     Weights associated to the delta field. None for no information
     """
     blinding = "none"
-    lambda_grid = None
-    lambda_rest_frame_grid = None
     log_lambda_grid = None
     log_lambda_rest_frame_grid = None
-    mask_fields = None
+    mask_fields = []
     wave_solution = None
 
     def __init__(self, **kwargs):
@@ -131,20 +103,11 @@ def __init__(self, **kwargs):
 
         Forest.class_variable_check()
 
-        if Forest.wave_solution == "log":
-            self.lambda_ = None
-            self.log_lambda = kwargs.get("log_lambda")
-            if self.log_lambda is None:
-                raise AstronomicalObjectError("Error constructing Forest. "
-                                              "Missing variable 'log_lambda'")
-            del kwargs["log_lambda"]
-        elif Forest.wave_solution == "lin":
-            self.log_lambda = None
-            self.lambda_ = kwargs.get("lambda")
-            if self.lambda_ is None:
-                raise AstronomicalObjectError("Error constructing Forest. "
-                                              "Missing variable 'lambda'")
-            del kwargs["lambda"]
+        self.log_lambda = kwargs.get("log_lambda")
+        if self.log_lambda is None:
+            raise AstronomicalObjectError("Error constructing Forest. "
+                                          "Missing variable 'log_lambda'")
+        del kwargs["log_lambda"]
 
         self.bad_continuum_reason = None
         self.continuum = kwargs.get("continuum")
@@ -185,9 +148,14 @@ def __init__(self, **kwargs):
     @classmethod
     def class_variable_check(cls):
         """Check that class variables have been correctly initialized"""
-        if cls.wave_solution is None:
+        if cls.log_lambda_grid is None:
             raise AstronomicalObjectError("Error constructing Forest. "
-                                          "Class variable 'wave_solution' "
+                                          "Class variable 'log_lambda_grid' "
+                                          "must be set prior to initialize "
+                                          "instances of this type")
+        if cls.log_lambda_rest_frame_grid is None:
+            raise AstronomicalObjectError("Error constructing Forest. "
+                                          "Class variable 'log_lambda_rest_frame_grid' "
                                           "must be set prior to initialize "
                                           "instances of this type")
         if cls.mask_fields is None:
@@ -196,74 +164,21 @@ def class_variable_check(cls):
                 "Expected list in class variable 'mask fields'. "
                 f"Found {cls.mask_fields}.")
 
-        if cls.wave_solution == "log":
-            if cls.log_lambda_grid is None:
-                raise AstronomicalObjectError("Error constructing Forest. "
-                                              "Class variable 'log_lambda_grid' "
-                                              "must be set prior to initialize "
-                                              "instances of this type")
-            if cls.log_lambda_rest_frame_grid is None:
-                raise AstronomicalObjectError("Error constructing Forest. "
-                                              "Class variable 'log_lambda_rest_frame_grid' "
-                                              "must be set prior to initialize "
-                                              "instances of this type")
-            if cls.lambda_grid is not None:
-                raise AstronomicalObjectError("Error constructing Forest. "
-                                              "Class variable 'lambda_grid' "
-                                              "should be None. Found: "
-                                              f"{cls.lambda_grid}")
-            if cls.lambda_rest_frame_grid is not None:
-                raise AstronomicalObjectError("Error constructing Forest. "
-                                              "Class variable 'lambda_rest_frame_grid' "
-                                              "should be None. Found: "
-                                              f"{cls.lambda_rest_frame_grid}")
-        elif cls.wave_solution == "lin":
-            if cls.log_lambda_grid is not None:
-                raise AstronomicalObjectError("Error constructing Forest. "
-                                              "Class variable 'log_lambda_grid' "
-                                              "should be None. Found: "
-                                              f"{cls.log_lambda_grid}")
-            if cls.log_lambda_rest_frame_grid is not None:
-                raise AstronomicalObjectError("Error constructing Forest. "
-                                              "Class variable 'log_lambda_rest_frame_grid' "
-                                              "should be None. Found: "
-                                              f"{cls.log_lambda_rest_frame_grid}")
-            if cls.lambda_grid is None:
-                raise AstronomicalObjectError("Error constructing Forest. "
-                                              "Class variable 'ambda_grid' "
-                                              "must be set prior to initialize "
-                                              "instances of this type")
-            if cls.lambda_rest_frame_grid is None:
-                raise AstronomicalObjectError("Error constructing Forest. "
-                                              "Class variable 'lambda_rest_frame_grid' "
-                                              "must be set prior to initialize "
-                                              "instances of this type")
-        else:
+        if cls.wave_solution is None:
             raise AstronomicalObjectError("Error constructing Forest. "
                                           "Class variable 'wave_solution' "
-                                          "must be either 'lin' or 'log'. "
-                                          f"Found: {cls.wave_solution}")
+                                          "must be set prior to initialize "
+                                          "instances of this type")
 
     def consistency_check(self):
         """Consistency checks after __init__"""
         if self.flux.size != self.ivar.size:
             raise AstronomicalObjectError("Error constructing Forest. 'flux', "
                                           "and 'ivar' don't have the same size")
-        if Forest.wave_solution == "log":
-            if self.log_lambda.size != self.flux.size:
-                raise AstronomicalObjectError("Error constructing Forest. "
-                                              "'flux'  and 'log_lambda' don't "
-                                              "have the same  size")
-        elif Forest.wave_solution == "lin":
-            if self.lambda_.size != self.flux.size:
-                raise AstronomicalObjectError("Error constructing Forest. "
-                                              "'flux' and 'lambda' don't have "
-                                              "the same size")
-        else:
+        if self.log_lambda.size != self.flux.size:
             raise AstronomicalObjectError("Error constructing Forest. "
-                                          "Class variable 'wave_solution' "
-                                          "must be either 'lin' or 'log'. "
-                                          f"Found: {Forest.wave_solution}")
+                                          "'flux'  and 'log_lambda' don't "
+                                          "have the same  size")
 
     def coadd(self, other):
         """Coadd the information of another forest.
@@ -292,15 +207,7 @@ def coadd(self, other):
                                           f"not happen. this.los_id={self.los_id}, "
                                           f"other.los_id={other.los_id}.")
 
-        if Forest.wave_solution == "log":
-            self.log_lambda = np.append(self.log_lambda, other.log_lambda)
-        elif Forest.wave_solution == "lin":
-            self.lambda_ = np.append(self.lambda_, other.lambda_)
-        else:
-            raise AstronomicalObjectError("Error in coadding Forest. "
-                                          "Class variable 'wave_solution' "
-                                          "must be either 'lin' or 'log'. "
-                                          f"Found: {Forest.wave_solution}")
+        self.log_lambda = np.append(self.log_lambda, other.log_lambda)
         self.flux = np.append(self.flux, other.flux)
         self.ivar = np.append(self.ivar, other.ivar)
         self.transmission_correction = np.append(self.transmission_correction,
@@ -346,11 +253,11 @@ def get_data(self):
             units += ["log Angstrom"]
             array_size = self.log_lambda.size
         elif Forest.wave_solution == "lin":
-            cols += [self.lambda_]
+            cols += [10**self.log_lambda]
             names += ["LAMBDA"]
             comments += ["Lambda"]
             units += ["Angstrom"]
-            array_size = self.lambda_.size
+            array_size = self.log_lambda.size
         else:
             raise AstronomicalObjectError("Error in getting data from Forest. "
                                           "Class variable 'wave_solution' "
@@ -409,9 +316,36 @@ def get_header(self):
                 'name': 'BLINDING',
                 'value': Forest.blinding,
                 'comment': "String specifying the blinding strategy"
-            }
+            },
+            {
+                'name': 'WAVE_SOLUTION',
+                'value': Forest.wave_solution,
+                'comment': "Chosen wavelength solution (linnear or logarithmic)"
+            },
         ]
 
+        if Forest.wave_solution == "log":
+            header += [
+                {
+                    'name': 'DELTA_LOG_LAMBDA',
+                    'value': Forest.log_lambda_grid[1] - Forest.log_lambda_grid[0],
+                    'comment': "Pixel step in log lambda [log(Angstrom)]"
+                },
+            ]
+        elif Forest.wave_solution == "lin":
+            header += [
+                {
+                    'name': 'DELTA_LAMBDA',
+                    'value': 10**Forest.log_lambda_grid[1] - 10**Forest.log_lambda_grid[0],
+                    'comment': "Pixel step in lambda [Angstrom]"
+                },
+            ]
+        else:
+            raise AstronomicalObjectError("Error in Forest.get_header(). "
+                                          "Class variable 'wave_solution' "
+                                          "must be either 'lin' or 'log'. "
+                                          f"Found: {Forest.wave_solution}")
+
         return header
 
     def rebin(self):
@@ -442,65 +376,58 @@ def rebin(self):
         AstronomicalObjectError if ivar only has zeros
         """
         orig_ivar = self.ivar.copy()
+        w1 = np.ones(self.log_lambda.size, dtype=bool)
+        pixel_step = np.nan
+
         # compute bins
         if Forest.wave_solution == "log":
-            delta_log_lambda = Forest.log_lambda_grid[1] - Forest.log_lambda_grid[0]
-            half_delta_log_lambda = delta_log_lambda / 2.
-
-            half_delta_log_lambda_rest_frame = (Forest.log_lambda_rest_frame_grid[1] -
-                                                Forest.log_lambda_rest_frame_grid[0]) / 2.
-
-            w1 = (self.log_lambda >= Forest.log_lambda_grid[0] - half_delta_log_lambda)
-            w1 = w1 & (self.log_lambda < Forest.log_lambda_grid[-1] + half_delta_log_lambda)
-            w1 = w1 & (self.log_lambda - np.log10(1. + self.z) >=
-                       Forest.log_lambda_rest_frame_grid[0] - half_delta_log_lambda_rest_frame)
-            w1 = w1 & (self.log_lambda - np.log10(1. + self.z) <
-                       Forest.log_lambda_rest_frame_grid[-1] + half_delta_log_lambda_rest_frame)
-            w1 = w1 & (self.ivar > 0.)
-            if w1.sum() == 0:
-                self.log_lambda = np.array([])
-                self.flux = np.array([])
-                self.ivar = np.array([])
-                self.transmission_correction = np.array([])
-                return [], [], [], [], []
-            self.log_lambda = self.log_lambda[w1]
-            self.flux = self.flux[w1]
-            self.ivar = self.ivar[w1]
-            self.transmission_correction = self.transmission_correction[w1]
-
-            bins = find_bins(self.log_lambda, Forest.log_lambda_grid)
-            self.log_lambda = Forest.log_lambda_grid[0] + bins * delta_log_lambda
+            pixel_step = Forest.log_lambda_grid[1] - Forest.log_lambda_grid[0]
+            half_pixel_step = pixel_step / 2.
+
+            half_pixel_step_rest_frame = (Forest.log_lambda_rest_frame_grid[1] -
+                                          Forest.log_lambda_rest_frame_grid[0]) / 2.
+
+            w1 &= (self.log_lambda >= Forest.log_lambda_grid[0] - half_pixel_step)
+            w1 &= (self.log_lambda < Forest.log_lambda_grid[-1] + half_pixel_step)
+            w1 &= (self.log_lambda - np.log10(1. + self.z) >=
+                       Forest.log_lambda_rest_frame_grid[0] - half_pixel_step_rest_frame)
+            w1 &= (self.log_lambda - np.log10(1. + self.z) <
+                       Forest.log_lambda_rest_frame_grid[-1] + half_pixel_step_rest_frame)
+            w1 &= (self.ivar > 0.)
 
         elif Forest.wave_solution == "lin":
-            delta_lambda = Forest.lambda_grid[1] - Forest.lambda_grid[0]
-            half_delta_lambda = delta_lambda / 2.
-
-            half_delta_lambda_rest_frame = (Forest.lambda_rest_frame_grid[1] -
-                                            Forest.lambda_rest_frame_grid[0]) / 2.
-            w1 = (self.lambda_ >= Forest.lambda_grid[0] - half_delta_lambda)
-            w1 = w1 & (self.lambda_ < Forest.lambda_grid[-1] + half_delta_lambda)
-            w1 = w1 & (self.lambda_ / (1. + self.z) >= Forest.lambda_rest_frame_grid[0] - half_delta_lambda_rest_frame)
-            w1 = w1 & (self.lambda_ / (1. + self.z) < Forest.lambda_rest_frame_grid[-1] + half_delta_lambda_rest_frame)
-            w1 = w1 & (self.ivar > 0.)
-            if w1.sum() == 0:
-                self.lambda_ = np.array([])
-                self.flux = np.array([])
-                self.ivar = np.array([])
-                self.transmission_correction = np.array([])
-                return [], [], [], [], []
-            self.lambda_ = self.lambda_[w1]
-            self.flux = self.flux[w1]
-            self.ivar = self.ivar[w1]
-            self.transmission_correction = self.transmission_correction[w1]
-
-            bins = find_bins(self.lambda_, Forest.lambda_grid)
-            self.lambda_ = Forest.lambda_grid[0] + bins * delta_lambda
+            pixel_step = 10**Forest.log_lambda_grid[1] - 10**Forest.log_lambda_grid[0]
+            half_pixel_step = pixel_step / 2.
+
+            half_pixel_step_rest_frame = (10**Forest.log_lambda_rest_frame_grid[1] -
+                                          10**Forest.log_lambda_rest_frame_grid[0]) / 2.
+            lambda_ = 10**self.log_lambda
+            w1 &= (lambda_ >= 10**Forest.log_lambda_grid[0] - half_pixel_step)
+            w1 &= (lambda_ < 10**Forest.log_lambda_grid[-1] + half_pixel_step)
+            w1 &= (lambda_ / (1. + self.z) >= 10**Forest.log_lambda_rest_frame_grid[0] - half_pixel_step_rest_frame)
+            w1 &= (lambda_ / (1. + self.z) < 10**Forest.log_lambda_rest_frame_grid[-1] + half_pixel_step_rest_frame)
+            w1 &= (self.ivar > 0.)
         else:
             raise AstronomicalObjectError("Error in rebinning Forest. "
                                           "Class variable 'wave_solution' "
                                           "must be either 'lin' or 'log'. "
                                           f"Found: {Forest.wave_solution}")
 
+        if w1.sum() == 0:
+            self.log_lambda = np.array([])
+            self.flux = np.array([])
+            self.ivar = np.array([])
+            self.transmission_correction = np.array([])
+            return [], [], [], [], []
+        self.log_lambda = self.log_lambda[w1]
+        self.flux = self.flux[w1]
+        self.ivar = self.ivar[w1]
+        self.transmission_correction = self.transmission_correction[w1]
+
+        bins = find_bins(self.log_lambda, Forest.log_lambda_grid,
+                         Forest.wave_solution)
+        self.log_lambda = Forest.log_lambda_grid[0] + bins * pixel_step
+
         # rebin flux, ivar and transmission_correction
         rebin_flux = np.zeros(bins.max() + 1)
         rebin_transmission_correction = np.zeros(bins.max() + 1)
@@ -528,13 +455,13 @@ def rebin(self):
 
         # then rebin wavelength
         if self.wave_solution == "log":
-            rebin_lambda = (Forest.log_lambda_grid[0] +
-                            np.arange(bins.max() + 1) * delta_log_lambda)
-            self.log_lambda = rebin_lambda[w2]
+            rebin_log_lambda = (Forest.log_lambda_grid[0] +
+                            np.arange(bins.max() + 1) * pixel_step)
+            self.log_lambda = rebin_log_lambda[w2]
         elif self.wave_solution == "lin":
-            rebin_lambda = (Forest.lambda_grid[0] +
-                            np.arange(bins.max() + 1) * delta_lambda)
-            self.lambda_ = rebin_lambda[w2]
+            rebin_lambda = (10**Forest.log_lambda_grid[0] +
+                            np.arange(bins.max() + 1) * pixel_step)
+            self.log_lambda = np.log10(rebin_lambda[w2])
         else:
             raise AstronomicalObjectError("wavelength solution must be either "
                                           "'log' or 'linear'")
@@ -551,22 +478,29 @@ def rebin(self):
     def set_class_variables(cls, lambda_min, lambda_max,
                             lambda_min_rest_frame,
                             lambda_max_rest_frame,
-                            pixel_step, pixel_step_rest_frame, wave_solution):
+                            pixel_step, pixel_step_rest_frame,
+                            wave_solution):
         """Set class variables
+
         Arguments
         ---------
         lambda_min: float
         Logarithm of the minimum wavelength (in Angs) to be considered in a forest.
+
         lambda_max: float
         Logarithm of the maximum wavelength (in Angs) to be considered in a forest.
+
         lambda_min_rest_frame: float or None
         As lambda_min but for rest-frame wavelength.
+
         lambda_max_rest_frame: float
         As lambda_max but for rest-frame wavelength.
+
         pixel_step: float
         Wavelength cahnge between two pixels. If pixel_step is "log" this is in
         units of the logarithm of the wavelength (in Angs). If pixel_step is "lin"
         this is in units of the wavelength (in Angs).
+
         wave_solution: "log" or "lin"
         Specifies whether we want to construct a wavelength grid that is evenly
         spaced on wavelength (lin) or on the logarithm of the wavelength (log)
@@ -580,17 +514,15 @@ def set_class_variables(cls, lambda_min, lambda_max,
                 np.log10(lambda_min_rest_frame) + pixel_step_rest_frame/2,
                 np.log10(lambda_max_rest_frame),
                 pixel_step_rest_frame)
-            cls.mask_fields = defaults.get("mask fields log").copy()
         elif wave_solution == "lin":
-            cls.lambda_grid = np.arange(
+            cls.log_lambda_grid = np.log10(np.arange(
                 lambda_min,
                 lambda_max + pixel_step/2,
-                pixel_step)
-            cls.lambda_rest_frame_grid = np.arange(
+                pixel_step))
+            cls.log_lambda_rest_frame_grid = np.log10(np.arange(
                 lambda_min_rest_frame + pixel_step_rest_frame/2,
                 lambda_max_rest_frame,
-                pixel_step_rest_frame)
-            cls.mask_fields = defaults.get("mask fields lin").copy()
+                pixel_step_rest_frame))
         else:
             raise AstronomicalObjectError("Error in setting Forest class "
                                           "variables. 'wave_solution' "
@@ -598,3 +530,5 @@ def set_class_variables(cls, lambda_min, lambda_max,
                                           f"Found: {wave_solution}")
 
         cls.wave_solution = wave_solution
+
+        cls.mask_fields = defaults.get("mask fields").copy()
diff --git a/py/picca/delta_extraction/astronomical_objects/pk1d_forest.py b/py/picca/delta_extraction/astronomical_objects/pk1d_forest.py
index f74f992ba..99923832a 100644
--- a/py/picca/delta_extraction/astronomical_objects/pk1d_forest.py
+++ b/py/picca/delta_extraction/astronomical_objects/pk1d_forest.py
@@ -21,83 +21,14 @@ class Pk1dForest(Forest):
 
     Class Attributes
     ----------------
-    blinding: str (from Forest)
-    Name of the blinding strategy used
-    
-    lambda_grid: array of float or None (from Forest)
-    Common grid in lambda_ based on the specified minimum and maximum
-    wavelengths, and delta_lambda.
-
-    lambda_rest_frame_grid: array of float or None (from Forest)
-    Same as lambda_grid but for rest-frame wavelengths
-
-    log_lambda_grid: array of float or None (from Forest)
-    Common grid in log_lambda based on the specified minimum and maximum
-    wavelengths, and delta_log_lambda.
-
-    log_lambda_rest_frame_grid: array of float or None (from Forest)
-    Same as log_lambda_grid but for rest-frame wavelengths.
-
-    mask_fields: list of str (from Forest)
-    Names of the fields that are affected by masking. In general it will
-    be "flux", "ivar", "transmission_correction", "exposures_diff", "reso" and
-    either "log_lambda" if Forest.wave_solution is "log" or "lambda_" if
-    Forests.wave_solution is "lin", but some child classes might add more.
-
-    wave_solution: "lin" or "log" (from Forest)
-    Determines whether the wavelength solution has linear spacing ("lin") or
-    logarithmic spacing ("log").
+    (see Forest in py/picca/delta_extraction/astronomical_objects/forest.py)
 
     lambda_abs_igm: float
     Wavelength of the IGM absorber
 
     Attributes
     ----------
-    dec: float (from AstronomicalObject)
-    Declination (in rad)
-
-    healpix: int (from AstronomicalObject)
-    Healpix number associated with (ra, dec)
-
-    los_id: longint (from AstronomicalObject)
-    Line-of-sight id. Same as thingid
-
-    ra: float (from AstronomicalObject)
-    Right ascention (in rad)
-
-    z: float (from AstronomicalObject)
-    Redshift
-
-    bad_continuum_reason: str or None (from Forest)
-    Reason as to why the continuum fit is not acceptable. None for acceptable
-    contiuum.
-
-    continuum: array of float or None (from Forest)
-    Quasar continuum. None for no information
-
-    deltas: array of float or None (from Forest)
-    Flux-transmission field (delta field). None for no information
-
-    flux: array of float (from Forest)
-    Flux
-
-    ivar: array of float (from Forest)
-    Inverse variance
-
-    lambda_: array of float (from Forest)
-    Wavelength (in Angstroms)
-
-    log_lambda: array of float or None
-    Logarithm of the wavelength (in Angstroms)
-
-    mean_snr: float (from Forest)
-    Mean signal-to-noise of the forest
-
-    transmission_correction: array of float (from Forest)
-    Transmission correction.
-
-    weights: array of float or None (from Forest)
-    Weights associated to the delta field. None for no information
+    (see Forest in py/picca/delta_extraction/astronomical_objects/forest.py)
 
     exposures_diff: array of floats
     Difference between exposures
@@ -148,21 +79,11 @@ def __init__(self, **kwargs):
 
         # compute mean quality variables
         self.mean_reso = self.reso.mean()
+        self.mean_z = (
+            (np.power(10., self.log_lambda[len(self.log_lambda) - 1]) +
+             np.power(10., self.log_lambda[0])) / 2. /
+            Pk1dForest.lambda_abs_igm - 1.0)
         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]) +
-                 np.power(10., self.log_lambda[0])) / 2. /
-                Pk1dForest.lambda_abs_igm - 1.0)
-        elif Forest.wave_solution == "lin":
-            self.mean_z = (
-                (self.lambda_[len(self.lambda_) - 1] + self.lambda_[0]) / 2. /
-                Pk1dForest.lambda_abs_igm - 1.0)
-        else:
-            raise AstronomicalObjectError("Error in constructing Pk1dForest. "
-                                          "Class variable 'wave_solution' "
-                                          "must be either 'lin' or 'log'. "
-                                          f"Found: '{Forest.wave_solution}'")
 
         self.consistency_check()
 
@@ -349,21 +270,11 @@ def rebin(self):
 
         # finally update control variables
         self.mean_reso = self.reso.mean()
+        self.mean_z = (
+            (np.power(10., self.log_lambda[len(self.log_lambda) - 1]) +
+             np.power(10., self.log_lambda[0])) / 2. /
+            Pk1dForest.lambda_abs_igm - 1.0)
         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]) +
-                 np.power(10., self.log_lambda[0])) / 2. /
-                Pk1dForest.lambda_abs_igm - 1.0)
-        elif Forest.wave_solution == "lin":
-            self.mean_z = (
-                (self.lambda_[len(self.lambda_) - 1] + self.lambda_[0]) / 2. /
-                Pk1dForest.lambda_abs_igm - 1.0)
-        else:
-            raise AstronomicalObjectError("Error in rebinning Pk1dForest. "
-                                          "Class variable 'wave_solution' "
-                                          "must be either 'lin' or 'log'. "
-                                          f"Found: {Forest.wave_solution}")
 
         # return weights and binning solution to be used by child classes if
         # required
diff --git a/py/picca/delta_extraction/astronomical_objects/sdss_forest.py b/py/picca/delta_extraction/astronomical_objects/sdss_forest.py
index 3e078d1a5..27c623812 100644
--- a/py/picca/delta_extraction/astronomical_objects/sdss_forest.py
+++ b/py/picca/delta_extraction/astronomical_objects/sdss_forest.py
@@ -19,78 +19,11 @@ class SdssForest(Forest):
 
     Class Attributes
     ----------------
-    blinding: str (from Forest)
-    Name of the blinding strategy used
-    
-    lambda_grid: array of float or None (from Forest)
-    Common grid in lambda_ based on the specified minimum and maximum
-    wavelengths, and delta_lambda.
-
-    lambda_rest_frame_grid: array of float or None (from Forest)
-    Same as lambda_grid but for rest-frame wavelengths
-
-    log_lambda_grid: array of float or None (from Forest)
-    Common grid in log_lambda based on the specified minimum and maximum
-    wavelengths, and delta_log_lambda.
-
-    log_lambda_rest_frame_grid: array of float or None (from Forest)
-    Same as log_lambda_grid but for rest-frame wavelengths.
-
-    mask_fields: list of str (from Forest)
-    Names of the fields that are affected by masking. In general it will
-    be "flux" and "ivar" but some child classes might add more.
-
-    wave_solution: "lin" or "log" (from Forest)
-    Determines whether the wavelength solution has linear spacing ("lin") or
-    logarithmic spacing ("log").
+    (see Forest in py/picca/delta_extraction/astronomical_objects/forest.py)
 
     Attributes
     ----------
-    dec: float (from AstronomicalObject)
-    Declination (in rad)
-
-    healpix: int (from AstronomicalObject)
-    Healpix number associated with (ra, dec)
-
-    los_id: longint (from AstronomicalObject)
-    Line-of-sight id. Same as thingid
-
-    ra: float (from AstronomicalObject)
-    Right ascention (in rad)
-
-    z: float (from AstronomicalObject)
-    Redshift
-
-    bad_continuum_reason: str or None
-    Reason as to why the continuum fit is not acceptable. None for acceptable
-    contiuum.
-
-    continuum: array of float or None (from Forest)
-    Quasar continuum. None for no information
-
-    deltas: array of float or None (from Forest)
-    Flux-transmission field (delta field). None for no information
-
-    flux: array of float (from Forest)
-    Flux
-
-    ivar: array of float (from Forest)
-    Inverse variance
-
-    lambda_: array of float or None (from Forest)
-    Wavelength (in Angstroms)
-
-    log_lambda: array of float or None (from Forest)
-    Logarithm of the wavelength (in Angstroms)
-
-    mean_snr: float (from Forest)
-    Mean signal-to-noise of the forest
-
-    transmission_correction: array of float (from Forest)
-    Transmission correction.
-
-    weights: array of float or None (from Forest)
-    Weights associated to the delta field. None for no information
+    (see Forest in py/picca/delta_extraction/astronomical_objects/forest.py)
 
     fiberid: list of int
     Fiberid of the observation
diff --git a/py/picca/delta_extraction/astronomical_objects/sdss_pk1d_forest.py b/py/picca/delta_extraction/astronomical_objects/sdss_pk1d_forest.py
index 2baac9546..aadcd7ffb 100644
--- a/py/picca/delta_extraction/astronomical_objects/sdss_pk1d_forest.py
+++ b/py/picca/delta_extraction/astronomical_objects/sdss_pk1d_forest.py
@@ -21,99 +21,13 @@ class SdssPk1dForest(SdssForest, Pk1dForest):
 
     Class Attributes
     ----------------
-    blinding: str (from Forest)
-    Name of the blinding strategy used
-    
-    lambda_grid: array of float or None (from Forest)
-    Common grid in lambda_ based on the specified minimum and maximum
-    wavelengths, and delta_lambda.
-
-    lambda_rest_frame_grid: array of float or None (from Forest)
-    Same as lambda_grid but for rest-frame wavelengths
-
-    log_lambda_grid: array of float or None (from Forest)
-    Common grid in log_lambda based on the specified minimum and maximum
-    wavelengths, and delta_log_lambda.
-
-    log_lambda_rest_frame_grid: array of float or None (from Forest)
-    Same as log_lambda_grid but for rest-frame wavelengths.
-
-    mask_fields: list of str (from Forest)
-    Names of the fields that are affected by masking. In general it will
-    be "flux" and "ivar" but some child classes might add more.
-
-    wave_solution: "lin" or "log" (from Forest)
-    Determines whether the wavelength solution has linear spacing ("lin") or
-    logarithmic spacing ("log").
+    (see SdssForest in py/picca/delta_extraction/astronomical_objects/sdss_forest.py)
+    (see Pk1dForest in py/picca/delta_extraction/astronomical_objects/pk1d_forest.py)
 
     Attributes
     ----------
-    dec: float (from AstronomicalObject)
-    Declination (in rad)
-
-    healpix: int (from AstronomicalObject)
-    Healpix number associated with (ra, dec)
-
-    los_id: longint (from AstronomicalObject)
-    Line-of-sight id. Same as thingid
-
-    ra: float (from AstronomicalObject)
-    Right ascention (in rad)
-
-    z: float (from AstronomicalObject)
-    Redshift
-
-    bad_continuum_reason: str or None
-    Reason as to why the continuum fit is not acceptable. None for acceptable
-    contiuum.
-
-    continuum: array of float or None (from Forest)
-    Quasar continuum. None for no information
-
-    deltas: array of float or None (from Forest)
-    Flux-transmission field (delta field). None for no information
-
-    flux: array of float (from Forest)
-    Flux
-
-    ivar: array of float (from Forest)
-    Inverse variance
-
-    lambda_: array of float or None (from Forest)
-    Wavelength (in Angstroms)
-
-    log_lambda: array of float or None (from Forest)
-    Logarithm of the wavelength (in Angstroms)
-
-    mean_snr: float (from Forest)
-    Mean signal-to-noise of the forest
-
-    transmission_correction: array of float (from Forest)
-    Transmission correction.
-
-    weights: array of float or None (from Forest)
-    Weights associated to the delta field. None for no information
-
-    fiberid: int (from SdssForest)
-    Fiberid of the observation
-
-    mjd: int (from SdssForest)
-    Modified Julian Date of the observation
-
-    plate: int (from SdssForest)
-    Plate of the observation
-
-    thingid: int (from SdssForest)
-    Thingid of the object
-
-    exposures_diff: array of floats (from Pk1dForest)
-    Difference between exposures
-
-    mean_z: float
-    Mean redshift of the forest (from Pk1dForest)
-
-    reso: array of floats or None (from Pk1dForest)
-    Resolution of the forest
+    (see SdssForest in py/picca/delta_extraction/astronomical_objects/sdss_forest.py)
+    (see Pk1dForest in py/picca/delta_extraction/astronomical_objects/pk1d_forest.py)
     """
     def __init__(self, **kwargs):
         """Initialize instance
diff --git a/py/picca/delta_extraction/corrections/calibration_correction.py b/py/picca/delta_extraction/corrections/calibration_correction.py
index f3cca0a4f..e92c85d68 100644
--- a/py/picca/delta_extraction/corrections/calibration_correction.py
+++ b/py/picca/delta_extraction/corrections/calibration_correction.py
@@ -95,20 +95,6 @@ def apply_correction(self, forest):
         CorrectionError if forest instance does not have the attribute
         'log_lambda'
         """
-        if self.wave_solution == "log":
-            if not hasattr(forest, "log_lambda"):
-                raise CorrectionError("Forest instance is missing "
-                                      "attribute 'log_lambda' required by "
-                                      "CalibrationCorrection")
-            correction = self.correct_flux(forest.log_lambda)
-        elif self.wave_solution == "lin":
-            if not hasattr(forest, "lambda_"):
-                raise CorrectionError("Forest instance is missing "
-                                      "attribute 'lambda_' required by "
-                                      "CalibrationCorrection")
-            correction = self.correct_flux(forest.lambda_)
-        else:
-            raise CorrectionError("In CalibrationCorrection wave_solution must "
-                                  "be either 'log' or 'lin'")
+        correction = self.correct_flux(forest.log_lambda)
         forest.flux /= correction
         forest.ivar *= correction**2
diff --git a/py/picca/delta_extraction/corrections/dust_correction.py b/py/picca/delta_extraction/corrections/dust_correction.py
index d5f8c0793..7eddd87ae 100644
--- a/py/picca/delta_extraction/corrections/dust_correction.py
+++ b/py/picca/delta_extraction/corrections/dust_correction.py
@@ -5,7 +5,6 @@
 import numpy as np
 from scipy import interpolate
 
-from picca.delta_extraction.astronomical_objects.forest import Forest
 from picca.delta_extraction.correction import Correction
 from picca.delta_extraction.errors import CorrectionError
 
@@ -45,22 +44,22 @@ def __init__(self, config):
         extinction_conversion_r = config.getfloat("extinction_conversion_r")
         if extinction_conversion_r is None:
             raise CorrectionError("Missing argument 'extinction_conversion_r' "
-                                  "required by SdssDustCorrection")
+                                  "required by DustCorrection")
 
         filename = config.get("filename")
         if filename is None:
             raise CorrectionError("Missing argument 'filename' "
-                                  "required by SdssDustCorrection")
+                                  "required by DustCorrection")
         try:
             hdu = fitsio.read(filename, ext="CATALOG")
             thingid = hdu['THING_ID']
             ext = hdu['EXTINCTION'][:, 1] / extinction_conversion_r
         except OSError:
-            raise CorrectionError("Error loading SdssDustCorrection. "
+            raise CorrectionError("Error loading DustCorrection. "
                                   f"File {filename} does not have extension "
                                   "'CATALOG'")
         except ValueError:
-            raise CorrectionError("Error loading SdssDustCorrection. "
+            raise CorrectionError("Error loading DustCorrection. "
                                   f"File {filename} does not have fields "
                                   "'THING_ID' and/or 'EXTINCTION' in HDU "
                                   "'CATALOG'")
@@ -88,13 +87,7 @@ def apply_correction(self, forest):
         if extinction is None:
             return
 
-        if Forest.wave_solution == "log":
-            correction = unred(10**forest.log_lambda, extinction)
-        elif Forsest.wave_solution == "lin":
-            correction = unred(forest.lambda_, extinction)
-        else:
-            raise CorrectionError("Forest.wave_solution must be either "
-                                  "'log' or 'lin'")
+        correction = unred(10**forest.log_lambda, extinction)
         forest.flux /= correction
         forest.ivar *= correction**2
         if hasattr(forest, "exposures_diff"):
diff --git a/py/picca/delta_extraction/corrections/ivar_correction.py b/py/picca/delta_extraction/corrections/ivar_correction.py
index c14673dbc..cf4bcd0fc 100644
--- a/py/picca/delta_extraction/corrections/ivar_correction.py
+++ b/py/picca/delta_extraction/corrections/ivar_correction.py
@@ -1,4 +1,6 @@
 """This module defines the abstract class IvarCorrection"""
+import logging
+
 import fitsio
 from scipy.interpolate import interp1d
 
@@ -36,6 +38,8 @@ def __init__(self, config):
         CorrectionError if input file does not have fields loglam and/or eta
         in extension VAR_FUNC
         """
+        self.logger = logging.getLogger(__name__)
+
         filename = config.get("filename")
         if filename is None:
             raise CorrectionError("Missing argument 'filename' required by SdssIvarCorrection")
@@ -43,10 +47,11 @@ def __init__(self, config):
             hdu = fitsio.read(filename, ext="VAR_FUNC")
             if "loglam" in hdu.dtype.names:
                 log_lambda = hdu['loglam']
-                self.wave_solution = "log"
             elif "lambda" in hdu.dtype.names:
-                lambda_ = hdu['lambda']
-                self.wave_solution = "lin"
+                self.logger.warning("DeprecationWarning: Reading correction using 'lambda'. "
+                                    "Newer versions of picca always save 'log_lambda' and "
+                                    "so this option will be removed in the future.")
+                log_lambda = np.log10(hdu['lambda'])
             else:
                 raise CorrectionError("Error loading IvarCorrection. In "
                                       "extension 'VAR_FUNC' in file "
@@ -63,20 +68,11 @@ def __init__(self, config):
             raise CorrectionError("Error loading IvarCorrection. "
                                   f"File {filename} does not have fields "
                                   "'loglam' and/or 'eta' in HDU 'VAR_FUNC'")
-        if self.wave_solution == "log":
-            self.correct_ivar = interp1d(log_lambda,
-                                         eta,
-                                         fill_value="extrapolate",
-                                         kind="nearest")
-        elif self.wave_solution == "lin":
-            self.correct_ivar = interp1d(lambda_,
-                                         eta,
-                                         fill_value="extrapolate",
-                                         kind="nearest")
-        else:
-            raise CorrectionError("Forest.wave_solution must be either "
-                                  "'log' or 'lin'")
 
+        self.correct_ivar = interp1d(log_lambda,
+                                     eta,
+                                     fill_value="extrapolate",
+                                     kind="nearest")
 
     def apply_correction(self, forest):
         """Apply the correction. Correction is applied by dividing the
@@ -92,20 +88,5 @@ def apply_correction(self, forest):
         CorrectionError if forest instance does not have the attribute
         'log_lambda'
         """
-        if self.wave_solution == "log":
-            if not hasattr(forest, "log_lambda"):
-                raise CorrectionError("Forest instance is missing "
-                                      "attribute 'log_lambda' required by "
-                                      "IvarCorrection")
-            correction = self.correct_ivar(forest.log_lambda)
-        elif self.wave_solution == "lin":
-            if not hasattr(forest, "lambda_"):
-                raise CorrectionError("Forest instance is missing "
-                                      "attribute 'lambda_' required by "
-                                      "IvarCorrection")
-            correction = self.correct_ivar(forest.lambda_)
-        else:
-            raise CorrectionError("In IvarCorrection wave_solution must "
-                                  "be either 'log' or 'lin'")
-
+        correction = self.correct_ivar(forest.log_lambda)
         forest.ivar /= correction
diff --git a/py/picca/delta_extraction/corrections/optical_depth_correction.py b/py/picca/delta_extraction/corrections/optical_depth_correction.py
index 40268d09a..82bb6bd0d 100644
--- a/py/picca/delta_extraction/corrections/optical_depth_correction.py
+++ b/py/picca/delta_extraction/corrections/optical_depth_correction.py
@@ -3,7 +3,6 @@
 
 import numpy as np
 
-from picca.delta_extraction.astronomical_objects.forest import Forest
 from picca.delta_extraction.correction import Correction
 from picca.delta_extraction.errors import CorrectionError
 from picca.delta_extraction.utils import ABSORBER_IGM
@@ -91,24 +90,12 @@ def apply_correction(self, forest):
         forest: Forest
         A Forest instance to which the correction is applied
         """
-        if Forest.wave_solution == "log":
-            mean_optical_depth = np.ones(forest.log_lambda.size)
-            for tau, gamma, lambda_rest_frame in zip(self.tau_list, self.gamma_list,
-                                                     self.lambda_rest_frame_list):
-
-                w = 10.**forest.log_lambda / (1. + forest.z) <= lambda_rest_frame
-                z = 10.**forest.log_lambda / lambda_rest_frame - 1.
-                mean_optical_depth[w] *= np.exp(-tau * (1. + z[w])**gamma)
-        elif Forest.wave_solution == "log":
-            mean_optical_depth = np.ones(forest.lambda_.size)
-            for tau, gamma, lambda_rest_frame in zip(self.tau_list, self.gamma_list,
-                                                     self.lambda_rest_frame_list):
-
-                w = forest.lambda_ / (1. + forest.z) <= lambda_rest_frame
-                z = forest.lambda_ / lambda_rest_frame - 1.
-                mean_optical_depth[w] *= np.exp(-tau * (1. + z[w])**gamma)
-        else:
-            raise CorrectionError("Forest.wave_solution must be either "
-                                  "'log' or 'lin'")
+        mean_optical_depth = np.ones(forest.log_lambda.size)
+        for tau, gamma, lambda_rest_frame in zip(self.tau_list, self.gamma_list,
+                                                 self.lambda_rest_frame_list):
 
+            w = 10.**forest.log_lambda / (1. + forest.z) <= lambda_rest_frame
+            z = 10.**forest.log_lambda / lambda_rest_frame - 1.
+            mean_optical_depth[w] *= np.exp(-tau * (1. + z[w])**gamma)
+        
         forest.transmission_correction *= mean_optical_depth
diff --git a/py/picca/delta_extraction/data.py b/py/picca/delta_extraction/data.py
index e2ae3fafa..036ab7dc0 100644
--- a/py/picca/delta_extraction/data.py
+++ b/py/picca/delta_extraction/data.py
@@ -148,7 +148,7 @@ def __parse_config(self, config):
             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"
+                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 "
diff --git a/py/picca/delta_extraction/data_catalogues/desi_healpix.py b/py/picca/delta_extraction/data_catalogues/desi_healpix.py
index 2046fdf35..30c3a9942 100644
--- a/py/picca/delta_extraction/data_catalogues/desi_healpix.py
+++ b/py/picca/delta_extraction/data_catalogues/desi_healpix.py
@@ -275,14 +275,8 @@ def read_file(self, filename, catalogue, forests_by_targetid):
                     "dec": row['DEC'],
                     "z": row['Z'],
                 }
-                if Forest.wave_solution == "log":
-                    args["log_lambda"] = np.log10(spec['WAVELENGTH'])
-                elif Forest.wave_solution == "lin":
-                    args["lambda"] = spec['WAVELENGTH']
-                else:
-                    raise DataError("Forest.wave_solution must be either "
-                                    "'log' or 'lin'")
-
+                args["log_lambda"] = np.log10(spec['WAVELENGTH'])
+                
                 if self.analysis_type == "BAO 3D":
                     forest = DesiForest(**args)
                 elif self.analysis_type == "PK 1D":
diff --git a/py/picca/delta_extraction/data_catalogues/desi_tile.py b/py/picca/delta_extraction/data_catalogues/desi_tile.py
index b28eade9e..accde9c4d 100644
--- a/py/picca/delta_extraction/data_catalogues/desi_tile.py
+++ b/py/picca/delta_extraction/data_catalogues/desi_tile.py
@@ -314,14 +314,8 @@ def read_data(self):
                         "tile": entry["TILEID"],
                         "night": entry["NIGHT"],
                     }
-                    if Forest.wave_solution == "log":
-                        args["log_lambda"] = np.log10(spec['WAVELENGTH'])
-                    elif Forest.wave_solution == "lin":
-                        args["lambda"] = spec['WAVELENGTH']
-                    else:
-                        raise DataError("Forest.wave_solution must be either "
-                                        "'log' or 'lin'")
-
+                    args["log_lambda"] = np.log10(spec['WAVELENGTH'])
+                    
                     if self.analysis_type == "BAO 3D":
                         forest = DesiForest(**args)
                     elif self.analysis_type == "PK 1D":
diff --git a/py/picca/delta_extraction/data_catalogues/desisim_mocks.py b/py/picca/delta_extraction/data_catalogues/desisim_mocks.py
index b0efaadec..9b97c63f7 100644
--- a/py/picca/delta_extraction/data_catalogues/desisim_mocks.py
+++ b/py/picca/delta_extraction/data_catalogues/desisim_mocks.py
@@ -60,6 +60,7 @@ def __init__(self, config):
 
         self.logger = logging.getLogger(__name__)
         super().__init__(config)
+        
     def read_data(self):
         """Read the spectra and formats its data as Forest instances.
 
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 c48e755ec..eb46f562f 100644
--- a/py/picca/delta_extraction/expected_fluxes/dr16_expected_flux.py
+++ b/py/picca/delta_extraction/expected_fluxes/dr16_expected_flux.py
@@ -9,7 +9,7 @@
 
 from picca.delta_extraction.astronomical_objects.forest import Forest
 from picca.delta_extraction.astronomical_objects.pk1d_forest import Pk1dForest
-from picca.delta_extraction.errors import ExpectedFluxError
+from picca.delta_extraction.errors import ExpectedFluxError, AstronomicalObjectError
 from picca.delta_extraction.expected_flux import ExpectedFlux
 from picca.delta_extraction.utils import find_bins
 
@@ -39,31 +39,22 @@ class Dr16ExpectedFlux(ExpectedFlux):
     -------
     extract_deltas (from ExpectedFlux)
     __init__
-    _initialize_arrays_lin
-    _initialize_arrays_log
+    _initialize_variables
     __parse_config
     compute_continuum
-        get_cont_model
-        chi2
     compute_delta_stack
-    compute_mean_cont_lin
-    compute_mean_cont_log
+    compute_mean_cont
     compute_expected_flux
     compute_var_stats
         chi2
+    get_continuum_model
+    get_continuum_weights
+    populate_los_ids
+    save_iteration_step
 
     Attributes
     ----------
-    los_ids: dict (from ExpectedFlux)
-    A dictionary to store the mean expected flux fraction, the weights, and
-    the inverse variance for each line of sight. Keys are the identifier for the
-    line of sight and values are dictionaries with the keys "mean expected flux",
-    and "weights" pointing to the respective arrays. If the given Forests are
-    also Pk1dForests, then the key "ivar" must be available. Arrays have the same
-    size as the flux array for the corresponding line of sight forest instance.
-
-    out_dir: str (from ExpectedFlux)
-    Directory where logs will be saved.
+    (see ExpectedFlux in py/picca/delta_extraction/expected_flux.py)
 
     continuum_fit_parameters: dict
     A dictionary containing the continuum fit parameters for each line of sight.
@@ -114,9 +105,9 @@ class Dr16ExpectedFlux(ExpectedFlux):
     limit_var_lss: tuple of floats
     Limits on the pixel variance due to Large Scale Structure
 
-    log_lambda: array of float or None
-    Logarithm of the rest frame wavelengths where the variance functions and
-    statistics are computed. None (and unused) for a linear wavelength solution.
+    log_lambda_var_func_grid: array of float
+    Logarithm of the wavelengths where the variance functions and
+    statistics are computed.
 
     num_bins_variance: int
     Number of bins to be used to compute variance functions and statistics as
@@ -147,10 +138,6 @@ def __init__(self, config):
         ---------
         config: configparser.SectionProxy
         Parsed options to initialize class
-
-        Raise
-        -----
-        ExpectedFluxError if Forest.wave_solution is not 'lin' or 'log'
         """
         self.logger = logging.getLogger(__name__)
         super().__init__(config)
@@ -167,6 +154,7 @@ def __init__(self, config):
         self.use_ivar_as_weight = None
         self.__parse_config(config)
 
+        # initialize variables
         self.get_eta = None
         self.get_fudge = None
         self.get_mean_cont = None
@@ -174,100 +162,16 @@ def __init__(self, config):
         self.get_num_pixels = None
         self.get_valid_fit = None
         self.get_var_lss = None
-        self.lambda_ = None
-        self.log_lambda = None
-        if Forest.wave_solution == "log":
-            self._initialize_variables_log()
-        elif Forest.wave_solution == "lin":
-            self._initialize_variables_lin()
-        else:
-            raise ExpectedFluxError("Forest.wave_solution must be either "
-                                    "'log' or 'lin'")
+        self.log_lambda_var_func_grid = None
+        self._initialize_variables()
 
         self.continuum_fit_parameters = None
 
         self.get_stack_delta = None
         self.get_stack_delta_weights = None
 
-    def _initialize_variables_lin(self):
-        """Initialize useful variables assuming a linear wavelength solution.
-        The initialized arrays are:
-        - self.get_eta
-        - self.get_fudge
-        - self.get_mean_cont
-        - self.get_num_pixels
-        - self.get_var_lss
-        - self.lambda_
-        """
-        # initialize the mean quasar continuum
-        # TODO: maybe we can drop this and compute first the mean quasar
-        # continuum on compute_mean_expected_flux
-        self.get_mean_cont = interp1d(Forest.lambda_rest_frame_grid,
-                                      np.ones_like(Forest.lambda_rest_frame_grid),
-                                      fill_value="extrapolate")
-        self.get_mean_cont_weight = interp1d(Forest.lambda_rest_frame_grid,
-                                             np.zeros_like(
-                                                 Forest.lambda_rest_frame_grid),
-                                             fill_value="extrapolate")
-
-        # initialize the variance-related variables (see equation 4 of
-        # du Mas des Bourboux et al. 2020 for details on these variables)
-        self.lambda_ = (
-            Forest.lambda_grid[0] + (np.arange(self.num_bins_variance) + .5) *
-            (Forest.lambda_grid[-1] - Forest.lambda_grid[0]) /
-            self.num_bins_variance)
-        # if use_ivar_as_weight is set, eta, var_lss and fudge will be ignored
-        # print a message to inform the user
-        if self.use_ivar_as_weight:
-            self.logger.info(("using ivar as weights, ignoring eta, "
-                              "var_lss, fudge fits"))
-            eta = np.ones(self.num_bins_variance)
-            var_lss = np.zeros(self.num_bins_variance)
-            fudge = np.zeros(self.num_bins_variance)
-            num_pixels = np.zeros(self.num_bins_variance)
-            valid_fit = np.ones(self.num_bins_variance)
-        # if use_constant_weight is set then initialize eta, var_lss, and fudge
-        # with values to have constant weights
-        elif self.use_constant_weight:
-            self.logger.info(("using constant weights, ignoring eta, "
-                              "var_lss, fudge fits"))
-            eta = np.zeros(self.num_bins_variance)
-            var_lss = np.ones(self.num_bins_variance)
-            fudge = np.zeros(self.num_bins_variance)
-            num_pixels = np.zeros(self.num_bins_variance)
-            valid_fit = np.ones(self.num_bins_variance)
-        # normal initialization: eta, var_lss, and fudge are ignored in the
-        # first iteration
-        else:
-            eta = np.ones(self.num_bins_variance)
-            var_lss = np.zeros(self.num_bins_variance)
-            fudge = np.zeros(self.num_bins_variance)
-            num_pixels = np.zeros(self.num_bins_variance)
-            valid_fit = np.zeros(self.num_bins_variance)
-
-        self.get_eta = interp1d(self.lambda_,
-                                eta,
-                                fill_value='extrapolate',
-                                kind='nearest')
-        self.get_var_lss = interp1d(self.lambda_,
-                                    var_lss,
-                                    fill_value='extrapolate',
-                                    kind='nearest')
-        self.get_fudge = interp1d(self.lambda_,
-                                  fudge,
-                                  fill_value='extrapolate',
-                                  kind='nearest')
-        self.get_num_pixels = interp1d(self.lambda_,
-                                       num_pixels,
-                                       fill_value="extrapolate",
-                                       kind='nearest')
-        self.get_valid_fit = interp1d(self.lambda_,
-                                      valid_fit,
-                                      fill_value="extrapolate",
-                                      kind='nearest')
-
-    def _initialize_variables_log(self):
-        """Initialize useful variables assuming a log-linear wavelength solution.
+    def _initialize_variables(self):
+        """Initialize useful variables
         The initialized arrays are:
         - self.get_eta
         - self.get_fudge
@@ -275,11 +179,22 @@ def _initialize_variables_log(self):
         - self.get_num_pixels
         - self.get_valid_fit
         - self.get_var_lss
-        - self.log_lambda
+        - self.log_lambda_var_func_grid
+
+        Raise
+        -----
+        ExpectedFluxError if Forest class variables are not set
         """
+        # check that Forest class variables are set
+        try:
+            Forest.class_variable_check()
+        except AstronomicalObjectError:
+            raise ExpectedFluxError("Forest class variables need to be set "
+                                    "before initializing variables here." )
+
         # initialize the mean quasar continuum
         # TODO: maybe we can drop this and compute first the mean quasar
-        # continuum on compute_mean_expected_flux
+        # continuum on compute_expected_flux
         self.get_mean_cont = interp1d(Forest.log_lambda_rest_frame_grid,
                                       np.ones_like(Forest.log_lambda_rest_frame_grid),
                                       fill_value="extrapolate")
@@ -291,10 +206,29 @@ def _initialize_variables_log(self):
 
         # initialize the variance-related variables (see equation 4 of
         # du Mas des Bourboux et al. 2020 for details on these variables)
-        self.log_lambda = (Forest.log_lambda_grid[0] +
-                           (np.arange(self.num_bins_variance) + .5) *
-                           (Forest.log_lambda_grid[-1] - Forest.log_lambda_grid[0]) /
-                           self.num_bins_variance)
+        if Forest.wave_solution == "log":
+            self.log_lambda_var_func_grid = (
+                Forest.log_lambda_grid[0] + (np.arange(self.num_bins_variance) + .5) *
+                (Forest.log_lambda_grid[-1] - Forest.log_lambda_grid[0]) /
+                self.num_bins_variance)
+        # TODO: this is related with the todo in check the effect of finding
+        # the nearest bin in log_lambda space versus lambda space infunction
+        # find_bins in utils.py. Once we understand that we can remove
+        # the dependence from Forest from here too.
+        elif Forest.wave_solution == "lin":
+            self.log_lambda_var_func_grid = np.log10(
+                10**Forest.log_lambda_grid[0] + (np.arange(self.num_bins_variance) + .5) *
+                (10**Forest.log_lambda_grid[-1] - 10**Forest.log_lambda_grid[0]) /
+                self.num_bins_variance)
+
+        # TODO: Replace the if/else block above by something like the commented
+        # block below. We need to check the impact of doing this on the final
+        # deltas first (eta, var_lss and fudge will be differently sampled).
+        #start of commented block
+        #resize = len(Forest.log_lambda_grid)/self.num_bins_variance
+        #print(resize)
+        #self.log_lambda_var_func_grid = Forest.log_lambda_grid[::int(resize)]
+        #end of commented block
 
         # if use_ivar_as_weight is set, eta, var_lss and fudge will be ignored
         # print a message to inform the user
@@ -325,24 +259,24 @@ def _initialize_variables_log(self):
             num_pixels = np.zeros(self.num_bins_variance)
             valid_fit = np.zeros(self.num_bins_variance, dtype=bool)
 
-        self.get_eta = interp1d(self.log_lambda,
+        self.get_eta = interp1d(self.log_lambda_var_func_grid,
                                 eta,
                                 fill_value='extrapolate',
                                 kind='nearest')
-        self.get_var_lss = interp1d(self.log_lambda,
+        self.get_var_lss = interp1d(self.log_lambda_var_func_grid,
                                     var_lss,
                                     fill_value='extrapolate',
                                     kind='nearest')
-        self.get_fudge = interp1d(self.log_lambda,
+        self.get_fudge = interp1d(self.log_lambda_var_func_grid,
                                   fudge,
                                   fill_value='extrapolate',
                                   kind='nearest')
-        self.get_num_pixels = interp1d(self.log_lambda,
+        self.get_num_pixels = interp1d(self.log_lambda_var_func_grid,
                                        num_pixels,
                                        fill_value="extrapolate",
                                        kind='nearest')
-        self.get_valid_fit = interp1d(self.log_lambda,
-                                      num_pixels,
+        self.get_valid_fit = interp1d(self.log_lambda_var_func_grid,
+                                      valid_fit,
                                       fill_value="extrapolate",
                                       kind='nearest')
 
@@ -428,7 +362,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.
@@ -446,40 +380,24 @@ def compute_continuum(self, forest):
         ------
         forest: Forest
         The modified forest instance
-
-        Raise
-        -----
-        ExpectedFluxError if Forest.wave_solution is not 'lin' or 'log'
         """
         self.continuum_fit_parameters = {}
 
-        if Forest.wave_solution == "log":
-            # get mean continuum
-            mean_cont = self.get_mean_cont(forest.log_lambda -
-                                           np.log10(1 + forest.z))
-        elif Forest.wave_solution == "lin":
-            # get mean continuum
-            mean_cont = self.get_mean_cont(forest.lambda_ / (1 + forest.z))
-        else:
-            raise ExpectedFluxError("Forest.wave_solution must be either "
-                                    "'log' or 'lin'")
+        # get mean continuum
+        mean_cont = self.get_mean_cont(forest.log_lambda -
+                                       np.log10(1 + forest.z))
 
         # add transmission correction
         # (previously computed using method add_optical_depth)
         mean_cont *= forest.transmission_correction
-        mean_cont_kwargs = {"mean_cont": mean_cont}
-
-        if Forest.wave_solution == "log":
-            mean_cont_kwargs["log_lambda_max"] = (
-                Forest.log_lambda_rest_frame_grid[-1] + np.log10(1 + forest.z))
-            mean_cont_kwargs["log_lambda_min"] = (
-                Forest.log_lambda_rest_frame_grid[0] + np.log10(1 + forest.z))
 
-        elif Forest.wave_solution == "lin":
-            mean_cont_kwargs["lambda_max"] = (
-                Forest.lambda_rest_frame_grid[-1] * (1 + forest.z))
-            mean_cont_kwargs["lambda_min"] = (
-                Forest.lambda_rest_frame_grid[0] * (1 + forest.z))
+        mean_cont_kwargs = {"mean_cont": mean_cont}
+        # TODO: This can probably be replaced by forest.log_lambda[-1] and
+        # forest.log_lambda[0]
+        mean_cont_kwargs["log_lambda_max"] = (
+            Forest.log_lambda_rest_frame_grid[-1] + np.log10(1 + forest.z))
+        mean_cont_kwargs["log_lambda_min"] = (
+            Forest.log_lambda_rest_frame_grid[0] + np.log10(1 + forest.z))
 
         leasts_squares = LeastsSquaresContModel(
             forest=forest,
@@ -500,6 +418,7 @@ def compute_continuum(self, forest):
         minimizer.fixed["bq"] = self.order == 0
         minimizer.migrad()
 
+
         forest.bad_continuum_reason = None
         temp_cont_model = self.get_continuum_model(forest,
                                                    minimizer.values["aq"],
@@ -532,22 +451,10 @@ def compute_delta_stack(self, forests, stack_from_deltas=False):
 
         stack_from_deltas: bool - default: False
         Flag to determine whether to stack from deltas or compute them
-
-        Raise
-        -----
-        ExpectedFluxError if Forest.wave_solution is not 'lin' or 'log'
         """
-        # TODO: move this to _initialize_variables_lin and
-        # _initialize_variables_log (after tests are done)
-        if Forest.wave_solution == "log":
-            stack_delta = np.zeros_like(Forest.log_lambda_grid)
-            stack_weight = np.zeros_like(Forest.log_lambda_grid)
-        elif Forest.wave_solution == "lin":
-            stack_delta = np.zeros_like(Forest.lambda_grid)
-            stack_weight = np.zeros_like(Forest.lambda_grid)
-        else:
-            raise ExpectedFluxError("Forest.wave_solution must be either "
-                                    "'log' or 'linear'")
+        # TODO: move this to _initialize_variables (after tests are done)
+        stack_delta = np.zeros_like(Forest.log_lambda_grid)
+        stack_weight = np.zeros_like(Forest.log_lambda_grid)
 
         for forest in forests:
             if stack_from_deltas:
@@ -558,30 +465,15 @@ def compute_delta_stack(self, forests, stack_from_deltas=False):
                 if forest.continuum is None:
                     continue
                 delta = forest.flux / forest.continuum
-                if Forest.wave_solution == "log":
-                    var_lss = self.get_var_lss(forest.log_lambda)
-                    eta = self.get_eta(forest.log_lambda)
-                    fudge = self.get_fudge(forest.log_lambda)
-                elif Forest.wave_solution == "lin":
-                    var_lss = self.get_var_lss(forest.lambda_)
-                    eta = self.get_eta(forest.lambda_)
-                    fudge = self.get_fudge(forest.lambda_)
-                else:
-                    raise ExpectedFluxError(
-                        "Forest.wave_solution must be either "
-                        "'log' or 'linear'")
+                var_lss = self.get_var_lss(forest.log_lambda)
+                eta = self.get_eta(forest.log_lambda)
+                fudge = self.get_fudge(forest.log_lambda)
                 var = 1. / forest.ivar / forest.continuum**2
                 variance = eta * var + var_lss + fudge / var
                 weights = 1. / variance
 
-            if Forest.wave_solution == "log":
-                bins = find_bins(forest.log_lambda, Forest.log_lambda_grid)
-
-            elif Forest.wave_solution == "lin":
-                bins = find_bins(forest.lambda_, Forest.lambda_grid)
-            else:
-                raise ExpectedFluxError("Forest.wave_solution must be either "
-                                        "'log' or 'linear'")
+            bins = find_bins(forest.log_lambda, Forest.log_lambda_grid,
+                             Forest.wave_solution)
             rebin = np.bincount(bins, weights=delta * weights)
             stack_delta[:len(rebin)] += rebin
             rebin = np.bincount(bins, weights=weights)
@@ -590,89 +482,20 @@ def compute_delta_stack(self, forests, stack_from_deltas=False):
         w = stack_weight > 0
         stack_delta[w] /= stack_weight[w]
 
-        if Forest.wave_solution == "log":
-            self.get_stack_delta = interp1d(Forest.log_lambda_grid[stack_weight > 0.],
-                                            stack_delta[stack_weight > 0.],
-                                            kind="nearest",
-                                            fill_value="extrapolate")
-            self.get_stack_delta_weights = interp1d(
-                Forest.log_lambda_grid[stack_weight > 0.],
-                stack_weight[stack_weight > 0.],
-                kind="nearest",
-                fill_value=0.0,
-                bounds_error=False)
-        elif Forest.wave_solution == "lin":
-            self.get_stack_delta = interp1d(Forest.lambda_grid[stack_weight > 0.],
-                                            stack_delta[stack_weight > 0.],
-                                            kind="nearest",
-                                            fill_value="extrapolate")
-            self.get_stack_delta_weights = interp1d(
-                Forest.lambda_grid[stack_weight > 0.],
-                stack_weight[stack_weight > 0.],
-                kind="nearest",
-                fill_value=0.0,
-                bounds_error=False)
-        else:
-            raise ExpectedFluxError("Forest.wave_solution must be either "
-                                    "'log' or 'linear'")
-
-    def compute_mean_cont_lin(self, forests):
-        """Compute the mean quasar continuum over the whole sample assuming a
-        linear wavelength solution. Then updates the value of self.get_mean_cont
-        to contain it
-
-        Arguments
-        ---------
-        forests: List of Forest
-        A list of Forest from which to compute the deltas.
-        """
-        mean_cont = np.zeros_like(Forest.lambda_rest_frame_grid)
-        mean_cont_weight = np.zeros_like(Forest.lambda_rest_frame_grid)
-
-        # first compute <F/C> in bins. C=Cont_old*spectrum_dependent_fitting_fct
-        # (and Cont_old is constant for all spectra in a bin), thus we actually
-        # compute
-        #    1/Cont_old * <F/spectrum_dependent_fitting_function>
-        for forest in forests:
-            if forest.bad_continuum_reason is not None:
-                continue
-            bins = find_bins(
-                forest.lambda_ / (1 + forest.z),
-                Forest.lambda_rest_frame_grid
-            )
-
-            var_lss = self.get_var_lss(forest.lambda_)
-            eta = self.get_eta(forest.lambda_)
-            fudge = self.get_fudge(forest.lambda_)
-            var_pipe = 1. / forest.ivar / forest.continuum**2
-            variance = eta * var_pipe + var_lss + fudge / var_pipe
-            weights = 1 / variance
-            cont = np.bincount(bins,
-                               weights=forest.flux / forest.continuum * weights)
-            mean_cont[:len(cont)] += cont
-            cont = np.bincount(bins, weights=weights)
-            mean_cont_weight[:len(cont)] += cont
-
-        w = mean_cont_weight > 0
-        mean_cont[w] /= mean_cont_weight[w]
-        mean_cont /= mean_cont.mean()
-        lambda_cont = Forest.lambda_rest_frame_grid[w]
-
-        # the new mean continuum is multiplied by the previous one to recover
-        # <F/spectrum_dependent_fitting_function>
-        new_cont = self.get_mean_cont(lambda_cont) * mean_cont[w]
-        self.get_mean_cont = interp1d(lambda_cont,
-                                      new_cont,
-                                      fill_value="extrapolate")
-        self.get_mean_cont_weight = interp1d(lambda_cont,
-                                             mean_cont_weight[w],
-                                             fill_value=0.0,
-                                             bounds_error=False)
-
-    def compute_mean_cont_log(self, forests):
-        """Compute the mean quasar continuum over the whole sample assuming a
-        log-linear wavelength solution. Then updates the value of
-        self.get_mean_cont to contain it
+        self.get_stack_delta = interp1d(Forest.log_lambda_grid[stack_weight > 0.],
+                                        stack_delta[stack_weight > 0.],
+                                        kind="nearest",
+                                        fill_value="extrapolate")
+        self.get_stack_delta_weights = interp1d(
+            Forest.log_lambda_grid[stack_weight > 0.],
+            stack_weight[stack_weight > 0.],
+            kind="nearest",
+            fill_value=0.0,
+            bounds_error=False)
+
+    def compute_mean_cont(self, forests):
+        """Compute the mean quasar continuum over the whole sample.
+        Then updates the value of self.get_mean_cont to contain it
 
         Arguments
         ---------
@@ -691,7 +514,8 @@ def compute_mean_cont_log(self, forests):
                 continue
             bins = find_bins(
                 forest.log_lambda - np.log10(1 + forest.z),
-                Forest.log_lambda_rest_frame_grid
+                Forest.log_lambda_rest_frame_grid,
+                Forest.wave_solution
             )
 
             var_lss = self.get_var_lss(forest.log_lambda)
@@ -732,10 +556,6 @@ def compute_expected_flux(self, forests):
         ---------
         forests: List of Forest
         A list of Forest from which to compute the deltas.
-
-        Raise
-        -----
-        ExpectedFluxError if Forest.wave_solution is not 'lin' or 'log'
         """
         context = multiprocessing.get_context('fork')
         for iteration in range(self.num_iterations):
@@ -751,13 +571,7 @@ def compute_expected_flux(self, forests):
 
             if iteration < self.num_iterations - 1:
                 # Compute mean continuum (stack in rest-frame)
-                if Forest.wave_solution == "log":
-                    self.compute_mean_cont_log(forests)
-                elif Forest.wave_solution == "lin":
-                    self.compute_mean_cont_lin(forests)
-                else:
-                    raise ExpectedFluxError("Forest.wave_solution must be "
-                                            "either 'log' or 'linear'")
+                self.compute_mean_cont(forests)
 
                 # Compute observer-frame mean quantities (var_lss, eta, fudge)
                 if not (self.use_ivar_as_weight or self.use_constant_weight):
@@ -834,33 +648,18 @@ def compute_var_stats(self, forests):
 
             # select the pipeline variance bins
             var_pipe_bins = np.floor(
-                (np.log10(var_pipe) - var_pipe_min) /
+                (np.log10(var_pipe[w]) - var_pipe_min) /
                 (var_pipe_max - var_pipe_min) * num_var_bins).astype(int)
-            # filter the values with a pipeline variance out of range
-            var_pipe_bins = var_pipe_bins[w]
 
             # select the wavelength bins
-            if Forest.wave_solution == "log":
-                log_lambda_bins = find_bins(
-                    forest.log_lambda,
-                    self.log_lambda
-                    )
-                # filter the values with a pipeline variance out of range
-                log_lambda_bins = log_lambda_bins[w]
-                # compute overall bin
-                bins = var_pipe_bins + num_var_bins * log_lambda_bins
-            elif Forest.wave_solution == "lin":
-                lambda_bins = find_bins(
-                    forest.lambda_,
-                    self.lambda_,
-                    )
-                # filter the values with a pipeline variance out of range
-                lambda_bins = lambda_bins[w]
-                # compute overall bin
-                bins = var_pipe_bins + num_var_bins * lambda_bins
-            else:
-                raise ExpectedFluxError("Forest.wave_solution must be either "
-                                        "'log' or 'linear'")
+            log_lambda_bins = find_bins(
+                forest.log_lambda[w],
+                self.log_lambda_var_func_grid,
+                Forest.wave_solution
+                )
+
+            # compute overall bin
+            bins = var_pipe_bins + num_var_bins * log_lambda_bins
 
             # compute deltas
             delta = (forest.flux / forest.continuum - 1)
@@ -894,15 +693,8 @@ def compute_var_stats(self, forests):
         fudge_ref = 1e-7
 
         self.logger.progress(" Mean quantities in observer-frame")
-        if Forest.wave_solution == "log":
-            self.logger.progress(
-                " loglam    eta      var_lss  fudge    chi2     num_pix ")
-        elif Forest.wave_solution == "lin":
-            self.logger.progress(
-                " lam    eta      var_lss  fudge    chi2     num_pix ")
-        else:
-            raise ExpectedFluxError("Forest.wave_solution must be either "
-                                    "'log' or 'linear'")
+        self.logger.progress(
+            " loglam    eta      var_lss  fudge    chi2     num_pix valid_fit")
         for index in range(self.num_bins_variance):
             # pylint: disable-msg=cell-var-from-loop
             # this function is defined differntly at each step of the loop
@@ -995,64 +787,33 @@ def chi2(eta, var_lss, fudge):
                                       num_var_bins].sum()
             chi2_in_bin[index] = minimizer.fval
 
-            if Forest.wave_solution == "log":
-                self.logger.progress(
-                    f" {self.log_lambda[index]:.3e} "
-                    f"{eta[index]:.2e} {var_lss[index]:.2e} {fudge[index]:.2e} "
-                    + f"{chi2_in_bin[index]:.2e} {num_pixels[index]:.2e} ")
-            elif Forest.wave_solution == "lin":
-                self.logger.progress(
-                    f" {self.lambda_[index]:.3e} "
-                    f"{eta[index]:.2e} {var_lss[index]:.2e} {fudge[index]:.2e} "
-                    + f"{chi2_in_bin[index]:.2e} {num_pixels[index]:.2e} ")
+            self.logger.progress(
+                f" {self.log_lambda_var_func_grid[index]:.3e} "
+                f"{eta[index]:.2e} {var_lss[index]:.2e} {fudge[index]:.2e} "
+                + f"{chi2_in_bin[index]:.2e} {num_pixels[index]:.2e} {valid_fit[index]}")
 
         w = num_pixels > 0
 
-        if Forest.wave_solution == "log":
-            self.get_eta = interp1d(self.log_lambda[w],
-                                    eta[w],
+        self.get_eta = interp1d(self.log_lambda_var_func_grid[w],
+                                eta[w],
+                                fill_value="extrapolate",
+                                kind="nearest")
+        self.get_var_lss = interp1d(self.log_lambda_var_func_grid[w],
+                                    var_lss[w],
                                     fill_value="extrapolate",
                                     kind="nearest")
-            self.get_var_lss = interp1d(self.log_lambda[w],
-                                        var_lss[w],
-                                        fill_value="extrapolate",
-                                        kind="nearest")
-            self.get_fudge = interp1d(self.log_lambda[w],
-                                      fudge[w],
-                                      fill_value="extrapolate",
-                                      kind="nearest")
-            self.get_num_pixels = interp1d(self.log_lambda[w],
-                                      num_pixels[w],
-                                      fill_value="extrapolate",
-                                      kind="nearest")
-            self.get_valid_fit = interp1d(self.log_lambda[w],
+        self.get_fudge = interp1d(self.log_lambda_var_func_grid[w],
+                                  fudge[w],
+                                  fill_value="extrapolate",
+                                  kind="nearest")
+        self.get_num_pixels = interp1d(self.log_lambda_var_func_grid[w],
+                                       num_pixels[w],
+                                       fill_value="extrapolate",
+                                       kind="nearest")
+        self.get_valid_fit = interp1d(self.log_lambda_var_func_grid[w],
                                       valid_fit[w],
                                       fill_value="extrapolate",
                                       kind="nearest")
-        elif Forest.wave_solution == "lin":
-            self.get_eta = interp1d(self.lambda_[w],
-                                    eta[w],
-                                    fill_value="extrapolate",
-                                    kind="nearest")
-            self.get_var_lss = interp1d(self.lambda_[w],
-                                        var_lss[w],
-                                        fill_value="extrapolate",
-                                        kind="nearest")
-            self.get_fudge = interp1d(self.lambda_[w],
-                                      fudge[w],
-                                      fill_value="extrapolate",
-                                      kind="nearest")
-            self.get_num_pixels = interp1d(self.lambda_[w],
-                                      num_pixels[w],
-                                      fill_value="extrapolate",
-                                      kind="nearest")
-            self.get_valid_fit = interp1d(self.lambda_[w],
-                                      valid_fit[w],
-                                      fill_value=0,
-                                      kind="nearest")
-        else:
-            raise ExpectedFluxError("Forest.wave_solution must be either "
-                                    "'log' or 'linear'")
 
     # pylint: dissable=no-self-use
     # We expect this function to be changed by some child classes
@@ -1075,21 +836,11 @@ def get_continuum_model(self, forest, aq, bq, **kwargs):
         mean_cont: array of floats
         Mean continuum. Required.
 
-        lambda_max: float
-        Maximum lambda_ for this forest. Required only if Forest.wave_solution
-        is lin.
-
-        lambda_min: float
-        Minimum lambda_ for this forest. Required only if Forest.wave_solution
-        is lin.
-
         log_lambda_max: float
-        Maximum log_lambda for this forest. Required only if Forest.wave_solution
-        is log.
+        Maximum log_lambda for this forest.
 
         log_lambda_min: float
-        Minimum log_lambda for this forest. Required only if Forest.wave_solution
-        is log.
+        Minimum log_lambda for this forest.
 
         Return
         ------
@@ -1101,34 +852,15 @@ def get_continuum_model(self, forest, aq, bq, **kwargs):
             raise ExpectedFluxError("Function get_cont_model requires "
                                     f"'mean_cont' in the **kwargs dictionary")
         mean_cont = kwargs.get("mean_cont")
-        if Forest.wave_solution == "log":
-            for key in ["log_lambda_max", "log_lambda_min"]:
-                if key not in kwargs:
-                    raise ExpectedFluxError("Function get_cont_model requires "
-                                            f"'{key}' in the **kwargs dictionary")
-            log_lambda_max = kwargs.get("log_lambda_max")
-            log_lambda_min = kwargs.get("log_lambda_min")
-        elif Forest.wave_solution == "lin":
-            for key in ["lambda_max", "lambda_min"]:
-                if key not in kwargs:
-                    raise ExpectedFluxError("Function get_cont_model requires "
-                                            f"'{key}' in the **kwargs dictionary")
-            lambda_max = kwargs.get("lambda_max")
-            lambda_min = kwargs.get("lambda_min")
-        else:
-            raise ExpectedFluxError("Forest.wave_solution must be either "
-                                    "'log' or 'lin'")
-
+        for key in ["log_lambda_max", "log_lambda_min"]:
+            if key not in kwargs:
+                raise ExpectedFluxError("Function get_cont_model requires "
+                                        f"'{key}' in the **kwargs dictionary")
+        log_lambda_max = kwargs.get("log_lambda_max")
+        log_lambda_min = kwargs.get("log_lambda_min")
         # compute continuum
-        if Forest.wave_solution == "log":
-            line = (bq * (forest.log_lambda - log_lambda_min) /
-                    (log_lambda_max - log_lambda_min) + aq)
-        elif Forest.wave_solution == "lin":
-            line = (bq * (forest.lambda_ - lambda_min) /
-                    (lambda_max - lambda_min) + aq)
-        else:
-            raise ExpectedFluxError("Forest.wave_solution must be either "
-                                    "'log' or 'lin'")
+        line = (bq * (forest.log_lambda - log_lambda_min) /
+                (log_lambda_max - log_lambda_min) + aq)
 
         return line * mean_cont
 
@@ -1153,25 +885,13 @@ def get_continuum_weights(self, forest, cont_model, **kwargs):
         if self.use_constant_weight:
             weights = np.ones_like(forest.flux)
         else:
-            if Forest.wave_solution == "log":
-                # pixel variance due to the Large Scale Strucure
-                var_lss = self.get_var_lss(forest.log_lambda)
-                # correction factor to the contribution of the pipeline
-                # estimate of the instrumental noise to the variance.
-                eta = self.get_eta(forest.log_lambda)
-                # fudge contribution to the variance
-                fudge = self.get_fudge(forest.log_lambda)
-            elif Forest.wave_solution == "lin":
-                # pixel variance due to the Large Scale Strucure
-                var_lss = self.get_var_lss(forest.lambda_)
-                # correction factor to the contribution of the pipeline
-                # estimate of the instrumental noise to the variance.
-                eta = self.get_eta(forest.lambda_)
-                # fudge contribution to the variance
-                fudge = self.get_fudge(forest.lambda_)
-            else:
-                raise ExpectedFluxError("Forest.wave_solution must be either "
-                                        "'log' or 'lin'")
+            # pixel variance due to the Large Scale Strucure
+            var_lss = self.get_var_lss(forest.log_lambda)
+            # correction factor to the contribution of the pipeline
+            # estimate of the instrumental noise to the variance.
+            eta = self.get_eta(forest.log_lambda)
+            # fudge contribution to the variance
+            fudge = self.get_fudge(forest.log_lambda)
 
             var_pipe = 1. / forest.ivar / cont_model**2
             ## prep_del.variance is the variance of delta
@@ -1194,20 +914,10 @@ def populate_los_ids(self, forests):
             if forest.bad_continuum_reason is not None:
                 continue
             # get the variance functions and statistics
-            if Forest.wave_solution == "log":
-                stack_delta = self.get_stack_delta(forest.log_lambda)
-                var_lss = self.get_var_lss(forest.log_lambda)
-                eta = self.get_eta(forest.log_lambda)
-                fudge = self.get_fudge(forest.log_lambda)
-            elif Forest.wave_solution == "lin":
-                stack_delta = self.get_stack_delta(forest.lambda_)
-                var_lss = self.get_var_lss(forest.lambda_)
-                eta = self.get_eta(forest.lambda_)
-                fudge = self.get_fudge(forest.lambda_)
-            else:
-                raise ExpectedFluxError("Forest.wave_solution must be either "
-                                        "'log' or 'lin'")
-
+            stack_delta = self.get_stack_delta(forest.log_lambda)
+            var_lss = self.get_var_lss(forest.log_lambda)
+            eta = self.get_eta(forest.log_lambda)
+            fudge = self.get_fudge(forest.log_lambda)
 
             mean_expected_flux = forest.continuum * stack_delta
             var_pipe = 1. / forest.ivar / mean_expected_flux**2
@@ -1239,10 +949,6 @@ def save_iteration_step(self, iteration):
         ---------
         iteration: int
         Iteration number. -1 for final iteration
-
-        Raise
-        -----
-        ExpectedFluxError if Forest.wave_solution is not 'lin' or 'log'
         """
         if iteration == -1:
             iter_out_file = self.iter_out_prefix + ".fits.gz"
@@ -1253,70 +959,36 @@ def save_iteration_step(self, iteration):
                          clobber=True) as results:
             header = {}
             header["FITORDER"] = self.order
-            if Forest.wave_solution == "log":
-                # TODO: update this once the TODO in compute continua is fixed
-                results.write([
-                    Forest.log_lambda_grid,
-                    self.get_stack_delta(Forest.log_lambda_grid),
-                    self.get_stack_delta_weights(Forest.log_lambda_grid)
-                ],
-                              names=['loglam', 'stack', 'weight'],
-                              header=header,
-                              extname='STACK_DELTAS')
-
-                results.write([
-                    self.log_lambda,
-                    self.get_eta(self.log_lambda),
-                    self.get_var_lss(self.log_lambda),
-                    self.get_fudge(self.log_lambda),
-                    self.get_num_pixels(self.log_lambda),
-                    self.get_valid_fit(self.log_lambda)
-                ],
-                              names=['loglam', 'eta', 'var_lss', 'fudge',
-                                     'num_pixels', 'valid_fit'],
-                              extname='VAR_FUNC')
-
-                results.write([
-                    Forest.log_lambda_rest_frame_grid,
-                    self.get_mean_cont(Forest.log_lambda_rest_frame_grid),
-                    self.get_mean_cont_weight(Forest.log_lambda_rest_frame_grid),
-                ],
-                              names=['loglam_rest', 'mean_cont', 'weight'],
-                              extname='CONT')
-            elif Forest.wave_solution == "lin":
-                # TODO: update this once the TODO in compute continua is fixed
-                results.write([
-                    Forest.lambda_grid,
-                    self.get_stack_delta(Forest.lambda_grid),
-                    self.get_stack_delta_weights(Forest.lambda_grid)
-                ],
-                              names=['lambda', 'stack', 'weight'],
-                              header=header,
-                              extname='STACK_DELTAS')
-
-                results.write([
-                    self.lambda_,
-                    self.get_eta(self.lambda_),
-                    self.get_var_lss(self.lambda_),
-                    self.get_fudge(self.lambda_),
-                    self.get_num_pixels(self.lambda_),
-                    self.get_valid_fit(self.lambda_)
-                ],
-                              names=['lambda', 'eta', 'var_lss', 'fudge',
-                                     'num_pixels', 'valid_fit'],
-                              extname='VAR_FUNC')
-
-                results.write([
-                    Forest.lambda_rest_frame_grid,
-                    self.get_mean_cont(Forest.lambda_rest_frame_grid),
-                    self.get_mean_cont_weight(Forest.lambda_rest_frame_grid),
-                ],
-                              names=['lambda_rest_frame', 'mean_cont', 'weight'],
-                              extname='CONT')
 
-            else:
-                raise ExpectedFluxError("Forest.wave_solution must be either "
-                                        "'log' or 'lin'")
+            # TODO: update this once the TODO in compute continua is fixed
+            results.write([
+                Forest.log_lambda_grid,
+                self.get_stack_delta(Forest.log_lambda_grid),
+                self.get_stack_delta_weights(Forest.log_lambda_grid)
+            ],
+                          names=['loglam', 'stack', 'weight'],
+                          header=header,
+                          extname='STACK_DELTAS')
+
+            results.write([
+                self.log_lambda_var_func_grid,
+                self.get_eta(self.log_lambda_var_func_grid),
+                self.get_var_lss(self.log_lambda_var_func_grid),
+                self.get_fudge(self.log_lambda_var_func_grid),
+                self.get_num_pixels(self.log_lambda_var_func_grid),
+                self.get_valid_fit(self.log_lambda_var_func_grid)
+            ],
+                          names=['loglam', 'eta', 'var_lss', 'fudge',
+                                 'num_pixels', 'valid_fit'],
+                          extname='VAR_FUNC')
+
+            results.write([
+                Forest.log_lambda_rest_frame_grid,
+                self.get_mean_cont(Forest.log_lambda_rest_frame_grid),
+                self.get_mean_cont_weight(Forest.log_lambda_rest_frame_grid),
+            ],
+                          names=['loglam_rest', 'mean_cont', 'weight'],
+                          extname='CONT')
 
 class LeastsSquaresContModel:
     def __init__(self, forest, expected_flux,
diff --git a/py/picca/delta_extraction/expected_fluxes/true_continuum.py b/py/picca/delta_extraction/expected_fluxes/true_continuum.py
index ef69540c7..599df6ed0 100644
--- a/py/picca/delta_extraction/expected_fluxes/true_continuum.py
+++ b/py/picca/delta_extraction/expected_fluxes/true_continuum.py
@@ -1,4 +1,5 @@
 """This module defines the class TrueContinuum"""
+from asyncio.log import logger
 import logging
 import multiprocessing
 
@@ -13,14 +14,16 @@
 from picca.delta_extraction.astronomical_objects.pk1d_forest import Pk1dForest
 from picca.delta_extraction.errors import ExpectedFluxError
 from picca.delta_extraction.expected_flux import ExpectedFlux
+from picca.delta_extraction.utils import find_bins
 
 accepted_options = ["input directory", "iter out prefix",
                     "num processors", "out dir",
-                    "raw statistics file"]
+                    "raw statistics file","use constant weight"]
 
 defaults = {
     "iter out prefix": "delta_attributes",
     "raw statistics file": "",
+    "use constant weight": False,
 }
 
 
@@ -94,7 +97,6 @@ def __init__(self, config):
         self.num_processors = None
         self.__parse_config(config)
 
-
         # read large scale structure variance and mean flux
         self.get_var_lss = None
         self.get_mean_flux = None
@@ -131,7 +133,7 @@ def __parse_config(self, config):
                 f"Found: {self.iter_out_prefix}")
 
         self.num_processors = config.getint("num processors")
-
+        self.use_constant_weight = config.get("use constant weight")
         self.raw_statistics_filename = config.get("raw statistics file")
 
 
@@ -157,92 +159,46 @@ def compute_expected_flux(self, forests):
             forests = pool.map(self.read_true_continuum, forests)
             pool.close()
 
-        if Forest.wave_solution == "log":
-            self.compute_mean_cont_log(forests)
-        elif Forest.wave_solution == "lin":
-            self.compute_mean_cont_lin(forests)
-        else:
-            raise ExpectedFluxError("Forest.wave_solution must be "
-                                            "either 'log' or 'linear'")
+        self.compute_mean_cont(forests)
+
         # now loop over forests to populate los_ids
         self.populate_los_ids(forests)
         # Save delta atributes
         self.save_delta_attributes()
 
-    def compute_mean_cont_lin(self, forests):
-        """Compute the mean quasar continuum over the whole sample assuming a
-        linear wavelength solution. Then updates the value of self.get_mean_cont
-        to contain it
+    def compute_mean_cont(self, forests):
+        """Compute the mean quasar continuum over the whole sample.
+        Then updates the value of self.get_mean_cont to contain it
 
         Arguments
         ---------
         forests: List of Forest
         A list of Forest from which to compute the deltas.
         """
-        mean_cont = np.zeros_like(Forest.lambda_rest_frame_grid.size)
-        mean_cont_weight = np.zeros_like(Forest.lambda_rest_frame_grid.size)
-
-        for forest in forests:
-            if forest.bad_continuum_reason is not None:
-                continue
-            bins = find_bins(
-                forest.lambda_ / (1 + forest.z),
-                Forest.lambda_rest_frame_grid
-                )
-
-            var_lss = self.get_var_lss(forest.lambda_)
-            var_pipe = 1. / forest.ivar / forest.continuum**2
-            variance = var_lss + var_pipe
-            weights = 1 / variance
-            cont = np.bincount(bins,
-                               weights=forest.continuum * weights)
-            mean_cont[:len(cont)] += cont
-            cont = np.bincount(bins, weights=weights)
-            mean_cont_weight[:len(cont)] += cont
-
-        w = mean_cont_weight > 0
-        mean_cont[w] /= mean_cont_weight[w]
-        mean_cont /= mean_cont.mean()
-        lambda_cont = Forest.lambda_rest_frame_grid[w]
-
-        self.get_mean_cont = interp1d(lambda_cont,
-                                      mean_cont,
-                                      fill_value="extrapolate")
-        self.get_mean_cont_weight = interp1d(lambda_cont,
-                                             mean_cont_weight,
-                                             fill_value=0.0,
-                                             bounds_error=False)
-
-    def compute_mean_cont_log(self, forests):
-        """Compute the mean quasar continuum over the whole sample assuming a
-        log-linear wavelength solution. Then updates the value of
-        self.get_mean_cont to contain it
-
-        Arguments
-        ---------
-        forests: List of Forest
-        A list of Forest from which to compute the deltas.
-        """
-        mean_cont = np.zeros_like(Forest.log_lambda_rest_frame_grid.size)
-        mean_cont_weight = np.zeros_like(Forest.log_lambda_rest_frame_grid.size)
+        mean_cont = np.zeros_like(Forest.log_lambda_rest_frame_grid)
+        mean_cont_weight = np.zeros_like(Forest.log_lambda_rest_frame_grid)
 
         for forest in forests:
             if forest.bad_continuum_reason is not None:
                 continue
             bins = find_bins(
                 forest.log_lambda -  np.log10(1 + forest.z),
-                Forest.log_lambda_rest_frame_grid
+                Forest.log_lambda_rest_frame_grid,
+                Forest.wave_solution
             )
 
-            var_lss = self.get_var_lss(forest.log_lambda)
-            var_pipe = 1. / forest.ivar / forest.continuum**2
-            variance = var_lss + var_pipe
-            weights = 1 / variance
+            if self.use_constant_weight:
+                weights = np.ones_like(forest.log_lambda)
+            else:
+                var_lss = self.get_var_lss(forest.log_lambda)
+                var_pipe = 1. / forest.ivar / forest.continuum**2
+                variance = var_lss + var_pipe
+                weights = 1 / variance
             cont = np.bincount(bins,
                                weights= forest.continuum * weights)
             mean_cont[:len(cont)] += cont
-            cont = np.bincount(bins, weights=weights)
-            mean_cont_weight[:len(cont)] += cont
+            cont_weight = np.bincount(bins, weights=weights)
+            mean_cont_weight[:len(cont)] += cont_weight
 
         w = mean_cont_weight > 0
         mean_cont[w] /= mean_cont_weight[w]
@@ -290,16 +246,11 @@ def read_true_continuum(self, forest):
         indx = np.where(true_cont["TARGETID"]==forest.targetid)
         true_continuum = interp1d(lambda_, true_cont["TRUE_CONT"][indx])
 
-        if Forest.wave_solution == "log":
-            forest.continuum = true_continuum(10**forest.log_lambda)[0]
-            forest.continuum *= self.get_mean_flux(10**forest.log_lambda)
-
-        elif Forest.wave_solution == "lin":
-            forest.continuum = true_continuum(forest.lambda_)[0]
-            forest.continuum *= self.get_mean_flux(forest.lambda_)
-        else:
-            raise ExpectedFluxError("Forest.wave_solution must be either 'log' "
-                                    "or 'lin'")
+        forest.continuum = true_continuum(10**forest.log_lambda)[0]
+        mean_optical_depth = np.ones(forest.log_lambda.size)
+        tau, gamma, lambda_rest_frame = 0.0023, 3.64, 1215.67
+        w = 10.**forest.log_lambda / (1. + forest.z) <= lambda_rest_frame
+        z = 10.**forest.log_lambda / lambda_rest_frame - 1.
 
         return forest
 
@@ -308,15 +259,19 @@ def read_raw_statistics(self):
         """
         #files are only for lya so far, this will need to be updated so that regions other than Lya are available
 
+        if self.use_constant_weight:
+            logger.info("using constant weights, skipping raw statistics as not needed")
+            return
+        
         if self.raw_statistics_filename != "":
             filename = self.raw_statistics_filename
         else:
             filename = resource_filename('picca', 'delta_extraction') + '/expected_fluxes/raw_stats/'
             if Forest.wave_solution == "log":
                 filename += 'colore_v9_lya_log.fits.gz'
-            elif Forest.wave_solution == "lin" and Forest.delta_lambda == 2.4:
+            elif Forest.wave_solution == "lin" and (10**Forest.log_lambda_grid[1]-10**Forest.log_lambda_grid[0]) == 2.4:
                 filename += 'colore_v9_lya_lin_2.4.fits.gz'
-            elif Forest.wave_solution == "lin" and Forest.delta_lambda == 3.2:
+            elif Forest.wave_solution == "lin" and (10**Forest.log_lambda_grid[1]-10**Forest.log_lambda_grid[0]) == 3.2:
                 filename += 'colore_v9_lya_lin_3.2.fits.gz'
             else:
                 raise ExpectedFluxError("Couldn't find compatible raw satistics file. Provide a custom one using 'raw statistics file' field.")
@@ -331,30 +286,30 @@ def read_raw_statistics(self):
         if Forest.wave_solution == "log":
             if (
                 header['LINEAR']
-                or not np.isclose(header['L_MIN'], 10**Forest.log_lambda_min, rtol=1e-3)
-                or not np.isclose(header['L_MAX'], 10**Forest.log_lambda_max, rtol=1e-3)
-                or not np.isclose(header['LR_MIN'], 10**Forest.log_lambda_min_rest_frame, rtol=1e-3)
-                or not np.isclose(header['LR_MAX'], 10**Forest.log_lambda_max_rest_frame, rtol=1e-3)
-                or not np.isclose(header['DEL_LL'], Forest.delta_log_lambda, rtol=1e-3)
+                or not np.isclose(header['L_MIN'], 10**Forest.log_lambda_grid[0], rtol=1e-3)
+                or not np.isclose(header['L_MAX'], 10**Forest.log_lambda_grid[-1], rtol=1e-3)
+                or not np.isclose(header['LR_MIN'], 10**Forest.log_lambda_rest_frame_grid[0], rtol=1e-3)
+                or not np.isclose(header['LR_MAX'], 10**Forest.log_lambda_rest_frame_grid[-1], rtol=1e-3)
+                or not np.isclose(header['DEL_LL'], Forest.log_lambda_grid[1]-Forest.log_lambda_grid[0], rtol=1e-3)
             ):
                 raise ExpectedFluxError(f'''raw statistics file pixelization scheme does not match input pixelization scheme.
                 \t\tL_MIN\tL_MAX\tLR_MIN\tLR_MAX\tDEL_LL
                 raw\t{header['L_MIN']}\t{header['L_MAX']}\t{header['LR_MIN']}\t{header['LR_MAX']}\t{header['DEL_LL']}
-                input\t{10**Forest.log_lambda_min}\t{10**Forest.log_lambda_max}\t{10**Forest.log_lambda_min_rest_frame}\t{10**Forest.log_lambda_max_rest_frame}\t{Forest.delta_log_lambda}
+                input\t{10**Forest.log_lambda_grid[0]}\t{10**Forest.log_lambda_grid[-1]}\t{10**Forest.log_lambda_rest_frame_grid[0]}\t{10**Forest.log_lambda_rest_frame_grid[-1]}\t{Forest.delta_log_lambda}
                 provide a custom file in 'raw statistics file' field matching input pixelization scheme''')
         elif Forest.wave_solution == "lin":
             if (
                 not header['LINEAR']
-                or not np.isclose(header['L_MIN'], Forest.lambda_min , rtol=1e-3)
-                or not np.isclose(header['L_MAX'], Forest.lambda_max , rtol=1e-3)
-                or not np.isclose(header['LR_MIN'], Forest.lambda_min_rest_frame, rtol=1e-3)
-                or not np.isclose(header['LR_MAX'], Forest.lambda_max_rest_frame, rtol=1e-3)
-                or not np.isclose(header['DEL_L'], Forest.delta_lambda, rtol=1e-3)
+                or not np.isclose(header['L_MIN'], 10**Forest.log_lambda_grid[0] , rtol=1e-3)
+                or not np.isclose(header['L_MAX'], 10**Forest.log_lambda_grid[-1] , rtol=1e-3)
+                or not np.isclose(header['LR_MIN'], 10**Forest.log_lambda_rest_frame_grid[0], rtol=1e-3)
+                or not np.isclose(header['LR_MAX'], 10**Forest.log_lambda_rest_frame_grid[-1], rtol=1e-3)
+                or not np.isclose(header['DEL_LL'], 10**Forest.log_lambda_grid[1]-10**Forest.log_lambda_grid[0], rtol=1e-3)
             ):
                 raise ExpectedFluxError(f'''raw statistics file pixelization scheme does not match input pixelization scheme.
-                \t\tL_MIN\tL_MAX\tLR_MIN\tLR_MAX\tDEL_LL
+                \t\tL_MIN\tL_MAX\tLR_MIN\tLR_MAX\tDEL_L
                 raw\t{header['L_MIN']}\t{header['L_MAX']}\t{header['LR_MIN']}\t{header['LR_MAX']}\t{header['DEL_LL']}
-                input\t{10**Forest.log_lambda_min}\t{10**Forest.log_lambda_max}\t{10**Forest.log_lambda_min_rest_frame}\t{10**Forest.log_lambda_max_rest_frame}\t{Forest.delta_log_lambda}
+                input\t{10**Forest.log_lambda_grid[0]}\t{10**Forest.log_lambda_grid[-1]}\t{ 10**Forest.log_lambda_rest_frame_grid[0]}\t{ 10**Forest.log_lambda_rest_frame_grid[-1]}\t{10**Forest.log_lambda_grid[1]-10**Forest.log_lambda_grid[0]}
                 provide a custom file in 'raw statistics file' field matching input pixelization scheme''')
 
         lambda_ = hdul[1].data['LAMBDA']
@@ -387,18 +342,16 @@ def populate_los_ids(self, forests):
             if forest.bad_continuum_reason is not None:
                 continue
             # get the variance functions
-            if Forest.wave_solution == "log":
-                var_lss = self.get_var_lss(forest.log_lambda)
-            elif Forest.wave_solution == "lin":
-                var_lss = self.get_var_lss(forest.lambda_)
+            if self.use_constant_weight:
+                weights = np.ones_like(forest.log_lambda)
+                mean_expected_flux = forest.continuum
             else:
-                raise ExpectedFluxError("Forest.wave_solution must be either "
-                                        "'log' or 'lin'")
+                var_lss = self.get_var_lss(forest.log_lambda)
 
-            mean_expected_flux = forest.continuum
-            var_pipe = 1. / forest.ivar/ forest.continuum**2
-            variance =  var_lss + var_pipe
-            weights = 1. / variance
+                mean_expected_flux = forest.continuum
+                var_pipe = 1. / forest.ivar/ forest.continuum**2
+                variance =  var_lss + var_pipe
+                weights = 1. / variance
 
             if isinstance(forest, Pk1dForest):
                 ivar = forest.ivar / mean_expected_flux**2
@@ -429,29 +382,11 @@ def save_delta_attributes(self):
                          clobber=True) as results:
             header = {}
             header["FITORDER"] = -1
-            if Forest.wave_solution == "log":
-                num_bins = int((Forest.log_lambda_max - Forest.log_lambda_min) /
-                               Forest.delta_log_lambda) + 1
-
-                results.write([
-                    Forest.log_lambda_rest_frame_grid,
-                    self.get_mean_cont(Forest.log_lambda_rest_frame_grid),
-                    self.get_mean_cont_weight(Forest.log_lambda_rest_frame_grid),
-                ],
-                              names=['loglam_rest', 'mean_cont', 'weight'],
-                              extname='CONT')
-            elif Forest.wave_solution == "lin":
-                num_bins = int((Forest.lambda_max - Forest.lambda_min) /
-                               Forest.delta_lambda) + 1
-
-                results.write([
-                    Forest.lambda_rest_frame_grid,
-                    self.get_mean_cont(Forest.lambda_rest_frame_grid),
-                    self.get_mean_cont_weight(Forest.lambda_rest_frame_grid),
-                ],
-                              names=['lambda_rest_frame', 'mean_cont', 'weight'],
-                              extname='CONT')
-
-            else:
-                raise ExpectedFluxError("Forest.wave_solution must be either "
-                                        "'log' or 'lin'")
+            results.write([
+                Forest.log_lambda_rest_frame_grid,
+                self.get_mean_cont(Forest.log_lambda_rest_frame_grid),
+                self.get_mean_cont_weight(Forest.log_lambda_rest_frame_grid),
+            ],
+                          names=['loglam_rest', 'mean_cont', 'weight'],
+                          extname='CONT',
+                          header=header)
diff --git a/py/picca/delta_extraction/masks/absorber_mask.py b/py/picca/delta_extraction/masks/absorber_mask.py
index 21fe7535a..88dd235f8 100644
--- a/py/picca/delta_extraction/masks/absorber_mask.py
+++ b/py/picca/delta_extraction/masks/absorber_mask.py
@@ -111,21 +111,13 @@ def apply_mask(self, forest):
         -----
         CorrectionError if Forest.wave_solution is not 'log'
         """
-        if Forest.wave_solution == "log":
-            log_lambda = forest.log_lambda
-        elif Forest.wave_solution == "lin":
-            log_labmda = np.log10(forest.lambda_)
-        else:
-            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)
+            w = np.ones(forest.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 * (forest.log_lambda - np.log10(lambda_absorber))) >
+                      self.absorber_mask_width)
 
             # do the actual masking
             for param in Forest.mask_fields:
diff --git a/py/picca/delta_extraction/masks/bal_mask.py b/py/picca/delta_extraction/masks/bal_mask.py
index 901d1900d..d20e52646 100644
--- a/py/picca/delta_extraction/masks/bal_mask.py
+++ b/py/picca/delta_extraction/masks/bal_mask.py
@@ -208,24 +208,11 @@ def apply_mask(self, forest):
         if (mask_table is not None) and len(mask_table) > 0:
 
             # find out which pixels to mask
-            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']))
-            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)
-                    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'")
+            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']))
 
             # do the actual masking
             for param in Forest.mask_fields:
diff --git a/py/picca/delta_extraction/masks/dla_mask.py b/py/picca/delta_extraction/masks/dla_mask.py
index a7024dc9c..6ddc3e7b3 100644
--- a/py/picca/delta_extraction/masks/dla_mask.py
+++ b/py/picca/delta_extraction/masks/dla_mask.py
@@ -138,13 +138,7 @@ def apply_mask(self, forest):
         -----
         MaskError if Forest.wave_solution is not 'log'
         """
-        if Forest.wave_solution == "log":
-            lambda_ = 10**forest.log_lambda
-        elif Forest.wave_solution == "lin":
-            lambda_ = forest.lambda_
-        else:
-            raise MaskError(
-                "Forest.wave_solution must be either 'log' or 'lin'")
+        lambda_ = 10**forest.log_lambda
 
         # load DLAs
         if self.los_ids.get(forest.los_id) is not None:
diff --git a/py/picca/delta_extraction/masks/lines_mask.py b/py/picca/delta_extraction/masks/lines_mask.py
index db8d8f09b..ffb84d1df 100644
--- a/py/picca/delta_extraction/masks/lines_mask.py
+++ b/py/picca/delta_extraction/masks/lines_mask.py
@@ -77,29 +77,14 @@ def apply_mask(self, forest):
         CorrectionError if Forest.wave_solution is not 'lin' or 'log'
         """
         # find masking array
-        if Forest.wave_solution == "log":
-            w = np.ones(forest.log_lambda.size, dtype=bool)
-            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)
-            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']))
-        elif Forest.wave_solution == "lin":
-            w = np.ones(forest.lambda_.size, dtype=bool)
-            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)
-            for mask_range in self.mask_rest_frame:
-                w &= ((lambda_rest_frame < mask_range['wave_min']) |
-                      (lambda_rest_frame > mask_range['wave_max']))
-        else:
-            raise MaskError("Unable to apply SkyMask. Forest.wave_solution is "
-                            "required to be either 'log' or 'lin'. Found "
-                            f"{Forest.wave_solution}")
+        w = np.ones(forest.log_lambda.size, dtype=bool)
+        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)
+        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']))
 
         # do the actual masking
         for param in Forest.mask_fields:
diff --git a/py/picca/delta_extraction/utils.py b/py/picca/delta_extraction/utils.py
index a72e824fa..1054ccfef 100644
--- a/py/picca/delta_extraction/utils.py
+++ b/py/picca/delta_extraction/utils.py
@@ -127,7 +127,7 @@ def class_from_string(class_name, module_name):
         accepted_options = []
     return class_object, default_args, accepted_options
 
-def find_bins(original_array, grid_array):
+def find_bins(original_array, grid_array, wave_solution):
     """For each element in original_array, find the corresponding bin in grid_array
 
     Arguments
@@ -138,13 +138,24 @@ def find_bins(original_array, grid_array):
     grid_array: array of float
     Common array, e.g. Forest.log_lambda_grid
 
+    wave_solution: "log" or "lin"
+    Specifies whether we want to construct a wavelength grid that is evenly
+    spaced on wavelength (lin) or on the logarithm of the wavelength (log)
+
     Return
     ------
     found_bin: array of int
     An array of size original_array.size filled with values smaller than
     grid_array.size with the bins correspondance
     """
-    found_bin = (np.abs(grid_array - original_array[:,np.newaxis])).argmin(axis=1)
+    # TODO: check the effect of finding the nearest bin in log_lambda space
+    # versus lambda space. Once we understand this we can possibly remove
+    # the dependence from Forest from here
+    if wave_solution == "log":
+        found_bin = (np.abs(grid_array - original_array[:,np.newaxis])).argmin(axis=1)
+    elif wave_solution == "lin":
+        found_bin = (np.abs(10**grid_array - 10**original_array[:,np.newaxis])).argmin(axis=1)
+
     return found_bin
 
 PROGRESS_LEVEL_NUM = 15


Test patch:
diff --git a/py/picca/tests/delta_extraction/abstract_test.py b/py/picca/tests/delta_extraction/abstract_test.py
index 8c7047db3..6ce2f4ec2 100644
--- a/py/picca/tests/delta_extraction/abstract_test.py
+++ b/py/picca/tests/delta_extraction/abstract_test.py
@@ -55,8 +55,8 @@ def compare_ascii(self, orig_file, new_file):
 
                 if not orig_line == new_line:
                     print("Lines not equal")
-                    print("Original line" + orig_line)
-                    print("New line" + new_line)
+                    print("Original line " + orig_line)
+                    print("New line " + new_line)
                     self.assertTrue(orig_line == new_line)
         finally:
             orig.close()
@@ -91,12 +91,15 @@ def compare_fits(self, orig_file, new_file):
                         if orig_header[key] != new_header[key]:
                             print(f"\nOriginal file: {orig_file}")
                             print(f"New file: {new_file}")
+                            print(f"\n For header {orig_header['EXTNAME']}")
                             print(f"Different values found for key {key}: "
                                   f"orig: {orig_header[key]}, new: {new_header[key]}")
                         self.assertTrue((orig_header[key] == new_header[key]) or
                                         (np.isclose(orig_header[key], new_header[key])))
                 for key in new_header:
                     if key not in orig_header:
+                        print(f"\nOriginal file: {orig_file}")
+                        print(f"New file: {new_file}")
                         print(f"key {key} missing in orig header")
                         if key in ["MEANSNR", "BLINDING"]:
                             continue
diff --git a/py/picca/tests/delta_extraction/astronomical_object_tests.py b/py/picca/tests/delta_extraction/astronomical_object_tests.py
index bd960e909..3ea7dfe67 100644
--- a/py/picca/tests/delta_extraction/astronomical_object_tests.py
+++ b/py/picca/tests/delta_extraction/astronomical_object_tests.py
@@ -46,9 +46,10 @@
 ])
 LOG_LAMBDA_REBIN = np.array(
     [3.5563025, 3.5566025, 3.5569025, 3.5572025, 3.5575025])
-LAMBDA_ = np.array(
-    [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])
+LOG_LAMBDA_LIN = np.log10(np.array([
+3610.0, 3610.4, 3650.0, 3650.4, 3670.0, 3670.4, 3680.0, 3680.4, 3700.0, 3700.4
+]))
+LOG_LAMBDA_REBIN_LIN = np.log10(np.array([3610, 3650, 3670, 3680, 3700]))
 
 RESOMAT = np.ones([7, 10])
 RESOMAT2 = np.ones([7, 10])*3
@@ -120,12 +121,12 @@
 
 kwargs_forest_lin = kwargs_forest.copy()
 kwargs_forest_lin.update({
-    "lambda": LAMBDA_,
+    "log_lambda": LOG_LAMBDA_LIN,
 })
 
 kwargs_forest_lin2 = kwargs_forest2.copy()
 kwargs_forest_lin2.update({
-    "lambda": LAMBDA_,
+    "log_lambda": LOG_LAMBDA_LIN,
 })
 
 kwargs_forest_rebin = kwargs_astronomical_object.copy()
@@ -141,7 +142,7 @@
 
 kwargs_forest_lin_rebin = kwargs_forest_rebin.copy()
 kwargs_forest_lin_rebin.update({
-    "lambda": LAMBDA_REBIN,
+    "log_lambda": LOG_LAMBDA_REBIN_LIN,
 })
 
 kwargs_forest_coadd = kwargs_astronomical_object.copy()
@@ -157,7 +158,7 @@
 
 kwargs_forest_lin_coadd = kwargs_forest_coadd.copy()
 kwargs_forest_lin_coadd.update({
-    "lambda": LAMBDA_REBIN,
+    "log_lambda": LOG_LAMBDA_REBIN_LIN,
 })
 
 # define contructors for Pk1dForest
@@ -188,12 +189,12 @@
 
 kwargs_pk1d_forest_lin = kwargs_pk1d_forest.copy()
 kwargs_pk1d_forest_lin.update({
-    "lambda": LAMBDA_,
+    "log_lambda": LOG_LAMBDA_LIN,
 })
 
 kwargs_pk1d_forest_lin2 = kwargs_pk1d_forest2.copy()
 kwargs_pk1d_forest_lin2.update({
-    "lambda": LAMBDA_,
+    "log_lambda": LOG_LAMBDA_LIN,
 })
 
 kwargs_pk1d_forest_rebin = kwargs_forest_rebin.copy()
@@ -210,7 +211,7 @@
 
 kwargs_pk1d_forest_lin_rebin = kwargs_pk1d_forest_rebin.copy()
 kwargs_pk1d_forest_lin_rebin.update({
-    "lambda": LAMBDA_REBIN,
+    "log_lambda": LOG_LAMBDA_REBIN_LIN,
 })
 
 kwargs_pk1d_forest_coadd = kwargs_forest_coadd.copy()
@@ -227,7 +228,7 @@
 
 kwargs_pk1d_forest_lin_coadd = kwargs_pk1d_forest_coadd.copy()
 kwargs_pk1d_forest_lin_coadd.update({
-    "lambda": LAMBDA_REBIN,
+    "log_lambda": LOG_LAMBDA_REBIN_LIN,
 })
 
 # define contructors for DesiForest
@@ -333,11 +334,11 @@
 
 # define contructors for SdssPk1dForest
 kwargs_sdss_pk1d_forest = kwargs_sdss_forest.copy()
-kwargs_sdss_pk1d_forest.update(kwargs_pk1d_forest_lin)
+kwargs_sdss_pk1d_forest.update(kwargs_pk1d_forest)
 del kwargs_sdss_pk1d_forest["los_id"]
 
 kwargs_sdss_pk1d_forest2 = kwargs_sdss_forest2.copy()
-kwargs_sdss_pk1d_forest2.update(kwargs_pk1d_forest_lin2)
+kwargs_sdss_pk1d_forest2.update(kwargs_pk1d_forest2)
 del kwargs_sdss_pk1d_forest2["los_id"]
 
 kwargs_sdss_pk1d_forest_rebin = kwargs_pk1d_forest_rebin.copy()
@@ -433,13 +434,12 @@ def assert_forest_object(self, test_obj, kwargs):
             self.assertTrue(np.allclose(test_obj.deltas, kwargs.get("deltas")))
         else:
             self.assertTrue(test_obj.deltas is None)
-        if Forest.wave_solution == "log":
-            self.assertTrue(
-                np.allclose(test_obj.log_lambda, kwargs.get("log_lambda")))
-            self.assertTrue(test_obj.lambda_ is None)
-        elif Forest.wave_solution == "lin":
-            self.assertTrue(test_obj.log_lambda is None)
-            self.assertTrue(np.allclose(test_obj.lambda_, kwargs.get("lambda")))
+        if len(test_obj.log_lambda) != len(kwargs.get("log_lambda")):
+            print()
+            print(test_obj.log_lambda)
+            print(kwargs.get("log_lambda"))
+        self.assertTrue(
+            np.allclose(test_obj.log_lambda, kwargs.get("log_lambda")))
         flux = kwargs.get("flux")
         ivar = kwargs.get("ivar")
         self.assertTrue(np.allclose(test_obj.flux, flux))
@@ -454,10 +454,7 @@ def assert_forest_object(self, test_obj, kwargs):
         self.assertTrue(Forest.mask_fields[0] == "flux")
         self.assertTrue(Forest.mask_fields[1] == "ivar")
         self.assertTrue(Forest.mask_fields[2] == "transmission_correction")
-        if Forest.wave_solution == "log":
-            self.assertTrue(Forest.mask_fields[3] == "log_lambda")
-        elif Forest.wave_solution == "lin":
-            self.assertTrue(Forest.mask_fields[3] == "lambda_")
+        self.assertTrue(Forest.mask_fields[3] == "log_lambda")
         self.assertTrue(
             np.allclose(test_obj.transmission_correction, np.ones_like(flux)))
         mean_snr = (flux * np.sqrt(ivar)).mean()
@@ -489,15 +486,10 @@ def assert_forest_object(self, test_obj, kwargs):
             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]) +
-                           np.power(10., log_lambda[0])) / 2. /
-                          Pk1dForest.lambda_abs_igm - 1.0)
-            elif Forest.wave_solution == "lin":
-                lambda_ = kwargs.get("lambda")
-                mean_z = ((lambda_[len(lambda_) - 1] + lambda_[0]) / 2. /
-                          Pk1dForest.lambda_abs_igm - 1.0)
+            log_lambda = kwargs.get("log_lambda")
+            mean_z = ((np.power(10., log_lambda[len(log_lambda) - 1]) +
+                       np.power(10., log_lambda[0])) / 2. /
+                      Pk1dForest.lambda_abs_igm - 1.0)
             if not np.isclose(test_obj.mean_z, mean_z):
                 print(test_obj.mean_z, mean_z)
                 print(log_lambda, test_obj.log_lambda,
@@ -525,7 +517,7 @@ def assert_get_data(self, test_obj):
             self.assertTrue(comments[0] == "Log lambda")
         elif Forest.wave_solution == "lin":
             self.assertTrue(names[0] == "LAMBDA")
-            self.assertTrue(np.allclose(cols[0], test_obj.lambda_))
+            self.assertTrue(np.allclose(cols[0], 10**test_obj.log_lambda))
             self.assertTrue(units[0] == "Angstrom")
             self.assertTrue(comments[0] == "Lambda")
 
@@ -597,7 +589,20 @@ def assert_get_header(self, test_obj):
             self.assertTrue(header[index + 1].get("value") == test_obj.mean_snr)
             self.assertTrue(header[index + 2].get("name") == "BLINDING")
             self.assertTrue(header[index + 2].get("value") == "none")
-            index += 2
+            self.assertTrue(header[index + 3].get("name") == "WAVE_SOLUTION")
+            if Forest.wave_solution == "log":
+                self.assertTrue(header[index + 3].get("value") == "log")
+                self.assertTrue(header[index + 4].get("name") == "DELTA_LOG_LAMBDA")
+                self.assertTrue(np.isclose(header[index + 4].get("value"), 3e-4))
+            elif Forest.wave_solution == "lin":
+                self.assertTrue(header[index + 3].get("value") == "lin")
+                self.assertTrue(header[index + 4].get("name") == "DELTA_LAMBDA")
+                self.assertTrue(np.isclose(header[index + 4].get("value"), 1.0))
+            else:
+                print(f"Forest.wave_solution={Forest.wave_solution}, expected "
+                      "'log' or 'lin'")
+                self.assertTrue(False)
+            index += 4
         if isinstance(test_obj, Pk1dForest):
             self.assertTrue(header[index + 1].get("name") == "MEANZ")
             self.assertTrue(header[index + 1].get("value") == test_obj.mean_z)
diff --git a/py/picca/tests/delta_extraction/data/delta_extraction_calib/delta-28.fits.gz b/py/picca/tests/delta_extraction/data/delta_extraction_calib/delta-28.fits.gz
index 56fe110b9..6b67c8d68 100644
Binary files a/py/picca/tests/delta_extraction/data/delta_extraction_calib/delta-28.fits.gz and b/py/picca/tests/delta_extraction/data/delta_extraction_calib/delta-28.fits.gz differ
diff --git a/py/picca/tests/delta_extraction/data/delta_extraction_calib/delta-45.fits.gz b/py/picca/tests/delta_extraction/data/delta_extraction_calib/delta-45.fits.gz
index 859c886b3..088c6f361 100644
Binary files a/py/picca/tests/delta_extraction/data/delta_extraction_calib/delta-45.fits.gz and b/py/picca/tests/delta_extraction/data/delta_extraction_calib/delta-45.fits.gz differ
diff --git a/py/picca/tests/delta_extraction/data/delta_extraction_calib2/delta-28.fits.gz b/py/picca/tests/delta_extraction/data/delta_extraction_calib2/delta-28.fits.gz
index 369fe8807..c5fa91db9 100644
Binary files a/py/picca/tests/delta_extraction/data/delta_extraction_calib2/delta-28.fits.gz and b/py/picca/tests/delta_extraction/data/delta_extraction_calib2/delta-28.fits.gz differ
diff --git a/py/picca/tests/delta_extraction/data/delta_extraction_calib2/delta-45.fits.gz b/py/picca/tests/delta_extraction/data/delta_extraction_calib2/delta-45.fits.gz
index 582164475..2b330d237 100644
Binary files a/py/picca/tests/delta_extraction/data/delta_extraction_calib2/delta-45.fits.gz and b/py/picca/tests/delta_extraction/data/delta_extraction_calib2/delta-45.fits.gz differ
diff --git a/py/picca/tests/delta_extraction/data/delta_extraction_calib2_nomask/delta-28.fits.gz b/py/picca/tests/delta_extraction/data/delta_extraction_calib2_nomask/delta-28.fits.gz
index 8cf7e15cd..a2a92fa81 100644
Binary files a/py/picca/tests/delta_extraction/data/delta_extraction_calib2_nomask/delta-28.fits.gz and b/py/picca/tests/delta_extraction/data/delta_extraction_calib2_nomask/delta-28.fits.gz differ
diff --git a/py/picca/tests/delta_extraction/data/delta_extraction_calib2_nomask/delta-45.fits.gz b/py/picca/tests/delta_extraction/data/delta_extraction_calib2_nomask/delta-45.fits.gz
index dc2274a8f..d1328507d 100644
Binary files a/py/picca/tests/delta_extraction/data/delta_extraction_calib2_nomask/delta-45.fits.gz and b/py/picca/tests/delta_extraction/data/delta_extraction_calib2_nomask/delta-45.fits.gz differ
diff --git a/py/picca/tests/delta_extraction/data/delta_extraction_calib_nomask/delta-28.fits.gz b/py/picca/tests/delta_extraction/data/delta_extraction_calib_nomask/delta-28.fits.gz
index 0bcb63096..f999f52d6 100644
Binary files a/py/picca/tests/delta_extraction/data/delta_extraction_calib_nomask/delta-28.fits.gz and b/py/picca/tests/delta_extraction/data/delta_extraction_calib_nomask/delta-28.fits.gz differ
diff --git a/py/picca/tests/delta_extraction/data/delta_extraction_calib_nomask/delta-45.fits.gz b/py/picca/tests/delta_extraction/data/delta_extraction_calib_nomask/delta-45.fits.gz
index 34cd6cdbb..306a12c69 100644
Binary files a/py/picca/tests/delta_extraction/data/delta_extraction_calib_nomask/delta-45.fits.gz and b/py/picca/tests/delta_extraction/data/delta_extraction_calib_nomask/delta-45.fits.gz differ
diff --git a/py/picca/tests/delta_extraction/data/delta_extraction_lya_nodla/delta-28.fits.gz b/py/picca/tests/delta_extraction/data/delta_extraction_lya_nodla/delta-28.fits.gz
index c33bcb0b0..c631217f9 100644
Binary files a/py/picca/tests/delta_extraction/data/delta_extraction_lya_nodla/delta-28.fits.gz and b/py/picca/tests/delta_extraction/data/delta_extraction_lya_nodla/delta-28.fits.gz differ
diff --git a/py/picca/tests/delta_extraction/data/delta_extraction_lya_nodla/delta-45.fits.gz b/py/picca/tests/delta_extraction/data/delta_extraction_lya_nodla/delta-45.fits.gz
index e0a695550..2d4c252aa 100644
Binary files a/py/picca/tests/delta_extraction/data/delta_extraction_lya_nodla/delta-45.fits.gz and b/py/picca/tests/delta_extraction/data/delta_extraction_lya_nodla/delta-45.fits.gz differ
diff --git a/py/picca/tests/delta_extraction/data/delta_extraction_lya_nomask_nodla/delta-28.fits.gz b/py/picca/tests/delta_extraction/data/delta_extraction_lya_nomask_nodla/delta-28.fits.gz
index c23b4964a..8a3ab236c 100644
Binary files a/py/picca/tests/delta_extraction/data/delta_extraction_lya_nomask_nodla/delta-28.fits.gz and b/py/picca/tests/delta_extraction/data/delta_extraction_lya_nomask_nodla/delta-28.fits.gz differ
diff --git a/py/picca/tests/delta_extraction/data/delta_extraction_lya_nomask_nodla/delta-45.fits.gz b/py/picca/tests/delta_extraction/data/delta_extraction_lya_nomask_nodla/delta-45.fits.gz
index 617bdfee1..38c92db11 100644
Binary files a/py/picca/tests/delta_extraction/data/delta_extraction_lya_nomask_nodla/delta-45.fits.gz and b/py/picca/tests/delta_extraction/data/delta_extraction_lya_nomask_nodla/delta-45.fits.gz differ
diff --git a/py/picca/tests/delta_extraction/expected_flux_tests.py b/py/picca/tests/delta_extraction/expected_flux_tests.py
index 96f9e6a6a..8df41b59f 100644
--- a/py/picca/tests/delta_extraction/expected_flux_tests.py
+++ b/py/picca/tests/delta_extraction/expected_flux_tests.py
@@ -80,15 +80,14 @@ def test_dr16_expected_flux(self):
             expected_flux = Dr16ExpectedFlux(config["expected flux"])
 
         # setup Forest variables; case: logarithmic wavelength solution
-        setup_forest("log")
+        setup_forest("log", rebin=3)
         expected_flux = Dr16ExpectedFlux(config["expected flux"])
 
         self.assertTrue(isinstance(expected_flux.get_eta, interp1d))
         self.assertTrue(isinstance(expected_flux.get_fudge, interp1d))
         self.assertTrue(isinstance(expected_flux.get_mean_cont, interp1d))
         self.assertTrue(isinstance(expected_flux.get_var_lss, interp1d))
-        self.assertTrue(expected_flux.lambda_ is None)
-        self.assertTrue(isinstance(expected_flux.log_lambda, np.ndarray))
+        self.assertTrue(isinstance(expected_flux.log_lambda_var_func_grid, np.ndarray))
 
         # setup Forest variables; case: linear wavelength solution
         reset_forest()
@@ -99,13 +98,12 @@ def test_dr16_expected_flux(self):
         self.assertTrue(isinstance(expected_flux.get_fudge, interp1d))
         self.assertTrue(isinstance(expected_flux.get_mean_cont, interp1d))
         self.assertTrue(isinstance(expected_flux.get_var_lss, interp1d))
-        self.assertTrue(isinstance(expected_flux.lambda_, np.ndarray))
-        self.assertTrue(expected_flux.log_lambda is None)
+        self.assertTrue(isinstance(expected_flux.log_lambda_var_func_grid, np.ndarray))
 
     def test_dr16_expected_flux_compute_continuum(self):
         """Test method compute_continuum for class Dr16ExpectedFlux"""
         # setup Forest variables; case: logarithmic wavelength solution
-        setup_forest("log")
+        setup_forest("log", rebin=3)
 
         out_file = f"{THIS_DIR}/results/continua_log.txt"
         test_file = f"{THIS_DIR}/data/continua_log.txt"
@@ -178,7 +176,7 @@ def test_dr16_expected_flux_compute_continuum(self):
     def test_dr16_expected_flux_compute_delta_stack(self):
         """Test method compute_delta_stack for class Dr16ExpectedFlux"""
         # setup Forest variables; case: logarithmic wavelength solution
-        setup_forest("log")
+        setup_forest("log", rebin=3)
 
         out_file = f"{THIS_DIR}/results/delta_stack_log.txt"
         test_file = f"{THIS_DIR}/data/delta_stack_log.txt"
@@ -238,7 +236,7 @@ def test_dr16_expected_flux_compute_delta_stack(self):
     def test_dr16_expected_flux_compute_expected_flux(self):
         """Test method compute_var_stats for class Dr16ExpectedFlux"""
         # setup Forest variables; case: logarithmic wavelength solution
-        setup_forest("log")
+        setup_forest("log", rebin=3)
 
         out_file = f"{THIS_DIR}/results/Log/iter_out_prefix_compute_expected_flux_log.fits.gz"
         test_file = f"{THIS_DIR}/data/iter_out_prefix_compute_expected_flux_log.fits.gz"
@@ -288,7 +286,7 @@ def test_dr16_expected_flux_compute_mean_cont_lin(self):
     def test_dr16_expected_flux_compute_mean_cont_log(self):
         """Test method compute_mean_cont_log for class Dr16ExpectedFlux"""
         # setup Forest variables; case: logarithmic wavelength solution
-        setup_forest("log")
+        setup_forest("log", rebin=3)
 
         out_file = f"{THIS_DIR}/results/mean_cont_log.txt"
         test_file = f"{THIS_DIR}/data/mean_cont_log.txt"
@@ -316,7 +314,7 @@ def test_dr16_expected_flux_compute_mean_cont_log(self):
             expected_flux.compute_continuum(forest)
 
         # compute mean quasar continuum
-        expected_flux.compute_mean_cont_log(data.forests)
+        expected_flux.compute_mean_cont(data.forests)
 
         # save results
         f = open(out_file, "w")
@@ -347,7 +345,7 @@ def test_dr16_expected_flux_compute_mean_cont_log(self):
     def test_dr16_expected_flux_compute_var_stats(self):
         """Test method compute_var_stats for class Dr16ExpectedFlux"""
         # setup Forest variables; case: logarithmic wavelength solution
-        setup_forest("log")
+        setup_forest("log", rebin=3)
 
         test_file = f"{THIS_DIR}/data/eta_var_lss_fudge.txt"
 
@@ -395,7 +393,7 @@ def test_dr16_expected_flux_compute_var_stats(self):
     def test_dr16_expected_flux_populate_los_ids(self):
         """Test method populate_los_ids for class Dr16ExpectedFlux"""
         # setup Forest variables; case: logarithmic wavelength solution
-        setup_forest("log")
+        setup_forest("log", rebin=3)
 
         # initialize Data and Dr16ExpectedFlux instances
         config = ConfigParser()
@@ -432,7 +430,7 @@ def test_dr16_expected_flux_populate_los_ids(self):
     def test_dr16_expected_flux_save_iteration_step(self):
         """Test method save_iteration_step for class Dr16ExpectedFlux"""
         # setup Forest variables; case: logarithmic wavelength solution
-        setup_forest("log")
+        setup_forest("log", rebin=3)
 
         out_file = f"{THIS_DIR}/results/Log/iter_out_prefix_log_iteration1.fits.gz"
         out_file2 = f"{THIS_DIR}/results/Log/iter_out_prefix_log.fits.gz"
diff --git a/py/picca/tests/delta_extraction/test_utils.py b/py/picca/tests/delta_extraction/test_utils.py
index 33564b6f4..b045356ad 100644
--- a/py/picca/tests/delta_extraction/test_utils.py
+++ b/py/picca/tests/delta_extraction/test_utils.py
@@ -13,14 +13,11 @@
 def reset_forest():
     """Reset the class variables of Forest and Pk1dForest"""
     Forest.wave_solution = None
-    Forest.lambda_grid = None
-    Forest.lambda_rest_frame_grid = None
     Forest.log_lambda_grid = None
     Forest.log_lambda_rest_frame_grid = None
     Forest.mask_fields = []
     Pk1dForest.lambda_abs_igm = None
 
-
 # setup Forest class variables
 def setup_forest(wave_solution, rebin=1):
     """Set Forest class variables
@@ -39,7 +36,7 @@ def setup_forest(wave_solution, rebin=1):
     if wave_solution == "log":
         pixel_step = 1e-4 * rebin
     elif wave_solution == "lin":
-        pixel_step = 1.
+        pixel_step = 1.0
     pixel_step_rf = pixel_step
 
     Forest.set_class_variables(3600.0, 5500.0, 1040.0, 1200.0, pixel_step, pixel_step_rf,