From bb485605f1dfa156c6523e61294c0b4c8ce423fd Mon Sep 17 00:00:00 2001
From: Maximilian Sackel <maximilian.sackel@tu-dortmund.de>
Date: Fri, 2 Oct 2020 11:47:53 +0000
Subject: [PATCH] Consider the types of the different coordinate systems when
calculating the momentum. Not using scoped enumeration because static_cast is
required to specify the tuple position. Instead mark variable with an extra
letter 'e' and use it only very scoped.
---
Processes/Proposal/ContinuousProcess.cc | 18 +++++++++---------
Processes/Proposal/ContinuousProcess.h | 2 +-
Processes/Proposal/Interaction.cc | 19 ++++++++-----------
Processes/Proposal/Interaction.h | 2 +-
4 files changed, 19 insertions(+), 22 deletions(-)
diff --git a/Processes/Proposal/ContinuousProcess.cc b/Processes/Proposal/ContinuousProcess.cc
index d3fcd605e..7ceb75b47 100644
--- a/Processes/Proposal/ContinuousProcess.cc
+++ b/Processes/Proposal/ContinuousProcess.cc
@@ -56,7 +56,8 @@ namespace corsika::process::proposal {
auto c = GetCalculator(vP, calc);
// Cast corsika vector to proposal vector
- auto d = vP.GetDirection().GetComponents();
+ auto vP_dir = vP.GetDirection();
+ auto d = vP_dir.GetComponents();
auto direction = PROPOSAL::Vector3D(d.GetX().magnitude(), d.GetY().magnitude(),
d.GetZ().magnitude());
@@ -68,7 +69,7 @@ namespace corsika::process::proposal {
for (auto& it : rnd) it = distr(fRNG);
// calculate deflection based on particle energy, loss
- auto [mean_dir, final_dir] = get<SCATTERING>(c->second)->Scatter(
+ auto [mean_dir, final_dir] = get<eSCATTERING>(c->second)->Scatter(
grammage / 1_g * square(1_cm), vP.GetEnergy() / 1_MeV, E_f / 1_MeV, direction,
rnd);
@@ -79,9 +80,8 @@ namespace corsika::process::proposal {
// scattering
auto vec = corsika::geometry::QuantityVector(
final_dir.GetX() * E_f, final_dir.GetY() * E_f, final_dir.GetZ() * E_f);
- vP.SetMomentum(corsika::stack::MomentumVector(
- corsika::geometry::RootCoordinateSystem::GetInstance().GetRootCoordinateSystem(),
- vec));
+ vP.SetMomentum(
+ corsika::stack::MomentumVector(vP_dir.GetCoordinateSystem(), vec));
}
template <>
@@ -90,7 +90,7 @@ namespace corsika::process::proposal {
using namespace corsika::units::si; // required for operator::_MeV
if (!CanInteract(vP.GetPID())) return process::EProcessReturn::eOk;
- if (vT.GetLength()==0_m) return process::EProcessReturn::eOk;
+ if (vT.GetLength() == 0_m) return process::EProcessReturn::eOk;
// calculate passed grammage
auto dX = vP.GetNode()->GetModelProperties().IntegratedGrammage(vT, vT.GetLength());
@@ -98,7 +98,7 @@ namespace corsika::process::proposal {
// Get or build corresponding track integral calculator and solve the
// integral
auto c = GetCalculator(vP, calc);
- auto final_energy = get<DISPLACEMENT>(c->second)->UpperLimitTrackIntegral(
+ auto final_energy = get<eDISPLACEMENT>(c->second)->UpperLimitTrackIntegral(
vP.GetEnergy() / 1_MeV, dX / 1_g * 1_cm * 1_cm) *
1_MeV;
auto dE = vP.GetEnergy() - final_energy;
@@ -144,11 +144,11 @@ namespace corsika::process::proposal {
// important, the important fact is that its E_kin is zero
// afterwards.
//
- auto energy_lim = std::max(0.9 * vP.GetEnergy(), 0.99*emCut_);
+ auto energy_lim = std::max(0.9 * vP.GetEnergy(), 0.99 * emCut_);
// solving the track integral for giving energy lim
auto c = GetCalculator(vP, calc);
- auto grammage = get<DISPLACEMENT>(c->second)->SolveTrackIntegral(
+ auto grammage = get<eDISPLACEMENT>(c->second)->SolveTrackIntegral(
vP.GetEnergy() / 1_MeV, energy_lim / 1_MeV) *
1_g / square(1_cm);
diff --git a/Processes/Proposal/ContinuousProcess.h b/Processes/Proposal/ContinuousProcess.h
index aefa4ffc6..7d9a88d28 100644
--- a/Processes/Proposal/ContinuousProcess.h
+++ b/Processes/Proposal/ContinuousProcess.h
@@ -27,7 +27,7 @@ namespace corsika::process::proposal {
class ContinuousProcess : public process::ContinuousProcess<ContinuousProcess>,
ProposalProcessBase {
- enum { DISPLACEMENT, SCATTERING };
+ enum { eDISPLACEMENT, eSCATTERING };
using calc_t = std::tuple<std::unique_ptr<PROPOSAL::Displacement>,
std::unique_ptr<PROPOSAL::Scattering>>;
diff --git a/Processes/Proposal/Interaction.cc b/Processes/Proposal/Interaction.cc
index d1b64dea2..623556b90 100644
--- a/Processes/Proposal/Interaction.cc
+++ b/Processes/Proposal/Interaction.cc
@@ -61,23 +61,24 @@ namespace corsika::process::proposal {
std::uniform_real_distribution<double> distr(0., 1.);
// sample a interaction-type, loss and component
- auto rates = get<INTERACTION>(c->second)->Rates(vP.GetEnergy() / 1_MeV);
- auto [type, comp_ptr, v] = get<INTERACTION>(c->second)->SampleLoss(
+ auto rates = get<eINTERACTION>(c->second)->Rates(vP.GetEnergy() / 1_MeV);
+ auto [type, comp_ptr, v] = get<eINTERACTION>(c->second)->SampleLoss(
vP.GetEnergy() / 1_MeV, rates, distr(fRNG));
// Read how much random numbers are required to calculate the secondaries.
// Calculate the secondaries and deploy them on the corsika stack.
- auto rnd = vector<double>(get<SECONDARIES>(c->second)->RequiredRandomNumbers(type));
+ auto rnd = vector<double>(get<eSECONDARIES>(c->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 vP_dir = vP.GetDirection();
+ auto d = vP_dir.GetComponents();
auto direction = PROPOSAL::Vector3D(d.GetX().magnitude(), d.GetY().magnitude(),
d.GetZ().magnitude());
auto loss = make_tuple(static_cast<int>(type), point, direction,
v * vP.GetEnergy() / 1_MeV, 0.);
- auto sec = get<SECONDARIES>(c->second)->CalculateSecondaries(vP.GetEnergy() / 1_MeV,
+ auto sec = get<eSECONDARIES>(c->second)->CalculateSecondaries(vP.GetEnergy() / 1_MeV,
loss, *comp_ptr, rnd);
for (auto& s : sec) {
auto E = get<PROPOSAL::Loss::ENERGY>(s) * 1_MeV;
@@ -85,13 +86,9 @@ namespace corsika::process::proposal {
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);
+ auto p = corsika::stack::MomentumVector(vP_dir.GetCoordinateSystem(), vec);
auto sec_code = corsika::particles::ConvertFromPDG(
static_cast<particles::PDGCode>(get<PROPOSAL::Loss::TYPE>(s)));
- std::cout << " proposal secondary: " << sec_code << " " << E/1_GeV << std::endl;
vP.AddSecondary(make_tuple(sec_code, E, p, vP.GetPosition(), vP.GetTime()));
}
}
@@ -105,7 +102,7 @@ namespace corsika::process::proposal {
if (CanInteract(vP.GetPID())) {
auto c = GetCalculator(vP, calc);
- return get<INTERACTION>(c->second)->MeanFreePath(vP.GetEnergy() / 1_MeV) * 1_g /
+ return get<eINTERACTION>(c->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 007cff4a8..b506dc739 100644
--- a/Processes/Proposal/Interaction.h
+++ b/Processes/Proposal/Interaction.h
@@ -27,7 +27,7 @@ namespace corsika::process::proposal {
//!
class Interaction : public InteractionProcess<Interaction>, ProposalProcessBase {
- enum { SECONDARIES, INTERACTION };
+ enum { eSECONDARIES, eINTERACTION };
using calculator_t = tuple<unique_ptr<PROPOSAL::SecondariesCalculator>,
unique_ptr<PROPOSAL::Interaction>>;
--
GitLab