Magic-simplified ER microphysics modeling

appletree/components/er.py

@@ -18,6 +18,20 @@ def __init__(self, *args, **kwargs):
         self.register_all(apt.plugins.efficiency)
 
 
+class ERBandSR(ComponentSim):
+    norm_type = "on_pdf"
+    add_eps_to_hist = False
+
+    def __init__(self, *args, **kwargs):
+        super().__init__(*args, **kwargs)
+        self.register(apt.plugins.common.UniformEnergySpectra)
+        self.register(apt.plugins.common.PositionSpectra)
+        self.register_all(apt.plugins.er_sr_microphys)
+        self.register_all(apt.plugins.detector)
+        self.register_all(apt.plugins.reconstruction)
+        self.register_all(apt.plugins.efficiency)
+
+
 class ERPeak(ComponentSim):
     norm_type = "on_pdf"
     add_eps_to_hist = False

appletree/parameters/er_sr.json

@@ -0,0 +1,229 @@
+{
+    "w": {
+        "prior_type": "norm",
+        "prior_args": {
+            "mean": 0.0137,
+            "std": 0.0002
+        },
+        "allowed_range": [
+            0,
+            1.0
+        ],
+        "init_mean": 0.0137,
+        "init_std": 0.0002,
+        "unit": "keV",
+        "doc": "Mean energy to generate a quanta in liquid xenon"
+    },
+    "fano": {
+        "prior_type": "fixed",
+        "prior_args": {
+            "val": 0.059
+        },
+        "allowed_range": null,
+        "init_mean": null,
+        "init_std": null,
+        "unit": "1",
+        "doc": "Fano factor which describes the fluctuation of num of quanta"
+    },
+    "gas_gain": {
+        "prior_type": "fixed",
+        "prior_args": {
+            "val": 31.3
+        },
+        "allowed_range": null,
+        "init_mean": null,
+        "init_std": null,
+        "unit": "PE/electron",
+        "doc": "Gas gain"
+    },
+    "drift_velocity": {
+        "prior_type": "fixed",
+        "prior_args": {
+            "val": 0.0677
+        },
+        "allowed_range": null,
+        "init_mean": null,
+        "init_std": null,
+        "unit": "cm/us",
+        "doc": "Drift velocity"
+    },
+    "s2_threshold": {
+        "prior_type": "fixed",
+        "prior_args": {
+            "val": 500.0
+        },
+        "allowed_range": null,
+        "init_mean": null,
+        "init_std": null,
+        "unit": "PE",
+        "doc": "S2 threshold"
+    },
+    "g1": {
+        "prior_type": "norm",
+        "prior_args": {
+            "mean": 0.1515,
+            "std": 0.0014
+        },
+        "allowed_range": [
+            0,
+            1.0
+        ],
+        "init_mean": 0.1515,
+        "init_std": 0.0014,
+        "unit": "PE/photon",
+        "doc": "g1"
+    },
+    "g2": {
+        "prior_type": "norm",
+        "prior_args": {
+            "mean": 16.45,
+            "std": 0.64
+        },
+        "allowed_range": [
+            0,
+            100.0
+        ],
+        "init_mean": 16.45,
+        "init_std": 0.64,
+        "unit": "PE/electron",
+        "doc": "g2"
+    },
+    "p_dpe": {
+        "prior_type": "norm",
+        "prior_args": {
+            "mean": 0.227,
+            "std": 0.024
+        },
+        "allowed_range": [
+            0,
+            1.0
+        ],
+        "init_mean": 0.227,
+        "init_std": 0.024,
+        "unit": "1",
+        "doc": "DPE probability"
+    },
+    "elife_sigma": {
+        "prior_type": "norm",
+        "prior_args": {
+            "mean": 0.0,
+            "std": 0.02
+        },
+        "allowed_range": [
+            -1.0,
+            1.0
+        ],
+        "init_mean": 0.0,
+        "init_std": 0.02,
+        "unit": "1",
+        "doc": "elife shifter"
+    },
+    "s1_eff_3f_sigma": {
+        "prior_type": "norm",
+        "prior_args": {
+            "mean": 0.0,
+            "std": 1.0
+        },
+        "allowed_range": [
+            -5.0,
+            5.0
+        ],
+        "init_mean": 0.0,
+        "init_std": 1.0,
+        "unit": "1",
+        "doc": "S1 efficiency shifter"
+    },
+    "s1_cut_acc_sigma": {
+        "prior_type": "norm",
+        "prior_args": {
+            "mean": 0.0,
+            "std": 1.0
+        },
+        "allowed_range": [
+            -5.0,
+            5.0
+        ],
+        "init_mean": 0.0,
+        "init_std": 1.0,
+        "unit": "1",
+        "doc": "S1 cut acceptance shifter"
+    },
+    "s2_cut_acc_sigma": {
+        "prior_type": "norm",
+        "prior_args": {
+            "mean": 0.0,
+            "std": 1.0
+        },
+        "allowed_range": [
+            -5.0,
+            5.0
+        ],
+        "init_mean": 0.0,
+        "init_std": 1.0,
+        "unit": "1",
+        "doc": "S2 cut acceptance shifter"
+    },
+    "p_mu_0": {
+        "prior_type": "free",
+        "prior_args": {},
+        "allowed_range": [
+            0.9,
+            1.1
+        ],
+        "init_mean": 1,
+        "init_std": 0.02,
+        "unit": "au",
+        "doc": "Parameter for Gaussian emission model mean"
+    },
+    "p_mu_1": {
+        "prior_type": "free",
+        "prior_args": {},
+        "allowed_range": [
+            0.1,
+            0.7
+        ],
+        "init_mean": 0.4,
+        "init_std": 0.02,
+        "unit": "au",
+        "doc": "Parameter for Gaussian emission model mean"
+    },
+    "p_sigma_0": {
+        "prior_type": "free",
+        "prior_args": {},
+        "allowed_range": [
+            0.2,
+            0.6
+        ],
+        "init_mean": 0.4,
+        "init_std": 0.02,
+        "unit": "au",
+        "doc": "Parameter for Gaussian emission model std"
+    },
+    "ac_rate": {
+        "prior_type": "norm",
+        "prior_args": {
+            "mean": 10.0,
+            "std": 2.0
+        },
+        "allowed_range": [
+            0,
+            10000000000.0
+        ],
+        "init_mean": 11.2,
+        "init_std": 2.24,
+        "unit": "1",
+        "doc": "ac_rate"
+    },
+    "er_rate": {
+        "prior_type": "free",
+        "prior_args": {},
+        "allowed_range": [
+            0,
+            10000000000.0
+        ],
+        "init_mean": 2000,
+        "init_std": 100,
+        "unit": "1",
+        "doc": "er_rate"
+    }
+}

