Merge pull request #362 from tudo-astroparticlephysics/moliere_interpol

Create `MoliereInterpol` parametrization for multiple scattering

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. 

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<Parametrization> clone() const final
-        {
-            return std::make_unique<Moliere>(*this);
-        }
-
         int numComp_; // number of components in medium
         double ZSq_A_average_;
         std::vector<double> Zi_; // nuclear charge of different components
@@ -62,25 +59,55 @@ namespace multiple_scattering {
         double chiCSq_; // characteristic angle² in rad²
         std::vector<double> 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<double, 4>& rnd) override;
+
+        std::unique_ptr<Parametrization> clone() const override
+        {
+            return std::make_unique<Moliere>(*this);
+        }
     };
+
+    class MoliereInterpol : public Moliere {
+        using interpolant_t = cubic_splines::Interpolant<cubic_splines::CubicSplines<double>>;
+        using interpolant_ptr = std::shared_ptr<interpolant_t>;
+
+        // 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<Parametrization> clone() const override
+        {
+            return std::make_unique<MoliereInterpol>(*this);
+        }
+    };
+
 } // namespace multiple_scattering
 
 template <typename... Args> inline auto make_moliere(Args... args)
@@ -89,4 +116,10 @@ template <typename... Args> inline auto make_moliere(Args... args)
         new multiple_scattering::Moliere(std::forward<Args>(args)...));
 }
 
+template <typename... Args> inline auto make_moliereinterpol(Args... args)
+{
+    return std::unique_ptr<multiple_scattering::Parametrization>(
+        new multiple_scattering::MoliereInterpol(std::forward<Args>(args)...));
+}
+
 } // namespace PROPOSAL

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;

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<std::string, MultipleScatteringType>
     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<multiple_scattering::Parametrization>(nullptr);
     default:

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 <CubicInterpolation/Axis.h>
+#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<double>::Definition();
+    def_f1M.f = [&](double x) {return Moliere::f1M(x);};
+    def_f1M.axis = std::make_unique<cubic_splines::LinAxis<double>>(0.f, 20, size_t(100));
+    f1M_interpolant_ = std::make_shared<interpolant_t>(std::move(def_f1M));
+
+    auto def_f2M = cubic_splines::CubicSplines<double>::Definition();
+    def_f2M.f = [&](double x) {return Moliere::f2M(x);};
+    def_f2M.axis = std::make_unique<cubic_splines::LinAxis<double>>(0.f, 20, size_t(100));
+    f2M_interpolant_ = std::make_shared<interpolant_t>(std::move(def_f2M));
+
+    auto def_F1M = cubic_splines::CubicSplines<double>::Definition();
+    def_F1M.f = [&](double x) {return Moliere::F1M(x);};
+    def_F1M.axis = std::make_unique<cubic_splines::LinAxis<double>>(0.f, 20, size_t(100));
+    F1M_interpolant_ = std::make_shared<interpolant_t>(std::move(def_F1M));
+
+    auto def_F2M = cubic_splines::CubicSplines<double>::Definition();
+    def_F2M.f = [&](double x) {return Moliere::F2M(x);};
+    def_F2M.axis = std::make_unique<cubic_splines::LinAxis<double>>(0.f, 20, size_t(100));
+    F2M_interpolant_ = std::make_shared<interpolant_t>(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);
+}

src/pyPROPOSAL/detail/scattering.cxx

@@ -34,6 +34,8 @@ void init_scattering(py::module& m)
 
     py::class_<multiple_scattering::Highland, multiple_scattering::Parametrization,
         std::shared_ptr<multiple_scattering::Highland>>(m_sub, "Highland")
+        .def(py::init<const ParticleDef&, const Medium&>(),
+             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_<multiple_scattering::HighlandIntegral, multiple_scattering::Highland,
             std::shared_ptr<multiple_scattering::HighlandIntegral>>(m_sub, "HighlandIntegral");
 
+    py::class_<multiple_scattering::Moliere, multiple_scattering::Parametrization,
+            std::shared_ptr<multiple_scattering::Moliere>>(m_sub, "Moliere")
+            .def(py::init<const ParticleDef&, const Medium&>(),
+                 py::arg("particle_def"), py::arg("medium"));
+
+    py::class_<multiple_scattering::MoliereInterpol, multiple_scattering::Parametrization,
+            std::shared_ptr<multiple_scattering::MoliereInterpol>>(m_sub, "MoliereInterpol")
+            .def(py::init<const ParticleDef&, const Medium&>(),
+                 py::arg("particle_def"), py::arg("medium"));
+
     py::class_<multiple_scattering::ScatteringOffset>(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<multiple_scattering::ScatteringOffset>);
 
     m_sub.def("scatter_initial_direction", &multiple_scattering::ScatterInitialDirection);
 

tests/Scattering_TEST.cxx

@@ -215,9 +215,10 @@ TEST(Scattering, ZeroDisplacement){
     auto cuts = std::make_shared<EnergyCutSettings>(INF, 1, false);
     auto cross = GetCrossSections(MuMinusDef(), medium, cuts, true);
 
-    std::array<std::unique_ptr<multiple_scattering::Parametrization>, 3> scatter_list = {make_multiple_scattering("moliere", MuMinusDef(), medium),
+    std::array<std::unique_ptr<multiple_scattering::Parametrization>, 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<std::unique_ptr<multiple_scattering::Parametrization>, 3> scatter_list = {make_multiple_scattering("moliere", MuMinusDef(), medium),
+    std::array<std::unique_ptr<multiple_scattering::Parametrization>, 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<EnergyCutSettings>(INF, 1, false);
 
-    int statistics = 1e3;
+    int statistics = 1e5;
 
     double E_i = 1e14;
     std::array<double, 6> final_energies = {1e13, 1e11, 1e9, 1e7, 1e5, 1e3};
     auto cross = GetCrossSections(MuMinusDef(), medium, cuts, true);
 
-    std::array<std::unique_ptr<multiple_scattering::Parametrization>, 3> scatter_list = {make_multiple_scattering("moliere", MuMinusDef(), medium),
+    std::array<std::unique_ptr<multiple_scattering::Parametrization>, 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;
@@ -497,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<EnergyCutSettings>(INF, 1, false);
+
+    std::vector<ParticleDef> 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<double, 6>{2e5, 1e6, 1e7, 1e8, 1e9, 1e10};
+        for (auto E_i : energies) {
+            double grammage = displacement.SolveTrackIntegral(E_i, E_f);
+            auto rnd = std::array<double, 4>{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);