From 9397c1a97bd5d77d93f82027932a505ad59d6f44 Mon Sep 17 00:00:00 2001
From: Maximilian Sackel <maximilian.sackel@tu-dortmund.de>
Date: Sun, 21 Jun 2020 10:30:48 +0000
Subject: [PATCH] add proposal sec to corsika stack
---
Processes/Proposal/Interaction.cc | 49 ++++++++++---------------------
Processes/Proposal/Interaction.h | 44 +++++++++++++--------------
2 files changed, 36 insertions(+), 57 deletions(-)
diff --git a/Processes/Proposal/Interaction.cc b/Processes/Proposal/Interaction.cc
index a62e45ca9..a19b6bf5f 100644
--- a/Processes/Proposal/Interaction.cc
+++ b/Processes/Proposal/Interaction.cc
@@ -19,13 +19,11 @@ namespace corsika::process::proposal {
using namespace corsika::setup;
using namespace corsika::environment;
using namespace corsika::units::si;
- using namespace corsika;
-
using SetupParticle = corsika::setup::Stack::StackIterator;
template <>
std::unordered_map<particles::Code, PROPOSAL::ParticleDef>
- Interaction<SetupEnvironment>::particle_map{
+ Interaction<SetupEnvironment>::particles{
{particles::Code::Gamma, PROPOSAL::GammaDef()},
{particles::Code::Electron, PROPOSAL::EMinusDef()},
{particles::Code::Positron, PROPOSAL::EPlusDef()},
@@ -40,8 +38,10 @@ namespace corsika::process::proposal {
CORSIKA_ParticleCut const& e_cut)
: fEnvironment(env)
, cut(make_shared<const PROPOSAL::EnergyCutSettings>(e_cut.GetCutEnergy() / 1_GeV,
- 0, false)) {
+ 1, false)) {}
+ template <>
+ void Interaction<SetupEnvironment>::Init() {
auto all_compositions = std::vector<NuclearComposition>();
fEnvironment.GetUniverse()->walk([&](auto& vtn) {
if (vtn.HasModelProperties())
@@ -55,44 +55,24 @@ namespace corsika::process::proposal {
*frac_iter);
++frac_iter;
}
- medium_map[&ncarg] = PROPOSAL::Medium(
+ 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);
}
- std::cout << "PROPOSAL done." << std::endl;
- }
-
- template <>
- template <>
- auto Interaction<SetupEnvironment>::GetCalculator(SetupParticle const& vP) {
- auto& mediumComposition = vP.GetNode()->GetModelProperties().GetNuclearComposition();
- auto calc_it = calculators.find(&mediumComposition);
- if (calc_it != calculators.end()) return calc_it;
- auto const& p_def = particle_map[vP.GetPID()];
- auto& medium = medium_map[&mediumComposition];
- auto cross = PROPOSAL::GetStdCrossSections(p_def, medium, cut, true);
- auto inter_types = PROPOSAL::CrossSectionVector::GetInteractionTypes(cross);
- auto [insert_it, success] = calculators.insert(
- {&mediumComposition,
- make_tuple(PROPOSAL::SecondariesCalculator(inter_types, p_def, medium),
- PROPOSAL::make_interaction(cross, true),
- PROPOSAL::make_displacement(cross, true))});
- if (success) return insert_it;
- throw std::logic_error("insert failed.");
}
template <>
template <>
corsika::process::EProcessReturn Interaction<SetupEnvironment>::DoInteraction(
- SetupParticle const& vP) {
+ SetupParticle& vP) {
auto calc = GetCalculator(vP); // [CrossSections]
std::uniform_real_distribution<double> distr(0., 1.);
auto [type, comp_ptr, v] = std::get<INTERACTION>(calc->second)
->TypeInteraction(vP.GetEnergy() / 1_GeV, distr(fRNG));
auto rnd = std::vector<double>();
- for (auto i = 0; i < std::get<SECONDARIES>(calc->second).RequiredRandomNumbers(type);
- ++i)
+ 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_GeV;
auto point =
@@ -114,11 +94,14 @@ namespace corsika::process::proposal {
corsika::geometry::RootCoordinateSystem::GetInstance()
.GetRootCoordinateSystem(),
vec);
- // particles::Code, units::si::HEPEnergyType, stack::MomentumVector,
- // geometry::Point, units::si::TimeType
- auto sec = make_tuple(get<PROPOSAL::Loss::TYPE>(s), energy, momentum,
- vP.GetPosition(), vP.GetTime());
- vP.AddSecondary(sec);
+ particles::Code sec_code = corsika::particles::ConvertFromPDG(
+ static_cast<particles::PDGCode>(get<PROPOSAL::Loss::TYPE>(s)));
+ corsika::units::si::HEPEnergyType sec_energy = energy;
+ stack::MomentumVector sec_momentum = momentum;
+ geometry::Point sec_position = vP.GetPosition();
+ corsika::units::si::TimeType sec_time = vP.GetTime();
+ vP.AddSecondary(
+ make_tuple(sec_code, sec_energy, sec_momentum, sec_position, sec_time));
}
return process::EProcessReturn::eOk;
}
diff --git a/Processes/Proposal/Interaction.h b/Processes/Proposal/Interaction.h
index 6dfe30995..8692dd4e6 100644
--- a/Processes/Proposal/Interaction.h
+++ b/Processes/Proposal/Interaction.h
@@ -27,16 +27,6 @@ using CORSIKA_ParticleCut = corsika::process::particle_cut::ParticleCut;
namespace corsika::process::proposal {
- /* static std::unordered_map<particles::Code, PROPOSAL::ParticleDef> particle_map{ */
- /* {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 <class TEnvironment>
class Interaction
: public corsika::process::InteractionProcess<Interaction<TEnvironment>> {
@@ -45,26 +35,17 @@ namespace corsika::process::proposal {
TEnvironment const& fEnvironment;
shared_ptr<const PROPOSAL::EnergyCutSettings> cut;
- static std::unordered_map<particles::Code, PROPOSAL::ParticleDef> particle_map;
- std::unordered_map<const NuclearComposition*, PROPOSAL::Medium> medium_map;
+ static std::unordered_map<particles::Code, PROPOSAL::ParticleDef> particles;
+ std::unordered_map<const NuclearComposition*, PROPOSAL::Medium> media;
enum { SECONDARIES, INTERACTION, DISPLACEMENT };
- /* std::uniform_real_distribution<> rnd_uniform(0., 1.); */
-
- /* std::map<particles::Code, std::unique_ptr<PROPOSAL::UtilityInterpolantInteraction>>
- * corsika_particle_to_utility_map; */
- /* std::map<particles::Code,
- * std::unique_ptr<PROPOSAL::UtilityInterpolantDisplacement>>
- * corsika_particle_to_displacement_map; */
corsika::random::RNG& fRNG =
corsika::random::RNGManager::GetInstance().GetRandomStream("s_rndm");
- // fTrackedParticles[proposal_particle.GetName()] return particle::Code
-
auto IsTracked(particles::Code pcode) const noexcept {
- auto search = particle_map.find(pcode);
- if (search != particle_map.end()) return true;
+ auto search = particles.find(pcode);
+ if (search != particles.end()) return true;
return false;
};
@@ -74,11 +55,26 @@ namespace corsika::process::proposal {
std::unordered_map<const NuclearComposition*, calculator_t> calculators;
template <typename Particle>
- auto GetCalculator(Particle&);
+ auto GetCalculator(Particle& vP) {
+ auto& comp = vP.GetNode()->GetModelProperties().GetNuclearComposition();
+ auto calc_it = calculators.find(&comp);
+ if (calc_it != calculators.end()) return calc_it;
+ auto cross =
+ PROPOSAL::GetStdCrossSections(particles[vP.GetPID()], media[&comp], cut, true);
+ auto inter_types = PROPOSAL::CrossSectionVector::GetInteractionTypes(cross);
+ auto [insert_it, success] = calculators.insert(
+ {&comp, make_tuple(PROPOSAL::SecondariesCalculator(
+ inter_types, particles[vP.GetPID()], media[&comp]),
+ PROPOSAL::make_interaction(cross, true),
+ PROPOSAL::make_displacement(cross, true))});
+ return insert_it;
+ }
public:
Interaction(TEnvironment const& env, CORSIKA_ParticleCut const& cut);
+ void Init();
+
template <typename Particle>
corsika::process::EProcessReturn DoInteraction(Particle&);
--
GitLab