From 1ef646de461aeb738e4354585396073dc2e838c4 Mon Sep 17 00:00:00 2001 From: Jean-Marco Alameddine Date: Sat, 15 Apr 2023 15:54:44 +0200 Subject: [PATCH 1/3] Add implementation of MoliereInterpol, based on Moliere, but using interpolation tables to evaluate f1M, F1M, f2M, F2M Add print for multiple_scattering::ScatteringOffset Add warning if iterations of Moliere::CalculateRandomAngle are over 100, and stop iteration process Add direct pybindings for Moliere and add simple constructor for Highland --- docs/config_docu.md | 1 + .../scattering/multiple_scattering/Moliere.h | 59 +++++++++++++---- .../multiple_scattering/Parametrization.h | 4 ++ .../multiple_scattering/ScatteringFactory.h | 8 ++- .../multiple_scattering/Moliere.cxx | 64 ++++++++++++++++++- src/pyPROPOSAL/detail/scattering.cxx | 15 ++++- 6 files changed, 134 insertions(+), 17 deletions(-) diff --git a/docs/config_docu.md b/docs/config_docu.md index ce52c4fe6..98fd76534 100644 --- a/docs/config_docu.md +++ b/docs/config_docu.md @@ -167,6 +167,7 @@ PROPOSAL provides different multiple scattering models, which can be set with th * `HighlandIntegral`: Gaussian approximation of Molière theory, derived by [Highland](https://doi.org/10.1016/0029-554X(75)90743-0) and corrected by [Lynch/Dahl](https://doi.org/10.1016/0168-583X(91)95671-Y). * `Highland`: Same as `HighlandIntegral`, but with the approximation that the particle energy is constant during a propagation step. *This parametrization should only be used for small step sizes where this approximation is valid.* * `Moliere`: Scattering based on [Molière's theory](https://zfn.mpdl.mpg.de/data/Reihe_A/3/ZNA-1948-3a-0078.pdf). *Significantly slower compared to other scattering models.* +* `MoliereInterpol`: Same as `Moliere`, but using interpolation tables to increase performance Multiple scattering effects can be scaled with a constant factor using the option `multiple_scattering_multiplier`. This scales the sampled deflections, which are described by their displacement in cartesian coordinates, with a constant factor for each component. diff --git a/src/PROPOSAL/PROPOSAL/scattering/multiple_scattering/Moliere.h b/src/PROPOSAL/PROPOSAL/scattering/multiple_scattering/Moliere.h index 5c2f12e6d..1f217efed 100644 --- a/src/PROPOSAL/PROPOSAL/scattering/multiple_scattering/Moliere.h +++ b/src/PROPOSAL/PROPOSAL/scattering/multiple_scattering/Moliere.h @@ -33,6 +33,8 @@ #include "PROPOSAL/medium/Medium.h" #include "PROPOSAL/scattering/multiple_scattering/Parametrization.h" +#include "CubicInterpolation/CubicSplines.h" +#include "CubicInterpolation/Interpolant.h" namespace PROPOSAL { struct ParticleDef; @@ -43,11 +45,6 @@ namespace multiple_scattering { bool compare(const Parametrization&) const override; void print(std::ostream&) const override; - std::unique_ptr clone() const final - { - return std::make_unique(*this); - } - int numComp_; // number of components in medium double ZSq_A_average_; std::vector Zi_; // nuclear charge of different components @@ -62,25 +59,55 @@ namespace multiple_scattering { double chiCSq_; // characteristic angle² in rad² std::vector B_; - double f1M(double x); - double f2M(double x); - double f(double theta); - - double F1M(double x); - double F2M(double x); - double F(double theta); - double GetRandom(double pre_factor, double rnd); + protected: + virtual double f1M(double x); + virtual double f2M(double x); + + virtual double F1M(double x); + virtual double F2M(double x); + public: // constructor Moliere(const ParticleDef&, Medium const&); ScatteringOffset CalculateRandomAngle(double grammage, double ei, double ef, const std::array& rnd) override; + + std::unique_ptr clone() const override + { + return std::make_unique(*this); + } }; + + class MoliereInterpol : public Moliere { + using interpolant_t = cubic_splines::Interpolant>; + using interpolant_ptr = std::shared_ptr; + + // interpolation tables + interpolant_ptr f1M_interpolant_; + interpolant_ptr f2M_interpolant_; + interpolant_ptr F1M_interpolant_; + interpolant_ptr F2M_interpolant_; + + double F1M(double x) override; + double F2M(double x) override; + + double f1M(double x) override; + double f2M(double x) override; + + public: + MoliereInterpol(const ParticleDef&, Medium const&); + + std::unique_ptr clone() const override + { + return std::make_unique(*this); + } + }; + } // namespace multiple_scattering template inline auto make_moliere(Args... args) @@ -89,4 +116,10 @@ template inline auto make_moliere(Args... args) new multiple_scattering::Moliere(std::forward(args)...)); } +template inline auto make_moliereinterpol(Args... args) +{ + return std::unique_ptr( + new multiple_scattering::MoliereInterpol(std::forward(args)...)); +} + } // namespace PROPOSAL diff --git a/src/PROPOSAL/PROPOSAL/scattering/multiple_scattering/Parametrization.h b/src/PROPOSAL/PROPOSAL/scattering/multiple_scattering/Parametrization.h index 73273c555..7c02c02b6 100644 --- a/src/PROPOSAL/PROPOSAL/scattering/multiple_scattering/Parametrization.h +++ b/src/PROPOSAL/PROPOSAL/scattering/multiple_scattering/Parametrization.h @@ -35,6 +35,10 @@ namespace PROPOSAL { namespace multiple_scattering { struct ScatteringOffset { ScatteringOffset() : sx(0.), sy(0.), tx(0.), ty(0.) {}; + friend std::ostream& operator<<(std::ostream& os, const ScatteringOffset& offset) { + os << "sx: " << offset.sx << ", sy: " << offset.sy << ", tx: " << offset.tx << ", ty: " << offset.ty << '\n'; + return os; + }; double sx; double sy; double tx; diff --git a/src/PROPOSAL/PROPOSAL/scattering/multiple_scattering/ScatteringFactory.h b/src/PROPOSAL/PROPOSAL/scattering/multiple_scattering/ScatteringFactory.h index 75c11ef90..5fbc8e438 100644 --- a/src/PROPOSAL/PROPOSAL/scattering/multiple_scattering/ScatteringFactory.h +++ b/src/PROPOSAL/PROPOSAL/scattering/multiple_scattering/ScatteringFactory.h @@ -51,14 +51,16 @@ enum class MultipleScatteringType : int { NoScattering, Moliere, Highland, - HighlandIntegral + HighlandIntegral, + MoliereInterpol }; static const std::unordered_map MultipleScatteringTable = { { "moliere", MultipleScatteringType::Moliere }, { "highland", MultipleScatteringType::Highland }, { "highlandintegral", MultipleScatteringType::HighlandIntegral }, - { "noscattering", MultipleScatteringType::NoScattering } }; + { "noscattering", MultipleScatteringType::NoScattering }, + { "moliereinterpol", MultipleScatteringType::MoliereInterpol }}; inline auto make_multiple_scattering( MultipleScatteringType t, ParticleDef const& p, Medium const& m) @@ -68,6 +70,8 @@ inline auto make_multiple_scattering( return make_highland(p, m); case MultipleScatteringType::Moliere: return make_moliere(p, m); + case MultipleScatteringType::MoliereInterpol: + return make_moliereinterpol(p, m); case MultipleScatteringType::NoScattering: return std::unique_ptr(nullptr); default: diff --git a/src/PROPOSAL/detail/PROPOSAL/scattering/multiple_scattering/Moliere.cxx b/src/PROPOSAL/detail/PROPOSAL/scattering/multiple_scattering/Moliere.cxx index dd652777b..7432640af 100644 --- a/src/PROPOSAL/detail/PROPOSAL/scattering/multiple_scattering/Moliere.cxx +++ b/src/PROPOSAL/detail/PROPOSAL/scattering/multiple_scattering/Moliere.cxx @@ -7,6 +7,8 @@ #include "PROPOSAL/particle/ParticleDef.h" #include "PROPOSAL/scattering/multiple_scattering/Coefficients.h" #include "PROPOSAL/scattering/multiple_scattering/Moliere.h" +#include +#include "PROPOSAL/Logging.h" using namespace PROPOSAL::multiple_scattering; @@ -362,11 +364,71 @@ double Moliere::GetRandom(double pre_factor, double rnd) rnd = rnd - 0.5; // iterating until the number of correct digits is greater than 4 + unsigned int i = 0; do { + i++; theta_n = theta_np1; theta_np1 = theta_n - (F(theta_n) - rnd) / f(theta_n); - + if (i == 100) { + Logging::Get("proposal.scattering")->warn( + "Iteration in Moliere::GetRandom did not converge after 100 iterations. " + "Return current value theta = {} (previous iteration: theta = {}", theta_np1, theta_n); + return theta_np1; + } } while (std::abs((theta_n - theta_np1) / theta_np1) > 1e-4); return theta_np1; } + +MoliereInterpol::MoliereInterpol(const ParticleDef& p_def, const Medium& medium) : Moliere(p_def, medium) { + // initialize interpolation tables + Logging::Get("proposal.scattering")->debug("Initialize interpooation tables for MoliereInterpol."); + + auto def_f1M = cubic_splines::CubicSplines::Definition(); + def_f1M.f = [&](double x) {return Moliere::f1M(x);}; + def_f1M.axis = std::make_unique>(0.f, 20, size_t(100)); + f1M_interpolant_ = std::make_shared(std::move(def_f1M)); + + auto def_f2M = cubic_splines::CubicSplines::Definition(); + def_f2M.f = [&](double x) {return Moliere::f2M(x);}; + def_f2M.axis = std::make_unique>(0.f, 20, size_t(100)); + f2M_interpolant_ = std::make_shared(std::move(def_f2M)); + + auto def_F1M = cubic_splines::CubicSplines::Definition(); + def_F1M.f = [&](double x) {return Moliere::F1M(x);}; + def_F1M.axis = std::make_unique>(0.f, 20, size_t(100)); + F1M_interpolant_ = std::make_shared(std::move(def_F1M)); + + auto def_F2M = cubic_splines::CubicSplines::Definition(); + def_F2M.f = [&](double x) {return Moliere::F2M(x);}; + def_F2M.axis = std::make_unique>(0.f, 20, size_t(100)); + F2M_interpolant_ = std::make_shared(std::move(def_F2M)); +} + +double MoliereInterpol::f1M(double x) +{ + if (x > 20.) + return Moliere::f1M(x); // use analytical evaluation outside range of interpolation tables + return f1M_interpolant_->evaluate(x); +} + +double MoliereInterpol::f2M(double x) +{ + if (x > 20) + return Moliere::f2M(x); // use analytical evaluation outside range of interpolation tables + return f2M_interpolant_->evaluate(x); +} + +double MoliereInterpol::F1M(double x) +{ + if (x > 20.) + return Moliere::F1M(x); // use analytical evaluation outside range of interpolation tables + return F1M_interpolant_->evaluate(x); +} + +double MoliereInterpol::F2M(double x) +{ + if (x > 20) + return Moliere::F2M(x); // use analytical evaluation outside range of interpolation tables + return F2M_interpolant_->evaluate(x); +} \ No newline at end of file diff --git a/src/pyPROPOSAL/detail/scattering.cxx b/src/pyPROPOSAL/detail/scattering.cxx index 6e48e3d30..6b6cfe355 100644 --- a/src/pyPROPOSAL/detail/scattering.cxx +++ b/src/pyPROPOSAL/detail/scattering.cxx @@ -34,6 +34,8 @@ void init_scattering(py::module& m) py::class_>(m_sub, "Highland") + .def(py::init(), + py::arg("particle_def"), py::arg("medium")) .def("CalculateTheta0", &multiple_scattering::Highland::CalculateTheta0, py::arg("grammage"), py::arg("e_i"), py::arg("e_f"), R"pbdoc( @@ -49,11 +51,22 @@ void init_scattering(py::module& m) py::class_>(m_sub, "HighlandIntegral"); + py::class_>(m_sub, "Moliere") + .def(py::init(), + py::arg("particle_def"), py::arg("medium")); + + py::class_>(m_sub, "MoliereInterpol") + .def(py::init(), + py::arg("particle_def"), py::arg("medium")); + py::class_(m_sub, "scattering_offset") .def_readwrite("sx", &multiple_scattering::ScatteringOffset::sx) .def_readwrite("sy", &multiple_scattering::ScatteringOffset::sy) .def_readwrite("tx", &multiple_scattering::ScatteringOffset::tx) - .def_readwrite("ty", &multiple_scattering::ScatteringOffset::ty); + .def_readwrite("ty", &multiple_scattering::ScatteringOffset::ty) + .def("__str__", &py_print); m_sub.def("scatter_initial_direction", &multiple_scattering::ScatterInitialDirection); From 4bd2369b782562bd880aedd3ebceba8e8373985e Mon Sep 17 00:00:00 2001 From: Jean-Marco Alameddine Date: Sat, 15 Apr 2023 21:03:53 +0200 Subject: [PATCH 2/3] Add MoliereInterpol to scattering tests --- tests/Scattering_TEST.cxx | 17 ++++++++++------- 1 file changed, 10 insertions(+), 7 deletions(-) diff --git a/tests/Scattering_TEST.cxx b/tests/Scattering_TEST.cxx index 84d60606e..b80e72c0d 100644 --- a/tests/Scattering_TEST.cxx +++ b/tests/Scattering_TEST.cxx @@ -215,9 +215,10 @@ TEST(Scattering, ZeroDisplacement){ auto cuts = std::make_shared(INF, 1, false); auto cross = GetCrossSections(MuMinusDef(), medium, cuts, true); - std::array, 3> scatter_list = {make_multiple_scattering("moliere", MuMinusDef(), medium), + std::array, 4> scatter_list = {make_multiple_scattering("moliere", MuMinusDef(), medium), make_multiple_scattering("highland", MuMinusDef(), medium), - make_multiple_scattering("highlandintegral", MuMinusDef(), medium, cross)}; + make_multiple_scattering("highlandintegral", MuMinusDef(), medium, cross), + make_multiple_scattering("moliereinterpol", MuMinusDef(), medium)}; for(auto const& scatter: scatter_list){ auto offset = scatter->CalculateRandomAngle(0, 1e5, 1e5, {0.1, 0.2, 0.3, 0.4}); @@ -238,9 +239,10 @@ TEST(Scattering, FirstMomentum){ int statistics = 1e3; auto cross = GetCrossSections(MuMinusDef(), medium, cuts, true); - std::array, 3> scatter_list = {make_multiple_scattering("moliere", MuMinusDef(), medium), + std::array, 4> scatter_list = {make_multiple_scattering("moliere", MuMinusDef(), medium), make_multiple_scattering("highland", MuMinusDef(), medium), - make_multiple_scattering("highlandintegral", MuMinusDef(), medium, cross)}; + make_multiple_scattering("highlandintegral", MuMinusDef(), medium, cross), + make_multiple_scattering("moliereinterpol", MuMinusDef(), medium),}; Cartesian3D scatter_sum; Cartesian3D offset_sum; @@ -277,15 +279,16 @@ TEST(Scattering, SecondMomentum){ auto direction_init = Cartesian3D(0, 0, 1); auto cuts = std::make_shared(INF, 1, false); - int statistics = 1e3; + int statistics = 1e5; double E_i = 1e14; std::array final_energies = {1e13, 1e11, 1e9, 1e7, 1e5, 1e3}; auto cross = GetCrossSections(MuMinusDef(), medium, cuts, true); - std::array, 3> scatter_list = {make_multiple_scattering("moliere", MuMinusDef(), medium), + std::array, 4> scatter_list = {make_multiple_scattering("moliere", MuMinusDef(), medium), make_multiple_scattering("highland", MuMinusDef(), medium), - make_multiple_scattering("highlandintegral", MuMinusDef(), medium, cross)}; + make_multiple_scattering("highlandintegral", MuMinusDef(), medium, cross), + make_multiple_scattering("moliereinterpol", MuMinusDef(), medium, cross)}; double scatter_sum; double offset_sum; double displacement; From f04afcf2018dac7e69147199ce025c262cdbad6e Mon Sep 17 00:00:00 2001 From: Jean-Marco Alameddine Date: Mon, 17 Apr 2023 13:33:52 +0200 Subject: [PATCH 3/3] Add test to compare results of Moliere to MoliereInterpol --- tests/Scattering_TEST.cxx | 50 +++++++++++++++++++++++++++++++++++++++ 1 file changed, 50 insertions(+) diff --git a/tests/Scattering_TEST.cxx b/tests/Scattering_TEST.cxx index b80e72c0d..cafc43c01 100644 --- a/tests/Scattering_TEST.cxx +++ b/tests/Scattering_TEST.cxx @@ -500,6 +500,56 @@ TEST(Scattering, NoScattering) EXPECT_EQ(std::get<1>(new_dir), init_dir); } +TEST(MoliereInterpol, ComparionToMoliere) { + // MoliereInterpol is based on Moliere, using interpolation tables to speed up the evaluation. + // Therefore, we expect results to be similar + + RandomGenerator::Get().SetSeed(24601); + auto medium = StandardRock(); + auto position_init = Cartesian3D(0, 0, 0); + auto direction_init = Cartesian3D(0, 0, 1); + + auto cuts = std::make_shared(INF, 1, false); + + std::vector particles = {EMinusDef(), MuMinusDef()}; + + for (auto p : particles) { + //displacement calculator + auto cross = GetCrossSections(p, medium, cuts, false); + auto displacement = DisplacementBuilder(cross, std::false_type()); + + auto moliere = make_multiple_scattering("moliere", p, medium); + auto moliere_interpol = make_multiple_scattering("moliereinterpol", p, medium); + double E_f = 1e5; + auto energies = std::array{2e5, 1e6, 1e7, 1e8, 1e9, 1e10}; + for (auto E_i : energies) { + double grammage = displacement.SolveTrackIntegral(E_i, E_f); + auto rnd = std::array{RandomGenerator::Get().RandomDouble(), + RandomGenerator::Get().RandomDouble(), + RandomGenerator::Get().RandomDouble(), + RandomGenerator::Get().RandomDouble()}; + auto coords = moliere->CalculateRandomAngle(grammage, E_i, E_f, rnd); + auto coords_interpol = moliere_interpol->CalculateRandomAngle(grammage, E_i, E_f, rnd); + auto vec = multiple_scattering::ScatterInitialDirection( + direction_init, coords); + auto vec_interpol = multiple_scattering::ScatterInitialDirection( + direction_init, coords_interpol); + EXPECT_NEAR(std::get<0>(vec).GetX(), std::get<0>(vec_interpol).GetX(), + std::abs(std::get<0>(vec).GetX() * 1e-2)); + EXPECT_NEAR(std::get<0>(vec).GetY(), std::get<0>(vec_interpol).GetY(), + std::abs(std::get<0>(vec).GetY() * 1e-2)); + EXPECT_NEAR(std::get<0>(vec).GetZ(), std::get<0>(vec_interpol).GetZ(), + std::abs(std::get<0>(vec).GetZ() * 1e-2)); + EXPECT_NEAR(std::get<1>(vec).GetX(), std::get<1>(vec_interpol).GetX(), + std::abs(std::get<1>(vec).GetX() * 1e-2)); + EXPECT_NEAR(std::get<1>(vec).GetY(), std::get<1>(vec_interpol).GetY(), + std::abs(std::get<1>(vec).GetY() * 1e-2)); + EXPECT_NEAR(std::get<1>(vec).GetZ(), std::get<1>(vec_interpol).GetZ(), + std::abs(std::get<1>(vec).GetZ() * 1e-2)); + } + } +} + int main(int argc, char** argv) { ::testing::InitGoogleTest(&argc, argv);