diff --git a/.gitlab-ci.yml b/.gitlab-ci.yml
index 879e1ac58766a03b0afd06f012e3bb207c4e8db8..c83d5488bb8ccdb2ad14b619b74bf39328d88ab9 100644
--- a/.gitlab-ci.yml
+++ b/.gitlab-ci.yml
@@ -18,8 +18,10 @@ build:
- cmake ..
- cmake --build .
- ctest -j4 -V >& test.log
- - gzip -9 -S .gz test.log
+ after_script:
+ - cd build
- ls
+ - gzip -v -9 -S .gz test.log
- pwd
artifacts:
expire_in: 1 week
diff --git a/Documentation/Examples/cascade_example.cc b/Documentation/Examples/cascade_example.cc
index 772853db877c6f099370b3b8f35ee949b2129c45..67b6f4d9160b18cfea3ed42b5564d9f24b688098 100644
--- a/Documentation/Examples/cascade_example.cc
+++ b/Documentation/Examples/cascade_example.cc
@@ -71,12 +71,22 @@ public:
template <typename Particle>
bool isBelowEnergyCut(Particle& p) const {
- // FOR NOW: center-of-mass energy hard coded
- const HEPEnergyType Ecm = sqrt(2. * p.GetEnergy() * 0.93827_GeV);
- if (p.GetEnergy() < fECut || Ecm < 10_GeV)
- return true;
- else
- return false;
+ // nuclei
+ if (p.GetPID() == corsika::particles::Code::Nucleus) {
+ auto const ElabNuc = p.GetEnergy() / p.GetNuclearA();
+ auto const EcmNN = sqrt(2. * ElabNuc * 0.93827_GeV);
+ if (ElabNuc < fECut || EcmNN < 10_GeV)
+ return true;
+ else
+ return false;
+ } else {
+ // TODO: center-of-mass energy hard coded
+ const HEPEnergyType Ecm = sqrt(2. * p.GetEnergy() * 0.93827_GeV);
+ if (p.GetEnergy() < fECut || Ecm < 10_GeV)
+ return true;
+ else
+ return false;
+ }
}
bool isEmParticle(Code pCode) const {
@@ -244,7 +254,6 @@ int main() {
corsika::random::RNGManager::GetInstance().RegisterRandomStream("s_rndm");
corsika::process::sibyll::Interaction sibyll(env);
- // corsika::process::sibyll::NuclearInteraction sibyllNuc(env);
corsika::process::sibyll::NuclearInteraction sibyllNuc(env, sibyll);
corsika::process::sibyll::Decay decay;
ProcessCut cut(20_GeV);
@@ -265,9 +274,14 @@ int main() {
// setup particle stack, and add primary particle
setup::Stack stack;
stack.Clear();
- const Code beamCode = Code::Proton;
+ const Code beamCode = Code::Nucleus;
+ const int nuclA = 56;
+ const int nuclZ = int(nuclA / 2.15 + 0.7);
+ const HEPMassType mass = corsika::particles::Proton::GetMass() * nuclZ +
+ (nuclA - nuclZ) * corsika::particles::Neutron::GetMass();
const HEPEnergyType E0 =
- 100_TeV; // 1_PeV crashes with bad COMboost in second interaction (crash later)
+ nuclA *
+ 100_GeV; // 1_PeV crashes with bad COMboost in second interaction (crash later)
double theta = 0.;
double phi = 0.;
@@ -275,7 +289,7 @@ int main() {
auto elab2plab = [](HEPEnergyType Elab, HEPMassType m) {
return sqrt(Elab * Elab - m * m);
};
- HEPMomentumType P0 = elab2plab(E0, corsika::particles::GetMass(beamCode));
+ HEPMomentumType P0 = elab2plab(E0, mass);
auto momentumComponents = [](double theta, double phi, HEPMomentumType ptot) {
return std::make_tuple(ptot * sin(theta) * cos(phi), ptot * sin(theta) * sin(phi),
-ptot * cos(theta));
@@ -287,7 +301,7 @@ int main() {
cout << "input angles: theta=" << theta << " phi=" << phi << endl;
cout << "input momentum: " << plab.GetComponents() / 1_GeV << endl;
Point pos(rootCS, 0_m, 0_m, 0_m);
- stack.AddParticle(beamCode, E0, plab, pos, 0_ns);
+ stack.AddParticle(beamCode, E0, plab, pos, 0_ns, nuclA, nuclZ);
}
// define air shower object, run simulation
diff --git a/Framework/Particles/NuclearData.xml b/Framework/Particles/NuclearData.xml
index 0fff50c33839093fc8bfacb9a7632f36110a2492..73eda0ef2350aa6849df60c61a34231dd9fff869 100644
--- a/Framework/Particles/NuclearData.xml
+++ b/Framework/Particles/NuclearData.xml
@@ -18,19 +18,19 @@
<particle id="1000020030" name="helium3" A="3" Z="2" >
</particle>
-<particle id="1000020030" name="lithium" A="6" Z="3" >
+<particle id="1000030070" name="lithium7" A="7" Z="3" >
</particle>
-<particle id="1000020030" name="beryllium" A="9" Z="4" >
+<particle id="1000040090" name="beryllium9" A="9" Z="4" >
</particle>
-<particle id="1000020030" name="boron" A="10" Z="5" >
+<particle id="1000110050" name="boron11" A="11" Z="5" >
</particle>
<particle id="1000060120" name="carbon" A="12" Z="6" >
</particle>
-<particle id="1000060120" name="carbon13" A="13" Z="6" >
+<particle id="1000060130" name="carbon13" A="13" Z="6" >
</particle>
<particle id="1000070140" name="nitrogen" A="14" Z="7" >
@@ -39,10 +39,10 @@
<particle id="1000080160" name="oxygen" A="16" Z="8" >
</particle>
-<particle id="1000080160" name="fluor" A="18" Z="9" >
+<particle id="1000080180" name="fluor" A="18" Z="9" >
</particle>
-<particle id="1000100220" name="neon21" A="21" Z="10" >
+<particle id="1000100210" name="neon21" A="21" Z="10" >
</particle>
<particle id="1000100220" name="neon" A="22" Z="10" >
diff --git a/Framework/Particles/ParticleProperties.cc b/Framework/Particles/ParticleProperties.cc
index 4932492ddf33028b456670e9ca2da0a5b8d0efe2..b90a744bd38f5dc2755a2f10753af2e7c7386081 100644
--- a/Framework/Particles/ParticleProperties.cc
+++ b/Framework/Particles/ParticleProperties.cc
@@ -27,5 +27,16 @@ namespace corsika::particles {
std::ostream& operator<<(std::ostream& stream, corsika::particles::Code const p) {
return stream << corsika::particles::GetName(p);
}
+
+ Code ConvertFromPDG(PDGCode p) {
+ static_assert(detail::conversionArray.size() % 2 == 1);
+ auto constexpr maxPDG{(detail::conversionArray.size() - 1) >> 1};
+ auto k = static_cast<PDGCodeType>(p);
+ if (abs(k) <= maxPDG) {
+ return detail::conversionArray[k + maxPDG];
+ } else {
+ return detail::conversionMap.at(p);
+ }
+ }
} // namespace corsika::particles
diff --git a/Framework/Particles/ParticleProperties.h b/Framework/Particles/ParticleProperties.h
index ad1a14e620899089e94093722c630e218b8db4c2..7b2b6f7e5594efa6aa4c5014f002bcc78d328701 100644
--- a/Framework/Particles/ParticleProperties.h
+++ b/Framework/Particles/ParticleProperties.h
@@ -23,13 +23,12 @@
#include <iosfwd>
#include <type_traits>
-#include <corsika/units/PhysicalConstants.h>
#include <corsika/units/PhysicalUnits.h>
/**
*
* The properties of all elementary particles is stored here. The data
- * is taken from the Pythia ParticleData.xml file.
+ * are taken from the Pythia ParticleData.xml file.
*
*/
@@ -40,14 +39,15 @@ namespace corsika::particles {
* The Code enum is the actual place to define CORSIKA 8 particle codes.
*/
enum class Code : int16_t;
+ enum class PDGCode : int32_t;
using CodeIntType = std::underlying_type<Code>::type;
- using PDGCodeType = int32_t;
+ using PDGCodeType = std::underlying_type<PDGCode>::type;
// forward declarations to be used in GeneratedParticleProperties
int16_t constexpr GetElectricChargeNumber(Code const);
corsika::units::si::ElectricChargeType constexpr GetElectricCharge(Code const);
corsika::units::si::HEPMassType constexpr GetMass(Code const);
- PDGCodeType constexpr GetPDG(Code const);
+ PDGCode constexpr GetPDG(Code const);
constexpr std::string const& GetName(Code const);
corsika::units::si::TimeType constexpr GetLifetime(Code const);
@@ -58,46 +58,52 @@ namespace corsika::particles {
#include <corsika/particles/GeneratedParticleProperties.inc>
/*!
- * returns mass of particle
+ * returns mass of particle in natural units
*/
corsika::units::si::HEPMassType constexpr GetMass(Code const p) {
- return detail::masses[static_cast<CodeIntType const>(p)];
+ return detail::masses[static_cast<CodeIntType>(p)];
}
- PDGCodeType constexpr GetPDG(Code const p) {
- return detail::pdg_codes[static_cast<CodeIntType const>(p)];
+ /*!
+ * returns PDG id
+ */
+ PDGCode constexpr GetPDG(Code const p) {
+ return detail::pdg_codes[static_cast<CodeIntType>(p)];
}
/*!
- * returns electric charge of particle / (e/3).
+ * returns electric charge of particle / (e/3), e.g. return 3 for a proton.
*/
int16_t constexpr GetElectricChargeNumber(Code const p) {
- return detail::electric_charges[static_cast<CodeIntType const>(p)];
+ return detail::electric_charges[static_cast<CodeIntType>(p)];
}
+ /*!
+ * returns electric charge of particle, e.g. return 1.602e-19_C for a proton.
+ */
corsika::units::si::ElectricChargeType constexpr GetElectricCharge(Code const p) {
- return GetElectricChargeNumber(p) * (corsika::units::constants::e / 3.);
+ return GetElectricChargeNumber(p) * (corsika::units::constants::e * (1. / 3.));
}
constexpr std::string const& GetName(Code const p) {
- return detail::names[static_cast<CodeIntType const>(p)];
+ return detail::names[static_cast<CodeIntType>(p)];
}
corsika::units::si::TimeType constexpr GetLifetime(Code const p) {
- return detail::lifetime[static_cast<CodeIntType const>(p)] *
+ return detail::lifetime[static_cast<CodeIntType>(p)] *
corsika::units::si::second;
}
bool constexpr IsNucleus(Code const p) {
- return detail::isNucleus[static_cast<CodeIntType const>(p)];
+ return detail::isNucleus[static_cast<CodeIntType>(p)];
}
int constexpr GetNucleusA(Code const p) {
- return detail::nucleusA[static_cast<CodeIntType const>(p)];
+ return detail::nucleusA[static_cast<CodeIntType>(p)];
}
int constexpr GetNucleusZ(Code const p) {
- return detail::nucleusZ[static_cast<CodeIntType const>(p)];
+ return detail::nucleusZ[static_cast<CodeIntType>(p)];
}
/**
@@ -105,7 +111,8 @@ namespace corsika::particles {
**/
std::ostream& operator<<(std::ostream& stream, corsika::particles::Code const p);
-
+
+ Code ConvertFromPDG(PDGCode);
} // namespace corsika::particles
#endif
diff --git a/Framework/Particles/pdxml_reader.py b/Framework/Particles/pdxml_reader.py
index 99544cf5029c933ef9a43478dae73d28c2efbb56..946940f7dcf54551e2888e438c6eeeae0e1659e1 100755
--- a/Framework/Particles/pdxml_reader.py
+++ b/Framework/Particles/pdxml_reader.py
@@ -4,6 +4,7 @@ import sys, math, itertools, re, csv, pprint
import xml.etree.ElementTree as ET
from collections import OrderedDict
import pickle
+import io
GeVfm = 0.19732696312541853
c_speed_of_light = 29.9792458e10 # mm / s
@@ -39,7 +40,7 @@ def parsePythia(filename):
ctau = 0.
else:
print ("missing lifetime: " + str(pdg_id) + " " + str(name))
- sys.exit(0)
+ sys.exit(1)
yield (pdg_id, name, mass, electric_charge, antiName, ctau/c_speed_of_light)
@@ -93,42 +94,7 @@ def class_names(filename):
pdg_id = int(particle.attrib["pdgID"])
map[pdg_id] = name
- return map
-
-
-
-##############################################################
-#
-# Automatically produce a string qualifying as C++ class name
-#
-# This function produces names of type "DELTA_PLUS_PLUS"
-#
-def c_identifier(name):
- orig = name
- name = name.upper()
- for c in "() /":
- name = name.replace(c, "_")
-
- name = name.replace("BAR", "_BAR")
- name = name.replace("0", "_0")
- name = name.replace("*", "_STAR")
- name = name.replace("'", "_PRIME")
- name = name.replace("+", "_PLUS")
- name = name.replace("-", "_MINUS")
-
- while True:
- tmp = name.replace("__", "_")
- if tmp == name:
- break
- else:
- name = tmp
-
- pattern = re.compile(r'^[A-Z_][A-Z_0-9]*$')
- if pattern.match(name):
- return name.strip("_")
- else:
- raise Exception("could not generate C identifier for '{:s}'".format(orig))
-
+ return map
##############################################################
#
@@ -163,7 +129,7 @@ def c_identifier_camel(name):
name = tmp
name.strip("_")
- # remove all "_", if this does not by accident concatenate two number
+ # remove all "_", if this does not by accident concatenate two numbers
istart = 0
while True:
i = name.find('_', istart)
@@ -253,6 +219,47 @@ def read_nuclei_db(filename, particle_db, classnames):
return particle_db
+###############################################################
+#
+# build conversion table PDG -> ngc
+#
+def gen_conversion_PDG_ngc(particle_db):
+ # todo: find a optimum value, think about cache miss with array vs lookup time with map
+ P_MAX = 500 # the maximum PDG code that is still filled into the table
+
+ conversionDict = dict()
+ conversionTable = [None] * (2*P_MAX + 1)
+ for cId, p in particle_db.items():
+ pdg = p['pdg']
+
+ if abs(pdg) < P_MAX:
+ if conversionTable[pdg + P_MAX]:
+ raise Exception("table entry already occupied")
+ else:
+ conversionTable[pdg + P_MAX] = cId
+ else:
+ if pdg in conversionDict.keys():
+ raise Exception(f"map entry {pdg} already occupied")
+ else:
+ conversionDict[pdg] = cId
+
+ output = io.StringIO()
+ def oprint(*args, **kwargs):
+ print(*args, **kwargs, file=output)
+
+ oprint(f"static std::array<Code, {len(conversionTable)}> constexpr conversionArray {{")
+ for ngc in conversionTable:
+ oprint(" Code::{0},".format(ngc if ngc else "Unknown"))
+ oprint("};")
+ oprint()
+
+ oprint("static std::map<PDGCode, Code> const conversionMap {")
+ for ngc in conversionDict.values():
+ oprint(f" {{PDGCode::{ngc}, Code::{ngc}}},")
+ oprint("};")
+ oprint()
+
+ return output.getvalue()
###############################################################
@@ -260,7 +267,7 @@ def read_nuclei_db(filename, particle_db, classnames):
# return string with enum of all internal particle codes
#
def gen_internal_enum(particle_db):
- string = ("enum class Code : int16_t {\n"
+ string = ("enum class Code : CodeIntType {\n"
" FirstParticle = 1, // if you want to loop over particles, you want to start with \"1\" \n") # identifier for eventual loops...
@@ -277,6 +284,20 @@ def gen_internal_enum(particle_db):
return string
+###############################################################
+#
+# return string with enum of all PDG particle codes
+#
+def gen_pdg_enum(particle_db):
+ string = "enum class PDGCode : PDGCodeType {\n"
+
+ for cId in particle_db:
+ pdgCode = particle_db[cId]['pdg']
+ string += " {key:s} = {code:d},\n".format(key = cId, code = pdgCode)
+
+ string += " };\n"
+
+ return string
###############################################################
@@ -296,9 +317,9 @@ def gen_properties(particle_db):
string += "};\n\n"
# PDG code table
- string += "static constexpr std::array<PDGCodeType const, size> pdg_codes = {\n"
- for p in particle_db.values():
- string += " {pdg:d}, // {name:s}\n".format(pdg = p['pdg'], name = p['name'])
+ string += "static constexpr std::array<PDGCode, size> pdg_codes = {\n"
+ for p in particle_db.keys():
+ string += f" PDGCode::{p},\n"
string += "};\n"
# name string table
@@ -473,8 +494,10 @@ if __name__ == "__main__":
with open("GeneratedParticleProperties.inc", "w") as f:
print(inc_start(), file=f)
print(gen_internal_enum(particle_db), file=f)
+ print(gen_pdg_enum(particle_db), file=f)
print(detail_start(), file=f)
print(gen_properties(particle_db), file=f)
+ print(gen_conversion_PDG_ngc(particle_db), file=f)
print(detail_end(), file=f)
print(gen_classes(particle_db), file=f)
print(inc_end(), file=f)
diff --git a/Framework/Particles/testParticles.cc b/Framework/Particles/testParticles.cc
index 1ca8b2188b5004c22c3168b5ad90290aa10a8ffb..cba91b476f42e6ee540ebfbf64b5990c4e03ea81 100644
--- a/Framework/Particles/testParticles.cc
+++ b/Framework/Particles/testParticles.cc
@@ -34,6 +34,8 @@ TEST_CASE("ParticleProperties", "[Particles]") {
SECTION("Masses") {
REQUIRE(Electron::GetMass() / (511_keV) == Approx(1));
REQUIRE(Electron::GetMass() / GetMass(Code::Electron) == Approx(1));
+
+ REQUIRE((Proton::GetMass() + Neutron::GetMass()) / corsika::units::constants::nucleonMass == Approx(2));
}
SECTION("Charges") {
@@ -50,11 +52,23 @@ TEST_CASE("ParticleProperties", "[Particles]") {
}
SECTION("PDG") {
- REQUIRE(GetPDG(Code::PiPlus) == 211);
- REQUIRE(GetPDG(Code::DPlus) == 411);
- REQUIRE(GetPDG(Code::NuMu) == 14);
- REQUIRE(GetPDG(Code::NuE) == 12);
- REQUIRE(GetPDG(Code::MuMinus) == 13);
+ REQUIRE(GetPDG(Code::PiPlus) == PDGCode::PiPlus);
+ REQUIRE(GetPDG(Code::DPlus) == PDGCode::DPlus);
+ REQUIRE(GetPDG(Code::NuMu) == PDGCode::NuMu);
+ REQUIRE(GetPDG(Code::NuE) == PDGCode::NuE);
+ REQUIRE(GetPDG(Code::MuMinus) == PDGCode::MuMinus);
+
+ REQUIRE(static_cast<int>(GetPDG(Code::PiPlus)) == 211);
+ REQUIRE(static_cast<int>(GetPDG(Code::DPlus)) == 411);
+ REQUIRE(static_cast<int>(GetPDG(Code::NuMu)) == 14);
+ REQUIRE(static_cast<int>(GetPDG(Code::NuEBar)) == -12);
+ REQUIRE(static_cast<int>(GetPDG(Code::MuMinus)) == 13);
+ }
+
+ SECTION("Conversion PDG -> internal") {
+ REQUIRE(ConvertFromPDG(PDGCode::KStarMinus) == Code::KStarMinus);
+ REQUIRE(ConvertFromPDG(PDGCode::MuPlus) == Code::MuPlus);
+ REQUIRE(ConvertFromPDG(PDGCode::SigmaStarCMinusBar) == Code::SigmaStarCMinusBar);
}
SECTION("Lifetimes") {
@@ -70,12 +84,12 @@ TEST_CASE("ParticleProperties", "[Particles]") {
}
SECTION("Nuclei") {
- REQUIRE(IsNucleus(Code::Gamma) == false);
- REQUIRE(IsNucleus(Code::Argon) == true);
- REQUIRE(IsNucleus(Code::Proton) == false);
- REQUIRE(IsNucleus(Code::Hydrogen) == true);
- REQUIRE(Argon::IsNucleus() == true);
- REQUIRE(EtaC::IsNucleus() == false);
+ REQUIRE_FALSE(IsNucleus(Code::Gamma));
+ REQUIRE(IsNucleus(Code::Argon));
+ REQUIRE_FALSE(IsNucleus(Code::Proton));
+ REQUIRE(IsNucleus(Code::Hydrogen));
+ REQUIRE(Argon::IsNucleus());
+ REQUIRE_FALSE(EtaC::IsNucleus());
REQUIRE(GetNucleusA(Code::Hydrogen) == 1);
REQUIRE(GetNucleusA(Code::Tritium) == 3);
diff --git a/Framework/Units/PhysicalConstants.h b/Framework/Units/PhysicalConstants.h
index 38c59992fe928228754e3a16658043308e71ab68..1bf9872a03ed3abee9bcb03ef93fd24250f14a50 100644
--- a/Framework/Units/PhysicalConstants.h
+++ b/Framework/Units/PhysicalConstants.h
@@ -61,7 +61,9 @@ namespace corsika::units::constants {
// unified atomic mass unit
constexpr quantity<mass_d> u{Rep(1.6605402e-27L) * kilogram};
-
+
+ auto constexpr nucleonMass = 0.5 * (0.93827 + 0.93957) * 1e9 * electronvolt;
+
// etc.
} // namespace corsika::units::constants
diff --git a/Framework/Units/PhysicalUnits.h b/Framework/Units/PhysicalUnits.h
index 90b98e170f2a912018e1486cac26b7da1c85f3ff..af4dbb966aaaf5dfe5bfe7ace9619c9d299e94c7 100644
--- a/Framework/Units/PhysicalUnits.h
+++ b/Framework/Units/PhysicalUnits.h
@@ -18,28 +18,11 @@ namespace phys::units {
/**
* @file PhysicalUnits
*
- * Add new units and types we need. Units are compile-time. We support
- * different system of units in parallel. Literals are used for
- * optimal coding style.
+ * Add new units and types we need. Units are compile-time. Literals are
+ * used for optimal coding style.
*
*/
-/*
-namespace corsika::units::hep {
- using namespace phys::units;
- using namespace phys::units::literals;
- using namespace phys::units::io;
-
- /// defining HEP energy, mass, momentum
- using energy_hep_d = phys::units::energy_d;
- constexpr phys::units::quantity<energy_hep_d> GeV{corsika::units::constants::eV};
- // corsika::units::constants::e / phys::units::coulomb * phys::units::joule };
-
- using MassType = phys::units::quantity<energy_hep_d, double>;
- using MomentumType = phys::units::quantity<energy_hep_d, double>;
- using EnergyType = phys::units::quantity<energy_hep_d, double>;
-} // namespace corsika::units::hep
-*/
namespace corsika::units::si {
using namespace phys::units;
using namespace phys::units::literals;
diff --git a/Processes/Sibyll/Interaction.h b/Processes/Sibyll/Interaction.h
index 966c2f7f03afd4a60dd1dfc7dab66540d1f636b8..69755f00eb2e2923bf43e28935f0617b1ccce860 100644
--- a/Processes/Sibyll/Interaction.h
+++ b/Processes/Sibyll/Interaction.h
@@ -56,6 +56,11 @@ namespace corsika::process::sibyll {
bool WasInitialized() { return fInitialized; }
+ bool ValidCoMEnergy(corsika::units::si::HEPEnergyType ecm) {
+ using namespace corsika::units::si;
+ return (10_GeV < ecm) && (ecm < 1_PeV);
+ }
+
std::tuple<corsika::units::si::CrossSectionType, corsika::units::si::CrossSectionType,
int>
GetCrossSection(const corsika::particles::Code BeamId,
@@ -124,6 +129,7 @@ namespace corsika::process::sibyll {
<< " beam can interact:" << kInteraction << endl
<< " beam pid:" << p.GetPID() << endl;
+ // TODO: move limits into variables
if (kInteraction && Elab >= 8.5_GeV && ECoM >= 10_GeV) {
// get target from environment
diff --git a/Processes/Sibyll/NuclearInteraction.h b/Processes/Sibyll/NuclearInteraction.h
index 50384ddeba4657bcc6988dc50acd4a0614e22a6c..7b7492d4c335477e196116f003e405d3829177bd 100644
--- a/Processes/Sibyll/NuclearInteraction.h
+++ b/Processes/Sibyll/NuclearInteraction.h
@@ -47,7 +47,7 @@ namespace corsika::process::sibyll {
using std::endl;
// initialize hadronic interaction module
- // TODO: save to run multiple initializations?
+ // TODO: safe to run multiple initializations?
if (!fHadronicInteraction.WasInitialized()) fHadronicInteraction.Init();
// initialize nuclib
@@ -56,38 +56,40 @@ namespace corsika::process::sibyll {
}
// TODO: remove number of nucleons, avg target mass is available in environment
- std::tuple<corsika::units::si::CrossSectionType, corsika::units::si::CrossSectionType,
- int>
- GetCrossSection(const corsika::particles::Code BeamId,
- const corsika::particles::Code TargetId,
- const corsika::units::si::HEPEnergyType CoMenergy) {
+ template <typename Particle>
+ std::tuple<corsika::units::si::CrossSectionType, corsika::units::si::CrossSectionType>
+ GetCrossSection(Particle const& p, const corsika::particles::Code TargetId) {
using namespace corsika::units::si;
double sigProd;
+ auto const pCode = p.GetPID();
+ if (pCode != corsika::particles::Code::Nucleus)
+ throw std::runtime_error(
+ "NuclearInteraction: GetCrossSection: particle not a nucleus!");
- std::cout << "NuclearInteraction: GetCrossSection: called with: beamId= " << BeamId
- << " TargetId= " << TargetId << " CoMenergy= " << CoMenergy / 1_GeV
- << std::endl;
-
- if (!corsika::particles::IsNucleus(BeamId)) {
- // ask hadronic interaction model for hadron cross section
- return fHadronicInteraction.GetCrossSection(BeamId, TargetId, CoMenergy);
- }
+ auto const iBeam = p.GetNuclearA();
+ HEPEnergyType LabEnergyPerNuc = p.GetEnergy() / iBeam;
+ std::cout << "NuclearInteraction: GetCrossSection: called with: beamNuclA= "
+ << iBeam << " TargetId= " << TargetId
+ << " LabEnergyPerNuc= " << LabEnergyPerNuc / 1_GeV << std::endl;
// use nuclib to calc. nuclear cross sections
// TODO: for now assumes air with hard coded composition
// extend to arbitrary mixtures, requires smarter initialization
-
- const int iBeam = corsika::particles::GetNucleusA(BeamId);
+ // get nuclib projectile code: nucleon number
if (iBeam > 56 || iBeam < 2) {
- // std::cout<<"NuclearInteraction: beam nucleus outside allowed range for NUCLIB!"
- // << std::endl;
+ std::cout << "NuclearInteraction: beam nucleus outside allowed range for NUCLIB!"
+ << std::endl
+ << "A=" << iBeam << std::endl;
throw std::runtime_error(
"NuclearInteraction: GetCrossSection: beam nucleus outside allowed range for "
"NUCLIB!");
}
- const double dEcm = CoMenergy / 1_GeV;
- if (dEcm < 10.)
+ const double dElabNuc = LabEnergyPerNuc / 1_GeV;
+ // TODO: these limitations are still sibyll specific.
+ // available target nuclei depends on the hadronic interaction model and the
+ // initialization
+ if (dElabNuc < 10.)
throw std::runtime_error("NuclearInteraction: GetCrossSection: energy too low!");
// TODO: these limitations are still sibyll specific.
@@ -99,12 +101,17 @@ namespace corsika::process::sibyll {
throw std::runtime_error(
"Sibyll target outside range. Only nuclei with A<18 are allowed.");
std::cout << "NuclearInteraction: calling signuc.." << std::endl;
- signuc_(iBeam, dEcm, sigProd);
+ std::cout << "WARNING: using hard coded cross section for Nucleus-Air with "
+ "SIBYLL! (fix me!)"
+ << std::endl;
+ // TODO: target id is not used because cross section is still hard coded and fixed
+ // to air.
+ signuc_(iBeam, dElabNuc, sigProd);
std::cout << "cross section: " << sigProd << std::endl;
- return std::make_tuple(sigProd * 1_mbarn, 0_mbarn, iTarget);
+ return std::make_tuple(sigProd * 1_mbarn, 0_mbarn);
}
return std::make_tuple(std::numeric_limits<double>::infinity() * 1_mbarn,
- std::numeric_limits<double>::infinity() * 1_mbarn, 0);
+ std::numeric_limits<double>::infinity() * 1_mbarn);
}
template <typename Particle, typename Track>
@@ -121,11 +128,15 @@ namespace corsika::process::sibyll {
RootCoordinateSystem::GetInstance().GetRootCoordinateSystem();
const particles::Code corsikaBeamId = p.GetPID();
-
if (!corsika::particles::IsNucleus(corsikaBeamId)) {
// no nuclear interaction
return std::numeric_limits<double>::infinity() * 1_g / (1_cm * 1_cm);
}
+ // check if target-style nucleus (enum)
+ if (corsikaBeamId != corsika::particles::Code::Nucleus)
+ throw std::runtime_error(
+ "NuclearInteraction: GetInteractionLength: Wrong nucleus type. Nuclear "
+ "projectiles should use NuclearStackExtension!");
// read from cross section code table
const HEPMassType nucleon_mass = 0.5 * (corsika::particles::Proton::GetMass() +
@@ -136,6 +147,9 @@ namespace corsika::process::sibyll {
// total momentum and energy
HEPEnergyType Elab = p.GetEnergy() + nucleon_mass;
+ int const nuclA = p.GetNuclearA();
+ auto const ElabNuc = p.GetEnergy() / nuclA;
+
MomentumVector pTotLab(rootCS, {0.0_GeV, 0.0_GeV, 0.0_GeV});
pTotLab += p.GetMomentum();
pTotLab += pTarget;
@@ -143,13 +157,19 @@ namespace corsika::process::sibyll {
// calculate cm. energy
const HEPEnergyType ECoM = sqrt(
(Elab + pTotLabNorm) * (Elab - pTotLabNorm)); // binomial for numerical accuracy
-
+ auto const ECoMNN = sqrt(2. * ElabNuc * nucleon_mass);
std::cout << "NuclearInteraction: LambdaInt: \n"
<< " input energy: " << Elab / 1_GeV << endl
<< " input energy CoM: " << ECoM / 1_GeV << endl
- << " beam pid:" << corsikaBeamId << endl;
+ << " beam pid:" << corsikaBeamId << endl
+ << " beam A: " << nuclA << endl
+ << " input energy per nucleon: " << ElabNuc / 1_GeV << endl
+ << " input energy CoM per nucleon: " << ECoMNN / 1_GeV << endl;
+ // throw std::runtime_error("stop here");
- if (Elab >= 8.5_GeV && ECoM >= 10_GeV) {
+ // energy limits
+ // TODO: values depend on hadronic interaction model !! this is sibyll specific
+ if (ElabNuc >= 8.5_GeV && ECoMNN >= 10_GeV) {
// get target from environment
/*
@@ -164,7 +184,6 @@ namespace corsika::process::sibyll {
// determine average interaction length
// weighted sum
int i = -1;
- double avgTargetMassNumber = 0.;
si::CrossSectionType weightedProdCrossSection = 0_mbarn;
// get weights of components from environment/medium
const auto w = mediumComposition.GetFractions();
@@ -173,24 +192,23 @@ namespace corsika::process::sibyll {
i++;
cout << "NuclearInteraction: get interaction length for target: " << targetId
<< endl;
- // TODO: remove number of nucleons, avg target mass is available in environment
- auto const [productionCrossSection, elaCrossSection, numberOfNucleons] =
- GetCrossSection(corsikaBeamId, targetId, ECoM);
+ auto const [productionCrossSection, elaCrossSection] =
+ GetCrossSection(p, targetId);
[[maybe_unused]] auto elaCrossSectionCopy = elaCrossSection;
std::cout << "NuclearInteraction: "
<< "IntLength: nuclib return (mb): "
<< productionCrossSection / 1_mbarn << std::endl;
weightedProdCrossSection += w[i] * productionCrossSection;
- avgTargetMassNumber += w[i] * numberOfNucleons;
}
cout << "NuclearInteraction: "
<< "IntLength: weighted CrossSection (mb): "
<< weightedProdCrossSection / 1_mbarn << endl;
// calculate interaction length in medium
- GrammageType const int_length =
- avgTargetMassNumber * corsika::units::constants::u / weightedProdCrossSection;
+ GrammageType const int_length = mediumComposition.GetAverageMassNumber() *
+ corsika::units::constants::u /
+ weightedProdCrossSection;
std::cout << "NuclearInteraction: "
<< "interaction length (g/cm2): "
<< int_length / (0.001_kg) * 1_cm * 1_cm << std::endl;
@@ -207,9 +225,6 @@ namespace corsika::process::sibyll {
// this routine superimposes different nucleon-nucleon interactions
// in a nucleus-nucleus interaction, based the SIBYLL routine SIBNUC
- // this routine superimposes different nucleon-nucleon interactions
- // in a nucleus-nucleus interaction, based the SIBYLL routine SIBNUC
-
using namespace corsika::units;
using namespace corsika::utl;
using namespace corsika::units::si;
@@ -217,18 +232,29 @@ namespace corsika::process::sibyll {
using std::cout;
using std::endl;
- const auto corsikaProjId = p.GetPID();
- cout << "NuclearInteraction: DoInteraction: called with:" << corsikaProjId << endl
- << "NuclearInteraction: projectile mass: "
- << corsika::particles::GetMass(corsikaProjId) / 1_GeV << endl;
+ const auto ProjId = p.GetPID();
+ // TODO: calculate projectile mass in nuclearStackExtension
+ // const auto ProjMass = p.GetMass();
+ cout << "NuclearInteraction: DoInteraction: called with:" << ProjId << endl;
- if (!IsNucleus(corsikaProjId)) {
+ if (!IsNucleus(ProjId)) {
cout << "WARNING: non nuclear projectile in NUCLIB!" << endl;
// this should not happen
// throw std::runtime_error("Non nuclear projectile in NUCLIB!");
return process::EProcessReturn::eOk;
}
+ // check if target-style nucleus (enum)
+ if (ProjId != corsika::particles::Code::Nucleus)
+ throw std::runtime_error(
+ "NuclearInteraction: DoInteraction: Wrong nucleus type. Nuclear projectiles "
+ "should use NuclearStackExtension!");
+
+ auto const ProjMass =
+ p.GetNuclearZ() * corsika::particles::Proton::GetMass() +
+ (p.GetNuclearA() - p.GetNuclearZ()) * corsika::particles::Neutron::GetMass();
+ cout << "NuclearInteraction: projectile mass: " << ProjMass / 1_GeV << endl;
+
fCount++;
const CoordinateSystem& rootCS =
@@ -242,7 +268,7 @@ namespace corsika::process::sibyll {
cout << "Interaction: time: " << tOrig << endl;
// projectile nucleon number
- const int kAProj = GetNucleusA(corsikaProjId);
+ const int kAProj = p.GetNuclearA(); // GetNucleusA(ProjId);
if (kAProj > 56)
throw std::runtime_error("Projectile nucleus too large for NUCLIB!");
@@ -284,6 +310,13 @@ namespace corsika::process::sibyll {
HEPEnergyType EcmNN = PtotNN4.GetNorm();
cout << "NuclearInteraction: nuc-nuc cm energy: " << EcmNN / 1_GeV << endl;
+ if (!fHadronicInteraction.ValidCoMEnergy(EcmNN)) {
+ cout << "NuclearInteraction: nuc-nuc. CoM energy too low for hadronic "
+ "interaction model!"
+ << endl;
+ throw std::runtime_error("NuclearInteraction: DoInteraction: energy too low!");
+ }
+
// define boost to NUCLEON-NUCLEON frame
COMBoost const boost(PprojNucLab, nucleon_mass);
// boost projecticle
@@ -308,8 +341,9 @@ namespace corsika::process::sibyll {
const auto currentNode = fEnvironment.GetUniverse()->GetContainingNode(pOrig);
const auto& mediumComposition =
currentNode->GetModelProperties().GetNuclearComposition();
- cout << "get cross sections for target materials.." << endl;
+ cout << "get nucleon-nucleus cross sections for target materials.." << endl;
// get cross sections for target materials
+ // using nucleon-target-nucleus cross section!!!
/*
Here we read the cross section from the interaction model again,
should be passed from GetInteractionLength if possible
@@ -320,8 +354,9 @@ namespace corsika::process::sibyll {
for (size_t i = 0; i < compVec.size(); ++i) {
auto const targetId = compVec[i];
cout << "target component: " << targetId << endl;
- cout << "beam id: " << targetId << endl;
- const auto [sigProd, sigEla, nNuc] = GetCrossSection(beamId, targetId, EcmNN);
+ cout << "beam id: " << beamId << endl;
+ const auto [sigProd, sigEla, nNuc] =
+ fHadronicInteraction.GetCrossSection(beamId, targetId, EcmNN);
cross_section_of_components[i] = sigProd;
[[maybe_unused]] auto sigElaCopy = sigEla; // ONLY TO AVOID COMPILER WARNINGS
[[maybe_unused]] auto sigNucCopy = nNuc; // ONLY TO AVOID COMPILER WARNINGS
@@ -396,97 +431,90 @@ namespace corsika::process::sibyll {
<< " px=" << fragments_.ppp[j][0] << " py=" << fragments_.ppp[j][1]
<< " pz=" << fragments_.ppp[j][2] << endl;
- auto Nucleus = [](int iA) {
- // int znew = iA / 2.15 + 0.7;
- corsika::particles::Code pCode;
- switch (iA) {
- case 2:
- pCode = corsika::particles::Deuterium::GetCode();
- break;
-
- case 3:
- pCode = corsika::particles::Tritium::GetCode();
- break;
-
- case 4:
- pCode = corsika::particles::Helium::GetCode();
- break;
-
- case 12:
- pCode = corsika::particles::Carbon::GetCode();
- break;
-
- case 14:
- pCode = corsika::particles::Nitrogen::GetCode();
- break;
-
- case 16:
- pCode = corsika::particles::Oxygen::GetCode();
- break;
-
- case 33:
- pCode = corsika::particles::Sulphur::GetCode();
- break;
-
- case 40:
- pCode = corsika::particles::Argon::GetCode();
- break;
-
- default:
- pCode = corsika::particles::Proton::GetCode();
- }
- return pCode;
- };
-
+ cout << "adding nuclear fragments to particle stack.." << endl;
// put nuclear fragments on corsika stack
for (int j = 0; j < nFragments; ++j) {
- // here we need the nucleonNumber to corsika Id conversion
- auto pCode = Nucleus(AFragments[j]);
+ corsika::particles::Code specCode;
+ const auto nuclA = AFragments[j];
+ // get Z from stability line
+ const auto nuclZ = int(nuclA / 2.15 + 0.7);
+
+ // TODO: do we need to catch single nucleons??
+ if (nuclA == 1)
+ // TODO: sample neutron or proton
+ specCode = corsika::particles::Code::Proton;
+ else
+ specCode = corsika::particles::Code::Nucleus;
+
+ // TODO: mass of nuclei?
+ const HEPMassType mass =
+ corsika::particles::Proton::GetMass() * nuclZ +
+ (nuclA - nuclZ) *
+ corsika::particles::Neutron::GetMass(); // this neglects binding energy
+
+ cout << "NuclearInteraction: adding fragment: " << specCode << endl;
+ cout << "NuclearInteraction: A,Z: " << nuclA << "," << nuclZ << endl;
+ cout << "NuclearInteraction: mass: " << mass / 1_GeV << endl;
// CORSIKA 7 way
// spectators inherit momentum from original projectile
- const double mass_ratio = corsika::particles::GetMass(pCode) /
- corsika::particles::GetMass(corsikaProjId);
+ const double mass_ratio = mass / ProjMass;
+
+ cout << "NuclearInteraction: mass ratio " << mass_ratio << endl;
+
auto const Plab = PprojLab * mass_ratio;
- p.AddSecondary(pCode, Plab.GetTimeLikeComponent(), Plab.GetSpaceLikeComponents(),
- pOrig, tOrig);
+ cout << "NuclearInteraction: fragment momentum: "
+ << Plab.GetSpaceLikeComponents().GetComponents() / 1_GeV << endl;
+
+ if (nuclA == 1)
+ // add nucleon
+ p.AddSecondary(specCode, Plab.GetTimeLikeComponent(),
+ Plab.GetSpaceLikeComponents(), pOrig, tOrig);
+ else
+ // add nucleus
+ p.AddSecondary(specCode, Plab.GetTimeLikeComponent(),
+ Plab.GetSpaceLikeComponents(), pOrig, tOrig, nuclA, nuclZ);
}
// add elastic nucleons to corsika stack
// TODO: the elastic interaction could be external like the inelastic interaction,
// e.g. use existing ElasticModel
+ cout << "adding elastically scattered nucleons to particle stack.." << endl;
for (int j = 0; j < nElasticNucleons; ++j) {
// TODO: sample proton or neutron
- auto pCode = corsika::particles::Proton::GetCode();
+ auto const elaNucCode = corsika::particles::Code::Proton;
// CORSIKA 7 way
// elastic nucleons inherit momentum from original projectile
// neglecting momentum transfer in interaction
- const double mass_ratio = corsika::particles::GetMass(pCode) /
- corsika::particles::GetMass(corsikaProjId);
+ const double mass_ratio = corsika::particles::GetMass(elaNucCode) / ProjMass;
auto const Plab = PprojLab * mass_ratio;
- p.AddSecondary(pCode, Plab.GetTimeLikeComponent(), Plab.GetSpaceLikeComponents(),
- pOrig, tOrig);
+ p.AddSecondary(elaNucCode, Plab.GetTimeLikeComponent(),
+ Plab.GetSpaceLikeComponents(), pOrig, tOrig);
}
// add inelastic interactions
+ cout << "calculate inelastic nucleon-nucleon interactions.." << endl;
for (int j = 0; j < nInelNucleons; ++j) {
-
// TODO: sample neutron or proton
auto pCode = corsika::particles::Proton::GetCode();
// temporarily add to stack, will be removed after interaction in DoInteraction
+ cout << "inelastic interaction no. " << j << endl;
auto inelasticNucleon =
p.AddSecondary(pCode, PprojNucLab.GetTimeLikeComponent(),
PprojNucLab.GetSpaceLikeComponents(), pOrig, tOrig);
// create inelastic interaction
+ cout << "calling HadronicInteraction..." << endl;
fHadronicInteraction.DoInteraction(inelasticNucleon, s);
}
// delete parent particle
p.Delete();
+ cout << "NuclearInteraction: DoInteraction: done" << endl;
+
return process::EProcessReturn::eOk;
}
diff --git a/Processes/Sibyll/SibStack.h b/Processes/Sibyll/SibStack.h
index 9b168f6eb97957b369f87a5c2066b9c601d039f5..fd847124abbcd34b42f2acd5ddb1bc3e5cef2e79 100644
--- a/Processes/Sibyll/SibStack.h
+++ b/Processes/Sibyll/SibStack.h
@@ -27,6 +27,7 @@ namespace corsika::process::sibyll {
public:
void Init();
+ void Dump() const {}
void Clear() { s_plist_.np = 0; }
unsigned int GetSize() const { return s_plist_.np; }
@@ -65,8 +66,7 @@ namespace corsika::process::sibyll {
RootCoordinateSystem::GetInstance().GetRootCoordinateSystem();
QuantityVector<hepmomentum_d> components = {
s_plist_.p[0][i] * 1_GeV, s_plist_.p[1][i] * 1_GeV, s_plist_.p[2][i] * 1_GeV};
- MomentumVector v1(rootCS, components);
- return v1;
+ return MomentumVector(rootCS, components);
}
void Copy(const unsigned int i1, const unsigned int i2) {
diff --git a/Processes/Sibyll/testSibyll.cc b/Processes/Sibyll/testSibyll.cc
index 6251aab11ed5efc5d8b620824ba2ba45847cfff5..8f53c6e12d687dcbfede62a31c3cc43917704ae4 100644
--- a/Processes/Sibyll/testSibyll.cc
+++ b/Processes/Sibyll/testSibyll.cc
@@ -136,13 +136,13 @@ TEST_CASE("SibyllInterface", "[processes]") {
SECTION("NuclearInteractionInterface") {
setup::Stack stack;
- const HEPEnergyType E0 = 100_GeV;
+ const HEPEnergyType E0 = 400_GeV;
HEPMomentumType P0 =
sqrt(E0 * E0 - particles::Proton::GetMass() * particles::Proton::GetMass());
auto plab = stack::super_stupid::MomentumVector(cs, {0_GeV, 0_GeV, -P0});
geometry::Point pos(cs, 0_m, 0_m, 0_m);
- auto particle = stack.AddParticle(particles::Code::Proton, E0, plab, pos, 0_ns);
+ auto particle = stack.AddParticle(particles::Code::Nucleus, E0, plab, pos, 0_ns, 4, 2);
Interaction hmodel(env);
NuclearInteraction model(env, hmodel);
diff --git a/Processes/StackInspector/StackInspector.cc b/Processes/StackInspector/StackInspector.cc
index 312fe40ec82dfd4f1faf12870e397ef4fbe794bc..fb749c0b8c7e6430c7283f079055979d099ee454 100644
--- a/Processes/StackInspector/StackInspector.cc
+++ b/Processes/StackInspector/StackInspector.cc
@@ -51,7 +51,10 @@ process::EProcessReturn StackInspector<Stack>::DoContinuous(Particle&, setup::Tr
auto pos = iterP.GetPosition().GetCoordinates(rootCS);
cout << "StackInspector: i=" << setw(5) << fixed << (i++) << ", id=" << setw(30)
<< iterP.GetPID() << " E=" << setw(15) << scientific << (E / 1_GeV) << " GeV, "
- << " pos=" << pos << endl;
+ << " pos=" << pos;
+ // if (iterP.GetPID()==Code::Nucleus)
+ // cout << " nuc_ref=" << iterP.GetNucleusRef();
+ cout << endl;
}
fCountStep++;
cout << "StackInspector: nStep=" << fCountStep << " stackSize=" << s.GetSize()
diff --git a/Stack/NuclearStackExtension/NuclearStackExtension.h b/Stack/NuclearStackExtension/NuclearStackExtension.h
index 4b3fa4ded55584a4daaff8c7fec0e1525aca48fb..355e9821d806fb3138cb36913b988e07ede4b234 100644
--- a/Stack/NuclearStackExtension/NuclearStackExtension.h
+++ b/Stack/NuclearStackExtension/NuclearStackExtension.h
@@ -24,11 +24,28 @@
namespace corsika::stack {
+ /**
+ * @namespace nuclear_extension
+ *
+ * Add A and Z data to existing stack of particle properties.
+ *
+ * Only for Code::Nucleus particles A and Z are stored, not for all
+ * normal elementary particles.
+ *
+ * Thus in your code, make sure to always check <code>
+ * particle.GetPID()==Code::Nucleus </code> before attempting to
+ * read any nuclear information.
+ *
+ *
+ */
+
namespace nuclear_extension {
typedef corsika::geometry::Vector<corsika::units::si::hepmomentum_d> MomentumVector;
/**
+ * @class NuclearParticleInterface
+ *
* Define ParticleInterface for NuclearStackExtension Stack derived from
* ParticleInterface of Inner stack class
*/
@@ -58,12 +75,12 @@ namespace corsika::stack {
err << "NuclearStackExtension: no A and Z specified for new Nucleus!";
throw std::runtime_error(err.str());
}
- SetNuclearRef(
+ SetNucleusRef(
GetStackData().GetNucleusNextRef()); // store this nucleus data ref
SetNuclearA(vA);
SetNuclearZ(vZ);
} else {
- SetNuclearRef(-1); // this is not a nucleus
+ SetNucleusRef(-1); // this is not a nucleus
}
InnerParticleInterface<StackIteratorInterface>::SetParticleData(
vDataPID, vDataE, vMomentum, vPosition, vTime);
@@ -99,31 +116,40 @@ namespace corsika::stack {
int GetNuclearZ() const { return GetStackData().GetNuclearZ(GetIndex()); }
/// @}
+ // int GetNucleusRef() const { return GetStackData().GetNucleusRef(GetIndex()); }
+
protected:
- void SetNuclearRef(const int vR) { GetStackData().SetNuclearRef(GetIndex(), vR); }
+ void SetNucleusRef(const int vR) { GetStackData().SetNucleusRef(GetIndex(), vR); }
};
/**
- * @class NuclearStackExtensionImpl
+ * @class NuclearStackExtension
*
- * Memory implementation of the extension of particle stack of
- * type InnerStackImpl with nuclear data
+ * Memory implementation of adding nuclear inforamtion to the
+ * existing particle stack defined in class InnerStackImpl.
+ *
+ * Inside the NuclearStackExtension class there is a dedicated
+ * fNucleusRef index, where fNucleusRef[i] is referring to the
+ * correct A and Z for a specific particle index i. fNucleusRef[i]
+ * == -1 means that this is not a nucleus, and a subsequent call to
+ * GetNucleusA would produce an exception.
*/
template <typename InnerStackImpl>
class NuclearStackExtensionImpl : public InnerStackImpl {
public:
void Init() { InnerStackImpl::Init(); }
+ void Dump() { InnerStackImpl::Dump(); }
void Clear() {
InnerStackImpl::Clear();
- fNuclearRef.clear();
+ fNucleusRef.clear();
fNuclearA.clear();
fNuclearZ.clear();
}
- unsigned int GetSize() const { return fNuclearRef.size(); }
- unsigned int GetCapacity() const { return fNuclearRef.size(); }
+ unsigned int GetSize() const { return fNucleusRef.size(); }
+ unsigned int GetCapacity() const { return fNucleusRef.capacity(); }
void SetNuclearA(const unsigned int i, const unsigned short vA) {
fNuclearA[GetNucleusRef(i)] = vA;
@@ -131,7 +157,7 @@ namespace corsika::stack {
void SetNuclearZ(const unsigned int i, const unsigned short vZ) {
fNuclearZ[GetNucleusRef(i)] = vZ;
}
- void SetNuclearRef(const unsigned int i, const int v) { fNuclearRef[i] = v; }
+ void SetNucleusRef(const unsigned int i, const int v) { fNucleusRef[i] = v; }
int GetNuclearA(const unsigned int i) const { return fNuclearA[GetNucleusRef(i)]; }
int GetNuclearZ(const unsigned int i) const { return fNuclearZ[GetNucleusRef(i)]; }
@@ -144,7 +170,7 @@ namespace corsika::stack {
}
int GetNucleusRef(const unsigned int i) const {
- if (fNuclearRef[i] >= 0) return fNuclearRef[i];
+ if (fNucleusRef[i] >= 0) return fNucleusRef[i];
std::ostringstream err;
err << "NuclearStackExtension: no nucleus at ref=" << i;
throw std::runtime_error(err.str());
@@ -154,35 +180,41 @@ namespace corsika::stack {
* Function to copy particle at location i1 in stack to i2
*/
void Copy(const unsigned int i1, const unsigned int i2) {
+ // index range check
if (i1 >= GetSize() || i2 >= GetSize()) {
std::ostringstream err;
err << "NuclearStackExtension: trying to access data beyond size of stack!";
throw std::runtime_error(err.str());
}
+ // copy internal particle data p[i2] = p[i1]
InnerStackImpl::Copy(i1, i2);
- const int ref1 = fNuclearRef[i1];
- const int ref2 = fNuclearRef[i2];
+ // check if any of p[i1] or p[i2] was a Code::Nucleus
+ const int ref1 = fNucleusRef[i1];
+ const int ref2 = fNucleusRef[i2];
if (ref2 < 0) {
if (ref1 >= 0) {
// i1 is nucleus, i2 is not
- fNuclearRef[i2] = GetNucleusNextRef();
- fNuclearA[fNuclearRef[i2]] = fNuclearA[ref1];
- fNuclearZ[fNuclearRef[i2]] = fNuclearZ[ref1];
+ fNucleusRef[i2] = GetNucleusNextRef();
+ fNuclearA[fNucleusRef[i2]] = fNuclearA[ref1];
+ fNuclearZ[fNucleusRef[i2]] = fNuclearZ[ref1];
} else {
// neither i1 nor i2 are nuclei
}
} else {
if (ref1 >= 0) {
// both are nuclei, i2 is overwritten with nucleus i1
+ // fNucleusRef stays the same, but A and Z data is overwritten
fNuclearA[ref2] = fNuclearA[ref1];
fNuclearZ[ref2] = fNuclearZ[ref1];
} else {
// i2 is overwritten with non-nucleus i1
- fNuclearA.erase(fNuclearA.cbegin() + ref2);
- fNuclearZ.erase(fNuclearZ.cbegin() + ref2);
- const int n = fNuclearRef.size();
+ fNucleusRef[i2] = -1; // flag as non-nucleus
+ fNuclearA.erase(fNuclearA.cbegin() + ref2); // remove data for i2
+ fNuclearZ.erase(fNuclearZ.cbegin() + ref2); // remove data for i2
+ const int n = fNucleusRef.size(); // update fNucleusRef: indices above ref2
+ // must be decremented by 1
for (int i = 0; i < n; ++i) {
- if (fNuclearRef[i] >= ref2) { fNuclearRef[i] -= 1; }
+ if (fNucleusRef[i] > ref2) { fNucleusRef[i] -= 1; }
}
}
}
@@ -192,48 +224,34 @@ namespace corsika::stack {
* Function to copy particle at location i2 in stack to i1
*/
void Swap(const unsigned int i1, const unsigned int i2) {
+ // index range check
if (i1 >= GetSize() || i2 >= GetSize()) {
std::ostringstream err;
err << "NuclearStackExtension: trying to access data beyond size of stack!";
throw std::runtime_error(err.str());
}
+ // swap original particle data
InnerStackImpl::Swap(i1, i2);
- const int ref1 = fNuclearRef[i1];
- const int ref2 = fNuclearRef[i2];
- if (ref2 < 0) {
- if (ref1 >= 0) {
- // i1 is nucleus, i2 is not
- std::swap(fNuclearRef[i2], fNuclearRef[i1]);
- } else {
- // neither i1 nor i2 are nuclei
- }
- } else {
- if (ref1 >= 0) {
- // both are nuclei, i2 is overwritten with nucleus i1
- std::swap(fNuclearRef[i2], fNuclearRef[i1]);
- } else {
- // i2 is overwritten with non-nucleus i1
- std::swap(fNuclearRef[i2], fNuclearRef[i1]);
- }
- }
+ // swap corresponding nuclear reference data
+ std::swap(fNucleusRef[i2], fNucleusRef[i1]);
}
void IncrementSize() {
InnerStackImpl::IncrementSize();
- fNuclearRef.push_back(-1);
+ fNucleusRef.push_back(-1);
}
void DecrementSize() {
InnerStackImpl::DecrementSize();
- if (fNuclearRef.size() > 0) {
- const int ref = fNuclearRef.back();
- fNuclearRef.pop_back();
+ if (fNucleusRef.size() > 0) {
+ const int ref = fNucleusRef.back();
+ fNucleusRef.pop_back();
if (ref >= 0) {
fNuclearA.erase(fNuclearA.begin() + ref);
fNuclearZ.erase(fNuclearZ.begin() + ref);
- const int n = fNuclearRef.size();
+ const int n = fNucleusRef.size();
for (int i = 0; i < n; ++i) {
- if (fNuclearRef[i] >= ref) { fNuclearRef[i] -= 1; }
+ if (fNucleusRef[i] >= ref) { fNucleusRef[i] -= 1; }
}
}
}
@@ -242,7 +260,7 @@ namespace corsika::stack {
private:
/// the actual memory to store particle data
- std::vector<int> fNuclearRef;
+ std::vector<int> fNucleusRef;
std::vector<unsigned short> fNuclearA;
std::vector<unsigned short> fNuclearZ;
diff --git a/Stack/SuperStupidStack/SuperStupidStack.h b/Stack/SuperStupidStack/SuperStupidStack.h
index 61ef97a96d3a3305156894f13cfd69b1f59b5a0e..cb63b03360b9df5c4c5593468258de900218b1bb 100644
--- a/Stack/SuperStupidStack/SuperStupidStack.h
+++ b/Stack/SuperStupidStack/SuperStupidStack.h
@@ -114,6 +114,7 @@ namespace corsika::stack {
public:
void Init() {}
+ void Dump() const {}
void Clear() {
fDataPID.clear();