IAP GITLAB

Skip to content
Snippets Groups Projects

Compare revisions

Changes are shown as if the source revision was being merged into the target revision. Learn more about comparing revisions.

Source

Select target project
No results found

Target

Select target project
  • AirShowerPhysics/corsika
  • rulrich/corsika
  • AAAlvesJr/corsika
  • Andre/corsika
  • arrabito/corsika
  • Nikos/corsika
  • olheiser73/corsika
  • AirShowerPhysics/papers/corsika
  • pranav/corsika
9 results
Show changes
Hide whitespace changes
Inline Side-by-side
Showing
with 2087 additions and 1256 deletions
corsika/detail/modules/sibyll/Interaction.inl deleted 100644 → 0
View file @ cd6c888d
/*
* (c) Copyright 2018 CORSIKA Project, corsika-project@lists.kit.edu
*
* 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.
*/
#pragma once
#include <corsika/modules/sibyll/Interaction.hpp>
#include <corsika/media/Environment.hpp>
#include <corsika/media/NuclearComposition.hpp>
#include <corsika/framework/geometry/FourVector.hpp>
#include <corsika/modules/sibyll/ParticleConversion.hpp>
#include <corsika/modules/sibyll/SibStack.hpp>
#include <corsika/framework/utility/COMBoost.hpp>
#include <sibyll2.3d.hpp>
#include <tuple>
namespace corsika::sibyll {
inline Interaction::Interaction(const bool sibyll_printout_on)
: sibyll_listing_(sibyll_printout_on) {
// initialize Sibyll
static bool initialized = false;
if (!initialized) {
sibyll_ini_();
initialized = true;
}
}
inline Interaction::~Interaction() {
CORSIKA_LOG_DEBUG("Sibyll::Interaction n={}, Nnuc={}", count_, nucCount_);
}
inline std::tuple<corsika::CrossSectionType, corsika::CrossSectionType>
Interaction::getCrossSection(const corsika::Code BeamId, const corsika::Code TargetId,
const corsika::HEPEnergyType CoMenergy) const {
double sigProd, sigEla, dummy, dum1, dum3, dum4;
double dumdif[3];
const int iBeam = corsika::sibyll::getSibyllXSCode(
BeamId); // 0 (can not interact, 1: proton-like, 2: pion-like, 3:kaon-like)
if (!iBeam)
throw std::runtime_error(
fmt::format("Interaction of beam {} not defined in "
"Sibyll!",
BeamId));
if (!isValidCoMEnergy(CoMenergy)) {
throw std::runtime_error(
"Interaction: getCrossSection: CoM energy outside range for Sibyll!");
}
const double dEcm = CoMenergy / 1_GeV;
// single nucleon target (p,n, hydrogen) or 4<=A<=18
if (isValidTarget(TargetId)) {
// single nucleon target
if (TargetId == corsika::Code::Proton || TargetId == Code::Hydrogen ||
TargetId == Code::Neutron) {
sib_sigma_hp_(iBeam, dEcm, dum1, sigEla, sigProd, dumdif, dum3, dum4);
} else {
// nuclear target
const int iTarget = corsika::get_nucleus_A(TargetId);
sib_sigma_hnuc_(iBeam, iTarget, dEcm, sigProd, dummy, sigEla);
}
} else {
// throw std::runtime_error(
// "Sibyll nuclear target outside range. Only nuclei with 4<=A<18 are
// allowed.");
// no interaction in sibyll possible, return infinite cross section? or throw?
sigProd = std::numeric_limits<double>::infinity();
sigEla = std::numeric_limits<double>::infinity();
}
return std::make_tuple(sigProd * 1_mb, sigEla * 1_mb);
}
template <typename TParticle>
inline corsika::GrammageType Interaction::getInteractionLength(
TParticle const& projectile) const {
const corsika::Code corsikaBeamId = projectile.getPID();
// beam corsika for sibyll : 1, 2, 3 for p, pi, k
// read from cross section code table
const bool kInteraction = corsika::sibyll::canInteract(corsikaBeamId);
MomentumVector const& pLab = projectile.getMomentum();
CoordinateSystemPtr const& labCS = pLab.getCoordinateSystem();
// FOR NOW: assume target is at rest
MomentumVector pTarget(labCS, {0_GeV, 0_GeV, 0_GeV});
// total momentum and energy
HEPEnergyType Elab = projectile.getEnergy() + constants::nucleonMass;
MomentumVector pTotLab(labCS, {0_GeV, 0_GeV, 0_GeV});
pTotLab += pLab;
pTotLab += pTarget;
auto const pTotLabNorm = pTotLab.getNorm();
// calculate cm. energy
const HEPEnergyType ECoM = sqrt(
(Elab + pTotLabNorm) * (Elab - pTotLabNorm)); // binomial for numerical accuracy
CORSIKA_LOG_DEBUG(
"Interaction: LambdaInt: \n"
" input energy: {} GeV "
" beam can interact: {} "
" beam pid: {}",
projectile.getEnergy() / 1_GeV, kInteraction, projectile.getPID());
// TODO: move limits into variables
// FR: removed && Elab >= 8.5_GeV
if (kInteraction && isValidCoMEnergy(ECoM)) {
// get target from environment
/*
the target should be defined by the Environment,
ideally as full particle object so that the four momenta
and the boosts can be defined..
*/
auto const* currentNode = projectile.getNode();
const auto& mediumComposition =
currentNode->getModelProperties().getNuclearComposition();
si::CrossSectionType weightedProdCrossSection = mediumComposition.getWeightedSum(
[=](corsika::Code targetID) -> si::CrossSectionType {
return std::get<0>(this->getCrossSection(corsikaBeamId, targetID, ECoM));
});
CORSIKA_LOG_DEBUG(
"Interaction: "
"IntLength: weighted CrossSection (mb): {} ",
weightedProdCrossSection / 1_mb);
// calculate interaction length in medium
GrammageType const int_length = mediumComposition.getAverageMassNumber() *
constants::u / weightedProdCrossSection;
CORSIKA_LOG_DEBUG(
"Interaction: "
"interaction length (g/cm2): {} ",
int_length / (0.001_kg) * 1_cm * 1_cm);
return int_length;
}
return std::numeric_limits<double>::infinity() * 1_g / (1_cm * 1_cm);
}
/**
In this function SIBYLL is called to produce one event. The
event is copied (and boosted) into the shower lab frame.
*/
template <typename TSecondaryView>
inline void Interaction::doInteraction(TSecondaryView& view) {
auto const projectile = view.getProjectile();
const auto corsikaBeamId = projectile.getPID();
if (corsika::is_nucleus(corsikaBeamId)) {
// nuclei handled by different process, this should not happen
throw std::runtime_error("Nuclear projectile are not handled by SIBYLL!");
}
// position and time of interaction, not used in Sibyll
Point const pOrig = projectile.getPosition();
TimeType const tOrig = projectile.getTime();
// define projectile
HEPEnergyType const eProjectileLab = projectile.getEnergy();
auto const pProjectileLab = projectile.getMomentum();
CoordinateSystemPtr const& originalCS = pProjectileLab.getCoordinateSystem();
CORSIKA_LOG_DEBUG(
"ProcessSibyll: "
"DoInteraction: pid {} interaction ",
corsikaBeamId);
// define target
// for Sibyll is always a single nucleon
// FOR NOW: target is always at rest
const auto eTargetLab = 0_GeV + constants::nucleonMass;
const auto pTargetLab = MomentumVector(originalCS, 0_GeV, 0_GeV, 0_GeV);
const FourVector PtargLab(eTargetLab, pTargetLab);
CORSIKA_LOG_DEBUG(
"Interaction: ebeam lab: {} GeV"
"Interaction: pbeam lab: {} GeV",
eProjectileLab / 1_GeV, pProjectileLab.getComponents());
CORSIKA_LOG_DEBUG(
"Interaction: etarget lab: {} GeV "
"Interaction: ptarget lab: {} GeV",
eTargetLab / 1_GeV, pTargetLab.getComponents() / 1_GeV);
const FourVector PprojLab(eProjectileLab, pProjectileLab);
// define target kinematics in lab frame
// define boost to and from CoM frame
// CoM frame definition in Sibyll projectile: +z
COMBoost const boost(PprojLab, constants::nucleonMass);
auto const& csPrime = boost.getRotatedCS();
// just for show:
// boost projecticle
[[maybe_unused]] auto const PprojCoM = boost.toCoM(PprojLab);
// boost target
[[maybe_unused]] auto const PtargCoM = boost.toCoM(PtargLab);
CORSIKA_LOG_DEBUG(
"Interaction: ebeam CoM: {} GeV "
"Interaction: pbeam CoM: {} GeV ",
PprojCoM.getTimeLikeComponent() / 1_GeV,
PprojCoM.getSpaceLikeComponents().getComponents(csPrime) / 1_GeV);
CORSIKA_LOG_DEBUG(
"Interaction: etarget CoM: {} GeV "
"Interaction: ptarget CoM: {} GeV ",
PtargCoM.getTimeLikeComponent() / 1_GeV,
PtargCoM.getSpaceLikeComponents().getComponents(csPrime) / 1_GeV);
CORSIKA_LOG_DEBUG("Interaction: position of interaction: {} ",
pOrig.getCoordinates());
CORSIKA_LOG_DEBUG("Interaction: time: {} ", tOrig);
HEPEnergyType Etot = eProjectileLab + eTargetLab;
MomentumVector Ptot = projectile.getMomentum();
// invariant mass, i.e. cm. energy
HEPEnergyType Ecm = sqrt(Etot * Etot - Ptot.getSquaredNorm());
// sample target mass number
auto const* currentNode = projectile.getNode();
auto const& mediumComposition =
currentNode->getModelProperties().getNuclearComposition();
// get cross sections for target materials
/*
Here we read the cross section from the interaction model again,
should be passed from getInteractionLength if possible
*/
//#warning reading interaction cross section again, should not be necessary
auto const& compVec = mediumComposition.getComponents();
std::vector<CrossSectionType> cross_section_of_components(compVec.size());
for (size_t i = 0; i < compVec.size(); ++i) {
auto const targetId = compVec[i];
const auto [sigProd, sigEla] = getCrossSection(corsikaBeamId, targetId, Ecm);
[[maybe_unused]] const auto& dummy_sigEla = sigEla;
cross_section_of_components[i] = sigProd;
}
const auto targetCode =
mediumComposition.sampleTarget(cross_section_of_components, RNG_);
CORSIKA_LOG_DEBUG("Interaction: target selected: {} ", targetCode);
/*
FOR NOW: allow nuclei with A<18 or protons only.
when medium composition becomes more complex, approximations will have to be
allowed air in atmosphere also contains some Argon.
*/
int targetSibCode = -1;
if (is_nucleus(targetCode)) targetSibCode = get_nucleus_A(targetCode);
if (targetCode == Proton::code) targetSibCode = 1;
CORSIKA_LOG_DEBUG("Interaction: sibyll code: {}", targetSibCode);
if (targetSibCode > maxTargetMassNumber_ || targetSibCode < 1)
throw std::runtime_error(
"Sibyll target outside range. Only nuclei with A<18 or protons are "
"allowed.");
// beam id for sibyll
const int kBeam = corsika::sibyll::convertToSibyllRaw(corsikaBeamId);
CORSIKA_LOG_DEBUG(
"Interaction: "
" DoInteraction: E(GeV): {} "
" Ecm(GeV): {} ",
eProjectileLab / 1_GeV, Ecm / 1_GeV);
if (Ecm > getMaxEnergyCoM())
throw std::runtime_error("Interaction::DoInteraction: CoM energy too high!");
// FR: removed eProjectileLab < 8.5_GeV ||
if (Ecm < getMinEnergyCoM()) {
CORSIKA_LOG_DEBUG(
"Interaction: "
" DoInteraction: should have dropped particle.. "
"THIS IS AN ERROR");
throw std::runtime_error("energy too low for SIBYLL");
} else {
count_++;
// Sibyll does not know about units..
const double sqs = Ecm / 1_GeV;
// running sibyll, filling stack
sibyll_(kBeam, targetSibCode, sqs);
if (sibyll_listing_) {
// print final state
int print_unit = 6;
sib_list_(print_unit);
nucCount_ += get_nwounded() - 1;
}
// add particles from sibyll to stack
// link to sibyll stack
SibStack ss;
MomentumVector Plab_final(originalCS, {0.0_GeV, 0.0_GeV, 0.0_GeV});
HEPEnergyType Elab_final = 0_GeV, Ecm_final = 0_GeV;
for (auto& psib : ss) {
// abort on particles that have decayed in Sibyll. Should not happen!
if (psib.hasDecayed())
throw std::runtime_error("found particle that decayed in SIBYLL!");
// transform 4-momentum to lab. frame
// note that the momentum needs to be rotated back
auto const tmp = psib.getMomentum().getComponents();
auto const pCoM = Vector<hepmomentum_d>(csPrime, tmp);
HEPEnergyType const eCoM = psib.getEnergy();
auto const Plab = boost.fromCoM(FourVector(eCoM, pCoM));
auto const p3lab = Plab.getSpaceLikeComponents();
assert(p3lab.getCoordinateSystem() == originalCS); // just to be sure!
// add to corsika stack
auto pnew = view.addSecondary(std::make_tuple(
corsika::sibyll::convertFromSibyll(psib.getPID()), p3lab, pOrig, tOrig));
Plab_final += pnew.getMomentum();
Elab_final += pnew.getEnergy();
Ecm_final += psib.getEnergy();
}
CORSIKA_LOG_DEBUG(
"conservation (all GeV):"
"Ecm_initial(per nucleon)={}, Ecm_final(per nucleon)={}, "
"Elab_initial={}, Elab_final={}, "
"diff (%)={}, "
"E in nucleons={}, "
"Plab_initial={}, "
"Plab_final={} ",
Ecm / 1_GeV, Ecm_final * 2. / (get_nwounded() + 1) / 1_GeV, Etot / 1_GeV,
Elab_final / 1_GeV, (Elab_final / Etot / get_nwounded() - 1) * 100,
constants::nucleonMass * get_nwounded() / 1_GeV,
(pProjectileLab / 1_GeV).getComponents(), (Plab_final / 1_GeV).getComponents());
}
}
} // namespace corsika::sibyll
corsika/detail/modules/sibyll/InteractionModel.inl 0 → 100644
View file @ 00b7c3f3
/*
* (c) Copyright 2022 CORSIKA Project, corsika-project@lists.kit.edu
*
* This software is distributed under the terms of the 3-clause BSD license.
* See file LICENSE for a full version of the license.
*/
#pragma once
#include <corsika/framework/core/ParticleProperties.hpp>
#include <corsika/framework/core/PhysicalUnits.hpp>
#include <corsika/framework/geometry/FourVector.hpp>
#include <corsika/modules/sibyll/HadronInteractionModel.hpp>
#include <corsika/modules/sibyll/NuclearInteractionModel.hpp>
namespace corsika::sibyll {
inline InteractionModel::InteractionModel(std::set<Code> const& nuccomp,
std::set<Code> const& list)
: hadronSibyll_{list}
, nuclearSibyll_{hadronSibyll_, nuccomp} {}
inline HadronInteractionModel& InteractionModel::getHadronInteractionModel() {
return hadronSibyll_;
}
inline HadronInteractionModel const& InteractionModel::getHadronInteractionModel()
const {
return hadronSibyll_;
}
inline typename InteractionModel::nuclear_model_type&
InteractionModel::getNuclearInteractionModel() {
return nuclearSibyll_;
}
inline typename InteractionModel::nuclear_model_type const&
InteractionModel::getNuclearInteractionModel() const {
return nuclearSibyll_;
}
inline CrossSectionType InteractionModel::getCrossSection(
Code projCode, Code targetCode, FourMomentum const& proj4mom,
FourMomentum const& target4mom) const {
if (is_nucleus(projCode))
return getNuclearInteractionModel().getCrossSection(projCode, targetCode, proj4mom,
target4mom);
else
return getHadronInteractionModel().getCrossSection(projCode, targetCode, proj4mom,
target4mom);
}
inline std::tuple<CrossSectionType, CrossSectionType>
InteractionModel::getCrossSectionInelEla(Code projCode, Code targetCode,
FourMomentum const& proj4mom,
FourMomentum const& target4mom) const {
if (is_nucleus(projCode))
return {getNuclearInteractionModel().getCrossSection(projCode, targetCode, proj4mom,
target4mom),
CrossSectionType::zero()};
else
return getHadronInteractionModel().getCrossSectionInelEla(projCode, targetCode,
proj4mom, target4mom);
}
template <typename TSecondaries>
inline void InteractionModel::doInteraction(TSecondaries& view, Code projCode,
Code targetCode,
FourMomentum const& proj4mom,
FourMomentum const& target4mom) {
if (is_nucleus(projCode))
return getNuclearInteractionModel().doInteraction(view, projCode, targetCode,
proj4mom, target4mom);
else
return getHadronInteractionModel().doInteraction(view, projCode, targetCode,
proj4mom, target4mom);
}
} // namespace corsika::sibyll
corsika/detail/modules/sibyll/NuclearInteraction.inl deleted 100644 → 0
View file @ cd6c888d
/*
* (c) Copyright 2018 CORSIKA Project, corsika-project@lists.kit.edu
*
* 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.
*/
#pragma once
#include <corsika/modules/sibyll/Interaction.hpp>
#include <corsika/modules/sibyll/NuclearInteraction.hpp>
#include <corsika/media/Environment.hpp>
#include <corsika/media/NuclearComposition.hpp>
#include <corsika/framework/geometry/FourVector.hpp>
#include <corsika/framework/core/PhysicalUnits.hpp>
#include <corsika/framework/utility/COMBoost.hpp>
#include <corsika/framework/core/Logging.hpp>
#include <nuclib.hpp>
namespace corsika::sibyll {
template <typename TEnvironment>
inline NuclearInteraction<TEnvironment>::NuclearInteraction(sibyll::Interaction& hadint,
TEnvironment const& env)
: environment_(env)
, hadronicInteraction_(hadint) {
// initialize hadronic interaction module
// check compatibility of energy ranges, someone could try to use low-energy model..
if (!hadronicInteraction_.isValidCoMEnergy(getMinEnergyPerNucleonCoM()) ||
!hadronicInteraction_.isValidCoMEnergy(getMaxEnergyPerNucleonCoM()))
throw std::runtime_error(
"NuclearInteraction: hadronic interaction model incompatible!");
// initialize nuclib
// TODO: make sure this does not overlap with sibyll
nuc_nuc_ini_();
// initialize cross sections
initializeNuclearCrossSections();
}
template <typename TEnvironment>
inline NuclearInteraction<TEnvironment>::~NuclearInteraction() {
CORSIKA_LOG_DEBUG("Nuclib::NuclearInteraction n={} Nnuc={}", count_, nucCount_);
}
template <typename TEnvironment>
inline void NuclearInteraction<TEnvironment>::printCrossSectionTable(Code pCode) {
const int k = targetComponentsIndex_.at(pCode);
Code pNuclei[] = {Code::Helium, Code::Lithium7, Code::Oxygen,
Code::Neon, Code::Argon, Code::Iron};
std::ostringstream table;
table << "Nuclear CrossSectionTable pCode=" << pCode << " :\n en/A ";
for (auto& j : pNuclei) table << std::setw(9) << j;
table << "\n";
// loop over energy bins
for (unsigned int i = 0; i < getNEnergyBins(); ++i) {
table << " " << i << " ";
for (auto& n : pNuclei) {
auto const j = get_nucleus_A(n);
table << " " << std::setprecision(5) << std::setw(8)
<< cnucsignuc_.sigma[j - 1][k][i];
}
table << "\n";
}
CORSIKA_LOG_DEBUG(table.str());
}
template <typename TEnvironment>
inline void NuclearInteraction<TEnvironment>::initializeNuclearCrossSections() {
auto& universe = *(environment_.getUniverse());
auto const allElementsInUniverse = std::invoke([&]() {
std::set<Code> allElementsInUniverse;
auto collectElements = [&](auto& vtn) {
if (vtn.hasModelProperties()) {
auto const& comp =
vtn.getModelProperties().getNuclearComposition().getComponents();
for (auto const c : comp) allElementsInUniverse.insert(c);
}
};
universe.walk(collectElements);
return allElementsInUniverse;
});
CORSIKA_LOG_DEBUG("NuclearInteraction: initializing nuclear cross sections...");
// loop over target components, at most 4!!
int k = -1;
for (auto& ptarg : allElementsInUniverse) {
++k;
CORSIKA_LOG_DEBUG("NuclearInteraction: init target component: {}", ptarg);
const int ib = get_nucleus_A(ptarg);
if (!hadronicInteraction_.isValidTarget(ptarg)) {
CORSIKA_LOG_DEBUG(
"NuclearInteraction::InitializeNuclearCrossSections: target nucleus? id={}",
ptarg);
throw std::runtime_error(
" target can not be handled by hadronic interaction model! ");
}
targetComponentsIndex_.insert(std::pair<Code, int>(ptarg, k));
// loop over energies, fNEnBins log. energy bins
for (unsigned int i = 0; i < getNEnergyBins(); ++i) {
// hard coded energy grid, has to be aligned to definition in signuc2!!, no
// comment..
const HEPEnergyType Ecm = pow(10., 1. + 1. * i) * 1_GeV;
// get p-p cross sections
auto const protonId = Code::Proton;
auto const [siginel, sigela] =
hadronicInteraction_.getCrossSection(protonId, protonId, Ecm);
const double dsig = siginel / 1_mb;
const double dsigela = sigela / 1_mb;
// loop over projectiles, mass numbers from 2 to fMaxNucleusAProjectile
for (unsigned int j = 1; j < gMaxNucleusAProjectile_; ++j) {
const int jj = j + 1;
double sig_out, dsig_out, sigqe_out, dsigqe_out;
sigma_mc_(jj, ib, dsig, dsigela, gNSample_, sig_out, dsig_out, sigqe_out,
dsigqe_out);
// write to table
cnucsignuc_.sigma[j][k][i] = sig_out;
cnucsignuc_.sigqe[j][k][i] = sigqe_out;
}
}
}
CORSIKA_LOG_DEBUG(
"NuclearInteraction: cross sections for {} "
" components initialized!",
targetComponentsIndex_.size());
for (auto& ptarg : allElementsInUniverse) { printCrossSectionTable(ptarg); }
}
template <typename TEnvironment>
inline CrossSectionType NuclearInteraction<TEnvironment>::readCrossSectionTable(
const int ia, Code pTarget, HEPEnergyType elabnuc) {
const int ib = targetComponentsIndex_.at(pTarget) + 1; // table index in fortran
auto const ECoMNuc = sqrt(2. * constants::nucleonMass * elabnuc);
if (ECoMNuc < getMinEnergyPerNucleonCoM() || ECoMNuc > getMaxEnergyPerNucleonCoM())
throw std::runtime_error("NuclearInteraction: energy outside tabulated range!");
const double e0 = elabnuc / 1_GeV;
double sig;
CORSIKA_LOG_DEBUG("ReadCrossSectionTable: {} {} {}", ia, ib, e0);
signuc2_(ia, ib, e0, sig);
CORSIKA_LOG_DEBUG("ReadCrossSectionTable: sig={}", sig);
return sig * 1_mb;
}
// TODO: remove elastic cross section?
template <typename TEnvironment>
template <typename TParticle>
std::tuple<CrossSectionType, CrossSectionType> inline NuclearInteraction<
TEnvironment>::getCrossSection(TParticle const& projectile, Code const TargetId) {
if (projectile.getPID() != Code::Nucleus)
throw std::runtime_error(
"NuclearInteraction: getCrossSection: particle not a nucleus!");
unsigned int const iBeamA = projectile.getNuclearA();
HEPEnergyType LabEnergyPerNuc = projectile.getEnergy() / iBeamA;
CORSIKA_LOG_DEBUG(
"NuclearInteraction: getCrossSection: called with: beamNuclA={} "
" TargetId={} LabEnergyPerNuc={}GeV ",
iBeamA, TargetId, LabEnergyPerNuc / 1_GeV);
// use nuclib to calc. nuclear cross sections
// TODO: for now assumes air with hard coded composition
// extend to arbitrary mixtures, requires smarter initialization
// get nuclib projectile code: nucleon number
if (iBeamA > getMaxNucleusAProjectile() || iBeamA < 2) {
CORSIKA_LOG_DEBUG(
"NuclearInteraction: beam nucleus outside allowed range for NUCLIB!"
"A=" +
std::to_string(iBeamA));
throw std::runtime_error(
"NuclearInteraction: getCrossSection: beam nucleus outside allowed range for "
"NUCLIB!");
}
if (hadronicInteraction_.isValidTarget(TargetId)) {
auto const sigProd = readCrossSectionTable(iBeamA, TargetId, LabEnergyPerNuc);
CORSIKA_LOG_DEBUG("cross section (mb): " + std::to_string(sigProd / 1_mb));
return std::make_tuple(sigProd, 0_mb);
} else {
throw std::runtime_error("target outside range.");
}
return std::make_tuple(std::numeric_limits<double>::infinity() * 1_mb,
std::numeric_limits<double>::infinity() * 1_mb);
}
template <typename TEnvironment>
template <typename TParticle>
inline GrammageType NuclearInteraction<TEnvironment>::getInteractionLength(
TParticle const& projectile) {
// coordinate system, get global frame of reference
const Code corsikaBeamId = projectile.getPID();
if (corsikaBeamId != Code::Nucleus) {
// check if target-style nucleus (enum), these are not allowed as projectile
if (is_nucleus(corsikaBeamId)) {
throw std::runtime_error(
"NuclearInteraction: getInteractionLength: Wrong nucleus type. Nuclear "
"projectiles should use NuclearStackExtension!");
} else {
// no nuclear interaction
return std::numeric_limits<double>::infinity() * 1_g / (1_cm * 1_cm);
}
}
// read from cross section code table
MomentumVector pLab = projectile.getMomentum();
CoordinateSystemPtr const& labCS = pLab.getCoordinateSystem();
// FOR NOW: assume target is at rest
MomentumVector pTarget(labCS, {0.0_GeV, 0.0_GeV, 0.0_GeV});
// total momentum and energy
HEPEnergyType Elab = projectile.getEnergy() + constants::nucleonMass;
int const nuclA = projectile.getNuclearA();
auto const ElabNuc = projectile.getEnergy() / nuclA;
MomentumVector pTotLab(labCS, {0.0_GeV, 0.0_GeV, 0.0_GeV});
pTotLab += pLab;
pTotLab += pTarget;
auto const pTotLabNorm = pTotLab.getNorm();
// calculate cm. energy
[[maybe_unused]] HEPEnergyType const ECoM = sqrt(
(Elab + pTotLabNorm) * (Elab - pTotLabNorm)); // binomial for numerical accuracy
auto const ECoMNN = sqrt(2. * ElabNuc * constants::nucleonMass);
CORSIKA_LOG_DEBUG(
"NuclearInteraction: LambdaInt: \n"
" input energy: {}GeV\n"
" input energy CoM: {}GeV\n"
" beam pid: {}\n"
" beam A: {}\n"
" input energy per nucleon: {}GeV\n"
" input energy CoM per nucleon: {}GeV ",
Elab / 1_GeV, ECoM / 1_GeV, get_name(corsikaBeamId), nuclA, ElabNuc / 1_GeV,
ECoMNN / 1_GeV);
// energy limits
// TODO: values depend on hadronic interaction model !! this is sibyll specific
if (ElabNuc >= 8.5_GeV && ECoMNN >= gMinEnergyPerNucleonCoM_ &&
ECoMNN < gMaxEnergyPerNucleonCoM_) {
// get target from environment
/*
the target should be defined by the Environment,
ideally as full particle object so that the four momenta
and the boosts can be defined..
*/
auto const* const currentNode = projectile.getNode();
auto const& mediumComposition =
currentNode->getModelProperties().getNuclearComposition();
// determine average interaction length
// weighted sum
int i = -1;
CrossSectionType weightedProdCrossSection = 0_mb;
// get weights of components from environment/medium
const auto& w = mediumComposition.getFractions();
// loop over components in medium
for (auto const targetId : mediumComposition.getComponents()) {
i++;
CORSIKA_LOG_DEBUG("NuclearInteraction: get interaction length for target: {}",
get_name(targetId));
auto const [productionCrossSection, elaCrossSection] =
getCrossSection(projectile, targetId);
[[maybe_unused]] auto& dummy_elaCrossSection = elaCrossSection;
CORSIKA_LOG_DEBUG(
"NuclearInteraction: "
"IntLength: nuclib return (mb): " +
std::to_string(productionCrossSection / 1_mb));
weightedProdCrossSection += w[i] * productionCrossSection;
}
CORSIKA_LOG_DEBUG(
"NuclearInteraction: "
"IntLength: weighted CrossSection (mb): {} ",
weightedProdCrossSection / 1_mb);
// calculate interaction length in medium
GrammageType const int_length = mediumComposition.getAverageMassNumber() *
constants::u / weightedProdCrossSection;
CORSIKA_LOG_DEBUG(
"NuclearInteraction: "
"interaction length (g/cm2): {} ",
int_length * (1_cm * 1_cm / (0.001_kg)));
return int_length;
} else {
return std::numeric_limits<double>::infinity() * 1_g / (1_cm * 1_cm);
}
}
template <typename TEnvironment>
template <typename TSecondaryView>
inline void NuclearInteraction<TEnvironment>::doInteraction(TSecondaryView& view) {
auto projectile = view.getProjectile();
// this routine superimposes different nucleon-nucleon interactions
// in a nucleus-nucleus interaction, based the SIBYLL routine SIBNUC
const auto ProjId = projectile.getPID();
// TODO: calculate projectile mass in nuclearStackExtension
// const auto ProjMass = projectile.getMass();
CORSIKA_LOG_DEBUG("NuclearInteraction: DoInteraction: called with: {}",
get_name(ProjId));
// check if target-style nucleus (enum)
if (ProjId != Code::Nucleus)
throw std::runtime_error(
"NuclearInteraction: DoInteraction: Wrong nucleus type. Nuclear projectiles "
"should use NuclearStackExtension!");
auto const ProjMass =
projectile.getNuclearZ() * Proton::mass +
(projectile.getNuclearA() - projectile.getNuclearZ()) * Neutron::mass;
CORSIKA_LOG_DEBUG("NuclearInteraction: projectile mass: {} ", ProjMass / 1_GeV);
count_++;
// position and time of interaction, not used in NUCLIB
Point pOrig = projectile.getPosition();
TimeType tOrig = projectile.getTime();
CORSIKA_LOG_DEBUG("Interaction: position of interaction: {}", pOrig.getCoordinates());
CORSIKA_LOG_DEBUG("Interaction: time: {} ", tOrig / 1_s);
// projectile nucleon number
const unsigned int kAProj = projectile.getNuclearA();
if (kAProj > getMaxNucleusAProjectile())
throw std::runtime_error("Projectile nucleus too large for NUCLIB!");
// kinematics
// define projectile nucleus
HEPEnergyType const eProjectileLab = projectile.getEnergy();
MomentumVector const pProjectileLab = projectile.getMomentum();
FourVector const PprojLab(eProjectileLab, pProjectileLab);
CoordinateSystemPtr const& labCS = pProjectileLab.getCoordinateSystem();
CORSIKA_LOG_DEBUG(
"NuclearInteraction: eProj lab: {} "
"pProj lab: {} ",
eProjectileLab / 1_GeV, pProjectileLab.getComponents() / 1_GeV);
// define projectile nucleon
HEPEnergyType const eProjectileNucLab = eProjectileLab / kAProj;
MomentumVector const pProjectileNucLab = pProjectileLab / kAProj;
FourVector const PprojNucLab(eProjectileNucLab, pProjectileNucLab);
CORSIKA_LOG_DEBUG(
"NuclearInteraction: eProjNucleon lab (GeV): {} "
"pProjNucleon lab (GeV): {} ",
eProjectileNucLab / 1_GeV, pProjectileNucLab.getComponents() / 1_GeV);
// define target
// always a nucleon
// target is always at rest
auto const eTargetNucLab = 0_GeV + constants::nucleonMass;
auto const pTargetNucLab = MomentumVector(labCS, 0_GeV, 0_GeV, 0_GeV);
FourVector const PtargNucLab(eTargetNucLab, pTargetNucLab);
CORSIKA_LOG_DEBUG(
"NuclearInteraction: etarget lab(GeV): {} "
"NuclearInteraction: ptarget lab(GeV): {} ",
eTargetNucLab / 1_GeV, pTargetNucLab.getComponents() / 1_GeV);
// center-of-mass energy in nucleon-nucleon frame
auto const PtotNN4 = PtargNucLab + PprojNucLab;
HEPEnergyType EcmNN = PtotNN4.getNorm();
CORSIKA_LOG_DEBUG("NuclearInteraction: nuc-nuc cm energy: {}", EcmNN / 1_GeV);
if (!hadronicInteraction_.isValidCoMEnergy(EcmNN)) {
CORSIKA_LOG_DEBUG(
"NuclearInteraction: nuc-nuc. CoM energy too low for hadronic "
"interaction model!");
throw std::runtime_error("NuclearInteraction: DoInteraction: energy too low!");
}
// define boost to NUCLEON-NUCLEON frame
COMBoost const boost(PprojNucLab, constants::nucleonMass);
// boost projecticle
auto const PprojNucCoM = boost.toCoM(PprojNucLab);
// boost target
auto const PtargNucCoM = boost.toCoM(PtargNucLab);
CORSIKA_LOG_DEBUG(
"Interaction: ebeam CoM: {} "
", pbeam CoM: {} ",
PprojNucCoM.getTimeLikeComponent() / 1_GeV,
PprojNucCoM.getSpaceLikeComponents().getComponents() / 1_GeV);
CORSIKA_LOG_DEBUG(
"Interaction: etarget CoM: {}"
", ptarget CoM: {}",
PtargNucCoM.getTimeLikeComponent() / 1_GeV,
PtargNucCoM.getSpaceLikeComponents().getComponents() / 1_GeV);
// sample target nucleon number
//
// proton stand-in for nucleon
const auto beamId = Code::Proton;
auto const* const currentNode = projectile.getNode();
const auto& mediumComposition =
currentNode->getModelProperties().getNuclearComposition();
CORSIKA_LOG_DEBUG("get nucleon-nucleus cross sections for target materials..");
// 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
*/
auto const& compVec = mediumComposition.getComponents();
std::vector<CrossSectionType> cross_section_of_components(compVec.size());
for (size_t i = 0; i < compVec.size(); ++i) {
auto const targetId = compVec[i];
CORSIKA_LOG_DEBUG("target component: {}", get_name(targetId));
CORSIKA_LOG_DEBUG("beam id: {}", get_name(beamId));
const auto [sigProd, sigEla] =
hadronicInteraction_.getCrossSection(beamId, targetId, EcmNN);
cross_section_of_components[i] = sigProd;
[[maybe_unused]] auto sigElaCopy = sigEla; // ONLY TO AVOID COMPILER WARNINGS
}
const auto targetCode =
mediumComposition.sampleTarget(cross_section_of_components, RNG_);
CORSIKA_LOG_DEBUG("Interaction: target selected: {}", get_name(targetCode));
/*
FOR NOW: allow nuclei with A<18 or protons only.
when medium composition becomes more complex, approximations will have to be
allowed air in atmosphere also contains some Argon.
*/
int kATarget = -1;
if (is_nucleus(targetCode))
kATarget = get_nucleus_A(targetCode);
else if (targetCode == Code::Proton)
kATarget = 1;
CORSIKA_LOG_DEBUG("NuclearInteraction: nuclib target code: " +
std::to_string(kATarget));
if (!hadronicInteraction_.isValidTarget(targetCode))
throw std::runtime_error("target outside range. ");
// end of target sampling
// superposition
CORSIKA_LOG_DEBUG(
"NuclearInteraction: sampling nuc. multiple interaction structure.. ");
// get nucleon-nucleon cross section
// (needed to determine number of nucleon-nucleon scatterings)
const auto protonId = Code::Proton;
const auto [prodCrossSection, elaCrossSection] =
hadronicInteraction_.getCrossSection(protonId, protonId, EcmNN);
const double sigProd = prodCrossSection / 1_mb;
const double sigEla = elaCrossSection / 1_mb;
// sample number of interactions (only input variables, output in common cnucms)
// nuclear multiple scattering according to glauber (r.i.p.)
int_nuc_(kATarget, kAProj, sigProd, sigEla);
CORSIKA_LOG_DEBUG(
"number of nucleons in target : {}\n"
"number of wounded nucleons in target : {}\n"
"number of nucleons in projectile : {}\n"
"number of wounded nucleons in project. : {}\n"
"number of inel. nuc.-nuc. interactions : {}\n"
"number of elastic nucleons in target : {}\n"
"number of elastic nucleons in project. : {}\n"
"impact parameter: {}",
kATarget, cnucms_.na, kAProj, cnucms_.nb, cnucms_.ni, cnucms_.nael, cnucms_.nbel,
cnucms_.b);
// calculate fragmentation
CORSIKA_LOG_DEBUG("calculating nuclear fragments..");
// number of interactions
// include elastic
const int nElasticNucleons = cnucms_.nbel;
const int nInelNucleons = cnucms_.nb;
const int nIntProj = nInelNucleons + nElasticNucleons;
const double impactPar = cnucms_.b; // only needed to avoid passing common var.
int nFragments = 0;
// number of fragments is limited to 60
int AFragments[60];
// call fragmentation routine
// input: target A, projectile A, number of int. nucleons in projectile, impact
// parameter (fm) output: nFragments, AFragments in addition the momenta ar stored
// in pf in common fragments, neglected
fragm_(kATarget, kAProj, nIntProj, impactPar, nFragments, AFragments);
// this should not occur but well :) (LCOV_EXCL_START)
if (nFragments > (int)getMaxNFragments())
throw std::runtime_error("Number of nuclear fragments in NUCLIB exceeded!");
// (LCOV_EXCL_STOP)
CORSIKA_LOG_DEBUG("number of fragments: " + std::to_string(nFragments));
CORSIKA_LOG_DEBUG("adding nuclear fragments to particle stack..");
// put nuclear fragments on corsika stack
for (int j = 0; j < nFragments; ++j) {
CORSIKA_LOG_DEBUG("fragment {}: A={} px={} py={} pz={}", j, AFragments[j],
fragments_.ppp[j][0], fragments_.ppp[j][1], fragments_.ppp[j][2]);
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 = Code::Proton;
else
specCode = Code::Nucleus;
const HEPMassType mass = get_nucleus_mass(nuclA, nuclZ);
CORSIKA_LOG_DEBUG("NuclearInteraction: adding fragment: {}", get_name(specCode));
CORSIKA_LOG_DEBUG("NuclearInteraction: A,Z: {}, {}", nuclA, nuclZ);
CORSIKA_LOG_DEBUG("NuclearInteraction: mass: {} GeV", std::to_string(mass / 1_GeV));
// CORSIKA 7 way
// spectators inherit momentum from original projectile
const double mass_ratio = mass / ProjMass;
CORSIKA_LOG_DEBUG("NuclearInteraction: mass ratio " + std::to_string(mass_ratio));
auto const Plab = PprojLab * mass_ratio;
CORSIKA_LOG_DEBUG("NuclearInteraction: fragment momentum: {}",
Plab.getSpaceLikeComponents().getComponents() / 1_GeV);
if (nuclA == 1)
// add nucleon
projectile.addSecondary(
std::make_tuple(specCode, Plab.getSpaceLikeComponents(), pOrig, tOrig));
else
// add nucleus
projectile.addSecondary(std::make_tuple(specCode, 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
CORSIKA_LOG_DEBUG("adding elastically scattered nucleons to particle stack..");
for (int j = 0; j < nElasticNucleons; ++j) {
// TODO: sample proton or neutron
auto const elaNucCode = Code::Proton;
// CORSIKA 7 way
// elastic nucleons inherit momentum from original projectile
// neglecting momentum transfer in interaction
const double mass_ratio = get_mass(elaNucCode) / ProjMass;
auto const Plab = PprojLab * mass_ratio;
projectile.addSecondary(
std::make_tuple(elaNucCode, Plab.getSpaceLikeComponents(), pOrig, tOrig));
}
// add inelastic interactions
CORSIKA_LOG_DEBUG("calculate inelastic nucleon-nucleon interactions..");
for (int j = 0; j < nInelNucleons; ++j) {
// TODO: sample neutron or proton
auto pCode = Code::Proton;
// temporarily add to stack, will be removed after interaction in DoInteraction
CORSIKA_LOG_DEBUG("inelastic interaction no. {}", j);
typename TSecondaryView::inner_stack_value_type nucleonStack;
auto inelasticNucleon = nucleonStack.addParticle(
std::make_tuple(pCode, PprojNucLab.getSpaceLikeComponents(), pOrig, tOrig));
inelasticNucleon.setNode(projectile.getNode());
// create inelastic interaction for each nucleon
CORSIKA_LOG_TRACE("calling HadronicInteraction...");
// create new StackView for each of the nucleons
TSecondaryView nucleon_secondaries(inelasticNucleon);
// all inner hadronic event generator
hadronicInteraction_.doInteraction(nucleon_secondaries);
for (const auto& pSec : nucleon_secondaries) {
projectile.addSecondary(std::make_tuple(pSec.getPID(), pSec.getMomentum(),
pSec.getPosition(), pSec.getTime()));
}
}
CORSIKA_LOG_DEBUG("NuclearInteraction: DoInteraction: done");
}
} // namespace corsika::sibyll
corsika/detail/modules/sibyll/NuclearInteractionModel.inl 0 → 100644
View file @ 00b7c3f3
/*
* (c) Copyright 2018 CORSIKA Project, corsika-project@lists.kit.edu
*
* This software is distributed under the terms of the 3-clause BSD license.
* See file LICENSE for a full version of the license.
*/
#pragma once
#include <corsika/media/Environment.hpp>
#include <corsika/media/NuclearComposition.hpp>
#include <corsika/modules/Random.hpp>
#include <corsika/framework/core/PhysicalUnits.hpp>
#include <corsika/framework/core/EnergyMomentumOperations.hpp>
#include <corsika/framework/utility/COMBoost.hpp>
#include <corsika/framework/core/Logging.hpp>
#include <nuclib.hpp>
namespace corsika::sibyll {
template <typename TNucleonModel>
inline NuclearInteractionModel<TNucleonModel>::NuclearInteractionModel(
TNucleonModel& hadint, std::set<Code> const& nuccomp)
: hadronicInteraction_(hadint) {
// initialize nuclib
// TODO: make sure this does not overlap with sibyll
corsika::connect_random_stream("sibyll", ::sibyll::set_rng_function);
nuc_nuc_ini_();
// initialize cross sections
initializeNuclearCrossSections(nuccomp);
}
template <typename TNucleonModel>
inline NuclearInteractionModel<TNucleonModel>::~NuclearInteractionModel() {
CORSIKA_LOGGER_DEBUG(logger_, "nuclear interactions handled by Sibyll n={}", count_);
}
template <typename TNucleonModel>
inline bool constexpr NuclearInteractionModel<TNucleonModel>::isValid(
Code const projectileId, Code const targetId, HEPEnergyType const sqrtSnn) const {
// also depends on underlying model, for Proton/Neutron projectile
if (!hadronicInteraction_.isValid(Code::Proton, targetId, sqrtSnn)) { return false; }
// projectile limits:
if (!is_nucleus(projectileId)) { return false; }
unsigned int projectileA = get_nucleus_A(projectileId);
if (projectileA > getMaxNucleusAProjectile() || projectileA < 2) { return false; }
return true;
} // namespace corsika::sibyll
template <typename TNucleonModel>
inline void NuclearInteractionModel<TNucleonModel>::printCrossSectionTable(
Code const pCode) const {
if (!hadronicInteraction_.isValid(Code::Proton, pCode, 100_GeV)) { // LCOV_EXCL_START
CORSIKA_LOGGER_WARN(logger_, "Invalid target type {} for hadron interaction model.",
pCode);
return;
} // LCOV_EXCL_STOP
int const k = targetComponentsIndex_.at(pCode);
Code const pNuclei[] = {Code::Helium, Code::Lithium7, Code::Oxygen,
Code::Neon, Code::Argon, Code::Iron};
std::ostringstream table;
table << "Nuclear CrossSectionTable pCode=" << pCode << " :\n en/A ";
for (auto& j : pNuclei) table << std::setw(9) << j;
table << "\n";
// loop over energy bins
for (unsigned int i = 0; i < getNEnergyBins(); ++i) {
table << " " << i << " ";
for (auto& n : pNuclei) {
auto const j = get_nucleus_A(n);
table << " " << std::setprecision(5) << std::setw(8)
<< cnucsignuc_.sigma[j - 1][k][i];
}
table << "\n";
}
CORSIKA_LOGGER_DEBUG(logger_, table.str());
}
template <typename TNucleonModel>
inline void NuclearInteractionModel<TNucleonModel>::initializeNuclearCrossSections(
std::set<Code> const& allElementsInUniverse) {
CORSIKA_LOGGER_DEBUG(logger_, "initializing nuclear cross sections...");
// loop over target components, at most 4!!
int k = -1;
for (Code const ptarg : allElementsInUniverse) {
++k;
CORSIKA_LOGGER_DEBUG(logger_, "init target component: {} A={}", ptarg,
get_nucleus_A(ptarg));
int const ib = get_nucleus_A(ptarg);
if (!hadronicInteraction_.isValid(Code::Proton, ptarg, 100_GeV)) {
CORSIKA_LOGGER_WARN(
logger_, "Invalid target type {} for hadron interaction model.", ptarg);
continue;
}
targetComponentsIndex_.insert(std::pair<Code, int>(ptarg, k));
// loop over energies, fNEnBins log. energy bins
for (size_t i = 0; i < getNEnergyBins(); ++i) {
// hard coded energy grid, has to be aligned to definition in signuc2!!, no
// comment..
HEPEnergyType const Ecm = pow(10., 1. + 1. * i) * 1_GeV;
// head-on pp collision:
HEPEnergyType const EcmHalve = Ecm / 2;
HEPMomentumType const pcm =
sqrt(EcmHalve * EcmHalve - Proton::mass * Proton::mass);
CoordinateSystemPtr cs = get_root_CoordinateSystem();
FourMomentum projectileP4(EcmHalve, {cs, pcm, 0_eV, 0_eV});
FourMomentum targetP4(EcmHalve, {cs, -pcm, 0_eV, 0_eV});
// get p-p cross sections
if (!hadronicInteraction_.isValid(Code::Proton, Code::Proton, Ecm)) {
throw std::runtime_error("invalid (projectile,target,ecm) combination");
}
auto const [siginel, sigela] = hadronicInteraction_.getCrossSectionInelEla(
Code::Proton, Code::Proton, projectileP4, targetP4);
double const dsig = siginel / 1_mb;
double const dsigela = sigela / 1_mb;
// loop over projectiles, mass numbers from 2 to fMaxNucleusAProjectile
CORSIKA_LOGGER_TRACE(logger_, "Ecm={} siginel={} sigela={}", Ecm / 1_GeV, dsig,
dsigela);
for (size_t j = 1; j < gMaxNucleusAProjectile_; ++j) {
const int jj = j + 1;
double sig_out, dsig_out, sigqe_out, dsigqe_out;
sigma_mc_(jj, ib, dsig, dsigela, gNSample_, sig_out, dsig_out, sigqe_out,
dsigqe_out);
// write to table
cnucsignuc_.sigma[j][k][i] = sig_out;
cnucsignuc_.sigqe[j][k][i] = sigqe_out;
CORSIKA_LOGGER_TRACE(logger_, "nuc A={} sig={} qe={}", j, sig_out, sigqe_out);
}
}
}
CORSIKA_LOGGER_DEBUG(logger_, "cross sections for {} components initialized!",
targetComponentsIndex_.size());
for (auto& ptarg : allElementsInUniverse) { printCrossSectionTable(ptarg); }
}
template <typename TNucleonModel>
inline CrossSectionType NuclearInteractionModel<TNucleonModel>::readCrossSectionTable(
int const ia, Code const pTarget, HEPEnergyType const elabnuc) const {
CORSIKA_LOGGER_DEBUG(logger_, "ia={}, target={}, ElabNuc={} GeV", ia, pTarget,
elabnuc / 1_GeV);
int const ib = targetComponentsIndex_.at(pTarget) + 1; // table index in fortran
auto const ECoMNuc = sqrt(2. * constants::nucleonMass * elabnuc);
CORSIKA_LOGGER_DEBUG(logger_, "sqrtSnn= {} GeV", ECoMNuc / 1_GeV);
if (ECoMNuc < getMinEnergyPerNucleonCoM() || ECoMNuc > getMaxEnergyPerNucleonCoM()) {
CORSIKA_LOGGER_WARN(
logger_,
"nucleon-nucleon energy outside range! sqrtSnn={}GeV (limits: {} .. {} GeV)",
ECoMNuc / 1_GeV, getMinEnergyPerNucleonCoM() / 1_GeV,
getMaxEnergyPerNucleonCoM() / 1_GeV);
// throw std::runtime_error("energy outside tabulated range!");
}
double const e0 = elabnuc / 1_GeV;
double sig;
CORSIKA_LOGGER_DEBUG(logger_, "ReadCrossSectionTable: {} {} {}", ia, ib, e0);
signuc2_(ia, ib, e0, sig);
CORSIKA_LOGGER_DEBUG(logger_, "ReadCrossSectionTable: sig={}", sig);
return sig * 1_mb;
}
template <typename TNucleonModel>
CrossSectionType inline NuclearInteractionModel<TNucleonModel>::getCrossSection(
Code const projectileId, Code const targetId, FourMomentum const& projectileP4,
FourMomentum const& targetP4) const {
CORSIKA_LOGGER_DEBUG(logger_, "projectile: E={}, p3={} \n target: E={}, p3={}",
projectileP4.getTimeLikeComponent() / 1_GeV,
projectileP4.getSpaceLikeComponents() / 1_GeV,
targetP4.getTimeLikeComponent() / 1_GeV,
targetP4.getSpaceLikeComponents() / 1_GeV);
// check if projectile and target are nuclei!
if (!is_nucleus(projectileId) || !is_nucleus(targetId)) {
return CrossSectionType::zero();
}
// calculate sqrt(Snn) (only works if projectile and target are nuclei)
HEPEnergyType const sqrtSnn =
(projectileP4 / get_nucleus_A(projectileId) + targetP4 / get_nucleus_A(targetId))
.getNorm();
CORSIKA_LOG_DEBUG("proj={}, targ={}, sqrtSNN={}GeV", projectileId, targetId,
sqrtSnn / 1_GeV);
if (!isValid(projectileId, targetId, sqrtSnn)) { return CrossSectionType::zero(); }
// lab-frame energy per projectile nucleon as required by signuc2()
HEPEnergyType const LabEnergyPerNuc = calculate_lab_energy(
static_pow<2>(sqrtSnn), get_mass(projectileId) / get_nucleus_A(projectileId),
get_mass(targetId) / get_nucleus_A(targetId));
auto const sigProd =
readCrossSectionTable(get_nucleus_A(projectileId), targetId, LabEnergyPerNuc);
CORSIKA_LOGGER_DEBUG(logger_, "cross section (mb): sqrtSnn={} sig={}",
sqrtSnn / 1_GeV, sigProd / 1_mb);
return sigProd;
}
template <typename TNucleonModel>
template <typename TSecondaryView>
inline void NuclearInteractionModel<TNucleonModel>::doInteraction(
TSecondaryView& view, Code const projectileId, Code const targetId,
FourMomentum const& projectileP4, FourMomentum const& targetP4) {
// model is only designed for projectile nuclei. Collisions are broken down into
// "nucleon-target" collisions.
if (!is_nucleus(projectileId)) {
throw std::runtime_error("Can only handle nuclear projectiles.");
}
size_t const projectileA = get_nucleus_A(projectileId);
// this is center-of-mass for projectile_nucleon - target
FourMomentum const nucleonP4 = projectileP4 / projectileA;
HEPEnergyType const sqrtSnucleon = (nucleonP4 + targetP4).getNorm();
if (!isValid(projectileId, targetId, sqrtSnucleon)) {
throw std::runtime_error("Invalid projectile/target/energy combination.");
}
// projectile is always nucleus!
// Elab corresponding to sqrtSnucleon -> fixed target projectile
COMBoost const boost(nucleonP4, targetP4);
CORSIKA_LOGGER_DEBUG(logger_, "pId={} tId={} sqrtSnucleon={}GeV Aproj={}",
projectileId, targetId, sqrtSnucleon / 1_GeV, projectileA);
count_++;
// lab. momentum per projectile nucleon
HEPMomentumType const pNucleonLab = nucleonP4.getSpaceLikeComponents().getNorm();
// nucleon momentum in direction of CM motion (lab system)
MomentumVector const p3NucleonLab(boost.getRotatedCS(), {0_GeV, 0_GeV, pNucleonLab});
/*
FOR NOW: allow nuclei with A<18 or protons/nucleon only.
when medium composition becomes more complex, approximations will have to be
allowed air in atmosphere also contains some Argon.
*/
int kATarget = -1;
size_t targetA = 1;
if (is_nucleus(targetId)) {
kATarget = get_nucleus_A(targetId);
targetA = kATarget;
} else if (targetId == Code::Proton || targetId == Code::Neutron ||
targetId == Code::Hydrogen) {
kATarget = 1;
}
CORSIKA_LOGGER_DEBUG(logger_, "nuclib target code: {}", kATarget);
// end of target sampling
// superposition
CORSIKA_LOGGER_DEBUG(logger_, "sampling nuc. multiple interaction structure.. ");
// get nucleon-nucleon cross section
// (needed to determine number of nucleon-nucleon scatterings)
auto const protonId = Code::Proton;
auto const [prodCrossSection, elaCrossSection] =
hadronicInteraction_.getCrossSectionInelEla(
protonId, protonId, nucleonP4,
targetP4 / targetA); // todo check, wrong RU
double const sigProd = prodCrossSection / 1_mb;
double const sigEla = elaCrossSection / 1_mb;
// sample number of interactions (only input variables, output in common cnucms)
// nuclear multiple scattering according to glauber (r.i.p.)
int_nuc_(kATarget, projectileA, sigProd, sigEla);
CORSIKA_LOGGER_DEBUG(logger_,
"number of nucleons in target : {}\n"
"number of wounded nucleons in target : {}\n"
"number of nucleons in projectile : {}\n"
"number of wounded nucleons in project. : {}\n"
"number of inel. nuc.-nuc. interactions : {}\n"
"number of elastic nucleons in target : {}\n"
"number of elastic nucleons in project. : {}\n"
"impact parameter: {}",
kATarget, cnucms_.na, projectileA, cnucms_.nb, cnucms_.ni,
cnucms_.nael, cnucms_.nbel, cnucms_.b);
// calculate fragmentation
CORSIKA_LOGGER_DEBUG(logger_, "calculating nuclear fragments..");
// number of interactions
// include elastic
int const nElasticNucleons = cnucms_.nbel;
int const nInelNucleons = cnucms_.nb;
int const nIntProj = nInelNucleons + nElasticNucleons;
double const impactPar = cnucms_.b; // only needed to avoid passing common var.
int nFragments = 0;
// number of fragments is limited to 60
int AFragments[60];
// call fragmentation routine
// input: target A, projectile A, number of int. nucleons in projectile, impact
// parameter (fm) output: nFragments, AFragments in addition the momenta ar stored
// in pf in common fragments, neglected
fragm_(kATarget, projectileA, nIntProj, impactPar, nFragments, AFragments);
// this should not occur but well :) (LCOV_EXCL_START)
if (nFragments > (int)getMaxNFragments()) {
throw std::runtime_error("Number of nuclear fragments in NUCLIB exceeded!");
}
// (LCOV_EXCL_STOP)
CORSIKA_LOGGER_DEBUG(logger_, "number of fragments: {}", nFragments);
CORSIKA_LOGGER_DEBUG(logger_, "adding nuclear fragments to particle stack..");
// put nuclear fragments on corsika stack
for (int j = 0; j < nFragments; ++j) {
CORSIKA_LOGGER_DEBUG(logger_, "fragment {}: A={} px={} py={} pz={}", j,
AFragments[j], fragments_.ppp[j][0], fragments_.ppp[j][1],
fragments_.ppp[j][2]);
auto const nuclA = AFragments[j];
// get Z from stability line
auto const nuclZ = int(nuclA / 2.15 + 0.7);
// TODO: do we need to catch single nucleons??
Code const specCode = (nuclA == 1 ?
// TODO: sample neutron or proton
Code::Proton
: get_nucleus_code(nuclA, nuclZ));
HEPMassType const mass = get_mass(specCode);
CORSIKA_LOGGER_DEBUG(logger_, "adding fragment: {}", get_name(specCode));
CORSIKA_LOGGER_DEBUG(logger_, "A,Z: {}, {}", nuclA, nuclZ);
CORSIKA_LOGGER_DEBUG(logger_, "mass: {} GeV", mass / 1_GeV);
// CORSIKA 7 way
// spectators inherit momentum from original projectile
auto const p3lab = p3NucleonLab * nuclA;
HEPEnergyType const Ekin = sqrt(p3lab.getSquaredNorm() + mass * mass) - mass;
CORSIKA_LOGGER_DEBUG(logger_, "fragment momentum {}",
p3lab.getComponents() / 1_GeV);
view.addSecondary(std::make_tuple(specCode, Ekin, p3lab.normalized()));
}
// add elastic nucleons to corsika stack
// TODO: the elastic interaction could be external like the inelastic interaction,
// e.g. use existing ElasticModel
CORSIKA_LOGGER_DEBUG(logger_,
"adding elastically scattered nucleons to particle stack..");
for (int j = 0; j < nElasticNucleons; ++j) {
// TODO: sample proton or neutron
Code const elaNucCode = Code::Proton;
// CORSIKA 7 way
// elastic nucleons inherit momentum from original projectile
// neglecting momentum transfer in interaction
auto const p3lab = p3NucleonLab;
HEPEnergyType const mass = get_mass(elaNucCode);
HEPEnergyType const Ekin = sqrt(p3lab.getSquaredNorm() + mass * mass) - mass;
view.addSecondary(std::make_tuple(elaNucCode, Ekin, p3lab.normalized()));
}
// add inelastic interactions
CORSIKA_LOGGER_DEBUG(logger_, "calculate inelastic nucleon-nucleon interactions..");
for (int j = 0; j < nInelNucleons; ++j) {
// TODO: sample neutron or proton
auto const pCode = Code::Proton;
HEPEnergyType const mass = get_mass(pCode);
HEPEnergyType const Ekin = sqrt(p3NucleonLab.getSquaredNorm() + mass * mass) - mass;
// temporarily add to stack, will be removed after interaction in DoInteraction
CORSIKA_LOGGER_DEBUG(logger_, "inelastic interaction no. {}", j);
typename TSecondaryView::inner_stack_value_type nucleonStack;
Point const pDummy(boost.getOriginalCS(), {0_m, 0_m, 0_m});
TimeType const tDummy = 0_ns;
auto inelasticNucleon = nucleonStack.addParticle(
std::make_tuple(pCode, Ekin, p3NucleonLab.normalized(), pDummy, tDummy));
inelasticNucleon.setNode(view.getProjectile().getNode());
// create inelastic interaction for each nucleon
CORSIKA_LOGGER_TRACE(logger_, "calling HadronicInteraction...");
// create new StackView for each of the nucleons
TSecondaryView nucleon_secondaries(inelasticNucleon);
// all inner hadronic event generator
hadronicInteraction_.doInteraction(nucleon_secondaries, pCode, targetId, nucleonP4,
targetP4);
for (const auto& pSec : nucleon_secondaries) {
auto const p3lab = pSec.getMomentum();
Code const pid = pSec.getPID();
HEPEnergyType const mass = get_mass(pid);
HEPEnergyType const Ekin = sqrt(p3lab.getSquaredNorm() + mass * mass) - mass;
view.addSecondary(std::make_tuple(pid, Ekin, p3lab.normalized()));
}
}
}
} // namespace corsika::sibyll
corsika/detail/modules/sibyll/ParticleConversion.inl
View file @ 00b7c3f3
/*
* (c) Copyright 2018 CORSIKA Project, corsika-project@lists.kit.edu
*
* 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.
* This software is distributed under the terms of the 3-clause BSD license.
* See file LICENSE for a full version of the license.
*/
#pragma once
......@@ -15,8 +14,7 @@
namespace corsika::sibyll {
inline HEPMassType getSibyllMass(Code const pCode) {
if (pCode == corsika::Code::Nucleus)
throw std::runtime_error("Cannot getMass() of particle::Nucleus -> unspecified");
if (is_nucleus(pCode)) throw std::runtime_error("Not defined for nuclei.");
auto sCode = convertToSibyllRaw(pCode);
if (sCode == 0)
throw std::runtime_error("getSibyllMass: unknown particle!");
......
corsika/detail/modules/sophia/InteractionModel.inl 0 → 100644
View file @ 00b7c3f3
/*
* (c) Copyright 2022 CORSIKA Project, corsika-project@lists.kit.edu
*
* This software is distributed under the terms of the 3-clause BSD license.
* See file LICENSE for a full version of the license.
*/
#pragma once
#include <corsika/framework/geometry/Point.hpp>
#include <corsika/modules/sophia/ParticleConversion.hpp>
#include <corsika/framework/utility/COMBoost.hpp>
#include <corsika/modules/Random.hpp>
#include <corsika/modules/sophia/SophiaStack.hpp>
#include <corsika/framework/core/EnergyMomentumOperations.hpp>
#include <sophia.hpp>
namespace corsika::sophia {
inline void InteractionModel::setVerbose(bool const flag) { sophia_listing_ = flag; }
inline InteractionModel::InteractionModel()
: sophia_listing_(false) {
corsika::connect_random_stream(RNG_, ::sophia::set_rng_function);
// set all particles stable in SOPHIA
for (int i = 0; i < 49; ++i) so_csydec_.idb[i] = -abs(so_csydec_.idb[i]);
}
inline InteractionModel::~InteractionModel() {
CORSIKA_LOG_DEBUG("Sophia::Model n={}", count_);
}
inline bool constexpr InteractionModel::isValid(Code const projectileId,
Code const targetId,
HEPEnergyType const sqrtSnn) const {
if ((minEnergyCoM_ > sqrtSnn) || (sqrtSnn > maxEnergyCoM_)) { return false; }
if (!(targetId == Code::Proton || targetId == Code::Neutron ||
targetId == Code::Hydrogen))
return false;
if (projectileId != Code::Photon) return false;
return true;
}
template <typename TSecondaryView>
inline void InteractionModel::doInteraction(TSecondaryView& secondaries,
Code const projectileId,
Code const targetId,
FourMomentum const& projectileP4,
FourMomentum const& targetP4) {
CORSIKA_LOGGER_DEBUG(logger_, "projectile: Id={}, E={} GeV, p3={} GeV", projectileId,
projectileP4.getTimeLikeComponent() / 1_GeV,
projectileP4.getSpaceLikeComponents().getComponents() / 1_GeV);
CORSIKA_LOGGER_DEBUG(logger_, "target: Id={}, E={} GeV, p3={} GeV", targetId,
targetP4.getTimeLikeComponent() / 1_GeV,
targetP4.getSpaceLikeComponents().getComponents() / 1_GeV);
// sqrtS per target nucleon
HEPEnergyType const sqrtS = (projectileP4 + targetP4).getNorm();
CORSIKA_LOGGER_DEBUG(logger_, "sqrtS={}GeV", sqrtS / 1_GeV);
// accepts only photon-nucleon interactions
if (!isValid(projectileId, targetId, sqrtS)) {
CORSIKA_LOGGER_ERROR(logger_,
"Invalid target/projectile/energy combination: {},{},{} GeV",
projectileId, targetId, sqrtS / 1_GeV);
throw std::runtime_error("SOPHIA: Invalid target/projectile/energy combination");
}
COMBoost const boost(projectileP4, targetP4);
int nucleonSophiaCode = convertToSophiaRaw(targetId); // either proton or neutron
// initialize resonance spectrum
initial_(nucleonSophiaCode);
// Sophia does sqrt(1 - mass_sophia / mass_c8), so we need to make sure that E0 >=
// m_sophia
double Enucleon = std::max(corsika::sophia::getSophiaMass(targetId) / 1_GeV,
targetP4.getTimeLikeComponent() / 1_GeV);
double Ephoton = projectileP4.getTimeLikeComponent() / 1_GeV;
double theta = 0.0; // set nucleon at rest in collision
int Imode = -1; // overwritten inside SOPHIA
CORSIKA_LOGGER_DEBUG(logger_,
"calling SOPHIA eventgen with L0={}, E0={}, eps={},theta={}",
nucleonSophiaCode, Enucleon, Ephoton, theta);
count_++;
// call sophia
eventgen_(nucleonSophiaCode, Enucleon, Ephoton, theta, Imode);
if (sophia_listing_) {
int arg = 3;
print_event_(arg);
}
auto const& originalCS = boost.getOriginalCS();
// SOPHIA has photon along -z and nucleon along +z (GZK calc..)
COMBoost const boostInternal(targetP4, projectileP4);
auto const& csPrime = boost.getRotatedCS();
CoordinateSystemPtr csPrimePrime =
make_rotation(csPrime, QuantityVector<length_d>{1_m, 0_m, 0_m}, M_PI);
SophiaStack ss;
MomentumVector P_final(originalCS, {0.0_GeV, 0.0_GeV, 0.0_GeV});
HEPEnergyType E_final = 0_GeV;
for (auto& psop : ss) {
// abort on particles that have decayed in SOPHIA. Should not happen!
if (psop.hasDecayed()) { // LCOV_EXCL_START
throw std::runtime_error("found particle that decayed in SOPHIA!");
} // LCOV_EXCL_STOP
auto momentumSophia = psop.getMomentum(csPrimePrime);
momentumSophia.rebase(csPrime);
auto const energySophia = psop.getEnergy();
auto const P4com = boostInternal.toCoM(FourVector{energySophia, momentumSophia});
auto const P4lab = boost.fromCoM(P4com);
SophiaCode const pidSophia = psop.getPID();
Code const pid = convertFromSophia(pidSophia);
auto momentum = P4lab.getSpaceLikeComponents();
momentum.rebase(originalCS);
HEPEnergyType const Ekin =
calculate_kinetic_energy(momentum.getNorm(), get_mass(pid));
CORSIKA_LOGGER_TRACE(logger_, "SOPHIA: pid={}, p={} GeV", pidSophia,
momentumSophia.getComponents() / 1_GeV);
CORSIKA_LOGGER_TRACE(logger_, "CORSIKA: pid={}, p={} GeV", pid,
momentum.getComponents() / 1_GeV);
auto pnew =
secondaries.addSecondary(std::make_tuple(pid, Ekin, momentum.normalized()));
P_final += pnew.getMomentum();
E_final += pnew.getEnergy();
}
CORSIKA_LOGGER_TRACE(logger_, "Efinal={} GeV,Pfinal={} GeV", E_final / 1_GeV,
P_final.getComponents() / 1_GeV);
}
} // namespace corsika::sophia
corsika/detail/modules/sophia/ParticleConversion.inl 0 → 100644
View file @ 00b7c3f3
/*
* (c) Copyright 2022 CORSIKA Project, corsika-project@lists.kit.edu
*
* This software is distributed under the terms of the 3-clause BSD license.
* See file LICENSE for a full version of the license.
*/
#pragma once
#include <corsika/framework/core/ParticleProperties.hpp>
#include <sophia.hpp>
namespace corsika::sophia {
inline HEPMassType getSophiaMass(Code const pCode) {
if (is_nucleus(pCode)) throw std::runtime_error("Not defined for nuclei.");
auto sCode = convertToSophiaRaw(pCode);
if (sCode == 0)
throw std::runtime_error("getSophiaMass: unknown particle!");
else
return sqrt(get_sophia_mass2(sCode)) * 1_GeV;
}
} // namespace corsika::sophia
\ No newline at end of file
corsika/detail/modules/tauola/Decay.inl 0 → 100644
View file @ 00b7c3f3
/*
* (c) Copyright 2023 CORSIKA Project, corsika-project@lists.kit.edu
*
* This software is distributed under the terms of the 3-clause BSD license.
* See file LICENSE for a full version of the license.
*/
#include <Tauola/Tauola.h>
#include <corsika/framework/utility/COMBoost.hpp>
namespace Tpp = Tauolapp;
namespace corsika::tauola {
// default constructor
Decay::Decay(const Helicity helicity, const bool radiative, const bool spincorr)
: helicity_(helicity) {
// if TAUOLA has not already been initialized,
if (!Tpp::Tauola::getIsTauolaIni()) {
// set the default units explicitly.
// we use GeV from momentum and cm for length
Tpp::Tauola::setUnits(Tpp::Tauola::GEV, Tpp::Tauola::CM);
// we set the lifetime as zero so that the tau decays
// at the vertex point given by CORSIKA since
// CORSIKA handles proapgating the tau until the decay point.
Tpp::Tauola::setTauLifetime(0.0);
// don't decay any intermediaries - add these to the stack
// in particular, eta's, K0, Pi-0's.
// See Pg. 46 in the TAUOLA manual.
Tpp::Tauola::setEtaK0sPi(1, 1, 0);
// use the new currents (Pg. 40 in the user manual).
Tpp::Tauola::setNewCurrents(1);
// whether to use QCD radiative corrections in calculating
// leptonic tau decay. If this is enabled, TAUOLA generates
// an order of magnitude more photons than Pythia (for example)
Tpp::Tauola::setRadiation(radiative);
// we have to set the spin correlation state after
// we have initialized TAUOLA
Tpp::Tauola::spin_correlation.setAll(spincorr);
// initialize the TAUOLA library
Tpp::Tauola::initialize();
}
}
Decay::~Decay() {
CORSIKA_LOGGER_INFO(logger_, "Number of particles decayed using TAUOLA: {}", count_);
}
bool Decay::isDecayHandled(const Code code) {
// if known to tauola, it can decay
if (code == Code::WPlus || code == Code::WMinus || code == Code::Z0 ||
code == Code::TauPlus || code == Code::TauMinus || code == Code::Eta ||
code == Code::K0Short || code == Code::RhoPlus || code == Code::RhoMinus ||
code == Code::KStarPlus || code == Code::KStarMinus)
return true;
else
return false;
}
// perform the decay
template <typename TSecondaryView>
void Decay::doDecay(TSecondaryView& view) {
CORSIKA_LOGGER_DEBUG(logger_, "Decaying primary particle using TAUOLA...");
count_++;
// get the coordinate system of the decay
auto projectile = view.getProjectile();
const auto& cs{projectile.getMomentum().getCoordinateSystem()};
// get the TAUOLA particle from this projectile
auto particle{TauolaInterfaceParticle::from_projectile(projectile)};
CORSIKA_LOGGER_TRACE(
logger_, "[decay primary] code: {}, momentum: ({}, {}, {}) GeV, energy: {} GeV",
convert_from_PDG(static_cast<PDGCode>(particle->getPdgID())), particle->getPx(),
particle->getPy(), particle->getPz(), particle->getE());
// decay this particle using TAUOLA
decay(particle.get());
std::vector<TauolaInterfaceParticle*> daughters;
// we iterate over the daughter particles and check for
// short-lived secondaries that we have to decay again using TAUOLA.
// We only ever have to go one layer down in the hierarchy.
for (const auto& daughter : particle->getDaughters()) {
// check if this is decay can be handled by TAUOLA
if (abs(daughter->getPdgID()) == 20213 ||
isDecayHandled(convert_from_PDG(static_cast<PDGCode>(daughter->getPdgID())))) {
// decay this secondary
decay(daughter);
// and now add all the daughters of this particle to our list
for (const auto& child : daughter->getDaughters()) {
// check if this is further decayable (K0 and Eta's typically)
if (abs(child->getPdgID()) == 20213 ||
isDecayHandled(convert_from_PDG(static_cast<PDGCode>(child->getPdgID())))) {
// decay the cild
decay(child);
// and add it's daughters to the top-level daughters vector
for (const auto& grandchild : child->getDaughters()) {
daughters.push_back(dynamic_cast<TauolaInterfaceParticle*>(grandchild));
}
} else {
// if we are here, the original daughter is non-decayable so add it to our
// struct.
daughters.push_back(dynamic_cast<TauolaInterfaceParticle*>(child));
}
}
} else {
// if we get here, `daughter` is long-lived and doesn't require a decay.
daughters.push_back(dynamic_cast<TauolaInterfaceParticle*>(daughter));
}
}
// we now have a complete list of the daughter particles of this decay
// we now iterate over them and add them to CORSIKA stack.
for (const auto& daughter : daughters) {
// get particle's PDG ID
const auto pdg{daughter->getPdgID()};
// get the CORSIKA PID from our PDG ID
const auto pid{convert_from_PDG(static_cast<PDGCode>(pdg))};
// construct the momentum in the ROOT coordinate system
const MomentumVector P(cs, {daughter->getPx() * 1_GeV, daughter->getPy() * 1_GeV,
daughter->getPz() * 1_GeV});
HEPEnergyType const mass = get_mass(pid);
HEPEnergyType const Ekin = sqrt(P.getSquaredNorm() + mass * mass) - mass;
CORSIKA_LOGGER_TRACE(logger_,
"[decay result] code: {}, momentum: {} GeV, energy: {} GeV",
pid, P.getComponents() / 1_GeV, daughter->getE());
// and add the particle to the CORSIKA stack
projectile.addSecondary(std::make_tuple(pid, Ekin, P.normalized()));
}
}
template <typename InterfaceParticle>
void Decay::decay(InterfaceParticle* particle) {
// return if we get a nullptr here
// this should never happen (but to be safe!)
// LCOV_EXCL_START
if (particle == nullptr) {
CORSIKA_LOGGER_WARN(logger_, "Got nullptr in decay routine. Skipping decay.");
return;
}
// LCOV_EXCL_STOP
// get the particle's PDG
const auto pdg{static_cast<PDGCode>(particle->getPdgID())};
// if this is a tau- or tau+, use the helicity in the decay
if (pdg == PDGCode::TauMinus || pdg == PDGCode::TauPlus) {
// get helicity from enum class
double helic;
switch (helicity_) {
case Helicity::LeftHanded:
helic = -1.;
break;
case Helicity::RightHanded:
helic = 1.;
break;
case Helicity::Unpolarized:
helic = 0.;
break;
// NOT testing undefined
default: // LCOV_EXCL_START
CORSIKA_LOGGER_ERROR(logger_, "Unknown helicity setting: {}", helicity_);
throw std::runtime_error("Unknown helicity setting!");
// LCOV_EXCL_STOP
}
double Polx{0};
double Poly{0};
double Polz{-helic};
// and call the decay routine with the polarization
Tpp::Tauola::decayOne(particle, false, Polx, Poly, Polz);
} else { // this is not a tau, so just do a standard decay
Tpp::Tauola::decayOne(particle);
}
}
template <typename TParticle>
TimeType Decay::getLifetime(TParticle const& particle) {
// Can be moved into Cascade.inl, see Issue #271
auto const pid = particle.getPID();
if ((pid == Code::TauMinus) || (pid == Code::TauPlus)) {
HEPEnergyType E = particle.getEnergy();
HEPMassType m = particle.getMass();
double const gamma = E / m;
TimeType const t0 = get_lifetime(pid);
auto const lifetime = gamma * t0;
CORSIKA_LOGGER_TRACE(logger_,
"[decay time] code: {}, t0: {} ns, momentum: {} GeV, energy: "
"{} GeV, gamma: {}, tau: {} ns",
particle.getPID(), t0 / 1_ns,
particle.getMomentum().getComponents() / 1_GeV, E / 1_GeV,
gamma, lifetime / 1_ns);
return lifetime;
} else
return std::numeric_limits<double>::infinity() * 1_s;
}
void Decay::PrintSummary() const { Tpp::Tauola::summary(); }
} // namespace corsika::tauola
corsika/detail/modules/thinning/EMThinning.inl 0 → 100644
View file @ 00b7c3f3
/*
* (c) Copyright 2022 CORSIKA Project, corsika-project@lists.kit.edu
*
* This software is distributed under the terms of the 3-clause BSD license.
* See file LICENSE for a full version of the license.
*/
#pragma once
#include <corsika/framework/core/PhysicalUnits.hpp>
#include <corsika/framework/process/SecondariesProcess.hpp>
namespace corsika {
EMThinning::EMThinning(HEPEnergyType threshold, double maxWeight,
bool const eraseParticles)
: threshold_{threshold}
, maxWeight_{maxWeight}
, eraseParticles_{eraseParticles} {}
template <typename TStackView>
void EMThinning::doSecondaries(TStackView& view) {
if (view.getSize() != 2) return;
auto projectile = view.getProjectile();
if (!is_em(projectile.getPID())) return;
double const parentWeight = projectile.getWeight();
if (parentWeight >= maxWeight_) { return; }
auto particle1 = view.begin();
auto particle2 = ++(view.begin());
// skip non-EM interactions
if (!is_em(particle1.getPID()) || !is_em(particle2.getPID())) { return; }
/* skip high-energy interactions
* We consider only the primary energy. A more sophisticated algorithm might
* sort out above-threshold secondaries and apply thinning only on the subset of
* those below the threshold.
*/
if (projectile.getEnergy() > threshold_) { return; }
corsika::HEPEnergyType const Esum = particle1.getEnergy() + particle2.getEnergy();
double const p1 = particle1.getEnergy() / Esum;
double const p2 = particle2.getEnergy() / Esum;
// weight factors
auto const f1 = 1 / p1;
auto const f2 = 1 / p2;
// max. allowed weight factor
double const fMax = maxWeight_ / parentWeight;
if (f1 <= fMax && f2 <= fMax) {
/* cannot possibly reach wmax in this vertex; apply
Hillas thinning; select one of the two */
if (uniform_(rng_) <= p1) { // keep 1st with probability p1
particle2.setWeight(0);
particle1.setWeight(parentWeight * f1);
} else { // keep 2nd
particle1.setWeight(0);
particle2.setWeight(parentWeight * f2);
}
} else { // weight limitation kicks in, do statistical thinning
double const f1prime = std::min(f1, fMax);
double const f2prime = std::min(f2, fMax);
if (uniform_(rng_) * f1prime <= 1) { // keep 1st
particle1.setWeight(parentWeight * f1prime);
} else {
particle1.setWeight(0);
}
if (uniform_(rng_) * f2prime <= 1) { // keep 2nd
particle2.setWeight(parentWeight * f2prime);
} else {
particle2.setWeight(0);
}
}
// erase discared particles in case of multithinning
if (eraseParticles_) {
for (auto& p : view) {
if (auto const w = p.getWeight(); w == 0) { p.erase(); }
}
} else {
return;
}
}
} // namespace corsika
corsika/detail/modules/tracking/Intersect.inl
View file @ 00b7c3f3
/*
* (c) Copyright 2020 CORSIKA Project, corsika-project@lists.kit.edu
*
* 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.
* This software is distributed under the terms of the 3-clause BSD license.
* See file LICENSE for a full version of the license.
*/
#pragma once
......@@ -44,9 +43,10 @@ namespace corsika {
fmt::ptr(&volumeNode), fmt::ptr(volumeNode.getParent()),
time_intersections_curr.hasIntersections());
if (time_intersections_curr.hasIntersections()) {
CORSIKA_LOG_DEBUG("intersection times with currentLogicalVolumeNode: {} s and {} s",
time_intersections_curr.getEntry() / 1_s,
time_intersections_curr.getExit() / 1_s);
CORSIKA_LOG_DEBUG(
"intersection times with currentLogicalVolumeNode: entry={} s and exit={} s",
time_intersections_curr.getEntry() / 1_s,
time_intersections_curr.getExit() / 1_s);
if (time_intersections_curr.getExit() <= minTime) {
minTime =
time_intersections_curr.getExit(); // we exit currentLogicalVolumeNode here
......@@ -59,15 +59,12 @@ namespace corsika {
for (auto const& node : volumeNode.getChildNodes()) {
Intersections const time_intersections = TDerived::intersect(particle, *node);
CORSIKA_LOG_TRACE("intersection times with child volume {}", fmt::ptr(node));
CORSIKA_LOG_DEBUG("search intersection times with child volume {}:",
fmt::ptr(node));
if (!time_intersections.hasIntersections()) { continue; }
CORSIKA_LOG_TRACE(" : enter {} s, exit {} s",
time_intersections.getEntry() / 1_s,
time_intersections.getExit() / 1_s);
auto const t_entry = time_intersections.getEntry();
auto const t_exit = time_intersections.getExit();
CORSIKA_LOG_TRACE("children t-entry: {}, t-exit: {}, smaller? {} ", t_entry, t_exit,
CORSIKA_LOG_DEBUG("children t-entry: {}, t-exit: {}, smaller? {} ", t_entry, t_exit,
t_entry <= minTime);
// note, theoretically t can even be smaller than 0 since we
// KNOW we can't yet be in this volume yet, so we HAVE TO
......@@ -95,7 +92,7 @@ namespace corsika {
time_intersections.getExit() / 1_s);
auto const t_entry = time_intersections.getEntry();
auto const t_exit = time_intersections.getExit();
CORSIKA_LOG_TRACE("children t-entry: {}, t-exit: {}, smaller? {} ", t_entry, t_exit,
CORSIKA_LOG_DEBUG("children t-entry: {}, t-exit: {}, smaller? {} ", t_entry, t_exit,
t_entry <= minTime);
// note, theoretically t can even be smaller than 0 since we
// KNOW we can't yet be in this volume yet, so we HAVE TO
......@@ -105,7 +102,18 @@ namespace corsika {
minNode = node;
}
}
CORSIKA_LOG_TRACE("t-intersect: {}, node {} ", minTime, fmt::ptr(minNode));
CORSIKA_LOG_DEBUG("next time-intersect: {}, node {} ", minTime, fmt::ptr(minNode));
// this branch cannot be unit-testes. This is malfunction: LCOV_EXCL_START
if (minTime < 0_s) {
if (minTime < -1e-8_s) {
CORSIKA_LOG_DEBUG(
"There is a very negative time step detected: {}. This is not physical and "
"may "
"easily crash subsequent modules. Set to 0_s, but CHECK AND FIX.",
minTime);
}
minTime = 0_s; // LCOV_EXCL_STOP
}
return std::make_tuple(minTime, minNode);
}
......
corsika/detail/modules/tracking/TrackingLeapFrogCurved.inl
View file @ 00b7c3f3
/*
* (c) Copyright 2020 CORSIKA Project, corsika-project@lists.kit.edu
*
* 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.
* This software is distributed under the terms of the 3-clause BSD license.
* See file LICENSE for a full version of the license.
*/
#pragma once
......@@ -28,33 +27,6 @@ namespace corsika {
namespace tracking_leapfrog_curved {
template <typename TParticle>
inline auto make_LeapFrogStep(TParticle const& particle, LengthType steplength) {
if (particle.getMomentum().getNorm() == 0_GeV) {
return std::make_tuple(particle.getPosition(), particle.getMomentum() / 1_GeV,
double(0));
} // charge of the particle
ElectricChargeType const charge = particle.getCharge();
auto const* currentLogicalVolumeNode = particle.getNode();
MagneticFieldVector const& magneticfield =
currentLogicalVolumeNode->getModelProperties().getMagneticField(
particle.getPosition());
VelocityVector velocity = particle.getVelocity();
decltype(meter / (second * volt)) k =
charge * constants::cSquared * 1_eV /
(velocity.getNorm() * particle.getEnergy() * 1_V);
DirectionVector direction = velocity.normalized();
auto position = particle.getPosition(); // First Movement
// assuming magnetic field does not change during movement
position =
position + direction * steplength / 2; // Change of direction by magnetic field
direction =
direction + direction.cross(magneticfield) * steplength * k; // Second Movement
position = position + direction * steplength / 2;
auto const steplength_true = steplength * (1 + direction.getNorm()) / 2;
return std::make_tuple(position, direction.normalized(), steplength_true);
}
template <typename TParticle>
inline auto Tracking::getTrack(TParticle const& particle) {
VelocityVector const initialVelocity = particle.getVelocity();
......@@ -94,6 +66,8 @@ namespace corsika {
particle.getMomentum().getParallelProjectionOnto(magneticfield))
.getNorm();
CORSIKA_LOG_TRACE("p_perp={} eV", p_perp / 1_eV);
if (p_perp < 1_eV) {
// particle travel along, parallel to magnetic field. Rg is
// "0", but for purpose of step limit we return infinity here.
......@@ -104,8 +78,18 @@ namespace corsika {
LengthType const gyroradius = convert_HEP_to_SI<MassType::dimension_type>(p_perp) *
constants::c / (abs(charge) * magnitudeB);
double const maxRadians = 0.01; // maximal allowed deflection
LengthType const steplimit = 2 * cos(maxRadians) * sin(maxRadians) * gyroradius;
if (gyroradius > 1e9_m) {
// this cannot be really unit-tested. It is hidden. LCOV_EXCL_START
CORSIKA_LOG_TRACE(
"CurvedLeapFrog is not very stable for extremely high gyroradius steps. "
"Rg={} -> straight tracking.",
gyroradius);
return getLinearTrajectory(particle);
// LCOV_EXCL_STOP
}
LengthType const steplimit = 2 * cos(maxMagneticDeflectionAngle_) *
sin(maxMagneticDeflectionAngle_) * gyroradius;
TimeType const steplimit_time = steplimit / initialVelocity.getNorm();
CORSIKA_LOG_DEBUG("gyroradius {}, steplimit: {} = {}", gyroradius, steplimit,
steplimit_time);
......@@ -121,19 +105,20 @@ namespace corsika {
decltype(1 / (tesla * second)) const k =
charge / p_norm *
initialVelocity.getNorm(); // since we use steps in time and not length
// units: C * s / m / kg * m/s = 1 / (T*m) * m/s = 1/(T s)
// units: C * s / m / kg * m/s = 1 / (T*m) * m/s = 1 / (T*s)
return std::make_tuple(
LeapFrogTrajectory(position, initialVelocity, magneticfield, k,
minTime), // trajectory
minNode); // next volume node
minTime), // --> trajectory
minNode); // --> next volume node
}
template <typename TParticle>
inline Intersections Tracking::intersect(TParticle const& particle,
Sphere const& sphere) {
if (sphere.getRadius() == 1_km * std::numeric_limits<double>::infinity()) {
LengthType const radius = sphere.getRadius();
if (radius == 1_km * std::numeric_limits<double>::infinity()) {
return Intersections();
}
......@@ -152,39 +137,66 @@ namespace corsika {
CORSIKA_LOG_TRACE("projectedDirectionSqrNorm={} T^2",
projectedDirectionSqrNorm / square(1_T));
if ((charge == 0 * constants::e) || magneticfield.getNorm() == 0_T || isParallel) {
if (isParallel) {
// particle moves parallel to field -> no deflection
return tracking_line::Tracking::intersect<TParticle>(particle, sphere);
}
bool const numericallyInside = sphere.contains(particle.getPosition());
CORSIKA_LOG_TRACE("numericallyInside={}", numericallyInside);
Vector<length_d> const deltaPos = position - sphere.getCenter();
{ // check extreme cases we don't want to solve analytically explicit
HEPMomentumType const p_perp =
(particle.getMomentum() -
particle.getMomentum().getParallelProjectionOnto(magneticfield))
.getNorm();
LengthType const gyroradius =
(convert_HEP_to_SI<MassType::dimension_type>(p_perp) * constants::c /
(abs(charge) * magneticfield.getNorm()));
LengthType const trackDist = abs(deltaPos.getNorm() - radius);
if (trackDist > gyroradius) {
// there is never a solution
return Intersections();
}
if (gyroradius > 1000 * trackDist) {
// the bending is negligible, use straight intersections instead
return tracking_line::Tracking::intersect(particle, sphere);
}
}
SpeedType const absVelocity = velocity.getNorm();
auto const p_norm =
constants::c * convert_HEP_to_SI<MassType::dimension_type>(
particle.getMomentum().getNorm()); // km * m /s
// this is: k = q/|p|
auto const k = charge / p_norm;
decltype(1 / (tesla * meter)) const k = charge / p_norm;
MagneticFieldVector const direction_x_B = directionBefore.cross(magneticfield);
auto const denom = 4. / (direction_x_B.getSquaredNorm() * k * k);
Vector<length_d> const deltaPos = position - sphere.getCenter();
double const b = (direction_x_B.dot(deltaPos) * k + 1) * denom / (1_m * 1_m);
double const c = directionBefore.dot(deltaPos) * 2 * denom / (1_m * 1_m * 1_m);
LengthType const deltaPosLength = deltaPos.getNorm();
double const d = (deltaPosLength + sphere.getRadius()) *
(deltaPosLength - sphere.getRadius()) * denom /
double const d = (deltaPosLength + radius) * (deltaPosLength - radius) * denom /
(1_m * 1_m * 1_m * 1_m);
CORSIKA_LOG_TRACE("denom={}, b={}, c={}, d={}", denom, b, c, d);
// solutions of deltaL are obtained from quartic equation. Note, deltaL/2 is the
// length of each half step, however, the second half step is slightly longer
// because of the non-conservation of norm/velocity.
// The leap-frog length L is deltaL/2 * (1+|u_{n+1}|)
std::vector<double> solutions = solve_quartic_real(1, 0, b, c, d);
if (!solutions.size()) {
return tracking_line::Tracking::intersect<TParticle>(particle, sphere);
}
if (!solutions.size()) { return Intersections(); }
LengthType d_enter, d_exit;
int first = 0, first_entry = 0, first_exit = 0;
for (auto solution : solutions) {
for (auto const solution : solutions) {
LengthType const dist = solution * 1_m;
CORSIKA_LOG_TRACE("Solution (real) for current Volume: {} ", dist);
CORSIKA_LOG_TRACE(
"Solution (real) for current Volume: deltaL/2*2={} (deltaL/2*2/v={}) ", dist,
dist / absVelocity);
if (numericallyInside) {
// there must be an entry (negative) and exit (positive) solution
if (dist < 0.0001_m) { // security margin to assure
......@@ -242,6 +254,8 @@ namespace corsika {
first_entry, first_exit);
return Intersections();
}
// return in units of time
return Intersections(d_enter / absVelocity, d_exit / absVelocity);
}
......@@ -253,7 +267,7 @@ namespace corsika {
ElectricChargeType const charge = particle.getCharge();
if (charge != 0 * constants::e) {
if (charge != ElectricChargeType::zero()) {
auto const* currentLogicalVolumeNode = particle.getNode();
VelocityVector const velocity = particle.getVelocity();
......@@ -284,6 +298,10 @@ namespace corsika {
plane.getNormal().dot(position - plane.getCenter())) // unit: kg*m/s *m
/ (1_m * 1_m * 1_kg) * 1_s;
// deltaL from quadratic solution return half-step length deltaL/2 for leap-frog
// algorithmus. Note, the leap-frog length L is longer by (1+|u_{n_1}|)/2 because
// the direction norm of the second half step is >1.
std::vector<double> const deltaLs = solve_quadratic_real(denom, p, q);
CORSIKA_LOG_TRACE("deltaLs=[{}]", fmt::join(deltaLs, ", "));
......@@ -310,15 +328,11 @@ namespace corsika {
return Intersections(std::numeric_limits<double>::infinity() * 1_s);
}
CORSIKA_LOG_TRACE("maxStepLength={}", maxStepLength);
CORSIKA_LOG_TRACE("maxStepLength={} s", maxStepLength / 1_s);
// with final length correction
auto const corr =
(direction + direction_x_B * maxStepLength * charge / (p_norm * 2))
.getNorm() /
absVelocity; // unit: s/m
// with final length correction, |direction| becomes >1 during step
return Intersections(maxStepLength * corr); // unit: s
return Intersections(maxStepLength / absVelocity); // unit: s
} // no charge
......@@ -329,13 +343,23 @@ namespace corsika {
return tracking_line::Tracking::intersect(particle, plane);
}
template <typename TParticle>
inline Intersections Tracking::intersect(TParticle const& particle,
SeparationPlane const& sepPlane) {
return intersect(particle, sepPlane.getPlane());
}
template <typename TParticle, typename TBaseNodeType>
inline Intersections Tracking::intersect(TParticle const& particle,
TBaseNodeType const& volumeNode) {
Sphere const* sphere = dynamic_cast<Sphere const*>(&volumeNode.getVolume());
if (sphere) { return intersect(particle, *sphere); }
SeparationPlane const* sepPlane =
dynamic_cast<SeparationPlane const*>(&volumeNode.getVolume());
if (sepPlane) { return intersect(particle, *sepPlane); }
throw std::runtime_error(
"The Volume type provided is not supported in intersect(particle, node)");
"The Volume type provided is not supported in "
"TrackingLeapFrogCurved::intersect(particle, node)");
}
template <typename TParticle>
......
corsika/detail/modules/tracking/TrackingLeapFrogStraight.inl
View file @ 00b7c3f3
/*
* (c) Copyright 2020 CORSIKA Project, corsika-project@lists.kit.edu
*
* 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.
* This software is distributed under the terms of the 3-clause BSD license.
* See file LICENSE for a full version of the license.
*/
#pragma once
......@@ -63,7 +62,7 @@ namespace corsika {
// check, where the first halve-step direction has geometric intersections
auto const [initialTrack, initialTrackNextVolume] =
tracking_line::Tracking::getTrack(particle);
{ [[maybe_unused]] auto& initialTrackNextVolume_dum = initialTrackNextVolume; }
//{ [[maybe_unused]] auto& initialTrackNextVolume_dum = initialTrackNextVolume; }
auto const initialTrackLength = initialTrack.getLength(1);
CORSIKA_LOG_DEBUG("initialTrack(0)={}, initialTrack(1)={}, initialTrackLength={}",
......@@ -82,8 +81,8 @@ namespace corsika {
.getNorm();
if (p_perp == 0_GeV) {
// particle travel along, parallel to magnetic field. Rg is
// "0", but for purpose of step limit we return infinity here.
// particle travels along, parallel to magnetic field. Rg is
// "0", return straight track here.
CORSIKA_LOG_TRACE("p_perp is 0_GeV --> parallel");
return std::make_tuple(initialTrack, initialTrackNextVolume);
}
......@@ -100,37 +99,38 @@ namespace corsika {
CORSIKA_LOG_DEBUG("gyroradius {}, Steplimit: {}", gyroradius, steplimit);
// calculate first halve step for "steplimit"
auto const initialMomentum = particle.getMomentum();
auto const absMomentum = initialMomentum.getNorm();
DirectionVector const direction = initialVelocity.normalized();
DirectionVector const initialDirection = particle.getDirection();
// avoid any intersections within first halve steplength
// aim for intersection in second halve step
LengthType const firstHalveSteplength =
std::min(steplimit, initialTrackLength * firstFraction_);
std::min(steplimit / 2, initialTrackLength * firstFraction_);
CORSIKA_LOG_DEBUG("first halve step length {}, steplimit={}, initialTrackLength={}",
firstHalveSteplength, steplimit, initialTrackLength);
// perform the first halve-step
Point const position_mid = initialPosition + direction * firstHalveSteplength;
Point const position_mid =
initialPosition + initialDirection * firstHalveSteplength;
auto const k = charge / (constants::c * convert_HEP_to_SI<MassType::dimension_type>(
particle.getMomentum().getNorm()));
DirectionVector const new_direction =
direction + direction.cross(magneticfield) * firstHalveSteplength * 2 * k;
initialDirection +
initialDirection.cross(magneticfield) * firstHalveSteplength * 2 * k;
auto const new_direction_norm = new_direction.getNorm(); // by design this is >1
CORSIKA_LOG_DEBUG(
"position_mid={}, new_direction={}, (new_direction_norm)={}, deflection={}",
position_mid.getCoordinates(), new_direction.getComponents(),
new_direction_norm,
acos(std::min(1.0, direction.dot(new_direction) / new_direction_norm)) * 180 /
M_PI);
acos(std::min(1.0, initialDirection.dot(new_direction) / new_direction_norm)) *
180 / M_PI);
// check, where the second halve-step direction has geometric intersections
particle.setPosition(position_mid);
particle.setMomentum(new_direction * absMomentum);
particle.setDirection(new_direction);
auto const [finalTrack, finalTrackNextVolume] =
tracking_line::Tracking::getTrack(particle);
particle.setPosition(initialPosition); // this is not nice...
particle.setMomentum(initialMomentum); // this is not nice...
particle.setPosition(initialPosition); // this is not nice...
particle.setDirection(initialDirection); // this is not nice...
LengthType const finalTrackLength = finalTrack.getLength(1);
LengthType const secondLeapFrogLength = firstHalveSteplength * new_direction_norm;
......
corsika/detail/modules/tracking/TrackingStraight.inl
View file @ 00b7c3f3
/*
* (c) Copyright 2020 CORSIKA Project, corsika-project@lists.kit.edu
*
* 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.
* This software is distributed under the terms of the 3-clause BSD license.
* See file LICENSE for a full version of the license.
*/
#pragma once
#include <corsika/framework/core/Logging.hpp>
#include <corsika/framework/core/PhysicalUnits.hpp>
#include <corsika/framework/geometry/Intersections.hpp>
#include <corsika/framework/geometry/Line.hpp>
#include <corsika/framework/geometry/Plane.hpp>
#include <corsika/framework/geometry/Sphere.hpp>
#include <corsika/framework/geometry/Vector.hpp>
#include <corsika/framework/geometry/StraightTrajectory.hpp>
#include <corsika/framework/geometry/Intersections.hpp>
#include <corsika/framework/core/Logging.hpp>
#include <corsika/framework/geometry/Vector.hpp>
#include <corsika/modules/tracking/Intersect.hpp>
#include <cmath>
#include <type_traits>
#include <utility>
#include <cmath>
namespace corsika::tracking_line {
......@@ -72,13 +71,68 @@ namespace corsika::tracking_line {
return Intersections();
}
template <typename TParticle>
inline Intersections Tracking::intersect(TParticle const& particle, Box const& box) {
Point const& position = particle.getPosition();
VelocityVector const velocity =
particle.getMomentum() / particle.getEnergy() * constants::c;
CoordinateSystemPtr const& cs = box.getCoordinateSystem();
LengthType x0 = position.getX(cs);
LengthType y0 = position.getY(cs);
LengthType z0 = position.getZ(cs);
SpeedType vx = velocity.getX(cs);
SpeedType vy = velocity.getY(cs);
SpeedType vz = velocity.getZ(cs);
CORSIKA_LOG_TRACE(
"particle in box coordinate: position: ({:.3f}, {:.3f}, "
"{:.3f}) m, veolocity: ({:.3f}, {:.3f}, {:.3f}) m/ns",
x0 / 1_m, y0 / 1_m, z0 / 1_m, vx / (1_m / 1_ns), vy / (1_m / 1_ns),
vz / (1_m / 1_ns));
auto get_intersect_min_max = [](LengthType x0, SpeedType v0, LengthType dx) {
auto t1 = (dx - x0) / v0;
auto t2 = (-dx - x0) / v0;
if (t1 > t2)
return std::make_pair(t1, t2);
else
return std::make_pair(t2, t1);
};
auto [tx_max, tx_min] = get_intersect_min_max(x0, vx, box.getX());
auto [ty_max, ty_min] = get_intersect_min_max(y0, vy, box.getY());
auto [tz_max, tz_min] = get_intersect_min_max(z0, vz, box.getZ());
TimeType t_exit = std::min(std::min(tx_max, ty_max), tz_max);
TimeType t_enter = std::max(std::max(tx_min, ty_min), tz_min);
CORSIKA_LOG_DEBUG("t_enter: {} ns, t_exit: {} ns", t_enter / 1_ns, t_exit / 1_ns);
if ((t_exit > t_enter)) {
if (t_enter < 0_s && t_exit > 0_s)
CORSIKA_LOG_DEBUG("numericallyInside={}", box.contains(position));
else if (t_enter < 0_s && t_exit < 0_s)
CORSIKA_LOG_DEBUG("oppisite direction");
return Intersections(std::move(t_enter), std::move(t_exit));
} else
return Intersections();
}
template <typename TParticle, typename TBaseNodeType>
inline Intersections Tracking::intersect(TParticle const& particle,
TBaseNodeType const& volumeNode) {
Sphere const* sphere = dynamic_cast<Sphere const*>(&volumeNode.getVolume());
if (sphere) { return Tracking::intersect<TParticle>(particle, *sphere); }
throw std::runtime_error(
"The Volume type provided is not supported in Intersect(particle, node)");
if (Sphere const* sphere = dynamic_cast<Sphere const*>(&volumeNode.getVolume());
sphere) {
return Tracking::intersect<TParticle>(particle, *sphere);
} else if (Box const* box = dynamic_cast<Box const*>(&volumeNode.getVolume()); box) {
return Tracking::intersect<TParticle>(particle, *box);
} else if (SeparationPlane const* sepPlane =
dynamic_cast<SeparationPlane const*>(&volumeNode.getVolume());
sepPlane) {
return Tracking::intersect<TParticle>(particle, *sepPlane);
} else {
throw std::runtime_error(
"The Volume type provided is not supported in "
"TrackingStraight::intersect(particle, node)");
}
}
template <typename TParticle>
......@@ -92,9 +146,16 @@ namespace corsika::tracking_line {
CORSIKA_LOG_TRACE("n_dot_v={}, delta={}, momentum={}", n_dot_v, delta,
particle.getMomentum());
return Intersections(n_dot_v.magnitude() == 0
? std::numeric_limits<TimeType::value_type>::infinity() * 1_s
: n.dot(delta) / n_dot_v);
if (n_dot_v.magnitude() == 0)
return Intersections();
else
return Intersections(n.dot(delta) / n_dot_v);
}
template <typename TParticle>
inline Intersections Tracking::intersect(TParticle const& particle,
SeparationPlane const& sepPlane) {
return intersect(particle, sepPlane.getPlane());
}
} // namespace corsika::tracking_line
corsika/detail/modules/urqmd/ParticleConversion.inl 0 → 100644
View file @ 00b7c3f3
/*
* (c) Copyright 2020 CORSIKA Project, corsika-project@lists.kit.edu
*
* This software is distributed under the terms of the 3-clause BSD license.
* See file LICENSE for a full version of the license.
*/
#pragma once
#include <corsika/modules/urqmd/ParticleConversion.hpp>
#include <corsika/framework/core/ParticleProperties.hpp>
#include <corsika/framework/core/PhysicalUnits.hpp>
#include <urqmd.hpp>
namespace corsika::urqmd {
inline bool canInteract(Code const vCode) {
// According to the manual, UrQMD can use all mesons, baryons and nucleons
// which are modeled also as input particles. I think it is safer to accept
// only the usual long-lived species as input.
// TODO: Charmed mesons should be added to the list, too
static std::array constexpr validProjectileCodes{
Code::Proton, Code::AntiProton, Code::Neutron, Code::AntiNeutron, Code::PiPlus,
Code::PiMinus, Code::KPlus, Code::KMinus, Code::K0Short, Code::K0Long};
return std::find(std::cbegin(validProjectileCodes), std::cend(validProjectileCodes),
vCode) != std::cend(validProjectileCodes);
}
inline std::pair<int, int> convertToUrQMD(Code const code) {
static const std::map<int, std::pair<int, int>> mapPDGToUrQMD{
// data mostly from github.com/afedynitch/ParticleDataTool
{22, {100, 0}}, // photon
{111, {101, 0}}, // pi0
{211, {101, 2}}, // pi+
{-211, {101, -2}}, // pi-
{321, {106, 1}}, // K+
{-321, {-106, -1}}, // K-
{311, {106, -1}}, // K0
{-311, {-106, 1}}, // K0bar
{2212, {1, 1}}, // p
{2112, {1, -1}}, // n
{-2212, {-1, -1}}, // pbar
{-2112, {-1, 1}}, // nbar
{221, {102, 0}}, // eta
{213, {104, 2}}, // rho+
{-213, {104, -2}}, // rho-
{113, {104, 0}}, // rho0
{323, {108, 2}}, // K*+
{-323, {108, -2}}, // K*-
{313, {108, 0}}, // K*0
{-313, {-108, 0}}, // K*0-bar
{223, {103, 0}}, // omega
{333, {109, 0}}, // phi
{3222, {40, 2}}, // Sigma+
{3212, {40, 0}}, // Sigma0
{3112, {40, -2}}, // Sigma-
{3322, {49, 0}}, // Xi0
{3312, {49, -1}}, // Xi-
{3122, {27, 0}}, // Lambda0
{2224, {17, 4}}, // Delta++
{2214, {17, 2}}, // Delta+
{2114, {17, 0}}, // Delta0
{1114, {17, -2}}, // Delta-
{3224, {41, 2}}, // Sigma*+
{3214, {41, 0}}, // Sigma*0
{3114, {41, -2}}, // Sigma*-
{3324, {50, 0}}, // Xi*0
{3314, {50, -1}}, // Xi*-
{3334, {55, 0}}, // Omega-
{411, {133, 2}}, // D+
{-411, {133, -2}}, // D-
{421, {133, 0}}, // D0
{-421, {-133, 0}}, // D0-bar
{441, {107, 0}}, // etaC
{431, {138, 1}}, // Ds+
{-431, {138, -1}}, // Ds-
{433, {139, 1}}, // Ds*+
{-433, {139, -1}}, // Ds*-
{413, {134, 1}}, // D*+
{-413, {134, -1}}, // D*-
{10421, {134, 0}}, // D*0
{-10421, {-134, 0}}, // D*0-bar
{443, {135, 0}}, // jpsi
};
return mapPDGToUrQMD.at(static_cast<int>(get_PDG(code)));
}
inline Code convertFromUrQMD(int vItyp, int vIso3) {
int const pdgInt =
::urqmd::pdgid_(vItyp, vIso3); // use the conversion function provided by UrQMD
if (pdgInt == 0) { // ::urqmd::pdgid_ returns 0 on error
throw std::runtime_error("UrQMD pdgid() returned 0");
}
auto const pdg = static_cast<PDGCode>(pdgInt);
return convert_from_PDG(pdg);
}
} // namespace corsika::urqmd
corsika/detail/modules/urqmd/UrQMD.inl
View file @ 00b7c3f3
/*
* (c) Copyright 2020 CORSIKA Project, corsika-project@lists.kit.edu
*
* 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.
* This software is distributed under the terms of the 3-clause BSD license.
* See file LICENSE for a full version of the license.
*/
#pragma once
#include <corsika/modules/Random.hpp>
#include <corsika/modules/urqmd/UrQMD.hpp>
#include <corsika/modules/urqmd/ParticleConversion.hpp>
#include <corsika/framework/core/ParticleProperties.hpp>
#include <corsika/framework/core/PhysicalUnits.hpp>
#include <corsika/framework/core/EnergyMomentumOperations.hpp>
#include <corsika/framework/geometry/QuantityVector.hpp>
#include <corsika/framework/geometry/Vector.hpp>
#include <corsika/framework/utility/COMBoost.hpp>
#include <boost/filesystem.hpp>
#include <boost/filesystem/fstream.hpp>
#include <boost/multi_array.hpp>
#include <algorithm>
......@@ -29,17 +33,28 @@
namespace corsika::urqmd {
inline UrQMD::UrQMD(boost::filesystem::path xs_file) {
inline UrQMD::UrQMD(boost::filesystem::path xs_file, int const retryFlag)
: iflb_(retryFlag) {
readXSFile(xs_file);
corsika::connect_random_stream(RNG_, ::urqmd::set_rng_function);
::urqmd::iniurqmdc8_();
}
inline CrossSectionType UrQMD::getTabulatedCrossSection(Code projectileCode,
Code targetCode,
HEPEnergyType labEnergy) const {
inline bool UrQMD::isValid(Code const projectileId, Code const targetId) const {
if (!is_hadron(projectileId) || !corsika::urqmd::canInteract(projectileId)) {
return false;
}
if (!is_nucleus(targetId)) { return false; }
return true;
}
inline CrossSectionType UrQMD::getTabulatedCrossSection(
Code const projectileId, Code const targetId, HEPEnergyType const labEnergy) const {
// translated to C++ from CORSIKA 7 subroutine cxtot_u
auto const kinEnergy = labEnergy - get_mass(projectileCode);
auto const kinEnergy = labEnergy - get_mass(projectileId);
assert(kinEnergy >= HEPEnergyType::zero());
......@@ -53,7 +68,7 @@ namespace corsika::urqmd {
w[2 - 1] = w[2 - 1] - 2 * w[3 - 1];
int projectileIndex;
switch (projectileCode) {
switch (projectileId) {
case Code::Proton:
projectileIndex = 0;
break;
......@@ -87,13 +102,14 @@ namespace corsika::urqmd {
projectileIndex = 8;
break;
default: { // LCOV_EXCL_START since this can never happen due to canInteract
CORSIKA_LOG_WARN("UrQMD cross-section not tabulated for {}", projectileCode);
return CrossSectionType::zero(); // LCOV_EXCL_STOP
CORSIKA_LOG_WARN("UrQMD cross-section not tabulated for {}", projectileId);
return CrossSectionType::zero();
// LCOV_EXCL_STOP
}
}
int targetIndex;
switch (targetCode) {
switch (targetId) {
case Code::Nitrogen:
targetIndex = 0;
break;
......@@ -105,7 +121,7 @@ namespace corsika::urqmd {
break;
default:
std::stringstream ss;
ss << "UrQMD cross-section not tabluated for target " << targetCode;
ss << "UrQMD cross-section not tabluated for target " << targetId;
throw std::runtime_error(ss.str().data());
}
......@@ -117,53 +133,17 @@ namespace corsika::urqmd {
CORSIKA_LOG_TRACE(
"UrQMD::GetTabulatedCrossSection proj={}, targ={}, E={} GeV, sigma={}",
get_name(projectileCode), get_name(targetCode), labEnergy / 1_GeV, result);
get_name(projectileId), get_name(targetId), labEnergy / 1_GeV, result);
return result;
}
inline CrossSectionType UrQMD::getCrossSection(Code vProjectileCode, Code vTargetCode,
HEPEnergyType vLabEnergy,
int vAProjectile = 1) {
// the following is a translation of ptsigtot() into C++
if (vProjectileCode != Code::Nucleus &&
!is_nucleus(vTargetCode)) { // both particles are "special"
auto const mProj = get_mass(vProjectileCode);
auto const mTar = get_mass(vTargetCode);
double sqrtS =
sqrt(static_pow<2>(mProj) + static_pow<2>(mTar) + 2 * vLabEnergy * mTar) *
(1 / 1_GeV);
inline CrossSectionType UrQMD::getCrossSection(Code const projectileId,
Code const targetId,
FourMomentum const& projectileP4,
FourMomentum const& targetP4) const {
// we must set some UrQMD globals first...
auto const [ityp, iso3] = convertToUrQMD(vProjectileCode);
::urqmd::inputs_.spityp[0] = ityp;
::urqmd::inputs_.spiso3[0] = iso3;
auto const [itypTar, iso3Tar] = convertToUrQMD(vTargetCode);
::urqmd::inputs_.spityp[1] = itypTar;
::urqmd::inputs_.spiso3[1] = iso3Tar;
int one = 1;
int two = 2;
return ::urqmd::sigtot_(one, two, sqrtS) * 1_mb;
} else {
int const Ap = vAProjectile;
int const At = is_nucleus(vTargetCode) ? get_nucleus_A(vTargetCode) : 1;
double const maxImpact = ::urqmd::nucrad_(Ap) + ::urqmd::nucrad_(At) +
2 * ::urqmd::options_.CTParam[30 - 1];
return 10_mb * M_PI * static_pow<2>(maxImpact);
// is a constant cross-section really reasonable?
}
}
template <typename TParticle>
inline CrossSectionType UrQMD::getCrossSection(TParticle const& projectile,
Code targetCode) const {
auto const projectileCode = projectile.getPID();
if (projectileCode == Code::Nucleus) {
if (!isValid(projectileId, targetId)) {
/*
* unfortunately unavoidable at the moment until we have tools to get the actual
* inealstic cross-section from UrQMD
......@@ -171,72 +151,61 @@ namespace corsika::urqmd {
return CrossSectionType::zero();
}
return getTabulatedCrossSection(projectileCode, targetCode, projectile.getEnergy());
}
inline bool UrQMD::canInteract(Code vCode) const {
// According to the manual, UrQMD can use all mesons, baryons and nucleons
// which are modeled also as input particles. I think it is safer to accept
// only the usual long-lived species as input.
// TODO: Charmed mesons should be added to the list, too
// define projectile, in lab frame
auto const sqrtS2 = (projectileP4 + targetP4).getNormSqr();
HEPEnergyType const Elab = (sqrtS2 - static_pow<2>(get_mass(projectileId)) -
static_pow<2>(get_mass(targetId))) /
(2 * get_mass(targetId));
static std::array constexpr validProjectileCodes{
Code::Proton, Code::AntiProton, Code::Neutron, Code::AntiNeutron, Code::PiPlus,
Code::PiMinus, Code::KPlus, Code::KMinus, Code::K0Short, Code::K0Long};
bool const tabulated = true;
if (tabulated) { return getTabulatedCrossSection(projectileId, targetId, Elab); }
return std::find(std::cbegin(validProjectileCodes), std::cend(validProjectileCodes),
vCode) != std::cend(validProjectileCodes);
}
// the following is a translation of ptsigtot() into C++
if (!is_nucleus(projectileId) &&
!is_nucleus(targetId)) { // both particles are "special"
template <typename TParticle>
inline GrammageType UrQMD::getInteractionLength(TParticle const& vParticle) const {
double sqrtS = sqrt(sqrtS2) / 1_GeV;
if (!canInteract(vParticle.getPID())) {
// we could do the canInteract check in getCrossSection, too but if
// we do it here we have the advantage of avoiding the loop
return std::numeric_limits<double>::infinity() * 1_g / (1_cm * 1_cm);
}
// we must set some UrQMD globals first...
auto const [ityp, iso3] = corsika::urqmd::convertToUrQMD(projectileId);
::urqmd::inputs_.spityp[0] = ityp;
::urqmd::inputs_.spiso3[0] = iso3;
auto const& mediumComposition =
vParticle.getNode()->getModelProperties().getNuclearComposition();
using namespace std::placeholders;
auto const [itypTar, iso3Tar] = corsika::urqmd::convertToUrQMD(targetId);
::urqmd::inputs_.spityp[1] = itypTar;
::urqmd::inputs_.spiso3[1] = iso3Tar;
CrossSectionType const weightedProdCrossSection = mediumComposition.getWeightedSum(
std::bind(&UrQMD::getCrossSection<decltype(vParticle)>, this, vParticle, _1));
int one = 1;
int two = 2;
return ::urqmd::sigtot_(one, two, sqrtS) * 1_mb;
}
return mediumComposition.getAverageMassNumber() * constants::u /
weightedProdCrossSection;
// at least one of them is a nucleus
int const Ap = is_nucleus(projectileId) ? get_nucleus_A(projectileId) : 1;
int const At = is_nucleus(targetId) ? get_nucleus_A(targetId) : 1;
double const maxImpact = ::urqmd::nucrad_(Ap) + ::urqmd::nucrad_(At) +
2 * ::urqmd::options_.CTParam[30 - 1];
return 10_mb * M_PI * static_pow<2>(maxImpact);
// is a constant cross-section really reasonable?
}
template <typename TView>
inline void UrQMD::doInteraction(TView& view) {
auto projectile = view.getProjectile();
auto projectileCode = projectile.getPID();
auto const projectileEnergyLab = projectile.getEnergy();
auto const& projectileMomentumLab = projectile.getMomentum();
auto const& projectilePosition = projectile.getPosition();
auto const projectileTime = projectile.getTime();
// sample target particle
auto const& mediumComposition =
projectile.getNode()->getModelProperties().getNuclearComposition();
auto const componentCrossSections = std::invoke([&]() {
auto const& components = mediumComposition.getComponents();
std::vector<CrossSectionType> crossSections;
crossSections.reserve(components.size());
for (auto const c : components) {
crossSections.push_back(getCrossSection(projectile, c));
}
return crossSections;
});
inline void UrQMD::doInteraction(TView& view, Code const projectileId,
Code const targetId, FourMomentum const& projectileP4,
FourMomentum const& targetP4) {
// define projectile, in lab frame
auto const sqrtS2 = (projectileP4 + targetP4).getNormSqr();
HEPEnergyType const Elab = (sqrtS2 - static_pow<2>(get_mass(projectileId)) -
static_pow<2>(get_mass(targetId))) /
(2 * get_mass(targetId));
if (!isValid(projectileId, targetId)) {
throw std::runtime_error("invalid target,projectile,energy combination");
}
auto const targetCode = mediumComposition.sampleTarget(componentCrossSections, RNG_);
auto const targetA = get_nucleus_A(targetCode);
auto const targetZ = get_nucleus_Z(targetCode);
size_t const targetA = get_nucleus_A(targetId);
size_t const targetZ = get_nucleus_Z(targetId);
::urqmd::inputs_.nevents = 1;
::urqmd::sys_.eos = 0; // could be configurable in principle
......@@ -244,14 +213,13 @@ namespace corsika::urqmd {
::urqmd::sys_.nsteps = 1;
// initialization regarding projectile
if (Code::Nucleus == projectileCode) {
if (is_nucleus(projectileId)) {
// is this everything?
::urqmd::inputs_.prspflg = 0;
::urqmd::sys_.Ap = projectile.getNuclearA();
::urqmd::sys_.Zp = projectile.getNuclearZ();
::urqmd::rsys_.ebeam = (projectileEnergyLab - projectile.getMass()) * (1 / 1_GeV) /
projectile.getNuclearA();
::urqmd::sys_.Ap = get_nucleus_A(projectileId);
::urqmd::sys_.Zp = get_nucleus_Z(projectileId);
::urqmd::rsys_.ebeam = (Elab - get_mass(projectileId)) / 1_GeV / ::urqmd::sys_.Ap;
::urqmd::rsys_.bdist = ::urqmd::nucrad_(targetA) +
::urqmd::nucrad_(::urqmd::sys_.Ap) +
......@@ -265,20 +233,19 @@ namespace corsika::urqmd {
1; // even for non-baryons this has to be set, see vanilla UrQMD.f
::urqmd::rsys_.bdist = ::urqmd::nucrad_(targetA) + ::urqmd::nucrad_(1) +
2 * ::urqmd::options_.CTParam[30 - 1];
::urqmd::rsys_.ebeam = (projectileEnergyLab - projectile.getMass()) * (1 / 1_GeV);
::urqmd::rsys_.ebeam = (Elab - get_mass(projectileId)) / 1_GeV;
if (projectileCode == Code::K0Long || projectileCode == Code::K0Short) {
projectileCode = booleanDist_(RNG_) ? Code::K0 : Code::K0Bar;
}
auto const [ityp, iso3] = convertToUrQMD(projectileCode);
auto const [ityp, iso3] = corsika::urqmd::convertToUrQMD(
(projectileId == Code::K0Long || projectileId == Code::K0Short)
? (booleanDist_(RNG_) ? Code::K0 : Code::K0Bar)
: projectileId);
// todo: conversion of K_long/short into strong eigenstates;
::urqmd::inputs_.spityp[0] = ityp;
::urqmd::inputs_.spiso3[0] = iso3;
}
// initilazation regarding target
if (is_nucleus(targetCode)) {
// initialization regarding target
if (is_nucleus(targetId)) {
::urqmd::sys_.Zt = targetZ;
::urqmd::sys_.At = targetA;
::urqmd::inputs_.trspflg = 0; // nucleus as target
......@@ -286,118 +253,57 @@ namespace corsika::urqmd {
::urqmd::cascinit_(::urqmd::sys_.Zt, ::urqmd::sys_.At, id);
} else {
::urqmd::inputs_.trspflg = 1; // special particle as target
auto const [ityp, iso3] = convertToUrQMD(targetCode);
auto const [ityp, iso3] = corsika::urqmd::convertToUrQMD(targetId);
::urqmd::inputs_.spityp[1] = ityp;
::urqmd::inputs_.spiso3[1] = iso3;
}
int iflb = 0; // flag for retrying interaction in case of empty event, 0 means retry
int iflb =
iflb_; // flag for retrying interaction in case of empty event, 0 means retry
::urqmd::urqmd_(iflb);
// now retrieve secondaries from UrQMD
auto const& originalCS = projectileMomentumLab.getCoordinateSystem();
CoordinateSystemPtr const& zAxisFrame =
make_rotationToZ(originalCS, projectileMomentumLab);
COMBoost const boost(projectileP4, targetP4);
auto const& originalCS = boost.getOriginalCS();
auto const& csPrime = boost.getRotatedCS();
for (int i = 0; i < ::urqmd::sys_.npart; ++i) {
auto code = convertFromUrQMD(::urqmd::isys_.ityp[i], ::urqmd::isys_.iso3[i]);
auto code = corsika::urqmd::convertFromUrQMD(::urqmd::isys_.ityp[i],
::urqmd::isys_.iso3[i]);
if (code == Code::K0 || code == Code::K0Bar) {
code = booleanDist_(RNG_) ? Code::K0Short : Code::K0Long;
}
// "coor_.p0[i] * 1_GeV" is likely off-shell as UrQMD doesn't preserve masses well
auto momentum =
Vector(zAxisFrame,
QuantityVector<dimensionless_d>{
::urqmd::coor_.px[i], ::urqmd::coor_.py[i], ::urqmd::coor_.pz[i]} *
1_GeV);
MomentumVector momentum{csPrime,
{::urqmd::coor_.px[i] * 1_GeV, ::urqmd::coor_.py[i] * 1_GeV,
::urqmd::coor_.pz[i] * 1_GeV}};
momentum.rebase(originalCS); // transform back into standard lab frame
CORSIKA_LOG_DEBUG(" {} {} {} ", i, code, momentum.getComponents());
projectile.addSecondary(
std::make_tuple(code, momentum, projectilePosition, projectileTime));
HEPEnergyType const mass = get_mass(code);
HEPEnergyType const Ekin = calculate_kinetic_energy(momentum.getNorm(), mass);
if (Ekin <= 0_GeV) {
if (Ekin < 0_GeV) {
CORSIKA_LOG_WARN("Negative kinetic energy {} {}. Skipping.", code, Ekin);
}
view.addSecondary(
std::make_tuple(code, 0_eV, DirectionVector{originalCS, {0, 0, 0}}));
} else {
view.addSecondary(std::make_tuple(code, Ekin, momentum.normalized()));
}
}
CORSIKA_LOG_DEBUG("UrQMD generated {} secondaries!", ::urqmd::sys_.npart);
}
inline Code convertFromUrQMD(int vItyp, int vIso3) {
int const pdgInt =
::urqmd::pdgid_(vItyp, vIso3); // use the conversion function provided by UrQMD
if (pdgInt == 0) { // ::urqmd::pdgid_ returns 0 on error
throw std::runtime_error("UrQMD pdgid() returned 0");
}
auto const pdg = static_cast<PDGCode>(pdgInt);
return convert_from_PDG(pdg);
}
inline std::pair<int, int> convertToUrQMD(Code code) {
static const std::map<int, std::pair<int, int>> mapPDGToUrQMD{
// data mostly from github.com/afedynitch/ParticleDataTool
{22, {100, 0}}, // photon
{111, {101, 0}}, // pi0
{211, {101, 2}}, // pi+
{-211, {101, -2}}, // pi-
{321, {106, 1}}, // K+
{-321, {-106, -1}}, // K-
{311, {106, -1}}, // K0
{-311, {-106, 1}}, // K0bar
{2212, {1, 1}}, // p
{2112, {1, -1}}, // n
{-2212, {-1, -1}}, // pbar
{-2112, {-1, 1}}, // nbar
{221, {102, 0}}, // eta
{213, {104, 2}}, // rho+
{-213, {104, -2}}, // rho-
{113, {104, 0}}, // rho0
{323, {108, 2}}, // K*+
{-323, {108, -2}}, // K*-
{313, {108, 0}}, // K*0
{-313, {-108, 0}}, // K*0-bar
{223, {103, 0}}, // omega
{333, {109, 0}}, // phi
{3222, {40, 2}}, // Sigma+
{3212, {40, 0}}, // Sigma0
{3112, {40, -2}}, // Sigma-
{3322, {49, 0}}, // Xi0
{3312, {49, -1}}, // Xi-
{3122, {27, 0}}, // Lambda0
{2224, {17, 4}}, // Delta++
{2214, {17, 2}}, // Delta+
{2114, {17, 0}}, // Delta0
{1114, {17, -2}}, // Delta-
{3224, {41, 2}}, // Sigma*+
{3214, {41, 0}}, // Sigma*0
{3114, {41, -2}}, // Sigma*-
{3324, {50, 0}}, // Xi*0
{3314, {50, -1}}, // Xi*-
{3334, {55, 0}}, // Omega-
{411, {133, 2}}, // D+
{-411, {133, -2}}, // D-
{421, {133, 0}}, // D0
{-421, {-133, 0}}, // D0-bar
{441, {107, 0}}, // etaC
{431, {138, 1}}, // Ds+
{-431, {138, -1}}, // Ds-
{433, {139, 1}}, // Ds*+
{-433, {139, -1}}, // Ds*-
{413, {134, 1}}, // D*+
{-413, {134, -1}}, // D*-
{10421, {134, 0}}, // D*0
{-10421, {-134, 0}}, // D*0-bar
{443, {135, 0}}, // jpsi
};
return mapPDGToUrQMD.at(static_cast<int>(get_PDG(code)));
}
inline void UrQMD::readXSFile(boost::filesystem::path const filename) {
boost::filesystem::ifstream file(filename, std::ios::in);
if (!file.is_open()) {
throw std::runtime_error(
filename.native() +
" could not be opened."); // LCOV_EXCL_LINE since this is pointless to test
// LCOV_EXCL_START since this is pointless to test
throw std::runtime_error(filename.native() + " could not be opened.");
// LCOV_EXCL_STOP
}
std::string line;
......@@ -405,7 +311,7 @@ namespace corsika::urqmd {
std::getline(file, line);
std::stringstream ss(line);
char dummy;
char dummy; // this is '#'
int nTargets, nProjectiles, nSupports;
ss >> dummy >> nTargets >> nProjectiles >> nSupports;
......
corsika/detail/modules/writers/EnergyLossWriter.inl 0 → 100644
View file @ 00b7c3f3
/*
* (c) Copyright 2021 CORSIKA Project, corsika-project@lists.kit.edu
*
* This software is distributed under the terms of the 3-clause BSD license.
* See file LICENSE for a full version of the license.
*/
#pragma once
#include <corsika/framework/core/ParticleProperties.hpp>
#include <corsika/framework/core/PhysicalUnits.hpp>
#include <corsika/framework/utility/FindXmax.hpp>
#include <corsika/media/ShowerAxis.hpp>
#include <exception>
namespace corsika {
template <typename TOutput>
inline EnergyLossWriter<TOutput>::EnergyLossWriter(ShowerAxis const& axis,
GrammageType dX,
GrammageType dX_threshold)
: EnergyLossWriter<TOutput>{axis,
static_cast<unsigned int>(axis.getMaximumX() / dX) + 1,
dX, dX_threshold} {}
template <typename TOutput>
inline EnergyLossWriter<TOutput>::EnergyLossWriter(ShowerAxis const& axis,
unsigned int const nBins,
GrammageType dX,
GrammageType dX_threshold)
: TOutput(dEdX_output::ProfileIndexNames)
, showerAxis_(axis)
, dX_(dX)
, nBins_(nBins)
, dX_threshold_(dX_threshold) {}
template <typename TOutput>
inline void EnergyLossWriter<TOutput>::startOfLibrary(
boost::filesystem::path const& directory) {
TOutput::startOfLibrary(directory);
}
template <typename TOutput>
inline void EnergyLossWriter<TOutput>::startOfShower(unsigned int const showerId) {
TOutput::startOfShower(showerId);
// reset profile
profile_.clear();
profile_.resize(nBins_);
}
template <typename TOutput>
inline void EnergyLossWriter<TOutput>::endOfLibrary() {
TOutput::endOfLibrary();
}
template <typename TOutput>
inline void EnergyLossWriter<TOutput>::endOfShower(unsigned int const showerId) {
int iRow{0};
for (dEdX_output::Profile const& row : profile_) {
// here: write to underlying writer (e.g. parquet)
TOutput::write(showerId, iRow * dX_, row);
iRow++;
}
TOutput::endOfShower(showerId);
}
template <typename TOutput>
inline void EnergyLossWriter<TOutput>::write(Point const& p0, Point const& p1,
Code const PID, HEPEnergyType const dE) {
GrammageType grammageStart = showerAxis_.getProjectedX(p0);
GrammageType grammageEnd = showerAxis_.getProjectedX(p1);
if (grammageStart > grammageEnd) { // particle going upstream
std::swap(grammageStart, grammageEnd);
}
GrammageType const deltaX = grammageEnd - grammageStart;
CORSIKA_LOGGER_TRACE(
TOutput::getLogger(),
"dE={} GeV, grammageStart={} g/cm2, End={}g /cm2, deltaX={} g/cm2", dE / 1_GeV,
grammageStart / 1_g * square(1_cm), grammageEnd / 1_g * square(1_cm),
deltaX / 1_g * square(1_cm));
if (deltaX < dX_threshold_) {
CORSIKA_LOGGER_TRACE(TOutput::getLogger(), "Point-like dE");
this->write(p0, PID, dE);
return;
}
// only register the range that is covered by the profile
int const maxBin = int(profile_.size() - 1);
int binStart = grammageStart / dX_;
if (binStart < 0) binStart = 0;
if (binStart > maxBin) binStart = maxBin;
int binEnd = grammageEnd / dX_;
if (binEnd < 0) binEnd = 0;
if (binEnd > maxBin) binEnd = maxBin;
CORSIKA_LOGGER_TRACE(TOutput::getLogger(), "maxBin={}, binStart={}, binEnd={}",
maxBin, binStart, binEnd);
auto energyCount = HEPEnergyType::zero();
auto const factor = dE / deltaX; // [ energy / grammage ]
auto fill = [&](int const bin, GrammageType const weight) {
auto const increment = factor * weight;
CORSIKA_LOGGER_TRACE(TOutput::getLogger(),
"filling bin={} with weight {} : dE={} GeV ", bin, weight,
increment / 1_GeV);
profile_[bin][static_cast<int>(dEdX_output::ProfileIndex::Total)] += increment;
energyCount += increment;
};
// fill longitudinal profile
if (binStart == binEnd) {
fill(binStart, deltaX);
} else {
fill(binStart, ((1 + binStart) * dX_ - grammageStart));
fill(binEnd, (grammageEnd - binEnd * dX_));
for (int bin = binStart + 1; bin < binEnd; ++bin) { fill(bin, dX_); }
}
CORSIKA_LOGGER_TRACE(TOutput::getLogger(), "total energy added to histogram: {} GeV ",
energyCount / 1_GeV);
}
template <typename TOutput>
inline void EnergyLossWriter<TOutput>::write(Point const& point, Code const,
HEPEnergyType const dE) {
GrammageType grammage = showerAxis_.getProjectedX(point);
int const maxBin = int(profile_.size() - 1);
int bin = grammage / dX_;
if (bin < 0) bin = 0;
if (bin > maxBin) bin = maxBin;
CORSIKA_LOGGER_TRACE(TOutput::getLogger(), "add local energy loss bin={} dE={} GeV ",
bin, dE / 1_GeV);
profile_[bin][static_cast<int>(dEdX_output::ProfileIndex::Total)] += dE;
}
template <typename TOutput>
inline void EnergyLossWriter<TOutput>::write(GrammageType const Xstart,
GrammageType const Xend, Code const,
HEPEnergyType const dE) {
double const bstart = Xstart / dX_;
double const bend = Xend / dX_;
if (abs(bstart - floor(bstart + 0.5)) > 1e-2 ||
abs(bend - floor(bend + 0.5)) > 1e-2 || abs(bend - bstart - 1) > 1e-2) {
CORSIKA_LOGGER_ERROR(
TOutput::getLogger(),
"CascadeEquation (CONEX) and Corsika8 dX grammage binning are not the same! "
"Xstart={} Xend={} dX={} g/cm2",
Xstart / 1_g * square(1_cm), Xend / 1_g * square(1_cm),
dX_ / 1_g * square(1_cm));
throw std::runtime_error(
"CONEX and Corsika8 dX grammage binning are not the same!");
}
size_t const bin = size_t((bend + bstart) / 2);
CORSIKA_LOGGER_TRACE(TOutput::getLogger(),
"add binned energy loss {} {} bin={} dE={} GeV ", bstart, bend,
bin, dE / 1_GeV);
if (bin >= profile_.size()) {
CORSIKA_LOGGER_WARN(TOutput::getLogger(),
"Grammage bin {} outside of profile {}. skipping.", bin,
profile_.size());
return;
}
profile_[bin][static_cast<int>(dEdX_output::ProfileIndex::Total)] += dE;
}
template <typename TOutput>
inline HEPEnergyType EnergyLossWriter<TOutput>::getEnergyLost() const {
HEPEnergyType tot = HEPEnergyType::zero();
for (dEdX_output::Profile const& row : profile_)
tot += row.at(static_cast<int>(dEdX_output::ProfileIndex::Total));
return tot;
}
template <typename TOutput>
inline YAML::Node EnergyLossWriter<TOutput>::getConfig() const {
YAML::Node node;
node["type"] = "EnergyLoss";
node["units"]["energy"] = "GeV";
node["units"]["grammage"] = "g/cm^2";
node["bin-size"] = dX_ / (1_g / square(1_cm));
node["nbins"] = nBins_;
node["grammage_threshold"] = dX_threshold_ / (1_g / square(1_cm));
return node;
}
template <typename TOutput>
inline YAML::Node EnergyLossWriter<TOutput>::getSummary() const {
// determined Xmax and dEdXmax from quadratic interpolation
double maximum = 0;
size_t iMaximum = 0;
for (size_t i = 0; i < profile_.size() - 3; ++i) {
double value =
(profile_[i + 0].at(static_cast<int>(dEdX_output::ProfileIndex::Total)) +
profile_[i + 1].at(static_cast<int>(dEdX_output::ProfileIndex::Total)) +
profile_[i + 2].at(static_cast<int>(dEdX_output::ProfileIndex::Total))) /
1_GeV;
if (value > maximum) {
maximum = value;
iMaximum = i;
}
}
double const dX = dX_ / 1_g * square(1_cm);
auto [Xmax, dEdXmax] = FindXmax::interpolateProfile(
dX * (0.5 + iMaximum), dX * (1.5 + iMaximum), dX * (2.5 + iMaximum),
profile_[iMaximum + 0].at(static_cast<int>(dEdX_output::ProfileIndex::Total)) /
1_GeV,
profile_[iMaximum + 1].at(static_cast<int>(dEdX_output::ProfileIndex::Total)) /
1_GeV,
profile_[iMaximum + 2].at(static_cast<int>(dEdX_output::ProfileIndex::Total)) /
1_GeV);
YAML::Node summary;
summary["sum_dEdX"] = getEnergyLost() / 1_GeV;
summary["Xmax"] = Xmax;
summary["dEdXmax"] = dEdXmax;
return summary;
}
} // namespace corsika
\ No newline at end of file
corsika/detail/modules/writers/EnergyLossWriterParquet.inl 0 → 100644
View file @ 00b7c3f3
/*
* (c) Copyright 2021 CORSIKA Project, corsika-project@lists.kit.edu
*
* This software is distributed under the terms of the 3-clause BSD license.
* See file LICENSE for a full version of the license.
*/
#pragma once
#include <corsika/framework/core/ParticleProperties.hpp>
#include <corsika/framework/core/PhysicalUnits.hpp>
#include <corsika/framework/utility/FindXmax.hpp>
#include <corsika/media/ShowerAxis.hpp>
#include <exception>
namespace corsika {
template <size_t NColumns>
inline EnergyLossWriterParquet<NColumns>::EnergyLossWriterParquet(
std::array<const char*, NColumns> const& colNames)
: columns_(colNames) {}
template <size_t NColumns>
inline void EnergyLossWriterParquet<NColumns>::startOfLibrary(
boost::filesystem::path const& directory) {
// setup the streamer
output_.initStreamer((directory / "dEdX.parquet").string());
// enable compression with the default level
output_.enableCompression();
// build the schema
output_.addField("X", parquet::Repetition::REQUIRED, parquet::Type::FLOAT,
parquet::ConvertedType::NONE);
for (auto const& col : columns_) {
output_.addField(col, parquet::Repetition::REQUIRED, parquet::Type::FLOAT,
parquet::ConvertedType::NONE);
}
// and build the streamer
output_.buildStreamer();
}
template <size_t NColumns>
inline void EnergyLossWriterParquet<NColumns>::write(
unsigned int const showerId, GrammageType const grammage,
std::array<HEPEnergyType, NColumns> const& data) {
// and write the data into the column
*(output_.getWriter()) << showerId
<< static_cast<float>(grammage / 1_g * square(1_cm));
for (HEPEnergyType const& dedx : data) {
*(output_.getWriter()) << static_cast<float>(dedx / 1_GeV);
}
*(output_.getWriter()) << parquet::EndRow;
}
template <size_t NColumns>
inline void EnergyLossWriterParquet<NColumns>::startOfShower(unsigned int const) {}
template <size_t NColumns>
inline void EnergyLossWriterParquet<NColumns>::endOfShower(unsigned int const) {}
template <size_t NColumns>
inline void EnergyLossWriterParquet<NColumns>::endOfLibrary() {
output_.closeStreamer();
}
} // namespace corsika
corsika/detail/modules/writers/InteractionWriter.inl 0 → 100644
View file @ 00b7c3f3
/*
* (c) Copyright 2024 CORSIKA Project, corsika-project@lists.kit.edu
*
* This software is distributed under the terms of the 3-clause BSD license.
* See file LICENSE for a full version of the license.
*/
#pragma once
#include <corsika/framework/core/Logging.hpp>
namespace corsika {
template <typename TTracking, typename TOutput>
inline InteractionWriter<TTracking, TOutput>::InteractionWriter(
ShowerAxis const& axis, ObservationPlane<TTracking, TOutput> const& obsPlane)
: obsPlane_(obsPlane)
, showerAxis_(axis)
, interactionCounter_(0)
, showerId_(0) {}
template <typename TTracking, typename TOutput>
template <typename TStackView>
inline void InteractionWriter<TTracking, TOutput>::doSecondaries(TStackView& vS) {
// Dumps the stack to the output parquet stream
// The primary and secondaries are all written along with the slant depth
if (interactionCounter_++) { return; } // Only run on the first interaction
auto primary = vS.getProjectile();
auto dX = showerAxis_.getProjectedX(primary.getPosition());
CORSIKA_LOG_INFO("First interaction at dX {}", dX);
CORSIKA_LOG_INFO("Primary: {}, E_kin {}", primary.getPID(),
primary.getKineticEnergy());
// get the location of the primary w.r.t. observation plane
Vector const displacement = primary.getPosition() - obsPlane_.getPlane().getCenter();
auto const x = displacement.dot(obsPlane_.getXAxis());
auto const y = displacement.dot(obsPlane_.getYAxis());
auto const z = displacement.dot(obsPlane_.getPlane().getNormal());
auto const nx = primary.getDirection().dot(obsPlane_.getXAxis());
auto const ny = primary.getDirection().dot(obsPlane_.getYAxis());
auto const nz = primary.getDirection().dot(obsPlane_.getPlane().getNormal());
auto const px = primary.getMomentum().dot(obsPlane_.getXAxis());
auto const py = primary.getMomentum().dot(obsPlane_.getYAxis());
auto const pz = primary.getMomentum().dot(obsPlane_.getPlane().getNormal());
summary_["shower_" + std::to_string(showerId_)]["pdg"] =
static_cast<int>(get_PDG(primary.getPID()));
summary_["shower_" + std::to_string(showerId_)]["name"] =
static_cast<std::string>(get_name(primary.getPID()));
summary_["shower_" + std::to_string(showerId_)]["total_energy"] =
(primary.getKineticEnergy() + get_mass(primary.getPID())) / 1_GeV;
summary_["shower_" + std::to_string(showerId_)]["kinetic_energy"] =
primary.getKineticEnergy() / 1_GeV;
summary_["shower_" + std::to_string(showerId_)]["x"] = x / 1_m;
summary_["shower_" + std::to_string(showerId_)]["y"] = y / 1_m;
summary_["shower_" + std::to_string(showerId_)]["z"] = z / 1_m;
summary_["shower_" + std::to_string(showerId_)]["nx"] = static_cast<double>(nx);
summary_["shower_" + std::to_string(showerId_)]["ny"] = static_cast<double>(ny);
summary_["shower_" + std::to_string(showerId_)]["nz"] = static_cast<double>(nz);
summary_["shower_" + std::to_string(showerId_)]["px"] =
static_cast<double>(px / 1_GeV);
summary_["shower_" + std::to_string(showerId_)]["py"] =
static_cast<double>(py / 1_GeV);
summary_["shower_" + std::to_string(showerId_)]["pz"] =
static_cast<double>(pz / 1_GeV);
summary_["shower_" + std::to_string(showerId_)]["time"] = primary.getTime() / 1_s;
summary_["shower_" + std::to_string(showerId_)]["slant_depth"] =
dX / (1_g / 1_cm / 1_cm);
uint nSecondaries = 0;
// Loop through secondaries
auto particle = vS.begin();
while (particle != vS.end()) {
// Get the momentum of the secondary, write w.r.t. observation plane
auto const p_2nd = particle.getMomentum();
*(output_.getWriter())
<< showerId_ << static_cast<int>(get_PDG(particle.getPID()))
<< static_cast<float>(p_2nd.dot(obsPlane_.getXAxis()) / 1_GeV)
<< static_cast<float>(p_2nd.dot(obsPlane_.getYAxis()) / 1_GeV)
<< static_cast<float>(p_2nd.dot(obsPlane_.getPlane().getNormal()) / 1_GeV)
<< parquet::EndRow;
CORSIKA_LOG_INFO(" 2ndary: {}, E_kin {}", particle.getPID(),
particle.getKineticEnergy());
++particle;
++nSecondaries;
}
summary_["shower_" + std::to_string(showerId_)]["n_secondaries"] = nSecondaries;
}
template <typename TTracking, typename TOutput>
inline void InteractionWriter<TTracking, TOutput>::startOfLibrary(
boost::filesystem::path const& directory) {
output_.initStreamer((directory / ("interactions.parquet")).string());
// enable compression with the default level
output_.enableCompression();
output_.addField("pdg", parquet::Repetition::REQUIRED, parquet::Type::INT32,
parquet::ConvertedType::INT_32);
output_.addField("px", parquet::Repetition::REQUIRED, parquet::Type::FLOAT,
parquet::ConvertedType::NONE);
output_.addField("py", parquet::Repetition::REQUIRED, parquet::Type::FLOAT,
parquet::ConvertedType::NONE);
output_.addField("pz", parquet::Repetition::REQUIRED, parquet::Type::FLOAT,
parquet::ConvertedType::NONE);
output_.buildStreamer();
showerId_ = 0;
interactionCounter_ = 0;
summary_ = YAML::Node();
}
template <typename TTracking, typename TOutput>
inline void InteractionWriter<TTracking, TOutput>::startOfShower(
unsigned int const showerId) {
showerId_ = showerId;
interactionCounter_ = 0;
}
template <typename TTracking, typename TOutput>
inline void InteractionWriter<TTracking, TOutput>::endOfShower(unsigned int const) {}
template <typename TTracking, typename TOutput>
inline void InteractionWriter<TTracking, TOutput>::endOfLibrary() {
output_.closeStreamer();
}
template <typename TTracking, typename TOutput>
inline YAML::Node InteractionWriter<TTracking, TOutput>::getConfig() const {
YAML::Node node;
node["type"] = "Interactions";
node["units"]["energy"] = "GeV";
node["units"]["length"] = "m";
node["units"]["time"] = "ns";
node["units"]["grammage"] = "g/cm^2";
return node;
}
template <typename TTracking, typename TOutput>
inline YAML::Node InteractionWriter<TTracking, TOutput>::getSummary() const {
return summary_;
}
} // namespace corsika
corsika/detail/modules/writers/LongitudinalProfileWriterParquet.inl 0 → 100644
View file @ 00b7c3f3
/*
* (c) Copyright 2021 CORSIKA Project, corsika-project@lists.kit.edu
*
* This software is distributed under the terms of the 3-clause BSD license.
* See file LICENSE for a full version of the license.
*/
#pragma once
#include <corsika/framework/core/ParticleProperties.hpp>
#include <corsika/framework/core/PhysicalUnits.hpp>
#include <corsika/framework/utility/FindXmax.hpp>
#include <corsika/media/ShowerAxis.hpp>
#include <string>
#include <exception>
namespace corsika {
template <size_t NColumns>
inline LongitudinalProfileWriterParquet<NColumns>::LongitudinalProfileWriterParquet(
std::array<const char*, NColumns> const& columns)
: columns_(columns) {}
template <size_t NColumns>
inline void LongitudinalProfileWriterParquet<NColumns>::startOfLibrary(
boost::filesystem::path const& directory) {
// setup the streamer
output_.initStreamer((directory / "profile.parquet").string());
// enable compression with the default level
output_.enableCompression();
// build the schema
output_.addField("X", parquet::Repetition::REQUIRED, parquet::Type::FLOAT,
parquet::ConvertedType::NONE);
for (auto const& col : columns_) {
output_.addField(col, parquet::Repetition::REQUIRED, parquet::Type::FLOAT,
parquet::ConvertedType::NONE);
}
// and build the streamer
output_.buildStreamer();
}
template <size_t NColumns>
inline void LongitudinalProfileWriterParquet<NColumns>::startOfShower(
unsigned int const) {}
template <size_t NColumns>
inline void LongitudinalProfileWriterParquet<NColumns>::endOfShower(
unsigned int const) {}
template <size_t NColumns>
inline void LongitudinalProfileWriterParquet<NColumns>::endOfLibrary() {
output_.closeStreamer();
}
template <size_t NColumns>
inline void LongitudinalProfileWriterParquet<NColumns>::write(
unsigned int const showerId, GrammageType const grammage,
std::array<double, NColumns> const& data) {
// and write the data into the column
*(output_.getWriter()) << showerId
<< static_cast<float>(grammage / 1_g * square(1_cm));
for (double const weight : data) {
*(output_.getWriter()) << static_cast<float>(weight);
}
*(output_.getWriter()) << parquet::EndRow;
}
} // namespace corsika
corsika/detail/modules/writers/LongitudinalWriter.inl 0 → 100644
View file @ 00b7c3f3
/*
* (c) Copyright 2021 CORSIKA Project, corsika-project@lists.kit.edu
*
* This software is distributed under the terms of the 3-clause BSD license.
* See file LICENSE for a full version of the license.
*/
#pragma once
#include <corsika/framework/core/ParticleProperties.hpp>
#include <corsika/framework/core/PhysicalUnits.hpp>
#include <corsika/framework/utility/FindXmax.hpp>
#include <corsika/media/ShowerAxis.hpp>
#include <exception>
#include <algorithm>
#include <iostream>
namespace corsika {
template <typename TOutput>
inline LongitudinalWriter<TOutput>::LongitudinalWriter(ShowerAxis const& axis,
GrammageType dX)
: LongitudinalWriter<TOutput>{
axis, static_cast<unsigned int>(axis.getMaximumX() / dX) + 1, dX} {}
template <typename TOutput>
inline LongitudinalWriter<TOutput>::LongitudinalWriter(ShowerAxis const& axis,
size_t nbins, GrammageType dX)
: TOutput(number_profile::ProfileIndexNames)
, showerAxis_(axis)
, dX_(dX)
, nBins_(nbins)
, profile_{nbins} {}
template <typename TOutput>
inline void LongitudinalWriter<TOutput>::startOfLibrary(
boost::filesystem::path const& directory) {
TOutput::startOfLibrary(directory);
}
template <typename TOutput>
inline void LongitudinalWriter<TOutput>::startOfShower(unsigned int const showerId) {
profile_.clear();
for (size_t i = 0; i < nBins_; ++i) { profile_.emplace_back(); }
TOutput::startOfShower(showerId);
}
template <typename TOutput>
inline void LongitudinalWriter<TOutput>::endOfLibrary() {
TOutput::endOfLibrary();
}
template <typename TOutput>
inline void LongitudinalWriter<TOutput>::endOfShower(unsigned int const showerId) {
int iRow{0};
for (number_profile::ProfileData const& row : profile_) {
// here: write to underlying writer (e.g. parquet)
TOutput::write(showerId, iRow * dX_, row);
iRow++;
}
TOutput::endOfShower(showerId);
}
template <typename TOutput>
inline void LongitudinalWriter<TOutput>::write(Point const& p0, Point const& p1,
Code const pid, double const weight) {
GrammageType const grammageStart = showerAxis_.getProjectedX(p0);
GrammageType const grammageEnd = showerAxis_.getProjectedX(p1);
// Avoid over counting in first bin when backscattered particle goes beyond the
// injection point.
if (grammageStart == grammageEnd) { return; }
// Note: particle may go also "upward", thus, grammageEnd<grammageStart
size_t const binStart = std::ceil(grammageStart / dX_);
size_t const binEnd = std::floor(grammageEnd / dX_);
CORSIKA_LOGGER_TRACE(TOutput::getLogger(),
"grammageStart={} End={} binStart={}, end={}",
grammageStart / 1_g * square(1_cm),
grammageEnd / 1_g * square(1_cm), binStart, binEnd);
for (size_t bin = binStart; bin <= std::min(binEnd, profile_.size() - 1); ++bin) {
if (pid == Code::Photon) {
profile_.at(bin)[static_cast<int>(number_profile::ProfileIndex::Photon)] +=
weight;
} else if (pid == Code::Positron) {
profile_.at(bin)[static_cast<int>(number_profile::ProfileIndex::Positron)] +=
weight;
} else if (pid == Code::Electron) {
profile_.at(bin)[static_cast<int>(number_profile::ProfileIndex::Electron)] +=
weight;
} else if (pid == Code::MuPlus) {
profile_.at(bin)[static_cast<int>(number_profile::ProfileIndex::MuPlus)] +=
weight;
} else if (pid == Code::MuMinus) {
profile_.at(bin)[static_cast<int>(number_profile::ProfileIndex::MuMinus)] +=
weight;
} else if (is_hadron(pid)) {
profile_.at(bin)[static_cast<int>(number_profile::ProfileIndex::Hadron)] +=
weight;
}
if (is_charged(pid)) {
profile_[bin][static_cast<int>(number_profile::ProfileIndex::Charged)] += weight;
}
}
}
template <typename TOutput>
inline void LongitudinalWriter<TOutput>::write(GrammageType const Xstart,
GrammageType const Xend, Code const pid,
double const weight) {
double const bstart = Xstart / dX_;
double const bend = Xend / dX_;
if (abs(bstart - floor(bstart + 0.5)) > 1e-2 ||
abs(bend - floor(bend + 0.5)) > 1e-2 || abs(bend - bstart - 1) > 1e-2) {
CORSIKA_LOGGER_ERROR(TOutput::getLogger(),
"CONEX and Corsika8 dX grammage binning are not the same! "
"Xstart={} Xend={} dX={}",
Xstart / 1_g * square(1_cm), Xend / 1_g * square(1_cm),
dX_ / 1_g * square(1_cm));
throw std::runtime_error(
"CONEX and Corsika8 dX grammage binning are not the same!");
}
size_t const bin = size_t((bend + bstart) / 2);
if (bin >= profile_.size()) {
CORSIKA_LOGGER_WARN(TOutput::getLogger(),
"Grammage bin {} outside of profile {}. skipping.", bin,
profile_.size());
return;
}
if (pid == Code::Photon) {
profile_.at(bin)[static_cast<int>(number_profile::ProfileIndex::Photon)] += weight;
} else if (pid == Code::Positron) {
profile_.at(bin)[static_cast<int>(number_profile::ProfileIndex::Positron)] +=
weight;
} else if (pid == Code::Electron) {
profile_.at(bin)[static_cast<int>(number_profile::ProfileIndex::Electron)] +=
weight;
} else if (pid == Code::MuPlus) {
profile_.at(bin)[static_cast<int>(number_profile::ProfileIndex::MuPlus)] += weight;
} else if (pid == Code::MuMinus) {
profile_.at(bin)[static_cast<int>(number_profile::ProfileIndex::MuMinus)] += weight;
} else if (is_hadron(pid)) {
profile_.at(bin)[static_cast<int>(number_profile::ProfileIndex::Hadron)] += weight;
}
if (is_charged(pid)) {
profile_[bin][static_cast<int>(number_profile::ProfileIndex::Charged)] += weight;
}
}
template <typename TOutput>
inline YAML::Node LongitudinalWriter<TOutput>::getConfig() const {
YAML::Node node;
node["type"] = "LongitudinalProfile";
node["units"]["grammage"] = "g/cm^2";
node["bin-size"] = dX_ / (1_g / square(1_cm));
node["nbins"] = nBins_;
return node;
}
template <typename TOutput>
inline YAML::Node LongitudinalWriter<TOutput>::getSummary() const {
YAML::Node summary;
return summary;
}
} // namespace corsika
\ No newline at end of file

data privacy - impressum