diff --git a/corsika/detail/framework/geometry/FourVector.inl b/corsika/detail/framework/geometry/FourVector.inl
index 7219a920b9d81b41687016cb3202eb2a837ed400..5d871430dad5390aaa7b80d8fb83d7e16c165f03 100644
--- a/corsika/detail/framework/geometry/FourVector.inl
+++ b/corsika/detail/framework/geometry/FourVector.inl
@@ -63,39 +63,39 @@ namespace corsika {
}
template <typename TTimeType, typename TSpaceVecType>
- inline FourVector<TTimeType, TSpaceVecType>&
- FourVector<TTimeType, TSpaceVecType>::operator-=(FourVector const& b) {
+ inline FourVector<TTimeType, TSpaceVecType>& FourVector<TTimeType, TSpaceVecType>::
+ operator-=(FourVector const& b) {
timeLike_ -= b.timeLike_;
spaceLike_ -= b.spaceLike_;
return *this;
}
template <typename TTimeType, typename TSpaceVecType>
- inline FourVector<TTimeType, TSpaceVecType>&
- FourVector<TTimeType, TSpaceVecType>::operator*=(double const b) {
+ inline FourVector<TTimeType, TSpaceVecType>& FourVector<TTimeType, TSpaceVecType>::
+ operator*=(double const b) {
timeLike_ *= b;
spaceLike_ *= b;
return *this;
}
template <typename TTimeType, typename TSpaceVecType>
- inline FourVector<TTimeType, TSpaceVecType>&
- FourVector<TTimeType, TSpaceVecType>::operator/=(double const b) {
+ inline FourVector<TTimeType, TSpaceVecType>& FourVector<TTimeType, TSpaceVecType>::
+ operator/=(double const b) {
timeLike_ /= b;
spaceLike_.getComponents() /= b;
return *this;
}
template <typename TTimeType, typename TSpaceVecType>
- inline FourVector<TTimeType, TSpaceVecType>&
- FourVector<TTimeType, TSpaceVecType>::operator/(double const b) {
+ inline FourVector<TTimeType, TSpaceVecType>& FourVector<TTimeType, TSpaceVecType>::
+ operator/(double const b) {
*this /= b;
return *this;
}
template <typename TTimeType, typename TSpaceVecType>
inline typename FourVector<TTimeType, TSpaceVecType>::norm_type
- FourVector<TTimeType, TSpaceVecType>::operator*(FourVector const& b) {
+ FourVector<TTimeType, TSpaceVecType>::operator*(FourVector const& b) {
if constexpr (std::is_same<time_type, decltype(std::declval<space_type>() / meter *
second)>::value)
return timeLike_ * b.timeLike_ * constants::cSquared - spaceLike_.norm();
diff --git a/corsika/detail/framework/geometry/QuantityVector.inl b/corsika/detail/framework/geometry/QuantityVector.inl
index c7f9457b98dc0935c710f6d93fef6c53aab69f56..a932e4de8b596217bde32b1bb672e69866f05dc6 100644
--- a/corsika/detail/framework/geometry/QuantityVector.inl
+++ b/corsika/detail/framework/geometry/QuantityVector.inl
@@ -17,8 +17,8 @@
namespace corsika {
template <typename TDimension>
- inline typename QuantityVector<TDimension>::quantity_type
- QuantityVector<TDimension>::operator[](size_t const index) const {
+ inline typename QuantityVector<TDimension>::quantity_type QuantityVector<TDimension>::
+ operator[](size_t const index) const {
return quantity_type(phys::units::detail::magnitude_tag, eigenVector_[index]);
}
diff --git a/corsika/detail/media/WeightProvider.inl b/corsika/detail/media/WeightProvider.inl
index 2bbc58dd3b994977cb8dfcc8e373b7e9e10cb053..b9cf1be7ae0e94d0a99b66366373ae62861fffc6 100644
--- a/corsika/detail/media/WeightProvider.inl
+++ b/corsika/detail/media/WeightProvider.inl
@@ -20,7 +20,7 @@ namespace corsika {
template <class AConstIterator, class BConstIterator>
inline typename WeightProviderIterator<AConstIterator, BConstIterator>::value_type
- WeightProviderIterator<AConstIterator, BConstIterator>::operator*() const {
+ WeightProviderIterator<AConstIterator, BConstIterator>::operator*() const {
return ((*aIter_) * (*bIter_)).magnitude();
}
diff --git a/corsika/detail/modules/TimeCut.inl b/corsika/detail/modules/TimeCut.inl
index 83cf4a57faaf19910dc611c70fc9592177e892b7..3dbdecb7f28dd0c7804308f2841b5e7ae78ebed1 100644
--- a/corsika/detail/modules/TimeCut.inl
+++ b/corsika/detail/modules/TimeCut.inl
@@ -16,8 +16,7 @@ namespace corsika {
: time_(time) {}
template <typename Particle>
- inline ProcessReturn TimeCut::doContinuous(Step<Particle> const& step,
- bool const) {
+ inline ProcessReturn TimeCut::doContinuous(Step<Particle> const& step, bool const) {
CORSIKA_LOG_TRACE("TimeCut::doContinuous");
if (step.getTimePost() >= time_) {
CORSIKA_LOG_TRACE("stopping continuous process");
diff --git a/corsika/detail/modules/proposal/ProposalProcessBase.inl b/corsika/detail/modules/proposal/ProposalProcessBase.inl
index 936798ea1cd05340ba2f0c4b397d5af9c47a068e..c7c38c10423ab6fe270c0851f7b77ca21d8ab73c 100644
--- a/corsika/detail/modules/proposal/ProposalProcessBase.inl
+++ b/corsika/detail/modules/proposal/ProposalProcessBase.inl
@@ -56,8 +56,8 @@ namespace corsika::proposal {
PROPOSAL::InterpolationSettings::TABLES_PATH = corsika_data("PROPOSAL").c_str();
}
- inline size_t ProposalProcessBase::hash::operator()(
- const calc_key_t& p) const noexcept {
+ inline size_t ProposalProcessBase::hash::operator()(const calc_key_t& p) const
+ noexcept {
return p.first ^ std::hash<Code>{}(p.second);
}
diff --git a/corsika/detail/modules/radio/CoREAS.inl b/corsika/detail/modules/radio/CoREAS.inl
index fc474ea4ab1a859e487c6e120586b851a3e7f4f9..9f25ae97265f97cddf955f19ad2eee0c199c01b8 100644
--- a/corsika/detail/modules/radio/CoREAS.inl
+++ b/corsika/detail/modules/radio/CoREAS.inl
@@ -18,13 +18,14 @@ namespace corsika {
Step<Particle> const& step) {
// get the global simulation time for that track.
- auto const startTime_{step.getTimePre()}; // time at the start point of the track. I
- // should use something similar to fCoreHitTime (?)
+ auto const startTime_{
+ step.getTimePre()}; // time at the start point of the track. I
+ // should use something similar to fCoreHitTime (?)
auto const endTime_{step.getTimePost()}; // time at end point of track.
- // LCOV_EXCL_START
+ // LCOV_EXCL_START
if (startTime_ == endTime_) {
- CORSIKA_LOG_ERROR("Time at the start and end of the track coincides! - radio");
+ CORSIKA_LOG_ERROR("Time at the start and end of the track coincides! - radio");
return ProcessReturn::Ok;
// LCOV_EXCL_STOP
} else {
@@ -33,7 +34,7 @@ namespace corsika {
Point const startPoint_{step.getPositionPre()};
Point const endPoint_{step.getPositionPost()};
// get the coordinate system of the startpoint and hence the track
- auto const cs_ {startPoint_.getCoordinateSystem()};
+ auto const cs_{startPoint_.getCoordinateSystem()};
auto const currDirection{(endPoint_ - startPoint_).normalized()};
// calculate the track length
@@ -41,7 +42,7 @@ namespace corsika {
// beta is velocity / speed of light. Start & end should be the same in endpoints!
auto const corrBetaValue{(endPoint_ - startPoint_).getNorm() /
- (constants::c * (endTime_ - startTime_))};
+ (constants::c * (endTime_ - startTime_))};
auto const beta_{currDirection * corrBetaValue};
// get particle charge
@@ -144,7 +145,7 @@ namespace corsika {
if ((paths1[i].refractive_index_destination_ > 1) &&
((std::fabs(preDoppler_) < approxThreshold_) ||
(std::fabs(postDoppler_) < approxThreshold_))) {
- CORSIKA_LOG_DEBUG("Used ZHS-like approximation in CoREAS - radio");
+ CORSIKA_LOG_DEBUG("Used ZHS-like approximation in CoREAS - radio");
// clear the existing paths for this particle and track, since we don't need
// them anymore
diff --git a/corsika/detail/modules/radio/ZHS.inl b/corsika/detail/modules/radio/ZHS.inl
index 4fa7816460123f0534f5a369e8610079287020ad..5c9e4aacc824b06499b65b988883d9bc1f3cbb00 100644
--- a/corsika/detail/modules/radio/ZHS.inl
+++ b/corsika/detail/modules/radio/ZHS.inl
@@ -19,9 +19,9 @@ namespace corsika {
auto const startTime{step.getTimePre()};
auto const endTime{step.getTimePost()};
- // LCOV_EXCL_START
+ // LCOV_EXCL_START
if (startTime == endTime) {
- CORSIKA_LOG_ERROR("Time at the start and end of the track coincides! - radio");
+ CORSIKA_LOG_ERROR("Time at the start and end of the track coincides! - radio");
return ProcessReturn::Ok;
// LCOV_EXCL_STOP
} else {
@@ -95,9 +95,11 @@ namespace corsika {
} // end if statement for time structure
double const startBin{std::floor(
- (detectionTime1 - antenna.getStartTime()) * antenna.getSampleRate() + 0.5)};
- double const endBin{std::floor(
- (detectionTime2 - antenna.getStartTime()) * antenna.getSampleRate() + 0.5)};
+ (detectionTime1 - antenna.getStartTime()) * antenna.getSampleRate() +
+ 0.5)};
+ double const endBin{std::floor((detectionTime2 - antenna.getStartTime()) *
+ antenna.getSampleRate() +
+ 0.5)};
auto const betaPerp{
newMidPaths[k].emit_.cross(beta.cross(newMidPaths[k].emit_))};
@@ -132,9 +134,9 @@ namespace corsika {
// intermidiate contributions
for (int it{1}; it < numberOfBins; ++it) {
Vp = betaPerp * constants / denominator / newMidPaths[k].R_distance_;
- antenna.receive(
- detectionTime1 + static_cast<double>(it) / antenna.getSampleRate(),
- betaPerp, Vp);
+ antenna.receive(detectionTime1 +
+ static_cast<double>(it) / antenna.getSampleRate(),
+ betaPerp, Vp);
} // end loop over bins in which potential vector is not zero
// final contribution// f +0.5 from new antenna rounding
f = std::fabs((detectionTime2 - antenna.getStartTime()) *
@@ -173,10 +175,12 @@ namespace corsika {
sign = -1.;
} // end if statement for time structure
- double const startBin{std::floor(
- (detectionTime1 - antenna.getStartTime()) * antenna.getSampleRate() + 0.5)};
- double const endBin{std::floor(
- (detectionTime2 - antenna.getStartTime()) * antenna.getSampleRate() + 0.5)};
+ double const startBin{std::floor((detectionTime1 - antenna.getStartTime()) *
+ antenna.getSampleRate() +
+ 0.5)};
+ double const endBin{std::floor((detectionTime2 - antenna.getStartTime()) *
+ antenna.getSampleRate() +
+ 0.5)};
auto const betaPerp{midPaths[i].emit_.cross(beta.cross(midPaths[i].emit_))};
double const denominator{1. -
diff --git a/corsika/detail/modules/radio/antennas/Antenna.inl b/corsika/detail/modules/radio/antennas/Antenna.inl
index 26310172fb3fdd7b7a2246da29cf0eceaf9bb72e..0b97560aaefc129779fdf1e39dfdfe8f68e9e61b 100644
--- a/corsika/detail/modules/radio/antennas/Antenna.inl
+++ b/corsika/detail/modules/radio/antennas/Antenna.inl
@@ -12,10 +12,11 @@
namespace corsika {
template <typename TAntennaImpl>
- inline Antenna<TAntennaImpl>::Antenna(std::string const& name, Point const& location, CoordinateSystemPtr const& coordinateSystem)
+ inline Antenna<TAntennaImpl>::Antenna(std::string const& name, Point const& location,
+ CoordinateSystemPtr const& coordinateSystem)
: name_(name)
, location_(location)
- , coordinateSystem_(coordinateSystem) {};
+ , coordinateSystem_(coordinateSystem){};
template <typename TAntennaImpl>
inline Point const& Antenna<TAntennaImpl>::getLocation() const {
@@ -97,9 +98,12 @@ namespace corsika {
} 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");
+ 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");
}
}
diff --git a/corsika/detail/modules/radio/antennas/TimeDomainAntenna.inl b/corsika/detail/modules/radio/antennas/TimeDomainAntenna.inl
index 3416caf3cc1214adb88485c8f7d99fe07d0f98fb..920b303b9d63281e708ad59b278b8089a4a98cf8 100644
--- a/corsika/detail/modules/radio/antennas/TimeDomainAntenna.inl
+++ b/corsika/detail/modules/radio/antennas/TimeDomainAntenna.inl
@@ -27,7 +27,7 @@ namespace corsika {
, waveformEX_(num_bins_, 0)
, waveformEY_(num_bins_, 0)
, waveformEZ_(num_bins_, 0)
- , time_axis_(createTimeAxis()) {};
+ , time_axis_(createTimeAxis()){};
inline void TimeDomainAntenna::receive(const TimeType time,
const Vector<dimensionless_d>& receive_vector,
@@ -46,7 +46,8 @@ namespace corsika {
waveformEX_.at(timebin_) += (Electric_field_components.getX() * (1_m / 1_V));
waveformEY_.at(timebin_) += (Electric_field_components.getY() * (1_m / 1_V));
waveformEZ_.at(timebin_) += (Electric_field_components.getZ() * (1_m / 1_V));
- // TODO: Check how they are stored in memory, row-wise or column-wise? Probably use a 3D object
+ // TODO: Check how they are stored in memory, row-wise or column-wise? Probably use
+ // a 3D object
}
}
@@ -64,10 +65,14 @@ namespace corsika {
// store the x,y,z electric field components.
auto const& Vector_potential_components{vectorP.getComponents(coordinateSystem_)};
- waveformEX_.at(timebin_) += (Vector_potential_components.getX() * (1_m / (1_V * 1_s)));
- waveformEY_.at(timebin_) += (Vector_potential_components.getY() * (1_m / (1_V * 1_s)));
- waveformEZ_.at(timebin_) += (Vector_potential_components.getZ() * (1_m / (1_V * 1_s)));
- // TODO: Check how they are stored in memory, row-wise or column-wise? Probably use a 3D object
+ waveformEX_.at(timebin_) +=
+ (Vector_potential_components.getX() * (1_m / (1_V * 1_s)));
+ waveformEY_.at(timebin_) +=
+ (Vector_potential_components.getY() * (1_m / (1_V * 1_s)));
+ waveformEZ_.at(timebin_) +=
+ (Vector_potential_components.getZ() * (1_m / (1_V * 1_s)));
+ // TODO: Check how they are stored in memory, row-wise or column-wise? Probably use
+ // a 3D object
}
}
@@ -97,7 +102,9 @@ namespace corsika {
inline auto const& TimeDomainAntenna::getAxis() const { return time_axis_; }
- inline InverseTimeType const& TimeDomainAntenna::getSampleRate() const { return sample_rate_; }
+ inline InverseTimeType const& TimeDomainAntenna::getSampleRate() const {
+ return sample_rate_;
+ }
inline TimeType const& TimeDomainAntenna::getStartTime() const { return start_time_; }
diff --git a/corsika/detail/modules/radio/detectors/AntennaCollection.inl b/corsika/detail/modules/radio/detectors/AntennaCollection.inl
index e3a790865933600fae38097df2c628da8c19b0aa..d985bb98027ccd9f89b202c8352355dccd3da35f 100644
--- a/corsika/detail/modules/radio/detectors/AntennaCollection.inl
+++ b/corsika/detail/modules/radio/detectors/AntennaCollection.inl
@@ -22,12 +22,15 @@ namespace corsika {
}
template <typename TAntennaImpl>
- inline TAntennaImpl const& AntennaCollection<TAntennaImpl>::at(std::size_t const i) const {
+ inline TAntennaImpl const& AntennaCollection<TAntennaImpl>::at(
+ std::size_t const i) const {
return antennas_.at(i);
}
template <typename TAntennaImpl>
- inline int AntennaCollection<TAntennaImpl>::size() const { return antennas_.size(); }
+ inline int AntennaCollection<TAntennaImpl>::size() const {
+ return antennas_.size();
+ }
template <typename TAntennaImpl>
inline std::vector<TAntennaImpl>& AntennaCollection<TAntennaImpl>::getAntennas() {
diff --git a/corsika/modules/TimeCut.hpp b/corsika/modules/TimeCut.hpp
index c17ef6f1b538fe1af99d727b64bdbc948ca5f8ea..2c846ad0ca5dbb92266321469a46b79789ec79c7 100644
--- a/corsika/modules/TimeCut.hpp
+++ b/corsika/modules/TimeCut.hpp
@@ -27,7 +27,8 @@ namespace corsika {
TimeCut(TimeType const time);
template <typename Particle>
- ProcessReturn doContinuous(Step<Particle> const& step,
+ ProcessReturn doContinuous(
+ Step<Particle> const& step,
const bool limitFlag = false); // this is not used for TimeCut
template <typename Particle, typename Track>
LengthType getMaxStepLength(Particle const&, Track const&) {
diff --git a/corsika/modules/radio/antennas/Antenna.hpp b/corsika/modules/radio/antennas/Antenna.hpp
index 226c37eaa9c47ef788e1768c5665c80958586067..b5aeb8fb99be1dd69e0ce57536a2e0e731921d09 100644
--- a/corsika/modules/radio/antennas/Antenna.hpp
+++ b/corsika/modules/radio/antennas/Antenna.hpp
@@ -25,8 +25,8 @@ namespace corsika {
class Antenna {
protected:
- std::string const name_; ///< The name/identifier of this antenna.
- Point const location_; ///< The location of this antenna.
+ 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.
@@ -40,7 +40,8 @@ namespace corsika {
* @param location The location of this antenna.
*
*/
- Antenna(std::string const& name, Point const& location, CoordinateSystemPtr const& coordinateSystem);
+ Antenna(std::string const& name, Point const& location,
+ CoordinateSystemPtr const& coordinateSystem);
/**
* Receive a signal at this antenna.
diff --git a/corsika/modules/radio/antennas/TimeDomainAntenna.hpp b/corsika/modules/radio/antennas/TimeDomainAntenna.hpp
index a1699b9233826e0d8fadf11a7171e330d39569c2..801c529b87a035a021ab1c3aeb52cd86adf9e14d 100644
--- a/corsika/modules/radio/antennas/TimeDomainAntenna.hpp
+++ b/corsika/modules/radio/antennas/TimeDomainAntenna.hpp
@@ -19,16 +19,16 @@ namespace corsika {
*/
class TimeDomainAntenna : public Antenna<TimeDomainAntenna> {
-
- TimeType const start_time_; ///< The start time of this waveform.
- TimeType const duration_; ///< The duration of this waveform.
- InverseTimeType const sample_rate_; ///< The sampling rate of this antenna.
- int const num_bins_; ///< The number of bins used.
- std::vector<double> waveformEX_; ///< EX polarization.
- std::vector<double> waveformEY_; ///< EY polarization.
- std::vector<double> waveformEZ_; ///< EZ polarization.
- TimeType const ground_hit_time_; ///< The time the primary particle hits the ground.
- std::vector<long double> const time_axis_; ///< The time axis corresponding to the electric field.
+ TimeType const start_time_; ///< The start time of this waveform.
+ TimeType const duration_; ///< The duration of this waveform.
+ InverseTimeType const sample_rate_; ///< The sampling rate of this antenna.
+ int const num_bins_; ///< The number of bins used.
+ std::vector<double> waveformEX_; ///< EX polarization.
+ std::vector<double> waveformEY_; ///< EY polarization.
+ std::vector<double> waveformEZ_; ///< EZ polarization.
+ TimeType const ground_hit_time_; ///< The time the primary particle hits the ground.
+ std::vector<long double> const
+ time_axis_; ///< The time axis corresponding to the electric field.
public:
// import the methods from the antenna
@@ -48,12 +48,13 @@ namespace corsika {
* @param start_time The starting time of this waveform.
* @param duration The duration of this waveform.
* @param sample_rate The sample rate of this waveform.
- * @param ground_hit_time The time the primary particle hits the ground on a straight vertical line.
+ * @param ground_hit_time The time the primary particle hits the ground on a
+ * straight vertical line.
*
*/
- TimeDomainAntenna(std::string const& name, Point const& location, CoordinateSystemPtr coordinateSystem,
- TimeType const& start_time, TimeType const& duration,
- InverseTimeType const& sample_rate,
+ TimeDomainAntenna(std::string const& name, Point const& location,
+ CoordinateSystemPtr coordinateSystem, TimeType const& start_time,
+ TimeType const& duration, InverseTimeType const& sample_rate,
TimeType const ground_hit_time);
/**
diff --git a/corsika/modules/radio/propagators/SignalPath.hpp b/corsika/modules/radio/propagators/SignalPath.hpp
index 59ce63e0d616c406b73e8b919f22d673405f52dd..9a1239223838b5c56dd49384d9152bb42ee9032d 100644
--- a/corsika/modules/radio/propagators/SignalPath.hpp
+++ b/corsika/modules/radio/propagators/SignalPath.hpp
@@ -25,14 +25,14 @@ namespace corsika {
double const average_refractive_index_; ///< The average refractive index.
double const refractive_index_source_; ///< The refractive index at the source.
double const
- refractive_index_destination_; ///< The refractive index at the destination point.
+ refractive_index_destination_; ///< The refractive index at the destination point.
Vector<dimensionless_d> const emit_; ///< The (unit-length) emission vector.
Vector<dimensionless_d> const receive_; ///< The (unit-length) receive vector.
std::deque<Point> const
- points_; ///< A collection of points that make up the geometrical path.
+ points_; ///< A collection of points that make up the geometrical path.
LengthType const
- R_distance_; ///< The distance from the point of emission to an observer. TODO:
- ///< optical path, not geometrical! (probably)
+ R_distance_; ///< The distance from the point of emission to an observer. TODO:
+ ///< optical path, not geometrical! (probably)
/**
* Create a new SignalPath instance.
@@ -40,8 +40,9 @@ namespace corsika {
SignalPath(TimeType const propagation_time, double const average_refractive_index,
double const refractive_index_source,
double const refractive_index_destination,
- Vector<dimensionless_d> const& emit, Vector<dimensionless_d> const& receive,
- LengthType const R_distance, std::deque<Point> const& points);
+ Vector<dimensionless_d> const& emit,
+ Vector<dimensionless_d> const& receive, LengthType const R_distance,
+ std::deque<Point> const& points);
}; // END: class SignalPath final
diff --git a/examples/clover_leaf.cpp b/examples/clover_leaf.cpp
index cb44e9d876ef7fd5112c26900c428eb491d4b05a..ccb8c7e79aed238d9bfeabe5994dac0efe41cbff 100644
--- a/examples/clover_leaf.cpp
+++ b/examples/clover_leaf.cpp
@@ -39,8 +39,7 @@
#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/TrackWriter.hpp>
#include <corsika/modules/StackInspector.hpp>
#include <corsika/modules/ParticleCut.hpp>
@@ -73,150 +72,157 @@ using namespace std;
//
int main() {
- logging::set_level(logging::level::warn);
- corsika_logger->set_pattern("[%n:%^%-8l%$] custom pattern: %v");
-
- CORSIKA_LOG_INFO("Synchrotron radiation");
-
- feenableexcept(FE_INVALID);
- RNGManager<>::getInstance().registerRandomStream("cascade");
- std::random_device rd;
- auto seed = rd();
- RNGManager<>::getInstance().setSeed(seed);
-
- OutputManager output("1 electron - 1 positron");
-
- // create a suitable environment
- using IModelInterface =
- IRefractiveIndexModel<IMediumPropertyModel<IMagneticFieldModel<IMediumModel>>>;
- using AtmModel = UniformRefractiveIndex<
- MediumPropertyModel<UniformMagneticField<HomogeneousMedium<IModelInterface>>>>;
- using EnvType = Environment<AtmModel>;
- EnvType env;
- CoordinateSystemPtr const& rootCS = env.getCoordinateSystem();
- Point const center{rootCS, 0_m, 0_m, 0_m};
- // a refractive index for the vacuum
- const double ri_{1.};
- // the composition we use for the homogeneous medium
- NuclearComposition const Composition({Code::Nitrogen}, {1.});
- // create magnetic field vector
- auto const Bmag {0.00005_T};
- MagneticFieldVector B(rootCS, 0_T, Bmag, 0_T);
- // create a Sphere for the medium
- auto world = EnvType::createNode<Sphere>(center, 150_km);
- // set the environment properties
- world->setModelProperties<AtmModel>(ri_, Medium::AirDry1Atm, B,
- 1_kg / (1_m * 1_m * 1_m), Composition);
- // bind things together
- env.getUniverse()->addChild(std::move(world));
-
-
- Point injectionPos(rootCS, 0_m, 0_m, 5_km);
-
- auto const centerX{center.getCoordinates().getX()};
- auto const centerY{center.getCoordinates().getY()};
- auto const centerZ{center.getCoordinates().getZ()};
- auto const injectionPosX_ {injectionPos.getCoordinates().getX()};
- auto const injectionPosY_ {injectionPos.getCoordinates().getY()};
- auto const injectionPosZ_ {injectionPos.getCoordinates().getZ()};
- auto const triggerpoint_ {Point(rootCS, injectionPosX_, injectionPosY_, injectionPosZ_)};
- std::cout << "Trigger Point is: " << triggerpoint_ << std::endl;
-
- // the antenna characteristics
- const TimeType duration_{2e-6_s}; //0.8e-4_s
- const InverseTimeType sampleRate_{1e+11_Hz}; // 1e+9_Hz
-
- // the detectors
- AntennaCollection<TimeDomainAntenna> detectorCoREAS;
- // AntennaCollection<TimeDomainAntenna> detectorZHS;
-
- std::string name_center = "CoREAS_R=0_m--Phi=0degrees";
- auto triggertime_center {((triggerpoint_ - center).getNorm() / constants::c) - 500_ns};
- TimeDomainAntenna antenna_center(name_center, center, rootCS, triggertime_center, duration_, sampleRate_, triggertime_center);
- detectorCoREAS.addAntenna(antenna_center);
-
- for (auto radius_1 = 25_m; radius_1 <= 1000_m; radius_1 += 25_m) {
- for (auto phi_1 = 0; phi_1 <= 315; phi_1 += 45) {
- auto phiRad_1 = phi_1 / 180. * M_PI;
- auto rr_1 = static_cast<int>(radius_1 / 1_m);
- auto const point_1 {Point(rootCS, centerX + radius_1 * cos(phiRad_1), centerY + radius_1 * sin(phiRad_1), centerZ)};
- std::cout << "Antenna point: " << point_1 << std::endl;
- auto triggertime_1 {((triggerpoint_ - point_1).getNorm() / constants::c) - 500_ns};
- std::string name_1 = "CoREAS_R=" + std::to_string(rr_1) + "_m--Phi=" + std::to_string(phi_1) + "degrees";
- TimeDomainAntenna antenna_1(name_1, point_1, rootCS, triggertime_1, duration_, sampleRate_, triggertime_1);
- detectorCoREAS.addAntenna(antenna_1);
- }
+ logging::set_level(logging::level::warn);
+ corsika_logger->set_pattern("[%n:%^%-8l%$] custom pattern: %v");
+
+ CORSIKA_LOG_INFO("Synchrotron radiation");
+
+ feenableexcept(FE_INVALID);
+ RNGManager<>::getInstance().registerRandomStream("cascade");
+ std::random_device rd;
+ auto seed = rd();
+ RNGManager<>::getInstance().setSeed(seed);
+
+ OutputManager output("1 electron - 1 positron");
+
+ // create a suitable environment
+ using IModelInterface =
+ IRefractiveIndexModel<IMediumPropertyModel<IMagneticFieldModel<IMediumModel>>>;
+ using AtmModel = UniformRefractiveIndex<
+ MediumPropertyModel<UniformMagneticField<HomogeneousMedium<IModelInterface>>>>;
+ using EnvType = Environment<AtmModel>;
+ EnvType env;
+ CoordinateSystemPtr const& rootCS = env.getCoordinateSystem();
+ Point const center{rootCS, 0_m, 0_m, 0_m};
+ // a refractive index for the vacuum
+ const double ri_{1.};
+ // the composition we use for the homogeneous medium
+ NuclearComposition const Composition({Code::Nitrogen}, {1.});
+ // create magnetic field vector
+ auto const Bmag{0.00005_T};
+ MagneticFieldVector B(rootCS, 0_T, Bmag, 0_T);
+ // create a Sphere for the medium
+ auto world = EnvType::createNode<Sphere>(center, 150_km);
+ // set the environment properties
+ world->setModelProperties<AtmModel>(ri_, Medium::AirDry1Atm, B,
+ 1_kg / (1_m * 1_m * 1_m), Composition);
+ // bind things together
+ env.getUniverse()->addChild(std::move(world));
+
+ Point injectionPos(rootCS, 0_m, 0_m, 5_km);
+
+ auto const centerX{center.getCoordinates().getX()};
+ auto const centerY{center.getCoordinates().getY()};
+ auto const centerZ{center.getCoordinates().getZ()};
+ auto const injectionPosX_{injectionPos.getCoordinates().getX()};
+ auto const injectionPosY_{injectionPos.getCoordinates().getY()};
+ auto const injectionPosZ_{injectionPos.getCoordinates().getZ()};
+ auto const triggerpoint_{Point(rootCS, injectionPosX_, injectionPosY_, injectionPosZ_)};
+ std::cout << "Trigger Point is: " << triggerpoint_ << std::endl;
+
+ // the antenna characteristics
+ const TimeType duration_{2e-6_s}; // 0.8e-4_s
+ const InverseTimeType sampleRate_{1e+11_Hz}; // 1e+9_Hz
+
+ // the detectors
+ AntennaCollection<TimeDomainAntenna> detectorCoREAS;
+ // AntennaCollection<TimeDomainAntenna> detectorZHS;
+
+ std::string name_center = "CoREAS_R=0_m--Phi=0degrees";
+ auto triggertime_center{((triggerpoint_ - center).getNorm() / constants::c) - 500_ns};
+ TimeDomainAntenna antenna_center(name_center, center, rootCS, triggertime_center,
+ duration_, sampleRate_, triggertime_center);
+ detectorCoREAS.addAntenna(antenna_center);
+
+ for (auto radius_1 = 25_m; radius_1 <= 1000_m; radius_1 += 25_m) {
+ for (auto phi_1 = 0; phi_1 <= 315; phi_1 += 45) {
+ auto phiRad_1 = phi_1 / 180. * M_PI;
+ auto rr_1 = static_cast<int>(radius_1 / 1_m);
+ auto const point_1{Point(rootCS, centerX + radius_1 * cos(phiRad_1),
+ centerY + radius_1 * sin(phiRad_1), centerZ)};
+ std::cout << "Antenna point: " << point_1 << std::endl;
+ auto triggertime_1{((triggerpoint_ - point_1).getNorm() / constants::c) - 500_ns};
+ std::string name_1 = "CoREAS_R=" + std::to_string(rr_1) +
+ "_m--Phi=" + std::to_string(phi_1) + "degrees";
+ TimeDomainAntenna antenna_1(name_1, point_1, rootCS, triggertime_1, duration_,
+ sampleRate_, triggertime_1);
+ detectorCoREAS.addAntenna(antenna_1);
}
-
-
- // setup particle stack, and add primary particle
- setup::Stack<EnvType> stack;
- stack.clear();
-
- // electron
- const Code beamCode = Code::Electron;
- auto const charge = get_charge(beamCode);
- auto const mass = get_mass(beamCode);
- auto const gyroradius = 100_m;
- auto const pLabMag = convert_SI_to_HEP(charge * Bmag * gyroradius);
- auto const omega_inv = convert_HEP_to_SI<MassType::dimension_type>(mass) / (abs(charge) * Bmag);
- MomentumVector const plab{rootCS, 0_MeV, 0_MeV, -10_MeV};
- auto const Elab = sqrt(plab.getSquaredNorm() + static_pow<2>(mass));
- auto gamma = Elab / mass;
- TimeType const period = 2 * M_PI * omega_inv * gamma;
-
- // positron
- const Code beamCodeP = Code::Positron;
- auto const chargeP = get_charge(beamCodeP);
- auto const massP = get_mass(beamCodeP);
- // auto const gyroradius = 100_m;
- auto const pLabMagP = convert_SI_to_HEP(chargeP * Bmag * gyroradius);
- auto const omega_invP = convert_HEP_to_SI<MassType::dimension_type>(massP) / (abs(chargeP) * Bmag);
- MomentumVector const plabP{rootCS, 0_MeV, 0_MeV, -10_MeV};
- auto const ElabP = sqrt(plabP.getSquaredNorm() + static_pow<2>(massP));
- auto gammaP = ElabP / massP;
- TimeType const periodP = 2 * M_PI * omega_invP * gammaP;
-
-
-
- stack.addParticle(std::make_tuple(beamCode, calculate_kinetic_energy(plab.getNorm(), mass), plab.normalized(), injectionPos, 0_ns));
- stack.addParticle(std::make_tuple(beamCodeP, calculate_kinetic_energy(plabP.getNorm(), massP), plabP.normalized(), injectionPos, 0_ns));
-
- // setup relevant processes
- setup::Tracking tracking;
-
- // put radio processes here
- RadioProcess<decltype(detectorCoREAS), CoREAS<decltype(detectorCoREAS),
- decltype(SimplePropagator(env))>, decltype(SimplePropagator(env))>
- coreas(detectorCoREAS, env);
- output.add("CoREAS", coreas);
-
- // RadioProcess<decltype(detectorZHS), ZHS<decltype(detectorZHS),
- // decltype(SimplePropagator(env))>, decltype(SimplePropagator(env))>
- // zhs(detectorZHS, env);
- // output.add("ZHS", zhs);
-
- TimeCut cut(period / 4);
-
- // TrackWriter trackWriter;
- // output.add("tracks", trackWriter); // register TrackWriter
-
- // assemble all processes into an ordered process list
- auto sequence = make_sequence(coreas, cut);
-
- // 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);
- CORSIKA_LOG_INFO("gamma: {}", gamma);
-
- CORSIKA_LOG_INFO("|p| positron = {} and E positron = {}",plabP.getNorm(), ElabP);
- CORSIKA_LOG_INFO("period: {}", periodP);
- CORSIKA_LOG_INFO("gamma: {}", gammaP);
-
- output.endOfLibrary();
+ }
+
+ // setup particle stack, and add primary particle
+ setup::Stack<EnvType> stack;
+ stack.clear();
+
+ // electron
+ const Code beamCode = Code::Electron;
+ auto const charge = get_charge(beamCode);
+ auto const mass = get_mass(beamCode);
+ auto const gyroradius = 100_m;
+ auto const pLabMag = convert_SI_to_HEP(charge * Bmag * gyroradius);
+ auto const omega_inv =
+ convert_HEP_to_SI<MassType::dimension_type>(mass) / (abs(charge) * Bmag);
+ MomentumVector const plab{rootCS, 0_MeV, 0_MeV, -10_MeV};
+ auto const Elab = sqrt(plab.getSquaredNorm() + static_pow<2>(mass));
+ auto gamma = Elab / mass;
+ TimeType const period = 2 * M_PI * omega_inv * gamma;
+
+ // positron
+ const Code beamCodeP = Code::Positron;
+ auto const chargeP = get_charge(beamCodeP);
+ auto const massP = get_mass(beamCodeP);
+ // auto const gyroradius = 100_m;
+ auto const pLabMagP = convert_SI_to_HEP(chargeP * Bmag * gyroradius);
+ auto const omega_invP =
+ convert_HEP_to_SI<MassType::dimension_type>(massP) / (abs(chargeP) * Bmag);
+ MomentumVector const plabP{rootCS, 0_MeV, 0_MeV, -10_MeV};
+ auto const ElabP = sqrt(plabP.getSquaredNorm() + static_pow<2>(massP));
+ auto gammaP = ElabP / massP;
+ TimeType const periodP = 2 * M_PI * omega_invP * gammaP;
+
+ stack.addParticle(std::make_tuple(beamCode,
+ calculate_kinetic_energy(plab.getNorm(), mass),
+ plab.normalized(), injectionPos, 0_ns));
+ stack.addParticle(std::make_tuple(beamCodeP,
+ calculate_kinetic_energy(plabP.getNorm(), massP),
+ plabP.normalized(), injectionPos, 0_ns));
+
+ // setup relevant processes
+ setup::Tracking tracking;
+
+ // put radio processes here
+ RadioProcess<decltype(detectorCoREAS),
+ CoREAS<decltype(detectorCoREAS), decltype(SimplePropagator(env))>,
+ decltype(SimplePropagator(env))>
+ coreas(detectorCoREAS, env);
+ output.add("CoREAS", coreas);
+
+ // RadioProcess<decltype(detectorZHS), ZHS<decltype(detectorZHS),
+ // decltype(SimplePropagator(env))>, decltype(SimplePropagator(env))>
+ // zhs(detectorZHS, env);
+ // output.add("ZHS", zhs);
+
+ TimeCut cut(period / 4);
+
+ // TrackWriter trackWriter;
+ // output.add("tracks", trackWriter); // register TrackWriter
+
+ // assemble all processes into an ordered process list
+ auto sequence = make_sequence(coreas, cut);
+
+ // 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);
+ CORSIKA_LOG_INFO("gamma: {}", gamma);
+
+ CORSIKA_LOG_INFO("|p| positron = {} and E positron = {}", plabP.getNorm(), ElabP);
+ CORSIKA_LOG_INFO("period: {}", periodP);
+ CORSIKA_LOG_INFO("gamma: {}", gammaP);
+
+ output.endOfLibrary();
}
diff --git a/examples/radio_em_shower.cpp b/examples/radio_em_shower.cpp
index 39147fffa4381e0812339e7de442a19aeeec94d9..f151b1e893d0444c47b54365f6145852f53dc0ff 100644
--- a/examples/radio_em_shower.cpp
+++ b/examples/radio_em_shower.cpp
@@ -1,10 +1,10 @@
/*
-* (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
-* the license.
-*/
+ * (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
+ * the license.
+ */
#include <corsika/framework/process/ProcessSequence.hpp>
#include <corsika/framework/process/SwitchProcessSequence.hpp>
@@ -76,230 +76,230 @@ using namespace corsika;
using namespace std;
void registerRandomStreams(int seed) {
- RNGManager<>::getInstance().registerRandomStream("cascade");
- RNGManager<>::getInstance().registerRandomStream("proposal");
- RNGManager<>::getInstance().registerRandomStream("sibyll");
- if (seed == 0) {
- std::random_device rd;
- seed = rd();
- cout << "new random seed (auto) " << seed << endl;
- }
- RNGManager<>::getInstance().setSeed(seed);
+ RNGManager<>::getInstance().registerRandomStream("cascade");
+ RNGManager<>::getInstance().registerRandomStream("proposal");
+ RNGManager<>::getInstance().registerRandomStream("sibyll");
+ if (seed == 0) {
+ std::random_device rd;
+ seed = rd();
+ cout << "new random seed (auto) " << seed << endl;
+ }
+ RNGManager<>::getInstance().setSeed(seed);
}
template <typename TInterface>
using MyExtraEnv =
- UniformRefractiveIndex<MediumPropertyModel<UniformMagneticField<TInterface>>>;
+ UniformRefractiveIndex<MediumPropertyModel<UniformMagneticField<TInterface>>>;
int main(int argc, char** argv) {
- logging::set_level(logging::level::warn);
-
- if (argc != 3) {
- std::cerr
- << "usage: radio_em_shower <energy/GeV> <seed> - put seed=0 to use random seed"
- << std::endl;
- return 1;
- }
-
- int seed{static_cast<int>(std::stof(std::string(argv[2])))};
- std::cout << "Seed: " << seed << std::endl;
- feenableexcept(FE_INVALID);
- // initialize random number sequence(s)
- registerRandomStreams(seed);
-
- // setup environment, geometry
- using EnvironmentInterface =
- IRefractiveIndexModel<IMediumPropertyModel<IMagneticFieldModel<IMediumModel>>>;
- using EnvType = Environment<EnvironmentInterface>;
- EnvType env;
- CoordinateSystemPtr const& rootCS = env.getCoordinateSystem();
- Point const center{rootCS, 0_m, 0_m, 0_m};
-
- create_5layer_atmosphere<EnvironmentInterface, MyExtraEnv>(
- env, AtmosphereId::LinsleyUSStd, center, 1.000327, Medium::AirDry1Atm,
- MagneticFieldVector{rootCS, 50_uT, 0_T, 0_T});
-
- std::unordered_map<Code, HEPEnergyType> energy_resolution = {
- {Code::Electron, 5_MeV},
- {Code::Positron, 5_MeV},
- {Code::Photon, 5_MeV},
- };
- for (auto [pcode, energy] : energy_resolution)
- set_energy_production_threshold(pcode, energy);
-
- // setup particle stack, and add primary particle
- setup::Stack<EnvType> stack;
- stack.clear();
- const Code beamCode = Code::Electron;
- auto const mass = get_mass(beamCode);
- const HEPEnergyType E0 = 1_GeV * std::stof(std::string(argv[1]));
- double theta = 0.;
- auto const thetaRad = theta / 180. * M_PI;
-
- HEPMomentumType P0 = calculate_momentum(E0, mass);
- auto momentumComponents = [](double thetaRad, HEPMomentumType ptot) {
- return std::make_tuple(ptot * sin(thetaRad), 0_eV, -ptot * cos(thetaRad));
- };
-
- auto const [px, py, pz] = momentumComponents(thetaRad, P0);
- auto plab = MomentumVector(rootCS, {px, py, pz});
- cout << "input particle: " << beamCode << endl;
- cout << "input angles: theta=" << theta << endl;
- cout << "input momentum: " << plab.getComponents() / 1_GeV
- << ", norm = " << plab.getNorm() << endl;
-
- auto const observationHeight = 1.4_km + constants::EarthRadius::Mean;
- auto const injectionHeight = 112.75_km + constants::EarthRadius::Mean;
- auto const t = -observationHeight * cos(thetaRad) +
- sqrt(-static_pow<2>(sin(thetaRad) * observationHeight) +
- static_pow<2>(injectionHeight));
- Point const showerCore{rootCS, 0_m, 0_m, observationHeight};
- Point const injectionPos =
- showerCore + DirectionVector{rootCS, {-sin(thetaRad), 0, cos(thetaRad)}} * t;
-
- std::cout << "point of injection: " << injectionPos.getCoordinates() << std::endl;
-
- stack.addParticle(std::make_tuple(
- beamCode, calculate_kinetic_energy(plab.getNorm(), get_mass(beamCode)),
- plab.normalized(), injectionPos, 0_ns));
-
- CORSIKA_LOG_INFO("shower axis length: {} ",
- (showerCore - injectionPos).getNorm() * 1.02);
-
- ShowerAxis const showerAxis{injectionPos, (showerCore - injectionPos) * 1.02, env,
- false, 1000};
-
- TimeType const groundHitTime{(showerCore - injectionPos).getNorm() / constants::c};
-
- std::string outname_ = "radio_em_shower_outputs"; // + std::to_string(rr_);
- OutputManager output(outname_);
-
- // Radio antennas and relevant information
- // the antenna time variables
- const TimeType duration_{1e-6_s};
- const InverseTimeType sampleRate_{1e+9_Hz};
-
- // the detector (aka antenna collection) for CoREAS and ZHS
- AntennaCollection<TimeDomainAntenna> detectorCoREAS;
- AntennaCollection<TimeDomainAntenna> detectorZHS;
-
- auto const showerCoreX_{showerCore.getCoordinates().getX()};
- auto const showerCoreY_{showerCore.getCoordinates().getY()};
- auto const injectionPosX_{injectionPos.getCoordinates().getX()};
- auto const injectionPosY_{injectionPos.getCoordinates().getY()};
- auto const injectionPosZ_{injectionPos.getCoordinates().getZ()};
- auto const triggerpoint_{Point(rootCS, injectionPosX_, injectionPosY_, injectionPosZ_)};
- std::cout << "Trigger Point is: " << triggerpoint_ << std::endl;
-
- // // setup CoREAS antennas - use the for loop for star shape pattern
-// for (auto radius_1 = 25_m; radius_1 <= 500_m; radius_1 += 25_m) {
-// for (auto phi_1 = 0; phi_1 <= 315; phi_1 += 45) {
- auto radius_1 = 200_m;
- auto phi_1 = 45;
- auto phiRad_1 = phi_1 / 180. * M_PI;
- auto rr_1 = static_cast<int>(radius_1 / 1_m);
- auto const point_1{Point(rootCS, showerCoreX_ + radius_1 * cos(phiRad_1),
- showerCoreY_ + radius_1 * sin(phiRad_1),
- constants::EarthRadius::Mean)};
- std::cout << "Antenna point: " << point_1 << std::endl;
- auto triggertime_1{(triggerpoint_ - point_1).getNorm() / constants::c};
- std::string name_1 = "CoREAS_R=" + std::to_string(rr_1) +
- "_m--Phi=" + std::to_string(phi_1) + "degrees";
- TimeDomainAntenna antenna_1(name_1, point_1, rootCS, triggertime_1, duration_, sampleRate_,
- triggertime_1);
- detectorCoREAS.addAntenna(antenna_1);
-// }
-// }
-
- // // setup ZHS antennas - use the for loop for star shape pattern
-// for (auto radius_ = 25_m; radius_ <= 500_m; radius_ += 25_m) {
-// for (auto phi_ = 0; phi_ <= 315; phi_ += 45) {
- auto radius_ = 200_m;
- auto phi_ = 45;
- auto phiRad_ = phi_ / 180. * M_PI;
- auto rr_ = static_cast<int>(radius_ / 1_m);
- auto const point_{Point(rootCS, showerCoreX_ + radius_ * cos(phiRad_),
- showerCoreY_ + radius_ * sin(phiRad_),
- constants::EarthRadius::Mean)};
- auto triggertime_{(triggerpoint_ - point_).getNorm() / constants::c};
- std::string name_ =
- "ZHS_R=" + std::to_string(rr_) + "_m--Phi=" + std::to_string(phi_) + "degrees";
- TimeDomainAntenna antenna_(name_, point_, rootCS, triggertime_, duration_, sampleRate_,
- triggertime_);
- detectorZHS.addAntenna(antenna_);
-// }
-// }
-
- // setup processes, decays and interactions
-
- EnergyLossWriter dEdX{showerAxis, 10_g / square(1_cm), 200};
- // register energy losses as output
- output.add("dEdX", dEdX);
-
- ParticleCut<SubWriter<decltype(dEdX)>> cut(5_MeV, 5_MeV, 100_GeV, 100_GeV, true, dEdX);
-
- corsika::sibyll::Interaction sibyll{env};
- HEPEnergyType heThresholdNN = 80_GeV;
- corsika::proposal::Interaction emCascade(env, sibyll.getHadronInteractionModel(),
- heThresholdNN);
- corsika::proposal::ContinuousProcess<SubWriter<decltype(dEdX)>> emContinuous(env, dEdX);
- // BetheBlochPDG<SubWriter<decltype(dEdX)>> emContinuous{dEdX};
-
- // NOT possible right now, due to interface differenc in PROPOSAL
- // InteractionCounter emCascadeCounted(emCascade);
-
- TrackWriter tracks;
- output.add("tracks", tracks);
-
- // long. profile
- LongitudinalWriter profile{showerAxis, 10_g / square(1_cm)};
- output.add("profile", profile);
- LongitudinalProfile<SubWriter<decltype(profile)>> longprof{profile};
-
- // initiate CoREAS
- RadioProcess<decltype(detectorCoREAS),
- CoREAS<decltype(detectorCoREAS), decltype(SimplePropagator(env))>,
- decltype(SimplePropagator(env))>
- coreas(detectorCoREAS, env);
-
- // register CoREAS with the output manager
- output.add("CoREAS", coreas);
-
- // initiate ZHS
- RadioProcess<decltype(detectorZHS),
- ZHS<decltype(detectorZHS), decltype(SimplePropagator(env))>,
- decltype(SimplePropagator(env))>
- zhs(detectorZHS, env);
-
- // // register ZHS with the output manager
- output.add("ZHS", zhs);
-
- Plane const obsPlane(showerCore, DirectionVector(rootCS, {0., 0., 1.}));
- ObservationPlane<setup::Tracking, ParticleWriterParquet> observationLevel{
- obsPlane, DirectionVector(rootCS, {1., 0., 0.})};
- output.add("particles", observationLevel);
-
- // auto sequence = make_sequence(emCascade, emContinuous, longprof, cut, coreas, zhs);
- auto sequence = make_sequence(emCascade, emContinuous, longprof, cut, coreas, zhs,
- observationLevel, tracks);
- // define air shower object, run simulation
- setup::Tracking tracking;
-
- output.startOfLibrary();
- Cascade EAS(env, tracking, sequence, output, stack);
-
- // to fix the point of first interaction, uncomment the following two lines:
- // EAS.forceInteraction();
-
- EAS.run();
-
- HEPEnergyType const Efinal = dEdX.getEnergyLost() + observationLevel.getEnergyGround();
-
- CORSIKA_LOG_INFO(
- "total energy budget (GeV): {}, "
- "relative difference (%): {}",
- Efinal / 1_GeV, (Efinal / E0 - 1) * 100);
-
- output.endOfLibrary();
+ logging::set_level(logging::level::warn);
+
+ if (argc != 3) {
+ std::cerr
+ << "usage: radio_em_shower <energy/GeV> <seed> - put seed=0 to use random seed"
+ << std::endl;
+ return 1;
+ }
+
+ int seed{static_cast<int>(std::stof(std::string(argv[2])))};
+ std::cout << "Seed: " << seed << std::endl;
+ feenableexcept(FE_INVALID);
+ // initialize random number sequence(s)
+ registerRandomStreams(seed);
+
+ // setup environment, geometry
+ using EnvironmentInterface =
+ IRefractiveIndexModel<IMediumPropertyModel<IMagneticFieldModel<IMediumModel>>>;
+ using EnvType = Environment<EnvironmentInterface>;
+ EnvType env;
+ CoordinateSystemPtr const& rootCS = env.getCoordinateSystem();
+ Point const center{rootCS, 0_m, 0_m, 0_m};
+
+ create_5layer_atmosphere<EnvironmentInterface, MyExtraEnv>(
+ env, AtmosphereId::LinsleyUSStd, center, 1.000327, Medium::AirDry1Atm,
+ MagneticFieldVector{rootCS, 50_uT, 0_T, 0_T});
+
+ std::unordered_map<Code, HEPEnergyType> energy_resolution = {
+ {Code::Electron, 5_MeV},
+ {Code::Positron, 5_MeV},
+ {Code::Photon, 5_MeV},
+ };
+ for (auto [pcode, energy] : energy_resolution)
+ set_energy_production_threshold(pcode, energy);
+
+ // setup particle stack, and add primary particle
+ setup::Stack<EnvType> stack;
+ stack.clear();
+ const Code beamCode = Code::Electron;
+ auto const mass = get_mass(beamCode);
+ const HEPEnergyType E0 = 1_GeV * std::stof(std::string(argv[1]));
+ double theta = 0.;
+ auto const thetaRad = theta / 180. * M_PI;
+
+ HEPMomentumType P0 = calculate_momentum(E0, mass);
+ auto momentumComponents = [](double thetaRad, HEPMomentumType ptot) {
+ return std::make_tuple(ptot * sin(thetaRad), 0_eV, -ptot * cos(thetaRad));
+ };
+
+ auto const [px, py, pz] = momentumComponents(thetaRad, P0);
+ auto plab = MomentumVector(rootCS, {px, py, pz});
+ cout << "input particle: " << beamCode << endl;
+ cout << "input angles: theta=" << theta << endl;
+ cout << "input momentum: " << plab.getComponents() / 1_GeV
+ << ", norm = " << plab.getNorm() << endl;
+
+ auto const observationHeight = 1.4_km + constants::EarthRadius::Mean;
+ auto const injectionHeight = 112.75_km + constants::EarthRadius::Mean;
+ auto const t = -observationHeight * cos(thetaRad) +
+ sqrt(-static_pow<2>(sin(thetaRad) * observationHeight) +
+ static_pow<2>(injectionHeight));
+ Point const showerCore{rootCS, 0_m, 0_m, observationHeight};
+ Point const injectionPos =
+ showerCore + DirectionVector{rootCS, {-sin(thetaRad), 0, cos(thetaRad)}} * t;
+
+ std::cout << "point of injection: " << injectionPos.getCoordinates() << std::endl;
+
+ stack.addParticle(std::make_tuple(
+ beamCode, calculate_kinetic_energy(plab.getNorm(), get_mass(beamCode)),
+ plab.normalized(), injectionPos, 0_ns));
+
+ CORSIKA_LOG_INFO("shower axis length: {} ",
+ (showerCore - injectionPos).getNorm() * 1.02);
+
+ ShowerAxis const showerAxis{injectionPos, (showerCore - injectionPos) * 1.02, env,
+ false, 1000};
+
+ TimeType const groundHitTime{(showerCore - injectionPos).getNorm() / constants::c};
+
+ std::string outname_ = "radio_em_shower_outputs"; // + std::to_string(rr_);
+ OutputManager output(outname_);
+
+ // Radio antennas and relevant information
+ // the antenna time variables
+ const TimeType duration_{1e-6_s};
+ const InverseTimeType sampleRate_{1e+9_Hz};
+
+ // the detector (aka antenna collection) for CoREAS and ZHS
+ AntennaCollection<TimeDomainAntenna> detectorCoREAS;
+ AntennaCollection<TimeDomainAntenna> detectorZHS;
+
+ auto const showerCoreX_{showerCore.getCoordinates().getX()};
+ auto const showerCoreY_{showerCore.getCoordinates().getY()};
+ auto const injectionPosX_{injectionPos.getCoordinates().getX()};
+ auto const injectionPosY_{injectionPos.getCoordinates().getY()};
+ auto const injectionPosZ_{injectionPos.getCoordinates().getZ()};
+ auto const triggerpoint_{Point(rootCS, injectionPosX_, injectionPosY_, injectionPosZ_)};
+ std::cout << "Trigger Point is: " << triggerpoint_ << std::endl;
+
+ // // setup CoREAS antennas - use the for loop for star shape pattern
+ // for (auto radius_1 = 25_m; radius_1 <= 500_m; radius_1 += 25_m) {
+ // for (auto phi_1 = 0; phi_1 <= 315; phi_1 += 45) {
+ auto radius_1 = 200_m;
+ auto phi_1 = 45;
+ auto phiRad_1 = phi_1 / 180. * M_PI;
+ auto rr_1 = static_cast<int>(radius_1 / 1_m);
+ auto const point_1{Point(rootCS, showerCoreX_ + radius_1 * cos(phiRad_1),
+ showerCoreY_ + radius_1 * sin(phiRad_1),
+ constants::EarthRadius::Mean)};
+ std::cout << "Antenna point: " << point_1 << std::endl;
+ auto triggertime_1{(triggerpoint_ - point_1).getNorm() / constants::c};
+ std::string name_1 =
+ "CoREAS_R=" + std::to_string(rr_1) + "_m--Phi=" + std::to_string(phi_1) + "degrees";
+ TimeDomainAntenna antenna_1(name_1, point_1, rootCS, triggertime_1, duration_,
+ sampleRate_, triggertime_1);
+ detectorCoREAS.addAntenna(antenna_1);
+ // }
+ // }
+
+ // // setup ZHS antennas - use the for loop for star shape pattern
+ // for (auto radius_ = 25_m; radius_ <= 500_m; radius_ += 25_m) {
+ // for (auto phi_ = 0; phi_ <= 315; phi_ += 45) {
+ auto radius_ = 200_m;
+ auto phi_ = 45;
+ auto phiRad_ = phi_ / 180. * M_PI;
+ auto rr_ = static_cast<int>(radius_ / 1_m);
+ auto const point_{Point(rootCS, showerCoreX_ + radius_ * cos(phiRad_),
+ showerCoreY_ + radius_ * sin(phiRad_),
+ constants::EarthRadius::Mean)};
+ auto triggertime_{(triggerpoint_ - point_).getNorm() / constants::c};
+ std::string name_ =
+ "ZHS_R=" + std::to_string(rr_) + "_m--Phi=" + std::to_string(phi_) + "degrees";
+ TimeDomainAntenna antenna_(name_, point_, rootCS, triggertime_, duration_, sampleRate_,
+ triggertime_);
+ detectorZHS.addAntenna(antenna_);
+ // }
+ // }
+
+ // setup processes, decays and interactions
+
+ EnergyLossWriter dEdX{showerAxis, 10_g / square(1_cm), 200};
+ // register energy losses as output
+ output.add("dEdX", dEdX);
+
+ ParticleCut<SubWriter<decltype(dEdX)>> cut(5_MeV, 5_MeV, 100_GeV, 100_GeV, true, dEdX);
+
+ corsika::sibyll::Interaction sibyll{env};
+ HEPEnergyType heThresholdNN = 80_GeV;
+ corsika::proposal::Interaction emCascade(env, sibyll.getHadronInteractionModel(),
+ heThresholdNN);
+ corsika::proposal::ContinuousProcess<SubWriter<decltype(dEdX)>> emContinuous(env, dEdX);
+ // BetheBlochPDG<SubWriter<decltype(dEdX)>> emContinuous{dEdX};
+
+ // NOT possible right now, due to interface differenc in PROPOSAL
+ // InteractionCounter emCascadeCounted(emCascade);
+
+ TrackWriter tracks;
+ output.add("tracks", tracks);
+
+ // long. profile
+ LongitudinalWriter profile{showerAxis, 10_g / square(1_cm)};
+ output.add("profile", profile);
+ LongitudinalProfile<SubWriter<decltype(profile)>> longprof{profile};
+
+ // initiate CoREAS
+ RadioProcess<decltype(detectorCoREAS),
+ CoREAS<decltype(detectorCoREAS), decltype(SimplePropagator(env))>,
+ decltype(SimplePropagator(env))>
+ coreas(detectorCoREAS, env);
+
+ // register CoREAS with the output manager
+ output.add("CoREAS", coreas);
+
+ // initiate ZHS
+ RadioProcess<decltype(detectorZHS),
+ ZHS<decltype(detectorZHS), decltype(SimplePropagator(env))>,
+ decltype(SimplePropagator(env))>
+ zhs(detectorZHS, env);
+
+ // // register ZHS with the output manager
+ output.add("ZHS", zhs);
+
+ Plane const obsPlane(showerCore, DirectionVector(rootCS, {0., 0., 1.}));
+ ObservationPlane<setup::Tracking, ParticleWriterParquet> observationLevel{
+ obsPlane, DirectionVector(rootCS, {1., 0., 0.})};
+ output.add("particles", observationLevel);
+
+ // auto sequence = make_sequence(emCascade, emContinuous, longprof, cut, coreas, zhs);
+ auto sequence = make_sequence(emCascade, emContinuous, longprof, cut, coreas, zhs,
+ observationLevel, tracks);
+ // define air shower object, run simulation
+ setup::Tracking tracking;
+
+ output.startOfLibrary();
+ Cascade EAS(env, tracking, sequence, output, stack);
+
+ // to fix the point of first interaction, uncomment the following two lines:
+ // EAS.forceInteraction();
+
+ EAS.run();
+
+ HEPEnergyType const Efinal = dEdX.getEnergyLost() + observationLevel.getEnergyGround();
+
+ CORSIKA_LOG_INFO(
+ "total energy budget (GeV): {}, "
+ "relative difference (%): {}",
+ Efinal / 1_GeV, (Efinal / E0 - 1) * 100);
+
+ output.endOfLibrary();
}
\ No newline at end of file
diff --git a/examples/synchrotron_test_C8tracking.cpp b/examples/synchrotron_test_C8tracking.cpp
index 944888383c47296ac6d167302b69d84fb47db5d8..cf1792188e4e463c33c2f5a9078b8121bdfb9ef3 100644
--- a/examples/synchrotron_test_C8tracking.cpp
+++ b/examples/synchrotron_test_C8tracking.cpp
@@ -76,7 +76,8 @@ int main() {
OutputManager output("synchrotron_radiation_C8tracking-output");
// set up the environment
- using EnvironmentInterface = IRefractiveIndexModel<IMediumPropertyModel<IMagneticFieldModel<IMediumModel>>>;
+ using EnvironmentInterface =
+ IRefractiveIndexModel<IMediumPropertyModel<IMagneticFieldModel<IMediumModel>>>;
using EnvType = Environment<EnvironmentInterface>;
EnvType env;
auto& universe = *(env.getUniverse());
@@ -84,8 +85,8 @@ int main() {
auto world = EnvType::createNode<Sphere>(Point{rootCS, 0_m, 0_m, 0_m}, 150_km);
- using MyHomogeneousModel = UniformRefractiveIndex<MediumPropertyModel<
- UniformMagneticField<HomogeneousMedium<EnvironmentInterface>>>>;
+ using MyHomogeneousModel = UniformRefractiveIndex<
+ MediumPropertyModel<UniformMagneticField<HomogeneousMedium<EnvironmentInterface>>>>;
auto const Bmag{0.0003809_T};
MagneticFieldVector B{rootCS, 0_T, 0_T, Bmag};
diff --git a/examples/synchrotron_test_manual_tracking.cpp b/examples/synchrotron_test_manual_tracking.cpp
index e2bcb3e1d4bf6b7e702524ad4d9305f4a0193923..8921329f5d365d7f0e54003b7ed51e9417435315 100644
--- a/examples/synchrotron_test_manual_tracking.cpp
+++ b/examples/synchrotron_test_manual_tracking.cpp
@@ -95,8 +95,10 @@ int main() {
std::cout << "number of points in time: " << duration * sampleRate_ << std::endl;
// create 2 antennas
- TimeDomainAntenna ant1("CoREAS antenna", point1, rootCS, start, duration, sampleRate_, start);
- TimeDomainAntenna ant2("ZHS antenna", point1, rootCS, start, duration, sampleRate_, start);
+ TimeDomainAntenna ant1("CoREAS antenna", point1, rootCS, start, duration, sampleRate_,
+ start);
+ TimeDomainAntenna ant2("ZHS antenna", point1, rootCS, start, duration, sampleRate_,
+ start);
// construct a radio detector instance to store our antennas
AntennaCollection<TimeDomainAntenna> detectorCoREAS;
diff --git a/tests/modules/testRadio.cpp b/tests/modules/testRadio.cpp
index f57fec8d91537fedda9698a82079a46fe2077992..4d019d9ecad783ee3e501d5b22a8130ca73026d2 100644
--- a/tests/modules/testRadio.cpp
+++ b/tests/modules/testRadio.cpp
@@ -259,245 +259,248 @@ TEST_CASE("Radio", "[processes]") {
} // END: SECTION("ZHS process")
- SECTION("Radio extreme cases") {
-
- // Environment
- using EnvironmentInterface =
- IRefractiveIndexModel<IMediumPropertyModel<IMagneticFieldModel<IMediumModel>>>;
-
- using EnvType = Environment<EnvironmentInterface>;
- EnvType envRadio;
- CoordinateSystemPtr const& rootCSRadio = envRadio.getCoordinateSystem();
- Point const center{rootCSRadio, 0_m, 0_m, 0_m};
-
- create_5layer_atmosphere<EnvironmentInterface, MyExtraEnv>(
- envRadio, AtmosphereId::LinsleyUSStd, center, 1.415, Medium::AirDry1Atm,
- MagneticFieldVector{rootCSRadio, 0_T, 50_uT, 0_T});
-
- // now create antennas and detectors
- // the antennas location
- const auto point1{Point(envRadio.getCoordinateSystem(), 0_m, 0_m, 0_m)};
-
- // track points
- 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});
-
- // create times for the antenna
- const TimeType start{0_s};
- const TimeType duration{100_ns};
- const InverseTimeType sample{1e+12_Hz};
- const TimeType duration_dummy{2_s};
- const InverseTimeType sample_dummy{1_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);
- // construct a radio detector instance to store our antennas
- AntennaCollection<TimeDomainAntenna> detector;
- AntennaCollection<TimeDomainAntenna> detector_dummy;
- // add the antennas to the detector
- detector.addAntenna(ant1);
- detector_dummy.addAntenna(ant2);
-
- // create a new stack for each trial
- setup::Stack<EnvType> stack;
- // create a particle
- const Code particle{Code::Electron};
- const Code particle2{Code::Proton};
-
- const auto pmass{get_mass(particle)};
- const auto pmass2{get_mass(particle2)};
- // construct an energy
- const HEPEnergyType E0{1_TeV};
- // compute the necessary momentumn
- const HEPMomentumType P0{sqrt(E0 * E0 - pmass * pmass)};
- // and create the momentum vector
- const auto plab{MomentumVector(rootCSRadio, {P0, 0_GeV, 0_GeV})};
- // add the particle to the stack
- auto const particle_stack{stack.addParticle(std::make_tuple(
- particle, calculate_kinetic_energy(plab.getNorm(), get_mass(particle)),
- plab.normalized(), point_1, 0_ns))};
-
- // particle stack with proton
- auto const particle_stack_proton{stack.addParticle(std::make_tuple(
- particle2, calculate_kinetic_energy(plab.getNorm(), get_mass(particle2)),
- plab.normalized(), point_1, 0_ns))};
-
- // feed radio with a proton track to check that it skips that track.
- TimeType tp{(point_2 - point_1).getNorm() / (0.999 * constants::c)};
- VelocityVector vp{(point_2 - point_1) / tp};
- Line lp{point_1, vp};
- StraightTrajectory track_p{lp, tp};
- Step step_proton(particle_stack_proton, track_p);
-
- // feed radio with a track that triggers zhs like approx in coreas and fraunhofer limit check for zhs
- TimeType th{(point_4 - point1).getNorm() / constants::c};
- VelocityVector vh{(point_4 - point1) / th};
- Line lh{point1, vh};
- StraightTrajectory track_h{lh, th};
- Step step_h(particle_stack, track_h);
-
- // create radio process instances
- RadioProcess<decltype(detector),
- CoREAS<decltype(detector), decltype(SimplePropagator(envRadio))>,
- decltype(SimplePropagator(envRadio))>
- coreas(detector, envRadio);
-
- RadioProcess<decltype(detector),
- 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);
-
- // 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))>
- coreas_dummy(detector_dummy, envRadio);
-
- RadioProcess<decltype(detector_dummy),
- ZHS<decltype(detector_dummy), decltype(SimplePropagator(envRadio))>,
- decltype(SimplePropagator(envRadio))>
- zhs_dummy(detector_dummy, envRadio);
-
- coreas_dummy.doContinuous(step_h, true);
- zhs_dummy.doContinuous(step_h, true);
-
- } // END: SECTION("Radio extreme cases")
+ SECTION("Radio extreme cases") {
+
+ // Environment
+ using EnvironmentInterface =
+ IRefractiveIndexModel<IMediumPropertyModel<IMagneticFieldModel<IMediumModel>>>;
+
+ using EnvType = Environment<EnvironmentInterface>;
+ EnvType envRadio;
+ CoordinateSystemPtr const& rootCSRadio = envRadio.getCoordinateSystem();
+ Point const center{rootCSRadio, 0_m, 0_m, 0_m};
+
+ create_5layer_atmosphere<EnvironmentInterface, MyExtraEnv>(
+ envRadio, AtmosphereId::LinsleyUSStd, center, 1.415, Medium::AirDry1Atm,
+ MagneticFieldVector{rootCSRadio, 0_T, 50_uT, 0_T});
+
+ // now create antennas and detectors
+ // the antennas location
+ const auto point1{Point(envRadio.getCoordinateSystem(), 0_m, 0_m, 0_m)};
+
+ // track points
+ 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});
+
+ // create times for the antenna
+ const TimeType start{0_s};
+ const TimeType duration{100_ns};
+ const InverseTimeType sample{1e+12_Hz};
+ const TimeType duration_dummy{2_s};
+ const InverseTimeType sample_dummy{1_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);
+ // construct a radio detector instance to store our antennas
+ AntennaCollection<TimeDomainAntenna> detector;
+ AntennaCollection<TimeDomainAntenna> detector_dummy;
+ // add the antennas to the detector
+ detector.addAntenna(ant1);
+ detector_dummy.addAntenna(ant2);
+
+ // create a new stack for each trial
+ setup::Stack<EnvType> stack;
+ // create a particle
+ const Code particle{Code::Electron};
+ const Code particle2{Code::Proton};
+
+ const auto pmass{get_mass(particle)};
+ const auto pmass2{get_mass(particle2)};
+ // construct an energy
+ const HEPEnergyType E0{1_TeV};
+ // compute the necessary momentumn
+ const HEPMomentumType P0{sqrt(E0 * E0 - pmass * pmass)};
+ // and create the momentum vector
+ const auto plab{MomentumVector(rootCSRadio, {P0, 0_GeV, 0_GeV})};
+ // add the particle to the stack
+ auto const particle_stack{stack.addParticle(std::make_tuple(
+ particle, calculate_kinetic_energy(plab.getNorm(), get_mass(particle)),
+ plab.normalized(), point_1, 0_ns))};
+
+ // particle stack with proton
+ auto const particle_stack_proton{stack.addParticle(std::make_tuple(
+ particle2, calculate_kinetic_energy(plab.getNorm(), get_mass(particle2)),
+ plab.normalized(), point_1, 0_ns))};
+
+ // feed radio with a proton track to check that it skips that track.
+ TimeType tp{(point_2 - point_1).getNorm() / (0.999 * constants::c)};
+ VelocityVector vp{(point_2 - point_1) / tp};
+ Line lp{point_1, vp};
+ StraightTrajectory track_p{lp, tp};
+ Step step_proton(particle_stack_proton, track_p);
+
+ // feed radio with a track that triggers zhs like approx in coreas and fraunhofer
+ // limit check for zhs
+ TimeType th{(point_4 - point1).getNorm() / constants::c};
+ VelocityVector vh{(point_4 - point1) / th};
+ Line lh{point1, vh};
+ StraightTrajectory track_h{lh, th};
+ Step step_h(particle_stack, track_h);
+
+ // create radio process instances
+ RadioProcess<decltype(detector),
+ CoREAS<decltype(detector), decltype(SimplePropagator(envRadio))>,
+ decltype(SimplePropagator(envRadio))>
+ coreas(detector, envRadio);
+
+ RadioProcess<decltype(detector),
+ 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);
+
+ // 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))>
+ coreas_dummy(detector_dummy, envRadio);
+
+ RadioProcess<decltype(detector_dummy),
+ ZHS<decltype(detector_dummy), decltype(SimplePropagator(envRadio))>,
+ decltype(SimplePropagator(envRadio))>
+ zhs_dummy(detector_dummy, envRadio);
+
+ coreas_dummy.doContinuous(step_h, true);
+ zhs_dummy.doContinuous(step_h, true);
+
+ } // END: SECTION("Radio extreme cases")
} // END: TEST_CASE("Radio", "[processes]")
TEST_CASE("Antennas") {
- SECTION("TimeDomainAntenna") {
+ SECTION("TimeDomainAntenna") {
- // create an environment so we can get a coordinate system
- using EnvType = Environment<IRefractiveIndexModel<IMediumModel>>;
- EnvType env6;
+ // create an environment so we can get a coordinate system
+ using EnvType = Environment<IRefractiveIndexModel<IMediumModel>>;
+ EnvType env6;
- using UniRIndex =
+ using UniRIndex =
UniformRefractiveIndex<HomogeneousMedium<IRefractiveIndexModel<IMediumModel>>>;
- // the antenna location
- const auto point1{Point(env6.getCoordinateSystem(), 1_m, 2_m, 3_m)};
- const auto point2{Point(env6.getCoordinateSystem(), 4_m, 5_m, 6_m)};
-
- // get a coordinate system
- const CoordinateSystemPtr rootCS6 = env6.getCoordinateSystem();
-
- auto Medium6 = EnvType::createNode<Sphere>(
- Point{rootCS6, 0_m, 0_m, 0_m}, 1_km * std::numeric_limits<double>::infinity());
-
- auto const props6 = Medium6->setModelProperties<UniRIndex>(
- 1, 1_kg / (1_m * 1_m * 1_m), NuclearComposition({Code::Nitrogen}, {1.}));
-
- env6.getUniverse()->addChild(std::move(Medium6));
-
- // create times for the antenna
- const TimeType t1{10_s};
- const TimeType t2{10_s};
- const InverseTimeType t3{1 / 1_s};
- const TimeType t4{11_s};
-
- // check that I can create an antenna at (1, 2, 3)
- TimeDomainAntenna ant1("antenna_name", point1, rootCS6, t1, t2, t3, t1);
- TimeDomainAntenna ant2("antenna_name2", point2, rootCS6, t4, t2, t3, t4);
-
- // assert that the antenna name is correct
- REQUIRE(ant1.getName() == "antenna_name");
- REQUIRE(ant2.getName() == "antenna_name2");
-
- // and check that the antenna is at the right location
- REQUIRE((ant1.getLocation() - point1).getNorm() < 1e-12 * 1_m);
- REQUIRE((ant2.getLocation() - point2).getNorm() < 1e-12 * 1_m);
-
- // construct a radio detector instance to store our antennas
- AntennaCollection<TimeDomainAntenna> detector;
-
- // add this antenna to the process
- detector.addAntenna(ant1);
- detector.addAntenna(ant2);
- CHECK(detector.size() == 2);
-
- // get a unit vector
- Vector<dimensionless_d> v1(rootCS6, {0, 0, 1});
- Vector<ElectricFieldType::dimension_type> v11(rootCS6,
- {10_V / 1_m, 10_V / 1_m, 10_V / 1_m});
-
- Vector<dimensionless_d> v2(rootCS6, {0, 1, 0});
- Vector<ElectricFieldType::dimension_type> v22(rootCS6,
- {20_V / 1_m, 20_V / 1_m, 20_V / 1_m});
-
- // use receive methods
- ant1.receive(15_s, v1, v11);
- ant2.receive(16_s, v2, v22);
-
- // use getDataX,Y,Z() and getAxis() methods
- auto Ex = ant1.getWaveformX();
- CHECK(Ex[5] - 10 == 0);
- auto tx = ant1.getAxis();
- CHECK(tx[5] - 5 * 1_s / 1_ns == Approx(0.0));
- auto Ey = ant1.getWaveformY();
- CHECK(Ey[5] - 10 == 0);
- auto Ez = ant1.getWaveformZ();
- CHECK(Ez[5] - 10 == 0);
- auto ty = ant1.getAxis();
- auto tz = ant1.getAxis();
- CHECK(tx[5] - ty[5] == 0);
- CHECK(ty[5] - tz[5] == 0);
- auto Ex2 = ant2.getWaveformX();
- CHECK(Ex2[5] - 20 == 0);
- auto Ey2 = ant2.getWaveformY();
- CHECK(Ey2[5] - 20 == 0);
- auto Ez2 = ant2.getWaveformZ();
- CHECK(Ez2[5] - 20 == 0);
-
- // the following creates a star-shaped pattern of antennas in the ground
- AntennaCollection<TimeDomainAntenna> detector__;
- const auto point11{Point(env6.getCoordinateSystem(), 1000_m, 20_m, 30_m)};
- const TimeType t2222{1e-6_s};
- const InverseTimeType t3333{1e+9_Hz};
-
- std::vector<std::string> antenna_names;
- std::vector<Point> antenna_locations;
- for (auto radius_ = 100_m; radius_ <= 200_m; radius_ += 100_m) {
- for (auto phi_ = 0; phi_ <= 315; phi_ += 45) {
- auto phiRad_ = phi_ / 180. * M_PI;
- auto const point_{Point(env6.getCoordinateSystem(), radius_ * cos(phiRad_),
- radius_ * sin(phiRad_), 0_m)};
- antenna_locations.push_back(point_);
- auto time__{(point11 - point_).getNorm() / constants::c};
- const int rr_ = static_cast<int>(radius_ / 1_m);
- std::string var_ = "antenna_R=" + std::to_string(rr_) +
- "_m-Phi=" + std::to_string(phi_) + "degrees";
- antenna_names.push_back(var_);
- TimeDomainAntenna ant111(var_, point_, rootCS6, time__, t2222, t3333, time__);
- detector__.addAntenna(ant111);
- }
- }
-
- CHECK(detector__.size() == 16);
- CHECK(detector__.getAntennas().size() == 16);
- int i = 0;
- // this prints out the antenna names and locations
- for (auto const antenna: detector__.getAntennas()) {
- CHECK(antenna.getName() == antenna_names[i]);
- CHECK(distance(antenna.getLocation(), antenna_locations[i]) / 1_m == 0);
- i++;
- }
-
- // Check the .at() method for radio detectors
- for (size_t i = 0; i <= (detector__.size()-1); i++) {
- CHECK(detector__.at(i).getName() == antenna_names[i]);
- CHECK(distance(detector__.at(i).getLocation(), antenna_locations[i]) / 1_m == 0);
- }
-
- } // END: SECTION("TimeDomainAntenna")
+ // the antenna location
+ const auto point1{Point(env6.getCoordinateSystem(), 1_m, 2_m, 3_m)};
+ const auto point2{Point(env6.getCoordinateSystem(), 4_m, 5_m, 6_m)};
+
+ // get a coordinate system
+ const CoordinateSystemPtr rootCS6 = env6.getCoordinateSystem();
+
+ auto Medium6 = EnvType::createNode<Sphere>(
+ Point{rootCS6, 0_m, 0_m, 0_m}, 1_km * std::numeric_limits<double>::infinity());
+
+ auto const props6 = Medium6->setModelProperties<UniRIndex>(
+ 1, 1_kg / (1_m * 1_m * 1_m), NuclearComposition({Code::Nitrogen}, {1.}));
+
+ env6.getUniverse()->addChild(std::move(Medium6));
+
+ // create times for the antenna
+ const TimeType t1{10_s};
+ const TimeType t2{10_s};
+ const InverseTimeType t3{1 / 1_s};
+ const TimeType t4{11_s};
+
+ // check that I can create an antenna at (1, 2, 3)
+ TimeDomainAntenna ant1("antenna_name", point1, rootCS6, t1, t2, t3, t1);
+ TimeDomainAntenna ant2("antenna_name2", point2, rootCS6, t4, t2, t3, t4);
+
+ // assert that the antenna name is correct
+ REQUIRE(ant1.getName() == "antenna_name");
+ REQUIRE(ant2.getName() == "antenna_name2");
+
+ // and check that the antenna is at the right location
+ REQUIRE((ant1.getLocation() - point1).getNorm() < 1e-12 * 1_m);
+ REQUIRE((ant2.getLocation() - point2).getNorm() < 1e-12 * 1_m);
+
+ // construct a radio detector instance to store our antennas
+ AntennaCollection<TimeDomainAntenna> detector;
+
+ // add this antenna to the process
+ detector.addAntenna(ant1);
+ detector.addAntenna(ant2);
+ CHECK(detector.size() == 2);
+
+ // get a unit vector
+ Vector<dimensionless_d> v1(rootCS6, {0, 0, 1});
+ Vector<ElectricFieldType::dimension_type> v11(rootCS6,
+ {10_V / 1_m, 10_V / 1_m, 10_V / 1_m});
+
+ Vector<dimensionless_d> v2(rootCS6, {0, 1, 0});
+ Vector<ElectricFieldType::dimension_type> v22(rootCS6,
+ {20_V / 1_m, 20_V / 1_m, 20_V / 1_m});
+
+ // use receive methods
+ ant1.receive(15_s, v1, v11);
+ ant2.receive(16_s, v2, v22);
+
+ // use getDataX,Y,Z() and getAxis() methods
+ auto Ex = ant1.getWaveformX();
+ CHECK(Ex[5] - 10 == 0);
+ auto tx = ant1.getAxis();
+ CHECK(tx[5] - 5 * 1_s / 1_ns == Approx(0.0));
+ auto Ey = ant1.getWaveformY();
+ CHECK(Ey[5] - 10 == 0);
+ auto Ez = ant1.getWaveformZ();
+ CHECK(Ez[5] - 10 == 0);
+ auto ty = ant1.getAxis();
+ auto tz = ant1.getAxis();
+ CHECK(tx[5] - ty[5] == 0);
+ CHECK(ty[5] - tz[5] == 0);
+ auto Ex2 = ant2.getWaveformX();
+ CHECK(Ex2[5] - 20 == 0);
+ auto Ey2 = ant2.getWaveformY();
+ CHECK(Ey2[5] - 20 == 0);
+ auto Ez2 = ant2.getWaveformZ();
+ CHECK(Ez2[5] - 20 == 0);
+
+ // the following creates a star-shaped pattern of antennas in the ground
+ AntennaCollection<TimeDomainAntenna> detector__;
+ const auto point11{Point(env6.getCoordinateSystem(), 1000_m, 20_m, 30_m)};
+ const TimeType t2222{1e-6_s};
+ const InverseTimeType t3333{1e+9_Hz};
+
+ std::vector<std::string> antenna_names;
+ std::vector<Point> antenna_locations;
+ for (auto radius_ = 100_m; radius_ <= 200_m; radius_ += 100_m) {
+ for (auto phi_ = 0; phi_ <= 315; phi_ += 45) {
+ auto phiRad_ = phi_ / 180. * M_PI;
+ auto const point_{Point(env6.getCoordinateSystem(), radius_ * cos(phiRad_),
+ radius_ * sin(phiRad_), 0_m)};
+ antenna_locations.push_back(point_);
+ auto time__{(point11 - point_).getNorm() / constants::c};
+ const int rr_ = static_cast<int>(radius_ / 1_m);
+ std::string var_ = "antenna_R=" + std::to_string(rr_) +
+ "_m-Phi=" + std::to_string(phi_) + "degrees";
+ antenna_names.push_back(var_);
+ TimeDomainAntenna ant111(var_, point_, rootCS6, time__, t2222, t3333, time__);
+ detector__.addAntenna(ant111);
+ }
+ }
+
+ CHECK(detector__.size() == 16);
+ CHECK(detector__.getAntennas().size() == 16);
+ int i = 0;
+ // this prints out the antenna names and locations
+ for (auto const antenna : detector__.getAntennas()) {
+ CHECK(antenna.getName() == antenna_names[i]);
+ CHECK(distance(antenna.getLocation(), antenna_locations[i]) / 1_m == 0);
+ i++;
+ }
+
+ // Check the .at() method for radio detectors
+ for (size_t i = 0; i <= (detector__.size() - 1); i++) {
+ CHECK(detector__.at(i).getName() == antenna_names[i]);
+ CHECK(distance(detector__.at(i).getLocation(), antenna_locations[i]) / 1_m == 0);
+ }
+
+ } // END: SECTION("TimeDomainAntenna")
} // END: TEST_CASE("Antennas")