From 9275f8f5f87b8a5ef4e5ee1dabc9b7e24220b68b Mon Sep 17 00:00:00 2001
From: Nikos Karastathis <n.karastathis@kit.edu>
Date: Thu, 27 Apr 2023 16:43:04 +0000
Subject: [PATCH] Radio output
---
corsika/detail/modules/radio/RadioProcess.inl | 78 ++++++++--
corsika/detail/modules/radio/ZHS.inl | 2 +-
.../detail/modules/radio/antennas/Antenna.inl | 79 -----------
.../radio/antennas/TimeDomainAntenna.inl | 4 +-
corsika/modules/radio/RadioProcess.hpp | 9 +-
corsika/modules/radio/antennas/Antenna.hpp | 16 +--
.../radio/antennas/TimeDomainAntenna.hpp | 14 +-
.../modules/writers/ParticleWriterParquet.hpp | 2 +-
examples/clover_leaf.cpp | 23 ++-
examples/radio_em_shower.cpp | 2 -
examples/synchrotron_test_C8tracking.cpp | 3 +-
tests/modules/testRadio.cpp | 134 +++++++++++++-----
12 files changed, 192 insertions(+), 174 deletions(-)
diff --git a/corsika/detail/modules/radio/RadioProcess.inl b/corsika/detail/modules/radio/RadioProcess.inl
index 0548b71a4..ff8e57010 100644
--- a/corsika/detail/modules/radio/RadioProcess.inl
+++ b/corsika/detail/modules/radio/RadioProcess.inl
@@ -36,13 +36,13 @@ namespace corsika {
TPropagator>::doContinuous(const Step<Particle>& step,
const bool) {
// we want the following particles:
- // Code::Electron & Code::Positron & Code::Gamma
+ // Code::Electron & Code::Positron
// we wrap Simulate() in doContinuous as the plan is to add particle level
// filtering or thinning for calculation of the radio emission. This is
// important for controlling the runtime of radio (by ignoring particles
// that aren't going to contribute i.e. heavy hadrons)
- // if (valid(particle, track)) {
+ // if (valid(step)) {
auto const particleID_{step.getParticlePre().getPID()};
if ((particleID_ == Code::Electron) || (particleID_ == Code::Positron)) {
CORSIKA_LOG_DEBUG("Particle for radio calculation: {} ", particleID_);
@@ -62,17 +62,28 @@ namespace corsika {
return meter * std::numeric_limits<double>::infinity();
}
- // this should all be moved at a separate radio output function
- // LCOV_EXCL_START
template <typename TAntennaCollection, typename TRadioImpl, typename TPropagator>
inline void RadioProcess<TAntennaCollection, TRadioImpl, TPropagator>::startOfLibrary(
const boost::filesystem::path& directory) {
- // loop over every antenna and set the initial path
- // this also writes the time-bins to disk.
- for (auto& antenna : antennas_.getAntennas()) {
- antenna.startOfLibrary(directory, this->implementation().algorithm);
- }
+ // setup the streamer
+ output_.initStreamer((directory / ("antennas.parquet")).string());
+ // LCOV_EXCL_START
+ // build the schema
+ output_.addField("Time", parquet::Repetition::REQUIRED, parquet::Type::DOUBLE,
+ parquet::ConvertedType::NONE);
+
+ output_.addField("Ex", parquet::Repetition::REQUIRED, parquet::Type::DOUBLE,
+ parquet::ConvertedType::NONE);
+
+ output_.addField("Ey", parquet::Repetition::REQUIRED, parquet::Type::DOUBLE,
+ parquet::ConvertedType::NONE);
+
+ output_.addField("Ez", parquet::Repetition::REQUIRED, parquet::Type::DOUBLE,
+ parquet::ConvertedType::NONE);
+ // LCOV_EXCL_STOP
+ // and build the streamer
+ output_.buildStreamer();
}
template <typename TAntennaCollection, typename TRadioImpl, typename TPropagator>
@@ -83,13 +94,55 @@ namespace corsika {
// flush data to disk, and then reset the antenna
// before the next event
for (auto& antenna : antennas_.getAntennas()) {
- antenna.endOfShower(event_, this->implementation().algorithm,
- antenna.getSampleRate() * 1_s);
+
+ auto const sampleRate = antenna.getSampleRate() * 1_s;
+ auto const radioImplementation =
+ static_cast<std::string>(this->implementation().algorithm);
+
+ // get the axis labels for this antenna and write the first row.
+ axistype axis = antenna.implementation().getAxis();
+
+ // get the copy of the waveform data for this event
+ std::vector<double> const& dataX = antenna.implementation().getWaveformX();
+ std::vector<double> const& dataY = antenna.implementation().getWaveformY();
+ std::vector<double> const& dataZ = antenna.implementation().getWaveformZ();
+
+ // check for the axis name
+ std::string label = "Unknown";
+ if (antenna.getDomainLabel() == "Time") {
+ label = "Time";
+ }
+ // LCOV_EXCL_START
+ else if (antenna.getDomainLabel() == "Frequency") {
+ label = "Frequency";
+ }
+ // LCOV_EXCL_STOP
+ if (radioImplementation == "ZHS" && label == "Time") {
+ for (size_t i = 0; i < axis.size() - 1; i++) {
+ auto time = (axis.at(i + 1) + axis.at(i)) / 2.;
+ auto Ex = -(dataX.at(i + 1) - dataX.at(i)) * sampleRate;
+ auto Ey = -(dataY.at(i + 1) - dataY.at(i)) * sampleRate;
+ auto Ez = -(dataZ.at(i + 1) - dataZ.at(i)) * sampleRate;
+
+ *(output_.getWriter())
+ << showerId_ << static_cast<double>(time) << static_cast<double>(Ex)
+ << static_cast<double>(Ey) << static_cast<double>(Ez) << parquet::EndRow;
+ }
+ } else if (radioImplementation == "CoREAS" && label == "Time") {
+ for (size_t i = 0; i < axis.size() - 1; i++) {
+ *(output_.getWriter())
+ << showerId_ << static_cast<double>(axis[i])
+ << static_cast<double>(dataX[i]) << static_cast<double>(dataY[i])
+ << static_cast<double>(dataZ[i]) << parquet::EndRow;
+ }
+ }
+
antenna.reset();
}
+ output_.closeStreamer();
// increment our event counter
- event_++;
+ showerId_++;
}
template <typename TAntennaCollection, typename TRadioImpl, typename TPropagator>
@@ -123,6 +176,5 @@ namespace corsika {
return config;
}
- // LCOV_EXCL_STOP
} // namespace corsika
\ No newline at end of file
diff --git a/corsika/detail/modules/radio/ZHS.inl b/corsika/detail/modules/radio/ZHS.inl
index 5c9e4aacc..6b8b4907a 100644
--- a/corsika/detail/modules/radio/ZHS.inl
+++ b/corsika/detail/modules/radio/ZHS.inl
@@ -214,7 +214,7 @@ namespace corsika {
VectorPotential Vp =
betaPerp * sign * constants * f / denominator / midPaths[i].R_distance_;
antenna.receive(detectionTime1, betaPerp, Vp);
- // intermidiate contributions
+ // intermediate contributions
for (int it{1}; it < numberOfBins; ++it) {
Vp = betaPerp * sign * constants / denominator / midPaths[i].R_distance_;
antenna.receive(
diff --git a/corsika/detail/modules/radio/antennas/Antenna.inl b/corsika/detail/modules/radio/antennas/Antenna.inl
index 89d511257..cb2967f63 100644
--- a/corsika/detail/modules/radio/antennas/Antenna.inl
+++ b/corsika/detail/modules/radio/antennas/Antenna.inl
@@ -28,85 +28,6 @@ namespace corsika {
return name_;
}
- template <typename TAntennaImpl>
- inline void Antenna<TAntennaImpl>::startOfLibrary(
- const boost::filesystem::path& directory, const std::string radioImplementation) {
-
- // calculate and save our filename
- filename_ = (directory / this->getName()).string() + ".npz";
-
- // get the axis labels for this antenna and write the first row.
- axistype axis = this->implementation().getAxis();
-
- // check for the axis name
- std::string label = "Unknown";
- if constexpr (TAntennaImpl::is_time_domain) {
- label = "Time";
- } else if constexpr (TAntennaImpl::is_freq_domain) {
- label = "Frequency";
- }
-
- if (radioImplementation == "ZHS" && TAntennaImpl::is_time_domain) {
- for (size_t i = 0; i < axis.size() - 1; i++) {
- axis.at(i) = (axis.at(i + 1) + axis.at(i)) / 2.;
- }
- // explicitly convert the arrays to the needed type for cnpy
- axis.pop_back();
- long double const* raw_data = axis.data();
- std::vector<size_t> N = {
- axis.size()}; // cnpy needs a vector here -- this should be axis.size() - 1 --
- // write the labels to the first row of the NumPy file
- cnpy::npz_save(filename_, label, raw_data, N, "w");
- } else {
- // explicitly convert the arrays to the needed type for cnpy
- long double const* raw_data = axis.data();
- std::vector<size_t> N = {axis.size()}; // cnpy needs a vector here
- // write the labels to the first row of the NumPy file
- cnpy::npz_save(filename_, label, raw_data, N, "w");
- }
- }
-
- template <typename TAntennaImpl>
- inline void Antenna<TAntennaImpl>::endOfShower(const int event,
- std::string const& radioImplementation,
- const double sampleRate) {
-
- // get the copy of the waveform data for this event
- // we transpose it so that we can match dimensions with the
- // time array that is already in the output file
- std::vector<double> const& dataX = this->implementation().getWaveformX();
- std::vector<double> const& dataY = this->implementation().getWaveformY();
- std::vector<double> const& dataZ = this->implementation().getWaveformZ();
-
- if (radioImplementation == "ZHS") {
- std::vector<double> electricFieldX(dataX.size() - 1,
- 0); // num_bins_, std::vector<double>(3, 0)
- std::vector<double> electricFieldY(dataY.size() - 1, 0);
- std::vector<double> electricFieldZ(dataZ.size() - 1, 0);
- for (size_t i = 0; i < electricFieldX.size(); i++) {
- electricFieldX.at(i) = -(dataX.at(i + 1) - dataX.at(i)) * sampleRate;
- electricFieldY.at(i) = -(dataY.at(i + 1) - dataY.at(i)) * sampleRate;
- electricFieldZ.at(i) = -(dataZ.at(i + 1) - dataZ.at(i)) * sampleRate;
- }
- // cnpy needs a vector for the shape
- cnpy::npz_save(filename_, std::to_string(event) + "X", electricFieldX.data(),
- {electricFieldX.size()}, "a");
- cnpy::npz_save(filename_, std::to_string(event) + "Y", electricFieldY.data(),
- {electricFieldY.size()}, "a");
- cnpy::npz_save(filename_, std::to_string(event) + "Z", electricFieldZ.data(),
- {electricFieldZ.size()}, "a");
- } else {
- // cnpy needs a vector for the shape
- // and write this event to the .npz archive
- cnpy::npz_save(filename_, std::to_string(event) + "X", dataX.data(), {dataX.size()},
- "a");
- cnpy::npz_save(filename_, std::to_string(event) + "Y", dataY.data(), {dataY.size()},
- "a");
- cnpy::npz_save(filename_, std::to_string(event) + "Z", dataZ.data(), {dataZ.size()},
- "a");
- }
- }
-
template <typename TAntennaImpl>
inline TAntennaImpl& Antenna<TAntennaImpl>::implementation() {
return static_cast<TAntennaImpl&>(*this);
diff --git a/corsika/detail/modules/radio/antennas/TimeDomainAntenna.inl b/corsika/detail/modules/radio/antennas/TimeDomainAntenna.inl
index c32f51657..0fa54f21e 100644
--- a/corsika/detail/modules/radio/antennas/TimeDomainAntenna.inl
+++ b/corsika/detail/modules/radio/antennas/TimeDomainAntenna.inl
@@ -97,6 +97,8 @@ namespace corsika {
inline auto const& TimeDomainAntenna::getWaveformZ() const { return waveformEZ_; }
+ inline std::string const TimeDomainAntenna::getDomainLabel() { return "Time"; }
+
inline std::vector<long double> TimeDomainAntenna::createTimeAxis() const {
// create a 1-D xtensor to store time values so we can print them later.
@@ -115,7 +117,7 @@ namespace corsika {
return times;
}
- inline auto const& TimeDomainAntenna::getAxis() const { return time_axis_; }
+ inline auto const TimeDomainAntenna::getAxis() const { return time_axis_; }
inline InverseTimeType const& TimeDomainAntenna::getSampleRate() const {
return sample_rate_;
diff --git a/corsika/modules/radio/RadioProcess.hpp b/corsika/modules/radio/RadioProcess.hpp
index 9ae9087bd..43b7a3453 100644
--- a/corsika/modules/radio/RadioProcess.hpp
+++ b/corsika/modules/radio/RadioProcess.hpp
@@ -7,11 +7,8 @@
*/
#pragma once
-#include <istream>
-#include <fstream>
-#include <iostream>
-#include <string>
#include <corsika/output/BaseOutput.hpp>
+#include <corsika/output/ParquetStreamer.hpp>
#include <corsika/framework/process/ContinuousProcess.hpp>
#include <corsika/framework/core/Step.hpp>
#include <corsika/setup/SetupStack.hpp>
@@ -49,9 +46,11 @@ namespace corsika {
protected:
TAntennaCollection& antennas_; ///< The radio antennas we store into.
TPropagator propagator_; ///< The propagator implementation.
- int event_{0}; ///< The current event ID.
+ unsigned int showerId_{0}; ///< The current event ID.
+ ParquetStreamer output_; //!< The parquet streamer for this process.
public:
+ using axistype = std::vector<long double>;
/**
* Construct a new RadioProcess.
*/
diff --git a/corsika/modules/radio/antennas/Antenna.hpp b/corsika/modules/radio/antennas/Antenna.hpp
index b5aeb8fb9..3ae084e71 100644
--- a/corsika/modules/radio/antennas/Antenna.hpp
+++ b/corsika/modules/radio/antennas/Antenna.hpp
@@ -7,10 +7,10 @@
*/
#pragma once
-#include <cnpy.hpp>
#include <boost/filesystem.hpp>
#include <corsika/framework/geometry/Point.hpp>
#include <corsika/framework/core/PhysicalGeometry.hpp>
+#include <corsika/output/ParquetStreamer.hpp>
namespace corsika {
@@ -28,7 +28,6 @@ namespace corsika {
std::string const name_; ///< The name/identifier of this antenna.
Point const location_; ///< The location of this antenna.
CoordinateSystemPtr const coordinateSystem_; ///< The coordinate system of the antenna
- std::string filename_ = ""; ///< The filename for the output file for this antenna.
public:
using axistype = std::vector<long double>;
@@ -102,19 +101,6 @@ namespace corsika {
*/
std::vector<double> const& getWaveformZ() const;
- /**
- * Prepare for the start of the library.
- */
- void startOfLibrary(boost::filesystem::path const& directory,
- std::string const radioImplementation);
-
- /**
- * Flush the data from this shower to disk.
- */
- void endOfShower(int const event, std::string const& radioImplementation,
- double const sampleRate);
-
- protected:
/**
* Get a reference to the underlying radio implementation.
*/
diff --git a/corsika/modules/radio/antennas/TimeDomainAntenna.hpp b/corsika/modules/radio/antennas/TimeDomainAntenna.hpp
index 85fd9dd10..8c4ee2147 100644
--- a/corsika/modules/radio/antennas/TimeDomainAntenna.hpp
+++ b/corsika/modules/radio/antennas/TimeDomainAntenna.hpp
@@ -33,10 +33,6 @@ namespace corsika {
public:
// import the methods from the antenna
-
- // label this as a time-domain antenna.
- static constexpr bool is_time_domain{true};
-
using Antenna<TimeDomainAntenna>::getName;
using Antenna<TimeDomainAntenna>::getLocation;
@@ -98,6 +94,14 @@ namespace corsika {
*/
auto const& getWaveformZ() const;
+ /**
+ * Return a label that indicates that this is a time
+ * domain antenna
+ *
+ * This returns the string "Time".
+ */
+ std::string const getDomainLabel();
+
/**
* Creates time-units of each waveform.
*
@@ -110,7 +114,7 @@ namespace corsika {
*
* This returns them in nanoseconds for ease of use.
*/
- auto const& getAxis() const;
+ auto const getAxis() const;
/**
* Returns the sampling rate of the time domain antenna.
diff --git a/corsika/modules/writers/ParticleWriterParquet.hpp b/corsika/modules/writers/ParticleWriterParquet.hpp
index d28abdd1c..2a431bc09 100644
--- a/corsika/modules/writers/ParticleWriterParquet.hpp
+++ b/corsika/modules/writers/ParticleWriterParquet.hpp
@@ -71,7 +71,7 @@ namespace corsika {
double countHadrons_ = 0; ///< count hadrons hitting plane
double countMuons_ = 0; ///< count muons hitting plane
double countEM_ = 0; ///< count EM particles hitting plane.
- double countOthers_ = 0; ///< count othe types of particles hitting plane
+ double countOthers_ = 0; ///< count other types of particles hitting plane
HEPEnergyType totalEnergy_; ///< energy absorbed in ground.
diff --git a/examples/clover_leaf.cpp b/examples/clover_leaf.cpp
index ccb8c7e79..0a93e1d9e 100644
--- a/examples/clover_leaf.cpp
+++ b/examples/clover_leaf.cpp
@@ -1,5 +1,5 @@
/*
- * (c) Copyright 2018 CORSIKA Project, corsika-project@lists.kit.edu
+ * (c) Copyright 2022 CORSIKA Project, corsika-project@lists.kit.edu
*
* This software is distributed under the terms of the GNU General Public
* Licence version 3 (GPL Version 3). See file LICENSE for a full version of
@@ -37,8 +37,6 @@
#include <corsika/modules/radio/detectors/AntennaCollection.hpp>
#include <corsika/modules/radio/propagators/StraightPropagator.hpp>
#include <corsika/modules/radio/propagators/SimplePropagator.hpp>
-#include <corsika/modules/radio/propagators/SignalPath.hpp>
-#include <corsika/modules/radio/propagators/RadioPropagator.hpp>
#include <corsika/modules/TrackWriter.hpp>
#include <corsika/modules/StackInspector.hpp>
@@ -47,13 +45,13 @@
//#include <corsika/modules/TrackWriter.hpp>
/*
- NOTE, WARNING, ATTENTION
+ NOTE, WARNING, ATTENTION
- The .../Random.hpppp implement the hooks of external modules to the C8 random
- number generator. It has to occur excatly ONCE per linked
- executable. If you include the header below multiple times and
- link this together, it will fail.
- */
+ The .../Random.hpppp implement the hooks of external modules to the C8 random
+ number generator. It has to occur excatly ONCE per linked
+ executable. If you include the header below multiple times and
+ link this together, it will fail.
+*/
#include <corsika/modules/sibyll/Random.hpp>
#include <corsika/modules/urqmd/Random.hpp>
@@ -72,7 +70,7 @@ using namespace std;
//
int main() {
- logging::set_level(logging::level::warn);
+ logging::set_level(logging::level::info);
corsika_logger->set_pattern("[%n:%^%-8l%$] custom pattern: %v");
CORSIKA_LOG_INFO("Synchrotron radiation");
@@ -210,11 +208,10 @@ int main() {
// assemble all processes into an ordered process list
auto sequence = make_sequence(coreas, cut);
+ output.startOfLibrary();
// define air shower object, run simulation
Cascade EAS(env, tracking, sequence, output, stack);
- output.startOfShower();
EAS.run();
- output.endOfShower();
CORSIKA_LOG_INFO("|p| electron = {} and E electron = {}", plab.getNorm(), Elab);
CORSIKA_LOG_INFO("period: {}", period);
@@ -225,4 +222,4 @@ int main() {
CORSIKA_LOG_INFO("gamma: {}", gammaP);
output.endOfLibrary();
-}
+}
\ No newline at end of file
diff --git a/examples/radio_em_shower.cpp b/examples/radio_em_shower.cpp
index 93eb14208..a088288de 100644
--- a/examples/radio_em_shower.cpp
+++ b/examples/radio_em_shower.cpp
@@ -51,8 +51,6 @@
#include <corsika/modules/radio/detectors/AntennaCollection.hpp>
#include <corsika/modules/radio/propagators/StraightPropagator.hpp>
#include <corsika/modules/radio/propagators/SimplePropagator.hpp>
-#include <corsika/modules/radio/propagators/SignalPath.hpp>
-#include <corsika/modules/radio/propagators/RadioPropagator.hpp>
#include <corsika/setup/SetupStack.hpp>
#include <corsika/setup/SetupTrajectory.hpp>
diff --git a/examples/synchrotron_test_C8tracking.cpp b/examples/synchrotron_test_C8tracking.cpp
index cf1792188..9654d6ef2 100644
--- a/examples/synchrotron_test_C8tracking.cpp
+++ b/examples/synchrotron_test_C8tracking.cpp
@@ -158,11 +158,10 @@ int main() {
// assemble all processes into an ordered process list
auto sequence = make_sequence(coreas, zhs, cut);
+ output.startOfLibrary();
// define air shower object, run simulation
Cascade EAS(env, tracking, sequence, output, stack);
- output.startOfShower();
EAS.run();
- output.endOfShower();
CORSIKA_LOG_INFO("|p| = {} and E = {}", plab.getNorm(), Elab);
CORSIKA_LOG_INFO("period: {}", period);
diff --git a/tests/modules/testRadio.cpp b/tests/modules/testRadio.cpp
index e31434fd4..58742fb83 100644
--- a/tests/modules/testRadio.cpp
+++ b/tests/modules/testRadio.cpp
@@ -7,6 +7,7 @@
*/
#include <catch2/catch.hpp>
+#include <corsika/modules/radio/RadioProcess.hpp>
#include <corsika/modules/radio/ZHS.hpp>
#include <corsika/modules/radio/CoREAS.hpp>
#include <corsika/modules/radio/antennas/TimeDomainAntenna.hpp>
@@ -65,9 +66,10 @@ TEST_CASE("Radio", "[processes]") {
SECTION("CoREAS process") {
- // This serves as a compiler test for any changes in the CoREAS algorithm
- // Environment
+ // This serves as a compiler test for any changes in the CoREAS algorithm and
+ // check the radio process output
+ // Environment
using EnvironmentInterface =
IRefractiveIndexModel<IMediumPropertyModel<IMagneticFieldModel<IMediumModel>>>;
@@ -178,6 +180,23 @@ TEST_CASE("Radio", "[processes]") {
REQUIRE(total < (1e-12 * ant.getWaveformX().size()));
}
+ // coreas output check
+ std::string const implemencoreas{"CoREAS"};
+
+ boost::filesystem::path const tempPathC{boost::filesystem::temp_directory_path() /
+ ("test_corsika_radio_" + implemencoreas)};
+ if (boost::filesystem::exists(tempPathC)) {
+ boost::filesystem::remove_all(tempPathC);
+ }
+
+ boost::filesystem::create_directory(tempPathC);
+ coreas.startOfLibrary(tempPathC);
+ auto const outputFileC = tempPathC / ("antennas.parquet");
+ CHECK(boost::filesystem::exists(outputFileC));
+ // run end of shower and make sure that something extra was added
+ auto const fileSizeC = boost::filesystem::file_size(outputFileC);
+ coreas.endOfShower(0);
+ CHECK(boost::filesystem::file_size(outputFileC) > fileSizeC);
coreas.endOfLibrary();
} // END: SECTION("CoREAS process")
@@ -374,6 +393,25 @@ TEST_CASE("Radio", "[processes]") {
// check doContinuous and simulate methods
zhs.doContinuous(step, true);
+ // zhs output check
+ std::string const implemenzhs{"ZHS"};
+
+ boost::filesystem::path const tempPathZ{boost::filesystem::temp_directory_path() /
+ ("test_corsika_radio_" + implemenzhs)};
+ if (boost::filesystem::exists(tempPathZ)) {
+ boost::filesystem::remove_all(tempPathZ);
+ }
+
+ boost::filesystem::create_directory(tempPathZ);
+ zhs.startOfLibrary(tempPathZ);
+ auto const outputFileZ = tempPathZ / ("antennas.parquet");
+ CHECK(boost::filesystem::exists(outputFileZ));
+ // run end of shower and make sure that something extra was added
+ auto const fileSizeZ = boost::filesystem::file_size(outputFileZ);
+ zhs.endOfShower(0);
+ CHECK(boost::filesystem::file_size(outputFileZ) > fileSizeZ);
+ zhs.endOfLibrary();
+
} // END: SECTION("ZHS process")
SECTION("Radio extreme cases") {
@@ -399,6 +437,7 @@ TEST_CASE("Radio", "[processes]") {
Point const point_1(rootCSRadio, {1_m, 1_m, 0_m});
Point const point_2(rootCSRadio, {2_km, 1_km, 0_m});
Point const point_4(rootCSRadio, {0_m, 1_m, 0_m});
+ const auto point_b{Point(rootCSRadio, 30000_m, 0_m, 0_m)};
// create times for the antenna
const TimeType start{0_s};
@@ -407,17 +446,26 @@ TEST_CASE("Radio", "[processes]") {
const TimeType duration_dummy{2_s};
const InverseTimeType sample_dummy{1_Hz};
+ // create specific times for antenna to do timebin check
+ const TimeType start_b{0.994e-4_s};
+ const TimeType duration_b{1.07e-4_s - 0.994e-4_s};
+ const InverseTimeType sampleRate_b{5e+11_Hz};
+
// check that I can create an antenna at (1, 2, 3)
TimeDomainAntenna ant1("antenna_name", point1, rootCSRadio, start, duration, sample,
start);
TimeDomainAntenna ant2("dummy", point1, rootCSRadio, start, duration_dummy,
sample_dummy, start);
+ TimeDomainAntenna ant_b("timebin", point_b, rootCSRadio, start_b, duration_b,
+ sampleRate_b, start_b);
// construct a radio detector instance to store our antennas
AntennaCollection<TimeDomainAntenna> detector;
AntennaCollection<TimeDomainAntenna> detector_dummy;
+ AntennaCollection<TimeDomainAntenna> detector_b;
// add the antennas to the detector
detector.addAntenna(ant1);
detector_dummy.addAntenna(ant2);
+ detector_b.addAntenna(ant_b);
// create a new stack for each trial
setup::Stack<EnvType> stack;
@@ -428,7 +476,7 @@ TEST_CASE("Radio", "[processes]") {
const auto pmass{get_mass(particle)};
// construct an energy
const HEPEnergyType E0{1_TeV};
- // compute the necessary momentumn
+ // compute the necessary momentum
const HEPMomentumType P0{sqrt(E0 * E0 - pmass * pmass)};
// and create the momentum vector
const auto plab{MomentumVector(rootCSRadio, {P0, 0_GeV, 0_GeV})};
@@ -455,7 +503,18 @@ TEST_CASE("Radio", "[processes]") {
VelocityVector vh{(point_4 - point1) / th};
Line lh{point1, vh};
StraightTrajectory track_h{lh, th};
+ StraightTrajectory track_h_neg_time{lh, -th};
Step step_h(particle_stack, track_h);
+ Step step_h_neg_time(particle_stack, track_h_neg_time);
+
+ // feed radio with an electron track that ends in a different antenna bin.
+ Point const point_start(rootCSRadio, {100_m, 0_m, 0_m});
+ Point const point_end(rootCSRadio, {100_m, 0.00628319_m, 0_m});
+ TimeType tb{(point_end - point_start).getNorm() / (0.999 * constants::c)};
+ VelocityVector vb{(point_end - point_start) / tb};
+ Line lb{point_start, vb};
+ StraightTrajectory track_b{lb, tb};
+ Step step_b(particle_stack, track_b);
// create radio process instances
RadioProcess<decltype(detector),
@@ -467,13 +526,12 @@ TEST_CASE("Radio", "[processes]") {
ZHS<decltype(detector), decltype(SimplePropagator(envRadio))>,
decltype(SimplePropagator(envRadio))>
zhs(detector, envRadio);
-
coreas.doContinuous(step_proton, true);
zhs.doContinuous(step_proton, true);
coreas.doContinuous(step_h, true);
zhs.doContinuous(step_h, true);
+ zhs.doContinuous(step_h_neg_time, true);
- // create radio processes with "dummy" antenna to trigger extreme time-binning
RadioProcess<decltype(detector_dummy),
CoREAS<decltype(detector_dummy), decltype(SimplePropagator(envRadio))>,
decltype(SimplePropagator(envRadio))>
@@ -483,10 +541,24 @@ TEST_CASE("Radio", "[processes]") {
ZHS<decltype(detector_dummy), decltype(SimplePropagator(envRadio))>,
decltype(SimplePropagator(envRadio))>
zhs_dummy(detector_dummy, envRadio);
-
+ coreas_dummy.doContinuous(step_proton, true);
+ zhs_dummy.doContinuous(step_proton, true);
coreas_dummy.doContinuous(step_h, true);
zhs_dummy.doContinuous(step_h, true);
+ // create radio process instances
+ RadioProcess<decltype(detector_b),
+ CoREAS<decltype(detector_b), decltype(SimplePropagator(envRadio))>,
+ decltype(SimplePropagator(envRadio))>
+ coreas_b(detector_b, envRadio);
+
+ RadioProcess<decltype(detector_b),
+ ZHS<decltype(detector_b), decltype(SimplePropagator(envRadio))>,
+ decltype(SimplePropagator(envRadio))>
+ zhs_b(detector_b, envRadio);
+ coreas_b.doContinuous(step_b, true);
+ zhs_b.doContinuous(step_b, true);
+
} // END: SECTION("Radio extreme cases")
SECTION("Process Library") {
@@ -574,7 +646,7 @@ TEST_CASE("Radio", "[processes]") {
CHECK(config["antennas"]["antenna_name2"]["location"][0].as<double>() == 4);
CHECK(config["antennas"]["antenna_name2"]["location"][1].as<double>() == 80);
CHECK(config["antennas"]["antenna_name2"]["location"][2].as<double>() == 6);
- }
+ } // END: SECTION("Process Library")
} // END: TEST_CASE("Radio", "[processes]")
@@ -827,7 +899,7 @@ TEST_CASE("Antennas") {
} // END: SECTION("TimeDomainAntenna AntennaCollection")
- SECTION("TimeDomainAntenna Library") {
+ SECTION("TimeDomainAntenna Config File") {
// Runs checks that the file readers are working properly
Environment<IRefractiveIndexModel<IMediumModel>> env;
const auto rootCS = env.getCoordinateSystem();
@@ -837,36 +909,24 @@ TEST_CASE("Antennas") {
InverseTimeType const sampleRate(1_GHz);
TimeType const groundHitTime(1e3_ns);
- TimeDomainAntenna antenna("test_antenna", antPos, rootCS, tStart, duration,
- sampleRate, groundHitTime);
-
- // Run the start of lib and end of shower functions on the antenna
- std::vector<std::string> implementations{"CoREAS", "ZHS"};
- for (auto const implemen : implementations) {
- // For each implementation type, save a file
- boost::filesystem::path const tempPath{boost::filesystem::temp_directory_path() /
- ("test_corsika_radio_" + implemen)};
- if (boost::filesystem::exists(tempPath)) {
- boost::filesystem::remove_all(tempPath);
- }
- boost::filesystem::create_directory(tempPath);
- antenna.startOfLibrary(tempPath, implemen);
- auto const outputFile = tempPath / (antenna.getName() + ".npz");
- CHECK(boost::filesystem::exists(outputFile));
-
- // run end of shower and make sure that something extra was added
- auto const fileSize = boost::filesystem::file_size(outputFile);
- antenna.endOfShower(0, implemen, sampleRate / 1_Hz);
- CHECK(boost::filesystem::file_size(outputFile) > fileSize);
- }
+ TimeDomainAntenna antennaC("test_antennaCoREAS", antPos, rootCS, tStart, duration,
+ sampleRate, groundHitTime);
+ TimeDomainAntenna antennaZ("test_antennaZHS", antPos, rootCS, tStart, duration,
+ sampleRate, groundHitTime);
// Check the YAML file output
- auto const config = antenna.getConfig();
- CHECK(config["type"].as<std::string>() == "TimeDomainAntenna");
- CHECK(config["start_time"].as<double>() == tStart / 1_ns);
- CHECK(config["duration"].as<double>() == duration / 1_ns);
- CHECK(config["sample_rate"].as<double>() == sampleRate / 1_GHz);
- } // END: SECTION("TimeDomainAntenna Library")
+ auto const configC = antennaC.getConfig();
+ CHECK(configC["type"].as<std::string>() == "TimeDomainAntenna");
+ CHECK(configC["start_time"].as<double>() == tStart / 1_ns);
+ CHECK(configC["duration"].as<double>() == duration / 1_ns);
+ CHECK(configC["sample_rate"].as<double>() == sampleRate / 1_GHz);
+
+ auto const configZ = antennaZ.getConfig();
+ CHECK(configZ["type"].as<std::string>() == "TimeDomainAntenna");
+ CHECK(configZ["start_time"].as<double>() == tStart / 1_ns);
+ CHECK(configZ["duration"].as<double>() == duration / 1_ns);
+ CHECK(configZ["sample_rate"].as<double>() == sampleRate / 1_GHz);
+ } // END: SECTION("TimeDomainAntenna Config File")
} // END: TEST_CASE("Antennas")
@@ -1156,7 +1216,7 @@ TEST_CASE("Propagators") {
Approx(0).margin(absMargin));
CHECK(path.average_refractive_index_ == Approx(0.210275935));
CHECK(path.refractive_index_source_ == Approx(2));
- // CHECK(path.refractive_index_destination_ == Approx(0.0000000041));
+ CHECK(path.refractive_index_destination_ == Approx(4.12231e-09));
CHECK(path.emit_.getComponents() == vvv1.getComponents());
CHECK(path.receive_.getComponents() == vvv2.getComponents());
CHECK(path.R_distance_ == 10_m);
--
GitLab