From c82b8207e17183d485456a380978f53e9801ed0a Mon Sep 17 00:00:00 2001 From: Jean-Marco Alameddine Date: Wed, 6 Dec 2023 15:28:00 +0100 Subject: [PATCH] Adapt transformation function for PhotoPairProduction interpolation tables to increase numerical efficiency. Add constructor to use Integral instead if Interpolation Crosssections for secondaries::PhotoPairProductionInterpolant classes --- .../CrossSectionDNDXInterpolant.h | 19 +++++++++++-- .../PhotoPairProductionInterpolant.h | 4 +-- .../PhotoPairProductionKochMotz.h | 8 +++--- .../PhotoPairProductionTsai.h | 12 ++++---- .../CrossSectionDNDXInterpolant.cxx | 18 ++++++++++++ .../parametrization/PhotoPairProduction.cxx | 2 +- src/pyPROPOSAL/detail/secondaries.cxx | 28 +++++++++++++++---- 7 files changed, 71 insertions(+), 20 deletions(-) diff --git a/src/PROPOSAL/PROPOSAL/crosssection/CrossSectionDNDX/CrossSectionDNDXInterpolant.h b/src/PROPOSAL/PROPOSAL/crosssection/CrossSectionDNDX/CrossSectionDNDXInterpolant.h index 80bf7eeb7..97ab8b252 100644 --- 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 ::value || std::is_base_of::value), bool> = true> +template ::value || std::is_base_of::value || std::is_base_of::value), bool> = true> double transform_loss(double v_cut, double v_max, double v) { return transform_relative_loss(v_cut, v_max, v); } -template ::value || std::is_base_of::value), bool> = true> +template ::value || std::is_base_of::value || std::is_base_of::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 ::value, bool> = true> + double transform_loss(double v_cut, double v_max, double v) { + return transform_loss_linear(v_cut, v_max, v); + } + + template ::value, bool> = true> + double retransform_loss(double v_cut, double v_max, double v) { + return retransform_loss_linear(v_cut, v_max, v); + } + template auto build_dndx_def(T1 const& param, ParticleDef const& p, Args... args) { diff --git a/src/PROPOSAL/PROPOSAL/secondaries/parametrization/photopairproduction/PhotoPairProductionInterpolant.h b/src/PROPOSAL/PROPOSAL/secondaries/parametrization/photopairproduction/PhotoPairProductionInterpolant.h index 783a63df1..4759ae9e1 100644 --- a/src/PROPOSAL/PROPOSAL/secondaries/parametrization/photopairproduction/PhotoPairProductionInterpolant.h +++ b/src/PROPOSAL/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/src/PROPOSAL/PROPOSAL/secondaries/parametrization/photopairproduction/PhotoPairProductionKochMotz.h b/src/PROPOSAL/PROPOSAL/secondaries/parametrization/photopairproduction/PhotoPairProductionKochMotz.h index 0fa2752a0..dce5b63f3 100644 --- a/src/PROPOSAL/PROPOSAL/secondaries/parametrization/photopairproduction/PhotoPairProductionKochMotz.h +++ b/src/PROPOSAL/PROPOSAL/secondaries/parametrization/photopairproduction/PhotoPairProductionKochMotz.h @@ -10,8 +10,8 @@ namespace PROPOSAL { public: PhotoPairProductionKochMotzForwardPeaked() = default; - PhotoPairProductionKochMotzForwardPeaked(ParticleDef p, Medium m) - : PhotoPairProductionInterpolant(p, m) + PhotoPairProductionKochMotzForwardPeaked(ParticleDef p, Medium m, bool interpol = true) + : PhotoPairProductionInterpolant(p, m, interpol) {} }; @@ -20,8 +20,8 @@ namespace PROPOSAL { public DefaultSecondaries { 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 CalculateDirections( diff --git a/src/PROPOSAL/PROPOSAL/secondaries/parametrization/photopairproduction/PhotoPairProductionTsai.h b/src/PROPOSAL/PROPOSAL/secondaries/parametrization/photopairproduction/PhotoPairProductionTsai.h index 715e2fc69..bef94efed 100644 --- a/src/PROPOSAL/PROPOSAL/secondaries/parametrization/photopairproduction/PhotoPairProductionTsai.h +++ b/src/PROPOSAL/PROPOSAL/secondaries/parametrization/photopairproduction/PhotoPairProductionTsai.h @@ -11,8 +11,8 @@ namespace secondaries { public: PhotoPairProductionTsaiForwardPeaked() = default; - PhotoPairProductionTsaiForwardPeaked(ParticleDef p, Medium m) - : PhotoPairProductionInterpolant(p, m) + PhotoPairProductionTsaiForwardPeaked(ParticleDef p, Medium m, bool interpol = true) + : PhotoPairProductionInterpolant(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 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 CalculateDirections( diff --git a/src/PROPOSAL/detail/PROPOSAL/crosssection/CrossSectionDNDX/CrossSectionDNDXInterpolant.cxx b/src/PROPOSAL/detail/PROPOSAL/crosssection/CrossSectionDNDX/CrossSectionDNDXInterpolant.cxx index 868280357..a95ae4080 100644 --- 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 index 20a407334..d31334a73 100644 --- 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 index 497cd3adb..5ec828608 100644 --- a/src/pyPROPOSAL/detail/secondaries.cxx +++ b/src/pyPROPOSAL/detail/secondaries.cxx @@ -54,6 +54,24 @@ template struct SecondariesBuilder { :math:`\frac{E_-}{E_{\gamma}}`. )pbdoc"); } + + template auto decl_rho_param_interpol(Args... args) + { + py::class_>(args...) + .def(py::init(), 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 {} - .decl_rho_param(m_sub, "PhotoTsai"); + .decl_rho_param_interpol(m_sub, "PhotoTsai"); SecondariesBuilder {} - .decl_rho_param(m_sub, "PhotoTsaiForwardPeaked"); + .decl_rho_param_interpol(m_sub, "PhotoTsaiForwardPeaked"); SecondariesBuilder {} - .decl_rho_param(m_sub, "PhotoTsaiSauter"); + .decl_rho_param_interpol(m_sub, "PhotoTsaiSauter"); SecondariesBuilder {} - .decl_rho_param(m_sub, "PhotoKochMotzForwardPeaked"); + .decl_rho_param_interpol(m_sub, "PhotoKochMotzForwardPeaked"); SecondariesBuilder {} - .decl_rho_param(m_sub, "PhotoKochMotzSauter"); + .decl_rho_param_interpol(m_sub, "PhotoKochMotzSauter"); SecondariesBuilder {} .decl_param(m_sub, "PhotoeffectNoDeflection");