From 17b8be7cdaa795248ff615b1f8fbae777b75e30b Mon Sep 17 00:00:00 2001
From: Jean-Marco Alameddine <jean-marco.alameddine@tu-dortmund.de>
Date: Wed, 1 Jul 2020 15:19:38 +0200
Subject: [PATCH] Add mapping of CORSIKA particles to PROPOSAL calculators
---
Documentation/Examples/proposal_example.cc | 6 ++---
Processes/Proposal/Interaction.cc | 4 ++--
Processes/Proposal/Interaction.h | 28 +++++++++++++++-------
3 files changed, 24 insertions(+), 14 deletions(-)
diff --git a/Documentation/Examples/proposal_example.cc b/Documentation/Examples/proposal_example.cc
index 44ccf38e9..3fb41d169 100644
--- a/Documentation/Examples/proposal_example.cc
+++ b/Documentation/Examples/proposal_example.cc
@@ -27,7 +27,7 @@
#include <corsika/setup/SetupTrajectory.h>
#include <corsika/units/PhysicalUnits.h>
#include <corsika/utl/CorsikaFenv.h>
-
+#include <corsika/process/track_writer/TrackWriter.h>
#include <corsika/process/proposal/Interaction.h>
#include <iomanip>
@@ -138,6 +138,7 @@ int main(int argc, char** argv) {
process::particle_cut::ParticleCut cut(10_GeV);
process::proposal::Interaction proposal(env, cut);
process::interaction_counter::InteractionCounter proposalCounted(proposal);
+ process::track_writer::TrackWriter trackWriter("tracks.dat");
// energy cut; n.b. ParticleCut needs to be modified not to discard EM particles!
/* process::particle_cut::ParticleCut cut{60_GeV}; */
@@ -149,8 +150,7 @@ int main(int argc, char** argv) {
process::observation_plane::ObservationPlane observationLevel(obsPlane,
"particles.dat");
- auto sequence = proposalCounted << longprof << proposal << cut << observationLevel;
-
+ auto sequence = proposalCounted << longprof << proposal << cut << observationLevel << trackWriter;
// define air shower object, run simulation
tracking_line::TrackingLine tracking;
cascade::Cascade EAS(env, tracking, sequence, stack);
diff --git a/Processes/Proposal/Interaction.cc b/Processes/Proposal/Interaction.cc
index 0e6234cf6..65d53c557 100644
--- a/Processes/Proposal/Interaction.cc
+++ b/Processes/Proposal/Interaction.cc
@@ -76,7 +76,7 @@ namespace corsika::process::proposal {
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)
+ i < std::get<SECONDARIES>(calc->second)->RequiredRandomNumbers(type); ++i)
rnd.push_back(distr(fRNG));
double primary_energy = vP.GetEnergy() / 1_MeV;
auto point = PROPOSAL::Vector3D(vP.GetPosition().GetX() / 1_cm,
@@ -89,7 +89,7 @@ namespace corsika::process::proposal {
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);
+ ->CalculateSecondaries(primary_energy, loss, *comp_ptr, rnd);
for (auto& s : sec) {
auto energy = get<PROPOSAL::Loss::ENERGY>(s) * 1_MeV;
auto vec = corsika::geometry::QuantityVector(
diff --git a/Processes/Proposal/Interaction.h b/Processes/Proposal/Interaction.h
index be6835069..48f3fbeb1 100644
--- a/Processes/Proposal/Interaction.h
+++ b/Processes/Proposal/Interaction.h
@@ -39,9 +39,19 @@ namespace corsika::process::proposal {
bool CanInteract(particles::Code pcode) const noexcept;
using calculator_t =
- tuple<PROPOSAL::SecondariesCalculator, unique_ptr<PROPOSAL::Interaction>,
+ tuple<unique_ptr<PROPOSAL::SecondariesCalculator>, unique_ptr<PROPOSAL::Interaction>,
unique_ptr<PROPOSAL::Displacement>>;
- std::unordered_map<const NuclearComposition*, calculator_t> calculators;
+
+ 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;
+ }
+ };
+
+ std::unordered_map<std::pair<const NuclearComposition*, particles::Code>, calculator_t, interaction_hash> calculators;
auto BuildCalculator(particles::Code corsika_code, NuclearComposition const& comp) {
@@ -52,8 +62,8 @@ namespace corsika::process::proposal {
GetStdCrossSections(PROPOSAL::GammaDef(), media.at(&comp), cut, true);
auto inter_types = PROPOSAL::CrossSectionVector::GetInteractionTypes(cross);
auto [insert_it, success] = calculators.insert(
- {&comp, make_tuple(PROPOSAL::SecondariesCalculator(
- inter_types, PROPOSAL::GammaDef(), media[&comp]),
+ {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;
@@ -64,8 +74,8 @@ namespace corsika::process::proposal {
GetStdCrossSections(PROPOSAL::EMinusDef(), media.at(&comp), cut, true);
auto inter_types = PROPOSAL::CrossSectionVector::GetInteractionTypes(cross);
auto [insert_it, success] = calculators.insert(
- {&comp, make_tuple(PROPOSAL::SecondariesCalculator(
- inter_types, PROPOSAL::EMinusDef(), media[&comp]),
+ {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;
@@ -76,8 +86,8 @@ namespace corsika::process::proposal {
GetStdCrossSections(PROPOSAL::EPlusDef(), media.at(&comp), cut, true);
auto inter_types = PROPOSAL::CrossSectionVector::GetInteractionTypes(cross);
auto [insert_it, success] = calculators.insert(
- {&comp, make_tuple(PROPOSAL::SecondariesCalculator(
- inter_types, PROPOSAL::EPlusDef(), media[&comp]),
+ {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;
@@ -88,7 +98,7 @@ namespace corsika::process::proposal {
template <typename Particle>
auto GetCalculator(Particle& vP) {
auto& comp = vP.GetNode()->GetModelProperties().GetNuclearComposition();
- auto calc_it = calculators.find(&comp);
+ auto calc_it = calculators.find(std::make_pair(&comp, vP.GetPID()));
if (calc_it != calculators.end()) return calc_it;
return BuildCalculator(vP.GetPID(), comp);
}
--
GitLab