From 7ba0023b3b7698a2d12f3531a297e2b364ca22e1 Mon Sep 17 00:00:00 2001
From: ralfulrich <ralf.ulrich@kit.edu>
Date: Wed, 16 Jan 2019 09:15:19 +0100
Subject: [PATCH] debuggin crash in COM at high energies
---
Documentation/Examples/cascade_example.cc | 2 +-
Framework/Geometry/FourVector.h | 16 ++---
Framework/Utilities/COMBoost.cc | 75 ++++++++++++++++-------
Framework/Utilities/COMBoost.h | 24 ++++----
Framework/Utilities/testCOMBoost.cc | 29 +++++----
Processes/Sibyll/Interaction.h | 60 ++++++++++--------
6 files changed, 128 insertions(+), 78 deletions(-)
diff --git a/Documentation/Examples/cascade_example.cc b/Documentation/Examples/cascade_example.cc
index e978949c..e49a41e3 100644
--- a/Documentation/Examples/cascade_example.cc
+++ b/Documentation/Examples/cascade_example.cc
@@ -245,7 +245,7 @@ int main() {
// setup particle stack, and add primary particle
setup::Stack stack;
stack.Clear();
- const HEPEnergyType E0 = 100_TeV;
+ const HEPEnergyType E0 = 1000_TeV;
double theta = 0.;
double phi = 0.;
{
diff --git a/Framework/Geometry/FourVector.h b/Framework/Geometry/FourVector.h
index bb848758..93f29dcf 100644
--- a/Framework/Geometry/FourVector.h
+++ b/Framework/Geometry/FourVector.h
@@ -33,7 +33,7 @@ namespace corsika::geometry {
public:
using SpaceType = typename std::decay<SpaceVecType>::type::Quantity;
- /// check the types and the physical units here:
+ //! check the types and the physical units here:
static_assert(
std::is_same<typename std::decay<TimeType>::type, SpaceType>::value ||
std::is_same<typename std::decay<TimeType>::type,
@@ -47,7 +47,9 @@ namespace corsika::geometry {
: fTimeLike(eT)
, fSpaceLike(eS) {}
- TimeType GetTime() { return fTimeLike; }
+ TimeType GetTimeLikeComponent() const { return fTimeLike; }
+ SpaceVecType& GetSpaceLikeComponents() { return fSpaceLike; }
+ const SpaceVecType& GetSpaceLikeComponents() const { return fSpaceLike; }
auto GetNormSqr() const { return GetTimeSquared() - fSpaceLike.squaredNorm(); }
@@ -55,15 +57,15 @@ namespace corsika::geometry {
bool IsTimelike() const {
return GetTimeSquared() < fSpaceLike.squaredNorm();
- } /// Norm2 < 0
+ } //! Norm2 < 0
bool IsSpacelike() const {
return GetTimeSquared() > fSpaceLike.squaredNorm();
- } /// Norm2 > 0
+ } //! Norm2 > 0
bool IsPhotonlike() const {
return GetTimeSquared() == fSpaceLike.squaredNorm();
- } /// Norm2 == 0
+ } //! Norm2 == 0
FourVector& operator+=(const FourVector& b) {
fTimeLike += b.fTimeLike;
@@ -129,11 +131,11 @@ namespace corsika::geometry {
}
protected:
- /// the data members
+ //! the data members
TimeType fTimeLike;
SpaceVecType fSpaceLike;
- /// the friends: math operators
+ //! the friends: math operators
template <typename T, typename U>
friend FourVector<typename std::decay<T>::type, typename std::decay<U>::type>
operator+(const FourVector<T, U>&, const FourVector<T, U>&);
diff --git a/Framework/Utilities/COMBoost.cc b/Framework/Utilities/COMBoost.cc
index 1bf8f839..3e9aafb8 100644
--- a/Framework/Utilities/COMBoost.cc
+++ b/Framework/Utilities/COMBoost.cc
@@ -8,16 +8,20 @@
* the license.
*/
+#include <corsika/geometry/CoordinateSystem.h>
+#include <corsika/geometry/Vector.h>
+#include <corsika/units/PhysicalUnits.h>
#include <corsika/utl/COMBoost.h>
using namespace corsika::utl;
using namespace corsika::units::si;
-COMBoost::COMBoost(HEPEnergyType eProjectile, COMBoost::MomentumVector const& pProjectile,
- HEPMassType mTarget)
+template <typename FourVector>
+COMBoost<FourVector>::COMBoost(const FourVector& Pprojectile, const FourVector& Ptarget)
: fRotation(Eigen::Matrix3d::Identity())
- , fCS(pProjectile.GetCoordinateSystem()) {
+ , fCS(Pprojectile.GetSpaceLikeComponents().GetCoordinateSystem()) {
// calculate matrix for rotating pProjectile to z-axis first
+ auto const pProjectile = Pprojectile.GetSpaceLikeComponents();
auto const pProjNorm = pProjectile.norm();
auto const a = (pProjectile / pProjNorm).GetComponents().eVector;
@@ -37,47 +41,76 @@ COMBoost::COMBoost(HEPEnergyType eProjectile, COMBoost::MomentumVector const& pP
}
// calculate boost
- double const x = pProjNorm / (eProjectile + mTarget);
+ double const beta =
+ pProjNorm / (Pprojectile.GetTimeLikeComponent() + Ptarget.GetTimeLikeComponent());
- /* Accurracy matters here, x = 1 - epsilon for ultra-relativistic boosts */
- double const coshEta = 1 / std::sqrt((1 + x) * (1 - x));
- //~ double const coshEta = 1 / std::sqrt((1-x*x));
- double const sinhEta = -x * coshEta;
+ /* Accurracy matters here, beta = 1 - epsilon for ultra-relativistic boosts */
+ double const coshEta = 1 / std::sqrt((1 + beta) * (1 - beta));
+ //~ double const coshEta = 1 / std::sqrt((1-beta*beta));
+ double const sinhEta = -beta * coshEta;
+ std::cout << "COMBoost (1-beta)=" << 1-beta << " gamma=" << coshEta << std::endl;
+
fBoost << coshEta, sinhEta, sinhEta, coshEta;
fInverseBoost << coshEta, -sinhEta, -sinhEta, coshEta;
}
-std::tuple<HEPEnergyType, corsika::geometry::QuantityVector<hepmomentum_d>>
-COMBoost::toCoM(HEPEnergyType E, COMBoost::MomentumVector p) const {
- corsika::geometry::QuantityVector<hepmomentum_d> pComponents = p.GetComponents(fCS);
+template <typename FourVector>
+FourVector COMBoost<FourVector>::toCoM(const FourVector& p) const {
+ auto pComponents = p.GetSpaceLikeComponents().GetComponents(fCS);
Eigen::Vector3d eVecRotated = fRotation * pComponents.eVector;
Eigen::Vector2d lab;
- lab << (E * (1 / 1_GeV)), (eVecRotated(2) * (1 / 1_GeV).magnitude());
+ lab << (p.GetTimeLikeComponent() * (1 / 1_GeV)),
+ (eVecRotated(2) * (1 / 1_GeV).magnitude());
auto const boostedZ = fBoost * lab;
auto const E_CoM = boostedZ(0) * 1_GeV;
eVecRotated(2) = boostedZ(1) * (1_GeV).magnitude();
- return std::make_tuple(E_CoM,
- corsika::geometry::QuantityVector<hepmomentum_d>{eVecRotated});
+ return FourVector(E_CoM, corsika::geometry::Vector<hepmomentum_d>(fCS, eVecRotated));
}
-std::tuple<HEPEnergyType, COMBoost::MomentumVector> COMBoost::fromCoM(
- HEPEnergyType E,
- corsika::geometry::QuantityVector<units::si::hepmomentum_d> pCoM) const {
+template <typename FourVector>
+FourVector COMBoost<FourVector>::fromCoM(const FourVector& p) const {
Eigen::Vector2d com;
- com << (E * (1 / 1_GeV)), (pCoM.eVector(2) * (1 / 1_GeV).magnitude());
+ com << (p.GetTimeLikeComponent() * (1 / 1_GeV)),
+ (p.GetSpaceLikeComponents().GetComponents().eVector(2) * (1 / 1_GeV).magnitude());
+ std::cout << "COMBoost::fromCoM Ecm=" << p.GetTimeLikeComponent() / 1_GeV << "GeV, "
+ << " pcm=" << p.GetSpaceLikeComponents().GetComponents() / 1_GeV << "GeV" << std::endl;
+
auto const boostedZ = fInverseBoost * com;
- auto const E_CoM = boostedZ(0) * 1_GeV;
+ auto const E_lab = boostedZ(0) * 1_GeV;
- auto pLab = pCoM;
+ auto pLab = p.GetSpaceLikeComponents().GetComponents();
pLab.eVector(2) = boostedZ(1) * (1_GeV).magnitude();
pLab.eVector = fRotation.transpose() * pLab.eVector;
- return std::make_tuple(E_CoM, MomentumVector(fCS, pLab));
+ std::cout << "COMBoost::fromCoM --> Elab=" << E_lab / 1_GeV << "GeV, "
+ << " pcm=" << pLab / 1_GeV << "GeV" << std::endl;
+
+ return FourVector(E_lab, corsika::geometry::Vector(fCS, pLab));
}
+
+/*
+ Here we instantiate all physically meaningful versions of COMBoost
+ */
+
+#include <corsika/geometry/FourVector.h>
+
+namespace corsika::utl {
+
+ using corsika::geometry::FourVector;
+ using corsika::geometry::Vector;
+ using namespace corsika::units;
+
+ // for normal HEP energy/momentum units in [GeV]
+ template class COMBoost<FourVector<HEPEnergyType, Vector<hepmomentum_d>>>;
+ // template class COMBoost<FourVector<HEPEnergyType&, Vector<hepmomentum_d>&>>;
+
+ // template class COMBoost<FourVector<TimeType, Vector<length_d>>>;
+ // template class COMBoost<FourVector<TimeType&, Vector<length_d>&>>;
+} // namespace corsika::utl
diff --git a/Framework/Utilities/COMBoost.h b/Framework/Utilities/COMBoost.h
index fa57d8f3..08dfe75a 100644
--- a/Framework/Utilities/COMBoost.h
+++ b/Framework/Utilities/COMBoost.h
@@ -12,33 +12,31 @@
#define _include_corsika_utilties_comboost_h_
#include <corsika/geometry/CoordinateSystem.h>
-#include <corsika/geometry/QuantityVector.h>
-#include <corsika/geometry/Vector.h>
-#include <corsika/units/PhysicalUnits.h>
+
#include <Eigen/Dense>
-#include <tuple>
namespace corsika::utl {
+
+ /**
+ This utility class handles Lorentz boost between different
+ referenence frames, using FourVectors.
+ */
+
+ template <typename FourVector>
class COMBoost {
Eigen::Matrix3d fRotation;
Eigen::Matrix2d fBoost, fInverseBoost;
corsika::geometry::CoordinateSystem const& fCS;
- using MomentumVector = corsika::geometry::Vector<units::si::hepmomentum_d>;
public:
//! construct a COMBoost given energy and momentum of projectile and mass of target
- COMBoost(units::si::HEPEnergyType eProjectile, MomentumVector const& pProjectile,
- units::si::HEPMassType mTarget);
+ COMBoost(const FourVector& Pprojectile, const FourVector& Ptarget);
//! transforms a 4-momentum from lab frame to the center-of-mass frame
- std::tuple<units::si::HEPEnergyType,
- geometry::QuantityVector<units::si::hepmomentum_d>>
- toCoM(units::si::HEPEnergyType E, MomentumVector p) const;
+ FourVector toCoM(const FourVector& p) const;
//! transforms a 4-momentum from the center-of-mass frame back to lab frame
- std::tuple<units::si::HEPEnergyType, MomentumVector> fromCoM(
- units::si::HEPEnergyType E,
- geometry::QuantityVector<units::si::hepmomentum_d> pCoM) const;
+ FourVector fromCoM(const FourVector& pCoM) const;
};
} // namespace corsika::utl
diff --git a/Framework/Utilities/testCOMBoost.cc b/Framework/Utilities/testCOMBoost.cc
index 6b3b9505..bbc6c295 100644
--- a/Framework/Utilities/testCOMBoost.cc
+++ b/Framework/Utilities/testCOMBoost.cc
@@ -12,7 +12,9 @@
// cpp file
#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>
@@ -43,6 +45,7 @@ TEST_CASE("boosts") {
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 target kinematics in lab frame
HEPMassType const targetMass = 1_GeV;
@@ -50,32 +53,34 @@ TEST_CASE("boosts") {
HEPEnergyType const eTargetLab = energy(targetMass, pTargetLab);
// define boost to com frame
- COMBoost boost(eProjectileLab, pProjectileLab, targetMass);
+ COMBoost boost(PprojLab, FourVector(targetMass, pTargetLab));
// boost projecticle
- auto const [eProjectileCoM, pProjectileCoM] =
- boost.toCoM(eProjectileLab, pProjectileLab);
+ auto const PprojCoM = boost.toCoM(PprojLab);
// boost target
- auto const [eTargetCoM, pTargetCoM] = boost.toCoM(eTargetLab, pTargetLab);
+ auto const PtargCoM = boost.toCoM(FourVector(targetMass, pTargetLab));
// sum of momenta in CoM, should be 0
- auto const sumPCoM = pProjectileCoM + pTargetCoM;
- CHECK(sumPCoM[2] / 1_GeV == Approx(0).margin(absMargin));
+ 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(eProjectileCoM + eTargetCoM, pProjectileCoM + pTargetCoM) / 1_GeV /
- 1_GeV));
+ Approx(s(PprojCoM.GetTimeLikeComponent() + PtargCoM.GetTimeLikeComponent(),
+ PprojCoM.GetSpaceLikeComponents().GetComponents() +
+ PtargCoM.GetSpaceLikeComponents().GetComponents()) /
+ 1_GeV / 1_GeV));
// boost back...
- auto const [eProjectileBack, pProjectileBack] =
- boost.fromCoM(eProjectileCoM, pProjectileCoM);
+ auto const PprojBack = boost.fromCoM(PprojCoM);
// ...should yield original values before the boosts
- CHECK(eProjectileBack / eProjectileLab == Approx(1));
- CHECK((pProjectileBack - pProjectileLab).norm() / pProjectileLab.norm() ==
+ CHECK(PprojBack.GetTimeLikeComponent() / PprojLab.GetTimeLikeComponent() == Approx(1));
+ CHECK((PprojBack.GetSpaceLikeComponents() - PprojLab.GetSpaceLikeComponents()).norm() /
+ PprojLab.GetSpaceLikeComponents().norm() ==
Approx(0).margin(absMargin));
}
diff --git a/Processes/Sibyll/Interaction.h b/Processes/Sibyll/Interaction.h
index f35b244f..91702692 100644
--- a/Processes/Sibyll/Interaction.h
+++ b/Processes/Sibyll/Interaction.h
@@ -16,6 +16,7 @@
#include <corsika/environment/Environment.h>
#include <corsika/environment/NuclearComposition.h>
+#include <corsika/geometry/FourVector.h>
#include <corsika/particles/ParticleProperties.h>
#include <corsika/process/sibyll/ParticleConversion.h>
#include <corsika/process/sibyll/SibStack.h>
@@ -97,11 +98,11 @@ namespace corsika::process::sibyll {
corsika::particles::Neutron::GetMass());
// FOR NOW: assume target is at rest
- MomentumVector pTarget(rootCS, {0.0_GeV, 0.0_GeV, 0.0_GeV});
+ MomentumVector pTarget(rootCS, {0_GeV, 0_GeV, 0_GeV});
// total momentum and energy
HEPEnergyType Elab = p.GetEnergy() + nucleon_mass;
- MomentumVector pTotLab(rootCS, {0.0_GeV, 0.0_GeV, 0.0_GeV});
+ MomentumVector pTotLab(rootCS, {0_GeV, 0_GeV, 0_GeV});
pTotLab += p.GetMomentum();
pTotLab += pTarget;
auto const pTotLabNorm = pTotLab.norm();
@@ -165,6 +166,11 @@ namespace corsika::process::sibyll {
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 Particle, typename Stack>
corsika::process::EProcessReturn DoInteraction(Particle& p, Stack& s) {
@@ -178,7 +184,8 @@ namespace corsika::process::sibyll {
const auto corsikaBeamId = p.GetPID();
cout << "ProcessSibyll: "
<< "DoInteraction: " << corsikaBeamId << " interaction? "
- << process::sibyll::CanInteract(corsikaBeamId) << endl;
+ << process::sibyll::CanInteract(corsikaBeamId)
+ << endl;
if (process::sibyll::CanInteract(corsikaBeamId)) {
const CoordinateSystem& rootCS =
RootCoordinateSystem::GetInstance().GetRootCoordinateSystem();
@@ -194,16 +201,12 @@ namespace corsika::process::sibyll {
// FOR NOW: target is always at rest
const auto eTargetLab = 0_GeV + nucleon_mass;
const auto pTargetLab = MomentumVector(rootCS, 0_GeV, 0_GeV, 0_GeV);
-
+ const FourVector PtargLab(eTargetLab, pTargetLab);
+
// define projectile
- auto const pProjectileLab = p.GetMomentum();
HEPEnergyType const eProjectileLab = p.GetEnergy();
-
- // define target kinematics in lab frame
- HEPMassType const targetMass = nucleon_mass;
- // define boost to and from CoM frame
- // CoM frame definition in Sibyll projectile: +z
- COMBoost const boost(eProjectileLab, pProjectileLab, targetMass);
+ auto const pProjectileLab = p.GetMomentum();
+ const FourVector PprojLab(eProjectileLab, pProjectileLab);
cout << "Interaction: ebeam lab: " << eProjectileLab / 1_GeV << endl
<< "Interaction: pbeam lab: " << pProjectileLab.GetComponents() / 1_GeV
@@ -212,18 +215,26 @@ namespace corsika::process::sibyll {
<< "Interaction: ptarget lab: " << pTargetLab.GetComponents() / 1_GeV
<< endl;
+ // define target kinematics in lab frame
+ // define boost to and from CoM frame
+ // CoM frame definition in Sibyll projectile: +z
+ COMBoost const boost(PprojLab, PtargLab);
+
// just for show:
// boost projecticle
- auto const [eProjectileCoM, pProjectileCoM] =
- boost.toCoM(eProjectileLab, pProjectileLab);
+ auto const PprojCoM = boost.toCoM(PprojLab);
// boost target
- auto const [eTargetCoM, pTargetCoM] = boost.toCoM(eTargetLab, pTargetLab);
+ auto const PtargCoM = boost.toCoM(PtargLab);
- cout << "Interaction: ebeam CoM: " << eProjectileCoM / 1_GeV << endl
- << "Interaction: pbeam CoM: " << pProjectileCoM / 1_GeV << endl;
- cout << "Interaction: etarget CoM: " << eTargetCoM / 1_GeV << endl
- << "Interaction: ptarget CoM: " << pTargetCoM / 1_GeV << endl;
+ cout << "Interaction: ebeam CoM: " << PprojCoM.GetTimeLikeComponent() / 1_GeV
+ << endl
+ << "Interaction: pbeam CoM: "
+ << PprojCoM.GetSpaceLikeComponents().GetComponents() / 1_GeV << endl;
+ cout << "Interaction: etarget CoM: " << PtargCoM.GetTimeLikeComponent() / 1_GeV
+ << endl
+ << "Interaction: ptarget CoM: "
+ << PtargCoM.GetSpaceLikeComponents().GetComponents() / 1_GeV << endl;
cout << "Interaction: position of interaction: " << pOrig.GetCoordinates()
<< endl;
@@ -246,7 +257,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(
- mediumComposition.GetComponents().size());
+ compVec.size());
for (size_t i = 0; i < compVec.size(); ++i) {
auto const targetId = compVec[i];
@@ -279,7 +290,8 @@ namespace corsika::process::sibyll {
<< " Ecm(GeV): " << Ecm / 1_GeV << std::endl;
if (eProjectileLab < 8.5_GeV || Ecm < 10_GeV) {
std::cout << "Interaction: "
- << " DoInteraction: should have dropped particle.." << std::endl;
+ << " DoInteraction: should have dropped particle.. "
+ << "THIS IS AN ERROR" << std::endl;
// p.Delete(); delete later... different process
} else {
fCount++;
@@ -310,15 +322,15 @@ namespace corsika::process::sibyll {
if (psib.HasDecayed()) continue;
// transform energy to lab. frame
- auto const pCoM = psib.GetMomentum().GetComponents();
+ auto const pCoM = psib.GetMomentum();
HEPEnergyType const eCoM = psib.GetEnergy();
- auto const [eLab, pLab] = boost.fromCoM(eCoM, pCoM);
+ auto const Plab = boost.fromCoM(FourVector(eCoM, pCoM));
// add to corsika stack
auto pnew = s.NewParticle();
pnew.SetPID(process::sibyll::ConvertFromSibyll(psib.GetPID()));
- pnew.SetEnergy(eLab);
- pnew.SetMomentum(pLab);
+ pnew.SetEnergy(Plab.GetTimeLikeComponent());
+ pnew.SetMomentum(Plab.GetSpaceLikeComponents());
pnew.SetPosition(pOrig);
pnew.SetTime(tOrig);
--
GitLab