diff --git a/Framework/Cascade/Cascade.h b/Framework/Cascade/Cascade.h index 20bfa7348b2517df8892634e1e780726c5b092aa..41e20931fd25bab27dce1474b8693ef9e4192476 100644 --- a/Framework/Cascade/Cascade.h +++ b/Framework/Cascade/Cascade.h @@ -212,75 +212,32 @@ 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 momentum after adding magnetic field - int chargeNumber; - if (corsika::particles::IsNucleus(vParticle.GetPID())) { - chargeNumber = vParticle.GetNuclearZ(); - } else { - chargeNumber = corsika::particles::GetChargeNumber(vParticle.GetPID()); - } - geometry::Vector<SpeedType::dimension_type> const velocity = - vParticle.GetMomentum() / vParticle.GetEnergy() * corsika::units::constants::c; - geometry::Vector<dimensionless_d> const directionBefore = - velocity / velocity.GetNorm(); - auto magMaxLength = 1_m / 0; - auto directionAfter = directionBefore; - if (chargeNumber != 0) { - auto magneticfield = currentLogicalNode->GetModelProperties(). - GetMagneticField(vParticle.GetPosition()); - geometry::Vector<SpeedType::dimension_type> const velocityVerticalMag = - velocity - velocity.parallelProjectionOnto(magneticfield); - LengthType const gyroradius = - vParticle.GetEnergy() * velocityVerticalMag.GetNorm() * 1_V / - (corsika::units::constants::cSquared * abs(chargeNumber) * - magneticfield.GetNorm() * 1_eV); - // steplength depending on how exact it should - LengthType const Steplength = 0.01 * gyroradius; - // First Movement - auto position = vParticle.GetPosition() + directionBefore * Steplength / 2; - // Change of direction by magnetic field at position - magneticfield = currentLogicalNode->GetModelProperties.GetMagneticField(position); - directionAfter = directionBefore + - velocity.cross(magneticfield) * chargeNumber * Steplength * - corsika::units::constants::cSquared * 1_eV / - (vParticle.GetEnergy() * velocity.GetSquaredNorm() * 1_V); - // Second Movement - position = position + directionAfter * Steplength / 2; - magMaxLength = (position - vParticle.GetPosition()).GetNorm(); - geometry::Vector<dimensionless_d> const direction = - (position - vParticle.GetPosition()) / magMaxLength; - vParticle.SetMomentum(direction * vParticle.GetMomentum().GetNorm()); - std::cout << "New direction because of magnetic field: " - << direction.GetComponents() << std::endl; - } - + // determine geometric tracking - auto [step, geomMaxLength, nextVol] = fTracking.GetTrack(vParticle); + auto [step, geomMaxLength, nextVol, magMaxLength, directionBefore, directionAfter] = 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, magMaxLength}); + auto const min_distance = + std::min({distance_interact, distance_decay, distance_max, geomMaxLength, magMaxLength}); C8LOG_DEBUG("transport particle by : {} m", min_distance / 1_m); // here the particle is actually moved along the trajectory to new position: // std::visit(setup::ParticleUpdate<Particle>{vParticle}, step); vParticle.SetPosition(step.PositionFromArclength(min_distance)); - // .... also update time, momentum, direction, ... - vParticle.SetMomentum((directionBefore * (1 - min_distance / magMaxLength) + - directionAfter * min_distance / magMaxLength) * - vParticle.GetMomentum().GetNorm()); + // .... also update time, momentum, direction, ... + vParticle.SetMomentum((directionBefore * (1 - min_distance / magMaxLength) + + directionAfter * min_distance /magMaxLength) * vParticle.GetMomentum().GetNorm()); vParticle.SetTime(vParticle.GetTime() + min_distance / units::constants::c); step.LimitEndTo(min_distance);