Add mapping of CORSIKA particles to PROPOSAL calculators

Documentation/Examples/proposal_example.cc

@@ -27,7 +27,7 @@
 #include <corsika/setup/SetupTrajectory.h>
 #include <corsika/units/PhysicalUnits.h>
 #include <corsika/utl/CorsikaFenv.h>
-
+#include <corsika/process/track_writer/TrackWriter.h>
 #include <corsika/process/proposal/Interaction.h>
 
 #include <iomanip>
@@ -138,6 +138,7 @@ int main(int argc, char** argv) {
   process::particle_cut::ParticleCut cut(10_GeV);
   process::proposal::Interaction proposal(env, cut);
   process::interaction_counter::InteractionCounter proposalCounted(proposal);
+  process::track_writer::TrackWriter trackWriter("tracks.dat");
 
   // energy cut; n.b. ParticleCut needs to be modified not to discard EM particles!
   /* process::particle_cut::ParticleCut cut{60_GeV}; */
@@ -149,8 +150,7 @@ int main(int argc, char** argv) {
   process::observation_plane::ObservationPlane observationLevel(obsPlane,
                                                                 "particles.dat");
 
-  auto sequence = proposalCounted << longprof << proposal << cut << observationLevel;
-
+  auto sequence = proposalCounted << longprof << proposal << cut << observationLevel << trackWriter;
   // define air shower object, run simulation
   tracking_line::TrackingLine tracking;
   cascade::Cascade EAS(env, tracking, sequence, stack);

Processes/Proposal/Interaction.cc

@@ -76,7 +76,7 @@ namespace corsika::process::proposal {
       std::cout << "InteractionType: "<< static_cast<int>(type) << std::endl;
       auto rnd = std::vector<double>();
       for (size_t i = 0;
-           i < std::get<SECONDARIES>(calc->second).RequiredRandomNumbers(type); ++i)
+           i < std::get<SECONDARIES>(calc->second)->RequiredRandomNumbers(type); ++i)
         rnd.push_back(distr(fRNG));
       double primary_energy = vP.GetEnergy() / 1_MeV;
       auto point = PROPOSAL::Vector3D(vP.GetPosition().GetX() / 1_cm,
@@ -89,7 +89,7 @@ namespace corsika::process::proposal {
       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);
+                     ->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(

Processes/Proposal/Interaction.h

@@ -39,9 +39,19 @@ namespace corsika::process::proposal {
     bool CanInteract(particles::Code pcode) const noexcept;
 
     using calculator_t =
-        tuple<PROPOSAL::SecondariesCalculator, unique_ptr<PROPOSAL::Interaction>,
+        tuple<unique_ptr<PROPOSAL::SecondariesCalculator>, unique_ptr<PROPOSAL::Interaction>,
               unique_ptr<PROPOSAL::Displacement>>;
-    std::unordered_map<const NuclearComposition*, calculator_t> calculators;
+
+    struct interaction_hash {
+        size_t operator()(const std::pair<const NuclearComposition*, particles::Code>& p) const
+        {
+            auto hash1 = std::hash<const NuclearComposition*>{}(p.first);
+            auto hash2 = std::hash<particles::Code>{}(p.second);
+            return hash1 ^ hash2;
+        }
+    };
+
+    std::unordered_map<std::pair<const NuclearComposition*, particles::Code>, calculator_t, interaction_hash> calculators;
 
 
     auto BuildCalculator(particles::Code corsika_code, NuclearComposition const& comp) {
@@ -52,8 +62,8 @@ namespace corsika::process::proposal {
             GetStdCrossSections(PROPOSAL::GammaDef(), media.at(&comp), cut, true);
         auto inter_types = PROPOSAL::CrossSectionVector::GetInteractionTypes(cross);
         auto [insert_it, success] = calculators.insert(
-            {&comp, make_tuple(PROPOSAL::SecondariesCalculator(
-                                   inter_types, PROPOSAL::GammaDef(), media[&comp]),
+            {std::make_pair(&comp, corsika_code), make_tuple(PROPOSAL::make_secondaries(
+                                   inter_types, PROPOSAL::GammaDef(), media.at(&comp)),
                                PROPOSAL::make_interaction(cross, true),
                                PROPOSAL::make_displacement(cross, true))});
         return insert_it;
@@ -64,8 +74,8 @@ namespace corsika::process::proposal {
             GetStdCrossSections(PROPOSAL::EMinusDef(), media.at(&comp), cut, true);
         auto inter_types = PROPOSAL::CrossSectionVector::GetInteractionTypes(cross);
         auto [insert_it, success] = calculators.insert(
-            {&comp, make_tuple(PROPOSAL::SecondariesCalculator(
-                                   inter_types, PROPOSAL::EMinusDef(), media[&comp]),
+            {std::make_pair(&comp, corsika_code), make_tuple(PROPOSAL::make_secondaries(
+                                   inter_types, PROPOSAL::EMinusDef(), media.at(&comp)),
                                PROPOSAL::make_interaction(cross, true),
                                PROPOSAL::make_displacement(cross, true))});
         return insert_it;
@@ -76,8 +86,8 @@ namespace corsika::process::proposal {
             GetStdCrossSections(PROPOSAL::EPlusDef(), media.at(&comp), cut, true);
         auto inter_types = PROPOSAL::CrossSectionVector::GetInteractionTypes(cross);
         auto [insert_it, success] = calculators.insert(
-            {&comp, make_tuple(PROPOSAL::SecondariesCalculator(
-                                   inter_types, PROPOSAL::EPlusDef(), media[&comp]),
+            {std::make_pair(&comp, corsika_code), make_tuple(PROPOSAL::make_secondaries(
+                                   inter_types, PROPOSAL::EPlusDef(), media.at(&comp)),
                                PROPOSAL::make_interaction(cross, true),
                                PROPOSAL::make_displacement(cross, true))});
         return insert_it;
@@ -88,7 +98,7 @@ namespace corsika::process::proposal {
     template <typename Particle>
     auto GetCalculator(Particle& vP) {
         auto& comp = vP.GetNode()->GetModelProperties().GetNuclearComposition();
-        auto calc_it = calculators.find(&comp);
+        auto calc_it = calculators.find(std::make_pair(&comp, vP.GetPID()));
         if (calc_it != calculators.end()) return calc_it;
         return BuildCalculator(vP.GetPID(), comp);
     }