IAP GITLAB

Skip to content
Snippets Groups Projects
Commit afeb2e86 authored by Maximilian Reininghaus's avatar Maximilian Reininghaus :vulcan: Committed by ralfulrich
Browse files

Stack and Sibyll

parent 6a460fac
No related branches found
No related tags found
1 merge request!280Refactory 2020
Hide whitespace changes
Inline Side-by-side
Showing
with 374 additions and 410 deletions
corsika/detail/modules/sibyll/Decay.inl
+ 2
− 2
View file @ afeb2e86
...@@ -90,7 +90,7 @@ namespace corsika::sibyll { ...@@ -90,7 +90,7 @@ namespace corsika::sibyll {
const double gamma = E / m; const double gamma = E / m;
const TimeType t0 = corsika::GetLifetime(vP.GetPID()); const TimeType t0 = corsika::lifetime(vP.GetPID());
auto const lifetime = gamma * t0; auto const lifetime = gamma * t0;
const auto mkin = const auto mkin =
...@@ -123,7 +123,7 @@ namespace corsika::sibyll { ...@@ -123,7 +123,7 @@ namespace corsika::sibyll {
vP.GetMomentum(), vP.GetMomentum(),
// setting particle mass with Corsika values, may be inconsistent // setting particle mass with Corsika values, may be inconsistent
// with sibyll internal values // with sibyll internal values
corsika::GetMass(pCode)); corsika::mass(pCode));
// remember position // remember position
Point const decayPoint = vP.GetPosition(); Point const decayPoint = vP.GetPosition();
TimeType const t0 = vP.GetTime(); TimeType const t0 = vP.GetTime();
......
corsika/detail/modules/sibyll/Interaction.inl
+ 6
− 6
View file @ afeb2e86
...@@ -87,13 +87,13 @@ namespace corsika::sibyll { ...@@ -87,13 +87,13 @@ namespace corsika::sibyll {
"Interaction: GetCrossSection: CoM energy outside range for Sibyll!"); "Interaction: GetCrossSection: CoM energy outside range for Sibyll!");
} }
const double dEcm = CoMenergy / 1_GeV; const double dEcm = CoMenergy / 1_GeV;
if (corsika::IsNucleus(TargetId)) { if (corsika::is_nucleus(TargetId)) {
const int iTarget = corsika::GetNucleusA(TargetId); const int iTarget = corsika::nucleus_A(TargetId);
if (iTarget > maxTargetMassNumber_ || iTarget == 0) if (iTarget > maxTargetMassNumber_ || iTarget == 0)
throw std::runtime_error( throw std::runtime_error(
"Sibyll target outside range. Only nuclei with A<18 are allowed."); "Sibyll target outside range. Only nuclei with A<18 are allowed.");
sib_sigma_hnuc_(iBeam, iTarget, dEcm, sigProd, dummy, sigEla); sib_sigma_hnuc_(iBeam, iTarget, dEcm, sigProd, dummy, sigEla);
} else if (TargetId == corsika::Proton::GetCode()) { } else if (TargetId == corsika::Code::Proton) {
sib_sigma_hp_(iBeam, dEcm, dum1, sigEla, sigProd, dumdif, dum3, dum4); sib_sigma_hp_(iBeam, dEcm, dum1, sigEla, sigProd, dumdif, dum3, dum4);
} else { } else {
// no interaction in sibyll possible, return infinite cross section? or throw? // no interaction in sibyll possible, return infinite cross section? or throw?
...@@ -186,7 +186,7 @@ namespace corsika::sibyll { ...@@ -186,7 +186,7 @@ namespace corsika::sibyll {
<< "DoInteraction: " << corsikaBeamId << " interaction? " << "DoInteraction: " << corsikaBeamId << " interaction? "
<< corsika::sibyll::CanInteract(corsikaBeamId) << std::endl; << corsika::sibyll::CanInteract(corsikaBeamId) << std::endl;
if (corsika::IsNucleus(corsikaBeamId)) { if (corsika::is_nucleus(corsikaBeamId)) {
// nuclei handled by different process, this should not happen // nuclei handled by different process, this should not happen
throw std::runtime_error("Nuclear projectile are not handled by SIBYLL!"); throw std::runtime_error("Nuclear projectile are not handled by SIBYLL!");
} }
...@@ -276,8 +276,8 @@ namespace corsika::sibyll { ...@@ -276,8 +276,8 @@ namespace corsika::sibyll {
allowed air in atmosphere also contains some Argon. allowed air in atmosphere also contains some Argon.
*/ */
int targetSibCode = -1; int targetSibCode = -1;
if (IsNucleus(targetCode)) targetSibCode = GetNucleusA(targetCode); if (is_nucleus(targetCode)) targetSibCode = nucleus_A(targetCode);
if (targetCode == corsika::Proton::GetCode()) targetSibCode = 1; if (targetCode == corsika::Code::Proton) targetSibCode = 1;
std::cout << "Interaction: sibyll code: " << targetSibCode << std::endl; std::cout << "Interaction: sibyll code: " << targetSibCode << std::endl;
if (targetSibCode > maxTargetMassNumber_ || targetSibCode < 1) if (targetSibCode > maxTargetMassNumber_ || targetSibCode < 1)
throw std::runtime_error( throw std::runtime_error(
......
corsika/detail/modules/sibyll/NuclearInteraction.inl
+ 11
− 16
View file @ afeb2e86
...@@ -50,7 +50,7 @@ namespace corsika::sibyll { ...@@ -50,7 +50,7 @@ namespace corsika::sibyll {
for (unsigned int i = 0; i < GetNEnergyBins(); ++i) { for (unsigned int i = 0; i < GetNEnergyBins(); ++i) {
std::cout << " " << i << " "; std::cout << " " << i << " ";
for (auto& n : pNuclei) { for (auto& n : pNuclei) {
auto const j = GetNucleusA(n); auto const j = corsika::nucleus_A(n);
std::cout << " " << std::setprecision(5) << std::setw(8) std::cout << " " << std::setprecision(5) << std::setw(8)
<< cnucsignuc_.sigma[j - 1][k][i]; << cnucsignuc_.sigma[j - 1][k][i];
} }
...@@ -85,7 +85,7 @@ namespace corsika::sibyll { ...@@ -85,7 +85,7 @@ namespace corsika::sibyll {
for (auto& ptarg : allElementsInUniverse) { for (auto& ptarg : allElementsInUniverse) {
++k; ++k;
std::cout << "NuclearInteraction: init target component: " << ptarg << std::endl; std::cout << "NuclearInteraction: init target component: " << ptarg << std::endl;
const int ib = GetNucleusA(ptarg); const int ib = corsika::nucleus_A(ptarg);
if (!hadronicInteraction_.IsValidTarget(ptarg)) { if (!hadronicInteraction_.IsValidTarget(ptarg)) {
std::cout std::cout
<< "NuclearInteraction::InitializeNuclearCrossSections: target nucleus? id=" << "NuclearInteraction::InitializeNuclearCrossSections: target nucleus? id="
...@@ -217,7 +217,7 @@ namespace corsika::sibyll { ...@@ -217,7 +217,7 @@ namespace corsika::sibyll {
if (corsikaBeamId != corsika::Code::Nucleus) { if (corsikaBeamId != corsika::Code::Nucleus) {
// check if target-style nucleus (enum), these are not allowed as projectile // check if target-style nucleus (enum), these are not allowed as projectile
if (corsika::IsNucleus(corsikaBeamId)) if (corsika::is_nucleus(corsikaBeamId))
throw std::runtime_error( throw std::runtime_error(
"NuclearInteraction: GetInteractionLength: Wrong nucleus type. Nuclear " "NuclearInteraction: GetInteractionLength: Wrong nucleus type. Nuclear "
"projectiles should use NuclearStackExtension!"); "projectiles should use NuclearStackExtension!");
...@@ -326,9 +326,7 @@ namespace corsika::sibyll { ...@@ -326,9 +326,7 @@ namespace corsika::sibyll {
"NuclearInteraction: DoInteraction: Wrong nucleus type. Nuclear projectiles " "NuclearInteraction: DoInteraction: Wrong nucleus type. Nuclear projectiles "
"should use NuclearStackExtension!"); "should use NuclearStackExtension!");
auto const ProjMass = auto const ProjMass = corsika::nucleus_mass(vP.GetNuclearA(), vP.GetNuclearZ());
vP.GetNuclearZ() * corsika::Proton::GetMass() +
(vP.GetNuclearA() - vP.GetNuclearZ()) * corsika::Neutron::GetMass();
std::cout << "NuclearInteraction: projectile mass: " << ProjMass / 1_GeV << std::endl; std::cout << "NuclearInteraction: projectile mass: " << ProjMass / 1_GeV << std::endl;
count_++; count_++;
...@@ -414,7 +412,7 @@ namespace corsika::sibyll { ...@@ -414,7 +412,7 @@ namespace corsika::sibyll {
// sample target nucleon number // sample target nucleon number
// //
// proton stand-in for nucleon // proton stand-in for nucleon
const auto beamId = corsika::Proton::GetCode(); const auto beamId = corsika::Code::Proton;
auto const* const currentNode = vP.GetNode(); auto const* const currentNode = vP.GetNode();
const auto& mediumComposition = const auto& mediumComposition =
currentNode->GetModelProperties().GetNuclearComposition(); currentNode->GetModelProperties().GetNuclearComposition();
...@@ -447,8 +445,8 @@ namespace corsika::sibyll { ...@@ -447,8 +445,8 @@ namespace corsika::sibyll {
allowed air in atmosphere also contains some Argon. allowed air in atmosphere also contains some Argon.
*/ */
int kATarget = -1; int kATarget = -1;
if (IsNucleus(targetCode)) kATarget = GetNucleusA(targetCode); if (corsika::is_nucleus(targetCode)) kATarget = corsika::nucleus_A(targetCode);
if (targetCode == corsika::Proton::GetCode()) kATarget = 1; else if (targetCode == corsika::Code::Proton) kATarget = 1;
std::cout << "NuclearInteraction: nuclib target code: " << kATarget << std::endl; std::cout << "NuclearInteraction: nuclib target code: " << kATarget << std::endl;
if (!hadronicInteraction_.IsValidTarget(targetCode)) if (!hadronicInteraction_.IsValidTarget(targetCode))
throw std::runtime_error("target outside range. "); throw std::runtime_error("target outside range. ");
...@@ -459,7 +457,7 @@ namespace corsika::sibyll { ...@@ -459,7 +457,7 @@ namespace corsika::sibyll {
<< std::endl; << std::endl;
// get nucleon-nucleon cross section // get nucleon-nucleon cross section
// (needed to determine number of nucleon-nucleon scatterings) // (needed to determine number of nucleon-nucleon scatterings)
const auto protonId = corsika::Proton::GetCode(); const auto protonId = corsika::Code::Proton;
const auto [prodCrossSection, elaCrossSection] = const auto [prodCrossSection, elaCrossSection] =
hadronicInteraction_.GetCrossSection(protonId, protonId, EcmNN); hadronicInteraction_.GetCrossSection(protonId, protonId, EcmNN);
const double sigProd = prodCrossSection / 1_mb; const double sigProd = prodCrossSection / 1_mb;
...@@ -519,10 +517,7 @@ namespace corsika::sibyll { ...@@ -519,10 +517,7 @@ namespace corsika::sibyll {
else else
specCode = corsika::Code::Nucleus; specCode = corsika::Code::Nucleus;
// TODO: mass of nuclei? const HEPMassType mass = corsika::nucleus_mass(nuclA, nuclZ);
const HEPMassType mass =
corsika::Proton::GetMass() * nuclZ +
(nuclA - nuclZ) * corsika::Neutron::GetMass(); // this neglects binding energy
std::cout << "NuclearInteraction: adding fragment: " << specCode << std::endl; std::cout << "NuclearInteraction: adding fragment: " << specCode << std::endl;
std::cout << "NuclearInteraction: A,Z: " << nuclA << "," << nuclZ << std::endl; std::cout << "NuclearInteraction: A,Z: " << nuclA << "," << nuclZ << std::endl;
...@@ -565,7 +560,7 @@ namespace corsika::sibyll { ...@@ -565,7 +560,7 @@ namespace corsika::sibyll {
// CORSIKA 7 way // CORSIKA 7 way
// elastic nucleons inherit momentum from original projectile // elastic nucleons inherit momentum from original projectile
// neglecting momentum transfer in interaction // neglecting momentum transfer in interaction
const double mass_ratio = corsika::GetMass(elaNucCode) / ProjMass; const double mass_ratio = corsika::mass(elaNucCode) / ProjMass;
auto const Plab = PprojLab * mass_ratio; auto const Plab = PprojLab * mass_ratio;
vP.AddSecondary(std::tuple<corsika::Code, si::HEPEnergyType, vP.AddSecondary(std::tuple<corsika::Code, si::HEPEnergyType,
...@@ -578,7 +573,7 @@ namespace corsika::sibyll { ...@@ -578,7 +573,7 @@ namespace corsika::sibyll {
std::cout << "calculate inelastic nucleon-nucleon interactions.." << std::endl; std::cout << "calculate inelastic nucleon-nucleon interactions.." << std::endl;
for (int j = 0; j < nInelNucleons; ++j) { for (int j = 0; j < nInelNucleons; ++j) {
// TODO: sample neutron or proton // TODO: sample neutron or proton
auto pCode = corsika::Proton::GetCode(); auto pCode = corsika::Code::Proton;
// temporarily add to stack, will be removed after interaction in DoInteraction // temporarily add to stack, will be removed after interaction in DoInteraction
std::cout << "inelastic interaction no. " << j << std::endl; std::cout << "inelastic interaction no. " << j << std::endl;
auto inelasticNucleon = vP.AddSecondary( auto inelasticNucleon = vP.AddSecondary(
......
corsika/media/NuclearComposition.hpp
+ 3
− 3
View file @ afeb2e86
...@@ -72,10 +72,10 @@ namespace corsika { ...@@ -72,10 +72,10 @@ namespace corsika {
, fAvgMassNumber(std::inner_product( , fAvgMassNumber(std::inner_product(
pComponents.cbegin(), pComponents.cend(), pFractions.cbegin(), 0., pComponents.cbegin(), pComponents.cend(), pFractions.cbegin(), 0.,
std::plus<double>(), [](auto const compID, auto const fraction) -> double { std::plus<double>(), [](auto const compID, auto const fraction) -> double {
if (IsNucleus(compID)) { if (is_nucleus(compID)) {
return GetNucleusA(compID) * fraction; return nucleus_A(compID) * fraction;
} else { } else {
return GetMass(compID) / ConvertSIToHEP(constants::u) * fraction; return mass(compID) / ConvertSIToHEP(constants::u) * fraction;
} }
})) { })) {
assert(pComponents.size() == pFractions.size()); assert(pComponents.size() == pFractions.size());
......
corsika/modules/sibyll/Interaction.hpp
+ 1
− 2
View file @ afeb2e86
...@@ -46,8 +46,7 @@ namespace corsika::sibyll { ...@@ -46,8 +46,7 @@ namespace corsika::sibyll {
HEPEnergyType GetMinEnergyCoM() const { return minEnergyCoM_; } HEPEnergyType GetMinEnergyCoM() const { return minEnergyCoM_; }
HEPEnergyType GetMaxEnergyCoM() const { return maxEnergyCoM_; } HEPEnergyType GetMaxEnergyCoM() const { return maxEnergyCoM_; }
bool IsValidTarget(corsika::Code TargetId) const { bool IsValidTarget(corsika::Code TargetId) const {
return (corsika::GetNucleusA(TargetId) < maxTargetMassNumber_) && return corsika::is_nucleus(TargetId) && (corsika::nucleus_A(TargetId) < maxTargetMassNumber_);
corsika::IsNucleus(TargetId);
} }
std::tuple<CrossSectionType, CrossSectionType> GetCrossSection( std::tuple<CrossSectionType, CrossSectionType> GetCrossSection(
......
corsika/stack/NuclearStackExtension.hpp
+ 338
− 368
View file @ afeb2e86
...@@ -23,372 +23,342 @@ ...@@ -23,372 +23,342 @@
namespace corsika { namespace corsika {
/** /**
* @namespace nuclear_extension * @namespace nuclear_extension
* *
* Add A and Z data to existing stack (currently SuperStupidStack) of particle * Add A and Z data to existing stack of particle properties.
* properties. This is done via inheritance, not via CombinedStack since the nuclear *
* data is stored ONLY when needed (for nuclei) and not for all particles. Thus, this is * Only for Code::Nucleus particles A and Z are stored, not for all
* a new, derived Stack object. * normal elementary particles.
* *
* Only for Code::Nucleus particles A and Z are stored, not for all * Thus in your code, make sure to always check <code>
* normal elementary particles. * particle.GetPID()==Code::Nucleus </code> before attempting to
* * read any nuclear information.
* Thus in your code, make sure to always check <code> *
* particle.GetPID()==Code::Nucleus </code> before attempting to *
* read any nuclear information. */
*
* typedef corsika::Vector<hepmomentum_d> MomentumVector;
*/
namespace nuclear_extension {
/** /**
* @class NuclearParticleInterface * @class NuclearParticleInterface
* *
* Define ParticleInterface for NuclearStackExtension Stack derived from * Define ParticleInterface for NuclearStackExtension Stack derived from
* ParticleInterface of Inner stack class * ParticleInterface of Inner stack class
*/ */
template < template <typename> class InnerParticleInterface, typename StackIteratorInterface> template <template <typename> typename InnerParticleInterface,
struct NuclearParticleInterface : public InnerParticleInterface<StackIteratorInterface> { typename StackIteratorInterface>
class NuclearParticleInterface
typedef InnerParticleInterface<StackIteratorInterface> super_type; : public InnerParticleInterface<StackIteratorInterface> {
protected:
using InnerParticleInterface<StackIteratorInterface>::GetStackData;
public: using InnerParticleInterface<StackIteratorInterface>::GetIndex;
public:
typedef std::tuple< void SetParticleData(
corsika::Code, corsika::units::si::HEPEnergyType, const std::tuple<corsika::Code, HEPEnergyType, corsika::MomentumVector,
momentum_type, corsika::Point, corsika::Point, TimeType>& v) {
corsika::units::si::TimeType> particle_data_type; if (std::get<0>(v) == corsika::Code::Nucleus) {
std::ostringstream err;
typedef std::tuple< err << "NuclearStackExtension: no A and Z specified for new Nucleus!";
corsika::Code, corsika::units::si::HEPEnergyType, throw std::runtime_error(err.str());
momentum_type, corsika::Point, }
corsika::units::si::TimeType, InnerParticleInterface<StackIteratorInterface>::SetParticleData(v);
unsigned short, unsigned short> altenative_particle_data_type; SetNucleusRef(-1); // this is not a nucleus
}
typedef corsika::Vector<corsika::units::si::hepmomentum_d> momentum_type; void SetParticleData(
const std::tuple<corsika::Code, HEPEnergyType, corsika::MomentumVector,
void setParticleData(particle_data_type const& v) { corsika::Point, TimeType, unsigned short, unsigned short>& v) {
const unsigned short A = std::get<5>(v);
if (std::get<0>(v) == corsika::Code::Nucleus) { const unsigned short Z = std::get<6>(v);
std::ostringstream err; if (std::get<0>(v) != corsika::Code::Nucleus || A == 0 || Z == 0) {
err << "NuclearStackExtension: no A and Z specified for new Nucleus!"; std::ostringstream err;
throw std::runtime_error(err.str()); err << "NuclearStackExtension: no A and Z specified for new Nucleus!";
} throw std::runtime_error(err.str());
}
super_type::setParticleData(v); SetNucleusRef(GetStackData().GetNucleusNextRef()); // store this nucleus data ref
setNucleusRef(-1); // this is not a nucleus SetNuclearA(A);
} SetNuclearZ(Z);
InnerParticleInterface<StackIteratorInterface>::SetParticleData(
void setParticleData( altenative_particle_data_type const& v) std::tuple<corsika::Code, HEPEnergyType, corsika::MomentumVector,
{ corsika::Point, TimeType>{std::get<0>(v), std::get<1>(v),
const unsigned short A = std::get<5>(v); std::get<2>(v), std::get<3>(v),
const unsigned short Z = std::get<6>(v); std::get<4>(v)});
if (std::get<0>(v) != corsika::Code::Nucleus || A == 0 || Z == 0) { }
std::ostringstream err;
err << "NuclearStackExtension: no A and Z specified for new Nucleus!"; void SetParticleData(
throw std::runtime_error(err.str()); InnerParticleInterface<StackIteratorInterface>& p,
} const std::tuple<corsika::Code, HEPEnergyType, corsika::MomentumVector,
setNucleusRef(super_type::GetStackData().getNucleusNextRef()); // store this nucleus data ref corsika::Point, TimeType>& v) {
setNuclearA(A); if (std::get<0>(v) == corsika::Code::Nucleus) {
setNuclearZ(Z); std::ostringstream err;
super_type::setParticleData(particle_data_type{std::get<0>(v), std::get<1>(v), err << "NuclearStackExtension: no A and Z specified for new Nucleus!";
std::get<2>(v), std::get<3>(v), std::get<4>(v)}); throw std::runtime_error(err.str());
} }
InnerParticleInterface<StackIteratorInterface>::SetParticleData(
void setParticleData( super_type& p, particle_data_type const& v) p, std::tuple<corsika::Code, HEPEnergyType, corsika::MomentumVector,
{ corsika::Point, TimeType>{std::get<0>(v), std::get<1>(v),
if (std::get<0>(v) == corsika::Code::Nucleus) { std::get<2>(v), std::get<3>(v),
std::ostringstream err; std::get<4>(v)});
err << "NuclearStackExtension: no A and Z specified for new Nucleus!"; SetNucleusRef(-1); // this is not a nucleus
throw std::runtime_error(err.str()); }
}
void SetParticleData(
super_type::setParticleData(p, particle_data_type{std::get<0>(v), std::get<1>(v), InnerParticleInterface<StackIteratorInterface>& p,
std::get<2>(v), std::get<3>(v), std::get<4>(v)}); const std::tuple<corsika::Code, HEPEnergyType, corsika::MomentumVector,
corsika::Point, TimeType, unsigned short, unsigned short>& v) {
setNucleusRef(-1); // this is not a nucleus const unsigned short A = std::get<5>(v);
} const unsigned short Z = std::get<6>(v);
if (std::get<0>(v) != corsika::Code::Nucleus || A == 0 || Z == 0) {
void setParticleData( super_type& p, altenative_particle_data_type const& v) { std::ostringstream err;
err << "NuclearStackExtension: no A and Z specified for new Nucleus!";
const unsigned short A = std::get<5>(v); throw std::runtime_error(err.str());
const unsigned short Z = std::get<6>(v); }
SetNucleusRef(GetStackData().GetNucleusNextRef()); // store this nucleus data ref
if (std::get<0>(v) != corsika::Code::Nucleus || A == 0 || Z == 0) { SetNuclearA(A);
std::ostringstream err; SetNuclearZ(Z);
err << "NuclearStackExtension: no A and Z specified for new Nucleus!"; InnerParticleInterface<StackIteratorInterface>::SetParticleData(
throw std::runtime_error(err.str()); p, std::tuple<corsika::Code, HEPEnergyType, corsika::MomentumVector,
} corsika::Point, TimeType>{std::get<0>(v), std::get<1>(v),
std::get<2>(v), std::get<3>(v),
setNucleusRef(super_type::GetStackData().getNucleusNextRef()); // store this nucleus data ref std::get<4>(v)});
setNuclearA(A); }
setNuclearZ(Z);
super_type::setParticleData(p, particle_data_type{std::get<0>(v), std::get<1>(v), /**
std::get<2>(v), std::get<3>(v), * @name individual setters
std::get<4>(v)}); * @{
} */
void SetNuclearA(const unsigned short vA) {
std::string as_string() const { GetStackData().SetNuclearA(GetIndex(), vA);
return fmt::format( }
"{}, nuc({})", super_type::as_string(), void SetNuclearZ(const unsigned short vZ) {
(isNucleus() ? fmt::format("A={}, Z={}", getNuclearA(), getNuclearZ()) GetStackData().SetNuclearZ(GetIndex(), vZ);
: "n/a")); }
} /// @}
/** /**
* @name individual setters * @name individual getters
* @{ * @{
*/ */
void setNuclearA(const unsigned short vA) { int GetNuclearA() const { return GetStackData().GetNuclearA(GetIndex()); }
super_type::GetStackData().setNuclearA(super_type::GetIndex(), vA); int GetNuclearZ() const { return GetStackData().GetNuclearZ(GetIndex()); }
} /// @}
void setNuclearZ(const unsigned short vZ) {
super_type::GetStackData().setNuclearZ(super_type::GetIndex(), vZ); /**
} * Overwrite normal GetParticleMass function with nuclear version
/// @} */
HEPMassType GetMass() const {
/** if (InnerParticleInterface<StackIteratorInterface>::GetPID() ==
* @name individual getters corsika::Code::Nucleus)
* @{ return corsika::nucleus_mass(GetNuclearA(), GetNuclearZ());
*/ return InnerParticleInterface<StackIteratorInterface>::GetMass();
int getNuclearA() const { return super_type::GetStackData().getNuclearA(super_type::GetIndex()); } }
int getNuclearZ() const { return super_type::GetStackData().getNuclearZ(super_type::GetIndex()); } /**
/// @} * Overwirte normal GetChargeNumber function with nuclear version
**/
/** int16_t GetChargeNumber() const {
* Overwrite normal GetParticleMass function with nuclear version if (InnerParticleInterface<StackIteratorInterface>::GetPID() ==
*/ corsika::Code::Nucleus)
corsika::units::si::HEPMassType getMass() const { return GetNuclearZ();
if (super_type::GetPID() == return InnerParticleInterface<StackIteratorInterface>::GetChargeNumber();
corsika::Code::Nucleus) }
return corsika::GetNucleusMass(getNuclearA(), getNuclearZ());
return super_type::getMass(); int GetNucleusRef() const { return GetStackData().GetNucleusRef(GetIndex()); }
}
/** protected:
* Overwirte normal GetChargeNumber function with nuclear version void SetNucleusRef(const int vR) { GetStackData().SetNucleusRef(GetIndex(), vR); }
**/ };
int16_t getChargeNumber() const {
if (super_type::GetPID() == /**
corsika::Code::Nucleus) * @class NuclearStackExtension
return getNuclearZ(); *
return super_type::getChargeNumber(); * Memory implementation of adding nuclear inforamtion to the
} * existing particle stack defined in class InnerStackImpl.
*
int getNucleusRef() const { * Inside the NuclearStackExtension class there is a dedicated
return super_type::GetStackData().getNucleusRef(GetIndex()); * fNucleusRef index, where fNucleusRef[i] is referring to the
} // LCOV_EXCL_LINE * 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
protected: * GetNucleusA would produce an exception.
*/
void setNucleusRef(const int vR) { template <typename InnerStackImpl>
super_type::GetStackData().setNucleusRef(super_type::GetIndex(), vR); class NuclearStackExtensionImpl : public InnerStackImpl {
}
public:
bool isNucleus() const { void Init() { InnerStackImpl::Init(); }
return super_type::GetStackData().isNucleus(super_type::GetIndex()); void Dump() { InnerStackImpl::Dump(); }
}
}; void Clear() {
InnerStackImpl::Clear();
/** fNucleusRef.clear();
* @class NuclearStackExtension fNuclearA.clear();
* fNuclearZ.clear();
* Memory implementation of adding nuclear inforamtion to the }
* existing particle stack defined in class InnerStackImpl.
* unsigned int GetSize() const { return fNucleusRef.size(); }
* Inside the NuclearStackExtension class there is a dedicated unsigned int GetCapacity() const { return fNucleusRef.capacity(); }
* fNucleusRef index, where fNucleusRef[i] is referring to the
* correct A and Z for a specific particle index i. fNucleusRef[i] void SetNuclearA(const unsigned int i, const unsigned short vA) {
* == -1 means that this is not a nucleus, and a subsequent call to fNuclearA[GetNucleusRef(i)] = vA;
* GetNucleusA would produce an exception. }
*/ void SetNuclearZ(const unsigned int i, const unsigned short vZ) {
template <typename InnerStackImpl> fNuclearZ[GetNucleusRef(i)] = vZ;
class NuclearStackExtensionImpl : public InnerStackImpl { }
void SetNucleusRef(const unsigned int i, const int v) { fNucleusRef[i] = v; }
typedef InnerStackImpl super_type;
int GetNuclearA(const unsigned int i) const { return fNuclearA[GetNucleusRef(i)]; }
public: int GetNuclearZ(const unsigned int i) const { return fNuclearZ[GetNucleusRef(i)]; }
// this function will create new storage for Nuclear Properties, and return the
typedef std::vector<int> nucleus_ref_type; // reference to it
typedef std::vector<unsigned short> nuclear_a_type; int GetNucleusNextRef() {
typedef std::vector<unsigned short> nuclear_z_type; fNuclearA.push_back(0);
fNuclearZ.push_back(0);
return fNuclearA.size() - 1;
NuclearStackExtensionImpl()= default; }
NuclearStackExtensionImpl( NuclearStackExtensionImpl<InnerStackImpl> const&)= default; int GetNucleusRef(const unsigned int i) const {
if (fNucleusRef[i] >= 0) return fNucleusRef[i];
NuclearStackExtensionImpl( NuclearStackExtensionImpl<InnerStackImpl> &&)= default; std::ostringstream err;
err << "NuclearStackExtension: no nucleus at ref=" << i;
NuclearStackExtensionImpl<InnerStackImpl>& throw std::runtime_error(err.str());
operator=( NuclearStackExtensionImpl<InnerStackImpl> const&)= default; }
NuclearStackExtensionImpl<InnerStackImpl>& /**
operator=( NuclearStackExtensionImpl<InnerStackImpl> &&)= default; * Function to copy particle at location i1 in stack to i2
*/
void init() { void Copy(const unsigned int i1, const unsigned int i2) {
super_type::init(); // index range check
} if (i1 >= GetSize() || i2 >= GetSize()) {
std::ostringstream err;
void dump() { err << "NuclearStackExtension: trying to access data beyond size of stack!";
super_type::dump(); throw std::runtime_error(err.str());
} }
// copy internal particle data p[i2] = p[i1]
void clear() { InnerStackImpl::Copy(i1, i2);
super_type::clear(); // check if any of p[i1] or p[i2] was a Code::Nucleus
nucleusRef_.clear(); const int ref1 = fNucleusRef[i1];
nuclearA_.clear(); const int ref2 = fNucleusRef[i2];
nuclearZ_.clear(); if (ref2 < 0) {
} if (ref1 >= 0) {
// i1 is nucleus, i2 is not
unsigned int getSize() const { fNucleusRef[i2] = GetNucleusNextRef();
return nucleusRef_.size(); fNuclearA[fNucleusRef[i2]] = fNuclearA[ref1];
} fNuclearZ[fNucleusRef[i2]] = fNuclearZ[ref1];
} else {
unsigned int getCapacity() const { // neither i1 nor i2 are nuclei
return nucleusRef_.capacity(); }
} } else {
if (ref1 >= 0) {
void setNuclearA(const unsigned int i, const unsigned short vA) { // both are nuclei, i2 is overwritten with nucleus i1
nuclearA_[getNucleusRef(i)] = vA; // fNucleusRef stays the same, but A and Z data is overwritten
} fNuclearA[ref2] = fNuclearA[ref1];
fNuclearZ[ref2] = fNuclearZ[ref1];
void setNuclearZ(const unsigned int i, const unsigned short vZ) { } else {
nuclearZ_[getNucleusRef(i)] = vZ; // i2 is overwritten with non-nucleus i1
} fNucleusRef[i2] = -1; // flag as non-nucleus
fNuclearA.erase(fNuclearA.cbegin() + ref2); // remove data for i2
void setNucleusRef(const unsigned int i, const int v) { fNuclearZ.erase(fNuclearZ.cbegin() + ref2); // remove data for i2
nucleusRef_[i] = v; const int n = fNucleusRef.size(); // update fNucleusRef: indices above ref2
} // must be decremented by 1
for (int i = 0; i < n; ++i) {
int getNuclearA(const unsigned int i) const { if (fNucleusRef[i] > ref2) { fNucleusRef[i] -= 1; }
return nuclearA_[getNucleusRef(i)]; }
} }
}
int getNuclearZ(const unsigned int i) const { }
return nuclearZ_[getNucleusRef(i)];
} /**
// this function will create new storage for Nuclear Properties, and return the * Function to copy particle at location i2 in stack to i1
// reference to it */
int getNucleusNextRef() { void Swap(const unsigned int i1, const unsigned int i2) {
nuclearA_.push_back(0); // index range check
nuclearZ_.push_back(0); if (i1 >= GetSize() || i2 >= GetSize()) {
return nuclearA_.size() - 1; std::ostringstream err;
} err << "NuclearStackExtension: trying to access data beyond size of stack!";
throw std::runtime_error(err.str());
int getNucleusRef(const unsigned int i) const { }
if (nucleusRef_[i] >= 0) return nucleusRef_[i]; // swap original particle data
std::ostringstream err; InnerStackImpl::Swap(i1, i2);
err << "NuclearStackExtension: no nucleus at ref=" << i; // swap corresponding nuclear reference data
throw std::runtime_error(err.str()); std::swap(fNucleusRef[i2], fNucleusRef[i1]);
} }
bool isNucleus(const unsigned int i) const { return nucleusRef_[i] >= 0; } void IncrementSize() {
InnerStackImpl::IncrementSize();
/** fNucleusRef.push_back(-1);
* Function to copy particle at location i1 in stack to i2 }
*/
void copy(const unsigned int i1, const unsigned int i2) { void DecrementSize() {
// index range check InnerStackImpl::DecrementSize();
if (i1 >= getSize() || i2 >= getSize()) { if (fNucleusRef.size() > 0) {
std::ostringstream err; const int ref = fNucleusRef.back();
err << "NuclearStackExtension: trying to access data beyond size of stack!"; fNucleusRef.pop_back();
throw std::runtime_error(err.str()); if (ref >= 0) {
} fNuclearA.erase(fNuclearA.begin() + ref);
// copy internal particle data p[i2] = p[i1] fNuclearZ.erase(fNuclearZ.begin() + ref);
super_type::copy(i1, i2); const int n = fNucleusRef.size();
// check if any of p[i1] or p[i2] was a Code::Nucleus for (int i = 0; i < n; ++i) {
const int ref1 = nucleusRef_[i1]; if (fNucleusRef[i] >= ref) { fNucleusRef[i] -= 1; }
const int ref2 = nucleusRef_[i2]; }
if (ref2 < 0) { }
if (ref1 >= 0) { }
// i1 is nucleus, i2 is not }
nucleusRef_[i2] = getNucleusNextRef();
nuclearA_[nucleusRef_[i2]] = nuclearA_[ref1]; private:
nuclearZ_[nucleusRef_[i2]] = nuclearZ_[ref1]; /// the actual memory to store particle data
} else {
// neither i1 nor i2 are nuclei std::vector<int> fNucleusRef;
} std::vector<unsigned short> fNuclearA;
} else { std::vector<unsigned short> fNuclearZ;
if (ref1 >= 0) {
// both are nuclei, i2 is overwritten with nucleus i1 }; // end class NuclearStackExtensionImpl
// fNucleusRef stays the same, but A and Z data is overwritten
nuclearA_[ref2] = nuclearA_[ref1]; // template<typename StackIteratorInterface>
nuclearZ_[ref2] = nuclearZ_[ref1]; // using NuclearParticleInterfaceType<StackIteratorInterface> =
} else { // NuclearParticleInterface< ,StackIteratorInterface>
// i2 is overwritten with non-nucleus i1
nucleusRef_[i2] = -1; // flag as non-nucleus // works, but requires stupd _PI class
nuclearA_.erase(nuclearA_.cbegin() + ref2); // remove data for i2 // template<typename SS> using TEST =
nuclearZ_.erase(nuclearZ_.cbegin() + ref2); // remove data for i2 // NuclearParticleInterface<corsika::super_stupid::SuperStupidStack::PIType,
const int n = nucleusRef_.size(); // update fNucleusRef: indices above ref2 // SS>;
// must be decremented by 1 template <typename InnerStack, template <typename> typename _PI>
for (int i = 0; i < n; ++i) { using NuclearStackExtension =
if (nucleusRef_[i] > ref2) { nucleusRef_[i] -= 1; } Stack<NuclearStackExtensionImpl<typename InnerStack::StackImpl>, _PI>;
}
} // ----
}
} // I'm dont't manage to do this properly.......
/*
/** template<typename TT, typename SS> using TESTi = typename
* Function to copy particle at location i2 in stack to i1 NuclearParticleInterface<TT::template PIType, SS>::ExtendedParticleInterface;
*/ template<typename TT, typename SS> using TEST1 = TESTi<TT, SS>;
void swap(const unsigned int i1, const unsigned int i2) { template<typename SS> using TEST2 = TEST1<typename
// index range check corsika::super_stupid::SuperStupidStack, SS>;
if (i1 >= getSize() || i2 >= getSize()) {
std::ostringstream err; using NuclearStackExtension = Stack<NuclearStackExtensionImpl<typename
err << "NuclearStackExtension: trying to access data beyond size of stack!"; InnerStack::StackImpl>, TEST2>;
throw std::runtime_error(err.str()); */
} /*
// swap original particle data // .... this should be fixed ....
super_type::swap(i1, i2);
// swap corresponding nuclear reference data template <typename InnerStack, typename SS=StackIteratorInterface>
std::swap(nucleusRef_[i2], nucleusRef_[i1]); //using NuclearStackExtension = Stack<NuclearStackExtensionImpl<typename
} InnerStack::StackImpl>, NuclearParticleInterface<typename InnerStack::template PIType,
StackIteratorInterface>::ExtendedParticleInterface>; using NuclearStackExtension =
void incrementSize() { Stack<NuclearStackExtensionImpl<typename InnerStack::StackImpl>, TEST1<typename
super_type::incrementSize(); corsika::super_stupid::SuperStupidStack, SS> >;
nucleusRef_.push_back(-1);
} //template <typename InnerStack>
// using NuclearStackExtension = Stack<NuclearStackExtensionImpl<typename
void decrementSize() { InnerStack::StackImpl>, TEST<typename
super_type::decrementSize(); corsika::super_stupid::SuperStupidStack::PIType>>;
if (nucleusRef_.size() > 0) { //using NuclearStackExtension = Stack<NuclearStackExtensionImpl<typename
const int ref = nucleusRef_.back(); InnerStack::StackImpl>, TEST>;
nucleusRef_.pop_back(); */
if (ref >= 0) {
nuclearA_.erase(nuclearA_.begin() + ref); } // namespace nuclear_extension
nuclearZ_.erase(nuclearZ_.begin() + ref);
const int n = nucleusRef_.size();
for (int i = 0; i < n; ++i) {
if (nucleusRef_[i] >= ref) { nucleusRef_[i] -= 1; }
}
}
}
}
private:
/// the actual memory to store particle data
nucleus_ref_type nucleusRef_;
nuclear_a_type nuclearA_;
nuclear_z_type nuclearZ_;
}; // end class NuclearStackExtensionImpl
template <typename InnerStack, template <typename> typename _PI>
using NuclearStackExtension =
Stack<NuclearStackExtensionImpl<typename InnerStack::StackImpl>, _PI> ;
//
template <typename StackIter>
using ExtendedParticleInterfaceType = NuclearParticleInterface<SuperStupidStack, StackIter>;
// the particle data stack with extra nuclear information:
using ParticleDataStack = NuclearStackExtension<SuperStupidStack, ExtendedParticleInterfaceType>;
} // namespace corsika } // namespace corsika
corsika/stack/SuperStupidStack.hpp
+ 2
− 2
View file @ afeb2e86
...@@ -108,8 +108,8 @@ namespace corsika { ...@@ -108,8 +108,8 @@ namespace corsika {
corsika::Vector<dimensionless_d> GetDirection() const { corsika::Vector<dimensionless_d> GetDirection() const {
return GetMomentum() / GetEnergy(); return GetMomentum() / GetEnergy();
} }
HEPMassType GetMass() const { return corsika::GetMass(GetPID()); } HEPMassType GetMass() const { return corsika::mass(GetPID()); }
int16_t GetChargeNumber() const { return corsika::GetChargeNumber(GetPID()); } int16_t GetChargeNumber() const { return corsika::charge_number(GetPID()); }
///@} ///@}
}; };
......
tests/modules/testSibyll.cpp
+ 11
− 11
View file @ afeb2e86
...@@ -23,26 +23,26 @@ using namespace corsika::sibyll; ...@@ -23,26 +23,26 @@ using namespace corsika::sibyll;
TEST_CASE("Sibyll", "[processes]") { TEST_CASE("Sibyll", "[processes]") {
SECTION("Sibyll -> Corsika") { SECTION("Sibyll -> Corsika") {
REQUIRE(Electron::GetCode() == REQUIRE(Code::Electron ==
corsika::sibyll::ConvertFromSibyll(corsika::sibyll::SibyllCode::Electron)); corsika::sibyll::ConvertFromSibyll(corsika::sibyll::SibyllCode::Electron));
} }
SECTION("Corsika -> Sibyll") { SECTION("Corsika -> Sibyll") {
REQUIRE(corsika::sibyll::ConvertToSibyll(Electron::GetCode()) == REQUIRE(corsika::sibyll::ConvertToSibyll(Code::Electron) ==
corsika::sibyll::SibyllCode::Electron); corsika::sibyll::SibyllCode::Electron);
REQUIRE(corsika::sibyll::ConvertToSibyllRaw(Proton::GetCode()) == 13); REQUIRE(corsika::sibyll::ConvertToSibyllRaw(Code::Proton) == 13);
} }
SECTION("canInteractInSibyll") { SECTION("canInteractInSibyll") {
REQUIRE(corsika::sibyll::CanInteract(Proton::GetCode())); REQUIRE(corsika::sibyll::CanInteract(Code::Proton));
REQUIRE(corsika::sibyll::CanInteract(Code::XiCPlus)); REQUIRE(corsika::sibyll::CanInteract(Code::XiCPlus));
REQUIRE_FALSE(corsika::sibyll::CanInteract(Electron::GetCode())); REQUIRE_FALSE(corsika::sibyll::CanInteract(Code::Electron));
REQUIRE_FALSE(corsika::sibyll::CanInteract(SigmaC0::GetCode())); REQUIRE_FALSE(corsika::sibyll::CanInteract(Code::SigmaC0));
REQUIRE_FALSE(corsika::sibyll::CanInteract(Nucleus::GetCode())); REQUIRE_FALSE(corsika::sibyll::CanInteract(Code::Nucleus));
REQUIRE_FALSE(corsika::sibyll::CanInteract(Helium::GetCode())); REQUIRE_FALSE(corsika::sibyll::CanInteract(Code::Helium));
} }
SECTION("cross-section type") { SECTION("cross-section type") {
...@@ -101,7 +101,7 @@ TEST_CASE("SibyllInterface", "[processes]") { ...@@ -101,7 +101,7 @@ TEST_CASE("SibyllInterface", "[processes]") {
corsika::setup::Stack stack; corsika::setup::Stack stack;
const HEPEnergyType E0 = 100_GeV; const HEPEnergyType E0 = 100_GeV;
HEPMomentumType P0 = sqrt(E0 * E0 - Proton::GetMass() * Proton::GetMass()); HEPMomentumType P0 = sqrt(E0 * E0 - Proton::mass() * Proton::mass());
auto plab = corsika::MomentumVector(cs, {0_GeV, 0_GeV, -P0}); auto plab = corsika::MomentumVector(cs, {0_GeV, 0_GeV, -P0});
corsika::Point pos(cs, 0_m, 0_m, 0_m); corsika::Point pos(cs, 0_m, 0_m, 0_m);
auto particle = stack.AddParticle( auto particle = stack.AddParticle(
...@@ -123,7 +123,7 @@ TEST_CASE("SibyllInterface", "[processes]") { ...@@ -123,7 +123,7 @@ TEST_CASE("SibyllInterface", "[processes]") {
setup::Stack stack; setup::Stack stack;
const HEPEnergyType E0 = 400_GeV; const HEPEnergyType E0 = 400_GeV;
HEPMomentumType P0 = sqrt(E0 * E0 - Proton::GetMass() * Proton::GetMass()); HEPMomentumType P0 = sqrt(E0 * E0 - Proton::mass() * Proton::mass());
auto plab = corsika::MomentumVector(cs, {0_GeV, 0_GeV, -P0}); auto plab = corsika::MomentumVector(cs, {0_GeV, 0_GeV, -P0});
corsika::Point pos(cs, 0_m, 0_m, 0_m); corsika::Point pos(cs, 0_m, 0_m, 0_m);
...@@ -147,7 +147,7 @@ TEST_CASE("SibyllInterface", "[processes]") { ...@@ -147,7 +147,7 @@ TEST_CASE("SibyllInterface", "[processes]") {
setup::Stack stack; setup::Stack stack;
const HEPEnergyType E0 = 10_GeV; const HEPEnergyType E0 = 10_GeV;
HEPMomentumType P0 = sqrt(E0 * E0 - Proton::GetMass() * Proton::GetMass()); HEPMomentumType P0 = sqrt(E0 * E0 - Proton::mass() * Proton::mass());
auto plab = corsika::MomentumVector(cs, {0_GeV, 0_GeV, -P0}); auto plab = corsika::MomentumVector(cs, {0_GeV, 0_GeV, -P0});
corsika::Point pos(cs, 0_m, 0_m, 0_m); corsika::Point pos(cs, 0_m, 0_m, 0_m);
auto particle = stack.AddParticle( auto particle = stack.AddParticle(
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment

data privacy - impressum