From 172cc93b2ef4e46e76c6955a77a4594ad1c4200b Mon Sep 17 00:00:00 2001
From: Felix Riehn <friehn@lip.pt>
Date: Tue, 5 Feb 2019 16:46:47 +0100
Subject: [PATCH] Resolve "use nuclear stack extension in nuclear interaction"
---
.gitlab-ci.yml | 4 +-
Documentation/Examples/cascade_example.cc | 36 ++-
.../StackInterface/testStackInterface.cc | 1 +
Processes/Sibyll/Interaction.h | 6 +
Processes/Sibyll/NuclearInteraction.h | 228 ++++++++++--------
Processes/Sibyll/SibStack.h | 4 +-
Processes/Sibyll/testSibyll.cc | 4 +-
Processes/StackInspector/StackInspector.cc | 5 +-
.../NuclearStackExtension.h | 107 ++++----
Stack/SuperStupidStack/SuperStupidStack.h | 1 +
10 files changed, 236 insertions(+), 160 deletions(-)
diff --git a/.gitlab-ci.yml b/.gitlab-ci.yml
index 879e1ac5..c83d5488 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 772853db..67b6f4d9 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/StackInterface/testStackInterface.cc b/Framework/StackInterface/testStackInterface.cc
index 293aa591..83ab16e1 100644
--- a/Framework/StackInterface/testStackInterface.cc
+++ b/Framework/StackInterface/testStackInterface.cc
@@ -38,6 +38,7 @@ class TestStackData {
public:
// these functions are needed for the Stack interface
void Init() {}
+ void Dump() const {}
void Clear() { fData.clear(); }
unsigned int GetSize() const { return fData.size(); }
unsigned int GetCapacity() const { return fData.size(); }
diff --git a/Processes/Sibyll/Interaction.h b/Processes/Sibyll/Interaction.h
index 966c2f7f..69755f00 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 50384dde..7b7492d4 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 9b168f6e..fd847124 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 6251aab1..8f53c6e1 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 312fe40e..fb749c0b 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 a08dad1d..355e9821 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,28 +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); }
};
/**
- * Memory implementation of the most simple (stupid) particle stack object.
+ * @class NuclearStackExtension
+ *
+ * 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;
@@ -128,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)]; }
@@ -141,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());
@@ -151,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; }
}
}
}
@@ -189,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; }
}
}
}
@@ -239,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 61ef97a9..cb63b033 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();
--
GitLab