diff --git a/JugDigi/src/components/CalorimeterHitDigi.cpp b/JugDigi/src/components/CalorimeterHitDigi.cpp
index 96b6a8fc3e004811a2fdde9a9fc33050280ea763..15fdd999b07a6fea529cab3c25f83512c042d9c4 100644
--- a/JugDigi/src/components/CalorimeterHitDigi.cpp
+++ b/JugDigi/src/components/CalorimeterHitDigi.cpp
@@ -212,7 +212,6 @@ namespace Jug::Digi {
}
// signal sum
- int nhits = 0;
for (auto &[id, hits] : merge_map) {
double edep = hits[0].getEnergy();
double time = hits[0].getContributions(0).getTime();
diff --git a/JugFast/src/components/ParticlesWithTruthPID.cpp b/JugFast/src/components/ParticlesWithTruthPID.cpp
index c78cf0f79e8cf65aee27e0978fb3c4588c6f899e..3886509da8f39a6dd44ce28b3843c02a5e6a9b5a 100644
--- a/JugFast/src/components/ParticlesWithTruthPID.cpp
+++ b/JugFast/src/components/ParticlesWithTruthPID.cpp
@@ -21,12 +21,11 @@ namespace Jug::Fast {
class ParticlesWithTruthPID : public GaudiAlgorithm {
public:
- DataHandle<edm4hep::MCParticleCollection> m_inputTruthCollection{"inputMCParticles", Gaudi::DataHandle::Reader,
- this};
+ DataHandle<edm4hep::MCParticleCollection> m_inputTruthCollection{"inputMCParticles", Gaudi::DataHandle::Reader, this};
DataHandle<eic::TrackParametersCollection> m_inputTrackCollection{"inputTrackParameters", Gaudi::DataHandle::Reader,
this};
- DataHandle<eic::ReconstructedParticleCollection> m_outputParticleCollection{
- "ReconstructedParticles", Gaudi::DataHandle::Writer, this};
+ DataHandle<eic::ReconstructedParticleCollection> m_outputParticleCollection{"ReconstructedParticles",
+ Gaudi::DataHandle::Writer, this};
// Matching momentum tolerance requires 10% by default;
Gaudi::Property<double> m_pRelativeTolerance{this, "pRelativeTolerance", {0.1}};
@@ -35,8 +34,7 @@ public:
// Matchin eta tolerance of 0.1
Gaudi::Property<double> m_etaTolerance{this, "etaTolerance", {0.2}};
- ParticlesWithTruthPID(const std::string& name, ISvcLocator* svcLoc)
- : GaudiAlgorithm(name, svcLoc) {
+ ParticlesWithTruthPID(const std::string& name, ISvcLocator* svcLoc) : GaudiAlgorithm(name, svcLoc) {
declareProperty("inputMCParticles", m_inputTruthCollection, "MCParticles");
declareProperty("inputTrackParameters", m_inputTrackCollection, "outputTrackParameters");
declareProperty("outputParticles", m_outputParticleCollection, "ReconstructedParticles");
@@ -56,28 +54,29 @@ public:
const double sinPhiOver2Tolerance = sin(0.5 * m_phiTolerance);
std::vector<bool> consumed(mc.size(), false);
for (const auto& trk : tracks) {
- const eic::VectorXYZ mom = eicd::sphericalToVector(1.0 / std::abs(trk.qOverP()), trk.theta(), trk.phi());
- const auto charge_rec = std::copysign(1., trk.qOverP());
+ const auto mom = eicd::sphericalToVector(1.0 / std::abs(trk.qOverP()), trk.theta(), trk.phi());
+ const auto charge_rec = trk.charge();
// utility variables for matching
int best_match = -1;
double best_delta = std::numeric_limits<double>::max();
for (size_t ip = 0; ip < mc.size(); ++ip) {
const auto& mcpart = mc[ip];
- if (consumed[ip] || mcpart.getGeneratorStatus() > 1 || mcpart.getCharge() == 0 || mcpart.getCharge() * charge_rec < 0) {
+ if (consumed[ip] || mcpart.getGeneratorStatus() > 1 || mcpart.getCharge() == 0 ||
+ mcpart.getCharge() * charge_rec < 0) {
if (msgLevel(MSG::DEBUG)) {
debug() << "ignoring non-primary/neutral/opposite charge particle" << endmsg;
}
continue;
}
- const auto& p = mcpart.getMomentum();
- const auto p_mag = std::hypot(p.x, p.y, p.z);
- const auto p_phi = std::atan2(p.y, p.x);
- const auto p_eta = std::atanh(p.z / p_mag);
- const double dp_rel = std::abs((mom.mag() - p_mag) / p_mag);
+ const auto& p = mcpart.getMomentum();
+ const auto p_mag = std::hypot(p.x, p.y, p.z);
+ const auto p_phi = std::atan2(p.y, p.x);
+ const auto p_eta = std::atanh(p.z / p_mag);
+ const double dp_rel = std::abs((eicd::magnitude(mom) - p_mag) / p_mag);
// check the tolerance for sin(dphi/2) to avoid the hemisphere problem and allow
// for phi rollovers
- const double dsphi = std::abs(sin(0.5 * (mom.phi() - p_phi)));
- const double deta = std::abs((mom.eta() - p_eta));
+ const double dsphi = std::abs(sin(0.5 * (eicd::angleAzimuthal(mom) - p_phi)));
+ const double deta = std::abs((eicd::eta(mom) - p_eta));
if (dp_rel < m_pRelativeTolerance && deta < m_etaTolerance && dsphi < sinPhiOver2Tolerance) {
const double delta =
@@ -88,21 +87,21 @@ public:
}
}
}
+ auto rec_part = part.create();
int32_t best_pid = 0;
- eic::VectorXYZ vertex;
- float time = 0;
- float mass = 0;
+ auto vertex = rec_part.v();
+ float time = 0;
+ float mass = 0;
eic::Index mcID;
if (best_match >= 0) {
consumed[best_match] = true;
const auto& mcpart = mc[best_match];
best_pid = mcpart.getPDG();
- vertex = {mcpart.getVertex().x, mcpart.getVertex().y, mcpart.getVertex().z};
+ vertex = mcpart.getVertex();
time = mcpart.getTime();
mass = mcpart.getMass();
- //mcID = {mcpart.ID()};
+ // mcID = {mcpart.ID()};
}
- auto rec_part = part.create();
rec_part.p(mom);
rec_part.v(vertex);
rec_part.time(time);
@@ -110,10 +109,7 @@ public:
rec_part.status(static_cast<int16_t>(best_match >= 0 ? 0 : -1));
rec_part.charge(static_cast<int16_t>(charge_rec));
rec_part.weight(1.);
- rec_part.theta(mom.theta());
- rec_part.phi(mom.phi());
- rec_part.momentum(mom.mag());
- rec_part.energy(std::hypot(mom.mag(), mass));
+ rec_part.energy(std::hypot(eicd::magnitude(mom), mass));
rec_part.mass(mass);
rec_part.mcID(mcID);
@@ -121,20 +117,20 @@ public:
if (best_match > 0) {
const auto& mcpart = mc[best_match];
debug() << fmt::format("Matched track with MC particle {}\n", best_match) << endmsg;
- debug() << fmt::format(" - Track: (mom: {}, theta: {}, phi: {}, charge: {})", mom.mag(), mom.theta(),
- mom.phi(), charge_rec)
+ debug() << fmt::format(" - Track: (mom: {}, theta: {}, phi: {}, charge: {})", eicd::magnitude(mom),
+ eicd::anglePolar(mom), eicd::angleAzimuthal(mom), charge_rec)
<< endmsg;
- const auto& p = mcpart.getMomentum();
- const auto p_mag = std::hypot(p.x, p.y, p.z);
- const auto p_phi = std::atan2(p.y, p.x);
- const auto p_theta = std::atan2(std::hypot(p.x, p.y), p.z);
- debug() << fmt::format(" - MC particle: (mom: {}, theta: {}, phi: {}, charge: {}, type: {}",
- p_mag, p_theta, p_phi, mcpart.getCharge(), mcpart.getPDG())
+ const auto& p = mcpart.getMomentum();
+ const auto p_mag = eicd::magnitude(p);
+ const auto p_phi = eicd::angleAzimuthal(p);
+ const auto p_theta = eicd::anglePolar(p);
+ debug() << fmt::format(" - MC particle: (mom: {}, theta: {}, phi: {}, charge: {}, type: {}", p_mag, p_theta,
+ p_phi, mcpart.getCharge(), mcpart.getPDG())
<< endmsg;
} else {
debug() << fmt::format("Did not find a good match for track \n") << endmsg;
- debug() << fmt::format(" - Track: (mom: {}, theta: {}, phi: {}, charge: {})", mom.mag(), mom.theta(),
- mom.phi(), charge_rec)
+ debug() << fmt::format(" - Track: (mom: {}, theta: {}, phi: {}, charge: {})", eicd::magnitude(mom),
+ eicd::anglePolar(mom), eicd::angleAzimuthal(mom), charge_rec)
<< endmsg;
}
}
diff --git a/JugFast/src/components/SmearedFarForwardParticles.cpp b/JugFast/src/components/SmearedFarForwardParticles.cpp
index 542c1ccb6de73b101111b861cc2c07de9db81ef5..6e961030294729e01e322bd952868834731c65a3 100644
--- a/JugFast/src/components/SmearedFarForwardParticles.cpp
+++ b/JugFast/src/components/SmearedFarForwardParticles.cpp
@@ -13,7 +13,7 @@
// Event Model related classes
#include "edm4hep/MCParticleCollection.h"
#include "eicd/ReconstructedParticleCollection.h"
-#include "eicd/VectorXYZ.h"
+#include "eicd/vector_utils.h"
namespace {
enum DetectorTags {
@@ -219,7 +219,7 @@ private:
debug()
<< fmt::format(
"Found ZDC particle: {}, Etrue: {}, Emeas: {}, ptrue: {}, pmeas: {}, theta_true: {}, theta_meas: {}",
- part.getPDG(), E, rec_part.energy(), part_p_mag, rec_part.p().mag(), th, rec_part.p().theta())
+ part.getPDG(), E, rec_part.energy(), part_p_mag, eicd::magnitude(rec_part.p()), th, eicd::anglePolar(rec_part.p()))
<< endmsg;
}
}
@@ -264,7 +264,7 @@ private:
debug() << fmt::format("Found B0 particle: {}, ptrue: {}, pmeas: {}, pttrue: {}, ptmeas: {}, theta_true: {}, "
"theta_meas: {}",
part.getPDG(), part_p_mag, rec_part.momentum(), part_p_pt,
- std::hypot(rec_part.p().x, rec_part.p().y), part_p_theta, rec_part.p().theta())
+ std::hypot(rec_part.p().x, rec_part.p().y), part_p_theta, eicd::anglePolar(rec_part.p()))
<< endmsg;
}
}
@@ -303,7 +303,7 @@ private:
debug() << fmt::format("Found RP particle: {}, ptrue: {}, pmeas: {}, pttrue: {}, ptmeas: {}, theta_true: {}, "
"theta_meas: {}",
part.getPDG(), part_p_mag, rec_part.momentum(), part_p_pt,
- std::hypot(rec_part.p().x, rec_part.p().y), part_p_theta, rec_part.p().theta())
+ std::hypot(rec_part.p().x, rec_part.p().y), part_p_theta, eicd::anglePolar(rec_part.p()))
<< endmsg;
}
}
@@ -348,7 +348,7 @@ private:
debug() << fmt::format("Found OMD particle: {}, ptrue: {}, pmeas: {}, pttrue: {}, ptmeas: {}, theta_true: {}, "
"theta_meas: {}",
part.getPDG(), part_p_mag, rec_part.momentum(), part_p_pt,
- std::hypot(rec_part.p().x, rec_part.p().y), part_p_theta, rec_part.p().theta())
+ std::hypot(rec_part.p().x, rec_part.p().y), part_p_theta, eicd::anglePolar(rec_part.p()))
<< endmsg;
}
}
diff --git a/JugReco/src/components/CalorimeterHitsEtaPhiProjector.cpp b/JugReco/src/components/CalorimeterHitsEtaPhiProjector.cpp
index 6fdf007b381ab978f79e1e391543e267e130681c..746d0fb0edf278c65e07b68abd6cdc3d4c553edb 100644
--- a/JugReco/src/components/CalorimeterHitsEtaPhiProjector.cpp
+++ b/JugReco/src/components/CalorimeterHitsEtaPhiProjector.cpp
@@ -98,8 +98,8 @@ namespace Jug::Reco {
std::unordered_map<std::pair<int64_t, int64_t>, std::vector<eic::ConstCalorimeterHit>, pair_hash> merged_hits;
for (const auto h : *m_inputHitCollection.get()) {
- auto bins = std::make_pair(static_cast<int64_t>(pos2bin(h.position().eta(), gridSizes[0], 0.)),
- static_cast<int64_t>(pos2bin(h.position().phi(), gridSizes[1], 0.)));
+ auto bins = std::make_pair(static_cast<int64_t>(pos2bin(eicd::eta(h.position()), gridSizes[0], 0.)),
+ static_cast<int64_t>(pos2bin(eicd::angleAzimuthal(h.position()), gridSizes[1], 0.)));
merged_hits[bins].push_back(h);
}
@@ -109,11 +109,11 @@ namespace Jug::Reco {
hit.cellID(ref.cellID());
// TODO, we can do timing cut to reject noises
hit.time(ref.time());
- double r = ref.position().mag();
+ double r = eicd::magnitude(ref.position());
double eta = bin2pos(bins.first, gridSizes[0], 0.);
double phi = bin2pos(bins.second, gridSizes[1], 1.);
hit.position(eicd::sphericalToVector(r, eicd::etaToAngle(eta), phi));
- hit.dimension({gridSizes[0], gridSizes[1], 0.});
+ hit.dimension({static_cast<float>(gridSizes[0]), static_cast<float>(gridSizes[1]), 0.});
// merge energy
hit.energy(0.);
for (const auto &h : hits) {
diff --git a/JugReco/src/components/CalorimeterIslandCluster.cpp b/JugReco/src/components/CalorimeterIslandCluster.cpp
index 5da860f7b0d23f5b5b1467549f767302affab934..220bc6a2eaf982e69c3da70b0f4be6de0c58441f 100644
--- a/JugReco/src/components/CalorimeterIslandCluster.cpp
+++ b/JugReco/src/components/CalorimeterIslandCluster.cpp
@@ -33,8 +33,8 @@
#include "eicd/CalorimeterHitCollection.h"
#include "eicd/ClusterCollection.h"
#include "eicd/ProtoClusterCollection.h"
-#include "eicd/VectorXY.h"
-#include "eicd/VectorXYZ.h"
+#include "eicd/Vector2f.h"
+#include "eicd/Vector3f.h"
#include "eicd/vector_utils.h"
using namespace Gaudi::Units;
@@ -45,35 +45,35 @@ using CaloHit = eic::ConstCalorimeterHit;
using CaloHitCollection = eic::CalorimeterHitCollection;
// helper functions to get distance between hits
-static eic::VectorXY localDistXY(const CaloHit& h1, const CaloHit& h2) {
+static eicd::Vector2f localDistXY(const CaloHit& h1, const CaloHit& h2) {
const auto delta = h1.local() - h2.local();
return {delta.x, delta.y};
}
-static eic::VectorXY localDistXZ(const CaloHit& h1, const CaloHit& h2) {
+static eicd::Vector2f localDistXZ(const CaloHit& h1, const CaloHit& h2) {
const auto delta = h1.local() - h2.local();
return {delta.x, delta.z};
}
-static eic::VectorXY localDistYZ(const CaloHit& h1, const CaloHit& h2) {
+static eicd::Vector2f localDistYZ(const CaloHit& h1, const CaloHit& h2) {
const auto delta = h1.local() - h2.local();
return {delta.y, delta.z};
}
-static eic::VectorXY dimScaledLocalDistXY(const CaloHit& h1, const CaloHit& h2) {
+static eicd::Vector2f dimScaledLocalDistXY(const CaloHit& h1, const CaloHit& h2) {
const auto delta = h1.local() - h2.local();
const auto dimsum = h1.dimension() + h2.dimension();
return {2 * delta.x / dimsum.x, 2 * delta.y / dimsum.y};
}
-static eic::VectorXY globalDistRPhi(const CaloHit& h1, const CaloHit& h2) {
+static eicd::Vector2f globalDistRPhi(const CaloHit& h1, const CaloHit& h2) {
return {eicd::magnitude(h1.position()) - eicd::magnitude(h2.position()),
eicd::angleAzimuthal(h1.position()) - eicd::angleAzimuthal(h2.position())};
}
-static eic::VectorXY globalDistEtaPhi(const CaloHit& h1,
+static eicd::Vector2f globalDistEtaPhi(const CaloHit& h1,
const CaloHit& h2) {
return {eicd::eta(h1.position()) - eicd::eta(h2.position()),
eicd::angleAzimuthal(h1.position()) - eicd::angleAzimuthal(h2.position())};
}
// name: {method, units}
static std::map<std::string,
- std::tuple<std::function<eic::VectorXY(const CaloHit&, const CaloHit&)>, std::vector<double>>>
+ std::tuple<std::function<eicd::Vector2f(const CaloHit&, const CaloHit&)>, std::vector<double>>>
distMethods{
{"localDistXY", {localDistXY, {mm, mm}}}, {"localDistXZ", {localDistXZ, {mm, mm}}},
{"localDistYZ", {localDistYZ, {mm, mm}}}, {"dimScaledLocalDistXY", {dimScaledLocalDistXY, {1., 1.}}},
@@ -113,7 +113,7 @@ public:
Gaudi::Property<std::vector<double>> u_globalDistEtaPhi{this, "globalDistEtaPhi", {}};
Gaudi::Property<std::vector<double>> u_dimScaledLocalDistXY{this, "dimScaledLocalDistXY", {1.8, 1.8}};
// neighbor checking function
- std::function<eic::VectorXY(const CaloHit&, const CaloHit&)> hitsDist;
+ std::function<eicd::Vector2f(const CaloHit&, const CaloHit&)> hitsDist;
// unitless counterparts of the input parameters
double minClusterHitEdep, minClusterCenterEdep, sectorDist;
@@ -231,7 +231,7 @@ private:
// in the same sector
if (h1.sector() == h2.sector()) {
auto dist = hitsDist(h1, h2);
- return (dist.x <= neighbourDist[0]) && (dist.y <= neighbourDist[1]);
+ return (dist.a <= neighbourDist[0]) && (dist.b <= neighbourDist[1]);
// different sector, local coordinates do not work, using global coordinates
} else {
// sector may have rotation (barrel), so z is included
@@ -354,7 +354,7 @@ private:
for (const auto& chit : maxima) {
double dist_ref = chit.dimension().x;
double energy = chit.energy();
- double dist = hitsDist(chit, hit).mag();
+ double dist = eicd::magnitude(hitsDist(chit, hit));
weights[j] = std::exp(-dist / dist_ref) * energy;
j += 1;
}
diff --git a/JugReco/src/components/ClusterRecoCoG.cpp b/JugReco/src/components/ClusterRecoCoG.cpp
index ec93e131095203a3d09761940f9cf95e81665129..f9652f3eee49e9bae063cde7480ae8969307c0d2 100644
--- a/JugReco/src/components/ClusterRecoCoG.cpp
+++ b/JugReco/src/components/ClusterRecoCoG.cpp
@@ -186,18 +186,18 @@ private:
// center of gravity with logarithmic weighting
float tw = 0.;
- eic::VectorXYZ v;
+ auto v = cl.position();
for (unsigned i = 0; i < pcl.hits().size(); ++i) {
const auto& hit = pcl.hits()[i];
const auto weight = pcl.weights()[i];
float w = weightFunc(hit.energy() * weight, totalE, m_logWeightBase.value(), 0);
tw += w;
- v = v.add(hit.position().scale(w));
+ v = v + (hit.position() * w);
}
if (tw == 0.) {
warning() << "zero total weights encountered, you may want to adjust your weighting parameter." << endmsg;
}
- cl.position(v.scale(1 / tw));
+ cl.position(v / tw);
cl.positionError({}); // @TODO: Covariance matrix
// Optionally constrain the cluster to the hit eta values
diff --git a/JugReco/src/components/EnergyPositionClusterMerger.cpp b/JugReco/src/components/EnergyPositionClusterMerger.cpp
index 7b56962071824b58cb5a511958cd9b4b82fe4360..76063ceaa177944e92d57231b20c37b7f8db6a64 100644
--- a/JugReco/src/components/EnergyPositionClusterMerger.cpp
+++ b/JugReco/src/components/EnergyPositionClusterMerger.cpp
@@ -13,6 +13,7 @@
// Event Model related classes
#include "eicd/ClusterCollection.h"
+#include "eicd/vector_utils.h"
using namespace Gaudi::Units;
@@ -46,8 +47,7 @@ public:
double m_phiTolerance;
public:
- EnergyPositionClusterMerger(const std::string& name, ISvcLocator* svcLoc)
- : GaudiAlgorithm(name, svcLoc) {
+ EnergyPositionClusterMerger(const std::string& name, ISvcLocator* svcLoc) : GaudiAlgorithm(name, svcLoc) {
declareProperty("energyClusters", m_energyClusters, "Cluster collection with good energy precision");
declareProperty("positionClusters", m_positionClusters, "Cluster collection with good position precision");
declareProperty("outputClusters", m_outputClusters, "");
@@ -64,7 +64,7 @@ public:
const auto& e_clus = *(m_energyClusters.get());
const auto& pos_clus = *(m_positionClusters.get());
// output
- auto& merged = *(m_outputClusters.createAndPut());
+ auto& merged = *(m_outputClusters.createAndPut());
std::vector<bool> consumed(e_clus.size(), false);
@@ -80,7 +80,7 @@ public:
const auto& ec = e_clus[ie];
// 1. stop if not within tolerance
// (make sure to handle rollover of phi properly)
- double dphi = pc.position().phi() - ec.position().phi();
+ double dphi = eicd::angleAzimuthal(pc.position()) - eicd::angleAzimuthal(ec.position());
if (std::abs(dphi) > M_PI) {
dphi = std::abs(dphi) - M_PI;
}
@@ -116,14 +116,14 @@ public:
consumed[best_match] = true;
if (msgLevel(MSG::DEBUG)) {
debug() << fmt::format("Matched position cluster {} with energy cluster {}\n", pc.id(), ec.id()) << endmsg;
- debug() << fmt::format(" - Position cluster: (E: {}, phi: {}, z: {})", pc.energy(), pc.position().phi(),
- pc.position().z)
+ debug() << fmt::format(" - Position cluster: (E: {}, phi: {}, z: {})", pc.energy(),
+ eicd::angleAzimuthal(pc.position()), pc.position().z)
<< endmsg;
- debug() << fmt::format(" - Energy cluster: (E: {}, phi: {}, z: {})", ec.energy(), ec.position().phi(),
- ec.position().z)
+ debug() << fmt::format(" - Energy cluster: (E: {}, phi: {}, z: {})", ec.energy(),
+ eicd::angleAzimuthal(ec.position()), ec.position().z)
<< endmsg;
debug() << fmt::format(" ---> Merged cluster: (E: {}, phi: {}, z: {})", new_clus.energy(),
- new_clus.position().phi(), new_clus.position().z)
+ eicd::angleAzimuthal(new_clus.position()), new_clus.position().z)
<< endmsg;
}
}
diff --git a/JugReco/src/components/FarForwardParticles.cpp b/JugReco/src/components/FarForwardParticles.cpp
index eb21a3985e083a44a5839732dc285abd0657ce54..8ebbd83cdcc4250513f40ea60c242b26ad7b2a6b 100644
--- a/JugReco/src/components/FarForwardParticles.cpp
+++ b/JugReco/src/components/FarForwardParticles.cpp
@@ -18,6 +18,7 @@
// Event Model related classes
#include "eicd/ReconstructedParticleCollection.h"
#include "eicd/TrackerHitCollection.h"
+#include <eicd/vector_utils.h>
namespace Jug::Reco {
@@ -39,16 +40,14 @@ class FarForwardParticles : public GaudiAlgorithm {
Gaudi::Property<std::string> m_readout{this, "readoutClass", ""};
Gaudi::Property<std::string> m_layerField{this, "layerField", ""};
Gaudi::Property<std::string> m_sectorField{this, "sectorField", ""};
- SmartIF<IGeoSvc> m_geoSvc;
- dd4hep::BitFieldCoder* id_dec = nullptr;
- size_t sector_idx, layer_idx;
+ SmartIF<IGeoSvc> m_geoSvc;
+ dd4hep::BitFieldCoder* id_dec = nullptr;
+ size_t sector_idx, layer_idx;
- Gaudi::Property<std::string> m_localDetElement{this, "localDetElement", ""};
+ Gaudi::Property<std::string> m_localDetElement{this, "localDetElement", ""};
Gaudi::Property<std::vector<std::string>> u_localDetFields{this, "localDetFields", {}};
- dd4hep::DetElement local;
- size_t local_mask = ~0;
-
-
+ dd4hep::DetElement local;
+ size_t local_mask = ~0;
const double aXRP[2][2] = {{2.102403743, 29.11067626}, {0.186640381, 0.192604619}};
const double aYRP[2][2] = {{0.0000159900, 3.94082098}, {0.0000079946, -0.1402995}};
@@ -57,28 +56,26 @@ class FarForwardParticles : public GaudiAlgorithm {
double aYRPinv[2][2] = {0.0, 0.0};
public:
- FarForwardParticles(const std::string& name, ISvcLocator* svcLoc)
- : GaudiAlgorithm(name, svcLoc) {
+ FarForwardParticles(const std::string& name, ISvcLocator* svcLoc) : GaudiAlgorithm(name, svcLoc) {
declareProperty("inputCollection", m_inputHitCollection, "FarForwardTrackerHits");
declareProperty("outputCollection", m_outputParticles, "ReconstructedParticles");
}
// See Wouter's example to extract local coordinates CalorimeterHitReco.cpp
- //includes DDRec/CellIDPositionConverter.here
- //See tutorial
+ // includes DDRec/CellIDPositionConverter.here
+ // See tutorial
// auto converter = m_GeoSvc ....
- //https://eicweb.phy.anl.gov/EIC/juggler/-/blob/master/JugReco/src/components/CalorimeterHitReco.cpp - line 200
- //include the Eigen libraries, used in ACTS, for the linear algebra.
+ // https://eicweb.phy.anl.gov/EIC/juggler/-/blob/master/JugReco/src/components/CalorimeterHitReco.cpp - line 200
+ // include the Eigen libraries, used in ACTS, for the linear algebra.
StatusCode initialize() override {
- if (GaudiAlgorithm::initialize().isFailure()){
+ if (GaudiAlgorithm::initialize().isFailure()) {
return StatusCode::FAILURE;
}
m_geoSvc = service(m_geoSvcName);
if (!m_geoSvc) {
error() << "Unable to locate Geometry Service. "
- << "Make sure you have GeoSvc and SimSvc in the right order in the configuration."
- << endmsg;
+ << "Make sure you have GeoSvc and SimSvc in the right order in the configuration." << endmsg;
return StatusCode::FAILURE;
}
@@ -107,14 +104,12 @@ public:
if (m_localDetElement.value().size()) {
try {
local = m_geoSvc->detector()->detector(m_localDetElement.value());
- info() << "Local coordinate system from DetElement " << m_localDetElement.value()
- << endmsg;
+ info() << "Local coordinate system from DetElement " << m_localDetElement.value() << endmsg;
} catch (...) {
- error() << "Failed to locate local coordinate system from DetElement "
- << m_localDetElement.value() << endmsg;
+ error() << "Failed to locate local coordinate system from DetElement " << m_localDetElement.value() << endmsg;
return StatusCode::FAILURE;
}
- // or get from fields
+ // or get from fields
} else {
std::vector<std::pair<std::string, int>> fields;
for (auto& f : u_localDetFields.value()) {
@@ -125,11 +120,11 @@ public:
if (fields.empty()) {
local_mask = ~0;
}
- //info() << fmt::format("Local DetElement mask {:#064b} from fields [{}]", local_mask,
+ // info() << fmt::format("Local DetElement mask {:#064b} from fields [{}]", local_mask,
// fmt::join(fields, ", "))
// << endmsg;
}
-
+
double det = aXRP[0][0] * aXRP[1][1] - aXRP[0][1] * aXRP[1][0];
if (det == 0) {
@@ -175,32 +170,30 @@ public:
for (const auto& h : *rawhits) {
- auto cellID = h.cellID();
+ auto cellID = h.cellID();
// The actual hit position in Global Coordinates
- //auto pos0 = h.position();
-
+ // auto pos0 = h.position();
+
auto gpos = converter->position(cellID);
- // local positions
- if (m_localDetElement.value().empty()) {
- auto volman = m_geoSvc->detector()->volumeManager();
- local = volman.lookupDetElement(cellID);
- }
- auto pos0 = local.nominal().worldToLocal(dd4hep::Position(gpos.x(), gpos.y(), gpos.z())); //hit position in local coordinates
+ // local positions
+ if (m_localDetElement.value().empty()) {
+ auto volman = m_geoSvc->detector()->volumeManager();
+ local = volman.lookupDetElement(cellID);
+ }
+ auto pos0 = local.nominal().worldToLocal(
+ dd4hep::Position(gpos.x(), gpos.y(), gpos.z())); // hit position in local coordinates
-
// auto mom0 = h.momentum;
// auto pidCode = h.g4ID;
auto eDep = h.edep();
-
-
if (eDep < 0.00001) {
continue;
}
if (eventReset < 2) {
- hitx.push_back(pos0.x());// - local_x_offset_station_2);
- } // use station 2 for both offsets since it is used for the reference orbit
+ hitx.push_back(pos0.x()); // - local_x_offset_station_2);
+ } // use station 2 for both offsets since it is used for the reference orbit
else {
hitx.push_back(pos0.x()); // - local_x_offset_station_2);
}
@@ -226,7 +219,7 @@ public:
if (base == 0) {
warning() << "Detector separation = 0!"
- << "Cannot calculate slope!" << endmsg;
+ << "Cannot calculate slope!" << endmsg;
return StatusCode::SUCCESS;
}
@@ -254,49 +247,23 @@ public:
double p = nomMomentum * (1 + 0.01 * Xip[0]);
double norm = std::sqrt(1.0 + rsx * rsx + rsy * rsy);
- double prec[3];
- prec[0] = p * rsx / norm;
- prec[1] = p * rsy / norm;
- prec[2] = p / norm;
-
- // double pT_reco = std::sqrt(prec[0]*prec[0] + prec[1]*prec[1]);
- float p_reco = std::sqrt(prec[0] * prec[0] + prec[1] * prec[1] + prec[2] * prec[2]);
+ const float prec[3] = {static_cast<float>(p * rsx / norm), static_cast<float>(p * rsy / norm),
+ static_cast<float>(p / norm)};
//----- end RP reconstruction code ------
- eic::VectorXYZ recoMom(prec[0], prec[1], prec[2]);
- eic::VectorXYZ primVtx(0, 0, 0);
- eic::Weight wgt(1);
-
- // ReconstructedParticle (eic::Index ID, eic::VectorXYZ p, eic::VectorXYZ v, float time, std::int32_t pid,
- // std::int16_t status, std::int16_t charge, eic::Weight weight, eic::Direction direction, float momentum,
- // float energy, float mass)
-
- // eic::ReconstructedParticle rec_part{{part.ID(), algorithmID()},
- // mom3s,
- //{part.vs().x, part.vs().y, part.vs().z},
- // static_cast<float>(part.time()),
- // part.pdgID(),
- // 0,
- // static_cast<int16_t>(part.charge()),
- // 1.,
- //{mom3s.theta(), mom3s.phi()},
- // static_cast<float>(moms),
- // static_cast<float>(Es),
- // static_cast<float>(part.mass())};
-
eic::ReconstructedParticle rpTrack;
- rpTrack.p(recoMom);
- rpTrack.v(primVtx);
+ rpTrack.p({prec});
+ rpTrack.v({0, 0, 0});
rpTrack.time(0);
rpTrack.pid(2122);
rpTrack.status(0);
rpTrack.charge(1);
rpTrack.weight(1.);
- rpTrack.theta(recoMom.theta());
- rpTrack.phi(recoMom.phi());
- rpTrack.momentum(p_reco);
- rpTrack.energy(std::hypot(p_reco, .938272));
+ rpTrack.theta(eicd::anglePolar(rpTrack.p()));
+ rpTrack.phi(eicd::angleAzimuthal(rpTrack.p()));
+ rpTrack.momentum(eicd::magnitude(rpTrack.p()));
+ rpTrack.energy(std::hypot(rpTrack.momentum(), .938272));
rpTrack.mass(.938272);
rc->push_back(rpTrack);
diff --git a/JugReco/src/components/FarForwardParticlesOMD.cpp b/JugReco/src/components/FarForwardParticlesOMD.cpp
index 6539799cdd53e796fee9f5c145453ccb9ed713a4..5da4f762e87e0831882a12018640b2c0452fc6cd 100644
--- a/JugReco/src/components/FarForwardParticlesOMD.cpp
+++ b/JugReco/src/components/FarForwardParticlesOMD.cpp
@@ -13,6 +13,7 @@
// Event Model related classes
#include "eicd/ReconstructedParticleCollection.h"
#include "eicd/TrackerHitCollection.h"
+#include <eicd/vector_utils.h>
namespace Jug::Reco {
@@ -153,49 +154,23 @@ public:
double p = nomMomentum * (1 + 0.01 * Xip[0]);
double norm = std::sqrt(1.0 + rsx * rsx + rsy * rsy);
- double prec[3];
- prec[0] = p * rsx / norm;
- prec[1] = p * rsy / norm;
- prec[2] = p / norm;
-
- // double pT_reco = std::sqrt(prec[0]*prec[0] + prec[1]*prec[1]);
- float p_reco = std::sqrt(prec[0] * prec[0] + prec[1] * prec[1] + prec[2] * prec[2]);
+ const float prec[3] = {static_cast<float>(p * rsx / norm), static_cast<float>(p * rsy / norm),
+ static_cast<float>(p / norm)};
//----- end RP reconstruction code ------
- eic::VectorXYZ recoMom(prec[0], prec[1], prec[2]);
- eic::VectorXYZ primVtx(0, 0, 0);
- eic::Weight wgt(1);
-
- // ReconstructedParticle (eic::Index ID, eic::VectorXYZ p, eic::VectorXYZ v, float time, std::int32_t pid,
- // std::int16_t status, std::int16_t charge, eic::Weight weight, eic::Direction direction, float momentum, float
- // energy, float mass)
-
- // eic::ReconstructedParticle rec_part{{part.ID(), algorithmID()},
- // mom3s,
- //{part.vs().x, part.vs().y, part.vs().z},
- // static_cast<float>(part.time()),
- // part.pdgID(),
- // 0,
- // static_cast<int16_t>(part.charge()),
- // 1.,
- //{mom3s.theta(), mom3s.phi()},
- // static_cast<float>(moms),
- // static_cast<float>(Es),
- // static_cast<float>(part.mass())};
-
eic::ReconstructedParticle rpTrack;
- rpTrack.p(recoMom);
- rpTrack.v(primVtx);
+ rpTrack.p({prec});
+ rpTrack.v({0, 0, 0});
rpTrack.time(0);
rpTrack.pid(2122);
rpTrack.status(0);
rpTrack.charge(1);
rpTrack.weight(1.);
- rpTrack.theta(recoMom.theta());
- rpTrack.phi(recoMom.phi());
- rpTrack.momentum(p_reco);
- rpTrack.energy(std::hypot(p_reco, .938272));
+ rpTrack.theta(eicd::anglePolar(rpTrack.p()));
+ rpTrack.phi(eicd::angleAzimuthal(rpTrack.p()));
+ rpTrack.momentum(eicd::magnitude(rpTrack.p()));
+ rpTrack.energy(std::hypot(rpTrack.momentum(), .938272));
rpTrack.mass(.938272);
rc->push_back(rpTrack);
diff --git a/JugReco/src/components/ImagingClusterReco.cpp b/JugReco/src/components/ImagingClusterReco.cpp
index 7075eddc6507c78dab07d4d5bd59b21a704e6212..baa683d7290dc97637ed514615cacb15aa38a5b8 100644
--- a/JugReco/src/components/ImagingClusterReco.cpp
+++ b/JugReco/src/components/ImagingClusterReco.cpp
@@ -134,7 +134,7 @@ private:
const auto& hits = pcl.hits();
const auto& weights = pcl.weights();
// using map to have hits sorted by layer
- std::map<int, std::vector<std::pair<eic::ConstCalorimeterHit, eic::Weight>>> layer_map;
+ std::map<int, std::vector<std::pair<eic::ConstCalorimeterHit, float>>> layer_map;
for (unsigned i = 0; i < hits.size(); ++i) {
const auto hit = hits[i];
auto lid = hit.layer();
@@ -153,7 +153,7 @@ private:
return cl_layers;
}
- eic::Cluster reconstruct_layer(const std::vector<std::pair<eic::ConstCalorimeterHit, eic::Weight>>& hits) const {
+ eic::Cluster reconstruct_layer(const std::vector<std::pair<eic::ConstCalorimeterHit, float>>& hits) const {
eic::Cluster layer;
layer.type(ClusterType::kClusterSlice);
// Calculate averages
@@ -162,7 +162,7 @@ private:
double time;
double timeError;
double sumOfWeights = 0;
- eic::VectorXYZ pos;
+ auto pos = layer.position();
for (const auto& [hit, weight] : hits) {
energy += hit.energy() * weight;
energyError += std::pow(hit.energyError() * weight, 2);
@@ -199,13 +199,13 @@ private:
const auto& weights = pcl.weights();
cluster.type(ClusterType::kCluster3D);
- double energy = 0.;
- double energyError = 0.;
- double time = 0.;
- double timeError = 0.;
- double meta = 0.;
- double mphi = 0.;
- double r = 9999 * cm;
+ double energy = 0.;
+ double energyError = 0.;
+ double time = 0.;
+ double timeError = 0.;
+ double meta = 0.;
+ double mphi = 0.;
+ double r = 9999 * cm;
for (unsigned i = 0; i < hits.size(); ++i) {
const auto& hit = hits[i];
const auto weight = weights[i];
@@ -230,8 +230,8 @@ private:
// shower radius estimate (eta-phi plane)
double radius = 0.;
for (const auto& hit : hits) {
- radius +=
- pow2(hit.position().eta() - cluster.position().eta()) + pow2(hit.position().phi() - cluster.position().phi());
+ radius += pow2(eicd::eta(hit.position()) - eicd::eta(cluster.position())) +
+ pow2(eicd::angleAzimuthal(hit.position()) - eicd::angleAzimuthal(cluster.position()));
}
cluster.addshapeParameters(std::sqrt(radius / cluster.nhits()));
// Skewedness not calculated TODO
@@ -248,7 +248,7 @@ private:
std::pair<double /* polar */, double /* azimuthal */> fit_track(const std::vector<eic::Cluster>& layers) const {
int nrows = 0;
- eic::VectorXYZ mean_pos{0, 0, 0};
+ eicd::Vector3f mean_pos{0, 0, 0};
for (const auto& layer : layers) {
if ((layer.nhits() > 0) && (layer.hits(0).layer() <= m_trackStopLayer)) {
mean_pos = mean_pos + layer.position();
diff --git a/JugReco/src/components/ImagingTopoCluster.cpp b/JugReco/src/components/ImagingTopoCluster.cpp
index 528df7c73e149bb52da8f9a185a22af589366e8f..0ae7dd3b32affeffc9bc8f8efabc521db0ccb650 100644
--- a/JugReco/src/components/ImagingTopoCluster.cpp
+++ b/JugReco/src/components/ImagingTopoCluster.cpp
@@ -27,208 +27,202 @@
#include "JugReco/ClusterTypes.h"
// Event Model related classes
-#include "eicd/ProtoClusterCollection.h"
#include "eicd/CalorimeterHitCollection.h"
+#include "eicd/ProtoClusterCollection.h"
+#include "eicd/vector_utils.h"
using namespace Gaudi::Units;
namespace Jug::Reco {
- /** Topological Cell Clustering Algorithm.
- *
- * Topological Cell Clustering Algorithm for Imaging Calorimetry
- * 1. group all the adjacent pixels
- *
- * Author: Chao Peng (ANL), 06/02/2021
- * References: https://arxiv.org/pdf/1603.02934.pdf
- *
- * \ingroup reco
- */
- class ImagingTopoCluster : public GaudiAlgorithm {
- public:
- // maximum difference in layer numbers that can be considered as neighbours
- Gaudi::Property<int> m_neighbourLayersRange{this, "neighbourLayersRange", 1};
- // maximum distance of local (x, y) to be considered as neighbors at the same layer
- Gaudi::Property<std::vector<double>> u_localDistXY{this, "localDistXY", {1.0 * mm, 1.0 * mm}};
- // maximum distance of global (eta, phi) to be considered as neighbors at different layers
- Gaudi::Property<std::vector<double>> u_layerDistEtaPhi{this, "layerDistEtaPhi", {0.01, 0.01}};
- // maximum global distance to be considered as neighbors in different sectors
- Gaudi::Property<double> m_sectorDist{this, "sectorDist", 1.0 * cm};
-
- // minimum hit energy to participate clustering
- Gaudi::Property<double> m_minClusterHitEdep{this, "minClusterHitEdep", 0.};
- // minimum cluster center energy (to be considered as a seed for cluster)
- Gaudi::Property<double> m_minClusterCenterEdep{this, "minClusterCenterEdep", 0.};
- // minimum cluster energy (to save this cluster)
- Gaudi::Property<double> m_minClusterEdep{this, "minClusterEdep", 0.5 * MeV};
- // minimum number of hits (to save this cluster)
- Gaudi::Property<int> m_minClusterNhits{this, "minClusterNhits", 10};
- // input hits collection
- DataHandle<eic::CalorimeterHitCollection> m_inputHitCollection{"inputHitCollection",
- Gaudi::DataHandle::Reader, this};
- // output clustered hits
- DataHandle<eic::ProtoClusterCollection> m_outputProtoClusterCollection{"outputProtoClusterCollection",
- Gaudi::DataHandle::Writer, this};
-
- // unitless counterparts of the input parameters
- double localDistXY[2], layerDistEtaPhi[2], sectorDist;
- double minClusterHitEdep, minClusterCenterEdep, minClusterEdep, minClusterNhits;
-
- ImagingTopoCluster(const std::string& name, ISvcLocator* svcLoc)
- : GaudiAlgorithm(name, svcLoc)
- {
- declareProperty("inputHitCollection", m_inputHitCollection, "");
- declareProperty("outputProtoClusterCollection", m_outputProtoClusterCollection, "");
+/** Topological Cell Clustering Algorithm.
+ *
+ * Topological Cell Clustering Algorithm for Imaging Calorimetry
+ * 1. group all the adjacent pixels
+ *
+ * Author: Chao Peng (ANL), 06/02/2021
+ * References: https://arxiv.org/pdf/1603.02934.pdf
+ *
+ * \ingroup reco
+ */
+class ImagingTopoCluster : public GaudiAlgorithm {
+public:
+ // maximum difference in layer numbers that can be considered as neighbours
+ Gaudi::Property<int> m_neighbourLayersRange{this, "neighbourLayersRange", 1};
+ // maximum distance of local (x, y) to be considered as neighbors at the same layer
+ Gaudi::Property<std::vector<double>> u_localDistXY{this, "localDistXY", {1.0 * mm, 1.0 * mm}};
+ // maximum distance of global (eta, phi) to be considered as neighbors at different layers
+ Gaudi::Property<std::vector<double>> u_layerDistEtaPhi{this, "layerDistEtaPhi", {0.01, 0.01}};
+ // maximum global distance to be considered as neighbors in different sectors
+ Gaudi::Property<double> m_sectorDist{this, "sectorDist", 1.0 * cm};
+
+ // minimum hit energy to participate clustering
+ Gaudi::Property<double> m_minClusterHitEdep{this, "minClusterHitEdep", 0.};
+ // minimum cluster center energy (to be considered as a seed for cluster)
+ Gaudi::Property<double> m_minClusterCenterEdep{this, "minClusterCenterEdep", 0.};
+ // minimum cluster energy (to save this cluster)
+ Gaudi::Property<double> m_minClusterEdep{this, "minClusterEdep", 0.5 * MeV};
+ // minimum number of hits (to save this cluster)
+ Gaudi::Property<int> m_minClusterNhits{this, "minClusterNhits", 10};
+ // input hits collection
+ DataHandle<eic::CalorimeterHitCollection> m_inputHitCollection{"inputHitCollection", Gaudi::DataHandle::Reader, this};
+ // output clustered hits
+ DataHandle<eic::ProtoClusterCollection> m_outputProtoClusterCollection{"outputProtoClusterCollection",
+ Gaudi::DataHandle::Writer, this};
+
+ // unitless counterparts of the input parameters
+ double localDistXY[2], layerDistEtaPhi[2], sectorDist;
+ double minClusterHitEdep, minClusterCenterEdep, minClusterEdep, minClusterNhits;
+
+ ImagingTopoCluster(const std::string& name, ISvcLocator* svcLoc) : GaudiAlgorithm(name, svcLoc) {
+ declareProperty("inputHitCollection", m_inputHitCollection, "");
+ declareProperty("outputProtoClusterCollection", m_outputProtoClusterCollection, "");
+ }
+
+ StatusCode initialize() override {
+ if (GaudiAlgorithm::initialize().isFailure()) {
+ return StatusCode::FAILURE;
}
- StatusCode initialize() override
- {
- if (GaudiAlgorithm::initialize().isFailure()) {
- return StatusCode::FAILURE;
- }
+ // unitless conversion
+ // sanity checks
+ if (u_localDistXY.size() != 2) {
+ error() << "Expected 2 values (x_dist, y_dist) for localDistXY" << endmsg;
+ return StatusCode::FAILURE;
+ }
+ if (u_layerDistEtaPhi.size() != 2) {
+ error() << "Expected 2 values (eta_dist, phi_dist) for layerDistEtaPhi" << endmsg;
+ return StatusCode::FAILURE;
+ }
- // unitless conversion
- // sanity checks
- if (u_localDistXY.size() != 2) {
- error() << "Expected 2 values (x_dist, y_dist) for localDistXY" << endmsg;
- return StatusCode::FAILURE;
+ // using juggler internal units (GeV, mm, ns, rad)
+ localDistXY[0] = u_localDistXY.value()[0] / mm;
+ localDistXY[1] = u_localDistXY.value()[1] / mm;
+ layerDistEtaPhi[0] = u_layerDistEtaPhi.value()[0];
+ layerDistEtaPhi[1] = u_layerDistEtaPhi.value()[1] / rad;
+ sectorDist = m_sectorDist.value() / mm;
+ minClusterHitEdep = m_minClusterHitEdep.value() / GeV;
+ minClusterCenterEdep = m_minClusterCenterEdep.value() / GeV;
+ minClusterEdep = m_minClusterEdep.value() / GeV;
+
+ // summarize the clustering parameters
+ info() << fmt::format("Local clustering (same sector and same layer): "
+ "Local [x, y] distance between hits <= [{:.4f} mm, {:.4f} mm].",
+ localDistXY[0], localDistXY[1])
+ << endmsg;
+ info() << fmt::format("Neighbour layers clustering (same sector and layer id within +- {:d}: "
+ "Global [eta, phi] distance between hits <= [{:.4f}, {:.4f} rad].",
+ m_neighbourLayersRange.value(), layerDistEtaPhi[0], layerDistEtaPhi[1])
+ << endmsg;
+ info() << fmt::format("Neighbour sectors clustering (different sector): "
+ "Global distance between hits <= {:.4f} mm.",
+ sectorDist)
+ << endmsg;
+
+ return StatusCode::SUCCESS;
+ }
+
+ StatusCode execute() override {
+ // input collections
+ const auto& hits = *m_inputHitCollection.get();
+ // Create output collections
+ auto& proto = *m_outputProtoClusterCollection.createAndPut();
+
+ // group neighboring hits
+ std::vector<bool> visits(hits.size(), false);
+ std::vector<std::vector<std::pair<uint32_t, eic::ConstCalorimeterHit>>> groups;
+ for (size_t i = 0; i < hits.size(); ++i) {
+ if (msgLevel(MSG::DEBUG)) {
+ debug() << fmt::format("hit {:d}: local position = ({}, {}, {}), global position = ({}, {}, {})", i + 1,
+ hits[i].local().x, hits[i].local().y, hits[i].position().z, hits[i].position().x,
+ hits[i].position().y, hits[i].position().z)
+ << endmsg;
}
- if (u_layerDistEtaPhi.size() != 2) {
- error() << "Expected 2 values (eta_dist, phi_dist) for layerDistEtaPhi" << endmsg;
- return StatusCode::FAILURE;
+ // already in a group, or not energetic enough to form a cluster
+ if (visits[i] || hits[i].energy() < minClusterCenterEdep) {
+ continue;
+ }
+ groups.emplace_back();
+ // create a new group, and group all the neighboring hits
+ dfs_group(groups.back(), i, hits, visits);
+ }
+ if (msgLevel(MSG::DEBUG)) {
+ debug() << "found " << groups.size() << " potential clusters (groups of hits)" << endmsg;
+ for (size_t i = 0; i < groups.size(); ++i) {
+ debug() << fmt::format("group {}: {} hits", i, groups[i].size()) << endmsg;
}
-
- // using juggler internal units (GeV, mm, ns, rad)
- localDistXY[0] = u_localDistXY.value()[0]/mm;
- localDistXY[1] = u_localDistXY.value()[1]/mm;
- layerDistEtaPhi[0] = u_layerDistEtaPhi.value()[0];
- layerDistEtaPhi[1] = u_layerDistEtaPhi.value()[1]/rad;
- sectorDist = m_sectorDist.value()/mm;
- minClusterHitEdep = m_minClusterHitEdep.value()/GeV;
- minClusterCenterEdep = m_minClusterCenterEdep.value()/GeV;
- minClusterEdep = m_minClusterEdep.value()/GeV;
-
- // summarize the clustering parameters
- info() << fmt::format("Local clustering (same sector and same layer): "
- "Local [x, y] distance between hits <= [{:.4f} mm, {:.4f} mm].",
- localDistXY[0], localDistXY[1]) << endmsg;
- info() << fmt::format("Neighbour layers clustering (same sector and layer id within +- {:d}: "
- "Global [eta, phi] distance between hits <= [{:.4f}, {:.4f} rad].",
- m_neighbourLayersRange.value(), layerDistEtaPhi[0], layerDistEtaPhi[1]) << endmsg;
- info() << fmt::format("Neighbour sectors clustering (different sector): "
- "Global distance between hits <= {:.4f} mm.",
- sectorDist) << endmsg;
-
- return StatusCode::SUCCESS;
}
- StatusCode execute() override
- {
- // input collections
- const auto& hits = *m_inputHitCollection.get();
- // Create output collections
- auto& proto = *m_outputProtoClusterCollection.createAndPut();
-
- // group neighboring hits
- std::vector<bool> visits(hits.size(), false);
- std::vector<std::vector<std::pair<uint32_t, eic::ConstCalorimeterHit>>> groups;
- for (size_t i = 0; i < hits.size(); ++i) {
- if (msgLevel(MSG::DEBUG)) {
- debug() << fmt::format("hit {:d}: local position = ({}, {}, {}), global position = ({}, {}, {})", i + 1,
- hits[i].local().x, hits[i].local().y, hits[i].position().z, hits[i].position().x,
- hits[i].position().y, hits[i].position().z)
- << endmsg;
- }
- // already in a group, or not energetic enough to form a cluster
- if (visits[i] || hits[i].energy() < minClusterCenterEdep) {
- continue;
- }
- groups.emplace_back();
- // create a new group, and group all the neighboring hits
- dfs_group(groups.back(), i, hits, visits);
+ // form clusters
+ for (const auto& group : groups) {
+ if (static_cast<int>(group.size()) < m_minClusterNhits.value()) {
+ continue;
}
- if (msgLevel(MSG::DEBUG)) {
- debug() << "found " << groups.size() << " potential clusters (groups of hits)" << endmsg;
- for (size_t i = 0; i < groups.size(); ++i) {
- debug() << fmt::format("group {}: {} hits", i, groups[i].size()) << endmsg;
- }
+ double energy = 0.;
+ for (const auto& [idx, hit] : group) {
+ energy += hit.energy();
}
-
- // form clusters
- for (const auto& group : groups) {
- if (static_cast<int>(group.size()) < m_minClusterNhits.value()) {
- continue;
- }
- double energy = 0.;
- for (const auto& [idx, hit] : group) {
- energy += hit.energy();
- }
- if (energy < minClusterEdep) {
- continue;
- }
- auto pcl = proto.create();
- for (const auto& [idx, hit] : group) {
- pcl.addhits(hit);
- }
+ if (energy < minClusterEdep) {
+ continue;
+ }
+ auto pcl = proto.create();
+ for (const auto& [idx, hit] : group) {
+ pcl.addhits(hit);
}
-
- return StatusCode::SUCCESS;
}
- private:
- template <typename T>
- static inline T pow2(const T& x)
- {
- return x * x;
- }
+ return StatusCode::SUCCESS;
+ }
- // helper function to group hits
- bool is_neighbor(const eic::ConstCalorimeterHit& h1, const eic::ConstCalorimeterHit& h2) const
- {
- // different sectors, simple distance check
- if (h1.sector() != h2.sector()) {
- return std::sqrt(pow2(h1.position().x - h2.position().x) + pow2(h1.position().y - h2.position().y) + pow2(h1.position().z - h2.position().z)) <=
- sectorDist;
- }
+private:
+ template <typename T> static inline T pow2(const T& x) { return x * x; }
- // layer check
- int ldiff = std::abs(h1.layer() - h2.layer());
- // same layer, check local positions
- if (!ldiff) {
- return (std::abs(h1.local().x - h2.local().x) <= localDistXY[0]) &&
- (std::abs(h1.local().y - h2.local().y) <= localDistXY[1]);
- } else if (ldiff <= m_neighbourLayersRange) {
- return (std::abs(h1.position().eta() - h2.position().eta()) <= layerDistEtaPhi[0]) &&
- (std::abs(h1.position().phi() - h2.position().phi()) <= layerDistEtaPhi[1]);
- }
+ // helper function to group hits
+ bool is_neighbor(const eic::ConstCalorimeterHit& h1, const eic::ConstCalorimeterHit& h2) const {
+ // different sectors, simple distance check
+ if (h1.sector() != h2.sector()) {
+ return std::sqrt(pow2(h1.position().x - h2.position().x) + pow2(h1.position().y - h2.position().y) +
+ pow2(h1.position().z - h2.position().z)) <= sectorDist;
+ }
- // not in adjacent layers
- return false;
+ // layer check
+ int ldiff = std::abs(h1.layer() - h2.layer());
+ // same layer, check local positions
+ if (!ldiff) {
+ return (std::abs(h1.local().x - h2.local().x) <= localDistXY[0]) &&
+ (std::abs(h1.local().y - h2.local().y) <= localDistXY[1]);
+ } else if (ldiff <= m_neighbourLayersRange) {
+ return (std::abs(eicd::eta(h1.position()) - eicd::eta(h2.position())) <= layerDistEtaPhi[0]) &&
+ (std::abs(eicd::angleAzimuthal(h1.position()) - eicd::angleAzimuthal(h2.position())) <=
+ layerDistEtaPhi[1]);
}
- // grouping function with Depth-First Search
- void dfs_group(std::vector<std::pair<uint32_t, eic::ConstCalorimeterHit>>& group, int idx,
- const eic::CalorimeterHitCollection& hits, std::vector<bool>& visits) const
- {
- // not a qualified hit to participate in clustering, stop here
- if (hits[idx].energy() < minClusterHitEdep) {
- visits[idx] = true;
- return;
- }
+ // not in adjacent layers
+ return false;
+ }
- group.push_back({idx, hits[idx]});
+ // grouping function with Depth-First Search
+ void dfs_group(std::vector<std::pair<uint32_t, eic::ConstCalorimeterHit>>& group, int idx,
+ const eic::CalorimeterHitCollection& hits, std::vector<bool>& visits) const {
+ // not a qualified hit to participate in clustering, stop here
+ if (hits[idx].energy() < minClusterHitEdep) {
visits[idx] = true;
- for (size_t i = 0; i < hits.size(); ++i) {
- // visited, or not a neighbor
- if (visits[i] || !is_neighbor(hits[idx], hits[i])) {
- continue;
- }
- dfs_group(group, i, hits, visits);
+ return;
+ }
+
+ group.push_back({idx, hits[idx]});
+ visits[idx] = true;
+ for (size_t i = 0; i < hits.size(); ++i) {
+ // visited, or not a neighbor
+ if (visits[i] || !is_neighbor(hits[idx], hits[i])) {
+ continue;
}
+ dfs_group(group, i, hits, visits);
}
- }; // namespace Jug::Reco
+ }
+}; // namespace Jug::Reco
- DECLARE_COMPONENT(ImagingTopoCluster)
+DECLARE_COMPONENT(ImagingTopoCluster)
} // namespace Jug::Reco
diff --git a/JugReco/src/components/SimpleClustering.cpp b/JugReco/src/components/SimpleClustering.cpp
index 671eac3bd79cb10d2626bbbc889224d732859e87..cd14b39cccabb534f76394e384984c0a214cf505 100644
--- a/JugReco/src/components/SimpleClustering.cpp
+++ b/JugReco/src/components/SimpleClustering.cpp
@@ -21,6 +21,7 @@
#include "eicd/ClusterCollection.h"
#include "eicd/ProtoClusterCollection.h"
#include "eicd/RawCalorimeterHitCollection.h"
+#include "eicd/vector_utils.h"
using namespace Gaudi::Units;
@@ -122,7 +123,7 @@ namespace Jug::Reco {
std::vector<std::pair<uint32_t, eic::ConstCalorimeterHit>> cluster_hits;
for (const auto& [idx, h] : the_hits) {
- if (h.position().subtract(ref_hit.position()).mag() < max_dist) {
+ if (eicd::magnitude(h.position() - ref_hit.position()) < max_dist) {
cluster_hits.push_back({idx, h});
} else {
remaining_hits.push_back({idx, h});
@@ -142,7 +143,7 @@ namespace Jug::Reco {
auto pcl = proto.create();
for (const auto& [idx, h] : cluster_hits) {
cl.energy(cl.energy() + h.energy());
- cl.position(cl.position().add(h.position().scale(h.energy()/total_energy)));
+ cl.position(cl.position() + (h.position() * h.energy() / total_energy));
pcl.addhits(h);
}
// Optionally store the MC truth associated with the first hit in this cluster
diff --git a/JugTrack/src/components/GenFitTrackFitter.cpp b/JugTrack/src/components/GenFitTrackFitter.cpp
index c31d2a91976c02839ffa10dba83ed3f3b2ac9a41..3169faa6928c127198967a760a2c3829f4cde9ce 100644
--- a/JugTrack/src/components/GenFitTrackFitter.cpp
+++ b/JugTrack/src/components/GenFitTrackFitter.cpp
@@ -20,6 +20,7 @@
#include "JugTrack/Track.hpp"
#include "eicd/TrackerHitCollection.h"
+#include "eicd/vector_utils.h"
#include <functional>
#include <random>
@@ -129,7 +130,7 @@ namespace Jug::Reco {
ROOT::Math::XYZVector tp(track_param.momentum()[0], track_param.momentum()[1], track_param.momentum()[2]);
auto first_hit = (*hits)[proto_track[0]];
- auto first_hit_phi = first_hit.position().phi();
+ auto first_hit_phi = eicd::angleAzimuthal(first_hit.position());
auto track_param_phi = tp.phi();
if (msgLevel(MSG::DEBUG)) {
debug() << " first hit phi: " << first_hit_phi << endmsg;
diff --git a/JugTrack/src/components/ParticlesFromTrackFit.cpp b/JugTrack/src/components/ParticlesFromTrackFit.cpp
index 49803460ddaea7ff80a98a25e7c570acdc8c92e7..02fc59100ef927f238ff15dc8397ab0cb22177b6 100644
--- a/JugTrack/src/components/ParticlesFromTrackFit.cpp
+++ b/JugTrack/src/components/ParticlesFromTrackFit.cpp
@@ -116,14 +116,14 @@ namespace Jug::Reco {
debug() << " chi2 = " << trajState.chi2Sum << endmsg;
}
- const eic::CovXYZ covMomentum {
+ const eicd::Cov3f covMomentum {
static_cast<float>(covariance(Acts::eBoundTheta, Acts::eBoundTheta)),
static_cast<float>(covariance(Acts::eBoundPhi, Acts::eBoundPhi)),
static_cast<float>(covariance(Acts::eBoundQOverP, Acts::eBoundQOverP)),
static_cast<float>(covariance(Acts::eBoundTheta, Acts::eBoundPhi)),
static_cast<float>(covariance(Acts::eBoundTheta, Acts::eBoundQOverP)),
static_cast<float>(covariance(Acts::eBoundPhi, Acts::eBoundQOverP))};
- const eic::CovXY covPos {
+ const eicd::Cov2f covPos {
static_cast<float>(covariance(Acts::eBoundLoc0, Acts::eBoundLoc0)),
static_cast<float>(covariance(Acts::eBoundLoc1, Acts::eBoundLoc1)),
static_cast<float>(covariance(Acts::eBoundLoc0, Acts::eBoundLoc1))};
diff --git a/JugTrack/src/components/TrackParamImagingClusterInit.cpp b/JugTrack/src/components/TrackParamImagingClusterInit.cpp
index f6eb347c997c720948ada935b497fc25bb0942e0..dd35c7689cc31b9aa87932c8f91aa1bb38b5c91f 100644
--- a/JugTrack/src/components/TrackParamImagingClusterInit.cpp
+++ b/JugTrack/src/components/TrackParamImagingClusterInit.cpp
@@ -15,8 +15,8 @@
#include "eicd/TrackerHitCollection.h"
#include "eicd/ClusterCollection.h"
+#include "eicd/vector_utils.h"
-#include "Math/Vector3D.h"
#include "Acts/Surfaces/PerigeeSurface.hpp"
///// (Reconstructed) track parameters e.g. close to the vertex.
@@ -76,28 +76,19 @@ namespace Jug::Reco {
using Acts::UnitConstants::mm;
using Acts::UnitConstants::ns;
- double p = c.energy()*GeV;
+ const double p = c.energy()*GeV;
// FIXME hardcoded value
if( p < 0.1*GeV) {
continue;
}
- double len = c.position().mag();
- ROOT::Math::XYZVector momentum(c.position().x * p / len, c.position().y * p / len, c.position().z * p / len);
-
- // build some track cov matrix
- //Acts::BoundSymMatrix cov = Acts::BoundSymMatrix::Zero();
- //cov(Acts::eBoundLoc0, Acts::eBoundLoc0) = 1.0 * mm*1.0 * mm;
- //cov(Acts::eBoundLoc1, Acts::eBoundLoc1) = 1.0 * mm*1.0 * mm;
- //cov(Acts::eBoundPhi, Acts::eBoundPhi) = M_PI / 180.0;
- //cov(Acts::eBoundTheta, Acts::eBoundTheta) = M_PI / 180.0;
- //cov(Acts::eBoundQOverP, Acts::eBoundQOverP) = 1.0 / (p*p);
- //cov(Acts::eBoundTime, Acts::eBoundTime) = Acts::UnitConstants::ns;
+ const double theta = eicd::anglePolar(c.position());
+ const double phi = eicd::angleAzimuthal(c.position());
Acts::BoundVector params;
params(Acts::eBoundLoc0) = 0.0 * mm ;
params(Acts::eBoundLoc1) = 0.0 * mm ;
- params(Acts::eBoundPhi) = momentum.Phi();
- params(Acts::eBoundTheta) = momentum.Theta();
+ params(Acts::eBoundPhi) = phi;
+ params(Acts::eBoundTheta) = theta;
params(Acts::eBoundQOverP) = 1/p;
params(Acts::eBoundTime) = 0 * ns;
@@ -111,24 +102,12 @@ namespace Jug::Reco {
Acts::BoundVector params2;
params2(Acts::eBoundLoc0) = 0.0 * mm ;
params2(Acts::eBoundLoc1) = 0.0 * mm ;
- params2(Acts::eBoundPhi) = momentum.Phi();
- params2(Acts::eBoundTheta) = momentum.Theta();
+ params2(Acts::eBoundPhi) = phi;
+ params2(Acts::eBoundTheta) = theta;
params2(Acts::eBoundQOverP) = -1/p;
params2(Acts::eBoundTime) = 0 * ns;
init_trk_params->push_back({pSurface, params2, -1});
- // acts v1.2.0:
- //init_trk_params->emplace_back(Acts::Vector4(0 * mm, 0 * mm, 0 * mm, 0),
- // Acts::Vector3(c.x() * p / len, c.y() * p / len, c.z() * p / len), p, -1,
- // std::make_optional(cov));
- //debug() << init_trk_params->back() << endmsg;
- //init_trk_params->emplace_back(Acts::Vector4(0 * mm, 0 * mm, 0 * mm, 0),
- // Acts::Vector3(c.x() * p / len, c.y() * p / len, c.z() * p / len), p, 1,
- // std::make_optional(cov));
- ////debug() << init_trk_params->back() << endmsg;
- //init_trk_params->emplace_back(Acts::Vector4(0 * mm, 0 * mm, 0 * mm, 0),
- // Acts::Vector3(c.x() * p / len, c.y() * p / len, c.z() * p / len), p, 0,
- // std::make_optional(cov));
}
return StatusCode::SUCCESS;
}