diff --git a/bxa/sherpa/background/pca.py b/bxa/sherpa/background/pca.py
index fbf07f2..84d3cc7 100644
--- a/bxa/sherpa/background/pca.py
+++ b/bxa/sherpa/background/pca.py
@@ -14,7 +14,7 @@
 	from sherpa.models.parameter import Parameter
 	from sherpa.models import ArithmeticModel, CompositeModel
 	from sherpa.astro.ui import *
-	from sherpa.astro.instrument import RSPModelNoPHA, RMFModelNoPHA
+	from sherpa.astro.instrument import RSPModelNoPHA, RMFModelNoPHA, create_delta_rmf, create_arf
 else:
 	# mock objects when sherpa doc is built
 	ArithmeticModel = object
@@ -103,13 +103,11 @@ def calc(self, p, x, xhi=None, *args, **kwargs):
 
 
 def get_identity_response(i):
-	n = get_bkg(i).counts.size
 	rmf = get_rmf(i)
-	try:
-		arf = get_arf(i)
-		return lambda model: IdentityResponse(n, model, arf=arf, rmf=rmf)
-	except:
-		return lambda model: IdentityRMF(n, model, rmf=rmf)
+	delta_rmf = create_delta_rmf(rmf.e_min, rmf.e_max, offset=1, e_min=rmf.e_min, e_max=rmf.e_max, ethresh=rmf.ethresh)
+	flat_arf = create_arf(rmf.e_min, rmf.e_max, ethresh=rmf.ethresh)
+
+	return lambda model: RSPModelNoPHA(flat_arf, delta_rmf, model)
 
 
 
@@ -250,7 +248,7 @@ def calc_bkg_stat(self):
 		assert len(ss) == 1
 		return ss[0].statval
 
-	def fit(self):
+	def fit(self, max_lines=10):
 		# try a PCA decomposition of this spectrum
 		logf.info('fitting background of ID=%s using PCA method' % (self.id))
 		initial = self.decompose()
@@ -347,9 +345,9 @@ def fit(self):
 			p.val = v
 
 		last_model = convbkgmodel
-		for i in range(10):
+		for i in range(1, max_lines + 1):
 			print()
-			print('Adding Gaussian#%d' % (i+1))
+			print('Adding Gaussian#%d' % (i))
 			# find largest discrepancy
 			set_analysis(id, "ener", "rate")
 			m = get_bkg_fit_plot(id)
@@ -362,12 +360,12 @@ def fit(self):
 			y = m.dataplot.y.cumsum()
 			z = m.modelplot.y.cumsum()
 			diff = numpy.abs(y - z)
-			i = numpy.argmax(diff)
+			imax = numpy.argmax(diff)
 			energies = x
-			e = x[i]
-			print('largest remaining discrepancy at %.3fkeV[%d], need %d counts' % (x[i], i, diff[i]))
-			#e = x[i]
-			power = diff_rate[i]
+			e = x[imax]
+			print('largest remaining discrepancy at %.3fkeV[%d], need %d counts' % (x[imax], imax, diff[imax]))
+			power = diff_rate[imax]
+			# power = diff[imax]
 			# lets try to inject a gaussian there
 
 			g = xsgaussian('g_%d_%d' % (id, i))
@@ -376,12 +374,13 @@ def fit(self):
 			g.LineE.min = energies[0]
 			g.LineE.max = energies[-1]
 			g.LineE.val = e
-			if i > len(diff) - 2:
-				i = len(diff) - 2
-			if i < 2:
-				i = 2
-			g.Sigma = (x[i + 1] - x[i - 1])
-			g.Sigma.min = (x[i + 1] - x[i - 1])/3
+			if imax > len(diff) - 2:
+				imax = len(diff) - 2
+			if imax < 2:
+				imax = 2
+			g.Sigma = (x[imax + 1] - x[imax - 1])
+			g.Sigma.min = (x[imax + 1] - x[imax - 1])/3
+			g.Sigma.hard_max = 1e10
 			g.Sigma.max = x[-1] - x[0]
 			g.norm.min = power * 1e-6
 			g.norm.val = power
@@ -405,7 +404,7 @@ def fit(self):
 					p.val = v
 				break
 
-def auto_background(id):
+def auto_background(id, max_lines=10):
 	"""Automatically fits background *id* based on PCA-based templates,
 	and additional gaussian lines as needed by AIC."""
 	bkgmodel = PCAFitter(id)
@@ -413,7 +412,7 @@ def auto_background(id):
 	prev_level = log_sherpa.level
 	try:
 		log_sherpa.setLevel(logging.WARN)
