From 4516ac8ce57c4e0dbfa3b3c47055495fe80a03c8 Mon Sep 17 00:00:00 2001
From: Jean-Marco Alameddine <jean-marco.alameddine@tu-dortmund.de>
Date: Thu, 2 Jul 2020 16:43:24 +0200
Subject: [PATCH] Include particle mapping for ContinuousProcess as well as
some unit fixing
---
Processes/Proposal/ContinuousProcess.cc | 20 ++++++-----
Processes/Proposal/ContinuousProcess.h | 45 ++++++++++++++++++++-----
2 files changed, 48 insertions(+), 17 deletions(-)
diff --git a/Processes/Proposal/ContinuousProcess.cc b/Processes/Proposal/ContinuousProcess.cc
index 523c9400a..5914d5296 100644
--- a/Processes/Proposal/ContinuousProcess.cc
+++ b/Processes/Proposal/ContinuousProcess.cc
@@ -36,22 +36,22 @@ namespace corsika::process::proposal {
template <>
ContinuousProcess::ContinuousProcess(SetupEnvironment const& _env,
CORSIKA_ParticleCut const& _cut)
- : cut(make_shared<const PROPOSAL::EnergyCutSettings>(_cut.GetCutEnergy() / 1_GeV, 1,
+ : cut(make_shared<const PROPOSAL::EnergyCutSettings>(_cut.GetECut() / 1_MeV, 1,
false)) {
- auto all_compositions = std::vector<NuclearComposition>();
+ auto all_compositions = std::vector<const NuclearComposition*>();
_env.GetUniverse()->walk([&](auto& vtn) {
if (vtn.HasModelProperties())
- all_compositions.push_back(vtn.GetModelProperties().GetNuclearComposition());
+ all_compositions.push_back(&vtn.GetModelProperties().GetNuclearComposition());
});
for (auto& ncarg : all_compositions) {
auto comp_vec = std::vector<PROPOSAL::Components::Component>();
- auto frac_iter = ncarg.GetFractions().cbegin();
- for (auto& pcode : ncarg.GetComponents()) {
+ auto frac_iter = ncarg->GetFractions().cbegin();
+ for (auto& pcode : ncarg->GetComponents()) {
comp_vec.emplace_back(GetName(pcode), GetNucleusZ(pcode), GetNucleusA(pcode),
*frac_iter);
++frac_iter;
}
- media[&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);
@@ -62,14 +62,16 @@ namespace corsika::process::proposal {
GrammageType const vDX) {
auto calc_ptr = GetCalculator(vP);
auto upper_energy = calc_ptr->second->UpperLimitTrackIntegral(
- vP.GetEnergy() / 1_GeV, vDX / 1_g * 1_cm * 1_cm);
- return upper_energy * 1_GeV;
+ 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;
}
+ void ContinuousProcess::Init() {}
template <>
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);
diff --git a/Processes/Proposal/ContinuousProcess.h b/Processes/Proposal/ContinuousProcess.h
index ef0782c70..e123ee938 100644
--- a/Processes/Proposal/ContinuousProcess.h
+++ b/Processes/Proposal/ContinuousProcess.h
@@ -42,18 +42,45 @@ namespace corsika::process::proposal {
bool CanInteract(particles::Code pcode) const noexcept;
- unordered_map<const NuclearComposition*, unique_ptr<PROPOSAL::Displacement>> calc;
-
+ 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;
+ }
+ };
+
+ unordered_map<std::pair<const NuclearComposition*, particles::Code>, unique_ptr<PROPOSAL::Displacement>, interaction_hash> calc;
+
+ 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;
+ }
+ }
+
template <typename Particle>
auto GetCalculator(Particle& vP) {
auto& comp = vP.GetNode()->GetModelProperties().GetNuclearComposition();
- auto calc_it = calc.find(&comp);
+ auto calc_it = calc.find(std::make_pair(&comp, vP.GetPID()));
if (calc_it != calc.end()) return calc_it;
- auto cross =
- PROPOSAL::GetStdCrossSections(particles[vP.GetPID()], media[&comp], cut, true);
- auto [insert_it, success] =
- calc.insert({&comp, PROPOSAL::make_displacement(cross, true)});
- return insert_it;
+ return BuildCalculator(vP.GetPID(), comp);
}
units::si::HEPEnergyType TotalEnergyLoss(setup::Stack::ParticleType const&,
@@ -63,6 +90,8 @@ namespace corsika::process::proposal {
template <typename TEnvironment>
ContinuousProcess(TEnvironment const& env, CORSIKA_ParticleCut const& cut);
+ void Init();
+
template <typename Particle, typename Track>
EProcessReturn DoContinuous(Particle&, Track const&) ;
--
GitLab