diff --git a/Documentation/Examples/proposal_example.cc b/Documentation/Examples/proposal_example.cc
index 3fb41d1698f02a1cd361088a2b44acf28def9b54..de16fc65db2f361cea0ccb74aa28d8cbe089bf3d 100644
--- a/Documentation/Examples/proposal_example.cc
+++ b/Documentation/Examples/proposal_example.cc
@@ -29,6 +29,7 @@
#include <corsika/utl/CorsikaFenv.h>
#include <corsika/process/track_writer/TrackWriter.h>
#include <corsika/process/proposal/Interaction.h>
+#include <corsika/process/proposal/ContinuousProcess.h>
#include <iomanip>
#include <iostream>
@@ -137,6 +138,7 @@ int main(int argc, char** argv) {
process::particle_cut::ParticleCut cut(10_GeV);
process::proposal::Interaction proposal(env, cut);
+ process::proposal::ContinuousProcess em_continuous(env, cut);
process::interaction_counter::InteractionCounter proposalCounted(proposal);
process::track_writer::TrackWriter trackWriter("tracks.dat");
@@ -150,7 +152,12 @@ int main(int argc, char** argv) {
process::observation_plane::ObservationPlane observationLevel(obsPlane,
"particles.dat");
- auto sequence = proposalCounted << longprof << proposal << cut << observationLevel << trackWriter;
+ auto sequence = proposalCounted
+ << cut
+ << em_continuous
+ << longprof
+ << observationLevel
+ << trackWriter;
// define air shower object, run simulation
tracking_line::TrackingLine tracking;
cascade::Cascade EAS(env, tracking, sequence, stack);
diff --git a/Processes/Proposal/ContinuousProcess.cc b/Processes/Proposal/ContinuousProcess.cc
index 5914d529606bf264da310ad46eed8717ddbe6e7e..d222a80b121e21c7de1de50e03fd1df820f026f4 100644
--- a/Processes/Proposal/ContinuousProcess.cc
+++ b/Processes/Proposal/ContinuousProcess.cc
@@ -1,6 +1,7 @@
#include <PROPOSAL/PROPOSAL.h>
#include <corsika/environment/IMediumModel.h>
#include <corsika/environment/NuclearComposition.h>
+#include <corsika/process/particle_cut/ParticleCut.h>
#include <corsika/process/proposal/ContinuousProcess.h>
#include <corsika/process/proposal/Interaction.h>
#include <corsika/setup/SetupEnvironment.h>
@@ -28,20 +29,20 @@ namespace corsika::process::proposal {
};
bool ContinuousProcess::CanInteract(particles::Code pcode) const noexcept {
- auto search = particles.find(pcode);
- if (search != particles.end()) return true;
+ if (particles.find(pcode) != particles.end()) return true;
return false;
}
template <>
ContinuousProcess::ContinuousProcess(SetupEnvironment const& _env,
- CORSIKA_ParticleCut const& _cut)
- : cut(make_shared<const PROPOSAL::EnergyCutSettings>(_cut.GetECut() / 1_MeV, 1,
- false)) {
+ CORSIKA_ParticleCut& _cut)
+ : cut(_cut)
+ , fRNG(corsika::random::RNGManager::GetInstance().GetRandomStream("proposal")) {
auto all_compositions = std::vector<const NuclearComposition*>();
_env.GetUniverse()->walk([&](auto& vtn) {
- if (vtn.HasModelProperties())
+ if (vtn.HasModelProperties()) {
all_compositions.push_back(&vtn.GetModelProperties().GetNuclearComposition());
+ }
});
for (auto& ncarg : all_compositions) {
auto comp_vec = std::vector<PROPOSAL::Components::Component>();
@@ -60,47 +61,32 @@ namespace corsika::process::proposal {
HEPEnergyType ContinuousProcess::TotalEnergyLoss(SetupParticle const& vP,
GrammageType const vDX) {
- auto calc_ptr = GetCalculator(vP);
- auto upper_energy = calc_ptr->second->UpperLimitTrackIntegral(
- 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;
+ 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) *
+ 1_MeV;
}
- void ContinuousProcess::Init() {}
template <>
- EProcessReturn ContinuousProcess::DoContinuous(SetupParticle& vP, SetupTrack const& vT) {
+ 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);
- 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;
+ auto dX = vP.GetNode()->GetModelProperties().IntegratedGrammage(vT, vT.GetLength());
+ vP.SetEnergy(vP.GetEnergy() - TotalEnergyLoss(vP, dX));
+ if (vP.GetEnergy() < cut.GetECut()) return process::EProcessReturn::eParticleAbsorbed;
+ vP.SetMomentum(vP.GetMomentum() * vP.GetEnergy() / vP.GetMomentum().GetNorm());
+ return process::EProcessReturn::eOk;
}
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
+ auto energy_lim = 0.9 * vP.GetEnergy() / 1_MeV;
+ auto grammage = GetCalculator(vP)->second->SolveTrackIntegral(vP.GetEnergy() / 1_MeV,
+ energy_lim) *
+ 1_g / square(1_cm);
+ return vP.GetNode()->GetModelProperties().ArclengthFromGrammage(vT, grammage) *
+ 1.0001;
}
} // namespace corsika::process::proposal
diff --git a/Processes/Proposal/ContinuousProcess.h b/Processes/Proposal/ContinuousProcess.h
index 58b19ec2910a61d84ef7965a6cdbe50a4d6ea58f..e0cfe6e2679fda8389c44a230160c571a81b0b9b 100644
--- a/Processes/Proposal/ContinuousProcess.h
+++ b/Processes/Proposal/ContinuousProcess.h
@@ -24,6 +24,7 @@
using std::unordered_map;
using namespace corsika::environment;
+using namespace corsika::units::si;
using CORSIKA_ParticleCut = corsika::process::particle_cut::ParticleCut;
@@ -31,75 +32,56 @@ namespace corsika::process::proposal {
class ContinuousProcess
: public corsika::process::ContinuousProcess<ContinuousProcess> {
- private:
- shared_ptr<const PROPOSAL::EnergyCutSettings> cut;
-
+ CORSIKA_ParticleCut& cut;
+ corsika::random::RNG& fRNG;
static unordered_map<particles::Code, PROPOSAL::ParticleDef> particles;
unordered_map<const NuclearComposition*, PROPOSAL::Medium> media;
- corsika::random::RNG& fRNG =
- corsika::random::RNGManager::GetInstance().GetRandomStream("proposal");
bool CanInteract(particles::Code pcode) const noexcept;
- 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;
- }
+ using calc_key_t = std::pair<const NuclearComposition*, particles::Code>;
+
+ struct disp_hash {
+ size_t operator()(const calc_key_t& p) const {
+ return std::hash<const NuclearComposition*>{}(p.first) ^
+ std::hash<particles::Code>{}(p.second);
+ }
};
- unordered_map<std::pair<const NuclearComposition*, particles::Code>, unique_ptr<PROPOSAL::Displacement>, interaction_hash> calc;
+ unordered_map<calc_key_t, unique_ptr<PROPOSAL::Displacement>, 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)});
+ return insert_it;
+ }
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;
- }
- if (corsika_code == particles::Code::MuMinus) {
- auto cross = GetStdCrossSections(PROPOSAL::MuMinusDef(), 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::MuPlus) {
- auto cross = GetStdCrossSections(PROPOSAL::MuPlusDef(), 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::TauMinus) {
- auto cross = GetStdCrossSections(PROPOSAL::TauMinusDef(), 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::TauPlus) {
- auto cross = GetStdCrossSections(PROPOSAL::TauPlusDef(), 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;
- }
- throw std::runtime_error("PROPOSAL could not find corresponding builder");
+ if (corsika_code == particles::Code::Gamma)
+ return BuildCalculator(particles::Code::Gamma, PROPOSAL::GammaDef(), comp);
+ if (corsika_code == particles::Code::Electron)
+ return BuildCalculator(particles::Code::Electron, PROPOSAL::EMinusDef(), comp);
+ if (corsika_code == particles::Code::Positron)
+ return BuildCalculator(particles::Code::Positron, PROPOSAL::EPlusDef(), comp);
+ if (corsika_code == particles::Code::MuMinus)
+ return BuildCalculator(particles::Code::MuMinus, PROPOSAL::MuMinusDef(), comp);
+ if (corsika_code == particles::Code::MuPlus)
+ return BuildCalculator(particles::Code::MuPlus, PROPOSAL::MuPlusDef(), comp);
+ if (corsika_code == particles::Code::TauMinus)
+ return BuildCalculator(particles::Code::TauMinus, PROPOSAL::TauMinusDef(), comp);
+ if (corsika_code == particles::Code::TauPlus)
+ return BuildCalculator(particles::Code::TauPlus, PROPOSAL::TauPlusDef(), comp);
+ throw std::runtime_error("PROPOSAL could not find corresponding builder");
}
-
+
template <typename Particle>
auto GetCalculator(Particle& vP) {
auto& comp = vP.GetNode()->GetModelProperties().GetNuclearComposition();
@@ -113,15 +95,15 @@ namespace corsika::process::proposal {
public:
template <typename TEnvironment>
- ContinuousProcess(TEnvironment const& env, CORSIKA_ParticleCut const& cut);
+ ContinuousProcess(TEnvironment const& env, CORSIKA_ParticleCut& cut);
- void Init();
+ void Init(){};
template <typename Particle, typename Track>
- EProcessReturn DoContinuous(Particle&, Track const&) ;
+ EProcessReturn DoContinuous(Particle&, Track const&);
template <typename Particle, typename Track>
- units::si::LengthType MaxStepLength(Particle const& p, Track const& track) ;
+ units::si::LengthType MaxStepLength(Particle const& p, Track const& track);
};
} // namespace corsika::process::proposal
diff --git a/Processes/Proposal/Interaction.cc b/Processes/Proposal/Interaction.cc
index 4394ae868f89dfce4871cc8c077f6eb2fb00cb34..60ae3c81f64fdeeafed295838272e35357f58ccd 100644
--- a/Processes/Proposal/Interaction.cc
+++ b/Processes/Proposal/Interaction.cc
@@ -35,9 +35,8 @@ namespace corsika::process::proposal {
}
template <>
- Interaction::Interaction(SetupEnvironment const& _env, CORSIKA_ParticleCut const& _cut)
- : cut(make_shared<const PROPOSAL::EnergyCutSettings>(_cut.GetECut() / 1_MeV, 1,
- false))
+ Interaction::Interaction(SetupEnvironment const& _env, CORSIKA_ParticleCut& _cut)
+ : cut(_cut)
, fRNG(corsika::random::RNGManager::GetInstance().GetRandomStream("proposal")) {
auto all_compositions = std::vector<const NuclearComposition*>();
_env.GetUniverse()->walk([&](auto& vtn) {
diff --git a/Processes/Proposal/Interaction.h b/Processes/Proposal/Interaction.h
index 54597a656717fe9dca37c5a3e58bf196c620a8bf..6d0ad2f37c8b0183598cebed10fce5a37cd5d749 100644
--- a/Processes/Proposal/Interaction.h
+++ b/Processes/Proposal/Interaction.h
@@ -21,6 +21,7 @@
#include "PROPOSAL/PROPOSAL.h"
using namespace corsika::environment;
+using namespace corsika::units::si;
using CORSIKA_ParticleCut = corsika::process::particle_cut::ParticleCut;
using std::make_pair;
@@ -29,14 +30,11 @@ using std::make_tuple;
namespace corsika::process::proposal {
class Interaction : public corsika::process::InteractionProcess<Interaction> {
-
- shared_ptr<const PROPOSAL::EnergyCutSettings> cut;
-
+ CORSIKA_ParticleCut& cut;
+ corsika::random::RNG& fRNG;
static std::unordered_map<particles::Code, PROPOSAL::ParticleDef> particles;
std::unordered_map<const NuclearComposition*, PROPOSAL::Medium> media;
- corsika::random::RNG& fRNG;
-
bool CanInteract(particles::Code pcode) const noexcept;
using calculator_t = tuple<unique_ptr<PROPOSAL::SecondariesCalculator>,
@@ -53,12 +51,15 @@ namespace corsika::process::proposal {
std::unordered_map<calc_key_t, calculator_t, interaction_hash> calculators;
template <typename Particle>
- auto BuildCalculator(particles::Code corsika_code, Particle p_def,
+ auto BuildCalculator(particles::Code code, Particle p_def,
NuclearComposition const& comp) {
- auto cross = GetStdCrossSections(p_def, media.at(&comp), cut, true);
+ auto cross = GetStdCrossSections(
+ p_def, media.at(&comp),
+ make_shared<const PROPOSAL::EnergyCutSettings>(cut.GetECut() / 1_MeV, 1, false),
+ true);
auto inter_types = PROPOSAL::CrossSectionVector::GetInteractionTypes(cross);
auto [insert_it, success] = calculators.insert(
- {make_pair(&comp, corsika_code),
+ {make_pair(&comp, code),
make_tuple(PROPOSAL::make_secondaries(inter_types, p_def, media.at(&comp)),
PROPOSAL::make_interaction(cross, true))});
return insert_it;
@@ -93,15 +94,15 @@ namespace corsika::process::proposal {
public:
template <typename TEnvironment>
- Interaction(TEnvironment const& env, CORSIKA_ParticleCut const& cut);
+ Interaction(TEnvironment const& env, CORSIKA_ParticleCut& cut);
- void Init() {};
+ void Init(){};
template <typename Particle>
corsika::process::EProcessReturn DoInteraction(Particle&);
template <typename TParticle>
corsika::units::si::GrammageType GetInteractionLength(TParticle const& p);
- }; // namespace corsika::process::proposal
+ };
} // namespace corsika::process::proposal
#endif