From 7a6af43617a4a48a217b30500453dae32f40fad3 Mon Sep 17 00:00:00 2001
From: Maximilian Reininghaus <maximilian.reininghaus@kit.edu>
Date: Thu, 24 Sep 2020 17:01:09 +0200
Subject: [PATCH] use SecondaryView in DoInteraction() instead of projectile
---
Framework/Cascade/Cascade.h | 26 ++++-----
.../HadronicElasticModel.cc | 5 +-
.../HadronicElasticModel.h | 4 +-
.../InteractionCounter/InteractionCounter.h | 7 +--
Processes/Pythia/Interaction.cc | 22 ++++----
Processes/Pythia/Interaction.h | 4 +-
Processes/QGSJetII/Interaction.cc | 26 ++++-----
Processes/QGSJetII/Interaction.h | 4 +-
Processes/Sibyll/Decay.cc | 2 +-
Processes/Sibyll/Interaction.cc | 22 ++++----
Processes/Sibyll/Interaction.h | 4 +-
Processes/Sibyll/NuclearInteraction.cc | 53 ++++++++++---------
Processes/Sibyll/NuclearInteraction.h | 4 +-
Processes/UrQMD/UrQMD.cc | 8 +--
Processes/UrQMD/UrQMD.h | 3 +-
15 files changed, 101 insertions(+), 93 deletions(-)
diff --git a/Framework/Cascade/Cascade.h b/Framework/Cascade/Cascade.h
index 490b04086..d83738c5e 100644
--- a/Framework/Cascade/Cascade.h
+++ b/Framework/Cascade/Cascade.h
@@ -98,9 +98,6 @@ namespace corsika::cascade {
, count_(0) {
std::cout << c8_ascii_ << std::endl;
-
-
-
}
/**
@@ -124,9 +121,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);
@@ -264,7 +261,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);
@@ -302,8 +299,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);
@@ -312,12 +308,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 =
@@ -327,13 +322,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
@@ -343,7 +338,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/Processes/HadronicElasticModel/HadronicElasticModel.cc b/Processes/HadronicElasticModel/HadronicElasticModel.cc
index d95c69d65..ece6447e8 100644
--- a/Processes/HadronicElasticModel/HadronicElasticModel.cc
+++ b/Processes/HadronicElasticModel/HadronicElasticModel.cc
@@ -21,6 +21,7 @@
using namespace corsika::setup;
using SetupParticle = corsika::setup::Stack::ParticleType;
+using SetupView = corsika::setup::Stack::StackView;
namespace corsika::process::HadronicElasticModel {
@@ -72,12 +73,14 @@ namespace corsika::process::HadronicElasticModel {
}
template <>
- process::EProcessReturn HadronicElasticInteraction::DoInteraction(SetupParticle& p) {
+ 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();
diff --git a/Processes/HadronicElasticModel/HadronicElasticModel.h b/Processes/HadronicElasticModel/HadronicElasticModel.h
index b719a3081..3273bfec1 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 2ca7db507..f91aa5101 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/Pythia/Interaction.cc b/Processes/Pythia/Interaction.cc
index f6b2a0140..fe641a626 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 d6e3862bf..4b06aeb6d 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/QGSJetII/Interaction.cc b/Processes/QGSJetII/Interaction.cc
index 4f2f1ca34..40aac437d 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});
diff --git a/Processes/QGSJetII/Interaction.h b/Processes/QGSJetII/Interaction.h
index fc44b8f88..bef66e3ce 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/Sibyll/Decay.cc b/Processes/Sibyll/Decay.cc
index 49ec39ca4..805c0f76a 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 d738dc81c..50da5e0f9 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 fcb3b94e5..cccc971bd 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 451a5c643..1a814e0dd 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 ba9c710bb..3f5c4aeb6 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/UrQMD/UrQMD.cc b/Processes/UrQMD/UrQMD.cc
index cc9e02f61..cec9beb3e 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 3461c8ca3..937b6d702 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;
--
GitLab