From 25580181e2db196ffdf2810e821265f668abf8fd Mon Sep 17 00:00:00 2001 From: Maximilian Reininghaus <maximilian.reininghaus@kit.edu> Date: Fri, 25 Sep 2020 15:52:27 +0200 Subject: [PATCH] adapted tests to DoInteraction(SecondaryView) --- Documentation/Examples/vertical_EAS.cc | 5 +- .../HadronicElasticModel.cc | 307 +++++++++--------- .../testInteractionCounter.cc | 6 +- Processes/Pythia/testPythia8.cc | 6 +- Processes/QGSJetII/testQGSJetII.cc | 2 +- Processes/Sibyll/testSibyll.cc | 12 +- Processes/UrQMD/testUrQMD.cc | 6 +- 7 files changed, 164 insertions(+), 180 deletions(-) diff --git a/Documentation/Examples/vertical_EAS.cc b/Documentation/Examples/vertical_EAS.cc index 88e9f172b..77c02ada9 100644 --- a/Documentation/Examples/vertical_EAS.cc +++ b/Documentation/Examples/vertical_EAS.cc @@ -225,10 +225,7 @@ int main(int argc, char** argv) { // define air shower object, run simulation tracking_line::TrackingLine tracking; - cascade::Cascade<decltype(tracking), decltype(sequence), decltype(stack) - //~ , corsika::history::HistorySecondaryView<corsika::setup::StackView> - > - EAS(env, tracking, sequence, stack); + cascade::Cascade EAS(env, tracking, sequence, stack); // to fix the point of first interaction, uncomment the following two lines: // EAS.SetNodes(); diff --git a/Processes/HadronicElasticModel/HadronicElasticModel.cc b/Processes/HadronicElasticModel/HadronicElasticModel.cc index ece6447e8..c8f3eed3d 100644 --- a/Processes/HadronicElasticModel/HadronicElasticModel.cc +++ b/Processes/HadronicElasticModel/HadronicElasticModel.cc @@ -6,11 +6,10 @@ * the license. */ -#include <corsika/process/hadronic_elastic_model/HadronicElasticModel.h> - #include <corsika/environment/Environment.h> #include <corsika/environment/NuclearComposition.h> #include <corsika/geometry/FourVector.h> +#include <corsika/process/hadronic_elastic_model/HadronicElasticModel.h> #include <corsika/random/ExponentialDistribution.h> #include <corsika/utl/COMBoost.h> @@ -19,181 +18,177 @@ #include <iomanip> #include <iostream> -using namespace corsika::setup; using SetupParticle = corsika::setup::Stack::ParticleType; -using SetupView = corsika::setup::Stack::StackView; - -namespace corsika::process::HadronicElasticModel { - - HadronicElasticInteraction::HadronicElasticInteraction(units::si::CrossSectionType x, - units::si::CrossSectionType y) - : fX(x) - , fY(y) {} - - template <> - units::si::GrammageType HadronicElasticInteraction::GetInteractionLength( - SetupParticle const& p) { - using namespace units::si; - if (p.GetPID() == particles::Code::Proton) { - auto const* currentNode = p.GetNode(); - auto const& mediumComposition = - currentNode->GetModelProperties().GetNuclearComposition(); - - auto const& components = mediumComposition.GetComponents(); - auto const& fractions = mediumComposition.GetFractions(); +using SetupView = corsika::setup::StackView; - auto const projectileMomentum = p.GetMomentum(); - auto const projectileMomentumSquaredNorm = projectileMomentum.squaredNorm(); - auto const projectileEnergy = p.GetEnergy(); +using namespace corsika::process::HadronicElasticModel; +using namespace corsika; - auto const avgCrossSection = [&]() { - CrossSectionType avgCrossSection = 0_b; +HadronicElasticInteraction::HadronicElasticInteraction(units::si::CrossSectionType x, + units::si::CrossSectionType y) + : fX(x) + , fY(y) {} - for (size_t i = 0; i < fractions.size(); ++i) { - auto const targetMass = particles::GetMass(components[i]); - auto const s = detail::static_pow<2>(projectileEnergy + targetMass) - - projectileMomentumSquaredNorm; - avgCrossSection += CrossSection(s) * fractions[i]; - } - - std::cout << "avgCrossSection: " << avgCrossSection / 1_mb << " mb" << std::endl; - - return avgCrossSection; - }(); - - auto const avgTargetMassNumber = mediumComposition.GetAverageMassNumber(); - - GrammageType const interactionLength = - avgTargetMassNumber * units::constants::u / avgCrossSection; - - return interactionLength; - } else { - return std::numeric_limits<double>::infinity() * 1_g / (1_cm * 1_cm); - } - } +template <> +units::si::GrammageType HadronicElasticInteraction::GetInteractionLength( + SetupParticle const& p) { + using namespace units::si; + if (p.GetPID() == particles::Code::Proton) { + auto const* currentNode = p.GetNode(); + auto const& mediumComposition = + currentNode->GetModelProperties().GetNuclearComposition(); - template <> - process::EProcessReturn HadronicElasticInteraction::DoInteraction(SetupView& view) { - if (p.GetPID() != particles::Code::Proton) { return process::EProcessReturn::eOk; } - - using namespace units::si; - using namespace units::constants; - - auto const p = view.GetProjectile(); - - const auto* currentNode = p.GetNode(); - const auto& composition = currentNode->GetModelProperties().GetNuclearComposition(); - const auto& components = composition.GetComponents(); - - std::vector<units::si::CrossSectionType> cross_section_of_components( - composition.GetComponents().size()); + auto const& components = mediumComposition.GetComponents(); + auto const& fractions = mediumComposition.GetFractions(); auto const projectileMomentum = p.GetMomentum(); auto const projectileMomentumSquaredNorm = projectileMomentum.squaredNorm(); auto const projectileEnergy = p.GetEnergy(); - for (size_t i = 0; i < components.size(); ++i) { - auto const targetMass = particles::GetMass(components[i]); - auto const s = units::si::detail::static_pow<2>(projectileEnergy + targetMass) - - projectileMomentumSquaredNorm; - cross_section_of_components[i] = CrossSection(s); - } - - const auto targetCode = composition.SampleTarget(cross_section_of_components, fRNG); - - auto const targetMass = particles::GetMass(targetCode); + auto const avgCrossSection = [&]() { + CrossSectionType avgCrossSection = 0_b; - std::uniform_real_distribution phiDist(0., 2 * M_PI); + for (size_t i = 0; i < fractions.size(); ++i) { + auto const targetMass = particles::GetMass(components[i]); + auto const s = detail::static_pow<2>(projectileEnergy + targetMass) - + projectileMomentumSquaredNorm; + avgCrossSection += CrossSection(s) * fractions[i]; + } - geometry::FourVector const projectileLab(projectileEnergy, projectileMomentum); - geometry::FourVector const targetLab( - targetMass, geometry::Vector<units::si::hepmomentum_d>( - projectileMomentum.GetCoordinateSystem(), {0_eV, 0_eV, 0_eV})); - utl::COMBoost const boost(projectileLab, targetMass); + std::cout << "avgCrossSection: " << avgCrossSection / 1_mb << " mb" << std::endl; - auto const projectileCoM = boost.toCoM(projectileLab); - auto const targetCoM = boost.toCoM(targetLab); - - auto const pProjectileCoMSqNorm = - projectileCoM.GetSpaceLikeComponents().squaredNorm(); - auto const pProjectileCoMNorm = sqrt(pProjectileCoMSqNorm); + return avgCrossSection; + }(); - auto const eProjectileCoM = projectileCoM.GetTimeLikeComponent(); - auto const eTargetCoM = targetCoM.GetTimeLikeComponent(); + auto const avgTargetMassNumber = mediumComposition.GetAverageMassNumber(); - auto const sqrtS = eProjectileCoM + eTargetCoM; - auto const s = units::si::detail::static_pow<2>(sqrtS); + GrammageType const interactionLength = + avgTargetMassNumber * units::constants::u / avgCrossSection; - auto const B = this->B(s); - std::cout << B << std::endl; + return interactionLength; + } else { + return std::numeric_limits<double>::infinity() * 1_g / (1_cm * 1_cm); + } +} - random::ExponentialDistribution tDist(1 / B); - auto const absT = [&]() { - decltype(tDist(fRNG)) absT; - auto const maxT = 4 * pProjectileCoMSqNorm; +template <> +process::EProcessReturn HadronicElasticInteraction::DoInteraction(SetupView& view) { + using namespace units::si; + using namespace units::constants; - do { - // |t| cannot become arbitrarily large, max. given by GER eq. (4.16), so we just - // throw again until we have an acceptable value. Note that the formula holds in - // any frame despite of what is stated in the book. - absT = tDist(fRNG); - } while (absT >= maxT); + auto p = view.GetProjectile(); + if (p.GetPID() != particles::Code::Proton) { return process::EProcessReturn::eOk; } - return absT; - }(); + const auto* currentNode = p.GetNode(); + const auto& composition = currentNode->GetModelProperties().GetNuclearComposition(); + const auto& components = composition.GetComponents(); - std::cout << "HadronicElasticInteraction: s = " << s * invGeVsq - << " GeV²; absT = " << absT * invGeVsq - << " GeV² (max./GeV² = " << 4 * invGeVsq * projectileMomentumSquaredNorm - << ')' << std::endl; - - auto const theta = 2 * asin(sqrt(absT / (4 * pProjectileCoMSqNorm))); - auto const phi = phiDist(fRNG); - - auto const projectileScatteredLab = boost.fromCoM( - geometry::FourVector<HEPEnergyType, geometry::Vector<hepmomentum_d>>( - eProjectileCoM, - geometry::Vector<hepmomentum_d>(projectileMomentum.GetCoordinateSystem(), - {pProjectileCoMNorm * sin(theta) * cos(phi), - pProjectileCoMNorm * sin(theta) * sin(phi), - pProjectileCoMNorm * cos(theta)}))); - - p.SetMomentum(projectileScatteredLab.GetSpaceLikeComponents()); - p.SetEnergy( - sqrt(projectileScatteredLab.GetSpaceLikeComponents().squaredNorm() + - units::si::detail::static_pow<2>(particles::GetMass( - p.GetPID())))); // Don't use energy from boost. It can be smaller than - // the momentum due to limited numerical accuracy. - - return process::EProcessReturn::eOk; - } + std::vector<units::si::CrossSectionType> cross_section_of_components( + composition.GetComponents().size()); - HadronicElasticInteraction::inveV2 HadronicElasticInteraction::B(eV2 s) const { - using namespace units::constants; - auto constexpr b_p = 2.3; - auto const result = - (2 * b_p + 2 * b_p + 4 * pow(s * invGeVsq, gfEpsilon) - 4.2) * invGeVsq; - std::cout << "B(" << s << ") = " << result / invGeVsq << " GeV¯²" << std::endl; - return result; - } + auto const projectileMomentum = p.GetMomentum(); + auto const projectileMomentumSquaredNorm = projectileMomentum.squaredNorm(); + auto const projectileEnergy = p.GetEnergy(); - units::si::CrossSectionType HadronicElasticInteraction::CrossSection( - SquaredHEPEnergyType s) const { - using namespace units::si; - using namespace units::constants; - // assuming every target behaves like a proton, fX and fY are universal - CrossSectionType const sigmaTotal = - fX * pow(s * invGeVsq, gfEpsilon) + fY * pow(s * invGeVsq, -gfEta); - - // according to Schuler & Sjöstrand, PRD 49, 2257 (1994) - // (we ignore rho because rho^2 is just ~2 %) - auto const sigmaElastic = - units::si::detail::static_pow<2>(sigmaTotal) / - (16 * M_PI * ConvertHEPToSI<CrossSectionType::dimension_type>(B(s))); - - std::cout << "HEM sigmaTot = " << sigmaTotal / 1_mb << " mb" << std::endl; - std::cout << "HEM sigmaElastic = " << sigmaElastic / 1_mb << " mb" << std::endl; - return sigmaElastic; + for (size_t i = 0; i < components.size(); ++i) { + auto const targetMass = particles::GetMass(components[i]); + auto const s = units::si::detail::static_pow<2>(projectileEnergy + targetMass) - + projectileMomentumSquaredNorm; + cross_section_of_components[i] = CrossSection(s); } -} // namespace corsika::process::HadronicElasticModel + const auto targetCode = composition.SampleTarget(cross_section_of_components, fRNG); + + auto const targetMass = particles::GetMass(targetCode); + + std::uniform_real_distribution phiDist(0., 2 * M_PI); + + geometry::FourVector const projectileLab(projectileEnergy, projectileMomentum); + geometry::FourVector const targetLab( + targetMass, geometry::Vector<units::si::hepmomentum_d>( + projectileMomentum.GetCoordinateSystem(), {0_eV, 0_eV, 0_eV})); + utl::COMBoost const boost(projectileLab, targetMass); + + auto const projectileCoM = boost.toCoM(projectileLab); + auto const targetCoM = boost.toCoM(targetLab); + + auto const pProjectileCoMSqNorm = projectileCoM.GetSpaceLikeComponents().squaredNorm(); + auto const pProjectileCoMNorm = sqrt(pProjectileCoMSqNorm); + + auto const eProjectileCoM = projectileCoM.GetTimeLikeComponent(); + auto const eTargetCoM = targetCoM.GetTimeLikeComponent(); + + auto const sqrtS = eProjectileCoM + eTargetCoM; + auto const s = units::si::detail::static_pow<2>(sqrtS); + + auto const B = this->B(s); + std::cout << B << std::endl; + + random::ExponentialDistribution tDist(1 / B); + auto const absT = [&]() { + decltype(tDist(fRNG)) absT; + auto const maxT = 4 * pProjectileCoMSqNorm; + + do { + // |t| cannot become arbitrarily large, max. given by GER eq. (4.16), so we just + // throw again until we have an acceptable value. Note that the formula holds in + // any frame despite of what is stated in the book. + absT = tDist(fRNG); + } while (absT >= maxT); + + return absT; + }(); + + std::cout << "HadronicElasticInteraction: s = " << s * invGeVsq + << " GeV²; absT = " << absT * invGeVsq + << " GeV² (max./GeV² = " << 4 * invGeVsq * projectileMomentumSquaredNorm + << ')' << std::endl; + + auto const theta = 2 * asin(sqrt(absT / (4 * pProjectileCoMSqNorm))); + auto const phi = phiDist(fRNG); + + auto const projectileScatteredLab = + boost.fromCoM(geometry::FourVector<HEPEnergyType, geometry::Vector<hepmomentum_d>>( + eProjectileCoM, + geometry::Vector<hepmomentum_d>(projectileMomentum.GetCoordinateSystem(), + {pProjectileCoMNorm * sin(theta) * cos(phi), + pProjectileCoMNorm * sin(theta) * sin(phi), + pProjectileCoMNorm * cos(theta)}))); + + p.SetMomentum(projectileScatteredLab.GetSpaceLikeComponents()); + p.SetEnergy( + sqrt(projectileScatteredLab.GetSpaceLikeComponents().squaredNorm() + + units::si::detail::static_pow<2>(particles::GetMass( + p.GetPID())))); // Don't use energy from boost. It can be smaller than + // the momentum due to limited numerical accuracy. + + return process::EProcessReturn::eOk; +} + +HadronicElasticInteraction::inveV2 HadronicElasticInteraction::B(eV2 s) const { + using namespace units::constants; + auto constexpr b_p = 2.3; + auto const result = + (2 * b_p + 2 * b_p + 4 * pow(s * invGeVsq, gfEpsilon) - 4.2) * invGeVsq; + std::cout << "B(" << s << ") = " << result / invGeVsq << " GeV¯²" << std::endl; + return result; +} + +units::si::CrossSectionType HadronicElasticInteraction::CrossSection( + SquaredHEPEnergyType s) const { + using namespace units::si; + using namespace units::constants; + // assuming every target behaves like a proton, fX and fY are universal + CrossSectionType const sigmaTotal = + fX * pow(s * invGeVsq, gfEpsilon) + fY * pow(s * invGeVsq, -gfEta); + + // according to Schuler & Sjöstrand, PRD 49, 2257 (1994) + // (we ignore rho because rho^2 is just ~2 %) + auto const sigmaElastic = + units::si::detail::static_pow<2>(sigmaTotal) / + (16 * M_PI * ConvertHEPToSI<CrossSectionType::dimension_type>(B(s))); + + std::cout << "HEM sigmaTot = " << sigmaTotal / 1_mb << " mb" << std::endl; + std::cout << "HEM sigmaElastic = " << sigmaElastic / 1_mb << " mb" << std::endl; + return sigmaElastic; +} diff --git a/Processes/InteractionCounter/testInteractionCounter.cc b/Processes/InteractionCounter/testInteractionCounter.cc index 13bf34c17..c63579e88 100644 --- a/Processes/InteractionCounter/testInteractionCounter.cc +++ b/Processes/InteractionCounter/testInteractionCounter.cc @@ -124,8 +124,7 @@ TEST_CASE("InteractionCounter") { REQUIRE(stackPtr->getEntries() == 1); REQUIRE(secViewPtr->getEntries() == 0); - auto projectile = secViewPtr->GetProjectile(); - auto const ret = countedProcess.DoInteraction(projectile); + auto const ret = countedProcess.DoInteraction(*secViewPtr); REQUIRE(ret == nullptr); auto const& h = countedProcess.GetHistogram().labHists().second.at(1'000'070'140); @@ -144,8 +143,7 @@ TEST_CASE("InteractionCounter") { REQUIRE(stackPtr->getEntries() == 1); REQUIRE(secViewPtr->getEntries() == 0); - auto projectile = secViewPtr->GetProjectile(); - auto const ret = countedProcess.DoInteraction(projectile); + auto const ret = countedProcess.DoInteraction(*secViewPtr); REQUIRE(ret == nullptr); auto const& h = countedProcess.GetHistogram().labHists().first; diff --git a/Processes/Pythia/testPythia8.cc b/Processes/Pythia/testPythia8.cc index 680b8cade..4b1da6d83 100644 --- a/Processes/Pythia/testPythia8.cc +++ b/Processes/Pythia/testPythia8.cc @@ -131,12 +131,11 @@ TEST_CASE("pythia process") { random::RNGManager::GetInstance().RegisterRandomStream("pythia"); corsika::stack::SecondaryView view(particle); - auto projectile = view.GetProjectile(); process::pythia::Decay model; [[maybe_unused]] const TimeType time = model.GetLifetime(particle); - model.DoDecay(projectile); + model.DoDecay(view); CHECK(stack.getEntries() == 3); auto const pSum = sumMomentum(view, cs); CHECK((pSum - plab).norm() / 1_GeV == Approx(0).margin(1e-4)); @@ -181,11 +180,10 @@ TEST_CASE("pythia process") { particles::Code::PiPlus, E0, plab, pos, 0_ns}); particle.SetNode(nodePtr); corsika::stack::SecondaryView view(particle); - auto projectile = view.GetProjectile(); process::pythia::Interaction model; - [[maybe_unused]] const process::EProcessReturn ret = model.DoInteraction(projectile); + [[maybe_unused]] const process::EProcessReturn ret = model.DoInteraction(view); [[maybe_unused]] const GrammageType length = model.GetInteractionLength(particle); } } diff --git a/Processes/QGSJetII/testQGSJetII.cc b/Processes/QGSJetII/testQGSJetII.cc index 6daa09512..0e6fddae9 100644 --- a/Processes/QGSJetII/testQGSJetII.cc +++ b/Processes/QGSJetII/testQGSJetII.cc @@ -139,7 +139,7 @@ TEST_CASE("QgsjetIIInterface", "[processes]") { Interaction model; - [[maybe_unused]] const process::EProcessReturn ret = model.DoInteraction(projectile); + [[maybe_unused]] const process::EProcessReturn ret = model.DoInteraction(view); [[maybe_unused]] const GrammageType length = model.GetInteractionLength(particle); CHECK(length / (1_g / square(1_cm)) == Approx(93.47).margin(0.1)); diff --git a/Processes/Sibyll/testSibyll.cc b/Processes/Sibyll/testSibyll.cc index bd7454856..e99400670 100644 --- a/Processes/Sibyll/testSibyll.cc +++ b/Processes/Sibyll/testSibyll.cc @@ -131,11 +131,10 @@ TEST_CASE("SibyllInterface", "[processes]") { particles::Code::Proton, E0, plab, pos, 0_ns}); particle.SetNode(nodePtr); corsika::stack::SecondaryView view(particle); - auto projectile = view.GetProjectile(); Interaction model; - [[maybe_unused]] const process::EProcessReturn ret = model.DoInteraction(projectile); + [[maybe_unused]] const process::EProcessReturn ret = model.DoInteraction(view); auto const pSum = sumMomentum(view, cs); /* @@ -217,11 +216,10 @@ TEST_CASE("SibyllInterface", "[processes]") { particles::Code::Proton, E0, plab, pos, 0_ns}); particle.SetNode(nodePtr); corsika::stack::SecondaryView view(particle); - auto projectile = view.GetProjectile(); Interaction model; - [[maybe_unused]] const process::EProcessReturn ret = model.DoInteraction(projectile); + [[maybe_unused]] const process::EProcessReturn ret = model.DoInteraction(view); auto const pSum = sumMomentum(view, cs); CHECK(pSum.GetComponents(cs).GetX() / P0 == Approx(1).margin(0.001)); CHECK(pSum.GetComponents(cs).GetY() / 1_GeV == Approx(0).margin(1e-4)); @@ -248,12 +246,11 @@ TEST_CASE("SibyllInterface", "[processes]") { particles::Code::Nucleus, E0, plab, pos, 0_ns, 4, 2}); particle.SetNode(nodePtr); corsika::stack::SecondaryView view(particle); - auto projectile = view.GetProjectile(); Interaction hmodel; NuclearInteraction model(hmodel, env); - [[maybe_unused]] const process::EProcessReturn ret = model.DoInteraction(projectile); + [[maybe_unused]] const process::EProcessReturn ret = model.DoInteraction(view); [[maybe_unused]] const GrammageType length = model.GetInteractionLength(particle); } @@ -270,7 +267,6 @@ TEST_CASE("SibyllInterface", "[processes]") { corsika::stack::MomentumVector, geometry::Point, units::si::TimeType>{ particles::Code::Lambda0, E0, plab, pos, 0_ns}); corsika::stack::SecondaryView view(particle); - auto projectile = view.GetProjectile(); Decay model; @@ -278,7 +274,7 @@ TEST_CASE("SibyllInterface", "[processes]") { [[maybe_unused]] const TimeType time = model.GetLifetime(particle); - /*[[maybe_unused]] const process::EProcessReturn ret =*/model.DoDecay(projectile); + /*[[maybe_unused]] const process::EProcessReturn ret =*/model.DoDecay(view); // run checks // lambda decays into proton and pi- or neutron and pi+ diff --git a/Processes/UrQMD/testUrQMD.cc b/Processes/UrQMD/testUrQMD.cc index 11f11adfd..f7d36c891 100644 --- a/Processes/UrQMD/testUrQMD.cc +++ b/Processes/UrQMD/testUrQMD.cc @@ -168,7 +168,7 @@ TEST_CASE("UrQMD") { // must be assigned to variable, cannot be used as rvalue?! auto projectile = secViewPtr->GetProjectile(); auto const projectileMomentum = projectile.GetMomentum(); - [[maybe_unused]] process::EProcessReturn const ret = urqmd.DoInteraction(projectile); + [[maybe_unused]] process::EProcessReturn const ret = urqmd.DoInteraction(*secViewPtr); REQUIRE(sumCharge(*secViewPtr) == Z + particles::GetChargeNumber(particles::Code::Oxygen)); @@ -193,7 +193,7 @@ TEST_CASE("UrQMD") { auto projectile = secViewPtr->GetProjectile(); auto const projectileMomentum = projectile.GetMomentum(); - [[maybe_unused]] process::EProcessReturn const ret = urqmd.DoInteraction(projectile); + [[maybe_unused]] process::EProcessReturn const ret = urqmd.DoInteraction(*secViewPtr); REQUIRE(sumCharge(*secViewPtr) == particles::GetChargeNumber(particles::Code::PiPlus) + @@ -219,7 +219,7 @@ TEST_CASE("UrQMD") { auto projectile = secViewPtr->GetProjectile(); auto const projectileMomentum = projectile.GetMomentum(); - [[maybe_unused]] process::EProcessReturn const ret = urqmd.DoInteraction(projectile); + [[maybe_unused]] process::EProcessReturn const ret = urqmd.DoInteraction(*secViewPtr); REQUIRE(sumCharge(*secViewPtr) == particles::GetChargeNumber(particles::Code::K0Long) + -- GitLab