resolves #103 by updating config, eas and cphotang

setup.cfg

@@ -56,6 +56,8 @@ console_scripts =
     nuspacesim = nuspacesim.apps.cli:cli
 
 [options.package_data]
+nuspacesim.dat.CONEX_table =
+    dumpGH_conex_pi_E17_95deg_0km_eposlhc_1394249052_211.dat
 nuspacesim.data.cloud_maps =
     nss_map_CloudTopPressure_01.v0.fits
     nss_map_CloudTopPressure_02.v0.fits

src/nuspacesim/config.py

@@ -316,7 +316,21 @@ def serialize_rad(self, x: float) -> str:
     """ Maximum Azimuthal Angle (Radians). """
     angle_from_limb: float = np.radians(7)
     """ Angle From Limb. Default (Radians). """
-    cherenkov_light_engine: Literal["Default"] = "Default"  # , "CHASM" , "EASCherSim"
+    cherenkov_light_engine: Literal[
+        "Greisen",
+        "Gaisser-Hillas Parameterized",
+        "Gaisser-Hillas Fluctuated",
+        "Default",
+    ] = "Greisen"  # "CHASM", "EASCherSim"
+    """Cherenkov Light Engine model: Default = 'Greisen'"""
+
+    @field_validator("cherenkov_light_engine", mode="before")
+    @classmethod
+    def validate_cherenkov_light_engine(cls, value: str) -> str:
+        if value == "Default":
+            return "Greisen"
+        return value
+
     ionosphere: Optional[Ionosphere] = Ionosphere()
     tau_shower: NuPyPropShower = NuPyPropShower()
     """ Tau Shower Generator. """
@@ -375,16 +389,22 @@ def config_from_fits(filename: str) -> NssConfig:
     def v(key: str):
         fullkey = "Config " + key
         if fullkey not in h:
-            raise KeyError(fullkey)
+            raise KeyError(f"Missing required key '{fullkey}' in FITS header.")
         return h[fullkey]
 
     # header (d)etector config value assocciated with partial key string.
     def d(key: str):
-        return v("detector " + key)
+        try:
+            return v("detector " + key)
+        except KeyError as e:
+            raise KeyError(f"Detector configuration key error: {e}")
 
     # header (s)etector config value assocciated with partial key string.
     def s(key: str):
