IAP GITLAB

Commit 9dece868 authored by Maximilian Reininghaus's avatar Maximilian Reininghaus 🖖 Committed by Maximilian Reininghaus

refactored interaction interface without track argument

parent eaa5461b
......@@ -153,7 +153,7 @@ namespace corsika::cascade {
// determine combined total interaction length (inverse)
InverseGrammageType const total_inv_lambda =
fProcessSequence.GetTotalInverseInteractionLength(vParticle, step);
fProcessSequence.GetTotalInverseInteractionLength(vParticle);
// sample random exponential step length in grammage
corsika::random::ExponentialDistribution expDist(1 / total_inv_lambda);
......@@ -242,13 +242,14 @@ namespace corsika::cascade {
if (min_distance == distance_interact) {
std::cout << "collide" << std::endl;
InverseGrammageType const actual_inv_length =
fProcessSequence.GetTotalInverseInteractionLength(vParticle, step);
InverseGrammageType const current_inv_length =
fProcessSequence.GetTotalInverseInteractionLength(vParticle);
random::UniformRealDistribution<InverseGrammageType> uniDist(actual_inv_length);
random::UniformRealDistribution<InverseGrammageType> uniDist(
current_inv_length);
const auto sample_process = uniDist(fRNG);
InverseGrammageType inv_lambda_count = 0. * meter * meter / gram;
fProcessSequence.SelectInteraction(vParticle, projectile, step, sample_process,
fProcessSequence.SelectInteraction(vParticle, projectile, sample_process,
inv_lambda_count);
} else if (min_distance == distance_decay) {
std::cout << "decay" << std::endl;
......
......@@ -69,8 +69,8 @@ public:
ProcessSplit(GrammageType const X0)
: fX0(X0) {}
template <typename Particle, typename Track>
corsika::units::si::GrammageType GetInteractionLength(Particle&, Track&) const {
template <typename Particle>
corsika::units::si::GrammageType GetInteractionLength(Particle const&) const {
return fX0;
}
......
......@@ -38,13 +38,12 @@ namespace corsika::process {
template <typename Particle>
EProcessReturn DoInteraction(Particle&);
template <typename Particle, typename Track>
corsika::units::si::GrammageType GetInteractionLength(Particle& p, Track& t);
template <typename TParticle>
corsika::units::si::GrammageType GetInteractionLength(TParticle& p);
template <typename Particle, typename Track>
corsika::units::si::InverseGrammageType GetInverseInteractionLength(Particle& p,
Track& t) {
return 1. / GetRef().GetInteractionLength(p, t);
template <typename TParticle>
corsika::units::si::InverseGrammageType GetInverseInteractionLength(TParticle& p) {
return 1. / GetRef().GetInteractionLength(p);
}
};
......
......@@ -165,51 +165,49 @@ namespace corsika::process {
return max_length;
}
template <typename TParticle, typename TTrack>
corsika::units::si::GrammageType GetTotalInteractionLength(TParticle& vP,
TTrack& vT) {
return 1. / GetInverseInteractionLength(vP, vT);
template <typename TParticle>
corsika::units::si::GrammageType GetTotalInteractionLength(TParticle& vP) {
return 1. / GetInverseInteractionLength(vP);
}
template <typename TParticle, typename TTrack>
template <typename TParticle>
corsika::units::si::InverseGrammageType GetTotalInverseInteractionLength(
TParticle& vP, TTrack& vT) {
return GetInverseInteractionLength(vP, vT);
TParticle& vP) {
return GetInverseInteractionLength(vP);
}
template <typename TParticle, typename TTrack>
corsika::units::si::InverseGrammageType GetInverseInteractionLength(TParticle& vP,
TTrack& vT) {
template <typename TParticle>
corsika::units::si::InverseGrammageType GetInverseInteractionLength(TParticle& vP) {
using namespace corsika::units::si;
InverseGrammageType tot = 0 * meter * meter / gram;
if constexpr (std::is_base_of<InteractionProcess<T1type>, T1type>::value ||
is_process_sequence<T1>::value) {
tot += A.GetInverseInteractionLength(vP, vT);
if constexpr (std::is_base_of_v<InteractionProcess<T1type>, T1type> ||
is_process_sequence_v<T1>) {
tot += A.GetInverseInteractionLength(vP);
}
if constexpr (std::is_base_of<InteractionProcess<T2type>, T2type>::value ||
is_process_sequence<T2>::value) {
tot += B.GetInverseInteractionLength(vP, vT);
if constexpr (std::is_base_of_v<InteractionProcess<T2type>, T2type> ||
is_process_sequence_v<T2>) {
tot += B.GetInverseInteractionLength(vP);
}
return tot;
}
template <typename TParticle, typename TSecondaries, typename TTrack>
template <typename TParticle, typename TSecondaries>
EProcessReturn SelectInteraction(
TParticle& vP, TSecondaries& vS, TTrack& vT,
TParticle& vP, TSecondaries& vS,
[[maybe_unused]] corsika::units::si::InverseGrammageType lambda_select,
corsika::units::si::InverseGrammageType& lambda_inv_count) {
if constexpr (is_process_sequence<T1type>::value) {
// if A is a process sequence --> check inside
const EProcessReturn ret =
A.SelectInteraction(vP, vS, vT, lambda_select, lambda_inv_count);
A.SelectInteraction(vP, vS, lambda_select, lambda_inv_count);
// if A did succeed, stop routine
if (ret != EProcessReturn::eOk) { return ret; }
} else if constexpr (std::is_base_of<InteractionProcess<T1type>, T1type>::value) {
// if this is not a ContinuousProcess --> evaluate probability
lambda_inv_count += A.GetInverseInteractionLength(vP, vT);
lambda_inv_count += A.GetInverseInteractionLength(vP);
// check if we should execute THIS process and then EXIT
if (lambda_select < lambda_inv_count) {
A.DoInteraction(vS);
......@@ -220,12 +218,12 @@ namespace corsika::process {
if constexpr (is_process_sequence<T2>::value) {
// if A is a process sequence --> check inside
const EProcessReturn ret =
B.SelectInteraction(vP, vS, vT, lambda_select, lambda_inv_count);
B.SelectInteraction(vP, vS, lambda_select, lambda_inv_count);
// if A did succeed, stop routine
if (ret != EProcessReturn::eOk) { return ret; }
} else if constexpr (std::is_base_of<InteractionProcess<T2type>, T2type>::value) {
// if this is not a ContinuousProcess --> evaluate probability
lambda_inv_count += B.GetInverseInteractionLength(vP, vT);
lambda_inv_count += B.GetInverseInteractionLength(vP);
// check if we should execute THIS process and then EXIT
if (lambda_select < lambda_inv_count) {
B.DoInteraction(vS);
......
......@@ -100,8 +100,8 @@ public:
cout << "Process2::DoInteraction" << endl;
return EProcessReturn::eOk;
}
template <typename Particle, typename Track>
GrammageType GetInteractionLength(Particle&, Track&) const {
template <typename Particle>
GrammageType GetInteractionLength(Particle&) const {
cout << "Process2::GetInteractionLength" << endl;
return 3_g / (1_cm * 1_cm);
}
......@@ -123,8 +123,8 @@ public:
cout << "Process3::DoInteraction" << endl;
return EProcessReturn::eOk;
}
template <typename Particle, typename Track>
GrammageType GetInteractionLength(Particle&, Track&) const {
template <typename Particle>
GrammageType GetInteractionLength(Particle&) const {
cout << "Process3::GetInteractionLength" << endl;
return 1_g / (1_cm * 1_cm);
}
......@@ -220,12 +220,11 @@ TEST_CASE("Process Sequence", "[Process Sequence]") {
Process3 m3(2);
DummyData particle;
DummyTrajectory track;
auto sequence2 = cp1 << m2 << m3;
GrammageType const tot = sequence2.GetTotalInteractionLength(particle, track);
GrammageType const tot = sequence2.GetTotalInteractionLength(particle);
InverseGrammageType const tot_inv =
sequence2.GetTotalInverseInteractionLength(particle, track);
sequence2.GetTotalInverseInteractionLength(particle);
cout << "lambda_tot=" << tot << "; lambda_tot_inv=" << tot_inv << endl;
}
......
......@@ -18,14 +18,12 @@
#include <corsika/utl/COMBoost.h>
#include <corsika/setup/SetupStack.h>
#include <corsika/setup/SetupTrajectory.h>
#include <iomanip>
#include <iostream>
using namespace corsika::setup;
using Particle = corsika::setup::Stack::ParticleType;
using Track = corsika::setup::Trajectory;
using SetupParticle = corsika::setup::Stack::ParticleType;
namespace corsika::process::HadronicElasticModel {
......@@ -38,7 +36,7 @@ namespace corsika::process::HadronicElasticModel {
template <>
units::si::GrammageType HadronicElasticInteraction::GetInteractionLength(
Particle const& p, Track&) {
SetupParticle const& p) {
using namespace units::si;
if (p.GetPID() == particles::Code::Proton) {
auto const* currentNode = p.GetNode();
......@@ -79,7 +77,7 @@ namespace corsika::process::HadronicElasticModel {
}
template <>
process::EProcessReturn HadronicElasticInteraction::DoInteraction(Particle& p) {
process::EProcessReturn HadronicElasticInteraction::DoInteraction(SetupParticle& p) {
if (p.GetPID() != particles::Code::Proton) { return process::EProcessReturn::eOk; }
using namespace units::si;
......
......@@ -56,8 +56,8 @@ namespace corsika::process::HadronicElasticModel {
units::si::CrossSectionType y = 0.05608 * units::si::barn);
void Init();
template <typename Particle, typename Track>
corsika::units::si::GrammageType GetInteractionLength(Particle const& p, Track&);
template <typename Particle>
corsika::units::si::GrammageType GetInteractionLength(Particle const& p);
template <typename Particle>
corsika::process::EProcessReturn DoInteraction(Particle&);
......
......@@ -15,7 +15,6 @@
#include <corsika/environment/NuclearComposition.h>
#include <corsika/geometry/FourVector.h>
#include <corsika/setup/SetupStack.h>
#include <corsika/setup/SetupTrajectory.h>
#include <corsika/utl/COMBoost.h>
#include <tuple>
......@@ -28,7 +27,6 @@ using namespace corsika;
using namespace corsika::setup;
using Projectile = corsika::setup::StackView::ParticleType;
using Particle = corsika::setup::Stack::ParticleType;
using Track = Trajectory;
namespace corsika::process::pythia {
......@@ -168,7 +166,7 @@ namespace corsika::process::pythia {
}
template <>
units::si::GrammageType Interaction::GetInteractionLength(Particle& p, Track&) {
units::si::GrammageType Interaction::GetInteractionLength(Particle& p) {
using namespace units;
using namespace units::si;
......
......@@ -51,16 +51,16 @@ namespace corsika::process::pythia {
const corsika::particles::Code TargetId,
const corsika::units::si::HEPEnergyType CoMenergy);
template <typename TParticle, typename TTrack>
corsika::units::si::GrammageType GetInteractionLength(TParticle&, TTrack&);
template <typename TParticle>
corsika::units::si::GrammageType GetInteractionLength(TParticle&);
/**
In this function PYTHIA is called to produce one event. The
event is copied (and boosted) into the shower lab frame.
*/
template <typename TProjctile>
corsika::process::EProcessReturn DoInteraction(TProjctile&);
template <typename TProjectile>
corsika::process::EProcessReturn DoInteraction(TProjectile&);
private:
corsika::random::RNG& fRNG =
......
......@@ -118,14 +118,6 @@ TEST_CASE("pythia process") {
auto const* nodePtr = theMedium.get(); // save the medium for later use before moving it
universe.AddChild(std::move(theMedium));
geometry::Point const origin(cs, {0_m, 0_m, 0_m});
geometry::Vector<units::si::SpeedType::dimension_type> v(cs, 0_m / second, 0_m / second,
1_m / second);
geometry::Line line(origin, v);
geometry::Trajectory<geometry::Line> track(line, 10_s);
random::RNGManager::GetInstance().RegisterRandomStream("s_rndm");
SECTION("pythia decay") {
setup::Stack stack;
......@@ -173,7 +165,6 @@ TEST_CASE("pythia process") {
process::pythia::Interaction model;
model.Init();
[[maybe_unused]] const process::EProcessReturn ret = model.DoInteraction(projectile);
[[maybe_unused]] const GrammageType length =
model.GetInteractionLength(particle, track);
[[maybe_unused]] const GrammageType length = model.GetInteractionLength(particle);
}
}
......@@ -118,8 +118,8 @@ namespace corsika::process::sibyll {
}
template <>
units::si::GrammageType Interaction::GetInteractionLength(SetupParticle& vP,
Track&) const {
units::si::GrammageType Interaction::GetInteractionLength(
SetupParticle const& vP) const {
using namespace units;
using namespace units::si;
......
......@@ -52,8 +52,8 @@ namespace corsika::process::sibyll {
GetCrossSection(const corsika::particles::Code, const corsika::particles::Code,
const corsika::units::si::HEPEnergyType) const;
template <typename TParticle, typename TTrack>
corsika::units::si::GrammageType GetInteractionLength(TParticle&, TTrack&) const;
template <typename TParticle>
corsika::units::si::GrammageType GetInteractionLength(TParticle const&) const;
/**
In this function SIBYLL is called to produce one event. The
......
......@@ -217,7 +217,7 @@ namespace corsika::process::sibyll {
template <>
template <>
units::si::GrammageType NuclearInteraction<SetupEnvironment>::GetInteractionLength(
Particle& vP, Track&) {
Particle& vP) {
using namespace units;
using namespace units::si;
......
......@@ -53,11 +53,11 @@ namespace corsika::process::sibyll {
std::tuple<corsika::units::si::CrossSectionType, corsika::units::si::CrossSectionType>
GetCrossSection(Particle& p, const corsika::particles::Code TargetId);
template <typename Particle, typename Track>
corsika::units::si::GrammageType GetInteractionLength(Particle&, Track&);
template <typename Particle>
corsika::units::si::GrammageType GetInteractionLength(Particle&);
template <typename Projectle>
corsika::process::EProcessReturn DoInteraction(Projectle&);
template <typename Projectile>
corsika::process::EProcessReturn DoInteraction(Projectile&);
private:
TEnvironment const& fEnvironment;
......
......@@ -98,12 +98,6 @@ TEST_CASE("SibyllInterface", "[processes]") {
const geometry::CoordinateSystem& cs = env.GetCoordinateSystem();
geometry::Point const origin(cs, {0_m, 0_m, 0_m});
geometry::Vector<units::si::SpeedType::dimension_type> v(cs, 0_m / second, 0_m / second,
1_m / second);
geometry::Line line(origin, v);
geometry::Trajectory<geometry::Line> track(line, 10_s);
random::RNGManager::GetInstance().RegisterRandomStream("s_rndm");
SECTION("InteractionInterface") {
......@@ -126,8 +120,7 @@ TEST_CASE("SibyllInterface", "[processes]") {
model.Init();
[[maybe_unused]] const process::EProcessReturn ret = model.DoInteraction(projectile);
[[maybe_unused]] const GrammageType length =
model.GetInteractionLength(particle, track);
[[maybe_unused]] const GrammageType length = model.GetInteractionLength(particle);
}
SECTION("NuclearInteractionInterface") {
......@@ -153,8 +146,7 @@ TEST_CASE("SibyllInterface", "[processes]") {
model.Init();
[[maybe_unused]] const process::EProcessReturn ret = model.DoInteraction(projectile);
[[maybe_unused]] const GrammageType length =
model.GetInteractionLength(particle, track);
[[maybe_unused]] const GrammageType length = model.GetInteractionLength(particle);
}
SECTION("DecayInterface") {
......
......@@ -2,9 +2,6 @@
#include <corsika/geometry/Vector.h>
#include <corsika/particles/ParticleProperties.h>
#include <corsika/process/urqmd/UrQMD.h>
#include <corsika/random/RNGManager.h>
#include <corsika/setup/SetupStack.h>
#include <corsika/setup/SetupTrajectory.h>
#include <corsika/units/PhysicalUnits.h>
#include <algorithm>
......@@ -20,7 +17,6 @@ UrQMD::UrQMD() { iniurqmd_(); }
using SetupStack = corsika::setup::Stack;
using SetupParticle = corsika::setup::Stack::StackIterator;
using SetupProjectile = corsika::setup::StackView::StackIterator;
using SetupTrack = corsika::setup::Trajectory;
template <typename TParticle> // need template here, as this is called both with
// SetupParticle as well as SetupProjectile
......@@ -82,7 +78,7 @@ bool UrQMD::CanInteract(particles::Code vCode) const {
vCode) != std::cend(validProjectileCodes);
}
GrammageType UrQMD::GetInteractionLength(SetupParticle& vParticle, SetupTrack&) const {
GrammageType UrQMD::GetInteractionLength(SetupParticle& vParticle) const {
if (!CanInteract(vParticle.GetPID())) {
// we could do the canInteract check in GetCrossSection, too but if
// we do it here we have the advantage of avoiding the loop
......
......@@ -5,7 +5,6 @@
#include <corsika/process/InteractionProcess.h>
#include <corsika/random/RNGManager.h>
#include <corsika/setup/SetupStack.h>
#include <corsika/setup/SetupTrajectory.h>
#include <corsika/units/PhysicalUnits.h>
#include <array>
......@@ -17,7 +16,7 @@ namespace corsika::process::UrQMD {
UrQMD();
void Init() {}
corsika::units::si::GrammageType GetInteractionLength(
corsika::setup::Stack::StackIterator&, corsika::setup::Trajectory&) const;
corsika::setup::Stack::StackIterator&) const;
template <typename TParticle>
corsika::units::si::CrossSectionType GetCrossSection(TParticle const&,
......
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment