IAP GITLAB

Commit 74d04301 authored by Andre Schmidt's avatar Andre Schmidt

Test

parent 8cf7ec86
Pipeline #1908 canceled with stages
......@@ -226,7 +226,7 @@ int main(int argc, char** argv) {
decaySibyll.PrintDecayConfig();
process::particle_cut::ParticleCut cut{60_GeV};
process::particle_cut::ParticleCut cut{50_GeV};
process::energy_loss::EnergyLoss eLoss(showerAxis);
process::longitudinal_profile::LongitudinalProfile longprof{showerAxis};
......
......@@ -8,8 +8,7 @@
* the license.
*/
#ifndef _include_COORDINATESYSTEM_H_
#define _include_COORDINATESYSTEM_H_
#pragma once
#include <corsika/geometry/QuantityVector.h>
#include <corsika/units/PhysicalUnits.h>
......@@ -118,8 +117,12 @@ namespace corsika::geometry {
auto const* GetReference() const { return reference; }
auto const& GetTransform() const { return transf; }
bool operator==(CoordinateSystem const& cs) const {
return reference == cs.reference && transf.matrix() == cs.transf.matrix();
}
bool operator!=(CoordinateSystem const& cs) const { return !(cs == *this); }
};
} // namespace corsika::geometry
#endif
/*
* (c) Copyright 2018 CORSIKA Project, corsika-project@lists.kit.edu
*
* See file AUTHORS for a list of contributors.
*
* This software is distributed under the terms of the GNU General Public
* Licence version 3 (GPL Version 3). See file LICENSE for a full version of
* the license.
*/
#ifndef _include_COORDINATESYSTEM_H_
#define _include_COORDINATESYSTEM_H_
#include <corsika/geometry/QuantityVector.h>
#include <corsika/units/PhysicalUnits.h>
#include <corsika/utl/sgn.h>
#include <Eigen/Dense>
#include <stdexcept>
typedef Eigen::Transform<double, 3, Eigen::Affine> EigenTransform;
typedef Eigen::Translation<double, 3> EigenTranslation;
namespace corsika::geometry {
class RootCoordinateSystem;
template <typename T>
class Vector;
using corsika::units::si::length_d;
class CoordinateSystem {
CoordinateSystem const* reference = nullptr;
EigenTransform transf;
CoordinateSystem(CoordinateSystem const& reference, EigenTransform const& transf)
: reference(&reference)
, transf(transf) {}
CoordinateSystem()
: // for creating the root CS
transf(EigenTransform::Identity()) {}
protected:
static auto CreateCS() { return CoordinateSystem(); }
friend corsika::geometry::RootCoordinateSystem; /// this is the only class that can
/// create ONE unique root CS
public:
static EigenTransform GetTransformation(CoordinateSystem const& c1,
CoordinateSystem const& c2);
auto& operator=(const CoordinateSystem& pCS) {
reference = pCS.reference;
transf = pCS.transf;
return *this;
}
auto translate(QuantityVector<length_d> vector) const {
EigenTransform const translation{EigenTranslation(vector.eVector)};
return CoordinateSystem(*this, translation);
}
template <typename TDim>
auto RotateToZ(Vector<TDim> vVec) const {
auto const a = vVec.normalized().GetComponents(*this).eVector;
auto const a1 = a(0), a2 = a(1);
auto const s = utl::sgn(a(2));
auto const c = 1 / (1 + s * a(2));
Eigen::Matrix3d A, B;
if (s > 0) {
A << 1, 0, a1, // comment to prevent clang-format
0, 1, a2, // .
-a1, -a2, 1; // .
B << -a1 * a1 * c, -a1 * a2 * c, 0, // .
-a1 * a2 * c, -a2 * a2 * c, 0, // .
0, 0, -(a1 * a1 + a2 * a2) * c; // .
} else {
A << 1, 0, a1, // .
0, -1, a2, // .
a1, -a2, -1; // .
B << -a1 * a1 * c, +a1 * a2 * c, 0, // .
-a1 * a2 * c, +a2 * a2 * c, 0, // .
0, 0, (a1 * a1 + a2 * a2) * c; // .
}
return CoordinateSystem(*this, EigenTransform(A + B));
}
template <typename TDim>
auto rotate(QuantityVector<TDim> axis, double angle) const {
if (axis.eVector.isZero()) {
throw std::runtime_error("null-vector given as axis parameter");
}
EigenTransform const rotation{Eigen::AngleAxisd(angle, axis.eVector.normalized())};
return CoordinateSystem(*this, rotation);
}
template <typename TDim>
auto translateAndRotate(QuantityVector<phys::units::length_d> translation,
QuantityVector<TDim> axis, double angle) {
if (axis.eVector.isZero()) {
throw std::runtime_error("null-vector given as axis parameter");
}
EigenTransform const transf{Eigen::AngleAxisd(angle, axis.eVector.normalized()) *
EigenTranslation(translation.eVector)};
return CoordinateSystem(*this, transf);
}
auto const* GetReference() const { return reference; }
auto const& GetTransform() const { return transf; }
};
} // namespace corsika::geometry
#endif
/*
* (c) Copyright 2018 CORSIKA Project, corsika-project@lists.kit.edu
*
* See file AUTHORS for a list of contributors.
*
* This software is distributed under the terms of the GNU General Public
* Licence version 3 (GPL Version 3). See file LICENSE for a full version of
* the license.
*/
#include <catch2/catch.hpp>
#include <corsika/geometry/FourVector.h>
#include <corsika/geometry/RootCoordinateSystem.h>
#include <corsika/geometry/Vector.h>
#include <corsika/units/PhysicalUnits.h>
#include <corsika/utl/COMBoost.h>
#include <iostream>
using namespace corsika::geometry;
using namespace corsika::utl;
using namespace corsika::units::si;
using corsika::units::constants::c;
using corsika::units::constants::cSquared;
double constexpr absMargin = 1e-6;
CoordinateSystem const& rootCS =
RootCoordinateSystem::GetInstance().GetRootCoordinateSystem();
// helper function for energy-momentum
// relativistic energy
auto const energy = [](HEPMassType m, Vector<hepmomentum_d> const& p) {
return sqrt(m * m + p.squaredNorm());
};
// helper function for mandelstam-s
auto const s = [](HEPEnergyType E, QuantityVector<hepmomentum_d> const& p) {
return E * E - p.squaredNorm();
};
TEST_CASE("boosts") {
// define target kinematics in lab frame
HEPMassType const targetMass = 1_GeV;
Vector<hepmomentum_d> pTargetLab{rootCS, {0_eV, 0_eV, 0_eV}};
HEPEnergyType const eTargetLab = energy(targetMass, pTargetLab);
/*
General tests check the interface and basic operation
*/
SECTION("General tests") {
// define projectile kinematics in lab frame
HEPMassType const projectileMass = 1._GeV;
Vector<hepmomentum_d> pProjectileLab{rootCS, {0_GeV, 1_PeV, 0_GeV}};
HEPEnergyType const eProjectileLab = energy(projectileMass, pProjectileLab);
const FourVector PprojLab(eProjectileLab, pProjectileLab);
// define boost to com frame
COMBoost boost(PprojLab, targetMass);
// boost projecticle
auto const PprojCoM = boost.toCoM(PprojLab);
// boost target
auto const PtargCoM = boost.toCoM(FourVector(targetMass, pTargetLab));
// sum of momenta in CoM, should be 0
auto const sumPCoM =
PprojCoM.GetSpaceLikeComponents() + PtargCoM.GetSpaceLikeComponents();
CHECK(sumPCoM.norm() / 1_GeV == Approx(0).margin(absMargin));
// mandelstam-s should be invariant under transformation
CHECK(s(eProjectileLab + eTargetLab,
pProjectileLab.GetComponents() + pTargetLab.GetComponents()) /
1_GeV / 1_GeV ==
Approx(s(PprojCoM.GetTimeLikeComponent() + PtargCoM.GetTimeLikeComponent(),
PprojCoM.GetSpaceLikeComponents().GetComponents() +
PtargCoM.GetSpaceLikeComponents().GetComponents()) /
1_GeV / 1_GeV));
// boost back...
auto const PprojBack = boost.fromCoM(PprojCoM);
// ...should yield original values before the boosts
CHECK(PprojBack.GetTimeLikeComponent() / PprojLab.GetTimeLikeComponent() ==
Approx(1));
CHECK(
(PprojBack.GetSpaceLikeComponents() - PprojLab.GetSpaceLikeComponents()).norm() /
PprojLab.GetSpaceLikeComponents().norm() ==
Approx(0).margin(absMargin));
}
/*
special case: projectile along -z
*/
SECTION("Test boost along z-axis") {
// define projectile kinematics in lab frame
HEPMassType const projectileMass = 1_GeV;
Vector<hepmomentum_d> pProjectileLab{rootCS, {0_GeV, 0_PeV, -1_PeV}};
HEPEnergyType const eProjectileLab = energy(projectileMass, pProjectileLab);
const FourVector PprojLab(eProjectileLab, pProjectileLab);
// define boost to com frame
COMBoost boost(PprojLab, targetMass);
// boost projecticle
auto const PprojCoM = boost.toCoM(PprojLab);
// boost target
auto const PtargCoM = boost.toCoM(FourVector(targetMass, pTargetLab));
// sum of momenta in CoM, should be 0
auto const sumPCoM =
PprojCoM.GetSpaceLikeComponents() + PtargCoM.GetSpaceLikeComponents();
CHECK(sumPCoM.norm() / 1_GeV == Approx(0).margin(absMargin));
}
/*
special case: projectile with arbitrary direction
*/
SECTION("Test boost along tilted axis") {
const HEPMomentumType P0 = 1_PeV;
double theta = 33.;
double phi = -10.;
auto momentumComponents = [](double theta, double phi, HEPMomentumType ptot) {
return std::make_tuple(ptot * sin(theta) * cos(phi), ptot * sin(theta) * sin(phi),
-ptot * cos(theta));
};
auto const [px, py, pz] =
momentumComponents(theta / 180. * M_PI, phi / 180. * M_PI, P0);
// define projectile kinematics in lab frame
HEPMassType const projectileMass = 1_GeV;
Vector<hepmomentum_d> pProjectileLab(rootCS, {px, py, pz});
HEPEnergyType const eProjectileLab = energy(projectileMass, pProjectileLab);
const FourVector PprojLab(eProjectileLab, pProjectileLab);
// define boost to com frame
COMBoost boost(PprojLab, targetMass);
// boost projecticle
auto const PprojCoM = boost.toCoM(PprojLab);
// boost target
auto const PtargCoM = boost.toCoM(FourVector(targetMass, pTargetLab));
// sum of momenta in CoM, should be 0
auto const sumPCoM =
PprojCoM.GetSpaceLikeComponents() + PtargCoM.GetSpaceLikeComponents();
CHECK(sumPCoM.norm() / 1_GeV == Approx(0).margin(absMargin));
}
/*
test the ultra-high energy behaviour: E=ZeV
*/
SECTION("High energy") {
// define projectile kinematics in lab frame
HEPMassType const projectileMass = 1_GeV;
HEPMomentumType P0 = 1_ZeV;
Vector<hepmomentum_d> pProjectileLab{rootCS, {0_GeV, 0_PeV, -P0}};
HEPEnergyType const eProjectileLab = energy(projectileMass, pProjectileLab);
const FourVector PprojLab(eProjectileLab, pProjectileLab);
// define boost to com frame
COMBoost boost(PprojLab, targetMass);
// boost projecticle
auto const PprojCoM = boost.toCoM(PprojLab);
// boost target
auto const PtargCoM = boost.toCoM(FourVector(targetMass, pTargetLab));
// sum of momenta in CoM, should be 0
auto const sumPCoM =
PprojCoM.GetSpaceLikeComponents() + PtargCoM.GetSpaceLikeComponents();
CHECK(sumPCoM.norm() / P0 == Approx(0).margin(absMargin)); // MAKE RELATIVE CHECK
}
SECTION("rest frame") {
HEPMassType const projectileMass = 1_GeV;
HEPMomentumType const P0 = 1_TeV;
Vector<hepmomentum_d> pProjectileLab{rootCS, {0_GeV, P0, 0_GeV}};
HEPEnergyType const eProjectileLab = energy(projectileMass, pProjectileLab);
const FourVector PprojLab(eProjectileLab, pProjectileLab);
COMBoost boostRest(pProjectileLab, projectileMass);
auto const& csPrime = boostRest.GetRotatedCS();
FourVector const rest4Mom = boostRest.toCoM(PprojLab);
CHECK(rest4Mom.GetTimeLikeComponent() / 1_GeV == Approx(projectileMass / 1_GeV));
CHECK(rest4Mom.GetSpaceLikeComponents().norm() / 1_GeV ==
Approx(0).margin(absMargin));
FourVector const a{0_eV, Vector{csPrime, 0_eV, 5_GeV, 0_eV}};
FourVector const b{0_eV, Vector{rootCS, 3_GeV, 0_eV, 0_eV}};
auto const aLab = boostRest.fromCoM(a);
auto const bLab = boostRest.fromCoM(b);
CHECK(aLab.GetNorm() / a.GetNorm() == Approx(1));
CHECK(aLab.GetSpaceLikeComponents().GetComponents(csPrime)[1].magnitude() ==
Approx((5_GeV).magnitude()));
CHECK(bLab.GetSpaceLikeComponents().GetComponents(rootCS)[0].magnitude() ==
Approx((3_GeV).magnitude()));
}
}
......@@ -35,6 +35,8 @@ auto const energy = [](HEPMassType m, Vector<hepmomentum_d> const& p) {
return sqrt(m * m + p.squaredNorm());
};
auto const momentum = [](HEPEnergyType E, HEPMassType m) { return sqrt(E * E - m * m); };
// helper function for mandelstam-s
auto const s = [](HEPEnergyType E, QuantityVector<hepmomentum_d> const& p) {
return E * E - p.squaredNorm();
......@@ -54,7 +56,7 @@ TEST_CASE("boosts") {
// define projectile kinematics in lab frame
HEPMassType const projectileMass = 1._GeV;
Vector<hepmomentum_d> pProjectileLab{rootCS, {0_GeV, 1_PeV, 0_GeV}};
Vector<hepmomentum_d> pProjectileLab{rootCS, {0_GeV, 20_GeV, 0_GeV}};
HEPEnergyType const eProjectileLab = energy(projectileMass, pProjectileLab);
const FourVector PprojLab(eProjectileLab, pProjectileLab);
......@@ -101,18 +103,26 @@ TEST_CASE("boosts") {
// define projectile kinematics in lab frame
HEPMassType const projectileMass = 1_GeV;
Vector<hepmomentum_d> pProjectileLab{rootCS, {0_GeV, 0_PeV, -1_PeV}};
Vector<hepmomentum_d> pProjectileLab{rootCS, {0_GeV, 0_GeV, -20_GeV}};
HEPEnergyType const eProjectileLab = energy(projectileMass, pProjectileLab);
const FourVector PprojLab(eProjectileLab, pProjectileLab);
auto const sqrt_s_lab =
sqrt(s(eProjectileLab + targetMass, pProjectileLab.GetComponents(rootCS)));
// define boost to com frame
COMBoost boost(PprojLab, targetMass);
// boost projecticle
auto const PprojCoM = boost.toCoM(PprojLab);
auto const a = PprojCoM.GetSpaceLikeComponents().GetComponents(boost.GetRotatedCS());
CHECK(a.GetX() / 1_GeV == Approx(0));
CHECK(a.GetY() / 1_GeV == Approx(0));
CHECK(a.GetZ() / (momentum(sqrt_s_lab / 2, projectileMass)) == Approx(1));
// boost target
auto const PtargCoM = boost.toCoM(FourVector(targetMass, pTargetLab));
CHECK(PtargCoM.GetTimeLikeComponent() / sqrt_s_lab == Approx(.5));
// sum of momenta in CoM, should be 0
auto const sumPCoM =
......@@ -183,31 +193,30 @@ TEST_CASE("boosts") {
PprojCoM.GetSpaceLikeComponents() + PtargCoM.GetSpaceLikeComponents();
CHECK(sumPCoM.norm() / P0 == Approx(0).margin(absMargin)); // MAKE RELATIVE CHECK
}
}
SECTION("rest frame") {
HEPMassType const projectileMass = 1_GeV;
HEPMomentumType const P0 = 1_TeV;
Vector<hepmomentum_d> pProjectileLab{rootCS, {0_GeV, P0, 0_GeV}};
HEPEnergyType const eProjectileLab = energy(projectileMass, pProjectileLab);
const FourVector PprojLab(eProjectileLab, pProjectileLab);
COMBoost boostRest(pProjectileLab, projectileMass);
auto const& csPrime = boostRest.GetRotatedCS();
FourVector const rest4Mom = boostRest.toCoM(PprojLab);
CHECK(rest4Mom.GetTimeLikeComponent() / 1_GeV == Approx(projectileMass / 1_GeV));
CHECK(rest4Mom.GetSpaceLikeComponents().norm() / 1_GeV ==
Approx(0).margin(absMargin));
FourVector const a{0_eV, Vector{csPrime, 0_eV, 5_GeV, 0_eV}};
FourVector const b{0_eV, Vector{rootCS, 3_GeV, 0_eV, 0_eV}};
auto const aLab = boostRest.fromCoM(a);
auto const bLab = boostRest.fromCoM(b);
CHECK(aLab.GetNorm() / a.GetNorm() == Approx(1));
CHECK(aLab.GetSpaceLikeComponents().GetComponents(csPrime)[1].magnitude() ==
Approx((5_GeV).magnitude()));
CHECK(bLab.GetSpaceLikeComponents().GetComponents(rootCS)[0].magnitude() ==
Approx((3_GeV).magnitude()));
}
TEST_CASE("rest frame") {
HEPMassType const projectileMass = 1_GeV;
HEPMomentumType const P0 = 1_TeV;
Vector<hepmomentum_d> pProjectileLab{rootCS, {0_GeV, P0, 0_GeV}};
HEPEnergyType const eProjectileLab = energy(projectileMass, pProjectileLab);
const FourVector PprojLab(eProjectileLab, pProjectileLab);
COMBoost boostRest(pProjectileLab, projectileMass);
auto const& csPrime = boostRest.GetRotatedCS();
FourVector const rest4Mom = boostRest.toCoM(PprojLab);
CHECK(rest4Mom.GetTimeLikeComponent() / 1_GeV == Approx(projectileMass / 1_GeV));
CHECK(rest4Mom.GetSpaceLikeComponents().norm() / 1_GeV == Approx(0).margin(absMargin));
FourVector const a{0_eV, Vector{csPrime, 0_eV, 5_GeV, 0_eV}};
FourVector const b{0_eV, Vector{rootCS, 3_GeV, 0_eV, 0_eV}};
auto const aLab = boostRest.fromCoM(a);
auto const bLab = boostRest.fromCoM(b);
CHECK(aLab.GetNorm() / a.GetNorm() == Approx(1));
CHECK(aLab.GetSpaceLikeComponents().GetComponents(csPrime)[1].magnitude() ==
Approx((5_GeV).magnitude()));
CHECK(bLab.GetSpaceLikeComponents().GetComponents(rootCS)[0].magnitude() ==
Approx((3_GeV).magnitude()));
}
......@@ -164,9 +164,8 @@ namespace corsika::process::sibyll {
template <>
process::EProcessReturn Interaction::DoInteraction(SetupProjectile& vP) {
using namespace units;
using namespace utl;
using namespace units;
using namespace units::si;
using namespace geometry;
......@@ -181,24 +180,22 @@ namespace corsika::process::sibyll {
}
if (process::sibyll::CanInteract(corsikaBeamId)) {
const CoordinateSystem& rootCS =
RootCoordinateSystem::GetInstance().GetRootCoordinateSystem();
// position and time of interaction, not used in Sibyll
Point pOrig = vP.GetPosition();
TimeType tOrig = vP.GetTime();
Point const pOrig = vP.GetPosition();
TimeType const tOrig = vP.GetTime();
// define projectile
HEPEnergyType const eProjectileLab = vP.GetEnergy();
auto const pProjectileLab = vP.GetMomentum();
const CoordinateSystem& originalCS = pProjectileLab.GetCoordinateSystem();
// 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(rootCS, 0_GeV, 0_GeV, 0_GeV);
const auto pTargetLab = MomentumVector(originalCS, 0_GeV, 0_GeV, 0_GeV);
const FourVector PtargLab(eTargetLab, pTargetLab);
// define projectile
HEPEnergyType const eProjectileLab = vP.GetEnergy();
auto const pProjectileLab = vP.GetMomentum();
cout << "Interaction: ebeam lab: " << eProjectileLab / 1_GeV << endl
<< "Interaction: pbeam lab: " << pProjectileLab.GetComponents() / 1_GeV
<< endl;
......@@ -211,6 +208,7 @@ namespace corsika::process::sibyll {
// 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
......@@ -222,11 +220,11 @@ namespace corsika::process::sibyll {
cout << "Interaction: ebeam CoM: " << PprojCoM.GetTimeLikeComponent() / 1_GeV
<< endl
<< "Interaction: pbeam CoM: "
<< PprojCoM.GetSpaceLikeComponents().GetComponents() / 1_GeV << endl;
<< PprojCoM.GetSpaceLikeComponents().GetComponents(csPrime) / 1_GeV << endl;
cout << "Interaction: etarget CoM: " << PtargCoM.GetTimeLikeComponent() / 1_GeV
<< endl
<< "Interaction: ptarget CoM: "
<< PtargCoM.GetSpaceLikeComponents().GetComponents() / 1_GeV << endl;
<< PtargCoM.GetSpaceLikeComponents().GetComponents(csPrime) / 1_GeV << endl;
cout << "Interaction: position of interaction: " << pOrig.GetCoordinates() << endl;
cout << "Interaction: time: " << tOrig << endl;
......@@ -247,7 +245,7 @@ namespace corsika::process::sibyll {
*/
//#warning reading interaction cross section again, should not be necessary
auto const& compVec = mediumComposition.GetComponents();
std::vector<si::CrossSectionType> cross_section_of_components(compVec.size());
std::vector<CrossSectionType> cross_section_of_components(compVec.size());
for (size_t i = 0; i < compVec.size(); ++i) {
auto const targetId = compVec[i];
......@@ -303,7 +301,7 @@ namespace corsika::process::sibyll {
// link to sibyll stack
SibStack ss;
MomentumVector Plab_final(rootCS, {0.0_GeV, 0.0_GeV, 0.0_GeV});
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) {
......@@ -311,25 +309,34 @@ namespace corsika::process::sibyll {
if (psib.HasDecayed())
throw std::runtime_error("found particle that decayed in SIBYLL!");
// transform energy to lab. frame
auto const pCoM = psib.GetMomentum();
// 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 = vP.AddSecondary(
tuple<particles::Code, units::si::HEPEnergyType, stack::MomentumVector,
geometry::Point, units::si::TimeType>{
process::sibyll::ConvertFromSibyll(psib.GetPID()),
Plab.GetTimeLikeComponent(), Plab.GetSpaceLikeComponents(), pOrig,
tOrig});
Plab.GetTimeLikeComponent(), p3lab, pOrig, tOrig});
Plab_final += pnew.GetMomentum();
Elab_final += pnew.GetEnergy();
Ecm_final += psib.GetEnergy();
}
cout << "conservation (all GeV): Ecm_final=" << Ecm_final / 1_GeV << endl
<< "Elab_final=" << Elab_final / 1_GeV
cout << "conservation (all GeV):" << endl
<< "Ecm_initial=" << Ecm / 1_GeV << " Ecm_final=" << Ecm_final / 1_GeV
<< endl
<< "Elab_initial=" << eProjectileLab / 1_GeV
<< " Elab_final=" << Elab_final / 1_GeV
<< " diff (%)=" << (Elab_final / eProjectileLab / get_nwounded() - 1) * 100
<< endl
<< "Plab_initial=" << (pProjectileLab / 1_GeV).GetComponents()
<< ", Plab_final=" << (Plab_final / 1_GeV).GetComponents() << endl;
}
}
......
This diff is collapsed.
/*
* (c) Copyright 2019 CORSIKA Project, corsika-project@lists.kit.edu
*
* See file AUTHORS for a list of contributors.
*
* This software is distributed under the terms of the GNU General Public
* Licence version 3 (GPL Version 3). See file LICENSE for a full version of
* the license.
*/
#include <corsika/process/sibyll/Decay.h>
#include <corsika/process/sibyll/Interaction.h>
#include <corsika/process/sibyll/NuclearInteraction.h>
#include <corsika/process/sibyll/ParticleConversion.h>
#include <corsika/random/RNGManager.h>