IAP GITLAB

Commit 13a9d066 authored by Ralf Ulrich's avatar Ralf Ulrich

more robust Stack::ParticleType, preventing dangling objects and crashes

parent 2a549e44
......@@ -201,7 +201,7 @@ namespace corsika::cascade {
actual_inv_length);
const auto sample_process = uniDist(fRNG);
InverseGrammageType inv_lambda_count = 0. * meter * meter / gram;
fProcessSequence.SelectInteraction(particle, fStack, sample_process,
fProcessSequence.SelectInteraction(particle, step, fStack, sample_process,
inv_lambda_count);
} else {
std::cout << "decay" << std::endl;
......
......@@ -132,24 +132,24 @@ namespace corsika::process {
return tot;
}
template <typename Particle, typename Stack>
template <typename Particle, typename Track, typename Stack>
EProcessReturn SelectInteraction(
Particle& p, Stack& s,
Particle& vP, Track& vT, Stack& vS,
[[maybe_unused]] corsika::units::si::InverseGrammageType lambda_select,
corsika::units::si::InverseGrammageType& lambda_inv_count) {
if constexpr (is_process_sequence<T1type>::value) {
// if A is a process sequence --> check inside
const EProcessReturn ret =
A.SelectInteraction(p, s, lambda_select, lambda_inv_count);
A.SelectInteraction(vP, vT, vS, lambda_select, lambda_inv_count);
// if A did succeed, stop routine
if (ret != EProcessReturn::eOk) { return ret; }
} else if constexpr (std::is_base_of<InteractionProcess<T1type>, T1type>::value) {
// if this is not a ContinuousProcess --> evaluate probability
lambda_inv_count += A.GetInverseInteractionLength(p, s);
lambda_inv_count += A.GetInverseInteractionLength(vP, vT);
// check if we should execute THIS process and then EXIT
if (lambda_select < lambda_inv_count) {
A.DoInteraction(p, s);
A.DoInteraction(vP, vS);
return EProcessReturn::eInteracted;
}
} // end branch A
......@@ -157,15 +157,15 @@ namespace corsika::process {
if constexpr (is_process_sequence<T2>::value) {
// if A is a process sequence --> check inside
const EProcessReturn ret =
B.SelectInteraction(p, s, lambda_select, lambda_inv_count);
B.SelectInteraction(vP, vT, vS, lambda_select, lambda_inv_count);
// if A did succeed, stop routine
if (ret != EProcessReturn::eOk) { return ret; }
} else if constexpr (std::is_base_of<InteractionProcess<T2type>, T2type>::value) {
// if this is not a ContinuousProcess --> evaluate probability
lambda_inv_count += B.GetInverseInteractionLength(p, s);
lambda_inv_count += B.GetInverseInteractionLength(vP, vT);
// check if we should execute THIS process and then EXIT
if (lambda_select < lambda_inv_count) {
B.DoInteraction(p, s);
B.DoInteraction(vP, vS);
return EProcessReturn::eInteracted;
}
} // end branch A
......
......@@ -48,6 +48,7 @@ namespace corsika::stack {
class ParticleBase {
public:
using StackIteratorType = StackIterator;
ParticleBase() = default;
private:
......
......@@ -19,22 +19,6 @@
#include <corsika/stack/SecondaryView.h>
// SFINAE test
template <typename T>
class HasGetIndexFromIterator {
private:
typedef char YesType[1];
typedef char NoType[2];
template <typename C>
static YesType& test(decltype(&C::GetIndexFromIterator));
template <typename C>
static NoType& test(...);
public:
enum { value = sizeof(test<T>(0)) == sizeof(YesType) };
};
/**
All classes around management of particles on a stack.
*/
......@@ -123,9 +107,14 @@ namespace corsika::stack {
ParticleInterface, StackType>;
/**
* this is the full type of the declared ParticleInterface: typedef typename
* this is the full type of the user-declared ParticleInterface
*/
using ParticleInterfaceType = typename StackIterator::ParticleInterfaceType;
/**
* In all programming context, the object to access, copy, and
* transport particle data is via the StackIterator
*/
typedef typename StackIterator::ParticleInterfaceType ParticleType;
using ParticleType = StackIterator;
// friends are needed since they need access to protected members
friend class StackIteratorInterface<
......@@ -216,7 +205,7 @@ namespace corsika::stack {
/**
* delete this particle
*/
void Delete(ParticleType p) { Delete(p.GetIterator()); }
void Delete(ParticleInterfaceType p) { Delete(p.GetIterator()); }
/**
* delete last particle on stack by decrementing stack size
......
......@@ -89,7 +89,6 @@ using CombinedTestInterfaceType =
TestParticleInterface2, StackIter>;
using StackTest = CombinedStack<TestStackData, TestStackData2, CombinedTestInterfaceType>;
typedef StackTest::ParticleType Particle;
TEST_CASE("Combined Stack", "[stack]") {
......@@ -161,10 +160,9 @@ TEST_CASE("Combined Stack", "[stack]") {
StackTest s;
REQUIRE(s.GetSize() == 0);
auto iter = s.AddParticle(std::tuple{9.9});
Particle& p = *iter; // also this is valid to access particle data
p.SetData2(2);
iter.SetData2(2);
REQUIRE(s.GetSize() == 1);
p.AddSecondary(std::tuple{4.4});
iter.AddSecondary(std::tuple{4.4});
REQUIRE(s.GetSize() == 2);
// p.AddSecondary(3.3, 2.2, 1.);
// REQUIRE(s.GetSize() == 3);
......@@ -259,7 +257,6 @@ using CombinedTestInterfaceType2 =
using StackTest2 = CombinedStack<typename StackTest::StackImpl, TestStackData3,
CombinedTestInterfaceType2>;
typedef StackTest2::ParticleType Particle2;
TEST_CASE("Combined Stack - multi", "[stack]") {
......
......@@ -32,7 +32,6 @@ using namespace corsika::stack;
using namespace std;
typedef Stack<TestStackData, TestParticleInterface> StackTest;
typedef StackTest::ParticleType Particle;
TEST_CASE("SecondaryStack", "[stack]") {
......
......@@ -33,7 +33,7 @@ using namespace corsika::stack;
using namespace std;
typedef Stack<TestStackData, TestParticleInterface> StackTest;
typedef StackTest::ParticleType Particle;
typedef StackTest::ParticleInterfaceType Particle;
TEST_CASE("Stack", "[Stack]") {
......
......@@ -18,6 +18,9 @@ add_custom_command (
set (
MODEL_SOURCES
ParticleConversion.cc
Interaction.cc
Decay.cc
NuclearInteraction.cc
sibyll2.3c.f
nuclib.f
signuc.f
......
/*
* (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/sibyll/Decay.h>
#include <corsika/process/sibyll/ParticleConversion.h>
#include <corsika/process/sibyll/SibStack.h>
#include <corsika/setup/SetupStack.h>
#include <corsika/setup/SetupTrajectory.h>
using std::cout;
using std::endl;
using std::tuple;
using std::vector;
using namespace corsika;
using namespace corsika::setup;
using Particle = Stack::StackIterator; // ParticleType;
using Track = Trajectory;
namespace corsika::process::sibyll {
Decay::Decay() {}
Decay::~Decay() { cout << "Sibyll::Decay n=" << fCount << endl; }
void Decay::Init() {
setHadronsUnstable();
setTrackedParticlesStable();
}
void Decay::setTrackedParticlesStable() {
/*
Sibyll is hadronic generator
only hadrons decay
*/
// set particles unstable
setHadronsUnstable();
// make tracked particles stable
cout << "Interaction: setting tracked hadrons stable.." << endl;
const vector<particles::Code> particleList = {
particles::Code::PiPlus, particles::Code::PiMinus, particles::Code::KPlus,
particles::Code::KMinus, particles::Code::K0Long, particles::Code::K0Short};
for (auto p : particleList) {
// set particle stable by setting table value negative
const int sibid = process::sibyll::ConvertToSibyllRaw(p);
s_csydec_.idb[sibid - 1] = (-1) * abs(s_csydec_.idb[sibid - 1]);
}
}
void Decay::setUnstable(const particles::Code pCode) {
int s_id = process::sibyll::ConvertToSibyllRaw(pCode);
s_csydec_.idb[s_id - 1] = abs(s_csydec_.idb[s_id - 1]);
}
void Decay::setStable(const particles::Code pCode) {
int s_id = process::sibyll::ConvertToSibyllRaw(pCode);
s_csydec_.idb[s_id - 1] = (-1) * abs(s_csydec_.idb[s_id - 1]);
}
void Decay::setAllStable() {
// name? also makes EM particles stable
cout << "Decay: setting all particles stable.." << endl;
// loop over all particles in sibyll
// should be changed to loop over human readable list
// i.e. particles::ListOfParticles()
for (auto& p : corsika2sibyll) {
// cout << (int)p << endl;
const int sibCode = static_cast<int>(p);
// skip unknown and antiparticles
if (sibCode < 1) continue;
s_csydec_.idb[sibCode - 1] = -1 * abs(s_csydec_.idb[sibCode - 1]);
}
}
void Decay::setHadronsUnstable() {
// name? also makes EM particles stable
// loop over all particles in sibyll
// should be changed to loop over human readable list
// i.e. particles::ListOfParticles()
cout << "Sibyll: setting hadrons unstable.." << endl;
// make ALL particles unstable, then set EM stable
for (int sibCode : corsika2sibyll) {
if (sibCode < 1) continue;
s_csydec_.idb[sibCode - 1] = abs(s_csydec_.idb[sibCode - 1]);
}
// set Leptons and Proton and Neutron stable
// use stack to loop over particles
constexpr particles::Code particleList[] = {
particles::Code::Proton, particles::Code::Neutron, particles::Code::Electron,
particles::Code::Positron, particles::Code::NuE, particles::Code::NuEBar,
particles::Code::MuMinus, particles::Code::MuPlus, particles::Code::NuMu,
particles::Code::NuMuBar};
for (auto p : particleList) {
const int sibid = process::sibyll::ConvertToSibyllRaw(p);
s_csydec_.idb[sibid - 1] = (-1) * abs(s_csydec_.idb[sibid - 1]);
}
}
template <>
units::si::TimeType Decay::GetLifetime(Particle const& p) {
using namespace units::si;
HEPEnergyType E = p.GetEnergy();
HEPMassType m = particles::GetMass(p.GetPID());
const double gamma = E / m;
const TimeType t0 = 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;
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;
return lifetime;
}
template <>
void Decay::DoDecay(Particle& p, Stack&) {
using geometry::Point;
using namespace units::si;
fCount++;
SibStack ss;
ss.Clear();
const particles::Code pCode = p.GetPID();
// copy particle to sibyll stack
ss.AddParticle(process::sibyll::ConvertToSibyllRaw(pCode), p.GetEnergy(),
p.GetMomentum(),
// setting particle mass with Corsika values, may be inconsistent
// with sibyll internal values
particles::GetMass(pCode));
// remember position
Point const decayPoint = p.GetPosition();
TimeType const t0 = p.GetTime();
// remove original particle from corsika stack
p.Delete();
// set all particles/hadrons unstable
// setHadronsUnstable();
setUnstable(pCode);
// call sibyll decay
cout << "Decay: calling Sibyll decay routine.." << endl;
decsib_();
// reset to stable
setStable(pCode);
// print output
int print_unit = 6;
sib_list_(print_unit);
// copy particles from sibyll stack to corsika
for (auto& psib : ss) {
// FOR NOW: skip particles that have decayed in Sibyll, move to iterator?
if (psib.HasDecayed()) continue;
// add to corsika stack
p.AddSecondary(
tuple<particles::Code, units::si::HEPEnergyType, corsika::stack::MomentumVector,
geometry::Point, units::si::TimeType>{
process::sibyll::ConvertFromSibyll(psib.GetPID()), psib.GetEnergy(),
psib.GetMomentum(), decayPoint, t0});
}
// empty sibyll stack
ss.Clear();
}
} // namespace corsika::process::sibyll
......@@ -12,14 +12,8 @@
#ifndef _include_corsika_process_sibyll_decay_h_
#define _include_corsika_process_sibyll_decay_h_
#include <corsika/process/DecayProcess.h>
#include <corsika/process/sibyll/ParticleConversion.h>
#include <corsika/process/sibyll/SibStack.h>
#include <corsika/setup/SetupStack.h>
#include <corsika/setup/SetupTrajectory.h>
#include <corsika/particles/ParticleProperties.h>
#include <corsika/process/DecayProcess.h>
namespace corsika::process {
......@@ -29,191 +23,21 @@ namespace corsika::process {
int fCount = 0;
public:
Decay() {}
~Decay() { std::cout << "Sibyll::Decay n=" << fCount << std::endl; }
void Init() {
setHadronsUnstable();
setTrackedParticlesStable();
}
void setTrackedParticlesStable() {
/*
Sibyll is hadronic generator
only hadrons decay
*/
// set particles unstable
setHadronsUnstable();
// make tracked particles stable
std::cout << "Interaction: setting tracked hadrons stable.." << std::endl;
const std::vector<corsika::particles::Code> particleList = {
particles::Code::PiPlus, particles::Code::PiMinus, particles::Code::KPlus,
particles::Code::KMinus, particles::Code::K0Long, particles::Code::K0Short};
for (auto p : particleList) {
// set particle stable by setting table value negative
const int sibid = process::sibyll::ConvertToSibyllRaw(p);
s_csydec_.idb[sibid - 1] = (-1) * abs(s_csydec_.idb[sibid - 1]);
}
}
void setUnstable(const corsika::particles::Code pCode) {
int s_id = process::sibyll::ConvertToSibyllRaw(pCode);
s_csydec_.idb[s_id - 1] = abs(s_csydec_.idb[s_id - 1]);
}
void setStable(const corsika::particles::Code pCode) {
int s_id = process::sibyll::ConvertToSibyllRaw(pCode);
s_csydec_.idb[s_id - 1] = (-1) * abs(s_csydec_.idb[s_id - 1]);
}
void setAllStable() {
// name? also makes EM particles stable
using std::cout;
using std::endl;
std::cout << "Decay: setting all particles stable.." << std::endl;
// loop over all particles in sibyll
// should be changed to loop over human readable list
// i.e. corsika::particles::ListOfParticles()
for (auto& p : corsika2sibyll) {
// std::cout << (int)p << std::endl;
const int sibCode = static_cast<int>(p);
// skip unknown and antiparticles
if (sibCode < 1) continue;
s_csydec_.idb[sibCode - 1] = -1 * abs(s_csydec_.idb[sibCode - 1]);
}
}
Decay();
~Decay();
void Init();
void setHadronsUnstable() {
using std::cout;
using std::endl;
using namespace corsika::units::si;
// name? also makes EM particles stable
// loop over all particles in sibyll
// should be changed to loop over human readable list
// i.e. corsika::particles::ListOfParticles()
std::cout << "Sibyll: setting hadrons unstable.." << std::endl;
// make ALL particles unstable, then set EM stable
for (int sibCode : corsika2sibyll) {
// std::cout << (int)p << std::endl;
// skip unknown and antiparticles
if (sibCode < 1) continue;
// std::cout << "Sibyll: Decay: setting " << ConvertFromSibyll(
// static_cast<SibyllCode> ( sibCode ) ) << " unstable" << std::endl;
s_csydec_.idb[sibCode - 1] = abs(s_csydec_.idb[sibCode - 1]);
// std::cout << "decay table value: " << s_csydec_.idb[ sibCode - 1 ] <<
// std::endl;
}
// set Leptons and Proton and Neutron stable
// use stack to loop over particles
constexpr corsika::particles::Code particleList[] = {
corsika::particles::Code::Proton, corsika::particles::Code::Neutron,
corsika::particles::Code::Electron, corsika::particles::Code::Positron,
corsika::particles::Code::NuE, corsika::particles::Code::NuEBar,
corsika::particles::Code::MuMinus, corsika::particles::Code::MuPlus,
corsika::particles::Code::NuMu, corsika::particles::Code::NuMuBar};
for (auto p : particleList) {
// set particle stable by setting table value negative
// cout << "Sibyll: setting " << p.GetPID() << "(" << s_id << ")"
// << " stable in Sibyll .." << endl;
const int sibid = process::sibyll::ConvertToSibyllRaw(p);
s_csydec_.idb[sibid - 1] = (-1) * abs(s_csydec_.idb[sibid - 1]);
}
}
void setTrackedParticlesStable();
void setUnstable(const corsika::particles::Code pCode);
void setStable(const corsika::particles::Code pCode);
void setAllStable();
void setHadronsUnstable();
template <typename Particle>
corsika::units::si::TimeType GetLifetime(Particle const& p) {
using std::cout;
using std::endl;
using namespace corsika::units::si;
corsika::units::si::HEPEnergyType E = p.GetEnergy();
corsika::units::si::HEPMassType m = corsika::particles::GetMass(p.GetPID());
const double gamma = E / m;
const corsika::units::si::TimeType t0 =
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;
return lifetime;
}
corsika::units::si::TimeType GetLifetime(Particle const& p);
template <typename Particle, typename Stack>
void DoDecay(Particle& p, Stack&) {
using corsika::geometry::Point;
using namespace corsika::units::si;
fCount++;
SibStack ss;
ss.Clear();
const corsika::particles::Code pCode = p.GetPID();
// copy particle to sibyll stack
ss.AddParticle(
// std::tuple<corsika::particles::Code, corsika::units::si::HEPEnergyType,
// corsika::stack::MomentumVector, corsika::geometry::Point,
// corsika::units::si::TimeType>{
corsika::process::sibyll::ConvertToSibyllRaw(pCode), p.GetEnergy(),
p.GetMomentum(),
// setting particle mass with Corsika values, may be inconsistent
// with sibyll internal values
// TODO: #warning setting particle mass with Corsika values, may be
// inconsistent with sibyll internal values
corsika::particles::GetMass(pCode));
// remember position
Point const decayPoint = p.GetPosition();
TimeType const t0 = p.GetTime();
// remove original particle from corsika stack
p.Delete();
// set all particles/hadrons unstable
// setHadronsUnstable();
setUnstable(pCode);
// call sibyll decay
std::cout << "Decay: calling Sibyll decay routine.." << std::endl;
decsib_();
// reset to stable
setStable(pCode);
// print output
int print_unit = 6;
sib_list_(print_unit);
// copy particles from sibyll stack to corsika
for (auto& psib : ss) {
// FOR NOW: skip particles that have decayed in Sibyll, move to iterator?
if (psib.HasDecayed()) continue;
// add to corsika stack
p.AddSecondary(
std::tuple<corsika::particles::Code, corsika::units::si::HEPEnergyType,
corsika::stack::MomentumVector, corsika::geometry::Point,
corsika::units::si::TimeType>{
corsika::process::sibyll::ConvertFromSibyll(psib.GetPID()),
psib.GetEnergy(), psib.GetMomentum(), decayPoint, t0});
}
// empty sibyll stack
ss.Clear();
}
void DoDecay(Particle& p, Stack&);
};
} // namespace sibyll
} // namespace corsika::process
......
This diff is collapsed.
This diff is collapsed.
This diff is collapsed.
This diff is collapsed.
/*
* (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/tracking_line/TrackingLine.h>
#include <corsika/environment/Environment.h>
......@@ -126,8 +137,6 @@ namespace corsika::process::tracking_line {
#include <corsika/setup/SetupStack.h>
#include <corsika/setup/SetupTrajectory.h>
using namespace corsika::setup;
using Particle = Stack::ParticleType;
// using Track = Trajectory;
template class corsika::process::tracking_line::TrackingLine<setup::Stack,
setup::Trajectory>;
......
......@@ -31,7 +31,7 @@ namespace corsika::process {
template <typename Stack, typename Trajectory>
class TrackingLine { //
using Particle = typename Stack::ParticleType;
using Particle = typename Stack::StackIterator;
corsika::environment::Environment const& fEnvironment;
......
/*
* (c) Copyright 2018 CORSIKA Project, corsika-project@lists.kit.edu
*
* See file AUTHORS for a list of contributors.
*