IAP GITLAB

Skip to content
Snippets Groups Projects
Commit 8463b28a authored by ralfulrich's avatar ralfulrich Committed by Ralf Ulrich
Browse files

rebase to master

parent 2e83a9cf
No related branches found
No related tags found
1 merge request!271Resolve "Assumption of no continuous between two stochastic losses"
Pipeline #4144 failed
......@@ -33,7 +33,6 @@ namespace corsika {
long double a2 = a * a;
long double q = (a2 - 3 * b) / 9;
long double r = (a * (2 * a2 - 9 * b) + 27 * c) / 54;
long double r2 = r * r;
long double q3 = q * q * q;
// disc = q**3 + r**2
......@@ -56,8 +55,6 @@ namespace corsika {
// sum all terms into final result
long double const disc = -(((e + uf) + au) + ab) + v;
CORSIKA_LOG_TRACE("solveP3 disc={} r2-r3={}", disc, r2 - q3);
if (disc >= 0) {
long double t = r / std::sqrt(q3);
if (t < -1) t = -1;
......
......@@ -61,9 +61,9 @@ namespace corsika {
}
} else {
if (Det < 0) return {};
long double sqDet = sqrt(Det);
q1 = (y + sqDet) * 0.5;
q2 = (y - sqDet) * 0.5;
long double sqDet1 = sqrt(Det);
q1 = (y + sqDet1) * 0.5;
q2 = (y - sqDet1) * 0.5;
// g1+g2 = b && g1*h2 + g2*h1 = c ( && g === p ) Krammer
p1 = (b * q1 - d) / (q1 - q2);
p2 = (d - b * q2) / (q1 - q2);
......@@ -71,32 +71,30 @@ namespace corsika {
// solving quadratic eqs. x^2 + p1*x + q1 = 0
// x^2 + p2*x + q2 = 0
#ifdef DEBUG
{
long double Det1 = p1 * p1 - 4 * q1;
long double Det2 = p2 * p2 - 4 * q2;
[[maybe_unused]] long double Det1 = p1 * p1 - 4 * q1;
[[maybe_unused]] long double Det2 = p2 * p2 - 4 * q2;
CORSIKA_LOG_TRACE("Det1={} Det2={}", Det1, Det2);
if (Det1 >= 0 && Det2 >= 0) {
long double sqDet1 = sqrt(Det1);
long double sqDet2 = sqrt(Det2);
// return {double(-p1 + sqDet1) * 0.5, double(-p1 - sqDet1) * 0.5, double(-p2 +
// sqDet2) * 0.5, double(-p2 - sqDet2) * 0.5};
[[maybe_unused]] long double sqDet1 = sqrt(Det1);
[[maybe_unused]] long double sqDet2 = sqrt(Det2);
CORSIKA_LOG_TRACE("check1 {} {} {} {}", double(-p1 + sqDet1) * 0.5,
double(-p1 - sqDet1) * 0.5, double(-p2 + sqDet2) * 0.5,
double(-p2 - sqDet2) * 0.5);
}
if (Det1 >= 0) {
long double sqDet1 = sqrt(Det1);
// return {double(-p1 + sqDet1) * 0.5, double(-p1 - sqDet1) * 0.5};
[[maybe_unused]] long double sqDet1 = sqrt(Det1);
CORSIKA_LOG_TRACE("check2 {} {} ", double(-p1 + sqDet1) * 0.5,
double(-p1 - sqDet1) * 0.5);
}
if (Det2 >= 0) {
long double sqDet2 = sqrt(Det2);
// return {double(-p2 + sqDet2) * 0.5, double(-p2 - sqDet2) * 0.5};
[[maybe_unused]] long double sqDet2 = sqrt(Det2);
CORSIKA_LOG_TRACE("check3 {} {}", double(-p1 + sqDet2) * 0.5,
double(-p1 - sqDet2) * 0.5);
}
}
#endif
std::vector<double> quad1 = solve_quadratic_real(1, p1, q1);
std::vector<double> quad2 = solve_quadratic_real(1, p2, q2);
......
......@@ -65,8 +65,7 @@ namespace corsika::tracking_line {
auto const sqDisc = sqrt(discriminant);
auto const invDenom = 1 / vSqNorm;
bool const numericallyInside = sphere.contains(position);
CORSIKA_LOG_TRACE("numericallyInside={}", numericallyInside);
CORSIKA_LOG_TRACE("numericallyInside={}", sphere.contains(position));
return Intersections((-vDotDelta - sqDisc) * invDenom,
(-vDotDelta + sqDisc) * invDenom);
}
......
......@@ -103,7 +103,7 @@ int main(int argc, char** argv) {
CORSIKA_LOG_INFO("vertical_EAS");
if (argc < 5) {
std::cerr << "usage: vertical_EAS <A> <Z> <energy/GeV> <Nevt> [seed] [output-dir] \n"
std::cerr << "usage: vertical_EAS <A> <Z> <energy/GeV> <Nevt> [seed] \n"
" if A=0, Z is interpreted as PDG code \n"
" if no seed is given, a random seed is chosen \n"
<< std::endl;
......@@ -116,8 +116,6 @@ int main(int argc, char** argv) {
if (argc > 5) { seed = std::stoi(std::string(argv[5])); }
CORSIKA_LOG_INFO("output_dir={} seed={}", output_dir, seed);
// initialize random number sequence(s)
registerRandomStreams(seed);
......@@ -400,7 +398,7 @@ int main(int argc, char** argv) {
BetheBlochPDG emContinuous(showerAxis);
OnShellCheck reset_particle_mass(1.e-3, 1.e-1, false);
TrackWriter trackWriter(tracks_dir);
TrackWriter trackWriter(tracks_file);
LongitudinalProfile longprof{showerAxis};
......@@ -433,8 +431,6 @@ int main(int argc, char** argv) {
cut.reset();
// emContinuous.reset();
longprof.save(longprof_dat);
auto const hists = sibyllCounted.getHistogram() + sibyllNucCounted.getHistogram() +
urqmdCounted.getHistogram();
......
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment