diff --git a/Framework/Cascade/Cascade.h b/Framework/Cascade/Cascade.h index 8779ed513a9fa10b3e0b9b077f373858573a1313..3a6df6da491912d30ecaa0a74494cb9c15c2beac 100644 --- a/Framework/Cascade/Cascade.h +++ b/Framework/Cascade/Cascade.h @@ -212,70 +212,77 @@ 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())) { + 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 = corsika::geometry::Vector - (fEnvironment.GetCoordinateSystem(), 0_uT, 50_uT, 0_uT); - geometry::Vector<SpeedType::dimension_type> const velocityVerticalMag = - velocity - magneticfield * velocity.dot(magneticfield) / (magneticfield.GetSquaredNorm()); - 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 = corsika::geometry::Vector(fEnvironment.GetCoordinateSystem(), 0_uT, 50_uT, 0_uT); - 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; - } - + } 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 = corsika::geometry::Vector(fEnvironment.GetCoordinateSystem(), + 0_uT, 50_uT, 0_uT); + geometry::Vector<SpeedType::dimension_type> const velocityVerticalMag = + velocity - magneticfield * velocity.dot(magneticfield) / + (magneticfield.GetSquaredNorm()); + 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 = corsika::geometry::Vector(fEnvironment.GetCoordinateSystem(), + 0_uT, 50_uT, 0_uT); + 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); [[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);