diff --git a/corsika/modules/radio/CoREAS.hpp b/corsika/modules/radio/CoREAS.hpp
index ef22fcd07768d4dc31faf17f1be4a2e2aeb7ea1d..b56da0f1db5cc3615515be31cf3cb6307155fb7c 100755
--- a/corsika/modules/radio/CoREAS.hpp
+++ b/corsika/modules/radio/CoREAS.hpp
@@ -53,27 +53,35 @@ namespace corsika {
// get the global simulation time for that track. (best guess for now)
auto startTime_{particle.getTime()}; // time at the start point of the track hopefully. I should use something similar to fCoreHitTime (?)
- std::cout << "startTime_: " << startTime_ << std::endl;
+ std::cout << "Particle time at start of track: " << startTime_ << std::endl;
auto endTime_{particle.getTime() + track.getDuration()};
- std::cout << "endTime_: " << endTime_ << std::endl;
- std::cout << "track.getDuration(): " << track.getDuration() << std::endl;
+ std::cout << "Particle time at end of track: " << endTime_ << std::endl;
+ std::cout << "Duration of track: " << track.getDuration() << std::endl;
+ // TODO: this should be fixed with the continuous processes new design, so we can get the energy at start and end of track
// gamma factor is calculated using beta
// auto startGamma_ {1. / sqrt(1. - (startBeta_ * startBeta_))};
// auto endGamma_ {1. / sqrt(1. - (endBeta_ * endBeta_))};
// get start and end position of the track
auto startPoint_{track.getPosition(0)};
- std::cout << "STARTPOINT : " << startPoint_ << std::endl;
+ std::cout << "STARTPOINT position: " << startPoint_ << std::endl;
auto endPoint_{track.getPosition(1)};
-// track.getVelocity(1); TODO: use this for velocity weight factors
- std::cout << "ENDPOINT : " << endPoint_ << std::endl;
+ std::cout << "ENDPOINT position: " << endPoint_ << std::endl;
// beta is velocity / speed of light. Start & end should be the same!
// auto beta_ {static_cast<double>((endPoint_.distance_to(startPoint_) / 1_m ) / (endTime_/ 1_s - startTime_/ 1_s) )};
// auto beta_ {((endPoint_.distance_to(startPoint_)) / (endTime_ - startTime_)).normalized()};
- auto beta_{((endPoint_ - startPoint_) / (constants::c * (endTime_ - startTime_))).normalized()};
- std::cout << "BETA_: " << beta_ << std::endl;
+ auto beta_ {((endPoint_ - startPoint_) / (constants::c * (endTime_ - startTime_))).normalized()};
+ std::cout << "BETA v/c: " << beta_ << std::endl;
+
+ // ideally these variables should use the improved getEnergy() method once it's there
+// auto preBeta_ {(track.getVelocity(0) / constants::c)};
+// std::cout << "PreBeta: " << preBeta_ << std::endl;
+// auto postBeta_ {(track.getVelocity(1) / constants::c)};
+// std::cout << "PostBeta: " << postBeta_ << std::endl;
+// auto velocityWeight_ {0.5 * (preBeta_ + postBeta_) / beta_};
+// std::cout << "Velocity Weight: " << velocityWeight_ << std::endl;
// get particle charge
auto const charge_{get_charge(particle.getPID())};
@@ -84,33 +92,69 @@ namespace corsika {
// loop over each antenna in the antenna collection (detector)
for (auto& antenna : detector_.getAntennas()) {
+ // check with which antenna we work in this loop
std::cout << "ANTENNA: " << antenna.getName() << std::endl;
// get the SignalPathCollection (path1) from the start "endpoint" to the antenna.
- auto paths1{this->propagator_.propagate(startPoint_, antenna.getLocation(), 1_m)}; // TODO: Need to add the stepsize to .propagate()!!!!
+ auto paths1{this->propagator_.propagate(startPoint_, antenna.getLocation(), 1_m)}; // TODO: Add the stepsize to .propagate() at some point
// get the SignalPathCollection (path2) from the end "endpoint" to the antenna.
auto paths2{this->propagator_.propagate(endPoint_, antenna.getLocation(), 1_m)};
+ // throw an exception if path sizes don't match
+ try {
// loop over both paths at once and directly compare 'start' and 'end' attributes
- for (size_t i = (paths1.size() == paths2.size()) ? 0 : paths1.size(); // TODO: throw an exception if the sizes don't match
+ for (size_t i = (paths1.size() == paths2.size()) ? 0 : throw (i = paths1.size());
(i < paths1.size() && i < paths2.size()); i++) {
// calculate preDoppler factor
double preDoppler_{1. - paths1[i].refractive_index_source_ *
- beta_.dot(paths1[i].emit_)}; // are you sure this is path.emit and not path.receive?
+ beta_.dot(paths1[i].emit_)}; // are you sure this is path.emit and not path.receive?
std::cout << "***** preDoppler: " << preDoppler_ << std::endl;
- std::cout << "***** preEmmit_: " << paths1[i].emit_ << std::endl;
- std::cout << "***** preRefractive_index: " << paths1[i].refractive_index_source_ << std::endl;
+ std::cout << "***** preEmmit_vector: " << paths1[i].emit_ << std::endl;
+ std::cout << "***** preRefractive_Index_Source: " << paths1[i].refractive_index_source_
+ << std::endl;
+
+ // check if preDoppler has become zero in case of refractive index of unity because of numerical limitations
+ if (preDoppler_ == 0) {
+ // redo calculation with higher precision
+ long double indexL_ {paths1[i].refractive_index_source_};
+ long double betaX_ {static_cast<double>(beta_.getComponents().getX())};
+ long double betaY_ {static_cast<double>(beta_.getComponents().getY())};
+ long double betaZ_ {static_cast<double>(beta_.getComponents().getZ())};
+ long double startX_ {static_cast<double>(paths1[i].emit_.getComponents().getX())};
+ long double startY_ {static_cast<double>(paths1[i].emit_.getComponents().getY())};
+ long double startZ_ {static_cast<double>(paths1[i].emit_.getComponents().getZ())};
+ long double doppler = 1.0l - indexL_ * (betaX_ * startX_ +
+ betaY_ * startY_ + betaZ_ * startZ_);
+ preDoppler_ = doppler;
+ }
// calculate postDoppler factor
- double postDoppler_{1. - paths2[i].refractive_index_source_ *
- beta_.dot(paths2[i].emit_)}; // maybe this is path.receive_ (?)
+ double postDoppler_{
+ 1. - paths2[i].refractive_index_source_ *
+ beta_.dot(paths2[i].emit_)}; // maybe this is path.receive_ (?)
std::cout << "***** postDoppler: " << postDoppler_ << std::endl;
- std::cout << "***** postEmmit_: " << paths2[i].emit_ << std::endl;
- std::cout << "***** postRefractive_index: " << paths2[i].refractive_index_source_ << std::endl;
-
- // calculate receive time for startpoint
+ std::cout << "***** postEmmit_vector: " << paths2[i].emit_ << std::endl;
+ std::cout << "***** postRefractive_Index_Source: "
+ << paths2[i].refractive_index_source_ << std::endl;
+
+ // check if postDoppler has become zero in case of refractive index of unity because of numerical limitations
+ if (postDoppler_ == 0) {
+ // redo calculation with higher precision
+ long double indexL_ {paths2[i].refractive_index_source_};
+ long double betaX_ {static_cast<double>(beta_.getComponents().getX())};
+ long double betaY_ {static_cast<double>(beta_.getComponents().getY())};
+ long double betaZ_ {static_cast<double>(beta_.getComponents().getZ())};
+ long double endX_ {static_cast<double>(paths2[i].emit_.getComponents().getX())};
+ long double endY_ {static_cast<double>(paths2[i].emit_.getComponents().getY())};
+ long double endZ_ {static_cast<double>(paths2[i].emit_.getComponents().getZ())};
+ long double doppler = 1.0l - indexL_ * (betaX_ * endX_ +
+ betaY_ * endY_ + betaZ_ * endZ_);
+ postDoppler_ = doppler;
+ }
+
+ // calculate receive time for startpoint (aka time delay)
auto startPointReceiveTime_{paths1[i].propagation_time_ + startTime_};
std::cout << "STARTPOINT RECEIVE TIME BEFORE IFS: " << startPointReceiveTime_ << std::endl; // TODO: time 0 is when the imaginary primary hits the ground
@@ -119,32 +163,32 @@ namespace corsika {
std::cout << "ENDPOINT RECEIVE TIME BEFORE IFS: " << endPointReceiveTime_ << std::endl;
// get receive unit vector for startpoint
- auto ReceiveVectorStart_ {paths1[i].receive_};
+ auto ReceiveVectorStart_{paths1[i].receive_};
std::cout << "start receive unit vector before ifs: " << ReceiveVectorStart_ << std::endl;
// get receive unit vector for endpoint
- auto ReceiveVectorEnd_ {paths2[i].receive_};
+ auto ReceiveVectorEnd_{paths2[i].receive_};
std::cout << "end receive unit vector before ifs: " << ReceiveVectorEnd_ << std::endl;
-
////////////////////////////////////////////////////////////////////////////////
- // start comparing stuff
+
// perform ZHS-like calculation close to Cherenkov angle and for refractive index at antenna location greater than 1
- if ( (paths1[i].refractive_index_destination_ > 1) &&
- (std::fabs(preDoppler_) <= approxThreshold_ || std::fabs(postDoppler_) <= approxThreshold_) ) {
+ if ((paths1[i].refractive_index_destination_ > 1) &&
+ (std::fabs(preDoppler_) <= approxThreshold_ || std::fabs(postDoppler_) <= approxThreshold_)) {
// clear the existing paths for this particle and track
paths1.clear();
paths2.clear();
// get global simulation time for the middle point of that track.
- TimeType midTime_ {((startPoint_ - endPoint_).getNorm()/2) / track.getVelocity(0).getNorm()};
-// auto midTime_{particle.getTime() - (track.getDuration() / 2)}; // this is not the geometrical calculation
+ TimeType midTime_{((startPoint_ - endPoint_).getNorm() / 2) /
+ track.getVelocity(0).getNorm()};
+ // auto midTime_{particle.getTime() - (track.getDuration() / 2)}; // this is not the geometrical calculation
// get "mid" position of the track geometrically
- auto midVector_ { (startPoint_ - endPoint_) / 2 };
- auto midPoint_ {Point(midVector_.getCoordinateSystem(), midVector_.getComponents().getX(),
- midVector_.getComponents().getY(), midVector_.getComponents().getZ())};
+ auto midVector_{(startPoint_ - endPoint_) / 2};
+ auto midPoint_{Point(midVector_.getCoordinateSystem(), midVector_.getComponents().getX(),
+ midVector_.getComponents().getY(), midVector_.getComponents().getZ())};
// get the SignalPathCollection (path3) from the middle "endpoint" to the antenna.
auto paths3{this->propagator_.propagate(midPoint_, antenna.getLocation(), 1_m)};
@@ -153,7 +197,22 @@ namespace corsika {
for (auto const& path : paths3) {
auto const midPointReceiveTime_{path.propagation_time_ + midTime_};
- auto midDoppler_{1. - path.refractive_index_source_ * beta_.dot(path.emit_)};
+ double midDoppler_{1. - path.refractive_index_source_ * beta_.dot(path.emit_)};
+
+ // check if midDoppler has become zero because of numerical limitations
+ if (midDoppler_ == 0) {
+ // redo calculation with higher precision
+ long double indexL_ {path.refractive_index_source_};
+ long double betaX_ {static_cast<double>(beta_.getComponents().getX())};
+ long double betaY_ {static_cast<double>(beta_.getComponents().getY())};
+ long double betaZ_ {static_cast<double>(beta_.getComponents().getZ())};
+ long double midX_ {static_cast<double>(path.emit_.getComponents().getX())};
+ long double midY_ {static_cast<double>(path.emit_.getComponents().getY())};
+ long double midZ_ {static_cast<double>(path.emit_.getComponents().getZ())};
+ long double doppler = 1.0l - indexL_ * (betaX_ * midX_ +
+ betaY_ * midY_ + betaZ_ * midZ_);
+ midDoppler_ = doppler;
+ }
// change the values of the receive unit vectors of start and end
ReceiveVectorStart_ = path.receive_;
@@ -162,19 +221,22 @@ namespace corsika {
// CoREAS calculation -> get ElectricFieldVector for "midPoint"
ElectricFieldVector EVmid_ =
path.receive_.cross(path.receive_.cross(beta_)).getComponents() /
- (path.R_distance_ * midDoppler_) * (antenna.sample_rate_ / (4 * M_PI)) * ((1 / constants::epsilonZero) * (1 / constants::c)) * charge_;
+ (path.R_distance_ * midDoppler_) * (antenna.sample_rate_ / (4 * M_PI)) *
+ ((1 / constants::epsilonZero) * (1 / constants::c)) * charge_;
- ElectricFieldVector EV1_ {EVmid_};
- ElectricFieldVector EV2_ {-EVmid_};
+ ElectricFieldVector EV1_{EVmid_};
+ ElectricFieldVector EV2_{-EVmid_};
- auto deltaT_{(endPoint_ - startPoint_).getNorm() / (constants::c * beta_.getNorm() * std::fabs(midDoppler_))}; // TODO: Caution with this!
+ auto deltaT_{(endPoint_ - startPoint_).getNorm() /
+ (constants::c * beta_.getNorm() *
+ std::fabs(midDoppler_))}; // TODO: Caution with this!
- if (startPointReceiveTime_ < endPointReceiveTime_) // EVstart_ arrives earlier
+ if (startPointReceiveTime_ <
+ endPointReceiveTime_) // EVstart_ arrives earlier
{
startPointReceiveTime_ = midPointReceiveTime_ - 0.5 * deltaT_;
endPointReceiveTime_ = midPointReceiveTime_ + 0.5 * deltaT_;
- }
- else // EVend_ arrives earlier
+ } else // EVend_ arrives earlier
{
startPointReceiveTime_ = midPointReceiveTime_ + 0.5 * deltaT_;
endPointReceiveTime_ = midPointReceiveTime_ - 0.5 * deltaT_;
@@ -189,66 +251,81 @@ namespace corsika {
EV1_ = EV1_ * std::fabs((deltaT_ / 1_s) / gridResolution_);
EV2_ = EV2_ * std::fabs((deltaT_ / 1_s) / gridResolution_);
- const long startBin = static_cast<long>(std::floor((startPointReceiveTime_ / 1_s)/gridResolution_+0.5l));
- const long endBin = static_cast<long>(std::floor((endPointReceiveTime_ / 1_s) /gridResolution_+0.5l));
- const double startBinFraction = ((startPointReceiveTime_ / 1_s)/gridResolution_)-std::floor((startPointReceiveTime_ / 1_s)/gridResolution_);
- const double endBinFraction = ((endPointReceiveTime_ / 1_s)/gridResolution_)-std::floor((endPointReceiveTime_ / 1_s)/gridResolution_);
+ const long startBin = static_cast<long>(
+ std::floor((startPointReceiveTime_ / 1_s) / gridResolution_ + 0.5l));
+ const long endBin = static_cast<long>(
+ std::floor((endPointReceiveTime_ / 1_s) / gridResolution_ + 0.5l));
+ const double startBinFraction =
+ ((startPointReceiveTime_ / 1_s) / gridResolution_) -
+ std::floor((startPointReceiveTime_ / 1_s) / gridResolution_);
+ const double endBinFraction =
+ ((endPointReceiveTime_ / 1_s) / gridResolution_) -
+ std::floor((endPointReceiveTime_ / 1_s) / gridResolution_);
// only do timing modification if contributions would land in same bin
if (startBin == endBin) {
// if startE arrives before endE
if ((deltaT_ / 1_s) >= 0) {
- if ((startBinFraction >= 0.5) && (endBinFraction >= 0.5)) // both points left of bin center
+ if ((startBinFraction >= 0.5) &&
+ (endBinFraction >= 0.5)) // both points left of bin center
{
- startPointReceiveTime_ = startPointReceiveTime_ - gridResolution_ * 1_s; // shift EV1_ to previous gridpoint
- }
- else if ((startBinFraction < 0.5) && (endBinFraction < 0.5)) // both points right of bin center
+ startPointReceiveTime_ =
+ startPointReceiveTime_ -
+ gridResolution_ * 1_s; // shift EV1_ to previous gridpoint
+ } else if ((startBinFraction < 0.5) &&
+ (endBinFraction < 0.5)) // both points right of bin center
{
- endPointReceiveTime_ = endPointReceiveTime_ + gridResolution_ * 1_s; // shift EV2_ to next gridpoint
- }
- else // points on both sides of bin center
+ endPointReceiveTime_ =
+ endPointReceiveTime_ +
+ gridResolution_ * 1_s; // shift EV2_ to next gridpoint
+ } else // points on both sides of bin center
{
- const double leftDist = 1.0-startBinFraction;
+ const double leftDist = 1.0 - startBinFraction;
const double rightDist = endBinFraction;
// check if asymmetry to right or left
- if (rightDist >= leftDist)
- {
- endPointReceiveTime_ = endPointReceiveTime_ + gridResolution_ * 1_s; // shift EV2_ to next gridpoint
- }
- else
- {
- startPointReceiveTime_ = startPointReceiveTime_ - gridResolution_ * 1_s; // shift EV1_ to previous gridpoint
+ if (rightDist >= leftDist) {
+ endPointReceiveTime_ =
+ endPointReceiveTime_ +
+ gridResolution_ * 1_s; // shift EV2_ to next gridpoint
+ } else {
+ startPointReceiveTime_ =
+ startPointReceiveTime_ -
+ gridResolution_ * 1_s; // shift EV1_ to previous gridpoint
}
}
- }
- else // if endE arrives before startE
+ } else // if endE arrives before startE
{
- if ((startBinFraction >= 0.5) && (endBinFraction >= 0.5)) // both points left of bin center
+ if ((startBinFraction >= 0.5) &&
+ (endBinFraction >= 0.5)) // both points left of bin center
{
- endPointReceiveTime_ = endPointReceiveTime_ - gridResolution_ * 1_s; // shift EV2_ to previous gridpoint
- }
- else if ((startBinFraction < 0.5) && (endBinFraction < 0.5)) // both points right of bin center
+ endPointReceiveTime_ =
+ endPointReceiveTime_ -
+ gridResolution_ * 1_s; // shift EV2_ to previous gridpoint
+ } else if ((startBinFraction < 0.5) &&
+ (endBinFraction < 0.5)) // both points right of bin center
{
- startPointReceiveTime_ = startPointReceiveTime_ + gridResolution_ * 1_s; // shift EV1_ to next gridpoint
- }
- else // points on both sides of bin center
+ startPointReceiveTime_ =
+ startPointReceiveTime_ +
+ gridResolution_ * 1_s; // shift EV1_ to next gridpoint
+ } else // points on both sides of bin center
{
- const double leftDist = 1.0-endBinFraction;
+ const double leftDist = 1.0 - endBinFraction;
const double rightDist = startBinFraction;
// check if asymmetry to right or left
- if (rightDist >= leftDist)
- {
- startPointReceiveTime_ = startPointReceiveTime_ + gridResolution_ * 1_s; // shift EV1_ to next gridpoint
- }
- else
- {
- endPointReceiveTime_ = endPointReceiveTime_ - gridResolution_ * 1_s; // shift EV2_ to previous gridpoint
+ if (rightDist >= leftDist) {
+ startPointReceiveTime_ =
+ startPointReceiveTime_ +
+ gridResolution_ * 1_s; // shift EV1_ to next gridpoint
+ } else {
+ endPointReceiveTime_ =
+ endPointReceiveTime_ -
+ gridResolution_ * 1_s; // shift EV2_ to previous gridpoint
}
}
} // End of else statement
- } // End of if for startbin == endbin
- } // End of if deltaT < gridresolution
+ } // End of if for startbin == endbin
+ } // End of if deltaT < gridresolution
// TODO: Be very careful with this. Maybe the EVs should be fed after the for loop of paths3
std::cout << "---------- RECEIVE INCIDENT USING ZHS-like APPROXIMATION ----------" << std::endl;
@@ -284,80 +361,94 @@ namespace corsika {
std::cout << "CHECK EV2 VALUE : " << EV2_ << std::endl;
-
if ((preDoppler_ < 1.e-9) || (postDoppler_ < 1.e-9)) {
std::cout << "--- Gets into if doppler factors are less than 1.e-9 ---" << std::endl;
- auto gridResolution_ {antenna.duration_};
- auto deltaT_ { endPointReceiveTime_ - startPointReceiveTime_ };
+ auto gridResolution_{antenna.duration_}; //TODO: maybe this should be antenna.sample_rate? (check)
+ auto deltaT_{endPointReceiveTime_ - startPointReceiveTime_};
if (std::fabs(deltaT_ / 1_s) < gridResolution_ / 1_s) {
EV1_ = EV1_ * std::fabs(deltaT_ / gridResolution_);
EV2_ = EV2_ * std::fabs(deltaT_ / gridResolution_);
- const long startBin = static_cast<long>(std::floor(startPointReceiveTime_/gridResolution_+0.5l));
- const long endBin = static_cast<long>(std::floor(endPointReceiveTime_/gridResolution_+0.5l));
- const double startBinFraction = (startPointReceiveTime_/gridResolution_)-std::floor(startPointReceiveTime_/gridResolution_);
- const double endBinFraction = (endPointReceiveTime_/gridResolution_)-std::floor(endPointReceiveTime_/gridResolution_);
+ const long startBin = static_cast<long>(
+ std::floor(startPointReceiveTime_ / gridResolution_ + 0.5l));
+ const long endBin = static_cast<long>(
+ std::floor(endPointReceiveTime_ / gridResolution_ + 0.5l));
+ const double startBinFraction =
+ (startPointReceiveTime_ / gridResolution_) -
+ std::floor(startPointReceiveTime_ / gridResolution_);
+ const double endBinFraction =
+ (endPointReceiveTime_ / gridResolution_) -
+ std::floor(endPointReceiveTime_ / gridResolution_);
// only do timing modification if contributions would land in same bin
if (startBin == endBin) {
// if startE arrives before endE
if (deltaT_ / 1_s >= 0) {
- if ((startBinFraction >= 0.5) && (endBinFraction >= 0.5)) // both points left of bin center
+ if ((startBinFraction >= 0.5) &&
+ (endBinFraction >= 0.5)) // both points left of bin center
{
- startPointReceiveTime_ = startPointReceiveTime_ - gridResolution_; // shift EV1_ to previous gridpoint
- }
- else if ((startBinFraction < 0.5) && (endBinFraction < 0.5)) // both points right of bin center
+ startPointReceiveTime_ =
+ startPointReceiveTime_ -
+ gridResolution_; // shift EV1_ to previous gridpoint
+ } else if ((startBinFraction < 0.5) &&
+ (endBinFraction < 0.5)) // both points right of bin center
{
- endPointReceiveTime_ = endPointReceiveTime_ + gridResolution_; // shift EV2_ to next gridpoint
- }
- else // points on both sides of bin center
+ endPointReceiveTime_ =
+ endPointReceiveTime_ +
+ gridResolution_; // shift EV2_ to next gridpoint
+ } else // points on both sides of bin center
{
- const double leftDist = 1.0-startBinFraction;
+ const double leftDist = 1.0 - startBinFraction;
const double rightDist = endBinFraction;
// check if asymmetry to right or left
- if (rightDist >= leftDist)
- {
- endPointReceiveTime_ = endPointReceiveTime_ + gridResolution_; // shift EV2_ to next gridpoint
- }
- else
- {
- startPointReceiveTime_ = startPointReceiveTime_ - gridResolution_; // shift EV1_ to previous gridpoint
+ if (rightDist >= leftDist) {
+ endPointReceiveTime_ =
+ endPointReceiveTime_ +
+ gridResolution_; // shift EV2_ to next gridpoint
+ } else {
+ startPointReceiveTime_ =
+ startPointReceiveTime_ -
+ gridResolution_; // shift EV1_ to previous gridpoint
}
}
- }
- else // if endE arrives before startE
+ } else // if endE arrives before startE
{
- if ((startBinFraction >= 0.5) && (endBinFraction >= 0.5)) // both points left of bin center
+ if ((startBinFraction >= 0.5) &&
+ (endBinFraction >= 0.5)) // both points left of bin center
{
- endPointReceiveTime_ = endPointReceiveTime_ - gridResolution_; // shift EV2_ to previous gridpoint
- }
- else if ((startBinFraction < 0.5) && (endBinFraction < 0.5)) // both points right of bin center
+ endPointReceiveTime_ =
+ endPointReceiveTime_ -
+ gridResolution_; // shift EV2_ to previous gridpoint
+ } else if ((startBinFraction < 0.5) &&
+ (endBinFraction < 0.5)) // both points right of bin center
{
- startPointReceiveTime_ = startPointReceiveTime_ + gridResolution_; // shift EV1_ to next gridpoint
- }
- else // points on both sides of bin center
+ startPointReceiveTime_ =
+ startPointReceiveTime_ +
+ gridResolution_; // shift EV1_ to next gridpoint
+ } else // points on both sides of bin center
{
- const double leftDist = 1.0-endBinFraction;
+ const double leftDist = 1.0 - endBinFraction;
const double rightDist = startBinFraction;
// check if asymmetry to right or left
- if (rightDist >= leftDist)
- {
- startPointReceiveTime_ = startPointReceiveTime_ + gridResolution_; // shift EV1_ to next gridpoint
- }
- else
- {
- endPointReceiveTime_ = endPointReceiveTime_ - gridResolution_; // shift EV2_ to previous gridpoint
+ if (rightDist >= leftDist) {
+ startPointReceiveTime_ =
+ startPointReceiveTime_ +
+ gridResolution_; // shift EV1_ to next gridpoint
+ } else {
+ endPointReceiveTime_ =
+ endPointReceiveTime_ -
+ gridResolution_; // shift EV2_ to previous gridpoint
}
}
} // End of else statement
- } // End of if for startbin == endbin
- } // End of if deltaT < gridresolution
- } // End of if that checks small doppler factors
+ } // End of if for startbin == endbin
+ } // End of if deltaT < gridresolution
+ } // End of if that checks small doppler factors
std::cout << "---------- RECEIVE INCIDENT WITH NO ZHS-like APPROXIMATION ----------" << std::endl;
std::cout << "Start Point Receive Time: " << startPointReceiveTime_ << std::endl;
@@ -371,9 +462,13 @@ namespace corsika {
} // End of else that does not perform ZHS-like approximation
} // End of loop over both paths to get signal info
+ } // End of try block
+ catch (size_t i) {
+ std::cout << " --- Signal Paths do not have the same size! --- " << std::endl;
+ }
} // End of looping over antennas
- return ProcessReturn::Ok;
+ return ProcessReturn::Ok;
} // End of simulate method
}; // END: class CoREAS
diff --git a/tests/modules/testRadio.cpp b/tests/modules/testRadio.cpp
index 4ae4358e3f46146bd27330279169dfdd7b29aba3..c38dfd47cd855fa404d5be7a7dcdde7987325450 100644
--- a/tests/modules/testRadio.cpp
+++ b/tests/modules/testRadio.cpp
@@ -269,6 +269,111 @@ TEST_CASE("Radio", "[processes]") {
}
+ SECTION("ZHS process") {
+
+ //////////////////////////////////////////////////////////////////////////////////////
+ // Environment
+ using IModelInterface = IRefractiveIndexModel<IMediumPropertyModel<IMagneticFieldModel<IMediumModel>>>;
+ using AtmModel = UniformRefractiveIndex<MediumPropertyModel<UniformMagneticField<HomogeneousMedium
+ <IModelInterface>>>>;
+ using EnvType = Environment<AtmModel>;
+ EnvType envZHS;
+ CoordinateSystemPtr const& rootCSZHS = envZHS.getCoordinateSystem();
+ // get the center point
+ Point const center{rootCSZHS, 0_m, 0_m, 0_m};
+ // a refractive index
+ const double ri_{1.000327};
+
+ // the constant density
+ const auto density{19.2_g / cube(1_cm)};
+
+ // the composition we use for the homogeneous medium
+ NuclearComposition const protonComposition(std::vector<Code>{Code::Proton},
+ std::vector<float>{1.f});
+
+ // create magnetic field vector
+ Vector B1(rootCSZHS, 0_T, 0_T, 1_T);
+
+ auto Medium = EnvType::createNode<Sphere>(
+ center, 1_km * std::numeric_limits<double>::infinity());
+
+ auto const props = Medium->setModelProperties<AtmModel>(ri_, Medium::AirDry1Atm, B1, density, protonComposition);
+ envZHS.getUniverse()->addChild(std::move(Medium));
+
+ // the antennas location
+ const auto point1{Point(envZHS.getCoordinateSystem(), 100_m, 2_m, 3_m)};
+ const auto point2{Point(envZHS.getCoordinateSystem(), 4_m, 80_m, 6_m)};
+ const auto point3{Point(envZHS.getCoordinateSystem(), 7_m, 8_m, 9_m)};
+ const auto point4{Point(envZHS.getCoordinateSystem(), 5_m, 5_m, 10_m)};
+
+ // create times for the antenna
+ const TimeType t1{0_s};
+ const TimeType t2{10_s};
+ const InverseTimeType t3{1e+3_Hz};
+ const TimeType t4{11_s};
+
+ // check that I can create an antenna at (1, 2, 3)
+ TimeDomainAntenna ant1("antenna_zhs", point1, t1, t2, t3);
+ TimeDomainAntenna ant2("antenna_zhs2", point2, t1, t2, t3);
+// TimeDomainAntenna ant3("antenna1", point1, 0_s, 2_s, 1/1e-7_s);
+
+ // construct a radio detector instance to store our antennas
+ AntennaCollection<TimeDomainAntenna> detector;
+
+ // add the antennas to the detector
+ detector.addAntenna(ant1);
+ detector.addAntenna(ant2);
+// detector.addAntenna(ant3);
+
+ // create a particle
+ auto const particle{Code::Electron};
+// auto const particle{Code::Gamma};
+ const auto pmass{get_mass(particle)};
+
+ VelocityVector v0(rootCSZHS, {5e+2_m / second, 5e+2_m / second, 5e+2_m / second});
+
+ Vector B0(rootCSZHS, 5_T, 5_T, 5_T);
+
+ Line const line(point3, v0);
+
+ auto const k{1_m * ((1_m) / ((1_s * 1_s) * 1_V))};
+
+ auto const t = 1_s;
+ LeapFrogTrajectory base(point4, v0, B0, k, t);
+
+ // create a new stack for each trial
+ setup::Stack stack;
+
+ // 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(rootCSZHS, {0_GeV, 0_GeV, P0})};
+
+ // and create the location of the particle in this coordinate system
+ const Point pos(rootCSZHS, 50_m, 10_m, 80_m);
+
+ // add the particle to the stack
+ auto const particle1{stack.addParticle(std::make_tuple(particle, E0, plab, pos, 0_ns))};
+
+ auto const charge_ {get_charge(particle1.getPID())};
+
+ // create a radio process instance using CoREAS
+ RadioProcess<decltype(detector), ZHS<decltype(detector), decltype(StraightPropagator(envZHS))>, decltype(StraightPropagator(envZHS))>
+ zhs(detector, envZHS);
+
+ // check doContinuous and simulate methods
+ zhs.doContinuous(particle1, base, true);
+// zhs.simulate(particle1, base);
+
+ // check writeOutput method -> should produce 2 csv files for each antenna
+ zhs.writeOutput();
+ }
+
+
SECTION("Synchrotron radiation") {
// create a suitable environment ///////////////////////////////////////////////////
@@ -910,7 +1015,7 @@ TEST_CASE("Radio", "[processes]") {
auto plab {beta * pmass * gamma};
Line l {point_1,v};
StraightTrajectory track {l,t};
- auto particle1{stack.addParticle(std::make_tuple(particle, E0, plab, point_1, t))}; //TODO: plab is inconsistent
+ auto particle1{stack.addParticle(std::make_tuple(particle, E0, plab, point_1, t))};
coreas.doContinuous(particle1,track,true);
stack.clear();
}