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 6e8e8e6c1..23395560e 100644
--- a/py/picca/delta_extraction/expected_fluxes/dr16_expected_flux.py
+++ b/py/picca/delta_extraction/expected_fluxes/dr16_expected_flux.py
@@ -16,7 +16,7 @@
 accepted_options = [
     "iter out prefix", "limit eta", "limit var lss", "min num qso in fit",
     "num bins variance", "num iterations", "num processors", "order", "out dir",
-    "use constant weight", "use ivar as weight"
+    "use constant weight", "use ivar as weight", "force stack delta to zero"
 ]
 
 defaults = {
@@ -29,6 +29,7 @@
     "order": 1,
     "use constant weight": False,
     "use ivar as weight": False,
+    "force stack delta to zero": True
 }
 
 
@@ -133,6 +134,9 @@ class Dr16ExpectedFlux(ExpectedFlux):
 
     use_ivar_as_weight: boolean
     If "True", use ivar as weights (implemented as eta = 1, sigma_lss = fudge = 0).
+
+    force_stack_delta_to_zero: boolean
+    If "True", continuum is corrected by stack_delta.
     """
 
     def __init__(self, config):
@@ -156,6 +160,7 @@ def __init__(self, config):
         self.order = None
         self.use_constant_weight = None
         self.use_ivar_as_weight = None
+        self.force_stack_delta_to_zero = None
         self.__parse_config(config)
 
         # initialize variables
@@ -376,6 +381,12 @@ def __parse_config(self, config):
                 "Missing argument 'use ivar as weight' required by Dr16ExpectedFlux"
             )
 
+        self.force_stack_delta_to_zero = config.getboolean("force stack delta to zero")
+        if self.force_stack_delta_to_zero is None:
+            raise ExpectedFluxError(
+                "Missing argument 'force stack delta to zero' required by Dr16ExpectedFlux"
+            )
+
     # this should be a read-only function as it is called in a parallelized way
     # TODO: consider making this not a function to minimize future bugs
     def compute_continuum(self, forest):
@@ -955,10 +966,15 @@ def populate_los_ids(self, forests):
             if forest.bad_continuum_reason is not None:
                 continue
             # get the variance functions and statistics
-            stack_delta = self.get_stack_delta(forest.log_lambda)
             eta = self.get_eta(forest.log_lambda)
 
-            mean_expected_flux = forest.continuum * stack_delta
+            # assignment operator (=) creates a reference, such that
+            # mean_expected_flux points to forest.continuum and
+            # forest.continuum gets modified within if statement
+            mean_expected_flux = np.copy(forest.continuum)
+            if self.force_stack_delta_to_zero:
+                stack_delta = self.get_stack_delta(forest.log_lambda)
+                mean_expected_flux *= stack_delta
             weights = self.get_continuum_weights(forest, mean_expected_flux)
             # this is needed as the weights from get_continuum_weights are
             # divided by the continuum model squared, in this case mean_expected_flux
diff --git a/py/picca/tests/delta_extraction/data/.config.ini b/py/picca/tests/delta_extraction/data/.config.ini
index b3b30747f..a1765a7e6 100644
--- a/py/picca/tests/delta_extraction/data/.config.ini
+++ b/py/picca/tests/delta_extraction/data/.config.ini
@@ -56,6 +56,7 @@ use ivar as weight = False
 out dir = /Users/iperezra/software/picca/py/picca/tests/delta_extraction/results/config_tests/
 num processors = 1
 min num qso in fit = 100
+force stack delta to zero = True
 
 [correction arguments 0]
 filename = /Users/iperezra/software/picca/py/picca/tests/delta_extraction/data/delta_attributes.fits.gz
diff --git a/tutorials/delta_extraction/picca_delta_extraction_configuration_tutorial.ipynb b/tutorials/delta_extraction/picca_delta_extraction_configuration_tutorial.ipynb
index 4b3014c41..b613bc7c4 100644
--- a/tutorials/delta_extraction/picca_delta_extraction_configuration_tutorial.ipynb
+++ b/tutorials/delta_extraction/picca_delta_extraction_configuration_tutorial.ipynb
@@ -298,6 +298,7 @@
     "- `order`: Order of the polynomial for the continuum fit. **Type: int, Required: no, Default: 1**\n",
     "- `use constant weight`: If \"True\", set all the delta weights to one (implemented as eta = 0, sigma_lss = 1, fudge = 0) **Type: bool, Required: no, Default: False**\n",
     "- `use ivar as weight`: If \"True\", use ivar as weights (implemented as eta = 1, sigma_lss = fudge = 0) **Type: bool, Required: no, Default: False**\n",
+    "- `force stack delta to zero`: If \"True\" (Default), continuum is corrected by stack_delta. **Type: bool, Required: no, Default: True**\n",
     "\n",
     "### TrueContinuum\n",
     "- `input directory`: Directory to spectra files. **Type: str, Required: yes**\n",