add flag if PROPOSAL can track the particle

Processes/Proposal/CMakeLists.txt

@@ -15,6 +15,11 @@ SET_TARGET_PROPERTIES ( ProcessProposal PROPERTIES VERSION ${PROJECT_VERSION}
 
 TARGET_LINK_LIBRARIES (
   ProcessProposal
+  CORSIKAparticles
+  CORSIKAutilities
+  CORSIKAunits
+  CORSIKAthirdparty
+  CORSIKAgeometry
   CORSIKAenvironment
   ${PROPOSAL_LIBRARY}
   )

Processes/Proposal/Interaction.cc

@@ -2,6 +2,7 @@
 #include <corsika/environment/IMediumModel.h>
 #include <corsika/environment/NuclearComposition.h>
 #include <corsika/process/proposal/Interaction.h>
+#include <corsika/setup/SetupEnvironment.h>
 #include <corsika/setup/SetupStack.h>
 #include <corsika/setup/SetupTrajectory.h>
 #include <corsika/units/PhysicalUnits.h>
@@ -14,12 +15,11 @@
 #include <tuple>
 
 using Component_PROPOSAL = PROPOSAL::Components::Component;
-namespace corsika::process::proposal {
 
+namespace corsika::process::proposal {
   using namespace corsika::setup;
   using namespace corsika::environment;
   using namespace corsika::units::si;
-  using SetupParticle = corsika::setup::Stack::StackIterator;
 
   template <>
   std::unordered_map<particles::Code, PROPOSAL::ParticleDef>
@@ -65,43 +65,42 @@ namespace corsika::process::proposal {
   template <>
   template <>
   corsika::process::EProcessReturn Interaction<SetupEnvironment>::DoInteraction(
-      SetupParticle& vP) {
-    auto calc = GetCalculator(vP); // [CrossSections]
-    std::uniform_real_distribution<double> distr(0., 1.);
-    auto [type, comp_ptr, v] = std::get<INTERACTION>(calc->second)
-                                   ->TypeInteraction(vP.GetEnergy() / 1_GeV, distr(fRNG));
-    auto rnd = std::vector<double>();
-    for (size_t i = 0;
-         i < std::get<SECONDARIES>(calc->second).RequiredRandomNumbers(type); ++i)
-      rnd.push_back(distr(fRNG));
-    double primary_energy = vP.GetEnergy() / 1_GeV;
-    auto point =
-        PROPOSAL::Vector3D(vP.GetPosition().GetX() / 1_cm, vP.GetPosition().GetY() / 1_cm,
-                           vP.GetPosition().GetZ() / 1_cm);
-    auto p = vP.GetMomentum().GetComponents();
-    auto direction = PROPOSAL::Vector3D(p[0] / 1_GeV, p[1] / 1_GeV, p[2] / 1_GeV);
-    auto loss =
-        make_tuple(static_cast<int>(type), point, direction, v * primary_energy, 0.);
-    auto sec = std::get<SECONDARIES>(calc->second)
-                   .CalculateSecondaries(primary_energy, loss, *comp_ptr, rnd);
-    for (auto& s : sec) {
-      auto energy = get<PROPOSAL::Loss::ENERGY>(s) * 1_MeV;
-      auto vec = corsika::geometry::QuantityVector(
-          std::get<PROPOSAL::Loss::DIRECTION>(s).GetX() * energy,
-          std::get<PROPOSAL::Loss::DIRECTION>(s).GetY() * energy,
-          std::get<PROPOSAL::Loss::DIRECTION>(s).GetZ() * energy);
-      auto momentum = corsika::stack::MomentumVector(
-          corsika::geometry::RootCoordinateSystem::GetInstance()
-              .GetRootCoordinateSystem(),
-          vec);
-      particles::Code sec_code = corsika::particles::ConvertFromPDG(
-          static_cast<particles::PDGCode>(get<PROPOSAL::Loss::TYPE>(s)));
-      corsika::units::si::HEPEnergyType sec_energy = energy;
-      stack::MomentumVector sec_momentum = momentum;
-      geometry::Point sec_position = vP.GetPosition();
-      corsika::units::si::TimeType sec_time = vP.GetTime();
-      vP.AddSecondary(
-          make_tuple(sec_code, sec_energy, sec_momentum, sec_position, sec_time));
+      setup::StackView::StackIterator& vP) {
+    if (CanInteract(vP.GetPID())) {
+      auto calc = GetCalculator(vP); // [CrossSections]
+      std::uniform_real_distribution<double> distr(0., 1.);
+      auto [type, comp_ptr, v] =
+          std::get<INTERACTION>(calc->second)
+              ->TypeInteraction(vP.GetEnergy() / 1_GeV, distr(fRNG));
+      auto rnd = std::vector<double>();
+      for (size_t i = 0;
+           i < std::get<SECONDARIES>(calc->second).RequiredRandomNumbers(type); ++i)
+        rnd.push_back(distr(fRNG));
+      double primary_energy = vP.GetEnergy() / 1_GeV;
+      auto point = PROPOSAL::Vector3D(vP.GetPosition().GetX() / 1_cm,
+                                      vP.GetPosition().GetY() / 1_cm,
+                                      vP.GetPosition().GetZ() / 1_cm);
+      auto p = vP.GetMomentum().GetComponents();
+      auto direction = PROPOSAL::Vector3D(p[0] / 1_GeV, p[1] / 1_GeV, p[2] / 1_GeV);
+      auto loss =
+          make_tuple(static_cast<int>(type), point, direction, v * primary_energy, 0.);
+      auto sec = std::get<SECONDARIES>(calc->second)
+                     .CalculateSecondaries(primary_energy, loss, *comp_ptr, rnd);
+      for (auto& s : sec) {
+        auto energy = get<PROPOSAL::Loss::ENERGY>(s) * 1_MeV;
+        auto vec = corsika::geometry::QuantityVector(
+            std::get<PROPOSAL::Loss::DIRECTION>(s).GetX() * energy,
+            std::get<PROPOSAL::Loss::DIRECTION>(s).GetY() * energy,
+            std::get<PROPOSAL::Loss::DIRECTION>(s).GetZ() * energy);
+        auto momentum = corsika::stack::MomentumVector(
+            corsika::geometry::RootCoordinateSystem::GetInstance()
+                .GetRootCoordinateSystem(),
+            vec);
+        particles::Code sec_code = corsika::particles::ConvertFromPDG(
+            static_cast<particles::PDGCode>(get<PROPOSAL::Loss::TYPE>(s)));
+        vP.AddSecondary(
+            make_tuple(sec_code, energy, momentum, vP.GetPosition(), vP.GetTime()));
+      }
     }
     return process::EProcessReturn::eOk;
   }
@@ -109,13 +108,16 @@ namespace corsika::process::proposal {
   template <>
   template <>
   corsika::units::si::GrammageType Interaction<SetupEnvironment>::GetInteractionLength(
-      SetupParticle const& vP) {
-    auto calc = GetCalculator(vP); // [CrossSections]
-    std::uniform_real_distribution<double> distr(0., 1.);
-    auto energy = get<INTERACTION>(calc->second)
-                      ->EnergyInteraction(vP.GetEnergy() / 1_GeV, distr(fRNG));
-    return get<DISPLACEMENT>(calc->second)
-               ->SolveTrackIntegral(vP.GetEnergy() / 1_GeV, energy) *
-           1_g / 1_cm / 1_cm;
+      setup::Stack::StackIterator& vP) {
+    if (CanInteract(vP.GetPID())) {
+      auto calc = GetCalculator(vP); // [CrossSections]
+      std::uniform_real_distribution<double> distr(0., 1.);
+      auto energy = get<INTERACTION>(calc->second)
+                        ->EnergyInteraction(vP.GetEnergy() / 1_GeV, distr(fRNG));
+      return get<DISPLACEMENT>(calc->second)
+                 ->SolveTrackIntegral(vP.GetEnergy() / 1_GeV, energy) *
+             1_g / 1_cm / 1_cm;
+    }
+    return std::numeric_limits<double>::infinity() * 1_g / (1_cm * 1_cm);
   }
 } // namespace corsika::process::proposal

Processes/Proposal/Interaction.h

@@ -38,12 +38,10 @@ namespace corsika::process::proposal {
     static std::unordered_map<particles::Code, PROPOSAL::ParticleDef> particles;
     std::unordered_map<const NuclearComposition*, PROPOSAL::Medium> media;
 
-    enum { SECONDARIES, INTERACTION, DISPLACEMENT };
-
     corsika::random::RNG& fRNG =
-        corsika::random::RNGManager::GetInstance().GetRandomStream("s_rndm");
+        corsika::random::RNGManager::GetInstance().GetRandomStream("p_rndm");
 
-    auto IsTracked(particles::Code pcode) const noexcept {
+    bool CanInteract(particles::Code pcode) const noexcept {
       auto search = particles.find(pcode);
       if (search != particles.end()) return true;
       return false;
@@ -54,6 +52,7 @@ namespace corsika::process::proposal {
               unique_ptr<PROPOSAL::Displacement>>;
     std::unordered_map<const NuclearComposition*, calculator_t> calculators;
 
+    enum { SECONDARIES, INTERACTION, DISPLACEMENT };
     template <typename Particle>
     auto GetCalculator(Particle& vP) {
       auto& comp = vP.GetNode()->GetModelProperties().GetNuclearComposition();
@@ -81,6 +80,5 @@ namespace corsika::process::proposal {
     template <typename TParticle>
     corsika::units::si::GrammageType GetInteractionLength(TParticle& p);
   };
-
 } // namespace corsika::process::proposal
 #endif