IAP GITLAB

Skip to content
Snippets Groups Projects
Commit 98d72bf5 authored by Maximilian Reininghaus's avatar Maximilian Reininghaus :vulcan: Committed by ralfulrich
Browse files

clang-format

parent 3ad101ec
No related branches found
No related tags found
1 merge request!278Magnetic Tracking
Hide whitespace changes
Inline Side-by-side
Framework/Cascade/Cascade.h
+ 50
− 43
View file @ 98d72bf5
...@@ -212,70 +212,77 @@ namespace corsika::cascade { ...@@ -212,70 +212,77 @@ namespace corsika::cascade {
// convert next_decay from time to length [m] // convert next_decay from time to length [m]
LengthType const distance_decay = next_decay * vParticle.GetMomentum().norm() / LengthType const distance_decay = next_decay * vParticle.GetMomentum().norm() /
vParticle.GetEnergy() * units::constants::c; vParticle.GetEnergy() * units::constants::c;
// determine momentum after adding magnetic field // determine momentum after adding magnetic field
int chargeNumber; int chargeNumber;
if(corsika::particles::IsNucleus(vParticle.GetPID())) { if (corsika::particles::IsNucleus(vParticle.GetPID())) {
chargeNumber = vParticle.GetNuclearZ(); chargeNumber = vParticle.GetNuclearZ();
} else { } else {
chargeNumber = corsika::particles::GetChargeNumber(vParticle.GetPID()); chargeNumber = corsika::particles::GetChargeNumber(vParticle.GetPID());
} }
geometry::Vector<SpeedType::dimension_type> const velocity = geometry::Vector<SpeedType::dimension_type> const velocity =
vParticle.GetMomentum() / vParticle.GetEnergy() * corsika::units::constants::c; vParticle.GetMomentum() / vParticle.GetEnergy() * corsika::units::constants::c;
geometry::Vector<dimensionless_d> const directionBefore = velocity / velocity.GetNorm(); geometry::Vector<dimensionless_d> const directionBefore =
auto magMaxLength = 1_m/0; velocity / velocity.GetNorm();
auto directionAfter = directionBefore; auto magMaxLength = 1_m / 0;
if(chargeNumber != 0) { auto directionAfter = directionBefore;
auto magneticfield = corsika::geometry::Vector if (chargeNumber != 0) {
(fEnvironment.GetCoordinateSystem(), 0_uT, 50_uT, 0_uT); auto magneticfield = corsika::geometry::Vector(fEnvironment.GetCoordinateSystem(),
geometry::Vector<SpeedType::dimension_type> const velocityVerticalMag = 0_uT, 50_uT, 0_uT);
velocity - magneticfield * velocity.dot(magneticfield) / (magneticfield.GetSquaredNorm()); geometry::Vector<SpeedType::dimension_type> const velocityVerticalMag =
LengthType const gyroradius = vParticle.GetEnergy() * velocityVerticalMag.GetNorm() * 1_V / velocity - magneticfield * velocity.dot(magneticfield) /
(corsika::units::constants::cSquared * abs(chargeNumber) * (magneticfield.GetSquaredNorm());
magneticfield.GetNorm() * 1_eV); LengthType const gyroradius =
//steplength depending on how exact it should vParticle.GetEnergy() * velocityVerticalMag.GetNorm() * 1_V /
LengthType const Steplength = 0.01 * gyroradius; (corsika::units::constants::cSquared * abs(chargeNumber) *
// First Movement magneticfield.GetNorm() * 1_eV);
auto position = vParticle.GetPosition() + directionBefore * Steplength / 2; // steplength depending on how exact it should
// Change of direction by magnetic field at position LengthType const Steplength = 0.01 * gyroradius;
magneticfield = corsika::geometry::Vector(fEnvironment.GetCoordinateSystem(), 0_uT, 50_uT, 0_uT); // First Movement
directionAfter = directionBefore + velocity.cross(magneticfield) * chargeNumber * auto position = vParticle.GetPosition() + directionBefore * Steplength / 2;
Steplength * corsika::units::constants::cSquared * 1_eV / // Change of direction by magnetic field at position
(vParticle.GetEnergy() * velocity.GetSquaredNorm() * 1_V); magneticfield = corsika::geometry::Vector(fEnvironment.GetCoordinateSystem(),
// Second Movement 0_uT, 50_uT, 0_uT);
position = position + directionAfter * Steplength / 2; directionAfter = directionBefore +
magMaxLength = (position - vParticle.GetPosition()).GetNorm(); velocity.cross(magneticfield) * chargeNumber * Steplength *
geometry::Vector<dimensionless_d> const direction = (position - vParticle.GetPosition()) / corsika::units::constants::cSquared * 1_eV /
magMaxLength; (vParticle.GetEnergy() * velocity.GetSquaredNorm() * 1_V);
vParticle.SetMomentum( direction * vParticle.GetMomentum().GetNorm()); // Second Movement
std::cout << "New direction because of magnetic field: " << direction.GetComponents() << std::endl; 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 // determine geometric tracking
auto [step, geomMaxLength, nextVol] = fTracking.GetTrack(vParticle); auto [step, geomMaxLength, nextVol] = fTracking.GetTrack(vParticle);
[[maybe_unused]] auto const& dummy_nextVol = nextVol; [[maybe_unused]] auto const& dummy_nextVol = nextVol;
// convert next_step from grammage to length // convert next_step from grammage to length
LengthType const distance_interact = LengthType const distance_interact =
currentLogicalNode->GetModelProperties().ArclengthFromGrammage(step, currentLogicalNode->GetModelProperties().ArclengthFromGrammage(step,
next_interact); next_interact);
// determine the maximum geometric step length // determine the maximum geometric step length
LengthType const distance_max = fProcessSequence.MaxStepLength(vParticle, step); LengthType const distance_max = fProcessSequence.MaxStepLength(vParticle, step);
std::cout << "distance_max=" << distance_max << std::endl; std::cout << "distance_max=" << distance_max << std::endl;
// take minimum of geometry, interaction, decay for next step // take minimum of geometry, interaction, decay for next step
auto const min_distance = auto const min_distance = std::min(
std::min({distance_interact, distance_decay, distance_max, geomMaxLength, magMaxLength}); {distance_interact, distance_decay, distance_max, geomMaxLength, magMaxLength});
C8LOG_DEBUG("transport particle by : {} m", min_distance / 1_m); C8LOG_DEBUG("transport particle by : {} m", min_distance / 1_m);
// here the particle is actually moved along the trajectory to new position: // here the particle is actually moved along the trajectory to new position:
// std::visit(setup::ParticleUpdate<Particle>{vParticle}, step); // std::visit(setup::ParticleUpdate<Particle>{vParticle}, step);
vParticle.SetPosition(step.PositionFromArclength(min_distance)); vParticle.SetPosition(step.PositionFromArclength(min_distance));
// .... also update time, momentum, direction, ... // .... also update time, momentum, direction, ...
vParticle.SetMomentum((directionBefore * (1 - min_distance / magMaxLength) + vParticle.SetMomentum((directionBefore * (1 - min_distance / magMaxLength) +
directionAfter * min_distance /magMaxLength) * vParticle.GetMomentum().GetNorm()); directionAfter * min_distance / magMaxLength) *
vParticle.GetMomentum().GetNorm());
vParticle.SetTime(vParticle.GetTime() + min_distance / units::constants::c); vParticle.SetTime(vParticle.GetTime() + min_distance / units::constants::c);
step.LimitEndTo(min_distance); step.LimitEndTo(min_distance);
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment

data privacy - impressum