preliminary working boundary crossings

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);

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 {};
     }

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);

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,