add proposal sec to corsika stack

Processes/Proposal/Interaction.cc

@@ -19,13 +19,11 @@ namespace corsika::process::proposal {
   using namespace corsika::setup;
   using namespace corsika::environment;
   using namespace corsika::units::si;
-  using namespace corsika;
-
   using SetupParticle = corsika::setup::Stack::StackIterator;
 
   template <>
   std::unordered_map<particles::Code, PROPOSAL::ParticleDef>
-      Interaction<SetupEnvironment>::particle_map{
+      Interaction<SetupEnvironment>::particles{
           {particles::Code::Gamma, PROPOSAL::GammaDef()},
           {particles::Code::Electron, PROPOSAL::EMinusDef()},
           {particles::Code::Positron, PROPOSAL::EPlusDef()},
@@ -40,8 +38,10 @@ namespace corsika::process::proposal {
                                              CORSIKA_ParticleCut const& e_cut)
       : fEnvironment(env)
       , cut(make_shared<const PROPOSAL::EnergyCutSettings>(e_cut.GetCutEnergy() / 1_GeV,
-                                                           0, false)) {
+                                                           1, false)) {}
 
+  template <>
+  void Interaction<SetupEnvironment>::Init() {
     auto all_compositions = std::vector<NuclearComposition>();
     fEnvironment.GetUniverse()->walk([&](auto& vtn) {
       if (vtn.HasModelProperties())
@@ -55,44 +55,24 @@ namespace corsika::process::proposal {
                               *frac_iter);
         ++frac_iter;
       }
-      medium_map[&ncarg] = PROPOSAL::Medium(
+      media[&ncarg] = PROPOSAL::Medium(
           "Modified Air", 1., PROPOSAL::Air().GetI(), PROPOSAL::Air().GetC(),
           PROPOSAL::Air().GetA(), PROPOSAL::Air().GetM(), PROPOSAL::Air().GetX0(),
           PROPOSAL::Air().GetX1(), PROPOSAL::Air().GetD0(), 1.0, comp_vec);
     }
-    std::cout << "PROPOSAL done." << std::endl;
-  }
-
-  template <>
-  template <>
-  auto Interaction<SetupEnvironment>::GetCalculator(SetupParticle const& vP) {
-    auto& mediumComposition = vP.GetNode()->GetModelProperties().GetNuclearComposition();
-    auto calc_it = calculators.find(&mediumComposition);
-    if (calc_it != calculators.end()) return calc_it;
-    auto const& p_def = particle_map[vP.GetPID()];
-    auto& medium = medium_map[&mediumComposition];
-    auto cross = PROPOSAL::GetStdCrossSections(p_def, medium, cut, true);
-    auto inter_types = PROPOSAL::CrossSectionVector::GetInteractionTypes(cross);
-    auto [insert_it, success] = calculators.insert(
-        {&mediumComposition,
-         make_tuple(PROPOSAL::SecondariesCalculator(inter_types, p_def, medium),
-                    PROPOSAL::make_interaction(cross, true),
-                    PROPOSAL::make_displacement(cross, true))});
-    if (success) return insert_it;
-    throw std::logic_error("insert failed.");
   }
 
   template <>
   template <>
   corsika::process::EProcessReturn Interaction<SetupEnvironment>::DoInteraction(
-      SetupParticle const& vP) {
+      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 (auto i = 0; i < std::get<SECONDARIES>(calc->second).RequiredRandomNumbers(type);
-         ++i)
+    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 =
@@ -114,11 +94,14 @@ namespace corsika::process::proposal {
           corsika::geometry::RootCoordinateSystem::GetInstance()
               .GetRootCoordinateSystem(),
           vec);
-      // particles::Code, units::si::HEPEnergyType, stack::MomentumVector,
-      // geometry::Point, units::si::TimeType
-      auto sec = make_tuple(get<PROPOSAL::Loss::TYPE>(s), energy, momentum,
-                            vP.GetPosition(), vP.GetTime());
-      vP.AddSecondary(sec);
+      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));
     }
     return process::EProcessReturn::eOk;
   }

