diff --git a/corsika/detail/modules/tracking/TrackingLeapFrogCurved.inl b/corsika/detail/modules/tracking/TrackingLeapFrogCurved.inl
index 1691b04862996a915d2534bf2ce57e554c31a541..5be8e516ef9549e7f58551964a1b68a3cb0db963 100644
--- a/corsika/detail/modules/tracking/TrackingLeapFrogCurved.inl
+++ b/corsika/detail/modules/tracking/TrackingLeapFrogCurved.inl
@@ -225,81 +225,84 @@ namespace corsika {
return Intersections(d_enter / absVelocity, d_exit / absVelocity);
}
-
template <typename TParticle, typename TMedium>
Intersections Tracking::intersect(TParticle const& particle, Plane const& plane,
- TMedium const& medium) {
-
+ TMedium const& medium) {
+
int chargeNumber;
if (is_nucleus(particle.getPID())) {
- chargeNumber = particle.getNuclearZ();
+ chargeNumber = particle.getNuclearZ();
} else {
- chargeNumber = get_charge_number(particle.getPID());
+ chargeNumber = get_charge_number(particle.getPID());
}
auto const* currentLogicalVolumeNode = particle.getNode();
- auto magneticfield = currentLogicalVolumeNode->getModelProperties().getMagneticField(
- particle.getPosition());
+ auto magneticfield =
+ currentLogicalVolumeNode->getModelProperties().getMagneticField(
+ particle.getPosition());
VelocityVector const velocity =
- particle.getMomentum() / particle.getEnergy() * constants::c;
+ particle.getMomentum() / particle.getEnergy() * constants::c;
auto const absVelocity = velocity.getNorm();
DirectionVector const direction =
- velocity.normalized(); // determine steplength to next volume
+ velocity.normalized(); // determine steplength to next volume
Point const position = particle.getPosition();
-
+
if (chargeNumber != 0 &&
- abs(plane.getNormal().dot(velocity.cross(magneticfield))) *
- 1_s / 1_m / 1_T > 1e-6) {
-
- auto const* currentLogicalVolumeNode = particle.getNode();
- auto magneticfield = currentLogicalVolumeNode->getModelProperties().getMagneticField(
- particle.getPosition());
- auto k = chargeNumber * constants::c * 1_eV /
- (particle.getMomentum().getNorm() * 1_V);
-
- if (direction.dot(plane.getNormal()) * direction.dot(plane.getNormal()) -
- (plane.getNormal().dot(position - plane.getCenter()) *
- plane.getNormal().dot(direction.cross(magneticfield)) * 2 * k) <
- 0) {
- return Intersections(std::numeric_limits<double>::infinity() * 1_s);
- }
-
- LengthType const MaxStepLength1 =
- (sqrt(direction.dot(plane.getNormal()) * direction.dot(plane.getNormal()) -
- (plane.getNormal().dot(position - plane.getCenter()) *
- plane.getNormal().dot(direction.cross(magneticfield)) * 2 * k)) -
- direction.dot(plane.getNormal()) / direction.getNorm()) /
- (plane.getNormal().dot(direction.cross(magneticfield)) * k);
-
- LengthType const MaxStepLength2 =
- (-sqrt(
- direction.dot(plane.getNormal()) * direction.dot(plane.getNormal()) -
- (plane.getNormal().dot(position - plane.getCenter()) *
- plane.getNormal().dot(direction.cross(magneticfield)) * 2 * k)) -
- direction.dot(plane.getNormal()) / direction.getNorm()) /
- (plane.getNormal().dot(direction.cross(magneticfield)) * k);
-
- if (MaxStepLength1 <= 0_m && MaxStepLength2 <= 0_m) {
- return Intersections(std::numeric_limits<double>::infinity() * 1_s);
- } else if (MaxStepLength1 <= 0_m || MaxStepLength2 < MaxStepLength1) {
- CORSIKA_LOG_TRACE(" steplength to obs plane 2: {} ", MaxStepLength2);
- return Intersections(MaxStepLength2 *
- (direction + direction.cross(magneticfield) * MaxStepLength2 * k / 2)
- .getNorm() / absVelocity);
- } else if (MaxStepLength2 <= 0_m || MaxStepLength1 < MaxStepLength2) {
- CORSIKA_LOG_TRACE(" steplength to obs plane 2: {} ", MaxStepLength1);
- return Intersections(MaxStepLength1 *
- (direction + direction.cross(magneticfield) * MaxStepLength2 * k / 2)
- .getNorm() / absVelocity);
- }
-
- CORSIKA_LOG_WARN("Particle wasn't tracked with curved trajectory -> straight");
-
+ abs(plane.getNormal().dot(velocity.cross(magneticfield))) * 1_s / 1_m / 1_T >
+ 1e-6) {
+
+ auto const* currentLogicalVolumeNode = particle.getNode();
+ auto magneticfield =
+ currentLogicalVolumeNode->getModelProperties().getMagneticField(
+ particle.getPosition());
+ auto k =
+ chargeNumber * constants::c * 1_eV / (particle.getMomentum().getNorm() * 1_V);
+
+ if (direction.dot(plane.getNormal()) * direction.dot(plane.getNormal()) -
+ (plane.getNormal().dot(position - plane.getCenter()) *
+ plane.getNormal().dot(direction.cross(magneticfield)) * 2 * k) <
+ 0) {
+ return Intersections(std::numeric_limits<double>::infinity() * 1_s);
+ }
+
+ LengthType const MaxStepLength1 =
+ (sqrt(direction.dot(plane.getNormal()) * direction.dot(plane.getNormal()) -
+ (plane.getNormal().dot(position - plane.getCenter()) *
+ plane.getNormal().dot(direction.cross(magneticfield)) * 2 * k)) -
+ direction.dot(plane.getNormal()) / direction.getNorm()) /
+ (plane.getNormal().dot(direction.cross(magneticfield)) * k);
+
+ LengthType const MaxStepLength2 =
+ (-sqrt(direction.dot(plane.getNormal()) * direction.dot(plane.getNormal()) -
+ (plane.getNormal().dot(position - plane.getCenter()) *
+ plane.getNormal().dot(direction.cross(magneticfield)) * 2 * k)) -
+ direction.dot(plane.getNormal()) / direction.getNorm()) /
+ (plane.getNormal().dot(direction.cross(magneticfield)) * k);
+
+ if (MaxStepLength1 <= 0_m && MaxStepLength2 <= 0_m) {
+ return Intersections(std::numeric_limits<double>::infinity() * 1_s);
+ } else if (MaxStepLength1 <= 0_m || MaxStepLength2 < MaxStepLength1) {
+ CORSIKA_LOG_TRACE(" steplength to obs plane 2: {} ", MaxStepLength2);
+ return Intersections(
+ MaxStepLength2 *
+ (direction + direction.cross(magneticfield) * MaxStepLength2 * k / 2)
+ .getNorm() /
+ absVelocity);
+ } else if (MaxStepLength2 <= 0_m || MaxStepLength1 < MaxStepLength2) {
+ CORSIKA_LOG_TRACE(" steplength to obs plane 2: {} ", MaxStepLength1);
+ return Intersections(
+ MaxStepLength1 *
+ (direction + direction.cross(magneticfield) * MaxStepLength2 * k / 2)
+ .getNorm() /
+ absVelocity);
+ }
+
+ CORSIKA_LOG_WARN("Particle wasn't tracked with curved trajectory -> straight");
+
} // end if curved-tracking
return tracking_line::Tracking::intersect(particle, plane, medium);
}
-
template <typename TParticle, typename TBaseNodeType>
inline Intersections Tracking::intersect(const TParticle& particle,
const TBaseNodeType& volumeNode) {
diff --git a/corsika/framework/stack/SecondaryView.hpp b/corsika/framework/stack/SecondaryView.hpp
index e18531b3154e250ab44c33014a92355e0b9a877a..e15179c7e8ade3da16d62e268845bfa5cbb0c8b5 100644
--- a/corsika/framework/stack/SecondaryView.hpp
+++ b/corsika/framework/stack/SecondaryView.hpp
@@ -189,6 +189,16 @@ namespace corsika {
// NOTE: 0 is special marker here for PROJECTILE, see getIndexFromIterator
return stack_view_iterator(*this, 0);
}
+
+ /**
+ * This return a projectile of this SecondaryView, which can be
+ * used to modify the SecondaryView
+ */
+ inline const_stack_view_iterator getProjectile() const {
+ // NOTE: 0 is special marker here for PROJECTILE, see getIndexFromIterator
+ return const_stack_view_iterator(*this, 0);
+ }
+
/**
* Method to add a new secondary particle on this SecondaryView
*/
diff --git a/corsika/modules/tracking/TrackingLeapFrogCurved.hpp b/corsika/modules/tracking/TrackingLeapFrogCurved.hpp
index 6fe3af378ac7a6bcde93c425931b7a1dadc2dd2f..05073691eafacb54ce8150ce73ae8a940c3b83eb 100644
--- a/corsika/modules/tracking/TrackingLeapFrogCurved.hpp
+++ b/corsika/modules/tracking/TrackingLeapFrogCurved.hpp
@@ -70,7 +70,7 @@ namespace corsika {
template <typename TParticle, typename TMedium>
static Intersections intersect(TParticle const& particle, Plane const& plane,
- TMedium const& medium);
+ TMedium const& medium);
protected:
/**
diff --git a/src/framework/core/ParticleData.xml b/src/framework/core/ParticleData.xml
index 40cd64c8a4e9b5dc03a8f87e32e308f16e3a3b94..db89d04628abe4613e0c4a673575bddb550aed9e 100644
--- a/src/framework/core/ParticleData.xml
+++ b/src/framework/core/ParticleData.xml
@@ -166,7 +166,7 @@
<channel onMode="1" bRatio="0.1081660" products="-13 14"/>
<channel onMode="1" bRatio="0.1080870" products="-15 16"/>
</particle>
-
+ -->
<particle id="25" name="h0" spinType="1" chargeType="0" colType="0"
m0="125.00000" mWidth="0.00374" mMin="50.00000" mMax="0.00000">
<channel onMode="1" bRatio="0.0000009" products="1 -1"/>
@@ -246,7 +246,7 @@
<channel onMode="1" bRatio="0.0000000" meMode="103" products="1000016 -2000016"/>
<channel onMode="1" bRatio="0.0000000" meMode="103" products="-1000016 2000016"/>
</particle>
-
+<!--
<particle id="32" name="Z'0" spinType="3" chargeType="0" colType="0"
m0="500.00000" mWidth="14.54029" mMin="10.00000" mMax="0.00000">
<channel onMode="1" bRatio="0.1458350" products="1 -1"/>
diff --git a/tests/common/SetupTestTrajectory.hpp b/tests/common/SetupTestTrajectory.hpp
index a33be655e0d9d28255738cc4df83f131cbb66810..8272317193a1d1864b5c7928c22947dcc7243e2f 100644
--- a/tests/common/SetupTestTrajectory.hpp
+++ b/tests/common/SetupTestTrajectory.hpp
@@ -17,13 +17,14 @@
//#include <corsika/modules/TrackingLine.hpp>
//#include <corsika/modules/TrackingCurved.hpp> // simple leap-frog implementation
//#include <corsika/modules/TrackingLeapFrog.hpp> // more complete leap-frog
- // implementation
+// implementation
namespace corsika::setup::testing {
template <typename TTrack>
- TTrack make_track(Line const line,
- TimeType const tEnd = std::numeric_limits<TimeType::value_type>::infinity() * 1_s);
+ TTrack make_track(
+ Line const line,
+ TimeType const tEnd = std::numeric_limits<TimeType::value_type>::infinity() * 1_s);
template <>
inline StraightTrajectory make_track<StraightTrajectory>(Line const line,
@@ -33,13 +34,13 @@ namespace corsika::setup::testing {
template <>
inline LeapFrogTrajectory make_track<LeapFrogTrajectory>(Line const line,
- TimeType const tEnd) {
+ TimeType const tEnd) {
auto const k = square(0_m) / (square(1_s) * 1_V);
return LeapFrogTrajectory(
- line.getStartPoint(), line.getVelocity(),
- MagneticFieldVector{line.getStartPoint().getCoordinateSystem(), 0_T, 0_T, 0_T},
- k, tEnd);
- }
-
+ line.getStartPoint(), line.getVelocity(),
+ MagneticFieldVector{line.getStartPoint().getCoordinateSystem(), 0_T, 0_T, 0_T}, k,
+ tEnd);
+ }
+
} // namespace corsika::setup::testing
diff --git a/tests/framework/testCascade.cpp b/tests/framework/testCascade.cpp
index 9070f0d89dae82f04fb608fe43440021f95f1db0..633f30f47328970be8374552ecbb86411aaa6936 100644
--- a/tests/framework/testCascade.cpp
+++ b/tests/framework/testCascade.cpp
@@ -73,8 +73,8 @@ public:
Line const theLine = Line(particle.getPosition(), initialVelocity);
TimeType const tEnd = std::numeric_limits<TimeType::value_type>::infinity() * 1_s;
return std::make_tuple(
- corsika::setup::testing::make_track<setup::Trajectory>(theLine, tEnd),
- // trajectory: just go ahead forever
+ corsika::setup::testing::make_track<setup::Trajectory>(theLine, tEnd),
+ // trajectory: just go ahead forever
particle.getNode()); // next volume node
}
};
diff --git a/tests/media/testMedium.cpp b/tests/media/testMedium.cpp
index d303c820f50b0a65f9b6d9800dc4224a6972dc6d..edcba45fc6fe01a3b5554878fcca735627b7a9d7 100644
--- a/tests/media/testMedium.cpp
+++ b/tests/media/testMedium.cpp
@@ -101,7 +101,8 @@ TEST_CASE("MediumPropertyModel w/ Homogeneous") {
auto const tEnd = 1_s;
// and the associated trajectory
- setup::Trajectory const trajectory = corsika::setup::testing::make_track<setup::Trajectory>(line, tEnd);
+ setup::Trajectory const trajectory =
+ corsika::setup::testing::make_track<setup::Trajectory>(line, tEnd);
// and check the integrated grammage
CHECK((medium.getIntegratedGrammage(trajectory, 3_m) / (density * 3_m)) == Approx(1));