naiv implementation of continuous em energy loss

Processes/Proposal/ContinuousProcess.cc

-// #include <PROPOSAL/PROPOSAL.h>
-// #include <corsika/process/proposal/ContinuousProcess.h>
-// // #include <corsika/process/sibyll/Interaction.h>
+#include <PROPOSAL/PROPOSAL.h>
+#include <corsika/environment/IMediumModel.h>
+#include <corsika/environment/NuclearComposition.h>
+#include <corsika/process/proposal/ContinuousProcess.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>
+#include <corsika/utl/COMBoost.h>
 
+using SetupParticle = corsika::setup::Stack::ParticleType;
+using SetupTrack = corsika::setup::Trajectory;
 
-// namespace corsika::process::proposal {
+namespace corsika::process::proposal {
+  using namespace corsika::setup;
+  using namespace corsika::environment;
+  using namespace corsika::units::si;
 
-// template <>
-// EProcessReturn ContinuousProcess::DoContinuous(Particle&, Track const&) const {};
+  unordered_map<particles::Code, PROPOSAL::ParticleDef> ContinuousProcess::particles{
+      {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 <>
-// units::si::LengthType ContinuousProcess::MaxStepLength(Particle const& p, Track const& track) const {};
+  bool ContinuousProcess::CanInteract(particles::Code pcode) const noexcept {
+    auto search = particles.find(pcode);
+    if (search != particles.end()) return true;
+    return false;
+  }
 
-// } // namespace corsika::process
+  template <>
+  ContinuousProcess::ContinuousProcess(SetupEnvironment const& _env,
+                                       CORSIKA_ParticleCut const& _cut)
+      : cut(make_shared<const PROPOSAL::EnergyCutSettings>(_cut.GetCutEnergy() / 1_GeV, 1,
+                                                           false)) {
+    auto all_compositions = std::vector<NuclearComposition>();
+    _env.GetUniverse()->walk([&](auto& vtn) {
+      if (vtn.HasModelProperties())
+        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()) {
+        comp_vec.emplace_back(GetName(pcode), GetNucleusZ(pcode), GetNucleusA(pcode),
+                              *frac_iter);
+        ++frac_iter;
+      }
+      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);
+    }
+  }
+
+  HEPEnergyType ContinuousProcess::TotalEnergyLoss(SetupParticle const& vP,
+                                                   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;
+  }
+
+  template <>
+  EProcessReturn ContinuousProcess::DoContinuous(SetupParticle& vP, SetupTrack const& vT) {
+    if (vP.GetChargeNumber() == 0) return process::EProcessReturn::eOk;
+
+    GrammageType const dX =
+        vP.GetNode()->GetModelProperties().IntegratedGrammage(vT, vT.GetLength());
+    HEPEnergyType dE = TotalEnergyLoss(vP, dX);
+    auto E = vP.GetEnergy();
+    const auto Ekin = E - vP.GetMass();
+    auto Enew = E + dE;
+    auto status = process::EProcessReturn::eOk;
+    if (-dE > Ekin) {
+      dE = -Ekin;
+      Enew = vP.GetMass();
+      status = process::EProcessReturn::eParticleAbsorbed;
+    }
+    vP.SetEnergy(Enew);
+    auto pnew = vP.GetMomentum();
+    vP.SetMomentum(pnew * Enew / pnew.GetNorm());
+    return status;
+  }
+
+  template <>
+  units::si::LengthType ContinuousProcess::MaxStepLength(SetupParticle const& vP,
+                                                         SetupTrack const& vT) {
+    auto constexpr dX = 1_g / square(1_cm);
+    auto const dE = -TotalEnergyLoss(vP, dX); // dE > 0
+    auto const maxLoss = 0.01 * vP.GetEnergy();
+    auto const maxGrammage = maxLoss / dE * dX;
+
+    return vP.GetNode()->GetModelProperties().ArclengthFromGrammage(vT,
+                                                                    maxGrammage) *
+           1.0001; // to make sure particle gets absorbed when DoContinuous() is called
+  }
+
+} // namespace corsika::process::proposal

Processes/Proposal/ContinuousProcess.h

@@ -8,26 +8,67 @@
  * the license.
  */
 
-// #ifndef _corsika_process_proposal_interaction_h_
-// #define _corsika_process_proposalythia_interaction_h_
+#ifndef _corsika_process_proposal_interaction_h_
+#define _corsika_process_proposal_interaction_h_
 
-// #include <PROPOSAL/PROPOSAL.h>
-// #include <corsika/process/ContinuousProcess.h>
+#include <PROPOSAL/PROPOSAL.h>
+#include <corsika/environment/Environment.h>
+#include <corsika/particles/ParticleProperties.h>
+#include <corsika/process/ContinuousProcess.h>
+#include <corsika/process/particle_cut/ParticleCut.h>
+#include <corsika/random/RNGManager.h>
+#include <corsika/random/UniformRealDistribution.h>
+#include <unordered_map>
+#include "PROPOSAL/PROPOSAL.h"
 
+using std::unordered_map;
 
-// namespace corsika::process::proposal {
+using namespace corsika::environment;
 
-//   class ContinuousProcess : public corsika::process::ContinuousProcess<Continuous> {
-//   private:
+using CORSIKA_ParticleCut = corsika::process::particle_cut::ParticleCut;
 
+namespace corsika::process::proposal {
 
-//   public:
-//     template <typename Particle, typename Track>
-//     EProcessReturn DoContinuous(Particle&, Track const&) const;
+  class ContinuousProcess
+      : public corsika::process::ContinuousProcess<ContinuousProcess> {
+  private:
+    shared_ptr<const PROPOSAL::EnergyCutSettings> cut;
 
-//     template <typename Particle, typename Track>
-//     units::si::LengthType MaxStepLength(Particle const& p, Track const& track) const;
-//   }
-// } // namespace corsika::process
+    static unordered_map<particles::Code, PROPOSAL::ParticleDef> particles;
+    unordered_map<const NuclearComposition*, PROPOSAL::Medium> media;
 
-// #endif
+    corsika::random::RNG& fRNG =
+        corsika::random::RNGManager::GetInstance().GetRandomStream("proposal");
+
+    bool CanInteract(particles::Code pcode) const noexcept;
+
+    unordered_map<const NuclearComposition*, unique_ptr<PROPOSAL::Displacement>> calc;
+
+    template <typename Particle>
+    auto GetCalculator(Particle& vP) {
+      auto& comp = vP.GetNode()->GetModelProperties().GetNuclearComposition();
+      auto calc_it = calc.find(&comp);
+      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;
+    }
+
+    units::si::HEPEnergyType TotalEnergyLoss(setup::Stack::ParticleType const&,
+                                             const units::si::GrammageType);
+
+  public:
+    template <typename TEnvironment>
+    ContinuousProcess(TEnvironment const& env, CORSIKA_ParticleCut const& cut);
+
+    template <typename Particle, typename Track>
+    EProcessReturn DoContinuous(Particle&, Track const&) ;
+
+    template <typename Particle, typename Track>
+    units::si::LengthType MaxStepLength(Particle const& p, Track const& track) ;
+  };
+} // namespace corsika::process::proposal
+
+#endif