Fix controller and weighter bugs from PR #125 review

Weighter (C++):
- Implement SurvivalProbability in Weighter.tcc with [0,1] clamping
- Add i_inj bounds-checking in GetInteractionProbabilities/GetSurvivalProbabilities
- Use SurvivalProbability method instead of inline 1-InteractionProbability

SIREN_Controller (Python):
- Fix SetInteractions assertion to check physical process type, not injection
- Fix save_int_params to produce per-event nested dicts instead of flat lists
- Fix InputDarkNewsDecay to only compute min decay width for matching decays
- Fix InputDarkNewsModel to reuse existing secondary processes instead of discarding
- Replace exit(0) calls with ValueError/TypeError/RuntimeError
- Fix get_material_model_path (non-existent function) with correct path construction
- Fix typos: "interseted", "sedcondary"

Tests:
- Add Weighter_TEST.cxx for one_minus_exp_of_negative helper accuracy
- Add test_controller_fixes.py covering all Python and Weighter fixes

Author: austinschneider

projects/injection/CMakeLists.txt

@@ -38,6 +38,7 @@ install(DIRECTORY "${PROJECT_SOURCE_DIR}/projects/injection/public/"
 )
 
 #package_add_test(UnitTest_Injector ${PROJECT_SOURCE_DIR}/projects/injection/private/test/Injector_TEST.cxx)
+package_add_test(UnitTest_Weighter ${PROJECT_SOURCE_DIR}/projects/injection/private/test/Weighter_TEST.cxx)
 if(NOT ${CIBUILDWHEEL})
 package_add_test(UnitTest_CCM_HNL ${PROJECT_SOURCE_DIR}/projects/injection/private/test/CCM_HNL_TEST.cxx)
 target_link_libraries(UnitTest_CCM_HNL pybind11::embed)

projects/injection/private/Weighter.cxx

