Dev/HNL_DIS

.github/workflows/build_wheels.yml

@@ -10,7 +10,7 @@ jobs:
       matrix:
         # macos-13 is an intel runner, macos-14 is apple silicon
         os: [macos-13, macos-14]
-        py: [cp38, cp39, cp311, cp312]
+        py: [cp38, cp39, cp310, cp311, cp312]
 
     steps:
       - uses: actions/checkout@v4
@@ -40,9 +40,9 @@ jobs:
     strategy:
       matrix:
         os: [ubuntu-latest]
-        py: [cp38, cp39, cp311, cp312]
+        py: [cp38, cp39, cp310, cp311, cp312]
         image: [manylinux, musllinux]
-        arch: [x86_64, i686, aarch64, ppc64le, s390x]
+        arch: [x86_64, aarch64]
 
     steps:
       - uses: actions/checkout@v4

CMakeLists.txt

@@ -63,6 +63,7 @@ include(pybind11)
 
 # load project dependencies
 include(rk)
+include(crundec3)
 include(cereal)
 include(delabella)
 include(CFITSIO)

cmake/Packages/crundec3.cmake

@@ -0,0 +1 @@
+add_subdirectory(${PROJECT_SOURCE_DIR}/vendor/crundec3 EXCLUDE_FROM_ALL)

projects/dataclasses/private/InteractionRecord.cxx

@@ -207,6 +207,14 @@ void PrimaryDistributionRecord::SetHelicity(double helicity) {
     this->helicity = helicity;
 }
 