appletree/plugins/__init__.py

@@ -4,6 +4,9 @@
 from . import er_microphys
 from .er_microphys import *
 
+from . import er_sr_microphys
+from .er_sr_microphys import *
+
 from . import lyqy
 from .lyqy import *
 

appletree/plugins/er_sr_microphys.py

@@ -0,0 +1,102 @@
+from jax import numpy as jnp
+from jax import scipy as jsp
+from jax import jit
+from functools import partial
+
+from appletree import randgen
+from appletree.plugin import Plugin
+from appletree.utils import exporter
+
+
+export, __all__ = exporter(export_self=False)
+
+
+def n_photon_mean_model(num_quanta, p_mu_0=0.9982161, p_mu_1=0.418):
+    """Mean of Gaussian model bridging quanta number to photon number. \begin{aligned} \mu= &
+    p_{\mu_0}\left( \sqrt{x_0}+ \operatorname{asinh}^3\left(\sin ^3.
+
+    \left(\log \left(x_0+p_{\mu_0}\right)\right)\right)-1/0.249\right)^2 \\ &
+    -\operatorname{atan}^2\left(\operatorname{erfc}\left(\sin \left(\log
+    \left(x_0+p_{\mu_0}\right)\right)\right)\right) & x_0 \log \left(p_{\mu_1} \tanh \left(\sin
+    \left(\sin \left(\tan \left( \sinh \left(\sin \left(\log
+    \left(x_0+1\right)\right)\right)\right)\right) \right)\right)+1\right)^3 \end{aligned}
+
+    """
+    result = (
+        p_mu_0
+        * (
+            jnp.sqrt(num_quanta)
+            + jnp.arcsinh(jnp.sin(jnp.log(num_quanta + p_mu_0)) ** 3)
+            - 1 / 0.249
+        )
+        ** 2
+        - jnp.log(jsp.special.erfc(jnp.sin(jnp.log(num_quanta + p_mu_0))) + 1)
+        - (
+            num_quanta
+            * jnp.log(
+                p_mu_1
+                * jnp.tanh(jnp.sin(jnp.sin(jnp.tan(jnp.sinh(jnp.sin(jnp.log(num_quanta + 1)))))))
+                + 1
+            )
+            ** 3
+        )
+    )
+    return result
+
+
+def n_photon_std_model(num_quanta, p_sigma_0=0.4):
+    """Standard deviation of Gaussian model bridging quanta number to photon number.
+
+    \sigma=x_0^{p_{\sigma_0}}-\sinh \left(\cosh \left(\sin \left(
+    \operatorname{asinh}\left(x_0\right)\right)\right)\right)
+
+    """
+    result = num_quanta**p_sigma_0 - jnp.sinh(jnp.cosh(jnp.sin(jnp.arcsinh(num_quanta))))
+    return result
+
+
+@export
+class Quanta(Plugin):
+    depends_on = ["energy"]
+    provides = ["num_quanta"]
+    parameters = (
+        "w",
+        "fano",
+    )
+
+    @partial(jit, static_argnums=(0,))
+    def simulate(self, key, parameters, energy):
+        num_quanta_mean = energy / parameters["w"]
+        num_quanta_std = jnp.sqrt(num_quanta_mean * parameters["fano"])
+        key, num_quanta = randgen.truncate_normal(key, num_quanta_mean, num_quanta_std, vmin=0)
+        return key, num_quanta.round().astype(int)
+
+
+@export
+class Emission(Plugin):
+    """A symbolic-regression-based empirical simulation for microphysics aiming at parameter number
+    reduction.
+
+    Wrapped up IonizationER + mTI + RecombFluct + TrueRecombER + RecombinationER in traditional ER
+    microphysics.
+
+    """
+
+    depends_on = ["num_quanta"]
+    provides = ["num_photon", "num_electron"]
+    parameters = ("p_mu_0", "p_mu_1", "p_sigma_0")
+
+    @partial(jit, static_argnums=(0,))
+    def simulate(self, key, parameters, num_quanta):
+        n_photon_mean = n_photon_mean_model(
+            num_quanta,
+            parameters["p_mu_0"],
+            parameters["p_mu_1"],
+        )
+        n_photon_std = n_photon_std_model(num_quanta, parameters["p_sigma_0"])
+        key, _num_photon = randgen.truncate_normal(
+            key, n_photon_mean, n_photon_std, vmin=0, vmax=num_quanta
+        )
+        num_photon = _num_photon.round().astype(int)
+        num_electron = num_quanta - num_photon
+        return key, num_photon, num_electron