Harrison 3029 - PMCL/SWYN - CoDICE Lo direct events: Calculate normalization_per_event and add WIP metadata for it

Author: Menlo Innovations - CAVA Project

imap_l3_processing/cdf/config/imap_codice_l3a_lo-direct-events_variable_attrs.yaml

@@ -162,8 +162,27 @@ normalization:
   LABLAXIS: Normalization
   UNITS: " "
   SCALETYP: linear
+  VALIDMIN: 0
+  VALIDMAX: 1000
   LABL_PTR_2: energy_bin_label
   LABL_PTR_3: spin_angle_bin_label
+  VARIABLE_PURPOSE: ""
+
+normalization_per_event:
+  <<: *per_event
+  NAME: normalization_per_event
+  VAR_TYPE: data
+  CATDESC: Normalization per event
+  FIELDNAM: Normalization factor
+  LABLAXIS: Normalization
+  UNITS: " "
+  SCALETYP: linear
+  VALIDMIN: 0
+  VALIDMAX: 1000
+  SCALEMIN: 0
+  SCALEMAX: 10
+  VARIABLE_PURPOSE: ""
+  DATA_TYPE: CDF_FLOAT
 
 mass_per_charge:
   <<: *per_event

imap_l3_processing/codice/l3/lo/codice_lo_processor.py

@@ -18,7 +18,7 @@
     calculate_mass, calculate_mass_per_charge, \
     rebin_to_counts_by_species_elevation_and_spin_sector, normalize_counts, \
     combine_priorities_and_convert_to_rate, rebin_3d_distribution_azimuth_to_elevation, convert_count_rate_to_intensity, \
-    calculate_normalization_factor
+    calculate_normalization_factor, lookup_normalization_per_event
 from imap_l3_processing.data_utils import safe_divide
 from imap_l3_processing.models import InputMetadata
 from imap_l3_processing.processor import Processor
@@ -192,11 +192,14 @@ def process_l3a_direct_event_data_product(self, dependencies: CodiceLoL3aDirectE
                                                     codice_direct_events.tof)
         mass = calculate_mass(codice_direct_events.apd_energy, codice_direct_events.tof, mass_coefficient_lookup)
         stacked_priorities = np.stack(priority_counts, axis=1)