+void PrimaryDistributionRecord::SetInteractionParameters(std::map<std::string, double> const & parameters) {
+    interaction_parameters = parameters;
+}
+
+void PrimaryDistributionRecord::SetInteractionParameter(std::string const & name, double value) {
+    interaction_parameters[name] = value;
+}
+
 void PrimaryDistributionRecord::UpdateMass() const {
     if(mass_set)
         return;
@@ -334,6 +342,9 @@ void PrimaryDistributionRecord::Finalize(InteractionRecord & record) const {
     record.primary_mass = GetMass();
     record.primary_momentum = GetFourMomentum();
     record.primary_helicity = GetHelicity();
+    for (auto x : interaction_parameters) {
+        record.interaction_parameters[x.first] = x.second;
+    }
 }
 
 /////////////////////////////////////////
@@ -680,7 +691,10 @@ void CrossSectionDistributionRecord::Finalize(InteractionRecord & record) const
     record.target_mass = target_mass;
     record.target_helicity = target_helicity;
 
-    record.interaction_parameters = interaction_parameters;
+    //record.interaction_parameters = interaction_parameters;
+    for (auto x : interaction_parameters) {
+        record.interaction_parameters[x.first] = x.second;
+    }
 
     record.secondary_ids.resize(secondary_particles.size());
     record.secondary_masses.resize(secondary_particles.size());

projects/dataclasses/private/pybindings/dataclasses.cxx

@@ -124,6 +124,7 @@ PYBIND11_MODULE(dataclasses,m) {
           .def(init<>())
           .def("__str__", [](InteractionRecord const & r) { std::stringstream ss; ss << r; return ss.str(); })
           .def_readwrite("signature",&InteractionRecord::signature)
+          .def_readwrite("primary_initial_position",&InteractionRecord::primary_initial_position)
           .def_readwrite("primary_mass",&InteractionRecord::primary_mass)
           .def_readwrite("primary_momentum",&InteractionRecord::primary_momentum)
           .def_readwrite("primary_helicity",&InteractionRecord::primary_helicity)

projects/dataclasses/public/SIREN/dataclasses/InteractionRecord.h

@@ -67,6 +67,7 @@ class PrimaryDistributionRecord {
     mutable std::array<double, 3> initial_position;
     mutable std::array<double, 3> interaction_vertex;
     mutable double helicity = 0;
+    std::map<std::string, double> interaction_parameters;
 public:
     friend std::ostream& ::operator<<(std::ostream& os, PrimaryDistributionRecord const& record);
 
@@ -104,6 +105,8 @@ class PrimaryDistributionRecord {
     void SetInitialPosition(std::array<double, 3> initial_position);
     void SetInteractionVertex(std::array<double, 3> interaction_vertex);
     void SetHelicity(double helicity);
+    void SetInteractionParameters(std::map<std::string, double> const & parameters);
+    void SetInteractionParameter(std::string const & name, double value);
 
     void UpdateMass() const;
     void UpdateEnergy() const;

projects/dataclasses/public/SIREN/dataclasses/ParticleTypes.def

@@ -1,4 +1,16 @@
 X(unknown, 0)
+X(d,1)
+X(dBar,-1)
+X(u,2)
+X(uBar,-2)
+X(s,3)
+X(sBar,-3)
+X(c,4)
+X(cBar,-4)
+X(b,5)
+X(bBar,-5)
+X(t,6)
+X(tBar,-6)
 X(Gamma, 22)
 X(EPlus, -11)
 X(EMinus, 11)
@@ -13,8 +25,9 @@ X(RhoMinus, -213)
 X(K0_Long, 130)
 X(KPlus, 321)
 X(KMinus, -321)
-X(KStarPlus, 9000311)
-X(KStarMinus, -9000311)
+X(KPrime0, 9000311)
+X(KPrimePlus, 9000321)
+X(KPrimeMinus, -9000321)
 X(Neutron, 2112)
 X(PPlus, 2212)
 X(PMinus, -2212)
@@ -43,7 +56,7 @@ X(DMinus, -411)
 X(D0, 421)
 X(D0Bar, -421)
 X(DsPlus, 431)
-X(DsMinusBar, -431)
+X(DsMinus, -431)
 X(LambdacPlus, 4122)
 X(WPlus, 24)
 X(WMinus, -24)
@@ -167,8 +180,11 @@ X(Nu1, 5911)
 X(Nu2, 5912)
 X(Nu3, 5913)
 X(N4, 5914)
+X(N4Bar, -5914)
 X(N5, 5915)
+X(N5Bar, -5915)
 X(N6, 5916)
+X(N6Bar, -5916)
 X(ZPrime, 5921)
 X(HPrime, 5901)
 X(PhiPrime, 5902)

projects/detector/private/DetectorModel.cxx

@@ -44,6 +44,24 @@ namespace {
 void string_to_lower(std::string & data) {
     std::transform(data.begin(), data.end(), data.begin(), [](unsigned char c){ return std::tolower(c); });
 }
+double estimateDotError(const std::array<double, 3>& p0, const std::array<double, 3>& p1, const std::array<double, 3>& int_dir) {
+    constexpr double epsilon = std::numeric_limits<double>::epsilon();
+
+    double dot_product = 0.0;
+    double error_sum = 0.0;
+
+    for (int i = 0; i < 3; ++i) {
+        double diff = p1[i] - p0[i];
+        double prod = int_dir[i] * diff;
+        dot_product += prod;
+        error_sum += std::fabs(int_dir[i] * diff);
+    }
+
+    // Estimated error of the expression
+    double error_estimate = epsilon * error_sum;
+
+    return error_estimate;
+}
 
 template <class InIt>
 typename std::iterator_traits<InIt>::value_type accumulate(InIt begin, InIt end) {
@@ -244,6 +262,14 @@ DetectorSector DetectorModel::GetSector(int heirarchy) const {
     return sectors_[index];
 }
 
+DetectorSector DetectorModel::GetSector(std::string name) const {
+    for (auto sector : sectors_) {
+        if (sector.name==name) return sector;
+    }
+    std::cout << "Sector " << name << " not found, returning empty sector\n";
+    return DetectorSector();
+}
+
 void DetectorModel::ClearSectors() {
     sectors_.clear();
     sector_map_.clear();
@@ -688,6 +714,11 @@ double DetectorModel::GetInteractionDensity(Geometry::IntersectionList const & i
     } else {
         direction.normalize();
     }
+    // If we have only decays, avoid the sector loop
+    // total_decay_length is in m
+    if(targets.empty()) {
+        return 1.0 / total_decay_length;
+    }
     double dot = direction * intersections.direction;
     assert(std::abs(1.0 - std::abs(dot)) < 1e-6);
     double offset = (intersections.position - p0) * direction;
@@ -698,11 +729,6 @@ double DetectorModel::GetInteractionDensity(Geometry::IntersectionList const & i
         dot = 1;
     }
 
-    // If we have only decays, avoid the sector loop
-    // total_decay_length is in m
-    if(targets.empty()) {
-        return 1.0 / total_decay_length;
-    }
 
     double interaction_density = std::numeric_limits<double>::quiet_NaN();
 
@@ -983,11 +1009,22 @@ double DetectorModel::GetInteractionDepthInCGS(Geometry::IntersectionList const
     }
     Vector3D direction = p1 - p0;
     double distance = direction.magnitude();
+    // this dot error turned out to be much smaller than 1e-6 for event that failed the assertion below
+    // double dot_error = estimateDotError(std::array<double, 3>(p0.get()),
+    //                                     std::array<double, 3>(p1.get()),
+    //                                     std::array<double, 3>(intersections.direction));
+
+    // If we have only decays, avoid the sector loop
+    if(targets.empty()) {
+      return distance/total_decay_length; // m / m --> dimensionless
+    }
     if(distance == 0.0) {
         return 0.0;
     }
     direction.normalize();
 
+    // TODO: a better numerical precision check when the traversed distance is very small
+    // this functionally only happens for decays right now, so we just check for decays at the top
     double dot = intersections.direction * direction;
     assert(std::abs(1.0 - std::abs(dot)) < 1e-6);
     double offset = (intersections.position - p0) * direction;
@@ -998,11 +1035,6 @@ double DetectorModel::GetInteractionDepthInCGS(Geometry::IntersectionList const
         dot = 1;
     }
 
-    // If we have only decays, avoid the sector loop
-    if(targets.empty()) {
-      return distance/total_decay_length; // m / m --> dimensionless
-    }
-
     std::vector<double> interaction_depths(targets.size(), 0.0);
 
     std::function<bool(std::vector<Geometry::Intersection>::const_iterator, std::vector<Geometry::Intersection>::const_iterator, double)> callback =
@@ -1331,6 +1363,11 @@ double DetectorModel::DistanceForInteractionDepthFromPoint(Geometry::Intersectio
         interaction_depth *= -1;
         direction = -direction;
     }
+    // If we have only decays, avoid the sector loop
+    // decay length in m
+    if(targets.empty()) {
+      return interaction_depth * total_decay_length;
+    }
 
     double dot = intersections.direction * direction;
     assert(std::abs(1.0 - std::abs(dot)) < 1e-6);
@@ -1342,11 +1379,6 @@ double DetectorModel::DistanceForInteractionDepthFromPoint(Geometry::Intersectio
         dot = 1;
     }
 
-    // If we have only decays, avoid the sector loop
-    // decay length in m
-    if(targets.empty()) {
-      return interaction_depth * total_decay_length;
-    }
 
     // Recast decay length to cm for density integral
     double total_decay_length_cm = total_decay_length / siren::utilities::Constants::cm;

projects/detector/private/pybindings/DetectorModel.h

@@ -180,7 +180,14 @@ void register_DetectorModel(pybind11::module_ & m) {
         .def_property("DetectorOrigin", &DetectorModel::GetDetectorOrigin, &DetectorModel::SetDetectorOrigin)
         .def_property("DetectorRotation", &DetectorModel::GetDetectorRotation, &DetectorModel::SetDetectorRotation)
         .def("AddSector", &DetectorModel::AddSector)
-        .def("GetSector", &DetectorModel::GetSector)
+        .def("GetSector", (
+                    DetectorSector (DetectorModel::*)(
+                        int) const
+                    )(&DetectorModel::GetSector))
+        .def("GetSector", (
+                    DetectorSector (DetectorModel::*)(
+                        std::string) const
+                    )(&DetectorModel::GetSector))
         .def("ClearSectors", &DetectorModel::ClearSectors)
         .def("GetIntersections", (
                     siren::geometry::Geometry::IntersectionList (DetectorModel::*)(

projects/detector/public/SIREN/detector/DetectorModel.h

@@ -260,6 +260,7 @@ friend siren::detector::Path;
 
     void AddSector(DetectorSector sector);
     DetectorSector GetSector(int level) const;
+    DetectorSector GetSector(std::string name) const;
 
     void ClearSectors();
 

projects/distributions/CMakeLists.txt

@@ -3,6 +3,8 @@
 LIST (APPEND distributions_SOURCES
     ${PROJECT_SOURCE_DIR}/projects/distributions/private/Distributions.cxx
 
+    ${PROJECT_SOURCE_DIR}/projects/distributions/private/primary/PrimaryExternalDistribution.cxx
+
     ${PROJECT_SOURCE_DIR}/projects/distributions/private/primary/direction/PrimaryDirectionDistribution.cxx
     ${PROJECT_SOURCE_DIR}/projects/distributions/private/primary/direction/IsotropicDirection.cxx
     ${PROJECT_SOURCE_DIR}/projects/distributions/private/primary/direction/FixedDirection.cxx
@@ -14,6 +16,9 @@ LIST (APPEND distributions_SOURCES
     ${PROJECT_SOURCE_DIR}/projects/distributions/private/primary/energy/PowerLaw.cxx
     ${PROJECT_SOURCE_DIR}/projects/distributions/private/primary/energy/Monoenergetic.cxx
 
+    ${PROJECT_SOURCE_DIR}/projects/distributions/private/primary/energy_direction/PrimaryEnergyDirectionDistribution.cxx
+    ${PROJECT_SOURCE_DIR}/projects/distributions/private/primary/energy_direction/Tabulated2DFluxDistribution.cxx
+
     ${PROJECT_SOURCE_DIR}/projects/distributions/private/primary/helicity/PrimaryNeutrinoHelicityDistribution.cxx
 
     ${PROJECT_SOURCE_DIR}/projects/distributions/private/primary/mass/PrimaryMass.cxx
@@ -29,6 +34,8 @@ LIST (APPEND distributions_SOURCES
     ${PROJECT_SOURCE_DIR}/projects/distributions/private/primary/vertex/OrientedCylinderPositionDistribution.cxx
     ${PROJECT_SOURCE_DIR}/projects/distributions/private/primary/vertex/CylinderVolumePositionDistribution.cxx
     ${PROJECT_SOURCE_DIR}/projects/distributions/private/primary/vertex/ColumnDepthPositionDistribution.cxx
+    ${PROJECT_SOURCE_DIR}/projects/distributions/private/primary/vertex/PrimaryPhysicalVertexDistribution.cxx
+    ${PROJECT_SOURCE_DIR}/projects/distributions/private/primary/vertex/PrimaryBoundedVertexDistribution.cxx
 
     ${PROJECT_SOURCE_DIR}/projects/distributions/private/secondary/vertex/SecondaryVertexPositionDistribution.cxx
     ${PROJECT_SOURCE_DIR}/projects/distributions/private/secondary/vertex/SecondaryPhysicalVertexDistribution.cxx