diff --git a/Processes/Proposal/ContinuousProcess.cc b/Processes/Proposal/ContinuousProcess.cc
index d222a80b121e21c7de1de50e03fd1df820f026f4..53eddcc85e810cd6cd1c00b99ee1033931209d65 100644
--- a/Processes/Proposal/ContinuousProcess.cc
+++ b/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;
diff --git a/Processes/Proposal/ContinuousProcess.h b/Processes/Proposal/ContinuousProcess.h
index e0cfe6e2679fda8389c44a230160c571a81b0b9b..7a42c122430eac408bf6939152d71edeba48dc9f 100644
--- a/Processes/Proposal/ContinuousProcess.h
+++ b/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