update Magnetic movement

Documentation/Examples/vertical_EAS.cc

@@ -97,7 +97,7 @@ int main(int argc, char** argv) {
   node1->SetModelProperties<UniformMagneticField<HomogeneousMedium<IEnvModel>>>(
       B, 1_g / (1_cm * 1_cm * 1e9_cm), composition);
 
-  /*auto node2 = std::make_unique<VolumeTreeNode<IEnvModel>>(
+  auto node2 = std::make_unique<VolumeTreeNode<IEnvModel>>(
       std::make_unique<geometry::Sphere>(center, seaLevel + 100.0_km));
   node2->SetModelProperties<UniformMagneticField<SlidingPlanarExponential<IEnvModel>>>(
       B, center, 540.1778_g / (1_cm * 1_cm) / 772170.16_cm, -772170.16_cm, composition,
@@ -124,7 +124,7 @@ int main(int argc, char** argv) {
   node4->AddChild(std::move(node5));
   node3->AddChild(std::move(node4));
   node2->AddChild(std::move(node3));
-  node1->AddChild(std::move(node2));*/
+  node1->AddChild(std::move(node2));
   env.GetUniverse()->AddChild(std::move(node1));
 
   // setup particle stack, and add primary particle

Framework/Cascade/Cascade.h

@@ -209,6 +209,10 @@ namespace corsika::cascade {
       std::cout << "distance_max=" << distance_max << std::endl;
 
       // take minimum of geometry, interaction, decay for next step
+      std::cout << "Interaction: " << distance_interact << std::endl;
+      std::cout << "Decay: " << distance_decay << std::endl;
+      std::cout << "ObsPlane: " << distance_max << std::endl;
+      std::cout << "Transition: " << geomMaxLength << std::endl;
       auto const min_distance = std::min(
           {distance_interact, distance_decay, distance_max, geomMaxLength});
 
@@ -216,44 +220,43 @@ namespace corsika::cascade {
 
       // determine displacement by the magnetic field
 	    auto const* currentLogicalVolumeNode = vParticle.GetNode();
+      auto magneticfield = currentLogicalVolumeNode->GetModelProperties().GetMagneticField(vParticle.GetPosition());
+      geometry::Vector<SpeedType::dimension_type> velocity = vParticle.GetMomentum() / vParticle.GetEnergy() *
+                                                              corsika::units::constants::c;
+  		geometry::Vector<dimensionless_d> const directionBefore = velocity.normalized();
       int chargeNumber;
       if (corsika::particles::IsNucleus(vParticle.GetPID())) {
         chargeNumber = vParticle.GetNuclearZ();
       } else {
      	  chargeNumber = corsika::particles::GetChargeNumber(vParticle.GetPID());
       }
-      if (chargeNumber != 0) {
-  		  auto magneticfield = currentLogicalVolumeNode->GetModelProperties().GetMagneticField(vParticle.GetPosition());
-    		geometry::Vector<SpeedType::dimension_type> velocity = vParticle.GetMomentum() / vParticle.GetEnergy() *
-    			                                                     corsika::units::constants::c;
-    		geometry::Vector<dimensionless_d> const directionBefore = velocity.normalized();
-    		auto k = chargeNumber * corsika::units::constants::cSquared * 1_eV / 
-                (velocity.GetNorm() * vParticle.GetEnergy() * 1_V);
-    		// First Movement
-    		// assuming magnetic field does not change during movement
-    		auto position = vParticle.GetPosition() + directionBefore * min_distance / 2;
-    		// Change of direction by magnetic field
-    		geometry::Vector<dimensionless_d> const directionAfter = directionBefore + directionBefore.cross(magneticfield) *
-                                                                min_distance * k; 
-    		// Second Movement
-    		position = position + directionAfter * min_distance / 2;
-        // here the particle is actually moved along the trajectory to new position:
-        // std::visit(setup::ParticleUpdate<Particle>{vParticle}, step);
-        vParticle.SetMomentum(directionAfter.normalized() * vParticle.GetMomentum().GetNorm());
-        vParticle.SetPosition(position);
-      } else {
-        vParticle.SetPosition(stepWithoutB.PositionFromArclength(min_distance));
-      }
-      // .... also update time, momentum, direction, ...  
-      vParticle.SetTime(vParticle.GetTime() + min_distance / units::constants::c);
+      auto k = chargeNumber * corsika::units::constants::cSquared * 1_eV / 
+              (velocity.GetNorm() * vParticle.GetEnergy() * 1_V);
+		  
+  		// First Movement
+  		// assuming magnetic field does not change during movement
+  		auto position = vParticle.GetPosition() + directionBefore * min_distance / 2;
+  		// Change of direction by magnetic field
+  		geometry::Vector<dimensionless_d> const directionAfter = directionBefore + directionBefore.cross(magneticfield) *
+                                                              min_distance * k; 
+  		// Second Movement
+  		position = position + directionAfter * min_distance / 2;
+      auto distance = position - vParticle.GetPosition();
+      //distance.norm() != min_distance for distance_interact, distance_decay if q != 0
+      //small error can be neglected
+      velocity = distance.normalized() * velocity.norm();
+      
+      // here the particle is actually moved along the trajectory to new position:
+      // std::visit(setup::ParticleUpdate<Particle>{vParticle}, step);
+      vParticle.SetMomentum(directionAfter.normalized() * vParticle.GetMomentum().GetNorm());
+      geometry::Line line(vParticle.GetPosition(), velocity);
+      geometry::Trajectory<geometry::Line> stepNew(line, distance.norm() / velocity.GetNorm());
+      vParticle.SetPosition(position);
+      vParticle.SetTime(vParticle.GetTime() + distance.norm() / units::constants::c);
       std::cout << "New Position: " << vParticle.GetPosition().GetCoordinates() << std::endl;
 
-      // is this necessary?
-      stepWithoutB.LimitEndTo(min_distance);
-      stepWithB.LimitEndTo(min_distance);
-
       // apply all continuous processes on particle + track
-      process::EProcessReturn status = fProcessSequence.DoContinuous(vParticle, stepWithoutB);
+      process::EProcessReturn status = fProcessSequence.DoContinuous(vParticle, stepNew);
 
       if (status == process::EProcessReturn::eParticleAbsorbed) {
         std::cout << "Cascade: delete absorbed particle " << vParticle.GetPID() << " "

Processes/ObservationPlane/ObservationPlane.cc

@@ -82,25 +82,22 @@ LengthType ObservationPlane::MaxStepLength(setup::Stack::ParticleType const& vPa
       plane_.GetNormal().dot(velocity.cross(magneticfield)) * 2 * k)) - 
       velocity.dot(plane_.GetNormal()) / velocity.GetNorm() ) / 
       (plane_.GetNormal().dot(velocity.cross(magneticfield)) * k);
-    std::cout << "Test: " << MaxStepLength1 << " " << MaxStepLength2 << std::endl;
-    if (MaxStepLength1 < 0_m && MaxStepLength2 < 0_m) {
-      return std::numeric_limits<double>::infinity() * 1_m;
-    } else if (MaxStepLength1 < 0_m || MaxStepLength2 < MaxStepLength1) {
-      return MaxStepLength2;
-    } else if (MaxStepLength2 < 0_m || MaxStepLength1 < MaxStepLength2) {
-      return MaxStepLength1;
-    }
-  } else {
-    TimeType const timeOfIntersection =
-      (plane_.GetCenter() - trajectory.GetR0()).dot(plane_.GetNormal()) /
-      trajectory.GetV0().dot(plane_.GetNormal());
-
-    if (timeOfIntersection < TimeType::zero()) {
+    if (MaxStepLength1 <= 0_m && MaxStepLength2 <= 0_m) {
       return std::numeric_limits<double>::infinity() * 1_m;
+    } else if (MaxStepLength1 <= 0_m || MaxStepLength2 < MaxStepLength1) {
+      return MaxStepLength2 * 1.0001;
+    } else if (MaxStepLength2 <= 0_m || MaxStepLength1 < MaxStepLength2) {
+      return MaxStepLength1 * 1.0001;
     }
+  } 
+  TimeType const timeOfIntersection =
+    (plane_.GetCenter() - trajectory.GetR0()).dot(plane_.GetNormal()) /
+    trajectory.GetV0().dot(plane_.GetNormal());
 
-    auto const pointOfIntersection = trajectory.GetPosition(timeOfIntersection);
-    return (trajectory.GetR0() - pointOfIntersection).norm() * 1.0001;
-    //why is it *1.0001? should i do that too?
+  if (timeOfIntersection < TimeType::zero()) {
+    return std::numeric_limits<double>::infinity() * 1_m;
   }
+
+  auto const pointOfIntersection = trajectory.GetPosition(timeOfIntersection);
+  return (trajectory.GetR0() - pointOfIntersection).norm() * 1.0001;
 }