diff --git a/Framework/Cascade/Cascade.h b/Framework/Cascade/Cascade.h
index f4b8e92d2c35be48f95d99b294f48d508756c79a..df9ee0e1820891ce8a59fe2817815cdd6992c1a0 100644
--- a/Framework/Cascade/Cascade.h
+++ b/Framework/Cascade/Cascade.h
@@ -36,6 +36,26 @@
#include <boost/type_index.hpp>
using boost::typeindex::type_id_with_cvr;
+#include <fstream>
+#include <boost/histogram.hpp>
+#include <boost/histogram/ostream.hpp>
+#include <corsika/process/tracking_line/dump_bh.hpp>
+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 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 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"));
+
/**
* The cascade namespace assembles all objects needed to simulate full particles cascades.
*/
@@ -107,6 +127,72 @@ 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.close();
+
+ /*std::cout << histLBrel << std::endl;
+ std::cout << histLSrel << std::endl;*/
+
+
+ std::ofstream file1("histL2.json");
+ dump_bh(file1, histL2);
+ file1.close();
+ std::ofstream file2("histS2.json");
+ dump_bh(file2, histS2);
+ file2.close();
+ std::ofstream file3("histB2.json");
+ dump_bh(file3, histB2);
+ file3.close();
+ /*std::ofstream file4("histLB.json");
+ dump_bh(file4, histLB);
+ file4.close();
+ std::ofstream file5("histLS.json");
+ dump_bh(file5, histLS);
+ file5.close();
+ std::ofstream file6("histBS.json");
+ dump_bh(file6, histBS);
+ file6.close();
+ std::ofstream file7("histLBrel.json");
+ dump_bh(file7, histLBrel);
+ file7.close();
+ std::ofstream file8("histLSrel.json");
+ dump_bh(file8, histLSrel);
+ file8.close();
+
+ std::ofstream file10("histELSrel.json");
+ dump_bh(file10, histELSrel);
+ file10.close();*/
+ std::ofstream file11("histLmu2.json");
+ dump_bh(file11, histLmu2);
+ file11.close();
+ std::ofstream file12("histLpi2.json");
+ dump_bh(file12, histLpi2);
+ file12.close();
+ std::ofstream file13("histLp2.json");
+ dump_bh(file13, histLp2);
+ file13.close();
+
+ };
+
+ /**
+ * set the nodes for all particles on the stack according to their numerical
+ * position
+ */
+ void SetNodes() {
+ std::for_each(fStack.begin(), fStack.end(), [&](auto& p) {
+ auto const* numericalNode =
+ fEnvironment.GetUniverse()->GetContainingNode(p.GetPosition());
+ p.SetNode(numericalNode);
+ });
+ }
+
/**
* The Run function is the main simulation loop, which processes
* particles from the Stack until the Stack is empty.
@@ -214,7 +300,7 @@ namespace corsika::cascade {
vParticle.GetEnergy() * units::constants::c;
// determine geometric tracking
- auto [lineWithoutB, stepWithoutB, stepWithB, geomMaxLength, nextVol] = fTracking.GetTrack(vParticle);
+ auto [stepWithoutB, stepWithB, geomMaxLength, magMaxLength, nextVol] = fTracking.GetTrack(vParticle);
[[maybe_unused]] auto const& dummy_nextVol = nextVol;
// convert next_step from grammage to length
@@ -228,20 +314,59 @@ namespace corsika::cascade {
// take minimum of geometry, interaction, decay for next step
auto min_distance = std::min(
- {distance_interact, distance_decay, distance_max, geomMaxLength});
-
- C8LOG_DEBUG("transport particle by : {} m", min_distance / 1_m);
-
- // determine displacement by the magnetic field
- auto [line, position, directionAfter] = fTracking.MagneticStep(vParticle, lineWithoutB, min_distance);
+ {distance_interact, distance_decay, distance_max, geomMaxLength, magMaxLength});
+
+ C8LOG_DEBUG("transport particle by : {} m "
+ "Max Displacement after: {} m "
+ "Medium transition after: {} m "
+ "Decay after: {} m "
+ "Interaction after: {} m",
+ min_distance/1_m, magMaxLength/1_m, geomMaxLength/1_m, distance_decay/1_m, distance_interact/1_m);
+
+ // determine displacement by the magnetic field
+ /*
+ int chargeNumber;
+ if (corsika::particles::IsNucleus(vParticle.GetPID())) {
+ chargeNumber = vParticle.GetNuclearZ();
+ } else {
+ chargeNumber = corsika::particles::GetChargeNumber(vParticle.GetPID());
+ }
+ auto const* currentLogicalVolumeNode = vParticle.GetNode();
+ auto magneticfield = currentLogicalVolumeNode->GetModelProperties().GetMagneticField(vParticle.GetPosition());
+ auto k = chargeNumber * corsika::units::constants::cSquared * 1_eV /
+ ((vParticle.GetMomentum() / vParticle.GetEnergy() *
+ corsika::units::constants::c).norm() * vParticle.GetEnergy() * 1_V);
+ geometry::Vector<dimensionless_d> const directionBefore = vParticle.GetMomentum().normalized();
+ LengthType Steplength = min_distance;
+ if (chargeNumber != 0) {
+ Steplength = (-sqrt(2 * k * min_distance * (directionBefore.cross(magneticfield)).norm() + 1) - 1)
+ / ((directionBefore.cross(magneticfield)).norm() * k);
+ std::cout << "Steplength2 " << Steplength << std::endl;
+
+ Steplength = (sqrt(2 * k * min_distance * (directionBefore.cross(magneticfield)).norm() + 1) - 1)
+ / ((directionBefore.cross(magneticfield)).norm() * k);
+ std::cout << "Steplength1 " << Steplength << std::endl;
+ // not totally sure if it is always the positive solution
+ }
+ */
+ // 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 distance = position - vParticle.GetPosition();
- // distance.norm() != min_distance if q != 0
- // small error can be neglected
+
+ //Building Trajectory for Continuous processes
+ geometry::Vector<SpeedType::dimension_type> velocity =
+ vParticle.GetMomentum() / vParticle.GetEnergy() * corsika::units::constants::c;
+ if (distance.norm() != 0_m) {
+ velocity = distance.normalized() * velocity.norm();
+ }
+ geometry::Line line(vParticle.GetPosition(), velocity);
+ geometry::Trajectory<geometry::Line> stepNew(line, distance.norm() / line.GetV0().norm());
// here the particle is actually moved along the trajectory to new position:
// std::visit(setup::ParticleUpdate<Particle>{vParticle}, step);
- vParticle.SetMomentum(directionAfter.normalized() * vParticle.GetMomentum().norm());
- geometry::Trajectory<geometry::Line> stepNew(line, distance.norm() / line.GetV0().norm());
+ vParticle.SetMomentum(direction * vParticle.GetMomentum().norm());
vParticle.SetPosition(position);
vParticle.SetTime(vParticle.GetTime() + distance.norm() / units::constants::c);
std::cout << "New Position: " << vParticle.GetPosition().GetCoordinates() << std::endl;
@@ -268,7 +393,7 @@ namespace corsika::cascade {
TStackView secondaries(vParticle);
- if (min_distance != distance_max) {
+ if (min_distance != distance_max && min_distance != magMaxLength) {
/*
Create SecondaryView object on Stack. The data container
remains untouched and identical, and 'projectil' is identical
diff --git a/Processes/TrackWriter/TrackWriter.cc b/Processes/TrackWriter/TrackWriter.cc
index a3d3edbf39eb27f4aaaf66c6255f6cb879572a9a..1bef877c3a48d327e56de4a70b6f13df5154f117 100644
--- a/Processes/TrackWriter/TrackWriter.cc
+++ b/Processes/TrackWriter/TrackWriter.cc
@@ -28,7 +28,7 @@ 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 "s
+ fFile << "# PID, E / eV, start coordinates / m, displacement vector to end / m, steplength / m "s
<< '\n';
}
@@ -47,7 +47,9 @@ namespace corsika::process::track_writer {
<< std::setw(width) << std::scientific << std::setprecision(precision) << start[2] / 1_m
<< std::setw(width) << std::scientific << std::setprecision(precision) << delta[0] / 1_m
<< std::setw(width) << std::scientific << std::setprecision(precision) << delta[1] / 1_m
- << std::setw(width) << std::scientific << std::setprecision(precision) << delta[2] / 1_m << '\n';
+ << std::setw(width) << std::scientific << std::setprecision(precision) << delta[2] / 1_m
+ << std::setw(width) << std::scientific << std::setprecision(precision) << delta.norm() / 1_m
+ << '\n';
// clang-format on
return process::EProcessReturn::eOk;
diff --git a/Processes/TrackingLine/CMakeLists.txt b/Processes/TrackingLine/CMakeLists.txt
index dd626f8778c5a4139fe042fc651bb5e20eb315fe..36d53e30c3b1dceccc83e079a52499d5e73734fa 100644
--- a/Processes/TrackingLine/CMakeLists.txt
+++ b/Processes/TrackingLine/CMakeLists.txt
@@ -1,6 +1,7 @@
set (
MODEL_HEADERS
TrackingLine.h
+ dump_bh.hpp
)
set (
diff --git a/Processes/TrackingLine/TrackingLine.h b/Processes/TrackingLine/TrackingLine.h
index d68a96d65c80c9f308c7a4694dc9f6faf8a60fda..875569079fa4a65c86c0dd5041eeb0b7f0e2f7a9 100644
--- a/Processes/TrackingLine/TrackingLine.h
+++ b/Processes/TrackingLine/TrackingLine.h
@@ -20,6 +20,12 @@
#include <type_traits>
#include <utility>
+#include <fstream>
+#include <iostream>
+#include <boost/histogram.hpp>
+#include <boost/histogram/ostream.hpp>
+#include <corsika/process/tracking_line/dump_bh.hpp>
+
namespace corsika::environment {
template <typename IEnvironmentModel>
class Environment;
@@ -27,6 +33,23 @@ namespace corsika::environment {
class VolumeTreeNode;
} // 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 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 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;
+
namespace corsika::process {
namespace tracking_line {
@@ -36,11 +59,64 @@ namespace corsika::process {
corsika::units::si::TimeType TimeOfIntersection(geometry::Line const&,
geometry::Plane const&);
-
+
class TrackingLine {
public:
TrackingLine() = default;
+ ~TrackingLine(){
+ std::ofstream myfile;
+ myfile.open ("histograms.txt");
+ 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 file1("histL.json");
+ dump_bh(file1, histL);
+ file1.close();
+ std::ofstream file2("histS.json");
+ dump_bh(file2, histS);
+ file2.close();
+ std::ofstream file3("histB.json");
+ dump_bh(file3, histB);
+ file3.close();
+ std::ofstream file4("histLB.json");
+ dump_bh(file4, histLB);
+ file4.close();
+ std::ofstream file5("histLS.json");
+ dump_bh(file5, histLS);
+ file5.close();
+ std::ofstream file6("histBS.json");
+ dump_bh(file6, histBS);
+ file6.close();
+ std::ofstream file7("histLBrel.json");
+ dump_bh(file7, histLBrel);
+ file7.close();
+ std::ofstream file8("histLSrel.json");
+ dump_bh(file8, histLSrel);
+ file8.close();
+
+ std::ofstream file10("histELSrel.json");
+ dump_bh(file10, histELSrel);
+ file10.close();
+ std::ofstream file11("histLmu.json");
+ dump_bh(file11, histLmu);
+ file11.close();
+ std::ofstream file12("histLpi.json");
+ dump_bh(file12, histLpi);
+ file12.close();
+ std::ofstream file13("histLp.json");
+ dump_bh(file13, histLp);
+ file13.close();
+
+ };
template <typename Particle> // was Stack previously, and argument was
// Stack::StackIterator
@@ -88,7 +164,10 @@ namespace corsika::process {
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 Steplength = 10_m; // length irrelevant if q=0 and else it gets changed again
+ LengthType Steplength1 = 10_m; // length irrelevant if q=0 and else it gets changed again
+ LengthType Steplength2 = 10_m;
+ LengthType Steplimit = 1_m / 0;
+ //infinite kann man auch anders schreiben
// for entering from outside
auto addIfIntersects = [&](auto const& vtn) {
@@ -98,7 +177,7 @@ namespace corsika::process {
// everything is a sphere, crashes with exception if not
// creating Line with magnetic field
- if (chargeNumber != 0) {
+ if (chargeNumber != 0) {
// 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);
@@ -116,18 +195,24 @@ namespace corsika::process {
}
}
if (tmp.size() > 0) {
- Steplength = 1_m * *std::min_element(tmp.begin(),tmp.end());
- std::cout << "Steplength to next volume = " << Steplength << std::endl;
+ Steplength1 = 1_m * *std::min_element(tmp.begin(),tmp.end());
+ std::cout << "Steplength to next volume = " << Steplength1 << std::endl;
+ 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)
}
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, position, direction] = MagneticStep(p, line, Steplength);
- // new particle position and direction not needed in this case
+ 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
@@ -168,15 +253,15 @@ namespace corsika::process {
std::cout << "Solutions for current Volume: " << solutions[i].real() << std::endl;
}
}
- LengthType Steplength;
if (tmp.size() > 0) {
- Steplength = 1_m * *std::min_element(tmp.begin(),tmp.end());
- if (numericallyInside == false) {
- int p = std::min_element(tmp.begin(),tmp.end()) - tmp.begin();
- tmp.erase(tmp.begin() + p);
- Steplength = 1_m * *std::min_element(tmp.begin(),tmp.end());
- }
- std::cout << "steplength out of current volume = " << Steplength << std::endl;
+ Steplength2 = 1_m * *std::min_element(tmp.begin(),tmp.end());
+ if (numericallyInside == false) {
+ int p = std::min_element(tmp.begin(),tmp.end()) - tmp.begin();
+ tmp.erase(tmp.begin() + p);
+ Steplength2 = 1_m * *std::min_element(tmp.begin(),tmp.end());
+ }
+ std::cout << "steplength out of current volume = " << Steplength2 << std::endl;
+ std::cout << "Time out of current volume = " << Steplength2 / velocity.norm() << std::endl;
} else {
std::cout << "no intersection (2)!" << std::endl;
// what to do when this happens? (very unlikely)
@@ -184,8 +269,7 @@ namespace corsika::process {
delete [] solutions;
}
- auto [line2, position, direction] = MagneticStep(p, line, Steplength);
- // new particle position and direction not needed in this case
+ auto line2 = MagneticStep(p, line, Steplength2, false);
// instead of changing the line with magnetic field, the TimeOfIntersection() could be changed
// using line has some errors for huge steps
@@ -212,22 +296,55 @@ namespace corsika::process {
<< min
// << " " << minIter->second->GetModelProperties().GetName()
<< std::endl;
-
- auto [lineWithB, position, direction] = MagneticStep(p, line, velocity.norm() * min);
- // new particle position and direction not needed in this case
+
+ auto lineWithB = MagneticStep(p, line, velocity.norm() * min, true);
+
+ if(min * velocity.norm() < 1000000_m) {
+ histL(L);
+ histS(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);
+
+
+ LengthType const gyroradius = p.GetEnergy() * velocityVerticalMag.GetNorm() * 1_V /
+ (corsika::units::constants::cSquared * abs(chargeNumber) *
+ magneticfield.GetNorm() * 1_eV);;
+ LengthType B = 2 * gyroradius * asin ( lineWithB.ArcLength(0_s,min) / 2 / gyroradius);
+ 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);
+ }
+
+ 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(line, geometry::Trajectory<geometry::Line>(line, min),
+ return std::make_tuple(geometry::Trajectory<geometry::Line>(line, min),
geometry::Trajectory<geometry::Line>(lineWithB, min),
- velocity.norm() * min, minIter->second);
+ velocity.norm() * min, Steplimit, minIter->second);
}
template <typename Particle> // was Stack previously, and argument was
// Stack::StackIterator
// determine direction of the particle after adding magnetic field
- auto MagneticStep(Particle const& p, corsika::geometry::Line line, corsika::units::si::LengthType Steplength) {
+ corsika::geometry::Line MagneticStep(
+ Particle const& p, corsika::geometry::Line line, corsika::units::si::LengthType Steplength, bool i) {
using namespace corsika::units::si;
-
+
//charge of the particle
int chargeNumber;
if (corsika::particles::IsNucleus(p.GetPID())) {
@@ -235,26 +352,69 @@ namespace corsika::process {
} else {
chargeNumber = corsika::particles::GetChargeNumber(p.GetPID());
}
+
auto const* currentLogicalVolumeNode = p.GetNode();
auto magneticfield = currentLogicalVolumeNode->GetModelProperties().GetMagneticField(p.GetPosition());
auto k = chargeNumber * corsika::units::constants::cSquared * 1_eV / (line.GetV0().norm() * p.GetEnergy() * 1_V);
- geometry::Vector<dimensionless_d> const directionBefore = line.GetV0().normalized();
-
+ geometry::Vector<dimensionless_d> direction = line.GetV0().normalized();
+ auto position = p.GetPosition();
+
// First Movement
// assuming magnetic field does not change during movement
- auto position = p.GetPosition() + directionBefore * Steplength / 2;
+ position = position + direction * Steplength / 2;
// Change of direction by magnetic field
- geometry::Vector<dimensionless_d> const directionAfter = directionBefore + directionBefore.cross(magneticfield) *
- Steplength * k;
+ direction = direction + direction.cross(magneticfield) * Steplength * k;
// Second Movement
- position = position + directionAfter * Steplength / 2;
- auto distance = position - p.GetPosition();
+ position = position + direction * Steplength / 2;
+ if(i) {
+ //if(chargeNumber != 0) {
+ if(Steplength < 1000000_m) {
+ L = direction.norm() * Steplength / 2_m + Steplength / 2_m;
+ }
+ //}
+ }
+
+ geometry::Vector<LengthType::dimension_type> const distance = position - p.GetPosition();
if(distance.norm() == 0_m) {
- return std::make_tuple(line, position, directionAfter);
+ return line;
}
geometry::Line newLine(p.GetPosition(), distance.normalized() * line.GetV0().norm());
- return std::make_tuple(newLine, position, directionAfter);
+ return newLine;
+ }
+
+ template <typename Particle> // was Stack previously, and argument was
+ // Stack::StackIterator
+
+ // determine direction of the particle after adding magnetic field
+ auto MagneticStep(Particle const& p, corsika::units::si::LengthType Steplength) {
+ using namespace corsika::units::si;
+
+ //charge of the particle
+ int chargeNumber;
+ if (corsika::particles::IsNucleus(p.GetPID())) {
+ chargeNumber = p.GetNuclearZ();
+ } else {
+ chargeNumber = corsika::particles::GetChargeNumber(p.GetPID());
+ }
+
+ auto const* currentLogicalVolumeNode = p.GetNode();
+ auto magneticfield = currentLogicalVolumeNode->GetModelProperties().GetMagneticField(p.GetPosition());
+ geometry::Vector<SpeedType::dimension_type> velocity =
+ p.GetMomentum() / p.GetEnergy() * corsika::units::constants::c;
+ auto k = chargeNumber * corsika::units::constants::cSquared * 1_eV / (velocity.norm() * p.GetEnergy() * 1_V);
+ geometry::Vector<dimensionless_d> direction = velocity.normalized();
+ auto position = p.GetPosition();
+
+ // First Movement
+ // assuming magnetic field does not change during movement
+ position = position + direction * Steplength / 2;
+ // Change of direction by magnetic field
+ direction = direction + direction.cross(magneticfield) * Steplength * k;
+ // Second Movement
+ position = position + direction * Steplength / 2;
+ return std::make_tuple(position, direction.normalized());
+
}
};
diff --git a/Processes/TrackingLine/dump_bh.hpp b/Processes/TrackingLine/dump_bh.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..01e0dac20726b2d2b4787fad3bb3a9f0fff27b09
--- /dev/null
+++ b/Processes/TrackingLine/dump_bh.hpp
@@ -0,0 +1,186 @@
+#pragma once
+
+#include <boost/histogram.hpp>
+#include <sstream>
+#include <string>
+#include <vector>
+
+/*
+ To save one 1D boost::histogram in C++, include this file and run:
+
+ dump_bh_1d("file_name.json", hist);
+
+ Note that the file_name.json will be overwritten, so use a new file for each histogram! or use:
+
+ std::ofstream file1("test_hist.json");
+ dump_bh(file1, h1);
+ file << ",\n";
+ dump_bh(file1, h2);
+ file << ",\n";
+ dump_bh(file1, h3);
+ file1.close();
+
+
+ In python, just read the json and plot this histogram with matplotib
+
+ for 1D histogram:
+
+ import json
+ import matplotlib.pyplot as plt
+
+ h1=json.load(open("test_hist_direct.json"))
+ print (h1)
+ plt.hist(h1['bins'][:-1], bins=h1['bins'], weights=h1['data'])
+ plt.show()
+
+
+ for 2D histogram (needs some list/array gymnastics):
+
+ import json
+ import matplotlib.pyplot as plt
+ import numpy as np
+
+ h2=json.load(open("test_hist_2D_direct.json"))
+ print (h2)
+ xx,yy = np.meshgrid(h2['xbins'][:-1], h2['ybins'][:-1])
+ plt.hist2d([i for s in xx for i in s], [i for s in yy for i in s],
+ bins=[h2['xbins'],h2['ybins']], weights=h2['data'])
+ plt.show()
+
+
+ Do not forget to add axis titles, legend, etc.
+ */
+
+
+template<typename T>
+void dump_bh_2d(std::ostream& os, T& h)
+{
+ // determine number of axes (rank) and axes sizes
+ const int rank = 2;
+ std::vector<int> axis_size(rank);
+ for (unsigned int i=0; i<rank; ++i)
+ axis_size[i] = h.axis(i).size();
+
+ // start dump json object
+ os << "{\n";
+ os << " \"rank\" : " << rank << ",\n";
+ os << " \"size\" : [";
+ bool first = true;
+ for (auto s : axis_size) {
+ os << (first?"":", ") << s;
+ first = false;
+ }
+ os << "],\n";
+
+ // bins-x
+ os << " \"xbins\" : [";
+ double lastx = 0;
+ for (unsigned int i=0; i<h.axis(0).size(); ++i) {
+ os << h.axis(0).bin(i).lower() << ", ";
+ lastx = h.axis(0).bin(i).upper();
+ }
+ os << lastx << "],\n";
+
+ // bins-y
+ os << " \"ybins\" : [";
+ double lasty = 0;
+ for (unsigned int i=0; i<h.axis(1).size(); ++i) {
+ os << h.axis(1).bin(i).lower() << ", ";
+ lasty = h.axis(1).bin(i).upper();
+ }
+ os << lasty << "],\n";
+
+ // binned data
+ os << " \"data\" : [";
+ first = true;
+ for (auto x : boost::histogram::indexed(h)) {
+ os << (first?"":", ") << *x;
+ first = false;
+ }
+ os << " ]\n";
+ os << "}\n";
+}
+
+
+
+
+template<typename T>
+void dump_bh_1d(std::ostream& os, T& h)
+{
+ const int rank = 1;
+ // determine number of axes (rank) and axes sizes
+ std::vector<int> axis_size(rank);
+ for (unsigned int i=0; i<rank; ++i)
+ axis_size[i] = h.axis(i).size();
+
+ // start dump json object
+ os << "{\n";
+ os << " \"rank\" : 1,\n";
+ os << " \"size\" : [";
+ bool first = true;
+ for (auto s : axis_size) {
+ os << (first?"":", ") << s << "],\n";
+ first = false;
+ }
+
+ // underflow data
+ os << " \"underflow\" : [";
+ os << h.at(boost::histogram::axis::option::underflow);
+ os << "],\n";
+
+ // overflow data
+ os << " \"overflow\" : [";
+ os << h.at(boost::histogram::axis::option::overflow);
+ os << "],\n";
+
+ // bins
+ os << " \"bins\" : [";
+ double last = 0;
+ for (auto x : boost::histogram::indexed(h)) {
+ os << x.bin().lower() << ", ";
+ last = x.bin().upper();
+ }
+ os << last << "],\n";
+
+ // data of bins
+ os << " \"data\" : [";
+ first = true;
+ for (auto x : boost::histogram::indexed(h)) {
+ os << (first?"":", ") << *x;
+ first = false;
+ }
+ os << " ]\n";
+ os << "}\n";
+}
+
+
+
+
+template<typename T>
+void dump_bh(std::ostream& os, T& h)
+{
+ // determine number of axes (rank) and axes sizes
+ const int rank = h.rank();
+ switch (rank) {
+ case 1:
+ dump_bh_1d(os, h);
+ break;
+ case 2:
+ dump_bh_2d(os, h);
+ break;
+ default:
+ {
+ throw std::runtime_error("multi dim hist not implemented");
+ }
+ break;
+ }
+}
+
+template<typename T>
+void dump_bh(const char* fname, T& h)
+{
+ std::ofstream file(fname);
+ dump_bh(file, h);
+ file.close();
+}
+