diff --git a/Common/Constants/include/CommonConstants/LHCConstants.h b/Common/Constants/include/CommonConstants/LHCConstants.h
index 1582f2166ca0f..2cdd0b5118731 100644
--- a/Common/Constants/include/CommonConstants/LHCConstants.h
+++ b/Common/Constants/include/CommonConstants/LHCConstants.h
@@ -16,6 +16,8 @@
#ifndef ALICEO2_LHCCONSTANTS_H_
#define ALICEO2_LHCCONSTANTS_H_
+#include "GPUCommonDouble.h"
+
#include "GPUCommonDef.h"
namespace o2
@@ -31,12 +33,12 @@ enum BeamDirection : int { BeamA, // beamA = beam 0,
InteractingBC = -1 // as used in the BunchFilling class
};
GPUglobalconstexpr() int LHCMaxBunches = 3564; // max N bunches
-constexpr double LHCRFFreq = 400.789e6; // LHC RF frequency in Hz
-constexpr double LHCBunchSpacingNS = 10 * 1.e9 / LHCRFFreq; // bunch spacing in ns (10 RFbuckets)
-constexpr double LHCOrbitNS = LHCMaxBunches * LHCBunchSpacingNS; // orbit duration in ns
-constexpr double LHCRevFreq = 1.e9 / LHCOrbitNS; // revolution frequency
-constexpr double LHCBunchSpacingMUS = LHCBunchSpacingNS * 1e-3; // bunch spacing in \mus (10 RFbuckets)
-constexpr double LHCOrbitMUS = LHCOrbitNS * 1e-3; // orbit duration in \mus
+GPUglobalconstexpr() o2::gpu::GPUdoubleValue LHCRFFreq = 400.789e6; // LHC RF frequency in Hz
+GPUglobalconstexpr() o2::gpu::GPUdoubleValue LHCBunchSpacingNS = 10 * 1.e9 / LHCRFFreq; // bunch spacing in ns (10 RFbuckets)
+GPUglobalconstexpr() o2::gpu::GPUdoubleValue LHCOrbitNS = LHCMaxBunches * LHCBunchSpacingNS; // orbit duration in ns
+GPUglobalconstexpr() o2::gpu::GPUdoubleValue LHCRevFreq = 1.e9 / LHCOrbitNS; // revolution frequency
+GPUglobalconstexpr() o2::gpu::GPUdoubleValue LHCBunchSpacingMUS = LHCBunchSpacingNS * 1e-3; // bunch spacing in \mus (10 RFbuckets)
+GPUglobalconstexpr() o2::gpu::GPUdoubleValue LHCOrbitMUS = LHCOrbitNS * 1e-3; // orbit duration in \mus
GPUglobalconstexpr() unsigned int MaxNOrbits = 0xffffffff;
// Offsets of A, C beam bunches at P2
diff --git a/Common/Constants/include/CommonConstants/PhysicsConstants.h b/Common/Constants/include/CommonConstants/PhysicsConstants.h
index 051fb6d2a6e89..86e5162e2b654 100644
--- a/Common/Constants/include/CommonConstants/PhysicsConstants.h
+++ b/Common/Constants/include/CommonConstants/PhysicsConstants.h
@@ -18,8 +18,19 @@
#ifndef ALICEO2_PHYSICSCONSTANTS_H_
#define ALICEO2_PHYSICSCONSTANTS_H_
+#include "GPUCommonDef.h"
+
namespace o2::constants::physics
{
+#ifdef __METAL__
+// MSL has no double. The masses only ever serve as compile-time initialisers,
+// and the tables built from them (PID::sMasses) are float already, so nothing
+// is lost. Every other backend, device included, keeps double.
+using MassType = float;
+#else
+using MassType = double;
+#endif
+
// particles masses
// BEGINNING OF THE GENERATED BLOCK.
@@ -102,153 +113,153 @@ enum Pdg {
};
/// \brief Declarations of masses for additional particles
-constexpr double MassEta = 0.547862;
-constexpr double MassOmega = 0.78266;
-constexpr double MassEtaPrime = 0.95778;
-constexpr double MassB0 = 5.27966;
-constexpr double MassB0Bar = 5.27966;
-constexpr double MassBPlus = 5.27934;
-constexpr double MassBCPlus = 6.27447;
-constexpr double MassBS = 5.36692;
-constexpr double MassBSBar = 5.36692;
-constexpr double MassD0 = 1.86484;
-constexpr double MassD0Bar = 1.86484;
-constexpr double MassD0StarPlus = 2.272;
-constexpr double MassD0Star0 = 2.343;
-constexpr double MassD1Plus = 2.372;
-constexpr double MassD10 = 2.412;
-constexpr double MassD2StarPlus = 2.4601;
-constexpr double MassD2Star0 = 2.4611;
-constexpr double MassDMinus = 1.86966;
-constexpr double MassDPlus = 1.86966;
-constexpr double MassDS = 1.96835;
-constexpr double MassDSBar = 1.96835;
-constexpr double MassDSStar = 2.1122;
-constexpr double MassDS1 = 2.53511;
-constexpr double MassDS1Star2700 = 2.714;
-constexpr double MassDS1Star2860 = 2.859;
-constexpr double MassDS2Star = 2.5691;
-constexpr double MassDS3Star2860 = 2.86;
-constexpr double MassDStar = 2.01026;
-constexpr double MassDStar0 = 2.00685;
-constexpr double MassChiC1 = 3.51067;
-constexpr double MassJPsi = 3.0969;
-constexpr double MassLambdaB0 = 5.6196;
-constexpr double MassLambdaCPlus = 2.28646;
-constexpr double MassLambdaCPlus2860 = 2.8561;
-constexpr double MassLambdaCPlus2880 = 2.8816;
-constexpr double MassLambdaCPlus2940 = 2.9396;
-constexpr double MassOmegaC0 = 2.6952;
-constexpr double MassK0Star892 = 0.89555;
-constexpr double MassKPlusStar892 = 0.89167;
-constexpr double MassPhi = 1.019461;
-constexpr double MassSigmaC0 = 2.45375;
-constexpr double MassSigmaCPlusPlus = 2.45397;
-constexpr double MassSigmaCStar0 = 2.51848;
-constexpr double MassSigmaCStarPlusPlus = 2.51841;
-constexpr double MassX3872 = 3.87165;
-constexpr double MassXi0 = 1.31486;
-constexpr double MassXiB0 = 5.7919;
-constexpr double MassXiCCPlusPlus = 3.62155;
-constexpr double MassXiCPlus = 2.46771;
-constexpr double MassXiC0 = 2.47044;
-constexpr double MassXiC3055Plus = 3.0559;
-constexpr double MassXiC3080Plus = 3.0772;
-constexpr double MassXiC3055_0 = 3.059;
-constexpr double MassXiC3080_0 = 3.0799;
-constexpr double MassDeuteron = 1.87561294257;
-constexpr double MassTriton = 2.80892113298;
-constexpr double MassHelium3 = 2.80839160743;
-constexpr double MassAlpha = 3.7273794066;
-constexpr double MassLithium4 = 3.7513;
-constexpr double MassHyperTriton = 2.991134;
-constexpr double MassHyperHydrogen4 = 3.922434;
-constexpr double MassHyperHelium4 = 3.921728;
-constexpr double MassHyperHelium5 = 4.839961;
-constexpr double MassHyperHelium4Sigma = 3.995;
-constexpr double MassLambda1520_Py = 1.5195;
-constexpr double MassK1_1270_0 = 1.253;
-constexpr double MassK1_1270Plus = 1.272;
-constexpr double MassCDeuteron = 3.226;
+GPUglobalconstexpr() MassType MassEta = 0.547862;
+GPUglobalconstexpr() MassType MassOmega = 0.78266;
+GPUglobalconstexpr() MassType MassEtaPrime = 0.95778;
+GPUglobalconstexpr() MassType MassB0 = 5.27966;
+GPUglobalconstexpr() MassType MassB0Bar = 5.27966;
+GPUglobalconstexpr() MassType MassBPlus = 5.27934;
+GPUglobalconstexpr() MassType MassBCPlus = 6.27447;
+GPUglobalconstexpr() MassType MassBS = 5.36692;
+GPUglobalconstexpr() MassType MassBSBar = 5.36692;
+GPUglobalconstexpr() MassType MassD0 = 1.86484;
+GPUglobalconstexpr() MassType MassD0Bar = 1.86484;
+GPUglobalconstexpr() MassType MassD0StarPlus = 2.272;
+GPUglobalconstexpr() MassType MassD0Star0 = 2.343;
+GPUglobalconstexpr() MassType MassD1Plus = 2.372;
+GPUglobalconstexpr() MassType MassD10 = 2.412;
+GPUglobalconstexpr() MassType MassD2StarPlus = 2.4601;
+GPUglobalconstexpr() MassType MassD2Star0 = 2.4611;
+GPUglobalconstexpr() MassType MassDMinus = 1.86966;
+GPUglobalconstexpr() MassType MassDPlus = 1.86966;
+GPUglobalconstexpr() MassType MassDS = 1.96835;
+GPUglobalconstexpr() MassType MassDSBar = 1.96835;
+GPUglobalconstexpr() MassType MassDSStar = 2.1122;
+GPUglobalconstexpr() MassType MassDS1 = 2.53511;
+GPUglobalconstexpr() MassType MassDS1Star2700 = 2.714;
+GPUglobalconstexpr() MassType MassDS1Star2860 = 2.859;
+GPUglobalconstexpr() MassType MassDS2Star = 2.5691;
+GPUglobalconstexpr() MassType MassDS3Star2860 = 2.86;
+GPUglobalconstexpr() MassType MassDStar = 2.01026;
+GPUglobalconstexpr() MassType MassDStar0 = 2.00685;
+GPUglobalconstexpr() MassType MassChiC1 = 3.51067;
+GPUglobalconstexpr() MassType MassJPsi = 3.0969;
+GPUglobalconstexpr() MassType MassLambdaB0 = 5.6196;
+GPUglobalconstexpr() MassType MassLambdaCPlus = 2.28646;
+GPUglobalconstexpr() MassType MassLambdaCPlus2860 = 2.8561;
+GPUglobalconstexpr() MassType MassLambdaCPlus2880 = 2.8816;
+GPUglobalconstexpr() MassType MassLambdaCPlus2940 = 2.9396;
+GPUglobalconstexpr() MassType MassOmegaC0 = 2.6952;
+GPUglobalconstexpr() MassType MassK0Star892 = 0.89555;
+GPUglobalconstexpr() MassType MassKPlusStar892 = 0.89167;
+GPUglobalconstexpr() MassType MassPhi = 1.019461;
+GPUglobalconstexpr() MassType MassSigmaC0 = 2.45375;
+GPUglobalconstexpr() MassType MassSigmaCPlusPlus = 2.45397;
+GPUglobalconstexpr() MassType MassSigmaCStar0 = 2.51848;
+GPUglobalconstexpr() MassType MassSigmaCStarPlusPlus = 2.51841;
+GPUglobalconstexpr() MassType MassX3872 = 3.87165;
+GPUglobalconstexpr() MassType MassXi0 = 1.31486;
+GPUglobalconstexpr() MassType MassXiB0 = 5.7919;
+GPUglobalconstexpr() MassType MassXiCCPlusPlus = 3.62155;
+GPUglobalconstexpr() MassType MassXiCPlus = 2.46771;
+GPUglobalconstexpr() MassType MassXiC0 = 2.47044;
+GPUglobalconstexpr() MassType MassXiC3055Plus = 3.0559;
+GPUglobalconstexpr() MassType MassXiC3080Plus = 3.0772;
+GPUglobalconstexpr() MassType MassXiC3055_0 = 3.059;
+GPUglobalconstexpr() MassType MassXiC3080_0 = 3.0799;
+GPUglobalconstexpr() MassType MassDeuteron = 1.87561294257;
+GPUglobalconstexpr() MassType MassTriton = 2.80892113298;
+GPUglobalconstexpr() MassType MassHelium3 = 2.80839160743;
+GPUglobalconstexpr() MassType MassAlpha = 3.7273794066;
+GPUglobalconstexpr() MassType MassLithium4 = 3.7513;
+GPUglobalconstexpr() MassType MassHyperTriton = 2.991134;
+GPUglobalconstexpr() MassType MassHyperHydrogen4 = 3.922434;
+GPUglobalconstexpr() MassType MassHyperHelium4 = 3.921728;
+GPUglobalconstexpr() MassType MassHyperHelium5 = 4.839961;
+GPUglobalconstexpr() MassType MassHyperHelium4Sigma = 3.995;
+GPUglobalconstexpr() MassType MassLambda1520_Py = 1.5195;
+GPUglobalconstexpr() MassType MassK1_1270_0 = 1.253;
+GPUglobalconstexpr() MassType MassK1_1270Plus = 1.272;
+GPUglobalconstexpr() MassType MassCDeuteron = 3.226;
/// \brief Declarations of masses for particles in ROOT PDG_t
-constexpr double MassDown = 0.00467;
-constexpr double MassDownBar = 0.00467;
-constexpr double MassUp = 0.00216;
-constexpr double MassUpBar = 0.00216;
-constexpr double MassStrange = 0.0934;
-constexpr double MassStrangeBar = 0.0934;
-constexpr double MassCharm = 1.27;
-constexpr double MassCharmBar = 1.27;
-constexpr double MassBottom = 4.18;
-constexpr double MassBottomBar = 4.18;
-constexpr double MassTop = 172.5;
-constexpr double MassTopBar = 172.5;
-constexpr double MassGluon = 0.0;
-constexpr double MassElectron = 0.000510999;
-constexpr double MassPositron = 0.000510999;
-constexpr double MassNuE = 0.0;
-constexpr double MassNuEBar = 0.0;
-constexpr double MassMuonMinus = 0.1056584;
-constexpr double MassMuonPlus = 0.1056584;
-constexpr double MassNuMu = 0.0;
-constexpr double MassNuMuBar = 0.0;
-constexpr double MassTauMinus = 1.77686;
-constexpr double MassTauPlus = 1.77686;
-constexpr double MassNuTau = 0.0;
-constexpr double MassNuTauBar = 0.0;
-constexpr double MassGamma = 0.0;
-constexpr double MassZ0 = 91.1876;
-constexpr double MassWPlus = 80.377;
-constexpr double MassWMinus = 80.377;
-constexpr double MassPi0 = 0.1349768;
-constexpr double MassK0Long = 0.497611;
-constexpr double MassPiPlus = 0.1395704;
-constexpr double MassPiMinus = 0.1395704;
-constexpr double MassProton = 0.9382721;
-constexpr double MassProtonBar = 0.9382721;
-constexpr double MassNeutron = 0.9395654;
-constexpr double MassNeutronBar = 0.9395654;
-constexpr double MassK0Short = 0.497611;
-constexpr double MassK0 = 0.497611;
-constexpr double MassK0Bar = 0.497611;
-constexpr double MassKPlus = 0.493677;
-constexpr double MassKMinus = 0.493677;
-constexpr double MassLambda0 = 1.115683;
-constexpr double MassLambda0Bar = 1.115683;
-constexpr double MassLambda1520 = 1.519;
-constexpr double MassSigmaMinus = 1.197449;
-constexpr double MassSigmaBarPlus = 1.197449;
-constexpr double MassSigmaPlus = 1.18937;
-constexpr double MassSigmaBarMinus = 1.18937;
-constexpr double MassSigma0 = 1.192642;
-constexpr double MassSigma0Bar = 1.192642;
-constexpr double MassXiMinus = 1.32171;
-constexpr double MassXiPlusBar = 1.32171;
-constexpr double MassOmegaMinus = 1.67245;
-constexpr double MassOmegaPlusBar = 1.67245;
+GPUglobalconstexpr() MassType MassDown = 0.00467;
+GPUglobalconstexpr() MassType MassDownBar = 0.00467;
+GPUglobalconstexpr() MassType MassUp = 0.00216;
+GPUglobalconstexpr() MassType MassUpBar = 0.00216;
+GPUglobalconstexpr() MassType MassStrange = 0.0934;
+GPUglobalconstexpr() MassType MassStrangeBar = 0.0934;
+GPUglobalconstexpr() MassType MassCharm = 1.27;
+GPUglobalconstexpr() MassType MassCharmBar = 1.27;
+GPUglobalconstexpr() MassType MassBottom = 4.18;
+GPUglobalconstexpr() MassType MassBottomBar = 4.18;
+GPUglobalconstexpr() MassType MassTop = 172.5;
+GPUglobalconstexpr() MassType MassTopBar = 172.5;
+GPUglobalconstexpr() MassType MassGluon = 0.0;
+GPUglobalconstexpr() MassType MassElectron = 0.000510999;
+GPUglobalconstexpr() MassType MassPositron = 0.000510999;
+GPUglobalconstexpr() MassType MassNuE = 0.0;
+GPUglobalconstexpr() MassType MassNuEBar = 0.0;
+GPUglobalconstexpr() MassType MassMuonMinus = 0.1056584;
+GPUglobalconstexpr() MassType MassMuonPlus = 0.1056584;
+GPUglobalconstexpr() MassType MassNuMu = 0.0;
+GPUglobalconstexpr() MassType MassNuMuBar = 0.0;
+GPUglobalconstexpr() MassType MassTauMinus = 1.77686;
+GPUglobalconstexpr() MassType MassTauPlus = 1.77686;
+GPUglobalconstexpr() MassType MassNuTau = 0.0;
+GPUglobalconstexpr() MassType MassNuTauBar = 0.0;
+GPUglobalconstexpr() MassType MassGamma = 0.0;
+GPUglobalconstexpr() MassType MassZ0 = 91.1876;
+GPUglobalconstexpr() MassType MassWPlus = 80.377;
+GPUglobalconstexpr() MassType MassWMinus = 80.377;
+GPUglobalconstexpr() MassType MassPi0 = 0.1349768;
+GPUglobalconstexpr() MassType MassK0Long = 0.497611;
+GPUglobalconstexpr() MassType MassPiPlus = 0.1395704;
+GPUglobalconstexpr() MassType MassPiMinus = 0.1395704;
+GPUglobalconstexpr() MassType MassProton = 0.9382721;
+GPUglobalconstexpr() MassType MassProtonBar = 0.9382721;
+GPUglobalconstexpr() MassType MassNeutron = 0.9395654;
+GPUglobalconstexpr() MassType MassNeutronBar = 0.9395654;
+GPUglobalconstexpr() MassType MassK0Short = 0.497611;
+GPUglobalconstexpr() MassType MassK0 = 0.497611;
+GPUglobalconstexpr() MassType MassK0Bar = 0.497611;
+GPUglobalconstexpr() MassType MassKPlus = 0.493677;
+GPUglobalconstexpr() MassType MassKMinus = 0.493677;
+GPUglobalconstexpr() MassType MassLambda0 = 1.115683;
+GPUglobalconstexpr() MassType MassLambda0Bar = 1.115683;
+GPUglobalconstexpr() MassType MassLambda1520 = 1.519;
+GPUglobalconstexpr() MassType MassSigmaMinus = 1.197449;
+GPUglobalconstexpr() MassType MassSigmaBarPlus = 1.197449;
+GPUglobalconstexpr() MassType MassSigmaPlus = 1.18937;
+GPUglobalconstexpr() MassType MassSigmaBarMinus = 1.18937;
+GPUglobalconstexpr() MassType MassSigma0 = 1.192642;
+GPUglobalconstexpr() MassType MassSigma0Bar = 1.192642;
+GPUglobalconstexpr() MassType MassXiMinus = 1.32171;
+GPUglobalconstexpr() MassType MassXiPlusBar = 1.32171;
+GPUglobalconstexpr() MassType MassOmegaMinus = 1.67245;
+GPUglobalconstexpr() MassType MassOmegaPlusBar = 1.67245;
// END OF THE GENERATED BLOCK
// legacy names
-constexpr double MassPhoton = MassGamma;
-constexpr double MassMuon = MassMuonMinus;
-constexpr double MassPionCharged = MassPiPlus;
-constexpr double MassPionNeutral = MassPi0;
-constexpr double MassKaonCharged = MassKPlus;
-constexpr double MassKaonNeutral = MassK0;
-constexpr double MassLambda = MassLambda0;
-constexpr double MassHyperhydrog4 = MassHyperHydrogen4;
-constexpr double MassHyperhelium4 = MassHyperHelium4;
-constexpr double MassHyperhelium4sigma = MassHyperHelium4Sigma;
+GPUglobalconstexpr() MassType MassPhoton = MassGamma;
+GPUglobalconstexpr() MassType MassMuon = MassMuonMinus;
+GPUglobalconstexpr() MassType MassPionCharged = MassPiPlus;
+GPUglobalconstexpr() MassType MassPionNeutral = MassPi0;
+GPUglobalconstexpr() MassType MassKaonCharged = MassKPlus;
+GPUglobalconstexpr() MassType MassKaonNeutral = MassK0;
+GPUglobalconstexpr() MassType MassLambda = MassLambda0;
+GPUglobalconstexpr() MassType MassHyperhydrog4 = MassHyperHydrogen4;
+GPUglobalconstexpr() MassType MassHyperhelium4 = MassHyperHelium4;
+GPUglobalconstexpr() MassType MassHyperhelium4sigma = MassHyperHelium4Sigma;
// Light speed
-constexpr float LightSpeedCm2S = 299792458.e2; // C in cm/s
-constexpr float LightSpeedCm2NS = LightSpeedCm2S * 1e-9; // C in cm/ns
-constexpr float LightSpeedCm2PS = LightSpeedCm2S * 1e-12; // C in cm/ps
+GPUglobalconstexpr() float LightSpeedCm2S = 299792458.e2; // C in cm/s
+GPUglobalconstexpr() float LightSpeedCm2NS = LightSpeedCm2S * 1e-9; // C in cm/ns
+GPUglobalconstexpr() float LightSpeedCm2PS = LightSpeedCm2S * 1e-12; // C in cm/ps
// Light speed inverse
-constexpr float invLightSpeedCm2PS = 1. / LightSpeedCm2PS; // 1/C in ps/cm
+GPUglobalconstexpr() float invLightSpeedCm2PS = 1. / LightSpeedCm2PS; // 1/C in ps/cm
} // namespace o2::constants::physics
diff --git a/Common/Constants/include/CommonConstants/make_pdg_header.py b/Common/Constants/include/CommonConstants/make_pdg_header.py
index b2dac688fd098..512f3f9223a52 100755
--- a/Common/Constants/include/CommonConstants/make_pdg_header.py
+++ b/Common/Constants/include/CommonConstants/make_pdg_header.py
@@ -169,9 +169,9 @@ def mass(code):
return dbPdg.Mass(code, success)
-def declare_mass(pdg, mass_type="double") -> str:
+def declare_mass(pdg, mass_type="MassType") -> str:
"""Returns a C++ declaration of a particle mass constant."""
- return f"constexpr {mass_type} Mass{pdg.name[1:]} = {mass(pdg.value)};"
+ return f"GPUglobalconstexpr() {mass_type} Mass{pdg.name[1:]} = {mass(pdg.value)};"
def main():
diff --git a/Common/ML/include/ML/3rdparty/GPUORTFloat16.h b/Common/ML/include/ML/3rdparty/GPUORTFloat16.h
index 75e146d872cd1..38ee5f6f7b5ba 100644
--- a/Common/ML/include/ML/3rdparty/GPUORTFloat16.h
+++ b/Common/ML/include/ML/3rdparty/GPUORTFloat16.h
@@ -58,19 +58,19 @@ struct Float16Impl {
///
///
///
- GPUd() constexpr static uint16_t ToUint16Impl(float v) noexcept;
+ GPUd() constexpr static uint16_t ToUint16Impl(float v) GPUnoexcept();
///
/// Converts float16 to float
///
/// float representation of float16 value
- GPUd() float ToFloatImpl() const noexcept;
+ GPUd() float ToFloatImpl() const GPUnoexcept();
///
/// Creates an instance that represents absolute value.
///
/// Absolute value
- GPUd() uint16_t AbsImpl() const noexcept
+ GPUd() uint16_t AbsImpl() const GPUnoexcept()
{
return static_cast(val & ~kSignMask);
}
@@ -79,24 +79,24 @@ struct Float16Impl {
/// Creates a new instance with the sign flipped.
///
/// Flipped sign instance
- GPUd() uint16_t NegateImpl() const noexcept
+ GPUd() uint16_t NegateImpl() const GPUnoexcept()
{
return IsNaN() ? val : static_cast(val ^ kSignMask);
}
public:
// uint16_t special values
- static constexpr uint16_t kSignMask = 0x8000U;
- static constexpr uint16_t kBiasedExponentMask = 0x7C00U;
- static constexpr uint16_t kPositiveInfinityBits = 0x7C00U;
- static constexpr uint16_t kNegativeInfinityBits = 0xFC00U;
- static constexpr uint16_t kPositiveQNaNBits = 0x7E00U;
- static constexpr uint16_t kNegativeQNaNBits = 0xFE00U;
- static constexpr uint16_t kEpsilonBits = 0x4170U;
- static constexpr uint16_t kMinValueBits = 0xFBFFU; // Minimum normal number
- static constexpr uint16_t kMaxValueBits = 0x7BFFU; // Largest normal number
- static constexpr uint16_t kOneBits = 0x3C00U;
- static constexpr uint16_t kMinusOneBits = 0xBC00U;
+ static GPUglobalconstexpr() uint16_t kSignMask = 0x8000U;
+ static GPUglobalconstexpr() uint16_t kBiasedExponentMask = 0x7C00U;
+ static GPUglobalconstexpr() uint16_t kPositiveInfinityBits = 0x7C00U;
+ static GPUglobalconstexpr() uint16_t kNegativeInfinityBits = 0xFC00U;
+ static GPUglobalconstexpr() uint16_t kPositiveQNaNBits = 0x7E00U;
+ static GPUglobalconstexpr() uint16_t kNegativeQNaNBits = 0xFE00U;
+ static GPUglobalconstexpr() uint16_t kEpsilonBits = 0x4170U;
+ static GPUglobalconstexpr() uint16_t kMinValueBits = 0xFBFFU; // Minimum normal number
+ static GPUglobalconstexpr() uint16_t kMaxValueBits = 0x7BFFU; // Largest normal number
+ static GPUglobalconstexpr() uint16_t kOneBits = 0x3C00U;
+ static GPUglobalconstexpr() uint16_t kMinusOneBits = 0xBC00U;
uint16_t val{0};
@@ -106,7 +106,7 @@ struct Float16Impl {
/// Checks if the value is negative
///
/// true if negative
- GPUd() bool IsNegative() const noexcept
+ GPUd() bool IsNegative() const GPUnoexcept()
{
return static_cast(val) < 0;
}
@@ -115,7 +115,7 @@ struct Float16Impl {
/// Tests if the value is NaN
///
/// true if NaN
- GPUd() bool IsNaN() const noexcept
+ GPUd() bool IsNaN() const GPUnoexcept()
{
return AbsImpl() > kPositiveInfinityBits;
}
@@ -124,7 +124,7 @@ struct Float16Impl {
/// Tests if the value is finite
///
/// true if finite
- GPUd() bool IsFinite() const noexcept
+ GPUd() bool IsFinite() const GPUnoexcept()
{
return AbsImpl() < kPositiveInfinityBits;
}
@@ -133,7 +133,7 @@ struct Float16Impl {
/// Tests if the value represents positive infinity.
///
/// true if positive infinity
- GPUd() bool IsPositiveInfinity() const noexcept
+ GPUd() bool IsPositiveInfinity() const GPUnoexcept()
{
return val == kPositiveInfinityBits;
}
@@ -142,7 +142,7 @@ struct Float16Impl {
/// Tests if the value represents negative infinity
///
/// true if negative infinity
- GPUd() bool IsNegativeInfinity() const noexcept
+ GPUd() bool IsNegativeInfinity() const GPUnoexcept()
{
return val == kNegativeInfinityBits;
}
@@ -151,7 +151,7 @@ struct Float16Impl {
/// Tests if the value is either positive or negative infinity.
///
/// True if absolute value is infinity
- GPUd() bool IsInfinity() const noexcept
+ GPUd() bool IsInfinity() const GPUnoexcept()
{
return AbsImpl() == kPositiveInfinityBits;
}
@@ -160,7 +160,7 @@ struct Float16Impl {
/// Tests if the value is NaN or zero. Useful for comparisons.
///
/// True if NaN or zero.
- GPUd() bool IsNaNOrZero() const noexcept
+ GPUd() bool IsNaNOrZero() const GPUnoexcept()
{
auto abs = AbsImpl();
return (abs == 0 || abs > kPositiveInfinityBits);
@@ -170,7 +170,7 @@ struct Float16Impl {
/// Tests if the value is normal (not zero, subnormal, infinite, or NaN).
///
/// True if so
- GPUd() bool IsNormal() const noexcept
+ GPUd() bool IsNormal() const GPUnoexcept()
{
auto abs = AbsImpl();
return (abs < kPositiveInfinityBits) // is finite
@@ -182,7 +182,7 @@ struct Float16Impl {
/// Tests if the value is subnormal (denormal).
///
/// True if so
- GPUd() bool IsSubnormal() const noexcept
+ GPUd() bool IsSubnormal() const GPUnoexcept()
{
auto abs = AbsImpl();
return (abs < kPositiveInfinityBits) // is finite
@@ -194,13 +194,13 @@ struct Float16Impl {
/// Creates an instance that represents absolute value.
///
/// Absolute value
- GPUd() Derived Abs() const noexcept { return Derived::FromBits(AbsImpl()); }
+ GPUd() Derived Abs() const GPUnoexcept() { return Derived::FromBits(AbsImpl()); }
///
/// Creates a new instance with the sign flipped.
///
/// Flipped sign instance
- GPUd() Derived Negate() const noexcept { return Derived::FromBits(NegateImpl()); }
+ GPUd() Derived Negate() const GPUnoexcept() { return Derived::FromBits(NegateImpl()); }
///
/// IEEE defines that positive and negative zero are equal, this gives us a quick equality check
@@ -210,12 +210,12 @@ struct Float16Impl {
/// first value
/// second value
/// True if both arguments represent zero
- GPUd() static bool AreZero(const Float16Impl& lhs, const Float16Impl& rhs) noexcept
+ GPUd() static bool AreZero(const Float16Impl& lhs, const Float16Impl& rhs) GPUnoexcept()
{
return static_cast((lhs.val | rhs.val) & ~kSignMask) == 0;
}
- GPUd() bool operator==(const Float16Impl& rhs) const noexcept
+ GPUd() bool operator==(const Float16Impl& rhs) const GPUnoexcept()
{
if (IsNaN() || rhs.IsNaN()) {
// IEEE defines that NaN is not equal to anything, including itself.
@@ -224,9 +224,9 @@ struct Float16Impl {
return val == rhs.val;
}
- GPUd() bool operator!=(const Float16Impl& rhs) const noexcept { return !(*this == rhs); }
+ GPUd() bool operator!=(const Float16Impl& rhs) const GPUnoexcept() { return !(*this == rhs); }
- GPUd() bool operator<(const Float16Impl& rhs) const noexcept
+ GPUd() bool operator<(const Float16Impl& rhs) const GPUnoexcept()
{
if (IsNaN() || rhs.IsNaN()) {
// IEEE defines that NaN is unordered with respect to everything, including itself.
@@ -275,7 +275,7 @@ union float32_bits {
}; // namespace detail
template
-GPUdi() constexpr uint16_t Float16Impl::ToUint16Impl(float v) noexcept
+GPUdi() constexpr uint16_t Float16Impl::ToUint16Impl(float v) GPUnoexcept()
{
detail::float32_bits f{};
f.f = v;
@@ -324,7 +324,7 @@ GPUdi() constexpr uint16_t Float16Impl::ToUint16Impl(float v) noexcept
}
template
-GPUdi() float Float16Impl::ToFloatImpl() const noexcept
+GPUdi() float Float16Impl::ToFloatImpl() const GPUnoexcept()
{
constexpr detail::float32_bits magic = {113 << 23};
constexpr unsigned int shifted_exp = 0x7c00 << 13; // exponent mask after shift
@@ -364,19 +364,19 @@ struct BFloat16Impl {
///
///
///
- GPUd() static uint16_t ToUint16Impl(float v) noexcept;
+ GPUd() static uint16_t ToUint16Impl(float v) GPUnoexcept();
///
/// Converts bfloat16 to float
///
/// float representation of bfloat16 value
- GPUd() float ToFloatImpl() const noexcept;
+ GPUd() float ToFloatImpl() const GPUnoexcept();
///
/// Creates an instance that represents absolute value.
///
/// Absolute value
- GPUd() uint16_t AbsImpl() const noexcept
+ GPUd() uint16_t AbsImpl() const GPUnoexcept()
{
return static_cast(val & ~kSignMask);
}
@@ -385,26 +385,26 @@ struct BFloat16Impl {
/// Creates a new instance with the sign flipped.
///
/// Flipped sign instance
- GPUd() uint16_t NegateImpl() const noexcept
+ GPUd() uint16_t NegateImpl() const GPUnoexcept()
{
return IsNaN() ? val : static_cast(val ^ kSignMask);
}
public:
// uint16_t special values
- static constexpr uint16_t kSignMask = 0x8000U;
- static constexpr uint16_t kBiasedExponentMask = 0x7F80U;
- static constexpr uint16_t kPositiveInfinityBits = 0x7F80U;
- static constexpr uint16_t kNegativeInfinityBits = 0xFF80U;
- static constexpr uint16_t kPositiveQNaNBits = 0x7FC1U;
- static constexpr uint16_t kNegativeQNaNBits = 0xFFC1U;
- static constexpr uint16_t kSignaling_NaNBits = 0x7F80U;
- static constexpr uint16_t kEpsilonBits = 0x0080U;
- static constexpr uint16_t kMinValueBits = 0xFF7FU;
- static constexpr uint16_t kMaxValueBits = 0x7F7FU;
- static constexpr uint16_t kRoundToNearest = 0x7FFFU;
- static constexpr uint16_t kOneBits = 0x3F80U;
- static constexpr uint16_t kMinusOneBits = 0xBF80U;
+ static GPUglobalconstexpr() uint16_t kSignMask = 0x8000U;
+ static GPUglobalconstexpr() uint16_t kBiasedExponentMask = 0x7F80U;
+ static GPUglobalconstexpr() uint16_t kPositiveInfinityBits = 0x7F80U;
+ static GPUglobalconstexpr() uint16_t kNegativeInfinityBits = 0xFF80U;
+ static GPUglobalconstexpr() uint16_t kPositiveQNaNBits = 0x7FC1U;
+ static GPUglobalconstexpr() uint16_t kNegativeQNaNBits = 0xFFC1U;
+ static GPUglobalconstexpr() uint16_t kSignaling_NaNBits = 0x7F80U;
+ static GPUglobalconstexpr() uint16_t kEpsilonBits = 0x0080U;
+ static GPUglobalconstexpr() uint16_t kMinValueBits = 0xFF7FU;
+ static GPUglobalconstexpr() uint16_t kMaxValueBits = 0x7F7FU;
+ static GPUglobalconstexpr() uint16_t kRoundToNearest = 0x7FFFU;
+ static GPUglobalconstexpr() uint16_t kOneBits = 0x3F80U;
+ static GPUglobalconstexpr() uint16_t kMinusOneBits = 0xBF80U;
uint16_t val{0};
@@ -414,7 +414,7 @@ struct BFloat16Impl {
/// Checks if the value is negative
///
/// true if negative
- GPUd() bool IsNegative() const noexcept
+ GPUd() bool IsNegative() const GPUnoexcept()
{
return static_cast(val) < 0;
}
@@ -423,7 +423,7 @@ struct BFloat16Impl {
/// Tests if the value is NaN
///
/// true if NaN
- GPUd() bool IsNaN() const noexcept
+ GPUd() bool IsNaN() const GPUnoexcept()
{
return AbsImpl() > kPositiveInfinityBits;
}
@@ -432,7 +432,7 @@ struct BFloat16Impl {
/// Tests if the value is finite
///
/// true if finite
- GPUd() bool IsFinite() const noexcept
+ GPUd() bool IsFinite() const GPUnoexcept()
{
return AbsImpl() < kPositiveInfinityBits;
}
@@ -441,7 +441,7 @@ struct BFloat16Impl {
/// Tests if the value represents positive infinity.
///
/// true if positive infinity
- GPUd() bool IsPositiveInfinity() const noexcept
+ GPUd() bool IsPositiveInfinity() const GPUnoexcept()
{
return val == kPositiveInfinityBits;
}
@@ -450,7 +450,7 @@ struct BFloat16Impl {
/// Tests if the value represents negative infinity
///
/// true if negative infinity
- GPUd() bool IsNegativeInfinity() const noexcept
+ GPUd() bool IsNegativeInfinity() const GPUnoexcept()
{
return val == kNegativeInfinityBits;
}
@@ -459,7 +459,7 @@ struct BFloat16Impl {
/// Tests if the value is either positive or negative infinity.
///
/// True if absolute value is infinity
- GPUd() bool IsInfinity() const noexcept
+ GPUd() bool IsInfinity() const GPUnoexcept()
{
return AbsImpl() == kPositiveInfinityBits;
}
@@ -468,7 +468,7 @@ struct BFloat16Impl {
/// Tests if the value is NaN or zero. Useful for comparisons.
///
/// True if NaN or zero.
- GPUd() bool IsNaNOrZero() const noexcept
+ GPUd() bool IsNaNOrZero() const GPUnoexcept()
{
auto abs = AbsImpl();
return (abs == 0 || abs > kPositiveInfinityBits);
@@ -478,7 +478,7 @@ struct BFloat16Impl {
/// Tests if the value is normal (not zero, subnormal, infinite, or NaN).
///
/// True if so
- GPUd() bool IsNormal() const noexcept
+ GPUd() bool IsNormal() const GPUnoexcept()
{
auto abs = AbsImpl();
return (abs < kPositiveInfinityBits) // is finite
@@ -490,7 +490,7 @@ struct BFloat16Impl {
/// Tests if the value is subnormal (denormal).
///
/// True if so
- GPUd() bool IsSubnormal() const noexcept
+ GPUd() bool IsSubnormal() const GPUnoexcept()
{
auto abs = AbsImpl();
return (abs < kPositiveInfinityBits) // is finite
@@ -502,13 +502,13 @@ struct BFloat16Impl {
/// Creates an instance that represents absolute value.
///
/// Absolute value
- GPUd() Derived Abs() const noexcept { return Derived::FromBits(AbsImpl()); }
+ GPUd() Derived Abs() const GPUnoexcept() { return Derived::FromBits(AbsImpl()); }
///
/// Creates a new instance with the sign flipped.
///
/// Flipped sign instance
- GPUd() Derived Negate() const noexcept { return Derived::FromBits(NegateImpl()); }
+ GPUd() Derived Negate() const GPUnoexcept() { return Derived::FromBits(NegateImpl()); }
///
/// IEEE defines that positive and negative zero are equal, this gives us a quick equality check
@@ -518,7 +518,7 @@ struct BFloat16Impl {
/// first value
/// second value
/// True if both arguments represent zero
- GPUd() static bool AreZero(const BFloat16Impl& lhs, const BFloat16Impl& rhs) noexcept
+ GPUd() static bool AreZero(const BFloat16Impl& lhs, const BFloat16Impl& rhs) GPUnoexcept()
{
// IEEE defines that positive and negative zero are equal, this gives us a quick equality check
// for two values by or'ing the private bits together and stripping the sign. They are both zero,
@@ -528,7 +528,7 @@ struct BFloat16Impl {
};
template
-GPUdi() uint16_t BFloat16Impl::ToUint16Impl(float v) noexcept
+GPUdi() uint16_t BFloat16Impl::ToUint16Impl(float v) GPUnoexcept()
{
uint16_t result;
if (o2::gpu::CAMath::IsNaN(v)) {
@@ -566,7 +566,7 @@ GPUdi() uint16_t BFloat16Impl::ToUint16Impl(float v) noexcept
}
template
-GPUdi() float BFloat16Impl::ToFloatImpl() const noexcept
+GPUdi() float BFloat16Impl::ToFloatImpl() const GPUnoexcept()
{
#ifndef __FAST_MATH__
if (IsNaN()) {
@@ -621,7 +621,7 @@ struct Float16_t : OrtDataType::Float16Impl {
/// No conversion is done here.
///
/// 16-bit representation
- constexpr explicit Float16_t(uint16_t v) noexcept { val = v; }
+ constexpr explicit Float16_t(uint16_t v) GPUnoexcept() { val = v; }
public:
using Base = OrtDataType::Float16Impl;
@@ -636,19 +636,19 @@ struct Float16_t : OrtDataType::Float16Impl {
///
/// uint16_t bit representation of float16
/// new instance of Float16_t
- GPUd() constexpr static Float16_t FromBits(uint16_t v) noexcept { return Float16_t(v); }
+ GPUd() constexpr static Float16_t FromBits(uint16_t v) GPUnoexcept() { return Float16_t(v); }
///
/// __ctor from float. Float is converted into float16 16-bit representation.
///
/// float value
- GPUd() explicit Float16_t(float v) noexcept { val = Base::ToUint16Impl(v); }
+ GPUd() explicit Float16_t(float v) GPUnoexcept() { val = Base::ToUint16Impl(v); }
///
/// Converts float16 to float
///
/// float representation of float16 value
- GPUd() float ToFloat() const noexcept { return Base::ToFloatImpl(); }
+ GPUd() float ToFloat() const GPUnoexcept() { return Base::ToFloatImpl(); }
///
/// Checks if the value is negative
@@ -729,7 +729,7 @@ struct Float16_t : OrtDataType::Float16Impl {
///
/// User defined conversion operator. Converts Float16_t to float.
///
- GPUdi() explicit operator float() const noexcept { return ToFloat(); }
+ GPUdi() explicit operator float() const GPUnoexcept() { return ToFloat(); }
using Base::operator==;
using Base::operator!=;
@@ -765,7 +765,7 @@ struct BFloat16_t : OrtDataType::BFloat16Impl {
/// No conversion is done.
///
/// 16-bit bfloat16 value
- constexpr explicit BFloat16_t(uint16_t v) noexcept { val = v; }
+ constexpr explicit BFloat16_t(uint16_t v) GPUnoexcept() { val = v; }
public:
using Base = OrtDataType::BFloat16Impl;
@@ -777,19 +777,19 @@ struct BFloat16_t : OrtDataType::BFloat16Impl {
///
/// uint16_t bit representation of bfloat16
/// new instance of BFloat16_t
- GPUd() static constexpr BFloat16_t FromBits(uint16_t v) noexcept { return BFloat16_t(v); }
+ GPUd() static constexpr BFloat16_t FromBits(uint16_t v) GPUnoexcept() { return BFloat16_t(v); }
///
/// __ctor from float. Float is converted into bfloat16 16-bit representation.
///
/// float value
- GPUd() explicit BFloat16_t(float v) noexcept { val = Base::ToUint16Impl(v); }
+ GPUd() explicit BFloat16_t(float v) GPUnoexcept() { val = Base::ToUint16Impl(v); }
///
/// Converts bfloat16 to float
///
/// float representation of bfloat16 value
- GPUd() float ToFloat() const noexcept { return Base::ToFloatImpl(); }
+ GPUd() float ToFloat() const GPUnoexcept() { return Base::ToFloatImpl(); }
///
/// Checks if the value is negative
@@ -870,13 +870,13 @@ struct BFloat16_t : OrtDataType::BFloat16Impl {
///
/// User defined conversion operator. Converts BFloat16_t to float.
///
- GPUdi() explicit operator float() const noexcept { return ToFloat(); }
+ GPUdi() explicit operator float() const GPUnoexcept() { return ToFloat(); }
// We do not have an inherited impl for the below operators
// as the internal class implements them a little differently
- bool operator==(const BFloat16_t& rhs) const noexcept;
- bool operator!=(const BFloat16_t& rhs) const noexcept { return !(*this == rhs); }
- bool operator<(const BFloat16_t& rhs) const noexcept;
+ bool operator==(const BFloat16_t& rhs) const GPUnoexcept();
+ bool operator!=(const BFloat16_t& rhs) const GPUnoexcept() { return !(*this == rhs); }
+ bool operator<(const BFloat16_t& rhs) const GPUnoexcept();
};
static_assert(sizeof(BFloat16_t) == sizeof(uint16_t), "Sizes must match");
diff --git a/Common/MathUtils/include/MathUtils/Cartesian.h b/Common/MathUtils/include/MathUtils/Cartesian.h
index e61b10a7caee9..99b946069f089 100644
--- a/Common/MathUtils/include/MathUtils/Cartesian.h
+++ b/Common/MathUtils/include/MathUtils/Cartesian.h
@@ -152,7 +152,9 @@ class Rotation2D
};
using Rotation2Df_t = Rotation2D;
+#ifndef __METAL__
using Rotation2Dd_t = Rotation2D;
+#endif
#if (!defined(GPUCA_STANDALONE) || !defined(DGPUCA_NO_ROOT)) && !defined(GPUCA_GPUCODE) && !defined(GPUCOMMONRTYPES_H_ACTIVE)
diff --git a/Common/MathUtils/include/MathUtils/Primitive2D.h b/Common/MathUtils/include/MathUtils/Primitive2D.h
index 052926d111594..e654b9cfb1f13 100644
--- a/Common/MathUtils/include/MathUtils/Primitive2D.h
+++ b/Common/MathUtils/include/MathUtils/Primitive2D.h
@@ -28,17 +28,23 @@ namespace math_utils
template
using CircleXY = detail::CircleXY;
using CircleXYf_t = detail::CircleXY;
+#ifndef __METAL__
using CircleXYd_t = detail::CircleXY;
+#endif
template
using IntervalXY = detail::IntervalXY;
using IntervalXYf_t = detail::IntervalXY;
+#ifndef __METAL__
using IntervalXYd_t = detail::IntervalXY;
+#endif
template
using Bracket = detail::Bracket;
using Bracketf_t = detail::Bracket;
+#ifndef __METAL__
using Bracketd_t = detail::Bracket;
+#endif
} // namespace math_utils
} // namespace o2
diff --git a/Common/MathUtils/include/MathUtils/SMatrixGPU.h b/Common/MathUtils/include/MathUtils/SMatrixGPU.h
index 8158a93666a92..497a142a1902c 100644
--- a/Common/MathUtils/include/MathUtils/SMatrixGPU.h
+++ b/Common/MathUtils/include/MathUtils/SMatrixGPU.h
@@ -340,7 +340,11 @@ class MatRepSymGPU
static GPUdi() int off(int i)
{
+#ifdef __METAL__ // MSL rejects variables declared static at function scope
+ constexpr auto v = row_offsets_utils::make(off1);
+#else
static constexpr auto v = row_offsets_utils::make(off1);
+#endif
return v[i];
}
@@ -518,7 +522,7 @@ class SMatrixGPU
R mRep;
};
-#ifndef __OPENCL__ // TODO: current C++ for OpenCL 2021 is at C++17, so no concepts. But we don't need this trick for OpenCL anyway, so we can just hide it.
+#if !defined(__OPENCL__) && !defined(__METAL__) // TODO: current C++ for OpenCL 2021 and MSL 4.1 are both at C++17, so no concepts. But we don't need this trick there anyway, so we can just hide it.
template
requires(sizeof(typename X::traits_type::pos_type) != 0) // do not provide a template to fair::Logger, etc... (pos_type is a member type of all std::ostream classes)
GPUd() X& operator<<(Y& y, const SMatrixGPU&)
@@ -1429,14 +1433,18 @@ template
template
GPUdi() SMatrixGPU& SMatrixGPU::operator*=(const SMatrixGPU& rhs)
{
- return operator=(*this* rhs);
+ // the product is an expression evaluated element by element, and every element
+ // of it reads the whole of *this, so it has to be materialised first
+ const SMatrixGPU tmp(*this * rhs);
+ return operator=(tmp);
}
template
template
GPUdi() SMatrixGPU& SMatrixGPU::operator*=(const Expr& rhs)
{
- return operator=(*this* rhs);
+ const SMatrixGPU tmp(*this * rhs);
+ return operator=(tmp);
}
template
diff --git a/Common/MathUtils/include/MathUtils/Utils.h b/Common/MathUtils/include/MathUtils/Utils.h
index 3c51245fc6c29..4d17efb0fe6c2 100644
--- a/Common/MathUtils/include/MathUtils/Utils.h
+++ b/Common/MathUtils/include/MathUtils/Utils.h
@@ -32,111 +32,133 @@ GPUdi() float to02Pi(float phi)
return detail::to02Pi(phi);
}
+#ifndef __METAL__ // MSL has no double; every other backend keeps these
GPUdi() double to02Pid(double phi)
{
return detail::to02Pi(phi);
}
+#endif
GPUdi() void bringTo02Pi(float& phi)
{
detail::bringTo02Pi(phi);
}
+#ifndef __METAL__ // MSL has no double; every other backend keeps these
GPUdi() void bringTo02Pid(double& phi)
{
detail::bringTo02Pi(phi);
}
+#endif
inline float toPMPiGen(float phi)
{
return detail::toPMPiGen(phi);
}
+#ifndef __METAL__ // MSL has no double; every other backend keeps these
inline double toPMPiGend(double phi)
{
return detail::toPMPiGen(phi);
}
+#endif
inline void bringToPMPiGen(float& phi)
{
detail::bringToPMPiGen(phi);
}
+#ifndef __METAL__ // MSL has no double; every other backend keeps these
inline void bringToPMPiGend(double& phi)
{
detail::bringToPMPiGen(phi);
}
+#endif
inline float to02PiGen(float phi)
{
return detail::to02PiGen(phi);
}
+#ifndef __METAL__ // MSL has no double; every other backend keeps these
inline double to02PiGend(double phi)
{
return detail::to02PiGen(phi);
}
+#endif
inline void bringTo02PiGen(float& phi)
{
detail::bringTo02PiGen(phi);
}
+#ifndef __METAL__ // MSL has no double; every other backend keeps these
inline void bringTo02PiGend(double& phi)
{
detail::bringTo02PiGen(phi);
}
+#endif
inline float toPMPi(float phi)
{
return detail::toPMPi(phi);
}
+#ifndef __METAL__ // MSL has no double; every other backend keeps these
inline double toPMPid(double phi)
{
return detail::toPMPi(phi);
}
+#endif
inline void bringToPMPi(float& phi)
{
return detail::bringToPMPi(phi);
}
+#ifndef __METAL__ // MSL has no double; every other backend keeps these
inline void bringToPMPid(double& phi)
{
return detail::bringToPMPi(phi);
}
+#endif
GPUdi() void sincos(float ang, float& s, float& c)
{
detail::sincos(ang, s, c);
}
#ifndef __OPENCL__
+#ifndef __METAL__ // MSL has no double; every other backend keeps these
GPUdi() void sincosd(double ang, double& s, double& c)
{
detail::sincos(ang, s, c);
}
#endif
+#endif
GPUdi() void rotateZ(float xL, float yL, float& xG, float& yG, float snAlp, float csAlp)
{
return detail::rotateZ(xL, yL, xG, yG, snAlp, csAlp);
}
+#ifndef __METAL__ // MSL has no double; every other backend keeps these
GPUdi() void rotateZd(double xL, double yL, double& xG, double& yG, double snAlp, double csAlp)
{
return detail::rotateZ(xL, yL, xG, yG, snAlp, csAlp);
}
+#endif
GPUdi() void rotateZInv(float xG, float yG, float& xL, float& yL, float snAlp, float csAlp)
{
detail::rotateZInv(xG, yG, xL, yL, snAlp, csAlp);
}
+#ifndef __METAL__ // MSL has no double; every other backend keeps these
GPUdi() void rotateZInvd(double xG, double yG, double& xL, double& yL, double snAlp, double csAlp)
{
detail::rotateZInv(xG, yG, xL, yL, snAlp, csAlp);
}
+#endif
#ifndef GPUCA_GPUCODE_DEVICE
inline std::tuple rotateZInv(float xG, float yG, float snAlp, float csAlp)
@@ -185,40 +207,48 @@ inline int angle2Sector(float phi)
return detail::angle2Sector(phi);
}
+#ifndef __METAL__ // MSL has no double; every other backend keeps these
inline int angle2Sectord(double phi)
{
return detail::angle2Sector(phi);
}
+#endif
inline float sector2Angle(int sect)
{
return detail::sector2Angle(sect);
}
+#ifndef __METAL__ // MSL has no double; every other backend keeps these
inline double sector2Angled(int sect)
{
return detail::sector2Angle(sect);
}
+#endif
inline float angle2Alpha(float phi)
{
return detail::angle2Alpha(phi);
}
+#ifndef __METAL__ // MSL has no double; every other backend keeps these
inline double angle2Alphad(double phi)
{
return detail::angle2Alpha(phi);
}
+#endif
GPUhdi() constexpr float fastATan2(float y, float x)
{
return detail::fastATan2(y, x);
}
+#ifndef __METAL__ // MSL has no double; every other backend keeps these
GPUhdi() constexpr double fastATan2d(double y, double x)
{
return detail::fastATan2(y, x);
}
+#endif
template
GPUhdi() T min(const T x, const T y)
@@ -226,10 +256,12 @@ GPUhdi() T min(const T x, const T y)
return detail::min(x, y);
};
+#ifndef __METAL__ // MSL has no double; every other backend keeps these
GPUhdi() double mind(const double x, const double y)
{
return detail::min(x, y);
};
+#endif
template
GPUhdi() T max(const T x, const T y)
@@ -237,130 +269,156 @@ GPUhdi() T max(const T x, const T y)
return detail::max(x, y);
};
+#ifndef __METAL__ // MSL has no double; every other backend keeps these
GPUhdi() double maxd(const double x, const double y)
{
return detail::max(x, y);
};
+#endif
GPUhdi() float sqrt(float x)
{
return detail::sqrt(x);
};
+#ifndef __METAL__ // MSL has no double; every other backend keeps these
GPUhdi() double sqrtd(double x)
{
return detail::sqrt(x);
};
+#endif
GPUhdi() float abs(float x)
{
return detail::abs(x);
};
+#ifndef __METAL__ // MSL has no double; every other backend keeps these
GPUhdi() double absd(double x)
{
return detail::abs(x);
};
+#endif
GPUdi() float asin(float x)
{
return detail::asin(x);
};
+#ifndef __METAL__ // MSL has no double; every other backend keeps these
GPUdi() double asind(double x)
{
return detail::asin(x);
};
+#endif
GPUdi() float atan(float x)
{
return detail::atan(x);
};
+#ifndef __METAL__ // MSL has no double; every other backend keeps these
GPUdi() double atand(double x)
{
return detail::atan(x);
};
+#endif
GPUdi() float atan2(float y, float x)
{
return detail::atan2(y, x);
};
+#ifndef __METAL__ // MSL has no double; every other backend keeps these
GPUdi() double atan2d(double y, double x)
{
return detail::atan2(y, x);
};
+#endif
GPUdi() float sin(float x)
{
return detail::sin(x);
};
+#ifndef __METAL__ // MSL has no double; every other backend keeps these
GPUdi() double sind(double x)
{
return detail::sin(x);
};
+#endif
GPUdi() float cos(float x)
{
return detail::cos(x);
};
+#ifndef __METAL__ // MSL has no double; every other backend keeps these
GPUdi() double cosd(double x)
{
return detail::cos(x);
};
+#endif
GPUdi() float tan(float x)
{
return detail::tan(x);
};
+#ifndef __METAL__ // MSL has no double; every other backend keeps these
GPUdi() double tand(double x)
{
return detail::tan(x);
};
+#endif
GPUdi() float twoPi()
{
return detail::twoPi();
};
+#ifndef __METAL__ // MSL has no double; every other backend keeps these
GPUdi() double twoPid()
{
return detail::twoPi();
};
+#endif
GPUdi() float pi()
{
return detail::pi();
}
+#ifndef __METAL__ // MSL has no double; every other backend keeps these
GPUdi() double pid()
{
return detail::pi();
}
+#endif
GPUdi() int nint(float x)
{
return detail::nint(x);
};
+#ifndef __METAL__ // MSL has no double; every other backend keeps these
GPUdi() int nintd(double x)
{
return detail::nint(x);
};
+#endif
GPUdi() bool finite(float x)
{
return detail::finite(x);
}
+#ifndef __METAL__ // MSL has no double; every other backend keeps these
GPUdi() bool finited(double x)
{
return detail::finite(x);
}
+#endif
GPUdi() unsigned int clz(unsigned int val)
{
@@ -377,12 +435,16 @@ GPUdi() float log(float x)
return detail::log(x);
};
+#ifndef __METAL__ // MSL has no double; every other backend keeps these
GPUdi() double logd(double x)
{
return detail::log(x);
};
+#endif
+#ifndef __METAL__
using detail::StatAccumulator;
+#endif
using detail::bit2Mask;
using detail::numberOfBitsSet;
diff --git a/Common/MathUtils/include/MathUtils/detail/Bracket.h b/Common/MathUtils/include/MathUtils/detail/Bracket.h
index 2da6949c4a6f8..450174dc9737f 100644
--- a/Common/MathUtils/include/MathUtils/detail/Bracket.h
+++ b/Common/MathUtils/include/MathUtils/detail/Bracket.h
@@ -16,6 +16,8 @@
#ifndef ALICEO2_BRACKET_H
#define ALICEO2_BRACKET_H
+#include "GPUCommonDef.h"
+
#include
#ifndef GPUCA_GPUCODE_DEVICE
#include
@@ -53,9 +55,9 @@ class Bracket
bool operator==(const Bracket& other) const;
bool operator!=(const Bracket& other) const;
- void setMax(T v) noexcept;
- void setMin(T v) noexcept;
- void set(T minv, T maxv) noexcept;
+ void setMax(T v) GPUnoexcept();
+ void setMin(T v) GPUnoexcept();
+ void set(T minv, T maxv) GPUnoexcept();
T& getMax();
T& getMin();
@@ -129,19 +131,19 @@ inline bool Bracket::operator!=(const Bracket& rhs) const
}
template
-inline void Bracket::setMax(T v) noexcept
+inline void Bracket::setMax(T v) GPUnoexcept()
{
mMax = v;
}
template
-inline void Bracket::setMin(T v) noexcept
+inline void Bracket::setMin(T v) GPUnoexcept()
{
mMin = v;
}
template
-inline void Bracket::set(T minv, T maxv) noexcept
+inline void Bracket::set(T minv, T maxv) GPUnoexcept()
{
this->setMin(minv);
this->setMax(maxv);
diff --git a/Common/MathUtils/include/MathUtils/detail/StatAccumulator.h b/Common/MathUtils/include/MathUtils/detail/StatAccumulator.h
index abb8a716cc5ee..2d0eb94fd951a 100644
--- a/Common/MathUtils/include/MathUtils/detail/StatAccumulator.h
+++ b/Common/MathUtils/include/MathUtils/detail/StatAccumulator.h
@@ -27,6 +27,7 @@ namespace math_utils
namespace detail
{
+#ifndef __METAL__ // host-only accumulator, and MSL has no double
struct StatAccumulator {
// mean / RMS accumulator
double sum = 0.;
@@ -83,6 +84,7 @@ struct StatAccumulator {
n = 0;
}
};
+#endif
} // namespace detail
} // namespace math_utils
diff --git a/Common/MathUtils/include/MathUtils/detail/trigonometric.h b/Common/MathUtils/include/MathUtils/detail/trigonometric.h
index e13d965663dc9..e0bf675d696cc 100644
--- a/Common/MathUtils/include/MathUtils/detail/trigonometric.h
+++ b/Common/MathUtils/include/MathUtils/detail/trigonometric.h
@@ -48,7 +48,7 @@ GPUhdi() T to02Pi(T phi)
template
GPUhdi() void bringTo02Pi(T& phi)
{
- phi = to02Pi(phi);
+ phi = to02Pi(phi);
}
template
@@ -68,7 +68,7 @@ inline T to02PiGen(T phi)
template
inline void bringTo02PiGen(T& phi)
{
- phi = to02PiGen(phi);
+ phi = to02PiGen(phi);
}
template
@@ -87,7 +87,7 @@ GPUhdi() T toPMPi(T phi)
template
GPUhdi() void bringToPMPi(T& phi)
{
- phi = toPMPi(phi);
+ phi = toPMPi(phi);
}
template
@@ -107,10 +107,10 @@ inline T toPMPiGen(T phi)
template
inline void bringToPMPiGen(T& phi)
{
- phi = toPMPiGen(phi);
+ phi = toPMPiGen(phi);
}
-#ifdef __OPENCL__ // TODO: get rid of that stupid workaround for OpenCL template address spaces
+#if defined(__OPENCL__) || defined(__METAL__) // TODO: get rid of that stupid workaround for OpenCL template address spaces
template
GPUhdi() void sincos(T ang, S& s, U& c)
{
@@ -122,12 +122,14 @@ GPUhdi() void sincos(T ang, T& s, T& c)
{
return o2::gpu::GPUCommonMath::SinCos(ang, s, c);
}
+#ifndef __METAL__ // MSL has no double; the primary template still serves float
template <>
GPUhdi() void sincos(double ang, double& s, double& c)
{
return o2::gpu::GPUCommonMath::SinCosd(ang, s, c);
}
#endif
+#endif
#ifndef GPUCA_GPUCODE_DEVICE
@@ -358,11 +360,13 @@ GPUdi() T twoPi()
return o2::gpu::GPUCommonMath::TwoPi();
};
+#ifndef __METAL__ // MSL has no double; the primary template still serves float
template <>
GPUdi() double twoPi()
{
return o2::constants::math::TwoPI;
};
+#endif
template
GPUdi() T pi()
@@ -370,11 +374,13 @@ GPUdi() T pi()
return o2::gpu::GPUCommonMath::Pi();
}
+#ifndef __METAL__ // MSL has no double; the primary template still serves float
template <>
GPUdi() double pi()
{
return o2::constants::math::PI;
}
+#endif
#ifndef GPUCA_GPUCODE_DEVICE
template <>
diff --git a/DataFormats/Detectors/TPC/include/DataFormatsTPC/Defs.h b/DataFormats/Detectors/TPC/include/DataFormatsTPC/Defs.h
index a5be0da32f641..15f72d6e2f68e 100644
--- a/DataFormats/Detectors/TPC/include/DataFormatsTPC/Defs.h
+++ b/DataFormats/Detectors/TPC/include/DataFormatsTPC/Defs.h
@@ -19,6 +19,8 @@
#ifndef AliceO2_TPC_Defs_H
#define AliceO2_TPC_Defs_H
+#include "GPUCommonDouble.h"
+
#include "GPUCommonDef.h"
#ifndef GPUCA_GPUCODE_DEVICE
@@ -42,9 +44,9 @@ enum Side { A = 0,
GPUglobalconstexpr() unsigned char SECTORSPERSIDE = 18;
GPUglobalconstexpr() unsigned char SIDES = 2;
-constexpr double PI = 3.14159265358979323846;
-constexpr double TWOPI = 2. * PI;
-constexpr double SECPHIWIDTH = TWOPI / 18.;
+GPUglobalconstexpr() o2::gpu::GPUdoubleValue PI = 3.14159265358979323846;
+GPUglobalconstexpr() o2::gpu::GPUdoubleValue TWOPI = 2. * PI;
+GPUglobalconstexpr() o2::gpu::GPUdoubleValue SECPHIWIDTH = TWOPI / 18.;
/// TPC ROC types
enum RocType { IROC = 0,
diff --git a/DataFormats/Reconstruction/include/ReconstructionDataFormats/HelixHelper.h b/DataFormats/Reconstruction/include/ReconstructionDataFormats/HelixHelper.h
index 47de5457cea16..0785afa7553ac 100644
--- a/DataFormats/Reconstruction/include/ReconstructionDataFormats/HelixHelper.h
+++ b/DataFormats/Reconstruction/include/ReconstructionDataFormats/HelixHelper.h
@@ -17,6 +17,7 @@
#define _ALICEO2_HELIX_HELPER_
#include "CommonConstants/MathConstants.h"
+#include "GPUCommonDouble.h"
#include "MathUtils/Utils.h"
#include "MathUtils/Primitive2D.h"
@@ -247,8 +248,8 @@ struct CrossInfo {
auto tgp = trcL.getSnp() * cspi;
float kx = traxL.c - traxL.s * tgp;
float ky = traxL.s + traxL.c * tgp;
- double dk = dx * kx + dy * ky;
- double det = dk * dk - cspi2 * (dx * dx + dy * dy - traxH.rC * traxH.rC);
+ o2::gpu::GPUdoubleCalc dk = dx * kx + dy * ky;
+ o2::gpu::GPUdoubleCalc det = dk * dk - cspi2 * (dx * dx + dy * dy - traxH.rC * traxH.rC);
if (det > 0) { // 2 crossings
det = o2::gpu::GPUCommonMath::Sqrt(det);
float t0 = (-dk + det) * cspi2;
diff --git a/DataFormats/Reconstruction/include/ReconstructionDataFormats/Track.h b/DataFormats/Reconstruction/include/ReconstructionDataFormats/Track.h
index c4b17158f7dc5..2d4d67574a12a 100644
--- a/DataFormats/Reconstruction/include/ReconstructionDataFormats/Track.h
+++ b/DataFormats/Reconstruction/include/ReconstructionDataFormats/Track.h
@@ -25,11 +25,15 @@ namespace track
{
using TrackParF = TrackParametrization;
+#ifndef __METAL__
using TrackParD = TrackParametrization;
+#endif
using TrackPar = TrackParF;
using TrackParCovF = TrackParametrizationWithError;
+#ifndef __METAL__
using TrackParCovD = TrackParametrizationWithError;
+#endif
using TrackParCov = TrackParCovF;
} // namespace track
diff --git a/DataFormats/Reconstruction/include/ReconstructionDataFormats/TrackParametrizationWithError.h b/DataFormats/Reconstruction/include/ReconstructionDataFormats/TrackParametrizationWithError.h
index 81280d090be71..d120e4c15077f 100644
--- a/DataFormats/Reconstruction/include/ReconstructionDataFormats/TrackParametrizationWithError.h
+++ b/DataFormats/Reconstruction/include/ReconstructionDataFormats/TrackParametrizationWithError.h
@@ -17,6 +17,7 @@
#ifndef INCLUDE_RECONSTRUCTIONDATAFORMATS_TRACKPARAMETRIZATIONWITHERROR_H_
#define INCLUDE_RECONSTRUCTIONDATAFORMATS_TRACKPARAMETRIZATIONWITHERROR_H_
+#include "GPUCommonDouble.h"
#include "ReconstructionDataFormats/TrackParametrization.h"
#include
@@ -41,8 +42,8 @@ class TrackParametrizationWithError : public TrackParametrization
#endif
using covMat_t = std::array;
- using MatrixDSym5 = o2::math_utils::SMatrix>;
- using MatrixD5 = o2::math_utils::SMatrix>;
+ using MatrixDSym5 = o2::math_utils::SMatrix>;
+ using MatrixD5 = o2::math_utils::SMatrix>;
GPUhd() TrackParametrizationWithError();
GPUd() TrackParametrizationWithError(value_t x, value_t alpha, const params_t& par, const covMat_t& cov, int charge = 1, const PID pid = PID::Pion);
diff --git a/DataFormats/Reconstruction/include/ReconstructionDataFormats/TrackUtils.h b/DataFormats/Reconstruction/include/ReconstructionDataFormats/TrackUtils.h
index 6befcd0dfc898..1e96ab224a308 100644
--- a/DataFormats/Reconstruction/include/ReconstructionDataFormats/TrackUtils.h
+++ b/DataFormats/Reconstruction/include/ReconstructionDataFormats/TrackUtils.h
@@ -138,7 +138,11 @@ GPUd() value_T BetheBlochSolid(value_T bg, value_T rho, value_T kp1, value_T kp2
if (x > kp2) {
d2 = lhwI + x - value_T(0.5);
} else if (x > kp1) {
+#ifdef __METAL__ // MSL has no double
+ float r = (kp2 - x) / (kp2 - kp1);
+#else
double r = (kp2 - x) / (kp2 - kp1);
+#endif
d2 = lhwI + x - value_T(0.5) + (value_T(0.5) - lhwI - kp1) * r * r * r;
}
auto dedx = mK * meanZA / beta2 * (value_T(0.5) * gpu::CAMath::Log(value_T(2) * me * bg2 * maxT / (meanI * meanI)) - beta2 - d2);
diff --git a/DataFormats/Reconstruction/src/TrackParametrization.cxx b/DataFormats/Reconstruction/src/TrackParametrization.cxx
index fb398bbcf07cf..cdf8f35b7cfb0 100644
--- a/DataFormats/Reconstruction/src/TrackParametrization.cxx
+++ b/DataFormats/Reconstruction/src/TrackParametrization.cxx
@@ -15,6 +15,7 @@
/// @brief
#include "ReconstructionDataFormats/TrackParametrization.h"
+#include "GPUCommonDouble.h"
#include "ReconstructionDataFormats/Vertex.h"
#include "ReconstructionDataFormats/DCA.h"
#include
@@ -477,7 +478,7 @@ GPUd() bool TrackParametrization::getYZAt(value_t xk, value_t b, value_
if (gpu::CAMath::Abs(r2) < constants::math::Almost0) {
return false;
}
- double dy2dx = (f1 + f2) / (r1 + r2);
+ GPUdoubleCalc dy2dx = (f1 + f2) / (r1 + r2);
y += dx * dy2dx;
if (gpu::CAMath::Abs(x2r) < 0.05f) {
z += dx * (r2 + f2 * dy2dx) * getTgl();
@@ -692,7 +693,7 @@ GPUd() bool TrackParametrization::getXatLabR(value_t r, value_t& x, val
// DirOutward (==1) - go along the track (increasing mX)
// DirInward (==-1) - go backward (decreasing mX)
//
- const double fy = mP[0], sn = mP[2];
+ const GPUdoubleCalc fy = mP[0], sn = mP[2];
const value_t kEps = 1.e-6;
//
if (gpu::CAMath::Abs(getSnp()) > constants::math::Almost1) {
@@ -711,7 +712,7 @@ GPUd() bool TrackParametrization::getXatLabR(value_t r, value_t& x, val
if (r0 <= constants::math::Almost0) {
return false; // the track is concentric to circle
}
- double tR2r0 = 1., g = 0., tmp = 0.;
+ GPUdoubleCalc tR2r0 = 1., g = 0., tmp = 0.;
if (gpu::CAMath::Abs(circle.rC - r0) > kEps) {
tR2r0 = circle.rC / r0;
g = 0.5f * (r * r / (r0 * circle.rC) - tR2r0 - 1.f / tR2r0);
@@ -786,7 +787,7 @@ GPUd() bool TrackParametrization::getXatLabR(value_t r, value_t& x, val
}
// this is a straight track
if (gpu::CAMath::Abs(sn) >= constants::math::Almost1) { // || to Y axis
- double det = (r - mX) * (r + mX);
+ GPUdoubleCalc det = (r - mX) * (r + mX);
if (det < 0.f) {
return false; // does not reach raduis r
}
@@ -815,7 +816,7 @@ GPUd() bool TrackParametrization::getXatLabR(value_t r, value_t& x, val
}
}
} else if (gpu::CAMath::Abs(sn) <= constants::math::Almost0) { // || to X axis
- double det = (r - fy) * (r + fy);
+ GPUdoubleCalc det = (r - fy) * (r + fy);
if (det < 0.) {
return false; // does not reach raduis r
}
diff --git a/DataFormats/Reconstruction/src/TrackParametrizationWithError.cxx b/DataFormats/Reconstruction/src/TrackParametrizationWithError.cxx
index 748cb47094d26..7b19e7a83f867 100644
--- a/DataFormats/Reconstruction/src/TrackParametrizationWithError.cxx
+++ b/DataFormats/Reconstruction/src/TrackParametrizationWithError.cxx
@@ -10,6 +10,7 @@
// or submit itself to any jurisdiction.
#include "ReconstructionDataFormats/TrackParametrizationWithError.h"
+#include "GPUCommonDouble.h"
#include "ReconstructionDataFormats/Vertex.h"
#include "ReconstructionDataFormats/DCA.h"
#include "CommonConstants/MathConstants.h"
@@ -67,9 +68,9 @@ GPUd() bool TrackParametrizationWithError::propagateTo(value_t xk, valu
if (gpu::CAMath::Abs(r2) < constants::math::Almost0) {
return false;
}
- double r1pr2Inv = 1. / (r1 + r2);
- double dy2dx = (f1 + f2) * r1pr2Inv;
- const auto dy2dxF = static_cast(dy2dx); // the parameter update does not need the double
+ GPUdoubleCalc r1pr2Inv = 1. / (r1 + r2);
+ GPUdoubleCalc dy2dx = (f1 + f2) * r1pr2Inv;
+ const auto dy2dxF = static_cast(dy2dx); // the parameter update does not need the GPUdoubleCalc
bool arcz = gpu::CAMath::Abs(x2r) > 0.05f;
params_t dP{0.f};
if (arcz) {
@@ -110,33 +111,33 @@ GPUd() bool TrackParametrizationWithError::propagateTo(value_t xk, valu
// evaluate matrix in double prec.
value_t kb = bz * constants::math::B2C;
- double r2inv = 1. / r2, r1inv = 1. / r1;
- double dx2r1pr2 = dx * r1pr2Inv;
+ GPUdoubleCalc r2inv = 1. / r2, r1inv = 1. / r1;
+ GPUdoubleCalc dx2r1pr2 = dx * r1pr2Inv;
- double hh = dx2r1pr2 * r2inv * (1. + r1 * r2 + f1 * f2), jj = dx * (dy2dx - f2 * r2inv);
- double f02 = hh * r1inv;
- double f04 = hh * dx2r1pr2 * kb;
- double f24 = dx * kb; // x2r/mP[kQ2Pt];
- double f12 = this->getTgl() * (f02 * f2 + jj);
- double f13 = dx * (r2 + f2 * dy2dx);
- double f14 = this->getTgl() * (f04 * f2 + jj * f24);
+ GPUdoubleCalc hh = dx2r1pr2 * r2inv * (1. + r1 * r2 + f1 * f2), jj = dx * (dy2dx - f2 * r2inv);
+ GPUdoubleCalc f02 = hh * r1inv;
+ GPUdoubleCalc f04 = hh * dx2r1pr2 * kb;
+ GPUdoubleCalc f24 = dx * kb; // x2r/mP[kQ2Pt];
+ GPUdoubleCalc f12 = this->getTgl() * (f02 * f2 + jj);
+ GPUdoubleCalc f13 = dx * (r2 + f2 * dy2dx);
+ GPUdoubleCalc f14 = this->getTgl() * (f04 * f2 + jj * f24);
// b = C*ft
- double b00 = f02 * c20 + f04 * c40, b01 = f12 * c20 + f14 * c40 + f13 * c30;
- double b02 = f24 * c40;
- double b10 = f02 * c21 + f04 * c41, b11 = f12 * c21 + f14 * c41 + f13 * c31;
- double b12 = f24 * c41;
- double b20 = f02 * c22 + f04 * c42, b21 = f12 * c22 + f14 * c42 + f13 * c32;
- double b22 = f24 * c42;
- double b40 = f02 * c42 + f04 * c44, b41 = f12 * c42 + f14 * c44 + f13 * c43;
- double b42 = f24 * c44;
- double b30 = f02 * c32 + f04 * c43, b31 = f12 * c32 + f14 * c43 + f13 * c33;
- double b32 = f24 * c43;
+ GPUdoubleCalc b00 = f02 * c20 + f04 * c40, b01 = f12 * c20 + f14 * c40 + f13 * c30;
+ GPUdoubleCalc b02 = f24 * c40;
+ GPUdoubleCalc b10 = f02 * c21 + f04 * c41, b11 = f12 * c21 + f14 * c41 + f13 * c31;
+ GPUdoubleCalc b12 = f24 * c41;
+ GPUdoubleCalc b20 = f02 * c22 + f04 * c42, b21 = f12 * c22 + f14 * c42 + f13 * c32;
+ GPUdoubleCalc b22 = f24 * c42;
+ GPUdoubleCalc b40 = f02 * c42 + f04 * c44, b41 = f12 * c42 + f14 * c44 + f13 * c43;
+ GPUdoubleCalc b42 = f24 * c44;
+ GPUdoubleCalc b30 = f02 * c32 + f04 * c43, b31 = f12 * c32 + f14 * c43 + f13 * c33;
+ GPUdoubleCalc b32 = f24 * c43;
// a = f*b = f*C*ft
- double a00 = f02 * b20 + f04 * b40, a01 = f02 * b21 + f04 * b41, a02 = f02 * b22 + f04 * b42;
- double a11 = f12 * b21 + f14 * b41 + f13 * b31, a12 = f12 * b22 + f14 * b42 + f13 * b32;
- double a22 = f24 * b42;
+ GPUdoubleCalc a00 = f02 * b20 + f04 * b40, a01 = f02 * b21 + f04 * b41, a02 = f02 * b22 + f04 * b42;
+ GPUdoubleCalc a11 = f12 * b21 + f14 * b41 + f13 * b31, a12 = f12 * b22 + f14 * b42 + f13 * b32;
+ GPUdoubleCalc a22 = f24 * b42;
// F*C*Ft = C + (b + bt + a)
c00 += b00 + b00 + a00;
@@ -180,17 +181,17 @@ GPUd() bool TrackParametrizationWithError::propagateTo(value_t xk, Trac
}
value_t kb = bz * constants::math::B2C;
// evaluate in double prec.
- double snpRef0 = linRef0.getSnp(), cspRef0 = gpu::CAMath::Sqrt((1 - snpRef0) * (1 + snpRef0));
- double snpRef1 = linRef1.getSnp(), cspRef1 = gpu::CAMath::Sqrt((1 - snpRef1) * (1 + snpRef1));
- double cspRef0Inv = 1 / cspRef0, cspRef1Inv = 1 / cspRef1, cc = cspRef0 + cspRef1, ccInv = 1 / cc, dy2dx = (snpRef0 + snpRef1) * ccInv;
- double dxccInv = dx * ccInv, hh = dxccInv * cspRef1Inv * (1 + cspRef0 * cspRef1 + snpRef0 * snpRef1), jj = dx * (dy2dx - snpRef1 * cspRef1Inv);
-
- double f02 = hh * cspRef0Inv;
- double f04 = hh * dxccInv * kb;
- double f24 = dx * kb;
- double f12 = linRef0.getTgl() * (f02 * snpRef1 + jj);
- double f13 = dx * (cspRef1 + snpRef1 * dy2dx); // dS
- double f14 = linRef0.getTgl() * (f04 * snpRef1 + jj * f24);
+ GPUdoubleCalc snpRef0 = linRef0.getSnp(), cspRef0 = gpu::CAMath::Sqrt((1.f - snpRef0) * (1.f + snpRef0));
+ GPUdoubleCalc snpRef1 = linRef1.getSnp(), cspRef1 = gpu::CAMath::Sqrt((1.f - snpRef1) * (1.f + snpRef1));
+ GPUdoubleCalc cspRef0Inv = 1.f / cspRef0, cspRef1Inv = 1.f / cspRef1, cc = cspRef0 + cspRef1, ccInv = 1.f / cc, dy2dx = (snpRef0 + snpRef1) * ccInv;
+ GPUdoubleCalc dxccInv = dx * ccInv, hh = dxccInv * cspRef1Inv * (1.f + cspRef0 * cspRef1 + snpRef0 * snpRef1), jj = dx * (dy2dx - snpRef1 * cspRef1Inv);
+
+ GPUdoubleCalc f02 = hh * cspRef0Inv;
+ GPUdoubleCalc f04 = hh * dxccInv * kb;
+ GPUdoubleCalc f24 = dx * kb;
+ GPUdoubleCalc f12 = linRef0.getTgl() * (f02 * snpRef1 + jj);
+ GPUdoubleCalc f13 = dx * (cspRef1 + snpRef1 * dy2dx); // dS
+ GPUdoubleCalc f14 = linRef0.getTgl() * (f04 * snpRef1 + jj * f24);
// difference between the current and reference state
value_t diff[5];
@@ -215,21 +216,21 @@ GPUd() bool TrackParametrizationWithError::propagateTo(value_t xk, Trac
&c44 = mC[kSigQ2Pt2];
// b = C*ft
- double b00 = f02 * c20 + f04 * c40, b01 = f12 * c20 + f14 * c40 + f13 * c30;
- double b02 = f24 * c40;
- double b10 = f02 * c21 + f04 * c41, b11 = f12 * c21 + f14 * c41 + f13 * c31;
- double b12 = f24 * c41;
- double b20 = f02 * c22 + f04 * c42, b21 = f12 * c22 + f14 * c42 + f13 * c32;
- double b22 = f24 * c42;
- double b40 = f02 * c42 + f04 * c44, b41 = f12 * c42 + f14 * c44 + f13 * c43;
- double b42 = f24 * c44;
- double b30 = f02 * c32 + f04 * c43, b31 = f12 * c32 + f14 * c43 + f13 * c33;
- double b32 = f24 * c43;
+ GPUdoubleCalc b00 = f02 * c20 + f04 * c40, b01 = f12 * c20 + f14 * c40 + f13 * c30;
+ GPUdoubleCalc b02 = f24 * c40;
+ GPUdoubleCalc b10 = f02 * c21 + f04 * c41, b11 = f12 * c21 + f14 * c41 + f13 * c31;
+ GPUdoubleCalc b12 = f24 * c41;
+ GPUdoubleCalc b20 = f02 * c22 + f04 * c42, b21 = f12 * c22 + f14 * c42 + f13 * c32;
+ GPUdoubleCalc b22 = f24 * c42;
+ GPUdoubleCalc b40 = f02 * c42 + f04 * c44, b41 = f12 * c42 + f14 * c44 + f13 * c43;
+ GPUdoubleCalc b42 = f24 * c44;
+ GPUdoubleCalc b30 = f02 * c32 + f04 * c43, b31 = f12 * c32 + f14 * c43 + f13 * c33;
+ GPUdoubleCalc b32 = f24 * c43;
// a = f*b = f*C*ft
- double a00 = f02 * b20 + f04 * b40, a01 = f02 * b21 + f04 * b41, a02 = f02 * b22 + f04 * b42;
- double a11 = f12 * b21 + f14 * b41 + f13 * b31, a12 = f12 * b22 + f14 * b42 + f13 * b32;
- double a22 = f24 * b42;
+ GPUdoubleCalc a00 = f02 * b20 + f04 * b40, a01 = f02 * b21 + f04 * b41, a02 = f02 * b22 + f04 * b42;
+ GPUdoubleCalc a11 = f12 * b21 + f14 * b41 + f13 * b31, a12 = f12 * b22 + f14 * b42 + f13 * b32;
+ GPUdoubleCalc a22 = f24 * b42;
// F*C*Ft = C + (b + bt + a)
c00 += b00 + b00 + a00;
@@ -636,9 +637,9 @@ GPUd() bool TrackParametrizationWithError::propagateTo(value_t xk, cons
if (gpu::CAMath::Abs(r2) < constants::math::Almost0) {
return false;
}
- double r1pr2Inv = 1. / (r1 + r2), r2inv = 1. / r2, r1inv = 1. / r1;
- double dy2dx = (f1 + f2) * r1pr2Inv, dx2r1pr2 = dx * r1pr2Inv;
- value_t step = (gpu::CAMath::Abs(x2r) < 0.05f) ? dx * gpu::CAMath::Abs(r2 + f2 * dy2dx) // chord
+ GPUdoubleCalc r1pr2Inv = 1. / (r1 + r2), r2inv = 1. / r2, r1inv = 1. / r1;
+ GPUdoubleCalc dy2dx = (f1 + f2) * r1pr2Inv, dx2r1pr2 = dx * r1pr2Inv;
+ value_t step = (gpu::CAMath::Abs(x2r) < 0.05f) ? value_t(dx * gpu::CAMath::Abs(r2 + f2 * dy2dx)) // chord
: 2.f * gpu::CAMath::ASin(0.5f * dx * gpu::CAMath::Sqrt(1.f + dy2dx * dy2dx) * crv) / crv; // arc
step *= gpu::CAMath::Sqrt(1.f + this->getTgl() * this->getTgl());
//
@@ -656,30 +657,30 @@ GPUd() bool TrackParametrizationWithError::propagateTo(value_t xk, cons
// evaluate matrix in double prec.
value_t kb = b[2] * constants::math::B2C;
- double hh = dx2r1pr2 * r2inv * (1. + r1 * r2 + f1 * f2), jj = dx * (dy2dx - f2 * r2inv);
- double f02 = hh * r1inv;
- double f04 = hh * dx2r1pr2 * kb;
- double f24 = dx * kb; // x2r/mP[kQ2Pt];
- double f12 = this->getTgl() * (f02 * f2 + jj);
- double f13 = dx * (r2 + f2 * dy2dx);
- double f14 = this->getTgl() * (f04 * f2 + jj * f24);
+ GPUdoubleCalc hh = dx2r1pr2 * r2inv * (1. + r1 * r2 + f1 * f2), jj = dx * (dy2dx - f2 * r2inv);
+ GPUdoubleCalc f02 = hh * r1inv;
+ GPUdoubleCalc f04 = hh * dx2r1pr2 * kb;
+ GPUdoubleCalc f24 = dx * kb; // x2r/mP[kQ2Pt];
+ GPUdoubleCalc f12 = this->getTgl() * (f02 * f2 + jj);
+ GPUdoubleCalc f13 = dx * (r2 + f2 * dy2dx);
+ GPUdoubleCalc f14 = this->getTgl() * (f04 * f2 + jj * f24);
// b = C*ft
- double b00 = f02 * c20 + f04 * c40, b01 = f12 * c20 + f14 * c40 + f13 * c30;
- double b02 = f24 * c40;
- double b10 = f02 * c21 + f04 * c41, b11 = f12 * c21 + f14 * c41 + f13 * c31;
- double b12 = f24 * c41;
- double b20 = f02 * c22 + f04 * c42, b21 = f12 * c22 + f14 * c42 + f13 * c32;
- double b22 = f24 * c42;
- double b40 = f02 * c42 + f04 * c44, b41 = f12 * c42 + f14 * c44 + f13 * c43;
- double b42 = f24 * c44;
- double b30 = f02 * c32 + f04 * c43, b31 = f12 * c32 + f14 * c43 + f13 * c33;
- double b32 = f24 * c43;
+ GPUdoubleCalc b00 = f02 * c20 + f04 * c40, b01 = f12 * c20 + f14 * c40 + f13 * c30;
+ GPUdoubleCalc b02 = f24 * c40;
+ GPUdoubleCalc b10 = f02 * c21 + f04 * c41, b11 = f12 * c21 + f14 * c41 + f13 * c31;
+ GPUdoubleCalc b12 = f24 * c41;
+ GPUdoubleCalc b20 = f02 * c22 + f04 * c42, b21 = f12 * c22 + f14 * c42 + f13 * c32;
+ GPUdoubleCalc b22 = f24 * c42;
+ GPUdoubleCalc b40 = f02 * c42 + f04 * c44, b41 = f12 * c42 + f14 * c44 + f13 * c43;
+ GPUdoubleCalc b42 = f24 * c44;
+ GPUdoubleCalc b30 = f02 * c32 + f04 * c43, b31 = f12 * c32 + f14 * c43 + f13 * c33;
+ GPUdoubleCalc b32 = f24 * c43;
// a = f*b = f*C*ft
- double a00 = f02 * b20 + f04 * b40, a01 = f02 * b21 + f04 * b41, a02 = f02 * b22 + f04 * b42;
- double a11 = f12 * b21 + f14 * b41 + f13 * b31, a12 = f12 * b22 + f14 * b42 + f13 * b32;
- double a22 = f24 * b42;
+ GPUdoubleCalc a00 = f02 * b20 + f04 * b40, a01 = f02 * b21 + f04 * b41, a02 = f02 * b22 + f04 * b42;
+ GPUdoubleCalc a11 = f12 * b21 + f14 * b41 + f13 * b31, a12 = f12 * b22 + f14 * b42 + f13 * b32;
+ GPUdoubleCalc a22 = f24 * b42;
// F*C*Ft = C + (b + bt + a)
c00 += b00 + b00 + a00;
@@ -889,13 +890,13 @@ GPUd() bool TrackParametrizationWithError::propagateTo(value_t xk, Trac
cc = cspRef0 + cspRef1;
ccInv = value_t(1) / cc;
dy2dx = (snpRef0 + snpRef1) * ccInv;
- double dxccInv = dx * ccInv, hh = dxccInv * cspRef1Inv * (1 + cspRef0 * cspRef1 + snpRef0 * snpRef1), jj = dx * (dy2dx - snpRef1 * cspRef1Inv);
- double f02 = hh * cspRef0Inv;
- double f04 = hh * dxccInv * kb;
- double f24 = dx * kb;
- double f12 = linRef0.getTgl() * (f02 * snpRef1 + jj);
- double f13 = dx * (cspRef1 + snpRef1 * dy2dx); // dS
- double f14 = linRef0.getTgl() * (f04 * snpRef1 + jj * f24);
+ GPUdoubleCalc dxccInv = dx * ccInv, hh = dxccInv * cspRef1Inv * (1 + cspRef0 * cspRef1 + snpRef0 * snpRef1), jj = dx * (dy2dx - snpRef1 * cspRef1Inv);
+ GPUdoubleCalc f02 = hh * cspRef0Inv;
+ GPUdoubleCalc f04 = hh * dxccInv * kb;
+ GPUdoubleCalc f24 = dx * kb;
+ GPUdoubleCalc f12 = linRef0.getTgl() * (f02 * snpRef1 + jj);
+ GPUdoubleCalc f13 = dx * (cspRef1 + snpRef1 * dy2dx); // dS
+ GPUdoubleCalc f14 = linRef0.getTgl() * (f04 * snpRef1 + jj * f24);
// difference between the current and reference state
value_t diff[5];
@@ -922,21 +923,21 @@ GPUd() bool TrackParametrizationWithError::propagateTo(value_t xk, Trac
&c44 = mC[kSigQ2Pt2];
// b = C*ft
- double b00 = f02 * c20 + f04 * c40, b01 = f12 * c20 + f14 * c40 + f13 * c30;
- double b02 = f24 * c40;
- double b10 = f02 * c21 + f04 * c41, b11 = f12 * c21 + f14 * c41 + f13 * c31;
- double b12 = f24 * c41;
- double b20 = f02 * c22 + f04 * c42, b21 = f12 * c22 + f14 * c42 + f13 * c32;
- double b22 = f24 * c42;
- double b40 = f02 * c42 + f04 * c44, b41 = f12 * c42 + f14 * c44 + f13 * c43;
- double b42 = f24 * c44;
- double b30 = f02 * c32 + f04 * c43, b31 = f12 * c32 + f14 * c43 + f13 * c33;
- double b32 = f24 * c43;
+ GPUdoubleCalc b00 = f02 * c20 + f04 * c40, b01 = f12 * c20 + f14 * c40 + f13 * c30;
+ GPUdoubleCalc b02 = f24 * c40;
+ GPUdoubleCalc b10 = f02 * c21 + f04 * c41, b11 = f12 * c21 + f14 * c41 + f13 * c31;
+ GPUdoubleCalc b12 = f24 * c41;
+ GPUdoubleCalc b20 = f02 * c22 + f04 * c42, b21 = f12 * c22 + f14 * c42 + f13 * c32;
+ GPUdoubleCalc b22 = f24 * c42;
+ GPUdoubleCalc b40 = f02 * c42 + f04 * c44, b41 = f12 * c42 + f14 * c44 + f13 * c43;
+ GPUdoubleCalc b42 = f24 * c44;
+ GPUdoubleCalc b30 = f02 * c32 + f04 * c43, b31 = f12 * c32 + f14 * c43 + f13 * c33;
+ GPUdoubleCalc b32 = f24 * c43;
// a = f*b = f*C*ft
- double a00 = f02 * b20 + f04 * b40, a01 = f02 * b21 + f04 * b41, a02 = f02 * b22 + f04 * b42;
- double a11 = f12 * b21 + f14 * b41 + f13 * b31, a12 = f12 * b22 + f14 * b42 + f13 * b32;
- double a22 = f24 * b42;
+ GPUdoubleCalc a00 = f02 * b20 + f04 * b40, a01 = f02 * b21 + f04 * b41, a02 = f02 * b22 + f04 * b42;
+ GPUdoubleCalc a11 = f12 * b21 + f14 * b41 + f13 * b31, a12 = f12 * b22 + f14 * b42 + f13 * b32;
+ GPUdoubleCalc a22 = f24 * b42;
// F*C*Ft = C + (b + bt + a)
c00 += b00 + b00 + a00;
@@ -1034,7 +1035,7 @@ template
GPUd() void TrackParametrizationWithError::resetCovariance(value_t s2)
{
// Reset the covarince matrix to "something big"
- double d0(kCY2max), d1(kCZ2max), d2(kCSnp2max), d3(kCTgl2max), d4(kC1Pt2max);
+ GPUdoubleCalc d0(kCY2max), d1(kCZ2max), d2(kCSnp2max), d3(kCTgl2max), d4(kC1Pt2max);
if (s2 > constants::math::Almost0) {
d0 = getSigmaY2() * s2;
d1 = getSigmaZ2() * s2;
@@ -1072,9 +1073,9 @@ template
GPUd() auto TrackParametrizationWithError::getPredictedChi2(const value_t* p, const value_t* cov) const -> value_t
{
// Estimate the chi2 of the space point "p" with the cov. matrix "cov"
- auto sdd = static_cast(getSigmaY2()) + static_cast(cov[0]);
- auto sdz = static_cast(getSigmaZY()) + static_cast(cov[1]);
- auto szz = static_cast(getSigmaZ2()) + static_cast(cov[2]);
+ auto sdd = static_cast(getSigmaY2()) + static_cast(cov[0]);
+ auto sdz = static_cast(getSigmaZY()) + static_cast(cov[1]);
+ auto szz = static_cast(getSigmaZ2()) + static_cast(cov[2]);
auto det = sdd * szz - sdz * sdz;
if (gpu::CAMath::Abs(det) < constants::math::Almost0) {
@@ -1086,7 +1087,7 @@ GPUd() auto TrackParametrizationWithError::getPredictedChi2(const value
auto chi2 = (d * (szz * d - sdz * z) + z * (sdd * z - d * sdz)) / det;
if (chi2 < 0.) {
#ifndef GPUCA_ALIGPUCODE
- LOGP(warning, "Negative chi2={}, Cluster: {} {} {} Dy:{} Dz:{} | sdd:{} sdz:{} szz:{} det:{}", chi2, cov[0], cov[1], cov[2], d, z, sdd, sdz, szz, det);
+ LOGP(warning, "Negative chi2={}, Cluster: {} {} {} Dy:{} Dz:{} | sdd:{} sdz:{} szz:{} det:{}", double(chi2), cov[0], cov[1], cov[2], d, z, double(sdd), double(sdz), double(szz), double(det));
LOGP(warning, "Track: {}", asString());
#endif
}
@@ -1098,9 +1099,9 @@ template
GPUd() auto TrackParametrizationWithError::getPredictedChi2Quiet(const value_t* p, const value_t* cov) const -> value_t
{
// Estimate the chi2 of the space point "p" with the cov. matrix "cov"
- auto sdd = static_cast(getSigmaY2()) + static_cast(cov[0]);
- auto sdz = static_cast(getSigmaZY()) + static_cast(cov[1]);
- auto szz = static_cast(getSigmaZ2()) + static_cast(cov[2]);
+ auto sdd = static_cast(getSigmaY2()) + static_cast(cov[0]);
+ auto sdz = static_cast(getSigmaZY()) + static_cast(cov[1]);
+ auto szz = static_cast(getSigmaZ2()) + static_cast(cov[2]);
auto det = sdd * szz - sdz * sdz;
if (gpu::CAMath::Abs(det) < constants::math::Almost0) {
@@ -1150,9 +1151,9 @@ GPUd() auto TrackParametrizationWithError::getPredictedChi2Fast(const T
// Factorize cov = L * D * L^T with L unit lower triangular. The strictly lower triangle of
// lmat holds L, its strictly upper triangle holds the transpose of L * D, so that the inner
// products below need no extra multiplication by D. dInv holds the inverted diagonal of D.
- double lmat[kNParams][kNParams], dInv[kNParams];
+ GPUdoubleCalc lmat[kNParams][kNParams], dInv[kNParams];
for (int j = 0; j < kNParams; j++) {
- double djj = cov(j, j);
+ GPUdoubleCalc djj = cov(j, j);
for (int k = 0; k < j; k++) {
djj -= lmat[j][k] * lmat[k][j];
}
@@ -1161,7 +1162,7 @@ GPUd() auto TrackParametrizationWithError::getPredictedChi2Fast(const T
}
dInv[j] = 1. / djj;
for (int i = j + 1; i < kNParams; i++) {
- double s = cov(i, j);
+ GPUdoubleCalc s = cov(i, j);
for (int k = 0; k < j; k++) {
s -= lmat[i][k] * lmat[k][j];
}
@@ -1171,9 +1172,9 @@ GPUd() auto TrackParametrizationWithError::getPredictedChi2Fast(const T
}
// chi2 = d^T C^-1 d = sum_i y_i^2 / D_i with y from the forward substitution L y = d
- double chi2 = 0., y[kNParams];
+ GPUdoubleCalc chi2 = 0., y[kNParams];
for (int i = 0; i < kNParams; i++) {
- double s = double(this->getParam(i)) - double(rhs.getParam(i));
+ GPUdoubleCalc s = GPUdoubleCalc(this->getParam(i)) - GPUdoubleCalc(rhs.getParam(i));
for (int k = 0; k < i; k++) {
s -= lmat[i][k] * y[k];
}
@@ -1188,21 +1189,21 @@ template
GPUd() void TrackParametrizationWithError::buildCombinedCovMatrix(const TrackParametrizationWithError& rhs, MatrixDSym5& cov) const
{
// fill combined cov.matrix (NOT inverted)
- cov(kY, kY) = static_cast(getSigmaY2()) + static_cast(rhs.getSigmaY2());
- cov(kZ, kY) = static_cast(getSigmaZY()) + static_cast(rhs.getSigmaZY());
- cov(kZ, kZ) = static_cast(getSigmaZ2()) + static_cast(rhs.getSigmaZ2());
- cov(kSnp, kY) = static_cast(getSigmaSnpY()) + static_cast(rhs.getSigmaSnpY());
- cov(kSnp, kZ) = static_cast(getSigmaSnpZ()) + static_cast(rhs.getSigmaSnpZ());
- cov(kSnp, kSnp) = static_cast(getSigmaSnp2()) + static_cast(rhs.getSigmaSnp2());
- cov(kTgl, kY) = static_cast(getSigmaTglY()) + static_cast(rhs.getSigmaTglY());
- cov(kTgl, kZ) = static_cast(getSigmaTglZ()) + static_cast(rhs.getSigmaTglZ());
- cov(kTgl, kSnp) = static_cast(getSigmaTglSnp()) + static_cast(rhs.getSigmaTglSnp());
- cov(kTgl, kTgl) = static_cast(getSigmaTgl2()) + static_cast(rhs.getSigmaTgl2());
- cov(kQ2Pt, kY) = static_cast(getSigma1PtY()) + static_cast(rhs.getSigma1PtY());
- cov(kQ2Pt, kZ) = static_cast(getSigma1PtZ()) + static_cast(rhs.getSigma1PtZ());
- cov(kQ2Pt, kSnp) = static_cast(getSigma1PtSnp()) + static_cast(rhs.getSigma1PtSnp());
- cov(kQ2Pt, kTgl) = static_cast(getSigma1PtTgl()) + static_cast(rhs.getSigma1PtTgl());
- cov(kQ2Pt, kQ2Pt) = static_cast(getSigma1Pt2()) + static_cast(rhs.getSigma1Pt2());
+ cov(kY, kY) = static_cast(getSigmaY2()) + static_cast(rhs.getSigmaY2());
+ cov(kZ, kY) = static_cast(getSigmaZY()) + static_cast(rhs.getSigmaZY());
+ cov(kZ, kZ) = static_cast(getSigmaZ2()) + static_cast(rhs.getSigmaZ2());
+ cov(kSnp, kY) = static_cast(getSigmaSnpY()) + static_cast(rhs.getSigmaSnpY());
+ cov(kSnp, kZ) = static_cast(getSigmaSnpZ()) + static_cast(rhs.getSigmaSnpZ());
+ cov(kSnp, kSnp) = static_cast(getSigmaSnp2()) + static_cast(rhs.getSigmaSnp2());
+ cov(kTgl, kY) = static_cast(getSigmaTglY()) + static_cast(rhs.getSigmaTglY());
+ cov(kTgl, kZ) = static_cast(getSigmaTglZ()) + static_cast(rhs.getSigmaTglZ());
+ cov(kTgl, kSnp) = static_cast(getSigmaTglSnp()) + static_cast(rhs.getSigmaTglSnp());
+ cov(kTgl, kTgl) = static_cast(getSigmaTgl2()) + static_cast(rhs.getSigmaTgl2());
+ cov(kQ2Pt, kY) = static_cast(getSigma1PtY()) + static_cast(rhs.getSigma1PtY());
+ cov(kQ2Pt, kZ) = static_cast(getSigma1PtZ()) + static_cast(rhs.getSigma1PtZ());
+ cov(kQ2Pt, kSnp) = static_cast(getSigma1PtSnp()) + static_cast(rhs.getSigma1PtSnp());
+ cov(kQ2Pt, kTgl) = static_cast