diff --git a/Documentation/Examples/vertical_EAS.cc b/Documentation/Examples/vertical_EAS.cc
index 4e8ca5cfd47d75d706a82c3bd320feffa57e1a3b..a19c8bd89b0d7cf32884ed5589d08e11b299e720 100644
--- a/Documentation/Examples/vertical_EAS.cc
+++ b/Documentation/Examples/vertical_EAS.cc
@@ -226,9 +226,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/Framework/Cascade/Cascade.h b/Framework/Cascade/Cascade.h
index fc59e6e8b480a1bb1535179db865df168c0f39f2..de033d190b037c5e30b713be5d9d6b39eedfeef5 100644
--- a/Framework/Cascade/Cascade.h
+++ b/Framework/Cascade/Cascade.h
@@ -96,7 +96,10 @@ namespace corsika::cascade {
, fTracking(tr)
, fProcessSequence(pl)
, fStack(stack)
- , energy_cut_(0 * corsika::units::si::electronvolt) {}
+ , energy_cut_(0 * corsika::units::si::electronvolt) {
+
+ std::cout << c8_ascii_ << std::endl;
+ }
corsika::units::si::HEPEnergyType GetEnergyCut() const { return energy_cut_; }
@@ -121,9 +124,9 @@ namespace corsika::cascade {
while (!fStack.IsEmpty()) {
while (!fStack.IsEmpty()) {
- count_++;
+ count_++;
auto pNext = fStack.GetNextParticle();
- std::cout << "========= next: count=" << count_ << ", pid=" << pNext.GetPID()
+ std::cout << "========= next: count=" << count_ << ", pid=" << pNext.GetPID()
<< ", stack entries=" << fStack.getEntries()
<< ", stack deleted=" << fStack.getDeleted() << std::endl;
Step(pNext);
@@ -145,8 +148,7 @@ namespace corsika::cascade {
std::cout << "forced interaction!" << std::endl;
auto vParticle = fStack.GetNextParticle();
TStackView secondaries(vParticle);
- auto projectile = secondaries.GetProjectile();
- interaction(vParticle, projectile);
+ interaction(vParticle, secondaries);
fProcessSequence.DoSecondaries(secondaries);
vParticle.Delete(); // todo: this should be reviewed, see below
}
@@ -262,7 +264,7 @@ namespace corsika::cascade {
[[maybe_unused]] auto projectile = secondaries.GetProjectile();
if (min_distance == distance_interact) {
- interaction(vParticle, projectile);
+ interaction(vParticle, secondaries);
} else {
assert(min_distance == distance_decay);
decay(vParticle, secondaries);
@@ -313,8 +315,7 @@ namespace corsika::cascade {
}
}
- auto decay(Particle& particle,
- TStackView& view) {
+ auto decay(Particle& particle, TStackView& view) {
std::cout << "decay" << std::endl;
units::si::InverseTimeType const actual_decay_time =
fProcessSequence.GetTotalInverseLifetime(particle);
@@ -323,12 +324,11 @@ namespace corsika::cascade {
actual_decay_time);
const auto sample_process = uniDist(fRNG);
units::si::InverseTimeType inv_decay_count = units::si::InverseTimeType::zero();
- return fProcessSequence.SelectDecay(particle, projectile, sample_process,
+ return fProcessSequence.SelectDecay(particle, view, sample_process,
inv_decay_count);
}
- auto interaction(Particle& particle,
- decltype(std::declval<TStackView>().GetProjectile()) projectile) {
+ auto interaction(Particle& particle, TStackView& view) {
std::cout << "collide" << std::endl;
units::si::InverseGrammageType const current_inv_length =
@@ -338,13 +338,13 @@ namespace corsika::cascade {
current_inv_length);
const auto sample_process = uniDist(fRNG);
auto inv_lambda_count = units::si::InverseGrammageType::zero();
- return fProcessSequence.SelectInteraction(particle, projectile, sample_process,
+ return fProcessSequence.SelectInteraction(particle, view, sample_process,
inv_lambda_count);
}
// but this here temporarily. Should go into dedicated file later:
- const char* c8_ascii_ =
-R"V0G0N(
+ const char* c8_ascii_ =
+ R"V0G0N(
,ad8888ba, ,ad8888ba, 88888888ba ad88888ba 88 88 a8P db ad88888ba
d8"' `"8b d8"' `"8b 88 "8b d8" "8b 88 88 ,88' d88b d8" "8b
d8' d8' `8b 88 ,8P Y8, 88 88 ,88" d8'`8b Y8a a8P
@@ -354,7 +354,6 @@ Y8, Y8, ,8P 88 `8b `8b 88 88P Y8b d8"
Y8a. .a8P Y8a. .a8P 88 `8b Y8a a8P 88 88 "88, d8' `8b Y8a a8P
`"Y8888Y"' `"Y8888Y"' 88 `8b "Y88888P" 88 88 Y8b d8' `8b "Y88888P"
)V0G0N";
-
- };
+ };
} // namespace corsika::cascade
diff --git a/Framework/Cascade/testCascade.cc b/Framework/Cascade/testCascade.cc
index 57f715ef5b0e6a71f541390e2fc96d49ebca9826..9db2c73ccba340df465f7ac75c518d163fef6cd2 100644
--- a/Framework/Cascade/testCascade.cc
+++ b/Framework/Cascade/testCascade.cc
@@ -69,18 +69,19 @@ public:
return fX0;
}
- template <typename TProjectile>
- corsika::process::EProcessReturn DoInteraction(TProjectile& vP) {
+ template <typename TSecondaryView>
+ corsika::process::EProcessReturn DoInteraction(TSecondaryView& view) {
fCalls++;
- const HEPEnergyType E = vP.GetEnergy();
- vP.AddSecondary(
+ auto const projectile = view.GetProjectile();
+ const HEPEnergyType E = projectile.GetEnergy();
+ view.AddSecondary(
std::tuple<particles::Code, units::si::HEPEnergyType,
corsika::stack::MomentumVector, geometry::Point, units::si::TimeType>{
- vP.GetPID(), E / 2, vP.GetMomentum(), vP.GetPosition(), vP.GetTime()});
- vP.AddSecondary(
+ projectile.GetPID(), E / 2, projectile.GetMomentum(), projectile.GetPosition(), projectile.GetTime()});
+ view.AddSecondary(
std::tuple<particles::Code, units::si::HEPEnergyType,
corsika::stack::MomentumVector, geometry::Point, units::si::TimeType>{
- vP.GetPID(), E / 2, vP.GetMomentum(), vP.GetPosition(), vP.GetTime()});
+ projectile.GetPID(), E / 2, projectile.GetMomentum(), projectile.GetPosition(), projectile.GetTime()});
return EProcessReturn::eInteracted;
}
diff --git a/Processes/HadronicElasticModel/HadronicElasticModel.cc b/Processes/HadronicElasticModel/HadronicElasticModel.cc
index d95c69d657cb36709ad25ae6febb3e293709e9dc..c8f3eed3dd74a39d1de39e40be9112c5dab5f54e 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,178 +18,177 @@
#include <iomanip>
#include <iostream>
-using namespace corsika::setup;
using SetupParticle = corsika::setup::Stack::ParticleType;
+using SetupView = corsika::setup::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();
-
- auto const projectileMomentum = p.GetMomentum();
- auto const projectileMomentumSquaredNorm = projectileMomentum.squaredNorm();
- auto const projectileEnergy = p.GetEnergy();
-
- auto const avgCrossSection = [&]() {
- CrossSectionType avgCrossSection = 0_b;
-
- 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();
+using namespace corsika::process::HadronicElasticModel;
+using namespace corsika;
- GrammageType const interactionLength =
- avgTargetMassNumber * units::constants::u / avgCrossSection;
+HadronicElasticInteraction::HadronicElasticInteraction(units::si::CrossSectionType x,
+ units::si::CrossSectionType y)
+ : fX(x)
+ , fY(y) {}
- return interactionLength;
- } else {
- return std::numeric_limits<double>::infinity() * 1_g / (1_cm * 1_cm);
- }
- }
-
- template <>
- process::EProcessReturn HadronicElasticInteraction::DoInteraction(SetupParticle& p) {
- if (p.GetPID() != particles::Code::Proton) { return process::EProcessReturn::eOk; }
-
- using namespace units::si;
- using namespace units::constants;
+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();
- 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/HadronicElasticModel/HadronicElasticModel.h b/Processes/HadronicElasticModel/HadronicElasticModel.h
index b719a30819c38e0057821e0af19208c118cbff4e..3273bfec15aad2a47128a169e1031790f0dda3c6 100644
--- a/Processes/HadronicElasticModel/HadronicElasticModel.h
+++ b/Processes/HadronicElasticModel/HadronicElasticModel.h
@@ -54,8 +54,8 @@ namespace corsika::process::HadronicElasticModel {
template <typename Particle>
corsika::units::si::GrammageType GetInteractionLength(Particle const& p);
- template <typename Particle>
- corsika::process::EProcessReturn DoInteraction(Particle&);
+ template <typename TStackView>
+ corsika::process::EProcessReturn DoInteraction(TStackView&);
};
} // namespace corsika::process::HadronicElasticModel
diff --git a/Processes/InteractionCounter/InteractionCounter.h b/Processes/InteractionCounter/InteractionCounter.h
index 2ca7db507620b2143991e0d6830225a0bcb8373b..f91aa51013db377eed82e2c3f3a48ea83915a916 100644
--- a/Processes/InteractionCounter/InteractionCounter.h
+++ b/Processes/InteractionCounter/InteractionCounter.h
@@ -30,8 +30,9 @@ namespace corsika::process::interaction_counter {
InteractionCounter(TCountedProcess& process)
: process_(process) {}
- template <typename TProjectile>
- auto DoInteraction(TProjectile& projectile) {
+ template <typename TSecondaryView>
+ auto DoInteraction(TSecondaryView& view) {
+ auto const projectile = view.GetProjectile();
auto const massNumber = projectile.GetNode()
->GetModelProperties()
.GetNuclearComposition()
@@ -46,7 +47,7 @@ namespace corsika::process::interaction_counter {
} else {
histogram_.fill(projectile_id, projectile.GetEnergy(), massTarget);
}
- return process_.DoInteraction(projectile);
+ return process_.DoInteraction(view);
}
template <typename TParticle>
diff --git a/Processes/InteractionCounter/testInteractionCounter.cc b/Processes/InteractionCounter/testInteractionCounter.cc
index 13bf34c1714c0554363b87296cf7f8584fa90127..c63579e8819920a298dd39ddc90ffa8973f31454 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/Interaction.cc b/Processes/Pythia/Interaction.cc
index f6b2a01403ce6bc3462c6d476fb94f1299c6bbf3..fe641a62629a5bfc2b75f4ce9d0ada6c984322ce 100644
--- a/Processes/Pythia/Interaction.cc
+++ b/Processes/Pythia/Interaction.cc
@@ -22,7 +22,7 @@ using std::tuple;
using namespace corsika;
using namespace corsika::setup;
-using Projectile = corsika::setup::StackView::ParticleType;
+using SecondaryView = corsika::setup::StackView;
using Particle = corsika::setup::Stack::ParticleType;
namespace corsika::process::pythia {
@@ -241,14 +241,16 @@ namespace corsika::process::pythia {
*/
template <>
- process::EProcessReturn Interaction::DoInteraction(Projectile& vP) {
+ process::EProcessReturn Interaction::DoInteraction(SecondaryView& view) {
using namespace units;
using namespace utl;
using namespace units::si;
using namespace geometry;
- const auto corsikaBeamId = vP.GetPID();
+ auto const projectile = view.GetProjectile();
+
+ const auto corsikaBeamId = projectile.GetPID();
cout << "Pythia::Interaction: "
<< "DoInteraction: " << corsikaBeamId << " interaction? "
<< process::pythia::Interaction::CanInteract(corsikaBeamId) << endl;
@@ -264,8 +266,8 @@ namespace corsika::process::pythia {
RootCoordinateSystem::GetInstance().GetRootCoordinateSystem();
// position and time of interaction, not used in Sibyll
- Point pOrig = vP.GetPosition();
- TimeType tOrig = vP.GetTime();
+ Point pOrig = projectile.GetPosition();
+ TimeType tOrig = projectile.GetTime();
// define target
// FOR NOW: target is always at rest
@@ -274,8 +276,8 @@ namespace corsika::process::pythia {
const FourVector PtargLab(eTargetLab, pTargetLab);
// define projectile
- HEPEnergyType const eProjectileLab = vP.GetEnergy();
- auto const pProjectileLab = vP.GetMomentum();
+ HEPEnergyType const eProjectileLab = projectile.GetEnergy();
+ auto const pProjectileLab = projectile.GetMomentum();
cout << "Interaction: ebeam lab: " << eProjectileLab / 1_GeV << endl
<< "Interaction: pbeam lab: " << pProjectileLab.GetComponents() / 1_GeV
@@ -310,12 +312,12 @@ namespace corsika::process::pythia {
cout << "Interaction: time: " << tOrig << endl;
HEPEnergyType Etot = eProjectileLab + eTargetLab;
- MomentumVector Ptot = vP.GetMomentum();
+ MomentumVector Ptot = projectile.GetMomentum();
// invariant mass, i.e. cm. energy
HEPEnergyType Ecm = sqrt(Etot * Etot - Ptot.squaredNorm());
// sample target mass number
- const auto* currentNode = vP.GetNode();
+ const auto* currentNode = projectile.GetNode();
const auto& mediumComposition =
currentNode->GetModelProperties().GetNuclearComposition();
// get cross sections for target materials
@@ -381,7 +383,7 @@ namespace corsika::process::pythia {
HEPEnergyType const pyEn = p8p.e() * 1_GeV;
// add to corsika stack
- auto pnew = vP.AddSecondary(
+ auto pnew = view.AddSecondary(
tuple<particles::Code, units::si::HEPEnergyType, stack::MomentumVector,
geometry::Point, units::si::TimeType>{pyId, pyEn, pyPlab, pOrig,
tOrig});
diff --git a/Processes/Pythia/Interaction.h b/Processes/Pythia/Interaction.h
index d6e3862bf72c3b1728c75931ee7bffb509adeed2..4b06aeb6d58f6bf034b49b054376e540fa279628 100644
--- a/Processes/Pythia/Interaction.h
+++ b/Processes/Pythia/Interaction.h
@@ -54,8 +54,8 @@ namespace corsika::process::pythia {
event is copied (and boosted) into the shower lab frame.
*/
- template <typename TProjectile>
- corsika::process::EProcessReturn DoInteraction(TProjectile&);
+ template <typename TSecondaryView>
+ corsika::process::EProcessReturn DoInteraction(TSecondaryView&);
private:
corsika::random::RNG& fRNG =
diff --git a/Processes/Pythia/testPythia8.cc b/Processes/Pythia/testPythia8.cc
index 680b8cadeb838a4f5e500d273a50b5aef185732c..4b1da6d837551d1bba4aded4d0c43cd8f315f1ca 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/Interaction.cc b/Processes/QGSJetII/Interaction.cc
index 4f2f1ca34ecbc6bbcc175040e9e99387f942e1ff..c16233af582c5908261f88e6de89af58748b0e09 100644
--- a/Processes/QGSJetII/Interaction.cc
+++ b/Processes/QGSJetII/Interaction.cc
@@ -33,7 +33,7 @@ using std::tuple;
using namespace corsika;
using namespace corsika::setup;
using SetupParticle = setup::Stack::StackIterator;
-using SetupProjectile = setup::StackView::StackIterator;
+using SetupView = setup::StackView;
using Track = Trajectory;
namespace corsika::process::qgsjetII {
@@ -169,14 +169,16 @@ namespace corsika::process::qgsjetII {
*/
template <>
- process::EProcessReturn Interaction::DoInteraction(SetupProjectile& vP) {
+ process::EProcessReturn Interaction::DoInteraction(SetupView& view) {
using namespace units;
using namespace utl;
using namespace units::si;
using namespace geometry;
- const auto corsikaBeamId = vP.GetPID();
+ auto const projectile = view.GetProjectile();
+
+ const auto corsikaBeamId = projectile.GetPID();
cout << "ProcessQgsjetII: "
<< "DoInteraction: " << corsikaBeamId << " interaction? "
<< process::qgsjetII::CanInteract(corsikaBeamId) << endl;
@@ -187,8 +189,8 @@ namespace corsika::process::qgsjetII {
RootCoordinateSystem::GetInstance().GetRootCoordinateSystem();
// position and time of interaction, not used in QgsjetII
- Point pOrig = vP.GetPosition();
- TimeType tOrig = vP.GetTime();
+ Point pOrig = projectile.GetPosition();
+ TimeType tOrig = projectile.GetTime();
// define target
// for QgsjetII is always a single nucleon
@@ -198,11 +200,11 @@ namespace corsika::process::qgsjetII {
const FourVector PtargLab(targetEnergyLab, targetMomentumLab);
// define projectile
- HEPEnergyType const projectileEnergyLab = vP.GetEnergy();
- auto const projectileMomentumLab = vP.GetMomentum();
+ HEPEnergyType const projectileEnergyLab = projectile.GetEnergy();
+ auto const projectileMomentumLab = projectile.GetMomentum();
int beamA = 1;
- if (particles::IsNucleus(corsikaBeamId)) beamA = vP.GetNuclearA();
+ if (particles::IsNucleus(corsikaBeamId)) beamA = projectile.GetNuclearA();
const HEPEnergyType projectileEnergyLabPerNucleon = projectileEnergyLab / beamA;
@@ -217,7 +219,7 @@ namespace corsika::process::qgsjetII {
cout << "Interaction: time: " << tOrig << endl;
// sample target mass number
- auto const* currentNode = vP.GetNode();
+ auto const* currentNode = projectile.GetNode();
auto const& mediumComposition =
currentNode->GetModelProperties().GetNuclearComposition();
// get cross sections for target materials
@@ -257,7 +259,7 @@ namespace corsika::process::qgsjetII {
QgsjetIIHadronType qgsjet_hadron_type =
process::qgsjetII::GetQgsjetIIHadronType(corsikaBeamId);
if (qgsjet_hadron_type == QgsjetIIHadronType::NucleusType) {
- projectileMassNumber = vP.GetNuclearA();
+ projectileMassNumber = projectile.GetNuclearA();
if (projectileMassNumber > maxMassNumber_)
throw std::runtime_error("QgsjetII projectile mass outside range.");
std::array<QgsjetIIHadronType, 2> constexpr nucleons = {
@@ -325,7 +327,7 @@ namespace corsika::process::qgsjetII {
momentum.rebase(originalCS); // transform back into standard lab frame
std::cout << "secondary fragment> id=" << idFragm
<< " p=" << momentum.GetComponents() << std::endl;
- auto pnew = vP.AddSecondary(
+ auto pnew = view.AddSecondary(
tuple<particles::Code, units::si::HEPEnergyType, stack::MomentumVector,
geometry::Point, units::si::TimeType>{idFragm, energy, momentum,
pOrig, tOrig});
@@ -361,7 +363,7 @@ namespace corsika::process::qgsjetII {
std::cout << "secondary fragment> id=" << idFragm
<< " p=" << momentum.GetComponents() << " A=" << A << " Z=" << Z
<< std::endl;
- auto pnew = vP.AddSecondary(
+ auto pnew = view.AddSecondary(
tuple<particles::Code, units::si::HEPEnergyType, stack::MomentumVector,
geometry::Point, units::si::TimeType, unsigned short, unsigned short>{
idFragm, energy, momentum, pOrig, tOrig, A, Z});
@@ -381,7 +383,7 @@ namespace corsika::process::qgsjetII {
std::cout << "secondary fragment> id="
<< process::qgsjetII::ConvertFromQgsjetII(psec.GetPID())
<< " p=" << momentum.GetComponents() << std::endl;
- auto pnew = vP.AddSecondary(
+ auto pnew = view.AddSecondary(
tuple<particles::Code, units::si::HEPEnergyType, stack::MomentumVector,
geometry::Point, units::si::TimeType>{
process::qgsjetII::ConvertFromQgsjetII(psec.GetPID()), energy, momentum,
diff --git a/Processes/QGSJetII/Interaction.h b/Processes/QGSJetII/Interaction.h
index fc44b8f884e487ce1ed904c47d7af3be2b51a177..bef66e3ce3ac61ea8975aaf89ed4ae3645b4d811 100644
--- a/Processes/QGSJetII/Interaction.h
+++ b/Processes/QGSJetII/Interaction.h
@@ -50,8 +50,8 @@ namespace corsika::process::qgsjetII {
event is copied (and boosted) into the shower lab frame.
*/
- template <typename TProjectile>
- corsika::process::EProcessReturn DoInteraction(TProjectile&);
+ template <typename TSecondaryView>
+ corsika::process::EProcessReturn DoInteraction(TSecondaryView&);
private:
corsika::random::RNG& rng_ =
diff --git a/Processes/QGSJetII/testQGSJetII.cc b/Processes/QGSJetII/testQGSJetII.cc
index 5acb717a51b77ef05d50b7eaf6b04e9bee674f2b..4f09f12841a185cf3295d6c592a9e617f65a4a13 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/Decay.cc b/Processes/Sibyll/Decay.cc
index 49ec39ca477cb23bff30654ab879fdb259500beb..805c0f76aea0ad8f49e39da08531c0b0b64b9b1c 100644
--- a/Processes/Sibyll/Decay.cc
+++ b/Processes/Sibyll/Decay.cc
@@ -182,7 +182,7 @@ namespace corsika::process::sibyll {
auto const projectile = view.GetProjectile();
- const particles::Code pCode = vP.GetPID();
+ const particles::Code pCode = projectile.GetPID();
// check if sibyll is configured to handle this decay!
if (!IsDecayHandled(pCode))
throw std::runtime_error("STOP! Sibyll not configured to execute this decay!");
diff --git a/Processes/Sibyll/Interaction.cc b/Processes/Sibyll/Interaction.cc
index d738dc81caa233416e74b1ca56d5e6e3d5662473..50da5e0f90729b0c56e7d0b45cd472aa37530842 100644
--- a/Processes/Sibyll/Interaction.cc
+++ b/Processes/Sibyll/Interaction.cc
@@ -27,7 +27,7 @@ using std::tuple;
using namespace corsika;
using namespace corsika::setup;
using SetupParticle = setup::Stack::StackIterator;
-using SetupProjectile = setup::StackView::StackIterator;
+using SetupView = setup::StackView;
using Track = Trajectory;
namespace corsika::process::sibyll {
@@ -160,13 +160,15 @@ namespace corsika::process::sibyll {
*/
template <>
- process::EProcessReturn Interaction::DoInteraction(SetupProjectile& vP) {
+ process::EProcessReturn Interaction::DoInteraction(SetupView& view) {
using namespace utl;
using namespace units;
using namespace units::si;
using namespace geometry;
- const auto corsikaBeamId = vP.GetPID();
+ auto const projectile = view.GetProjectile();
+
+ const auto corsikaBeamId = projectile.GetPID();
cout << "ProcessSibyll: "
<< "DoInteraction: " << corsikaBeamId << " interaction? "
<< process::sibyll::CanInteract(corsikaBeamId) << endl;
@@ -178,12 +180,12 @@ namespace corsika::process::sibyll {
if (process::sibyll::CanInteract(corsikaBeamId)) {
// position and time of interaction, not used in Sibyll
- Point const pOrig = vP.GetPosition();
- TimeType const tOrig = vP.GetTime();
+ Point const pOrig = projectile.GetPosition();
+ TimeType const tOrig = projectile.GetTime();
// define projectile
- HEPEnergyType const eProjectileLab = vP.GetEnergy();
- auto const pProjectileLab = vP.GetMomentum();
+ HEPEnergyType const eProjectileLab = projectile.GetEnergy();
+ auto const pProjectileLab = projectile.GetMomentum();
const CoordinateSystem& originalCS = pProjectileLab.GetCoordinateSystem();
// define target
@@ -227,12 +229,12 @@ namespace corsika::process::sibyll {
cout << "Interaction: time: " << tOrig << endl;
HEPEnergyType Etot = eProjectileLab + eTargetLab;
- MomentumVector Ptot = vP.GetMomentum();
+ MomentumVector Ptot = projectile.GetMomentum();
// invariant mass, i.e. cm. energy
HEPEnergyType Ecm = sqrt(Etot * Etot - Ptot.squaredNorm());
// sample target mass number
- auto const* currentNode = vP.GetNode();
+ auto const* currentNode = projectile.GetNode();
auto const& mediumComposition =
currentNode->GetModelProperties().GetNuclearComposition();
// get cross sections for target materials
@@ -316,7 +318,7 @@ namespace corsika::process::sibyll {
assert(p3lab.GetCoordinateSystem() == originalCS); // just to be sure!
// add to corsika stack
- auto pnew = vP.AddSecondary(
+ auto pnew = view.AddSecondary(
tuple<particles::Code, units::si::HEPEnergyType, stack::MomentumVector,
geometry::Point, units::si::TimeType>{
process::sibyll::ConvertFromSibyll(psib.GetPID()),
diff --git a/Processes/Sibyll/Interaction.h b/Processes/Sibyll/Interaction.h
index fcb3b94e5f342134b16289286e66ee977d157abf..cccc971bd36517f955f23da88d6b7a588fd4341d 100644
--- a/Processes/Sibyll/Interaction.h
+++ b/Processes/Sibyll/Interaction.h
@@ -52,8 +52,8 @@ namespace corsika::process::sibyll {
event is copied (and boosted) into the shower lab frame.
*/
- template <typename TProjectile>
- corsika::process::EProcessReturn DoInteraction(TProjectile&);
+ template <typename TSecondaryView>
+ corsika::process::EProcessReturn DoInteraction(TSecondaryView&);
private:
corsika::random::RNG& RNG_ =
diff --git a/Processes/Sibyll/NuclearInteraction.cc b/Processes/Sibyll/NuclearInteraction.cc
index 451a5c643478a54db7199a0ef919ddc464cacb5c..1a814e0dd00f0684e228450541496e6c4dcdf85a 100644
--- a/Processes/Sibyll/NuclearInteraction.cc
+++ b/Processes/Sibyll/NuclearInteraction.cc
@@ -28,8 +28,8 @@ using std::vector;
using namespace corsika;
using namespace corsika::setup;
-using Particle = Stack::ParticleType; // StackIterator; // ParticleType;
-using Projectile = StackView::StackIterator; // StackView::ParticleType;
+using Particle = corsika::setup::Stack::ParticleType; // StackIterator; // ParticleType;
+using View = corsika::setup::StackView; // StackView::ParticleType;
using Track = Trajectory;
namespace corsika::process::sibyll {
@@ -293,7 +293,7 @@ namespace corsika::process::sibyll {
template <>
template <>
process::EProcessReturn NuclearInteraction<SetupEnvironment>::DoInteraction(
- Projectile& vP) {
+ View& view) {
// this routine superimposes different nucleon-nucleon interactions
// in a nucleus-nucleus interaction, based the SIBYLL routine SIBNUC
@@ -303,9 +303,11 @@ namespace corsika::process::sibyll {
using namespace units::si;
using namespace geometry;
- const auto ProjId = vP.GetPID();
+ auto projectile = view.GetProjectile();
+
+ const auto ProjId = projectile.GetPID();
// TODO: calculate projectile mass in nuclearStackExtension
- // const auto ProjMass = vP.GetMass();
+ // const auto ProjMass = projectile.GetMass();
cout << "NuclearInteraction: DoInteraction: called with:" << ProjId << endl;
// check if target-style nucleus (enum)
@@ -314,9 +316,9 @@ namespace corsika::process::sibyll {
"NuclearInteraction: DoInteraction: Wrong nucleus type. Nuclear projectiles "
"should use NuclearStackExtension!");
- auto const ProjMass =
- vP.GetNuclearZ() * particles::Proton::GetMass() +
- (vP.GetNuclearA() - vP.GetNuclearZ()) * particles::Neutron::GetMass();
+ auto const ProjMass = projectile.GetNuclearZ() * particles::Proton::GetMass() +
+ (projectile.GetNuclearA() - projectile.GetNuclearZ()) *
+ particles::Neutron::GetMass();
cout << "NuclearInteraction: projectile mass: " << ProjMass / 1_GeV << endl;
count_++;
@@ -325,21 +327,21 @@ namespace corsika::process::sibyll {
RootCoordinateSystem::GetInstance().GetRootCoordinateSystem();
// position and time of interaction, not used in NUCLIB
- Point pOrig = vP.GetPosition();
- TimeType tOrig = vP.GetTime();
+ Point pOrig = projectile.GetPosition();
+ TimeType tOrig = projectile.GetTime();
cout << "Interaction: position of interaction: " << pOrig.GetCoordinates() << endl;
cout << "Interaction: time: " << tOrig << endl;
// projectile nucleon number
- const int kAProj = vP.GetNuclearA();
+ const int kAProj = projectile.GetNuclearA();
if (kAProj > GetMaxNucleusAProjectile())
throw std::runtime_error("Projectile nucleus too large for NUCLIB!");
// kinematics
// define projectile nucleus
- HEPEnergyType const eProjectileLab = vP.GetEnergy();
- auto const pProjectileLab = vP.GetMomentum();
+ HEPEnergyType const eProjectileLab = projectile.GetEnergy();
+ auto const pProjectileLab = projectile.GetMomentum();
const FourVector PprojLab(eProjectileLab, pProjectileLab);
cout << "NuclearInteraction: eProj lab: " << eProjectileLab / 1_GeV << endl
@@ -347,8 +349,8 @@ namespace corsika::process::sibyll {
<< endl;
// define projectile nucleon
- HEPEnergyType const eProjectileNucLab = vP.GetEnergy() / kAProj;
- auto const pProjectileNucLab = vP.GetMomentum() / kAProj;
+ HEPEnergyType const eProjectileNucLab = projectile.GetEnergy() / kAProj;
+ auto const pProjectileNucLab = projectile.GetMomentum() / kAProj;
const FourVector PprojNucLab(eProjectileNucLab, pProjectileNucLab);
cout << "NuclearInteraction: eProjNucleon lab: " << eProjectileNucLab / 1_GeV << endl
@@ -400,7 +402,7 @@ namespace corsika::process::sibyll {
//
// proton stand-in for nucleon
const auto beamId = particles::Proton::GetCode();
- auto const* const currentNode = vP.GetNode();
+ auto const* const currentNode = projectile.GetNode();
const auto& mediumComposition =
currentNode->GetModelProperties().GetNuclearComposition();
cout << "get nucleon-nucleus cross sections for target materials.." << endl;
@@ -525,18 +527,19 @@ namespace corsika::process::sibyll {
if (nuclA == 1)
// add nucleon
- vP.AddSecondary(
+ projectile.AddSecondary(
tuple<particles::Code, units::si::HEPEnergyType, stack::MomentumVector,
geometry::Point, units::si::TimeType>{
specCode, Plab.GetTimeLikeComponent(), Plab.GetSpaceLikeComponents(),
pOrig, tOrig});
else
// add nucleus
- vP.AddSecondary(tuple<particles::Code, units::si::HEPEnergyType,
- corsika::stack::MomentumVector, geometry::Point,
- units::si::TimeType, unsigned short, unsigned short>{
- specCode, Plab.GetTimeLikeComponent(), Plab.GetSpaceLikeComponents(), pOrig,
- tOrig, nuclA, nuclZ});
+ projectile.AddSecondary(
+ tuple<particles::Code, units::si::HEPEnergyType,
+ corsika::stack::MomentumVector, geometry::Point, units::si::TimeType,
+ unsigned short, unsigned short>{specCode, Plab.GetTimeLikeComponent(),
+ Plab.GetSpaceLikeComponents(), pOrig,
+ tOrig, nuclA, nuclZ});
}
// add elastic nucleons to corsika stack
@@ -553,7 +556,7 @@ namespace corsika::process::sibyll {
const double mass_ratio = particles::GetMass(elaNucCode) / ProjMass;
auto const Plab = PprojLab * mass_ratio;
- vP.AddSecondary(
+ projectile.AddSecondary(
tuple<particles::Code, units::si::HEPEnergyType, corsika::stack::MomentumVector,
geometry::Point, units::si::TimeType>{
elaNucCode, Plab.GetTimeLikeComponent(), Plab.GetSpaceLikeComponents(),
@@ -567,14 +570,14 @@ namespace corsika::process::sibyll {
auto pCode = particles::Proton::GetCode();
// temporarily add to stack, will be removed after interaction in DoInteraction
cout << "inelastic interaction no. " << j << endl;
- auto inelasticNucleon = vP.AddSecondary(
+ auto inelasticNucleon = projectile.AddSecondary(
tuple<particles::Code, units::si::HEPEnergyType, corsika::stack::MomentumVector,
geometry::Point, units::si::TimeType>{
pCode, PprojNucLab.GetTimeLikeComponent(),
PprojNucLab.GetSpaceLikeComponents(), pOrig, tOrig});
// create inelastic interaction
cout << "calling HadronicInteraction..." << endl;
- hadronicInteraction_.DoInteraction(inelasticNucleon);
+ //~ hadronicInteraction_.DoInteraction(inelasticNucleon);
}
cout << "NuclearInteraction: DoInteraction: done" << endl;
diff --git a/Processes/Sibyll/NuclearInteraction.h b/Processes/Sibyll/NuclearInteraction.h
index ba9c710bb09363da5ee9f82d15bfeb6f41f397ab..3f5c4aeb6df7aedb7586d1a9247190e0e55e9698 100644
--- a/Processes/Sibyll/NuclearInteraction.h
+++ b/Processes/Sibyll/NuclearInteraction.h
@@ -52,8 +52,8 @@ namespace corsika::process::sibyll {
template <typename Particle>
corsika::units::si::GrammageType GetInteractionLength(Particle const&);
- template <typename Projectile>
- corsika::process::EProcessReturn DoInteraction(Projectile&);
+ template <typename TSecondaryView>
+ corsika::process::EProcessReturn DoInteraction(TSecondaryView&);
private:
TEnvironment const& environment_;
diff --git a/Processes/Sibyll/testSibyll.cc b/Processes/Sibyll/testSibyll.cc
index bd7454856560ba6b6bccda22d820faaff311a0a9..e9940067054987ca8eb1108fb0df895669daa004 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/UrQMD.cc b/Processes/UrQMD/UrQMD.cc
index cc9e02f61d2f7c45144260e6871a17e63aa42c60..cec9beb3ee2608a4d0a0314ba857e8ff47cd926b 100644
--- a/Processes/UrQMD/UrQMD.cc
+++ b/Processes/UrQMD/UrQMD.cc
@@ -27,7 +27,7 @@ using namespace corsika::units::si;
using SetupStack = corsika::setup::Stack;
using SetupParticle = corsika::setup::Stack::StackIterator;
-using SetupProjectile = corsika::setup::StackView::StackIterator;
+using SetupView = corsika::setup::StackView;
UrQMD::UrQMD(std::string const& xs_file) {
readXSFile(xs_file);
@@ -240,9 +240,11 @@ GrammageType UrQMD::GetInteractionLength(SetupParticle const& particle) const {
weightedProdCrossSection;
}
-corsika::process::EProcessReturn UrQMD::DoInteraction(SetupProjectile& projectile) {
+corsika::process::EProcessReturn UrQMD::DoInteraction(SetupView& view) {
using namespace units::si;
+ auto const projectile = view.GetProjectile();
+
auto projectileCode = projectile.GetPID();
auto const projectileEnergyLab = projectile.GetEnergy();
auto const& projectileMomentumLab = projectile.GetMomentum();
@@ -343,7 +345,7 @@ corsika::process::EProcessReturn UrQMD::DoInteraction(SetupProjectile& projectil
momentum.rebase(originalCS); // transform back into standard lab frame
std::cout << i << " " << code << " " << momentum.GetComponents() << std::endl;
- projectile.AddSecondary(
+ view.AddSecondary(
std::tuple<particles::Code, HEPEnergyType, stack::MomentumVector, geometry::Point,
TimeType>{code, energy, momentum, projectilePosition, projectileTime});
}
diff --git a/Processes/UrQMD/UrQMD.h b/Processes/UrQMD/UrQMD.h
index 3461c8ca36cd082eb0afc4731cdb44177c883201..937b6d702f63443fe0e221060e86739010ffcae8 100644
--- a/Processes/UrQMD/UrQMD.h
+++ b/Processes/UrQMD/UrQMD.h
@@ -40,8 +40,7 @@ namespace corsika::process::UrQMD {
corsika::units::si::CrossSectionType GetTabulatedCrossSection(
particles::Code, particles::Code, corsika::units::si::HEPEnergyType) const;
- corsika::process::EProcessReturn DoInteraction(
- corsika::setup::StackView::StackIterator&);
+ corsika::process::EProcessReturn DoInteraction(corsika::setup::StackView&);
bool CanInteract(particles::Code) const;
diff --git a/Processes/UrQMD/testUrQMD.cc b/Processes/UrQMD/testUrQMD.cc
index 11f11adfd621894c95a50e8ea4cf234ba0ba66c3..f7d36c891dc07b942186213a794184598cc12506 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) +