From a5d290bc5296218575032fb01f333395eee4972e Mon Sep 17 00:00:00 2001
From: Maximilian Reininghaus <maximilian.reininghaus@kit.edu>
Date: Fri, 17 Jul 2020 17:21:12 +0200
Subject: [PATCH] x,y coordinates for ObservationPlane
---
Documentation/Examples/vertical_EAS.cc | 20 +++++++++++++++----
.../ObservationPlane/ObservationPlane.cc | 13 +++++++++---
Processes/ObservationPlane/ObservationPlane.h | 5 ++++-
3 files changed, 30 insertions(+), 8 deletions(-)
diff --git a/Documentation/Examples/vertical_EAS.cc b/Documentation/Examples/vertical_EAS.cc
index 388255bfc..0f4d57c87 100644
--- a/Documentation/Examples/vertical_EAS.cc
+++ b/Documentation/Examples/vertical_EAS.cc
@@ -162,8 +162,20 @@ int main(int argc, char** argv) {
stack.AddParticle(std::make_tuple(beamCode, E0, plab, injectionPos, 0_ns, A, Z));
} else {
- stack.AddParticle(
- std::make_tuple(particles::Code::Proton, E0, plab, injectionPos, 0_ns));
+ if (Z == 1) {
+ stack.AddParticle(std::tuple<particles::Code, units::si::HEPEnergyType,
+ corsika::stack::MomentumVector, geometry::Point,
+ units::si::TimeType>{particles::Code::Proton, E0, plab,
+ injectionPos, 0_ns});
+ } else if (Z == 0) {
+ stack.AddParticle(std::tuple<particles::Code, units::si::HEPEnergyType,
+ corsika::stack::MomentumVector, geometry::Point,
+ units::si::TimeType>{particles::Code::Neutron, E0,
+ plab, injectionPos, 0_ns});
+ } else {
+ std::cerr << "illegal parameters" << std::endl;
+ return EXIT_FAILURE;
+ }
}
// we make the axis much longer than the inj-core distance since the
@@ -217,8 +229,8 @@ int main(int argc, char** argv) {
process::longitudinal_profile::LongitudinalProfile longprof{showerAxis};
Plane const obsPlane(showerCore, Vector<dimensionless_d>(rootCS, {0., 0., 1.}));
- process::observation_plane::ObservationPlane observationLevel(obsPlane,
- "particles.dat");
+ process::observation_plane::ObservationPlane observationLevel(
+ obsPlane, Vector<dimensionless_d>(rootCS, {1., 0., 0.}), "particles.dat");
process::UrQMD::UrQMD urqmd;
process::interaction_counter::InteractionCounter urqmdCounted{urqmd};
diff --git a/Processes/ObservationPlane/ObservationPlane.cc b/Processes/ObservationPlane/ObservationPlane.cc
index ab3a80232..7abff3513 100644
--- a/Processes/ObservationPlane/ObservationPlane.cc
+++ b/Processes/ObservationPlane/ObservationPlane.cc
@@ -14,14 +14,18 @@
using namespace corsika::process::observation_plane;
using namespace corsika::units::si;
-ObservationPlane::ObservationPlane(geometry::Plane const& obsPlane,
- std::string const& filename, bool deleteOnHit)
+ObservationPlane::ObservationPlane(
+ geometry::Plane const& obsPlane,
+ geometry::Vector<units::si::dimensionless_d> const& x_axis,
+ std::string const& filename, bool deleteOnHit)
: plane_(obsPlane)
, outputStream_(filename)
, deleteOnHit_(deleteOnHit)
, energy_ground_(0_GeV)
, count_ground_(0) {
- outputStream_ << "#PDG code, energy / eV, distance to center / m" << std::endl;
+ , xAxis_(x_axis.normalized())
+ , yAxis_(obsPlane.GetNormal().cross(xAxis_)) {
+ outputStream_ << "#PDG code, energy / eV, x distance / m, y distance / m" << std::endl;
}
corsika::process::EProcessReturn ObservationPlane::DoContinuous(
@@ -37,8 +41,11 @@ corsika::process::EProcessReturn ObservationPlane::DoContinuous(
}
const auto energy = particle.GetEnergy();
+ auto const displacement = trajectory.GetPosition(1) - plane_.GetCenter();
+
outputStream_ << static_cast<int>(particles::GetPDG(particle.GetPID())) << ' '
<< energy / 1_eV << ' '
+ << displacement.dot(xAxis_) / 1_m << ' ' << displacement.dot(yAxis_) / 1_m
<< (trajectory.GetPosition(1) - plane_.GetCenter()).norm() / 1_m
<< std::endl;
diff --git a/Processes/ObservationPlane/ObservationPlane.h b/Processes/ObservationPlane/ObservationPlane.h
index 86c5ba0f2..3b530f621 100644
--- a/Processes/ObservationPlane/ObservationPlane.h
+++ b/Processes/ObservationPlane/ObservationPlane.h
@@ -26,7 +26,9 @@ namespace corsika::process::observation_plane {
class ObservationPlane : public corsika::process::ContinuousProcess<ObservationPlane> {
public:
- ObservationPlane(geometry::Plane const&, std::string const&, bool = true);
+ ObservationPlane(geometry::Plane const&,
+ geometry::Vector<units::si::dimensionless_d> const&,
+ std::string const&, bool = true);
corsika::process::EProcessReturn DoContinuous(
corsika::setup::Stack::ParticleType& vParticle,
@@ -47,5 +49,6 @@ namespace corsika::process::observation_plane {
units::si::HEPEnergyType energy_ground_ = 0 * units::si::electronvolt;
unsigned int count_ground_ = 0;
+ geometry::Vector<units::si::dimensionless_d> const xAxis_, yAxis_;
};
} // namespace corsika::process::observation_plane
--
GitLab