testCOMBoost.cc

/*
 * (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.
 */

#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());
};

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();
};

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, 20_GeV, 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