this is part of the production

plenoirf/analysis/__init__.py

@@ -5,4 +5,3 @@
 from . import energy
 from . import reweight
 from . import pulsar_timing
-from . import particles_on_ground

plenoirf/production/inspect_particle_pool.py

@@ -2,7 +2,6 @@
 import numpy as np
 import sparse_numeric_table as spt
 
-from .. import analysis
 from .. import bookkeeping
 
 
@@ -43,7 +42,7 @@ def run_job(job, logger):
                 aaa["num_air_cherenkov_on_aperture"] = 0
 
             for parblock in parreader:
-                cer = analysis.particles_on_ground.mask_cherenkov_emission(
+                cer = mask_cherenkov_emission(
                     corsika_particles=parblock,
                     corsika_particle_zoo=corsika_particle_zoo,
                 )
@@ -53,10 +52,12 @@ def run_job(job, logger):
 
                 if core:
                     subparblock = parblock[cer["media"]["air"]]
-                    subparblock_r_m = analysis.particles_on_ground.distances_to_point_on_observation_level_m(
-                        corsika_particles=subparblock,
-                        x_m=core["core_x_m"],
-                        y_m=core["core_y_m"],
+                    subparblock_r_m = (
+                        distances_to_point_on_observation_level_m(
+                            corsika_particles=subparblock,
+                            x_m=core["core_x_m"],
+                            y_m=core["core_y_m"],
+                        )
                     )
                     mask_on_aperture = subparblock_r_m <= bin_radius_m
                     aaa["num_air_cherenkov_on_aperture"] += np.sum(
@@ -85,3 +86,56 @@ def record_by_uid(dynamicsizerecarray, uid):
             for name in dynamicsizerecarray._recarray.dtype.names:
                 out[name] = dynamicsizerecarray._recarray[name][i]
             return out
+
+
+def mask_cherenkov_emission(corsika_particles, corsika_particle_zoo):
+    num_particles = corsika_particles.shape[0]
+    media = corsika_particle_zoo.media_cherenkov_threshold_lorentz_factor
+    out = {}
+
+    out["unknown"] = np.zeros(num_particles, dtype=int)
+    out["media"] = {}
+
+    for medium_key in media:
+        out["media"][medium_key] = np.zeros(num_particles, dtype=int)
+
+    for i in range(num_particles):
+        particle = corsika_particles[i]
+
+        corsika_particle_id = cpw.particles.decode_particle_id(
+            code=particle[cpw.I.PARTICLE.CODE]
+        )
+
+        if corsika_particle_zoo.has(corsika_id=corsika_particle_id):
+            momentum_GeV = np.array(
+                [
+                    particle[cpw.I.PARTICLE.PX],
+                    particle[cpw.I.PARTICLE.PY],
+                    particle[cpw.I.PARTICLE.PZ],
+                ]
+            )
+
+            for medium_key in media:
+                if corsika_particle_zoo.cherenkov_emission(
+                    corsika_id=corsika_particle_id,
+                    momentum_GeV=momentum_GeV,
+                    medium_key=medium_key,
+                ):
+                    out["media"][medium_key][i] = 1
+        else:
+            out["unknown"][i] = 1
+
+    return out
+
+
+def distances_to_point_on_observation_level_m(corsika_particles, x_m, y_m):
+    num_particles = corsika_particles.shape[0]
+    distances = np.zeros(num_particles)
+    for i in range(num_particles):
+        par = corsika_particles[i]
+        par_x_m = par[cpw.I.PARTICLE.X] * cpw.CM2M
+        par_y_m = par[cpw.I.PARTICLE.Y] * cpw.CM2M
+        dx = par_x_m - x_m
+        dy = par_y_m - y_m
+        distances[i] = np.hypot(dx, dy)
+    return distances