Processes/Proposal/Interaction.h

@@ -27,16 +27,6 @@ using CORSIKA_ParticleCut = corsika::process::particle_cut::ParticleCut;
 
 namespace corsika::process::proposal {
 
-  /* static std::unordered_map<particles::Code, PROPOSAL::ParticleDef> particle_map{ */
-  /*     {particles::Code::Gamma, PROPOSAL::GammaDef()}, */
-  /*     {particles::Code::Electron, PROPOSAL::EMinusDef()}, */
-  /*     {particles::Code::Positron, PROPOSAL::EPlusDef()}, */
-  /*     {particles::Code::MuMinus, PROPOSAL::MuMinusDef()}, */
-  /*     {particles::Code::MuPlus, PROPOSAL::MuPlusDef()}, */
-  /*     {particles::Code::TauPlus, PROPOSAL::TauPlusDef()}, */
-  /*     {particles::Code::TauMinus, PROPOSAL::TauMinusDef()}, */
-  /* }; */
-
   template <class TEnvironment>
   class Interaction
       : public corsika::process::InteractionProcess<Interaction<TEnvironment>> {
@@ -45,26 +35,17 @@ namespace corsika::process::proposal {
     TEnvironment const& fEnvironment;
     shared_ptr<const PROPOSAL::EnergyCutSettings> cut;
 
-    static std::unordered_map<particles::Code, PROPOSAL::ParticleDef> particle_map;
-    std::unordered_map<const NuclearComposition*, PROPOSAL::Medium> medium_map;
+    static std::unordered_map<particles::Code, PROPOSAL::ParticleDef> particles;
+    std::unordered_map<const NuclearComposition*, PROPOSAL::Medium> media;
 
     enum { SECONDARIES, INTERACTION, DISPLACEMENT };
-    /* std::uniform_real_distribution<> rnd_uniform(0., 1.); */
-
-    /* std::map<particles::Code, std::unique_ptr<PROPOSAL::UtilityInterpolantInteraction>>
-     * corsika_particle_to_utility_map; */
-    /* std::map<particles::Code,
-     * std::unique_ptr<PROPOSAL::UtilityInterpolantDisplacement>>
-     * corsika_particle_to_displacement_map; */
 
     corsika::random::RNG& fRNG =
         corsika::random::RNGManager::GetInstance().GetRandomStream("s_rndm");
 
-    // fTrackedParticles[proposal_particle.GetName()] return  particle::Code
-
     auto IsTracked(particles::Code pcode) const noexcept {
-      auto search = particle_map.find(pcode);
-      if (search != particle_map.end()) return true;
+      auto search = particles.find(pcode);
+      if (search != particles.end()) return true;
       return false;
     };
 
@@ -74,11 +55,26 @@ namespace corsika::process::proposal {
     std::unordered_map<const NuclearComposition*, calculator_t> calculators;
 
     template <typename Particle>
-    auto GetCalculator(Particle&);
+    auto GetCalculator(Particle& vP) {
+      auto& comp = vP.GetNode()->GetModelProperties().GetNuclearComposition();
+      auto calc_it = calculators.find(&comp);
+      if (calc_it != calculators.end()) return calc_it;
+      auto cross =
+          PROPOSAL::GetStdCrossSections(particles[vP.GetPID()], media[&comp], cut, true);
+      auto inter_types = PROPOSAL::CrossSectionVector::GetInteractionTypes(cross);
+      auto [insert_it, success] = calculators.insert(
+          {&comp, make_tuple(PROPOSAL::SecondariesCalculator(
+                                 inter_types, particles[vP.GetPID()], media[&comp]),
+                             PROPOSAL::make_interaction(cross, true),
+                             PROPOSAL::make_displacement(cross, true))});
+      return insert_it;
+    }
 
   public:
     Interaction(TEnvironment const& env, CORSIKA_ParticleCut const& cut);
 
+    void Init();
+
     template <typename Particle>
     corsika::process::EProcessReturn DoInteraction(Particle&);