Merge pull request #128 from tudo-astroparticlephysics/photonucl

Implementation of new Photonuclear Parametrizations
tudo-astroparticlephysics · Mar 16, 2021 · 37059e4 · 37059e4
2 parents 2fa6770 + 290ada5
commit 37059e4
Show file tree

Hide file tree

Showing 8 changed files with 259 additions and 12 deletions.
diff --git a/.travis.yml b/.travis.yml
@@ -67,13 +67,13 @@ jobs:
         - CXX=clang++-3.6
 
     - os: linux
-      dist: xenial
+      dist: bionic
       name: Linux clang 8
       addons:
         apt:
           sources:
             - ubuntu-toolchain-r-test
-            - llvm-toolchain-xenial-8
+            - llvm-toolchain-bionic-8
           packages:
             - clang-8
             - g++-8

diff --git a/CMakeLists.txt b/CMakeLists.txt
@@ -1,6 +1,6 @@
 cmake_minimum_required(VERSION 3.8)
 
-project(PROPOSAL VERSION 6.1.5 LANGUAGES CXX)
+project(PROPOSAL VERSION 6.1.6 LANGUAGES CXX)
 
 
 IF(I3_PROJECTS)

diff --git a/private/Interface/python/parametrization.cxx b/private/Interface/python/parametrization.cxx
@@ -787,6 +787,10 @@ void init_parametrization(py::module& m) {
 
             * RenoSarcevicSu
 
+            * AbtFT
+
+            * BlockDurandHa
+
             * AbramowiczLevinLevyMaor91Interpolant
 
             * AbramowiczLevinLevyMaor97Interpolant
@@ -795,6 +799,10 @@ void init_parametrization(py::module& m) {
 
             * RenoSarcevicSuInterpolant
 
+            * AbtFTInterpolant
+
+            * BlockDurandHaInterpolant
+
             The parametrization with "Interpolant" as a suffix creates an interpolation table for the :math:`Q^2` integration, which improved the perfomance.
 
             Example:
@@ -818,11 +826,15 @@ void init_parametrization(py::module& m) {
     PHOTO_Q2_DEF(m_sub_photo, AbramowiczLevinLevyMaor97)
     PHOTO_Q2_DEF(m_sub_photo, ButkevichMikhailov)
     PHOTO_Q2_DEF(m_sub_photo, RenoSarcevicSu)
+    PHOTO_Q2_DEF(m_sub_photo, AbtFT)
+    PHOTO_Q2_DEF(m_sub_photo, BlockDurandHa)
 
     PHOTO_Q2_INTERPOL_DEF(m_sub_photo, AbramowiczLevinLevyMaor91)
     PHOTO_Q2_INTERPOL_DEF(m_sub_photo, AbramowiczLevinLevyMaor97)
     PHOTO_Q2_INTERPOL_DEF(m_sub_photo, ButkevichMikhailov)
     PHOTO_Q2_INTERPOL_DEF(m_sub_photo, RenoSarcevicSu)
+    PHOTO_Q2_INTERPOL_DEF(m_sub_photo, AbtFT)
+    PHOTO_Q2_INTERPOL_DEF(m_sub_photo, BlockDurandHa)
 
     py::enum_<PhotonuclearFactory::Enum>(m_sub_photo, "PhotoParametrization")
         .value("Zeus", PhotonuclearFactory::Zeus)
@@ -835,6 +847,8 @@ void init_parametrization(py::module& m) {
                PhotonuclearFactory::AbramowiczLevinLevyMaor97)
         .value("ButkevichMikhailov", PhotonuclearFactory::ButkevichMikhailov)
         .value("RenoSarcevicSu", PhotonuclearFactory::RenoSarcevicSu)
+        .value("AbtFT", PhotonuclearFactory::AbtFT)
+        .value("BlockDurandHa", PhotonuclearFactory::BlockDurandHa)
         .value("None", PhotonuclearFactory::None);
 
     py::enum_<PhotonuclearFactory::Shadow>(m_sub_photo, "PhotoShadow")

diff --git a/private/PROPOSAL/crossection/factories/PhotonuclearFactory.cxx b/private/PROPOSAL/crossection/factories/PhotonuclearFactory.cxx
@@ -49,6 +49,12 @@ PhotonuclearFactory::PhotonuclearFactory()
     RegisterQ2("photorenosarcevicsu",
                RenoSarcevicSu,
                std::make_pair(&PhotoRenoSarcevicSu::create, &PhotoQ2Interpolant<PhotoRenoSarcevicSu>::create));
+    RegisterQ2("photoabtft",
+               AbtFT,
+               std::make_pair(&PhotoAbtFT::create, &PhotoQ2Interpolant<PhotoAbtFT>::create));
+    RegisterQ2("photoblockdurandha",
+               BlockDurandHa,
+               std::make_pair(&PhotoBlockDurandHa::create, &PhotoQ2Interpolant<PhotoBlockDurandHa>::create));
 }
 
 PhotonuclearFactory::~PhotonuclearFactory()

diff --git a/private/PROPOSAL/crossection/parametrization/PhotoQ2Integration.cxx b/private/PROPOSAL/crossection/parametrization/PhotoQ2Integration.cxx
@@ -139,6 +139,8 @@ Q2_PHOTO_PARAM_INTEGRAL_IMPL(AbramowiczLevinLevyMaor91)
 Q2_PHOTO_PARAM_INTEGRAL_IMPL(AbramowiczLevinLevyMaor97)
 Q2_PHOTO_PARAM_INTEGRAL_IMPL(ButkevichMikhailov)
 Q2_PHOTO_PARAM_INTEGRAL_IMPL(RenoSarcevicSu)
+Q2_PHOTO_PARAM_INTEGRAL_IMPL(AbtFT)
+Q2_PHOTO_PARAM_INTEGRAL_IMPL(BlockDurandHa)
 
 // ------------------------------------------------------------------------- //
 // Abramowicz Levin Levy Maor 91
@@ -259,8 +261,7 @@ double PhotoAbramowiczLevinLevyMaor91::FunctionToQ2Integral(double energy, doubl
     //     ( 1 - v + \frac{M x v}{2E}
     //     - \frac{v^2}{2} (1 - \frac{2m_{particle}}{Q^2})
     //          \frac{1 + \frac{4 M^2 x^2}{Q^2}}{1 + R})
-    double result = ME * RE / Q2;
-    result *= result * 4 * PI * structure_function_nucleus / v *
+    double result = std::pow(ME * RE / Q2, 2) * 4 * PI * structure_function_nucleus / v *
               (1 - v - mass_nucleus * bjorken_x * v / (2 * energy) +
                (1 - 2 * particle_def_.mass * particle_def_.mass / Q2) * v * v *
                    (1 + 4 * mass_nucleus * mass_nucleus * bjorken_x * bjorken_x / Q2) / (2 * (1 + R)));
@@ -387,8 +388,7 @@ double PhotoAbramowiczLevinLevyMaor97::FunctionToQ2Integral(double energy, doubl
     //     ( 1 - v + \frac{M x v}{2E}
     //     - \frac{v^2}{2} (1 - \frac{2m_{particle}}{Q^2})
     //          \frac{1 + \frac{4 M^2 x^2}{Q^2}}{1 + R})
-    double result = ME * RE / Q2;
-    result *= result * 4 * PI * structure_function_nucleus / v *
+    double result = std::pow(ME * RE / Q2, 2) * 4 * PI * structure_function_nucleus / v *
               (1 - v - mass_nucleus * bjorken_x * v / (2 * energy) +
                (1 - 2 * particle_def_.mass * particle_def_.mass / Q2) * v * v *
                    (1 + 4 * mass_nucleus * mass_nucleus * bjorken_x * bjorken_x / Q2) / (2 * (1 + R)));
@@ -476,8 +476,7 @@ double PhotoButkevichMikhailov::FunctionToQ2Integral(double energy, double v, do
     //     ( 1 - v + \frac{M x v}{2E}
     //     - \frac{v^2}{2} (1 - \frac{2m_{particle}}{Q^2})
     //          \frac{1 + \frac{4 M^2 x^2}{Q^2}}{1 + R})
-    double result = ME * RE / Q2;
-    result *= result * 4 * PI * structure_function_nucleus / v *
+    double result = std::pow(ME * RE / Q2, 2) * 4 * PI * structure_function_nucleus / v *
               (1 - v - mass_nucleus * bjorken_x * v / (2 * energy) +
                (1 - 2 * particle_def_.mass * particle_def_.mass / Q2) * v * v *
                    (1 + 4 * mass_nucleus * mass_nucleus * bjorken_x * bjorken_x / Q2) / (2 * (1 + R)));
@@ -607,13 +606,231 @@ double PhotoRenoSarcevicSu::FunctionToQ2Integral(double energy, double v, double
     //     ( 1 - v + \frac{v^2}{4}
     //     - \frac{v^2}{4} (1 - \frac{4m_{particle}}{Q^2})
     //          \frac{1 + \frac{4M^2x^2}{Q^2}}{1 + R})
-    double result = ME * RE / Q2;
-    result *= result * 4 * PI * structure_function_nucleus / v *
+    double result = std::pow(ME * RE / Q2, 2) * 4 * PI * structure_function_nucleus / v *
               (1 - v + 0.25 * v * v -
                (1 + 4 * particle_def_.mass * particle_def_.mass / Q2) * 0.25 * v * v *
                    (1 + 4 * mass_nucleus * mass_nucleus * bjorken_x * bjorken_x / Q2) / (1 + R));
 
     return result;
 }
 
+// ------------------------------------------------------------------------- //
+// Abt et al. HHT-ALLM-FT
+// Phys. Rev. D 96 (2017) 014001
+// ------------------------------------------------------------------------- //
+double PhotoAbtFT::FunctionToQ2Integral(double energy, double v, double Q2)
+{
+    Components::Component component = components_[component_index_];
+
+    double mass_nucleus = component.GetAverageNucleonWeight();
+
+    // Bjorken x = \frac{Q^2}{2pq}
+    double bjorken_x = Q2 / (2 * mass_nucleus * v * energy);
+
+    // --------------------------------------------------------------------- //
+    // Evaluate form factor F_2 for the nucleus same like ALLM91
+    // but using new parameters
+    // --------------------------------------------------------------------- //
+
+    // Parameter for Pomeron and Reggeon
+    const double a1_pomeron = -0.075;
+    const double a2_pomeron = -0.470;
+    const double a3_pomeron = 9.2;
+
+    const double b1_pomeron = -0.477;
+    const double b2_pomeron = 54.0;
+    const double b3_pomeron = 0.073;
+
+    const double c1_pomeron = 0.356;
+    const double c2_pomeron = 0.171;
+    const double c3_pomeron = 18.6;
+
+    const double a1_reggeon = 0.882;
+    const double a2_reggeon = 0.082;
+    const double a3_reggeon = -8.5;
+
+    const double b1_reggeon = 0.339;
+    const double b2_reggeon = 3.38;
+    const double b3_reggeon = 1.07;
+
+    const double c1_reggeon = -0.636;
+    const double c2_reggeon = 3.37;
+    const double c3_reggeon = -0.660;
+
+    // parameter with conversion from GeV^2 to Mev^2
+    const double mass_photon_eff = 0.388 * 1e6;
+    const double mass_reggeon    = 0.838 * 1e6;
+    const double mass_pomeron    = 50.8 * 1e6;
+    const double scaleParameter  = 4.4e-9 * 1e6;
+    const double Q20_free_param  = 1.87e-5 * 1e6;
+
+    // R(x, Q^2) is approximated to 0
+    // relation between structure functions F_1 and F_2
+    const double R = 0;
+
+    // t = log( \frac{log \frac{Q^2 + Q_0^2}{\Lambda^2}}{log \frac{Q_0^2}{\Lambda^2}} )
+    // eq. 22
+    double t = std::log(std::log((Q2 + Q20_free_param) / scaleParameter) / std::log(Q20_free_param / scaleParameter));
+
+    if (t < 0)
+        t = 0;
+
+    // parameter, that increase with Q^2
+    // eq. 23
+    // f(t) = f_1 + f_2 t^{f_3}
+    double b_reggeon = b1_reggeon + b2_reggeon * std::pow(t, b3_reggeon);
+    double b_pomeron = b1_pomeron + b2_pomeron * std::pow(t, b3_pomeron);
+    double a_reggeon;
+    double c_reggeon;
+    // handle the case if t = 0
+    // since a3_reggeon and c3_reggeon are < 0
+    // resulting in a_reggeon and c_reggeon being inf
+    if (t>0)
+    {
+        a_reggeon = a1_reggeon + a2_reggeon * std::pow(t, a3_reggeon);
+        c_reggeon = c1_reggeon + c2_reggeon * std::pow(t, c3_reggeon);
+    } else {
+        a_reggeon = a1_reggeon;
+        c_reggeon = c1_reggeon;
+    }
+
+    // parameter, that decrease with Q^2
+    // eq. 24
+    // f'(t) = f_1 + (f_1 - f_2) (\frac{1}{1 + t^{f_3}} - 1)
+    double a_pomeron = a1_pomeron + (a1_pomeron - a2_pomeron) * (1 / (1 + std::pow(t, a3_pomeron)) - 1);
+    double c_pomeron = c1_pomeron + (c1_pomeron - c2_pomeron) * (1 / (1 + std::pow(t, c3_pomeron)) - 1);
+
+    // invariant mass of nucleus and virtual photon
+    // W^2 = (p + q)^2 = M^2 + 2MEv - Q^2
+    double W2 = mass_nucleus * mass_nucleus + 2 * mass_nucleus * energy * v - Q2;
+
+    // Relation between structure function of proton and structure function of neutron
+    // from the BCDMS Collaboration
+    // Phys. Lett. B 237 (1990), 599
+    // eq. 4
+    // P(x) = 1 - 1.85x + 2.45x^2 - 2.35x^3 + x^4
+    double relation_proton_neutron = bjorken_x * bjorken_x;
+    relation_proton_neutron        = 1 - 1.85 * bjorken_x + 2.45 * relation_proton_neutron -
+                              2.35 * relation_proton_neutron * bjorken_x +
+                              relation_proton_neutron * relation_proton_neutron;
+
+    // eq. 17 and 18
+    // x_i = \frac{Q^2 + m_i}{Q^2 + m_i + W^2 - M^2}
+    double bjorken_x_pomeron = (Q2 + mass_pomeron) / (Q2 + mass_pomeron + W2 - mass_nucleus * mass_nucleus);
+    double bjorken_x_reggeon = (Q2 + mass_reggeon) / (Q2 + mass_reggeon + W2 - mass_nucleus * mass_nucleus);
+    // eq. 20 and 21
+    // F_{2, i}(x, Q^2) = c_i(t) x^{a_i(t)} (1 - x)^{b_i(t)}
+    double pomeron_contribution = c_pomeron * std::pow(bjorken_x_pomeron, a_pomeron) * std::pow(1 - bjorken_x, b_pomeron);
+    double reggeon_controbution = c_reggeon * std::pow(bjorken_x_reggeon, a_reggeon) * std::pow(1 - bjorken_x, b_reggeon);
+    // ALLM97 eq. 1
+    // F_{2, proton}(x, Q^2) = \frac{Q^2}{Q^2 + m_0^2} (F_{2, Pomeron} + F_{2, Reggeon})
+    double structure_function_proton = Q2 / (Q2 + mass_photon_eff) * (pomeron_contribution + reggeon_controbution);
+    // structure function from Dutta et al
+    // Phys. Rev. D 63 (2001), 094020
+    // eq. 3.11
+    // F_{2, nucleus} = G(x) (Z + (A - Z)P(x)) F_{2, Proton}
+    double structure_function_nucleus =
+        structure_function_proton * shadow_effect_->CalculateShadowEffect(component, bjorken_x, v * energy) *
+        (component.GetNucCharge() + (component.GetAtomicNum() - component.GetNucCharge()) * relation_proton_neutron);
+
+    // differential cross section from Dutta et al.
+    // Phys. Rev. D 63 (2001), 094020
+    // eq. 3.11
+    // eq. 3.4
+    // \frac{d^2 \sigma}{dQ^2 dx} = \frac{4\pi\alpha^2}{Q^4} \frac{F_{2, Proton}}{v}
+    //     ( 1 - v + \frac{M x v}{2E}
+    //     - \frac{v^2}{2} (1 - \frac{2m_{particle}}{Q^2})
+    //          \frac{1 + \frac{4 M^2 x^2}{Q^2}}{1 + R})
+    double result = std::pow(ME * RE / Q2, 2) * 4 * PI * structure_function_nucleus / v *
+              (1 - v - mass_nucleus * bjorken_x * v / (2 * energy) +
+               (1 - 2 * particle_def_.mass * particle_def_.mass / Q2) * v * v *
+                   (1 + 4 * mass_nucleus * mass_nucleus * bjorken_x * bjorken_x / Q2) / (2 * (1 + R)));
+
+    return result;
+}
+
+// ------------------------------------------------------------------------- //
+// Martin M. Block, Loyal Durand, Phuoc Ha
+// Phys. Rev. D 89 (2014) 094027
+// ------------------------------------------------------------------------- //
+double PhotoBlockDurandHa::FunctionToQ2Integral(double energy, double v,
+  double Q2)
+{
+
+    Components::Component component = components_[component_index_];
+
+    double mass_nucleus = component.GetAverageNucleonWeight();
+
+    // Bjorken x = \frac{Q^2}{2pq}
+    double bjorken_x = Q2 / (2 * mass_nucleus * v * energy);
+
+    // Parameters of BDH fit from Table I
+    // Dimensionless parameter values
+    const double a0 = 8.205e-4;
+    const double a1 = -5.148e-2;
+    const double a2 = -4.725e-3;
+
+    const double b0 = 2.217e-3;
+    const double b1 = 1.244e-2;
+    const double b2 = 5.958e-4;
+
+    const double c0 = 0.255;
+    const double c1 = 1.475e-1;
+
+    const double lambda = 2.430;
+    const double n      = 11.49;
+
+    // parameters with Conversion from GeV^2 to MeV^2
+    const double M2  = 0.753 * 1e6;
+    const double mu2 = 2.82  * 1e6;
+
+    // R(x, Q^2) is approximated to 0
+    // relation between structure functions F_1 and F_2
+    const double R = 0;
+
+    // eq. 8, 10
+    double aux = std::log(1 + Q2/mu2);
+    double A   = a0 + a1 * aux + a2 * aux * aux;
+    double B   = b0 + b1 * aux + b2 * aux * aux;
+    double C   = c0 + c1 * aux;
+
+    double D = Q2 * (Q2 + lambda * M2)/std::pow((Q2 + M2), 2);
+
+    double aux1 = std::log(Q2/bjorken_x/(Q2 + mu2));
+    double structure_function_proton = D * std::pow(1 - bjorken_x, n)
+      * (C + A * aux1 + B * aux1 * aux1);
+
+    // Relation between structure function of proton and structure function of neutron
+    // from the BCDMS Collaboration
+    // Phys. Lett. B 237 (1990), 599
+    // eq. 4
+    // P(x) = 1 - 1.85x + 2.45x^2 - 2.35x^3 + x^4
+    double relation_proton_neutron = bjorken_x * bjorken_x;
+    relation_proton_neutron        = 1 - 1.85 * bjorken_x + 2.45 * relation_proton_neutron -
+                              2.35 * relation_proton_neutron * bjorken_x +
+                              relation_proton_neutron * relation_proton_neutron;
+
+    // structure function from Dutta et al
+    // Phys. Rev. D 63 (2001), 094020
+    // eq. 3.11
+    // F_{2, nucleus} = G(x) (Z + (A - Z)P(x)) F_{2, Proton}
+    double structure_function_nucleus =
+        structure_function_proton * shadow_effect_->CalculateShadowEffect(component, bjorken_x, v * energy) *
+        (component.GetNucCharge() + (component.GetAtomicNum() - component.GetNucCharge()) * relation_proton_neutron);
+
+    // differential cross section from Dutta et al.
+    // Phys. Rev. D 63 (2001), 094020
+    // eq. 3.11
+    // eq. 3.4
+    // \frac{d^2 \sigma}{dQ^2 dx} = \frac{4\pi\alpha^2}{Q^4} \frac{F_{2, Proton}}{v}
+    //     ( 1 - v + \frac{M x v}{2E}
+    //     - \frac{v^2}{2} (1 - \frac{2m_{particle}}{Q^2})
+    //          \frac{1 + \frac{4 M^2 x^2}{Q^2}}{1 + R})
+    double result = std::pow(ME * RE / Q2, 2) * 4 * PI * structure_function_nucleus / v *
+              (1 - v - mass_nucleus * bjorken_x * v / (2 * energy) +
+               (1 - 2 * particle_def_.mass * particle_def_.mass / Q2) * v * v *
+                   (1 + 4 * mass_nucleus * mass_nucleus * bjorken_x * bjorken_x / Q2) / (2 * (1 + R)));
+
+    return result;
+}
 #undef Q2_PHOTO_PARAM_INTEGRAL_IMPL
diff --git a/public/PROPOSAL/crossection/factories/PhotonuclearFactory.h b/public/PROPOSAL/crossection/factories/PhotonuclearFactory.h
@@ -68,6 +68,8 @@ class PhotonuclearFactory
         AbramowiczLevinLevyMaor97,
         ButkevichMikhailov,
         RenoSarcevicSu,
+        AbtFT,
+        BlockDurandHa,
     };
 
     enum Shadow

diff --git a/public/PROPOSAL/crossection/parametrization/PhotoQ2Integration.h b/public/PROPOSAL/crossection/parametrization/PhotoQ2Integration.h
@@ -120,6 +120,8 @@ Q2_PHOTO_PARAM_INTEGRAL_DEC(AbramowiczLevinLevyMaor91)
 Q2_PHOTO_PARAM_INTEGRAL_DEC(AbramowiczLevinLevyMaor97)
 Q2_PHOTO_PARAM_INTEGRAL_DEC(ButkevichMikhailov)
 Q2_PHOTO_PARAM_INTEGRAL_DEC(RenoSarcevicSu)
+Q2_PHOTO_PARAM_INTEGRAL_DEC(AbtFT)
+Q2_PHOTO_PARAM_INTEGRAL_DEC(BlockDurandHa)
 
 /******************************************************************************
  *                    Declare Interpolant Parametrizations                    *

diff --git a/resources/examples/standalone/plot_photo.py b/resources/examples/standalone/plot_photo.py
@@ -43,6 +43,12 @@
         ),
         parametrization.photonuclear.ButkevichMikhailovInterpolant(
             *param_defs
+        ),
+        parametrization.photonuclear.AbtFTInterpolant(
+            *param_defs
+        ),
+        parametrization.photonuclear.BlockDurandHaInterpolant(
+            *param_defs
         )
     ]
 
@@ -58,7 +64,7 @@
             interpolation_def
         ))
 
-    print(crosssections[0])
+    # print(crosssections[0])
 
     # =========================================================
     # 	Calculate DE/dx at the given energies