-        return v("simulation " + key)
+        try:
+            return v("simulation " + key)
+        except KeyError as e:
+            raise KeyError(f"Simulation configuration key error: {e}")
 
     c = {
         "detector": {

src/nuspacesim/simulation/eas_optical/cphotang.py

@@ -43,10 +43,15 @@
 
 import dask.bag as db
 import numpy as np
-from dask.diagnostics import ProgressBar
+from dask.diagnostics.progress import ProgressBar
 from numpy.polynomial import Polynomial
 
 from .detector_geometry import distance_to_detector
+from .shower_properties import (
+    greisen_particle_count,
+    particle_count_fluctuated_gaisser_hillas,
+    particle_count_parameterized_gaisser_hillas,
+)
 
 # Wrapped in try-catch block as a hack to enable sphinx documentation to be generated
 # on ReadTheDocs without pre-compiling.
@@ -55,13 +60,18 @@
 except ImportError:
     pass
 
+try:
+    from importlib.resources import as_file, files
+except ImportError:
+    from importlib_resources import as_file, files
+
 __all__ = ["CphotAng"]
 
 
 class CphotAng:
     r"""Cherenkov Photon Angle"""
 
-    def __init__(self, detector_altitude):
+    def __init__(self, detector_altitude, longitudinal_profile_func="Greisen"):
         r"""CphotAng: Cherenkov photon density and angle determination class.
 
         Iterative summation of cherenkov radiation reimplemented in numpy and
@@ -244,8 +254,26 @@ def __init__(self, detector_altitude):
         self.zMaxZ = self.dtype(65.0)
         self.RadE = self.dtype(6378.14)
 
-        Zair = self.dtype(7.4)
-        self.ecrit = self.dtype(0.710 / (Zair + 0.96))
+        # Longitudinal Profile Funciton selection
+        if longitudinal_profile_func == "Greisen":
+            self.particle_count = greisen_particle_count
+        elif longitudinal_profile_func == "Gaisser-Hillas Parameterized":
+            self.particle_count = (
+                lambda *args, **kwargs: particle_count_parameterized_gaisser_hillas(
+                    *args, **kwargs
+                )
+            )
+        elif longitudinal_profile_func == "Gaisser-Hillas Fluctuated":
+            with as_file(
+                files("nuspacesim.data.CONEX_table")
+                / "dumpGH_conex_pi_E17_95deg_0km_eposlhc_1394249052_211.dat"
+            ) as file:
+                CONEX_table = np.loadtxt(file, usecols=(4, 5, 6, 7, 8, 9))
+                self.particle_count = (
+                    lambda *args, **kwargs: particle_count_fluctuated_gaisser_hillas(
+                        CONEX_table, *args, **kwargs
+                    )
+                )
 
     def theta_view(self, ThetProp):
         """
@@ -512,19 +540,19 @@ def valid_arrays(self, zsave, delgram, gramsum, gramz, ZonZ, ThetPrpA, Eshow):
         t = np.zeros_like(zsave, dtype=self.dtype)
         t[mask] = gramsum[mask] / self.dtype(36.66)
 
-        greisen_beta = self.dtype(np.log(np.float64(Eshow) / np.float64(self.ecrit)))
+        Zair = self.dtype(7.4)
+        ecrit = self.dtype(0.710 / (Zair + 0.96))
+        greisen_beta = self.dtype(np.log(np.float64(Eshow) / np.float64(ecrit)))
         s = np.zeros_like(zsave, dtype=self.dtype)
         s[mask] = self.dtype(3) * t[mask] / (t[mask] + self.dtype(2) * greisen_beta)
 
-        RN = np.zeros_like(zsave, dtype=self.dtype)
-        RN[mask] = (
-            self.dtype(0.31)
-            / np.sqrt(greisen_beta, dtype=self.dtype)
-            * np.exp(
-                t[mask] * (1 - self.dtype(3 / 2) * np.log(s[mask], dtype=self.dtype)),
-                dtype=self.dtype,
-            )
+        # Greisen  or  Gaisser-Hillas  Longitudinal Paramaterization
+        rn, mask = self.particle_count(
+            mask=mask, t=t, s=s, greisen_beta=greisen_beta, gramsum=gramsum, Eshow=Eshow
         )
+
+        RN = np.zeros_like(zsave, dtype=self.dtype)
+        RN[mask] = rn
         RN[RN < 0] = self.dtype(0)
 
         mask &= ~((RN < 1) & (s > 1))
@@ -616,7 +644,7 @@ def run(self, betaE, alt, Eshow100PeV, lat, long, cloudf=None):
         cloud_top_height = cloudf(lat, long) if cloudf else -np.inf
 
         # early return check, ensure at least 2 segments above cloud-top height.
-        if zs[-2] < cloud_top_height:
+        if len(zs) <= 1 or zs[-2] < cloud_top_height:
             return self.dtype(0), self.dtype(0)
 
         # Cloud mask. No particles will be considered if generated below the clouds.

src/nuspacesim/simulation/eas_optical/eas.py

@@ -53,7 +53,10 @@ class EAS:
 
     def __init__(self, config: NssConfig):
         self.config = config
-        self.CphotAng = CphotAng(self.config.detector.initial_position.altitude)
+        self.CphotAng = CphotAng(
+            self.config.detector.initial_position.altitude,
+            self.config.simulation.cherenkov_light_engine,
+        )
 
     @decorators.nss_result_store("altDec", "lenDec")
     def altDec(self, beta, tauBeta, tauLorentz, u=None, *args, **kwargs):
@@ -88,7 +91,7 @@ def __call__(
         init_long,
         *args,
         cloudf=None,
-        **kwargs
+        **kwargs,
     ):
         """
         Electromagnetic Air Shower operation.

src/nuspacesim/simulation/eas_optical/shower_properties.py

@@ -97,7 +97,7 @@ def shower_age(T):
     return 3.0 * T / (T + 41.77325895999150334743982471)
 
 
-def greisen_particle_count(T, s):
+def greisen_particle_count(t, s, greisen_beta, mask, *args, dtype=np.float32, **kwargs):
     r"""Particle count as a function of radiation length from atmospheric depth
 
     Hillas 1461 eqn (6)
@@ -106,12 +106,93 @@ def greisen_particle_count(T, s):
     (0.31 / sqrt (10^8 / (0.710 / 8.36)))
     = 0.0678308895484773316048795658058110209448440898800928880798622962...
     """
+
     # , param_beta=np.log(10 ** 8 / (0.710 / 8.36))
     # N_e = (0.31 / np.sqrt(param_beta)) * np.exp(T * (1.0 - 1.5 * np.log(s)))
     # N_e[N_e < 0] = 0.0
-    N_e = 0.067830889548477331 * np.exp(T * (1.0 - 1.5 * np.log(s)))
-    N_e[N_e < 0] = 0.0
-    return N_e
+    alpha = dtype(0.31) / np.sqrt(greisen_beta, dtype=dtype)
+
+    RN = alpha * np.exp(
+        t * (dtype(1) - dtype(1.5) * np.log(s, dtype=dtype)), dtype=dtype
+    )
+    RN[RN < 0] = 0.0
+    return RN, mask
+
+
+def gaisser_hillas_particle_count_exp_form(
+    gramsum, X0, Xmax, Nmax, gh_lam, *args, dtype=np.float32, **kwargs
+):
+    # Parametric Form Parameters
+    x = (gramsum - X0) / gh_lam
+    m = (Xmax - X0) / gh_lam
+    return Nmax * np.exp(m * (np.log(x) - np.log(m)) - (x - m))
+
+
+def particle_count_parameterized_gaisser_hillas(
+    gramsum, Eshow, *args, mask=None, dtype=np.float32, **kwargs
+):
+    """
+    Shower particle count from Gaisser Hillas formula with static parameters.
+    """
+
+    # Nuclear Collision length in Air from PDG.
+    # From https://pdg.lbl.gov/2024/AtomicNuclearProperties/HTML/air_dry_1_atm.html
+    X0 = 61.3
+    Xm = 739.0
+    # 65.17 g/cm^2 value obtained by evaluating lambda at Xmax for 1000 upward pion 10^17 eV EAS at 5 deg Earth-emergence angle and starting at sea level
+    gh_lam = 65.12
+    Nmax = 0.045 * (1.0 + 0.0217 * np.log(Eshow / 1.0e5)) * Eshow / 0.074
+    XmaxOff = 58.0 * np.log10(Eshow / 1.0e8)
+    Xmax = Xm + XmaxOff
+
+    particle_count = gaisser_hillas_particle_count_exp_form(
+        gramsum, X0, Xmax, Nmax, gh_lam
+    )
+
+    return particle_count, mask
+
+
+def particle_count_fluctuated_gaisser_hillas(
+    CONEX_table, gramsum, Eshow, mask, *args, dtype=np.float32, **kwargs
+):
+    """
+    Shower particle count from Gaisser Hillas formula with fluctuated parameters.
+    """
+    # Gaisser Hillas Values from CONEX File
+    idx = np.random.randint(low=0, high=CONEX_table.shape[0])
+    Nm, Xm, X0, p1, p2, p3 = CONEX_table[idx]
+    X0 = 0.0 if X0 < 0.0 else X0
+
+    # Masking Gramsum
+    gramsum_mask = gramsum > X0
+    mask &= gramsum_mask
+    Xmask = gramsum[gramsum_mask]
+
+    # JFK : put in form from Tom Gaisser's book, pg: 238 - 239
+    Nmax100 = 6.99e7
+    NmaxE = 0.045 * (1.0 + 0.0217 * np.log(Eshow / 1.0e5)) * Eshow / 0.074
+    Nmax = Nm * NmaxE / Nmax100
+
+    # the following from Gaisser leads to ~80 g/cm^2/decade elongation rate
+    # DOI:10.1103/RevModPhys.83.907, Letessier-Selvon & Stanev
+    #  gives in Egn 3 ~ 85 g/cm^2/decade is for EM showers
+    #    Xmax = 36. * np.log(Eshow/0.074)
+    # HiRes Measurement:  R. U. Abbasi et al 2005 ApJ 622 910 : gives ~ 56 g/cm^2, Auger ~60 g/cm^2
+    # DOI:10.1103/RevModPhys.83.907, Letessier-Selvon & Stanev gives in Egn 7 ~ 62 g/cm^2/decade
+    #   use a single 58 g/cm^2 per decade energy addition/subtraction,
+    #   assume using only 100 PeV energy file for this to be correct
+    XmaxOff = 58.0 * np.log10(Eshow / 1.0e8)
+    Xmax = Xm + XmaxOff
+
+    gh_lam = p1 + p2 * Xmask + p3 * Xmask * Xmask
+    gh_lam[gh_lam > 100.0] = 100.0
+    gh_lam[gh_lam < 1.0e-5] = 1.0e-5
+
+    particle_count = gaisser_hillas_particle_count_exp_form(
+        Xmask, X0, Xmax, Nmax, gh_lam
+    )
+
+    return particle_count, mask
 
 
 def shower_age_of_greisen_particle_count(target_count, x0=2):
@@ -130,13 +211,6 @@ def rns(s):
     return newton(rns, x0)
 
 
-def gaisser_hillas_particle_count(X, Nmax, X0, Xmax, invlam):
-    # return ((X - X0) / (Xmax - X0)) ** xmax * np.exp((Xmax - X) * invlam)
-    xmax = (Xmax - X0) * invlam
-    x = (X - X0) * invlam
-    return Nmax * (x / xmax) ** xmax * np.exp(xmax - x)
-
-
 def slant_depth_trig_approx(z_lo, z_hi, theta_tr, z_max=100.0):
     rho = us_std_atm_density
     r0 = rho(0)

test/core/test_config.py

@@ -290,7 +290,7 @@ def test_default_simulation():
     assert a.max_cherenkov_angle == np.radians(3.0)
     assert a.max_azimuth_angle == np.radians(360.0)
     assert a.angle_from_limb == np.radians(7.0)
-    assert a.cherenkov_light_engine == "Default"
+    assert a.cherenkov_light_engine == "Greisen"
     assert a.ionosphere is not None
     assert a.ionosphere.total_electron_content == 10.0
     assert a.ionosphere.total_electron_error == 0.1
@@ -306,7 +306,7 @@ def test_default_simulation():
         "max_cherenkov_angle": "3.0000000000000004 deg",
         "max_azimuth_angle": "360.0 deg",
         "angle_from_limb": "7.0 deg",
-        "cherenkov_light_engine": "Default",
+        "cherenkov_light_engine": "Greisen",
         "ionosphere": {
             "enable": True,
             "total_electron_content": 10.0,
@@ -437,3 +437,41 @@ def test_config_serialization():
         loaded_config = config_from_toml(tmpfile_name)
 
     assert loaded_config.model_dump() == a.model_dump()
+
+
+def test_cherenkov_light_engine_default_replacement():
+    sim_default = Simulation(cherenkov_light_engine="Default")
+    assert sim_default.cherenkov_light_engine == "Greisen"
+
+
+def test_cherenkov_light_engine_invalid_value():
+    with pytest.raises(ValueError):
+        Simulation(cherenkov_light_engine="InvalidValue")
+
+
+def test_simulation_cherenkov_light_engine_default():
+    sim = Simulation()
+    assert sim.cherenkov_light_engine == "Greisen"
+
+
+def test_simulation_cherenkov_light_engine_default_replacement():
+    sim = Simulation(cherenkov_light_engine="Default")
+    assert sim.cherenkov_light_engine == "Greisen"
+
+
+def test_simulation_cherenkov_light_engine_valid_values():
+    sim_greisen = Simulation(cherenkov_light_engine="Greisen")
+    assert sim_greisen.cherenkov_light_engine == "Greisen"
+
+    sim_gaisser_hillas = Simulation(
+        cherenkov_light_engine="Gaisser-Hillas Parameterized"
+    )
+    assert sim_gaisser_hillas.cherenkov_light_engine == "Gaisser-Hillas Parameterized"
+
+    sim_gaisser_hillas = Simulation(cherenkov_light_engine="Gaisser-Hillas Fluctuated")
+    assert sim_gaisser_hillas.cherenkov_light_engine == "Gaisser-Hillas Fluctuated"
+
+
+def test_simulation_cherenkov_light_engine_invalid_value():
+    with pytest.raises(ValidationError):
+        Simulation(cherenkov_light_engine="InvalidValue")