IAP GITLAB

Commit c10b0ac9 authored by Ralf Ulrich's avatar Ralf Ulrich

fixed all of Remy's comments

parent e3a74d98
......@@ -232,7 +232,6 @@ add_subdirectory (documentation)
#+++++++++++++++++++++++++++++
# modules
#
set (NO_BOOST 1) # do not use libboost_iostream in CorsikaData
add_subdirectory (modules/data)
add_subdirectory (modules/pythia)
add_subdirectory (modules/sibyll)
......
......@@ -26,7 +26,7 @@ namespace corsika {
}
template <typename TDimension>
QuantityVector<TDimension>& BaseVector<TDimension>::quantityVector() {
QuantityVector<TDimension>& BaseVector<TDimension>::getQuantityVector() {
return quantityVector_;
}
......
......@@ -20,76 +20,7 @@
namespace corsika {
inline CoordinateSystemPtr CoordinateSystem::translate(
QuantityVector<length_d> vector) const {
EigenTransform const translation{EigenTranslation(vector.eigenVector_)};
return std::make_shared<CoordinateSystem const>(*(new CoordinateSystem(
std::make_shared<CoordinateSystem const>(*this), translation)));
}
template <typename TDim>
CoordinateSystemPtr CoordinateSystem::rotateToZ(Vector<TDim> vVec) const {
auto const a = vVec.normalized()
.getComponents(std::make_shared<CoordinateSystem const>(*this))
.getEigenVector();
auto const a1 = a(0), a2 = a(1), a3 = a(2);
Eigen::Matrix3d A, B;
if (a3 > 0) {
auto const c = 1 / (1 + a3);
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 {
auto const c = 1 / (1 - a3);
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 std::make_shared<CoordinateSystem const>(*(new CoordinateSystem(
std::make_shared<CoordinateSystem const>(*this), EigenTransform(A + B))));
}
template <typename TDim>
CoordinateSystemPtr CoordinateSystem::rotate(QuantityVector<TDim> axis,
double angle) const {
if (axis.eigenVector_.isZero()) {
throw std::runtime_error("null-vector given as axis parameter");
}
EigenTransform const rotation{
Eigen::AngleAxisd(angle, axis.eigenVector_.normalized())};
return std::make_shared<CoordinateSystem const>(
CoordinateSystem(std::make_shared<CoordinateSystem const>(*this), rotation));
}
template <typename TDim>
CoordinateSystemPtr CoordinateSystem::translateAndRotate(
QuantityVector<length_d> translation, QuantityVector<TDim> axis, double angle) {
if (axis.eigenVector_.isZero()) {
throw std::runtime_error("null-vector given as axis parameter");
}
EigenTransform const transf{Eigen::AngleAxisd(angle, axis.eigenVector_.normalized()) *
EigenTranslation(translation.eigenVector_)};
return std::make_shared<CoordinateSystem const>(CoordinateSystem(*this, transf));
}
CoordinateSystemPtr CoordinateSystem::getReferenceCS() const {
return referenceCS_; //*(referenceCS_.get());
}
CoordinateSystemPtr CoordinateSystem::getReferenceCS() const { return referenceCS_; }
EigenTransform const& CoordinateSystem::getTransform() const { return transf_; }
......@@ -101,30 +32,27 @@ namespace corsika {
return !(cs == *this);
}
/**
* returns the transformation matrix necessary to transform primitives with coordinates
* in \a pFrom to \a pTo, e.g.
* \f$ \vec{v}^{\text{(to)}} = \mathcal{M} \vec{v}^{\text{(from)}} \f$
* (\f$ \vec{v}^{(.)} \f$ denotes the coordinates/components of the component in
* the indicated CoordinateSystem).
*/
inline EigenTransform getTransformation(CoordinateSystemPtr pFrom,
CoordinateSystemPtr pTo) {
/// find transformation between two CS, using most optimal common base
inline EigenTransform get_transformation(CoordinateSystemPtr const& pFrom,
CoordinateSystemPtr const& pTo) {
CoordinateSystemPtr a{pFrom};
CoordinateSystemPtr b{pTo};
CoordinateSystemPtr commonBase{nullptr};
while (a != b && b != nullptr) {
a = pFrom;
while (a != b && b) {
while (a != b && a != nullptr) { a = a->getReferenceCS(); }
// traverse pFrom
a = pFrom;
while (a != b && a) {
a = a->getReferenceCS();
}
if (a == b) break;
b = b->getReferenceCS();
}
if (a == b && a != nullptr) {
if (a == b && a) {
commonBase = a;
} else {
......@@ -149,4 +77,70 @@ namespace corsika {
return t;
}
inline CoordinateSystemPtr make_translation(CoordinateSystemPtr const& cs,
QuantityVector<length_d> const& vector) {
EigenTransform const translation{EigenTranslation(vector.getEigenVector())};
return std::make_shared<CoordinateSystem const>(CoordinateSystem(cs, translation));
}
template <typename TDim>
inline CoordinateSystemPtr make_rotationToZ(CoordinateSystemPtr const& cs,
Vector<TDim> const& vVec) {
auto const a = vVec.normalized().getComponents(cs).getEigenVector();
auto const a1 = a(0), a2 = a(1), a3 = a(2);
Eigen::Matrix3d A, B;
if (a3 > 0) {
auto const c = 1 / (1 + a3);
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 {
auto const c = 1 / (1 - a3);
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 std::make_shared<CoordinateSystem const>(
CoordinateSystem(cs, EigenTransform(A + B)));
}
template <typename TDim>
inline CoordinateSystemPtr make_rotation(CoordinateSystemPtr const& cs,
QuantityVector<TDim> const& axis,
double const angle) {
if (axis.getEigenVector().isZero()) {
throw std::runtime_error("null-vector given as axis parameter");
}
EigenTransform const rotation{
Eigen::AngleAxisd(angle, axis.getEigenVector().normalized())};
return std::make_shared<CoordinateSystem const>(CoordinateSystem(cs, rotation));
}
template <typename TDim>
inline CoordinateSystemPtr make_translationAndRotation(
CoordinateSystemPtr const& cs, QuantityVector<length_d> const& translation,
QuantityVector<TDim> const& axis, double const angle) {
if (axis.getEigenVector().isZero()) {
throw std::runtime_error("null-vector given as axis parameter");
}
EigenTransform const transf{
Eigen::AngleAxisd(angle, axis.getEigenVector().normalized()) *
EigenTranslation(translation.getEigenVector())};
return std::make_shared<CoordinateSystem const>(CoordinateSystem(cs, transf));
}
} // namespace corsika
/*
* (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.
......@@ -115,4 +113,13 @@ namespace corsika {
return timeLike_ * timeLike_;
}
template <typename TTimeType, typename TSpaceVecType>
inline std::ostream& operator<<(std::ostream& os,
corsika::FourVector<TTimeType, TSpaceVecType> const qv) {
os << '(' << qv.timeLike_ << ", " << qv.spaceLike_ << ") ";
return os;
}
} // namespace corsika
......@@ -11,29 +11,29 @@
#pragma once
#include <corsika/framework/geometry/Point.hpp>
#include <cmath>
#include <corsika/framework/core/PhysicalUnits.hpp>
#include <corsika/framework/geometry/Vector.hpp>
#include <cmath>
namespace corsika {
LengthType Helix::getRadius() const { return radius_; }
Point Helix::getPosition(TimeType t) const {
Point Helix::getPosition(TimeType const t) const {
return r0_ + vPar_ * t +
(vPerp_ * (std::cos(omegaC_ * t) - 1) + uPerp_ * std::sin(omegaC_ * t)) /
omegaC_;
}
Point Helix::getPositionFromArclength(LengthType l) const {
Point Helix::getPositionFromArclength(LengthType const l) const {
return getPosition(getTimeFromArclength(l));
}
LengthType Helix::getArcLength(TimeType t1, TimeType t2) const {
LengthType Helix::getArcLength(TimeType const t1, TimeType const t2) const {
return (vPar_ + vPerp_).getNorm() * (t2 - t1);
}
TimeType Helix::getTimeFromArclength(LengthType l) const {
TimeType Helix::getTimeFromArclength(LengthType const l) const {
return l / (vPar_ + vPerp_).getNorm();
}
......
/*
* (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.
......@@ -16,17 +14,17 @@
namespace corsika {
Point Line::getPosition(TimeType t) const { return start_point_ + velocity_ * t; }
Point Line::getPosition(TimeType const t) const { return start_point_ + velocity_ * t; }
Point Line::getPositionFromArclength(LengthType l) const {
Point Line::getPositionFromArclength(LengthType const l) const {
return start_point_ + velocity_.normalized() * l;
}
LengthType Line::getArcLength(TimeType t1, TimeType t2) const {
LengthType Line::getArcLength(TimeType const t1, TimeType const t2) const {
return velocity_.getNorm() * (t2 - t1);
}
TimeType Line::getTimeFromArclength(LengthType t) const {
TimeType Line::getTimeFromArclength(LengthType const t) const {
return t / velocity_.getNorm();
}
......
......@@ -17,72 +17,82 @@
namespace corsika {
auto Point::getCoordinates() const { return BaseVector<length_d>::getQuantityVector(); }
QuantityVector<length_d> const& Point::getCoordinates() const {
return BaseVector<length_d>::getQuantityVector();
}
QuantityVector<length_d>& Point::getCoordinates() {
return BaseVector<length_d>::getQuantityVector();
}
inline LengthType Point::getX(CoordinateSystemPtr cs) const {
if (*cs == *BaseVector<length_d>::getCoordinateSystem()) {
inline LengthType Point::getX(CoordinateSystemPtr const& pCS) const {
CoordinateSystemPtr const& cs = BaseVector<length_d>::getCoordinateSystem();
if (*pCS == *cs) {
return BaseVector<length_d>::getQuantityVector().getX();
} else {
return QuantityVector<length_d>(
getTransformation(BaseVector<length_d>::getCoordinateSystem(), cs) *
get_transformation(cs, pCS) *
BaseVector<length_d>::getQuantityVector().eigenVector_)
.getX();
}
}
inline LengthType Point::getY(CoordinateSystemPtr cs) const {
if (*cs == *BaseVector<length_d>::getCoordinateSystem()) {
inline LengthType Point::getY(CoordinateSystemPtr const& pCS) const {
CoordinateSystemPtr const& cs = BaseVector<length_d>::getCoordinateSystem();
if (*pCS == *cs) {
return BaseVector<length_d>::getQuantityVector().getY();
} else {
return QuantityVector<length_d>(
getTransformation(BaseVector<length_d>::getCoordinateSystem(), cs) *
get_transformation(cs, pCS) *
BaseVector<length_d>::getQuantityVector().eigenVector_)
.getY();
}
}
inline LengthType Point::getZ(CoordinateSystemPtr cs) const {
if (*cs == *BaseVector<length_d>::getCoordinateSystem()) {
inline LengthType Point::getZ(CoordinateSystemPtr const& pCS) const {
CoordinateSystemPtr const& cs = BaseVector<length_d>::getCoordinateSystem();
if (*pCS == *cs) {
return BaseVector<length_d>::getQuantityVector().getZ();
} else {
return QuantityVector<length_d>(
getTransformation(BaseVector<length_d>::getCoordinateSystem(), cs) *
get_transformation(cs, pCS) *
BaseVector<length_d>::getQuantityVector().eigenVector_)
.getZ();
}
}
/// this always returns a QuantityVector as triple
auto Point::getCoordinates(CoordinateSystemPtr pCS) const {
if (*pCS == *BaseVector<length_d>::getCoordinateSystem()) {
QuantityVector<length_d> Point::getCoordinates(CoordinateSystemPtr const& pCS) const {
CoordinateSystemPtr const& cs = BaseVector<length_d>::getCoordinateSystem();
if (*pCS == *cs) {
return BaseVector<length_d>::getQuantityVector();
} else {
return QuantityVector<length_d>(
getTransformation(BaseVector<length_d>::getCoordinateSystem(), pCS) *
get_transformation(cs, pCS) *
BaseVector<length_d>::getQuantityVector().eigenVector_);
}
}
/*!
* transforms the Point into another CoordinateSystem by changing its
* coordinates interally
*/
void Point::rebase(CoordinateSystemPtr pCS) {
BaseVector<length_d>::quantityVector() = getCoordinates(pCS);
/// this always returns a QuantityVector as triple
QuantityVector<length_d>& Point::getCoordinates(CoordinateSystemPtr const& pCS) {
if (*pCS != *BaseVector<length_d>::getCoordinateSystem()) { rebase(pCS); }
return BaseVector<length_d>::getQuantityVector();
}
void Point::rebase(CoordinateSystemPtr const& pCS) {
BaseVector<length_d>::setQuantityVector(QuantityVector<length_d>(
get_transformation(BaseVector<length_d>::getCoordinateSystem(), pCS) *
BaseVector<length_d>::getQuantityVector().eigenVector_));
BaseVector<length_d>::setCoordinateSystem(pCS);
}
Point Point::operator+(Vector<length_d> const& pVec) const {
return Point(BaseVector<length_d>::getCoordinateSystem(),
getCoordinates() +
pVec.getComponents(BaseVector<length_d>::getCoordinateSystem()));
CoordinateSystemPtr const& cs = BaseVector<length_d>::getCoordinateSystem();
return Point(cs, getCoordinates() + pVec.getComponents(cs));
}
/*!
* returns the distance Vector between two points
*/
Vector<length_d> Point::operator-(Point const& pB) const {
CoordinateSystemPtr cs = BaseVector<length_d>::getCoordinateSystem();
CoordinateSystemPtr const& cs = BaseVector<length_d>::getCoordinateSystem();
return Vector<length_d>(cs, getCoordinates() - pB.getCoordinates(cs));
}
......
......@@ -18,58 +18,63 @@
namespace corsika {
template <typename dim>
inline typename QuantityVector<dim>::quantity_type QuantityVector<dim>::operator[](
size_t index) const {
template <typename TDimension>
inline typename QuantityVector<TDimension>::quantity_type QuantityVector<TDimension>::
operator[](size_t const index) const {
return quantity_type(phys::units::detail::magnitude_tag, eigenVector_[index]);
}
template <typename dim>
inline typename QuantityVector<dim>::quantity_type QuantityVector<dim>::getX() const {
template <typename TDimension>
inline typename QuantityVector<TDimension>::quantity_type
QuantityVector<TDimension>::getX() const {
return (*this)[0];
}
template <typename dim>
inline typename QuantityVector<dim>::quantity_type QuantityVector<dim>::getY() const {
template <typename TDimension>
inline typename QuantityVector<TDimension>::quantity_type
QuantityVector<TDimension>::getY() const {
return (*this)[1];
}
template <typename dim>
inline typename QuantityVector<dim>::quantity_type QuantityVector<dim>::getZ() const {
template <typename TDimension>
inline typename QuantityVector<TDimension>::quantity_type
QuantityVector<TDimension>::getZ() const {
return (*this)[2];
}
template <typename dim>
inline typename QuantityVector<dim>::quantity_type QuantityVector<dim>::getNorm()
const {
template <typename TDimension>
inline typename QuantityVector<TDimension>::quantity_type
QuantityVector<TDimension>::getNorm() const {
return quantity_type(phys::units::detail::magnitude_tag, eigenVector_.norm());
}
template <typename dim>
inline auto QuantityVector<dim>::getSquaredNorm() const {
template <typename TDimension>
inline typename QuantityVector<TDimension>::quantity_square_type
QuantityVector<TDimension>::getSquaredNorm() const {
using QuantitySquared =
decltype(std::declval<quantity_type>() * std::declval<quantity_type>());
return QuantitySquared(phys::units::detail::magnitude_tag,
eigenVector_.squaredNorm());
}
template <typename dim>
inline QuantityVector<dim> QuantityVector<dim>::operator+(
QuantityVector<dim> const& pQVec) const {
return QuantityVector<dim>(eigenVector_ + pQVec.eigenVector_);
template <typename TDimension>
inline QuantityVector<TDimension> QuantityVector<TDimension>::operator+(
QuantityVector<TDimension> const& pQVec) const {
return QuantityVector<TDimension>(eigenVector_ + pQVec.eigenVector_);
}
template <typename dim>
inline QuantityVector<dim> QuantityVector<dim>::operator-(
QuantityVector<dim> const& pQVec) const {
return QuantityVector<dim>(eigenVector_ - pQVec.eigenVector_);
template <typename TDimension>
inline QuantityVector<TDimension> QuantityVector<TDimension>::operator-(
QuantityVector<TDimension> const& pQVec) const {
return QuantityVector<TDimension>(eigenVector_ - pQVec.eigenVector_);
}
template <typename dim>
template <typename ScalarDim>
inline auto QuantityVector<dim>::operator*(
phys::units::quantity<ScalarDim, double> const p) const {
using ResQuantity = phys::units::detail::Product<ScalarDim, dim, double, double>;
template <typename TDimension>
template <typename TScalarDim>
inline auto QuantityVector<TDimension>::operator*(
phys::units::quantity<TScalarDim, double> const p) const {
using ResQuantity =
phys::units::detail::Product<TScalarDim, TDimension, double, double>;
if constexpr (std::is_same<ResQuantity, double>::value) // result dimensionless, not
// a "quantity_type" anymore
......@@ -81,69 +86,73 @@ namespace corsika {
}
}
template <typename dim>
template <typename ScalarDim>
inline auto QuantityVector<dim>::operator/(
phys::units::quantity<ScalarDim, double> const p) const {
template <typename TDimension>
template <typename TScalarDim>
inline auto QuantityVector<TDimension>::operator/(
phys::units::quantity<TScalarDim, double> const p) const {
return (*this) * (1 / p);
}
template <typename dim>
inline auto QuantityVector<dim>::operator*(double const p) const {
return QuantityVector<dim>(eigenVector_ * p);
template <typename TDimension>
inline auto QuantityVector<TDimension>::operator*(double const p) const {
return QuantityVector<TDimension>(eigenVector_ * p);
}
template <typename dim>
inline auto QuantityVector<dim>::operator/(double const p) const {
return QuantityVector<dim>(eigenVector_ / p);
template <typename TDimension>
inline auto QuantityVector<TDimension>::operator/(double const p) const {
return QuantityVector<TDimension>(eigenVector_ / p);
}
template <typename dim>
inline auto& QuantityVector<dim>::operator/=(double const p) {
template <typename TDimension>
inline auto& QuantityVector<TDimension>::operator/=(double const p) {
eigenVector_ /= p;
return *this;
}
template <typename dim>
inline auto& QuantityVector<dim>::operator*=(double const p) {
template <typename TDimension>
inline auto& QuantityVector<TDimension>::operator*=(double const p) {
eigenVector_ *= p;
return *this;
}
template <typename dim>
inline auto& QuantityVector<dim>::operator+=(QuantityVector<dim> const& pQVec) {
template <typename TDimension>
inline auto& QuantityVector<TDimension>::operator+=(
QuantityVector<TDimension> const& pQVec) {
eigenVector_ += pQVec.eigenVector_;
return *this;
}
template <typename dim>
inline auto& QuantityVector<dim>::operator-=(QuantityVector<dim> const& pQVec) {
template <typename TDimension>
inline auto& QuantityVector<TDimension>::operator-=(
QuantityVector<TDimension> const& pQVec) {
eigenVector_ -= pQVec.eigenVector_;
return *this;
}
template <typename dim>
inline auto& QuantityVector<dim>::operator-() const {
return QuantityVector<dim>(-eigenVector_);
template <typename TDimension>
inline auto& QuantityVector<TDimension>::operator-() const {
return QuantityVector<TDimension>(-eigenVector_);
}
template <typename dim>
inline auto QuantityVector<dim>::normalized() const {
return QuantityVector<dim>(eigenVector_.normalized());
template <typename TDimension>
inline auto QuantityVector<TDimension>::normalized() const {
return QuantityVector<TDimension>(eigenVector_.normalized());
}
template <typename dim>
inline auto QuantityVector<dim>::operator==(QuantityVector<dim> const& p) const {
template <typename TDimension>
inline auto QuantityVector<TDimension>::operator==(
QuantityVector<TDimension> const& p) const {
return eigenVector_ == p.eigenVector_;
}
template <typename dim>
inline std::ostream& operator<<(std::ostream& os, corsika::QuantityVector<dim> qv) {
using quantity_type = phys::units::quantity<dim, double>;
template <typename TDimension>
inline std::ostream& operator<<(std::ostream& os,
corsika::QuantityVector<TDimension> const qv) {
using quantity_type = phys::units::quantity<TDimension, double>;
os << '(' << qv.eigenVector_(0) << ' ' << qv.eigenVector_(1) << ' '
<< qv.eigenVector_(2) << ") "
<< phys::units::to_unit_symbol<dim, double>(
<< phys::units::to_unit_symbol<TDimension, double>(
quantity_type(phys::units::detail::magnitude_tag, 1));
return os;
}
......
......@@ -17,7 +17,7 @@
namespace corsika {
template <typename TType>
Point Trajectory<TType>::getPosition(double u) const {
Point Trajectory<TType>::getPosition(double const u) const {
return TType::getPosition(timeLength_ * u);
}
......@@ -32,14 +32,14 @@ namespace corsika {
}
template <typename TType>
LengthType Trajectory<TType>::getDistance(TimeType t) const {
LengthType Trajectory<TType>::getDistance(TimeType const t) const {
assert(t <= timeLength_);
assert(t >= 0 * second);
return TType::getArcLength(0 * second, t);
}
template <typename TType>
void Trajectory<TType>::getLimitEndTo(LengthType limit) {
void Trajectory<TType>::getLimitEndTo(LengthType const limit) {
timeLength_ = TType::getTimeFromArclength(limit);
}
......
......@@ -26,7 +26,7 @@ namespace corsika {
COMBoost::COMBoost(FourVector<HEPEnergyType, Vector<hepmomentum_d>> const& Pprojectile,
HEPMassType const massTarget)
: originalCS_{Pprojectile.getSpaceLikeComponents().getCoordinateSystem()}
, rotatedCS_{originalCS_->rotateToZ(Pprojectile.getSpaceLikeComponents())} {
, rotatedCS_{make_rotationToZ(originalCS_, Pprojectile.getSpaceLikeComponents())} {
auto const pProjectile = Pprojectile.getSpaceLikeComponents();
auto const pProjNormSquared = pProjectile.getSquaredNorm();
auto const pProjNorm = sqrt(pProjNormSquared);
......@@ -48,7 +48,7 @@ namespace corsika {
COMBoost::COMBoost(Vector<hepmomentum_d> const& momentum, HEPEnergyType mass)
: originalCS_{momentum.getCoordinateSystem()}
, rotatedCS_{originalCS_->rotateToZ(momentum)} {
, rotatedCS_{make_rotationToZ(originalCS_, momentum)} {
auto const squaredNorm = momentum.getSquaredNorm();
auto const norm = sqrt(squaredNorm);
auto const sinhEta = -norm / mass;
......@@ -89,7 +89,7 @@ namespace corsika {
auto const boostedZ = inverseBoost_ * com;
auto const E_lab = boostedZ(0) * 1_GeV;
pCM.eigenVector()(2) = boostedZ(1) * (1_GeV).magnitude();
pCM.getEigenVector()[2] = boostedZ(1) * (1_GeV).magnitude();