Include particle mapping for ContinuousProcess as well as some unit fixing

Processes/Proposal/ContinuousProcess.cc

@@ -36,22 +36,22 @@ namespace corsika::process::proposal {
   template <>
   ContinuousProcess::ContinuousProcess(SetupEnvironment const& _env,
                                        CORSIKA_ParticleCut const& _cut)
-      : cut(make_shared<const PROPOSAL::EnergyCutSettings>(_cut.GetCutEnergy() / 1_GeV, 1,
+      : cut(make_shared<const PROPOSAL::EnergyCutSettings>(_cut.GetECut() / 1_MeV, 1,
                                                            false)) {
-    auto all_compositions = std::vector<NuclearComposition>();
+    auto all_compositions = std::vector<const NuclearComposition*>();
     _env.GetUniverse()->walk([&](auto& vtn) {
       if (vtn.HasModelProperties())
-        all_compositions.push_back(vtn.GetModelProperties().GetNuclearComposition());
+        all_compositions.push_back(&vtn.GetModelProperties().GetNuclearComposition());
     });
     for (auto& ncarg : all_compositions) {
       auto comp_vec = std::vector<PROPOSAL::Components::Component>();
-      auto frac_iter = ncarg.GetFractions().cbegin();
-      for (auto& pcode : ncarg.GetComponents()) {
+      auto frac_iter = ncarg->GetFractions().cbegin();
+      for (auto& pcode : ncarg->GetComponents()) {
         comp_vec.emplace_back(GetName(pcode), GetNucleusZ(pcode), GetNucleusA(pcode),
                               *frac_iter);
         ++frac_iter;
       }
-      media[&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);
@@ -62,14 +62,16 @@ namespace corsika::process::proposal {
                                                    GrammageType const vDX) {
     auto calc_ptr = GetCalculator(vP);
     auto upper_energy = calc_ptr->second->UpperLimitTrackIntegral(
-        vP.GetEnergy() / 1_GeV, vDX / 1_g * 1_cm * 1_cm);
-    return upper_energy * 1_GeV;
+        vP.GetEnergy() / 1_MeV, vDX / 1_g * 1_cm * 1_cm);
+    std::cout << "upper_energy: " << upper_energy << "MeV" << std::endl;
+    return upper_energy * 1_MeV;
   }
 
+   void ContinuousProcess::Init() {}
   template <>
   EProcessReturn ContinuousProcess::DoContinuous(SetupParticle& vP, SetupTrack const& vT) {
     if (vP.GetChargeNumber() == 0) return process::EProcessReturn::eOk;
-
+    std::cout << "DoContinuous..." << std::endl;
     GrammageType const dX =
         vP.GetNode()->GetModelProperties().IntegratedGrammage(vT, vT.GetLength());
     HEPEnergyType dE = TotalEnergyLoss(vP, dX);

Processes/Proposal/ContinuousProcess.h

@@ -42,18 +42,45 @@ namespace corsika::process::proposal {
 
     bool CanInteract(particles::Code pcode) const noexcept;
 
-    unordered_map<const NuclearComposition*, unique_ptr<PROPOSAL::Displacement>> calc;
-
+    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;
+        }
+    };
+
+    unordered_map<std::pair<const NuclearComposition*, particles::Code>, unique_ptr<PROPOSAL::Displacement>, interaction_hash> calc;
+
+    auto BuildCalculator(particles::Code corsika_code, NuclearComposition const& comp) {
+        auto medium = media.at(&comp);
+        if (corsika_code == particles::Code::Gamma) {
+            auto cross = GetStdCrossSections(PROPOSAL::GammaDef(), media.at(&comp), cut, true);
+            auto [insert_it, success] =
+                    calc.insert({std::make_pair(&comp, corsika_code), PROPOSAL::make_displacement(cross, true)});
+            return insert_it;
+        }
+        if (corsika_code == particles::Code::Electron) {
+            auto cross = GetStdCrossSections(PROPOSAL::EMinusDef(), media.at(&comp), cut, true);
+            auto [insert_it, success] =
+            calc.insert({std::make_pair(&comp, corsika_code), PROPOSAL::make_displacement(cross, true)});
+            return insert_it;
+        }
+        if (corsika_code == particles::Code::Positron) {
+            auto cross = GetStdCrossSections(PROPOSAL::EPlusDef(), media.at(&comp), cut, true);
+            auto [insert_it, success] =
+            calc.insert({std::make_pair(&comp, corsika_code), PROPOSAL::make_displacement(cross, true)});
+            return insert_it;
+        }
+    }
+    
     template <typename Particle>
     auto GetCalculator(Particle& vP) {
       auto& comp = vP.GetNode()->GetModelProperties().GetNuclearComposition();
-      auto calc_it = calc.find(&comp);
+      auto calc_it = calc.find(std::make_pair(&comp, vP.GetPID()));
       if (calc_it != calc.end()) return calc_it;
-      auto cross =
-          PROPOSAL::GetStdCrossSections(particles[vP.GetPID()], media[&comp], cut, true);
-      auto [insert_it, success] =
-          calc.insert({&comp, PROPOSAL::make_displacement(cross, true)});
-      return insert_it;
+      return BuildCalculator(vP.GetPID(), comp);
     }
 
     units::si::HEPEnergyType TotalEnergyLoss(setup::Stack::ParticleType const&,
@@ -63,6 +90,8 @@ namespace corsika::process::proposal {
     template <typename TEnvironment>
     ContinuousProcess(TEnvironment const& env, CORSIKA_ParticleCut const& cut);
 
+    void Init();
+
     template <typename Particle, typename Track>
     EProcessReturn DoContinuous(Particle&, Track const&) ;