diff --git a/corsika/modules/radio/CoREAS.hpp b/corsika/modules/radio/CoREAS.hpp index 032b41c9e0e3d06cfba88a4fac5f205c6ffe9289..1cbd1e5ffb6fff4a7db2d774cd5335f3d5e33df4 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 2abeedaf351dc82f46446eece9712f174faef336..7d5ff6e31786ea5d686f47a44129bc0f034d50e2 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; - }