Fork/alexwenym (#38)

* atlas notebook, added some pybindings * updates to SNe Ranged injection * incremental changes for atlas injection, columndepth debug ongoing * added some pybindings to help debug classes Path and ColumnDepthPositionDistribution * Fix column depth unit conversions. Fix unsigned integer overflow in bounds clipping. Fix spacing. * Special case for constructing 180deg rotation Quaternion for a rotation between opposing directions. * Include what you use * Convert from detector coordinates to earth coordinates before checking bounds * Removed extra files, some may need to be added to resources later * TreeWeighter constructor w/o secondary processes * Delete Doxyfile Removing unnecessary file * small changes to print error messages if needed, M-H sampling in log space for rapidly falling distributions * added inverse CDF sampling to replace the M-H sampling in TablulatedFluxDistribution class. testing is ongoing * Copy photospline's CFITSIO cmake file into our cmake Packages * added units argument to DIS constructors in pybindings * removed comments * remove unncecessary comments * remove comments * removed comments * remove ATLAS directory * remove old legacy analysis * remove old ATLAS geometries * sneaky whitespace * remove whitespace * just kidding, now all the whitespace is gone * add WeightingUtils include --------- Co-authored-by: Austin Schneider <[email protected]> Co-authored-by: Philip Weigel <[email protected]> Co-authored-by: Alex Wen <[email protected]>
Harvard-Neutrino · Jan 10, 2024 · 46841b0 · 46841b0
1 parent c14285f
commit 46841b0
Show file tree

Hide file tree

Showing 52 changed files with 271 additions and 24,568 deletions.
diff --git a/CMakeLists.txt b/CMakeLists.txt
@@ -65,6 +65,7 @@ include(pybind11)
 include(rk)
 include(cereal)
 include(delabella)
+include(CFITSIO)
 include(photospline)
 include(googletest)
 

diff --git a/cmake/Packages/CFITSIO.cmake b/cmake/Packages/CFITSIO.cmake
@@ -0,0 +1,66 @@
+################################################################################
+# Module to find cfitsio                                                       #
+#                                                                              #
+# This module defines:                                                         #
+#                                                                              #
+#   CFITSIO_FOUND                                                              #
+#   CFITSIO_VERSION                                                            #
+#   CFITSIO_LIBRARIES                                                          #
+#   CFITSIO_INCLUDE_DIR                                                        #
+#   CFITSIO_LIB_DIR                                                            #
+#   CFITSIO_CPPFLAGS                                                           #
+#   CFITSIO_LDFLAGS                                                            #
+################################################################################
+
+SET (CFITSIO_FIND_QUIETLY TRUE)
+SET (CFITSIO_FIND_REQUIRED TRUE)
+
+IF (NOT CFITSIO_FOUND)
+
+  # Search user environment for headers, then default paths; extract version
+  FIND_PATH (CFITSIO_INCLUDE_DIR fitsio.h
+    PATHS $ENV{CFITSIOROOT}/include
+    NO_DEFAULT_PATH)
+  FIND_PATH (CFITSIO_INCLUDE_DIR fitsio.h)
+  if(CFITSIO_INCLUDE_DIR)
+    GET_FILENAME_COMPONENT (CFITSIOROOT ${CFITSIO_INCLUDE_DIR} PATH)
+  else()
+    FIND_PATH (CFITSIO_INCLUDE_DIR cfitsio/fitsio.h
+      PATHS $ENV{CFITSIOROOT}/include
+      )
+    SET(CFITSIO_INCLUDE_DIR "${CFITSIO_INCLUDE_DIR}/cfitsio" CACHE PATH "Path to cfitsio headers" FORCE)
+  endif()
+
+  SET (CFITSIO_VERSION 0)
+  IF (CFITSIO_INCLUDE_DIR)
+    FILE (READ "${CFITSIO_INCLUDE_DIR}/fitsio.h" _cfitsio_VERSION)
+    STRING (REGEX REPLACE ".*define CFITSIO_VERSION ([0-9]+\\.[0-9]+).*" "\\1"
+            CFITSIO_VERSION "${_cfitsio_VERSION}")
+  ENDIF (CFITSIO_INCLUDE_DIR)
+
+  # Search user environment for libraries, then default paths
+  FIND_LIBRARY (CFITSIO_LIBRARIES NAMES cfitsio
+    PATHS $ENV{CFITSIOROOT}/lib
+    NO_DEFAULT_PATH)
+  FIND_LIBRARY (CFITSIO_LIBRARIES NAMES cfitsio)
+  GET_FILENAME_COMPONENT (CFITSIO_LIB_DIR ${CFITSIO_LIBRARIES} PATH)
+
+  # Set CFITSIO_FOUND and error out if cfitsio is not found
+  INCLUDE (FindPackageHandleStandardArgs)
+  FIND_PACKAGE_HANDLE_STANDARD_ARGS (CFITSIO
+    DEFAULT_MSG CFITSIO_LIBRARIES CFITSIO_INCLUDE_DIR)
+  ADD_DEFINITIONS ("-I${CFITSIO_INCLUDE_DIR}")
+
+  IF (CFITSIO_FOUND)
+    # Set flags and print a status message
+    MESSAGE (STATUS "CFITSIO version ${CFITSIO_VERSION} found:")
+
+    SET (CFITSIO_CPPFLAGS "-I${CFITSIO_INCLUDE_DIR}")
+    SET (CFITSIO_LDFLAGS "${CFITSIO_LIBRARIES}")
+
+    MESSAGE (STATUS "  * includes: ${CFITSIO_INCLUDE_DIR}")
+    MESSAGE (STATUS "  * libs:     ${CFITSIO_LIBRARIES}")
+  ENDIF (CFITSIO_FOUND)
+
+ENDIF (NOT CFITSIO_FOUND)
+
diff --git a/projects/crosssections/private/DISFromSpline.cxx b/projects/crosssections/private/DISFromSpline.cxx
@@ -54,36 +54,54 @@ bool kinematicallyAllowed(double x, double y, double E, double M, double m) {
 
 DISFromSpline::DISFromSpline() {}
 
-DISFromSpline::DISFromSpline(std::vector<char> differential_data, std::vector<char> total_data, int interaction, double target_mass, double minimum_Q2, std::set<LI::dataclasses::Particle::ParticleType> primary_types, std::set<LI::dataclasses::Particle::ParticleType> target_types) : primary_types_(primary_types), target_types_(target_types), interaction_type_(interaction), target_mass_(target_mass), minimum_Q2_(minimum_Q2) {
+DISFromSpline::DISFromSpline(std::vector<char> differential_data, std::vector<char> total_data, int interaction, double target_mass, double minimum_Q2, std::set<LI::dataclasses::Particle::ParticleType> primary_types, std::set<LI::dataclasses::Particle::ParticleType> target_types, std::string units) : primary_types_(primary_types), target_types_(target_types), interaction_type_(interaction), target_mass_(target_mass), minimum_Q2_(minimum_Q2) {
     LoadFromMemory(differential_data, total_data);
     InitializeSignatures();
+    SetUnits(units);
 }
 
-DISFromSpline::DISFromSpline(std::vector<char> differential_data, std::vector<char> total_data, int interaction, double target_mass, double minimum_Q2, std::vector<LI::dataclasses::Particle::ParticleType> primary_types, std::vector<LI::dataclasses::Particle::ParticleType> target_types) : primary_types_(primary_types.begin(), primary_types.end()), target_types_(target_types.begin(), target_types.end()), interaction_type_(interaction), target_mass_(target_mass), minimum_Q2_(minimum_Q2) {
+DISFromSpline::DISFromSpline(std::vector<char> differential_data, std::vector<char> total_data, int interaction, double target_mass, double minimum_Q2, std::vector<LI::dataclasses::Particle::ParticleType> primary_types, std::vector<LI::dataclasses::Particle::ParticleType> target_types, std::string units) : primary_types_(primary_types.begin(), primary_types.end()), target_types_(target_types.begin(), target_types.end()), interaction_type_(interaction), target_mass_(target_mass), minimum_Q2_(minimum_Q2) {
     LoadFromMemory(differential_data, total_data);
     InitializeSignatures();
+    SetUnits(units);
 }
 
-DISFromSpline::DISFromSpline(std::string differential_filename, std::string total_filename, int interaction, double target_mass, double minimum_Q2, std::set<LI::dataclasses::Particle::ParticleType> primary_types, std::set<LI::dataclasses::Particle::ParticleType> target_types) : primary_types_(primary_types), target_types_(target_types), interaction_type_(interaction), target_mass_(target_mass), minimum_Q2_(minimum_Q2) {
+DISFromSpline::DISFromSpline(std::string differential_filename, std::string total_filename, int interaction, double target_mass, double minimum_Q2, std::set<LI::dataclasses::Particle::ParticleType> primary_types, std::set<LI::dataclasses::Particle::ParticleType> target_types, std::string units) : primary_types_(primary_types), target_types_(target_types), interaction_type_(interaction), target_mass_(target_mass), minimum_Q2_(minimum_Q2) {
     LoadFromFile(differential_filename, total_filename);
     InitializeSignatures();
+    SetUnits(units);
 }
 
-DISFromSpline::DISFromSpline(std::string differential_filename, std::string total_filename, std::set<LI::dataclasses::Particle::ParticleType> primary_types, std::set<LI::dataclasses::Particle::ParticleType> target_types) : primary_types_(primary_types), target_types_(target_types) {
+DISFromSpline::DISFromSpline(std::string differential_filename, std::string total_filename, std::set<LI::dataclasses::Particle::ParticleType> primary_types, std::set<LI::dataclasses::Particle::ParticleType> target_types, std::string units) : primary_types_(primary_types), target_types_(target_types) {
     LoadFromFile(differential_filename, total_filename);
     ReadParamsFromSplineTable();
     InitializeSignatures();
+    SetUnits(units);
 }
 
-DISFromSpline::DISFromSpline(std::string differential_filename, std::string total_filename, int interaction, double target_mass, double minimum_Q2, std::vector<LI::dataclasses::Particle::ParticleType> primary_types, std::vector<LI::dataclasses::Particle::ParticleType> target_types) : primary_types_(primary_types.begin(), primary_types.end()), target_types_(target_types.begin(), target_types.end()), interaction_type_(interaction), target_mass_(target_mass), minimum_Q2_(minimum_Q2) {
+DISFromSpline::DISFromSpline(std::string differential_filename, std::string total_filename, int interaction, double target_mass, double minimum_Q2, std::vector<LI::dataclasses::Particle::ParticleType> primary_types, std::vector<LI::dataclasses::Particle::ParticleType> target_types, std::string units) : primary_types_(primary_types.begin(), primary_types.end()), target_types_(target_types.begin(), target_types.end()), interaction_type_(interaction), target_mass_(target_mass), minimum_Q2_(minimum_Q2) {
     LoadFromFile(differential_filename, total_filename);
     InitializeSignatures();
+    SetUnits(units);
 }
 
-DISFromSpline::DISFromSpline(std::string differential_filename, std::string total_filename, std::vector<LI::dataclasses::Particle::ParticleType> primary_types, std::vector<LI::dataclasses::Particle::ParticleType> target_types) : primary_types_(primary_types.begin(), primary_types.end()), target_types_(target_types.begin(), target_types.end()) {
+DISFromSpline::DISFromSpline(std::string differential_filename, std::string total_filename, std::vector<LI::dataclasses::Particle::ParticleType> primary_types, std::vector<LI::dataclasses::Particle::ParticleType> target_types, std::string units) : primary_types_(primary_types.begin(), primary_types.end()), target_types_(target_types.begin(), target_types.end()) {
     LoadFromFile(differential_filename, total_filename);
     ReadParamsFromSplineTable();
     InitializeSignatures();
+    SetUnits(units);
+}
+
+void DISFromSpline::SetUnits(std::string units) {
+    std::transform(units.begin(), units.end(), units.begin(),
+        [](unsigned char c){ return std::tolower(c); });
+    if(units == "cm") {
+        unit = 1.0;
+    } else if(units == "m") {
+        unit = 10000.0;
+    } else {
+        throw std::runtime_error("Cross section units not supported!");
+    }
 }
 
 bool DISFromSpline::equal(CrossSection const & other) const {
@@ -263,7 +281,7 @@ double DISFromSpline::TotalCrossSection(LI::dataclasses::Particle::ParticleType
     total_cross_section_.searchcenters(&log_energy, &center);
     double log_xs = total_cross_section_.ndsplineeval(&log_energy, &center, 0);
 
-    return std::pow(10.0, log_xs);
+    return unit * std::pow(10.0, log_xs);
 }
 
 // No implementation for DIS yet, just use non-target function
@@ -336,7 +354,7 @@ double DISFromSpline::DifferentialCrossSection(double energy, double x, double y
         return 0;
     double result = pow(10., differential_cross_section_.ndsplineeval(coordinates.data(), centers.data(), 0));
     assert(result >= 0);
-    return result;
+    return unit * result;
 }
 
 double DISFromSpline::InteractionThreshold(dataclasses::InteractionRecord const & interaction) const {

diff --git a/projects/crosssections/private/pybindings/DISFromSpline.h b/projects/crosssections/private/pybindings/DISFromSpline.h
@@ -22,12 +22,12 @@ void register_DISFromSpline(pybind11::module_ & m) {
 
     disfromspline
         .def(init<>())
-        .def(init<std::vector<char>, std::vector<char>, int, double, double, std::set<LI::dataclasses::Particle::ParticleType>, std::set<LI::dataclasses::Particle::ParticleType>>())
-        .def(init<std::vector<char>, std::vector<char>, int, double, double, std::vector<LI::dataclasses::Particle::ParticleType>, std::vector<LI::dataclasses::Particle::ParticleType>>())
-        .def(init<std::string, std::string, int, double, double, std::set<LI::dataclasses::Particle::ParticleType>, std::set<LI::dataclasses::Particle::ParticleType>>())
-        .def(init<std::string, std::string, std::set<LI::dataclasses::Particle::ParticleType>, std::set<LI::dataclasses::Particle::ParticleType>>())
-        .def(init<std::string, std::string, int, double, double, std::vector<LI::dataclasses::Particle::ParticleType>, std::vector<LI::dataclasses::Particle::ParticleType>>())
-        .def(init<std::string, std::string, std::vector<LI::dataclasses::Particle::ParticleType>, std::vector<LI::dataclasses::Particle::ParticleType>>())
+        .def(init<std::vector<char>, std::vector<char>, int, double, double, std::set<LI::dataclasses::Particle::ParticleType>, std::set<LI::dataclasses::Particle::ParticleType>, std::string>())
+        .def(init<std::vector<char>, std::vector<char>, int, double, double, std::vector<LI::dataclasses::Particle::ParticleType>, std::vector<LI::dataclasses::Particle::ParticleType>, std::string>())
+        .def(init<std::string, std::string, int, double, double, std::set<LI::dataclasses::Particle::ParticleType>, std::set<LI::dataclasses::Particle::ParticleType>, std::string>())
+        .def(init<std::string, std::string, std::set<LI::dataclasses::Particle::ParticleType>, std::set<LI::dataclasses::Particle::ParticleType>, std::string>())
+        .def(init<std::string, std::string, int, double, double, std::vector<LI::dataclasses::Particle::ParticleType>, std::vector<LI::dataclasses::Particle::ParticleType>, std::string>())
+        .def(init<std::string, std::string, std::vector<LI::dataclasses::Particle::ParticleType>, std::vector<LI::dataclasses::Particle::ParticleType>, std::string>())
         .def(self == self)
         .def("TotalCrossSection",overload_cast<LI::dataclasses::InteractionRecord const &>(&DISFromSpline::TotalCrossSection, const_))
         .def("TotalCrossSection",overload_cast<LI::dataclasses::Particle::ParticleType, double>(&DISFromSpline::TotalCrossSection, const_))

diff --git a/projects/crosssections/private/pybindings/crosssections.cxx b/projects/crosssections/private/pybindings/crosssections.cxx
@@ -6,7 +6,6 @@
 #include "../../public/LeptonInjector/crosssections/CrossSectionCollection.h"
 #include "../../public/LeptonInjector/crosssections/DISFromSpline.h"
 #include "../../public/LeptonInjector/crosssections/HNLFromSpline.h"
-#include "../../public/LeptonInjector/crosssections/CrossSectionCollection.h"
 #include "../../public/LeptonInjector/crosssections/DipoleFromTable.h"
 #include "../../public/LeptonInjector/crosssections/DarkNewsCrossSection.h"
 #include "../../public/LeptonInjector/crosssections/DarkNewsDecay.h"

diff --git a/projects/crosssections/public/LeptonInjector/crosssections/DISFromSpline.h b/projects/crosssections/public/LeptonInjector/crosssections/DISFromSpline.h
@@ -49,15 +49,19 @@ friend cereal::access;
     int interaction_type_;
     double target_mass_;
     double minimum_Q2_;
+
+    double unit;
 
 public:
     DISFromSpline();
-    DISFromSpline(std::vector<char> differential_data, std::vector<char> total_data, int interaction, double target_mass, double minumum_Q2, std::set<LI::dataclasses::Particle::ParticleType> primary_types, std::set<LI::dataclasses::Particle::ParticleType> target_types);
-    DISFromSpline(std::vector<char> differential_data, std::vector<char> total_data, int interaction, double target_mass, double minumum_Q2, std::vector<LI::dataclasses::Particle::ParticleType> primary_types, std::vector<LI::dataclasses::Particle::ParticleType> target_types);
-    DISFromSpline(std::string differential_filename, std::string total_filename, int interaction, double target_mass, double minumum_Q2, std::set<LI::dataclasses::Particle::ParticleType> primary_types, std::set<LI::dataclasses::Particle::ParticleType> target_types);
-    DISFromSpline(std::string differential_filename, std::string total_filename, std::set<LI::dataclasses::Particle::ParticleType> primary_types, std::set<LI::dataclasses::Particle::ParticleType> target_types);
-    DISFromSpline(std::string differential_filename, std::string total_filename, int interaction, double target_mass, double minumum_Q2, std::vector<LI::dataclasses::Particle::ParticleType> primary_types, std::vector<LI::dataclasses::Particle::ParticleType> target_types);
-    DISFromSpline(std::string differential_filename, std::string total_filename, std::vector<LI::dataclasses::Particle::ParticleType> primary_types, std::vector<LI::dataclasses::Particle::ParticleType> target_types);
+    DISFromSpline(std::vector<char> differential_data, std::vector<char> total_data, int interaction, double target_mass, double minumum_Q2, std::set<LI::dataclasses::Particle::ParticleType> primary_types, std::set<LI::dataclasses::Particle::ParticleType> target_types, std::string units = "cm");
+    DISFromSpline(std::vector<char> differential_data, std::vector<char> total_data, int interaction, double target_mass, double minumum_Q2, std::vector<LI::dataclasses::Particle::ParticleType> primary_types, std::vector<LI::dataclasses::Particle::ParticleType> target_types, std::string units = "cm");
+    DISFromSpline(std::string differential_filename, std::string total_filename, int interaction, double target_mass, double minumum_Q2, std::set<LI::dataclasses::Particle::ParticleType> primary_types, std::set<LI::dataclasses::Particle::ParticleType> target_types, std::string units = "cm");
+    DISFromSpline(std::string differential_filename, std::string total_filename, std::set<LI::dataclasses::Particle::ParticleType> primary_types, std::set<LI::dataclasses::Particle::ParticleType> target_types, std::string units = "cm");
+    DISFromSpline(std::string differential_filename, std::string total_filename, int interaction, double target_mass, double minumum_Q2, std::vector<LI::dataclasses::Particle::ParticleType> primary_types, std::vector<LI::dataclasses::Particle::ParticleType> target_types, std::string units = "cm");
+    DISFromSpline(std::string differential_filename, std::string total_filename, std::vector<LI::dataclasses::Particle::ParticleType> primary_types, std::vector<LI::dataclasses::Particle::ParticleType> target_types, std::string units = "cm");
+
+    void SetUnits(std::string units);
 
     virtual bool equal(CrossSection const & other) const override;
 

diff --git a/projects/detector/private/EarthModel.cxx b/projects/detector/private/EarthModel.cxx
@@ -697,6 +697,7 @@ double EarthModel::DistanceForColumnDepthFromPoint(Geometry::IntersectionList co
             EarthSector sector = GetSector(current_intersection->hierarchy);
             double target = column_depth - total_column_depth;
             double distance = sector.density->InverseIntegral(p0+start_point*direction, direction, target, segment_length);
+
             done = distance >= 0;
             double integral = sector.density->Integral(p0+start_point*direction, direction, segment_length);
             total_column_depth += integral;
@@ -873,6 +874,7 @@ double EarthModel::GetInteractionDepthInCGS(Geometry::IntersectionList const & i
             double segment_length = end_point - start_point;
             EarthSector sector = GetSector(current_intersection->hierarchy);
             double integral = sector.density->Integral(p0+start_point*direction, direction, segment_length);
+
             std::vector<double> particle_fractions = materials_.GetTargetParticleFraction(sector.material_id, targets.begin(), targets.end());
             for(unsigned int i=0; i<targets.size(); ++i) {
                 interaction_depths[i] += (integral * 100) * particle_fractions[i]; // cm^-3 * m --> cm^-2