Adapt transformation function for PhotoPairProduction interpolation t…

…ables to increase numerical efficiency. Add constructor to use Integral instead if Interpolation Crosssections for secondaries::PhotoPairProductionInterpolant classes
tudo-astroparticlephysics · Dec 6, 2023 · c82b820 · c82b820
1 parent c3b01b7
commit c82b820
Show file tree

Hide file tree

Showing 7 changed files with 71 additions and 20 deletions.
diff --git a/src/PROPOSAL/PROPOSAL/crosssection/CrossSectionDNDX/CrossSectionDNDXInterpolant.h b/src/PROPOSAL/PROPOSAL/crosssection/CrossSectionDNDX/CrossSectionDNDXInterpolant.h
@@ -16,21 +16,24 @@ double transform_relative_loss(double v_cut, double v_max, double v, double c =
 double retransform_relative_loss(double v_cut, double v_max, double v, double c = 1);
 double transform_loss_log(double v_cut, double v_max, double v);
 double retransform_loss_log(double v_cut, double v_max, double v);
+double transform_loss_linear(double v_cut, double v_max, double v);
+double retransform_loss_linear(double v_cut, double v_max, double v);
 
 namespace crosssection {
     struct Compton;
     struct Ionization;
+    struct PhotoPairProduction;
 }
 
 // general transformation rules for v nodes
 // TODO: in C++17, this section can be rewritten much cleaner using if constexpr
 
-template <typename Param, std::enable_if_t<!(std::is_base_of<crosssection::Ionization, Param>::value || std::is_base_of<crosssection::Compton, Param>::value), bool> = true>
+template <typename Param, std::enable_if_t<!(std::is_base_of<crosssection::Ionization, Param>::value || std::is_base_of<crosssection::Compton, Param>::value || std::is_base_of<crosssection::PhotoPairProduction, Param>::value), bool> = true>
 double transform_loss(double v_cut, double v_max, double v) {
     return transform_relative_loss(v_cut, v_max, v);
 }
 
-template <typename Param, std::enable_if_t<!(std::is_base_of<crosssection::Ionization, Param>::value || std::is_base_of<crosssection::Compton, Param>::value), bool> = true>
+template <typename Param, std::enable_if_t<!(std::is_base_of<crosssection::Ionization, Param>::value || std::is_base_of<crosssection::Compton, Param>::value || std::is_base_of<crosssection::PhotoPairProduction, Param>::value), bool> = true>
 double retransform_loss(double v_cut, double v_max, double v) {
     return retransform_relative_loss(v_cut, v_max, v);
 }
@@ -59,6 +62,18 @@ double retransform_loss(double v_cut, double v_max, double v) {
     return retransform_loss_log(v_cut, v_max, v);
 }
 
+// specifications for PhotoPairProduction
+
+    template <typename Param, std::enable_if_t<std::is_base_of<crosssection::PhotoPairProduction, Param>::value, bool> = true>
+    double transform_loss(double v_cut, double v_max, double v) {
+        return transform_loss_linear(v_cut, v_max, v);
+    }
+
+    template <typename Param, std::enable_if_t<std::is_base_of<crosssection::PhotoPairProduction, Param>::value, bool> = true>
+    double retransform_loss(double v_cut, double v_max, double v) {
+        return retransform_loss_linear(v_cut, v_max, v);
+    }
+
 template <typename T1, typename... Args>
 auto build_dndx_def(T1 const& param, ParticleDef const& p, Args... args)
 {

diff --git a/...PROPOSAL/secondaries/parametrization/photopairproduction/PhotoPairProductionInterpolant.h b/...PROPOSAL/secondaries/parametrization/photopairproduction/PhotoPairProductionInterpolant.h
@@ -28,9 +28,9 @@ namespace secondaries {
         static constexpr int n_rnd = 4;
 
         PhotoPairProductionInterpolant() = default;
-        PhotoPairProductionInterpolant(ParticleDef p, Medium m)
+        PhotoPairProductionInterpolant(ParticleDef p, Medium m, bool interpol = true)
             : medium(m)
-            , dndx(detail::build_dndx(std::true_type {}, true,
+            , dndx(detail::build_dndx(std::true_type {}, interpol,
                                       Param(), p, medium, nullptr))
         {}
 

diff --git a/...AL/PROPOSAL/secondaries/parametrization/photopairproduction/PhotoPairProductionKochMotz.h b/...AL/PROPOSAL/secondaries/parametrization/photopairproduction/PhotoPairProductionKochMotz.h
@@ -10,8 +10,8 @@ namespace PROPOSAL {
 
         public:
             PhotoPairProductionKochMotzForwardPeaked() = default;
-            PhotoPairProductionKochMotzForwardPeaked(ParticleDef p, Medium m)
-                : PhotoPairProductionInterpolant<crosssection::PhotoPairKochMotz>(p, m)
+            PhotoPairProductionKochMotzForwardPeaked(ParticleDef p, Medium m, bool interpol = true)
+                : PhotoPairProductionInterpolant<crosssection::PhotoPairKochMotz>(p, m, interpol)
             {}
 
         };
@@ -20,8 +20,8 @@ namespace PROPOSAL {
                                                   public DefaultSecondaries<PhotoPairProductionKochMotzSauter> {
         public:
             PhotoPairProductionKochMotzSauter() = default;
-            PhotoPairProductionKochMotzSauter(ParticleDef p, Medium m)
-                : PhotoPairProductionKochMotzForwardPeaked(p, m)
+            PhotoPairProductionKochMotzSauter(ParticleDef p, Medium m, bool interpol = true)
+                : PhotoPairProductionKochMotzForwardPeaked(p, m, interpol)
             {}
 
             std::tuple<Cartesian3D, Cartesian3D> CalculateDirections(

diff --git a/...OPOSAL/PROPOSAL/secondaries/parametrization/photopairproduction/PhotoPairProductionTsai.h b/...OPOSAL/PROPOSAL/secondaries/parametrization/photopairproduction/PhotoPairProductionTsai.h
@@ -11,8 +11,8 @@ namespace secondaries {
 
     public:
         PhotoPairProductionTsaiForwardPeaked() = default;
-        PhotoPairProductionTsaiForwardPeaked(ParticleDef p, Medium m)
-            : PhotoPairProductionInterpolant<crosssection::PhotoPairTsai>(p, m)
+        PhotoPairProductionTsaiForwardPeaked(ParticleDef p, Medium m, bool interpol = true)
+            : PhotoPairProductionInterpolant<crosssection::PhotoPairTsai>(p, m, interpol)
             {}
 
     };
@@ -24,8 +24,8 @@ namespace secondaries {
 
     public:
         PhotoPairProductionTsai() = default;
-        PhotoPairProductionTsai(ParticleDef p, Medium m)
-            : PhotoPairProductionTsaiForwardPeaked(p, m) {}
+        PhotoPairProductionTsai(ParticleDef p, Medium m, bool interpol = true)
+            : PhotoPairProductionTsaiForwardPeaked(p, m, interpol) {}
 
         std::tuple<Cartesian3D, Cartesian3D> CalculateDirections(
                 const Vector3D&, double, double, const Component&,
@@ -35,8 +35,8 @@ namespace secondaries {
     class PhotoPairProductionTsaiSauter : public PhotoPairProductionTsaiForwardPeaked, public SauterSampling {
     public:
         PhotoPairProductionTsaiSauter() = default;
-        PhotoPairProductionTsaiSauter(ParticleDef p, Medium m)
-            : PhotoPairProductionTsaiForwardPeaked(p, m)
+        PhotoPairProductionTsaiSauter(ParticleDef p, Medium m, bool interpol = true)
+            : PhotoPairProductionTsaiForwardPeaked(p, m, interpol)
         {}
 
         std::tuple<Cartesian3D, Cartesian3D> CalculateDirections(

diff --git a/src/PROPOSAL/detail/PROPOSAL/crosssection/CrossSectionDNDX/CrossSectionDNDXInterpolant.cxx b/src/PROPOSAL/detail/PROPOSAL/crosssection/CrossSectionDNDX/CrossSectionDNDXInterpolant.cxx
@@ -39,6 +39,24 @@ namespace PROPOSAL {
         auto xi = std::log((1. - v_cut)/(1 - v_max));
         return 1. - std::log((1. - v)/(1. - v_max)) / xi;
     }
+
+    double transform_loss_linear(double v_cut, double v_max, double v)
+    {
+        if (v < 0 || v_max == 0)
+            return v_cut;
+        if (v >= 1)
+            return v_max;
+        return v_cut + (v_max - v_cut) * v;
+    }
+
+    double retransform_loss_linear(double v_cut, double v_max, double v)
+    {
+        if (v <= v_cut)
+            return 0;
+        if (v >= v_max)
+            return 1;
+        return (v - v_cut) / (v_max - v_cut);
+    }
 }
 
 std::string CrossSectionDNDXInterpolant::gen_path() const

diff --git a/src/PROPOSAL/detail/PROPOSAL/crosssection/parametrization/PhotoPairProduction.cxx b/src/PROPOSAL/detail/PROPOSAL/crosssection/parametrization/PhotoPairProduction.cxx
@@ -126,7 +126,7 @@ double crosssection::PhotoPairKochMotz::DifferentialCrossSection(
     };
     auto dNdx_nocorrection = i.Integrate(limits.v_min, limits.v_max, integrand, 3);
     auto A = interpolant_->InterpolateArray(comp.GetNucCharge(), energy) / dNdx_nocorrection;
-    return A * DifferentialCrossSectionWithoutA(p, comp, energy, v);
+    return std::max(A * DifferentialCrossSectionWithoutA(p, comp, energy, v), 0.);
 }
 
 double crosssection::PhotoPairKochMotz::DifferentialCrossSectionWithoutA(

diff --git a/src/pyPROPOSAL/detail/secondaries.cxx b/src/pyPROPOSAL/detail/secondaries.cxx
@@ -54,6 +54,24 @@ template <typename T, typename BaseT> struct SecondariesBuilder {
                 :math:`\frac{E_-}{E_{\gamma}}`.
             )pbdoc");
     }
+
+    template <typename... Args> auto decl_rho_param_interpol(Args... args)
+    {
+        py::class_<T, BaseT, std::shared_ptr<T>>(args...)
+                .def(py::init<ParticleDef, Medium, bool>(), py::arg("particle_def"),
+                     py::arg("medium"), py::arg("interpolate") = true)
+                .def("calculate_rho", &T::CalculateRho,
+                     R"pbdoc(
+                Samples the asymmetry factor for interactions where two particles
+                are created. For EpairProduction and MupairProduction, this
+                factor is defined as :math:`\frac{E_+ - E_-}{E_+ + E_-}`, where
+                :math:`E_-` is the energy of the created particle and :math:`E_+`
+                the energy of the created antiparticle. For annihilation and
+                photopairproduction, this factor is defined as
+                :math:`\frac{E_{\gamma,1}}{E_+ + m_e}`, respectively
+                :math:`\frac{E_-}{E_{\gamma}}`.
+            )pbdoc");
+    }
 };
 
 void init_secondaries(py::module& m)
@@ -147,19 +165,19 @@ void init_secondaries(py::module& m)
 
     SecondariesBuilder<secondaries::PhotoPairProductionTsai,
         secondaries::PhotoPairProduction> {}
-        .decl_rho_param(m_sub, "PhotoTsai");
+        .decl_rho_param_interpol(m_sub, "PhotoTsai");
     SecondariesBuilder<secondaries::PhotoPairProductionTsaiForwardPeaked,
         secondaries::PhotoPairProduction> {}
-        .decl_rho_param(m_sub, "PhotoTsaiForwardPeaked");
+        .decl_rho_param_interpol(m_sub, "PhotoTsaiForwardPeaked");
     SecondariesBuilder<secondaries::PhotoPairProductionTsaiSauter,
         secondaries::PhotoPairProduction> {}
-        .decl_rho_param(m_sub, "PhotoTsaiSauter");
+        .decl_rho_param_interpol(m_sub, "PhotoTsaiSauter");
     SecondariesBuilder<secondaries::PhotoPairProductionKochMotzForwardPeaked,
         secondaries::PhotoPairProduction> {}
-        .decl_rho_param(m_sub, "PhotoKochMotzForwardPeaked");
+        .decl_rho_param_interpol(m_sub, "PhotoKochMotzForwardPeaked");
     SecondariesBuilder<secondaries::PhotoPairProductionKochMotzSauter,
         secondaries::PhotoPairProduction> {}
-        .decl_rho_param(m_sub, "PhotoKochMotzSauter");
+        .decl_rho_param_interpol(m_sub, "PhotoKochMotzSauter");
 
     SecondariesBuilder<secondaries::PhotoeffectNoDeflection, secondaries::Photoeffect> {}
         .decl_param(m_sub, "PhotoeffectNoDeflection");