@@ -136,10 +136,9 @@ double Weighter::EventWeight(siren::dataclasses::InteractionTree const & tree) c
 }
 
 std::vector<double> Weighter::GetInteractionProbabilities(siren::dataclasses::InteractionTree const & tree, int i_inj) const {
-    // Returns the vector of interaction physical probabilities for each process
-    // Since we are concerned only with the physical probability, we use the first injector since physical processes are the same for all injectors
-    // HOWEVER the injection bounds will change based on the injector
-    // so, we allow the user to specify which injector they are interseted in
+    if(i_inj < 0 || static_cast<size_t>(i_inj) >= injectors.size()) {
+        throw std::out_of_range("i_inj index out of range in GetInteractionProbabilities");
+    }
 
     std::vector<double> int_probs;
     for(auto const & datum : tree.tree) {
@@ -162,22 +161,23 @@ std::vector<double> Weighter::GetInteractionProbabilities(siren::dataclasses::In
 }
 
 std::vector<double> Weighter::GetSurvivalProbabilities(siren::dataclasses::InteractionTree const & tree, int i_inj) const {
-    // This allows the user to get the survival probabilities for each interaction
-    // Useful in the case that secondary interactions are restricted to fiducial volumes
+    if(i_inj < 0 || static_cast<size_t>(i_inj) >= injectors.size()) {
+        throw std::out_of_range("i_inj index out of range in GetSurvivalProbabilities");
+    }
 
     std::vector<double> survival_probs;
     for(auto const & datum : tree.tree) {
         std::tuple<siren::math::Vector3D, siren::math::Vector3D> bounds;
         if(datum->depth() == 0) {
-            std::get<0>(bounds) = datum->record.primary_initial_position; // start location
-            std::get<1>(bounds) = std::get<0>(injectors[i_inj]->PrimaryInjectionBounds(datum->record)); // start of injection bounds
-            survival_probs.push_back(1 - primary_process_weighters[i_inj]->InteractionProbability(bounds, datum->record));
+            std::get<0>(bounds) = datum->record.primary_initial_position;
+            std::get<1>(bounds) = std::get<0>(injectors[i_inj]->PrimaryInjectionBounds(datum->record));
+            survival_probs.push_back(primary_process_weighters[i_inj]->SurvivalProbability(bounds, datum->record));
         }
         else {
             try {
-                std::get<0>(bounds) = datum->record.primary_initial_position; // start location
-                std::get<1>(bounds) = std::get<0>(injectors[i_inj]->SecondaryInjectionBounds(datum->record)); // start of injection bounds
-                survival_probs.push_back(1 - secondary_process_weighter_maps[i_inj].at(datum->record.signature.primary_type)->InteractionProbability(bounds, datum->record));
+                std::get<0>(bounds) = datum->record.primary_initial_position;
+                std::get<1>(bounds) = std::get<0>(injectors[i_inj]->SecondaryInjectionBounds(datum->record));
+                survival_probs.push_back(secondary_process_weighter_maps[i_inj].at(datum->record.signature.primary_type)->SurvivalProbability(bounds, datum->record));
             } catch(const std::out_of_range& oor) {
                 std::cout << "Out of Range error: " << oor.what() << '\n';
                 return {};

projects/injection/private/test/Weighter_TEST.cxx

@@ -0,0 +1,111 @@
+#include <cmath>
+#include <stdexcept>
+
+#include <gtest/gtest.h>
+
+// Forward-declare the helper functions from Weighter.tcc (linked via SIREN lib)
+// to avoid ODR violations from including the full header.
+namespace siren { namespace injection {
+    double one_minus_exp_of_negative(double x);
+    double log_one_minus_exp_of_negative(double x);
+}}
+
+using namespace siren::injection;
+
+// ---------------------------------------------------------------------------
+// one_minus_exp_of_negative
+// ---------------------------------------------------------------------------
+
+TEST(WeighterHelpers, OneMinusExpSmallX) {
+    double x = 1e-5;
+    double result = one_minus_exp_of_negative(x);
+    double exact = 1.0 - std::exp(-x);
+    EXPECT_NEAR(result, exact, 1e-15);
+}
+
+TEST(WeighterHelpers, OneMinusExpMediumX) {
+    double x = 0.5;
+    double result = one_minus_exp_of_negative(x);
+    double exact = 1.0 - std::exp(-x);
+    EXPECT_NEAR(result, exact, 1e-12);
+}
+
+TEST(WeighterHelpers, OneMinusExpLargeX) {
+    double x = 5.0;
+    double result = one_minus_exp_of_negative(x);
+    double exact = 1.0 - std::exp(-x);
+    EXPECT_NEAR(result, exact, 1e-12);
+}
+
+TEST(WeighterHelpers, OneMinusExpVerySmallX) {
+    // Exercises the Taylor expansion branch
+    double x = 1e-8;
+    double result = one_minus_exp_of_negative(x);
+    double exact = 1.0 - std::exp(-x);
+    EXPECT_NEAR(result, exact, 1e-15);
+}
+
+TEST(WeighterHelpers, OneMinusExpAtBranchPoint) {
+    // Near the 0.1 branch boundary
+    double x = 0.099;
+    double result = one_minus_exp_of_negative(x);
+    double exact = 1.0 - std::exp(-x);
+    EXPECT_NEAR(result, exact, 1e-12);
+
+    x = 0.101;
+    result = one_minus_exp_of_negative(x);
+    exact = 1.0 - std::exp(-x);
+    EXPECT_NEAR(result, exact, 1e-12);
+}
+
+TEST(WeighterHelpers, OneMinusExpResultBounded) {
+    // Result should always be in [0, 1) for non-negative x
+    for(double x = 0.0; x <= 20.0; x += 0.1) {
+        double result = one_minus_exp_of_negative(x);
+        EXPECT_GE(result, 0.0) << "Failed at x=" << x;
+        EXPECT_LT(result, 1.0) << "Failed at x=" << x;
+    }
+}
+
+// ---------------------------------------------------------------------------
+// log_one_minus_exp_of_negative
+// ---------------------------------------------------------------------------
+
+TEST(WeighterHelpers, LogOneMinusExpSmallX) {
+    double x = 1e-5;
+    double result = log_one_minus_exp_of_negative(x);
+    double exact = std::log(1.0 - std::exp(-x));
+    EXPECT_NEAR(result, exact, 1e-10);
+}
+
+TEST(WeighterHelpers, LogOneMinusExpMidX) {
+    double x = 1.5;
+    double result = log_one_minus_exp_of_negative(x);
+    double exact = std::log(1.0 - std::exp(-x));
+    EXPECT_NEAR(result, exact, 1e-12);
+}
+
+TEST(WeighterHelpers, LogOneMinusExpLargeX) {
+    // Exercises the exp-series branch (x > 3)
+    double x = 5.0;
+    double result = log_one_minus_exp_of_negative(x);
+    double exact = std::log(1.0 - std::exp(-x));
+    EXPECT_NEAR(result, exact, 1e-12);
+}
+
+TEST(WeighterHelpers, LogOneMinusExpAtBranchPoints) {
+    // Test near the 0.1 and 3.0 branch boundaries
+    for(double x : {0.09, 0.11, 2.99, 3.01}) {
+        double result = log_one_minus_exp_of_negative(x);
+        double exact = std::log(1.0 - std::exp(-x));
+        EXPECT_NEAR(result, exact, 1e-9) << "Failed at x=" << x;
+    }
+}
+
+TEST(WeighterHelpers, LogOneMinusExpIsNegative) {
+    // log(1 - exp(-x)) is always negative for x > 0
+    for(double x = 0.001; x <= 20.0; x += 0.1) {
+        double result = log_one_minus_exp_of_negative(x);
+        EXPECT_LT(result, 0.0) << "Failed at x=" << x;
+    }
+}

projects/injection/public/SIREN/injection/Weighter.tcc

@@ -142,6 +142,15 @@ double ProcessWeighter<ProcessType>::InteractionProbability(std::tuple<siren::ma
     return interaction_probability;
 }
 
+template<typename ProcessType>
+double ProcessWeighter<ProcessType>::SurvivalProbability(std::tuple<siren::math::Vector3D, siren::math::Vector3D> const & bounds, siren::dataclasses::InteractionRecord const & record) const {
+    double interaction_probability = InteractionProbability(bounds, record);
+    double survival = 1.0 - interaction_probability;
+    if(survival < 0.0) survival = 0.0;
+    if(survival > 1.0) survival = 1.0;
+    return survival;
+}
+
 template<typename ProcessType>
 double ProcessWeighter<ProcessType>::NormalizedPositionProbability(std::tuple<siren::math::Vector3D, siren::math::Vector3D> const & bounds, siren::dataclasses::InteractionRecord const & record) const {
     using siren::detector::DetectorPosition;

python/SIREN_Controller.py

@@ -69,11 +69,11 @@ def __init__(self, events_to_inject, experiment=None, detector_model_file=None,
         self.materials_model_file = materials_model_file
         if experiment is not None:
             # Find the density and materials files
-            self.materials_model_file = _util.get_material_model_path(experiment)
-            self.detector_model_file = _util.get_detector_model_path(experiment)
+            detector_dir = _util.get_detector_model_path(experiment)
+            self.materials_model_file = os.path.join(detector_dir, "materials.dat")
+            self.detector_model_file = os.path.join(detector_dir, "densities.dat")
         elif (self.detector_model_file is None or self.materials_model_file is None):
-            print("Must provide either an experiment name or both a detector model file and materials model file. Exiting")
-            exit(0)
+            raise ValueError("Must provide either an experiment name or both a detector model file and materials model file")
 
         self.detector_model = _detector.DetectorModel()
         self.detector_model.LoadMaterialModel(self.materials_model_file)
@@ -277,20 +277,23 @@ def InputDarkNewsModel(self, primary_type, table_dir, upscattering=True, decay=T
         secondary_interaction_collections = []
         for secondary_type, decay_list in secondary_decays.items():
 
-            # Define a sedcondary injection distribution if necessary
+            # Define a secondary injection distribution if necessary
             inj_sec_defined = False
             phys_sec_defined = False
+            existing_inj_process = None
             for secondary_injection_process in self.secondary_injection_processes:
                 if secondary_injection_process.primary_type == secondary_type:
                     inj_sec_defined = True
+                    existing_inj_process = secondary_injection_process
             for secondary_physical_process in self.secondary_physical_processes:
                 if secondary_physical_process.primary_type == secondary_type:
                     phys_sec_defined = True
 
-            secondary_injection_process = _injection.SecondaryInjectionProcess()
-            secondary_physical_process = _injection.PhysicalProcess()
-            secondary_injection_process.primary_type = secondary_type
-            secondary_physical_process.primary_type = secondary_type
+            if inj_sec_defined:
+                secondary_injection_process = existing_inj_process
+            else:
+                secondary_injection_process = _injection.SecondaryInjectionProcess()
+                secondary_injection_process.primary_type = secondary_type
 
             # Add the secondary position distribution
             if fid_vol_secondary and self.fid_vol is not None:
@@ -302,8 +305,12 @@ def InputDarkNewsModel(self, primary_type, table_dir, upscattering=True, decay=T
                     _distributions.SecondaryPhysicalVertexDistribution()
                 )
 
-            if not inj_sec_defined: self.secondary_injection_processes.append(secondary_injection_process)
-            if not phys_sec_defined: self.secondary_physical_processes.append(secondary_physical_process)
+            if not inj_sec_defined:
+                self.secondary_injection_processes.append(secondary_injection_process)
+            if not phys_sec_defined:
+                secondary_physical_process = _injection.PhysicalProcess()
+                secondary_physical_process.primary_type = secondary_type
+                self.secondary_physical_processes.append(secondary_physical_process)
 
             secondary_interaction_collections.append(
                 _interactions.InteractionCollection(secondary_type, decay_list)
@@ -343,9 +350,9 @@ def InputDarkNewsDecay(self, primary_type, table_dir, **kwargs):
                 decay.dec_case.nu_parent.pdgid
             ):
                 primary_decays.append(decay)
-            total_decay_width = decay.TotalDecayWidth(primary_type)
-            if total_decay_width < self.DN_min_decay_width:
-                self.DN_min_decay_width = total_decay_width
+                total_decay_width = decay.TotalDecayWidth(primary_type)
+                if total_decay_width > 0 and total_decay_width < self.DN_min_decay_width:
+                    self.DN_min_decay_width = total_decay_width
         primary_interaction_collection = _interactions.InteractionCollection(
             primary_type, primary_decays
         )
@@ -377,8 +384,7 @@ def GetFiducialVolume(self):
     def GetVolumePositionDistributionFromSector(self, sector_name):
         geo = self.GetDetectorSectorGeometry(sector_name)
         if geo is None:
-            print("Sector %s not found. Exiting"%sector_name)
-            exit(0)
+            raise ValueError("Sector %s not found" % sector_name)
         # the position is in geometry coordinates
         # must update to detector coordintes
         det_position = self.detector_model.GeoPositionToDetPosition(_detector.GeometryPosition(geo.placement.Position))
@@ -391,8 +397,7 @@ def GetVolumePositionDistributionFromSector(self, sector_name):
             sphere = _geometry.Sphere(det_placement,geo.Radius,geo.InnerRadius)
             return _distributions.SphereVolumePositionDistribution(sphere)
         else:
-            print("Geometry type %s not supported for position distribution. Exiting"%str(type(geo)))
-            exit(0)
+            raise TypeError("Geometry type %s not supported for position distribution" % str(type(geo)))
 
     def GetDetectorModelTargets(self):
         """
@@ -449,7 +454,7 @@ def SetInteractions(
                 if self.primary_physical_process.primary_type == _dataclasses.Particle.ParticleType.unknown:
                     self.primary_physical_process.primary_type = primary_interaction_collection.GetPrimaryType()
                 else:
-                    assert(self.primary_injection_process.primary_type == primary_interaction_collection.GetPrimaryType())
+                    assert(self.primary_physical_process.primary_type == primary_interaction_collection.GetPrimaryType())
                 if self.primary_physical_process.interactions is None:
                     self.primary_physical_process.interactions = primary_interaction_collection
                 else:
@@ -484,11 +489,10 @@ def SetInteractions(
                                                                                     [sec_phys.interactions, sec_ints])
                         found_collection = True
                 if not found_collection and(sec_inj.interactions is None or sec_phys.interactions is None):
-                    print(
-                        "Couldn't find cross section collection for secondary particle %s; Exiting"
+                    raise RuntimeError(
+                        "Couldn't find cross section collection for secondary particle %s"
                         % record.signature.primary_type
                     )
-                    exit(0)
 
     # set the stopping condition of the injector with a python function
     # must accept two arguments, assumes first is datum and the second is the index of the secondary particle
@@ -636,12 +640,16 @@ def SaveEvents(self, filename, fill_tables_at_exit=True,
                       "parent_idx",
                       "num_daughters"]:
                 datasets[k].append([])
+            if save_int_params:
+                datasets.setdefault("int_params", [])
+                datasets["int_params"].append({})
             # loop over interactions
             for id, datum in enumerate(event.tree):
                 if save_int_params:
-                    for param_name,param_value in datum.record.interaction_parameters.items():
-                        if ie==0: datasets[param_name] = []
-                        datasets[param_name].append(param_value)
+                    for param_name, param_value in datum.record.interaction_parameters.items():
+                        if param_name not in datasets["int_params"][-1]:
+                            datasets["int_params"][-1][param_name] = []
+                        datasets["int_params"][-1][param_name].append(param_value)
                 datasets["vertex"][-1].append(np.array(datum.record.interaction_vertex,dtype=float))
                 datasets["primary_initial_position"][-1].append(np.array(datum.record.primary_initial_position,dtype=float))