diff --git a/Documentation/Examples/cascade_example.cc b/Documentation/Examples/cascade_example.cc
index 0d4dc8c0f72b61dd15391ac4012196021a2e8e15..c0dcbdb7323704c658d897f57ca207a397fb4192 100644
--- a/Documentation/Examples/cascade_example.cc
+++ b/Documentation/Examples/cascade_example.cc
@@ -220,10 +220,7 @@ struct MyBoundaryCrossingProcess
: public BoundaryCrossingProcess<MyBoundaryCrossingProcess> {
//~ environment::BaseNodeType const& fA, fB;
- MyBoundaryCrossingProcess() {}
-
- //~ MyBoundaryCrossingProcess(environment::BaseNodeType const& a,
- // environment::BaseNodeType const& b) : fA(a), fB(b) {}
+ MyBoundaryCrossingProcess(std::string const& filename) {fFile.open(filename);}
template <typename Particle>
EProcessReturn DoBoundaryCrossing(Particle& p,
@@ -232,11 +229,18 @@ struct MyBoundaryCrossingProcess
std::cout << "boundary crossing! from:" << &from << "; to: " << &to << std::endl;
//~ if ((&fA == &from && &fB == &to) || (&fA == &to && &fB == &from)) {
- p.Delete();
+ //~ p.Delete();
//~ }
+
+ auto const& name = particles::GetName(p.GetPID());
+ auto const start = p.GetPosition().GetCoordinates();
+
+ fFile << name << " " << start[0] / 1_m << ' ' << start[1] / 1_m << ' ' << start[2] / 1_m << '\n';
return EProcessReturn::eOk;
}
+
+ std::ofstream fFile;
void Init() {}
};
@@ -282,7 +286,7 @@ int main() {
std::vector<float>{(float)1. - fox, fox}));
auto innerMedium = EnvType::CreateNode<Sphere>(
- Point{env.GetCoordinateSystem(), 0_m, 0_m, 0_m}, 2000_m);
+ Point{env.GetCoordinateSystem(), 0_m, 0_m, -500_m}, 1000_m);
innerMedium->SetModelProperties<
TheNameModel<environment::HomogeneousMedium<setup::IEnvironmentModel>>>(
@@ -312,12 +316,13 @@ int main() {
process::HadronicElasticModel::HadronicElasticInteraction hadronicElastic;
process::TrackWriter::TrackWriter trackWriter("tracks.dat");
+ MyBoundaryCrossingProcess boundaryCrossing("crossings.dat");
// assemble all processes into an ordered process list
// auto sequence = p0 << sibyll << decay << hadronicElastic << cut << trackWriter;
auto sequence =
- p0 /*<< sibyll << sibyllNuc << decay << cut*/ << hadronicElastic
- << MyBoundaryCrossingProcess()
+ /*p0 <<*/ sibyll << sibyllNuc << decay << cut /* << hadronicElastic */
+ << boundaryCrossing
<< trackWriter;
// setup particle stack, and add primary particle
@@ -328,10 +333,9 @@ int main() {
const int nuclZ = int(nuclA / 2.15 + 0.7);
const HEPMassType mass = particles::Proton::GetMass() * nuclZ +
(nuclA - nuclZ) * particles::Neutron::GetMass();
- const HEPEnergyType E0 = nuclA * 100_TeV;
- double theta = 0;
+ const HEPEnergyType E0 = nuclA * 10_TeV;
double phi = 0;
-
+ for (double theta = 0; theta < 180; theta += 20)
{
auto elab2plab = [](HEPEnergyType Elab, HEPMassType m) {
return sqrt(Elab * Elab - m * m);
diff --git a/Processes/TrackingLine/TrackingLine.cc b/Processes/TrackingLine/TrackingLine.cc
index 6b119c890aa573b4d62fea9b693472384e1bdfaf..319030f91d5347d47af58cddbdea607684a16677 100644
--- a/Processes/TrackingLine/TrackingLine.cc
+++ b/Processes/TrackingLine/TrackingLine.cc
@@ -44,8 +44,8 @@ namespace corsika::process::tracking_line {
if (discriminant.magnitude() > 0) {
auto const sqDisc = sqrt(discriminant);
auto const invDenom = 1 / vSqNorm;
- return std::make_pair((vDotDelta - sqDisc) * invDenom,
- (vDotDelta + sqDisc) * invDenom);
+ return std::make_pair((-vDotDelta - sqDisc) * invDenom,
+ (-vDotDelta + sqDisc) * invDenom);
} else {
return {};
}
diff --git a/Processes/TrackingLine/TrackingLine.h b/Processes/TrackingLine/TrackingLine.h
index ff3b178aa4ee597e3a13598e3514fb1c2e7340e7..001297004ed185d12094b7d0127ba61d6b2ff4cf 100644
--- a/Processes/TrackingLine/TrackingLine.h
+++ b/Processes/TrackingLine/TrackingLine.h
@@ -12,8 +12,6 @@
#ifndef _include_corsika_processes_TrackingLine_h_
#define _include_corsika_processes_TrackingLine_h_
-//~ #include <corsika/environment/Environment.h>
-//~ #include <corsika/environment/VolumeTreeNode.h>
#include <corsika/units/PhysicalUnits.h>
#include <corsika/geometry/Vector.h>
#include <corsika/geometry/Trajectory.h>
@@ -27,10 +25,6 @@ namespace corsika::environment {
template <typename IEnvironmentModel> class Environment;
template <typename IEnvironmentModel> class VolumeTreeNode;
}
-namespace corsika::geometry {
- class Line;
- class Sphere;
-} // namespace corsika::geometry
namespace corsika::process {
@@ -70,14 +64,16 @@ namespace corsika::process {
//~ fEnvironment.GetUniverse()->GetContainingNode(currentPosition);
auto const numericallyInside =
currentLogicalVolumeNode->GetVolume().Contains(currentPosition);
+
+ std::cout << "numericallyInside = " << (numericallyInside ? "true":"false");
auto const& children = currentLogicalVolumeNode->GetChildNodes();
auto const& excluded = currentLogicalVolumeNode->GetExcludedNodes();
std::vector<std::pair<TimeType, decltype(p.GetNode())>> intersections;
- auto addIfIntersects = [&](auto const& vtn, auto const& nextNode) {
- static_assert(std::is_same_v<decltype(vtn), decltype(nextNode)>);
+ // for entering from outside
+ auto addIfIntersects = [&](auto const& vtn) {
auto const& volume = vtn.GetVolume();
auto const& sphere = dynamic_cast<geometry::Sphere const&>(
volume); // for the moment we are a bit bold here and assume
@@ -88,14 +84,14 @@ namespace corsika::process {
std::cout << "intersection times: " << t1 / 1_s << "; " << t2 / 1_s
<< std::endl;
if (t1.magnitude() > 0)
- intersections.emplace_back(t1, &nextNode);
+ intersections.emplace_back(t1, &vtn);
else if (t2.magnitude() > 0)
throw std::runtime_error("inside other volume");
}
};
- for (auto const& child : children) { addIfIntersects(*child, *child); }
- for (auto const* ex : excluded) { addIfIntersects(*ex, *ex); }
+ for (auto const& child : children) { addIfIntersects(*child); }
+ for (auto const* ex : excluded) { addIfIntersects(*ex); }
{
auto const& sphere = dynamic_cast<geometry::Sphere const&>(
@@ -120,7 +116,8 @@ namespace corsika::process {
min = minIter->first;
}
- std::cout << " t-intersect: " << min << std::endl;
+ std::cout << " t-intersect: " << min << " " << minIter->second->GetModelProperties().GetName() << std::endl;
+ std::cout << "point of intersection: " << line.GetPosition(min).GetCoordinates() << std::endl;
return std::make_tuple(geometry::Trajectory<geometry::Line>(line, min), velocity.norm() * min,
minIter->second);
diff --git a/Processes/TrackingLine/testTrackingLine.cc b/Processes/TrackingLine/testTrackingLine.cc
index 26dce27d178bae5bb7c7cb4f48b0f4d2c58df060..fb213f9d8e3729015a8ad2f19c092e3d476bf9dc 100644
--- a/Processes/TrackingLine/testTrackingLine.cc
+++ b/Processes/TrackingLine/testTrackingLine.cc
@@ -41,7 +41,7 @@ TEST_CASE("TrackingLine") {
tracking_line::TrackingLine tracking;
SECTION("intersection with sphere") {
- Point const origin(cs, {0_m, 0_m, 0_m});
+ Point const origin(cs, {0_m, 0_m, -5_m});
Point const center(cs, {0_m, 0_m, 10_m});
Sphere const sphere(center, 1_m);
Vector<corsika::units::si::SpeedType::dimension_type> v(cs, 0_m / second,
@@ -55,8 +55,8 @@ TEST_CASE("TrackingLine") {
REQUIRE(opt.has_value());
auto [t1, t2] = opt.value();
- REQUIRE(t1 / 9_s == Approx(1));
- REQUIRE(t2 / 11_s == Approx(1));
+ REQUIRE(t1 / 14_s == Approx(1));
+ REQUIRE(t2 / 16_s == Approx(1));
auto const optNoIntersection =
tracking_line::TimeOfIntersection(traj, Sphere(Point(cs, {5_m, 0_m, 10_m}), 1_m));
@@ -70,6 +70,7 @@ TEST_CASE("TrackingLine") {
auto theMedium = environment::Environment<environment::Empty>::CreateNode<Sphere>(
Point{env.GetCoordinateSystem(), 0_m, 0_m, 0_m}, radius);
+ auto const* theMediumPtr = theMedium.get();
universe.AddChild(std::move(theMedium));
TestTrackingLineStack stack;
@@ -82,6 +83,7 @@ TEST_CASE("TrackingLine") {
{cs, {0_m, 0_m, 0_km}},
0_ns});
auto p = stack.GetNextParticle();
+ p.SetNode(theMediumPtr);
Point const origin(cs, {0_m, 0_m, 0_m});
Vector<corsika::units::si::SpeedType::dimension_type> v(cs, 0_m / second,