diff --git a/eic_data.yaml b/eic_data.yaml
index 7932b84e9c481ee339af4b108f32b1fc646626e0..3d64601e6146b74adc32941f2ecb34405aa703ac 100644
--- a/eic_data.yaml
+++ b/eic_data.yaml
@@ -45,39 +45,6 @@ components:
explicit operator bool() const {return valid();}
"
- ## simple weight that defaults to 1 if not set
- eic::Weight:
- Members:
- - float value
- ExtraCode:
- declaration: "
- Weight() : value{1.} {}\n
- Weight(double w) : value {static_cast<float>(w)} {}\n
- Weight& operator=(double w) {value = static_cast<float>(w); return *this;}\n
- operator float() const {return value;}
- "
-
- eic::VectorXY:
- Members:
- - float x // [mm] or [GeV]
- - float y //
- ExtraCode:
- includes: "#include <cmath>\n"
- declaration: "
- VectorXY() : x{0}, y{0} {}\n
- VectorXY(double xx, double yy) : x{static_cast<float>(xx)}, y{static_cast<float>(yy)} {}\n
- float& operator[](unsigned i) {return *(&x + i);}\n
- const float& operator[](unsigned i) const {return *(&x + i);}\n
- float mag() const {return std::hypot(x, y);}\n
- float r() const {return mag();}\n
- float phi() const {return atan2(y, x);}\n
- operator std::pair<float, float>() const {return {x, y};}\n
- float dot(const VectorXY& rhs) const {return x*rhs.x + y*rhs.y;}\n
- VectorXY add(const VectorXY& rhs) const {return {x+rhs.x, y+rhs.y};}\n
- VectorXY subtract(const VectorXY& rhs) const {return {x-rhs.x, y-rhs.y};}\n
- VectorXY scale(double f) const {return {f*x, f*y};}\n
- "
-
# Copy from EDM4hep with conversion operators added
# TODO: uncomment Vector2f conversions when we migrate to new EDM4hep version
eicd::Vector2f :
@@ -117,94 +84,47 @@ components:
includes: "
#include <edm4hep/Vector3f.h>\n
"
-
- eic::VectorXYZ:
- Members:
- - float x // [mm] or [GeV]
- - float y //
- - float z //
- ExtraCode:
- includes: "#include <cmath>\n#include<tuple>"
- declaration: "
- VectorXYZ() : x{0}, y{0}, z{0} {}\n
- VectorXYZ(double xx, double yy, double zz) : x{static_cast<float>(xx)}, y{static_cast<float>(yy)}, z{static_cast<float>(zz)} {}\n
- template <class V>
- VectorXYZ(const V& v) : x{v.x}, y{v.y}, z{v.z} {}\n
- float& operator[](unsigned i) {return *(&x + i);}\n
- const float& operator[](unsigned i) const {return *(&x + i);}\n
- float mag() const {return std::hypot(x, y, z);}\n
- float r() const {return mag();}\n
- float theta() const {return acos(z/mag());}\n
- float phi() const {return atan2(y,x);}\n
- float eta() const {return -log(tan(0.5*theta()));}\n
- operator std::tuple<float, float, float>() {return {x, y, z};}\n
- float dot(const VectorXYZ& rhs) const {return x*rhs.x + y*rhs.y + z*rhs.z;}\n
- VectorXYZ add(const VectorXYZ& rhs) const {return {x+rhs.x, y+rhs.y, z+rhs.z};}\n
- VectorXYZ subtract(const VectorXYZ& rhs) const {return {x-rhs.x, y-rhs.y, z-rhs.z};}\n
- VectorXYZ scale(double f) const {return {f*x, f*y, f*z};}\n
- "
-
- eic::VectorXYZT:
+ # Copy from EDM4hep with conversion operators added
+ eicd::Vector3d :
Members:
- - double x // [mm] or [GeV]
- - double y //
- - double z //
- - double t // [ns] or [GeV]
+ - double x
+ - double y
+ - double z
ExtraCode:
- includes: "#include <cmath>\n#include <tuple>"
declaration: "
- VectorXYZT() : x{0}, y{0}, z{0}, t{0} {}\n
- VectorXYZT(double xx, double yy, double zz, double tt) : x{xx}, y{yy}, z{zz}, t{tt} {}\n
- double& operator[](unsigned i) {return *(&x + i);}\n
- const double& operator[](unsigned i) const {return *(&x + i);}\n
- double mag() const {return std::hypot(x, y, z);}\n
- double r() const {return mag();}\n
- double theta() const {return acos(z/mag());}\n
- double phi() const {return atan2(y,x);}\n
- double eta() const {return -log(tan(0.5*theta()));}
- double energy() const {return t;}\n
- double mass() const {return sqrt(t*t - x*x - y*y - z*z);}\n
- operator std::tuple<double, double, double, double>() {return {x, y, z, t};}\n
- double dot(const VectorXYZT& rhs) const {return t*rhs.t - x*rhs.x - y*rhs.y - z*rhs.z;}\n
- VectorXYZT add(const VectorXYZT& rhs) const {return {x+rhs.x, y+rhs.y, z+rhs.z, t+rhs.t};}\n
- VectorXYZT subtract(const VectorXYZT& rhs) const {return {x-rhs.x, y-rhs.y, z-rhs.z, t-rhs.t};}\n
- VectorXYZT scale(double f) const {return {f*x, f*y, f*z, f*t};}\n
+ Vector3d() : x(0),y(0),z(0) {}\n
+ Vector3d(double xx, double yy, double zz) : x(xx),y(yy),z(zz) {}\n
+ Vector3d(const double* v) : x(v[0]),y(v[1]),z(v[2]) {}\n
+ bool operator==(const Vector3d& v) const { return (x==v.x&&y==v.y&&z==v.z) ; }\n
+ double operator[](unsigned i) const { return *( &x + i ) ; }\n
+ Vector3d(const edm4hep::Vector3d& v) : x{v.x}, y{v.y}, z{v.z} {}
+ operator edm4hep::Vector3d() const {return {x, y, z};}\n
+ "
+ includes: "
+ #include <edm4hep/Vector3d.h>\n
"
- eic::CovDiagXYZ:
+ eicd::CovDiag3f:
Members:
- float xx
- float yy
- float zz
ExtraCode:
declaration: "
- CovDiagXYZ() : xx{0}, yy{0}, zz{0} {}\n
- CovDiagXYZ(double x, double y, double z) : xx{static_cast<float>(x)}, yy{static_cast<float>(y)}, zz{static_cast<float>(z)} {}\n
+ CovDiag3f() : xx{0}, yy{0}, zz{0} {}\n
+ CovDiag3f(double x, double y, double z) : xx{static_cast<float>(x)}, yy{static_cast<float>(y)}, zz{static_cast<float>(z)} {}\n
float operator()(unsigned i, unsigned j) const {return (i == j) ? *(&xx + i) : 0.;}\n
"
- eic::CovDiagXYZT:
- Members:
- - double xx
- - double yy
- - double zz
- - double tt
- ExtraCode:
- declaration: "
- CovDiagXYZT() : xx{0}, yy{0}, zz{0}, tt{0} {}\n
- CovDiagXYZT(double x, double y, double z, double t) : xx{x}, yy{y}, zz{z}, tt{t} {}\n
- double operator()(unsigned i, unsigned j) const {return (i == j) ? *(&xx + i) : 0.;}\n
- "
-
- eic::CovXY:
+ eicd::Cov2f:
Members:
- float xx
- float yy
- float xy
ExtraCode:
declaration: "
- CovXY() : xx{0}, yy{0}, xy{0} {}\n
- CovXY(double vx, double vy, double vxy = 0)\n
+ Cov2f() : xx{0}, yy{0}, xy{0} {}\n
+ Cov2f(double vx, double vy, double vxy = 0)\n
: xx{static_cast<float>(vx)}, yy{static_cast<float>(vy)}, xy{static_cast<float>(vxy)} {}\n
float operator()(unsigned i, unsigned j) const {\n
// diagonal\n
@@ -222,9 +142,9 @@ components:
Members:
- eic::Index ID // ID of the hit
- uint32_t index // Raw index of the hit in the relevant array
- - eic::Weight weight // weight of the hit
+ - float weight // weight of the hit
- eic::CovXYZ:
+ eicd::Cov3f:
Members:
- float xx
- float yy
@@ -234,12 +154,12 @@ components:
- float yz
ExtraCode:
declaration: "
- CovXYZ() : xx{0}, yy{0}, zz{0}, xy{0}, xz{0}, yz{0} {}\n
- CovXYZ(double vx, double vy, double vz, double vxy, double vxz, double vyz)\n
+ Cov3f() : xx{0}, yy{0}, zz{0}, xy{0}, xz{0}, yz{0} {}\n
+ Cov3f(double vx, double vy, double vz, double vxy, double vxz, double vyz)\n
: xx{static_cast<float>(vx)}, yy{static_cast<float>(vy)}, zz{static_cast<float>(vz)},\n
xy{static_cast<float>(vxy)}, xz{static_cast<float>(vxz)}, yz{static_cast<float>(vyz)} {}\n
- CovXYZ(double vx, double vy, double vz)\n
- : CovXYZ(vx, vy, vz, 0, 0, 0) {} \n
+ Cov3f(double vx, double vy, double vz)\n
+ : Cov3f(vx, vy, vz, 0, 0, 0) {} \n
float operator()(unsigned i, unsigned j) const {\n
// diagonal\n
if (i == j) {\n
@@ -255,15 +175,15 @@ components:
## A point along a track
eic::TrackPoint:
Members:
- - eic::VectorXYZ position // Position of the trajectory point [mm]
- - eic::CovXYZ positionError // Error on the position
- - eic::VectorXYZ momentum // 3-momentum at the point [GeV]
- - eic::CovXYZ momentumError // Error on the 3-momentum
+ - eicd::Vector3f position // Position of the trajectory point [mm]
+ - eicd::Cov3f positionError // Error on the position
+ - eicd::Vector3f momentum // 3-momentum at the point [GeV]
+ - eicd::Cov3f momentumError // Error on the 3-momentum
- float time // Time at this point [ns]
- float timeError // Error on the time at this point
- float theta // polar direction of the track at the surface [rad]
- float phi // azimuthal direction of the track at the surface [rad]
- - eic::CovXY directionError // Error on the polar and azimuthal angles
+ - eicd::Cov2f directionError // Error on the polar and azimuthal angles
- float pathlength // Pathlength from the origin to this point
- float pathlengthError // Error on the pathlenght
@@ -282,7 +202,7 @@ datatypes:
- int32_t type // Event type identifier (TBD).
- int32_t proc // Process identifier (TBD).
- int32_t source // Source/identifier (TBD)
- - eic::Weight weight // Optional event weight (useful for MC)
+ - float weight // Optional event weight (useful for MC)
## ==========================================================================
@@ -295,25 +215,25 @@ datatypes:
Members:
- eic::Index ID // Unique particle index
- eicd::Vector3f p // Momentum [GeV]
- - eicd::Vector3f v // Vertex [mm]
+ - eicd::Vector3d v // Vertex [mm]
- float time // Time in [ns]
- int32_t pid // Particle PDG code
- int16_t status // Status code
- int16_t charge // Particle charge (or sign)
- - eic::Weight weight // Particle weight, e.g. from PID algorithm [0-1]
+ - float weight // Particle weight, e.g. from PID algorithm [0-1]
eic::ReconstructedParticle:
Description: "EIC Reconstructed Particle"
Author: "W. Armstrong, S. Joosten"
Members:
- eic::Index ID // Unique particle index
- - eic::VectorXYZ p // Momentum vector [GeV]
- - eic::VectorXYZ v // Vertex [mm]
+ - eicd::Vector3f p // Momentum vector [GeV]
+ - eicd::Vector3d v // Vertex [mm]
- float time // Time in [ns]
- int32_t pid // PID of reconstructed particle.
- int16_t status // Status code
- int16_t charge // Particle charge (or sign)
- - eic::Weight weight // Particle weight, e.g. from PID algorithm [0-1]
+ - float weight // Particle weight, e.g. from PID algorithm [0-1]
- float theta // Polar angle of this particle [rad]
- float phi // Azimuthal angle of this particle [rad]
- float momentum // particle 3-momentum magnitude [GeV]
@@ -347,11 +267,11 @@ datatypes:
- float energyError // Error on energy [GeV].
- float time // The time of the hit in [ns].
- float timeError // Error on the time
- - eic::VectorXYZ position // The global position of the hit in world coordinates [mm].
- - eic::VectorXYZ dimension // The dimension information of the cell [mm].
+ - eicd::Vector3f position // The global position of the hit in world coordinates [mm].
+ - eicd::Vector3f dimension // The dimension information of the cell [mm].
- int32_t sector // Sector that this hit occured in
- int32_t layer // Layer that the hit occured in
- - eic::VectorXYZ local // The local coordinates of the hit in the detector segment [mm].
+ - eicd::Vector3f local // The local coordinates of the hit in the detector segment [mm].
## ==========================================================================
## Clustering
@@ -376,11 +296,11 @@ datatypes:
- float time // [ns]
- float timeError // Error on the cluster time
- uint32_t nhits // Number of hits in the cluster.
- - eic::VectorXYZ position // Global position of the cluster [mm].
- - eic::CovXYZ positionError // Covariance matrix of the position (6 Parameters).
+ - eicd::Vector3f position // Global position of the cluster [mm].
+ - eicd::Cov3f positionError // Covariance matrix of the position (6 Parameters).
- float intrinsicTheta // Intrinsic cluster propagation direction polar angle [rad]
- float intrinsicPhi // Intrinsic cluster propagation direction azimuthal angle [rad]
- - eic::CovXY intrinsicDirectionError // Error on the intrinsic cluster propagation direction
+ - eicd::Cov2f intrinsicDirectionError // Error on the intrinsic cluster propagation direction
VectorMembers:
- float shapeParameters // Should be set in metadata, for now radius/skewness
- float hitContributions // Energy contributions of the hits. Runs parallel to ::hits()
@@ -412,9 +332,9 @@ datatypes:
- float npe // Estimated number of photo-electrons [#]
- float time // Time [ns]
- float timeError // Error on the time [ns]
- - eic::VectorXYZ position // PMT hit position [mm]
- - eic::VectorXYZ local // The local position of the hit in detector coordinates [mm]
- - eic::VectorXYZ dimension // The dimension information of the pixel [mm].
+ - eicd::Vector3f position // PMT hit position [mm]
+ - eicd::Vector3f local // The local position of the hit in detector coordinates [mm]
+ - eicd::Vector3f dimension // The dimension information of the pixel [mm].
eic::RingImage:
Description: "EIC Ring Image Cluster"
@@ -422,8 +342,8 @@ datatypes:
Members:
- eic::Index ID // Unique cluster ID
- float npe // Number of photo-electrons [#]
- - eic::VectorXYZ position // Global position of the cluster [mm]
- - eic::VectorXYZ positionError // Error on the position
+ - eicd::Vector3f position // Global position of the cluster [mm]
+ - eicd::Vector3f positionError // Error on the position
- float theta // Opening angle of the ring [rad, 0->pi]
- float thetaError // Error on the opening angle
- float radius // Radius of the best fit ring [mm]
@@ -448,8 +368,8 @@ datatypes:
Members:
- eic::Index ID // Unique hit ID, same as one of the involved raw hits.
- uint64_t cellID // The detector specific (geometrical) cell id.
- - eic::VectorXYZ position // Hit (cell) position and time [mm, ns]
- - eic::CovDiagXYZ positionError // Covariance Matrix
+ - eicd::Vector3f position // Hit (cell) position and time [mm, ns]
+ - eicd::CovDiag3f positionError // Covariance Matrix
- float time // Hit time
- float timeError // Error on the time
- float edep // Energy deposit in this hit [GeV]
@@ -464,7 +384,7 @@ datatypes:
## Members:
## eic::Index ID // Unique identifier
## eic::Index seedID // Index of corresponding initial track parameters
- ## eic::Weight weight // prototrack weight, in case we share pixels [0-1]
+ ## float weight // prototrack weight, in case we share pixels [0-1]
## VectorMembers:
## int32_t hits // tracker hit indicies
@@ -492,12 +412,12 @@ datatypes:
Author: "W. Armstrong, S. Joosten"
Members:
- int32_t type // Type of track parameters (-1/seed, 0/head, ...)
- - eic::VectorXY loc // 2D location on surface
- - eic::CovXY locError // Covariance on loc
+ - eicd::Vector2f loc // 2D location on surface
+ - eicd::Cov2f locError // Covariance on loc
- float theta // Track polar angle [rad]
- float phi // Track azimuthal angle [rad]
- float qOverP // [e/GeV]
- - eic::CovXYZ momentumError // Covariance on theta, phi and qOverP
+ - eicd::Cov3f momentumError // Covariance on theta, phi and qOverP
- float time // Track time [ns]
- float timeError // Error on the time
- float charge // Particle charge
@@ -511,8 +431,8 @@ datatypes:
- int32_t type // Flag that defines the type of track
- float chi2 // Total chi2 (sum) of the track fit
- int32_t ndf // Numbers of degrees of freedom of the track fit
- - eic::VectorXYZ momentum // Track 3-momentum at the vertex [GeV]
- - eic::CovXYZ momentumError // Covariance matrix on the momentum
+ - eicd::Vector3f momentum // Track 3-momentum at the vertex [GeV]
+ - eicd::Cov3f momentumError // Covariance matrix on the momentum
- float time // Track time at the vertex [ns]
- float timeError // Error on the track vertex time
- float charge // Particle charge
@@ -541,7 +461,7 @@ datatypes:
Author: "W. Armstrong, S. Joosten"
Members:
- eic::Index ID // Unique vertex ID
- - eic::VectorXYZ position // Postion of vertex [mm]
+ - eicd::Vector3f position // Postion of vertex [mm]
- float time // Time of vertex [ns]
- float chi2 // Chi squared of the vertex fit.
- float probability // Probability of the vertex fit
diff --git a/utils/include/eicd/vector_utils.h b/utils/include/eicd/vector_utils.h
index 1c816282ece9f3b1c51447c2f3aeb87103747fb5..96c80abf7b994e201d8c48033dd7d173bdd4f190 100644
--- a/utils/include/eicd/vector_utils.h
+++ b/utils/include/eicd/vector_utils.h
@@ -21,14 +21,21 @@ namespace eicd {
template <class V> concept VectorHasX = requires(V v) { v.x; };
template <class V> concept VectorHasY = requires(V v) { v.y; };
template <class V> concept VectorHasZ = requires(V v) { v.z; };
+template <class V> concept VectorHasA = requires(V v) { v.a; };
+template <class V> concept VectorHasB = requires(V v) { v.b; };
template <class V> concept ClassVector = requires(V v) { v.x(); };
template <class V>
-concept Vector2D =
+concept Vector2D_XY =
VectorHasX<V>&& VectorHasY<V> && !VectorHasZ<V> && !ClassVector<V>;
template <class V>
+concept Vector2D_AB =
+ VectorHasA<V>&& VectorHasB<V> && !VectorHasZ<V> && !ClassVector<V>;
+template <class V> concept Vector2D = Vector2D_XY<V> || Vector2D_AB<V>;
+template <class V>
concept Vector3D =
VectorHasX<V>&& VectorHasY<V>&& VectorHasZ<V> && !ClassVector<V>;
template <class V> concept VectorND = Vector2D<V> || Vector3D<V>;
+template <class V> concept VectorND_XYZ = Vector2D_XY<V> || Vector3D<V>;
inline double etaToAngle(const double eta) {
return std::atan(std::exp(-eta)) * 2.;
@@ -37,13 +44,20 @@ inline double angleToEta(const double theta) {
return -std::log(std::tan(0.5 * theta));
}
+// Utility getters to accomodate different vector types
+template <VectorND_XYZ V> auto vector_x(const V& v) { return v.x; }
+template <VectorND_XYZ V> auto vector_y(const V& v) { return v.y; }
+template <Vector3D V> auto vector_z(const V& v) { return v.z; }
+template <Vector2D_AB V> auto vector_x(const V& v) { return v.a; }
+template <Vector2D_AB V> auto vector_y(const V& v) { return v.b; }
+
// inline edm4hep::Vector2f VectorFromPolar(const double r, const double theta)
// {
// return {r * sin(theta), r * cos(theta)};
//}
template <Vector3D V = edm4hep::Vector3f>
V sphericalToVector(const double r, const double theta, const double phi) {
- using FloatType = decltype(V().x);
+ using FloatType = decltype(vector_x(V()));
const double sth = sin(theta);
const double cth = cos(theta);
const double sph = sin(phi);
@@ -55,24 +69,26 @@ V sphericalToVector(const double r, const double theta, const double phi) {
}
template <Vector3D V> double anglePolar(const V& v) {
- return std::atan2(std::hypot(v.x, v.y), v.z);
+ return std::atan2(std::hypot(vector_x(v), vector_y(v)), vector_z(v));
}
template <VectorND V> double angleAzimuthal(const V& v) {
- return std::atan2(v.y, v.x);
+ return std::atan2(vector_y(v), vector_x(v));
}
template <Vector3D V> double eta(const V& v) {
return angleToEta(anglePolar(v));
}
template <Vector2D V> double magnitude(const V& v) {
- return std::hypot(v.x, v.y);
+ return std::hypot(vector_x(v), vector_y(v));
}
template <Vector3D V> double magnitude(const V& v) {
- return std::hypot(v.x, v.y, v.z);
+ return std::hypot(vector_x(v), vector_y(v), vector_z(v));
}
template <Vector3D V> double magnitudeTransverse(const V& v) {
- return std::hypot(v.x, v.y);
+ return std::hypot(vector_x(v), vector_y(v));
+}
+template <Vector3D V> double magnitudeLongitudinal(const V& v) {
+ return vector_z(v);
}
-template <Vector3D V> double magnitudeLongitudinal(const V& v) { return v.z; }
template <VectorND V> V normalizeVector(const V& v, double norm = 1.) {
const double old = magnitude(v);
if (old == 0) {
@@ -80,8 +96,12 @@ template <VectorND V> V normalizeVector(const V& v, double norm = 1.) {
}
return (norm > 0) ? v * norm / old : v * 0;
}
-template <Vector3D V> V vectorTransverse(const V& v) { return {v.x, v.y, 0}; }
-template <Vector3D V> V vectorLongitudinal(const V& v) { return {0, 0, v.z}; }
+template <Vector3D V> V vectorTransverse(const V& v) {
+ return {vector_x(v), vector_y(v), 0};
+}
+template <Vector3D V> V vectorLongitudinal(const V& v) {
+ return {0, 0, vector_z(v)};
+}
// Two vector functions
template <VectorND V> double angleBetween(const V& v1, const V& v2) {
const double dot = v1 * v2;
@@ -101,28 +121,30 @@ template <Vector3D V> double projection(const V& v, const V& v1) {
} // namespace eicd
template <eicd::Vector2D V> V operator+(const V& v1, const V& v2) {
- return {v1.x + v2.x, v1.y + v2.y};
+ return {vector_x(v1) + vector_x(v2), vector_y(v1) + vector_y(v2)};
}
template <eicd::Vector3D V> V operator+(const V& v1, const V& v2) {
- return {v1.x + v2.x, v1.y + v2.y, v1.z + v2.z};
+ return {vector_x(v1) + vector_x(v2), vector_y(v1) + vector_y(v2),
+ vector_z(v1) + vector_z(v2)};
}
template <eicd::Vector2D V> double operator*(const V& v1, const V& v2) {
- return v1.x * v2.x + v1.y * v2.y;
+ return vector_x(v1) * vector_x(v2) + vector_y(v1) * vector_y(v2);
}
template <eicd::Vector3D V> double operator*(const V& v1, const V& v2) {
- return v1.x * v2.x + v1.y * v2.y + v1.z * v2.z;
+ return vector_x(v1) * vector_x(v2) + vector_y(v1) * vector_y(v2) +
+ vector_z(v1) * vector_z(v2);
}
template <eicd::Vector2D V> V operator*(const double d, const V& v) {
- return {d * v.x, d * v.y};
+ return {d * vector_x(v), d * vector_y(v)};
}
template <eicd::Vector3D V> V operator*(const double d, const V& v) {
- return {d * v.x, d * v.y, d * v.z};
+ return {d * vector_x(v), d * vector_y(v), d * vector_z(v)};
}
template <eicd::Vector2D V> V operator*(const V& v, const double d) {
- return {d * v.x, d * v.y};
+ return {d * vector_x(v), d * vector_y(v)};
}
template <eicd::Vector3D V> V operator*(const V& v, const double d) {
- return {d * v.x, d * v.y, d * v.z};
+ return {d * vector_x(v), d * vector_y(v), d * vector_z(v)};
}
template <eicd::VectorND V> V operator-(const V& v1, const V& v2) {
return v1 + (-1. * v2);