-		bkgmodel.fit()
+		bkgmodel.fit(max_lines)
 	finally:
 		log_sherpa.setLevel(prev_level)
 	m = get_bkg_fit_plot(id)
diff --git a/examples/sherpa/chandra/179_pcabkg.json b/examples/sherpa/chandra/179_pcabkg.json
new file mode 100644
index 0000000..fb6394d
--- /dev/null
+++ b/examples/sherpa/chandra/179_pcabkg.json
@@ -0,0 +1,13 @@
+{
+  "lognorm": 3.9371033003777676,
+  "PC1": 0.00792083863859961,
+  "PC2": 0.013572359127513619,
+  "PC3": 0.003391010057753655,
+  "PC4": 0.0033947454403729987,
+  "PC5": -0.013415970033646883,
+  "PC6": 0.005511275689439748,
+  "PC7": 0.00018728120427328259,
+  "PC8": 0.0006131251304566246,
+  "PC9": -0.0007613235551938485,
+  "PC10": -0.0006913655223552859
+}
\ No newline at end of file
diff --git a/examples/sherpa/swift/interval0pc_pcabkg.json b/examples/sherpa/swift/interval0pc_pcabkg.json
new file mode 100644
index 0000000..9aaa19d
--- /dev/null
+++ b/examples/sherpa/swift/interval0pc_pcabkg.json
@@ -0,0 +1,13 @@
+{
+  "lognorm": 2.3985433250317354,
+  "PC1": -0.07305826667314888,
+  "PC2": -0.08753541499291051,
+  "PC3": 0.0,
+  "PC4": 0.0,
+  "PC5": 0.0,
+  "PC6": 0.0,
+  "PC7": 0.0,
+  "PC8": 0.0,
+  "PC9": 0.0,
+  "PC10": 0.0
+}
\ No newline at end of file
diff --git a/examples/sherpa/test_auto_background.py b/examples/sherpa/test_auto_background.py
new file mode 100644
index 0000000..19788b4
--- /dev/null
+++ b/examples/sherpa/test_auto_background.py
@@ -0,0 +1,271 @@
+import json
+import numpy as np
+import sherpa.astro.ui as shp
+
+from bxa.sherpa.background.pca import auto_background
+
+
+def test_chandra_pcabkg():
+    _reset_sherpa()
+    _load_chandra_data(emin=0.5, emax=8)
+
+    filename = "chandra/179_pcabkg.json"
+    bkgmodel = auto_background(1, max_lines=0)
+    _test_pca_params(filename, bkgmodel)
+    # _save_pca_params(bkgmodel, filename)
+
+def test_chandra_detect_line():
+    _reset_sherpa()
+    bkgdata = _load_chandra_data(emin=0.5, emax=8)
+    channel = _inject_count_excess(bkgdata)
+    bkgmodel = auto_background(1, max_lines=1)
+
+    _check_injected_line(bkgdata, bkgmodel, channel)
+
+def _load_chandra_data(emin, emax):
+    shp.load_pha("chandra/179.pi")
+    _ungroup_and_ignore_bkg(emin, emax)
+
+    return shp.get_bkg()
+
+
+def test_swift_pcabkg():
+    _reset_sherpa()
+    _load_swift_data(emin=0.5, emax=5.0)
+
+    filename = "swift/interval0pc_pcabkg.json"
+    bkgmodel = auto_background(1, max_lines=0)
+    _test_pca_params(filename, bkgmodel)
+    # _save_pca_params(bkgmodel, filename)
+
+def test_swift_detect_line():
+    _reset_sherpa()
+    bkgdata = _load_swift_data(emin=0.5, emax=5.0)
+    channel = _inject_count_excess(bkgdata)
+    bkgmodel = auto_background(1, max_lines=1)
+
+    _check_injected_line(bkgdata, bkgmodel, channel)
+
+def _load_swift_data(emin, emax):
+    shp.load_pha("swift/interval0pc.pi")
+    _ungroup_and_ignore_bkg(emin, emax)
+
+    return shp.get_bkg()
+
+
+def test_xmmpn_pcabkg():
+    _reset_sherpa()
+    _load_xmmpn_data(emin=0.2, emax=10.0)
+
+    filename = "xmm/pn_pcabkg.json"
+    bkgmodel = auto_background(1, max_lines=0)
+    _test_pca_params(filename, bkgmodel)
+    # _save_pca_params(bkgmodel, filename)
+
+def test_xmmpn_detect_line():
+    _reset_sherpa()
+    bkgdata = _load_xmmpn_data(emin=0.2, emax=10.0)
+    channel = _inject_count_excess(bkgdata)
+    bkgmodel = auto_background(1, max_lines=1)
+
+    _check_injected_line(bkgdata, bkgmodel, channel)
+
+def _load_xmmpn_data(emin, emax):
+    shp.load_pha("xmm/pn_spec.fits")
+    shp.load_bkg("xmm/pn_backspec.fits")
+    shp.load_rmf("xmm/pn.rmf")
+    shp.load_arf("xmm/pn.arf")
+    shp.load_bkg_rmf("xmm/pn.rmf")
+    shp.load_bkg_arf("xmm/pn.arf")
+    _ungroup_and_ignore_bkg(emin, emax)
+
+    return shp.get_bkg()
+
+
+def test_xmmmos_pcabkg():
+    _reset_sherpa()
+    _load_xmmmos_data(emin=0.2, emax=10.0)
+
+    filename = "xmm/mos_pcabkg.json"
+    bkgmodel = auto_background(1, max_lines=0)
+    _test_pca_params(filename, bkgmodel)
+    # _save_pca_params(bkgmodel, filename)
+
+def test_xmmmos_detect_line():
+    _reset_sherpa()
+    bkgdata = _load_xmmmos_data(emin=0.2, emax=10.0)
+    channel = _inject_count_excess(bkgdata)
+    bkgmodel = auto_background(1, max_lines=1)
+
+    _check_injected_line(bkgdata, bkgmodel, channel)
+
+def _load_xmmmos_data(emin, emax):
+    shp.load_pha("xmm/mos_spec.fits")
+    shp.load_bkg("xmm/mos_backspec.fits")
+    shp.load_rmf("xmm/mos.rmf")
+    shp.load_arf("xmm/mos.arf")
+    shp.load_bkg_rmf("xmm/mos.rmf")
+    shp.load_bkg_arf("xmm/mos.arf")
+    _ungroup_and_ignore_bkg(emin, emax)
+
+    return shp.get_bkg()
+
+
+def test_nustar_pcabkg():
+    _reset_sherpa()
+    _load_nustar_data(emin=5.0, emax=77.0)
+
+    filename = "nustar/nu60360003002A01_pcabkg.json"
+    bkgmodel = auto_background(1, max_lines=0)
+    _test_pca_params(filename, bkgmodel)
+    # _save_pca_params(bkgmodel, filename)
+
+def test_nustar_detect_line():
+    _reset_sherpa()
+    bkgdata = _load_nustar_data(emin=5.0, emax=77.0)
+    channel = _inject_count_excess(bkgdata)
+    bkgmodel = auto_background(1, max_lines=1)
+
+    _check_injected_line(bkgdata, bkgmodel, channel)
+
+def _load_nustar_data(emin, emax):
+    shp.load_pha("nustar/nu60360003002A01_sr_grp.pha")
+    _ungroup_and_ignore_bkg(emin, emax)
+
+    return shp.get_bkg()
+
+
+def test_erosita_pcabkg():
+    _reset_sherpa()
+    _load_erosita_data(emin=0.2, emax=8.0)
+
+    filename = "erosita/em01_182117_020_pcabkg.json"
+    bkgmodel = auto_background(1, max_lines=0)
+    _test_pca_params(filename, bkgmodel)
+    # _save_pca_params(bkgmodel, filename)
+
+def test_erosita_detect_line():
+    _reset_sherpa()
+    bkgdata = _load_erosita_data(emin=0.2, emax=8.0)
+    channel = _inject_count_excess(bkgdata)
+    bkgmodel = auto_background(1, max_lines=1)
+
+    _check_injected_line(bkgdata, bkgmodel, channel)
+
+def _load_erosita_data(emin, emax):
+    shp.load_pha("erosita/em01_182117_020_SourceSpec_00005_c010.fits")
+    _ungroup_and_ignore_bkg(emin, emax)
+
+    return shp.get_bkg()
+
+
+def test_suzaku_pcabkg():
+    _reset_sherpa()
+    _load_suzaku_data(emin=0.2, emax=8.0)
+
+    filename = "suzaku/ae900001010xi0_0_3x3n000a_pcabkg.json"
+    bkgmodel = auto_background(1, max_lines=0)
+    _test_pca_params(filename, bkgmodel)
+    # _save_pca_params(bkgmodel, filename)
+
+def test_suzaku_detect_line():
+    _reset_sherpa()
+    bkgdata = _load_suzaku_data(emin=0.2, emax=8.0)
+    channel = _inject_count_excess(bkgdata)
+    bkgmodel = auto_background(1, max_lines=1)
+
+    _check_injected_line(bkgdata, bkgmodel, channel)
+
+def _load_suzaku_data(emin, emax):
+    shp.load_pha("suzaku/ae900001010xi0_0_3x3n000a_sr.pi")
+    _ungroup_and_ignore_bkg(emin, emax)
+
+    return shp.get_bkg()
+
+
+def _reset_sherpa():
+    shp.clean()
+    shp.set_stat("cstat")
+
+def _ungroup_and_ignore_bkg(emin, emax):
+    shp.ungroup(bkg_id=1)
+    shp.notice(bkg_id=1)
+    shp.ignore(lo=None, hi=emin, bkg_id=1)
+    shp.ignore(lo=emax, hi=None, bkg_id=1)
+    
+    # Set dummy source model for auto_background
+    shp.set_model("powlaw1d.dummy")
+
+
+def _save_pca_params(bkgmodel, filename):
+    pca_params = {p.name: p.val for p in bkgmodel.pars}
+
+    with open(filename, "w") as fp:
+        json.dump(pca_params, fp, indent=2)
+
+def _load_pca_params(filename):
+    with open(filename, "r") as fp:
+        params = json.load(fp)
+
+    return params
+
+
+def _test_pca_params(filename, bkgmodel):
+    test_pca_params = _load_pca_params(filename)
+    bkgdata = shp.get_bkg()
+
+    assert len(bkgdata.channel) == bkgmodel.rmf.detchans
+    assert len(bkgmodel.pars) == len(test_pca_params)
+    
+    for p in bkgmodel.pars:
+        assert p.name in test_pca_params
+        assert np.isclose(p.val, test_pca_params[p.name])
+
+
+def _inject_count_excess(bkgdata):
+    # Insert a high excess of counts in the middle
+    # of the noticed energy range
+    cmin = bkgdata.channel[bkgdata.mask][0]
+    cmax = bkgdata.channel[bkgdata.mask][-1]
+    c = int((cmin + cmax) / 2)
+    c0 = int(bkgdata.channel[0])
+    bkgdata.counts[c - c0] = 100 * bkgdata.counts.max()
+    print(f"Injecting {bkgdata.counts[c - c0]:.0f} counts at {bkgdata.x[c - c0]:.2f} keV [channel {c}]")
+
+    return c
+
+def _check_injected_line(bkgdata, bkgmodel, channel_injected):
+    for p in bkgmodel.pars:
+        if p.name == "LineE":
+            first_line_energy = p.val
+            break
+
+    first_line_channel = bkgdata.channel[0] + np.argmin(np.absolute(bkgdata.x - first_line_energy))
+    assert int(first_line_channel) == channel_injected
+
+
+
+def run_tests():
+    instruments = [
+        "chandra",
+        "swift",
+        "xmmpn",
+        "xmmmos",
+        # "nustar",
+        # "erosita",
+        # "suzaku",
+        # "rxtepca",
+        # "rxtehexa",
+        # "rxtehexb",
+    ]
+
+    for instrument in instruments:
+        test_continuum = globals()[f"test_{instrument}_pcabkg"]
+        test_continuum()
+
+        test_lines = globals()[f"test_{instrument}_detect_line"]
+        test_lines()
+
+
+if __name__ == "__main__":
+    run_tests()
diff --git a/examples/sherpa/xmm/mos_pcabkg.json b/examples/sherpa/xmm/mos_pcabkg.json
new file mode 100644
index 0000000..a0beb9b
--- /dev/null
+++ b/examples/sherpa/xmm/mos_pcabkg.json
@@ -0,0 +1,8 @@
+{
+  "lognorm": 2.4530946692283666,
+  "PC1": -0.0038685903834438174,
+  "PC2": 0.0016748138563737666,
+  "PC3": -0.00028983876139867214,
+  "PC4": -0.0037276618440139284,
+  "PC5": 0.00329618452813025
+}
\ No newline at end of file
diff --git a/examples/sherpa/xmm/pn_pcabkg.json b/examples/sherpa/xmm/pn_pcabkg.json
new file mode 100644
index 0000000..a509f28
--- /dev/null
+++ b/examples/sherpa/xmm/pn_pcabkg.json
@@ -0,0 +1,8 @@
+{
+  "lognorm": 3.1704175359092646,
+  "PC1": -0.0026390052954809416,
+  "PC2": 0.0017207083291831816,
+  "PC3": 0.0021540150804401806,
+  "PC4": 0.0,
+  "PC5": 0.0
+}
\ No newline at end of file