From 2e94dec739bb9d3480cb9b3cf916dbc76f8d7323 Mon Sep 17 00:00:00 2001
From: Nikos Karastathis <n.karastathis@kit.edu>
Date: Wed, 14 Apr 2021 15:36:21 +0200
Subject: [PATCH] Synchrotron radiation fixed momentum
---
corsika/modules/radio/CoREAS.hpp | 6 +++---
tests/modules/testRadio.cpp | 30 ++++++++++++++----------------
2 files changed, 17 insertions(+), 19 deletions(-)
diff --git a/corsika/modules/radio/CoREAS.hpp b/corsika/modules/radio/CoREAS.hpp
index 032b41c9e..1cbd1e5ff 100755
--- a/corsika/modules/radio/CoREAS.hpp
+++ b/corsika/modules/radio/CoREAS.hpp
@@ -96,13 +96,13 @@ namespace corsika {
// calculate preDoppler factor
double preDoppler_{1. - paths1[i].refractive_index_source_ *
- beta_.dot(paths1[i].receive_)};
+ beta_.dot(paths1[i].emit_)}; // TODO: are you sure this is path.receive and not emit?
std::cout << "***** preDoppler: " << preDoppler_ << std::endl;
// calculate postDoppler factor
double postDoppler_{
1. - paths2[i].refractive_index_source_ *
- beta_.dot(paths2[i].receive_)}; // maybe this is path.receive_ (?)
+ beta_.dot(paths2[i].emit_)}; // maybe this is path.receive_ (?)
std::cout << "***** postDoppler: " << postDoppler_ << std::endl;
// calculate receive time for startpoint
@@ -147,7 +147,7 @@ namespace corsika {
for (auto const& path : paths3) {
auto const midPointReceiveTime_{path.propagation_time_ + midTime_};
- auto midDoppler_{1. - path.refractive_index_source_ * beta_.dot(path.receive_)};
+ auto midDoppler_{1. - path.refractive_index_source_ * beta_.dot(path.emit_)};
// change the values of the receive unit vectors of start and end
ReceiveVectorStart_ = path.receive_;
diff --git a/tests/modules/testRadio.cpp b/tests/modules/testRadio.cpp
index 2abeedaf3..7d5ff6e31 100644
--- a/tests/modules/testRadio.cpp
+++ b/tests/modules/testRadio.cpp
@@ -776,22 +776,18 @@ TEST_CASE("Radio", "[processes]") {
std::vector<TimeType> times_;
- // create a particle /////////////////////////////////////////////////////////////
- auto const particle{Code::Electron};
- const auto pmass{get_mass(particle)};
+ //////////////////////////////////////////////////////////////////////////////////
// create a new stack for each trial
setup::Stack stack;
+ stack.clear();
+
+ const Code particle{Code::Electron};
+ const HEPMassType pmass{get_mass(particle)};
// construct an energy // move in the for loop
const HEPEnergyType E0{11.4_MeV};
- // compute the necessary momentum // move in the for loop
- const HEPMomentumType P0{sqrt(E0 * E0 - pmass * pmass)};
-
- // and create the momentum vector
- const auto plab{MomentumVector(rootCS, {0_GeV, 0_GeV, P0})};
-
// create a radio process instance using CoREAS
RadioProcess<decltype(detector), CoREAS<decltype(detector), decltype(StraightPropagator(env))>, decltype(StraightPropagator(env))>
coreas(detector, env);
@@ -800,29 +796,31 @@ TEST_CASE("Radio", "[processes]") {
for (size_t i = 1; i <= 399; i++) {
TimeType t {(points_[i] - points_[i-1]).getNorm() / (0.999 * constants::c)};
VelocityVector v { (points_[i] - points_[i-1]) / t };
+ auto beta {v / constants::c};
+ auto gamma {E0/pmass};
+ auto plab {beta * pmass * gamma};
Line l {points_[i-1],v};
StraightTrajectory track {l,t};
- auto const particle1{stack.addParticle(std::make_tuple(particle, E0, plab, points_[i-1], t))}; //TODO: plab is inconsistent
+ auto particle1{stack.addParticle(std::make_tuple(particle, E0, plab, points_[i-1], t))}; //TODO: plab is inconsistent
coreas.doContinuous(particle1,track,true);
}
+ // ToDo: use just one electron, this way I have 400! (although it shouldn't affect the physics)
// get the last track
TimeType t {(points_[0] - points_[399]).getNorm() / (0.999 * constants::c)};
VelocityVector v { (points_[0] - points_[399]) / t };
+ auto beta {v / constants::c};
+ auto gamma {E0/pmass};
+ auto plab {beta * pmass * gamma};
Line l {points_[399],v};
StraightTrajectory track {l,t};
- auto const particle1{stack.addParticle(std::make_tuple(particle, E0, plab, points_[399], t))};
+ auto particle1{stack.addParticle(std::make_tuple(particle, E0, plab, points_[399], t))};
coreas.doContinuous(particle1,track,true);
// get the output
coreas.writeOutput();
- // VelocityVector v0 {(p1 - p0) / 1_s}; //v = ((x1 - x2) + (y1 - y2)) / dt
-// Line l1{p0,v0};
-// StraightTrajectory st0 {l1,1_s};
-// std::cout << times_.size() << std::endl;
-
}
--
GitLab