From e5cdd8bc989024aaa39cf0b73dc0a06b8b1f064f Mon Sep 17 00:00:00 2001 From: Austin Schneider Date: Mon, 15 Jan 2024 17:50:18 -0700 Subject: [PATCH 01/42] Add NamedType as dependency --- .gitmodules | 3 +++ vendor/NamedType | 1 + 2 files changed, 4 insertions(+) create mode 160000 vendor/NamedType diff --git a/.gitmodules b/.gitmodules index 50c8c3f5c..fd9eb6b5f 100644 --- a/.gitmodules +++ b/.gitmodules @@ -14,3 +14,6 @@ [submodule "vendor/photospline"] path = vendor/photospline url = https://github.com/austinschneider/photospline.git +[submodule "vendor/NamedType"] + path = vendor/NamedType + url = https://github.com/joboccara/NamedType.git diff --git a/vendor/NamedType b/vendor/NamedType new file mode 160000 index 000000000..76668abe0 --- /dev/null +++ b/vendor/NamedType @@ -0,0 +1 @@ +Subproject commit 76668abe09807f92a695ee5e868f9719e888e65f From 3a183cc781c4cbf31f140dbf90d8f8fb8467fd3b Mon Sep 17 00:00:00 2001 From: Austin Schneider Date: Mon, 15 Jan 2024 17:50:41 -0700 Subject: [PATCH 02/42] Add NamedType to cmake --- CMakeLists.txt | 1 + cmake/Packages/NamedType.cmake | 23 +++++++++++++++++++++++ 2 files changed, 24 insertions(+) create mode 100644 cmake/Packages/NamedType.cmake diff --git a/CMakeLists.txt b/CMakeLists.txt index 1a085ab82..0dfafd561 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -68,6 +68,7 @@ include(delabella) include(CFITSIO) include(photospline) include(googletest) +include(NamedType) # load macros for googletest include(testing) diff --git a/cmake/Packages/NamedType.cmake b/cmake/Packages/NamedType.cmake new file mode 100644 index 000000000..14d5e4c71 --- /dev/null +++ b/cmake/Packages/NamedType.cmake @@ -0,0 +1,23 @@ +find_package(Git QUIET) +if(GIT_FOUND AND EXISTS "${PROJECT_SOURCE_DIR}/.git") +# Update submodules as needed + option(GIT_SUBMODULE "Check submodules during build" ON) + if(GIT_SUBMODULE) + if(NOT EXISTS "${PROJECT_SOURCE_DIR}/vendor/NamedType/CMakeLists.txt") + message(STATUS "Submodule update") + execute_process(COMMAND ${GIT_EXECUTABLE} submodule update --init ./ + WORKING_DIRECTORY ${CMAKE_CURRENT_SOURCE_DIR}/vendor/NamedType + RESULT_VARIABLE GIT_SUBMOD_RESULT) + if(NOT GIT_SUBMOD_RESULT EQUAL "0") + message(FATAL_ERROR "git submodule update --init failed with ${GIT_SUBMOD_RESULT}, please checkout submodules") + endif() + endif() + endif() +endif() + +if(NOT EXISTS "${PROJECT_SOURCE_DIR}/vendor/NamedType/CMakeLists.txt") + message(FATAL_ERROR "The NamedType submodule was not downloaded! GIT_SUBMODULE was turned off or failed. Please update submodules and try again.") +endif() + +add_subdirectory("${PROJECT_SOURCE_DIR}/vendor/NamedType" "extern/NamedType") +include_directories("${PROJECT_SOURCE_DIR}/vendor/NamedType/include") From 001a388eb63d45a3046d87394b24c719be95ae54 Mon Sep 17 00:00:00 2001 From: Austin Schneider Date: Mon, 15 Jan 2024 17:51:14 -0700 Subject: [PATCH 03/42] Define special named types for detector and geometry --- .../LeptonInjector/detector/Coordinates.h | 20 +++++++++++++++++++ 1 file changed, 20 insertions(+) create mode 100644 projects/detector/public/LeptonInjector/detector/Coordinates.h diff --git a/projects/detector/public/LeptonInjector/detector/Coordinates.h b/projects/detector/public/LeptonInjector/detector/Coordinates.h new file mode 100644 index 000000000..81d212245 --- /dev/null +++ b/projects/detector/public/LeptonInjector/detector/Coordinates.h @@ -0,0 +1,20 @@ +#pragma once +#ifndef LI_Coordinates_H +#define LI_Coordinates_H + +#include "LeptonInjector/math/Vector3D.h" // for Vector3D +#include + + +namespace LI { +namespace detector { + +using GeometryPosition = fluent::NamedType; +using GeometryDirection = fluent::NamedType; +using DetectorPosition = fluent::NamedType; +using DetectorDirection = fluent::NamedType; + +} // namespace detector +} // namespace LI + +#endif // LI_Coordinates_H From c1d339c1cff9813293024da40bd0c097343c5605 Mon Sep 17 00:00:00 2001 From: Austin Schneider Date: Mon, 15 Jan 2024 17:51:37 -0700 Subject: [PATCH 04/42] Sketch the new interface for DetectorModel --- projects/detector/private/DetectorModel.cxx | 240 +++++++++++++++--- .../LeptonInjector/detector/DetectorModel.h | 162 +++++++++--- .../LeptonInjector/detector/DetectorModel.tcc | 32 ++- 3 files changed, 350 insertions(+), 84 deletions(-) diff --git a/projects/detector/private/DetectorModel.cxx b/projects/detector/private/DetectorModel.cxx index 045d0aad1..3849cbb4d 100644 --- a/projects/detector/private/DetectorModel.cxx +++ b/projects/detector/private/DetectorModel.cxx @@ -118,6 +118,38 @@ std::ostream & DetectorSector::Print(std::ostream& oss) const { return oss; } +DetectorPosition DetectorModel::ToDet(GeometryPosition const & pos) const { + return DetectorPosition(detector_rotation_.rotate(pos - detector_origin_, false)); +} + +DetectorDirection DetectorModel::ToDet(GeometryDirection const & dir) const { + return DetectorDirection(detector_rotation_.rotate(dir, false)); +} + +DetectorPosition DetectorModel::ToDet(GeometryPosition && pos) const { + return DetectorPosition(detector_rotation_.rotate(pos - detector_origin_, false)); +} + +DetectorDirection DetectorModel::ToDet(GeometryDirection && dir) const { + return DetectorDirection(detector_rotation_.rotate(dir, false)); +} + +GeometryPosition DetectorModel::ToGeo(DetectorPosition const & pos) const { + return GeometryPosition(detector_rotation_..rotate(pos, true) + detector_origin_); +} + +GeometryDirection DetectorModel::ToGeo(DetectorDirection const & dir) const { + return GeometryDirection(detector_rotation_.rotate(dir, true)); +} + +GeometryPosition DetectorModel::ToGeo(DetectorPosition && pos) const { + return GeometryPosition(detector_rotation_.rotate(pos, true) + detector_origin_); +} + +GeometryDirection DetectorModel::ToGeo(DetectorDirection && dir) const { + return GeometryDirection(detector_rotation_.rotate(dir, true)); +} + std::ostream& operator<<(std::ostream& oss, DetectorSector const & bcm) { return(bcm.Print(oss)); } @@ -176,11 +208,11 @@ void DetectorModel::SetSectors(std::vector const & sectors) { sectors_ = sectors; } -Vector3D DetectorModel::GetDetectorOrigin() const { +GeometryPosition DetectorModel::GetDetectorOrigin() const { return detector_origin_; } -void DetectorModel::SetDetectorOrigin(Vector3D const & detector_origin) { +void DetectorModel::SetDetectorOrigin(GeometryPosition const & detector_origin) { detector_origin_ = detector_origin; } @@ -444,7 +476,7 @@ void DetectorModel::LoadMaterialModel(std::string const & material_model) { } -double DetectorModel::GetMassDensity(Geometry::IntersectionList const & intersections, Vector3D const & p0) const { +double DetectorModel::GetMassDensity(Geometry::IntersectionList const & intersections, GeometryPosition const & p0) const { Vector3D direction = p0 - intersections.position; if(direction.magnitude() == 0) { direction = intersections.direction; @@ -484,13 +516,13 @@ double DetectorModel::GetMassDensity(Geometry::IntersectionList const & intersec return density; } -double DetectorModel::GetMassDensity(Vector3D const & p0) const { +double DetectorModel::GetMassDensity(GeometryPosition const & p0) const { Vector3D direction(1,0,0); // Any direction will work for determining the sector heirarchy Geometry::IntersectionList intersections = GetIntersections(p0, direction); return GetMassDensity(intersections, p0); } -double DetectorModel::GetParticleDensity(Geometry::IntersectionList const & intersections, Vector3D const & p0, LI::dataclasses::Particle::ParticleType target) const { +double DetectorModel::GetParticleDensity(Geometry::IntersectionList const & intersections, GeometryPosition const & p0, LI::dataclasses::Particle::ParticleType target) const { Vector3D direction = p0 - intersections.position; if(direction.magnitude() == 0) { direction = intersections.direction; @@ -531,13 +563,13 @@ double DetectorModel::GetParticleDensity(Geometry::IntersectionList const & inte return density; } -double DetectorModel::GetParticleDensity(Vector3D const & p0, LI::dataclasses::Particle::ParticleType target) const { +double DetectorModel::GetParticleDensity(GeometryPosition const & p0, LI::dataclasses::Particle::ParticleType target) const { Vector3D direction(1,0,0); // Any direction will work for determining the sector heirarchy Geometry::IntersectionList intersections = GetIntersections(p0, direction); return GetParticleDensity(intersections, p0, target); } -double DetectorModel::GetInteractionDensity(Geometry::IntersectionList const & intersections, Vector3D const & p0, +double DetectorModel::GetInteractionDensity(Geometry::IntersectionList const & intersections, GeometryPosition const & p0, std::vector const & targets, std::vector const & total_cross_sections, double const & total_decay_length) const { @@ -595,7 +627,7 @@ double DetectorModel::GetInteractionDensity(Geometry::IntersectionList const & i return interaction_density; } -double DetectorModel::GetInteractionDensity(Vector3D const & p0, +double DetectorModel::GetInteractionDensity(GeometryPosition const & p0, std::vector const & targets, std::vector const & total_cross_sections, double const & total_decay_length) const { @@ -604,7 +636,7 @@ double DetectorModel::GetInteractionDensity(Vector3D const & p0, return GetInteractionDensity(intersections, p0, targets, total_cross_sections, total_decay_length); } -double DetectorModel::GetColumnDepthInCGS(Geometry::IntersectionList const & intersections, Vector3D const & p0, Vector3D const & p1) const { +double DetectorModel::GetColumnDepthInCGS(Geometry::IntersectionList const & intersections, GeometryPosition const & p0, GeometryPosition const & p1) const { if(p0 == p1) { return 0.0; } @@ -649,7 +681,7 @@ double DetectorModel::GetColumnDepthInCGS(Geometry::IntersectionList const & int return column_depth * 100; } -double DetectorModel::GetColumnDepthInCGS(Vector3D const & p0, Vector3D const & p1) const { +double DetectorModel::GetColumnDepthInCGS(GeometryPosition const & p0, GeometryPosition const & p1) const { if(p0 == p1) { return 0.0; } @@ -664,7 +696,7 @@ double DetectorModel::GetColumnDepthInCGS(Vector3D const & p0, Vector3D const & return GetColumnDepthInCGS(intersections, p0, p1); } -double DetectorModel::DistanceForColumnDepthFromPoint(Geometry::IntersectionList const & intersections, Vector3D const & p0, Vector3D const & dir, double column_depth) const { +double DetectorModel::DistanceForColumnDepthFromPoint(Geometry::IntersectionList const & intersections, GeometryPosition const & p0, GeometryDirection const & dir, double column_depth) const { Vector3D direction = dir; column_depth /= 100; bool flip = column_depth < 0; @@ -720,20 +752,20 @@ double DetectorModel::DistanceForColumnDepthFromPoint(Geometry::IntersectionList return total_distance; } -double DetectorModel::DistanceForColumnDepthFromPoint(Vector3D const & p0, Vector3D const & direction, double column_depth) const { +double DetectorModel::DistanceForColumnDepthFromPoint(GeometryPosition const & p0, GeometryDirection const & direction, double column_depth) const { Geometry::IntersectionList intersections = GetIntersections(p0, direction); return DistanceForColumnDepthFromPoint(intersections, p0, direction, column_depth); } -double DetectorModel::DistanceForColumnDepthToPoint(Geometry::IntersectionList const & intersections, Vector3D const & p0, Vector3D const & direction, double column_depth) const { +double DetectorModel::DistanceForColumnDepthToPoint(Geometry::IntersectionList const & intersections, GeometryPosition const & p0, GeometryDirection const & direction, double column_depth) const { return DistanceForColumnDepthFromPoint(intersections, p0, -direction, column_depth); } -double DetectorModel::DistanceForColumnDepthToPoint(Vector3D const & p0, Vector3D const & direction, double column_depth) const { +double DetectorModel::DistanceForColumnDepthToPoint(GeometryPosition const & p0, GeometryDirection const & direction, double column_depth) const { return DistanceForColumnDepthFromPoint(p0, -direction, column_depth); } -double DetectorModel::GetMassDensity(Geometry::IntersectionList const & intersections, Vector3D const & p0, std::set targets) const { +double DetectorModel::GetMassDensity(Geometry::IntersectionList const & intersections, GeometryPosition const & p0, std::set targets) const { Vector3D direction = p0 - intersections.position; if(direction.magnitude() == 0) { direction = intersections.direction; @@ -775,13 +807,13 @@ double DetectorModel::GetMassDensity(Geometry::IntersectionList const & intersec return density; } -double DetectorModel::GetMassDensity(Vector3D const & p0, std::set targets) const { +double DetectorModel::GetMassDensity(GeometryPosition const & p0, std::set targets) const { Vector3D direction(1,0,0); // Any direction will work for determining the sector heirarchy Geometry::IntersectionList intersections = GetIntersections(p0, direction); return GetMassDensity(intersections, p0, targets); } -std::vector DetectorModel::GetParticleDensity(Geometry::IntersectionList const & intersections, Vector3D const & p0, std::set targets) const { +std::vector DetectorModel::GetParticleDensity(Geometry::IntersectionList const & intersections, GeometryPosition const & p0, std::set targets) const { Vector3D direction = p0 - intersections.position; if(direction.magnitude() == 0) { direction = intersections.direction; @@ -827,13 +859,13 @@ std::vector DetectorModel::GetParticleDensity(Geometry::IntersectionList return particle_fractions; } -std::vector DetectorModel::GetParticleDensity(Vector3D const & p0, std::set targets) const { +std::vector DetectorModel::GetParticleDensity(GeometryPosition const & p0, std::set targets) const { Vector3D direction(1,0,0); // Any direction will work for determining the sector heirarchy Geometry::IntersectionList intersections = GetIntersections(p0, direction); return GetParticleDensity(intersections, p0, targets); } -double DetectorModel::GetInteractionDepthInCGS(Geometry::IntersectionList const & intersections, Vector3D const & p0, Vector3D const & p1, +double DetectorModel::GetInteractionDepthInCGS(Geometry::IntersectionList const & intersections, GeometryPosition const & p0, GeometryPosition const & p1, std::vector const & targets, std::vector const & total_cross_sections, double const & total_decay_length) const { @@ -898,7 +930,7 @@ double DetectorModel::GetInteractionDepthInCGS(Geometry::IntersectionList const return interaction_depth; } -std::vector DetectorModel::GetParticleColumnDepth(Geometry::IntersectionList const & intersections, Vector3D const & p0, Vector3D const & p1, std::vector const & targets) const { +std::vector DetectorModel::GetParticleColumnDepth(Geometry::IntersectionList const & intersections, GeometryPosition const & p0, GeometryPosition const & p1, std::vector const & targets) const { if(p0 == p1) { return std::vector(targets.size(), 0.0); } @@ -947,7 +979,7 @@ std::vector DetectorModel::GetParticleColumnDepth(Geometry::Intersection return target_counts; } -double DetectorModel::GetInteractionDepthInCGS(Vector3D const & p0, Vector3D const & p1, +double DetectorModel::GetInteractionDepthInCGS(GeometryPosition const & p0, GeometryPosition const & p1, std::vector const & targets, std::vector const & total_cross_sections, double const & total_decay_length) const { @@ -965,7 +997,7 @@ double DetectorModel::GetInteractionDepthInCGS(Vector3D const & p0, Vector3D con return GetInteractionDepthInCGS(intersections, p0, p1, targets, total_cross_sections, total_decay_length); } -DetectorSector DetectorModel::GetContainingSector(Geometry::IntersectionList const & intersections, Vector3D const & p0) const { +DetectorSector DetectorModel::GetContainingSector(Geometry::IntersectionList const & intersections, GeometryPosition const & p0) const { Vector3D direction = intersections.direction; double offset = (intersections.position - p0) * direction; @@ -996,13 +1028,13 @@ DetectorSector DetectorModel::GetContainingSector(Geometry::IntersectionList con return sector; } -DetectorSector DetectorModel::GetContainingSector(Vector3D const & p0) const { +DetectorSector DetectorModel::GetContainingSector(GeometryPosition const & p0) const { Vector3D direction(0, 0, 1); Geometry::IntersectionList intersections = GetIntersections(p0, direction); return GetContainingSector(intersections, p0); } -Geometry::IntersectionList DetectorModel::GetIntersections(Vector3D const & p0, Vector3D const & direction) const { +Geometry::IntersectionList DetectorModel::GetIntersections(GeometryPosition const & p0, GeometryPosition const & direction) const { Geometry::IntersectionList intersections; intersections.position = p0; intersections.direction = direction; @@ -1078,17 +1110,17 @@ Geometry::IntersectionList DetectorModel::GetOuterBounds(Geometry::IntersectionL return result; } -Geometry::IntersectionList DetectorModel::GetOuterBounds(Vector3D const & p0, Vector3D const & direction) const { +Geometry::IntersectionList DetectorModel::GetOuterBounds(GeometryPosition const & p0, GeometryDirection const & direction) const { Geometry::IntersectionList intersections = GetIntersections(p0, direction); return GetOuterBounds(intersections); } -std::set DetectorModel::GetAvailableTargets(std::array const & vertex) const { +std::set DetectorModel::GetAvailableTargets(GeometryPosition const & vertex) const { Geometry::IntersectionList intersections = GetIntersections(vertex, math::Vector3D(0,0,1)); return GetAvailableTargets(intersections, vertex); } -std::set DetectorModel::GetAvailableTargets(geometry::Geometry::IntersectionList const & intersections, std::array const & vertex) const { +std::set DetectorModel::GetAvailableTargets(geometry::Geometry::IntersectionList const & intersections, GeometryPosition const & vertex) const { int matID = GetContainingSector(intersections, Vector3D(vertex[0],vertex[1],vertex[2])).material_id; std::vector particles = materials_.GetMaterialConstituents(matID); return std::set(particles.begin(), particles.end()); @@ -1180,7 +1212,7 @@ void DetectorModel::SectorLoop(std::function const & targets, std::vector const & total_cross_sections, double const & total_decay_length) const { @@ -1264,7 +1296,7 @@ double DetectorModel::DistanceForInteractionDepthFromPoint(Geometry::Intersectio return total_distance; } -double DetectorModel::DistanceForInteractionDepthFromPoint(Vector3D const & p0, Vector3D const & direction, double interaction_depth, +double DetectorModel::DistanceForInteractionDepthFromPoint(GeometryPosition const & p0, GeometryDirection const & direction, double interaction_depth, std::vector const & targets, std::vector const & total_cross_sections, double const & total_decay_length) const { @@ -1272,36 +1304,166 @@ double DetectorModel::DistanceForInteractionDepthFromPoint(Vector3D const & p0, return DistanceForInteractionDepthFromPoint(intersections, p0, direction, interaction_depth, targets, total_cross_sections, total_decay_length); } -double DetectorModel::DistanceForInteractionDepthToPoint(Geometry::IntersectionList const & intersections, Vector3D const & p0, Vector3D const & direction, double interaction_depth, +double DetectorModel::DistanceForInteractionDepthToPoint(Geometry::IntersectionList const & intersections, GeometryPosition const & p0, GeometryDirection const & direction, double interaction_depth, std::vector const & targets, std::vector const & total_cross_sections, double const & total_decay_length) const { return DistanceForInteractionDepthFromPoint(intersections, p0, -direction, interaction_depth, targets, total_cross_sections, total_decay_length); } -double DetectorModel::DistanceForInteractionDepthToPoint(Vector3D const & p0, Vector3D const & direction, double interaction_depth, +double DetectorModel::DistanceForInteractionDepthToPoint(GeometryPosition const & p0, GeometryDirection const & direction, double interaction_depth, std::vector const & targets, std::vector const & total_cross_sections, double const & total_decay_length) const { return DistanceForInteractionDepthFromPoint(p0, -direction, interaction_depth, targets, total_cross_sections, total_decay_length); } -Vector3D DetectorModel::GetEarthCoordPosFromDetCoordPos(Vector3D const & point) const { - return point + detector_origin_; +////////////////////////////////////////////////////////// + + +double DetectorModel::GetMassDensity(Geometry::IntersectionList const & intersections, DetectorPosition const & p0) const { + return GetMassDensity(intersections, ToGeo(p0)); +} + +double DetectorModel::GetMassDensity(DetectorPosition const & p0) const { + return GetMassDensity(ToGeo(p0)); +} + +double DetectorModel::GetParticleDensity(Geometry::IntersectionList const & intersections, DetectorPosition const & p0, LI::dataclasses::Particle::ParticleType target) const { + return GetParticleDensity(intersections, ToGeo(p0), target); +} + +double DetectorModel::GetParticleDensity(DetectorPosition const & p0, LI::dataclasses::Particle::ParticleType target) const { + return GetParticleDensity(ToGeo(p0), target); +} + +double DetectorModel::GetInteractionDensity(Geometry::IntersectionList const & intersections, DetectorPosition const & p0, + std::vector const & targets, + std::vector const & total_cross_sections, + double const & total_decay_length) const { + return GetInteractionDensity(intersections, ToGeo(p0), targets, total_cross_sections, total_decay_length); +} + +double DetectorModel::GetInteractionDensity(DetectorPosition const & p0, + std::vector const & targets, + std::vector const & total_cross_sections, + double const & total_decay_length) const { + return GetInteractionDensity(ToGeo(p0), targets, total_cross_sections, total_decay_length); +} + +double DetectorModel::GetColumnDepthInCGS(Geometry::IntersectionList const & intersections, DetectorPosition const & p0, DetectorPosition const & p1) const { + return GetColumnDepthInCGS(intersections, ToGeo(p0), ToGeo(p1)); +} + +double DetectorModel::GetColumnDepthInCGS(DetectorPosition const & p0, DetectorPosition const & p1) const { + return GetColumnDepthInCGS(ToGeo(p0), ToGeo(p1)); +} + +double DetectorModel::DistanceForColumnDepthFromPoint(Geometry::IntersectionList const & intersections, DetectorPosition const & p0, DetectorDirection const & dir, double column_depth) const { + return DistanceForColumnDepthFromPoint(intersections, ToGeo(p0), ToGeo(dir), column_depth); +} + +double DetectorModel::DistanceForColumnDepthFromPoint(DetectorPosition const & p0, DetectorDirection const & direction, double column_depth) const { + return DistanceForColumnDepthFromPoint(ToGeo(p0), ToGeo(direction), column_depth); +} + +double DetectorModel::DistanceForColumnDepthToPoint(Geometry::IntersectionList const & intersections, DetectorPosition const & p0, DetectorDirection const & direction, double column_depth) const { + return DistanceForColumnDepthFromPoint(intersections, ToGeo(p0), ToGeo(direction), column_depth); +} + +double DetectorModel::DistanceForColumnDepthToPoint(DetectorPosition const & p0, DetectorDirection const & direction, double column_depth) const { + return DistanceForColumnDepthFromPoint(ToGeo(p0), ToGeo(direction), column_depth); +} + +double DetectorModel::GetMassDensity(Geometry::IntersectionList const & intersections, DetectorPosition const & p0, std::set targets) const { + return GetMassDensity(intersections, ToGeo(p0), targets); +} + +double DetectorModel::GetMassDensity(DetectorPosition const & p0, std::set targets) const { + return GetMassDensity(ToGeo(p0), targets); +} + +std::vector DetectorModel::GetParticleDensity(Geometry::IntersectionList const & intersections, DetectorPosition const & p0, std::set targets) const { + return GetParticleDensity(intersections, ToGeo(p0), targets); +} + +std::vector DetectorModel::GetParticleDensity(DetectorPosition const & p0, std::set targets) const { + return GetParticleDensity(ToGeo(p0), targets); +} + +double DetectorModel::GetInteractionDepthInCGS(Geometry::IntersectionList const & intersections, DetectorPosition const & p0, DetectorPosition const & p1, + std::vector const & targets, + std::vector const & total_cross_sections, + double const & total_decay_length) const { + return GetInteractionDepthInCGS(intersections, ToGeo(p0), ToGeo(p1), targets, total_cross_sections, total_decay_length); +} + +std::vector DetectorModel::GetParticleColumnDepth(Geometry::IntersectionList const & intersections, DetectorPosition const & p0, DetectorPosition const & p1, std::vector const & targets) const { + return GetParticleColumnDepth(intersections, ToGeo(p0), ToGeo(p1), targets); +} + +double DetectorModel::GetInteractionDepthInCGS(DetectorPosition const & p0, DetectorPosition const & p1, + std::vector const & targets, + std::vector const & total_cross_sections, + double const & total_decay_length) const { + return GetInteractionDepthInCGS(ToGeo(p0), ToGeo(p1), targets, total_cross_sections, total_decay_length); +} + +DetectorSector DetectorModel::GetContainingSector(Geometry::IntersectionList const & intersections, DetectorPosition const & p0) const { + return GetContainingSector(intersections, ToGeo(p0)); } -Vector3D DetectorModel::GetEarthCoordDirFromDetCoordDir(Vector3D const & direction) const { - return direction; +DetectorSector DetectorModel::GetContainingSector(DetectorPosition const & p0) const { + return GetContainingSector(ToGeo(p0)); } -Vector3D DetectorModel::GetDetCoordPosFromEarthCoordPos(Vector3D const & point) const { - return point - detector_origin_; +Geometry::IntersectionList DetectorModel::GetIntersections(DetectorPosition const & p0, DetectorPosition const & direction) const { + return GetIntersections(ToGeo(p0), ToGeo(direction)); } -Vector3D DetectorModel::GetDetCoordDirFromEarthCoordDir(Vector3D const & direction) const { - return direction; +Geometry::IntersectionList DetectorModel::GetOuterBounds(DetectorPosition const & p0, DetectorDirection const & direction) const { + return GetOuterBounds(ToGeo(p0), ToGeo(direction)); } +std::set DetectorModel::GetAvailableTargets(DetectorPosition const & vertex) const { + return GetAvailableTargets(ToGeo(vertex)); +} + +std::set DetectorModel::GetAvailableTargets(geometry::Geometry::IntersectionList const & intersections, DetectorPosition const & vertex) const { + return GetAvailableTargets(intersections, ToGeo(vertex)); +} + + +double DetectorModel::DistanceForInteractionDepthFromPoint(Geometry::IntersectionList const & intersections, DetectorPosition const & p0, DetectorDirection const & dir, double interaction_depth, + std::vector const & targets, + std::vector const & total_cross_sections, + double const & total_decay_length) const { + return DistanceForInteractionDepthFromPoint(intersections, ToGeo(p0), ToGeo(dir), interaction_depth, targets, total_cross_sections, total_decay_length); +} + +double DetectorModel::DistanceForInteractionDepthFromPoint(DetectorPosition const & p0, DetectorDirection const & direction, double interaction_depth, + std::vector const & targets, + std::vector const & total_cross_sections, + double const & total_decay_length) const { + return DistanceForInteractionDepthFromPoint(ToGeo(p0), ToGeo(direction), interaction_depth, targets, total_cross_sections, total_decay_length); +} + +double DetectorModel::DistanceForInteractionDepthToPoint(Geometry::IntersectionList const & intersections, DetectorPosition const & p0, DetectorDirection const & direction, double interaction_depth, + std::vector const & targets, + std::vector const & total_cross_sections, + double const & total_decay_length) const { + return DistanceForInteractionDepthFromPoint(intersections, ToGeo(p0), ToGeo(direction), interaction_depth, targets, total_cross_sections, total_decay_length); +} + +double DetectorModel::DistanceForInteractionDepthToPoint(DetectorPosition const & p0, DetectorDirection const & direction, double interaction_depth, + std::vector const & targets, + std::vector const & total_cross_sections, + double const & total_decay_length) const { + return DistanceForInteractionDepthFromPoint(ToGeo(p0), ToGeo(direction), interaction_depth, targets, total_cross_sections, total_decay_length); +} + +////////////////////////////////////////////////////////// + void DetectorModel::LoadConcentricShellsFromLegacyFile(std::string model_fname, double detector_depth, double ice_cap_angle) { if(model_fname.empty()) throw(std::runtime_error("Received empty detector model filename!")); diff --git a/projects/detector/public/LeptonInjector/detector/DetectorModel.h b/projects/detector/public/LeptonInjector/detector/DetectorModel.h index 34cea86d2..a44de0388 100644 --- a/projects/detector/public/LeptonInjector/detector/DetectorModel.h +++ b/projects/detector/public/LeptonInjector/detector/DetectorModel.h @@ -28,9 +28,12 @@ #include "LeptonInjector/detector/MaterialModel.h" // for MaterialModel #include "LeptonInjector/geometry/Geometry.h" // for Geometry #include "LeptonInjector/math/Vector3D.h" // for Vector3D +#include "LeptonInjector/math/Quaternion.h" // for Quaternion +#include "LeptonInjector/detector/Coordinates.h" #include "LeptonInjector/dataclasses/Particle.h" +namespace LI { namespace detector { class Path; } } namespace LI { namespace detector { class DensityDistribution; } } namespace LI { @@ -60,11 +63,13 @@ struct DetectorSector { class DetectorModel { private: +friend LI::detector::Path; std::string path_; MaterialModel materials_; std::vector sectors_; std::map sector_map_; math::Vector3D detector_origin_; + math::Quaternion detector_rotation_; public: DetectorModel(); DetectorModel(std::string const & detector_model, std::string const & material_model); @@ -87,74 +92,157 @@ class DetectorModel { void LoadDetectorModel(std::string const & detector_model); void LoadMaterialModel(std::string const & material_model); +private: + double GetMassDensity(geometry::Geometry::IntersectionList const & intersections, GeometryPosition const & p0) const; + double GetMassDensity(GeometryPosition const & p0) const; + double GetParticleDensity(geometry::Geometry::IntersectionList const & intersections, GeometryPosition const & p0, LI::dataclasses::Particle::ParticleType target) const; + double GetParticleDensity(GeometryPosition const & p0, LI::dataclasses::Particle::ParticleType target) const; + double GetInteractionDensity(geometry::Geometry::IntersectionList const & intersections, GeometryPosition const & p0, + std::vector const & targets, + std::vector const & total_cross_sections, + double const & total_decay_length) const; + double GetInteractionDensity(GeometryPosition const & p0, + std::vector const & targets, + std::vector const & total_cross_sections, + double const & total_decay_length) const; + + double GetColumnDepthInCGS(geometry::Geometry::IntersectionList const & intersections, GeometryPosition const & p0, GeometryPosition const & p1) const; + double GetColumnDepthInCGS(GeometryPosition const & p0, GeometryPosition const & p1) const; + double DistanceForColumnDepthFromPoint(geometry::Geometry::IntersectionList const & intersections, GeometryPosition const & end_point, GeometryDirection const & direction, double column_depth) const; + double DistanceForColumnDepthFromPoint(GeometryPosition const & end_point, GeometryDirection const & direction, double column_depth) const; + double DistanceForColumnDepthToPoint(geometry::Geometry::IntersectionList const & intersections, GeometryPosition const & end_point, GeometryDirection const & direction, double column_depth) const; + double DistanceForColumnDepthToPoint(GeometryPosition const & end_point, GeometryDirection const & direction, double column_depth) const; + + // Density/CD calculations with general target list, not just nucleon/electron + double GetMassDensity(geometry::Geometry::IntersectionList const & intersections, GeometryPosition const & p0, std::set targets) const; + double GetMassDensity(GeometryPosition const & p0, std::set targets) const; + template::value, Iterator>::type> + double GetMassDensity(geometry::Geometry::IntersectionList const & intersections, GeometryPosition const & p0, Iterator begin, Iterator end) const; + template::value, Iterator>::type> + double GetMassDensity(GeometryPosition const & p0, Iterator begin, Iterator end) const; + std::vector GetParticleDensity(geometry::Geometry::IntersectionList const & intersections, GeometryPosition const & p0, std::set targets) const; + std::vector GetParticleDensity(GeometryPosition const & p0, std::set targets) const; + template::value, Iterator>::type> + std::vector GetParticleDensity(geometry::Geometry::IntersectionList const & intersections, GeometryPosition const & p0, Iterator begin, Iterator end) const; + template::value, Iterator>::type> + std::vector GetParticleDensity(GeometryPosition const & p0, Iterator begin, Iterator end) const; + + double GetInteractionDepthInCGS(geometry::Geometry::IntersectionList const & intersections, GeometryPosition const & p0, GeometryPosition const & p1, + std::vector const & targets, + std::vector const & total_cross_sections, + double const & total_decay_length) const; + double GetInteractionDepthInCGS(GeometryPosition const & p0, GeometryPosition const & p1, + std::vector const & targets, + std::vector const & total_cross_sections, + double const & total_decay_length) const; + double DistanceForInteractionDepthFromPoint(geometry::Geometry::IntersectionList const & intersections, GeometryPosition const & end_point, GeometryDirection const & direction, double interaction_depth, + std::vector const & targets, + std::vector const & total_cross_sections, + double const & total_decay_length) const; + double DistanceForInteractionDepthFromPoint(GeometryPosition const & end_point, GeometryDirection const & direction, double interaction_depth, + std::vector const & targets, + std::vector const & total_cross_sections, + double const & total_decay_length) const; + double DistanceForInteractionDepthToPoint(geometry::Geometry::IntersectionList const & intersections, GeometryPosition const & end_point, GeometryDirection const & direction, double interaction_depth, + std::vector const & targets, + std::vector const & total_cross_sections, + double const & total_decay_length) const; + double DistanceForInteractionDepthToPoint(GeometryPosition const & end_point, GeometryDirection const & direction, double interaction_depth, + std::vector const & targets, + std::vector const & total_cross_sections, + double const & total_decay_length) const; + + std::vector GetParticleColumnDepth(geometry::Geometry::IntersectionList const & intersections, GeometryPosition const & p0, GeometryPosition const & p1, std::vector const & targets) const; + + DetectorSector GetContainingSector(geometry::Geometry::IntersectionList const & intersections, GeometryPosition const & p0) const; + DetectorSector GetContainingSector(GeometryPosition const & p0) const; + + geometry::Geometry::IntersectionList GetIntersections(GeometryPosition const & p0, GeometryDirection const & direction) const; + geometry::Geometry::IntersectionList GetOuterBounds(GeometryPosition const & p0, GeometryDirection const & direction) const; + + std::set GetAvailableTargets(GeometryPosition const & vertex) const; + std::set GetAvailableTargets(geometry::Geometry::IntersectionList const & intersections, GeometryPosition const & vertex) const; - double GetMassDensity(geometry::Geometry::IntersectionList const & intersections, math::Vector3D const & p0) const; - double GetMassDensity(math::Vector3D const & p0) const; - double GetParticleDensity(geometry::Geometry::IntersectionList const & intersections, math::Vector3D const & p0, LI::dataclasses::Particle::ParticleType target) const; - double GetParticleDensity(math::Vector3D const & p0, LI::dataclasses::Particle::ParticleType target) const; - double GetInteractionDensity(geometry::Geometry::IntersectionList const & intersections, math::Vector3D const & p0, +public: + double GetMassDensity(geometry::Geometry::IntersectionList const & intersections, DetectorPosition const & p0) const; + double GetMassDensity(DetectorPosition const & p0) const; + double GetParticleDensity(geometry::Geometry::IntersectionList const & intersections, DetectorPosition const & p0, LI::dataclasses::Particle::ParticleType target) const; + double GetParticleDensity(DetectorPosition const & p0, LI::dataclasses::Particle::ParticleType target) const; + double GetInteractionDensity(geometry::Geometry::IntersectionList const & intersections, DetectorPosition const & p0, std::vector const & targets, std::vector const & total_cross_sections, double const & total_decay_length) const; - double GetInteractionDensity(math::Vector3D const & p0, + double GetInteractionDensity(DetectorPosition const & p0, std::vector const & targets, std::vector const & total_cross_sections, double const & total_decay_length) const; - double GetColumnDepthInCGS(geometry::Geometry::IntersectionList const & intersections, math::Vector3D const & p0, math::Vector3D const & p1) const; - double GetColumnDepthInCGS(math::Vector3D const & p0, math::Vector3D const & p1) const; - double DistanceForColumnDepthFromPoint(geometry::Geometry::IntersectionList const & intersections, math::Vector3D const & end_point, math::Vector3D const & direction, double column_depth) const; - double DistanceForColumnDepthFromPoint(math::Vector3D const & end_point, math::Vector3D const & direction, double column_depth) const; - double DistanceForColumnDepthToPoint(geometry::Geometry::IntersectionList const & intersections, math::Vector3D const & end_point, math::Vector3D const & direction, double column_depth) const; - double DistanceForColumnDepthToPoint(math::Vector3D const & end_point, math::Vector3D const & direction, double column_depth) const; + double GetColumnDepthInCGS(geometry::Geometry::IntersectionList const & intersections, DetectorPosition const & p0, DetectorPosition const & p1) const; + double GetColumnDepthInCGS(DetectorPosition const & p0, DetectorPosition const & p1) const; + double DistanceForColumnDepthFromPoint(geometry::Geometry::IntersectionList const & intersections, DetectorPosition const & end_point, DetectorDirection const & direction, double column_depth) const; + double DistanceForColumnDepthFromPoint(DetectorPosition const & end_point, DetectorDirection const & direction, double column_depth) const; + double DistanceForColumnDepthToPoint(geometry::Geometry::IntersectionList const & intersections, DetectorPosition const & end_point, DetectorDirection const & direction, double column_depth) const; + double DistanceForColumnDepthToPoint(DetectorPosition const & end_point, DetectorDirection const & direction, double column_depth) const; // Density/CD calculations with general target list, not just nucleon/electron - double GetMassDensity(geometry::Geometry::IntersectionList const & intersections, math::Vector3D const & p0, std::set targets) const; - double GetMassDensity(math::Vector3D const & p0, std::set targets) const; + double GetMassDensity(geometry::Geometry::IntersectionList const & intersections, DetectorPosition const & p0, std::set targets) const; + double GetMassDensity(DetectorPosition const & p0, std::set targets) const; template::value, Iterator>::type> - double GetMassDensity(geometry::Geometry::IntersectionList const & intersections, math::Vector3D const & p0, Iterator begin, Iterator end) const; + double GetMassDensity(geometry::Geometry::IntersectionList const & intersections, DetectorPosition const & p0, Iterator begin, Iterator end) const; template::value, Iterator>::type> - double GetMassDensity(math::Vector3D const & p0, Iterator begin, Iterator end) const; - std::vector GetParticleDensity(geometry::Geometry::IntersectionList const & intersections, math::Vector3D const & p0, std::set targets) const; - std::vector GetParticleDensity(math::Vector3D const & p0, std::set targets) const; + double GetMassDensity(DetectorPosition const & p0, Iterator begin, Iterator end) const; + std::vector GetParticleDensity(geometry::Geometry::IntersectionList const & intersections, DetectorPosition const & p0, std::set targets) const; + std::vector GetParticleDensity(DetectorPosition const & p0, std::set targets) const; template::value, Iterator>::type> - std::vector GetParticleDensity(geometry::Geometry::IntersectionList const & intersections, math::Vector3D const & p0, Iterator begin, Iterator end) const; + std::vector GetParticleDensity(geometry::Geometry::IntersectionList const & intersections, DetectorPosition const & p0, Iterator begin, Iterator end) const; template::value, Iterator>::type> - std::vector GetParticleDensity(math::Vector3D const & p0, Iterator begin, Iterator end) const; + std::vector GetParticleDensity(DetectorPosition const & p0, Iterator begin, Iterator end) const; - double GetInteractionDepthInCGS(geometry::Geometry::IntersectionList const & intersections, math::Vector3D const & p0, math::Vector3D const & p1, + double GetInteractionDepthInCGS(geometry::Geometry::IntersectionList const & intersections, DetectorPosition const & p0, DetectorPosition const & p1, std::vector const & targets, std::vector const & total_cross_sections, double const & total_decay_length) const; - double GetInteractionDepthInCGS(math::Vector3D const & p0, math::Vector3D const & p1, + double GetInteractionDepthInCGS(DetectorPosition const & p0, DetectorPosition const & p1, std::vector const & targets, std::vector const & total_cross_sections, double const & total_decay_length) const; - double DistanceForInteractionDepthFromPoint(geometry::Geometry::IntersectionList const & intersections, math::Vector3D const & end_point, math::Vector3D const & direction, double interaction_depth, + double DistanceForInteractionDepthFromPoint(geometry::Geometry::IntersectionList const & intersections, DetectorPosition const & end_point, DetectorDirection const & direction, double interaction_depth, std::vector const & targets, std::vector const & total_cross_sections, double const & total_decay_length) const; - double DistanceForInteractionDepthFromPoint(math::Vector3D const & end_point, math::Vector3D const & direction, double interaction_depth, + double DistanceForInteractionDepthFromPoint(DetectorPosition const & end_point, DetectorDirection const & direction, double interaction_depth, std::vector const & targets, std::vector const & total_cross_sections, double const & total_decay_length) const; - double DistanceForInteractionDepthToPoint(geometry::Geometry::IntersectionList const & intersections, math::Vector3D const & end_point, math::Vector3D const & direction, double interaction_depth, + double DistanceForInteractionDepthToPoint(geometry::Geometry::IntersectionList const & intersections, DetectorPosition const & end_point, DetectorDirection const & direction, double interaction_depth, std::vector const & targets, std::vector const & total_cross_sections, double const & total_decay_length) const; - double DistanceForInteractionDepthToPoint(math::Vector3D const & end_point, math::Vector3D const & direction, double interaction_depth, + double DistanceForInteractionDepthToPoint(DetectorPosition const & end_point, DetectorDirection const & direction, double interaction_depth, std::vector const & targets, std::vector const & total_cross_sections, double const & total_decay_length) const; - std::vector GetParticleColumnDepth(geometry::Geometry::IntersectionList const & intersections, math::Vector3D const & p0, math::Vector3D const & p1, std::vector const & targets) const; + std::vector GetParticleColumnDepth(geometry::Geometry::IntersectionList const & intersections, DetectorPosition const & p0, DetectorPosition const & p1, std::vector const & targets) const; + + DetectorSector GetContainingSector(geometry::Geometry::IntersectionList const & intersections, DetectorPosition const & p0) const; + DetectorSector GetContainingSector(DetectorPosition const & p0) const; + + geometry::Geometry::IntersectionList GetIntersections(DetectorPosition const & p0, DetectorDirection const & direction) const; + geometry::Geometry::IntersectionList GetOuterBounds(DetectorPosition const & p0, DetectorDirection const & direction) const; + + std::set GetAvailableTargets(DetectorPosition const & vertex) const; + std::set GetAvailableTargets(geometry::Geometry::IntersectionList const & intersections, DetectorPosition const & vertex) const; + + DetectorPosition ToDet(GeometryPosition const &) const; + DetectorDirection ToDet(GeometryDirection const &) const; + DetectorPosition ToDet(GeometryPosition &&) const; + DetectorDirection ToDet(GeometryDirection &&) const; + GeometryPosition ToGeo(DetectorPosition const &) const; + GeometryDirection ToGeo(DetectorDirection const &) const; + GeometryPosition ToGeo(DetectorPosition &&) const; + GeometryDirection ToGeo(DetectorDirection &&) const; - DetectorSector GetContainingSector(geometry::Geometry::IntersectionList const & intersections, math::Vector3D const & p0) const; - DetectorSector GetContainingSector(math::Vector3D const & p0) const; - math::Vector3D GetEarthCoordPosFromDetCoordPos(math::Vector3D const & point) const; - math::Vector3D GetEarthCoordDirFromDetCoordDir(math::Vector3D const & direction) const; - math::Vector3D GetDetCoordPosFromEarthCoordPos(math::Vector3D const & point) const; - math::Vector3D GetDetCoordDirFromEarthCoordDir(math::Vector3D const & direction) const; std::string GetPath() const; void SetPath(std::string const & path); @@ -165,23 +253,19 @@ class DetectorModel { std::vector const & GetSectors() const; void SetSectors(std::vector const & sectors); - math::Vector3D GetDetectorOrigin() const; - void SetDetectorOrigin(math::Vector3D const & detector_origin); + GeometryPosition GetDetectorOrigin() const; + void SetDetectorOrigin(GeometryPosition const & detector_origin); void AddSector(DetectorSector sector); DetectorSector GetSector(int level) const; void ClearSectors(); - geometry::Geometry::IntersectionList GetIntersections(math::Vector3D const & p0, math::Vector3D const & direction) const; static void SortIntersections(geometry::Geometry::IntersectionList & intersections); static void SortIntersections(std::vector & intersections); static void SectorLoop(std::function::const_iterator, std::vector::const_iterator, double)> callback, geometry::Geometry::IntersectionList const & intersections, bool reverse=false); static geometry::Geometry::IntersectionList GetOuterBounds(geometry::Geometry::IntersectionList const & intersections); - geometry::Geometry::IntersectionList GetOuterBounds(math::Vector3D const & p0, math::Vector3D const & direction) const; - std::set GetAvailableTargets(std::array const & vertex) const; - std::set GetAvailableTargets(geometry::Geometry::IntersectionList const & intersections, std::array const & vertex) const; private: void LoadDefaultMaterials(); diff --git a/projects/detector/public/LeptonInjector/detector/DetectorModel.tcc b/projects/detector/public/LeptonInjector/detector/DetectorModel.tcc index 9be16ed40..6e8f2d1d4 100644 --- a/projects/detector/public/LeptonInjector/detector/DetectorModel.tcc +++ b/projects/detector/public/LeptonInjector/detector/DetectorModel.tcc @@ -5,6 +5,7 @@ #include #include +#include "LeptonInjector/detector/Coordinates.h" #include "LeptonInjector/detector/DetectorModel.h" #include "LeptonInjector/detector/DensityDistribution.h" @@ -12,7 +13,7 @@ namespace LI { namespace detector { template -double DetectorModel::GetMassDensity(geometry::Geometry::IntersectionList const & intersections, math::Vector3D const & p0, Iterator begin, Iterator end) const { +double DetectorModel::GetMassDensity(geometry::Geometry::IntersectionList const & intersections, GeometryPosition const & p0, Iterator begin, Iterator end) const { math::Vector3D direction = p0 - intersections.position; if(direction.magnitude() == 0) { direction = intersections.direction; @@ -55,15 +56,15 @@ double DetectorModel::GetMassDensity(geometry::Geometry::IntersectionList const } template -double DetectorModel::GetMassDensity(math::Vector3D const & p0, Iterator begin, Iterator end) const { +double DetectorModel::GetMassDensity(GeometryPosition const & p0, Iterator begin, Iterator end) const { math::Vector3D direction(1,0,0); // Any direction will work for determining the sector heirarchy - geometry::Geometry::IntersectionList intersections = GetIntersections(p0, direction); + geometry::Geometry::IntersectionList intersections = GetIntersections(p0, GeometryDirection(direction)); return GetMassDensity(intersections, p0, begin, end); } template -std::vector DetectorModel::GetParticleDensity(geometry::Geometry::IntersectionList const & intersections, math::Vector3D const & p0, Iterator begin, Iterator end) const { +std::vector DetectorModel::GetParticleDensity(geometry::Geometry::IntersectionList const & intersections, GeometryPosition const & p0, Iterator begin, Iterator end) const { math::Vector3D direction = p0 - intersections.position; if(direction.magnitude() == 0) { direction = intersections.direction; @@ -110,12 +111,31 @@ std::vector DetectorModel::GetParticleDensity(geometry::Geometry::Inters } template -std::vector DetectorModel::GetParticleDensity(math::Vector3D const & p0, Iterator begin, Iterator end) const { +std::vector DetectorModel::GetParticleDensity(GeometryPosition const & p0, Iterator begin, Iterator end) const { math::Vector3D direction(1,0,0); // Any direction will work for determining the sector heirarchy - geometry::Geometry::IntersectionList intersections = GetIntersections(p0, direction); + geometry::Geometry::IntersectionList intersections = GetIntersections(p0, GeometryDirection(direction)); return GetParticleDensity(intersections, p0, begin, end); } +template +double DetectorModel::GetMassDensity(geometry::Geometry::IntersectionList const & intersections, DetectorPosition const & p0, Iterator begin, Iterator end) const { + return GetMassDensity(intersections, ToGeo(p0), begin, end); +} + +template +double DetectorModel::GetMassDensity(DetectorPosition const & p0, Iterator begin, Iterator end) const { + return GetMassDensity(ToGeo(p0), begin, end); +} + +template +std::vector DetectorModel::GetParticleDensity(geometry::Geometry::IntersectionList const & intersections, DetectorPosition const & p0, Iterator begin, Iterator end) const { + return GetParticleDensity(intersections, ToGeo(p0), begin, end); +} + +template +std::vector DetectorModel::GetParticleDensity(DetectorPosition const & p0, Iterator begin, Iterator end) const { + return GetParticleDensity(ToGeo(p0), begin, end); +} } // namespace detector } // namespace LI From c12fd1de3a73497283e41e4f8bc8a966f785a30f Mon Sep 17 00:00:00 2001 From: Austin Schneider Date: Mon, 15 Jan 2024 18:40:11 -0700 Subject: [PATCH 05/42] Sketch new interface for Path --- projects/detector/private/Path.cxx | 136 ++++++++++++++++-- .../public/LeptonInjector/detector/Path.h | 41 ++++-- 2 files changed, 150 insertions(+), 27 deletions(-) diff --git a/projects/detector/private/Path.cxx b/projects/detector/private/Path.cxx index d50d4e01e..9c05d993d 100644 --- a/projects/detector/private/Path.cxx +++ b/projects/detector/private/Path.cxx @@ -14,6 +14,20 @@ namespace LI { namespace detector { +void Path::UpdatePoints() { + if((not set_points_) and set_det_points_ and set_detector_model_) { + first_point_ = detector_model_->ToGeo(first_point_det_); + last_point_ = detector_model_->ToGeo(last_point_det_); + direction_ = detector_model_->ToGeo(direction_det_); + set_points_ = true; + } else if ((not set_det_points_) and set_points_ and set_detector_model_) { + first_point_det_ = detector_model_->ToDet(first_point_); + last_point_det_ = detector_model_->ToDet(last_point_); + direction_det_ = detector_model_->ToDet(direction_); + set_det_points_ = true; + } +} + Path::Path() { } @@ -52,20 +66,38 @@ std::shared_ptr Path::GetDetectorModel() { return detector_model_; } -math::Vector3D const & Path::GetFirstPoint() { +DetectorPosition const & Path::GetFirstPoint() { + UpdatePoints(); + return first_point_det_; +} + +DetectorPosition const & Path::GetLastPoint() { + UpdatePoints(); + return last_point_det_; +} + +DetectorPosition const & Path::GetDirection() { + UpdatePoints(); + return direction_det_; +} + +GeometryPosition const & Path::GetGeoFirstPoint() { + UpdatePoints(); return first_point_; } -math::Vector3D const & Path::GetLastPoint() { +GeometryPosition const & Path::GetGeoLastPoint() { + UpdatePoints(); return last_point_; } -math::Vector3D const & Path::GetDirection() { +GeometryPosition const & Path::GetGeoDirection() { + UpdatePoints(); return direction_; } double Path::GetDistance() { - return distance_; + return distance__; } geometry::Geometry::IntersectionList const & Path::GetIntersections() { @@ -73,8 +105,12 @@ geometry::Geometry::IntersectionList const & Path::GetIntersections() { } void Path::SetDetectorModel(std::shared_ptr detector_model) { + if(set_detector_model_ and set_det_points_) { + set_points_ = false; + } detector_model_ = detector_model; set_detector_model_ = true; + UpdatePoints(); } void Path::EnsureDetectorModel() { @@ -83,18 +119,33 @@ void Path::EnsureDetectorModel() { } } -void Path::SetPoints(math::Vector3D first_point, math::Vector3D last_point) { +void Path::SetPoints(GeometryPosition first_point, GeometryPosition last_point) { first_point_ = first_point; last_point_ = last_point; direction_ = last_point_ - first_point_; distance_ = direction_.magnitude(); direction_.normalize(); set_points_ = true; + set_det_points_ = false; + set_intersections_ = false; + set_column_depth_ = false; + UpdatePoints(); +} + +void Path::SetPoints(DetectorPosition first_point, DetectorPosition last_point) { + first_point_det_ = first_point; + last_point_det_ = last_point; + direction_det_ = last_point_det_ - first_point_det_; + distance_ = direction_det_.magnitude(); + direction_det_.normalize(); + set_points_ = false; + set_det_points_ = true; set_intersections_ = false; set_column_depth_ = false; + UpdatePoints(); } -void Path::SetPointsWithRay(math::Vector3D first_point, math::Vector3D direction, double distance) { +void Path::SetPointsWithRay(GeometryPosition first_point, GeometryPosition direction, double distance) { first_point_ = first_point; direction_ = direction; direction_.normalize(); @@ -103,11 +154,29 @@ void Path::SetPointsWithRay(math::Vector3D first_point, math::Vector3D direction distance_ = distance; last_point_ = first_point + direction * distance; set_points_ = true; + set_det_points_ = false; + set_intersections_ = false; + set_column_depth_ = false; + UpdatePoints(); +} + +void Path::SetPointsWithRay(DetectorPosition first_point, DetectorPosition direction, double distance) { + first_point_det_ = first_point; + direction_det_ = direction; + direction_det_.normalize(); + //double dif = std::abs(direction_.magnitude() - direction.magnitude()) / std::max(direction_.magnitude(), direction.magnitude()); + //if(not std::isnan(dif)) assert(dif < 1e-12); + distance_det_ = distance; + last_point_det_ = first_point + direction * distance; + set_points_ = false; + set_det_points_ = true; set_intersections_ = false; set_column_depth_ = false; + UpdatePoints(); } void Path::EnsurePoints() { + UpdatePoints(); if(not set_points_) { throw(std::runtime_error("Points not set!")); } @@ -159,6 +228,7 @@ void Path::ClipToOuterBounds() { distance_ = (last_point_ - first_point_).magnitude(); set_column_depth_ = false; } + set_det_points_ = false; } else { return; } @@ -168,6 +238,7 @@ void Path::Flip() { EnsurePoints(); std::swap(first_point_, last_point_); direction_ *= -1; + set_det_points_ = false; } @@ -183,6 +254,7 @@ void Path::ExtendFromEndByDistance(double distance) { last_point_ = first_point_; } set_column_depth_ = false; + set_det_points_ = false; } void Path::ExtendFromStartByDistance(double distance) { @@ -194,6 +266,7 @@ void Path::ExtendFromStartByDistance(double distance) { first_point_ = last_point_; } set_column_depth_ = false; + set_det_points_ = false; } void Path::ShrinkFromEndByDistance(double distance) { @@ -651,17 +724,50 @@ double Path::GetDistanceFromEndInReverse(double interaction_depth, } -bool Path::IsWithinBounds(math::Vector3D point) { - EnsurePoints(); - double d0 = LI::math::scalar_product(direction_, first_point_ - point); - double d1 = LI::math::scalar_product(direction_, last_point_ - point); - return d0 <= 0 and d1 >= 0; +bool Path::IsWithinBounds(GeometryPosition point) { + UpdatePoints(); + if(set_points_) { + double d0 = LI::math::scalar_product(direction_, first_point_ - point); + double d1 = LI::math::scalar_product(direction_, last_point_ - point); + return d0 <= 0 and d1 >= 0; + } else { + EnsurePoints(); + } } -double Path::GetDistanceFromStartInBounds(math::Vector3D point) { - EnsurePoints(); - double d0 = LI::math::scalar_product(direction_, point - first_point_); - return std::max(0.0, d0); +double Path::GetDistanceFromStartInBounds(GeometryPosition point) { + UpdatePoints(); + if(set_points_) { + double d0 = LI::math::scalar_product(direction_, point - first_point_); + return std::max(0.0, d0); + } else { + EnsurePoints(); + } +} + +bool Path::IsWithinBounds(DetectorPosition point) { + UpdatePoints(); + if(set_det_points_) { + double d0 = LI::math::scalar_product(direction_det_, first_point_det_ - point); + double d1 = LI::math::scalar_product(direction_det_, last_point_det_ - point); + return d0 <= 0 and d1 >= 0; + } else if(set_points_ and set_detector_model_) { + return IsWithinBounds(detector_model_->ToGeo(point)); + } else { + throw(std::runtime_error("Detector points not set!")); + } +} + +double Path::GetDistanceFromStartInBounds(DetectorPosition point) { + UpdatePoints(); + if(set_det_points) { + double d0 = LI::math::scalar_product(direction_det_, point - first_point_det_); + return std::max(0.0, d0); + else if (set_points_ and set_detector_model_) { + return GetDistanceFromStartInBounds(detector_model_->ToGeo(point)); + } else { + throw(std::runtime_error("Detector points not set!")); + } } } // namespace detector diff --git a/projects/detector/public/LeptonInjector/detector/Path.h b/projects/detector/public/LeptonInjector/detector/Path.h index b26a17b63..0ea1f7a2c 100644 --- a/projects/detector/public/LeptonInjector/detector/Path.h +++ b/projects/detector/public/LeptonInjector/detector/Path.h @@ -23,6 +23,7 @@ #include "LeptonInjector/dataclasses/Particle.h" // for Particle #include "LeptonInjector/geometry/Geometry.h" // for Geometry #include "LeptonInjector/detector/DetectorModel.h" // for DetectorModel +#include "LeptonInjector/detector/Coordinates.h" #include "LeptonInjector/math/Vector3D.h" // for Vector3D namespace LI { @@ -33,22 +34,31 @@ class Path { std::shared_ptr detector_model_; bool set_detector_model_ = false; - math::Vector3D first_point_; - math::Vector3D last_point_; - math::Vector3D direction_; + GeometryPosition first_point_; + GeometryPosition last_point_; + GeometryDirection direction_; double distance_ = 0; bool set_points_ = false; + DetectorPosition first_point_det_; + DetectorPosition last_point_det_; + DetectorDirection direction_det_; + bool set_det_points_ = false; + double column_depth_cached_; bool set_column_depth_ = false; geometry::Geometry::IntersectionList intersections_; bool set_intersections_ = false; + + void UpdatePoints(); public: Path(); Path(std::shared_ptr detector_model); - Path(std::shared_ptr detector_model, math::Vector3D const & first_point, math::Vector3D const & last_point); - Path(std::shared_ptr detector_model, math::Vector3D const & first_point, math::Vector3D const & direction, double distance); + Path(std::shared_ptr detector_model, GeometryPosition const & first_point, GeometryPosition const & last_point); + Path(std::shared_ptr detector_model, GeometryPosition const & first_point, GeometryDirection const & direction, double distance); + Path(std::shared_ptr detector_model, DetectorPosition const & first_point, DetectorPosition const & last_point); + Path(std::shared_ptr detector_model, DetectorPosition const & first_point, DetectorDirection const & direction, double distance); template void serialize(Archive & archive, std::uint32_t const version) { @@ -70,17 +80,22 @@ class Path { bool HasColumnDepth(); std::shared_ptr GetDetectorModel(); - math::Vector3D const & GetFirstPoint(); - math::Vector3D const & GetLastPoint(); - math::Vector3D const & GetDirection(); + DetectorPosition const & GetFirstPoint(); + DetectorPosition const & GetLastPoint(); + DetectorPosition const & GetDirection(); + GeometryPosition const & GetGeoFirstPoint(); + GeometryPosition const & GetGeoLastPoint(); + GeometryPosition const & GetGeoDirection(); double GetDistance(); geometry::Geometry::IntersectionList const & GetIntersections(); void SetDetectorModel(std::shared_ptr detector_model); void EnsureDetectorModel(); - void SetPoints(math::Vector3D first_point, math::Vector3D last_point); - void SetPointsWithRay(math::Vector3D first_point, math::Vector3D direction, double distance); + void SetPoints(GeometryPosition first_point, GeometryPosition last_point); + void SetPointsWithRay(GeometryPosition first_point, GeometryDirection direction, double distance); + void SetPoints(DetectorPosition first_point, DetectorPosition last_point); + void SetPointsWithRay(DetectorPosition first_point, DetectorDirection direction, double distance); void EnsurePoints(); void SetIntersections(geometry::Geometry::IntersectionList const & intersections); @@ -224,8 +239,10 @@ class Path { double const & total_decay_length); // - bool IsWithinBounds(math::Vector3D point); - double GetDistanceFromStartInBounds(math::Vector3D point); + bool IsWithinBounds(GeometryPosition point); + double GetDistanceFromStartInBounds(GeometryPosition point); + bool IsWithinBounds(DetectorPosition point); + double GetDistanceFromStartInBounds(DetectorPosition point); }; } // namespace detector From 17941fabff898369d72c7d1e9091d7641e017fee Mon Sep 17 00:00:00 2001 From: Austin Schneider Date: Mon, 15 Jan 2024 18:57:00 -0700 Subject: [PATCH 06/42] Fix compilation issues in DetectorModel --- projects/detector/private/DetectorModel.cxx | 28 +++++++++---------- .../LeptonInjector/detector/Coordinates.h | 8 +++--- 2 files changed, 18 insertions(+), 18 deletions(-) diff --git a/projects/detector/private/DetectorModel.cxx b/projects/detector/private/DetectorModel.cxx index 3849cbb4d..2a130b76a 100644 --- a/projects/detector/private/DetectorModel.cxx +++ b/projects/detector/private/DetectorModel.cxx @@ -135,7 +135,7 @@ DetectorDirection DetectorModel::ToDet(GeometryDirection && dir) const { } GeometryPosition DetectorModel::ToGeo(DetectorPosition const & pos) const { - return GeometryPosition(detector_rotation_..rotate(pos, true) + detector_origin_); + return GeometryPosition(detector_rotation_.rotate(pos, true) + detector_origin_); } GeometryDirection DetectorModel::ToGeo(DetectorDirection const & dir) const { @@ -209,7 +209,7 @@ void DetectorModel::SetSectors(std::vector const & sectors) { } GeometryPosition DetectorModel::GetDetectorOrigin() const { - return detector_origin_; + return GeometryPosition(detector_origin_); } void DetectorModel::SetDetectorOrigin(GeometryPosition const & detector_origin) { @@ -518,7 +518,7 @@ double DetectorModel::GetMassDensity(Geometry::IntersectionList const & intersec double DetectorModel::GetMassDensity(GeometryPosition const & p0) const { Vector3D direction(1,0,0); // Any direction will work for determining the sector heirarchy - Geometry::IntersectionList intersections = GetIntersections(p0, direction); + Geometry::IntersectionList intersections = GetIntersections(p0, GeometryDirection(direction)); return GetMassDensity(intersections, p0); } @@ -565,7 +565,7 @@ double DetectorModel::GetParticleDensity(Geometry::IntersectionList const & inte double DetectorModel::GetParticleDensity(GeometryPosition const & p0, LI::dataclasses::Particle::ParticleType target) const { Vector3D direction(1,0,0); // Any direction will work for determining the sector heirarchy - Geometry::IntersectionList intersections = GetIntersections(p0, direction); + Geometry::IntersectionList intersections = GetIntersections(p0, GeometryDirection(direction)); return GetParticleDensity(intersections, p0, target); } @@ -632,7 +632,7 @@ double DetectorModel::GetInteractionDensity(GeometryPosition const & p0, std::vector const & total_cross_sections, double const & total_decay_length) const { Vector3D direction(1,0,0); // Any direction will work for determining the sector heirarchy - Geometry::IntersectionList intersections = GetIntersections(p0, direction); + Geometry::IntersectionList intersections = GetIntersections(p0, GeometryDirection(direction)); return GetInteractionDensity(intersections, p0, targets, total_cross_sections, total_decay_length); } @@ -692,7 +692,7 @@ double DetectorModel::GetColumnDepthInCGS(GeometryPosition const & p0, GeometryP } direction.normalize(); - Geometry::IntersectionList intersections = GetIntersections(p0, direction); + Geometry::IntersectionList intersections = GetIntersections(p0, GeometryDirection(direction)); return GetColumnDepthInCGS(intersections, p0, p1); } @@ -809,7 +809,7 @@ double DetectorModel::GetMassDensity(Geometry::IntersectionList const & intersec double DetectorModel::GetMassDensity(GeometryPosition const & p0, std::set targets) const { Vector3D direction(1,0,0); // Any direction will work for determining the sector heirarchy - Geometry::IntersectionList intersections = GetIntersections(p0, direction); + Geometry::IntersectionList intersections = GetIntersections(p0, GeometryDirection(direction)); return GetMassDensity(intersections, p0, targets); } @@ -861,7 +861,7 @@ std::vector DetectorModel::GetParticleDensity(Geometry::IntersectionList std::vector DetectorModel::GetParticleDensity(GeometryPosition const & p0, std::set targets) const { Vector3D direction(1,0,0); // Any direction will work for determining the sector heirarchy - Geometry::IntersectionList intersections = GetIntersections(p0, direction); + Geometry::IntersectionList intersections = GetIntersections(p0, GeometryDirection(direction)); return GetParticleDensity(intersections, p0, targets); } @@ -993,7 +993,7 @@ double DetectorModel::GetInteractionDepthInCGS(GeometryPosition const & p0, Geom } direction.normalize(); - Geometry::IntersectionList intersections = GetIntersections(p0, direction); + Geometry::IntersectionList intersections = GetIntersections(p0, GeometryDirection(direction)); return GetInteractionDepthInCGS(intersections, p0, p1, targets, total_cross_sections, total_decay_length); } @@ -1030,11 +1030,11 @@ DetectorSector DetectorModel::GetContainingSector(Geometry::IntersectionList con DetectorSector DetectorModel::GetContainingSector(GeometryPosition const & p0) const { Vector3D direction(0, 0, 1); - Geometry::IntersectionList intersections = GetIntersections(p0, direction); + Geometry::IntersectionList intersections = GetIntersections(p0, GeometryDirection(direction)); return GetContainingSector(intersections, p0); } -Geometry::IntersectionList DetectorModel::GetIntersections(GeometryPosition const & p0, GeometryPosition const & direction) const { +Geometry::IntersectionList DetectorModel::GetIntersections(GeometryPosition const & p0, GeometryDirection const & direction) const { Geometry::IntersectionList intersections; intersections.position = p0; intersections.direction = direction; @@ -1116,12 +1116,12 @@ Geometry::IntersectionList DetectorModel::GetOuterBounds(GeometryPosition const } std::set DetectorModel::GetAvailableTargets(GeometryPosition const & vertex) const { - Geometry::IntersectionList intersections = GetIntersections(vertex, math::Vector3D(0,0,1)); + Geometry::IntersectionList intersections = GetIntersections(vertex, GeometryDirection(math::Vector3D(0,0,1))); return GetAvailableTargets(intersections, vertex); } std::set DetectorModel::GetAvailableTargets(geometry::Geometry::IntersectionList const & intersections, GeometryPosition const & vertex) const { - int matID = GetContainingSector(intersections, Vector3D(vertex[0],vertex[1],vertex[2])).material_id; + int matID = GetContainingSector(intersections, vertex).material_id; std::vector particles = materials_.GetMaterialConstituents(matID); return std::set(particles.begin(), particles.end()); } @@ -1417,7 +1417,7 @@ DetectorSector DetectorModel::GetContainingSector(DetectorPosition const & p0) c return GetContainingSector(ToGeo(p0)); } -Geometry::IntersectionList DetectorModel::GetIntersections(DetectorPosition const & p0, DetectorPosition const & direction) const { +Geometry::IntersectionList DetectorModel::GetIntersections(DetectorPosition const & p0, DetectorDirection const & direction) const { return GetIntersections(ToGeo(p0), ToGeo(direction)); } diff --git a/projects/detector/public/LeptonInjector/detector/Coordinates.h b/projects/detector/public/LeptonInjector/detector/Coordinates.h index 81d212245..bab3578de 100644 --- a/projects/detector/public/LeptonInjector/detector/Coordinates.h +++ b/projects/detector/public/LeptonInjector/detector/Coordinates.h @@ -9,10 +9,10 @@ namespace LI { namespace detector { -using GeometryPosition = fluent::NamedType; -using GeometryDirection = fluent::NamedType; -using DetectorPosition = fluent::NamedType; -using DetectorDirection = fluent::NamedType; +using GeometryPosition = fluent::NamedType; +using GeometryDirection = fluent::NamedType; +using DetectorPosition = fluent::NamedType; +using DetectorDirection = fluent::NamedType; } // namespace detector } // namespace LI From 90cc1b21c8c1813d70e8b4a9d46ed832a01d4cc1 Mon Sep 17 00:00:00 2001 From: Austin Schneider Date: Mon, 15 Jan 2024 19:11:20 -0700 Subject: [PATCH 07/42] Fix compilation issues in Path --- projects/detector/private/Path.cxx | 80 +++++++++++-------- .../public/LeptonInjector/detector/Path.h | 4 +- 2 files changed, 47 insertions(+), 37 deletions(-) diff --git a/projects/detector/private/Path.cxx b/projects/detector/private/Path.cxx index 9c05d993d..b1238e814 100644 --- a/projects/detector/private/Path.cxx +++ b/projects/detector/private/Path.cxx @@ -36,12 +36,22 @@ Path::Path(std::shared_ptr detector_model) { SetDetectorModel(detector_model); } -Path::Path(std::shared_ptr detector_model, math::Vector3D const & first_point, math::Vector3D const & last_point) { +Path::Path(std::shared_ptr detector_model, GeometryPosition const & first_point, GeometryPosition const & last_point) { SetDetectorModel(detector_model); SetPoints(first_point, last_point); } -Path::Path(std::shared_ptr detector_model, math::Vector3D const & first_point, math::Vector3D const & direction, double distance) { +Path::Path(std::shared_ptr detector_model, GeometryPosition const & first_point, GeometryDirection const & direction, double distance) { + SetDetectorModel(detector_model); + SetPointsWithRay(first_point, direction, distance); +} + +Path::Path(std::shared_ptr detector_model, DetectorPosition const & first_point, DetectorPosition const & last_point) { + SetDetectorModel(detector_model); + SetPoints(first_point, last_point); +} + +Path::Path(std::shared_ptr detector_model, DetectorPosition const & first_point, DetectorDirection const & direction, double distance) { SetDetectorModel(detector_model); SetPointsWithRay(first_point, direction, distance); } @@ -76,7 +86,7 @@ DetectorPosition const & Path::GetLastPoint() { return last_point_det_; } -DetectorPosition const & Path::GetDirection() { +DetectorDirection const & Path::GetDirection() { UpdatePoints(); return direction_det_; } @@ -91,13 +101,13 @@ GeometryPosition const & Path::GetGeoLastPoint() { return last_point_; } -GeometryPosition const & Path::GetGeoDirection() { +GeometryDirection const & Path::GetGeoDirection() { UpdatePoints(); return direction_; } double Path::GetDistance() { - return distance__; + return distance_; } geometry::Geometry::IntersectionList const & Path::GetIntersections() { @@ -122,9 +132,9 @@ void Path::EnsureDetectorModel() { void Path::SetPoints(GeometryPosition first_point, GeometryPosition last_point) { first_point_ = first_point; last_point_ = last_point; - direction_ = last_point_ - first_point_; - distance_ = direction_.magnitude(); - direction_.normalize(); + direction_ = GeometryDirection(last_point_ - first_point_); + distance_ = direction_->magnitude(); + direction_->normalize(); set_points_ = true; set_det_points_ = false; set_intersections_ = false; @@ -135,9 +145,9 @@ void Path::SetPoints(GeometryPosition first_point, GeometryPosition last_point) void Path::SetPoints(DetectorPosition first_point, DetectorPosition last_point) { first_point_det_ = first_point; last_point_det_ = last_point; - direction_det_ = last_point_det_ - first_point_det_; - distance_ = direction_det_.magnitude(); - direction_det_.normalize(); + direction_det_ = DetectorDirection(last_point_det_ - first_point_det_); + distance_ = direction_det_->magnitude(); + direction_det_->normalize(); set_points_ = false; set_det_points_ = true; set_intersections_ = false; @@ -145,14 +155,14 @@ void Path::SetPoints(DetectorPosition first_point, DetectorPosition last_point) UpdatePoints(); } -void Path::SetPointsWithRay(GeometryPosition first_point, GeometryPosition direction, double distance) { +void Path::SetPointsWithRay(GeometryPosition first_point, GeometryDirection direction, double distance) { first_point_ = first_point; direction_ = direction; - direction_.normalize(); + direction_->normalize(); //double dif = std::abs(direction_.magnitude() - direction.magnitude()) / std::max(direction_.magnitude(), direction.magnitude()); //if(not std::isnan(dif)) assert(dif < 1e-12); distance_ = distance; - last_point_ = first_point + direction * distance; + last_point_ = first_point + GeometryPosition(direction * distance); set_points_ = true; set_det_points_ = false; set_intersections_ = false; @@ -160,14 +170,14 @@ void Path::SetPointsWithRay(GeometryPosition first_point, GeometryPosition direc UpdatePoints(); } -void Path::SetPointsWithRay(DetectorPosition first_point, DetectorPosition direction, double distance) { +void Path::SetPointsWithRay(DetectorPosition first_point, DetectorDirection direction, double distance) { first_point_det_ = first_point; direction_det_ = direction; - direction_det_.normalize(); + direction_det_->normalize(); //double dif = std::abs(direction_.magnitude() - direction.magnitude()) / std::max(direction_.magnitude(), direction.magnitude()); //if(not std::isnan(dif)) assert(dif < 1e-12); - distance_det_ = distance; - last_point_det_ = first_point + direction * distance; + distance_ = distance; + last_point_det_ = first_point + DetectorPosition(direction * distance); set_points_ = false; set_det_points_ = true; set_intersections_ = false; @@ -219,13 +229,13 @@ void Path::ClipToOuterBounds() { bool clip_1 = (p1 - last_point_) * direction_ < 0; bool clip = clip_0 or clip_1; if(clip_0) { - first_point_ = p0; + first_point_ = GeometryPosition(p0); } if(clip_1) { - last_point_ = p1; + last_point_ = GeometryPosition(p1); } if(clip) { - distance_ = (last_point_ - first_point_).magnitude(); + distance_ = (last_point_ - first_point_)->magnitude(); set_column_depth_ = false; } set_det_points_ = false; @@ -487,7 +497,7 @@ double Path::GetColumnDepthFromStartInBounds(double distance) { } EnsureIntersections(); EnsurePoints(); - return detector_model_->GetColumnDepthInCGS(intersections_, first_point_, first_point_ + direction_ * distance); + return detector_model_->GetColumnDepthInCGS(intersections_, first_point_, GeometryPosition(first_point_ + direction_ * distance)); } double Path::GetColumnDepthFromEndInBounds(double distance) { @@ -498,31 +508,31 @@ double Path::GetColumnDepthFromEndInBounds(double distance) { } EnsureIntersections(); EnsurePoints(); - return detector_model_->GetColumnDepthInCGS(intersections_, last_point_, last_point_ + direction_ * -distance); + return detector_model_->GetColumnDepthInCGS(intersections_, last_point_, GeometryPosition(last_point_ + direction_ * -distance)); } double Path::GetColumnDepthFromStartAlongPath(double distance) { EnsureIntersections(); EnsurePoints(); - return std::copysign(detector_model_->GetColumnDepthInCGS(intersections_, first_point_, first_point_ + direction_ * distance), distance); + return std::copysign(detector_model_->GetColumnDepthInCGS(intersections_, first_point_, GeometryPosition(first_point_ + direction_ * distance)), distance); } double Path::GetColumnDepthFromEndAlongPath(double distance) { EnsureIntersections(); EnsurePoints(); - return std::copysign(detector_model_->GetColumnDepthInCGS(intersections_, last_point_, last_point_ + direction_ * distance), distance); + return std::copysign(detector_model_->GetColumnDepthInCGS(intersections_, last_point_, GeometryPosition(last_point_ + direction_ * distance)), distance); } double Path::GetColumnDepthFromStartInReverse(double distance) { EnsureIntersections(); EnsurePoints(); - return std::copysign(detector_model_->GetColumnDepthInCGS(intersections_, first_point_, first_point_ + direction_ * -distance), distance); + return std::copysign(detector_model_->GetColumnDepthInCGS(intersections_, first_point_, GeometryPosition(first_point_ + direction_ * -distance)), distance); } double Path::GetColumnDepthFromEndInReverse(double distance) { EnsureIntersections(); EnsurePoints(); - return std::copysign(detector_model_->GetColumnDepthInCGS(intersections_, last_point_, last_point_ + direction_ * -distance), distance); + return std::copysign(detector_model_->GetColumnDepthInCGS(intersections_, last_point_, GeometryPosition(last_point_ + direction_ * -distance)), distance); } @@ -540,7 +550,7 @@ double Path::GetInteractionDepthFromStartInBounds(double distance, } EnsureIntersections(); EnsurePoints(); - return detector_model_->GetInteractionDepthInCGS(intersections_, first_point_, first_point_ + direction_ * distance, targets, total_cross_sections, total_decay_length); + return detector_model_->GetInteractionDepthInCGS(intersections_, first_point_, GeometryPosition(first_point_ + direction_ * distance), targets, total_cross_sections, total_decay_length); } double Path::GetInteractionDepthFromEndInBounds(double distance, @@ -554,7 +564,7 @@ double Path::GetInteractionDepthFromEndInBounds(double distance, } EnsureIntersections(); EnsurePoints(); - return detector_model_->GetInteractionDepthInCGS(intersections_, last_point_, last_point_ + direction_ * -distance, targets, total_cross_sections, total_decay_length); + return detector_model_->GetInteractionDepthInCGS(intersections_, last_point_, GeometryPosition(last_point_ + direction_ * -distance), targets, total_cross_sections, total_decay_length); } double Path::GetInteractionDepthFromStartAlongPath(double distance, @@ -563,7 +573,7 @@ double Path::GetInteractionDepthFromStartAlongPath(double distance, double const & total_decay_length) { EnsureIntersections(); EnsurePoints(); - return std::copysign(detector_model_->GetInteractionDepthInCGS(intersections_, first_point_, first_point_ + direction_ * distance, targets, total_cross_sections, total_decay_length), distance); + return std::copysign(detector_model_->GetInteractionDepthInCGS(intersections_, first_point_, GeometryPosition(first_point_ + direction_ * distance), targets, total_cross_sections, total_decay_length), distance); } double Path::GetInteractionDepthFromEndAlongPath(double distance, @@ -572,7 +582,7 @@ double Path::GetInteractionDepthFromEndAlongPath(double distance, double const & total_decay_length) { EnsureIntersections(); EnsurePoints(); - return std::copysign(detector_model_->GetInteractionDepthInCGS(intersections_, last_point_, last_point_ + direction_ * distance, targets, total_cross_sections, total_decay_length), distance); + return std::copysign(detector_model_->GetInteractionDepthInCGS(intersections_, last_point_, GeometryPosition(last_point_ + direction_ * distance), targets, total_cross_sections, total_decay_length), distance); } double Path::GetInteractionDepthFromStartInReverse(double distance, @@ -581,7 +591,7 @@ double Path::GetInteractionDepthFromStartInReverse(double distance, double const & total_decay_length) { EnsureIntersections(); EnsurePoints(); - return std::copysign(detector_model_->GetInteractionDepthInCGS(intersections_, first_point_, first_point_ + direction_ * -distance, targets, total_cross_sections, total_decay_length), distance); + return std::copysign(detector_model_->GetInteractionDepthInCGS(intersections_, first_point_, GeometryPosition(first_point_ + direction_ * -distance), targets, total_cross_sections, total_decay_length), distance); } double Path::GetInteractionDepthFromEndInReverse(double distance, @@ -590,7 +600,7 @@ double Path::GetInteractionDepthFromEndInReverse(double distance, double const & total_decay_length) { EnsureIntersections(); EnsurePoints(); - return std::copysign(detector_model_->GetInteractionDepthInCGS(intersections_, last_point_, last_point_ + direction_ * -distance, targets, total_cross_sections, total_decay_length), distance); + return std::copysign(detector_model_->GetInteractionDepthInCGS(intersections_, last_point_, GeometryPosition(last_point_ + direction_ * -distance), targets, total_cross_sections, total_decay_length), distance); } @@ -760,10 +770,10 @@ bool Path::IsWithinBounds(DetectorPosition point) { double Path::GetDistanceFromStartInBounds(DetectorPosition point) { UpdatePoints(); - if(set_det_points) { + if(set_det_points_) { double d0 = LI::math::scalar_product(direction_det_, point - first_point_det_); return std::max(0.0, d0); - else if (set_points_ and set_detector_model_) { + } else if (set_points_ and set_detector_model_) { return GetDistanceFromStartInBounds(detector_model_->ToGeo(point)); } else { throw(std::runtime_error("Detector points not set!")); diff --git a/projects/detector/public/LeptonInjector/detector/Path.h b/projects/detector/public/LeptonInjector/detector/Path.h index 0ea1f7a2c..8944a3531 100644 --- a/projects/detector/public/LeptonInjector/detector/Path.h +++ b/projects/detector/public/LeptonInjector/detector/Path.h @@ -82,10 +82,10 @@ class Path { std::shared_ptr GetDetectorModel(); DetectorPosition const & GetFirstPoint(); DetectorPosition const & GetLastPoint(); - DetectorPosition const & GetDirection(); + DetectorDirection const & GetDirection(); GeometryPosition const & GetGeoFirstPoint(); GeometryPosition const & GetGeoLastPoint(); - GeometryPosition const & GetGeoDirection(); + GeometryDirection const & GetGeoDirection(); double GetDistance(); geometry::Geometry::IntersectionList const & GetIntersections(); From 26e90690d88635cda6effcc8812ccda2a4b34322 Mon Sep 17 00:00:00 2001 From: Austin Schneider Date: Mon, 15 Jan 2024 23:49:06 -0700 Subject: [PATCH 08/42] Use new path interface for PointSourcePositionDistribution --- .../PointSourcePositionDistribution.cxx | 21 +++++++++++-------- 1 file changed, 12 insertions(+), 9 deletions(-) diff --git a/projects/distributions/private/primary/vertex/PointSourcePositionDistribution.cxx b/projects/distributions/private/primary/vertex/PointSourcePositionDistribution.cxx index c6167ee91..fd38c9e31 100644 --- a/projects/distributions/private/primary/vertex/PointSourcePositionDistribution.cxx +++ b/projects/distributions/private/primary/vertex/PointSourcePositionDistribution.cxx @@ -13,6 +13,7 @@ #include "LeptonInjector/dataclasses/Particle.h" // for Par... #include "LeptonInjector/detector/DetectorModel.h" // for Ear... #include "LeptonInjector/detector/Path.h" // for Path +#include "LeptonInjector/detector/Coordinates.h" #include "LeptonInjector/distributions/Distributions.h" // for Inj... #include "LeptonInjector/math/Vector3D.h" // for Vec... #include "LeptonInjector/utilities/Errors.h" // for Inj... @@ -21,6 +22,9 @@ namespace LI { namespace distributions { +using detector::DetectorPosition; +using detector::DetectorDirection; + namespace { double log_one_minus_exp_of_negative(double x) { if(x < 1e-1) { @@ -47,11 +51,10 @@ LI::math::Vector3D PointSourcePositionDistribution::SamplePosition(std::shared_p LI::math::Vector3D dir(record.primary_momentum[1], record.primary_momentum[2], record.primary_momentum[3]); dir.normalize(); - LI::math::Vector3D endcap_0 = origin; LI::math::Vector3D endcap_1 = origin + max_distance * dir; - LI::detector::Path path(detector_model, detector_model->GetEarthCoordPosFromDetCoordPos(endcap_0), detector_model->GetEarthCoordDirFromDetCoordDir(dir), max_distance); + LI::detector::Path path(detector_model, DetectorPosition(endcap_0), DetectorDirection(dir), max_distance); path.ClipToOuterBounds(); std::set const & possible_targets = interactions->TargetTypes(); @@ -84,7 +87,7 @@ LI::math::Vector3D PointSourcePositionDistribution::SamplePosition(std::shared_p } double dist = path.GetDistanceFromStartAlongPath(traversed_interaction_depth, targets, total_cross_sections, total_decay_length); - LI::math::Vector3D vertex = detector_model->GetDetCoordPosFromEarthCoordPos(path.GetFirstPoint() + dist * path.GetDirection()); + LI::math::Vector3D vertex = path.GetFirstPoint() + dist * path.GetDirection(); return vertex; } @@ -97,10 +100,10 @@ double PointSourcePositionDistribution::GenerationProbability(std::shared_ptr
  • GetEarthCoordPosFromDetCoordPos(endcap_0), detector_model->GetEarthCoordDirFromDetCoordDir(dir), max_distance); + LI::detector::Path path(detector_model, DetectorPosition(endcap_0), DetectorDirection(dir), max_distance); path.ClipToOuterBounds(); - if(not path.IsWithinBounds(detector_model->GetEarthCoordPosFromDetCoordPos(vertex))) + if(not path.IsWithinBounds(DetectorPosition(vertex))) return 0.0; @@ -121,11 +124,11 @@ double PointSourcePositionDistribution::GenerationProbability(std::shared_ptr
  • GetEarthCoordPosFromDetCoordPos(vertex))); + path.SetPointsWithRay(path.GetFirstPoint(), path.GetDirection(), path.GetDistanceFromStartInBounds(DetectorPosition(vertex))); double traversed_interaction_depth = path.GetInteractionDepthInBounds(targets, total_cross_sections, total_decay_length); - double interaction_density = detector_model->GetInteractionDensity(path.GetIntersections(), detector_model->GetEarthCoordPosFromDetCoordPos(vertex), targets, total_cross_sections, total_decay_length); + double interaction_density = detector_model->GetInteractionDensity(path.GetIntersections(), DetectorPosition(vertex), targets, total_cross_sections, total_decay_length); double prob_density; if(total_interaction_depth < 1e-6) { @@ -157,10 +160,10 @@ std::pair PointSourcePositionDistributio LI::math::Vector3D endcap_0 = origin; LI::math::Vector3D endcap_1 = origin + max_distance * dir; - LI::detector::Path path(detector_model, detector_model->GetEarthCoordPosFromDetCoordPos(endcap_0), detector_model->GetEarthCoordDirFromDetCoordDir(dir), max_distance); + LI::detector::Path path(detector_model, DetectorPosition(endcap_0), DetectorDirection(dir), max_distance); path.ClipToOuterBounds(); - if(not path.IsWithinBounds(vertex)) + if(not path.IsWithinBounds(DetectorPosition(vertex))) return std::pair(LI::math::Vector3D(0, 0, 0), LI::math::Vector3D(0, 0, 0)); return std::pair(path.GetFirstPoint(), path.GetLastPoint()); } From df8ea5da9bea090455575c7acd096c666117d091 Mon Sep 17 00:00:00 2001 From: Austin Schneider Date: Mon, 15 Jan 2024 23:50:16 -0700 Subject: [PATCH 09/42] Use new path interface for SecondaryPositionDistribution --- .../vertex/SecondaryPositionDistribution.cxx | 32 +++++++++++-------- 1 file changed, 18 insertions(+), 14 deletions(-) diff --git a/projects/distributions/private/primary/vertex/SecondaryPositionDistribution.cxx b/projects/distributions/private/primary/vertex/SecondaryPositionDistribution.cxx index 70a508e1f..493ceb032 100644 --- a/projects/distributions/private/primary/vertex/SecondaryPositionDistribution.cxx +++ b/projects/distributions/private/primary/vertex/SecondaryPositionDistribution.cxx @@ -14,6 +14,7 @@ #include "LeptonInjector/dataclasses/Particle.h" // for Par... #include "LeptonInjector/detector/DetectorModel.h" // for Ear... #include "LeptonInjector/detector/Path.h" // for Path +#include "LeptonInjector/detector/Coordinates.h" #include "LeptonInjector/distributions/Distributions.h" // for Inj... #include "LeptonInjector/geometry/Geometry.h" // for Geo... #include "LeptonInjector/math/Vector3D.h" // for Vec... @@ -23,6 +24,9 @@ namespace LI { namespace distributions { +using detector::DetectorPosition; +using detector::DetectorDirection; + namespace { double log_one_minus_exp_of_negative(double x) { if(x < 1e-1) { @@ -70,13 +74,13 @@ LI::math::Vector3D SecondaryPositionDistribution::SamplePosition(std::shared_ptr LI::math::Vector3D endcap_0 = LI::math::Vector3D(datum.parent->record.interaction_vertex); LI::math::Vector3D endcap_1 = endcap_0 + max_length * dir; - LI::detector::Path path(detector_model, detector_model->GetEarthCoordPosFromDetCoordPos(endcap_0), detector_model->GetEarthCoordDirFromDetCoordDir(dir), max_length); + LI::detector::Path path(detector_model, DetectorPosition(endcap_0), DetectorDirection(dir), max_length); path.ClipToOuterBounds(); // Check if fiducial volume is provided if(fiducial) { - std::vector fid_intersections = fiducial->Intersections(detector_model->GetEarthCoordPosFromDetCoordPos(endcap_0), - detector_model->GetEarthCoordDirFromDetCoordDir(dir)); + std::vector fid_intersections = fiducial->Intersections(DetectorPosition(endcap_0), + DetectorDirection(dir)); // If the path intersects the fiducial volume, restrict position to that volume if(!fid_intersections.empty()) { // make sure the first intersection happens before the maximum generation length @@ -86,7 +90,7 @@ LI::math::Vector3D SecondaryPositionDistribution::SamplePosition(std::shared_ptr if(update_path) { LI::math::Vector3D first_point = (fid_intersections.front().distance > 0) ? fid_intersections.front().position : endcap_0; LI::math::Vector3D last_point = (fid_intersections.back().distance < max_length) ? fid_intersections.back().position : endcap_1; - path.SetPoints(first_point,last_point); + path.SetPoints(DetectorPosition(first_point), DetectorPosition(last_point)); } } } @@ -121,7 +125,7 @@ LI::math::Vector3D SecondaryPositionDistribution::SamplePosition(std::shared_ptr } double dist = path.GetDistanceFromStartAlongPath(traversed_interaction_depth, targets, total_cross_sections, total_decay_length); - LI::math::Vector3D vertex = detector_model->GetDetCoordPosFromEarthCoordPos(path.GetFirstPoint() + dist * path.GetDirection()); + LI::math::Vector3D vertex = path.GetFirstPoint() + dist * path.GetDirection(); return vertex; } @@ -139,16 +143,16 @@ double SecondaryPositionDistribution::GenerationProbability(std::shared_ptrrecord.interaction_vertex); LI::math::Vector3D endcap_1 = endcap_0 + max_length * dir; - LI::detector::Path path(detector_model, detector_model->GetEarthCoordPosFromDetCoordPos(endcap_0), detector_model->GetEarthCoordDirFromDetCoordDir(dir), max_length); + LI::detector::Path path(detector_model, DetectorPosition(endcap_0), DetectorDirection(dir), max_length); path.ClipToOuterBounds(); - if(not path.IsWithinBounds(detector_model->GetEarthCoordPosFromDetCoordPos(vertex))) + if(not path.IsWithinBounds(DetectorPosition(vertex))) return 0.0; // Check if fiducial volume is provided if(fiducial) { - std::vector fid_intersections = fiducial->Intersections(detector_model->GetEarthCoordPosFromDetCoordPos(endcap_0), - detector_model->GetEarthCoordDirFromDetCoordDir(dir)); + std::vector fid_intersections = fiducial->Intersections(DetectorPosition(endcap_0), + DetectorDirection(dir)); // If the path intersects the fiducial volume, restrict position to that volume if(!fid_intersections.empty()) { // make sure the first intersection happens before the maximum generation length @@ -158,7 +162,7 @@ double SecondaryPositionDistribution::GenerationProbability(std::shared_ptr 0) ? fid_intersections.front().position : endcap_0; LI::math::Vector3D last_point = (fid_intersections.back().distance < max_length) ? fid_intersections.back().position : endcap_1; - path.SetPoints(first_point,last_point); + path.SetPoints(DetectorPosition(first_point), DetectorPosition(last_point)); } } } @@ -180,11 +184,11 @@ double SecondaryPositionDistribution::GenerationProbability(std::shared_ptrGetEarthCoordPosFromDetCoordPos(vertex))); + path.SetPointsWithRay(path.GetFirstPoint(), path.GetDirection(), path.GetDistanceFromStartInBounds(DetectorPosition(vertex))); double traversed_interaction_depth = path.GetInteractionDepthInBounds(targets, total_cross_sections, total_decay_length); - double interaction_density = detector_model->GetInteractionDensity(path.GetIntersections(), detector_model->GetEarthCoordPosFromDetCoordPos(vertex), targets, total_cross_sections, total_decay_length); + double interaction_density = detector_model->GetInteractionDensity(path.GetIntersections(), DetectorPosition(vertex), targets, total_cross_sections, total_decay_length); double prob_density; if(total_interaction_depth < 1e-6) { @@ -227,10 +231,10 @@ std::pair SecondaryPositionDistribution: LI::math::Vector3D endcap_0 = LI::math::Vector3D(datum.parent->record.interaction_vertex); LI::math::Vector3D endcap_1 = endcap_0 + max_length * dir; - LI::detector::Path path(detector_model, detector_model->GetEarthCoordPosFromDetCoordPos(endcap_0), detector_model->GetEarthCoordDirFromDetCoordDir(dir), max_length); + LI::detector::Path path(detector_model, DetectorPosition(endcap_0), DetectorDirection(dir), max_length); path.ClipToOuterBounds(); - if(not path.IsWithinBounds(vertex)) + if(not path.IsWithinBounds(DetectorPosition(vertex))) return std::pair(LI::math::Vector3D(0, 0, 0), LI::math::Vector3D(0, 0, 0)); return std::pair(path.GetFirstPoint(), path.GetLastPoint()); } From 293387a67b650965e943fe7d28241c7d228b9962 Mon Sep 17 00:00:00 2001 From: Austin Schneider Date: Mon, 15 Jan 2024 23:55:29 -0700 Subject: [PATCH 10/42] Use new path interface for DecayRangePositionDistribution --- .../vertex/DecayRangePositionDistribution.cxx | 22 +++++++++---------- 1 file changed, 11 insertions(+), 11 deletions(-) diff --git a/projects/distributions/private/primary/vertex/DecayRangePositionDistribution.cxx b/projects/distributions/private/primary/vertex/DecayRangePositionDistribution.cxx index 932f64356..e09aa0e63 100644 --- a/projects/distributions/private/primary/vertex/DecayRangePositionDistribution.cxx +++ b/projects/distributions/private/primary/vertex/DecayRangePositionDistribution.cxx @@ -9,6 +9,7 @@ #include "LeptonInjector/dataclasses/Particle.h" #include "LeptonInjector/detector/DetectorModel.h" #include "LeptonInjector/detector/Path.h" +#include "LeptonInjector/detector/Coordinates.h" #include "LeptonInjector/distributions/Distributions.h" #include "LeptonInjector/distributions/primary/vertex/DecayRangeFunction.h" #include "LeptonInjector/distributions/primary/vertex/RangeFunction.h" @@ -19,6 +20,9 @@ namespace LI { namespace distributions { +using detector::DetectorPosition; +using detector::DetectorDirection; + //--------------- // class DecayRangePositionDistribution : public VertexPositionDistribution //--------------- @@ -40,7 +44,7 @@ LI::math::Vector3D DecayRangePositionDistribution::SamplePosition(std::shared_pt LI::math::Vector3D endcap_0 = pca - endcap_length * dir; LI::math::Vector3D endcap_1 = pca + endcap_length * dir; - LI::detector::Path path(detector_model, detector_model->GetEarthCoordPosFromDetCoordPos(endcap_0), detector_model->GetEarthCoordDirFromDetCoordDir(dir), endcap_length*2); + LI::detector::Path path(detector_model, DetectorPosition(endcap_0), DetectorDirection(dir), endcap_length*2); path.ExtendFromStartByDistance(decay_length * range_function->Multiplier()); path.ClipToOuterBounds(); @@ -48,7 +52,7 @@ LI::math::Vector3D DecayRangePositionDistribution::SamplePosition(std::shared_pt double total_distance = path.GetDistance(); double dist = -decay_length * log(y * (exp(-total_distance/decay_length) - 1) + 1); - LI::math::Vector3D vertex = detector_model->GetDetCoordPosFromEarthCoordPos(path.GetFirstPoint() + dist * path.GetDirection()); + LI::math::Vector3D vertex = path.GetFirstPoint() + dist * path.GetDirection(); return vertex; } @@ -67,17 +71,15 @@ double DecayRangePositionDistribution::GenerationProbability(std::shared_ptrGetEarthCoordPosFromDetCoordPos(endcap_0), detector_model->GetEarthCoordDirFromDetCoordDir(dir), endcap_length*2); + LI::detector::Path path(detector_model, DetectorPosition(endcap_0), DetectorDirection(dir), endcap_length*2); path.ExtendFromStartByDistance(decay_length * range_function->Multiplier()); path.ClipToOuterBounds(); - LI::math::Vector3D earth_vertex = detector_model->GetEarthCoordPosFromDetCoordPos(vertex); - - if(not path.IsWithinBounds(earth_vertex)) + if(not path.IsWithinBounds(vertex)) return 0.0; double total_distance = path.GetDistance(); - double dist = LI::math::scalar_product(path.GetDirection(), earth_vertex - path.GetFirstPoint()); + double dist = LI::math::scalar_product(path.GetDirection(), vertex - path.GetFirstPoint()); double prob_density = exp(-dist / decay_length) / (decay_length * (1.0 - exp(-total_distance / decay_length))); // m^-1 prob_density /= (M_PI * radius * radius); // (m^-1 * m^-2) -> m^-3 @@ -110,13 +112,11 @@ std::pair DecayRangePositionDistribution LI::math::Vector3D endcap_0 = pca - endcap_length * dir; LI::math::Vector3D endcap_1 = pca + endcap_length * dir; - LI::detector::Path path(detector_model, detector_model->GetEarthCoordPosFromDetCoordPos(endcap_0), detector_model->GetEarthCoordDirFromDetCoordDir(dir), endcap_length*2); + LI::detector::Path path(detector_model, DetectorPosition(endcap_0), DetectorDirection(dir), endcap_length*2); path.ExtendFromStartByDistance(decay_length * range_function->Multiplier()); path.ClipToOuterBounds(); - LI::math::Vector3D earth_vertex = detector_model->GetEarthCoordPosFromDetCoordPos(vertex); - - if(not path.IsWithinBounds(earth_vertex)) + if(not path.IsWithinBounds(DetectorPosition(vertex))) return std::pair(LI::math::Vector3D(0, 0, 0), LI::math::Vector3D(0, 0, 0)); return std::pair(path.GetFirstPoint(), path.GetLastPoint()); From 65ad7aff5ed317c93c74ef8fd14129c174ac65af Mon Sep 17 00:00:00 2001 From: Austin Schneider Date: Mon, 15 Jan 2024 23:58:21 -0700 Subject: [PATCH 11/42] Use new path interface for RangePositionDistribution --- .../vertex/RangePositionDistribution.cxx | 20 +++++++++++-------- 1 file changed, 12 insertions(+), 8 deletions(-) diff --git a/projects/distributions/private/primary/vertex/RangePositionDistribution.cxx b/projects/distributions/private/primary/vertex/RangePositionDistribution.cxx index c9ba0fcc6..6d498b525 100644 --- a/projects/distributions/private/primary/vertex/RangePositionDistribution.cxx +++ b/projects/distributions/private/primary/vertex/RangePositionDistribution.cxx @@ -13,6 +13,7 @@ #include "LeptonInjector/dataclasses/Particle.h" #include "LeptonInjector/detector/DetectorModel.h" #include "LeptonInjector/detector/Path.h" +#include "LeptonInjector/detector/Coordinates.h" #include "LeptonInjector/distributions/Distributions.h" #include "LeptonInjector/distributions/primary/vertex/RangeFunction.h" #include "LeptonInjector/math/Quaternion.h" @@ -23,6 +24,9 @@ namespace LI { namespace distributions { +using detector::DetectorPosition; +using detector::DetectorDirection; + namespace { double log_one_minus_exp_of_negative(double x) { if(x < 1e-1) { @@ -62,7 +66,7 @@ LI::math::Vector3D RangePositionDistribution::SamplePosition(std::shared_ptrGetEarthCoordPosFromDetCoordPos(endcap_0), detector_model->GetEarthCoordDirFromDetCoordDir(dir), endcap_length*2); + LI::detector::Path path(detector_model, DetectorPosition(endcap_0), DetectorDirection(dir), endcap_length*2); path.ExtendFromStartByDistance(lepton_range); path.ClipToOuterBounds(); @@ -96,7 +100,7 @@ LI::math::Vector3D RangePositionDistribution::SamplePosition(std::shared_ptrGetDetCoordPosFromEarthCoordPos(path.GetFirstPoint() + dist * path.GetDirection()); + LI::math::Vector3D vertex = path.GetFirstPoint() + dist * path.GetDirection(); return vertex; } @@ -115,11 +119,11 @@ double RangePositionDistribution::GenerationProbability(std::shared_ptrGetEarthCoordPosFromDetCoordPos(endcap_0), detector_model->GetEarthCoordDirFromDetCoordDir(dir), endcap_length*2); + LI::detector::Path path(detector_model, DetectorPosition(endcap_0), DetectorDirection(dir), endcap_length*2); path.ExtendFromStartByDistance(lepton_range); path.ClipToOuterBounds(); - if(not path.IsWithinBounds(detector_model->GetEarthCoordPosFromDetCoordPos(vertex))) + if(not path.IsWithinBounds(DetectorPosition(vertex))) return 0.0; std::set const & possible_targets = interactions->TargetTypes(); @@ -139,11 +143,11 @@ double RangePositionDistribution::GenerationProbability(std::shared_ptrGetEarthCoordPosFromDetCoordPos(vertex))); + path.SetPointsWithRay(path.GetFirstPoint(), path.GetDirection(), path.GetDistanceFromStartInBounds(DetectorPosition(vertex))); double traversed_interaction_depth = path.GetInteractionDepthInBounds(targets, total_cross_sections, total_decay_length); - double interaction_density = detector_model->GetInteractionDensity(path.GetIntersections(), detector_model->GetEarthCoordPosFromDetCoordPos(vertex), targets, total_cross_sections, total_decay_length); + double interaction_density = detector_model->GetInteractionDensity(path.GetIntersections(), DetectorPosition(vertex), targets, total_cross_sections, total_decay_length); double prob_density; if(total_interaction_depth < 1e-6) { @@ -182,11 +186,11 @@ std::pair RangePositionDistribution::Inj LI::math::Vector3D endcap_0 = pca - endcap_length * dir; LI::math::Vector3D endcap_1 = pca + endcap_length * dir; - LI::detector::Path path(detector_model, detector_model->GetEarthCoordPosFromDetCoordPos(endcap_0), detector_model->GetEarthCoordDirFromDetCoordDir(dir), endcap_length*2); + LI::detector::Path path(detector_model, DetectorPosition(endcap_0), DetectorDirection(dir), endcap_length*2); path.ExtendFromStartByDistance(lepton_range); path.ClipToOuterBounds(); - if(not path.IsWithinBounds(vertex)) + if(not path.IsWithinBounds(DetectorPosition(vertex))) return std::pair(LI::math::Vector3D(0, 0, 0), LI::math::Vector3D(0, 0, 0)); return std::pair(path.GetFirstPoint(), path.GetLastPoint()); } From 808268addd742573123bdf3317a133449d5d2ccc Mon Sep 17 00:00:00 2001 From: Austin Schneider Date: Tue, 16 Jan 2024 00:02:11 -0700 Subject: [PATCH 12/42] Use new path interface for ColumnDepthPositionDistribution --- .../ColumnDepthPositionDistribution.cxx | 20 ++++++++++--------- 1 file changed, 11 insertions(+), 9 deletions(-) diff --git a/projects/distributions/private/primary/vertex/ColumnDepthPositionDistribution.cxx b/projects/distributions/private/primary/vertex/ColumnDepthPositionDistribution.cxx index 1d905a495..04d1fcca5 100644 --- a/projects/distributions/private/primary/vertex/ColumnDepthPositionDistribution.cxx +++ b/projects/distributions/private/primary/vertex/ColumnDepthPositionDistribution.cxx @@ -13,6 +13,7 @@ #include "LeptonInjector/dataclasses/Particle.h" #include "LeptonInjector/detector/DetectorModel.h" #include "LeptonInjector/detector/Path.h" +#include "LeptonInjector/detector/Coordinates.h" #include "LeptonInjector/distributions/Distributions.h" #include "LeptonInjector/distributions/primary/vertex/DepthFunction.h" #include "LeptonInjector/math/Quaternion.h" @@ -23,6 +24,9 @@ namespace LI { namespace distributions { +using detector::DetectorPosition; +using detector::DetectorDirection; + namespace { double log_one_minus_exp_of_negative(double x) { if(x < 1e-1) { @@ -62,7 +66,7 @@ LI::math::Vector3D ColumnDepthPositionDistribution::SamplePosition(std::shared_p LI::math::Vector3D endcap_0 = pca - endcap_length * dir; LI::math::Vector3D endcap_1 = pca + endcap_length * dir; - LI::detector::Path path(detector_model, detector_model->GetEarthCoordPosFromDetCoordPos(endcap_0), detector_model->GetEarthCoordDirFromDetCoordDir(dir), endcap_length*2); + LI::detector::Path path(detector_model, DetectorPosition(endcap_0), DetectorDirection(dir), endcap_length*2); path.ExtendFromStartByColumnDepth(lepton_depth); path.ClipToOuterBounds(); @@ -97,7 +101,7 @@ LI::math::Vector3D ColumnDepthPositionDistribution::SamplePosition(std::shared_p } double dist = path.GetDistanceFromStartAlongPath(traversed_interaction_depth, targets, total_cross_sections, total_decay_length); - LI::math::Vector3D vertex = detector_model->GetDetCoordPosFromEarthCoordPos(path.GetFirstPoint() + dist * path.GetDirection()); + LI::math::Vector3D vertex = path.GetFirstPoint() + dist * path.GetDirection(); return vertex; } @@ -125,13 +129,11 @@ double ColumnDepthPositionDistribution::GenerationProbability(std::shared_ptr
  • GetEarthCoordPosFromDetCoordPos(endcap_0), detector_model->GetEarthCoordDirFromDetCoordDir(dir), endcap_length*2); + LI::detector::Path path(detector_model, DetectorPosition(endcap_0), DetectorDirection(dir), endcap_length*2); path.ExtendFromStartByColumnDepth(lepton_depth); path.ClipToOuterBounds(); - LI::math::Vector3D earth_vertex = detector_model->GetEarthCoordPosFromDetCoordPos(vertex); - - if(not path.IsWithinBounds(earth_vertex)) + if(not path.IsWithinBounds(DetectorPosition(vertex))) return 0.0; std::set const & possible_targets = interactions->TargetTypes(); @@ -152,11 +154,11 @@ double ColumnDepthPositionDistribution::GenerationProbability(std::shared_ptr
  • GetInteractionDensity(path.GetIntersections(), earth_vertex, targets, total_cross_sections, total_decay_length); + double interaction_density = detector_model->GetInteractionDensity(path.GetIntersections(), DetectorPosition(vertex), targets, total_cross_sections, total_decay_length); double prob_density; if(total_interaction_depth < 1e-6) { @@ -194,7 +196,7 @@ std::pair ColumnDepthPositionDistributio LI::math::Vector3D endcap_0 = pca - endcap_length * dir; LI::math::Vector3D endcap_1 = pca + endcap_length * dir; - LI::detector::Path path(detector_model, detector_model->GetEarthCoordPosFromDetCoordPos(endcap_0), detector_model->GetEarthCoordDirFromDetCoordDir(dir), endcap_length*2); + LI::detector::Path path(detector_model, DetectorPosition(endcap_0), DetectorDirection(dir), endcap_length*2); path.ExtendFromStartByColumnDepth(lepton_depth); path.ClipToOuterBounds(); return std::pair(path.GetFirstPoint(), path.GetLastPoint()); From 67dbd1958360f24595b5cad0bee7ca52de54fd7a Mon Sep 17 00:00:00 2001 From: Austin Schneider Date: Tue, 16 Jan 2024 00:03:22 -0700 Subject: [PATCH 13/42] Missing one type tag --- .../private/primary/vertex/DecayRangePositionDistribution.cxx | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/projects/distributions/private/primary/vertex/DecayRangePositionDistribution.cxx b/projects/distributions/private/primary/vertex/DecayRangePositionDistribution.cxx index e09aa0e63..107b41673 100644 --- a/projects/distributions/private/primary/vertex/DecayRangePositionDistribution.cxx +++ b/projects/distributions/private/primary/vertex/DecayRangePositionDistribution.cxx @@ -75,7 +75,7 @@ double DecayRangePositionDistribution::GenerationProbability(std::shared_ptrMultiplier()); path.ClipToOuterBounds(); - if(not path.IsWithinBounds(vertex)) + if(not path.IsWithinBounds(DetectorPosition(vertex))) return 0.0; double total_distance = path.GetDistance(); From 74e6e06d51c6b48afb067975ef80b1d0e1dee220 Mon Sep 17 00:00:00 2001 From: Austin Schneider Date: Tue, 16 Jan 2024 00:14:16 -0700 Subject: [PATCH 14/42] Use new path interface for Injector. Removed SamplePairProduction. --- projects/injection/private/Injector.cxx | 77 ++----------------- .../injection/private/test/Injector_TEST.cxx | 2 +- .../LeptonInjector/injection/Injector.h | 1 - 3 files changed, 8 insertions(+), 72 deletions(-) diff --git a/projects/injection/private/Injector.cxx b/projects/injection/private/Injector.cxx index 95193048e..0fb695dfb 100644 --- a/projects/injection/private/Injector.cxx +++ b/projects/injection/private/Injector.cxx @@ -17,6 +17,7 @@ #include "LeptonInjector/detector/DetectorModel.h" #include "LeptonInjector/detector/MaterialModel.h" #include "LeptonInjector/detector/Path.h" +#include "LeptonInjector/detector/Coordinates.h" #include "LeptonInjector/distributions/Distributions.h" #include "LeptonInjector/distributions/primary/vertex/DecayRangeFunction.h" #include "LeptonInjector/distributions/primary/vertex/VertexPositionDistribution.h" @@ -31,6 +32,9 @@ namespace LI { namespace injection { +using detector::DetectorPosition; +using detector::DetectorDirection; + //--------------- // class Injector //--------------- @@ -144,8 +148,8 @@ void Injector::SampleCrossSection(LI::dataclasses::InteractionRecord & record, s primary_direction.normalize(); - LI::geometry::Geometry::IntersectionList intersections = detector_model->GetIntersections(interaction_vertex, primary_direction); - std::set available_targets = detector_model->GetAvailableTargets(intersections, record.interaction_vertex); + LI::geometry::Geometry::IntersectionList intersections = detector_model->GetIntersections(DetectorPosition(interaction_vertex), DetectorDirection(primary_direction)); + std::set available_targets = detector_model->GetAvailableTargets(intersections, DetectorPosition(record.interaction_vertex)); double total_prob = 0.0; double xsec_prob = 0.0; @@ -160,7 +164,7 @@ void Injector::SampleCrossSection(LI::dataclasses::InteractionRecord & record, s for(auto const target : available_targets) { if(possible_targets.find(target) != possible_targets.end()) { // Get target density - double target_density = detector_model->GetParticleDensity(intersections, interaction_vertex, target); + double target_density = detector_model->GetParticleDensity(intersections, DetectorPosition(interaction_vertex), target); // Loop over cross sections that have this target std::vector> const & target_cross_sections = interactions->GetCrossSectionsForTarget(target); for(auto const & cross_section : target_cross_sections) { @@ -316,73 +320,6 @@ void Injector::SampleNeutrissimoDecay(LI::dataclasses::InteractionRecord const & decay.decay_parameters[4] = b; } -void Injector::SamplePairProduction(LI::dataclasses::DecayRecord const & decay, LI::dataclasses::InteractionRecord & interaction) const { - // function for simulating the pair production of the photon created in HNL decay - // considers the different radiation lengths of materials in the detector - // Nick TODO: comment more - - LI::detector::MaterialModel const & mat_model = detector_model->GetMaterials(); - LI::math::Vector3D decay_vtx(decay.decay_vertex); - unsigned int gamma_index = 0; - LI::math::Vector3D decay_dir(decay.secondary_momenta[gamma_index][1], - decay.secondary_momenta[gamma_index][2], - decay.secondary_momenta[gamma_index][3]); - - interaction.signature.primary_type = decay.signature.secondary_types[gamma_index]; - interaction.primary_mass = decay.secondary_masses[gamma_index]; - interaction.primary_momentum = decay.secondary_momenta[gamma_index]; - interaction.primary_helicity = decay.secondary_helicity[gamma_index]; - - decay_dir.normalize(); - - LI::detector::Path path(detector_model, decay_vtx, decay_dir, 0); - path.ComputeIntersections(); - - std::vector X0; - std::vector P; - std::vector D; - double x0; double p; double density; - D.push_back(0.); - double N = 0; - double lnP_nopp = 0; // for calculating the probability that no pair production occurs - LI::geometry::Geometry::IntersectionList const & ilist = path.GetIntersections(); - LI::math::Vector3D density_point = decay_vtx; - unsigned int i = 0; - for(unsigned int j=0; j < ilist.intersections.size(); ++j) { - auto const & intersection = ilist.intersections[j]; - if(intersection.distance<0 || std::isinf(intersection.distance)) - continue; - D.push_back(intersection.distance); - x0 = (9./7.)*mat_model.GetMaterialRadiationLength(intersection.matID); // in g/cm^2 - density_point += 0.5*(intersection.position - density_point); - density = detector_model->GetMassDensity(density_point); - x0 *= 0.01/density; // in m - X0.push_back(x0); - p = std::exp(-D[i]/x0) - std::exp(-D[i+1]/x0); - P.push_back(p); - N += p; - lnP_nopp += -(D[i+1] - D[i])/x0; - density_point = intersection.position; - ++i; - } - - interaction.interaction_parameters.resize(1); - interaction.interaction_parameters[0] = std::exp(lnP_nopp); - - // sample the PDF by inverting the CDF - double X = random->Uniform(0, 1); - double C = 0; - if(P.size() > 0) { - unsigned int j=0; - for(; j < P.size(); ++j){ - C += P[j]/N; - if(C>X) {C -= P[j]/N; break;} - } - double pairprod_dist = -X0[j]*std::log(-N*(X - C) + std::exp(-D[j]/X0[j])); - interaction.interaction_vertex = LI::math::Vector3D(decay.decay_vertex) + pairprod_dist * decay_dir; - } -} - // Function to sample secondary processes // // Throws exception if no secondary process exists for the given particle diff --git a/projects/injection/private/test/Injector_TEST.cxx b/projects/injection/private/test/Injector_TEST.cxx index c1fcfccbe..620ee5e21 100644 --- a/projects/injection/private/test/Injector_TEST.cxx +++ b/projects/injection/private/test/Injector_TEST.cxx @@ -386,7 +386,7 @@ TEST(Injector, Generation) if(event.signature.target_type != Particle::ParticleType::unknown) { if(miniboone) injector->SampleSecondaryDecay(event, decay, HNL_decay_width, 1, 0, &MiniBooNE_fiducial, 0.1); else injector->SampleSecondaryDecay(event, decay, HNL_decay_width, 1, 0, &MINERvA_fiducial, 1.0); - injector->SamplePairProduction(decay, pair_prod); + //injector->SamplePairProduction(decay, pair_prod); //basic_weight = weighter.EventWeight(event); simplified_weight = weighter.SimplifiedEventWeight(event); interaction_lengths = ComputeInteractionLengths(detector_model, injector->GetInteractions(), injector->InjectionBounds(event), event); diff --git a/projects/injection/public/LeptonInjector/injection/Injector.h b/projects/injection/public/LeptonInjector/injection/Injector.h index 0a10ef402..85372fba9 100644 --- a/projects/injection/public/LeptonInjector/injection/Injector.h +++ b/projects/injection/public/LeptonInjector/injection/Injector.h @@ -81,7 +81,6 @@ friend cereal::access; virtual void SampleCrossSection(LI::dataclasses::InteractionRecord & record, std::shared_ptr interactions) const; virtual void SampleNeutrissimoDecay(LI::dataclasses::InteractionRecord const & interaction, LI::dataclasses::DecayRecord & decay, double width, double alpha_gen, double alpha_phys, LI::geometry::Geometry *fiducial, double buffer) const; - virtual void SamplePairProduction(LI::dataclasses::DecayRecord const & decay, LI::dataclasses::InteractionRecord & pairprod) const; bool SampleSecondaryProcess(unsigned int idx, std::shared_ptr parent, LI::dataclasses::InteractionTreeDatum & datum); From 1ed97c261678c6a0ebd4c805c697a57584313961 Mon Sep 17 00:00:00 2001 From: Austin Schneider Date: Tue, 16 Jan 2024 00:14:41 -0700 Subject: [PATCH 15/42] Use new path interface for WeightingUtils. --- projects/injection/private/WeightingUtils.cxx | 10 +++++++--- 1 file changed, 7 insertions(+), 3 deletions(-) diff --git a/projects/injection/private/WeightingUtils.cxx b/projects/injection/private/WeightingUtils.cxx index 8675969ab..b8faaecba 100644 --- a/projects/injection/private/WeightingUtils.cxx +++ b/projects/injection/private/WeightingUtils.cxx @@ -11,6 +11,7 @@ #include "LeptonInjector/dataclasses/InteractionSignature.h" // for Int... #include "LeptonInjector/dataclasses/Particle.h" // for Par... #include "LeptonInjector/detector/DetectorModel.h" // for Ear... +#include "LeptonInjector/detector/Coordinates.h" #include "LeptonInjector/geometry/Geometry.h" // for Geo... #include "LeptonInjector/math/Vector3D.h" // for Vec... #include "LeptonInjector/utilities/Constants.h" // for cm @@ -18,9 +19,12 @@ namespace LI { namespace injection { +using detector::DetectorPosition; +using detector::DetectorDirection; + double CrossSectionProbability(std::shared_ptr detector_model, std::shared_ptr interactions, LI::dataclasses::InteractionRecord const & record) { std::set const & possible_targets = interactions->TargetTypes(); - std::set available_targets_list = detector_model->GetAvailableTargets(detector_model->GetEarthCoordPosFromDetCoordPos(record.interaction_vertex)); + std::set available_targets_list = detector_model->GetAvailableTargets(DetectorPosition(record.interaction_vertex)); std::set available_targets(available_targets_list.begin(), available_targets_list.end()); LI::math::Vector3D interaction_vertex( @@ -34,7 +38,7 @@ double CrossSectionProbability(std::shared_ptrGetIntersections(detector_model->GetEarthCoordPosFromDetCoordPos(interaction_vertex), detector_model->GetEarthCoordDirFromDetCoordDir(primary_direction)); + LI::geometry::Geometry::IntersectionList intersections = detector_model->GetIntersections(DetectorPosition(interaction_vertex), DetectorDirection(primary_direction)); double total_prob = 0.0; // double selected_prob = 0.0; @@ -58,7 +62,7 @@ double CrossSectionProbability(std::shared_ptrGetParticleDensity(intersections, detector_model->GetEarthCoordPosFromDetCoordPos(interaction_vertex), target); + double target_density = detector_model->GetParticleDensity(intersections, DetectorPosition(interaction_vertex), target); // Loop over cross sections that have this target std::vector> const & target_cross_sections = interactions->GetCrossSectionsForTarget(target); for(auto const & cross_section : target_cross_sections) { From 3fa70d2f12a6db0c1f4dd61dbcecdbf1c507ad9f Mon Sep 17 00:00:00 2001 From: Austin Schneider Date: Tue, 16 Jan 2024 00:19:50 -0700 Subject: [PATCH 16/42] Use new path interface for Weighter. --- projects/injection/private/Weighter.cxx | 24 ++++++++++++++---------- 1 file changed, 14 insertions(+), 10 deletions(-) diff --git a/projects/injection/private/Weighter.cxx b/projects/injection/private/Weighter.cxx index a105172a3..c1519f0e6 100644 --- a/projects/injection/private/Weighter.cxx +++ b/projects/injection/private/Weighter.cxx @@ -17,6 +17,7 @@ #include "LeptonInjector/dataclasses/InteractionSignature.h" #include "LeptonInjector/dataclasses/Particle.h" #include "LeptonInjector/detector/DetectorModel.h" +#include "LeptonInjector/detector/Coordinates.h" #include "LeptonInjector/distributions/Distributions.h" #include "LeptonInjector/distributions/primary/vertex/VertexPositionDistribution.h" #include "LeptonInjector/geometry/Geometry.h" @@ -27,6 +28,9 @@ namespace LI { namespace injection { +using detector::DetectorPosition; +using detector::DetectorDirection; + namespace { template typename std::iterator_traits::value_type accumulate(InIt begin, InIt end) { @@ -99,7 +103,7 @@ double LeptonWeighter::InteractionProbability(std::pairGetIntersections(detector_model->GetEarthCoordPosFromDetCoordPos(interaction_vertex), detector_model->GetEarthCoordDirFromDetCoordDir(primary_direction)); + LI::geometry::Geometry::IntersectionList intersections = detector_model->GetIntersections(DetectorPosition(interaction_vertex), DetectorDirection(primary_direction)); std::map>> const & cross_sections_by_target = interactions->GetCrossSectionsByTarget(); std::vector targets; targets.reserve(cross_sections_by_target.size()); @@ -123,7 +127,7 @@ double LeptonWeighter::InteractionProbability(std::pairGetInteractionDepthInCGS(intersections, bounds.first, bounds.second, targets, total_cross_sections, total_decay_length); + double total_interaction_depth = detector_model->GetInteractionDepthInCGS(intersections, DetectorPosition(bounds.first), DetectorPosition(bounds.second), targets, total_cross_sections, total_decay_length); double interaction_probability; if(total_interaction_depth < 1e-6) { interaction_probability = total_interaction_depth; @@ -150,7 +154,7 @@ double LeptonWeighter::UnnormalizedPositionProbability(std::pairGetIntersections(detector_model->GetEarthCoordPosFromDetCoordPos(interaction_vertex), primary_direction); + LI::geometry::Geometry::IntersectionList intersections = detector_model->GetIntersections(DetectorPosition(interaction_vertex), DetectorDirection(primary_direction)); std::map>> const & cross_sections_by_target = interactions->GetCrossSectionsByTarget(); unsigned int n_targets = cross_sections_by_target.size(); @@ -176,9 +180,9 @@ double LeptonWeighter::UnnormalizedPositionProbability(std::pairGetInteractionDepthInCGS(intersections, bounds.first, bounds.second, targets, total_cross_sections, total_decay_length); - double traversed_interaction_depth = detector_model->GetInteractionDepthInCGS(intersections, bounds.first, detector_model->GetEarthCoordPosFromDetCoordPos(interaction_vertex), targets, total_cross_sections, total_decay_length); - double interaction_density = detector_model->GetInteractionDensity(intersections, detector_model->GetEarthCoordPosFromDetCoordPos(interaction_vertex), targets, total_cross_sections, total_decay_length); + double total_interaction_depth = detector_model->GetInteractionDepthInCGS(intersections, DetectorPosition(bounds.first), DetectorPosition(bounds.second), targets, total_cross_sections, total_decay_length); + double traversed_interaction_depth = detector_model->GetInteractionDepthInCGS(intersections, DetectorPosition(bounds.first), DetectorPosition(interaction_vertex), targets, total_cross_sections, total_decay_length); + double interaction_density = detector_model->GetInteractionDensity(intersections, DetectorPosition(interaction_vertex), targets, total_cross_sections, total_decay_length); double prob_density; if(total_interaction_depth < 1e-6) { @@ -202,7 +206,7 @@ double LeptonWeighter::NormalizedPositionProbability(std::pairGetIntersections(detector_model->GetEarthCoordPosFromDetCoordPos(interaction_vertex), primary_direction); + LI::geometry::Geometry::IntersectionList intersections = detector_model->GetIntersections(DetectorPosition(interaction_vertex), DetectorDirection(primary_direction)); std::map>> const & cross_sections_by_target = interactions->GetCrossSectionsByTarget(); unsigned int n_targets = cross_sections_by_target.size(); @@ -228,9 +232,9 @@ double LeptonWeighter::NormalizedPositionProbability(std::pairGetInteractionDepthInCGS(intersections, bounds.first, bounds.second, targets, total_cross_sections, total_decay_length); - double traversed_interaction_depth = detector_model->GetInteractionDepthInCGS(intersections, bounds.first, detector_model->GetEarthCoordPosFromDetCoordPos(interaction_vertex), targets, total_cross_sections, total_decay_length); - double interaction_density = detector_model->GetInteractionDensity(intersections, detector_model->GetEarthCoordPosFromDetCoordPos(interaction_vertex), targets, total_cross_sections, total_decay_length); + double total_interaction_depth = detector_model->GetInteractionDepthInCGS(intersections, DetectorPosition(bounds.first), DetectorPosition(bounds.second), targets, total_cross_sections, total_decay_length); + double traversed_interaction_depth = detector_model->GetInteractionDepthInCGS(intersections, DetectorPosition(bounds.first), DetectorPosition(interaction_vertex), targets, total_cross_sections, total_decay_length); + double interaction_density = detector_model->GetInteractionDensity(intersections, DetectorPosition(interaction_vertex), targets, total_cross_sections, total_decay_length); double prob_density; if(total_interaction_depth < 1e-6) { From 36f703b9f3597e7a7c8496337b497e5ea0e236e8 Mon Sep 17 00:00:00 2001 From: Austin Schneider Date: Tue, 16 Jan 2024 00:23:26 -0700 Subject: [PATCH 17/42] Use new path interface for TreeWeighter. --- projects/injection/private/TreeWeighter.cxx | 16 ++++++++++------ 1 file changed, 10 insertions(+), 6 deletions(-) diff --git a/projects/injection/private/TreeWeighter.cxx b/projects/injection/private/TreeWeighter.cxx index c208e46fd..5380afb45 100644 --- a/projects/injection/private/TreeWeighter.cxx +++ b/projects/injection/private/TreeWeighter.cxx @@ -13,6 +13,7 @@ #include "LeptonInjector/dataclasses/InteractionRecord.h" // for Int... #include "LeptonInjector/dataclasses/InteractionSignature.h" // for Int... #include "LeptonInjector/detector/DetectorModel.h" // for Ear... +#include "LeptonInjector/detector/Coordinates.h" #include "LeptonInjector/distributions/Distributions.h" // for Inj... #include "LeptonInjector/geometry/Geometry.h" // for Geo... #include "LeptonInjector/injection/Injector.h" // for Inj... @@ -39,6 +40,9 @@ namespace LI { namespace injection { +using detector::DetectorPosition; +using detector::DetectorDirection; + namespace { template typename std::iterator_traits::value_type accumulate(InIt begin, InIt end) { @@ -224,7 +228,7 @@ double LeptonProcessWeighter::InteractionProbability(std::pairGetIntersections(detector_model->GetEarthCoordPosFromDetCoordPos(interaction_vertex), detector_model->GetEarthCoordDirFromDetCoordDir(primary_direction)); + LI::geometry::Geometry::IntersectionList intersections = detector_model->GetIntersections(DetectorPosition(interaction_vertex), DetectorDirection(primary_direction)); std::map>> const & cross_sections_by_target = phys_process->GetInteractions()->GetCrossSectionsByTarget(); std::vector targets; targets.reserve(cross_sections_by_target.size()); @@ -249,7 +253,7 @@ double LeptonProcessWeighter::InteractionProbability(std::pairGetInteractionDepthInCGS(intersections, bounds.first, bounds.second, targets, total_cross_sections, total_decay_length); + double total_interaction_depth = detector_model->GetInteractionDepthInCGS(intersections, DetectorPosition(bounds.first), DetectorPosition(bounds.second), targets, total_cross_sections, total_decay_length); double interaction_probability; if(total_interaction_depth < 1e-6) { @@ -272,7 +276,7 @@ double LeptonProcessWeighter::NormalizedPositionProbability(std::pairGetIntersections(detector_model->GetEarthCoordPosFromDetCoordPos(interaction_vertex), primary_direction); + LI::geometry::Geometry::IntersectionList intersections = detector_model->GetIntersections(DetectorPosition(interaction_vertex), DetectorDirection(primary_direction)); std::map>> const & cross_sections_by_target = phys_process->GetInteractions()->GetCrossSectionsByTarget(); unsigned int n_targets = cross_sections_by_target.size(); @@ -298,9 +302,9 @@ double LeptonProcessWeighter::NormalizedPositionProbability(std::pairGetInteractionDepthInCGS(intersections, bounds.first, bounds.second, targets, total_cross_sections, total_decay_length); // unitless - double traversed_interaction_depth = detector_model->GetInteractionDepthInCGS(intersections, bounds.first, detector_model->GetEarthCoordPosFromDetCoordPos(interaction_vertex), targets, total_cross_sections, total_decay_length); - double interaction_density = detector_model->GetInteractionDensity(intersections, detector_model->GetEarthCoordPosFromDetCoordPos(interaction_vertex), targets, total_cross_sections, total_decay_length); //units of m^-1 + double total_interaction_depth = detector_model->GetInteractionDepthInCGS(intersections, DetectorPosition(bounds.first), DetectorPosition(bounds.second), targets, total_cross_sections, total_decay_length); // unitless + double traversed_interaction_depth = detector_model->GetInteractionDepthInCGS(intersections, DetectorPosition(bounds.first), DetectorPosition(interaction_vertex), targets, total_cross_sections, total_decay_length); + double interaction_density = detector_model->GetInteractionDensity(intersections, DetectorPosition(interaction_vertex), targets, total_cross_sections, total_decay_length); //units of m^-1 double prob_density; From 2b950b1d1dd2881da726e12891cb7de454e99047 Mon Sep 17 00:00:00 2001 From: Austin Schneider Date: Tue, 16 Jan 2024 11:26:36 -0700 Subject: [PATCH 18/42] Use new path interface for DetectorModel_TEST. --- .../private/test/DetectorModel_TEST.cxx | 55 ++++++++++--------- 1 file changed, 28 insertions(+), 27 deletions(-) diff --git a/projects/detector/private/test/DetectorModel_TEST.cxx b/projects/detector/private/test/DetectorModel_TEST.cxx index 125fc3372..b3dd74845 100644 --- a/projects/detector/private/test/DetectorModel_TEST.cxx +++ b/projects/detector/private/test/DetectorModel_TEST.cxx @@ -15,6 +15,7 @@ #include "LeptonInjector/geometry/Sphere.h" #include "LeptonInjector/detector/DetectorModel.h" #include "LeptonInjector/utilities/Constants.h" +#include "LeptonInjector/detector/Coordinates.h" #include "LeptonInjector/detector/DensityDistribution.h" #include "LeptonInjector/detector/Distribution1D.h" #include "LeptonInjector/detector/Axis1D.h" @@ -171,7 +172,7 @@ TEST(EdgeCases, ColumnDepthWithEqualPoints) unsigned int N_rand = 1000; for(unsigned int i=0; i const * density = dynamic_cast const *>(sector.density.get()); ASSERT_TRUE(density); double rho = density->Evaluate(Vector3D()); - double sum = A.GetColumnDepthInCGS(p0, p1); + double sum = A.GetColumnDepthInCGS(DetectorPosition(p0), DetectorPosition(p1)); EXPECT_DOUBLE_EQ((p1 - p0).magnitude() * 100 * rho, sum); } } @@ -440,7 +441,7 @@ TEST_F(FakeLegacyDetectorModelTest, LegacyFileConstantGetMassDensityInternal) DensityDistribution1D const * density_dist = dynamic_cast const *>(sector.density.get()); ASSERT_TRUE(density_dist); double rho = density_dist->Evaluate(Vector3D()); - double density = A.GetMassDensity(p0); + double density = A.GetMassDensity(DetectorPosition(p0)); EXPECT_DOUBLE_EQ(density, rho); } } @@ -483,7 +484,7 @@ TEST_F(FakeLegacyDetectorModelTest, LegacyFileConstantIntegralNested) bool in_0 = p0.magnitude() <= sphere_0->GetRadius(); bool in_1 = p1.magnitude() <= sphere_0->GetRadius(); double integral = 0.0; - double sum = A.GetColumnDepthInCGS(p0, p1); + double sum = A.GetColumnDepthInCGS(DetectorPosition(p0), DetectorPosition(p1)); if(in_0 and in_1) { integral = rho_0 * (p1-p0).magnitude(); ASSERT_DOUBLE_EQ(integral * 100, sum); @@ -581,7 +582,7 @@ TEST_F(FakeLegacyDetectorModelTest, LegacyFileConstantGetMassDensityNested) double rho_1 = density_1->Evaluate(Vector3D()); ASSERT_LE(p0.magnitude(), max_radius); bool in_0 = p0.magnitude() <= sphere_0->GetRadius(); - double density = A.GetMassDensity(p0); + double density = A.GetMassDensity(DetectorPosition(p0)); if(in_0) { ASSERT_DOUBLE_EQ(density, rho_0); } @@ -652,7 +653,7 @@ TEST_F(FakeLegacyDetectorModelTest, LegacyFileConstantIntegralIntersecting) double rho_upper = density_1->Evaluate(upper_center); double integral = 0.0; - double sum = A.GetColumnDepthInCGS(p0, p1); + double sum = A.GetColumnDepthInCGS(DetectorPosition(p0), DetectorPosition(p1)); std::vector lower_intersections = sphere_0->Intersections(p0, direction); std::vector upper_intersections = sphere_1->Intersections(p0, direction); @@ -941,7 +942,7 @@ TEST_F(FakeLegacyDetectorModelTest, LegacyFileConstantGetMassDensityIntersecting double rho_lower = density_0->Evaluate(lower_center); double rho_upper = density_1->Evaluate(upper_center); - double density = A.GetMassDensity(p0); + double density = A.GetMassDensity(DetectorPosition(p0)); bool p0_in_lower = (p0-lower_center).magnitude() < radius; bool p0_in_upper = (p0-upper_center).magnitude() < radius; @@ -1021,7 +1022,7 @@ TEST_F(FakeLegacyDetectorModelTest, LegacyFileConstantIntegralHidden) bool p0_in_upper = (p0 - upper_center).magnitude() < radius; bool p1_in_upper = (p1 - upper_center).magnitude() < radius; - double sum = A.GetColumnDepthInCGS(p0, p1); + double sum = A.GetColumnDepthInCGS(DetectorPosition(p0), DetectorPosition(p1)); double integral = 0.0; if(p0_in_upper) { if(p1_in_upper) { @@ -1117,7 +1118,7 @@ TEST_F(FakeLegacyDetectorModelTest, LegacyFileConstantGetMassDensityHidden) bool p0_in_upper = (p0 - upper_center).magnitude() < radius; - double density = A.GetMassDensity(p0); + double density = A.GetMassDensity(DetectorPosition(p0)); if(p0_in_upper) { EXPECT_DOUBLE_EQ(density, rho_upper); } else { @@ -1152,12 +1153,12 @@ TEST_F(FakeLegacyDetectorModelTest, LegacyFileConstantInverseIntegralInternal) Vector3D p0 = RandomVector(max_radius, min_radius); Vector3D p1 = RandomVector(max_radius, min_radius); - double sum = A.GetColumnDepthInCGS(p0, p1); + double sum = A.GetColumnDepthInCGS(DetectorPosition(p0), DetectorPosition(p1)); Vector3D direction = p1 - p0; double distance = direction.magnitude(); direction.normalize(); - double found_distance = A.DistanceForColumnDepthFromPoint(p0, direction, sum); + double found_distance = A.DistanceForColumnDepthFromPoint(DetectorPosition(p0), DetectorDirection(direction), sum); EXPECT_NEAR(distance, found_distance, std::max(std::abs(distance), std::abs(found_distance))*1e-12); } @@ -1228,9 +1229,9 @@ TEST_F(FakeLegacyDetectorModelTest, LegacyFileConstantInverseIntegralHidden) ASSERT_TRUE(density_0); ASSERT_TRUE(density_1); - double integral = A.GetColumnDepthInCGS(p0, p1); + double integral = A.GetColumnDepthInCGS(DetectorPosition(p0), DetectorPosition(p1)); - double found_distance = A.DistanceForColumnDepthFromPoint(p0, direction, integral); + double found_distance = A.DistanceForColumnDepthFromPoint(DetectorPosition(p0), DetectorDirection(direction), integral); EXPECT_NEAR(distance, found_distance, std::max(std::abs(distance), std::abs(found_distance))*1e-12); } } @@ -1301,9 +1302,9 @@ TEST_F(FakeLegacyDetectorModelTest, LegacyFileConstantInverseIntegralIntersectin ASSERT_TRUE(density_0); ASSERT_TRUE(density_1); - double integral = A.GetColumnDepthInCGS(p0, p1); + double integral = A.GetColumnDepthInCGS(DetectorPosition(p0), DetectorPosition(p1)); - double found_distance = A.DistanceForColumnDepthFromPoint(p0, direction, integral); + double found_distance = A.DistanceForColumnDepthFromPoint(DetectorPosition(p0), DetectorDirection(direction), integral); EXPECT_NEAR(distance, found_distance, std::max(std::abs(distance), std::abs(found_distance))*1e-12); } } @@ -1362,9 +1363,9 @@ TEST_F(FakeLegacyDetectorModelTest, LegacyFileConstantInverseIntegralNested) ASSERT_TRUE(density_1); ASSERT_LE(p0.magnitude(), max_radius); ASSERT_LE(p1.magnitude(), max_radius); - double integral = A.GetColumnDepthInCGS(p0, p1); + double integral = A.GetColumnDepthInCGS(DetectorPosition(p0), DetectorPosition(p1)); - double found_distance = A.DistanceForColumnDepthFromPoint(p0, direction, integral); + double found_distance = A.DistanceForColumnDepthFromPoint(DetectorPosition(p0), DetectorDirection(direction), integral); EXPECT_NEAR(distance, found_distance, std::max(std::abs(distance), std::abs(found_distance))*1e-12); } } From 37567bbac82263d2eb689990f7ae719ec494e8d3 Mon Sep 17 00:00:00 2001 From: Austin Schneider Date: Tue, 16 Jan 2024 12:07:20 -0700 Subject: [PATCH 19/42] Use new path interface for Path_TEST. --- projects/detector/private/test/Path_TEST.cxx | 249 ++++++++++--------- 1 file changed, 125 insertions(+), 124 deletions(-) diff --git a/projects/detector/private/test/Path_TEST.cxx b/projects/detector/private/test/Path_TEST.cxx index d3f9a42bd..565c47971 100644 --- a/projects/detector/private/test/Path_TEST.cxx +++ b/projects/detector/private/test/Path_TEST.cxx @@ -12,6 +12,7 @@ #include "LeptonInjector/geometry/Sphere.h" #include "LeptonInjector/math/Vector3D.h" #include "LeptonInjector/detector/DetectorModel.h" +#include "LeptonInjector/detector/Coordinates.h" #include "LeptonInjector/detector/DensityDistribution.h" #include "LeptonInjector/detector/DensityDistribution1D.h" #include "LeptonInjector/detector/Distribution1D.h" @@ -77,7 +78,7 @@ TEST(PointsConstructor, NoThrow) std::shared_ptr EMp(new DetectorModel()); Vector3D B(1, 2, 3); Vector3D C(4, 6, 8); - EXPECT_NO_THROW(Path A(EMp, B, C)); + EXPECT_NO_THROW(Path A(EMp, DetectorPosition(B), DetectorPosition(C))); } TEST(PointsConstructor, MemberValues) @@ -88,7 +89,7 @@ TEST(PointsConstructor, MemberValues) Vector3D direction = C - B; double distance = direction.magnitude(); direction.normalize(); - Path A(EMp, B, C); + Path A(EMp, DetectorPosition(B), DetectorPosition(C)); EXPECT_EQ(EMp, A.GetDetectorModel()); EXPECT_EQ(B, A.GetFirstPoint()); EXPECT_EQ(C, A.GetLastPoint()); @@ -105,7 +106,7 @@ TEST(RayConstructor, NoThrow) Vector3D direction = C - B; double distance = direction.magnitude(); direction.normalize(); - EXPECT_NO_THROW(Path A(EMp, B, direction, distance)); + EXPECT_NO_THROW(Path A(EMp, DetectorPosition(B), DetectorDirection(direction), distance)); } TEST(RayConstructor, MemberValues) @@ -121,7 +122,7 @@ TEST(RayConstructor, MemberValues) direction.normalize(); direction.normalize(); direction.normalize(); - Path A(EMp, B, direction, distance); + Path A(EMp, DetectorPosition(B), DetectorDirection(direction), distance); EXPECT_EQ(EMp, A.GetDetectorModel()); EXPECT_EQ(B, A.GetFirstPoint()); EXPECT_EQ(B + direction * distance, A.GetLastPoint()); @@ -155,7 +156,7 @@ TEST(Points, SetGet) { double distance = direction.magnitude(); direction.normalize(); Path C; - C.SetPoints(A, B); + C.SetPoints(DetectorPosition(A), DetectorPosition(B)); EXPECT_EQ(A, C.GetFirstPoint()); EXPECT_EQ(B, C.GetLastPoint()); EXPECT_EQ(direction, C.GetDirection()); @@ -174,7 +175,7 @@ TEST(Points, SetGetRay) { direction.normalize(); direction.normalize(); Path A; - A.SetPointsWithRay(B, direction, distance); + A.SetPointsWithRay(DetectorPosition(B), DetectorDirection(direction), distance); EXPECT_EQ(B, A.GetFirstPoint()); EXPECT_EQ(B + direction * distance, A.GetLastPoint()); EXPECT_EQ(direction, A.GetDirection()); @@ -192,7 +193,7 @@ TEST(EnsurePoints, NoThrow) { Vector3D direction = B - A; direction.normalize(); Path C; - C.SetPoints(A, B); + C.SetPoints(DetectorPosition(A), DetectorPosition(B)); EXPECT_NO_THROW(C.EnsurePoints()); } @@ -208,7 +209,7 @@ TEST(EnsurePoints, NoThrowRay) { direction.normalize(); direction.normalize(); Path A; - A.SetPointsWithRay(B, direction, distance); + A.SetPointsWithRay(DetectorPosition(B), DetectorDirection(direction), distance); EXPECT_NO_THROW(A.EnsurePoints()); } @@ -226,11 +227,11 @@ TEST(EnsureIntersections, Throw) { double distance = direction.magnitude(); direction.normalize(); A = Path(); - A.SetPoints(B, C); + A.SetPoints(DetectorPosition(B), DetectorPosition(C)); EXPECT_THROW(A.EnsureIntersections(), std::runtime_error); A = Path(); - A.SetPointsWithRay(B, direction, distance); + A.SetPointsWithRay(DetectorPosition(B), DetectorDirection(direction), distance); EXPECT_THROW(A.EnsureIntersections(), std::runtime_error); } @@ -246,26 +247,26 @@ TEST(EnsureIntersections, NoThrow) { A = Path(); A.SetDetectorModel(EMp); - A.SetPoints(B, C); + A.SetPoints(DetectorPosition(B), DetectorPosition(C)); EXPECT_NO_THROW(A.EnsureIntersections()); A = Path(); A.SetDetectorModel(EMp); - A.SetPointsWithRay(B, direction, distance); + A.SetPointsWithRay(DetectorPosition(B), DetectorDirection(direction), distance); EXPECT_NO_THROW(A.EnsureIntersections()); A = Path(EMp); - A.SetPoints(B, C); + A.SetPoints(DetectorPosition(B), DetectorPosition(C)); EXPECT_NO_THROW(A.EnsureIntersections()); A = Path(EMp); - A.SetPointsWithRay(B, direction, distance); + A.SetPointsWithRay(DetectorPosition(B), DetectorDirection(direction), distance); EXPECT_NO_THROW(A.EnsureIntersections()); - A = Path(EMp, B, C); + A = Path(EMp, DetectorPosition(B), DetectorPosition(C)); EXPECT_NO_THROW(A.EnsureIntersections()); - A = Path(EMp, B, direction, distance); + A = Path(EMp, DetectorPosition(B), DetectorDirection(direction), distance); EXPECT_NO_THROW(A.EnsureIntersections()); } @@ -279,7 +280,7 @@ TEST(PointManipulation, Flip) { direction.normalize(); reverse.normalize(); - Path A(EMp, B, C); + Path A(EMp, DetectorPosition(B), DetectorPosition(C)); EXPECT_EQ(B, A.GetFirstPoint()); EXPECT_EQ(C, A.GetLastPoint()); EXPECT_EQ(direction, A.GetDirection()); @@ -306,7 +307,7 @@ TEST(PointManipulation, ExtendFromEndByDistance) { double distance = direction.magnitude(); direction.normalize(); - Path A(EMp, B, C); + Path A(EMp, DetectorPosition(B), DetectorPosition(C)); double extra_distance = 101; A.ExtendFromEndByDistance(extra_distance); @@ -315,14 +316,14 @@ TEST(PointManipulation, ExtendFromEndByDistance) { EXPECT_EQ(end, A.GetLastPoint()); EXPECT_EQ(distance + extra_distance, A.GetDistance()); - A = Path(EMp, B, C); + A = Path(EMp, DetectorPosition(B), DetectorPosition(C)); extra_distance = -1; end = C + direction * extra_distance; A.ExtendFromEndByDistance(extra_distance); EXPECT_EQ(end, A.GetLastPoint()); EXPECT_EQ(distance + extra_distance, A.GetDistance()); - A = Path(EMp, B, C); + A = Path(EMp, DetectorPosition(B), DetectorPosition(C)); extra_distance = -101; end = C + direction * extra_distance; A.ExtendFromEndByDistance(extra_distance); @@ -362,44 +363,44 @@ TEST_F(FakeLegacyDetectorModelTest, ExtendFromEndByColumnDepth) { direction = inner_p1 - inner_p0; distance = direction.magnitude(); direction.normalize(); - Path P(A, inner_p0, inner_p1); + Path P(A, DetectorPosition(inner_p0), DetectorPosition(inner_p1)); DensityDistribution1D const * density = dynamic_cast const *>(sector.density.get()); ASSERT_TRUE(density); double rho = density->Evaluate(Vector3D()); double sum = P.GetColumnDepthInBounds(); ASSERT_DOUBLE_EQ((inner_p1 - inner_p0).magnitude() * rho * 100, sum); - P = Path(A, p0, p1); + P = Path(A, DetectorPosition(p0), DetectorPosition(p1)); sum = P.GetColumnDepthInBounds(); ASSERT_DOUBLE_EQ((p1 - p0).magnitude() * rho * 100, sum); - P = Path(A, inner_p0, inner_p1); + P = Path(A, DetectorPosition(inner_p0), DetectorPosition(inner_p1)); P.EnsureIntersections(); double extra_distance = distance/3.0; double extra_column_depth = extra_distance * rho * 100; P.ExtendFromEndByColumnDepth(extra_column_depth); - extra_distance = A->DistanceForColumnDepthFromPoint(inner_p1, direction, extra_column_depth); + extra_distance = A->DistanceForColumnDepthFromPoint(DetectorPosition(inner_p1), DetectorDirection(direction), extra_column_depth); EXPECT_DOUBLE_EQ(distance/3.0, extra_distance); Vector3D end = inner_p1 + direction * extra_distance; EXPECT_EQ(end, P.GetLastPoint()); EXPECT_EQ(distance + extra_distance, P.GetDistance()); - P = Path(A, inner_p0, inner_p1); + P = Path(A, DetectorPosition(inner_p0), DetectorPosition(inner_p1)); P.EnsureIntersections(); extra_distance = -distance/3.0; extra_column_depth = extra_distance * rho * 100; P.ExtendFromEndByColumnDepth(extra_column_depth); - extra_distance = A->DistanceForColumnDepthFromPoint(inner_p1, direction, extra_column_depth); + extra_distance = A->DistanceForColumnDepthFromPoint(DetectorPosition(inner_p1), DetectorDirection(direction), extra_column_depth); EXPECT_DOUBLE_EQ(-distance/3.0, extra_distance); end = inner_p1 + direction * extra_distance; EXPECT_EQ(end, P.GetLastPoint()); EXPECT_EQ(distance + extra_distance, P.GetDistance()); - P = Path(A, inner_p0, inner_p1); + P = Path(A, DetectorPosition(inner_p0), DetectorPosition(inner_p1)); P.EnsureIntersections(); extra_distance = -distance*1.5; extra_column_depth = extra_distance * rho * 100; P.ExtendFromEndByColumnDepth(extra_column_depth); - extra_distance = A->DistanceForColumnDepthFromPoint(inner_p1, direction, extra_column_depth); + extra_distance = A->DistanceForColumnDepthFromPoint(DetectorPosition(inner_p1), DetectorDirection(direction), extra_column_depth); EXPECT_DOUBLE_EQ(-distance*1.5, extra_distance); end = inner_p1 + direction * extra_distance; EXPECT_EQ(inner_p0, P.GetLastPoint()); @@ -417,7 +418,7 @@ TEST(PointManipulation, ExtendFromStartByDistance) { direction.normalize(); reverse.normalize(); - Path A(EMp, B, C); + Path A(EMp, DetectorPosition(B), DetectorPosition(C)); double extra_distance = 101; A.ExtendFromStartByDistance(extra_distance); @@ -425,14 +426,14 @@ TEST(PointManipulation, ExtendFromStartByDistance) { EXPECT_EQ(end, A.GetFirstPoint()); EXPECT_EQ(distance + extra_distance, A.GetDistance()); - A = Path(EMp, B, C); + A = Path(EMp, DetectorPosition(B), DetectorPosition(C)); extra_distance = -1; end = B + reverse * extra_distance; A.ExtendFromStartByDistance(extra_distance); EXPECT_EQ(end, A.GetFirstPoint()); EXPECT_EQ(distance + extra_distance, A.GetDistance()); - A = Path(EMp, B, C); + A = Path(EMp, DetectorPosition(B), DetectorPosition(C)); extra_distance = -101; end = B + reverse * extra_distance; A.ExtendFromStartByDistance(extra_distance); @@ -472,46 +473,46 @@ TEST_F(FakeLegacyDetectorModelTest, ExtendFromStartByColumnDepth) { direction = inner_p1 - inner_p0; distance = direction.magnitude(); direction.normalize(); - Path P(A, inner_p0, inner_p1); + Path P(A, DetectorPosition(inner_p0), DetectorPosition(inner_p1)); DensityDistribution1D const * density = dynamic_cast const *>(sector.density.get()); ASSERT_TRUE(density); double rho = density->Evaluate(Vector3D()); double sum = P.GetColumnDepthInBounds(); ASSERT_DOUBLE_EQ((inner_p1 - inner_p0).magnitude() * rho * 100, sum); - P = Path(A, p0, p1); + P = Path(A, DetectorPosition(p0), DetectorPosition(p1)); sum = P.GetColumnDepthInBounds(); ASSERT_DOUBLE_EQ((p1 - p0).magnitude() * rho * 100, sum); - P = Path(A, inner_p0, inner_p1); + P = Path(A, DetectorPosition(inner_p0), DetectorPosition(inner_p1)); P.EnsureIntersections(); double extra_distance = distance/3.0; double extra_column_depth = extra_distance * rho * 100; - extra_distance = A->DistanceForColumnDepthFromPoint(inner_p0, -direction, extra_column_depth); + extra_distance = A->DistanceForColumnDepthFromPoint(DetectorPosition(inner_p0), DetectorDirection(-direction), extra_column_depth); ASSERT_DOUBLE_EQ(distance/3.0, extra_distance); - ASSERT_DOUBLE_EQ(A->DistanceForColumnDepthFromPoint(P.GetFirstPoint(), -direction, extra_column_depth), A->DistanceForColumnDepthFromPoint(P.GetIntersections(), P.GetFirstPoint(), -direction, extra_column_depth)); + ASSERT_DOUBLE_EQ(A->DistanceForColumnDepthFromPoint(P.GetFirstPoint(), DetectorDirection(-direction), extra_column_depth), A->DistanceForColumnDepthFromPoint(P.GetIntersections(), P.GetFirstPoint(), DetectorDirection(-direction), extra_column_depth)); ASSERT_DOUBLE_EQ(extra_distance, P.GetDistanceFromStartInReverse(extra_column_depth)); P.ExtendFromStartByColumnDepth(extra_column_depth); Vector3D end = inner_p0 - direction * extra_distance; ASSERT_EQ(end, P.GetFirstPoint()); ASSERT_EQ(distance + extra_distance, P.GetDistance()); - P = Path(A, inner_p0, inner_p1); + P = Path(A, DetectorPosition(inner_p0), DetectorPosition(inner_p1)); P.EnsureIntersections(); extra_distance = -distance/3.0; extra_column_depth = extra_distance * rho * 100; P.ExtendFromStartByColumnDepth(extra_column_depth); - extra_distance = A->DistanceForColumnDepthFromPoint(inner_p0, -direction, extra_column_depth); + extra_distance = A->DistanceForColumnDepthFromPoint(DetectorPosition(inner_p0), DetectorDirection(-direction), extra_column_depth); EXPECT_DOUBLE_EQ(-distance/3.0, extra_distance); end = inner_p0 - direction * extra_distance; EXPECT_EQ(end, P.GetFirstPoint()); EXPECT_EQ(distance + extra_distance, P.GetDistance()); - P = Path(A, inner_p0, inner_p1); + P = Path(A, DetectorPosition(inner_p0), DetectorPosition(inner_p1)); P.EnsureIntersections(); extra_distance = -distance*1.5; extra_column_depth = extra_distance * rho * 100; P.ExtendFromStartByColumnDepth(extra_column_depth); - extra_distance = A->DistanceForColumnDepthFromPoint(inner_p0, -direction, extra_column_depth); + extra_distance = A->DistanceForColumnDepthFromPoint(DetectorPosition(inner_p0), DetectorDirection(-direction), extra_column_depth); EXPECT_DOUBLE_EQ(-distance*1.5, extra_distance); end = inner_p0 - direction * extra_distance; EXPECT_EQ(inner_p1, P.GetFirstPoint()); @@ -527,7 +528,7 @@ TEST(PointManipulation, ShrinkFromEndByDistance) { double distance = direction.magnitude(); direction.normalize(); - Path A(EMp, B, C); + Path A(EMp, DetectorPosition(B), DetectorPosition(C)); double extra_distance = -101; A.ShrinkFromEndByDistance(extra_distance); @@ -536,14 +537,14 @@ TEST(PointManipulation, ShrinkFromEndByDistance) { EXPECT_EQ(end, A.GetLastPoint()); EXPECT_EQ(distance - extra_distance, A.GetDistance()); - A = Path(EMp, B, C); + A = Path(EMp, DetectorPosition(B), DetectorPosition(C)); extra_distance = 1; end = C - direction * extra_distance; A.ShrinkFromEndByDistance(extra_distance); EXPECT_EQ(end, A.GetLastPoint()); EXPECT_EQ(distance - extra_distance, A.GetDistance()); - A = Path(EMp, B, C); + A = Path(EMp, DetectorPosition(B), DetectorPosition(C)); extra_distance = 101; end = C - direction * extra_distance; A.ShrinkFromEndByDistance(extra_distance); @@ -583,44 +584,44 @@ TEST_F(FakeLegacyDetectorModelTest, ShrinkFromEndByColumnDepth) { direction = inner_p1 - inner_p0; distance = direction.magnitude(); direction.normalize(); - Path P(A, inner_p0, inner_p1); + Path P(A, DetectorPosition(inner_p0), DetectorPosition(inner_p1)); DensityDistribution1D const * density = dynamic_cast const *>(sector.density.get()); ASSERT_TRUE(density); double rho = density->Evaluate(Vector3D()); double sum = P.GetColumnDepthInBounds(); ASSERT_DOUBLE_EQ((inner_p1 - inner_p0).magnitude() * rho * 100, sum); - P = Path(A, p0, p1); + P = Path(A, DetectorPosition(p0), DetectorPosition(p1)); sum = P.GetColumnDepthInBounds(); ASSERT_DOUBLE_EQ((p1 - p0).magnitude() * rho * 100, sum); - P = Path(A, inner_p0, inner_p1); + P = Path(A, DetectorPosition(inner_p0), DetectorPosition(inner_p1)); P.EnsureIntersections(); double extra_distance = distance/3.0; double extra_column_depth = extra_distance * rho * 100; P.ShrinkFromEndByColumnDepth(extra_column_depth); - extra_distance = A->DistanceForColumnDepthFromPoint(inner_p1, -direction, extra_column_depth); + extra_distance = A->DistanceForColumnDepthFromPoint(DetectorPosition(inner_p1), DetectorDirection(-direction), extra_column_depth); EXPECT_DOUBLE_EQ(distance/3.0, extra_distance); Vector3D end = inner_p1 - direction * extra_distance; EXPECT_EQ(end, P.GetLastPoint()); EXPECT_EQ(distance - extra_distance, P.GetDistance()); - P = Path(A, inner_p0, inner_p1); + P = Path(A, DetectorPosition(inner_p0), DetectorPosition(inner_p1)); P.EnsureIntersections(); extra_distance = -distance/3.0; extra_column_depth = extra_distance * rho * 100; P.ShrinkFromEndByColumnDepth(extra_column_depth); - extra_distance = A->DistanceForColumnDepthFromPoint(inner_p1, -direction, extra_column_depth); + extra_distance = A->DistanceForColumnDepthFromPoint(DetectorPosition(inner_p1), DetectorDirection(-direction), extra_column_depth); EXPECT_DOUBLE_EQ(-distance/3.0, extra_distance); end = inner_p1 - direction * extra_distance; EXPECT_EQ(end, P.GetLastPoint()); EXPECT_EQ(distance - extra_distance, P.GetDistance()); - P = Path(A, inner_p0, inner_p1); + P = Path(A, DetectorPosition(inner_p0), DetectorPosition(inner_p1)); P.EnsureIntersections(); extra_distance = distance*1.5; extra_column_depth = extra_distance * rho * 100; P.ShrinkFromEndByColumnDepth(extra_column_depth); - extra_distance = A->DistanceForColumnDepthFromPoint(inner_p1, -direction, extra_column_depth); + extra_distance = A->DistanceForColumnDepthFromPoint(DetectorPosition(inner_p1), DetectorDirection(-direction), extra_column_depth); EXPECT_DOUBLE_EQ(distance*1.5, extra_distance); end = inner_p1 - direction * extra_distance; EXPECT_EQ(inner_p0, P.GetLastPoint()); @@ -638,7 +639,7 @@ TEST(PointManipulation, ShrinkFromStartByDistance) { direction.normalize(); reverse.normalize(); - Path A(EMp, B, C); + Path A(EMp, DetectorPosition(B), DetectorPosition(C)); double extra_distance = -101; A.ShrinkFromStartByDistance(extra_distance); @@ -646,14 +647,14 @@ TEST(PointManipulation, ShrinkFromStartByDistance) { EXPECT_EQ(end, A.GetFirstPoint()); EXPECT_EQ(distance - extra_distance, A.GetDistance()); - A = Path(EMp, B, C); + A = Path(EMp, DetectorPosition(B), DetectorPosition(C)); extra_distance = 1; end = B - reverse * extra_distance; A.ShrinkFromStartByDistance(extra_distance); EXPECT_EQ(end, A.GetFirstPoint()); EXPECT_EQ(distance - extra_distance, A.GetDistance()); - A = Path(EMp, B, C); + A = Path(EMp, DetectorPosition(B), DetectorPosition(C)); extra_distance = 101; end = B - reverse * extra_distance; A.ShrinkFromStartByDistance(extra_distance); @@ -693,21 +694,21 @@ TEST_F(FakeLegacyDetectorModelTest, ShrinkFromStartByColumnDepth) { direction = inner_p1 - inner_p0; distance = direction.magnitude(); direction.normalize(); - Path P(A, inner_p0, inner_p1); + Path P(A, DetectorPosition(inner_p0), DetectorPosition(inner_p1)); DensityDistribution1D const * density = dynamic_cast const *>(sector.density.get()); ASSERT_TRUE(density); double rho = density->Evaluate(Vector3D()); double sum = P.GetColumnDepthInBounds(); ASSERT_DOUBLE_EQ((inner_p1 - inner_p0).magnitude() * rho * 100, sum); - P = Path(A, p0, p1); + P = Path(A, DetectorPosition(p0), DetectorPosition(p1)); sum = P.GetColumnDepthInBounds(); ASSERT_DOUBLE_EQ((p1 - p0).magnitude() * rho * 100, sum); - P = Path(A, inner_p0, inner_p1); + P = Path(A, DetectorPosition(inner_p0), DetectorPosition(inner_p1)); P.EnsureIntersections(); double extra_distance = distance/3.0; double extra_column_depth = extra_distance * rho * 100; - extra_distance = A->DistanceForColumnDepthFromPoint(inner_p0, direction, extra_column_depth); + extra_distance = A->DistanceForColumnDepthFromPoint(DetectorPosition(inner_p0), DetectorDirection(direction), extra_column_depth); ASSERT_DOUBLE_EQ(distance/3.0, extra_distance); ASSERT_DOUBLE_EQ(extra_distance, P.GetDistanceFromStartInReverse(extra_column_depth)); P.ShrinkFromStartByColumnDepth(extra_column_depth); @@ -715,23 +716,23 @@ TEST_F(FakeLegacyDetectorModelTest, ShrinkFromStartByColumnDepth) { ASSERT_EQ(end, P.GetFirstPoint()); ASSERT_EQ(distance - extra_distance, P.GetDistance()); - P = Path(A, inner_p0, inner_p1); + P = Path(A, DetectorPosition(inner_p0), DetectorPosition(inner_p1)); P.EnsureIntersections(); extra_distance = -distance/3.0; extra_column_depth = extra_distance * rho * 100; P.ShrinkFromStartByColumnDepth(extra_column_depth); - extra_distance = A->DistanceForColumnDepthFromPoint(inner_p0, direction, extra_column_depth); + extra_distance = A->DistanceForColumnDepthFromPoint(DetectorPosition(inner_p0), DetectorDirection(direction), extra_column_depth); EXPECT_DOUBLE_EQ(-distance/3.0, extra_distance); end = inner_p0 + direction * extra_distance; EXPECT_EQ(end, P.GetFirstPoint()); EXPECT_EQ(distance - extra_distance, P.GetDistance()); - P = Path(A, inner_p0, inner_p1); + P = Path(A, DetectorPosition(inner_p0), DetectorPosition(inner_p1)); P.EnsureIntersections(); extra_distance = distance*1.5; extra_column_depth = extra_distance * rho * 100; P.ShrinkFromStartByColumnDepth(extra_column_depth); - extra_distance = A->DistanceForColumnDepthFromPoint(inner_p0, direction, extra_column_depth); + extra_distance = A->DistanceForColumnDepthFromPoint(DetectorPosition(inner_p0), DetectorDirection(direction), extra_column_depth); EXPECT_DOUBLE_EQ(distance*1.5, extra_distance); end = inner_p0 + direction * extra_distance; EXPECT_EQ(inner_p1, P.GetFirstPoint()); @@ -747,7 +748,7 @@ TEST(PointManipulation, ExtendFromEndToDistance) { double distance = direction.magnitude(); direction.normalize(); - Path A(EMp, B, C); + Path A(EMp, DetectorPosition(B), DetectorPosition(C)); double target_distance = 101; double extra_distance = target_distance - distance; @@ -757,7 +758,7 @@ TEST(PointManipulation, ExtendFromEndToDistance) { EXPECT_EQ(end, A.GetLastPoint()); EXPECT_EQ(distance + extra_distance, A.GetDistance()); - A = Path(EMp, B, C); + A = Path(EMp, DetectorPosition(B), DetectorPosition(C)); target_distance = 1; extra_distance = target_distance - distance; end = C + direction * extra_distance; @@ -765,7 +766,7 @@ TEST(PointManipulation, ExtendFromEndToDistance) { EXPECT_EQ(C, A.GetLastPoint()); EXPECT_EQ(distance, A.GetDistance()); - A = Path(EMp, B, C); + A = Path(EMp, DetectorPosition(B), DetectorPosition(C)); target_distance = -101; extra_distance = target_distance - distance; end = C + direction * extra_distance; @@ -806,35 +807,35 @@ TEST_F(FakeLegacyDetectorModelTest, ExtendFromEndToColumnDepth) { direction = inner_p1 - inner_p0; distance = direction.magnitude(); direction.normalize(); - Path P(A, inner_p0, inner_p1); + Path P(A, DetectorPosition(inner_p0), DetectorPosition(inner_p1)); DensityDistribution1D const * density = dynamic_cast const *>(sector.density.get()); ASSERT_TRUE(density); double rho = density->Evaluate(Vector3D()); double sum = P.GetColumnDepthInBounds(); ASSERT_DOUBLE_EQ((inner_p1 - inner_p0).magnitude() * rho * 100, sum); - P = Path(A, p0, p1); + P = Path(A, DetectorPosition(p0), DetectorPosition(p1)); sum = P.GetColumnDepthInBounds(); ASSERT_DOUBLE_EQ((p1 - p0).magnitude() * rho * 100, sum); - P = Path(A, inner_p0, inner_p1); + P = Path(A, DetectorPosition(inner_p0), DetectorPosition(inner_p1)); P.EnsureIntersections(); double extra_distance = distance/3.0; double target_distance = distance + extra_distance; double target_column_depth = target_distance * rho * 100; - target_distance = A->DistanceForColumnDepthFromPoint(inner_p0, direction, target_column_depth); + target_distance = A->DistanceForColumnDepthFromPoint(DetectorPosition(inner_p0), DetectorDirection(direction), target_column_depth); ASSERT_DOUBLE_EQ(distance + distance/3.0, target_distance); ASSERT_DOUBLE_EQ(target_distance, P.GetDistanceFromStartAlongPath(target_column_depth)); P.ExtendFromEndToColumnDepth(target_column_depth); Vector3D end = inner_p0 + direction * target_distance; - ASSERT_TRUE((end - P.GetLastPoint()).magnitude() < 1e-6 * std::max(end.magnitude(), P.GetLastPoint().magnitude())); + ASSERT_TRUE((end - P.GetLastPoint()).magnitude() < 1e-6 * std::max(end.magnitude(), P.GetLastPoint()->magnitude())); ASSERT_DOUBLE_EQ(distance + extra_distance, P.GetDistance()); - P = Path(A, inner_p0, inner_p1); + P = Path(A, DetectorPosition(inner_p0), DetectorPosition(inner_p1)); P.EnsureIntersections(); extra_distance = -distance/3.0; target_distance = distance + extra_distance; target_column_depth = target_distance * rho * 100; - target_distance = A->DistanceForColumnDepthFromPoint(inner_p0, direction, target_column_depth); + target_distance = A->DistanceForColumnDepthFromPoint(DetectorPosition(inner_p0), DetectorDirection(direction), target_column_depth); P.ExtendFromEndToColumnDepth(target_column_depth); ASSERT_DOUBLE_EQ(distance - distance/3.0, target_distance); ASSERT_DOUBLE_EQ(target_distance, P.GetDistanceFromStartAlongPath(target_column_depth)); @@ -842,12 +843,12 @@ TEST_F(FakeLegacyDetectorModelTest, ExtendFromEndToColumnDepth) { ASSERT_EQ(inner_p1, P.GetLastPoint()); ASSERT_DOUBLE_EQ(distance, P.GetDistance()); - P = Path(A, inner_p0, inner_p1); + P = Path(A, DetectorPosition(inner_p0), DetectorPosition(inner_p1)); P.EnsureIntersections(); extra_distance = -distance*1.5; target_distance = distance + extra_distance; target_column_depth = target_distance * rho * 100; - target_distance = A->DistanceForColumnDepthFromPoint(inner_p0, direction, target_column_depth); + target_distance = A->DistanceForColumnDepthFromPoint(DetectorPosition(inner_p0), DetectorDirection(direction), target_column_depth); P.ExtendFromEndToColumnDepth(target_column_depth); ASSERT_DOUBLE_EQ(distance - distance*1.5, target_distance); ASSERT_DOUBLE_EQ(target_distance, P.GetDistanceFromStartAlongPath(target_column_depth)); @@ -868,7 +869,7 @@ TEST(PointManipulation, ExtendFromStartToDistance) { direction.normalize(); reverse.normalize(); - Path A(EMp, B, C); + Path A(EMp, DetectorPosition(B), DetectorPosition(C)); double target_distance = 101; double extra_distance = target_distance - distance; @@ -877,7 +878,7 @@ TEST(PointManipulation, ExtendFromStartToDistance) { EXPECT_EQ(end, A.GetFirstPoint()); EXPECT_EQ(distance + extra_distance, A.GetDistance()); - A = Path(EMp, B, C); + A = Path(EMp, DetectorPosition(B), DetectorPosition(C)); target_distance = 1; extra_distance = target_distance - distance; end = B + reverse * extra_distance; @@ -885,7 +886,7 @@ TEST(PointManipulation, ExtendFromStartToDistance) { EXPECT_EQ(B, A.GetFirstPoint()); EXPECT_EQ(distance , A.GetDistance()); - A = Path(EMp, B, C); + A = Path(EMp, DetectorPosition(B), DetectorPosition(C)); target_distance = -101; extra_distance = target_distance - distance; end = B + reverse * extra_distance; @@ -926,22 +927,22 @@ TEST_F(FakeLegacyDetectorModelTest, ExtendFromStartToColumnDepth) { direction = inner_p1 - inner_p0; distance = direction.magnitude(); direction.normalize(); - Path P(A, inner_p0, inner_p1); + Path P(A, DetectorPosition(inner_p0), DetectorPosition(inner_p1)); DensityDistribution1D const * density = dynamic_cast const *>(sector.density.get()); ASSERT_TRUE(density); double rho = density->Evaluate(Vector3D()); double sum = P.GetColumnDepthInBounds(); ASSERT_DOUBLE_EQ((inner_p1 - inner_p0).magnitude() * rho * 100, sum); - P = Path(A, p0, p1); + P = Path(A, DetectorPosition(p0), DetectorPosition(p1)); sum = P.GetColumnDepthInBounds(); ASSERT_DOUBLE_EQ((p1 - p0).magnitude() * rho * 100, sum); - P = Path(A, inner_p0, inner_p1); + P = Path(A, DetectorPosition(inner_p0), DetectorPosition(inner_p1)); P.EnsureIntersections(); double extra_distance = distance/3.0; double target_distance = distance + extra_distance; double target_column_depth = target_distance * rho * 100; - target_distance = A->DistanceForColumnDepthFromPoint(inner_p1, -direction, target_column_depth); + target_distance = A->DistanceForColumnDepthFromPoint(DetectorPosition(inner_p1), DetectorDirection(-direction), target_column_depth); ASSERT_DOUBLE_EQ(distance + distance/3.0, target_distance); ASSERT_DOUBLE_EQ(target_distance, P.GetDistanceFromEndInReverse(target_column_depth)); P.ExtendFromStartToColumnDepth(target_column_depth); @@ -949,14 +950,14 @@ TEST_F(FakeLegacyDetectorModelTest, ExtendFromStartToColumnDepth) { ASSERT_DOUBLE_EQ(target_column_depth, sum); Vector3D end = inner_p1 - direction * target_distance; ASSERT_DOUBLE_EQ(target_distance, P.GetDistance()); - ASSERT_TRUE((end - P.GetFirstPoint()).magnitude() < 1e-6 * std::max(end.magnitude(), P.GetFirstPoint().magnitude())); + ASSERT_TRUE((end - P.GetFirstPoint()).magnitude() < 1e-6 * std::max(end.magnitude(), P.GetFirstPoint()->magnitude())); - P = Path(A, inner_p0, inner_p1); + P = Path(A, DetectorPosition(inner_p0), DetectorPosition(inner_p1)); P.EnsureIntersections(); extra_distance = -distance/3.0; target_distance = distance + extra_distance; target_column_depth = target_distance * rho * 100; - target_distance = A->DistanceForColumnDepthFromPoint(inner_p1, -direction, target_column_depth); + target_distance = A->DistanceForColumnDepthFromPoint(DetectorPosition(inner_p1), DetectorDirection(-direction), target_column_depth); P.ExtendFromStartToColumnDepth(target_column_depth); ASSERT_DOUBLE_EQ(distance - distance/3.0, target_distance); ASSERT_DOUBLE_EQ(target_distance, P.GetDistanceFromEndInReverse(target_column_depth)); @@ -964,12 +965,12 @@ TEST_F(FakeLegacyDetectorModelTest, ExtendFromStartToColumnDepth) { ASSERT_EQ(inner_p0, P.GetFirstPoint()); ASSERT_DOUBLE_EQ(distance, P.GetDistance()); - P = Path(A, inner_p0, inner_p1); + P = Path(A, DetectorPosition(inner_p0), DetectorPosition(inner_p1)); P.EnsureIntersections(); extra_distance = -distance*1.5; target_distance = distance + extra_distance; target_column_depth = target_distance * rho * 100; - target_distance = A->DistanceForColumnDepthFromPoint(inner_p1, -direction, target_column_depth); + target_distance = A->DistanceForColumnDepthFromPoint(DetectorPosition(inner_p1), DetectorDirection(-direction), target_column_depth); P.ExtendFromStartToColumnDepth(target_column_depth); ASSERT_DOUBLE_EQ(distance - distance*1.5, target_distance); ASSERT_DOUBLE_EQ(target_distance, P.GetDistanceFromEndAlongPath(target_column_depth)); @@ -987,7 +988,7 @@ TEST(PointManipulation, ShrinkFromEndToDistance) { double distance = direction.magnitude(); direction.normalize(); - Path A(EMp, B, C); + Path A(EMp, DetectorPosition(B), DetectorPosition(C)); double target_distance = 101; double extra_distance = distance - target_distance; @@ -997,7 +998,7 @@ TEST(PointManipulation, ShrinkFromEndToDistance) { EXPECT_EQ(C, A.GetLastPoint()); EXPECT_EQ(distance, A.GetDistance()); - A = Path(EMp, B, C); + A = Path(EMp, DetectorPosition(B), DetectorPosition(C)); target_distance = 1; extra_distance = distance - target_distance; end = C - direction * extra_distance; @@ -1005,7 +1006,7 @@ TEST(PointManipulation, ShrinkFromEndToDistance) { EXPECT_EQ(end, A.GetLastPoint()); EXPECT_EQ(distance - extra_distance, A.GetDistance()); - A = Path(EMp, B, C); + A = Path(EMp, DetectorPosition(B), DetectorPosition(C)); target_distance = -101; extra_distance = distance - target_distance; end = C - direction * extra_distance; @@ -1046,50 +1047,50 @@ TEST_F(FakeLegacyDetectorModelTest, ShrinkFromEndToColumnDepth) { direction = inner_p1 - inner_p0; distance = direction.magnitude(); direction.normalize(); - Path P(A, inner_p0, inner_p1); + Path P(A, DetectorPosition(inner_p0), DetectorPosition(inner_p1)); DensityDistribution1D const * density = dynamic_cast const *>(sector.density.get()); ASSERT_TRUE(density); double rho = density->Evaluate(Vector3D()); double sum = P.GetColumnDepthInBounds(); ASSERT_DOUBLE_EQ((inner_p1 - inner_p0).magnitude() * rho * 100, sum); - P = Path(A, p0, p1); + P = Path(A, DetectorPosition(p0), DetectorPosition(p1)); sum = P.GetColumnDepthInBounds(); ASSERT_DOUBLE_EQ((p1 - p0).magnitude() * rho * 100, sum); - P = Path(A, inner_p0, inner_p1); + P = Path(A, DetectorPosition(inner_p0), DetectorPosition(inner_p1)); P.EnsureIntersections(); double extra_distance = distance/3.0; double extra_column_depth = extra_distance * rho * 100; double target_distance = distance - extra_distance; double target_column_depth = target_distance * rho * 100; P.ShrinkFromEndToColumnDepth(target_column_depth); - extra_distance = A->DistanceForColumnDepthFromPoint(inner_p1, -direction, extra_column_depth); + extra_distance = A->DistanceForColumnDepthFromPoint(DetectorPosition(inner_p1), DetectorDirection(-direction), extra_column_depth); EXPECT_DOUBLE_EQ(distance/3.0, extra_distance); Vector3D end = inner_p1 - direction * extra_distance; - ASSERT_TRUE((end - P.GetLastPoint()).magnitude() < 1e-6 * std::max(end.magnitude(), P.GetLastPoint().magnitude())); + ASSERT_TRUE((end - P.GetLastPoint()).magnitude() < 1e-6 * std::max(end.magnitude(), P.GetLastPoint()->magnitude())); EXPECT_DOUBLE_EQ(distance - extra_distance, P.GetDistance()); - P = Path(A, inner_p0, inner_p1); + P = Path(A, DetectorPosition(inner_p0), DetectorPosition(inner_p1)); P.EnsureIntersections(); extra_distance = -distance/3.0; extra_column_depth = extra_distance * rho * 100; target_distance = distance - extra_distance; target_column_depth = target_distance * rho * 100; P.ShrinkFromEndToColumnDepth(target_column_depth); - extra_distance = A->DistanceForColumnDepthFromPoint(inner_p1, -direction, extra_column_depth); + extra_distance = A->DistanceForColumnDepthFromPoint(DetectorPosition(inner_p1), DetectorDirection(-direction), extra_column_depth); EXPECT_DOUBLE_EQ(-distance/3.0, extra_distance); end = inner_p1 - direction * extra_distance; EXPECT_EQ(inner_p1, P.GetLastPoint()); EXPECT_DOUBLE_EQ(distance, P.GetDistance()); - P = Path(A, inner_p0, inner_p1); + P = Path(A, DetectorPosition(inner_p0), DetectorPosition(inner_p1)); P.EnsureIntersections(); extra_distance = distance*1.5; extra_column_depth = extra_distance * rho * 100; target_distance = distance - extra_distance; target_column_depth = target_distance * rho * 100; P.ShrinkFromEndToColumnDepth(target_column_depth); - extra_distance = A->DistanceForColumnDepthFromPoint(inner_p1, -direction, extra_column_depth); + extra_distance = A->DistanceForColumnDepthFromPoint(DetectorPosition(inner_p1), DetectorDirection(-direction), extra_column_depth); EXPECT_DOUBLE_EQ(distance*1.5, extra_distance); end = inner_p1 - direction * extra_distance; EXPECT_EQ(inner_p0, P.GetLastPoint()); @@ -1107,7 +1108,7 @@ TEST(PointManipulation, ShrinkFromStartToDistance) { direction.normalize(); reverse.normalize(); - Path A(EMp, B, C); + Path A(EMp, DetectorPosition(B), DetectorPosition(C)); double target_distance = 101; double extra_distance = distance - target_distance; @@ -1116,7 +1117,7 @@ TEST(PointManipulation, ShrinkFromStartToDistance) { EXPECT_EQ(B, A.GetFirstPoint()); EXPECT_EQ(distance, A.GetDistance()); - A = Path(EMp, B, C); + A = Path(EMp, DetectorPosition(B), DetectorPosition(C)); target_distance = 1; extra_distance = distance - target_distance; end = B - reverse * extra_distance; @@ -1124,7 +1125,7 @@ TEST(PointManipulation, ShrinkFromStartToDistance) { EXPECT_EQ(end, A.GetFirstPoint()); EXPECT_EQ(distance - extra_distance, A.GetDistance()); - A = Path(EMp, B, C); + A = Path(EMp, DetectorPosition(B), DetectorPosition(C)); target_distance = -101; extra_distance = distance - target_distance; end = B - reverse * extra_distance; @@ -1165,23 +1166,23 @@ TEST_F(FakeLegacyDetectorModelTest, ShrinkFromStartToColumnDepth) { direction = inner_p1 - inner_p0; distance = direction.magnitude(); direction.normalize(); - Path P(A, inner_p0, inner_p1); + Path P(A, DetectorPosition(inner_p0), DetectorPosition(inner_p1)); DensityDistribution1D const * density = dynamic_cast const *>(sector.density.get()); ASSERT_TRUE(density); double rho = density->Evaluate(Vector3D()); double sum = P.GetColumnDepthInBounds(); ASSERT_DOUBLE_EQ((inner_p1 - inner_p0).magnitude() * rho * 100, sum); - P = Path(A, p0, p1); + P = Path(A, DetectorPosition(p0), DetectorPosition(p1)); sum = P.GetColumnDepthInBounds(); ASSERT_DOUBLE_EQ((p1 - p0).magnitude() * rho * 100, sum); - P = Path(A, inner_p0, inner_p1); + P = Path(A, DetectorPosition(inner_p0), DetectorPosition(inner_p1)); P.EnsureIntersections(); double extra_distance = distance/3.0; double target_distance = distance + extra_distance; double extra_column_depth = extra_distance * rho * 100; double target_column_depth = target_distance * rho * 100; - extra_distance = A->DistanceForColumnDepthFromPoint(inner_p0, -direction, extra_column_depth); + extra_distance = A->DistanceForColumnDepthFromPoint(DetectorPosition(inner_p0), DetectorDirection(-direction), extra_column_depth); ASSERT_DOUBLE_EQ(distance/3.0, extra_distance); ASSERT_DOUBLE_EQ(extra_distance, P.GetDistanceFromStartInReverse(extra_column_depth)); P.ShrinkFromStartToColumnDepth(target_column_depth); @@ -1189,30 +1190,30 @@ TEST_F(FakeLegacyDetectorModelTest, ShrinkFromStartToColumnDepth) { ASSERT_EQ(inner_p0, P.GetFirstPoint()); ASSERT_EQ(distance, P.GetDistance()); - P = Path(A, inner_p0, inner_p1); + P = Path(A, DetectorPosition(inner_p0), DetectorPosition(inner_p1)); P.EnsureIntersections(); extra_distance = -distance/3.0; target_distance = distance + extra_distance; extra_column_depth = extra_distance * rho * 100; target_column_depth = target_distance * rho * 100; P.ShrinkFromStartToColumnDepth(target_column_depth); - extra_distance = A->DistanceForColumnDepthFromPoint(inner_p0, -direction, extra_column_depth); + extra_distance = A->DistanceForColumnDepthFromPoint(DetectorPosition(inner_p0), DetectorDirection(-direction), extra_column_depth); EXPECT_DOUBLE_EQ(-distance/3.0, extra_distance); end = inner_p0 - direction * extra_distance; - ASSERT_TRUE((end - P.GetFirstPoint()).magnitude() < 1e-6 * std::max(end.magnitude(), P.GetFirstPoint().magnitude())); + ASSERT_TRUE((end - P.GetFirstPoint()).magnitude() < 1e-6 * std::max(end.magnitude(), P.GetFirstPoint()->magnitude())); EXPECT_DOUBLE_EQ(distance + extra_distance, P.GetDistance()); - P = Path(A, inner_p0, inner_p1); + P = Path(A, DetectorPosition(inner_p0), DetectorPosition(inner_p1)); P.EnsureIntersections(); extra_distance = distance*1.5; target_distance = distance + extra_distance; extra_column_depth = extra_distance * rho * 100; target_column_depth = target_distance * rho * 100; P.ShrinkFromStartToColumnDepth(target_column_depth); - extra_distance = A->DistanceForColumnDepthFromPoint(inner_p0, direction, extra_column_depth); + extra_distance = A->DistanceForColumnDepthFromPoint(DetectorPosition(inner_p0), DetectorDirection(direction), extra_column_depth); EXPECT_DOUBLE_EQ(distance*1.5, extra_distance); end = inner_p0 + direction * extra_distance; - ASSERT_TRUE((inner_p0 - P.GetFirstPoint()).magnitude() < 1e-6 * std::max(inner_p0.magnitude(), P.GetFirstPoint().magnitude())); + ASSERT_TRUE((inner_p0 - P.GetFirstPoint()).magnitude() < 1e-6 * std::max(inner_p0.magnitude(), P.GetFirstPoint()->magnitude())); EXPECT_EQ(distance, P.GetDistance()); } } @@ -1243,7 +1244,7 @@ TEST_F(FakeLegacyDetectorModelTest, GetColumnDepthInBounds) direction.normalize(); Vector3D inner_p0 = p0 + direction * distance / 4.0; Vector3D inner_p1 = p1 - direction * distance / 4.0; - Path P(A, inner_p0, inner_p1); + Path P(A, DetectorPosition(inner_p0), DetectorPosition(inner_p1)); DensityDistribution1D const * density = dynamic_cast const *>(sector.density.get()); ASSERT_TRUE(density); double rho = density->Evaluate(Vector3D()); @@ -1281,7 +1282,7 @@ TEST_F(FakeLegacyDetectorModelTest, GetColumnDepthFromStartInBounds) direction = inner_p1 - inner_p0; distance = direction.magnitude(); direction.normalize(); - Path P(A, inner_p0, inner_p1); + Path P(A, DetectorPosition(inner_p0), DetectorPosition(inner_p1)); DensityDistribution1D const * density = dynamic_cast const *>(sector.density.get()); ASSERT_TRUE(density); double rho = density->Evaluate(Vector3D()); @@ -1332,7 +1333,7 @@ TEST_F(FakeLegacyDetectorModelTest, GetColumnDepthFromEndInBounds) direction = inner_p1 - inner_p0; distance = direction.magnitude(); direction.normalize(); - Path P(A, inner_p0, inner_p1); + Path P(A, DetectorPosition(inner_p0), DetectorPosition(inner_p1)); DensityDistribution1D const * density = dynamic_cast const *>(sector.density.get()); ASSERT_TRUE(density); double rho = density->Evaluate(Vector3D()); @@ -1383,7 +1384,7 @@ TEST_F(FakeLegacyDetectorModelTest, GetColumnDepthFromStartAlongPath) direction = inner_p1 - inner_p0; distance = direction.magnitude(); direction.normalize(); - Path P(A, inner_p0, inner_p1); + Path P(A, DetectorPosition(inner_p0), DetectorPosition(inner_p1)); DensityDistribution1D const * density = dynamic_cast const *>(sector.density.get()); ASSERT_TRUE(density); double rho = density->Evaluate(Vector3D()); @@ -1434,7 +1435,7 @@ TEST_F(FakeLegacyDetectorModelTest, GetColumnDepthFromEndAlongPath) direction = inner_p1 - inner_p0; distance = direction.magnitude(); direction.normalize(); - Path P(A, inner_p0, inner_p1); + Path P(A, DetectorPosition(inner_p0), DetectorPosition(inner_p1)); DensityDistribution1D const * density = dynamic_cast const *>(sector.density.get()); ASSERT_TRUE(density); double rho = density->Evaluate(Vector3D()); @@ -1485,7 +1486,7 @@ TEST_F(FakeLegacyDetectorModelTest, GetColumnDepthFromStartInReverse) direction = inner_p1 - inner_p0; distance = direction.magnitude(); direction.normalize(); - Path P(A, inner_p0, inner_p1); + Path P(A, DetectorPosition(inner_p0), DetectorPosition(inner_p1)); DensityDistribution1D const * density = dynamic_cast const *>(sector.density.get()); ASSERT_TRUE(density); double rho = density->Evaluate(Vector3D()); @@ -1536,7 +1537,7 @@ TEST_F(FakeLegacyDetectorModelTest, GetColumnDepthFromEndInReverse) direction = inner_p1 - inner_p0; distance = direction.magnitude(); direction.normalize(); - Path P(A, inner_p0, inner_p1); + Path P(A, DetectorPosition(inner_p0), DetectorPosition(inner_p1)); DensityDistribution1D const * density = dynamic_cast const *>(sector.density.get()); ASSERT_TRUE(density); double rho = density->Evaluate(Vector3D()); @@ -1589,7 +1590,7 @@ TEST_F(FakeLegacyDetectorModelTest, GetDistanceFromStartInBounds) direction = inner_p1 - inner_p0; distance = direction.magnitude(); direction.normalize(); - Path P(A, inner_p0, inner_p1); + Path P(A, DetectorPosition(inner_p0), DetectorPosition(inner_p1)); DensityDistribution1D const * density = dynamic_cast const *>(sector.density.get()); ASSERT_TRUE(density); double rho = density->Evaluate(Vector3D()); @@ -1640,7 +1641,7 @@ TEST_F(FakeLegacyDetectorModelTest, GetDistanceFromEndInBounds) direction = inner_p1 - inner_p0; distance = direction.magnitude(); direction.normalize(); - Path P(A, inner_p0, inner_p1); + Path P(A, DetectorPosition(inner_p0), DetectorPosition(inner_p1)); DensityDistribution1D const * density = dynamic_cast const *>(sector.density.get()); ASSERT_TRUE(density); double rho = density->Evaluate(Vector3D()); @@ -1691,7 +1692,7 @@ TEST_F(FakeLegacyDetectorModelTest, GetDistanceFromStartAlongPath) direction = inner_p1 - inner_p0; distance = direction.magnitude(); direction.normalize(); - Path P(A, inner_p0, inner_p1); + Path P(A, DetectorPosition(inner_p0), DetectorPosition(inner_p1)); DensityDistribution1D const * density = dynamic_cast const *>(sector.density.get()); ASSERT_TRUE(density); double rho = density->Evaluate(Vector3D()); @@ -1742,7 +1743,7 @@ TEST_F(FakeLegacyDetectorModelTest, GetDistanceFromEndAlongPath) direction = inner_p1 - inner_p0; distance = direction.magnitude(); direction.normalize(); - Path P(A, inner_p0, inner_p1); + Path P(A, DetectorPosition(inner_p0), DetectorPosition(inner_p1)); DensityDistribution1D const * density = dynamic_cast const *>(sector.density.get()); ASSERT_TRUE(density); double rho = density->Evaluate(Vector3D()); @@ -1793,7 +1794,7 @@ TEST_F(FakeLegacyDetectorModelTest, GetDistanceFromStartInReverse) direction = inner_p1 - inner_p0; distance = direction.magnitude(); direction.normalize(); - Path P(A, inner_p0, inner_p1); + Path P(A, DetectorPosition(inner_p0), DetectorPosition(inner_p1)); DensityDistribution1D const * density = dynamic_cast const *>(sector.density.get()); ASSERT_TRUE(density); double rho = density->Evaluate(Vector3D()); @@ -1844,7 +1845,7 @@ TEST_F(FakeLegacyDetectorModelTest, GetDistanceFromEndInReverse) direction = inner_p1 - inner_p0; distance = direction.magnitude(); direction.normalize(); - Path P(A, inner_p0, inner_p1); + Path P(A, DetectorPosition(inner_p0), DetectorPosition(inner_p1)); DensityDistribution1D const * density = dynamic_cast const *>(sector.density.get()); ASSERT_TRUE(density); double rho = density->Evaluate(Vector3D()); From 7b7267978d824ef7273c8b04d10def0c8f0ba120 Mon Sep 17 00:00:00 2001 From: Austin Schneider Date: Tue, 16 Jan 2024 13:08:14 -0700 Subject: [PATCH 20/42] Fix overloaded function casts in DetectorModel pybindings --- .../private/pybindings/DetectorModel.h | 107 ++++++++---------- 1 file changed, 48 insertions(+), 59 deletions(-) diff --git a/projects/detector/private/pybindings/DetectorModel.h b/projects/detector/private/pybindings/DetectorModel.h index 88ecf4600..0702c062d 100644 --- a/projects/detector/private/pybindings/DetectorModel.h +++ b/projects/detector/private/pybindings/DetectorModel.h @@ -6,6 +6,7 @@ #include #include "../../public/LeptonInjector/detector/DetectorModel.h" +#include "../../public/LeptonInjector/detector/Coordinates.h" #include "../../../geometry/public/LeptonInjector/geometry/Geometry.h" void register_DetectorModel(pybind11::module_ & m) { @@ -19,26 +20,26 @@ void register_DetectorModel(pybind11::module_ & m) { .def("LoadDetectorModel",&DetectorModel::LoadDetectorModel) .def("LoadMaterialModel",&DetectorModel::LoadMaterialModel) .def("GetMassDensity", ( - double (DetectorModel::*)(LI::geometry::Geometry::IntersectionList const &, LI::math::Vector3D const &) const + double (DetectorModel::*)(LI::geometry::Geometry::IntersectionList const &, DetectorPosition const &) const )(&DetectorModel::GetMassDensity)) .def("GetMassDensity", ( - double (DetectorModel::*)(LI::math::Vector3D const &) const + double (DetectorModel::*)(DetectorPosition const &) const )(&DetectorModel::GetMassDensity)) .def("GetParticleDensity", ( - double (DetectorModel::*)(LI::geometry::Geometry::IntersectionList const &, LI::math::Vector3D const &, LI::dataclasses::Particle::ParticleType) const + double (DetectorModel::*)(LI::geometry::Geometry::IntersectionList const &, DetectorPosition const &, LI::dataclasses::Particle::ParticleType) const )(&DetectorModel::GetParticleDensity)) .def("GetParticleDensity", ( - double (DetectorModel::*)(LI::math::Vector3D const &, LI::dataclasses::Particle::ParticleType) const + double (DetectorModel::*)(DetectorPosition const &, LI::dataclasses::Particle::ParticleType) const )(&DetectorModel::GetParticleDensity)) .def("GetInteractionDensity", ( double (DetectorModel::*)(LI::geometry::Geometry::IntersectionList const &, - LI::math::Vector3D const &, + DetectorPosition const &, std::vector const &, std::vector const &, double const &) const )(&DetectorModel::GetInteractionDensity)) .def("GetInteractionDensity", ( - double (DetectorModel::*)(LI::math::Vector3D const &, + double (DetectorModel::*)(DetectorPosition const &, std::vector const &, std::vector const &, double const &) const @@ -46,75 +47,75 @@ void register_DetectorModel(pybind11::module_ & m) { .def("GetColumnDepthInCGS", ( double (DetectorModel::*)( LI::geometry::Geometry::IntersectionList const &, - LI::math::Vector3D const &, - LI::math::Vector3D const &) const + DetectorPosition const &, + DetectorPosition const &) const )(&DetectorModel::GetColumnDepthInCGS)) .def("GetColumnDepthInCGS", ( double (DetectorModel::*)( - LI::math::Vector3D const &, - LI::math::Vector3D const &) const + DetectorPosition const &, + DetectorPosition const &) const )(&DetectorModel::GetColumnDepthInCGS)) .def("DistanceForColumnDepthFromPoint", ( double (DetectorModel::*)( LI::geometry::Geometry::IntersectionList const &, - LI::math::Vector3D const &, - LI::math::Vector3D const &, + DetectorPosition const &, + DetectorDirection const &, double) const )(&DetectorModel::DistanceForColumnDepthFromPoint)) .def("DistanceForColumnDepthFromPoint", ( double (DetectorModel::*)( - LI::math::Vector3D const &, - LI::math::Vector3D const &, + DetectorPosition const &, + DetectorDirection const &, double) const )(&DetectorModel::DistanceForColumnDepthFromPoint)) .def("DistanceForColumnDepthToPoint", ( double (DetectorModel::*)( LI::geometry::Geometry::IntersectionList const &, - LI::math::Vector3D const &, - LI::math::Vector3D const &, + DetectorPosition const &, + DetectorDirection const &, double) const )(&DetectorModel::DistanceForColumnDepthToPoint)) .def("DistanceForColumnDepthToPoint", ( double (DetectorModel::*)( - LI::math::Vector3D const &, - LI::math::Vector3D const &, + DetectorPosition const &, + DetectorDirection const &, double) const )(&DetectorModel::DistanceForColumnDepthToPoint)) .def("GetMassDensity", ( double (DetectorModel::*)( LI::geometry::Geometry::IntersectionList const &, - LI::math::Vector3D const &, + DetectorPosition const &, std::set) const )(&DetectorModel::GetMassDensity)) .def("GetMassDensity", ( double (DetectorModel::*)( - LI::math::Vector3D const &, + DetectorPosition const &, std::set) const )(&DetectorModel::GetMassDensity)) .def("GetParticleDensity", ( std::vector (DetectorModel::*)( LI::geometry::Geometry::IntersectionList const &, - LI::math::Vector3D const &, + DetectorPosition const &, std::set) const )(&DetectorModel::GetParticleDensity)) .def("GetParticleDensity", ( std::vector (DetectorModel::*)( - LI::math::Vector3D const &, + DetectorPosition const &, std::set) const )(&DetectorModel::GetParticleDensity)) .def("GetInteractionDepthInCGS", ( double (DetectorModel::*)( LI::geometry::Geometry::IntersectionList const &, - LI::math::Vector3D const &, - LI::math::Vector3D const &, + DetectorPosition const &, + DetectorPosition const &, std::vector const &, std::vector const &, double const &) const )(&DetectorModel::GetInteractionDepthInCGS)) .def("GetInteractionDepthInCGS", ( double (DetectorModel::*)( - LI::math::Vector3D const &, - LI::math::Vector3D const &, + DetectorPosition const &, + DetectorPosition const &, std::vector const &, std::vector const &, double const &) const @@ -122,8 +123,8 @@ void register_DetectorModel(pybind11::module_ & m) { .def("DistanceForInteractionDepthFromPoint", ( double (DetectorModel::*)( LI::geometry::Geometry::IntersectionList const &, - LI::math::Vector3D const &, - LI::math::Vector3D const &, + DetectorPosition const &, + DetectorDirection const &, double, std::vector const &, std::vector const &, @@ -131,8 +132,8 @@ void register_DetectorModel(pybind11::module_ & m) { )(&DetectorModel::DistanceForInteractionDepthFromPoint)) .def("DistanceForInteractionDepthFromPoint", ( double (DetectorModel::*)( - LI::math::Vector3D const &, - LI::math::Vector3D const &, + DetectorPosition const &, + DetectorDirection const &, double, std::vector const &, std::vector const &, @@ -141,8 +142,8 @@ void register_DetectorModel(pybind11::module_ & m) { .def("DistanceForInteractionDepthToPoint", ( double (DetectorModel::*)( LI::geometry::Geometry::IntersectionList const &, - LI::math::Vector3D const &, - LI::math::Vector3D const &, + DetectorPosition const &, + DetectorDirection const &, double, std::vector const &, std::vector const &, @@ -150,8 +151,8 @@ void register_DetectorModel(pybind11::module_ & m) { )(&DetectorModel::DistanceForInteractionDepthToPoint)) .def("DistanceForInteractionDepthToPoint", ( double (DetectorModel::*)( - LI::math::Vector3D const &, - LI::math::Vector3D const &, + DetectorPosition const &, + DetectorDirection const &, double, std::vector const &, std::vector const &, @@ -160,35 +161,19 @@ void register_DetectorModel(pybind11::module_ & m) { .def("GetParticleColumnDepth", ( std::vector (DetectorModel::*)( LI::geometry::Geometry::IntersectionList const &, - LI::math::Vector3D const &, - LI::math::Vector3D const &, + DetectorPosition const &, + DetectorPosition const &, std::vector const &) const )(&DetectorModel::GetParticleColumnDepth)) .def("GetContainingSector", ( DetectorSector (DetectorModel::*)( LI::geometry::Geometry::IntersectionList const &, - LI::math::Vector3D const & p0) const + DetectorPosition const & p0) const )(&DetectorModel::GetContainingSector)) .def("GetContainingSector", ( DetectorSector (DetectorModel::*)( - LI::math::Vector3D const & p0) const + DetectorPosition const & p0) const )(&DetectorModel::GetContainingSector)) - .def("GetEarthCoordPosFromDetCoordPos", ( - LI::math::Vector3D (DetectorModel::*)( - LI::math::Vector3D const &) const - )(&DetectorModel::GetEarthCoordPosFromDetCoordPos)) - .def("GetEarthCoordDirFromDetCoordDir", ( - LI::math::Vector3D (DetectorModel::*)( - LI::math::Vector3D const &) const - )(&DetectorModel::GetEarthCoordDirFromDetCoordDir)) - .def("GetDetCoordPosFromEarthCoordPos", ( - LI::math::Vector3D (DetectorModel::*)( - LI::math::Vector3D const &) const - )(&DetectorModel::GetDetCoordPosFromEarthCoordPos)) - .def("GetDetCoordDirFromEarthCoordDir", ( - LI::math::Vector3D (DetectorModel::*)( - LI::math::Vector3D const &) const - )(&DetectorModel::GetDetCoordDirFromEarthCoordDir)) .def_property("Path", &DetectorModel::GetPath, &DetectorModel::SetPath) .def_property("Materials", &DetectorModel::GetMaterials, &DetectorModel::SetMaterials) .def_property("Sectors", &DetectorModel::GetSectors, &DetectorModel::SetSectors) @@ -196,7 +181,11 @@ void register_DetectorModel(pybind11::module_ & m) { .def("AddSector", &DetectorModel::AddSector) .def("GetSector", &DetectorModel::GetSector) .def("ClearSectors", &DetectorModel::ClearSectors) - .def("GetIntersections", &DetectorModel::GetIntersections) + .def("GetIntersections", ( + LI::geometry::Geometry::IntersectionList (DetectorModel::*)( + DetectorPosition const &, + DetectorDirection const &) const + )(&DetectorModel::GetIntersections)) .def_static("SortIntersections", ( void (*)(LI::geometry::Geometry::IntersectionList &) )(&DetectorModel::SortIntersections)) @@ -213,17 +202,17 @@ void register_DetectorModel(pybind11::module_ & m) { }) .def("GetOuterBounds", ( LI::geometry::Geometry::IntersectionList (DetectorModel::*)( - LI::math::Vector3D const &, - LI::math::Vector3D const &) const + DetectorPosition const &, + DetectorDirection const &) const )(&DetectorModel::GetOuterBounds)) .def("GetAvailableTargets", ( std::set (DetectorModel::*)( LI::geometry::Geometry::IntersectionList const &, - std::array const &) const + DetectorPosition const &) const )(&DetectorModel::GetAvailableTargets)) .def("GetAvailableTargets", ( std::set (DetectorModel::*)( - std::array const &) const + DetectorPosition const &) const )(&DetectorModel::GetAvailableTargets)) .def("GetTargetMass", &DetectorModel::GetTargetMass) .def("GetMaterials",&DetectorModel::GetMaterials) From dd61e70646452713eae4b66281e67d6f07913558 Mon Sep 17 00:00:00 2001 From: Austin Schneider Date: Tue, 16 Jan 2024 16:15:25 -0700 Subject: [PATCH 21/42] Perform correct coordinate system conversions in DetectorModel_TEST --- .../private/test/DetectorModel_TEST.cxx | 41 ++++++++++++------- 1 file changed, 26 insertions(+), 15 deletions(-) diff --git a/projects/detector/private/test/DetectorModel_TEST.cxx b/projects/detector/private/test/DetectorModel_TEST.cxx index b3dd74845..3f1decc16 100644 --- a/projects/detector/private/test/DetectorModel_TEST.cxx +++ b/projects/detector/private/test/DetectorModel_TEST.cxx @@ -410,6 +410,8 @@ TEST_F(FakeLegacyDetectorModelTest, LegacyFileConstantIntegralInternal) double min_radius = sphere->GetInnerRadius(); Vector3D p0 = RandomVector(max_radius, min_radius); Vector3D p1 = RandomVector(max_radius, min_radius); + p0 = A.ToDet(GeometryPosition(p0)); + p1 = A.ToDet(GeometryPosition(p1)); DensityDistribution1D const * density = dynamic_cast const *>(sector.density.get()); ASSERT_TRUE(density); double rho = density->Evaluate(Vector3D()); @@ -438,6 +440,7 @@ TEST_F(FakeLegacyDetectorModelTest, LegacyFileConstantGetMassDensityInternal) double max_radius = sphere->GetRadius(); double min_radius = sphere->GetInnerRadius(); Vector3D p0 = RandomVector(max_radius, min_radius); + p0 = A.ToDet(GeometryPosition(p0)); DensityDistribution1D const * density_dist = dynamic_cast const *>(sector.density.get()); ASSERT_TRUE(density_dist); double rho = density_dist->Evaluate(Vector3D()); @@ -472,6 +475,10 @@ TEST_F(FakeLegacyDetectorModelTest, LegacyFileConstantIntegralNested) double min_radius = sphere_0->GetInnerRadius(); Vector3D p0 = RandomVector(max_radius, min_radius); Vector3D p1 = RandomVector(max_radius, min_radius); + GeometryPosition p0_geo(p0); + GeometryPosition p1_geo(p1); + DetectorPosition p0_det(A.ToDet(p0_geo)); + DetectorPosition p1_det(A.ToDet(p1_geo)); double distance = (p1-p0).magnitude(); DensityDistribution1D const * density_0 = dynamic_cast const *>(sector_0.density.get()); DensityDistribution1D const * density_1 = dynamic_cast const *>(sector_1.density.get()); @@ -484,15 +491,15 @@ TEST_F(FakeLegacyDetectorModelTest, LegacyFileConstantIntegralNested) bool in_0 = p0.magnitude() <= sphere_0->GetRadius(); bool in_1 = p1.magnitude() <= sphere_0->GetRadius(); double integral = 0.0; - double sum = A.GetColumnDepthInCGS(DetectorPosition(p0), DetectorPosition(p1)); + double sum = A.GetColumnDepthInCGS(p0_det, p1_det); if(in_0 and in_1) { integral = rho_0 * (p1-p0).magnitude(); - ASSERT_DOUBLE_EQ(integral * 100, sum); + ASSERT_NEAR(integral * 100, sum, std::min(std::abs(sum), std::abs(integral * 100)) * 1e-8); } else { - Vector3D direction = p1 - p0; + Vector3D direction = p1_geo - p0_geo; direction.normalize(); - std::vector intersections = sphere_0->Intersections(p0, direction); + std::vector intersections = sphere_0->Intersections(p0_geo, direction); if(intersections.size() > 1) { double dist_0 = intersections[0].distance; double dist_1 = intersections[1].distance; @@ -501,11 +508,11 @@ TEST_F(FakeLegacyDetectorModelTest, LegacyFileConstantIntegralNested) if((not in_0) and in_1) { integral += rho_1*near; integral += rho_0*((p1-p0).magnitude() - near); - ASSERT_DOUBLE_EQ(integral * 100, sum); + ASSERT_NEAR(integral * 100, sum, std::min(std::abs(integral * 100), std::abs(sum)) * 1e-12); } else if(in_0 and not in_1) { integral += rho_0*near; integral += rho_1*((p1-p0).magnitude() - near); - ASSERT_DOUBLE_EQ(integral * 100, sum); + ASSERT_NEAR(integral * 100, sum, std::min(std::abs(integral * 100), std::abs(sum)) * 1e-12); } else if((not in_0) and (not in_1)) { double far = std::max(std::max(dist_0, 0.0), std::max(dist_1, 0.0)); if(near < distance) { @@ -519,12 +526,12 @@ TEST_F(FakeLegacyDetectorModelTest, LegacyFileConstantIntegralNested) } else { integral += rho_1*distance; } - ASSERT_DOUBLE_EQ(integral * 100, sum); + ASSERT_NEAR(integral * 100, sum, std::min(std::abs(integral * 100), std::abs(sum)) * 1e-12); } } else if (dist_0 <= 0 and dist_1 <= 0) { integral = rho_1 * (p1-p0).magnitude(); - ASSERT_DOUBLE_EQ(integral * 100, sum); + ASSERT_NEAR(integral * 100, sum, std::min(std::abs(integral * 100), std::abs(sum)) * 1e-12); } else { double dist = std::max(dist_0, dist_1); @@ -533,7 +540,7 @@ TEST_F(FakeLegacyDetectorModelTest, LegacyFileConstantIntegralNested) } else if(in_0 and not in_1) { integral += rho_0*dist; integral += rho_1*((p1-p0).magnitude() - dist); - ASSERT_DOUBLE_EQ(integral * 100, sum); + ASSERT_NEAR(integral * 100, sum, std::min(std::abs(integral * 100), std::abs(sum)) * 1e-12); } else if((not in_0) and (not in_1)) { assert(false); } @@ -541,10 +548,10 @@ TEST_F(FakeLegacyDetectorModelTest, LegacyFileConstantIntegralNested) } else { integral = rho_1 * (p1-p0).magnitude(); - ASSERT_DOUBLE_EQ(integral * 100, sum); + ASSERT_NEAR(integral * 100, sum, std::min(std::abs(integral * 100), std::abs(sum)) * 1e-8); } } - ASSERT_DOUBLE_EQ(integral * 100, sum); + ASSERT_NEAR(integral * 100, sum, std::min(std::abs(integral * 100), std::abs(sum)) * 1e-12); } } @@ -573,16 +580,18 @@ TEST_F(FakeLegacyDetectorModelTest, LegacyFileConstantGetMassDensityNested) double max_radius = sphere_1->GetRadius(); double min_radius = sphere_0->GetInnerRadius(); Vector3D p0 = RandomVector(max_radius, min_radius); - Vector3D p1 = RandomVector(max_radius, min_radius); + GeometryPosition p0_geo(p0); + DetectorPosition p0_det(A.ToDet(p0_geo)); DensityDistribution1D const * density_0 = dynamic_cast const *>(sector_0.density.get()); DensityDistribution1D const * density_1 = dynamic_cast const *>(sector_1.density.get()); ASSERT_TRUE(density_0); ASSERT_TRUE(density_1); double rho_0 = density_0->Evaluate(Vector3D()); double rho_1 = density_1->Evaluate(Vector3D()); - ASSERT_LE(p0.magnitude(), max_radius); - bool in_0 = p0.magnitude() <= sphere_0->GetRadius(); - double density = A.GetMassDensity(DetectorPosition(p0)); + ASSERT_LE(min_radius, max_radius); + ASSERT_LE(p0_geo->magnitude(), max_radius); + bool in_0 = p0_geo->magnitude() <= sphere_0->GetRadius(); + double density = A.GetMassDensity(p0_det); if(in_0) { ASSERT_DOUBLE_EQ(density, rho_0); } @@ -1152,6 +1161,8 @@ TEST_F(FakeLegacyDetectorModelTest, LegacyFileConstantInverseIntegralInternal) Vector3D p0 = RandomVector(max_radius, min_radius); Vector3D p1 = RandomVector(max_radius, min_radius); + p0 = A.ToDet(GeometryPosition(p0)); + p1 = A.ToDet(GeometryPosition(p1)); double sum = A.GetColumnDepthInCGS(DetectorPosition(p0), DetectorPosition(p1)); Vector3D direction = p1 - p0; From bab7af2d851d98bb5a14bb50b1c92027a6187bd1 Mon Sep 17 00:00:00 2001 From: Austin Schneider Date: Tue, 16 Jan 2024 17:14:42 -0700 Subject: [PATCH 22/42] Add pybindings for Coordinates --- .../detector/private/pybindings/Coordinates.h | 31 +++++++++++++++++++ .../detector/private/pybindings/detector.cxx | 2 ++ 2 files changed, 33 insertions(+) create mode 100644 projects/detector/private/pybindings/Coordinates.h diff --git a/projects/detector/private/pybindings/Coordinates.h b/projects/detector/private/pybindings/Coordinates.h new file mode 100644 index 000000000..764f5ab62 --- /dev/null +++ b/projects/detector/private/pybindings/Coordinates.h @@ -0,0 +1,31 @@ +#include +#include + +#include +#include + +#include "../../../math/public/LeptonInjector/math/Vector3D.h" +#include "../../public/LeptonInjector/detector/DetectorModel.h" +#include "../../public/LeptonInjector/detector/Coordinates.h" + +void register_Coordinates(pybind11::module_ & m) { + using namespace pybind11; + using namespace LI::math; + using namespace LI::detector; + + // Register only `DetectorPosition` and `DetectorDirection` + // The python user does not need to know about or manipulate the Geometry coordinates + + class_>(m, "DetectorPosition") + .def(init()) + ; + + class_>(m, "DetectorDirection") + .def(init()) + ; + + // Allow implicit conversion from Vector3D to DetectorPosition and DetectorDirection + // But only on the python side! + implicitly_convertible(); + implicitly_convertible(); +} diff --git a/projects/detector/private/pybindings/detector.cxx b/projects/detector/private/pybindings/detector.cxx index 37a33cbc6..d37700979 100644 --- a/projects/detector/private/pybindings/detector.cxx +++ b/projects/detector/private/pybindings/detector.cxx @@ -1,6 +1,7 @@ #include #include +#include "./Coordinates.h" #include "./DetectorModel.h" #include "./DetectorSector.h" #include "./Path.h" @@ -20,6 +21,7 @@ using namespace pybind11; PYBIND11_MODULE(detector,m) { + register_Coordinates(m); register_DetectorModel(m); register_DetectorSector(m); register_Path(m); From 5d11f688057eacd18410dc982559c8a84caf86d4 Mon Sep 17 00:00:00 2001 From: Austin Schneider Date: Tue, 16 Jan 2024 17:15:16 -0700 Subject: [PATCH 23/42] Update pybindings for Path to use the new interface --- projects/detector/private/pybindings/Path.h | 204 +++++++++++++++----- 1 file changed, 156 insertions(+), 48 deletions(-) diff --git a/projects/detector/private/pybindings/Path.h b/projects/detector/private/pybindings/Path.h index 34d552067..a600c1296 100644 --- a/projects/detector/private/pybindings/Path.h +++ b/projects/detector/private/pybindings/Path.h @@ -5,6 +5,7 @@ #include #include "../../public/LeptonInjector/detector/DetectorModel.h" +#include "../../public/LeptonInjector/detector/Coordinates.h" #include "../../public/LeptonInjector/detector/Path.h" void register_Path(pybind11::module_ & m) { @@ -13,52 +14,159 @@ void register_Path(pybind11::module_ & m) { class_>(m, "Path") .def(init>()) - .def(init, LI::math::Vector3D const &, LI::math::Vector3D const &>()) - .def(init, LI::math::Vector3D const &, LI::math::Vector3D const &, double>()) - - .def("HasDetectorModel",&Path::HasDetectorModel) - .def("HasPoints",&Path::HasPoints) - .def("HasIntersections",&Path::HasIntersections) - .def("HasColumnDepth",&Path::HasColumnDepth) - .def("GetDetectorModel",&Path::GetDetectorModel) - .def("HasDetectorModel",&Path::HasDetectorModel) - .def("GetFirstPoint",&Path::GetFirstPoint) - .def("GetLastPoint",&Path::GetLastPoint) - .def("GetDirection",&Path::GetDirection) - .def("GetDistance",&Path::GetDistance) - - .def("ClipToOuterBounds",&Path::ClipToOuterBounds) - - .def("ExtendFromEndByDistance",&Path::ExtendFromEndByDistance) - .def("ExtendFromStartByDistance",&Path::ExtendFromStartByDistance) - .def("ShrinkFromEndByDistance",&Path::ShrinkFromEndByDistance) - .def("ShrinkFromStartByDistance",&Path::ShrinkFromStartByDistance) - .def("ExtendFromEndByColumnDepth",&Path::ExtendFromEndByColumnDepth) - .def("ExtendFromStartByColumnDepth",&Path::ExtendFromStartByColumnDepth) - .def("ShrinkFromEndByColumnDepth",&Path::ShrinkFromEndByColumnDepth) - .def("ShrinkFromStartByColumnDepth",&Path::ShrinkFromStartByColumnDepth) - .def("ExtendFromEndByInteractionDepth",&Path::ExtendFromEndByInteractionDepth) - .def("ExtendFromStartByInteractionDepth",&Path::ExtendFromStartByInteractionDepth) - .def("ShrinkFromEndByInteractionDepth",&Path::ShrinkFromEndByInteractionDepth) - .def("ShrinkFromStartByInteractionDepth",&Path::ShrinkFromStartByInteractionDepth) - .def("ExtendFromEndToDistance",&Path::ExtendFromEndToDistance) - .def("ExtendFromStartToDistance",&Path::ExtendFromStartToDistance) - .def("ShrinkFromEndToDistance",&Path::ShrinkFromEndToDistance) - .def("ShrinkFromStartToDistance",&Path::ShrinkFromStartToDistance) - .def("ExtendFromEndToColumnDepth",&Path::ExtendFromEndToColumnDepth) - .def("ExtendFromStartToColumnDepth",&Path::ExtendFromStartToColumnDepth) - .def("ShrinkFromEndToColumnDepth",&Path::ShrinkFromEndToColumnDepth) - .def("ShrinkFromStartToColumnDepth",&Path::ShrinkFromStartToColumnDepth) - .def("ExtendFromEndToInteractionDepth",&Path::ExtendFromEndToInteractionDepth) - .def("ExtendFromStartToInteractionDepth",&Path::ExtendFromStartToInteractionDepth) - .def("ShrinkFromEndToInteractionDepth",&Path::ShrinkFromEndToInteractionDepth) - .def("ShrinkFromStartToInteractionDepth",&Path::ShrinkFromStartToInteractionDepth) - .def("GetColumnDepthInBounds",&Path::GetColumnDepthInBounds) - .def("GetInteractionDepthInBounds",&Path::GetInteractionDepthInBounds) - - .def("GetDistanceFromStartAlongPath",overload_cast const &, std::vector< double > const &, double const &>(&Path::GetDistanceFromStartAlongPath)) - .def("GetDistanceFromStartAlongPath",overload_cast(&Path::GetDistanceFromStartAlongPath)) - - .def("GetDistanceFromStartInReverse",overload_cast const &, std::vector< double > const &, double const &>(&Path::GetDistanceFromStartInReverse)) - .def("GetDistanceFromStartInReverse",overload_cast(&Path::GetDistanceFromStartInReverse)); + .def(init, DetectorPosition const &, DetectorPosition const &>()) + .def(init, DetectorPosition const &, DetectorDirection const &, double>()) + + .def("HasDetectorModel", &Path::HasDetectorModel) + .def("HasPoints", &Path::HasPoints) + .def("HasIntersections", &Path::HasIntersections) + .def("HasColumnDepth", &Path::HasColumnDepth) + + .def("GetDetectorModel", &Path::GetDetectorModel) + .def("GetFirstPoint", &Path::GetFirstPoint) + .def("GetLastPoint", &Path::GetLastPoint) + .def("GetDirection", &Path::GetDirection) + .def("GetGeoFirstPoint", &Path::GetGeoFirstPoint) + .def("GetGeoLastPoint", &Path::GetGeoLastPoint) + .def("GetGeoDirection", &Path::GetGeoDirection) + .def("GetDistance", &Path::GetDistance) + .def("GetIntersections", &Path::GetIntersections) + + .def("SetDetectorModel", &Path::SetDetectorModel) + .def("EnsureDetectorModel", &Path::EnsureDetectorModel) + + .def("SetPoints", overload_cast< + DetectorPosition, + DetectorPosition + >(&Path::SetPoints)) + .def("SetPointsWithRay", overload_cast< + DetectorPosition, + DetectorDirection, + double + >(&Path::SetPointsWithRay)) + .def("EnsurePoints", &Path::EnsurePoints) + + .def("SetIntersections", &Path::SetIntersections) + .def("ComputeIntersections", &Path::ComputeIntersections) + .def("EnsureIntersections", &Path::EnsureIntersections) + + .def("ClipToOuterBounds", &Path::ClipToOuterBounds) + + .def("Flip", &Path::Flip) + + // Extend/Shrink By + .def("ExtendFromEndByDistance", &Path::ExtendFromEndByDistance) + .def("ExtendFromStartByDistance", &Path::ExtendFromStartByDistance) + .def("ShrinkFromEndByDistance", &Path::ShrinkFromEndByDistance) + .def("ShrinkFromStartByDistance", &Path::ShrinkFromStartByDistance) + + .def("ExtendFromEndByColumnDepth", &Path::ExtendFromEndByColumnDepth) + .def("ExtendFromStartByColumnDepth", &Path::ExtendFromStartByColumnDepth) + .def("ShrinkFromEndByColumnDepth", &Path::ShrinkFromEndByColumnDepth) + .def("ShrinkFromStartByColumnDepth", &Path::ShrinkFromStartByColumnDepth) + + .def("ExtendFromEndByInteractionDepth", &Path::ExtendFromEndByInteractionDepth) + .def("ExtendFromStartByInteractionDepth", &Path::ExtendFromStartByInteractionDepth) + .def("ShrinkFromEndByInteractionDepth", &Path::ShrinkFromEndByInteractionDepth) + .def("ShrinkFromStartByInteractionDepth", &Path::ShrinkFromStartByInteractionDepth) + + // Extend/Shrink To + .def("ExtendFromEndToDistance", &Path::ExtendFromEndToDistance) + .def("ExtendFromStartToDistance", &Path::ExtendFromStartToDistance) + .def("ShrinkFromEndToDistance", &Path::ShrinkFromEndToDistance) + .def("ShrinkFromStartToDistance", &Path::ShrinkFromStartToDistance) + + .def("ExtendFromEndToColumnDepth", &Path::ExtendFromEndToColumnDepth) + .def("ExtendFromStartToColumnDepth", &Path::ExtendFromStartToColumnDepth) + .def("ShrinkFromEndToColumnDepth", &Path::ShrinkFromEndToColumnDepth) + .def("ShrinkFromStartToColumnDepth", &Path::ShrinkFromStartToColumnDepth) + + .def("ExtendFromEndToInteractionDepth", &Path::ExtendFromEndToInteractionDepth) + .def("ExtendFromStartToInteractionDepth", &Path::ExtendFromStartToInteractionDepth) + .def("ShrinkFromEndToInteractionDepth", &Path::ShrinkFromEndToInteractionDepth) + .def("ShrinkFromStartToInteractionDepth", &Path::ShrinkFromStartToInteractionDepth) + // + + // Get + .def("GetColumnDepthInBounds", &Path::GetColumnDepthInBounds) + .def("GetInteractionDepthInBounds", &Path::GetInteractionDepthInBounds) + // + + // Get * From + .def("GetColumnDepthFromStartInBounds", &Path::GetColumnDepthFromStartInBounds) + .def("GetColumnDepthFromEndInBounds", &Path::GetColumnDepthFromEndInBounds) + .def("GetColumnDepthFromStartAlongPath", &Path::GetColumnDepthFromStartAlongPath) + .def("GetColumnDepthFromEndAlongPath", &Path::GetColumnDepthFromEndAlongPath) + .def("GetColumnDepthFromStartInReverse", &Path::GetColumnDepthFromStartInReverse) + .def("GetColumnDepthFromEndInReverse", &Path::GetColumnDepthFromEndInReverse) + + .def("GetInteractionDepthFromStartInBounds", &Path::GetInteractionDepthFromStartInBounds) + .def("GetInteractionDepthFromEndInBounds", &Path::GetInteractionDepthFromEndInBounds) + .def("GetInteractionDepthFromStartAlongPath", &Path::GetInteractionDepthFromStartAlongPath) + .def("GetInteractionDepthFromEndAlongPath", &Path::GetInteractionDepthFromEndAlongPath) + .def("GetInteractionDepthFromStartInReverse", &Path::GetInteractionDepthFromStartInReverse) + .def("GetInteractionDepthFromEndInReverse", &Path::GetInteractionDepthFromEndInReverse) + // + + // Get Distance From + .def("GetDistanceFromStartInBounds", + overload_cast + (&Path::GetDistanceFromStartInBounds)) + .def("GetDistanceFromEndInBounds", + overload_cast + (&Path::GetDistanceFromEndInBounds)) + .def("GetDistanceFromStartAlongPath", + overload_cast + (&Path::GetDistanceFromStartAlongPath)) + .def("GetDistanceFromEndAlongPath", + overload_cast + (&Path::GetDistanceFromEndAlongPath)) + .def("GetDistanceFromStartInReverse", + overload_cast + (&Path::GetDistanceFromStartInReverse)) + .def("GetDistanceFromEndInReverse", + overload_cast + (&Path::GetDistanceFromEndInReverse)) + + .def("GetDistanceFromStartInBounds", overload_cast< + double, + std::vector const &, + std::vector const &, + double const & + >(&Path::GetDistanceFromStartInBounds)) + .def("GetDistanceFromEndInBounds", overload_cast< + double, + std::vector const &, + std::vector const &, + double const & + >(&Path::GetDistanceFromEndInBounds)) + .def("GetDistanceFromStartAlongPath", overload_cast< + double, + std::vector const &, + std::vector const &, + double const & + >(&Path::GetDistanceFromStartAlongPath)) + .def("GetDistanceFromEndAlongPath", overload_cast< + double, + std::vector const &, + std::vector const &, + double const & + >(&Path::GetDistanceFromEndAlongPath)) + .def("GetDistanceFromStartInReverse", overload_cast< + double, + std::vector const &, + std::vector const &, + double const & + >(&Path::GetDistanceFromStartInReverse)) + .def("GetDistanceFromEndInReverse", overload_cast< + double, + std::vector const &, + std::vector const &, + double const & + >(&Path::GetDistanceFromEndInReverse)) + // + + .def("IsWithinBounds", overload_cast(&Path::IsWithinBounds)) + .def("GetDistanceFromStartInBounds", overload_cast(&Path::GetDistanceFromStartInBounds)) + ; } From 7d7d978ce10cd53414d77b306064b4c61b064515 Mon Sep 17 00:00:00 2001 From: Austin Schneider Date: Tue, 16 Jan 2024 18:25:16 -0700 Subject: [PATCH 24/42] Rework Flip() --- projects/detector/private/Path.cxx | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/projects/detector/private/Path.cxx b/projects/detector/private/Path.cxx index b1238e814..824835846 100644 --- a/projects/detector/private/Path.cxx +++ b/projects/detector/private/Path.cxx @@ -245,10 +245,10 @@ void Path::ClipToOuterBounds() { } void Path::Flip() { - EnsurePoints(); std::swap(first_point_, last_point_); + std::swap(first_point_det_, last_point_det_); direction_ *= -1; - set_det_points_ = false; + direction_det_ *= -1; } From 559c4c62e6aeae1d3838107693c4ece09dda1795 Mon Sep 17 00:00:00 2001 From: Austin Schneider Date: Tue, 16 Jan 2024 18:25:42 -0700 Subject: [PATCH 25/42] Start to fix Path_TEST --- projects/detector/private/test/Path_TEST.cxx | 97 ++++++++++++++++---- 1 file changed, 80 insertions(+), 17 deletions(-) diff --git a/projects/detector/private/test/Path_TEST.cxx b/projects/detector/private/test/Path_TEST.cxx index 565c47971..2dc4227b5 100644 --- a/projects/detector/private/test/Path_TEST.cxx +++ b/projects/detector/private/test/Path_TEST.cxx @@ -182,22 +182,32 @@ TEST(Points, SetGetRay) { EXPECT_EQ(distance, A.GetDistance()); } -TEST(EnsurePoints, Throw) { +TEST(EnsurePoints, ThrowDefault) { Path A; EXPECT_THROW(A.EnsurePoints(), std::runtime_error); } -TEST(EnsurePoints, NoThrow) { +TEST(EnsurePoints, ThrowDetectorPositions) { Vector3D A(1,2,3); Vector3D B(4,6,8); Vector3D direction = B - A; direction.normalize(); Path C; C.SetPoints(DetectorPosition(A), DetectorPosition(B)); + EXPECT_THROW(C.EnsurePoints(), std::runtime_error); +} + +TEST(EnsurePoints, NoThrowGeometryPositions) { + Vector3D A(1,2,3); + Vector3D B(4,6,8); + Vector3D direction = B - A; + direction.normalize(); + Path C; + C.SetPoints(GeometryPosition(A), GeometryPosition(B)); EXPECT_NO_THROW(C.EnsurePoints()); } -TEST(EnsurePoints, NoThrowRay) { +TEST(EnsurePoints, ThrowRayDetectorPositions) { Vector3D B(1, 2, 3); Vector3D C(4, 6, 8); Vector3D direction = C - B; @@ -210,6 +220,22 @@ TEST(EnsurePoints, NoThrowRay) { direction.normalize(); Path A; A.SetPointsWithRay(DetectorPosition(B), DetectorDirection(direction), distance); + EXPECT_THROW(A.EnsurePoints(), std::runtime_error); +} + +TEST(EnsurePoints, NoThrowRayGeometryPositions) { + Vector3D B(1, 2, 3); + Vector3D C(4, 6, 8); + Vector3D direction = C - B; + double distance = direction.magnitude(); + direction.normalize(); + direction.normalize(); + direction.normalize(); + direction.normalize(); + direction.normalize(); + direction.normalize(); + Path A; + A.SetPointsWithRay(GeometryPosition(B), GeometryDirection(direction), distance); EXPECT_NO_THROW(A.EnsurePoints()); } @@ -233,6 +259,14 @@ TEST(EnsureIntersections, Throw) { A = Path(); A.SetPointsWithRay(DetectorPosition(B), DetectorDirection(direction), distance); EXPECT_THROW(A.EnsureIntersections(), std::runtime_error); + + A = Path(); + A.SetPoints(GeometryPosition(B), GeometryPosition(C)); + EXPECT_THROW(A.EnsureIntersections(), std::runtime_error); + + A = Path(); + A.SetPointsWithRay(GeometryPosition(B), GeometryDirection(direction), distance); + EXPECT_THROW(A.EnsureIntersections(), std::runtime_error); } TEST(EnsureIntersections, NoThrow) { @@ -255,6 +289,16 @@ TEST(EnsureIntersections, NoThrow) { A.SetPointsWithRay(DetectorPosition(B), DetectorDirection(direction), distance); EXPECT_NO_THROW(A.EnsureIntersections()); + A = Path(); + A.SetDetectorModel(EMp); + A.SetPoints(GeometryPosition(B), GeometryPosition(C)); + EXPECT_NO_THROW(A.EnsureIntersections()); + + A = Path(); + A.SetDetectorModel(EMp); + A.SetPointsWithRay(GeometryPosition(B), GeometryDirection(direction), distance); + EXPECT_NO_THROW(A.EnsureIntersections()); + A = Path(EMp); A.SetPoints(DetectorPosition(B), DetectorPosition(C)); EXPECT_NO_THROW(A.EnsureIntersections()); @@ -263,11 +307,25 @@ TEST(EnsureIntersections, NoThrow) { A.SetPointsWithRay(DetectorPosition(B), DetectorDirection(direction), distance); EXPECT_NO_THROW(A.EnsureIntersections()); + A = Path(EMp); + A.SetPoints(GeometryPosition(B), GeometryPosition(C)); + EXPECT_NO_THROW(A.EnsureIntersections()); + + A = Path(EMp); + A.SetPointsWithRay(GeometryPosition(B), GeometryDirection(direction), distance); + EXPECT_NO_THROW(A.EnsureIntersections()); + A = Path(EMp, DetectorPosition(B), DetectorPosition(C)); EXPECT_NO_THROW(A.EnsureIntersections()); A = Path(EMp, DetectorPosition(B), DetectorDirection(direction), distance); EXPECT_NO_THROW(A.EnsureIntersections()); + + A = Path(EMp, GeometryPosition(B), GeometryPosition(C)); + EXPECT_NO_THROW(A.EnsureIntersections()); + + A = Path(EMp, GeometryPosition(B), GeometryDirection(direction), distance); + EXPECT_NO_THROW(A.EnsureIntersections()); } TEST(PointManipulation, Flip) { @@ -352,10 +410,14 @@ TEST_F(FakeLegacyDetectorModelTest, ExtendFromEndByColumnDepth) { Vector3D p0 = RandomVector(max_radius, min_radius); Vector3D p1 = RandomVector(max_radius, min_radius); Vector3D direction = p1 - p0; + DetectorPosition p0_det = A->ToDet(GeometryPosition(p0)); + DetectorPosition p1_det = A->ToDet(GeometryPosition(p1)); double distance = direction.magnitude(); direction.normalize(); Vector3D inner_p0 = p0 + direction * distance / 4.0; Vector3D inner_p1 = p1 - direction * distance / 4.0; + DetectorPosition inner_p0_det = A->ToDet(GeometryPosition(inner_p0)); + DetectorPosition inner_p1_det = A->ToDet(GeometryPosition(inner_p1)); ASSERT_TRUE(p0.magnitude() < max_radius); ASSERT_TRUE(p1.magnitude() < max_radius); ASSERT_TRUE(inner_p0.magnitude() < max_radius); @@ -363,47 +425,48 @@ TEST_F(FakeLegacyDetectorModelTest, ExtendFromEndByColumnDepth) { direction = inner_p1 - inner_p0; distance = direction.magnitude(); direction.normalize(); - Path P(A, DetectorPosition(inner_p0), DetectorPosition(inner_p1)); + DetectorDirection direction_det = A->ToDet(GeometryDirection(direction)); + Path P(A, GeometryPosition(inner_p0), GeometryPosition(inner_p1)); DensityDistribution1D const * density = dynamic_cast const *>(sector.density.get()); ASSERT_TRUE(density); double rho = density->Evaluate(Vector3D()); double sum = P.GetColumnDepthInBounds(); ASSERT_DOUBLE_EQ((inner_p1 - inner_p0).magnitude() * rho * 100, sum); - P = Path(A, DetectorPosition(p0), DetectorPosition(p1)); + P = Path(A, GeometryPosition(p0), GeometryPosition(p1)); sum = P.GetColumnDepthInBounds(); ASSERT_DOUBLE_EQ((p1 - p0).magnitude() * rho * 100, sum); - P = Path(A, DetectorPosition(inner_p0), DetectorPosition(inner_p1)); + P = Path(A, GeometryPosition(inner_p0), GeometryPosition(inner_p1)); P.EnsureIntersections(); double extra_distance = distance/3.0; double extra_column_depth = extra_distance * rho * 100; P.ExtendFromEndByColumnDepth(extra_column_depth); - extra_distance = A->DistanceForColumnDepthFromPoint(DetectorPosition(inner_p1), DetectorDirection(direction), extra_column_depth); + extra_distance = A->DistanceForColumnDepthFromPoint(inner_p1_det, direction_det, extra_column_depth); EXPECT_DOUBLE_EQ(distance/3.0, extra_distance); - Vector3D end = inner_p1 + direction * extra_distance; - EXPECT_EQ(end, P.GetLastPoint()); + Vector3D end = inner_p1_det + direction_det.get() * extra_distance; + EXPECT_NEAR((end - P.GetLastPoint()).magnitude(), 0, 1e-8); EXPECT_EQ(distance + extra_distance, P.GetDistance()); - P = Path(A, DetectorPosition(inner_p0), DetectorPosition(inner_p1)); + P = Path(A, GeometryPosition(inner_p0), GeometryPosition(inner_p1)); P.EnsureIntersections(); extra_distance = -distance/3.0; extra_column_depth = extra_distance * rho * 100; P.ExtendFromEndByColumnDepth(extra_column_depth); - extra_distance = A->DistanceForColumnDepthFromPoint(DetectorPosition(inner_p1), DetectorDirection(direction), extra_column_depth); + extra_distance = A->DistanceForColumnDepthFromPoint(inner_p1_det, direction_det, extra_column_depth); EXPECT_DOUBLE_EQ(-distance/3.0, extra_distance); - end = inner_p1 + direction * extra_distance; - EXPECT_EQ(end, P.GetLastPoint()); + end = inner_p1_det + direction_det.get() * extra_distance; + EXPECT_NEAR((end - P.GetLastPoint()).magnitude(), 0, 1e-8); EXPECT_EQ(distance + extra_distance, P.GetDistance()); - P = Path(A, DetectorPosition(inner_p0), DetectorPosition(inner_p1)); + P = Path(A, GeometryPosition(inner_p0), GeometryPosition(inner_p1)); P.EnsureIntersections(); extra_distance = -distance*1.5; extra_column_depth = extra_distance * rho * 100; P.ExtendFromEndByColumnDepth(extra_column_depth); - extra_distance = A->DistanceForColumnDepthFromPoint(DetectorPosition(inner_p1), DetectorDirection(direction), extra_column_depth); + extra_distance = A->DistanceForColumnDepthFromPoint(inner_p1_det, direction_det, extra_column_depth); EXPECT_DOUBLE_EQ(-distance*1.5, extra_distance); - end = inner_p1 + direction * extra_distance; - EXPECT_EQ(inner_p0, P.GetLastPoint()); + end = inner_p1_det + direction_det.get() * extra_distance; + EXPECT_NEAR((inner_p0_det - P.GetLastPoint())->magnitude(), 0, 1e-8); EXPECT_EQ(0, P.GetDistance()); } } From fabcd4f07b254f613d3b4d45ebae770161145884 Mon Sep 17 00:00:00 2001 From: Austin Schneider Date: Tue, 16 Jan 2024 18:33:18 -0700 Subject: [PATCH 26/42] Add build badge --- README.md | 2 ++ 1 file changed, 2 insertions(+) diff --git a/README.md b/README.md index 025bc0948..f14cacac5 100644 --- a/README.md +++ b/README.md @@ -1,3 +1,5 @@ +[![Build](https://github.com/Harvard-Neutrino/LeptonInjector/actions/workflows/build_wheels.yml/badge.svg)](https://github.com/Harvard-Neutrino/LeptonInjector/actions/workflows/build_wheels.yml) + # LeptonInjector LeptonInjector is a framework for injecting and weighting interaction final states of complex topology, with specific concern for the detector geometry. LeptonInjector is designed to support a wide variety of neutrino experimental setups, including atmospheric neutrinos, accelerator beam decay-in-flight neutrinos, and neutrinos from decay-at-rest sources. The original [LeptonInjector-v1](https://github.com/icecube/LeptonInjector) was developed within the IceCube collaboration to study atmospheric and astrophysical neutrino interactions in the IceCube detector. From 98c05617168273e8f2b220aed31f39c08724ffae Mon Sep 17 00:00:00 2001 From: Austin Schneider Date: Wed, 17 Jan 2024 13:40:38 -0700 Subject: [PATCH 27/42] Add arithmetic operators with doubles --- .../LeptonInjector/detector/Coordinates.h | 54 +++++++++++++++++-- 1 file changed, 50 insertions(+), 4 deletions(-) diff --git a/projects/detector/public/LeptonInjector/detector/Coordinates.h b/projects/detector/public/LeptonInjector/detector/Coordinates.h index bab3578de..5efb25564 100644 --- a/projects/detector/public/LeptonInjector/detector/Coordinates.h +++ b/projects/detector/public/LeptonInjector/detector/Coordinates.h @@ -5,14 +5,60 @@ #include "LeptonInjector/math/Vector3D.h" // for Vector3D #include +namespace fluent { + +template +struct AltArithmetic { + +template +struct MultiplicableBy : crtp +{ + FLUENT_NODISCARD constexpr T operator*(O const& other) const + { + return T(this->underlying().get() * other); + } + FLUENT_CONSTEXPR17 T& operator*=(O const& other) + { + this->underlying().get() *= other; + return this->underlying(); + } +}; + +template +struct DivisibleBy : crtp +{ + FLUENT_NODISCARD constexpr T operator/(O const& other) const + { + return T(this->underlying().get() / other); + } + FLUENT_CONSTEXPR17 T& operator/=(O const& other) + { + this->underlying().get() /= other; + return this->underlying(); + } +}; + +template +struct MultiplicableTo : crtp +{ + FLUENT_NODISCARD constexpr O operator*(T const& other) const + { + return O(this->underlying().get() * other.get()); + } +}; + +}; + +} // namespace fluent + namespace LI { namespace detector { -using GeometryPosition = fluent::NamedType; -using GeometryDirection = fluent::NamedType; -using DetectorPosition = fluent::NamedType; -using DetectorDirection = fluent::NamedType; +using GeometryPosition = fluent::NamedType::MultiplicableBy, fluent::AltArithmetic::DivisibleBy, fluent::AltArithmetic::MultiplicableTo>; +using GeometryDirection = fluent::NamedType::MultiplicableBy, fluent::AltArithmetic::DivisibleBy, fluent::AltArithmetic::MultiplicableTo>; +using DetectorPosition = fluent::NamedType::MultiplicableBy, fluent::AltArithmetic::DivisibleBy, fluent::AltArithmetic::MultiplicableTo>; +using DetectorDirection = fluent::NamedType::MultiplicableBy, fluent::AltArithmetic::DivisibleBy, fluent::AltArithmetic::MultiplicableTo>; } // namespace detector } // namespace LI From 39463aef00f37d431e8a0e203f2a09fbd7cf959c Mon Sep 17 00:00:00 2001 From: Austin Schneider Date: Wed, 17 Jan 2024 13:55:16 -0700 Subject: [PATCH 28/42] Allow Geometry coordinates on the python side. Only return Vector3D in Path --- .../detector/private/pybindings/Coordinates.h | 24 +++++++++++++++++-- .../private/pybindings/DetectorModel.h | 16 +++++++++++-- projects/detector/private/pybindings/Path.h | 7 +++--- 3 files changed, 40 insertions(+), 7 deletions(-) diff --git a/projects/detector/private/pybindings/Coordinates.h b/projects/detector/private/pybindings/Coordinates.h index 764f5ab62..046928ecd 100644 --- a/projects/detector/private/pybindings/Coordinates.h +++ b/projects/detector/private/pybindings/Coordinates.h @@ -13,19 +13,39 @@ void register_Coordinates(pybind11::module_ & m) { using namespace LI::math; using namespace LI::detector; - // Register only `DetectorPosition` and `DetectorDirection` - // The python user does not need to know about or manipulate the Geometry coordinates + // Register the NamedTypes class_>(m, "DetectorPosition") .def(init()) + .def("get", ( + Vector3D & (DetectorPosition::*)() + )(&DetectorPosition::get)) ; class_>(m, "DetectorDirection") .def(init()) + .def("get", ( + Vector3D & (DetectorDirection::*)() + )(&DetectorDirection::get)) + ; + + class_>(m, "GeometryPosition") + .def(init()) + .def("get", ( + Vector3D & (GeometryPosition::*)() + )(&GeometryPosition::get)) + ; + + class_>(m, "GeometryDirection") + .def(init()) + .def("get", ( + Vector3D & (GeometryDirection::*)() + )(&GeometryDirection::get)) ; // Allow implicit conversion from Vector3D to DetectorPosition and DetectorDirection // But only on the python side! + // This works because the interface exposed to python does not include methods that use GeometryPosition and GeometryDirection implicitly_convertible(); implicitly_convertible(); } diff --git a/projects/detector/private/pybindings/DetectorModel.h b/projects/detector/private/pybindings/DetectorModel.h index 0702c062d..dea1c1516 100644 --- a/projects/detector/private/pybindings/DetectorModel.h +++ b/projects/detector/private/pybindings/DetectorModel.h @@ -216,6 +216,18 @@ void register_DetectorModel(pybind11::module_ & m) { )(&DetectorModel::GetAvailableTargets)) .def("GetTargetMass", &DetectorModel::GetTargetMass) .def("GetMaterials",&DetectorModel::GetMaterials) - .def("GetDetectorOrigin",&DetectorModel::GetDetectorOrigin); - + .def("GetDetectorOrigin",&DetectorModel::GetDetectorOrigin) + .def("GeoPositionToDetPosition", ( + DetectorPosition (DetectorModel::*)(GeometryPosition const &) const + )(&DetectorModel::ToDet)) + .def("GeoDirectionToDetDirection", ( + DetectorDirection (DetectorModel::*)(GeometryDirection const &) const + )(&DetectorModel::ToDet)) + .def("DetPositionToGeoPosition", ( + GeometryPosition (DetectorModel::*)(DetectorPosition const &) const + )(&DetectorModel::ToGeo)) + .def("DetDirectionToGeoDirection", ( + GeometryDirection (DetectorModel::*)(DetectorDirection const &) const + )(&DetectorModel::ToGeo)) + ; } diff --git a/projects/detector/private/pybindings/Path.h b/projects/detector/private/pybindings/Path.h index a600c1296..ec1957907 100644 --- a/projects/detector/private/pybindings/Path.h +++ b/projects/detector/private/pybindings/Path.h @@ -11,6 +11,7 @@ void register_Path(pybind11::module_ & m) { using namespace pybind11; using namespace LI::detector; + using namespace LI::math; class_>(m, "Path") .def(init>()) @@ -23,9 +24,9 @@ void register_Path(pybind11::module_ & m) { .def("HasColumnDepth", &Path::HasColumnDepth) .def("GetDetectorModel", &Path::GetDetectorModel) - .def("GetFirstPoint", &Path::GetFirstPoint) - .def("GetLastPoint", &Path::GetLastPoint) - .def("GetDirection", &Path::GetDirection) + .def("GetFirstPoint", [](Path & p)->Vector3D{return p.GetFirstPoint().get(); }) + .def("GetLastPoint", [](Path & p)->Vector3D{return p.GetLastPoint().get(); }) + .def("GetDirection", [](Path & p)->Vector3D{return p.GetDirection().get(); }) .def("GetGeoFirstPoint", &Path::GetGeoFirstPoint) .def("GetGeoLastPoint", &Path::GetGeoLastPoint) .def("GetGeoDirection", &Path::GetGeoDirection) From 3e3c6e43f6754b18c79252dc2a065832bf376e16 Mon Sep 17 00:00:00 2001 From: Austin Schneider Date: Wed, 17 Jan 2024 18:54:01 -0700 Subject: [PATCH 29/42] Remove MultiplicableTo in an attempt to make operator overloads unambigous to gcc --- .../LeptonInjector/detector/Coordinates.h | 20 ++++++------------- 1 file changed, 6 insertions(+), 14 deletions(-) diff --git a/projects/detector/public/LeptonInjector/detector/Coordinates.h b/projects/detector/public/LeptonInjector/detector/Coordinates.h index 5efb25564..0f9d18d8c 100644 --- a/projects/detector/public/LeptonInjector/detector/Coordinates.h +++ b/projects/detector/public/LeptonInjector/detector/Coordinates.h @@ -2,13 +2,14 @@ #ifndef LI_Coordinates_H #define LI_Coordinates_H +#include #include "LeptonInjector/math/Vector3D.h" // for Vector3D #include namespace fluent { template -struct AltArithmetic { +struct A { template struct MultiplicableBy : crtp @@ -38,15 +39,6 @@ struct DivisibleBy : crtp } }; -template -struct MultiplicableTo : crtp -{ - FLUENT_NODISCARD constexpr O operator*(T const& other) const - { - return O(this->underlying().get() * other.get()); - } -}; - }; } // namespace fluent @@ -55,10 +47,10 @@ struct MultiplicableTo : crtp namespace LI { namespace detector { -using GeometryPosition = fluent::NamedType::MultiplicableBy, fluent::AltArithmetic::DivisibleBy, fluent::AltArithmetic::MultiplicableTo>; -using GeometryDirection = fluent::NamedType::MultiplicableBy, fluent::AltArithmetic::DivisibleBy, fluent::AltArithmetic::MultiplicableTo>; -using DetectorPosition = fluent::NamedType::MultiplicableBy, fluent::AltArithmetic::DivisibleBy, fluent::AltArithmetic::MultiplicableTo>; -using DetectorDirection = fluent::NamedType::MultiplicableBy, fluent::AltArithmetic::DivisibleBy, fluent::AltArithmetic::MultiplicableTo>; +using GeometryPosition = fluent::NamedType::MultiplicableBy, fluent::A::DivisibleBy>; +using GeometryDirection = fluent::NamedType::MultiplicableBy, fluent::A::DivisibleBy>; +using DetectorPosition = fluent::NamedType::MultiplicableBy, fluent::A::DivisibleBy>; +using DetectorDirection = fluent::NamedType::MultiplicableBy, fluent::A::DivisibleBy>; } // namespace detector } // namespace LI From e43bc10363ba63899ca4dbdead822cdd3477672a Mon Sep 17 00:00:00 2001 From: Austin Schneider Date: Thu, 18 Jan 2024 16:33:51 -0700 Subject: [PATCH 30/42] Missing coordinate tranform in volume sampling --- .../primary/vertex/CylinderVolumePositionDistribution.cxx | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/projects/distributions/private/primary/vertex/CylinderVolumePositionDistribution.cxx b/projects/distributions/private/primary/vertex/CylinderVolumePositionDistribution.cxx index aa393cb59..0ff9db99a 100644 --- a/projects/distributions/private/primary/vertex/CylinderVolumePositionDistribution.cxx +++ b/projects/distributions/private/primary/vertex/CylinderVolumePositionDistribution.cxx @@ -33,7 +33,7 @@ LI::math::Vector3D CylinderVolumePositionDistribution::SamplePosition(std::share } double CylinderVolumePositionDistribution::GenerationProbability(std::shared_ptr detector_model, std::shared_ptr interactions, LI::dataclasses::InteractionRecord const & record) const { - LI::math::Vector3D pos(record.interaction_vertex); + LI::math::Vector3D pos(cylinder.GlobalToLocalPosition(record.interaction_vertex)); double z = pos.GetZ(); double r = sqrt(pos.GetX() * pos.GetX() + pos.GetY() * pos.GetY()); if(abs(z) >= 0.5 * cylinder.GetZ() From f5a6020fa3bc52dd86037a29c147f201d35dcd81 Mon Sep 17 00:00:00 2001 From: Austin Schneider Date: Thu, 18 Jan 2024 18:36:35 -0700 Subject: [PATCH 31/42] Compute Q2 more robustly before checking bounds --- projects/interactions/private/DISFromSpline.cxx | 9 ++++++--- projects/interactions/private/pybindings/DISFromSpline.h | 2 +- .../public/LeptonInjector/interactions/DISFromSpline.h | 2 +- 3 files changed, 8 insertions(+), 5 deletions(-) diff --git a/projects/interactions/private/DISFromSpline.cxx b/projects/interactions/private/DISFromSpline.cxx index 67fe38c0d..bc393cb15 100644 --- a/projects/interactions/private/DISFromSpline.cxx +++ b/projects/interactions/private/DISFromSpline.cxx @@ -322,10 +322,10 @@ double DISFromSpline::DifferentialCrossSection(dataclasses::InteractionRecord co double x = Q2 / (2.0 * p2.dot(q)); double lepton_mass = particleMass(interaction.signature.secondary_types[lepton_index]); - return DifferentialCrossSection(primary_energy, x, y, lepton_mass); + return DifferentialCrossSection(primary_energy, x, y, lepton_mass, Q2); } -double DISFromSpline::DifferentialCrossSection(double energy, double x, double y, double secondary_lepton_mass) const { +double DISFromSpline::DifferentialCrossSection(double energy, double x, double y, double secondary_lepton_mass, double Q2) const { double log_energy = log10(energy); // check preconditions if(log_energy < differential_cross_section_.lower_extent(0) @@ -339,7 +339,9 @@ double DISFromSpline::DifferentialCrossSection(double energy, double x, double y // we assume that: // the target is stationary so its energy is just its mass // the incoming neutrino is massless, so its kinetic energy is its total energy - double Q2 = 2.0 * energy * target_mass_ * x * y; + if(std::isnan(Q2)) { + Q2 = 2.0 * energy * target_mass_ * x * y; + } if(Q2 < minimum_Q2_) // cross section not calculated, assumed to be zero return 0; @@ -402,6 +404,7 @@ void DISFromSpline::SampleFinalState(dataclasses::InteractionRecord& interaction double E1_lab = p1_lab.e(); double E2_lab = p2_lab.e(); + // The out-going particle always gets at least enough energy for its rest mass double yMax = 1 - m / primary_energy; double logYMax = log10(yMax); diff --git a/projects/interactions/private/pybindings/DISFromSpline.h b/projects/interactions/private/pybindings/DISFromSpline.h index db6965985..566d9edfc 100644 --- a/projects/interactions/private/pybindings/DISFromSpline.h +++ b/projects/interactions/private/pybindings/DISFromSpline.h @@ -33,7 +33,7 @@ void register_DISFromSpline(pybind11::module_ & m) { .def("TotalCrossSection",overload_cast(&DISFromSpline::TotalCrossSection, const_)) .def("TotalCrossSection",overload_cast(&DISFromSpline::TotalCrossSection, const_)) .def("DifferentialCrossSection",overload_cast(&DISFromSpline::DifferentialCrossSection, const_)) - .def("DifferentialCrossSection",overload_cast(&DISFromSpline::DifferentialCrossSection, const_)) + .def("DifferentialCrossSection",overload_cast(&DISFromSpline::DifferentialCrossSection, const_)) .def("InteractionThreshold",&DISFromSpline::InteractionThreshold) .def("GetPossibleTargets",&DISFromSpline::GetPossibleTargets) .def("GetPossibleTargetsFromPrimary",&DISFromSpline::GetPossibleTargetsFromPrimary) diff --git a/projects/interactions/public/LeptonInjector/interactions/DISFromSpline.h b/projects/interactions/public/LeptonInjector/interactions/DISFromSpline.h index 569568865..9d0995601 100644 --- a/projects/interactions/public/LeptonInjector/interactions/DISFromSpline.h +++ b/projects/interactions/public/LeptonInjector/interactions/DISFromSpline.h @@ -69,7 +69,7 @@ friend cereal::access; double TotalCrossSection(LI::dataclasses::Particle::ParticleType primary, double energy) const; double TotalCrossSection(LI::dataclasses::Particle::ParticleType primary, double energy, LI::dataclasses::Particle::ParticleType target) const override; double DifferentialCrossSection(dataclasses::InteractionRecord const &) const override; - double DifferentialCrossSection(double energy, double x, double y, double secondary_lepton_mass) const; + double DifferentialCrossSection(double energy, double x, double y, double secondary_lepton_mass, double Q2=std::numeric_limits::quiet_NaN()) const; double InteractionThreshold(dataclasses::InteractionRecord const &) const override; void SampleFinalState(dataclasses::InteractionRecord &, std::shared_ptr random) const override; From e63847ff8bb3a01767c6f7ad646f71f9817e673e Mon Sep 17 00:00:00 2001 From: Austin Schneider Date: Thu, 18 Jan 2024 18:37:07 -0700 Subject: [PATCH 32/42] Fix namespace of ostream operator --- .../public/LeptonInjector/dataclasses/InteractionRecord.h | 6 +++++- 1 file changed, 5 insertions(+), 1 deletion(-) diff --git a/projects/dataclasses/public/LeptonInjector/dataclasses/InteractionRecord.h b/projects/dataclasses/public/LeptonInjector/dataclasses/InteractionRecord.h index 7010dccb7..ae5327692 100644 --- a/projects/dataclasses/public/LeptonInjector/dataclasses/InteractionRecord.h +++ b/projects/dataclasses/public/LeptonInjector/dataclasses/InteractionRecord.h @@ -21,6 +21,10 @@ #include "LeptonInjector/dataclasses/InteractionSignature.h" +namespace LI { namespace dataclasses { struct InteractionRecord; } } + +std::ostream& operator<<(std::ostream& os, LI::dataclasses::InteractionRecord const& record); + namespace LI { namespace dataclasses { @@ -38,7 +42,7 @@ struct InteractionRecord { std::vector secondary_helicity; std::vector interaction_parameters; bool operator==(InteractionRecord const & other) const; - friend std::ostream& operator<<(std::ostream& os, InteractionRecord const& record); + friend std::ostream& ::operator<<(std::ostream& os, InteractionRecord const& record); template void serialize(Archive & archive, std::uint32_t const version) { if(version == 0) { From 3559bd776d5ff6ccb0f09497da61ecfabdc531af Mon Sep 17 00:00:00 2001 From: Austin Schneider Date: Sun, 21 Jan 2024 17:06:19 -0700 Subject: [PATCH 33/42] Add physically normalized property to primary energy distributions in pybindinds --- projects/distributions/private/pybindings/distributions.cxx | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/projects/distributions/private/pybindings/distributions.cxx b/projects/distributions/private/pybindings/distributions.cxx index cff0cd07c..1b122aaa2 100644 --- a/projects/distributions/private/pybindings/distributions.cxx +++ b/projects/distributions/private/pybindings/distributions.cxx @@ -81,7 +81,7 @@ PYBIND11_MODULE(distributions,m) { // Energy distributions - class_, InjectionDistribution>(m, "PrimaryEnergyDistribution") + class_, InjectionDistribution, PhysicallyNormalizedDistribution>(m, "PrimaryEnergyDistribution") .def("Sample",&PrimaryEnergyDistribution::Sample); class_, PrimaryEnergyDistribution>(m, "Monoenergetic") From d83f13b5e9ffa9b0df7267aa5eb0e44d175291ab Mon Sep 17 00:00:00 2001 From: Austin Schneider Date: Sun, 21 Jan 2024 17:06:40 -0700 Subject: [PATCH 34/42] Throw when adding duplicate distributions --- projects/injection/private/Process.cxx | 6 ++++-- 1 file changed, 4 insertions(+), 2 deletions(-) diff --git a/projects/injection/private/Process.cxx b/projects/injection/private/Process.cxx index ee2e0d27d..8a0a490d8 100644 --- a/projects/injection/private/Process.cxx +++ b/projects/injection/private/Process.cxx @@ -79,7 +79,8 @@ PhysicalProcess & PhysicalProcess::operator=(PhysicalProcess && other) { void PhysicalProcess::AddPhysicalDistribution(std::shared_ptr dist) { for(auto _dist: physical_distributions) { - if((*_dist) == (*dist)) return; + if((*_dist) == (*dist)) + throw std::runtime_error("Cannot add duplicate WeightableDistributions"); } physical_distributions.push_back(dist); } @@ -111,7 +112,8 @@ void InjectionProcess::AddPhysicalDistribution(std::shared_ptr dist) { for(auto _dist: injection_distributions) { - if((*_dist) == (*dist)) return; + if((*_dist) == (*dist)) + throw std::runtime_error("Cannot add duplicate InjectionDistributions"); } injection_distributions.push_back(dist); physical_distributions.push_back(std::static_pointer_cast(dist)); From 384055dabdcd59de30a22e1c900b148d9fc171bf Mon Sep 17 00:00:00 2001 From: Austin Schneider Date: Sun, 21 Jan 2024 17:09:07 -0700 Subject: [PATCH 35/42] Fix how weighters are cancelled so that iterators are not invalidated --- projects/injection/private/TreeWeighter.cxx | 332 ++++++++++---------- 1 file changed, 165 insertions(+), 167 deletions(-) diff --git a/projects/injection/private/TreeWeighter.cxx b/projects/injection/private/TreeWeighter.cxx index 5380afb45..41072d433 100644 --- a/projects/injection/private/TreeWeighter.cxx +++ b/projects/injection/private/TreeWeighter.cxx @@ -44,145 +44,146 @@ using detector::DetectorPosition; using detector::DetectorDirection; namespace { - template - typename std::iterator_traits::value_type accumulate(InIt begin, InIt end) { - typedef typename std::iterator_traits::value_type real; - real sum = real(0); - real running_error = real(0); - real temp; - real difference; - - for (; begin != end; ++begin) { - difference = *begin; - difference -= running_error; - temp = sum; - temp += difference; - running_error = temp; - running_error -= sum; - running_error -= difference; - sum = std::move(temp); - } - return sum; +template +typename std::iterator_traits::value_type accumulate(InIt begin, InIt end) { + typedef typename std::iterator_traits::value_type real; + real sum = real(0); + real running_error = real(0); + real temp; + real difference; + + for (; begin != end; ++begin) { + difference = *begin; + difference -= running_error; + temp = sum; + temp += difference; + running_error = temp; + running_error -= sum; + running_error -= difference; + sum = std::move(temp); } + return sum; +} - template - T accumulate(std::initializer_list list) { - return accumulate(list.begin(), list.end()); - } +template +T accumulate(std::initializer_list list) { + return accumulate(list.begin(), list.end()); +} - double one_minus_exp_of_negative(double x) { - if(x < 1e-1) { - return std::exp(std::log(x) - x/2.0 + x*x/24.0 - x*x*x*x/2880.0); - } else { - return 1.0 - std::exp(-x); - } +double one_minus_exp_of_negative(double x) { + if(x < 1e-1) { + return std::exp(std::log(x) - x/2.0 + x*x/24.0 - x*x*x*x/2880.0); + } else { + return 1.0 - std::exp(-x); } - double log_one_minus_exp_of_negative(double x) { - if(x < 1e-1) { - return std::log(x) - x/2.0 + x*x/24.0 - x*x*x*x/2880.0; - } else if(x > 3) { - double ex = std::exp(-x); - double ex2 = ex * ex; - double ex3 = ex2 * ex; - double ex4 = ex3 * ex; - double ex5 = ex4 * ex; - double ex6 = ex5 * ex; - return -(ex + ex2 / 2.0 + ex3 / 3.0 + ex4 / 4.0 + ex5 / 5.0 + ex6 / 6.0); - } else { - return std::log(1.0 - std::exp(-x)); - } +} + +double log_one_minus_exp_of_negative(double x) { + if(x < 1e-1) { + return std::log(x) - x/2.0 + x*x/24.0 - x*x*x*x/2880.0; + } else if(x > 3) { + double ex = std::exp(-x); + double ex2 = ex * ex; + double ex3 = ex2 * ex; + double ex4 = ex3 * ex; + double ex5 = ex4 * ex; + double ex6 = ex5 * ex; + return -(ex + ex2 / 2.0 + ex3 / 3.0 + ex4 / 4.0 + ex5 / 5.0 + ex6 / 6.0); + } else { + return std::log(1.0 - std::exp(-x)); } } +} // namespace //--------------- // class LeptonTreeWeighter //--------------- void LeptonTreeWeighter::Initialize() { - int i = 0; - primary_process_weighters.reserve(injectors.size()); - secondary_process_weighter_maps.reserve(injectors.size()); - for(auto const & injector : injectors) { - assert(primary_physical_process->MatchesHead(injector->GetPrimaryProcess())); - primary_process_weighters.push_back(std::make_shared(LeptonProcessWeighter(primary_physical_process,injector->GetPrimaryProcess(),detector_model))); - std::map - > injector_sec_process_weighter_map; - std::map - > injector_sec_process_map = injector->GetSecondaryProcessMap(); - for(auto const & sec_phys_process : secondary_physical_processes) { - try{ - std::shared_ptr sec_inj_process = injector_sec_process_map.at(sec_phys_process->GetPrimaryType()); - assert(sec_phys_process->MatchesHead(sec_inj_process)); // make sure cross section collection matches - injector_sec_process_weighter_map[sec_phys_process->GetPrimaryType()] = std::make_shared(LeptonProcessWeighter(sec_phys_process,sec_inj_process,detector_model)); - } catch(const std::out_of_range& oor) { - std::cout << "Out of Range error: " << oor.what() << '\n'; - std::cout << "Initialization Incomplete: Particle " << sec_phys_process->GetPrimaryType() << " does not exist in injector\n"; - return; - } - } - if(injector_sec_process_weighter_map.size() != injector_sec_process_map.size()) { - std::cout << "Initialization Incomplete: No one-to-one mapping between injection and physical distributions for injector " << i << "\n"; - return; + int i = 0; + primary_process_weighters.reserve(injectors.size()); + secondary_process_weighter_maps.reserve(injectors.size()); + for(auto const & injector : injectors) { + assert(primary_physical_process->MatchesHead(injector->GetPrimaryProcess())); + primary_process_weighters.push_back(std::make_shared(LeptonProcessWeighter(primary_physical_process, injector->GetPrimaryProcess(), detector_model))); + std::map + > injector_sec_process_weighter_map; + std::map + > injector_sec_process_map = injector->GetSecondaryProcessMap(); + for(auto const & sec_phys_process : secondary_physical_processes) { + try{ + std::shared_ptr sec_inj_process = injector_sec_process_map.at(sec_phys_process->GetPrimaryType()); + assert(sec_phys_process->MatchesHead(sec_inj_process)); // make sure cross section collection matches + injector_sec_process_weighter_map[sec_phys_process->GetPrimaryType()] = std::make_shared(LeptonProcessWeighter(sec_phys_process,sec_inj_process,detector_model)); + } catch(const std::out_of_range& oor) { + std::cout << "Out of Range error: " << oor.what() << '\n'; + std::cout << "Initialization Incomplete: Particle " << sec_phys_process->GetPrimaryType() << " does not exist in injector\n"; + return; + } + } + if(injector_sec_process_weighter_map.size() != injector_sec_process_map.size()) { + std::cout << "Initialization Incomplete: No one-to-one mapping between injection and physical distributions for injector " << i << "\n"; + return; + } + secondary_process_weighter_maps.push_back(injector_sec_process_weighter_map); } - secondary_process_weighter_maps.push_back(injector_sec_process_weighter_map); - } } double LeptonTreeWeighter::EventWeight(LI::dataclasses::InteractionTree const & tree) const { - // The weight is given by - // - // [sum_{injectors i} - // x prod_{tree datum d} - // x (prod_{generation dist j} p_gen^{idj}) - // / (prod_{physical dist j} p_phys^{idj}) ] ^-1 - // - double inv_weight = 0; - for(unsigned int idx = 0; idx < injectors.size(); ++idx) { - double physical_probability = 1.0; - double generation_probability = injectors[idx]->EventsToInject();//GenerationProbability(tree); - for(auto const & datum : tree.tree) { - std::pair bounds; - if(datum->depth()==0) { - bounds = injectors[idx]->InjectionBounds(datum->record); - physical_probability *= primary_process_weighters[idx]->PhysicalProbability(bounds, datum->record); - generation_probability *= primary_process_weighters[idx]->GenerationProbability(*datum); - } - else { - try { - bounds = injectors[idx]->InjectionBounds(*datum, datum->record.signature.primary_type); - double phys_prob = secondary_process_weighter_maps[idx].at(datum->record.signature.primary_type)->PhysicalProbability(bounds, datum->record); - double gen_prob = secondary_process_weighter_maps[idx].at(datum->record.signature.primary_type)->GenerationProbability(*datum); - physical_probability *= phys_prob; - generation_probability *= gen_prob; - } catch(const std::out_of_range& oor) { - std::cout << "Out of Range error: " << oor.what() << '\n'; - return 0; + // The weight is given by + // + // [sum_{injectors i} + // x prod_{tree datum d} + // x (prod_{generation dist j} p_gen^{idj}) + // / (prod_{physical dist j} p_phys^{idj}) ] ^-1 + // + double inv_weight = 0; + for(unsigned int idx = 0; idx < injectors.size(); ++idx) { + double physical_probability = 1.0; + double generation_probability = injectors[idx]->EventsToInject();//GenerationProbability(tree); + for(auto const & datum : tree.tree) { + std::pair bounds; + if(datum->depth()==0) { + bounds = injectors[idx]->InjectionBounds(datum->record); + physical_probability *= primary_process_weighters[idx]->PhysicalProbability(bounds, datum->record); + generation_probability *= primary_process_weighters[idx]->GenerationProbability(*datum); + } + else { + try { + bounds = injectors[idx]->InjectionBounds(*datum, datum->record.signature.primary_type); + double phys_prob = secondary_process_weighter_maps[idx].at(datum->record.signature.primary_type)->PhysicalProbability(bounds, datum->record); + double gen_prob = secondary_process_weighter_maps[idx].at(datum->record.signature.primary_type)->GenerationProbability(*datum); + physical_probability *= phys_prob; + generation_probability *= gen_prob; + } catch(const std::out_of_range& oor) { + std::cout << "Out of Range error: " << oor.what() << '\n'; + return 0; + } + } } - } + inv_weight += generation_probability / physical_probability; } - inv_weight += generation_probability / physical_probability; - } - return 1./inv_weight; + return 1./inv_weight; } LeptonTreeWeighter::LeptonTreeWeighter(std::vector> injectors, std::shared_ptr detector_model, std::shared_ptr primary_physical_process, std::vector> secondary_physical_processes) : injectors(injectors) - , detector_model(detector_model) - , primary_physical_process(primary_physical_process) - , secondary_physical_processes(secondary_physical_processes) + , detector_model(detector_model) + , primary_physical_process(primary_physical_process) + , secondary_physical_processes(secondary_physical_processes) { - Initialize(); + Initialize(); } LeptonTreeWeighter::LeptonTreeWeighter(std::vector> injectors, std::shared_ptr detector_model, std::shared_ptr primary_physical_process) : injectors(injectors) - , detector_model(detector_model) - , primary_physical_process(primary_physical_process) - , secondary_physical_processes(std::vector>()) + , detector_model(detector_model) + , primary_physical_process(primary_physical_process) + , secondary_physical_processes(std::vector>()) { - Initialize(); + Initialize(); } //--------------- @@ -190,30 +191,32 @@ LeptonTreeWeighter::LeptonTreeWeighter(std::vector> in //--------------- void LeptonProcessWeighter::Initialize() { - normalization = 1.0; - for(auto physical_dist : phys_process->GetPhysicalDistributions()) { - const LI::distributions::PhysicallyNormalizedDistribution* p = dynamic_cast(physical_dist.get()); - if(p) { - if(p->IsNormalizationSet()) { - normalization *= p->GetNormalization(); - } + normalization = 1.0; + for(auto physical_dist : phys_process->GetPhysicalDistributions()) { + const LI::distributions::PhysicallyNormalizedDistribution* p = dynamic_cast(physical_dist.get()); + if(p) { + if(p->IsNormalizationSet()) { + normalization *= p->GetNormalization(); + } + } } - } - unique_gen_distributions = inj_process->GetInjectionDistributions(); - unique_phys_distributions = phys_process->GetPhysicalDistributions(); - std::vector>::iterator gen_iterator = unique_gen_distributions.begin(); - while(gen_iterator != unique_gen_distributions.end()) { - std::vector>::iterator phys_iterator = unique_phys_distributions.begin(); - while(phys_iterator != unique_phys_distributions.end()) { - if((*gen_iterator) == (*phys_iterator)) { - unique_gen_distributions.erase(gen_iterator); - unique_phys_distributions.erase(phys_iterator); - break; - } - ++phys_iterator; + unique_gen_distributions = inj_process->GetInjectionDistributions(); + unique_phys_distributions = phys_process->GetPhysicalDistributions(); + std::cout << "Num gen distributions: " << unique_gen_distributions.size() << std::endl; + std::cout << "Num phys distributions: " << unique_phys_distributions.size() << std::endl; + for(std::vector>::reverse_iterator gen_it = unique_gen_distributions.rbegin(); + gen_it != unique_gen_distributions.rend(); ++gen_it) { + for(std::vector>::reverse_iterator phys_it = unique_phys_distributions.rbegin(); + phys_it != unique_phys_distributions.rend(); ++phys_it) { + if((*gen_it) == (*phys_it)) { + unique_gen_distributions.erase(std::next(gen_it).base()); + unique_phys_distributions.erase(std::next(phys_it).base()); + break; + } + } } - ++gen_iterator; - } + std::cout << "Num gen distributions: " << unique_gen_distributions.size() << std::endl; + std::cout << "Num phys distributions: " << unique_phys_distributions.size() << std::endl; } double LeptonProcessWeighter::InteractionProbability(std::pair const & bounds, LI::dataclasses::InteractionRecord const & record) const { @@ -234,7 +237,7 @@ double LeptonProcessWeighter::InteractionProbability(std::pair total_cross_sections; double total_decay_length = phys_process->GetInteractions()->TotalDecayLength(record); - + LI::dataclasses::InteractionRecord fake_record = record; for(auto const & target_xs : cross_sections_by_target) { targets.push_back(target_xs.first); @@ -254,7 +257,7 @@ double LeptonProcessWeighter::InteractionProbability(std::pairGetInteractionDepthInCGS(intersections, DetectorPosition(bounds.first), DetectorPosition(bounds.second), targets, total_cross_sections, total_decay_length); - + double interaction_probability; if(total_interaction_depth < 1e-6) { interaction_probability = total_interaction_depth; @@ -305,7 +308,6 @@ double LeptonProcessWeighter::NormalizedPositionProbability(std::pairGetInteractionDepthInCGS(intersections, DetectorPosition(bounds.first), DetectorPosition(bounds.second), targets, total_cross_sections, total_decay_length); // unitless double traversed_interaction_depth = detector_model->GetInteractionDepthInCGS(intersections, DetectorPosition(bounds.first), DetectorPosition(interaction_vertex), targets, total_cross_sections, total_decay_length); double interaction_density = detector_model->GetInteractionDensity(intersections, DetectorPosition(interaction_vertex), targets, total_cross_sections, total_decay_length); //units of m^-1 - double prob_density; if(total_interaction_depth < 1e-6) { @@ -318,49 +320,45 @@ double LeptonProcessWeighter::NormalizedPositionProbability(std::pair const & bounds, - LI::dataclasses::InteractionRecord const & record ) const { - - - double physical_probability = 1.0; - double prob = InteractionProbability(bounds, record); - physical_probability *= prob; - - - prob = NormalizedPositionProbability(bounds, record); - physical_probability *= prob; - - - prob = LI::injection::CrossSectionProbability(detector_model, phys_process->GetInteractions(), record); - physical_probability *= prob; - - + LI::dataclasses::InteractionRecord const & record ) const { + + double physical_probability = 1.0; + double prob = InteractionProbability(bounds, record); + physical_probability *= prob; + + prob = NormalizedPositionProbability(bounds, record); + physical_probability *= prob; + + prob = LI::injection::CrossSectionProbability(detector_model, phys_process->GetInteractions(), record); + physical_probability *= prob; + for(auto physical_dist : unique_phys_distributions) { - physical_probability *= physical_dist->GenerationProbability(detector_model, phys_process->GetInteractions(), record); - } - - return normalization * physical_probability; + physical_probability *= physical_dist->GenerationProbability(detector_model, phys_process->GetInteractions(), record); + } + + return normalization * physical_probability; } double LeptonProcessWeighter::GenerationProbability(LI::dataclasses::InteractionTreeDatum const & datum ) const { - double gen_probability = LI::injection::CrossSectionProbability(detector_model, phys_process->GetInteractions(), datum.record); - - for(auto gen_dist : unique_gen_distributions) { - gen_probability *= gen_dist->GenerationProbability(detector_model, phys_process->GetInteractions(), datum); - } - return gen_probability; + double gen_probability = LI::injection::CrossSectionProbability(detector_model, phys_process->GetInteractions(), datum.record); + + for(auto gen_dist : unique_gen_distributions) { + gen_probability *= gen_dist->GenerationProbability(detector_model, phys_process->GetInteractions(), datum); + } + return gen_probability; } double LeptonProcessWeighter::EventWeight(std::pair const & bounds, - LI::dataclasses::InteractionTreeDatum const & datum) const { - return PhysicalProbability(bounds,datum.record)/GenerationProbability(datum); + LI::dataclasses::InteractionTreeDatum const & datum) const { + return PhysicalProbability(bounds,datum.record)/GenerationProbability(datum); } LeptonProcessWeighter::LeptonProcessWeighter(std::shared_ptr phys_process, std::shared_ptr inj_process, std::shared_ptr detector_model) : phys_process(phys_process) - , inj_process(inj_process) - , detector_model(detector_model) + , inj_process(inj_process) + , detector_model(detector_model) { - Initialize(); + Initialize(); } } // namespace injection From cb37c67489532cad5bb0a913073c5ea4ca882976 Mon Sep 17 00:00:00 2001 From: Nicholas Kamp Date: Sat, 13 Jan 2024 23:06:15 -0500 Subject: [PATCH 36/42] updated tabulated flux distribution to handle zeros in input flux file and properly integrate by trapezoidal rule --- .../energy/TabulatedFluxDistribution.cxx | 43 +++++++++++++------ 1 file changed, 30 insertions(+), 13 deletions(-) diff --git a/projects/distributions/private/primary/energy/TabulatedFluxDistribution.cxx b/projects/distributions/private/primary/energy/TabulatedFluxDistribution.cxx index 9709d9cf6..4acfd8ccc 100644 --- a/projects/distributions/private/primary/energy/TabulatedFluxDistribution.cxx +++ b/projects/distributions/private/primary/energy/TabulatedFluxDistribution.cxx @@ -112,36 +112,53 @@ void TabulatedFluxDistribution::ComputeCDF() { energies.push_back(energyMax); // assign the energies vector so it's accessible outside of the function - cdf_energy_nodes = energies; - - // declare the cdf vector - std::vector cdfvector(energies.size()); + //cdf_energy_nodes = energies; + + // declare the cdf vectors + std::vector cdf_vector; + std::vector cdf_energy_nodes; + + cdf_vector.push_back(0); // start the CDF at zero + cdf_energy_nodes.push_back(energies[0]); + for(int i = 1; i < energies.size(); ++i) { + double p1 = pdf(energies[i-1]); + double p2 = pdf(energies[i]); + + // skip unsupported parts of pdf + if((p1+p2)<=0) continue; + + // Check if we have skipped previous energy node + if(cdf_energy_nodes.back() != energies[i-1]){ + // if so, we need to add small epsilon to CDF to + // ignore unsupported part of PDF during interpolation + cdf_energy_nodes.push_back(energies[i-1]); + cdf_vector.push_back(cdf_vector.back()+1e-12); + } - cdfvector[0] = 0; - for (std::size_t i = 1; i < energies.size(); ++i) { - double y = pdf(energies[i]); - double area = y * (energies[i]-energies[i-1]); - cdfvector[i] = cdfvector[i-1] + area; + // trapezoidal area + double area = 0.5 * (p1 + p2) * (energies[i]-energies[i-1]); + cdf_vector.push_back(cdf_vector.back() + area); + cdf_energy_nodes.push_back(energies[i]); } // find the max of CDF (should be 1 since we computed from normalized PDF) - auto max_it = std::max_element(cdfvector.begin(), cdfvector.end()); + auto max_it = std::max_element(cdf_vector.begin(), cdf_vector.end()); double max_cdf = *max_it; // should be normalized, but just to make sure in case energy nodes are too sparse - for (double& value : cdfvector) { + for (double& value : cdf_vector) { value *= 1/max_cdf; } // assign the cdf vector so it's accessible outside of the function - cdf = cdfvector; + cdf = cdf_vector; LI::utilities::TableData1D inverse_cdf_data; inverse_cdf_data.x = cdf; - inverse_cdf_data.f = energies; + inverse_cdf_data.f = cdf_energy_nodes; inverseCdfTable = LI::utilities::Interpolator1D(inverse_cdf_data); From 1b4efca2cdad215cb0d2474798e8c365123d1ca6 Mon Sep 17 00:00:00 2001 From: Nicholas Kamp Date: Sat, 13 Jan 2024 23:08:09 -0500 Subject: [PATCH 37/42] clean up names in python code and existing examples --- python/LIDarkNews.py | 5 +++-- resources/Examples/Example1/DIS_IceCube.py | 4 ++-- resources/Examples/Example2/DipolePortal_MINERvA.py | 3 +-- resources/Examples/Example2/DipolePortal_MiniBooNE.py | 4 +--- 4 files changed, 7 insertions(+), 9 deletions(-) diff --git a/python/LIDarkNews.py b/python/LIDarkNews.py index 6ff05750d..cb7012bf4 100644 --- a/python/LIDarkNews.py +++ b/python/LIDarkNews.py @@ -97,7 +97,7 @@ def __init__(self, self.decays.append(PyDarkNewsDecay(dec_case, table_dir = self.table_dir + table_subdirs)) - def SaveCrossSectionTables(self,fill_tables_at_exit=True): + def SaveCrossSectionTables(self,fill_tables_at_exit=False): if not fill_tables_at_exit: print('WARNING: Saving tables without filling PyDarkNewsCrossSection interpolation tables. Future updates to DarkNews can lead to inconsistent behavior if new entries are ever added to this table') for cross_section in self.cross_sections: @@ -198,7 +198,8 @@ def _interpolation_flags(self,inputs,mode): exit(0) # first check if we have saved table points already - if len(interp_table) == 0: return 0 + if len(interp_table) == 0: + return False,False,-1 # bools to keep track of whether to use a single point or interpolate UseSinglePoint = True diff --git a/resources/Examples/Example1/DIS_IceCube.py b/resources/Examples/Example1/DIS_IceCube.py index 924e6d98a..d53dc909a 100644 --- a/resources/Examples/Example1/DIS_IceCube.py +++ b/resources/Examples/Example1/DIS_IceCube.py @@ -25,12 +25,12 @@ xsfiledir = LI_SRC+'resources/CrossSectionTables/DISSplines/' target_type = LI.dataclasses.Particle.ParticleType.Nucleon -DIS_xs = LI.crosssections.DISFromSpline(xsfiledir+'test_xs.fits', +DIS_xs = LI.interactions.DISFromSpline(xsfiledir+'test_xs.fits', xsfiledir+'test_xs_total.fits', [primary_type], [target_type]) -primary_xs = LI.crosssections.CrossSectionCollection(primary_type, [DIS_xs]) +primary_xs = LI.interactions.CrossSectionCollection(primary_type, [DIS_xs]) controller.SetCrossSections(primary_xs) # Primary distributions diff --git a/resources/Examples/Example2/DipolePortal_MINERvA.py b/resources/Examples/Example2/DipolePortal_MINERvA.py index b13f7fbf8..7334fe668 100644 --- a/resources/Examples/Example2/DipolePortal_MINERvA.py +++ b/resources/Examples/Example2/DipolePortal_MINERvA.py @@ -1,4 +1,3 @@ -import numpy as np import os import leptoninjector as LI @@ -67,7 +66,7 @@ 240) position_distribution = LI.distributions.RangePositionDistribution(1.24, 5., decay_range_func, - set(controller.GetEarthModelTargets()[0])) + set(controller.GetDetectorModelTargets()[0])) primary_injection_distributions['position'] = position_distribution # SetProcesses diff --git a/resources/Examples/Example2/DipolePortal_MiniBooNE.py b/resources/Examples/Example2/DipolePortal_MiniBooNE.py index d175b5ab1..d13ede9b9 100644 --- a/resources/Examples/Example2/DipolePortal_MiniBooNE.py +++ b/resources/Examples/Example2/DipolePortal_MiniBooNE.py @@ -1,6 +1,4 @@ -import numpy as np import os -import matplotlib.pyplot as plt import leptoninjector as LI import sys @@ -68,7 +66,7 @@ 541) position_distribution = LI.distributions.RangePositionDistribution(6.2, 6.2, decay_range_func, - set(controller.GetEarthModelTargets()[0])) + set(controller.GetDetectorModelTargets()[0])) primary_injection_distributions['position'] = position_distribution # SetProcesses From 9ab22a0ba77b2402d2a2482eb2fd3fe605a06015 Mon Sep 17 00:00:00 2001 From: Austin Schneider Date: Mon, 22 Jan 2024 12:01:14 -0700 Subject: [PATCH 38/42] Add inf flags to Path (not yet used) --- projects/detector/private/Path.cxx | 24 +++++++++++++++++++ .../public/LeptonInjector/detector/Path.h | 7 ++++++ 2 files changed, 31 insertions(+) diff --git a/projects/detector/private/Path.cxx b/projects/detector/private/Path.cxx index 824835846..824c11860 100644 --- a/projects/detector/private/Path.cxx +++ b/projects/detector/private/Path.cxx @@ -28,6 +28,30 @@ void Path::UpdatePoints() { } } +bool Path::IsInfinite(LI::math::Vector3D const & vec) { + return std::isinf(vec.GetX()) or std::isinf(vec.GetY()) or std::isinf(vec.GetZ()); +} + +void Path::RequireFirstFinite() { + if(first_inf_) + throw std::runtime_error("First point is required to be finite here"); +} + +void Path::RequireLastFinite() { + if(last_inf_) + throw std::runtime_error("Last point is required to be finite here"); +} + +void Path::RequireBothFinite() { + if(first_inf_ or last_inf_) + throw std::runtime_error("Both points are required to be finite here"); +} + +void Path::RequireOneFinite() { + if(first_inf_ and last_inf_) + throw std::runtime_error("At least one point is required to be finite here"); +} + Path::Path() { } diff --git a/projects/detector/public/LeptonInjector/detector/Path.h b/projects/detector/public/LeptonInjector/detector/Path.h index 8944a3531..b047f4ed3 100644 --- a/projects/detector/public/LeptonInjector/detector/Path.h +++ b/projects/detector/public/LeptonInjector/detector/Path.h @@ -39,6 +39,8 @@ class Path { GeometryDirection direction_; double distance_ = 0; bool set_points_ = false; + bool first_inf_ = false; + bool last_inf_ = false; DetectorPosition first_point_det_; DetectorPosition last_point_det_; @@ -52,6 +54,11 @@ class Path { bool set_intersections_ = false; void UpdatePoints(); + static bool IsInfinite(LI::math::Vector3D const & vec); + void RequireFirstFinite(); + void RequireLastFinite(); + void RequireBothFinite(); + void RequireOneFinite(); public: Path(); Path(std::shared_ptr detector_model); From 9bd3020e4a0b8dec3f454c67623278682dcadb7c Mon Sep 17 00:00:00 2001 From: Nicholas Kamp Date: Sat, 13 Jan 2024 23:08:26 -0500 Subject: [PATCH 39/42] starting to implement dipole CCM example --- .../Examples/Example2/DipolePortal_CCM.py | 83 +++++++++++++++++++ 1 file changed, 83 insertions(+) create mode 100644 resources/Examples/Example2/DipolePortal_CCM.py diff --git a/resources/Examples/Example2/DipolePortal_CCM.py b/resources/Examples/Example2/DipolePortal_CCM.py new file mode 100644 index 000000000..0702c2197 --- /dev/null +++ b/resources/Examples/Example2/DipolePortal_CCM.py @@ -0,0 +1,83 @@ +import os +import numpy as np + +import leptoninjector as LI +import sys +sys.path.insert(1,'/n/holylfs05/LABS/arguelles_delgado_lab/Everyone/nkamp/LIV2/sources/LeptonInjector/python') +from LIController import LIController + +LI_SRC = os.environ.get('LEPTONINJECTOR_SRC') + +# Define a DarkNews model +model_kwargs = { + 'm4': 0.02, + 'mu_tr_mu4': 2.5e-6, #1e-6, # GeV^-1 + 'UD4': 0, + 'Umu4': 0, + 'epsilon': 0.0, + 'gD': 0.0, + 'decay_product':'photon', + 'noHC':True, + 'HNLtype':"dirac", +} + +# Number of events to inject +events_to_inject = 1000 + +# Expeirment to run +experiment = 'CCM' + +# Define the controller +controller = LIController(events_to_inject, + experiment) + +# Particle to inject +primary_type = LI.dataclasses.Particle.ParticleType.NuMu + +# Define DarkNews Model +table_dir = LI_SRC+'/resources/CrossSectionTables/DarkNewsTables/Dipole_M%2.2f_mu%2.2e/'%(model_kwargs['m4'],model_kwargs['mu_tr_mu4']) +controller.InputDarkNewsModel(primary_type, + table_dir, + model_kwargs) + +# Primary distributions +primary_injection_distributions = {} +primary_physical_distributions = {} + +# energy distribution +nu_energy = 0.02965 # from pi+ DAR +edist = LI.distributions.Monoenergetic(nu_energy) +primary_injection_distributions['energy'] = edist +primary_physical_distributions['energy'] = edist + +# Flux normalization: +# using the number quoted in 2105.14020, 4.74e9 nu/m^2/s / (6.2e14 POT/s) * 4*pi*20m^2 to get nu/POT +flux_units = LI.distributions.NormalizationConstant(3.76e-2) +primary_physical_distributions['flux_units'] = flux_units + +# direction distribution: cone from lower W target +opening_angle = np.arctan(12/23.); # slightly larger than CCM +lower_target_origin = LI.math.Vector3D(0, 0, -0.241) +detector_origin = LI.math.Vector3D(23, 0, -0.65) +lower_dir = detector_origin - lower_target_origin +lower_dir.normalize() +lower_inj_ddist = LI.distributions.Cone(lower_dir,opening_angle) +phys_ddist = LI.distributions.IsotropicDirection() # truly we are isotropicprimary_injection_distributions['direction'] = direction_distribution +primary_injection_distributions['direction'] = lower_inj_ddist +primary_physical_distributions['direction'] = phys_ddist + +# Position distribution: consider neutrinos from a point source +max_dist = 25 +lower_pos_dist = LI.distributions.PointSourcePositionDistribution(lower_target_origin, max_dist, set(controller.GetDetectorModelTargets()[0])) +primary_injection_distributions['position'] = lower_pos_dist + +# SetProcesses +controller.SetProcesses(primary_type, + primary_injection_distributions, + primary_physical_distributions) + +controller.Initialize() + +events = controller.GenerateEvents() + +controller.SaveEvents('output/CCM_Dipole_M%2.2f_mu%2.2e_example.hdf5'%(model_kwargs['m4'],model_kwargs['mu_tr_mu4'])) \ No newline at end of file From c150ee289b99955dd5d0ca834065332b856184df Mon Sep 17 00:00:00 2001 From: Nicholas Kamp Date: Sun, 14 Jan 2024 22:43:39 -0500 Subject: [PATCH 40/42] fix dark news cross section bugs, reject unphysical cross sections --- .../private/DarkNewsCrossSection.cxx | 25 ++++++++++++++++--- 1 file changed, 22 insertions(+), 3 deletions(-) diff --git a/projects/interactions/private/DarkNewsCrossSection.cxx b/projects/interactions/private/DarkNewsCrossSection.cxx index 67dc82dce..6e27794f0 100644 --- a/projects/interactions/private/DarkNewsCrossSection.cxx +++ b/projects/interactions/private/DarkNewsCrossSection.cxx @@ -172,6 +172,16 @@ void DarkNewsCrossSection::SampleFinalState(dataclasses::InteractionRecord & int double log_minQ2 = log10(minQ2); double log_maxQ2 = log10(maxQ2); + // Function for computing CosTheta in the lab frame + auto ComputeCosThetaLab = [&](double Q2) { + double E4Lab = (Q2 + m2*m2 + m4*m4) / (2*m2); + double E3Lab = primary_energy + p2_lab.e() - E4Lab; + double P1Lab = sqrt(primary_energy*primary_energy - m1*m1); + double P3Lab = sqrt(E3Lab*E3Lab - m3*m3); + double CosThetaLab = (primary_energy*E3Lab - 0.5*(Q2 + m1*m1 + m3*m3))/(P1Lab*P3Lab); + return CosThetaLab; + }; + bool accept; // kin_vars and its twin are 2-vectors containing [nu-energy, Q2] @@ -185,6 +195,10 @@ void DarkNewsCrossSection::SampleFinalState(dataclasses::InteractionRecord & int // sample an intial point kin_vars[1] = std::pow(10,random->Uniform(log_minQ2,log_maxQ2)); + // make sure it is physical + while(std::abs(ComputeCosThetaLab(kin_vars[1]))>1) { + kin_vars[1] = std::pow(10,random->Uniform(log_minQ2,log_maxQ2)); + } test_cross_section = DifferentialCrossSection(interaction.signature.primary_type, interaction.signature.target_type, primary_energy, kin_vars[1]); cross_section = test_cross_section; @@ -196,13 +210,18 @@ void DarkNewsCrossSection::SampleFinalState(dataclasses::InteractionRecord & int for(size_t j = 0; j <= burnin; j++) { // repeat the sampling from above to get a new valid point test_kin_vars[1] = std::pow(10,random->Uniform(log_minQ2,log_maxQ2)); - test_cross_section = DifferentialCrossSection(interaction.signature.primary_type, interaction.signature.target_type, primary_energy, kin_vars[1]); + test_cross_section = DifferentialCrossSection(interaction.signature.primary_type, interaction.signature.target_type, primary_energy, test_kin_vars[1]); double odds = (test_cross_section / cross_section); accept = (cross_section == 0 || (odds > 1.) || random->Uniform(0, 1) < odds); if(accept) { - kin_vars = test_kin_vars; - cross_section = test_cross_section; + // Make sure we have a physical point + double CosThetaLabTest = ComputeCosThetaLab(test_kin_vars[1]); + bool physical = (std::abs(CosThetaLabTest) <= 1); + if(physical) { + kin_vars = test_kin_vars; + cross_section = test_cross_section; + } } } double final_Q2 = kin_vars[1]; From b3a2c21365b37bb8a616a9cecee50891e133731a Mon Sep 17 00:00:00 2001 From: Nicholas Kamp Date: Mon, 22 Jan 2024 17:03:21 -0500 Subject: [PATCH 41/42] Require finite points in Path where necessary, make sure to properly clip infinite points to bounds --- projects/detector/private/Path.cxx | 51 ++++++++++++++++++++++++++++-- 1 file changed, 49 insertions(+), 2 deletions(-) diff --git a/projects/detector/private/Path.cxx b/projects/detector/private/Path.cxx index 824c11860..1098038bc 100644 --- a/projects/detector/private/Path.cxx +++ b/projects/detector/private/Path.cxx @@ -163,6 +163,9 @@ void Path::SetPoints(GeometryPosition first_point, GeometryPosition last_point) set_det_points_ = false; set_intersections_ = false; set_column_depth_ = false; + first_inf_ = IsInfinite(first_point); + last_inf_ = IsInfinite(last_point); + RequireBothFinite(); UpdatePoints(); } @@ -176,6 +179,9 @@ void Path::SetPoints(DetectorPosition first_point, DetectorPosition last_point) set_det_points_ = true; set_intersections_ = false; set_column_depth_ = false; + first_inf_ = IsInfinite(first_point); + last_inf_ = IsInfinite(last_point); + RequireBothFinite(); UpdatePoints(); } @@ -191,6 +197,9 @@ void Path::SetPointsWithRay(GeometryPosition first_point, GeometryDirection dire set_det_points_ = false; set_intersections_ = false; set_column_depth_ = false; + first_inf_ = IsInfinite(first_point_); // Set using GeometryPosition + last_inf_ = IsInfinite(last_point_); // Set using GeometryPosition + RequireFirstFinite(); UpdatePoints(); } @@ -206,6 +215,9 @@ void Path::SetPointsWithRay(DetectorPosition first_point, DetectorDirection dire set_det_points_ = true; set_intersections_ = false; set_column_depth_ = false; + first_inf_ = IsInfinite(first_point_det_); // Set using DetectorPosition + last_inf_ = IsInfinite(last_point_det_); // Set using DetectorPosition + RequireFirstFinite(); UpdatePoints(); } @@ -249,14 +261,16 @@ void Path::ClipToOuterBounds() { if(dot < 0) { p0.swap(p1); } - bool clip_0 = (p0 - first_point_) * direction_ > 0; - bool clip_1 = (p1 - last_point_) * direction_ < 0; + bool clip_0 = ((p0 - first_point_) * direction_ > 0) || (first_inf_); + bool clip_1 = ((p1 - last_point_) * direction_ < 0) || (last_inf_); bool clip = clip_0 or clip_1; if(clip_0) { first_point_ = GeometryPosition(p0); + first_inf_ = IsInfinite(first_point_); // re-check infinite } if(clip_1) { last_point_ = GeometryPosition(p1); + last_inf_ = IsInfinite(last_point_); // re-check infinite } if(clip) { distance_ = (last_point_ - first_point_)->magnitude(); @@ -271,6 +285,7 @@ void Path::ClipToOuterBounds() { void Path::Flip() { std::swap(first_point_, last_point_); std::swap(first_point_det_, last_point_det_); + std::swap(first_inf_, last_inf_); direction_ *= -1; direction_det_ *= -1; } @@ -281,6 +296,7 @@ void Path::Flip() { //// void Path::ExtendFromEndByDistance(double distance) { EnsurePoints(); + RequireLastFinite(); distance_ += distance; last_point_ += direction_ * distance; if(distance_ < 0) { @@ -293,6 +309,7 @@ void Path::ExtendFromEndByDistance(double distance) { void Path::ExtendFromStartByDistance(double distance) { EnsurePoints(); + RequireFirstFinite(); distance_ += distance; first_point_ += direction_ * -distance; if(distance_ < 0) { @@ -486,6 +503,7 @@ void Path::ShrinkFromStartToInteractionDepth(double interaction_depth, double Path::GetColumnDepthInBounds() { EnsureIntersections(); EnsurePoints(); + RequireBothFinite(); if(HasColumnDepth()) { return column_depth_cached_; } else { @@ -505,6 +523,7 @@ double Path::GetInteractionDepthInBounds( double const & total_decay_length) { EnsureIntersections(); EnsurePoints(); + RequireBothFinite(); double interaction_depth = detector_model_->GetInteractionDepthInCGS(intersections_, first_point_, last_point_, targets, total_cross_sections, total_decay_length); return interaction_depth; } @@ -521,6 +540,7 @@ double Path::GetColumnDepthFromStartInBounds(double distance) { } EnsureIntersections(); EnsurePoints(); + RequireFirstFinite(); return detector_model_->GetColumnDepthInCGS(intersections_, first_point_, GeometryPosition(first_point_ + direction_ * distance)); } @@ -532,30 +552,35 @@ double Path::GetColumnDepthFromEndInBounds(double distance) { } EnsureIntersections(); EnsurePoints(); + RequireLastFinite(); return detector_model_->GetColumnDepthInCGS(intersections_, last_point_, GeometryPosition(last_point_ + direction_ * -distance)); } double Path::GetColumnDepthFromStartAlongPath(double distance) { EnsureIntersections(); EnsurePoints(); + RequireFirstFinite(); return std::copysign(detector_model_->GetColumnDepthInCGS(intersections_, first_point_, GeometryPosition(first_point_ + direction_ * distance)), distance); } double Path::GetColumnDepthFromEndAlongPath(double distance) { EnsureIntersections(); EnsurePoints(); + RequireLastFinite(); return std::copysign(detector_model_->GetColumnDepthInCGS(intersections_, last_point_, GeometryPosition(last_point_ + direction_ * distance)), distance); } double Path::GetColumnDepthFromStartInReverse(double distance) { EnsureIntersections(); EnsurePoints(); + RequireFirstFinite(); return std::copysign(detector_model_->GetColumnDepthInCGS(intersections_, first_point_, GeometryPosition(first_point_ + direction_ * -distance)), distance); } double Path::GetColumnDepthFromEndInReverse(double distance) { EnsureIntersections(); EnsurePoints(); + RequireLastFinite(); return std::copysign(detector_model_->GetColumnDepthInCGS(intersections_, last_point_, GeometryPosition(last_point_ + direction_ * -distance)), distance); } @@ -574,6 +599,7 @@ double Path::GetInteractionDepthFromStartInBounds(double distance, } EnsureIntersections(); EnsurePoints(); + RequireFirstFinite(); return detector_model_->GetInteractionDepthInCGS(intersections_, first_point_, GeometryPosition(first_point_ + direction_ * distance), targets, total_cross_sections, total_decay_length); } @@ -588,6 +614,7 @@ double Path::GetInteractionDepthFromEndInBounds(double distance, } EnsureIntersections(); EnsurePoints(); + RequireLastFinite(); return detector_model_->GetInteractionDepthInCGS(intersections_, last_point_, GeometryPosition(last_point_ + direction_ * -distance), targets, total_cross_sections, total_decay_length); } @@ -597,6 +624,7 @@ double Path::GetInteractionDepthFromStartAlongPath(double distance, double const & total_decay_length) { EnsureIntersections(); EnsurePoints(); + RequireFirstFinite(); return std::copysign(detector_model_->GetInteractionDepthInCGS(intersections_, first_point_, GeometryPosition(first_point_ + direction_ * distance), targets, total_cross_sections, total_decay_length), distance); } @@ -606,6 +634,7 @@ double Path::GetInteractionDepthFromEndAlongPath(double distance, double const & total_decay_length) { EnsureIntersections(); EnsurePoints(); + RequireLastFinite(); return std::copysign(detector_model_->GetInteractionDepthInCGS(intersections_, last_point_, GeometryPosition(last_point_ + direction_ * distance), targets, total_cross_sections, total_decay_length), distance); } @@ -615,6 +644,7 @@ double Path::GetInteractionDepthFromStartInReverse(double distance, double const & total_decay_length) { EnsureIntersections(); EnsurePoints(); + RequireFirstFinite(); return std::copysign(detector_model_->GetInteractionDepthInCGS(intersections_, first_point_, GeometryPosition(first_point_ + direction_ * -distance), targets, total_cross_sections, total_decay_length), distance); } @@ -624,6 +654,7 @@ double Path::GetInteractionDepthFromEndInReverse(double distance, double const & total_decay_length) { EnsureIntersections(); EnsurePoints(); + RequireLastFinite(); return std::copysign(detector_model_->GetInteractionDepthInCGS(intersections_, last_point_, GeometryPosition(last_point_ + direction_ * -distance), targets, total_cross_sections, total_decay_length), distance); } @@ -634,6 +665,7 @@ double Path::GetInteractionDepthFromEndInReverse(double distance, double Path::GetDistanceFromStartInBounds(double column_depth) { EnsureIntersections(); EnsurePoints(); + RequireFirstFinite(); double distance = detector_model_->DistanceForColumnDepthFromPoint(intersections_, first_point_, direction_, column_depth); if(distance > distance_) { distance = distance_; @@ -646,6 +678,7 @@ double Path::GetDistanceFromStartInBounds(double column_depth) { double Path::GetDistanceFromEndInBounds(double column_depth) { EnsureIntersections(); EnsurePoints(); + RequireLastFinite(); double distance = detector_model_->DistanceForColumnDepthFromPoint(intersections_, last_point_, -direction_, column_depth); if(distance > distance_) { distance = distance_; @@ -658,6 +691,7 @@ double Path::GetDistanceFromEndInBounds(double column_depth) { double Path::GetDistanceFromStartAlongPath(double column_depth) { EnsureIntersections(); EnsurePoints(); + RequireFirstFinite(); double distance = detector_model_->DistanceForColumnDepthFromPoint(intersections_, first_point_, direction_, column_depth); return distance; } @@ -665,6 +699,7 @@ double Path::GetDistanceFromStartAlongPath(double column_depth) { double Path::GetDistanceFromEndAlongPath(double column_depth) { EnsureIntersections(); EnsurePoints(); + RequireLastFinite(); double distance = detector_model_->DistanceForColumnDepthFromPoint(intersections_, last_point_, direction_, column_depth); return distance; } @@ -672,6 +707,7 @@ double Path::GetDistanceFromEndAlongPath(double column_depth) { double Path::GetDistanceFromStartInReverse(double column_depth) { EnsureIntersections(); EnsurePoints(); + RequireFirstFinite(); double distance = detector_model_->DistanceForColumnDepthFromPoint(intersections_, first_point_, -direction_, column_depth); return distance; } @@ -679,6 +715,7 @@ double Path::GetDistanceFromStartInReverse(double column_depth) { double Path::GetDistanceFromEndInReverse(double column_depth) { EnsureIntersections(); EnsurePoints(); + RequireLastFinite(); double distance = detector_model_->DistanceForColumnDepthFromPoint(intersections_, last_point_, -direction_, column_depth); return distance; } @@ -693,6 +730,7 @@ double Path::GetDistanceFromStartInBounds(double interaction_depth, double const & total_decay_length) { EnsureIntersections(); EnsurePoints(); + RequireFirstFinite(); double distance = detector_model_->DistanceForInteractionDepthFromPoint(intersections_, first_point_, direction_, interaction_depth, targets, total_cross_sections, total_decay_length); if(distance > distance_) { distance = distance_; @@ -708,6 +746,7 @@ double Path::GetDistanceFromEndInBounds(double interaction_depth, double const & total_decay_length) { EnsureIntersections(); EnsurePoints(); + RequireLastFinite(); double distance = detector_model_->DistanceForInteractionDepthFromPoint(intersections_, last_point_, -direction_, interaction_depth, targets, total_cross_sections, total_decay_length); if(distance > distance_) { distance = distance_; @@ -723,6 +762,7 @@ double Path::GetDistanceFromStartAlongPath(double interaction_depth, double const & total_decay_length) { EnsureIntersections(); EnsurePoints(); + RequireFirstFinite(); double distance = detector_model_->DistanceForInteractionDepthFromPoint(intersections_, first_point_, direction_, interaction_depth, targets, total_cross_sections, total_decay_length); return distance; } @@ -733,6 +773,7 @@ double Path::GetDistanceFromEndAlongPath(double interaction_depth, double const & total_decay_length) { EnsureIntersections(); EnsurePoints(); + RequireLastFinite(); double distance = detector_model_->DistanceForInteractionDepthFromPoint(intersections_, last_point_, direction_, interaction_depth, targets, total_cross_sections, total_decay_length); return distance; } @@ -743,6 +784,7 @@ double Path::GetDistanceFromStartInReverse(double interaction_depth, double const & total_decay_length) { EnsureIntersections(); EnsurePoints(); + RequireFirstFinite(); double distance = detector_model_->DistanceForInteractionDepthFromPoint(intersections_, first_point_, -direction_, interaction_depth, targets, total_cross_sections, total_decay_length); return distance; } @@ -753,6 +795,7 @@ double Path::GetDistanceFromEndInReverse(double interaction_depth, double const & total_decay_length) { EnsureIntersections(); EnsurePoints(); + RequireLastFinite(); double distance = detector_model_->DistanceForInteractionDepthFromPoint(intersections_, last_point_, -direction_, interaction_depth, targets, total_cross_sections, total_decay_length); return distance; } @@ -760,6 +803,7 @@ double Path::GetDistanceFromEndInReverse(double interaction_depth, bool Path::IsWithinBounds(GeometryPosition point) { UpdatePoints(); + RequireBothFinite(); if(set_points_) { double d0 = LI::math::scalar_product(direction_, first_point_ - point); double d1 = LI::math::scalar_product(direction_, last_point_ - point); @@ -771,6 +815,7 @@ bool Path::IsWithinBounds(GeometryPosition point) { double Path::GetDistanceFromStartInBounds(GeometryPosition point) { UpdatePoints(); + RequireFirstFinite(); if(set_points_) { double d0 = LI::math::scalar_product(direction_, point - first_point_); return std::max(0.0, d0); @@ -781,6 +826,7 @@ double Path::GetDistanceFromStartInBounds(GeometryPosition point) { bool Path::IsWithinBounds(DetectorPosition point) { UpdatePoints(); + RequireBothFinite(); if(set_det_points_) { double d0 = LI::math::scalar_product(direction_det_, first_point_det_ - point); double d1 = LI::math::scalar_product(direction_det_, last_point_det_ - point); @@ -794,6 +840,7 @@ bool Path::IsWithinBounds(DetectorPosition point) { double Path::GetDistanceFromStartInBounds(DetectorPosition point) { UpdatePoints(); + RequireFirstFinite(); if(set_det_points_) { double d0 = LI::math::scalar_product(direction_det_, point - first_point_det_); return std::max(0.0, d0); From 2b409a6442d99d6ae70bc374e878610553ee5ece Mon Sep 17 00:00:00 2001 From: Nicholas Kamp Date: Mon, 22 Jan 2024 17:39:59 -0500 Subject: [PATCH 42/42] reverse infinite check order in clip_0/1 --- projects/detector/private/Path.cxx | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/projects/detector/private/Path.cxx b/projects/detector/private/Path.cxx index 1098038bc..c993e8b25 100644 --- a/projects/detector/private/Path.cxx +++ b/projects/detector/private/Path.cxx @@ -261,8 +261,8 @@ void Path::ClipToOuterBounds() { if(dot < 0) { p0.swap(p1); } - bool clip_0 = ((p0 - first_point_) * direction_ > 0) || (first_inf_); - bool clip_1 = ((p1 - last_point_) * direction_ < 0) || (last_inf_); + bool clip_0 = (first_inf_) || ((p0 - first_point_) * direction_ > 0); + bool clip_1 = (last_inf_) || ((p1 - last_point_) * direction_ < 0); bool clip = clip_0 or clip_1; if(clip_0) { first_point_ = GeometryPosition(p0);