diff --git a/Processes/Proposal/Interaction.cc b/Processes/Proposal/Interaction.cc
index 65d53c557a77e5b8b5204252fb1c4e5dcf4619db..4394ae868f89dfce4871cc8c077f6eb2fb00cb34 100644
--- a/Processes/Proposal/Interaction.cc
+++ b/Processes/Proposal/Interaction.cc
@@ -30,17 +30,15 @@ namespace corsika::process::proposal {
};
bool Interaction::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 <>
Interaction::Interaction(SetupEnvironment const& _env, CORSIKA_ParticleCut const& _cut)
: cut(make_shared<const PROPOSAL::EnergyCutSettings>(_cut.GetECut() / 1_MeV, 1,
- false)) {
- std::cout << "corsika set cut to " << _cut.GetECut() / 1_GeV << " [GeV]" << std::endl;
- std::cout << "proposal set ecut to " << cut->GetEcut() << " [MeV]" << std::endl;
+ false))
+ , fRNG(corsika::random::RNGManager::GetInstance().GetRandomStream("proposal")) {
auto all_compositions = std::vector<const NuclearComposition*>();
_env.GetUniverse()->walk([&](auto& vtn) {
if (vtn.HasModelProperties()) {
@@ -62,8 +60,6 @@ namespace corsika::process::proposal {
}
}
- void Interaction::Init() {}
-
template <>
corsika::process::EProcessReturn Interaction::DoInteraction(
setup::StackView::StackIterator& vP) {
@@ -71,39 +67,34 @@ namespace corsika::process::proposal {
auto calc = GetCalculator(vP); // [CrossSections]
std::uniform_real_distribution<double> distr(0., 1.);
auto [type, comp_ptr, v] =
- std::get<INTERACTION>(calc->second)
+ get<INTERACTION>(calc->second)
->TypeInteraction(vP.GetEnergy() / 1_MeV, distr(fRNG));
- std::cout << "InteractionType: "<< static_cast<int>(type) << std::endl;
- auto rnd = std::vector<double>();
- for (size_t i = 0;
- i < std::get<SECONDARIES>(calc->second)->RequiredRandomNumbers(type); ++i)
- rnd.push_back(distr(fRNG));
- double primary_energy = vP.GetEnergy() / 1_MeV;
+ auto rnd =
+ vector<double>(get<SECONDARIES>(calc->second)->RequiredRandomNumbers(type));
+ for (auto& it : rnd) it = distr(fRNG);
auto point = PROPOSAL::Vector3D(vP.GetPosition().GetX() / 1_cm,
vP.GetPosition().GetY() / 1_cm,
vP.GetPosition().GetZ() / 1_cm);
auto d = vP.GetDirection().GetComponents();
- auto direction = PROPOSAL::Vector3D(d.GetX().magnitude(),
- d.GetY().magnitude(),
+ auto direction = PROPOSAL::Vector3D(d.GetX().magnitude(), d.GetY().magnitude(),
d.GetZ().magnitude());
- auto loss =
- make_tuple(static_cast<int>(type), point, direction, v * primary_energy, 0.);
- auto sec = std::get<SECONDARIES>(calc->second)
- ->CalculateSecondaries(primary_energy, loss, *comp_ptr, rnd);
+ auto loss = make_tuple(static_cast<int>(type), point, direction,
+ v * vP.GetEnergy() / 1_MeV, 0.);
+ auto sec = get<SECONDARIES>(calc->second)
+ ->CalculateSecondaries(vP.GetEnergy() / 1_MeV, loss, *comp_ptr, rnd);
for (auto& s : sec) {
- auto energy = get<PROPOSAL::Loss::ENERGY>(s) * 1_MeV;
+ auto E = get<PROPOSAL::Loss::ENERGY>(s) * 1_MeV;
auto vec = corsika::geometry::QuantityVector(
- std::get<PROPOSAL::Loss::DIRECTION>(s).GetX() * energy,
- std::get<PROPOSAL::Loss::DIRECTION>(s).GetY() * energy,
- std::get<PROPOSAL::Loss::DIRECTION>(s).GetZ() * energy);
- auto momentum = corsika::stack::MomentumVector(
+ get<PROPOSAL::Loss::DIRECTION>(s).GetX() * E,
+ get<PROPOSAL::Loss::DIRECTION>(s).GetY() * E,
+ get<PROPOSAL::Loss::DIRECTION>(s).GetZ() * E);
+ auto p = corsika::stack::MomentumVector(
corsika::geometry::RootCoordinateSystem::GetInstance()
.GetRootCoordinateSystem(),
vec);
- particles::Code sec_code = corsika::particles::ConvertFromPDG(
+ auto sec_code = corsika::particles::ConvertFromPDG(
static_cast<particles::PDGCode>(get<PROPOSAL::Loss::TYPE>(s)));
- vP.AddSecondary(
- make_tuple(sec_code, energy, momentum, vP.GetPosition(), vP.GetTime()));
+ vP.AddSecondary(make_tuple(sec_code, E, p, vP.GetPosition(), vP.GetTime()));
}
}
return process::EProcessReturn::eOk;
@@ -113,14 +104,9 @@ namespace corsika::process::proposal {
corsika::units::si::GrammageType Interaction::GetInteractionLength(
setup::Stack::StackIterator const& vP) {
if (CanInteract(vP.GetPID())) {
- auto calc = GetCalculator(vP); // [CrossSections]
- std::uniform_real_distribution<double> distr(0., 1.);
- auto rnd = distr(fRNG);
- auto energy = get<INTERACTION>(calc->second)
- ->EnergyInteraction(vP.GetEnergy() / 1_MeV, rnd);
- return get<DISPLACEMENT>(calc->second)
- ->SolveTrackIntegral(vP.GetEnergy() / 1_MeV, energy) *
- 1_g / 1_cm / 1_cm;
+ auto calc = GetCalculator(vP);
+ return get<INTERACTION>(calc->second)->MeanFreePath(vP.GetEnergy() / 1_MeV) * 1_g /
+ (1_cm * 1_cm);
}
return std::numeric_limits<double>::infinity() * 1_g / (1_cm * 1_cm);
}
diff --git a/Processes/Proposal/Interaction.h b/Processes/Proposal/Interaction.h
index b95b6afdbe18133bf49f66c88adf69378d3d5fba..54597a656717fe9dca37c5a3e58bf196c620a8bf 100644
--- a/Processes/Proposal/Interaction.h
+++ b/Processes/Proposal/Interaction.h
@@ -23,6 +23,8 @@
using namespace corsika::environment;
using CORSIKA_ParticleCut = corsika::process::particle_cut::ParticleCut;
+using std::make_pair;
+using std::make_tuple;
namespace corsika::process::proposal {
@@ -33,130 +35,67 @@ namespace corsika::process::proposal {
static std::unordered_map<particles::Code, PROPOSAL::ParticleDef> particles;
std::unordered_map<const NuclearComposition*, PROPOSAL::Medium> media;
- corsika::random::RNG& fRNG =
- corsika::random::RNGManager::GetInstance().GetRandomStream("proposal");
+ corsika::random::RNG& fRNG;
bool CanInteract(particles::Code pcode) const noexcept;
- using calculator_t =
- tuple<unique_ptr<PROPOSAL::SecondariesCalculator>, unique_ptr<PROPOSAL::Interaction>,
- unique_ptr<PROPOSAL::Displacement>>;
+ using calculator_t = tuple<unique_ptr<PROPOSAL::SecondariesCalculator>,
+ unique_ptr<PROPOSAL::Interaction>>;
+ using calc_key_t = std::pair<const NuclearComposition*, particles::Code>;
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;
- }
+ size_t operator()(const calc_key_t& p) const {
+ return std::hash<const NuclearComposition*>{}(p.first) ^
+ std::hash<particles::Code>{}(p.second);
+ }
};
- std::unordered_map<std::pair<const NuclearComposition*, particles::Code>, calculator_t, interaction_hash> calculators;
+ std::unordered_map<calc_key_t, calculator_t, interaction_hash> calculators;
+ template <typename Particle>
+ auto BuildCalculator(particles::Code corsika_code, Particle p_def,
+ NuclearComposition const& comp) {
+ auto cross = GetStdCrossSections(p_def, media.at(&comp), cut, true);
+ auto inter_types = PROPOSAL::CrossSectionVector::GetInteractionTypes(cross);
+ auto [insert_it, success] = calculators.insert(
+ {make_pair(&comp, corsika_code),
+ make_tuple(PROPOSAL::make_secondaries(inter_types, p_def, media.at(&comp)),
+ PROPOSAL::make_interaction(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) {
- std::cout << "Build gamma tables" << std::endl;
- auto cross =
- GetStdCrossSections(PROPOSAL::GammaDef(), media.at(&comp), cut, true);
- auto inter_types = PROPOSAL::CrossSectionVector::GetInteractionTypes(cross);
- auto [insert_it, success] = calculators.insert(
- {std::make_pair(&comp, corsika_code), make_tuple(PROPOSAL::make_secondaries(
- inter_types, PROPOSAL::GammaDef(), media.at(&comp)),
- PROPOSAL::make_interaction(cross, true),
- PROPOSAL::make_displacement(cross, true))});
- return insert_it;
- }
- if (corsika_code == particles::Code::Electron) {
- std::cout << "Build electron tables" << std::endl;
- auto cross =
- GetStdCrossSections(PROPOSAL::EMinusDef(), media.at(&comp), cut, true);
- auto inter_types = PROPOSAL::CrossSectionVector::GetInteractionTypes(cross);
- auto [insert_it, success] = calculators.insert(
- {std::make_pair(&comp, corsika_code), make_tuple(PROPOSAL::make_secondaries(
- inter_types, PROPOSAL::EMinusDef(), media.at(&comp)),
- PROPOSAL::make_interaction(cross, true),
- PROPOSAL::make_displacement(cross, true))});
- return insert_it;
- }
- if (corsika_code == particles::Code::Positron) {
- std::cout << "Build positron tables" << std::endl;
- auto cross =
- GetStdCrossSections(PROPOSAL::EPlusDef(), media.at(&comp), cut, true);
- auto inter_types = PROPOSAL::CrossSectionVector::GetInteractionTypes(cross);
- auto [insert_it, success] = calculators.insert(
- {std::make_pair(&comp, corsika_code), make_tuple(PROPOSAL::make_secondaries(
- inter_types, PROPOSAL::EPlusDef(), media.at(&comp)),
- PROPOSAL::make_interaction(cross, true),
- PROPOSAL::make_displacement(cross, true))});
- return insert_it;
- }
- if (corsika_code == particles::Code::MuMinus) {
- std::cout << "Build MuMinus tables" << std::endl;
- auto cross =
- GetStdCrossSections(PROPOSAL::MuMinusDef(), media.at(&comp), cut, true);
- auto inter_types = PROPOSAL::CrossSectionVector::GetInteractionTypes(cross);
- auto [insert_it, success] = calculators.insert(
- {std::make_pair(&comp, corsika_code), make_tuple(PROPOSAL::make_secondaries(
- inter_types, PROPOSAL::MuMinusDef(), media.at(&comp)),
- PROPOSAL::make_interaction(cross, true),
- PROPOSAL::make_displacement(cross, true))});
- return insert_it;
- }
- if (corsika_code == particles::Code::MuPlus) {
- std::cout << "Build MuPlus tables" << std::endl;
- auto cross =
- GetStdCrossSections(PROPOSAL::MuPlusDef(), media.at(&comp), cut, true);
- auto inter_types = PROPOSAL::CrossSectionVector::GetInteractionTypes(cross);
- auto [insert_it, success] = calculators.insert(
- {std::make_pair(&comp, corsika_code), make_tuple(PROPOSAL::make_secondaries(
- inter_types, PROPOSAL::MuPlusDef(), media.at(&comp)),
- PROPOSAL::make_interaction(cross, true),
- PROPOSAL::make_displacement(cross, true))});
- return insert_it;
- }
- if (corsika_code == particles::Code::TauMinus) {
- std::cout << "Build TauMinus tables" << std::endl;
- auto cross =
- GetStdCrossSections(PROPOSAL::TauMinusDef(), media.at(&comp), cut, true);
- auto inter_types = PROPOSAL::CrossSectionVector::GetInteractionTypes(cross);
- auto [insert_it, success] = calculators.insert(
- {std::make_pair(&comp, corsika_code), make_tuple(PROPOSAL::make_secondaries(
- inter_types, PROPOSAL::TauMinusDef(), media.at(&comp)),
- PROPOSAL::make_interaction(cross, true),
- PROPOSAL::make_displacement(cross, true))});
- return insert_it;
- }
- if (corsika_code == particles::Code::TauPlus) {
- std::cout << "Build TauPlus tables" << std::endl;
- auto cross =
- GetStdCrossSections(PROPOSAL::TauPlusDef(), media.at(&comp), cut, true);
- auto inter_types = PROPOSAL::CrossSectionVector::GetInteractionTypes(cross);
- auto [insert_it, success] = calculators.insert(
- {std::make_pair(&comp, corsika_code), make_tuple(PROPOSAL::make_secondaries(
- inter_types, PROPOSAL::TauPlusDef(), media.at(&comp)),
- PROPOSAL::make_interaction(cross, true),
- PROPOSAL::make_displacement(cross, true))});
- return insert_it;
- }
+ 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");
- } // namespace corsika::process::proposal
+ }
- enum { SECONDARIES, INTERACTION, DISPLACEMENT };
+ enum { SECONDARIES, INTERACTION };
template <typename Particle>
auto GetCalculator(Particle& vP) {
- auto& comp = vP.GetNode()->GetModelProperties().GetNuclearComposition();
- auto calc_it = calculators.find(std::make_pair(&comp, vP.GetPID()));
- if (calc_it != calculators.end()) return calc_it;
- return BuildCalculator(vP.GetPID(), comp);
+ auto& comp = vP.GetNode()->GetModelProperties().GetNuclearComposition();
+ auto calc_it = calculators.find(make_pair(&comp, vP.GetPID()));
+ if (calc_it != calculators.end()) return calc_it;
+ return BuildCalculator(vP.GetPID(), comp);
}
public:
template <typename TEnvironment>
Interaction(TEnvironment const& env, CORSIKA_ParticleCut const& cut);
- void Init();
+ void Init() {};
template <typename Particle>
corsika::process::EProcessReturn DoInteraction(Particle&);