IAP GITLAB

Skip to content
Snippets Groups Projects
Commit 8bf837b8 authored by ralfulrich's avatar ralfulrich
Browse files

Max's comments, and a bug...

parent 8d07b2b5
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
Showing
with 108 additions and 126 deletions
corsika/detail/modules/ObservationPlane.inl
+ 2
− 6
View file @ 8bf837b8
......@@ -59,13 +59,9 @@ namespace corsika {
corsika::setup::Stack::particle_type const& particle,
corsika::setup::Trajectory const& trajectory) {
auto const& volumeNode = particle.getNode();
typedef typename std::remove_const_t<
std::remove_reference_t<decltype(volumeNode->getModelProperties())>>
medium_type;
Intersections const intersection =
setup::Tracking::intersect<corsika::setup::Stack::particle_type, medium_type>(
particle, plane_, volumeNode->getModelProperties());
setup::Tracking::intersect<corsika::setup::Stack::particle_type>(particle,
plane_);
TimeType const timeOfIntersection = intersection.getEntry();
CORSIKA_LOG_TRACE("particle={}, pos={}, dir={}, plane={}, timeOfIntersection={}",
particle.asString(), particle.getPosition(),
......
corsika/detail/modules/tracking/Intersect.inl
+ 8
− 8
View file @ 8bf837b8
......@@ -38,7 +38,7 @@ namespace corsika {
// first check, where we leave the current volume
// this assumes our convention, that all Volume primitives must be convex
// thus, the last entry is always the exit point
const Intersections time_intersections_curr =
Intersections const time_intersections_curr =
TDerived::intersect(particle, volumeNode);
CORSIKA_LOG_TRACE("curr node {}, parent node {}, hasIntersections={} ",
fmt::ptr(&volumeNode), fmt::ptr(volumeNode.getParent()),
......@@ -56,16 +56,16 @@ namespace corsika {
// where do we collide with any of the next-tree-level volumes
// entirely contained by currentLogicalVolumeNode
for (const auto& node : volumeNode.getChildNodes()) {
for (auto const& node : volumeNode.getChildNodes()) {
const Intersections time_intersections = TDerived::intersect(particle, *node);
Intersections const time_intersections = TDerived::intersect(particle, *node);
if (!time_intersections.hasIntersections()) { continue; }
CORSIKA_LOG_DEBUG("intersection times with child volume {} : enter {} s, exit {} s",
fmt::ptr(node), time_intersections.getEntry() / 1_s,
time_intersections.getExit() / 1_s);
const auto t_entry = time_intersections.getEntry();
const auto t_exit = time_intersections.getExit();
auto const t_entry = time_intersections.getEntry();
auto const t_exit = time_intersections.getExit();
CORSIKA_LOG_TRACE("children t-entry: {}, t-exit: {}, smaller? {} ", t_entry, t_exit,
t_entry <= minTime);
// note, theoretically t can even be smaller than 0 since we
......@@ -86,14 +86,14 @@ namespace corsika {
// current volume
for (node_type* node : volumeNode.getExcludedNodes()) {
const Intersections time_intersections = TDerived::intersect(particle, *node);
Intersections const time_intersections = TDerived::intersect(particle, *node);
if (!time_intersections.hasIntersections()) { continue; }
CORSIKA_LOG_DEBUG(
"intersection times with exclusion volume {} : enter {} s, exit {} s",
fmt::ptr(node), time_intersections.getEntry() / 1_s,
time_intersections.getExit() / 1_s);
const auto t_entry = time_intersections.getEntry();
const auto t_exit = time_intersections.getExit();
auto const t_entry = time_intersections.getEntry();
auto const t_exit = time_intersections.getExit();
CORSIKA_LOG_TRACE("children t-entry: {}, t-exit: {}, smaller? {} ", t_entry, t_exit,
t_entry <= minTime);
// note, theoretically t can even be smaller than 0 since we
......
corsika/detail/modules/tracking/TrackingLeapFrogCurved.inl
+ 59
− 66
View file @ 8bf837b8
......@@ -103,9 +103,9 @@ namespace corsika {
LengthType const gyroradius =
(pAlongB_delta * 1_V / (constants::c * abs(chargeNumber) * magnitudeB * 1_eV));
const double maxRadians = 0.01;
const LengthType steplimit = 2 * cos(maxRadians) * sin(maxRadians) * gyroradius;
const TimeType steplimit_time = steplimit / initialVelocity.getNorm();
double const maxRadians = 0.01;
LengthType const steplimit = 2 * cos(maxRadians) * sin(maxRadians) * gyroradius;
TimeType const steplimit_time = steplimit / initialVelocity.getNorm();
CORSIKA_LOG_DEBUG("gyroradius {}, steplimit: {} = {}", gyroradius, steplimit,
steplimit_time);
......@@ -113,7 +113,7 @@ namespace corsika {
// intersection
auto [minTime, minNode] = nextIntersect(particle, steplimit_time);
const auto k =
auto const k =
chargeNumber * constants::cSquared * 1_eV / (particle.getEnergy() * 1_V);
return std::make_tuple(
LeapFrogTrajectory(position, initialVelocity, magneticfield, k,
......@@ -121,10 +121,9 @@ namespace corsika {
minNode); // next volume node
}
template <typename TParticle, typename TMedium>
inline Intersections Tracking::intersect(const TParticle& particle,
const Sphere& sphere,
const TMedium& medium) {
template <typename TParticle>
inline Intersections Tracking::intersect(TParticle const& particle,
Sphere const& sphere) {
if (sphere.getRadius() == 1_km * std::numeric_limits<double>::infinity()) {
return Intersections();
......@@ -132,7 +131,9 @@ namespace corsika {
int const chargeNumber = particle.getChargeNumber();
auto const& position = particle.getPosition();
MagneticFieldVector const& magneticfield = medium.getMagneticField(position);
auto const* currentLogicalVolumeNode = particle.getNode();
MagneticFieldVector const& magneticfield =
currentLogicalVolumeNode->getModelProperties().getMagneticField(position);
VelocityVector const velocity =
particle.getMomentum() / particle.getEnergy() * constants::c;
......@@ -144,26 +145,23 @@ namespace corsika {
bool const isParallel = (projectedDirectionSqrNorm == 0 * square(1_T));
if (chargeNumber == 0 || magneticfield.getNorm() == 0_T || isParallel) {
return tracking_line::Tracking::intersect<TParticle, TMedium>(particle, sphere,
medium);
return tracking_line::Tracking::intersect<TParticle>(particle, sphere);
}
bool const numericallyInside = sphere.contains(particle.getPosition());
const auto absVelocity = velocity.getNorm();
auto const absVelocity = velocity.getNorm();
auto energy = particle.getEnergy();
auto k = chargeNumber * constants::cSquared * 1_eV / (absVelocity * energy * 1_V);
auto const denom = (directionBefore.cross(magneticfield)).getSquaredNorm() * k * k;
const double a =
((directionBefore.cross(magneticfield)).dot(position - sphere.getCenter()) * k +
1) *
4 / (1_m * 1_m * denom);
const double b = directionBefore.dot(position - sphere.getCenter()) * 8 /
(denom * 1_m * 1_m * 1_m);
const double c = ((position - sphere.getCenter()).getSquaredNorm() -
(sphere.getRadius() * sphere.getRadius())) *
4 / (denom * 1_m * 1_m * 1_m * 1_m);
auto const direction_B_perp = directionBefore.cross(magneticfield);
auto const denom = direction_B_perp.getSquaredNorm() * k * k;
Vector<length_d> const deltaPos = position - sphere.getCenter();
double const a = (direction_B_perp.dot(deltaPos) * k + 1) * 4 / (1_m * 1_m * denom);
double const b = directionBefore.dot(deltaPos) * 8 / (denom * 1_m * 1_m * 1_m);
double const c =
(deltaPos.getSquaredNorm() - (sphere.getRadius() * sphere.getRadius())) * 4 /
(denom * 1_m * 1_m * 1_m * 1_m);
CORSIKA_LOG_TRACE("denom={}, a={}, b={}, c={}", denom, a, b, c);
std::complex<double>* solutions = quartic_solver::solve_quartic(0, a, b, c);
LengthType d_enter, d_exit;
......@@ -225,9 +223,9 @@ namespace corsika {
return Intersections(d_enter / absVelocity, d_exit / absVelocity);
}
template <typename TParticle, typename TMedium>
template <typename TParticle>
inline Intersections Tracking::intersect(TParticle const& particle,
Plane const& plane, TMedium const& medium) {
Plane const& plane) {
int chargeNumber;
if (is_nucleus(particle.getPID())) {
......@@ -236,9 +234,6 @@ namespace corsika {
chargeNumber = get_charge_number(particle.getPID());
}
auto const* currentLogicalVolumeNode = particle.getNode();
auto magneticfield =
currentLogicalVolumeNode->getModelProperties().getMagneticField(
particle.getPosition());
VelocityVector const velocity =
particle.getMomentum() / particle.getEnergy() * constants::c;
auto const absVelocity = velocity.getNorm();
......@@ -246,70 +241,68 @@ namespace corsika {
velocity.normalized(); // determine steplength to next volume
Point const position = particle.getPosition();
if (chargeNumber != 0 &&
abs(plane.getNormal().dot(velocity.cross(magneticfield))) * 1_s / 1_m / 1_T >
1e-6) {
auto const magneticfield =
currentLogicalVolumeNode->getModelProperties().getMagneticField(position);
if (chargeNumber != 0 && abs(plane.getNormal().dot(velocity.cross(magneticfield))) >
1e-6_T * 1_m / 1_s) {
auto const* currentLogicalVolumeNode = particle.getNode();
auto magneticfield =
currentLogicalVolumeNode->getModelProperties().getMagneticField(
particle.getPosition());
auto k =
chargeNumber * constants::c * 1_eV / (particle.getMomentum().getNorm() * 1_V);
if (direction.dot(plane.getNormal()) * direction.dot(plane.getNormal()) -
(plane.getNormal().dot(position - plane.getCenter()) *
plane.getNormal().dot(direction.cross(magneticfield)) * 2 * k) <
0) {
auto const magneticfield =
currentLogicalVolumeNode->getModelProperties().getMagneticField(position);
auto const k =
chargeNumber * (constants::c * 1_eV / 1_V) / particle.getMomentum().getNorm();
auto const direction_B_perp = direction.cross(magneticfield);
auto const denom = plane.getNormal().dot(direction_B_perp) * k;
auto const sqrtArg =
direction.dot(plane.getNormal()) * direction.dot(plane.getNormal()) -
(plane.getNormal().dot(position - plane.getCenter()) * denom * 2);
if (sqrtArg < 0) {
return Intersections(std::numeric_limits<double>::infinity() * 1_s);
}
double const sqrSqrtArg = sqrt(sqrtArg);
auto const norm_projected =
direction.dot(plane.getNormal()) / direction.getNorm();
LengthType const MaxStepLength1 = (sqrSqrtArg - norm_projected) / denom;
LengthType const MaxStepLength2 = (-sqrSqrtArg - norm_projected) / denom;
LengthType const MaxStepLength1 =
(sqrt(direction.dot(plane.getNormal()) * direction.dot(plane.getNormal()) -
(plane.getNormal().dot(position - plane.getCenter()) *
plane.getNormal().dot(direction.cross(magneticfield)) * 2 * k)) -
direction.dot(plane.getNormal()) / direction.getNorm()) /
(plane.getNormal().dot(direction.cross(magneticfield)) * k);
LengthType const MaxStepLength2 =
(-sqrt(direction.dot(plane.getNormal()) * direction.dot(plane.getNormal()) -
(plane.getNormal().dot(position - plane.getCenter()) *
plane.getNormal().dot(direction.cross(magneticfield)) * 2 * k)) -
direction.dot(plane.getNormal()) / direction.getNorm()) /
(plane.getNormal().dot(direction.cross(magneticfield)) * k);
// check: both intersections in past
if (MaxStepLength1 <= 0_m && MaxStepLength2 <= 0_m) {
return Intersections(std::numeric_limits<double>::infinity() * 1_s);
// check: next intersection is MaxStepLength2
} else if (MaxStepLength1 <= 0_m || MaxStepLength2 < MaxStepLength1) {
CORSIKA_LOG_TRACE(" steplength to obs plane 2: {} ", MaxStepLength2);
return Intersections(
MaxStepLength2 *
(direction + direction.cross(magneticfield) * MaxStepLength2 * k / 2)
.getNorm() /
(direction + direction_B_perp * MaxStepLength2 * k / 2).getNorm() /
absVelocity);
// check: next intersections is MaxStepLength1
} else if (MaxStepLength2 <= 0_m || MaxStepLength1 < MaxStepLength2) {
CORSIKA_LOG_TRACE(" steplength to obs plane 2: {} ", MaxStepLength1);
return Intersections(
MaxStepLength1 *
(direction + direction.cross(magneticfield) * MaxStepLength2 * k / 2)
.getNorm() /
(direction + direction_B_perp * MaxStepLength1 * k / 2).getNorm() /
absVelocity);
}
CORSIKA_LOG_WARN("Particle wasn't tracked with curved trajectory -> straight");
CORSIKA_LOG_WARN(
"Particle wasn't tracked with curved trajectory -> straight (is this an "
"ERROR?)");
} // end if curved-tracking
return tracking_line::Tracking::intersect(particle, plane, medium);
return tracking_line::Tracking::intersect(particle, plane);
}
template <typename TParticle, typename TBaseNodeType>
inline Intersections Tracking::intersect(const TParticle& particle,
const TBaseNodeType& volumeNode) {
const Sphere* sphere = dynamic_cast<const Sphere*>(&volumeNode.getVolume());
if (sphere) {
return intersect(particle, *sphere, volumeNode.getModelProperties());
}
inline Intersections Tracking::intersect(TParticle const& particle,
TBaseNodeType const& volumeNode) {
Sphere const* sphere = dynamic_cast<Sphere const*>(&volumeNode.getVolume());
if (sphere) { return intersect(particle, *sphere); }
throw std::runtime_error(
"The Volume type provided is not supported in intersect(particle, node)");
}
......
corsika/detail/modules/tracking/TrackingLeapFrogStraight.inl
+ 16
− 16
View file @ 8bf837b8
......@@ -28,7 +28,7 @@ namespace corsika {
VelocityVector initialVelocity =
particle.getMomentum() / particle.getEnergy() * constants::c;
const Point initialPosition = particle.getPosition();
Point const initialPosition = particle.getPosition();
CORSIKA_LOG_DEBUG(
"TrackingB pid: {}"
" , E = {} GeV",
......@@ -53,17 +53,17 @@ namespace corsika {
const int chargeNumber = particle.getChargeNumber();
auto magneticfield =
volumeNode->getModelProperties().getMagneticField(initialPosition);
const auto magnitudeB = magneticfield.getNorm();
auto const magnitudeB = magneticfield.getNorm();
CORSIKA_LOG_DEBUG("field={} uT, chargeNumber={}, magnitudeB={} uT",
magneticfield.getComponents() / 1_uT, chargeNumber,
magnitudeB / 1_T);
bool const no_deflection = chargeNumber == 0 || magnitudeB == 0_T;
// check, where the first halve-step direction has geometric intersections
const auto [initialTrack, initialTrackNextVolume] =
auto const [initialTrack, initialTrackNextVolume] =
tracking_line::Tracking::getTrack(particle);
{ [[maybe_unused]] auto& initialTrackNextVolume_dum = initialTrackNextVolume; }
const auto initialTrackLength = initialTrack.getLength(1);
auto const initialTrackLength = initialTrack.getLength(1);
CORSIKA_LOG_DEBUG("initialTrack(0)={}, initialTrack(1)={}, initialTrackLength={}",
initialTrack.getPosition(0).getCoordinates(),
......@@ -94,8 +94,8 @@ namespace corsika {
// need to follow strongly curved trajectories segment-wise,
// at least if we don't employ concepts as "Helix
// Trajectories" or similar
const double maxRadians = 0.01;
const LengthType steplimit = 2 * cos(maxRadians) * sin(maxRadians) * gyroradius;
double const maxRadians = 0.01;
LengthType const steplimit = 2 * cos(maxRadians) * sin(maxRadians) * gyroradius;
CORSIKA_LOG_DEBUG("gyroradius {}, Steplimit: {}", gyroradius, steplimit);
// calculate first halve step for "steplimit"
......@@ -110,12 +110,12 @@ namespace corsika {
CORSIKA_LOG_DEBUG("first halve step length {}, steplimit={}, initialTrackLength={}",
firstHalveSteplength, steplimit, initialTrackLength);
// perform the first halve-step
const Point position_mid = initialPosition + direction * firstHalveSteplength;
const auto k =
chargeNumber * constants::c * 1_eV / (particle.getMomentum().getNorm() * 1_V);
const auto new_direction =
Point const position_mid = initialPosition + direction * firstHalveSteplength;
auto const k =
chargeNumber * (constants::c * 1_eV / 1_V) / particle.getMomentum().getNorm();
auto const new_direction =
direction + direction.cross(magneticfield) * firstHalveSteplength * 2 * k;
const auto new_direction_norm = new_direction.getNorm(); // by design this is >1
auto const new_direction_norm = new_direction.getNorm(); // by design this is >1
CORSIKA_LOG_DEBUG(
"position_mid={}, new_direction={}, (new_direction_norm)={}, deflection={}",
position_mid.getCoordinates(), new_direction.getComponents(),
......@@ -126,7 +126,7 @@ namespace corsika {
// check, where the second halve-step direction has geometric intersections
particle.setPosition(position_mid);
particle.setMomentum(new_direction * absMomentum);
const auto [finalTrack, finalTrackNextVolume] =
auto const [finalTrack, finalTrackNextVolume] =
tracking_line::Tracking::getTrack(particle);
particle.setPosition(initialPosition); // this is not nice...
particle.setMomentum(initialMomentum); // this is not nice...
......@@ -160,12 +160,12 @@ namespace corsika {
// perform the second halve-step
auto const new_direction_normalized = new_direction.normalized();
const Point finalPosition =
Point const finalPosition =
position_mid + new_direction_normalized * secondHalveStepLength;
const LengthType totalStep = firstHalveSteplength + secondHalveStepLength;
const auto delta_pos = finalPosition - initialPosition;
const auto distance = delta_pos.getNorm();
LengthType const totalStep = firstHalveSteplength + secondHalveStepLength;
auto const delta_pos = finalPosition - initialPosition;
auto const distance = delta_pos.getNorm();
return std::make_tuple(
StraightTrajectory(
......
corsika/detail/modules/tracking/TrackingStraight.inl
+ 6
− 12
View file @ 8bf837b8
......@@ -48,9 +48,9 @@ namespace corsika::tracking_line {
minNode); // next volume node
}
template <typename TParticle, typename TMedium>
template <typename TParticle>
inline Intersections Tracking::intersect(TParticle const& particle,
Sphere const& sphere, TMedium const&) {
Sphere const& sphere) {
auto const delta = particle.getPosition() - sphere.getCenter();
auto const velocity = particle.getMomentum() / particle.getEnergy() * constants::c;
auto const vSqNorm = velocity.getSquaredNorm();
......@@ -73,20 +73,14 @@ namespace corsika::tracking_line {
inline Intersections Tracking::intersect(TParticle const& particle,
TBaseNodeType const& volumeNode) {
Sphere const* sphere = dynamic_cast<Sphere const*>(&volumeNode.getVolume());
if (sphere) {
typedef typename std::remove_const_t<
std::remove_reference_t<decltype(volumeNode.getModelProperties())>>
medium_type;
return Tracking::intersect<TParticle, medium_type>(particle, *sphere,
volumeNode.getModelProperties());
}
if (sphere) { return Tracking::intersect<TParticle>(particle, *sphere); }
throw std::runtime_error(
"The Volume type provided is not supported in Intersect(particle, node)");
}
template <typename TParticle, typename TMedium>
inline Intersections Tracking::intersect(TParticle const& particle, Plane const& plane,
TMedium const&) {
template <typename TParticle>
inline Intersections Tracking::intersect(TParticle const& particle,
Plane const& plane) {
auto const delta = plane.getCenter() - particle.getPosition();
auto const velocity = particle.getMomentum() / particle.getEnergy() * constants::c;
auto const n = plane.getNormal();
......
corsika/modules/tracking/TrackingLeapFrogCurved.hpp
+ 8
− 7
View file @ 8bf837b8
......@@ -60,17 +60,18 @@ namespace corsika {
template <typename TParticle>
auto getTrack(TParticle const& particle);
template <typename TParticle, typename TMedium>
static Intersections intersect(TParticle const& particle, Sphere const& sphere,
TMedium const& medium);
//! find intersection of Sphere with Track
template <typename TParticle>
static Intersections intersect(TParticle const& particle, Sphere const& sphere);
//! find intersection of Volume node with Track of particle
template <typename TParticle, typename TBaseNodeType>
static Intersections intersect(TParticle const& particle,
TBaseNodeType const& volumeNode);
TBaseNodeType const& node);
template <typename TParticle, typename TMedium>
static Intersections intersect(TParticle const& particle, Plane const& plane,
TMedium const& medium);
//! find intersection of Plane with Track
template <typename TParticle>
static Intersections intersect(TParticle const& particle, Plane const& plane);
protected:
/**
......
corsika/modules/tracking/TrackingStraight.hpp
+ 6
− 9
View file @ 8bf837b8
......@@ -41,19 +41,16 @@ namespace corsika::tracking_line {
auto getTrack(TParticle const& particle);
//! find intersection of Sphere with Track
template <typename TParticle, typename TMedium>
static Intersections intersect(TParticle const& particle, Sphere const& sphere,
TMedium const&);
template <typename TParticle>
static Intersections intersect(TParticle const& particle, Sphere const& sphere);
//! find intersection of Volume with Track
//! find intersection of Volume node with Track of particle
template <typename TParticle, typename TBaseNodeType>
static Intersections intersect(TParticle const& particle,
TBaseNodeType const& volumeNode);
static Intersections intersect(TParticle const& particle, TBaseNodeType const& node);
//! find intersection of Plane with Track
template <typename TParticle, typename TMedium>
static Intersections intersect(TParticle const& particle, Plane const& plane,
TMedium const&);
template <typename TParticle>
static Intersections intersect(TParticle const& particle, Plane const& plane);
};
} // namespace corsika::tracking_line
......
do-clang-format.py
+ 1
− 0
View file @ 8bf837b8
......@@ -83,6 +83,7 @@ cmd += " -style=file"
version = subp.check_output(cmd.split() + ["--version"]).decode("utf-8")
print (version)
print ("Note: the clang-format version has an impact on the result. Make sure you are consistent with current CI. Consider \'--docker\' option.")
if args.apply:
for filename in filelist:
......
src/modules/qgsjetII/code_generator.py
+ 2
− 2
View file @ 8bf837b8
......@@ -93,7 +93,7 @@ def read_qgsjetII_codes(filename, particle_db):
if len(line)==0 or line[0] == '#':
continue
line = line.split('#')[0]
print (line)
print ('QGSJetII codes: ', line)
identifier, model_code, xsType = line.split()
try:
particle_db[identifier]["qgsjetII_code"] = int(model_code)
......@@ -202,7 +202,7 @@ if __name__ == "__main__":
particle_db = load_particledb(sys.argv[1])
read_qgsjetII_codes(sys.argv[2], particle_db)
set_default_qgsjetII_definition(particle_db)
with open("Generated.inc", "w") as f:
print("// this file is automatically generated\n// edit at your own risk!\n", file=f)
print(generate_qgsjetII_start(), file=f)
......
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