diff --git a/.gitlab-ci.yml b/.gitlab-ci.yml
index 877225617a805a9e7dc0f8c5225cbfe083c0bb8b..c89e989d5612562bd824949cbd8171fe93845a0d 100644
--- a/.gitlab-ci.yml
+++ b/.gitlab-ci.yml
@@ -17,7 +17,13 @@ build:
- cd build
- cmake ..
- cmake --build .
- - ctest
+ - ctest -V >& test.log
+ - gzip -9 test.log
+ artifacts:
+ paths:
+ - build/test.log.gz
+ when: on_failure
+
pages:
stage: build
diff --git a/CMakeLists.txt b/CMakeLists.txt
index f0a9f08988fefa78050cc9b625ebe2149df55a03..d780315d8939bee82fdb79fb09ecff9e4178d608 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -37,10 +37,14 @@ if(NOT CMAKE_BUILD_TYPE AND NOT CMAKE_CONFIGURATION_TYPES)
endif()
# enable warnings and disallow non-standard language
-set(CMAKE_CXX_FLAGS "-Wall -pedantic -Wextra -Wno-ignored-qualifiers -Wno-nonportable-include-path")
+set(CMAKE_CXX_FLAGS "-Wall -pedantic -Wextra -Wno-ignored-qualifiers")
set(CMAKE_CXX_FLAGS_DEBUG "${CMAKE_CXX_FLAGS} -O0 -g")
set(CMAKE_CXX_FLAGS_RELEASE "${CMAKE_CXX_FLAGS} -O3")
+# clang produces a lot of unecessary warnings without this:
+add_compile_options("$<$<CXX_COMPILER_ID:Clang>:-Wno-nonportable-include-path>")
+
+
# unit testing coverage, does not work yet
#include (CodeCoverage)
##set(COVERAGE_LCOV_EXCLUDES 'Documentation/*')
diff --git a/Documentation/Examples/cascade_example.cc b/Documentation/Examples/cascade_example.cc
index f3ebb57d541bf585312f47083d92b745bf82fc36..a76ce24b4b2dcbb1f70a583a296b00372e90ff96 100644
--- a/Documentation/Examples/cascade_example.cc
+++ b/Documentation/Examples/cascade_example.cc
@@ -80,6 +80,9 @@ public:
case Code::Electron:
is_em = true;
break;
+ case Code::Positron:
+ is_em = true;
+ break;
case Code::Gamma:
is_em = true;
break;
@@ -234,18 +237,27 @@ int main() {
// setup particle stack, and add primary particle
setup::Stack stack;
stack.Clear();
- const hep::EnergyType E0 = 10_TeV;
+ const hep::EnergyType E0 = 100_TeV;
+ double theta = 0.;
+ double phi = 0.;
{
auto particle = stack.NewParticle();
particle.SetPID(Code::Proton);
hep::MomentumType P0 = sqrt(E0 * E0 - Proton::GetMass() * Proton::GetMass());
- auto plab = stack::super_stupid::MomentumVector(rootCS, 0_GeV, 0_GeV, -P0);
+ auto momentumComponents = [](double theta, double phi, MomentumType& ptot) {
+ return std::make_tuple(ptot * sin(theta) * cos(phi), ptot * sin(theta) * sin(phi),
+ -ptot * cos(theta));
+ };
+ auto const [px, py, pz] =
+ momentumComponents(theta / 180. * M_PI, phi / 180. * M_PI, P0);
+ auto plab = stack::super_stupid::MomentumVector(rootCS, {px, py, pz});
+ cout << "input angles: theta=" << theta << " phi=" << phi << endl;
+ cout << "input momentum: " << plab.GetComponents() / 1_GeV << endl;
particle.SetEnergy(E0);
particle.SetMomentum(plab);
particle.SetTime(0_ns);
Point p(rootCS, 0_m, 0_m, 0_m);
particle.SetPosition(p);
- cout << particle.GetEnergy() / 1_GeV << endl;
}
// define air shower object, run simulation
diff --git a/Processes/Sibyll/Decay.h b/Processes/Sibyll/Decay.h
index 3d118d2180778de5225910e47d09ed4786c36f67..74bc02418b91f993c0d8ef3017c3dc91b52bd642 100644
--- a/Processes/Sibyll/Decay.h
+++ b/Processes/Sibyll/Decay.h
@@ -21,6 +21,8 @@
#include <corsika/particles/ParticleProperties.h>
+#include <fenv.h>
+
namespace corsika::process {
namespace sibyll {
@@ -143,9 +145,19 @@ namespace corsika::process {
corsika::particles::GetLifetime(p.GetPID());
auto const lifetime = gamma * t0;
+ const auto mkin =
+ (E * E - p.GetMomentum().squaredNorm()); // delta_mass(p.GetMomentum(), E, m);
cout << "Decay: code: " << p.GetPID() << endl;
cout << "Decay: MinStep: t0: " << t0 << endl;
cout << "Decay: MinStep: energy: " << E / 1_GeV << " GeV" << endl;
+ cout << "Decay: momentum: " << p.GetMomentum().GetComponents() / 1_GeV << " GeV"
+ << endl;
+ cout << "Decay: momentum: shell mass-kin. inv. mass " << mkin / 1_GeV / 1_GeV
+ << " " << m / 1_GeV * m / 1_GeV << endl;
+ // cout << "Decay: sib mass: " << s_mass1_.am2[
+ // process::sibyll::ConvertToSibyllRaw(p.GetPID()) ] << endl;
+ auto sib_id = process::sibyll::ConvertToSibyllRaw(p.GetPID());
+ cout << "Decay: sib mass: " << get_sibyll_mass2(sib_id) << endl;
cout << "Decay: MinStep: gamma: " << gamma << endl;
cout << "Decay: MinStep: tau: " << lifetime << endl;
@@ -156,6 +168,9 @@ namespace corsika::process {
void DoDecay(Particle& p, Stack& s) {
using corsika::geometry::Point;
using namespace corsika::units::si;
+
+ feenableexcept(FE_INVALID);
+
fCount++;
SibStack ss;
ss.Clear();
@@ -165,6 +180,10 @@ namespace corsika::process {
pin.SetPID(process::sibyll::ConvertToSibyllRaw(pCode));
pin.SetEnergy(p.GetEnergy());
pin.SetMomentum(p.GetMomentum());
+ // setting particle mass with Corsika values, may be inconsistent with sibyll
+ // internal values
+#warning setting particle mass with Corsika values, may be inconsistent with sibyll internal values
+ pin.SetMass(corsika::particles::GetMass(pCode));
// remember position
Point decayPoint = p.GetPosition();
TimeType t0 = p.GetTime();
@@ -181,15 +200,12 @@ namespace corsika::process {
// print output
int print_unit = 6;
sib_list_(print_unit);
+
// copy particles from sibyll stack to corsika
- int i = -1;
for (auto& psib : ss) {
- ++i;
// FOR NOW: skip particles that have decayed in Sibyll, move to iterator?
- if (abs(s_plist_.llist[i]) > 100) continue;
+ if (psib.HasDecayed()) continue;
// add to corsika stack
- // cout << "decay product: " << process::sibyll::ConvertFromSibyll(
- // psib.GetPID() ) << endl;
auto pnew = s.NewParticle();
pnew.SetEnergy(psib.GetEnergy());
pnew.SetPID(process::sibyll::ConvertFromSibyll(psib.GetPID()));
@@ -199,6 +215,8 @@ namespace corsika::process {
}
// empty sibyll stack
ss.Clear();
+
+ fedisableexcept(FE_INVALID);
}
};
} // namespace sibyll
diff --git a/Processes/Sibyll/Interaction.h b/Processes/Sibyll/Interaction.h
index 2b03b26009c1d6f0e4f9c9be75d9452d2715124f..e9edec807e1ba1df784f5a6ecd15517ddbde7425 100644
--- a/Processes/Sibyll/Interaction.h
+++ b/Processes/Sibyll/Interaction.h
@@ -27,10 +27,14 @@ namespace corsika::process::sibyll {
class Interaction : public corsika::process::InteractionProcess<Interaction> {
mutable int fCount = 0;
+ mutable int fNucCount = 0;
public:
Interaction() {}
- ~Interaction() { std::cout << "Sibyll::Interaction n=" << fCount << std::endl; }
+ ~Interaction() {
+ std::cout << "Sibyll::Interaction n=" << fCount << " Nnuc=" << fNucCount
+ << std::endl;
+ }
void Init() {
@@ -147,15 +151,28 @@ namespace corsika::process::sibyll {
// 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]);
+ QuantityVector<length_d> const yAxis{0_m, 1_m, 0_m};
+ auto rotation_angles = [](MomentumVector const& pin) {
+ const auto p = pin.GetComponents();
+ const auto th = acos(p[2] / p.norm());
+ const auto ph = atan2(
+ p[1] / 1_GeV, p[0] / 1_GeV); // acos( p[0] / sqrt(p[0]*p[0]+p[1]*p[1] ) );
+ return std::make_tuple(th, ph);
};
- 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);
+ // 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());
+ auto const [theta, phi] = rotation_angles(p.GetMomentum());
+ cout << "ProcessSibyll: zenith angle between sibyllCS and rootCS: "
+ << theta / M_PI * 180. << endl;
+ cout << "ProcessSibyll: azimuth angle between sibyllCS and rootCS: "
+ << phi / M_PI * 180. << endl;
+ // double phi = asin( p.GetMomentum().GetComponents()[0]/pt(p.GetMomentum() ) );
+ const CoordinateSystem tempCS = rootCS.rotate(zAxis, phi);
+ const CoordinateSystem sibyllCS = tempCS.rotate(yAxis, theta);
/*
the target should be defined by the Environment,
@@ -232,6 +249,7 @@ namespace corsika::process::sibyll {
// print final state
int print_unit = 6;
sib_list_(print_unit);
+ fNucCount += get_nwounded() - 1;
// delete current particle
p.Delete();
diff --git a/Processes/Sibyll/SibStack.h b/Processes/Sibyll/SibStack.h
index 3364e7deac87ca2fda683d188aeaf3495eaa1fb4..bc3dc4c0fba66fadba4e277c0e305f01b65b1131 100644
--- a/Processes/Sibyll/SibStack.h
+++ b/Processes/Sibyll/SibStack.h
@@ -39,6 +39,11 @@ namespace corsika::process::sibyll {
using namespace corsika::units::hep;
s_plist_.p[3][i] = v / 1_GeV;
}
+ void SetMass(const int i, const corsika::units::hep::MassType v) {
+ using namespace corsika::units::hep;
+ s_plist_.p[4][i] = v / 1_GeV;
+ }
+
void SetMomentum(const int i, const MomentumVector& v) {
using namespace corsika::units;
using namespace corsika::units::hep;
@@ -52,6 +57,10 @@ namespace corsika::process::sibyll {
using namespace corsika::units::hep;
return s_plist_.p[3][i] * 1_GeV;
}
+ corsika::units::hep::EnergyType GetMass(const int i) const {
+ using namespace corsika::units::hep;
+ return s_plist_.p[4][i] * 1_GeV;
+ }
MomentumVector GetMomentum(const int i) const {
using corsika::geometry::CoordinateSystem;
@@ -91,6 +100,14 @@ namespace corsika::process::sibyll {
corsika::units::hep::EnergyType GetEnergy() const {
return GetStackData().GetEnergy(GetIndex());
}
+
+ void SetMass(const corsika::units::hep::MassType v) {
+ GetStackData().SetMass(GetIndex(), v);
+ }
+ corsika::units::hep::EnergyType GetMass() const {
+ return GetStackData().GetMass(GetIndex());
+ }
+
bool HasDecayed() const {
return abs(GetStackData().GetId(GetIndex())) > 100 ? true : false;
}
diff --git a/Processes/Sibyll/sibyll2.3c.cc b/Processes/Sibyll/sibyll2.3c.cc
index 630b913c66d444bf5311e272b5640759eef9e2f9..f614b596bdcb376fc138af70a1af6539852b54f0 100644
--- a/Processes/Sibyll/sibyll2.3c.cc
+++ b/Processes/Sibyll/sibyll2.3c.cc
@@ -14,6 +14,9 @@
#include <corsika/random/RNGManager.h>
#include <random>
+int get_nwounded() { return s_chist_.nwd; }
+double get_sibyll_mass2(int& id) { return s_mass1_.am2[id]; }
+
double s_rndm_(int&) {
static corsika::random::RNG& rng =
corsika::random::RNGManager::GetInstance().GetRandomStream("s_rndm");
diff --git a/Processes/Sibyll/sibyll2.3c.h b/Processes/Sibyll/sibyll2.3c.h
index d915e71401ab397b3966f96b615ba0b73e2b7dee..104b55eb0bf2e6fc12277f6d62bf1adcf567acb4 100644
--- a/Processes/Sibyll/sibyll2.3c.h
+++ b/Processes/Sibyll/sibyll2.3c.h
@@ -30,6 +30,10 @@ extern struct {
int np;
} s_plist_;
+// additional information about interactions.
+// number of wounded nucleons, number of hard and soft scatterings etc.
+extern struct { int nnsof[20], nnjet[20], jdif[20], nwd, njet, nsof; } s_chist_;
+
extern struct {
double cbr[223 + 16 + 12 + 8];
int kdec[1338 + 6 * (16 + 12 + 8)];
@@ -96,6 +100,9 @@ void sib_sigma_hp_(const int&, const double&, double&, double&, double&, double*
double s_rndm_(int&);
+int get_nwounded();
+double get_sibyll_mass2(int&);
+
// phojet random generator setup
void pho_rndin_(int&, int&, int&, int&);
}