diff --git a/Framework/Cascade/Cascade.h b/Framework/Cascade/Cascade.h index 43a5275435480e45aaa3bb06ac824f2ba5636f40..db628600737f260e007e3acde8eca63ae2b2943e 100644 --- a/Framework/Cascade/Cascade.h +++ b/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);