-        normalization = calculate_normalization_factor(
-            stacked_priorities,
+        normalization = calculate_normalization_factor(stacked_priorities, codice_direct_events.num_events,
+                                                       codice_direct_events.energy_step,
+                                                       codice_direct_events.spin_sector)
+        normalization_per_event = lookup_normalization_per_event(
+            normalization,
             codice_direct_events.num_events,
-            codice_direct_events.spin_sector,
             codice_direct_events.energy_step,
+            codice_direct_events.spin_sector
         )
 
         return CodiceLoL3aDirectEventDataProduct(
@@ -227,7 +230,7 @@ def process_l3a_direct_event_data_product(self, dependencies: CodiceLoL3aDirectE
             rgfo_esa_step=codice_sw_priority_counts_l1a_data.rgfo_esa_step,
             nso_spin_sector=codice_sw_priority_counts_l1a_data.nso_spin_sector,
             nso_esa_step=codice_sw_priority_counts_l1a_data.nso_esa_step,
-            normalization_per_event=1,
+            normalization_per_event=normalization_per_event,
             spin_sector=codice_direct_events.spin_sector,
         )
 

imap_l3_processing/codice/l3/lo/science/codice_lo_calculations.py

@@ -83,7 +83,8 @@ def rebin_direct_events_by_energy_and_spin_sector(num_events: np.ndarray, spin_s
             np.add.at(rebinned_output[time_index, priority_index], (energy_indices, spin_sectors), 1)
     return rebinned_output
 
-def calculate_normalization_factor(priority_counts: np.ndarray, num_events: np.ndarray, spin_sectors: np.ndarray, energy_steps: np.ndarray) -> np.ndarray:
+def calculate_normalization_factor(priority_counts: np.ndarray, num_events: np.ndarray, energy_steps: np.ndarray,
+                                   spin_sectors: np.ndarray) -> np.ndarray:
     numerator = np.concatenate((priority_counts, priority_counts), axis=3)
     num_energies = priority_counts.shape[2]
     num_spin_sectors = 2 * priority_counts.shape[3]
@@ -98,6 +99,16 @@ def calculate_normalization_factor(priority_counts: np.ndarray, num_events: np.n
     output[(denominator==0)&(numerator!=0)] = np.nan
     return output
 
+def lookup_normalization_per_event(normalization: np.ndarray, num_events: np.ndarray, energy_steps: np.ndarray, spin_sectors: np.ndarray) -> np.ndarray:
+    results = np.full(spin_sectors.shape, np.nan)
+    for (epoch, priority), count in np.ma.ndenumerate(num_events, compressed=True):
+        energy_step = energy_steps[epoch, priority, :count]
+        spin_sector = spin_sectors[epoch, priority, :count]
+
+        assert 0 == np.ma.count_masked(energy_step) == np.ma.count_masked(spin_sector) , "Expected all events to have an energy_step and spin_sector!"
+        results[epoch, priority, :count] = normalization[epoch, priority, energy_step, spin_sector]
+    return results
+
 def rebin_to_counts_by_species_elevation_and_spin_sector(num_events: np.ndarray, mass: np.ndarray,
                                                          mass_per_charge: np.ndarray,
                                                          energy: np.ndarray,

tests/codice/l3/lo/science/test_codice_lo_calculations.py

@@ -15,7 +15,7 @@
     rebin_to_counts_by_species_elevation_and_spin_sector, rebin_direct_events_by_energy_and_spin_sector, \
     CODICE_LO_NUM_AZIMUTH_BINS, normalize_counts, combine_priorities_and_convert_to_rate, \
     rebin_3d_distribution_azimuth_to_elevation, convert_count_rate_to_intensity, rebin_direct_events_for_normalization, \
-    calculate_normalization_factor
+    calculate_normalization_factor, lookup_normalization_per_event
 from imap_l3_processing.constants import ONE_SECOND_IN_MICROSECONDS
 
 
@@ -336,7 +336,7 @@ def test_calculate_normalization_factor(self):
         spin_sectors[0, 0, :3] = [0, 2, 12]
         energy_steps[0, 0, :3] = [0, 0, 0]
 
-        result = calculate_normalization_factor(priority_counts, num_events, spin_sectors, energy_steps)
+        result = calculate_normalization_factor(priority_counts, num_events, energy_steps, spin_sectors)
 
         expected_shape = (num_epochs, num_priorities, num_esa_steps, num_l2_spin_sectors)
         expected = np.full(expected_shape, 0.0)
@@ -367,7 +367,7 @@ def test_calculate_normalization_factor_floor_is_one(self):
         spin_sectors[0, 0, :3] = [0, 12, 2]
         energy_steps[0, 0, :3] = [0, 0, 0]
 
-        result = calculate_normalization_factor(priority_counts, num_events, spin_sectors, energy_steps)
+        result = calculate_normalization_factor(priority_counts, num_events, energy_steps, spin_sectors)
 
         expected_shape = (num_epochs, num_priorities, num_esa_steps, num_l2_spin_sectors)
         expected = np.full(expected_shape, 0.0)
@@ -395,7 +395,7 @@ def test_normalization_factor_is_fill_when_priority_nonzero_and_direct_zero(self
         spin_sectors = np.zeros((num_epochs, num_priorities, event_buffer_len), dtype=int)
         energy_steps = np.zeros((num_epochs, num_priorities, event_buffer_len), dtype=int)
 
-        result = calculate_normalization_factor(priority_counts, num_events, spin_sectors, energy_steps)
+        result = calculate_normalization_factor(priority_counts, num_events, energy_steps, spin_sectors)
 
         expected_shape = (num_epochs, num_priorities, num_esa_steps, num_l2_spin_sectors)
         expected = np.full(expected_shape, 0.0)
@@ -407,6 +407,59 @@ def test_normalization_factor_is_fill_when_priority_nonzero_and_direct_zero(self
         expected[0, 0, 0, (0, 2, 12, 14)] = np.nan
         np.testing.assert_equal(result, expected)
 
+    def test_lookup_normalization_per_event(self):
+        num_epochs = 2
+        num_priorities = 7
+        num_esa_steps = 128
+        num_l1_spin_sectors = 12
+        num_l2_spin_sectors = 2*num_l1_spin_sectors
+
+        normalization = np.zeros((num_epochs, num_priorities, num_esa_steps, num_l2_spin_sectors))
+        normalization[0,0,0,0] = 123.
+        normalization[0,0,0,12] = 456.
+        normalization[0,0,1,2] = 789.
+
+        event_buffer_len = 15
+        num_events = np.zeros((num_epochs, num_priorities), dtype=int)
+        spin_sectors = np.ma.array(np.full((num_epochs, num_priorities, event_buffer_len), 255, dtype=int), mask=True)
+        energy_steps = np.ma.array(np.full((num_epochs, num_priorities, event_buffer_len), 255, dtype=int), mask=True)
+
+        num_events[0, 0] = 3
+        spin_sectors[0, 0, :3] = [0, 12, 2]
+        energy_steps[0, 0, :3] = [0, 0, 1]
+
+        normalization_per_event = lookup_normalization_per_event(normalization, num_events, energy_steps, spin_sectors)
+
+        np.testing.assert_array_equal(normalization_per_event[0, 0, :3], [123, 456, 789])
+
+        expected_normalization_per_event = np.full((num_epochs, num_priorities, event_buffer_len), np.nan)
+        expected_normalization_per_event[0, 0, :3] = [123, 456, 789]
+
+        np.testing.assert_array_equal(normalization_per_event, expected_normalization_per_event)
+
+
+    def test_lookup_normalization_per_event_skips_masked_num_events(self):
+        num_epochs = 2
+        num_priorities = 7
+        num_esa_steps = 128
+        num_l1_spin_sectors = 12
+        num_l2_spin_sectors = 2*num_l1_spin_sectors
+
+        normalization = np.zeros((num_epochs, num_priorities, num_esa_steps, num_l2_spin_sectors))
+        normalization[0,0,0,0] = 123.
+        normalization[0,0,0,12] = 456.
+        normalization[0,0,1,2] = 789.
+
+        event_buffer_len = 15
+        num_events = np.ma.array(np.full((num_epochs, num_priorities),65535), mask=True)
+        spin_sectors = np.ma.array(np.full((num_epochs, num_priorities, event_buffer_len), 255, dtype=int), mask=True)
+        energy_steps = np.ma.array(np.full((num_epochs, num_priorities, event_buffer_len), 255, dtype=int), mask=True)
+
+        normalization_per_event = lookup_normalization_per_event(normalization, num_events, energy_steps, spin_sectors)
+
+        expected_normalization_per_event = np.full((num_epochs, num_priorities, event_buffer_len), np.nan)
+        np.testing.assert_array_equal(normalization_per_event, expected_normalization_per_event)
+
     def test_normalize_counts(self):
         num_epochs = 2
         num_priorities = 3

tests/codice/l3/lo/test_codice_lo_processor.py

@@ -499,11 +499,13 @@ def test_process_l3a_partial_densities(self, mock_calculate_partial_densities):
         self.assertEqual(fe_hiq_partial_density, result_data.fe_hiq_partial_density),
 
     @patch('imap_l3_processing.codice.l3.lo.codice_lo_processor.SpinAngleLookup')
-    @patch('imap_l3_processing.codice.l3.lo.codice_lo_processor.calculate_normalization_factor')
-    @patch('imap_l3_processing.codice.l3.lo.codice_lo_processor.calculate_mass_per_charge')
-    @patch('imap_l3_processing.codice.l3.lo.codice_lo_processor.calculate_mass')
+    @patch('imap_l3_processing.codice.l3.lo.codice_lo_processor.lookup_normalization_per_event', autospec=True)
+    @patch('imap_l3_processing.codice.l3.lo.codice_lo_processor.calculate_normalization_factor', autospec=True)
+    @patch('imap_l3_processing.codice.l3.lo.codice_lo_processor.calculate_mass_per_charge', autospec=True)
+    @patch('imap_l3_processing.codice.l3.lo.codice_lo_processor.calculate_mass', autospec=True)
     def test_process_l3a_direct_events(self, mock_calculate_mass, mock_calculate_mass_per_charge,
                                        mock_calculate_normalization_factor,
+                                       mock_lookup_normalization_per_event,
                                        mock_spin_angle_lookup_class):
         rng = np.random.default_rng()
 
@@ -607,9 +609,14 @@ def test_process_l3a_direct_events(self, mock_calculate_mass, mock_calculate_mas
         np.testing.assert_equal(actual_stacked_priorities[:, 6], nsw_priority_rates.p6_hplus_heplusplus)
 
         np.testing.assert_equal(mock_calculate_normalization_factor.call_args.args[1], expected_num_events)
-        np.testing.assert_equal(mock_calculate_normalization_factor.call_args.args[2], expected_spin_sector)
-        np.testing.assert_equal(mock_calculate_normalization_factor.call_args.args[3], expected_energy_step)
+        np.testing.assert_equal(mock_calculate_normalization_factor.call_args.args[2], expected_energy_step)
+        np.testing.assert_equal(mock_calculate_normalization_factor.call_args.args[3], expected_spin_sector)
 
+        self.assertEqual(1, mock_lookup_normalization_per_event.call_count)
+        np.testing.assert_equal(mock_lookup_normalization_per_event.call_args.args[0], expected_normalization)
+        np.testing.assert_equal(mock_lookup_normalization_per_event.call_args.args[1], expected_num_events)
+        np.testing.assert_equal(mock_lookup_normalization_per_event.call_args.args[2], expected_energy_step)
+        np.testing.assert_equal(mock_lookup_normalization_per_event.call_args.args[3], expected_spin_sector)
 
         self.assertIsInstance(l3a_direct_event_data_product, CodiceLoL3aDirectEventDataProduct)
         self.assertEqual(input_metadata, l3a_direct_event_data_product.input_metadata)
@@ -624,6 +631,8 @@ def test_process_l3a_direct_events(self, mock_calculate_mass, mock_calculate_mas
 
         np.testing.assert_array_equal(l3a_direct_event_data_product.normalization,
                                       np.flip(expected_normalization, axis=2))
+        self.assertEqual(l3a_direct_event_data_product.normalization_per_event, mock_lookup_normalization_per_event.return_value)
+
 
         np.testing.assert_array_equal(expected_spin_angle, l3a_direct_event_data_product.spin_angle)
         np.testing.assert_array_equal(expected_spin_sector, l3a_direct_event_data_product.spin_sector)