IAP GITLAB

Commit 9a7a8ab5 authored by Ralf Ulrich's avatar Ralf Ulrich

more debug output better names

parent 2c0c4c24
Pipeline #3526 failed with stages
in 13 minutes and 7 seconds
...@@ -59,9 +59,10 @@ namespace corsika { ...@@ -59,9 +59,10 @@ namespace corsika {
for (auto const& node : volumeNode.getChildNodes()) { for (auto const& node : volumeNode.getChildNodes()) {
Intersections const time_intersections = TDerived::intersect(particle, *node); Intersections const time_intersections = TDerived::intersect(particle, *node);
CORSIKA_LOG_TRACE("intersection times with child volume {}", fmt::ptr(node));
if (!time_intersections.hasIntersections()) { continue; } if (!time_intersections.hasIntersections()) { continue; }
CORSIKA_LOG_DEBUG("intersection times with child volume {} : enter {} s, exit {} s", CORSIKA_LOG_TRACE(" : enter {} s, exit {} s",
fmt::ptr(node), time_intersections.getEntry() / 1_s, time_intersections.getEntry() / 1_s,
time_intersections.getExit() / 1_s); time_intersections.getExit() / 1_s);
auto const t_entry = time_intersections.getEntry(); auto const t_entry = time_intersections.getEntry();
......
...@@ -149,21 +149,25 @@ namespace corsika { ...@@ -149,21 +149,25 @@ namespace corsika {
} }
bool const numericallyInside = sphere.contains(particle.getPosition()); bool const numericallyInside = sphere.contains(particle.getPosition());
CORSIKA_LOG_TRACE("numericallyInside={}", numericallyInside);
auto const absVelocity = velocity.getNorm(); auto const absVelocity = velocity.getNorm();
auto energy = particle.getEnergy(); auto const energy = particle.getEnergy();
auto k = chargeNumber * constants::cSquared * 1_eV / (absVelocity * energy * 1_V); // this is: k = q/|p|
auto const k =
chargeNumber * constants::cSquared * 1_eV / (absVelocity * energy * 1_V);
auto const direction_B_perp = directionBefore.cross(magneticfield); auto const direction_B_perp = directionBefore.cross(magneticfield);
auto const denom = direction_B_perp.getSquaredNorm() * k * k; auto const denom = 4. / (direction_B_perp.getSquaredNorm() * k * k);
Vector<length_d> const deltaPos = position - sphere.getCenter(); 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 = (direction_B_perp.dot(deltaPos) * k + 1) * denom / (1_m * 1_m);
double const b = directionBefore.dot(deltaPos) * 8 / (denom * 1_m * 1_m * 1_m); double const c = directionBefore.dot(deltaPos) * 2 * denom / (1_m * 1_m * 1_m);
double const c = LengthType const deltaPosLength = deltaPos.getNorm();
(deltaPos.getSquaredNorm() - (sphere.getRadius() * sphere.getRadius())) * 4 / double const d = (deltaPosLength + sphere.getRadius()) *
(denom * 1_m * 1_m * 1_m * 1_m); (deltaPosLength - sphere.getRadius()) * denom /
CORSIKA_LOG_TRACE("denom={}, a={}, b={}, c={}", denom, a, b, c); (1_m * 1_m * 1_m * 1_m);
std::complex<double>* solutions = quartic_solver::solve_quartic(0, a, b, c); CORSIKA_LOG_TRACE("denom={}, b={}, c={}, d={}", denom, b, c, d);
std::complex<double> const* solutions = quartic_solver::solve_quartic(0, b, c, d);
LengthType d_enter, d_exit; LengthType d_enter, d_exit;
int first = 0, first_entry = 0, first_exit = 0; int first = 0, first_entry = 0, first_exit = 0;
for (int i = 0; i < 4; i++) { for (int i = 0; i < 4; i++) {
...@@ -172,8 +176,8 @@ namespace corsika { ...@@ -172,8 +176,8 @@ namespace corsika {
CORSIKA_LOG_TRACE("Solution (real) for current Volume: {} ", dist); CORSIKA_LOG_TRACE("Solution (real) for current Volume: {} ", dist);
if (numericallyInside) { if (numericallyInside) {
// there must be an entry (negative) and exit (positive) solution // there must be an entry (negative) and exit (positive) solution
if (dist < -0.0001_m) { // security margin to assure transfer to next if (dist < 0.0001_m) { // security margin to assure transfer to next
// logical volume // logical volume
if (first_entry == 0) { if (first_entry == 0) {
d_enter = dist; d_enter = dist;
} else { } else {
...@@ -181,7 +185,7 @@ namespace corsika { ...@@ -181,7 +185,7 @@ namespace corsika {
} }
first_entry++; first_entry++;
} else { // thus, dist >= -0.0001_m } else { // thus, dist > -0.0001_m
if (first_exit == 0) { if (first_exit == 0) {
d_exit = dist; d_exit = dist;
...@@ -196,7 +200,7 @@ namespace corsika { ...@@ -196,7 +200,7 @@ namespace corsika {
// both physical solutions (entry, exit) must be positive, and as small as // both physical solutions (entry, exit) must be positive, and as small as
// possible // possible
if (dist < -0.0001_m) { // need small numerical margin, to assure transport if (dist < 0.0001_m) { // need small numerical margin, to assure transport
// into next logical volume // into next logical volume
continue; continue;
} }
...@@ -268,6 +272,9 @@ namespace corsika { ...@@ -268,6 +272,9 @@ namespace corsika {
LengthType const MaxStepLength1 = (sqrSqrtArg - norm_projected) / denom; LengthType const MaxStepLength1 = (sqrSqrtArg - norm_projected) / denom;
LengthType const MaxStepLength2 = (-sqrSqrtArg - norm_projected) / denom; LengthType const MaxStepLength2 = (-sqrSqrtArg - norm_projected) / denom;
CORSIKA_LOG_TRACE("MaxStepLength1={}, MaxStepLength2={}", MaxStepLength1,
MaxStepLength2);
// check: both intersections in past // check: both intersections in past
if (MaxStepLength1 <= 0_m && MaxStepLength2 <= 0_m) { if (MaxStepLength1 <= 0_m && MaxStepLength2 <= 0_m) {
return Intersections(std::numeric_limits<double>::infinity() * 1_s); return Intersections(std::numeric_limits<double>::infinity() * 1_s);
...@@ -295,6 +302,9 @@ namespace corsika { ...@@ -295,6 +302,9 @@ namespace corsika {
} // end if curved-tracking } // end if curved-tracking
CORSIKA_LOG_TRACE("straight tracking with chargeNumber={}, B={}", chargeNumber,
magneticfield);
return tracking_line::Tracking::intersect(particle, plane); return tracking_line::Tracking::intersect(particle, plane);
} }
......
...@@ -84,13 +84,14 @@ namespace corsika::tracking_line { ...@@ -84,13 +84,14 @@ namespace corsika::tracking_line {
auto const delta = plane.getCenter() - particle.getPosition(); auto const delta = plane.getCenter() - particle.getPosition();
auto const velocity = particle.getMomentum() / particle.getEnergy() * constants::c; auto const velocity = particle.getMomentum() / particle.getEnergy() * constants::c;
auto const n = plane.getNormal(); auto const n = plane.getNormal();
auto const c = n.dot(velocity); auto const n_dot_v = n.dot(velocity);
CORSIKA_LOG_TRACE("c={}, delta={}, momentum={}", c, delta, particle.getMomentum()); CORSIKA_LOG_TRACE("n_dot_v={}, delta={}, momentum={}", n_dot_v, delta,
particle.getMomentum());
return Intersections(c.magnitude() == 0 return Intersections(n_dot_v.magnitude() == 0
? std::numeric_limits<TimeType::value_type>::infinity() * 1_s ? std::numeric_limits<TimeType::value_type>::infinity() * 1_s
: n.dot(delta) / c); : n.dot(delta) / n_dot_v);
} }
} // namespace corsika::tracking_line } // namespace corsika::tracking_line
...@@ -6,6 +6,8 @@ ...@@ -6,6 +6,8 @@
* the license. * the license.
*/ */
#define DEBUG 1
/* clang-format off */ /* clang-format off */
// InteractionCounter used boost/histogram, which // InteractionCounter used boost/histogram, which
// fails if boost/type_traits have been included before. Thus, we have // fails if boost/type_traits have been included before. Thus, we have
...@@ -129,9 +131,25 @@ int main(int argc, char** argv) { ...@@ -129,9 +131,25 @@ int main(int argc, char** argv) {
builder.addExponentialLayer(1144.9069_g / (1_cm * 1_cm), 878153.55_cm, 10_km); builder.addExponentialLayer(1144.9069_g / (1_cm * 1_cm), 878153.55_cm, 10_km);
builder.addExponentialLayer(1305.5948_g / (1_cm * 1_cm), 636143.04_cm, 40_km); builder.addExponentialLayer(1305.5948_g / (1_cm * 1_cm), 636143.04_cm, 40_km);
builder.addExponentialLayer(540.1778_g / (1_cm * 1_cm), 772170.16_cm, 100_km); builder.addExponentialLayer(540.1778_g / (1_cm * 1_cm), 772170.16_cm, 100_km);
builder.addLinearLayer(1e9_cm, 112.8_km); builder.addLinearLayer(1e9_cm, 112.8_km+constants::EarthRadius::Mean);
builder.assemble(env); builder.assemble(env);
CORSIKA_LOG_DEBUG(
"environment setup: universe={}, layer1={}, layer2={}, layer3={}, layer4={}, "
"layer5={}",
fmt::ptr(env.getUniverse()->getContainingNode(
Point(rootCS, {constants::EarthRadius::Mean + 130_km, 0_m, 0_m}))),
fmt::ptr(env.getUniverse()->getContainingNode(
Point(rootCS, {constants::EarthRadius::Mean + 110_km, 0_m, 0_m}))),
fmt::ptr(env.getUniverse()->getContainingNode(
Point(rootCS, {constants::EarthRadius::Mean + 50_km, 0_m, 0_m}))),
fmt::ptr(env.getUniverse()->getContainingNode(
Point(rootCS, {constants::EarthRadius::Mean + 20_km, 0_m, 0_m}))),
fmt::ptr(env.getUniverse()->getContainingNode(
Point(rootCS, {constants::EarthRadius::Mean + 5_km, 0_m, 0_m}))),
fmt::ptr(env.getUniverse()->getContainingNode(
Point(rootCS, {constants::EarthRadius::Mean + 2_km, 0_m, 0_m}))));
// setup particle stack, and add primary particle // setup particle stack, and add primary particle
setup::Stack stack; setup::Stack stack;
stack.clear(); stack.clear();
...@@ -159,7 +177,7 @@ int main(int argc, char** argv) { ...@@ -159,7 +177,7 @@ int main(int argc, char** argv) {
<< ", norm = " << plab.getNorm() << endl; << ", norm = " << plab.getNorm() << endl;
auto const observationHeight = 0_km + builder.getEarthRadius(); auto const observationHeight = 0_km + builder.getEarthRadius();
auto const injectionHeight = 112.75_km + builder.getEarthRadius(); auto const injectionHeight = 111.75_km + builder.getEarthRadius();
auto const t = -observationHeight * cos(thetaRad) + auto const t = -observationHeight * cos(thetaRad) +
sqrt(-static_pow<2>(sin(thetaRad) * observationHeight) + sqrt(-static_pow<2>(sin(thetaRad) * observationHeight) +
static_pow<2>(injectionHeight)); static_pow<2>(injectionHeight));
......
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