diff --git a/Documentation/Examples/vertical_EAS.cc b/Documentation/Examples/vertical_EAS.cc
index d5cb7baa024ca7eecbd0e8b88c97b27fe1fd96d9..235c5b3f5156224acef15a301b3cd84a1a306dca 100644
--- a/Documentation/Examples/vertical_EAS.cc
+++ b/Documentation/Examples/vertical_EAS.cc
@@ -227,7 +227,6 @@ int main(int argc, char** argv) {
decaySibyll.PrintDecayConfig();
-
process::on_shell_check::OnShellCheck reset_particle_mass(1.e-3, 1.e-1, false);
process::track_writer::TrackWriter trackWriter("tracks.dat");
diff --git a/Environment/BaseExponential.h b/Environment/BaseExponential.h
index c5ac3181f49d06d9e2c74231e3ce40580ecd265e..94b454bf686f4b0657b2e5692df8fb0e3b287816 100644
--- a/Environment/BaseExponential.h
+++ b/Environment/BaseExponential.h
@@ -80,8 +80,7 @@ namespace corsika::environment {
*/
// clang-format on
units::si::LengthType ArclengthFromGrammage(
- geometry::LineTrajectory const& vLine,
- units::si::GrammageType vGrammage,
+ geometry::LineTrajectory const& vLine, units::si::GrammageType vGrammage,
geometry::Vector<units::si::dimensionless_d> const& vAxis) const {
auto const uDotA = vLine.GetDirection(0).dot(vAxis).magnitude();
auto const rhoStart = GetImplementation().GetMassDensity(vLine.GetLine().GetR0());
diff --git a/Environment/FlatExponential.h b/Environment/FlatExponential.h
index c29dd91c7c0b3bc39e5e0f26d28e89f4e8bf1d55..31f3309943f3cc2c927b80dc4f29a32d8c61c82d 100644
--- a/Environment/FlatExponential.h
+++ b/Environment/FlatExponential.h
@@ -51,9 +51,8 @@ namespace corsika::environment {
NuclearComposition const& GetNuclearComposition() const override { return fNuclComp; }
- units::si::GrammageType IntegratedGrammage(
- geometry::LineTrajectory const& vLine,
- units::si::LengthType vTo) const override {
+ units::si::GrammageType IntegratedGrammage(geometry::LineTrajectory const& vLine,
+ units::si::LengthType vTo) const override {
return Base::IntegratedGrammage(vLine, vTo, fAxis);
}
diff --git a/Environment/LinearApproximationIntegrator.h b/Environment/LinearApproximationIntegrator.h
index 66bef41e2691b21152e230674c239315901e37c1..230d756fc9fcb7893f51b72605625945c8da8b6e 100644
--- a/Environment/LinearApproximationIntegrator.h
+++ b/Environment/LinearApproximationIntegrator.h
@@ -20,21 +20,19 @@ namespace corsika::environment {
auto const& GetImplementation() const { return *static_cast<TDerived const*>(this); }
public:
- auto IntegrateGrammage(
- corsika::geometry::LineTrajectory const& line,
- corsika::units::si::LengthType length) const {
+ auto IntegrateGrammage(corsika::geometry::LineTrajectory const& line,
+ corsika::units::si::LengthType length) const {
auto const c0 = GetImplementation().EvaluateAt(line.GetPosition(0));
- auto const c1 = GetImplementation().fRho.FirstDerivative(
- line.GetPosition(0), line.GetDirection(0));
+ auto const c1 = GetImplementation().fRho.FirstDerivative(line.GetPosition(0),
+ line.GetDirection(0));
return (c0 + 0.5 * c1 * length) * length;
}
- auto ArclengthFromGrammage(
- corsika::geometry::LineTrajectory const& line,
- corsika::units::si::GrammageType grammage) const {
+ auto ArclengthFromGrammage(corsika::geometry::LineTrajectory const& line,
+ corsika::units::si::GrammageType grammage) const {
auto const c0 = GetImplementation().fRho(line.GetPosition(0));
- auto const c1 = GetImplementation().fRho.FirstDerivative(
- line.GetPosition(0), line.GetDirection(0));
+ auto const c1 = GetImplementation().fRho.FirstDerivative(line.GetPosition(0),
+ line.GetDirection(0));
return (1 - 0.5 * grammage * c1 / (c0 * c0)) * grammage / c0;
}
diff --git a/Environment/NoMagneticField.h b/Environment/NoMagneticField.h
index 38e99347e7f474244483003d4dd3dc58ea7d1143..74cbd8cfc51edac6db9bd05269f33ddb37aa23b9 100644
--- a/Environment/NoMagneticField.h
+++ b/Environment/NoMagneticField.h
@@ -38,8 +38,7 @@ namespace corsika::environment {
*/
template <typename... Args>
NoMagneticField(Args&&... args)
- : T(std::forward<Args>(args)...)
- {}
+ : T(std::forward<Args>(args)...) {}
/**
* Evaluate the magnetic field at a given location.
@@ -50,7 +49,7 @@ namespace corsika::environment {
MagneticFieldVector GetMagneticField(
corsika::geometry::Point const&) const final override {
CoordinateSystem const& gCS =
- RootCoordinateSystem::GetInstance().GetRootCoordinateSystem();
+ RootCoordinateSystem::GetInstance().GetRootCoordinateSystem();
return MagneticFieldVector(gCS, {0_T, 0_T, 0_T});
}
diff --git a/Environment/VolumeTreeNode.h b/Environment/VolumeTreeNode.h
index 1cd5cf222ff1beb60c8ef896bc0dac2707087632..8cad83555833331e3d3c57912379f33ca3ac8eee 100644
--- a/Environment/VolumeTreeNode.h
+++ b/Environment/VolumeTreeNode.h
@@ -92,7 +92,7 @@ namespace corsika::environment {
void ExcludeOverlapWith(VTNUPtr const& pNode) {
fExcludedNodes.push_back(pNode.get());
}
-
+
const VTN_type* GetParent() const { return fParentNode; };
auto const& GetChildNodes() const { return fChildNodes; }
diff --git a/Framework/Cascade/Cascade.h b/Framework/Cascade/Cascade.h
index 9ae6724df21134a98c69a96834c69ea39f31ebef..3e81ef973e333f1f3d264810f0f1a88138e5acfa 100644
--- a/Framework/Cascade/Cascade.h
+++ b/Framework/Cascade/Cascade.h
@@ -38,7 +38,6 @@ using boost::typeindex::type_id_with_cvr;
#include <fstream>
-
/**
* The cascade namespace assembles all objects needed to simulate full particles cascades.
*/
@@ -98,10 +97,8 @@ namespace corsika::cascade {
C8LOG_INFO(" - With full cascade HISTORY.");
}
}
-
- ~Cascade(){
- };
+ ~Cascade(){};
/**
* The Run function is the main simulation loop, which processes
@@ -195,36 +192,39 @@ namespace corsika::cascade {
// convert next_decay from time to length [m]
LengthType const distance_decay = next_decay * vParticle.GetMomentum().norm() /
vParticle.GetEnergy() * units::constants::c;
-
+
// determine geometric tracking
auto [step, nextVol] = tracking_.GetTrack(vParticle);
auto geomMaxLength = step.GetLength(1);
-
+
// convert next_step from grammage to length
LengthType const distance_interact =
currentLogicalNode->GetModelProperties().ArclengthFromGrammage(step,
next_interact);
-
+
// determine the maximum geometric step length
LengthType const distance_max = process_sequence_.MaxStepLength(vParticle, step);
C8LOG_DEBUG("distance_max={} m", distance_max / 1_m);
// take minimum of geometry, interaction, decay for next step
- auto min_distance = std::min(
- {distance_interact, distance_decay, distance_max, geomMaxLength});
+ auto min_distance =
+ std::min({distance_interact, distance_decay, distance_max, geomMaxLength});
- C8LOG_DEBUG("transport particle by : {} m "
- "Medium transition after: {} m "
- "Decay after: {} m "
- "Interaction after: {} m",
- min_distance/1_m, geomMaxLength/1_m, distance_decay/1_m, distance_interact/1_m);
+ C8LOG_DEBUG(
+ "transport particle by : {} m "
+ "Medium transition after: {} m "
+ "Decay after: {} m "
+ "Interaction after: {} m",
+ min_distance / 1_m, geomMaxLength / 1_m, distance_decay / 1_m,
+ distance_interact / 1_m);
// here the particle is actually moved along the trajectory to new position:
step.SetLength(min_distance);
vParticle.SetPosition(step.GetPosition(1));
- vParticle.SetMomentum(step.GetDirection(1)*vParticle.GetMomentum().norm());
+ vParticle.SetMomentum(step.GetDirection(1) * vParticle.GetMomentum().norm());
vParticle.SetTime(vParticle.GetTime() + step.GetDuration());
- std::cout << "New Position: " << vParticle.GetPosition().GetCoordinates() << std::endl;
+ std::cout << "New Position: " << vParticle.GetPosition().GetCoordinates()
+ << std::endl;
// apply all continuous processes on particle + track
process::EProcessReturn status = process_sequence_.DoContinuous(vParticle, step);
@@ -380,7 +380,6 @@ Y8, Y8, ,8P 88 `8b `8b 88 88P Y8b d8"
corsika::random::RNGManager::GetInstance().GetRandomStream("cascade");
unsigned int count_ = 0;
-
}; // end class Cascade
} // namespace corsika::cascade
diff --git a/Framework/Cascade/Cascade_interpolation.h b/Framework/Cascade/Cascade_interpolation.h
index dccd3f8e705ca264fa27e7b637bdc4ab332b74b8..41808e5c5575e6816a66a85a5ce8f9e228c7139f 100644
--- a/Framework/Cascade/Cascade_interpolation.h
+++ b/Framework/Cascade/Cascade_interpolation.h
@@ -194,16 +194,17 @@ namespace corsika::cascade {
// convert next_decay from time to length [m]
LengthType const distance_decay = next_decay * vParticle.GetMomentum().norm() /
vParticle.GetEnergy() * units::constants::c;
-
+
// determine geometric tracking
- auto [step, geomMaxLength, nextVol, magMaxLength, directionBefore, directionAfter] = fTracking.GetTrack(vParticle);
+ auto [step, geomMaxLength, nextVol, magMaxLength, directionBefore, directionAfter] =
+ fTracking.GetTrack(vParticle);
[[maybe_unused]] auto const& dummy_nextVol = nextVol;
-
+
// convert next_step from grammage to length
LengthType const distance_interact =
currentLogicalNode->GetModelProperties().ArclengthFromGrammage(step,
next_interact);
-
+
// determine the maximum geometric step length
LengthType const distance_max = fProcessSequence.MaxStepLength(vParticle, step);
std::cout << "distance_max=" << distance_max << std::endl;
@@ -217,9 +218,10 @@ namespace corsika::cascade {
// here the particle is actually moved along the trajectory to new position:
// std::visit(setup::ParticleUpdate<Particle>{vParticle}, step);
vParticle.SetPosition(step.PositionFromArclength(min_distance));
- // .... also update time, momentum, direction, ...
- vParticle.SetMomentum((directionBefore * (1 - min_distance / magMaxLength) +
- directionAfter * min_distance /magMaxLength) * vParticle.GetMomentum().GetNorm());
+ // .... also update time, momentum, direction, ...
+ vParticle.SetMomentum((directionBefore * (1 - min_distance / magMaxLength) +
+ directionAfter * min_distance / magMaxLength) *
+ vParticle.GetMomentum().GetNorm());
vParticle.SetTime(vParticle.GetTime() + min_distance / units::constants::c);
step.LimitEndTo(min_distance);
diff --git a/Framework/Cascade/testCascade.h b/Framework/Cascade/testCascade.h
index 8c78f740c1eb0338ee05eb36766c736bf1b5cf7b..848c597ee3985e15aa2cfdad956b1dc621e2ca21 100644
--- a/Framework/Cascade/testCascade.h
+++ b/Framework/Cascade/testCascade.h
@@ -16,7 +16,8 @@
#include <corsika/stack/node/GeometryNodeStackExtension.h>
#include <corsika/stack/nuclear_extension/NuclearStackExtension.h>
-using TestEnvironmentInterface = corsika::environment::IMagneticFieldModel<corsika::environment::IMediumModel>;
+using TestEnvironmentInterface =
+ corsika::environment::IMagneticFieldModel<corsika::environment::IMediumModel>;
using TestEnvironmentType = corsika::environment::Environment<TestEnvironmentInterface>;
template <typename T>
diff --git a/Framework/Geometry/Line.h b/Framework/Geometry/Line.h
index 2d27e64df19415d0ab7bc1dfbb2705f58bef76db..089cc6120c4a54fbab5dbe3c507519afbc0bdea8 100644
--- a/Framework/Geometry/Line.h
+++ b/Framework/Geometry/Line.h
@@ -16,7 +16,7 @@ namespace corsika::geometry {
/**
* \class Line
- *
+ *
* A Line describes a movement in three dimensional space. It
* consists of a Point `$\vec{p_0}$` and and a speed-Vector
* `$\vec{v}$`, so that it can return GetPosition as
@@ -30,7 +30,7 @@ namespace corsika::geometry {
Point const r0;
VelocityVec const v0;
-
+
public:
Line() = delete;
Line(const Line&) = default;
diff --git a/Framework/Geometry/Sphere.h b/Framework/Geometry/Sphere.h
index 7d2ab39a360a090d0565d3c1fecbd579ab249b5a..a2a77efec78e2d6c3b063185c73be8423c265c5a 100644
--- a/Framework/Geometry/Sphere.h
+++ b/Framework/Geometry/Sphere.h
@@ -31,7 +31,6 @@ namespace corsika::geometry {
const Point& GetCenter() const { return fCenter; }
LengthType GetRadius() const { return fRadius; }
-
};
} // namespace corsika::geometry
diff --git a/Framework/Geometry/Trajectory.h b/Framework/Geometry/Trajectory.h
index 2c0e17afc340a8887553d21319ccfa5912d893ff..de49a7272d6c4525f7c98c6852aac7d9a668a1c3 100644
--- a/Framework/Geometry/Trajectory.h
+++ b/Framework/Geometry/Trajectory.h
@@ -162,14 +162,15 @@ namespace corsika::geometry {
return initialVelocity_ +
initialVelocity_.cross(magneticfield_) * timeStep_ * u * k_;
}
-
+
Vector<corsika::units::si::dimensionless_d> GetDirection(double u) const {
return GetVelocity(u).normalized();
}
///! duration along potentially bend trajectory
corsika::units::si::TimeType GetDuration(double u = 1) const {
- return u * timeStep_ * (double(GetVelocity(u).norm()/initialVelocity_.norm()) + 1.0) / 2;
+ return u * timeStep_ *
+ (double(GetVelocity(u).norm() / initialVelocity_.norm()) + 1.0) / 2;
}
///! total length along potentially bend trajectory
diff --git a/Framework/Geometry/testGeometry.cc b/Framework/Geometry/testGeometry.cc
index e522d0ab590ffd7d87967fe8cba37566ca54090d..b528bf1bf3d4d922092a0e517ba8dd1208184749 100644
--- a/Framework/Geometry/testGeometry.cc
+++ b/Framework/Geometry/testGeometry.cc
@@ -260,6 +260,5 @@ TEST_CASE("Trajectories") {
(helix.GetPosition(7_s) - helix.PositionFromArclength(helix.ArcLength(0_s, 7_s)))
.norm()
.magnitude() == Approx(0).margin(absMargin));
-
}
}
diff --git a/Framework/Utilities/quartic.cpp b/Framework/Utilities/quartic.cpp
index f11a8e083e887a3f11ac05dba394f232fdcedd3f..88da8cc905fc94307d9e89104a9d47add6747702 100644
--- a/Framework/Utilities/quartic.cpp
+++ b/Framework/Utilities/quartic.cpp
@@ -1,22 +1,10 @@
-/***************************************************************************
- * Copyright (C) 2016 by Саша Миленковић *
- * sasa.milenkovic.xyz@gmail.com *
- * *
- * This program is free software; you can redistribute it and/or modify *
- * it under the terms of the GNU General Public License as published by *
- * the Free Software Foundation; either version 2 of the License, or *
- * (at your option) any later version. *
- * This program is distributed in the hope that it will be useful, *
- * but WITHOUT ANY WARRANTY; without even the implied warranty of *
- * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
- * GNU General Public License for more details. *
- * ( http://www.gnu.org/licenses/gpl-3.0.en.html ) *
- * *
- * You should have received a copy of the GNU General Public License *
- * along with this program; if not, write to the *
- * Free Software Foundation, Inc., *
- * 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. *
- ***************************************************************************/
+/*
+ * (c) Copyright 2020 CORSIKA Project, corsika-project@lists.kit.edu
+ *
+ * This software is distributed under the terms of the GNU General Public
+ * Licence version 3 (GPL Version 3). See file LICENSE for a full version of
+ * the license.
+ */
#include <cmath>
#include "quartic.h"
@@ -27,126 +15,116 @@
// In case 3 real roots: => x[0], x[1], x[2], return 3
// 2 real roots: x[0], x[1], return 2
// 1 real root : x[0], x[1] ± i*x[2], return 1
-unsigned int solveP3(double *x,double a,double b,double c) {
- double a2 = a*a;
- double q = (a2 - 3*b)/9;
- double r = (a*(2*a2-9*b) + 27*c)/54;
- double r2 = r*r;
- double q3 = q*q*q;
- double A,B;
- if(r2<q3)
- {
- double t=r/sqrt(q3);
- if( t<-1) t=-1;
- if( t> 1) t= 1;
- t=acos(t);
- a/=3; q=-2*sqrt(q);
- x[0]=q*cos(t/3)-a;
- x[1]=q*cos((t+M_2PI)/3)-a;
- x[2]=q*cos((t-M_2PI)/3)-a;
- return 3;
- }
- else
- {
- A =-pow(fabs(r)+sqrt(r2-q3),1./3);
- if( r<0 ) A=-A;
- B = (0==A ? 0 : q/A);
-
- a/=3;
- x[0] =(A+B)-a;
- x[1] =-0.5*(A+B)-a;
- x[2] = 0.5*sqrt(3.)*(A-B);
- if(fabs(x[2])<eps) { x[2]=x[1]; return 2; }
-
- return 1;
- }
+unsigned int solveP3(double* x, double a, double b, double c) {
+ double a2 = a * a;
+ double q = (a2 - 3 * b) / 9;
+ double r = (a * (2 * a2 - 9 * b) + 27 * c) / 54;
+ double r2 = r * r;
+ double q3 = q * q * q;
+ double A, B;
+ if (r2 < q3) {
+ double t = r / sqrt(q3);
+ if (t < -1) t = -1;
+ if (t > 1) t = 1;
+ t = acos(t);
+ a /= 3;
+ q = -2 * sqrt(q);
+ x[0] = q * cos(t / 3) - a;
+ x[1] = q * cos((t + M_2PI) / 3) - a;
+ x[2] = q * cos((t - M_2PI) / 3) - a;
+ return 3;
+ } else {
+ A = -pow(fabs(r) + sqrt(r2 - q3), 1. / 3);
+ if (r < 0) A = -A;
+ B = (0 == A ? 0 : q / A);
+
+ a /= 3;
+ x[0] = (A + B) - a;
+ x[1] = -0.5 * (A + B) - a;
+ x[2] = 0.5 * sqrt(3.) * (A - B);
+ if (fabs(x[2]) < eps) {
+ x[2] = x[1];
+ return 2;
+ }
+
+ return 1;
+ }
}
//---------------------------------------------------------------------------
// solve quartic equation x^4 + a*x^3 + b*x^2 + c*x + d
-// Attention - this function returns dynamically allocated array. It has to be released afterwards.
-DComplex* solve_quartic(double a, double b, double c, double d)
-{
- double a3 = -b;
- double b3 = a*c -4.*d;
- double c3 = -a*a*d - c*c + 4.*b*d;
-
- // cubic resolvent
- // y^3 − b*y^2 + (ac−4d)*y − a^2*d−c^2+4*b*d = 0
-
- double x3[3];
- unsigned int iZeroes = solveP3(x3, a3, b3, c3);
-
- double q1, q2, p1, p2, D, sqD, y;
-
- y = x3[0];
- // The essence - choosing Y with maximal absolute value.
- if(iZeroes != 1)
- {
- if(fabs(x3[1]) > fabs(y)) y = x3[1];
- if(fabs(x3[2]) > fabs(y)) y = x3[2];
- }
-
- // h1+h2 = y && h1*h2 = d <=> h^2 -y*h + d = 0 (h === q)
-
- D = y*y - 4*d;
- if(fabs(D) < eps) //in other words - D==0
- {
- q1 = q2 = y * 0.5;
- // g1+g2 = a && g1+g2 = b-y <=> g^2 - a*g + b-y = 0 (p === g)
- D = a*a - 4*(b-y);
- if(fabs(D) < eps) //in other words - D==0
- p1 = p2 = a * 0.5;
-
- else
- {
- sqD = sqrt(D);
- p1 = (a + sqD) * 0.5;
- p2 = (a - sqD) * 0.5;
- }
- }
- else
- {
- sqD = sqrt(D);
- q1 = (y + sqD) * 0.5;
- q2 = (y - sqD) * 0.5;
- // g1+g2 = a && g1*h2 + g2*h1 = c ( && g === p ) Krammer
- p1 = (a*q1-c)/(q1-q2);
- p2 = (c-a*q2)/(q1-q2);
- }
-
- DComplex* retval = new DComplex[4];
-
- // solving quadratic eq. - x^2 + p1*x + q1 = 0
- D = p1*p1 - 4*q1;
- if(D < 0.0)
- {
- retval[0].real( -p1 * 0.5 );
- retval[0].imag( sqrt(-D) * 0.5 );
- retval[1] = std::conj(retval[0]);
- }
- else
- {
- sqD = sqrt(D);
- retval[0].real( (-p1 + sqD) * 0.5 );
- retval[1].real( (-p1 - sqD) * 0.5 );
- }
-
- // solving quadratic eq. - x^2 + p2*x + q2 = 0
- D = p2*p2 - 4*q2;
- if(D < 0.0)
- {
- retval[2].real( -p2 * 0.5 );
- retval[2].imag( sqrt(-D) * 0.5 );
- retval[3] = std::conj(retval[2]);
- }
- else
- {
- sqD = sqrt(D);
- retval[2].real( (-p2 + sqD) * 0.5 );
- retval[3].real( (-p2 - sqD) * 0.5 );
- }
-
- return retval;
-}
+// Attention - this function returns dynamically allocated array. It has to be released
+// afterwards.
+DComplex* solve_quartic(double a, double b, double c, double d) {
+ double a3 = -b;
+ double b3 = a * c - 4. * d;
+ double c3 = -a * a * d - c * c + 4. * b * d;
+
+ // cubic resolvent
+ // y^3 − b*y^2 + (ac−4d)*y − a^2*d−c^2+4*b*d = 0
+
+ double x3[3];
+ unsigned int iZeroes = solveP3(x3, a3, b3, c3);
+
+ double q1, q2, p1, p2, D, sqD, y;
+
+ y = x3[0];
+ // The essence - choosing Y with maximal absolute value.
+ if (iZeroes != 1) {
+ if (fabs(x3[1]) > fabs(y)) y = x3[1];
+ if (fabs(x3[2]) > fabs(y)) y = x3[2];
+ }
+ // h1+h2 = y && h1*h2 = d <=> h^2 -y*h + d = 0 (h === q)
+
+ D = y * y - 4 * d;
+ if (fabs(D) < eps) // in other words - D==0
+ {
+ q1 = q2 = y * 0.5;
+ // g1+g2 = a && g1+g2 = b-y <=> g^2 - a*g + b-y = 0 (p === g)
+ D = a * a - 4 * (b - y);
+ if (fabs(D) < eps) // in other words - D==0
+ p1 = p2 = a * 0.5;
+
+ else {
+ sqD = sqrt(D);
+ p1 = (a + sqD) * 0.5;
+ p2 = (a - sqD) * 0.5;
+ }
+ } else {
+ sqD = sqrt(D);
+ q1 = (y + sqD) * 0.5;
+ q2 = (y - sqD) * 0.5;
+ // g1+g2 = a && g1*h2 + g2*h1 = c ( && g === p ) Krammer
+ p1 = (a * q1 - c) / (q1 - q2);
+ p2 = (c - a * q2) / (q1 - q2);
+ }
+
+ DComplex* retval = new DComplex[4];
+
+ // solving quadratic eq. - x^2 + p1*x + q1 = 0
+ D = p1 * p1 - 4 * q1;
+ if (D < 0.0) {
+ retval[0].real(-p1 * 0.5);
+ retval[0].imag(sqrt(-D) * 0.5);
+ retval[1] = std::conj(retval[0]);
+ } else {
+ sqD = sqrt(D);
+ retval[0].real((-p1 + sqD) * 0.5);
+ retval[1].real((-p1 - sqD) * 0.5);
+ }
+
+ // solving quadratic eq. - x^2 + p2*x + q2 = 0
+ D = p2 * p2 - 4 * q2;
+ if (D < 0.0) {
+ retval[2].real(-p2 * 0.5);
+ retval[2].imag(sqrt(-D) * 0.5);
+ retval[3] = std::conj(retval[2]);
+ } else {
+ sqD = sqrt(D);
+ retval[2].real((-p2 + sqD) * 0.5);
+ retval[3].real((-p2 - sqD) * 0.5);
+ }
+
+ return retval;
+}
diff --git a/Framework/Utilities/quartic.h b/Framework/Utilities/quartic.h
index 968ebbbd79ca170d8a7ae2ec27737240203b57b5..ac96afb10e57cff019344d650a1df141ad69b223 100644
--- a/Framework/Utilities/quartic.h
+++ b/Framework/Utilities/quartic.h
@@ -17,41 +17,38 @@
* Free Software Foundation, Inc., *
* 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. *
***************************************************************************/
-
+
#ifndef QUARTIC_H_INCLUDED
#define QUARTIC_H_INCLUDED
#include <complex>
const double PI = 3.141592653589793238463L;
-const double M_2PI = 2*PI;
-const double eps=1e-12;
+const double M_2PI = 2 * PI;
+const double eps = 1e-12;
typedef std::complex<double> DComplex;
//---------------------------------------------------------------------------
// useful for testing
- inline DComplex polinom_2(DComplex x, double a, double b)
- {
- //Horner's scheme for x*x + a*x + b
- return x * (x + a) + b;
- }
+inline DComplex polinom_2(DComplex x, double a, double b) {
+ // Horner's scheme for x*x + a*x + b
+ return x * (x + a) + b;
+}
//---------------------------------------------------------------------------
// useful for testing
- inline DComplex polinom_3(DComplex x, double a, double b, double c)
- {
- //Horner's scheme for x*x*x + a*x*x + b*x + c;
- return x * (x * (x + a) + b) + c;
- }
+inline DComplex polinom_3(DComplex x, double a, double b, double c) {
+ // Horner's scheme for x*x*x + a*x*x + b*x + c;
+ return x * (x * (x + a) + b) + c;
+}
//---------------------------------------------------------------------------
// useful for testing
- inline DComplex polinom_4(DComplex x, double a, double b, double c, double d)
- {
- //Horner's scheme for x*x*x*x + a*x*x*x + b*x*x + c*x + d;
- return x * (x * (x * (x + a) + b) + c) + d;
- }
+inline DComplex polinom_4(DComplex x, double a, double b, double c, double d) {
+ // Horner's scheme for x*x*x*x + a*x*x*x + b*x*x + c*x + d;
+ return x * (x * (x * (x + a) + b) + c) + d;
+}
//---------------------------------------------------------------------------
// x - array of size 3
@@ -62,8 +59,8 @@ unsigned int solveP3(double* x, double a, double b, double c);
//---------------------------------------------------------------------------
// solve quartic equation x^4 + a*x^3 + b*x^2 + c*x + d
-// Attention - this function returns dynamically allocated array. It has to be released afterwards.
+// Attention - this function returns dynamically allocated array. It has to be released
+// afterwards.
DComplex* solve_quartic(double a, double b, double c, double d);
-
#endif // QUARTIC_H_INCLUDED
diff --git a/Processes/ObservationPlane/testObservationPlane.cc b/Processes/ObservationPlane/testObservationPlane.cc
index 3010d34b1faabbaedfe547ff5e98fe0d6cbc0d2b..ac5573e0b2df8552bd09ed246793e9452d75712b 100644
--- a/Processes/ObservationPlane/testObservationPlane.cc
+++ b/Processes/ObservationPlane/testObservationPlane.cc
@@ -43,10 +43,10 @@ TEST_CASE("ContinuousProcess interface", "[proccesses][observation_plane]") {
*/
auto [stack, viewPtr] =
- setup::testing::setupStack(particles::Code::NuE, 0, 0, 1_GeV, nodePtr, cs);
+ setup::testing::setupStack(particles::Code::NuE, 0, 0, 1_GeV, nodePtr, cs);
[[maybe_unused]] setup::StackView& view = *viewPtr;
- auto particle = stack->GetNextParticle();
-
+ auto particle = stack->GetNextParticle();
+
Point const start(cs, {0_m, 1_m, 10_m});
Vector<units::si::SpeedType::dimension_type> vec(cs, 0_m / second, 0_m / second,
-units::constants::c);
diff --git a/Processes/TrackWriter/TrackWriter.cc b/Processes/TrackWriter/TrackWriter.cc
index 1bef877c3a48d327e56de4a70b6f13df5154f117..2d9825d869f69a1da7956c7631e8ffdb6d85fea2 100644
--- a/Processes/TrackWriter/TrackWriter.cc
+++ b/Processes/TrackWriter/TrackWriter.cc
@@ -28,8 +28,9 @@ namespace corsika::process::track_writer {
using namespace std::string_literals;
fFile.open(fFilename);
- fFile << "# PID, E / eV, start coordinates / m, displacement vector to end / m, steplength / m "s
- << '\n';
+ fFile
+ << "# PID, E / eV, start coordinates / m, displacement vector to end / m, steplength / m "s
+ << '\n';
}
template <>
diff --git a/Processes/Tracking/Intersect.hpp b/Processes/Tracking/Intersect.hpp
index 19e4819bf704f5a2c4013ee8d0a5ae812707f9f0..277b66f0252931d82b468c781509938eba141d36 100644
--- a/Processes/Tracking/Intersect.hpp
+++ b/Processes/Tracking/Intersect.hpp
@@ -1,3 +1,11 @@
+/*
+ * (c) Copyright 2020 CORSIKA Project, corsika-project@lists.kit.edu
+ *
+ * This software is distributed under the terms of the GNU General Public
+ * Licence version 3 (GPL Version 3). See file LICENSE for a full version of
+ * the license.
+ */
+
#pragma once
#include <corsika/geometry/Point.h>
diff --git a/Processes/Tracking/testTracking.cc b/Processes/Tracking/testTracking.cc
index fffe1ba2a3410ae07981557f8acdd935e7552a87..aabd87603fdc888fbc4d5a2a3a39fde50a26adb1 100644
--- a/Processes/Tracking/testTracking.cc
+++ b/Processes/Tracking/testTracking.cc
@@ -62,8 +62,8 @@ TEMPLATE_TEST_CASE("TrackingLeapfrog_Curved", "tracking",
[]([[maybe_unused]] MagneticFluxType v) {
if constexpr (std::is_same_v<TestType, tracking_line::Tracking>)
return v == 0_uT;
- else
- return true;
+ else
+ return true;
},
values<MagneticFluxType>({50_uT, 0_uT, -50_uT})));
// particle --> (world) --> | --> (target)
diff --git a/Processes/TrackingLeapFrogCurved/Tracking.h b/Processes/TrackingLeapFrogCurved/Tracking.h
index f97ec44e099b77c80154300805c9b4a68b447d89..63e1dad56fc512b941c122b41cc09358e0039937 100644
--- a/Processes/TrackingLeapFrogCurved/Tracking.h
+++ b/Processes/TrackingLeapFrogCurved/Tracking.h
@@ -182,7 +182,7 @@ namespace corsika::process {
if (first != 2) {
C8LOG_DEBUG("no intersection! count={}", first);
- return geometry::Intersections();
+ return geometry::Intersections();
}
return geometry::Intersections(d_enter / absVelocity, d_exit / absVelocity);
}
diff --git a/Processes/TrackingLeapFrogCurved/testTrackingLeapFrogCurved.cc b/Processes/TrackingLeapFrogCurved/testTrackingLeapFrogCurved.cc
index cf45ccc48721a851e491eda2c0b8b5ef923c329f..2a76dbac8d58e154492a70e6d476e684a2912136 100644
--- a/Processes/TrackingLeapFrogCurved/testTrackingLeapFrogCurved.cc
+++ b/Processes/TrackingLeapFrogCurved/testTrackingLeapFrogCurved.cc
@@ -57,7 +57,7 @@ TEST_CASE("TrackingLeapfrog_Curved") {
if (chargeNumber != 0 and Bfield != 0_T) {
deflect = -sgn(chargeNumber) * sgn(Bfield / 1_T); // direction of deflection
LengthType const gyroradius =
- P0 * 1_V / (constants::c * abs(chargeNumber) * abs(Bfield) * 1_eV);
+ P0 * 1_V / (constants::c * abs(chargeNumber) * abs(Bfield) * 1_eV);
radius = gyroradius;
}
@@ -115,13 +115,14 @@ TEST_CASE("TrackingLeapfrog_Curved") {
particle.SetMomentum(traj2.GetDirection(1) * particle.GetMomentum().norm());
}
CHECK(nextVol == worldPtr);
-
+
Point pointCheck(cs, (deflect == 0 ? radius : 0_m), (deflect * radius), 0_m);
C8LOG_DEBUG("testTrackingLineStack: deflect={}, momentum={}, pos={}, pos_check={}",
deflect, particle.GetMomentum().GetComponents(),
particle.GetPosition().GetCoordinates(), pointCheck.GetCoordinates());
- CHECK((particle.GetPosition() - pointCheck).norm() / radius == Approx(0).margin(1e-3));
+ CHECK((particle.GetPosition() - pointCheck).norm() / radius ==
+ Approx(0).margin(1e-3));
}
}
diff --git a/Processes/TrackingLeapFrogCurved/testTrackingLeapFrogStraight.cc b/Processes/TrackingLeapFrogCurved/testTrackingLeapFrogStraight.cc
index c96c2ffc84ae3f126c0b74e2fae3ce8e2a534dd2..895a185866dbd9d9b0fc2e6efffbcd82413e349a 100644
--- a/Processes/TrackingLeapFrogCurved/testTrackingLeapFrogStraight.cc
+++ b/Processes/TrackingLeapFrogCurved/testTrackingLeapFrogStraight.cc
@@ -57,7 +57,7 @@ TEST_CASE("TrackingBField") {
if (chargeNumber != 0 and Bfield != 0_T) {
deflect = -sgn(chargeNumber) * sgn(Bfield / 1_T); // direction of deflection
LengthType const gyroradius =
- P0 * 1_V / (constants::c * abs(chargeNumber) * abs(Bfield) * 1_eV);
+ P0 * 1_V / (constants::c * abs(chargeNumber) * abs(Bfield) * 1_eV);
radius = gyroradius;
}
@@ -122,6 +122,7 @@ TEST_CASE("TrackingBField") {
deflect, particle.GetMomentum().GetComponents(),
particle.GetPosition().GetCoordinates(), pointCheck.GetCoordinates());
- CHECK((particle.GetPosition() - pointCheck).norm() / radius == Approx(0).margin(1e-3));
+ CHECK((particle.GetPosition() - pointCheck).norm() / radius ==
+ Approx(0).margin(1e-3));
}
}
diff --git a/Processes/TrackingLine/Tracking.cc b/Processes/TrackingLine/Tracking.cc
index 01cd6c134a1b3401a6e2204afb1400d1da6087c2..6529bf8b705594809f8ab1515ec4d2891a08aa66 100644
--- a/Processes/TrackingLine/Tracking.cc
+++ b/Processes/TrackingLine/Tracking.cc
@@ -21,6 +21,4 @@
using namespace corsika::geometry;
using namespace corsika::units::si;
-namespace corsika::process::tracking_line {
-
-} // namespace corsika::process::tracking_line
+namespace corsika::process::tracking_line {} // namespace corsika::process::tracking_line
diff --git a/Processes/TrackingLine/Tracking.h b/Processes/TrackingLine/Tracking.h
index 799f2194f9f0dd752dc42951b4095cc621aa0a1e..ab68080f9bfa6573c8f5185cb0d5b7aa0a18536d 100644
--- a/Processes/TrackingLine/Tracking.h
+++ b/Processes/TrackingLine/Tracking.h
@@ -54,9 +54,9 @@ namespace corsika::process {
particle.GetMomentum().GetComponents() / 1_GeV);
C8LOG_DEBUG("Tracking v: {} ", initialVelocity.GetComponents());
- // traverse the environment volume tree and find next
- // intersection
- auto [minTime, minNode] = tracking::Intersect<Tracking>::nextIntersect(particle);
+ // traverse the environment volume tree and find next
+ // intersection
+ auto [minTime, minNode] = tracking::Intersect<Tracking>::nextIntersect(particle);
return std::make_tuple(
geometry::LineTrajectory(geometry::Line(initialPosition, initialVelocity),
@@ -64,61 +64,58 @@ namespace corsika::process {
minNode); // next volume node
}
-
-
template <typename TParticle, typename TMedium>
static geometry::Intersections Intersect(const TParticle& particle,
- const corsika::geometry::Sphere& sphere,
- const TMedium&) {
- using namespace corsika::units::si;
- auto const delta = particle.GetPosition() - sphere.GetCenter();
- auto const velocity =
- particle.GetMomentum() / particle.GetEnergy() * corsika::units::constants::c;
- auto const vSqNorm = velocity.squaredNorm();
- auto const R = sphere.GetRadius();
-
- auto const vDotDelta = velocity.dot(delta);
- auto const discriminant =
- vDotDelta * vDotDelta - vSqNorm * (delta.squaredNorm() - R * R);
-
- if (discriminant.magnitude() > 0) {
- auto const sqDisc = sqrt(discriminant);
- auto const invDenom = 1 / vSqNorm;
- return geometry::Intersections((-vDotDelta - sqDisc) * invDenom,
- (-vDotDelta + sqDisc) * invDenom);
+ const corsika::geometry::Sphere& sphere,
+ const TMedium&) {
+ using namespace corsika::units::si;
+ auto const delta = particle.GetPosition() - sphere.GetCenter();
+ auto const velocity =
+ particle.GetMomentum() / particle.GetEnergy() * corsika::units::constants::c;
+ auto const vSqNorm = velocity.squaredNorm();
+ auto const R = sphere.GetRadius();
+
+ auto const vDotDelta = velocity.dot(delta);
+ auto const discriminant =
+ vDotDelta * vDotDelta - vSqNorm * (delta.squaredNorm() - R * R);
+
+ if (discriminant.magnitude() > 0) {
+ auto const sqDisc = sqrt(discriminant);
+ auto const invDenom = 1 / vSqNorm;
+ return geometry::Intersections((-vDotDelta - sqDisc) * invDenom,
+ (-vDotDelta + sqDisc) * invDenom);
+ }
+ return geometry::Intersections();
}
- return geometry::Intersections();
- }
-
- template <typename TParticle, typename TBaseNodeType>
- static geometry::Intersections Intersect(const TParticle& particle,
- const TBaseNodeType& volumeNode) {
- const geometry::Sphere* sphere =
- dynamic_cast<const geometry::Sphere*>(&volumeNode.GetVolume());
- if (sphere) {
- return Intersect(particle, *sphere, volumeNode.GetModelProperties());
+
+ template <typename TParticle, typename TBaseNodeType>
+ static geometry::Intersections Intersect(const TParticle& particle,
+ const TBaseNodeType& volumeNode) {
+ const geometry::Sphere* sphere =
+ dynamic_cast<const geometry::Sphere*>(&volumeNode.GetVolume());
+ if (sphere) {
+ return Intersect(particle, *sphere, volumeNode.GetModelProperties());
+ }
+ throw std::runtime_error(
+ "The Volume type provided is not supported in Intersect(particle, node)");
}
- throw std::runtime_error(
- "The Volume type provided is not supported in Intersect(particle, node)");
- }
-
-
- template <typename TParticle, typename TMedium>
- static geometry::Intersections Intersect(const TParticle& particle,
- const geometry::Plane& plane,
- const TMedium& medium) {
- using namespace corsika::units::si;
- auto const delta = plane.GetCenter() - particle.GetPosition();
- auto const velocity =
- particle.GetMomentum() / particle.GetEnergy() * corsika::units::constants::c;
- auto const n = plane.GetNormal();
- auto const c = n.dot(velocity);
-
- return Intersections(
- c.magnitude() == 0 ? std::numeric_limits<TimeType::value_type>::infinity() * 1_s
- : n.dot(delta) / c);
- }
+ template <typename TParticle, typename TMedium>
+ static geometry::Intersections Intersect(const TParticle& particle,
+ const geometry::Plane& plane,
+ const TMedium& medium) {
+ using namespace corsika::units::si;
+ auto const delta = plane.GetCenter() - particle.GetPosition();
+ auto const velocity =
+ particle.GetMomentum() / particle.GetEnergy() * corsika::units::constants::c;
+ auto const n = plane.GetNormal();
+ auto const c = n.dot(velocity);
+
+ return Intersections(c.magnitude() == 0
+ ? std::numeric_limits<TimeType::value_type>::infinity() *
+ 1_s
+ : n.dot(delta) / c);
+ }
};
} // namespace tracking_line
diff --git a/Stack/History/HistoryObservationPlane.cpp b/Stack/History/HistoryObservationPlane.cpp
index eb5092d1f77f0e59a605adbad7ec368a6b01014e..6a8c94ec7d9e33686383dadf479ef56e413d5522 100644
--- a/Stack/History/HistoryObservationPlane.cpp
+++ b/Stack/History/HistoryObservationPlane.cpp
@@ -28,12 +28,13 @@ HistoryObservationPlane::HistoryObservationPlane(setup::Stack const& stack,
corsika::process::EProcessReturn HistoryObservationPlane::DoContinuous(
setup::Stack::ParticleType const& particle, setup::Trajectory const& trajectory) {
TimeType const timeOfIntersection =
- (plane_.GetCenter() - trajectory.GetLine().GetR0()).dot(plane_.GetNormal()) /
+ (plane_.GetCenter() - trajectory.GetLine().GetR0()).dot(plane_.GetNormal()) /
trajectory.GetLine().GetV0().dot(plane_.GetNormal());
if (timeOfIntersection < TimeType::zero()) { return process::EProcessReturn::eOk; }
- if (plane_.IsAbove(trajectory.GetLine().GetR0()) == plane_.IsAbove(trajectory.GetPosition(1))) {
+ if (plane_.IsAbove(trajectory.GetLine().GetR0()) ==
+ plane_.IsAbove(trajectory.GetPosition(1))) {
return process::EProcessReturn::eOk;
}
diff --git a/do-copyright.py b/do-copyright.py
index 6b235afce2a7ccd03f342ce0597ed07c9ff354e5..adc278d1124eafdc943d65937ecdb9f552933a40 100755
--- a/do-copyright.py
+++ b/do-copyright.py
@@ -192,10 +192,9 @@ def next_file(dir_name, files, justCheck, forYear, updateMessage):
excludes if wished, process otherwise
"""
for check in excludeDirs :
- print (dir_name)
if check in dir_name:
- # if Debug>1:
- print ("exclude-dir: " + check)
+ if Debug>1:
+ print ("exclude-dir: " + check)
return True
for check in files :
if (os.path.isdir(check)):