Replace Cascade.h

Framework/Cascade/Cascade.h

@@ -209,10 +209,23 @@ namespace corsika::cascade {
       // convert next_decay from time to length [m]
       LengthType const distance_decay = next_decay * vParticle.GetMomentum().norm() /
                                         vParticle.GetEnergy() * units::constants::c;
+                                    
+      // determine geometric tracking
+      auto [step, geomMaxLength, nextVol] = fTracking.GetTrack(vParticle);
+      [[maybe_unused]] auto const& dummy_nextVol = nextVol;
+      
+      // convert next_step from grammage to length
+      LengthType const distance_interact =
+          currentLogicalNode->GetModelProperties().ArclengthFromGrammage(step,
+                                                                         next_interact);
+      
+      // determine the maximum geometric step length
+      LengthType const distance_max = fProcessSequence.MaxStepLength(vParticle, step);
+      std::cout << "distance_max=" << distance_max << std::endl;
 
       // take minimum of geometry, interaction, decay for next step
       auto const min_distance =
-          std::min({distance_interact, distance_decay, distance_max, geomMaxLength});
+          std::min({distance_interact, distance_decay/*, distance_max, geomMaxLength*/});
 
       C8LOG_DEBUG("transport particle by : {} m", min_distance / 1_m);
 
@@ -220,6 +233,19 @@ namespace corsika::cascade {
       // std::visit(setup::ParticleUpdate<Particle>{vParticle}, step);
       vParticle.SetPosition(step.PositionFromArclength(min_distance));
       // .... also update time, momentum, direction, ...
+
+// magneticfield/u have wrong units, not sure if corsika::stack is right
+	  auto magneticfield = corsika::stack::MomentumVector(fEnvironment.GetCoordinateSystem(),
+	        {0_eV, 50e-06_eV, 0_eV});
+	  geometry::Vector<SpeedType::dimension_type> const velocity =
+            vParticle.GetMomentum() / vParticle.GetEnergy() * corsika::units::constants::c;
+      geometry::Vector<SpeedType::dimension_type> const u = velocity * 1_m / (velocity.GetNorm() * 1_s) + 
+        velocity.cross(magneticfield) * corsika::particles::GetChargeNumber(vParticle.GetPID()) * 
+            min_distance * corsika::units::constants::cSquared / 
+                (vParticle.GetEnergy() * velocity.GetNorm() * velocity.GetNorm()* 1_m);
+	  vParticle.SetMomentum(u * 1_s / 1_m * vParticle.GetMomentum().GetNorm());
+//      std::cout << " Test: " << magneticfield.GetComponents() << std::endl;
+//      std::cout << " Test: " << vParticle.GetMomentum().cross(magneticfield).GetComponents() << std::endl;
       vParticle.SetTime(vParticle.GetTime() + min_distance / units::constants::c);
 
       step.LimitEndTo(min_distance);