diff --git a/Documentation/Examples/CMakeLists.txt b/Documentation/Examples/CMakeLists.txt
index 5a5bd65531b948d1b71536e83f7f88a905c6f4fa..cd7e5251321adc682422d25ecfe992a2d0b63b40 100644
--- a/Documentation/Examples/CMakeLists.txt
+++ b/Documentation/Examples/CMakeLists.txt
@@ -28,6 +28,7 @@ target_link_libraries (cascade_example SuperStupidStack CORSIKAunits CORSIKAlogg
ProcessSibyll
CORSIKAcascade
ProcessStackInspector
+ ProcessTrackWriter
CORSIKAprocesses
CORSIKAparticles
CORSIKAgeometry
diff --git a/Documentation/Examples/cascade_example.cc b/Documentation/Examples/cascade_example.cc
index e8377f80c39576f17323dd284f67fac4be6acf5e..f3ebb57d541bf585312f47083d92b745bf82fc36 100644
--- a/Documentation/Examples/cascade_example.cc
+++ b/Documentation/Examples/cascade_example.cc
@@ -25,6 +25,8 @@
#include <corsika/process/sibyll/Decay.h>
#include <corsika/process/sibyll/Interaction.h>
+#include <corsika/process/track_writer/TrackWriter.h>
+
#include <corsika/units/PhysicalUnits.h>
#include <corsika/random/RNGManager.h>
@@ -192,7 +194,6 @@ public:
// The example main program for a particle cascade
//
int main() {
-
// initialize random number sequence(s)
corsika::random::RNGManager::GetInstance().RegisterRandomStream("cascade");
@@ -222,9 +223,10 @@ int main() {
corsika::process::sibyll::Interaction sibyll;
corsika::process::sibyll::Decay decay;
ProcessCut cut(8_GeV);
+ corsika::process::TrackWriter::TrackWriter trackWriter("tracks.dat");
// assemble all processes into an ordered process list
- auto sequence = p0 << sibyll << decay << cut;
+ auto sequence = p0 << sibyll << decay << cut << trackWriter;
// cout << "decltype(sequence)=" << type_id_with_cvr<decltype(sequence)>().pretty_name()
// << "\n";
@@ -232,7 +234,7 @@ int main() {
// setup particle stack, and add primary particle
setup::Stack stack;
stack.Clear();
- const hep::EnergyType E0 = 1_TeV;
+ const hep::EnergyType E0 = 10_TeV;
{
auto particle = stack.NewParticle();
particle.SetPID(Code::Proton);
@@ -241,7 +243,7 @@ int main() {
particle.SetEnergy(E0);
particle.SetMomentum(plab);
particle.SetTime(0_ns);
- Point p(rootCS, 0_m, 0_m, 10_km);
+ Point p(rootCS, 0_m, 0_m, 0_m);
particle.SetPosition(p);
cout << particle.GetEnergy() / 1_GeV << endl;
}
diff --git a/Framework/Cascade/Cascade.h b/Framework/Cascade/Cascade.h
index 973c614391299349c0081e50fbb4b19831d1b69f..6aed37d644b868c1ed92d5ce87475db45245225f 100644
--- a/Framework/Cascade/Cascade.h
+++ b/Framework/Cascade/Cascade.h
@@ -57,13 +57,7 @@ namespace corsika::cascade {
}
void Step(Particle& particle) {
-
using namespace corsika::units::si;
- using std::cout;
- using std::endl;
- using std::log;
-
- cout << particle.GetEnergy() / 1_GeV << endl;
// determine geometric tracking
corsika::setup::Trajectory step = fTracking.GetTrack(particle);
@@ -80,11 +74,15 @@ namespace corsika::cascade {
<< ", next_interact=" << next_interact << std::endl;
// convert next_step from grammage to length
+ auto const* currentNode =
+ fEnvironment.GetUniverse()->GetContainingNode(particle.GetPosition());
+
+ if (currentNode == &*fEnvironment.GetUniverse()) {
+ throw std::runtime_error("particle entered void universe");
+ }
+
LengthType const distance_interact =
- fEnvironment.GetUniverse()
- ->GetContainingNode(particle.GetPosition())
- ->GetModelProperties()
- .ArclengthFromGrammage(step, next_interact);
+ currentNode->GetModelProperties().ArclengthFromGrammage(step, next_interact);
// determine the maximum geometric step length
LengthType const distance_max = fProcessSequence.MaxStepLength(particle, step);
@@ -116,6 +114,8 @@ namespace corsika::cascade {
particle.SetPosition(step.PositionFromArclength(min_distance));
// .... also update time, momentum, direction, ...
+ step.LimitEndTo(min_distance);
+
// apply all continuous processes on particle + track
corsika::process::EProcessReturn status =
fProcessSequence.DoContinuous(particle, step, fStack);
@@ -127,11 +127,11 @@ namespace corsika::cascade {
return;
}
- std::cout << "sth. happening before geometric limit ?"
+ std::cout << "sth. happening before geometric limit ? "
<< ((min_distance < distance_max) ? "yes" : "no") << std::endl;
if (min_distance < distance_max) { // interaction to happen within geometric limit
- // check weather decay or interaction limits this step
+ // check whether decay or interaction limits this step
if (min_distance == distance_interact) {
std::cout << "collide" << std::endl;
diff --git a/Framework/Geometry/Trajectory.h b/Framework/Geometry/Trajectory.h
index 5a9d1ed56a1e1232bb0591f2646a4f4dad1f6f85..bcaa45fba325367048a2a7b93c6ddb0752428300 100644
--- a/Framework/Geometry/Trajectory.h
+++ b/Framework/Geometry/Trajectory.h
@@ -43,6 +43,10 @@ namespace corsika::geometry {
assert(t >= 0 * corsika::units::si::second);
return T::ArcLength(0, t);
}
+
+ void LimitEndTo(corsika::units::si::LengthType limit) {
+ fTimeLength = T::TimeFromArclength(limit);
+ }
};
} // namespace corsika::geometry
diff --git a/Processes/CMakeLists.txt b/Processes/CMakeLists.txt
index 41a8277d9c736f387d68b255ff2af7415133c10c..d6b90e70ec3cc7d1efb5266626d41f50624cb019 100644
--- a/Processes/CMakeLists.txt
+++ b/Processes/CMakeLists.txt
@@ -2,6 +2,7 @@
add_subdirectory (NullModel)
add_subdirectory (Sibyll)
add_subdirectory (StackInspector)
+add_subdirectory (TrackWriter)
add_subdirectory (TrackingLine)
#add_custom_target(CORSIKAprocesses)
diff --git a/Processes/Sibyll/Interaction.h b/Processes/Sibyll/Interaction.h
index 101fbdd7b53aae2c4acbb1ebc071796449a6d299..2b03b26009c1d6f0e4f9c9be75d9452d2715124f 100644
--- a/Processes/Sibyll/Interaction.h
+++ b/Processes/Sibyll/Interaction.h
@@ -75,8 +75,9 @@ namespace corsika::process::sibyll {
// FOR NOW: assume target is oxygen
const int kTarget = corsika::particles::Oxygen::GetNucleusA();
- hep::EnergyType Etot =
- p.GetEnergy() + kTarget * corsika::particles::Proton::GetMass();
+ const hep::MassType nucleon_mass = 0.5 * (corsika::particles::Proton::GetMass() +
+ corsika::particles::Neutron::GetMass());
+ hep::EnergyType Etot = p.GetEnergy() + nucleon_mass;
MomentumVector Ptot(rootCS, {0.0_GeV, 0.0_GeV, 0.0_GeV});
// FOR NOW: assume target is at rest
MomentumVector pTarget(rootCS, {0.0_GeV, 0.0_GeV, 0.0_GeV});
@@ -143,6 +144,18 @@ namespace corsika::process::sibyll {
Point pOrig = p.GetPosition();
TimeType tOrig = p.GetTime();
+ // sibyll CS has z along particle momentum
+ // FOR NOW: hard coded z-axis for corsika frame
+ QuantityVector<length_d> const zAxis{0_m, 0_m, 1_m};
+ QuantityVector<length_d> const xAxis{1_m, 0_m, 0_m};
+ [[maybe_unused]] auto pt = [](MomentumVector p) {
+ return sqrt(p.GetComponents()[0] * p.GetComponents()[0] +
+ p.GetComponents()[1] * p.GetComponents()[1]);
+ };
+ double theta = acos(p.GetMomentum().GetComponents()[2] / p.GetMomentum().norm());
+ cout << "ProcessSibyll: polar angle between sibyllCS and rootCS: " << theta
+ << endl;
+ CoordinateSystem sibyllCS = rootCS.rotate(xAxis, theta);
/*
the target should be defined by the Environment,
@@ -157,11 +170,16 @@ namespace corsika::process::sibyll {
GetNucleusA(); // env.GetTargetParticle().GetPID();
// FOR NOW: target is always at rest
- const EnergyType Etarget = 0_GeV + corsika::particles::Proton::GetMass();
+ const hep::MassType nucleon_mass = 0.5 * (corsika::particles::Proton::GetMass() +
+ corsika::particles::Neutron::GetMass());
+ const EnergyType Etarget = 0_GeV + nucleon_mass;
const auto pTarget = MomentumVector(rootCS, 0_GeV, 0_GeV, 0_GeV);
cout << "target momentum (GeV/c): " << pTarget.GetComponents() / 1_GeV << endl;
cout << "beam momentum (GeV/c): " << p.GetMomentum().GetComponents() / 1_GeV
<< endl;
+ cout << "beam momentum in sibyll frame (GeV/c): "
+ << p.GetMomentum().GetComponents(sibyllCS) / 1_GeV << endl;
+
cout << "position of interaction: " << pOrig.GetCoordinates() << endl;
cout << "time: " << tOrig << endl;
@@ -220,23 +238,31 @@ namespace corsika::process::sibyll {
// add particles from sibyll to stack
// link to sibyll stack
+ // here we need to pass projectile momentum and energy to define the local
+ // sibyll frame and the boosts to the lab. frame
SibStack ss;
- // SibStack does not know about momentum yet so we need counter to access
// momentum array in Sibyll
- int i = -1;
MomentumVector Plab_final(rootCS, {0.0_GeV, 0.0_GeV, 0.0_GeV});
EnergyType E_final = 0_GeV, Ecm_final = 0_GeV;
for (auto& psib : ss) {
- ++i;
+
// skip particles that have decayed in Sibyll
- if (abs(s_plist_.llist[i]) > 100) continue;
+ if (psib.HasDecayed()) continue;
// transform energy to lab. frame, primitve
// compute beta_vec * p_vec
// arbitrary Lorentz transformation based on sibyll routines
const auto gammaBetaComponents = gambet.GetComponents();
- const auto pSibyllComponents = psib.GetMomentum().GetComponents();
+ // FOR NOW: fill vector in sibCS and then rotate into rootCS
+ // can be done in SibStack by passing sibCS
+ // get momentum vector in sibyllCS
+ const auto pSibyllComponentsSibCS = psib.GetMomentum().GetComponents();
+ // temporary vector in sibyllCS
+ auto SibVector = MomentumVector(sibyllCS, pSibyllComponentsSibCS);
+ // rotatate to rootCS
+ const auto pSibyllComponents = SibVector.GetComponents(rootCS);
+ // boost to lab. frame
hep::EnergyType en_lab = 0. * 1_GeV;
hep::MomentumType p_lab_components[3];
en_lab = psib.GetEnergy() * gamma;
@@ -255,7 +281,6 @@ namespace corsika::process::sibyll {
auto pnew = s.NewParticle();
pnew.SetEnergy(en_lab);
pnew.SetPID(process::sibyll::ConvertFromSibyll(psib.GetPID()));
-
corsika::geometry::QuantityVector<energy_hep_d> p_lab_c{
p_lab_components[0], p_lab_components[1], p_lab_components[2]};
pnew.SetMomentum(MomentumVector(rootCS, p_lab_c));
diff --git a/Processes/Sibyll/SibStack.h b/Processes/Sibyll/SibStack.h
index 1e44931b2cc79ea64233b3643f5ef64163d33a83..3364e7deac87ca2fda683d188aeaf3495eaa1fb4 100644
--- a/Processes/Sibyll/SibStack.h
+++ b/Processes/Sibyll/SibStack.h
@@ -31,7 +31,7 @@ namespace corsika::process::sibyll {
void Clear() { s_plist_.np = 0; }
int GetSize() const { return s_plist_.np; }
-#warning check actual capacity of sibyll stack
+
int GetCapacity() const { return 8000; }
void SetId(const int i, const int v) { s_plist_.llist[i] = v; }
@@ -91,6 +91,9 @@ namespace corsika::process::sibyll {
corsika::units::hep::EnergyType GetEnergy() const {
return GetStackData().GetEnergy(GetIndex());
}
+ bool HasDecayed() const {
+ return abs(GetStackData().GetId(GetIndex())) > 100 ? true : false;
+ }
void SetPID(const int v) { GetStackData().SetId(GetIndex(), v); }
corsika::process::sibyll::SibyllCode GetPID() const {
return static_cast<corsika::process::sibyll::SibyllCode>(
diff --git a/Processes/TrackWriter/CMakeLists.txt b/Processes/TrackWriter/CMakeLists.txt
new file mode 100644
index 0000000000000000000000000000000000000000..e9c73d68b7f0f072e07c0e007b0dcd993adbf0ca
--- /dev/null
+++ b/Processes/TrackWriter/CMakeLists.txt
@@ -0,0 +1,65 @@
+
+set (
+ MODEL_SOURCES
+ TrackWriter.cc
+ )
+
+set (
+ MODEL_HEADERS
+ TrackWriter.h
+ )
+
+set (
+ MODEL_NAMESPACE
+ corsika/process/track_writer
+ )
+
+add_library (ProcessTrackWriter STATIC ${MODEL_SOURCES})
+CORSIKA_COPY_HEADERS_TO_NAMESPACE (ProcessTrackWriter ${MODEL_NAMESPACE} ${MODEL_HEADERS})
+
+set_target_properties (
+ ProcessTrackWriter
+ PROPERTIES
+ VERSION ${PROJECT_VERSION}
+ SOVERSION 1
+# PUBLIC_HEADER "${MODEL_HEADERS}"
+ )
+
+# target dependencies on other libraries (also the header onlys)
+target_link_libraries (
+ ProcessTrackWriter
+ CORSIKAunits
+ CORSIKAparticles
+ CORSIKAgeometry
+ CORSIKAsetup
+ )
+
+target_include_directories (
+ ProcessTrackWriter
+ INTERFACE
+ $<BUILD_INTERFACE:${PROJECT_BINARY_DIR}/include>
+ $<INSTALL_INTERFACE:include/include>
+ )
+
+install (
+ TARGETS ProcessTrackWriter
+ LIBRARY DESTINATION lib
+ ARCHIVE DESTINATION lib
+# PUBLIC_HEADER DESTINATION include/${MODEL_NAMESPACE}
+ )
+
+
+# --------------------
+# code unit testing
+#add_executable (testNullModel testNullModel.cc)
+
+#target_link_libraries (
+# testNullModel ProcessNullModel
+# CORSIKAsetup
+# CORSIKAgeometry
+# CORSIKAunits
+# CORSIKAthirdparty # for catch2
+# )
+
+#add_test (NAME testNullModel COMMAND testNullModel)
+
diff --git a/Processes/TrackWriter/TrackWriter.cc b/Processes/TrackWriter/TrackWriter.cc
new file mode 100644
index 0000000000000000000000000000000000000000..ffc67cde8a478486cca12670a7e0bb7838d130fa
--- /dev/null
+++ b/Processes/TrackWriter/TrackWriter.cc
@@ -0,0 +1,20 @@
+/**
+ * (c) Copyright 2018 CORSIKA Project, corsika-project@lists.kit.edu
+ *
+ * See file AUTHORS for a list of contributors.
+ *
+ * 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 <corsika/process/track_writer/TrackWriter.h>
+#include <string>
+
+void corsika::process::TrackWriter::TrackWriter::Init() {
+ using namespace std::string_literals;
+
+ fFile.open(fFilename);
+ fFile << "# PID, E / eV, start coordinates / m, displacement vector to end / m "s
+ << '\n';
+}
diff --git a/Processes/TrackWriter/TrackWriter.h b/Processes/TrackWriter/TrackWriter.h
new file mode 100644
index 0000000000000000000000000000000000000000..10f2d574f355f72f46e9d681adcd853dac436c4c
--- /dev/null
+++ b/Processes/TrackWriter/TrackWriter.h
@@ -0,0 +1,58 @@
+/**
+ * (c) Copyright 2018 CORSIKA Project, corsika-project@lists.kit.edu
+ *
+ * See file AUTHORS for a list of contributors.
+ *
+ * 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.
+ */
+
+#ifndef _Processes_TrackWriter_h_
+#define _Processes_TrackWriter_h_
+
+#include <corsika/particles/ParticleProperties.h>
+#include <corsika/process/ContinuousProcess.h>
+#include <corsika/setup/SetupTrajectory.h>
+#include <corsika/units/PhysicalUnits.h>
+#include <fstream>
+#include <limits>
+#include <string>
+
+namespace corsika::process::TrackWriter {
+
+ class TrackWriter : public corsika::process::ContinuousProcess<TrackWriter> {
+
+ public:
+ TrackWriter(std::string const& filename)
+ : fFilename(filename) {}
+
+ void Init();
+
+ template <typename Particle, typename Track, typename Stack>
+ corsika::process::EProcessReturn DoContinuous(Particle& p, Track& t, Stack&) {
+ using namespace corsika::units::si;
+ auto const start = t.GetPosition(0).GetCoordinates();
+ auto const delta = t.GetPosition(1).GetCoordinates() - start;
+ auto const& name = corsika::particles::GetName(p.GetPID());
+
+ fFile << name << " " << p.GetEnergy() / 1_eV << ' ' << start[0] / 1_m << ' '
+ << start[1] / 1_m << ' ' << start[2] / 1_m << " " << delta[0] / 1_m << ' '
+ << delta[1] / 1_m << ' ' << delta[2] / 1_m << '\n';
+
+ return corsika::process::EProcessReturn::eOk;
+ }
+
+ template <typename Particle, typename Track>
+ corsika::units::si::LengthType MaxStepLength(Particle&, Track&) {
+ return corsika::units::si::meter * std::numeric_limits<double>::infinity();
+ }
+
+ private:
+ std::string const fFilename;
+ std::ofstream fFile;
+ };
+
+} // namespace corsika::process::TrackWriter
+
+#endif
diff --git a/Processes/TrackingLine/TrackingLine.h b/Processes/TrackingLine/TrackingLine.h
index c64522c40349f873f45693f19e37f978972cc13d..31966195751d0a1dfaee6d61a60421a0f80d4676 100644
--- a/Processes/TrackingLine/TrackingLine.h
+++ b/Processes/TrackingLine/TrackingLine.h
@@ -84,7 +84,8 @@ namespace corsika::process {
p.GetMomentum() / p.GetEnergy() * corsika::units::constants::c;
auto const currentPosition = p.GetPosition();
-
+ std::cout << "TrackingLine pid: " << p.GetPID()
+ << " , E = " << p.GetEnergy() / 1_GeV << " GeV" << std::endl;
std::cout << "TrackingLine pos: " << currentPosition.GetCoordinates()
<< std::endl;
std::cout << "TrackingLine E: " << p.GetEnergy() / 1_GeV << " GeV" << std::endl;
diff --git a/Processes/TrackingLine/testTrackingLine.cc b/Processes/TrackingLine/testTrackingLine.cc
index 737b536f2fad7ddae2403870bd94b9e3cab5c7a5..285f275c0ae53ab2b6b128be270860819326cac6 100644
--- a/Processes/TrackingLine/testTrackingLine.cc
+++ b/Processes/TrackingLine/testTrackingLine.cc
@@ -9,6 +9,7 @@
*/
#include <corsika/environment/Environment.h>
+#include <corsika/particles/ParticleProperties.h>
#include <corsika/process/tracking_line/TrackingLine.h>
@@ -47,6 +48,7 @@ struct DummyParticle {
auto GetEnergy() const { return fEnergy; }
auto GetMomentum() const { return fMomentum; }
auto GetPosition() const { return fPosition; }
+ auto GetPID() const { return corsika::particles::Code::Unknown; }
};
struct DummyStack {
diff --git a/Tools/plot_tracks.sh b/Tools/plot_tracks.sh
new file mode 100755
index 0000000000000000000000000000000000000000..135b410dc8d80523c8a5f340bf3c26634e35bb7f
--- /dev/null
+++ b/Tools/plot_tracks.sh
@@ -0,0 +1,39 @@
+#!/bin/sh
+
+# (c) Copyright 2018 CORSIKA Project, corsika-project@lists.kit.edu
+#
+# See file AUTHORS for a list of contributors.
+#
+# 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.
+
+# with this script you can plot an animation of output of TrackWriter
+
+track_dat=$1
+if [ -z "$track_dat" ]; then
+ echo "usage: $0 <track.dat> [output.gif]" >&2
+ exit 1
+fi
+
+output=$2
+if [ -z "$output" ]; then
+ output="$track_dat.gif"
+fi
+
+cat <<EOF | gnuplot
+set term gif animate size 600,600
+set output "$output"
+
+set xlabel "x / m"
+set ylabel "y / m"
+set zlabel "z / m"
+set title "CORSIKA 8 preliminary"
+
+do for [t=0:360:1] {
+ set view 60, t
+ splot "$track_dat" u 3:4:5:6:7:8 w vectors nohead t ""
+}
+EOF
+
+exit $?