From cc0c7d1b0cf436f56fd4ec46bd7f8ea00ad8778a Mon Sep 17 00:00:00 2001
From: Andre Schmidt <schmidt-a@iklx288.ikp.kit.edu>
Date: Sat, 5 Sep 2020 16:40:53 +0200
Subject: [PATCH] histograms finished
---
Framework/Cascade/Cascade.h | 38 +++-
Processes/TrackingLine/TrackingLine.h | 271 ++++++++++++++++++++------
2 files changed, 236 insertions(+), 73 deletions(-)
diff --git a/Framework/Cascade/Cascade.h b/Framework/Cascade/Cascade.h
index df9ee0e18..2fd8de8da 100644
--- a/Framework/Cascade/Cascade.h
+++ b/Framework/Cascade/Cascade.h
@@ -44,15 +44,18 @@ using namespace boost::histogram;
static auto histL2 = make_histogram(axis::regular<>(100, 0, 60000, "L'"));
static auto histS2 = make_histogram(axis::regular<>(100, 0, 60000, "S"));
static auto histB2 = make_histogram(axis::regular<>(100, 0, 60000, "Bogenlänge"));
+static auto histLlog2 = make_histogram(axis::regular<double, axis::transform::log>(100, 1, 1e7, "Leap-Frog-ength L'"));
+static auto histSlog2 = make_histogram(axis::regular<double, axis::transform::log>(100, 1, 1e7, "Direct Length S"));
+static auto histBlog2 = make_histogram(axis::regular<double, axis::transform::log>(100, 1, 1e7, "Arc Length B"));
static auto histLB2 = make_histogram(axis::regular<>(100, 0, 0.01, "L - B"));
static auto histLS2 = make_histogram(axis::regular<>(100, 0, 0.01, "L - S"));
static auto histLBrel2 = make_histogram(axis::regular<double, axis::transform::log> (20,1e-11,1e-6,"L/B -1"));
static auto histLSrel2 = make_histogram(axis::regular<double, axis::transform::log>(20,1e-11,1e-6, "L/S -1"));
static auto histELSrel2 = make_histogram(axis::regular<double, axis::transform::log>(20,1e-11,1e-6, "L/S -1"),axis::regular<double, axis::transform::log>(20, 0.1, 1e4, "E / GeV"));
static auto histBS2 = make_histogram(axis::regular<>(100, 0, 0.01, "B - S"));
-static auto histLp2 = make_histogram(axis::regular<>(100, 0, 60000, "L' für Protonen"));
-static auto histLpi2 = make_histogram(axis::regular<>(100, 0, 60000, "L' für Pionen"));
-static auto histLmu2 = make_histogram(axis::regular<>(100, 0, 60000, "L' für Myonen"));
+static auto histLp2 = make_histogram(axis::regular<double, axis::transform::log>(100, 1, 1e7, "L' für Protonen"));
+static auto histLpi2 = make_histogram(axis::regular<double, axis::transform::log>(100, 1, 1e7, "L' für Pionen"));
+static auto histLmu2 = make_histogram(axis::regular<double, axis::transform::log>(100, 1, 1e7, "L' für Myonen"));
//static auto histLe = make_histogram(axis::regular<>(100, 0, 60000, "L' für Elektronen"));
//static auto histLy = make_histogram(axis::regular<>(100, 0, 60000, "L' für Photonen"));
@@ -130,11 +133,10 @@ namespace corsika::cascade {
~Cascade(){
std::ofstream myfile;
myfile.open ("histograms2.txt");
- myfile << histLB2 << std::endl;
- myfile << histLBrel2 << std::endl;
- myfile << histLS2 << std::endl;
- myfile << histLSrel2 << std::endl;
- myfile << histELSrel2 << std::endl;
+ myfile << histLp2 << std::endl;
+ myfile << histLpi2 << std::endl;
+ myfile << histLmu2 << std::endl;
+ myfile << histL2 << std::endl;
myfile.close();
/*std::cout << histLBrel << std::endl;
@@ -178,6 +180,15 @@ namespace corsika::cascade {
std::ofstream file13("histLp2.json");
dump_bh(file13, histLp2);
file13.close();
+ std::ofstream file14("histLlog2.json");
+ dump_bh(file14, histLlog2);
+ file14.close();
+ std::ofstream file15("histBlog2.json");
+ dump_bh(file15, histBlog2);
+ file15.close();
+ std::ofstream file16("histSlog2.json");
+ dump_bh(file16, histSlog2);
+ file16.close();
};
@@ -352,7 +363,16 @@ namespace corsika::cascade {
// This formula has an error or doesnt work for specific conditions
// Steplength should not be min_distance
- auto [position, direction] = fTracking.MagneticStep(vParticle, min_distance);
+ auto [position, direction, L2] = fTracking.MagneticStep(vParticle, min_distance);
+ histL2(L2);
+ histLlog2(L2);
+ int pdg = static_cast<int>(particles::GetPDG(vParticle.GetPID()));
+ if (abs(pdg) == 13)
+ histLmu2(L2);
+ if (abs(pdg) == 211 || abs(pdg) == 111)
+ histLpi2(L2);
+ if (abs(pdg) == 2212 || abs(pdg) == 2112)
+ histLp2(L2);
auto distance = position - vParticle.GetPosition();
//Building Trajectory for Continuous processes
diff --git a/Processes/TrackingLine/TrackingLine.h b/Processes/TrackingLine/TrackingLine.h
index 875569079..cc283646b 100644
--- a/Processes/TrackingLine/TrackingLine.h
+++ b/Processes/TrackingLine/TrackingLine.h
@@ -34,26 +34,39 @@ namespace corsika::environment {
} // namespace corsika::environment
using namespace boost::histogram;
-static auto histL = make_histogram(axis::regular<>(100, 0, 60000, "L'"));
-static auto histS = make_histogram(axis::regular<>(100, 0, 60000, "S"));
-static auto histB = make_histogram(axis::regular<>(100, 0, 60000, "Bogenlänge"));
-static auto histLB = make_histogram(axis::regular<>(100, 0, 0.01, "L - B"));
-static auto histLS = make_histogram(axis::regular<>(100, 0, 0.01, "L - S"));
-static auto histLBrel = make_histogram(axis::regular<double, axis::transform::log> (20,1e-11,1e-6,"L/B -1"));
-static auto histLSrel = make_histogram(axis::regular<double, axis::transform::log>(20,1e-11,1e-6, "L/S -1"));
-static auto histELSrel = make_histogram(axis::regular<double, axis::transform::log>(20,1e-11,1e-6, "L/S -1"),axis::regular<double, axis::transform::log>(20, 0.1, 1e4, "E / GeV"));
+static auto histL = make_histogram(axis::regular<>(100, 0, 100000, "Leap-Frog-ength L'"));
+static auto histS = make_histogram(axis::regular<>(100, 0, 100000, "Direct Length S"));
+static auto histB = make_histogram(axis::regular<>(100, 0, 100000, "Arc Length B"));
+static auto histLlog = make_histogram(axis::regular<double, axis::transform::log>(100, 1, 1e7, "Leap-Frog-ength L'"));
+static auto histSlog = make_histogram(axis::regular<double, axis::transform::log>(100, 1, 1e7, "Direct Length S"));
+static auto histBlog = make_histogram(axis::regular<double, axis::transform::log>(100, 1, 1e7, "Arc Length B"));
+static auto histLB = make_histogram(axis::regular<>(100, 0, 100, "L - B"));
+static auto histLS = make_histogram(axis::regular<>(100, 0, 100, "L - S"));
+static auto histLBrel = make_histogram(axis::regular<double, axis::transform::log> (40,1e-12,1e-2,"L/B -1"));
+static auto histLSrel = make_histogram(axis::regular<double, axis::transform::log>(40,1e-12,1e-2, "L/S -1"));
+static auto histELSrel = make_histogram(axis::regular<double, axis::transform::log>(40,1e-12,1e-4, "L/S -1"),axis::regular<double, axis::transform::log>(30, 3, 3e1, "E / GeV"));
static auto histBS = make_histogram(axis::regular<>(100, 0, 0.01, "B - S"));
-static auto histLp = make_histogram(axis::regular<>(100, 0, 60000, "L' für Protonen"));
-static auto histLpi = make_histogram(axis::regular<>(100, 0, 60000, "L' für Pionen"));
-static auto histLmu = make_histogram(axis::regular<>(100, 0, 60000, "L' für Myonen"));
+static auto histLp = make_histogram(axis::regular<double, axis::transform::log>(100, 1, 1e7, "L' für Protonen"));
+static auto histLpi = make_histogram(axis::regular<double, axis::transform::log>(100, 1, 1e7, "L' für Pionen"));
+static auto histLmu = make_histogram(axis::regular<double, axis::transform::log>(100, 1, 1e7, "L' für Myonen"));
//static auto histLe = make_histogram(axis::regular<>(100, 0, 60000, "L' für Elektronen"));
//static auto histLy = make_histogram(axis::regular<>(100, 0, 60000, "L' für Photonen"));
static double L = 0;
+static std::vector<double> Lsmall;
+static std::vector<double> Bbig;
+static std::vector<double> Sbig;
+static std::vector<double> Lsmallmag;
+static std::vector<double> Bbigmag;
+static std::vector<double> Sbigmag;
+static std::vector<double> Emag;
+static std::vector<double> E;
+static std::vector<double> Steplimitmag;
+
namespace corsika::process {
namespace tracking_line {
-
+
std::optional<std::pair<corsika::units::si::TimeType, corsika::units::si::TimeType>>
TimeOfIntersection(geometry::Line const&, geometry::Sphere const&);
@@ -67,16 +80,53 @@ namespace corsika::process {
~TrackingLine(){
std::ofstream myfile;
myfile.open ("histograms.txt");
+ myfile << histLlog << std::endl;
+ myfile << histBlog << std::endl;
myfile << histLB << std::endl;
myfile << histLBrel << std::endl;
myfile << histLS << std::endl;
myfile << histLSrel << std::endl;
- myfile << histELSrel << std::endl;
myfile.close();
-
- /*std::cout << histLBrel << std::endl;
- std::cout << histLSrel << std::endl;*/
-
+
+ std::ofstream myfile2;
+ myfile2.open ("lengen.txt");
+ myfile2 << "L " << std::endl;
+ for(double n : Lsmall) {
+ myfile2 << n << '\n';
+ }
+ myfile2 << "B " << std::endl;
+ for(double n : Bbig) {
+ myfile2 << n << '\n';
+ }
+ myfile2 << "S " << std::endl;
+ for(double n : Sbig) {
+ myfile2 << n << '\n';
+ }
+ myfile2 << "E in GeV" << std::endl;
+ for(double n : E) {
+ myfile2 << n << '\n';
+ }
+ myfile2 << "L mag" << std::endl;
+ for(double n : Lsmallmag) {
+ myfile2 << n << '\n';
+ }
+ myfile2 << "B mag" << std::endl;
+ for(double n : Bbigmag) {
+ myfile2 << n << '\n';
+ }
+ myfile2 << "S mag" << std::endl;
+ for(double n : Sbigmag) {
+ myfile2 << n << '\n';
+ }
+ myfile2 << "E mag in GeV" << std::endl;
+ for(double n : Emag) {
+ myfile2 << n << '\n';
+ }
+ myfile2 << "Schrittlänge" << std::endl;
+ for(double n : Steplimitmag) {
+ myfile2 << n << '\n';
+ }
+ myfile2.close();
std::ofstream file1("histL.json");
dump_bh(file1, histL);
@@ -115,6 +165,15 @@ namespace corsika::process {
std::ofstream file13("histLp.json");
dump_bh(file13, histLp);
file13.close();
+ std::ofstream file14("histLlog.json");
+ dump_bh(file14, histLlog);
+ file14.close();
+ std::ofstream file15("histBlog.json");
+ dump_bh(file15, histBlog);
+ file15.close();
+ std::ofstream file16("histSlog.json");
+ dump_bh(file16, histSlog);
+ file16.close();
};
@@ -154,7 +213,7 @@ namespace corsika::process {
std::vector<std::pair<TimeType, decltype(p.GetNode())>> intersections;
//charge of the particle
- int chargeNumber;
+ int chargeNumber = 0;
if (corsika::particles::IsNucleus(p.GetPID())) {
chargeNumber = p.GetNuclearZ();
} else {
@@ -162,13 +221,12 @@ namespace corsika::process {
}
auto magneticfield = currentLogicalVolumeNode->GetModelProperties().GetMagneticField(currentPosition);
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();
- LengthType Steplength1 = 10_m; // length irrelevant if q=0 and else it gets changed again
- LengthType Steplength2 = 10_m;
+ LengthType Steplength1 = 0_m;
+ LengthType Steplength2 = 0_m;
LengthType Steplimit = 1_m / 0;
- //infinite kann man auch anders schreiben
-
+ //infinite kann man auch anders schreiben
+ bool intersection = true;
+
// for entering from outside
auto addIfIntersects = [&](auto const& vtn) {
auto const& volume = vtn.GetVolume();
@@ -177,7 +235,9 @@ namespace corsika::process {
// everything is a sphere, crashes with exception if not
// creating Line with magnetic field
- if (chargeNumber != 0) {
+ if (chargeNumber != 0) {
+ auto k = chargeNumber * corsika::units::constants::cSquared * 1_eV / (velocity.norm() * p.GetEnergy() * 1_V);
+ geometry::Vector<dimensionless_d> const directionBefore = velocity.normalized();
// determine steplength to next volume
double a = ((directionBefore.cross(magneticfield)).dot(currentPosition - sphere.GetCenter()) * k + 1) * 4 /
(1_m * 1_m * (directionBefore.cross(magneticfield)).GetSquaredNorm() * k * k);
@@ -200,29 +260,27 @@ namespace corsika::process {
std::cout << "Time to next volume = " << Steplength1 / velocity.norm() << std::endl;
} else {
std::cout << "no intersection (1)!" << std::endl;
- // what to do when this happens? (very unlikely)
+ intersection = false;
}
delete [] solutions;
-
- geometry::Vector<SpeedType::dimension_type> const velocityVerticalMag = velocity -
- velocity.parallelProjectionOnto(magneticfield);
- LengthType const gyroradius = p.GetEnergy() * velocityVerticalMag.GetNorm() * 1_V /
- (corsika::units::constants::cSquared * abs(chargeNumber) *
- magneticfield.GetNorm() * 1_eV);
- Steplimit = 2 * sin(0.1) * gyroradius;
}
-
- auto line1 = MagneticStep(p, line, Steplength1, false);
- // instead of changing the line with magnetic field, the TimeOfIntersection() could be changed
- // using line has some errors for huge steps
-
- if (auto opt = TimeOfIntersection(line1, sphere); opt.has_value()) {
- auto const [t1, t2] = *opt;
- C8LOG_DEBUG("intersection times: {} s; {} s", t1 / 1_s, t2 / 1_s);
- if (t1.magnitude() > 0)
- intersections.emplace_back(t1, &vtn);
- else if (t2.magnitude() > 0)
- C8LOG_DEBUG("inside other volume");
+
+ if (intersection) {
+ auto line1 = MagneticStep(p, line, Steplength1, false);
+ // instead of changing the line with magnetic field, the TimeOfIntersection() could be changed
+ // using line has some errors for huge steps
+
+ if (auto opt = TimeOfIntersection(line1, sphere); opt.has_value()) {
+ auto const [t1, t2] = *opt;
+ std::cout << "intersection times: " << t1 / 1_s << "; "
+ << t2 / 1_s
+ // << " " << vtn.GetModelProperties().GetName()
+ << std::endl;
+ if (t1.magnitude() > 0)
+ intersections.emplace_back(t1, &vtn);
+ else if (t2.magnitude() > 0)
+ std::cout << "inside other volume" << std::endl;
+ }
}
};
@@ -237,6 +295,8 @@ namespace corsika::process {
// creating Line with magnetic field
if (chargeNumber != 0) {
+ auto k = chargeNumber * corsika::units::constants::cSquared * 1_eV / (velocity.norm() * p.GetEnergy() * 1_V);
+ geometry::Vector<dimensionless_d> const directionBefore = velocity.normalized();
// determine steplength to next volume
double a = ((directionBefore.cross(magneticfield)).dot(currentPosition - sphere.GetCenter()) * k + 1) * 4 /
(1_m * 1_m * (directionBefore.cross(magneticfield)).GetSquaredNorm() * k * k);
@@ -267,6 +327,14 @@ namespace corsika::process {
// what to do when this happens? (very unlikely)
}
delete [] solutions;
+
+ geometry::Vector<SpeedType::dimension_type> const velocityVerticalMag = velocity -
+ velocity.parallelProjectionOnto(magneticfield);
+ LengthType const gyroradius = p.GetEnergy() * velocityVerticalMag.GetNorm() * 1_V /
+ (corsika::units::constants::cSquared * abs(chargeNumber) *
+ magneticfield.GetNorm() * 1_eV);
+ Steplimit = 2 * sin(0.1) * gyroradius;
+ std::cout << "Steplimit: " << Steplimit << std::endl;
}
auto line2 = MagneticStep(p, line, Steplength2, false);
@@ -296,30 +364,101 @@ namespace corsika::process {
<< min
// << " " << minIter->second->GetModelProperties().GetName()
<< std::endl;
-
+
+ //if the used Steplengths are too big, min gets false and we do another Step
+ std::cout << Steplength1 << Steplength2 << std::endl;
+ std::cout << intersection << " " << Steplimit << std::endl;
+ if ((Steplength1 > Steplimit || Steplength1 == 0_m) && Steplength2 > Steplimit) {
+ min = Steplimit / velocity.norm();
+ auto lineWithB = MagneticStep(p, line, Steplimit, true);
+
+ histL(L);
+ histLlog(L);
+ histS(lineWithB.ArcLength(0_s,min) / 1_m);
+ histSlog(lineWithB.ArcLength(0_s,min) / 1_m);
+ histLS(L-lineWithB.ArcLength(0_s,min) / 1_m);
+
+ if(chargeNumber != 0) {
+ if(L * 1_m/lineWithB.ArcLength(0_s,min) -1 > 0) {
+ histLSrel(L * 1_m/lineWithB.ArcLength(0_s,min) -1);
+ //histELSrel(L * 1_m/lineWithB.ArcLength(0_s,min) -1, p.GetEnergy() / 1_GeV);
+ } else {
+Lsmallmag.push_back(L); Sbigmag.push_back(lineWithB.ArcLength(0_s,min) / 1_m);
+Emag.push_back(p.GetEnergy() / 1_GeV);Steplimitmag.push_back(Steplimit / 1_m);
+ }
+ auto magneticfield = currentLogicalVolumeNode->GetModelProperties().GetMagneticField(currentPosition);
+ geometry::Vector<SpeedType::dimension_type> const velocityVerticalMag = velocity -
+ velocity.parallelProjectionOnto(magneticfield);
+ LengthType const gyroradius = p.GetEnergy() * velocityVerticalMag.GetNorm() * 1_V /
+ (corsika::units::constants::cSquared * abs(chargeNumber) *
+ magneticfield.GetNorm() * 1_eV);
+ std::cout << "Schrittlaenge " << velocity.norm() * min << " gyroradius " << gyroradius << std::endl;
+ LengthType B = 2 * gyroradius * asin ( lineWithB.ArcLength(0_s,min) / 2 / gyroradius);
+ std::cout << "Bogenlaenge" << B << std::endl;
+ histB(B / 1_m);
+ histBlog(B / 1_m);
+ if(L-B/1_m > 0) {
+ histLB(L-B/1_m);
+ histBS(B/1_m - lineWithB.ArcLength(0_s,min) / 1_m);
+ histLBrel(L * 1_m/B-1);
+ } else {
+Bbigmag.push_back(B / 1_m);
+ }
+ } else {
+ histB(lineWithB.ArcLength(0_s,min) / 1_m);
+ histBlog(lineWithB.ArcLength(0_s,min) / 1_m);
+ }
+
+ int pdg = static_cast<int>(particles::GetPDG(p.GetPID()));
+ if (abs(pdg) == 13)
+ histLmu(L);
+ if (abs(pdg) == 211 || abs(pdg) == 111)
+ histLpi(L);
+ if (abs(pdg) == 2212 || abs(pdg) == 2112)
+ histLp(L);
+
+ return std::make_tuple(geometry::Trajectory<geometry::Line>(line, min),
+ geometry::Trajectory<geometry::Line>(lineWithB, min),
+ Steplimit * 1.0001, Steplimit, minIter->second);
+ }
+
auto lineWithB = MagneticStep(p, line, velocity.norm() * min, true);
- if(min * velocity.norm() < 1000000_m) {
histL(L);
+ histLlog(L);
histS(lineWithB.ArcLength(0_s,min) / 1_m);
+ histSlog(lineWithB.ArcLength(0_s,min) / 1_m);
histLS(L-lineWithB.ArcLength(0_s,min) / 1_m);
if(chargeNumber != 0) {
- histLSrel(L * 1_m/lineWithB.ArcLength(0_s,min) -1);
- histELSrel(L * 1_m/lineWithB.ArcLength(0_s,min) -1, p.GetEnergy() / 1_GeV);
- auto magneticfield = currentLogicalVolumeNode->GetModelProperties().GetMagneticField(currentPosition);
- geometry::Vector<SpeedType::dimension_type> const velocityVerticalMag = velocity -
- velocity.parallelProjectionOnto(magneticfield);
-
-
+ if(L * 1_m/lineWithB.ArcLength(0_s,min) -1 > 0) {
+ histLSrel(L * 1_m/lineWithB.ArcLength(0_s,min) -1);
+ histELSrel(L * 1_m/lineWithB.ArcLength(0_s,min) -1, p.GetEnergy() / 1_GeV);
+ } else {
+Lsmall.push_back(L); Sbig.push_back(lineWithB.ArcLength(0_s,min) / 1_m);
+E.push_back(p.GetEnergy() / 1_GeV);
+ }
+ auto magneticfield = currentLogicalVolumeNode->GetModelProperties().GetMagneticField(currentPosition);
+ geometry::Vector<SpeedType::dimension_type> const velocityVerticalMag = velocity -
+ velocity.parallelProjectionOnto(magneticfield);
LengthType const gyroradius = p.GetEnergy() * velocityVerticalMag.GetNorm() * 1_V /
(corsika::units::constants::cSquared * abs(chargeNumber) *
- magneticfield.GetNorm() * 1_eV);;
+ magneticfield.GetNorm() * 1_eV);
+ std::cout << "Schrittlaenge " << velocity.norm() * min << " gyroradius " << gyroradius << std::endl;
LengthType B = 2 * gyroradius * asin ( lineWithB.ArcLength(0_s,min) / 2 / gyroradius);
+ std::cout << "Bogenlaenge" << B << std::endl;
histB(B / 1_m);
- histLB(L-B/1_m);
- histBS(B/1_m - lineWithB.ArcLength(0_s,min) / 1_m);
- histLBrel(L * 1_m/B-1);
+ histBlog(B / 1_m);
+ if(L-B/1_m > 0) {
+ histLB(L-B/1_m);
+ histBS(B/1_m - lineWithB.ArcLength(0_s,min) / 1_m);
+ histLBrel(L * 1_m/B-1);
+ } else {
+Bbig.push_back(B / 1_m);
+ }
+ } else {
+ histB(lineWithB.ArcLength(0_s,min) / 1_m);
+ histBlog(lineWithB.ArcLength(0_s,min) / 1_m);
}
int pdg = static_cast<int>(particles::GetPDG(p.GetPID()));
@@ -329,7 +468,6 @@ namespace corsika::process {
histLpi(L);
if (abs(pdg) == 2212 || abs(pdg) == 2112)
histLp(L);
- }
return std::make_tuple(geometry::Trajectory<geometry::Line>(line, min),
@@ -345,6 +483,9 @@ namespace corsika::process {
Particle const& p, corsika::geometry::Line line, corsika::units::si::LengthType Steplength, bool i) {
using namespace corsika::units::si;
+ if (line.GetV0().norm() * 1_s == 0_m)
+ return line;
+
//charge of the particle
int chargeNumber;
if (corsika::particles::IsNucleus(p.GetPID())) {
@@ -367,11 +508,9 @@ namespace corsika::process {
// Second Movement
position = position + direction * Steplength / 2;
- if(i) {
+ if(i == true) {
//if(chargeNumber != 0) {
- if(Steplength < 1000000_m) {
L = direction.norm() * Steplength / 2_m + Steplength / 2_m;
- }
//}
}
@@ -390,6 +529,10 @@ namespace corsika::process {
auto MagneticStep(Particle const& p, corsika::units::si::LengthType Steplength) {
using namespace corsika::units::si;
+ if (p.GetMomentum().norm() == 0_GeV) {
+ return std::make_tuple(p.GetPosition(), p.GetMomentum() / 1_GeV, double(0));
+ }
+
//charge of the particle
int chargeNumber;
if (corsika::particles::IsNucleus(p.GetPID())) {
@@ -413,8 +556,8 @@ namespace corsika::process {
direction = direction + direction.cross(magneticfield) * Steplength * k;
// Second Movement
position = position + direction * Steplength / 2;
- return std::make_tuple(position, direction.normalized());
-
+ auto L2 = direction.norm() * Steplength / 2_m + Steplength / 2_m;
+ return std::make_tuple(position, direction.normalized(), L2);
}
};
--
GitLab