IAP GITLAB

Commit 28713318 authored by André Schmidt's avatar André Schmidt

Merge branch 'andre-master' into 'master'

Test

See merge request !3
parents 8cf7ec86 b371555b
Pipeline #1913 failed 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
......@@ -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;
}
}
......
......@@ -86,6 +86,15 @@ TEST_CASE("Sibyll", "[processes]") {
using namespace corsika::units::si;
using namespace corsika::units;
template <typename TStackView>
auto sumMomentum(TStackView const& view, geometry::CoordinateSystem const& vCS) {
geometry::Vector<hepenergy_d> sum{vCS, 0_eV, 0_eV, 0_eV};
for (auto const& p : view) { sum += p.GetMomentum(); }
return sum;
}
TEST_CASE("SibyllInterface", "[processes]") {
// setup environment, geometry
......@@ -113,10 +122,10 @@ TEST_CASE("SibyllInterface", "[processes]") {
SECTION("InteractionInterface") {
setup::Stack stack;
const HEPEnergyType E0 = 100_GeV;
const HEPEnergyType E0 = 60_GeV;
HEPMomentumType P0 =
sqrt(E0 * E0 - particles::Proton::GetMass() * particles::Proton::GetMass());
auto plab = corsika::stack::MomentumVector(cs, {0_GeV, 0_GeV, -P0});
auto plab = corsika::stack::MomentumVector(cs, {P0, 0_eV, 0_eV});
geometry::Point pos(cs, 0_m, 0_m, 0_m);
auto particle = stack.AddParticle(
std::tuple<particles::Code, units::si::HEPEnergyType,
......@@ -130,6 +139,18 @@ TEST_CASE("SibyllInterface", "[processes]") {
model.Init();
[[maybe_unused]] const process::EProcessReturn ret = model.DoInteraction(projectile);
auto const pSum = sumMomentum(view, cs);
// for debugging!
std::cout << pSum.GetComponents(cs) << std::endl;
std::cout << plab.GetComponents(cs) << std::endl;
CHECK(pSum.GetComponents(cs).GetX() / P0 == Approx(1).margin(1e-4));
CHECK(pSum.GetComponents(cs).GetY() / 1_GeV == Approx(0).margin(1e-4));
CHECK(pSum.GetComponents(cs).GetZ() / 1_GeV == Approx(0).margin(1e-4));
CHECK((pSum - plab).norm() / 1_GeV == Approx(0).margin(1e-4));
CHECK(pSum.norm() / P0 == Approx(1).margin(1e-4));
[[maybe_unused]] const GrammageType length = model.GetInteractionLength(particle);
}
......
Markdown is supported
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