From 75c50f103fcc7a2b986806b797842340221d29d0 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.
---
Processes/Proposal/ContinuousProcess.cc | 12 ++++++------
Processes/Proposal/Interaction.cc | 9 +++------
2 files changed, 9 insertions(+), 12 deletions(-)
diff --git a/Processes/Proposal/ContinuousProcess.cc b/Processes/Proposal/ContinuousProcess.cc
index d515dd3a9..e552f49dd 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());
@@ -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());
@@ -138,7 +138,7 @@ 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);
diff --git a/Processes/Proposal/Interaction.cc b/Processes/Proposal/Interaction.cc
index d1b64dea2..4a1b767be 100644
--- a/Processes/Proposal/Interaction.cc
+++ b/Processes/Proposal/Interaction.cc
@@ -72,7 +72,8 @@ namespace corsika::process::proposal {
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,
@@ -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()));
}
}
--
GitLab