Fix: Rema has to be mirrored in both dimensions

README.md

@@ -61,7 +61,7 @@ Reconstruct incident_spectrum based on observed_spectrum using the supplied dete
 positional arguments:
   matrixfile            container file containing detector response matrix
   observed_spectrum     file containing the observed spectrum
-  incident_spectrum     write trace of incident spectrum to this path (.nc
+  incident_spectrum     write trace of incident spectrum to this path (.h5
                         file)
 
 options:

src/boris/app.py

@@ -30,6 +30,7 @@
 from boris.utils import (
     create_matrix,
     get_rema,
+    mult_rema,
     one_sigma,
     read_rebin_spectrum,
     QuantityExtractor,
@@ -206,12 +207,12 @@ def boris(
         if fit is None:
             from boris.core import fit
         trace = fit(
-            rema.values(),
+            rema.values()[::-1, ::-1],
             spectrum.values(),
             rema.axes[0].edges,
             background.values() if background else None,
             background_scale,
-            rema_alt.values() if rema_alt else None,
+            rema_alt.values()[::-1, ::-1] if rema_alt else None,
             fit_beam,
             **kwargs,
         )
@@ -320,7 +321,7 @@ def sirob(
 
     with do_step("Calculating observed (convoluted) spectrum"):
         observed = hist.Hist(rema.axes[0], storage=hist.storage.Double())
-        observed.values()[:] = incident.values() @ rema.values()
+        observed.values()[:] = mult_rema(rema, incident)
         if background is not None:
             observed.values()[:] += background_scale * background.values()
 

src/boris/boris_app.py

@@ -159,7 +159,7 @@ def boris_app():
     )
     parser.add_argument(
         "incident_spectrum",
-        help="write trace of incident spectrum to this path (.nc file)",
+        help="write trace of incident spectrum to this path (.h5 file)",
         type=Path,
     )
 

src/boris/core.py

@@ -39,7 +39,7 @@ def theuerkauf_norm(sigma, tl):
     Normalization factor for Theuerkauf peak shape
     """
     return 1 / (
-        (sigma ** 2) / tl * pt.exp(-(tl * tl) / (2.0 * sigma ** 2))
+        (sigma**2) / tl * pt.exp(-(tl * tl) / (2.0 * sigma**2))
         + pt.sqrt(np.pi / 2.0)
         * sigma
         * (1 + pt.erf(tl / (pt.sqrt(2.0) * sigma)))
@@ -58,8 +58,8 @@ def theuerkauf(x, pos, vol, sigma, tl):
     norm = theuerkauf_norm(sigma, tl)
     _x = pt.switch(
         dx < -tl,
-        tl / (sigma ** 2) * (dx + tl / 2.0),
-        -dx * dx / (2.0 * sigma ** 2),
+        tl / (sigma**2) * (dx + tl / 2.0),
+        -dx * dx / (2.0 * sigma**2),
     )
     return vol * norm * pt.exp(_x)
 
@@ -117,15 +117,18 @@ def fit(
             )
         else:
             incident = pm.HalfFlat("incident", shape=spectrum.shape[0])
+
         if rema_alt is None:
             rema_eff = rema
-            rema_eff_diag = np.diag(rema)
+            rema_eff_diag = np.diag(rema_eff)
         else:
             interpolation = pm.Uniform("interpolation", 0, 1, shape=())
             rema_eff = (1 - interpolation) * rema + interpolation * rema_alt
             rema_eff_diag = (1 - interpolation) * np.diag(
                 rema
             ) + interpolation * np.diag(rema_alt)
+
+        # rema_eff = pm.ConstantData("rema", rema_eff)
         folded = pm.Deterministic("folded", incident @ rema_eff)
         if background is None:
             folded_plus_bg = folded
@@ -211,7 +214,6 @@ def fit(
             ndraws,
             step=step,
             start=start,
-            return_inferencedata=True,
             idata_kwargs=dict(log_likelihood=True),
             **kwargs,
         )

src/boris/utils.py

@@ -486,3 +486,28 @@ def get_quantities(
         ).T
 
     return res
+
+
+def mult_rema(
+    rema: hist.hist.Hist, incident_spec: hist.hist.Hist
+) -> np.ndarray:
+    """
+    Apply a response matrix to an incident spectrum
+
+    :param rema: Detector response matrix
+    :param spec: Incident spectrum
+    """
+    return rema.values()[::-1, ::-1] @ incident_spec.values()
+
+
+def mult_rema_inv(
+    rema: hist.hist.Hist, observed_spec: hist.hist.Hist
+) -> np.ndarray:
+    """
+    Apply the inverse of a response matrix to an observed spectrum
+
+    :param rema: Detector response matrix
+    :param spec: Observed spectrum
+    """
+    rema_inv = np.linalg.inv(rema.values()[::-1, ::-1])
+    return observed_spec.values() @ rema_inv

tests/test_app.py

@@ -42,6 +42,7 @@
     sirob,
 )
 from boris.io import write_specs, read_spectrum
+from boris.utils import mult_rema
 
 from tests.helpers.utils import create_simulations
 
@@ -96,9 +97,7 @@ def test_sirob(rema, incident, tmp_path):
     assert observed.ndim == 1
     assert observed.shape[0] == 10
     assert np.isclose(observed.axes[0].edges, incident.axes[0].edges).all()
-    assert np.isclose(
-        incident.values() @ rema.values(), observed.values()
-    ).all()
+    assert np.isclose(mult_rema(rema, incident), observed.values()).all()
 
     background = hist.Hist.new.Regular(10, 2000.0, 2200.0).Int64(
         data=np.random.uniform(10, 100, size=10).astype(np.int64)
@@ -122,7 +121,7 @@ def test_sirob(rema, incident, tmp_path):
     assert observed_bg.shape[0] == 10
     assert np.isclose(observed_bg.axes[0].edges, incident.axes[0].edges).all()
     assert np.isclose(
-        incident.values() @ rema.values() + 2.0 * background.values(),
+        mult_rema(rema, incident) + 2.0 * background.values(),
         observed_bg.values(),
     ).all()
 
@@ -150,7 +149,7 @@ def test_do_step():
 
 def test_boris(rema, incident, tmp_path):
     observed = hist.Hist.new.Regular(10, 2000.0, 2200.0).Int64(
-        data=(incident.values() @ rema.values()).astype(np.int64)
+        data=mult_rema(rema, incident).astype(np.int64)
     )
 
     write_specs(tmp_path / "rema.npz", {"rema": rema})
@@ -173,7 +172,7 @@ def test_boris(rema, incident, tmp_path):
         data=np.random.uniform(10, 100, size=10).astype(np.int64)
     )
     observed_bg = hist.Hist.new.Regular(10, 2000.0, 2200.0).Int64(
-        data=(incident.values() @ rema.values() + background.values())
+        data=(mult_rema(rema, incident) + background.values())
     )
     write_specs(tmp_path / "observed_bg.npz", {"observed_bg": observed_bg})
     write_specs(tmp_path / "background.npz", {"background": background})

tests/test_boris_app.py

@@ -26,6 +26,8 @@
 from boris.io import write_specs
 from boris.boris_app import boris_app
 
+from boris.utils import mult_rema
+
 
 def test_boris_app(tmp_path):
     sys.argv = [
@@ -63,12 +65,10 @@ def test_boris_app(tmp_path):
         data=np.random.uniform(10, 100, size=10).astype(np.int64)
     )
     observed_wobg = hist.Hist.new.Regular(10, 2000, 2200).Int64(
-        data=(incident.values() @ rema.values()).astype(np.int64)
+        data=mult_rema(rema, incident).astype(np.int64)
     )
     observed = hist.Hist.new.Regular(10, 2000, 2200).Int64(
-        data=(background.values() + incident.values() @ rema.values()).astype(
-            np.int64
-        )
+        data=mult_rema(rema, background + incident).astype(np.int64)
     )
 
     write_specs(tmp_path / "rema.npz", {"rema": rema})