diff --git a/Framework/Cascade/Cascade.h b/Framework/Cascade/Cascade.h index 490b0408620982515036391ce539b81edc6ca537..d83738c5e01bfff9ecd37bc93c9ff6ea213fac99 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 d95c69d657cb36709ad25ae6febb3e293709e9dc..ece6447e8f862d46d73efe2e28ec9b8c67005b74 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 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/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/QGSJetII/Interaction.cc b/Processes/QGSJetII/Interaction.cc index 4f2f1ca34ecbc6bbcc175040e9e99387f942e1ff..40aac437d8d74c5b28ec3a958786c321af4a6248 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 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/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/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;