IAP GITLAB

Skip to content
Snippets Groups Projects
Commit d7ac9f8a authored by ralfulrich's avatar ralfulrich Committed by Ralf Ulrich
Browse files

qgsjetII new interface

parent 6b237f7c
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
Showing
with 514 additions and 90 deletions
corsika/detail/modules/qgsjetII/Interaction.inl → corsika/detail/modules/qgsjetII/InteractionModel.inl
+ 308
− 0
View file @ d7ac9f8a
......@@ -6,10 +6,8 @@
* the license.
*/
#include <corsika/modules/qgsjetII/Interaction.hpp>
#include <corsika/modules/qgsjetII/InteractionModel.hpp>
#include <corsika/media/Environment.hpp>
#include <corsika/media/NuclearComposition.hpp>
#include <corsika/framework/geometry/QuantityVector.hpp>
#include <corsika/framework/geometry/FourVector.hpp>
#include <corsika/modules/qgsjetII/ParticleConversion.hpp>
......@@ -18,14 +16,13 @@
#include <corsika/framework/utility/COMBoost.hpp>
#include <sstream>
#include <string>
#include <tuple>
#include <qgsjet-II-04.hpp>
namespace corsika::qgsjetII {
inline Interaction::Interaction(boost::filesystem::path dataPath) {
inline InteractionModel::InteractionModel(boost::filesystem::path const dataPath) {
CORSIKA_LOG_DEBUG("Reading QGSJetII data tables from {}", dataPath);
// initialize QgsjetII
......@@ -38,246 +35,121 @@ namespace corsika::qgsjetII {
}
}
inline Interaction::~Interaction() {
CORSIKA_LOG_DEBUG("QgsjetII::Interaction n= {}", count_);
inline InteractionModel::~InteractionModel() {
CORSIKA_LOG_DEBUG("QgsjetII::InteractionModel n= {}", count_);
}
inline CrossSectionType Interaction::getCrossSection(const Code beamId,
const Code targetId,
const HEPEnergyType Elab,
const unsigned int Abeam,
const unsigned int targetA) const {
double sigProd = std::numeric_limits<double>::infinity();
if (corsika::qgsjetII::canInteract(beamId)) {
int const iBeam = static_cast<QgsjetIIXSClassIntType>(
corsika::qgsjetII::getQgsjetIIXSCode(beamId));
int iTarget = 1;
if (is_nucleus(targetId)) {
iTarget = targetA;
if (iTarget > int(maxMassNumber_) || iTarget <= 0) {
std::ostringstream txt;
txt << "QgsjetII target outside range. Atarget=" << iTarget;
throw std::runtime_error(txt.str().c_str());
}
}
int iProjectile = 1;
if (is_nucleus(beamId)) {
iProjectile = Abeam;
if (iProjectile > int(maxMassNumber_) || iProjectile <= 0) {
std::ostringstream txt;
txt << "QgsjetII projectile outside range. Aprojectile=" << iProjectile;
throw std::runtime_error(txt.str().c_str());
}
inline void InteractionModel::isValid(Code const projectileId,
Code const targetId) const {
if (is_nucleus(targetId)) {
size_t iTarget = get_nucleus_A(targetId);
if (iTarget > int(maxMassNumber_) || iTarget <= 0) {
std::ostringstream txt;
txt << "QgsjetII target outside range. Atarget=" << iTarget;
throw std::runtime_error(txt.str().c_str());
}
CORSIKA_LOG_DEBUG(
"QgsjetII::getCrossSection Elab= {} GeV iBeam= {}"
" iProjectile= {} iTarget= {}",
Elab / 1_GeV, iBeam, iProjectile, iTarget);
sigProd = qgsect_(Elab / 1_GeV, iBeam, iProjectile, iTarget);
CORSIKA_LOG_DEBUG("QgsjetII::getCrossSection sigProd= {} mb", sigProd);
} else if (targetId != Proton::code) {
std::ostringstream txt;
txt << "QgsjetII not valid target=" << targetId;
throw std::runtime_error(txt.str().c_str());
}
return sigProd * 1_mb;
if (is_nucleus(projectileId)) {
size_t iProjectile = get_nucleus_A(projectileId);
if (iProjectile > int(maxMassNumber_) || iProjectile <= 0) {
std::ostringstream txt;
txt << "QgsjetII projectile outside range. Aprojectile=" << iProjectile;
throw std::runtime_error(txt.str().c_str());
}
} else if (!is_hadron(projectileId)) {
std::ostringstream txt;
txt << "QgsjetII projectile can only be hadrons. Is=" << projectileId;
throw std::runtime_error(txt.str().c_str());
}
}
template <typename TParticle>
inline GrammageType Interaction::getInteractionLength(const TParticle& particle) const {
inline CrossSectionType InteractionModel::getCrossSection(
Code const projectileId, Code const targetId, FourMomentum const& projectileP4,
FourMomentum const& targetP4) const {
// coordinate system, get global frame of reference
CoordinateSystemPtr const& rootCS = get_root_CoordinateSystem();
double sigProd = std::numeric_limits<double>::infinity();
const Code corsikaBeamId = particle.getPID();
if (!corsika::qgsjetII::canInteract(projectileId)) { return sigProd * 1_mb; }
// beam particles for qgsjetII : 1, 2, 3 for p, pi, k
// read from cross section code table
const bool kInteraction = corsika::qgsjetII::canInteract(corsikaBeamId);
isValid(projectileId, targetId); // throws
// FOR NOW: assume target is at rest
MomentumVector pTarget(rootCS, {0_GeV, 0_GeV, 0_GeV});
// 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));
// total momentum and energy
HEPEnergyType const Elab = particle.getEnergy();
int const iBeam = static_cast<QgsjetIIXSClassIntType>(
corsika::qgsjetII::getQgsjetIIXSCode(projectileId));
int iTarget = 1;
if (is_nucleus(targetId)) { iTarget = get_nucleus_A(targetId); }
int iProjectile = 1;
if (is_nucleus(projectileId)) { iProjectile = get_nucleus_A(projectileId); }
CORSIKA_LOG_DEBUG(
"Interaction: LambdaInt: \n"
" input energy: {} GeV"
" beam can interact: {}"
" beam pid: {}",
particle.getEnergy() / 1_GeV, kInteraction, corsikaBeamId);
if (kInteraction) {
int Abeam = 0;
if (is_nucleus(corsikaBeamId)) Abeam = get_nucleus_A(corsikaBeamId);
// 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 = particle.getNode();
const auto& mediumComposition =
currentNode->getModelProperties().getNuclearComposition();
CrossSectionType weightedProdCrossSection =
mediumComposition.getWeightedSum([=](Code targetID) -> CrossSectionType {
int targetA = 0;
if (is_nucleus(targetID)) targetA = get_nucleus_A(targetID);
return getCrossSection(corsikaBeamId, targetID, Elab, Abeam, targetA);
});
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;
}
"QgsjetII::getCrossSection Elab= {} GeV iBeam= {}"
" iProjectile= {} iTarget= {}",
Elab / 1_GeV, iBeam, iProjectile, iTarget);
sigProd = qgsect_(Elab / 1_GeV, iBeam, iProjectile, iTarget);
CORSIKA_LOG_DEBUG("QgsjetII::getCrossSection sigProd= {} mb", sigProd);
return std::numeric_limits<double>::infinity() * 1_g / (1_cm * 1_cm);
return sigProd * 1_mb;
}
/**
In this function QGSJETII is called to produce one event. The
event is copied (and boosted) into the shower lab frame.
*/
template <typename TSecondaries>
inline void InteractionModel::doInteraction(TSecondaries& view, Code const projectileId,
Code const targetId,
FourMomentum const& projectileP4,
FourMomentum const& targetP4) {
template <typename TView>
inline void Interaction::doInteraction(TView& view) {
auto const projectile = view.getProjectile();
auto const corsikaBeamId = projectile.getPID();
CORSIKA_LOG_DEBUG(
"ProcessQgsjetII: "
"doInteraction: {} interaction possible? {}",
corsikaBeamId, corsika::qgsjetII::canInteract(corsikaBeamId));
projectileId, corsika::qgsjetII::canInteract(projectileId));
if (!corsika::qgsjetII::canInteract(corsikaBeamId)) return;
if (!corsika::qgsjetII::canInteract(projectileId)) return;
CoordinateSystemPtr const& rootCS = get_root_CoordinateSystem();
// position and time of interaction, not used in QgsjetII
Point const pOrig = projectile.getPosition();
TimeType const tOrig = projectile.getTime();
isValid(projectileId, targetId); // throws
// define target
// for QgsjetII is always a single nucleon
// FOR NOW: target is always at rest
auto const targetEnergyLab = 0_GeV + constants::nucleonMass;
auto const targetMomentumLab = MomentumVector(rootCS, 0_GeV, 0_GeV, 0_GeV);
FourVector const PtargLab(targetEnergyLab, targetMomentumLab);
CoordinateSystemPtr const& rootCS = get_root_CoordinateSystem();
// define projectile
HEPEnergyType const projectileEnergyLab = projectile.getEnergy();
auto const projectileMomentumLab = projectile.getMomentum();
// 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));
int beamA = 0;
if (is_nucleus(corsikaBeamId)) beamA = get_nucleus_A(corsikaBeamId);
HEPEnergyType const projectileEnergyLabPerNucleon = projectileEnergyLab / beamA;
if (is_nucleus(projectileId)) { beamA = get_nucleus_A(projectileId); }
CORSIKA_LOG_DEBUG(
"ebeam lab: {} GeV "
"pbeam lab: {} GeV ",
projectileEnergyLab / 1_GeV, projectileMomentumLab.getComponents() / 1_GeV);
CORSIKA_LOG_DEBUG(
"etarget lab: {} GeV "
"ptarget lab: {} GeV ",
targetEnergyLab / 1_GeV, targetMomentumLab.getComponents() / 1_GeV);
CORSIKA_LOG_DEBUG("position of interaction: {}", pOrig.getCoordinates());
CORSIKA_LOG_DEBUG("time: {} ", tOrig);
// 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
*/
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];
int targetA = 0;
if (is_nucleus(targetId)) targetA = get_nucleus_A(targetId);
const auto sigProd =
getCrossSection(corsikaBeamId, targetId, projectileEnergyLab, beamA, targetA);
cross_section_of_components[i] = sigProd;
}
CORSIKA_LOG_DEBUG("ebeam lab: {} GeV ", Elab / 1_GeV);
const auto targetCode =
mediumComposition.sampleTarget(cross_section_of_components, rng_);
int targetMassNumber = 1; // proton
if (is_nucleus(targetCode)) { // nucleus
targetMassNumber = get_nucleus_A(targetCode);
if (targetMassNumber > int(maxMassNumber_))
throw std::runtime_error(
"QgsjetII target mass outside range."); // LCOV_EXCL_LINE there is no
// allowed path here
} else {
if (targetCode != Proton::code) // LCOV_EXCL_LINE there is no allowed path here
throw std::runtime_error(
"QgsjetII Taget not possible."); // LCOV_EXCL_LINE there is no allowed path
// here
int targetMassNumber = 1; // proton
if (is_nucleus(targetId)) { // nucleus
targetMassNumber = get_nucleus_A(targetId);
}
CORSIKA_LOG_DEBUG("target: {}, qgsjetII code/A: {}", targetCode, targetMassNumber);
CORSIKA_LOG_DEBUG("target: {}, qgsjetII code/A: {}", targetId, targetMassNumber);
// select QGSJetII internal projectile type
int projectileMassNumber = 1; // "1" means "hadron"
QgsjetIIHadronType qgsjet_hadron_type =
qgsjetII::getQgsjetIIHadronType(corsikaBeamId);
QgsjetIIHadronType qgsjet_hadron_type = qgsjetII::getQgsjetIIHadronType(projectileId);
if (qgsjet_hadron_type == QgsjetIIHadronType::NucleusType) {
projectileMassNumber = get_nucleus_A(corsikaBeamId);
if (projectileMassNumber > int(maxMassNumber_))
throw std::runtime_error(
"QgsjetII projectile mass outside range."); // LCOV_EXCL_LINE there is no
// allowed path here
projectileMassNumber = get_nucleus_A(projectileId);
std::array<QgsjetIIHadronType, 2> constexpr nucleons = {
QgsjetIIHadronType::ProtonType, QgsjetIIHadronType::NeutronType};
std::uniform_int_distribution select(0, 1);
qgsjet_hadron_type = nucleons[select(rng_)];
} else {
} else if (qgsjet_hadron_type == QgsjetIIHadronType::NeutralLightMesonType) {
// from conex: replace pi0 or rho0 with pi+/pi- in alternating sequence
if (qgsjet_hadron_type == QgsjetIIHadronType::NeutralLightMesonType) {
qgsjet_hadron_type = alternate_;
alternate_ = (alternate_ == QgsjetIIHadronType::PiPlusType
? QgsjetIIHadronType::PiMinusType
: QgsjetIIHadronType::PiPlusType);
}
}
// beam id for qgsjetII
int kBeam = 2; // default: proton Shouldn't we randomize neutron/proton for nuclei?
if (!is_nucleus(corsikaBeamId)) {
kBeam = corsika::qgsjetII::convertToQgsjetIIRaw(corsikaBeamId);
// from conex
if (kBeam == 0) { // replace pi0 or rho0 with pi+/pi-
static int select = 1;
kBeam = select;
select *= -1;
}
// replace lambda by neutron
if (kBeam == 6)
kBeam = 3;
else if (kBeam == -6)
kBeam = -3;
// else if (abs(kBeam)>6) -> throw
qgsjet_hadron_type = alternate_;
alternate_ =
(alternate_ == QgsjetIIHadronType::PiPlusType ? QgsjetIIHadronType::PiMinusType
: QgsjetIIHadronType::PiPlusType);
}
count_++;
......@@ -285,8 +157,7 @@ namespace corsika::qgsjetII {
CORSIKA_LOG_DEBUG(
"qgsjet_hadron_type_int={} projectileMassNumber={} targetMassNumber={}",
qgsjet_hadron_type_int, projectileMassNumber, targetMassNumber);
qgini_(projectileEnergyLab / 1_GeV, qgsjet_hadron_type_int, projectileMassNumber,
targetMassNumber);
qgini_(Elab / 1_GeV, qgsjet_hadron_type_int, projectileMassNumber, targetMassNumber);
qgconf_();
// bookkeeping
......@@ -296,9 +167,29 @@ namespace corsika::qgsjetII {
// to read the secondaries
// define rotation to and from CoM frame
// CoM frame definition in QgsjetII projectile: +z
auto const& originalCS = projectileMomentumLab.getCoordinateSystem();
CoordinateSystemPtr const zAxisFrame =
make_rotationToZ(originalCS, projectileMomentumLab);
// QGSJetII, both, in input and output only considers the lab frame with a target at
// rest.
// system of initial-state
COMBoost boost(projectileP4, targetP4);
auto const& originalCS = boost.getOriginalCS();
auto const& csPrime =
boost.getRotatedCS(); // z is along the CM motion (projectile, in Cascade)
MomentumVector const projectileMomentum = projectileP4.getSpaceLikeComponents();
HEPMomentumType const pLabMag =
sqrt((Elab - get_mass(projectileId)) * (Elab + get_mass(projectileId)));
MomentumVector pLab(csPrime, {0_eV, 0_eV, pLabMag});
// internal QGSJetII system
COMBoost boostInternal({Elab, pLab}, get_mass(targetId));
// position and time of interaction, not used in QgsjetII
auto const projectile = view.getProjectile();
Point const pOrig = projectile.getPosition();
TimeType const tOrig = projectile.getTime();
// fragments
QGSJetIIFragmentsStack qfs;
......@@ -310,19 +201,26 @@ namespace corsika::qgsjetII {
if (select(rng_) > 0.5) { idFragm = Code::Neutron; }
const HEPMassType nucleonMass = get_mass(idFragm);
// no pT, frgments just go forward
auto momentum =
Vector(zAxisFrame, corsika::QuantityVector<hepmomentum_d>{
0.0_GeV, 0.0_GeV,
sqrt((projectileEnergyLabPerNucleon + nucleonMass) *
(projectileEnergyLabPerNucleon - nucleonMass))});
// no pT, fragments just go forward
HEPEnergyType const projectileEnergyLabPerNucleon = Elab / beamA;
MomentumVector momentum{csPrime,
{0.0_GeV, 0.0_GeV,
sqrt((projectileEnergyLabPerNucleon + nucleonMass) *
(projectileEnergyLabPerNucleon - nucleonMass))}};
// this is not "CoM" here, but rather the system defined by projectile+target,
// which in Cascade-mode is already lab
auto const P4com =
boostInternal.toCoM(FourVector{projectileEnergyLabPerNucleon, momentum});
auto const P4output = boost.fromCoM(P4com);
auto p3output = P4output.getSpaceLikeComponents();
p3output.rebase(originalCS); // transform back into standard lab frame
momentum.rebase(originalCS); // transform back into standard lab frame
CORSIKA_LOG_DEBUG(
"secondary fragment> id= {}"
" p={}",
idFragm, momentum.getComponents());
auto pnew = view.addSecondary(std::make_tuple(idFragm, momentum, pOrig, tOrig));
idFragm, p3output.getComponents());
auto pnew = view.addSecondary(std::make_tuple(idFragm, p3output, pOrig, tOrig));
Plab_final += pnew.getMomentum();
Elab_final += pnew.getEnergy();
......@@ -347,22 +245,30 @@ namespace corsika::qgsjetII {
HEPMassType const nucleusMass = Proton::mass * Z + Neutron::mass * (A - Z);
// no pT, frgments just go forward
auto momentum = Vector(
zAxisFrame, QuantityVector<hepmomentum_d>{
0.0_GeV, 0.0_GeV,
sqrt((projectileEnergyLabPerNucleon * A + nucleusMass) *
(projectileEnergyLabPerNucleon * A - nucleusMass))});
HEPEnergyType const projectileEnergyLabPerNucleon = Elab / beamA;
MomentumVector momentum{
csPrime,
{0.0_GeV, 0.0_GeV,
sqrt((projectileEnergyLabPerNucleon * A + nucleusMass) *
(projectileEnergyLabPerNucleon * A - nucleusMass))}};
// this is not "CoM" here, but rather the system defined by projectile+target,
// which in Cascade-mode is already lab
auto const P4com =
boostInternal.toCoM(FourVector{projectileEnergyLabPerNucleon * A, momentum});
auto const P4output = boost.fromCoM(P4com);
auto p3output = P4output.getSpaceLikeComponents();
p3output.rebase(originalCS); // transform back into standard lab frame
momentum.rebase(originalCS); // transform back into standard lab frame
CORSIKA_LOG_DEBUG(
"secondary fragment> id={}"
" p={}"
" A={}"
" Z={}",
get_nucleus_code(A, Z), momentum.getComponents(), A, Z);
get_nucleus_code(A, Z), p3output.getComponents(), A, Z);
auto pnew = view.addSecondary(
std::make_tuple(get_nucleus_code(A, Z), momentum, pOrig, tOrig));
std::make_tuple(get_nucleus_code(A, Z), p3output, pOrig, tOrig));
Plab_final += pnew.getMomentum();
Elab_final += pnew.getEnergy();
}
......@@ -372,14 +278,19 @@ namespace corsika::qgsjetII {
QGSJetIIStack qs;
for (auto& psec : qs) {
auto momentum = psec.getMomentum(zAxisFrame);
auto momentum = psec.getMomentum(csPrime);
// this is not "CoM" here, but rather the system defined by projectile+target,
// which in Cascade-mode is already lab
auto const P4com = boostInternal.toCoM(FourVector{psec.getEnergy(), momentum});
auto const P4output = boost.fromCoM(P4com);
auto p3output = P4output.getSpaceLikeComponents();
p3output.rebase(originalCS); // transform back into standard lab frame
momentum.rebase(originalCS); // transform back into standard lab frame
CORSIKA_LOG_DEBUG("secondary> id= {}, p= {}",
corsika::qgsjetII::convertFromQgsjetII(psec.getPID()),
momentum.getComponents());
p3output.getComponents());
auto pnew = view.addSecondary(std::make_tuple(
corsika::qgsjetII::convertFromQgsjetII(psec.getPID()), momentum, pOrig, tOrig));
corsika::qgsjetII::convertFromQgsjetII(psec.getPID()), p3output, pOrig, tOrig));
Plab_final += pnew.getMomentum();
Elab_final += pnew.getEnergy();
}
......
corsika/detail/modules/sibyll/InteractionModel.inl
+ 18
− 17
View file @ d7ac9f8a
......@@ -156,8 +156,8 @@ namespace corsika::sibyll {
auto const tmp = psib.getMomentum().getComponents();
auto const pCoM = MomentumVector(csPrime, tmp);
HEPEnergyType const eCoM = psib.getEnergy();
auto const Plab = boost.fromCoM(FourVector{eCoM, pCoM});
auto const p3lab = Plab.getSpaceLikeComponents();
auto const P4lab = boost.fromCoM(FourVector{eCoM, pCoM});
auto const p3lab = P4lab.getSpaceLikeComponents();
// add to corsika stack
auto pnew = secondaries.addSecondary(std::make_tuple(
......@@ -167,19 +167,20 @@ namespace corsika::sibyll {
Elab_final += pnew.getEnergy();
Ecm_final += psib.getEnergy();
}
HEPEnergyType const Elab_initial =
static_pow<2>(sqrtSnn) / (2 * constants::nucleonMass);
CORSIKA_LOG_DEBUG(
"conservation (all GeV): "
"sqrtSnn={}, sqrtSnn_final={}, "
"Elab_initial={}, Elab_final={}, "
"diff(%)={}, "
"E in nucleons={}, "
"Plab_final={} ",
sqrtSnn / 1_GeV, Ecm_final * 2. / (get_nwounded() + 1) / 1_GeV, Elab_initial,
Elab_final / 1_GeV, (Elab_final - Elab_initial) / Elab_initial * 100,
constants::nucleonMass * get_nwounded() / 1_GeV,
(Plab_final / 1_GeV).getComponents());
} // namespace corsika::sibyll
{ // just output
HEPEnergyType const Elab_initial =
static_pow<2>(sqrtSnn) / (2 * constants::nucleonMass);
CORSIKA_LOG_DEBUG(
"conservation (all GeV): "
"sqrtSnn={}, sqrtSnn_final={}, "
"Elab_initial={}, Elab_final={}, "
"diff(%)={}, "
"E in nucleons={}, "
"Plab_final={} ",
sqrtSnn / 1_GeV, Ecm_final * 2. / (get_nwounded() + 1) / 1_GeV, Elab_initial,
Elab_final / 1_GeV, (Elab_final - Elab_initial) / Elab_initial * 100,
constants::nucleonMass * get_nwounded() / 1_GeV,
(Plab_final / 1_GeV).getComponents());
}
}
} // namespace corsika::sibyll
corsika/framework/utility/COMBoost.hpp
+ 21
− 11
View file @ d7ac9f8a
......@@ -25,11 +25,11 @@ namespace corsika {
*/
/**
* @class COMBoost
* @class COMBoost
* @ingroup Utilities
*
* This utility class handles Lorentz boost between different
* referenence frames, using FourVector.
* This utility class handles Lorentz boost (in one spatial direction)
* between different referenence frames, using FourVector.
*
* The class is initialized with projectile and optionally target
* energy/momentum data. During initialization, a rotation matrix is
......@@ -47,16 +47,26 @@ namespace corsika {
class COMBoost {
public:
//! construct a COMBoost given four-vector of projectile and mass of target (target at
//! rest)
/**
* Construct a COMBoost given four-vector of projectile and mass of target (target at
* rest).
*
* The FourMomentum and mass define the lab system.
*/
COMBoost(FourMomentum const& P4projectile, HEPEnergyType const massTarget);
//! construct a COMBoost given two four-vectors of projectile target
COMBoost(FourMomentum const& P4projectile, FourMomentum const& P4target);
//! construct a COMBoost to boost into the rest frame given a 3-momentum and mass
/**
* Construct a COMBoost to boost into the rest frame given a 3-momentum and mass.
*/
COMBoost(MomentumVector const& momentum, HEPEnergyType const mass);
/**
* Construct a COMBoost given two four-vectors of projectile target.
*
* The tow FourMomentum can define an arbitrary system.
*/
COMBoost(FourMomentum const& P4projectile, FourMomentum const& P4target);
//! transforms a 4-momentum from lab frame to the center-of-mass frame
template <typename FourVector>
FourVector toCoM(FourVector const& p4) const;
......@@ -65,10 +75,10 @@ namespace corsika {
template <typename FourVector>
FourVector fromCoM(FourVector const& p4) const;
//! returns the rotated coordinate system
//! returns the rotated coordinate system: +z is projectile direction
CoordinateSystemPtr getRotatedCS() const;
//! returns the original coordinate system
//! returns the original coordinate system of the projectile (lab)
CoordinateSystemPtr getOriginalCS() const;
protected:
......
corsika/modules/LongitudinalProfile.hpp
+ 1
− 2
View file @ d7ac9f8a
......@@ -31,8 +31,7 @@ namespace corsika {
* simulation into a projected grammage range and counts for
* different particle species when they cross dX (default: 10g/cm2)
* boundaries.
*
**/
*/
class LongitudinalProfile : public ContinuousProcess<LongitudinalProfile> {
......
corsika/modules/QGSJetII.hpp
+ 20
− 1
View file @ d7ac9f8a
......@@ -8,4 +8,23 @@
#pragma once
#include <corsika/modules/qgsjetII/Interaction.hpp>
#include <corsika/modules/qgsjetII/InteractionModel.hpp>
#include <corsika/framework/process/InteractionProcess.hpp>
/**
* @file QGSJetII.hpp
*
* Includes all the parts of the QGSJetII model. Defines the InteractionProcess<TModel>
* classes needed for the ProcessSequence.
*/
namespace corsika::qgsjetII {
/**
* @brief qgsjetII::Interaction is the process for ProcessSequence.
*
* The qgsjetII::InteractionModel is wrapped as an InteractionProcess here in order
* to provide all the functions for ProcessSequence.
*/
class Interaction : public InteractionModel, public InteractionProcess<Interaction> {};
} // namespace corsika::qgsjetII
corsika/modules/Sibyll.hpp
+ 2
− 2
View file @ d7ac9f8a
......@@ -18,7 +18,7 @@
/**
* @file Sibyll.hpp
*
* Includes all the parts of the Sibyll model. Defines the InteractoinProcess<TModel>
* Includes all the parts of the Sibyll model. Defines the InteractionProcess<TModel>
* classes needed for the ProcessSequence.
*/
......@@ -26,7 +26,7 @@ namespace corsika::sibyll {
/**
* @brief sibyll::Interaction is the process for ProcessSequence.
*
* The sibyll::Model is wrapped as an InteractionProcess here in order
* The sibyll::InteractionModel is wrapped as an InteractionProcess here in order
* to provide all the functions for ProcessSequence.
*/
class Interaction : public InteractionModel, public InteractionProcess<Interaction> {};
......
corsika/modules/qgsjetII/Interaction.hpp → corsika/modules/qgsjetII/InteractionModel.hpp
+ 77
− 0
View file @ d7ac9f8a
......@@ -8,58 +8,70 @@
#pragma once
#include <corsika/modules/qgsjetII/ParticleConversion.hpp>
#include <qgsjet-II-04.hpp>
#include <corsika/framework/core/ParticleProperties.hpp>
#include <corsika/framework/core/PhysicalUnits.hpp>
#include <corsika/framework/random/RNGManager.hpp>
#include <corsika/framework/process/InteractionProcess.hpp>
#include <corsika/modules/qgsjetII/ParticleConversion.hpp>
#include <corsika/framework/utility/COMBoost.hpp>
#include <corsika/framework/utility/CorsikaData.hpp>
#include <boost/filesystem/path.hpp>
#include <qgsjet-II-04.hpp>
#include <string>
namespace corsika::qgsjetII {
class Interaction : public corsika::InteractionProcess<Interaction> {
class InteractionModel {
public:
Interaction(boost::filesystem::path dataPath = corsika_data("QGSJetII"));
~Interaction();
InteractionModel(boost::filesystem::path dataPath = corsika_data("QGSJetII"));
~InteractionModel();
bool wasInitialized() { return initialized_; }
unsigned int getMaxTargetMassNumber() const { return maxMassNumber_; }
bool isValidTarget(corsika::Code TargetId) const {
return is_nucleus(TargetId) && (get_nucleus_A(TargetId) < maxMassNumber_);
}
CrossSectionType getCrossSection(const Code, const Code, const HEPEnergyType,
const unsigned int Abeam = 0,
const unsigned int Atarget = 0) const;
template <typename TParticle>
GrammageType getInteractionLength(TParticle const&) const;
/**
* Throws exception if invalid system is passed.
*
* @param beamId
* @param targetId
*/
void isValid(Code const beamId, Code const targetId) const;
/**
In this function QGSJETII is called to produce one event. The
event is copied (and boosted) into the shower lab frame.
* Return the QGSJETII inelastic/production cross section.
*
* This cross section must correspond to the process described in doInteraction.
* Allowed targets are: nuclei or single nucleons (p,n,hydrogen).
*
* @param projectile is the Code of the projectile
* @param target is the Code of the target
* @param sqrtSnn is the center-of-mass energy (per nucleon pair)
* @param Aprojectile is the mass number of the projectils, if it is a nucleus
* @param Atarget is the mass number of the target, if it is a nucleus
*
* @return inelastic cross section.
*/
CrossSectionType getCrossSection(Code const projectile, Code const target,
FourMomentum const& projectileP4,
FourMomentum const& targetP4) const;
template <typename TSecondaryView>
void doInteraction(TSecondaryView&);
/**
* In this function QGSJETII is called to produce one event.
*
* The event is copied (and boosted) into the shower lab frame.
*/
template <typename TSecondaries>
void doInteraction(TSecondaries&, Code const projectile, Code const target,
FourMomentum const& projectileP4, FourMomentum const& targetP4);
private:
int count_ = 0;
bool initialized_ = false;
QgsjetIIHadronType alternate_ =
QgsjetIIHadronType::PiPlusType; // for pi0, rho0 projectiles
corsika::default_prng_type& rng_ =
corsika::RNGManager<>::getInstance().getRandomStream("qgsjet");
static unsigned int constexpr maxMassNumber_ = 208;
static size_t constexpr maxMassNumber_ = 208;
};
} // namespace corsika::qgsjetII
#include <corsika/detail/modules/qgsjetII/Interaction.inl>
#include <corsika/detail/modules/qgsjetII/InteractionModel.inl>
corsika/modules/qgsjetII/ParticleConversion.hpp
+ 4
− 4
View file @ d7ac9f8a
......@@ -15,13 +15,13 @@
namespace corsika::qgsjetII {
/**
These are the possible secondaries produced by QGSJetII
* These are the possible secondaries produced by QGSJetII.
*/
enum class QgsjetIICode : int8_t;
using QgsjetIICodeIntType = std::underlying_type<QgsjetIICode>::type;
/**
These are the possible projectile for which QGSJetII knwos cross section
* These are the possible projectile for which QGSJetII knwos cross section.
*/
enum class QgsjetIIXSClass : int8_t {
CannotInteract = 0,
......@@ -32,7 +32,7 @@ namespace corsika::qgsjetII {
using QgsjetIIXSClassIntType = std::underlying_type<QgsjetIIXSClass>::type;
/**
These are the only possible projectile types in QGSJetII
* These are the only possible projectile types in QGSJetII.
*/
enum class QgsjetIIHadronType : int8_t {
UndefinedType = 0,
......@@ -58,7 +58,7 @@ namespace corsika::qgsjetII {
namespace corsika::qgsjetII {
QgsjetIICode constexpr convertToQgsjetII(Code pCode) {
QgsjetIICode constexpr convertToQgsjetII(Code const pCode) {
return corsika2qgsjetII[static_cast<CodeIntType>(pCode)];
}
......
corsika/modules/sibyll/InteractionModel.hpp
+ 2
− 2
View file @ d7ac9f8a
......@@ -48,7 +48,7 @@ namespace corsika::sibyll {
HEPEnergyType const sqrtSnn) const;
/**
* @brief returns inelastic AND elastic cross sections.
* Returns inelastic AND elastic cross sections.
*
* These cross sections must correspond to the process described in doInteraction
* AND elastic scattering (sigma_tot = sigma_inel + sigma_el). Allowed targets are:
......@@ -68,7 +68,7 @@ namespace corsika::sibyll {
FourMomentum const& targetP4) const;
/**
* @brief returns inelastic (production) cross section.
* Returns inelastic (production) cross section.
*
* This cross section must correspond to the process described in doInteraction.
* Allowed targets are: nuclei or single nucleons (p,n,hydrogen).
......
tests/modules/CMakeLists.txt
+ 1
− 1
View file @ d7ac9f8a
......@@ -4,7 +4,7 @@ set (test_modules_sources
testTracking.cpp
## testExecTime.cpp --> needs to be fixed, see #326
testObservationPlane.cpp
#testQGSJetII.cpp
testQGSJetII.cpp
#testPythia8.cpp
#testUrQMD.cpp
#testCONEX.cpp
......
tests/modules/testQGSJetII.cpp
+ 57
− 46
View file @ d7ac9f8a
......@@ -6,8 +6,7 @@
* the license.
*/
#include <corsika/modules/qgsjetII/Interaction.hpp>
#include <corsika/modules/qgsjetII/ParticleConversion.hpp>
#include <corsika/modules/QGSJetII.hpp>
#include <corsika/framework/core/ParticleProperties.hpp>
#include <corsika/framework/core/PhysicalUnits.hpp>
......@@ -49,6 +48,7 @@ auto sumMomentum(TStackView const& view, CoordinateSystemPtr const& vCS) {
TEST_CASE("QgsjetII", "[processes]") {
logging::set_level(logging::level::info);
RNGManager<>::getInstance().registerRandomStream("qgsjet");
SECTION("Corsika -> QgsjetII") {
CHECK(corsika::qgsjetII::convertToQgsjetII(PiMinus::code) ==
......@@ -91,6 +91,14 @@ TEST_CASE("QgsjetII", "[processes]") {
CHECK(corsika::qgsjetII::getQgsjetIIXSCode(Code::PiMinus) ==
corsika::qgsjetII::QgsjetIIXSClass::LightMesons);
}
SECTION("valid") {
corsika::qgsjetII::InteractionModel model;
CHECK_THROWS(model.isValid(Code::Electron, Code::Proton));
CHECK_THROWS(model.isValid(Code::Proton, Code::Electron));
}
}
#include <corsika/framework/geometry/Point.hpp>
......@@ -115,25 +123,23 @@ TEST_CASE("QgsjetIIInterface", "interaction,processes") {
logging::set_level(logging::level::info);
auto [env, csPtr, nodePtr] = setup::testing::setup_environment(Code::Oxygen);
auto const& cs = *csPtr;
[[maybe_unused]] auto const& env_dummy = env;
[[maybe_unused]] auto const& node_dummy = nodePtr;
RNGManager<>::getInstance().registerRandomStream("qgsjet");
SECTION("InteractionInterface") {
auto [stackPtr, secViewPtr] = setup::testing::setup_stack(
Code::Proton, 110_GeV, (setup::Environment::BaseNodeType* const)nodePtr, *csPtr);
setup::StackView& view = *(secViewPtr.get());
auto particle = stackPtr->first();
auto projectile = secViewPtr->getProjectile();
auto const projectileMomentum = projectile.getMomentum();
corsika::qgsjetII::Interaction model;
model.doInteraction(view);
[[maybe_unused]] const GrammageType length = model.getInteractionLength(particle);
CHECK(length / (1_g / square(1_cm)) == Approx(93.04).margin(0.1));
corsika::qgsjetII::InteractionModel model;
model.doInteraction(view, Code::Proton, Code::Oxygen,
{sqrt(static_pow<2>(110_GeV) + static_pow<2>(Proton::mass)),
MomentumVector{cs, 110_GeV, 0_GeV, 0_GeV}},
{Oxygen::mass, MomentumVector{cs, {0_eV, 0_eV, 0_eV}}});
/* **********************************
As it turned out already two times (#291 and #307) that the detailed output of
......@@ -152,24 +158,25 @@ TEST_CASE("QgsjetIIInterface", "interaction,processes") {
SECTION("InteractionInterface Nuclei") {
HEPEnergyType const P0 = 20100_GeV;
MomentumVector const plab = MomentumVector(cs, {P0, 0_eV, 0_eV});
Code const pid = get_nucleus_code(60, 30);
auto [stackPtr, secViewPtr] = setup::testing::setup_stack(
get_nucleus_code(60, 30), 20100_GeV,
(setup::Environment::BaseNodeType* const)nodePtr, *csPtr);
pid, P0, (setup::Environment::BaseNodeType* const)nodePtr, *csPtr);
setup::StackView& view = *(secViewPtr.get());
auto particle = stackPtr->first();
auto projectile = secViewPtr->getProjectile();
auto const projectileMomentum = projectile.getMomentum();
corsika::qgsjetII::Interaction model;
model.doInteraction(view); // this also should produce some fragments
HEPEnergyType const Elab = sqrt(static_pow<2>(P0) + static_pow<2>(get_mass(pid)));
FourMomentum const projectileP4(Elab, plab);
FourMomentum const targetP4(Oxygen::mass, MomentumVector(cs, {0_eV, 0_eV, 0_eV}));
view.clear();
corsika::qgsjetII::InteractionModel model;
model.doInteraction(view, pid, Code::Oxygen, projectileP4,
targetP4); // this also should produce some fragments
CHECK(view.getSize() == Approx(300).margin(150)); // this is not physics validation
int countFragments = 0;
for (auto const& sec : view) { countFragments += (is_nucleus(sec.getPID())); }
CHECK(countFragments == Approx(4).margin(2)); // this is not physics validation
[[maybe_unused]] const GrammageType length = model.getInteractionLength(particle);
CHECK(length / (1_g / square(1_cm)) ==
Approx(12).margin(2)); // this is not physics validation
}
SECTION("Heavy nuclei") {
......@@ -178,38 +185,33 @@ TEST_CASE("QgsjetIIInterface", "interaction,processes") {
get_nucleus_code(1000, 1000), 1100_GeV,
(setup::Environment::BaseNodeType* const)nodePtr, *csPtr);
setup::StackView& view = *(secViewPtr.get());
auto particle = stackPtr->first();
auto projectile = secViewPtr->getProjectile();
auto const projectileMomentum = projectile.getMomentum();
corsika::qgsjetII::Interaction model;
corsika::qgsjetII::InteractionModel model;
FourMomentum const aP4(100_GeV, {cs, 99_GeV, 0_GeV, 0_GeV});
FourMomentum const bP4(1_TeV, {cs, 0.9_TeV, 0_GeV, 0_GeV});
CHECK_THROWS(model.getCrossSection(get_nucleus_code(10, 5),
get_nucleus_code(1000, 500), aP4, bP4));
CHECK_THROWS(model.getCrossSection(Code::Nucleus, Code::Nucleus, aP4, bP4));
CHECK_THROWS(
model.getCrossSection(Code::Nucleus, Code::Nucleus, 100_GeV, 10., 1000.));
CHECK_THROWS(
model.getCrossSection(Code::Nucleus, Code::Nucleus, 100_GeV, 1000., 10.));
CHECK_THROWS(model.doInteraction(view));
CHECK_THROWS(model.getInteractionLength(particle));
model.doInteraction(view, get_nucleus_code(1000, 500), Code::Oxygen, aP4, bP4));
}
SECTION("Allowed Particles") {
{ // electron
auto [stackPtr, secViewPtr] = setup::testing::setup_stack(
Code::Electron, 100_GeV, (setup::Environment::BaseNodeType* const)nodePtr,
*csPtr);
[[maybe_unused]] setup::StackView& view = *(secViewPtr.get());
auto particle = stackPtr->first();
corsika::qgsjetII::Interaction model;
GrammageType const length = model.getInteractionLength(particle);
CHECK(length / (1_g / square(1_cm)) == std::numeric_limits<double>::infinity());
}
{ // pi0 is internally converted into pi+/pi-
auto [stackPtr, secViewPtr] = setup::testing::setup_stack(
Code::Pi0, 1000_GeV, (setup::Environment::BaseNodeType* const)nodePtr, *csPtr);
[[maybe_unused]] setup::StackView& view = *(secViewPtr.get());
[[maybe_unused]] auto particle = stackPtr->first();
corsika::qgsjetII::Interaction model;
model.doInteraction(view);
corsika::qgsjetII::InteractionModel model;
model.doInteraction(view, Code::Pi0, Code::Oxygen,
{sqrt(static_pow<2>(1_TeV) + static_pow<2>(Pi0::mass)),
MomentumVector{cs, 1_TeV, 0_GeV, 0_GeV}},
{Oxygen::mass, MomentumVector{cs, 0_eV, 0_eV, 0_eV}});
CHECK(view.getSize() == Approx(10).margin(8)); // this is not physics validation
}
{ // rho0 is internally converted into pi-/pi+
......@@ -217,8 +219,11 @@ TEST_CASE("QgsjetIIInterface", "interaction,processes") {
Code::Rho0, 1000_GeV, (setup::Environment::BaseNodeType* const)nodePtr, *csPtr);
[[maybe_unused]] setup::StackView& view = *(secViewPtr.get());
[[maybe_unused]] auto particle = stackPtr->first();
corsika::qgsjetII::Interaction model;
model.doInteraction(view);
corsika::qgsjetII::InteractionModel model;
model.doInteraction(view, Code::Rho0, Code::Oxygen,
{sqrt(static_pow<2>(1_TeV) + static_pow<2>(Rho0::mass)),
MomentumVector{cs, 1_TeV, 0_GeV, 0_GeV}},
{Oxygen::mass, MomentumVector{cs, 0_eV, 0_eV, 0_eV}});
CHECK(view.getSize() == Approx(25).margin(20)); // this is not physics validation
}
{ // Lambda is internally converted into neutron
......@@ -227,8 +232,11 @@ TEST_CASE("QgsjetIIInterface", "interaction,processes") {
*csPtr);
[[maybe_unused]] setup::StackView& view = *(secViewPtr.get());
[[maybe_unused]] auto particle = stackPtr->first();
corsika::qgsjetII::Interaction model;
model.doInteraction(view);
corsika::qgsjetII::InteractionModel model;
model.doInteraction(view, Code::Lambda0, Code::Oxygen,
{sqrt(static_pow<2>(100_GeV) + static_pow<2>(Lambda0::mass)),
MomentumVector{cs, 100_GeV, 0_GeV, 0_GeV}},
{Oxygen::mass, MomentumVector{cs, 0_eV, 0_eV, 0_eV}});
CHECK(view.getSize() == Approx(25).margin(20)); // this is not physics validation
}
{ // AntiLambda is internally converted into anti neutron
......@@ -237,8 +245,11 @@ TEST_CASE("QgsjetIIInterface", "interaction,processes") {
*csPtr);
[[maybe_unused]] setup::StackView& view = *(secViewPtr.get());
[[maybe_unused]] auto particle = stackPtr->first();
corsika::qgsjetII::Interaction model;
model.doInteraction(view);
corsika::qgsjetII::InteractionModel model;
model.doInteraction(view, Code::Lambda0Bar, Code::Oxygen,
{sqrt(static_pow<2>(1_TeV) + static_pow<2>(Lambda0Bar::mass)),
MomentumVector{cs, 1_TeV, 0_GeV, 0_GeV}},
{Oxygen::mass, MomentumVector{cs, 0_eV, 0_eV, 0_eV}});
CHECK(view.getSize() == Approx(70).margin(67)); // this is not physics validation
}
}
......
tests/modules/testSibyll.cpp
+ 3
− 4
View file @ d7ac9f8a
......@@ -88,7 +88,6 @@ TEST_CASE("Sibyll", "modules") {
#include <corsika/framework/geometry/Vector.hpp>
#include <corsika/framework/core/PhysicalUnits.hpp>
#include <corsika/framework/core/ParticleProperties.hpp>
#include <SetupTestEnvironment.hpp>
......@@ -186,9 +185,9 @@ TEST_CASE("SibyllInterface", "modules") {
model.setVerbose(true);
HEPEnergyType const Elab = sqrt(static_pow<2>(P0) + static_pow<2>(Proton::mass));
FourMomentum const projectileP4(Elab, plab);
FourMomentum const nucleonP4(Oxygen::mass, MomentumVector(cs, {0_eV, 0_eV, 0_eV}));
FourMomentum const nucleusP4(Oxygen::mass, MomentumVector(cs, {0_eV, 0_eV, 0_eV}));
view.clear();
model.doInteraction(view, Code::Proton, Code::Oxygen, projectileP4, nucleonP4);
model.doInteraction(view, Code::Proton, Code::Oxygen, projectileP4, nucleusP4);
auto const pSum = sumMomentum(view, cs);
/*
......@@ -255,7 +254,7 @@ TEST_CASE("SibyllInterface", "modules") {
Approx(0).margin(plab.getNorm() * 0.05 / 1_GeV));
CHECK(pSum.getNorm() / P0 == Approx(1).margin(0.05));
[[maybe_unused]] CrossSectionType const cx =
model.getCrossSection(Code::Proton, Code::Oxygen, projectileP4, nucleonP4);
model.getCrossSection(Code::Proton, Code::Oxygen, projectileP4, nucleusP4);
CHECK(cx / 1_mb == Approx(300).margin(1));
// CHECK(view.getEntries() == 9); //! \todo: this was 20 before refactory-2020: check
// "also sibyll not stable wrt. to compiler
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment

data privacy - impressum