Add deflection caused by scattering. Scattering along the track is not yet...

Add deflection caused by scattering. Scattering along the track is not yet possible because the track is a const object.

Processes/Proposal/ContinuousProcess.cc

@@ -59,20 +59,46 @@ namespace corsika::process::proposal {
     }
   }
 
+  template <>
   HEPEnergyType ContinuousProcess::TotalEnergyLoss(SetupParticle const& vP,
-                                                   GrammageType const vDX) {
-    std::cout << "distance: " << vDX / 1_g * 1_cm * 1_cm << std::endl;
-    return vP.GetEnergy() - GetCalculator(vP)->second->UpperLimitTrackIntegral(
-                                vP.GetEnergy() / 1_MeV, vDX / 1_g * 1_cm * 1_cm) *
+                                                   GrammageType const& vDX) {
+    auto calc = GetCalculator(vP);
+    return vP.GetEnergy() - get<DISPLACEMENT>(calc->second)
+                                    ->UpperLimitTrackIntegral(vP.GetEnergy() / 1_MeV,
+                                                              vDX / 1_g * 1_cm * 1_cm) *
                                 1_MeV;
   }
 
+  template <>
+  void ContinuousProcess::Scatter(SetupParticle& vP, HEPEnergyType const& loss,
+                                           GrammageType const& grammage) {
+    auto calc = GetCalculator(vP);
+    auto d = vP.GetDirection().GetComponents();
+    auto direction = PROPOSAL::Vector3D(d.GetX().magnitude(), d.GetY().magnitude(),
+                                        d.GetZ().magnitude());
+    auto E_f = vP.GetEnergy() - loss; // final energy
+    std::uniform_real_distribution<double> distr(0., 1.);
+    auto rnd = array<double, 4>();
+    for (auto& it : rnd) it = distr(fRNG);
+    auto [mean_dir, final_dir] =
+        get<SCATTERING>(calc->second)
+            ->Scatter(grammage / 1_g * square(1_cm), vP.GetEnergy() / 1_MeV, E_f / 1_MeV,
+                      direction, rnd);
+    auto vec = corsika::geometry::QuantityVector(
+        final_dir.GetX() * E_f, final_dir.GetY() * E_f, final_dir.GetZ() * E_f);
+    vP.SetMomentum(corsika::stack::MomentumVector(
+        corsika::geometry::RootCoordinateSystem::GetInstance().GetRootCoordinateSystem(),
+        vec));
+  }
+
   template <>
   EProcessReturn ContinuousProcess::DoContinuous(SetupParticle& vP,
                                                  SetupTrack const& vT) {
     if (vP.GetChargeNumber() == 0) return process::EProcessReturn::eOk;
     auto dX = vP.GetNode()->GetModelProperties().IntegratedGrammage(vT, vT.GetLength());
-    vP.SetEnergy(vP.GetEnergy() - TotalEnergyLoss(vP, dX));
+    auto energy_loss = TotalEnergyLoss(vP, dX);
+    Scatter(vP, energy_loss, dX);
+    vP.SetEnergy(vP.GetEnergy() - energy_loss);
     if (vP.GetEnergy() < cut.GetECut()) return process::EProcessReturn::eParticleAbsorbed;
     vP.SetMomentum(vP.GetMomentum() * vP.GetEnergy() / vP.GetMomentum().GetNorm());
     return process::EProcessReturn::eOk;
@@ -82,8 +108,9 @@ namespace corsika::process::proposal {
   units::si::LengthType ContinuousProcess::MaxStepLength(SetupParticle const& vP,
                                                          SetupTrack const& vT) {
     auto energy_lim = 0.9 * vP.GetEnergy() / 1_MeV;
-    auto grammage = GetCalculator(vP)->second->SolveTrackIntegral(vP.GetEnergy() / 1_MeV,
-                                                                  energy_lim) *
+    auto calc = GetCalculator(vP);
+    auto grammage = get<DISPLACEMENT>(calc->second)
+                        ->SolveTrackIntegral(vP.GetEnergy() / 1_MeV, energy_lim) *
                     1_g / square(1_cm);
     return vP.GetNode()->GetModelProperties().ArclengthFromGrammage(vT, grammage) *
            1.0001;

Processes/Proposal/ContinuousProcess.h

@@ -37,10 +37,11 @@ namespace corsika::process::proposal {
     static unordered_map<particles::Code, PROPOSAL::ParticleDef> particles;
     unordered_map<const NuclearComposition*, PROPOSAL::Medium> media;
 
-
     bool CanInteract(particles::Code pcode) const noexcept;
 
     using calc_key_t = std::pair<const NuclearComposition*, particles::Code>;
+    using calc_t =
+        tuple<unique_ptr<PROPOSAL::Displacement>, unique_ptr<PROPOSAL::Scattering>>;
 
     struct disp_hash {
       size_t operator()(const calc_key_t& p) const {
@@ -49,18 +50,21 @@ namespace corsika::process::proposal {
       }
     };
 
-    unordered_map<calc_key_t, unique_ptr<PROPOSAL::Displacement>, disp_hash> calc;
+    enum { DISPLACEMENT, SCATTERING };
+    unordered_map<calc_key_t, calc_t, disp_hash> calc;
 
     template <typename Particle>
     auto BuildCalculator(particles::Code code, Particle p_def,
                          NuclearComposition const& comp) {
-      auto medium = media.at(&comp);
       auto cross = GetStdCrossSections(
           p_def, media.at(&comp),
           make_shared<const PROPOSAL::EnergyCutSettings>(cut.GetECut() / 1_MeV, 1, false),
           true);
-      auto [insert_it, success] = calc.insert(
-          {std::make_pair(&comp, code), PROPOSAL::make_displacement(cross, true)});
+      auto [insert_it, success] =
+          calc.insert({std::make_pair(&comp, code),
+                       std::make_tuple(PROPOSAL::make_displacement(cross, true),
+                                       PROPOSAL::make_scattering("highland", p_def,
+                                                                 media.at(&comp)))});
       return insert_it;
     }
 
@@ -90,20 +94,25 @@ namespace corsika::process::proposal {
       return BuildCalculator(vP.GetPID(), comp);
     }
 
-    units::si::HEPEnergyType TotalEnergyLoss(setup::Stack::ParticleType const&,
-                                             const units::si::GrammageType);
-
   public:
     template <typename TEnvironment>
-    ContinuousProcess(TEnvironment const& env, CORSIKA_ParticleCut& cut);
+    ContinuousProcess(TEnvironment const&, CORSIKA_ParticleCut&);
 
     void Init(){};
 
+    template <typename Particle>
+    corsika::units::si::HEPEnergyType TotalEnergyLoss(
+        Particle const&, corsika::units::si::GrammageType const&);
+
+    template <typename Particle>
+    void Scatter(Particle&, corsika::units::si::HEPEnergyType const&,
+                 corsika::units::si::GrammageType const&);
+
     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);
+    corsika::units::si::LengthType MaxStepLength(Particle const&, Track const&);
   };
 } // namespace corsika::process::proposal