diff --git a/.travis.yml b/.travis.yml index 9eecda0b9..84bda5371 100644 --- 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 index ff475d9a2..3bdfb2e84 100644 --- 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 index 1ead04df9..80c33f6e6 100644 --- 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_(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_(m_sub_photo, "PhotoShadow") diff --git a/private/PROPOSAL/crossection/factories/PhotonuclearFactory.cxx b/private/PROPOSAL/crossection/factories/PhotonuclearFactory.cxx index a11e13e6f..e4e71b116 100644 --- 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::create)); + RegisterQ2("photoabtft", + AbtFT, + std::make_pair(&PhotoAbtFT::create, &PhotoQ2Interpolant::create)); + RegisterQ2("photoblockdurandha", + BlockDurandHa, + std::make_pair(&PhotoBlockDurandHa::create, &PhotoQ2Interpolant::create)); } PhotonuclearFactory::~PhotonuclearFactory() diff --git a/private/PROPOSAL/crossection/parametrization/PhotoQ2Integration.cxx b/private/PROPOSAL/crossection/parametrization/PhotoQ2Integration.cxx index 5c19e293b..c7eae9bed 100644 --- 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,8 +606,7 @@ 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)); @@ -616,4 +614,223 @@ double PhotoRenoSarcevicSu::FunctionToQ2Integral(double energy, double v, double 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 index 8e5a8dfff..0bc2f4842 100644 --- 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 index b580d6ba4..080382b14 100644 --- 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 index 511dbcd43..bf46cc5f5 100644 --- 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