IAP GITLAB

Commit 3f3225f1 authored by André Schmidt's avatar André Schmidt

Merge branch 'andre-master' into 'master'

added function for magnetic step

See merge request !13
parents a5f10ee3 84b1b737
Pipeline #2038 failed with stages
in 0 seconds
......@@ -196,7 +196,7 @@ namespace corsika::cascade {
vParticle.GetEnergy() * units::constants::c;
// determine geometric tracking
auto [stepWithoutB, stepWithB, geomMaxLength, nextVol] = fTracking.GetTrack(vParticle);
auto [lineWithoutB, stepWithoutB, stepWithB, geomMaxLength, nextVol] = fTracking.GetTrack(vParticle);
[[maybe_unused]] auto const& dummy_nextVol = nextVol;
// convert next_step from grammage to length
......@@ -212,45 +212,20 @@ namespace corsika::cascade {
auto min_distance = std::min(
{distance_interact, distance_decay, distance_max, geomMaxLength});
std::cout << "Decay after " << distance_decay << std::endl;
std::cout << "Interaction after " << distance_interact << std::endl;
std::cout << "Decay after " << distance_decay << std::endl;
std::cout << "Interaction after " << distance_interact << std::endl;
std::cout << " move particle by : " << min_distance << std::endl;
// determine displacement by the magnetic field
auto const* currentLogicalVolumeNode = vParticle.GetNode();
auto magneticfield = currentLogicalVolumeNode->GetModelProperties().GetMagneticField(vParticle.GetPosition());
geometry::Vector<SpeedType::dimension_type> velocity = vParticle.GetMomentum() / vParticle.GetEnergy() *
corsika::units::constants::c;
geometry::Vector<dimensionless_d> const directionBefore = velocity.normalized();
int chargeNumber;
if (corsika::particles::IsNucleus(vParticle.GetPID())) {
chargeNumber = vParticle.GetNuclearZ();
} else {
chargeNumber = corsika::particles::GetChargeNumber(vParticle.GetPID());
}
auto k = chargeNumber * corsika::units::constants::cSquared * 1_eV /
(velocity.GetNorm() * vParticle.GetEnergy() * 1_V);
// First Movement
// assuming magnetic field does not change during movement
auto position = vParticle.GetPosition() + directionBefore * min_distance / 2;
// Change of direction by magnetic field
geometry::Vector<dimensionless_d> const directionAfter = directionBefore + directionBefore.cross(magneticfield) *
min_distance * k;
// Second Movement
position = position + directionAfter * min_distance / 2;
// determine displacement by the magnetic field
auto [line, position, directionAfter] = fTracking.MagneticStep(vParticle, lineWithoutB, min_distance);
auto distance = position - vParticle.GetPosition();
// distance.norm() != min_distance if q != 0
// small error can be neglected
if (distance.norm() != 0_m) {
velocity = distance.normalized() * velocity.norm();
} // no velocity update for very small steps
// here the particle is actually moved along the trajectory to new position:
// std::visit(setup::ParticleUpdate<Particle>{vParticle}, step);
vParticle.SetMomentum(directionAfter.normalized() * vParticle.GetMomentum().norm());
geometry::Line line(vParticle.GetPosition(), velocity);
geometry::Trajectory<geometry::Line> stepNew(line, distance.norm() / velocity.norm());
geometry::Trajectory<geometry::Line> stepNew(line, distance.norm() / line.GetV0().norm());
vParticle.SetPosition(position);
vParticle.SetTime(vParticle.GetTime() + distance.norm() / units::constants::c);
std::cout << "New Position: " << vParticle.GetPosition().GetCoordinates() << std::endl;
......
......@@ -65,6 +65,8 @@ namespace corsika::process {
<< " GeV " << std::endl;
std::cout << "TrackingLine v: " << velocity.GetComponents() << std::endl;
geometry::Line line(currentPosition, velocity);
auto const* currentLogicalVolumeNode = p.GetNode();
//~ auto const* currentNumericalVolumeNode =
//~ fEnvironment.GetUniverse()->GetContainingNode(currentPosition);
......@@ -89,8 +91,7 @@ namespace corsika::process {
std::cout << " Magnetic Field: " << magneticfield.GetComponents() / 1_uT << " uT " << std::endl;
auto k = chargeNumber * corsika::units::constants::cSquared * 1_eV / (velocity.GetNorm() * p.GetEnergy() * 1_V);
geometry::Vector<dimensionless_d> const directionBefore = velocity.normalized();
geometry::Vector<SpeedType::dimension_type> velocity1 = velocity;
geometry::Vector<SpeedType::dimension_type> velocity2 = velocity;
LengthType Steplength = 10_m; // length irrelevant if q=0 and else it gets changed again
// for entering from outside
auto addIfIntersects = [&](auto const& vtn) {
......@@ -117,7 +118,6 @@ namespace corsika::process {
std::cout << "Solutions for next Volume: " << solutions[i].real() << std::endl;
}
}
LengthType Steplength;
if (tmp.size() > 0) {
Steplength = 1_m * *std::min_element(tmp.begin(),tmp.end());
std::cout << "Steplength to next volume = " << Steplength << std::endl;
......@@ -126,23 +126,15 @@ namespace corsika::process {
// what to do when this happens? (very unlikely)
}
delete [] solutions;
// First Movement
// assuming magnetic field does not change during movement
auto position = currentPosition + directionBefore * Steplength / 2;
// Change of direction by magnetic field
geometry::Vector<dimensionless_d> const directionAfter = directionBefore + directionBefore.cross(magneticfield) *
Steplength * k;
// Second Movement
position = position + directionAfter * Steplength / 2;
geometry::Vector<dimensionless_d> const direction = (position - currentPosition) /
(position - currentPosition).GetNorm();
velocity1 = direction * velocity.GetNorm();
} // instead of changing the line with magnetic field, the TimeOfIntersection() could be changed
}
auto [line1, position, direction] = MagneticStep(p, line, Steplength);
// new particle position and direction not needed in this case
// instead of changing the line with magnetic field, the TimeOfIntersection() could be changed
// using line has some errors for huge steps
geometry::Line line(currentPosition, velocity1);
if (auto opt = TimeOfIntersection(line, sphere); opt.has_value()) {
if (auto opt = TimeOfIntersection(line1, sphere); opt.has_value()) {
auto const [t1, t2] = *opt;
std::cout << "intersection times: " << t1 / 1_s << "; "
<< t2 / 1_s
......@@ -184,34 +176,26 @@ namespace corsika::process {
}
LengthType Steplength;
if (tmp.size() > 0) {
Steplength = 1_m * *std::min_element(tmp.begin(),tmp.end());
if (numericallyInside == false) {
int p = std::min_element(tmp.begin(),tmp.end()) - tmp.begin();
tmp.erase(tmp.begin() + p);
Steplength = 1_m * *std::min_element(tmp.begin(),tmp.end());
}
std::cout << "steplength out of current volume = " << Steplength << std::endl;
Steplength = 1_m * *std::min_element(tmp.begin(),tmp.end());
if (numericallyInside == false) {
int p = std::min_element(tmp.begin(),tmp.end()) - tmp.begin();
tmp.erase(tmp.begin() + p);
Steplength = 1_m * *std::min_element(tmp.begin(),tmp.end());
}
std::cout << "steplength out of current volume = " << Steplength << std::endl;
} else {
std::cout << "no intersection (2)!" << std::endl;
// what to do when this happens? (very unlikely)
}
delete [] solutions;
}
// First Movement
// assuming magnetic field does not change during movement
auto position = currentPosition + directionBefore * Steplength / 2;
// Change of direction by magnetic field
geometry::Vector<dimensionless_d> const directionAfter = directionBefore + directionBefore.cross(magneticfield) *
Steplength * k;
// Second Movement
position = position + directionAfter * Steplength / 2;
geometry::Vector<dimensionless_d> const direction = (position - currentPosition) /
(position - currentPosition).GetNorm();
velocity2 = direction * velocity.GetNorm();
} // instead of changing the line with magnetic field, the TimeOfIntersection() could be changed
geometry::Line line(currentPosition, velocity2);
auto [line2, position, direction] = MagneticStep(p, line, Steplength);
// new particle position and direction not needed in this case
// instead of changing the line with magnetic field, the TimeOfIntersection() could be changed
// using line has some errors for huge steps
[[maybe_unused]] auto const [t1, t2] = *TimeOfIntersection(line, sphere);
[[maybe_unused]] auto const [t1, t2] = *TimeOfIntersection(line2, sphere);
[[maybe_unused]] auto dummy_t1 = t1;
intersections.emplace_back(t2, currentLogicalVolumeNode->GetParent());
}
......@@ -234,26 +218,49 @@ namespace corsika::process {
<< min
// << " " << minIter->second->GetModelProperties().GetName()
<< std::endl;
auto [lineWithB, position, direction] = MagneticStep(p, line, velocity.norm() * min);
// new particle position and direction not needed in this case
return std::make_tuple(line, geometry::Trajectory<geometry::Line>(line, min),
geometry::Trajectory<geometry::Line>(lineWithB, min),
velocity.norm() * min, minIter->second);
}
template <typename Particle> // was Stack previously, and argument was
// Stack::StackIterator
// determine direction of the particle after adding magnetic field
auto MagneticStep(Particle const& p, corsika::geometry::Line line, corsika::units::si::LengthType Steplength) {
using namespace corsika::units::si;
//charge of the particle
int chargeNumber;
if (corsika::particles::IsNucleus(p.GetPID())) {
chargeNumber = p.GetNuclearZ();
} else {
chargeNumber = corsika::particles::GetChargeNumber(p.GetPID());
}
auto const* currentLogicalVolumeNode = p.GetNode();
auto magneticfield = currentLogicalVolumeNode->GetModelProperties().GetMagneticField(p.GetPosition());
auto k = chargeNumber * corsika::units::constants::cSquared * 1_eV / (line.GetV0().norm() * p.GetEnergy() * 1_V);
geometry::Vector<dimensionless_d> const directionBefore = line.GetV0().normalized();
geometry::Line lineWithoutB(currentPosition, velocity);
// determine direction of the particle after adding magnetic field
// assuming magnetic field does not change during movement
// First Movement
auto position = currentPosition + velocity * min / 2;
// assuming magnetic field does not change during movement
auto position = p.GetPosition() + directionBefore * Steplength / 2;
// Change of direction by magnetic field
geometry::Vector<dimensionless_d> const directionAfter = directionBefore + velocity.cross(magneticfield) *
min * k;
geometry::Vector<dimensionless_d> const directionAfter = directionBefore + directionBefore.cross(magneticfield) *
Steplength * k;
// Second Movement
position = position + directionAfter * velocity.norm() * min / 2;
if ((position - currentPosition).GetNorm() != 0_m) {
geometry::Vector<dimensionless_d> const direction = (position - currentPosition).normalized();
velocity = direction * velocity.norm();
} // no velocity update for very small steps
geometry::Line lineWithB(currentPosition, velocity);
return std::make_tuple(geometry::Trajectory<geometry::Line>(lineWithoutB, min),
geometry::Trajectory<geometry::Line>(lineWithB, min),
velocity.norm() * min, minIter->second);
position = position + directionAfter * Steplength / 2;
auto distance = position - p.GetPosition();
if(distance.norm() == 0_m) {
return std::make_tuple(line, position, directionAfter);
}
geometry::Line newLine(p.GetPosition(), distance.normalized() * line.GetV0().norm());
return std::make_tuple(newLine, position, directionAfter);
}
};
......
Markdown is supported
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