Skip to content
GitLab
Commit 2289a999 authored by Sylvester Joosten's avatar Sylvester Joosten
Browse files

Update for new EICD components

parent b0b64dc6
Hide whitespace changes
Inline Side-by-side
JugReco/src/components/CalorimeterHitsEtaPhiProjector.cpp
View file @ 2289a999
......@@ -31,6 +31,7 @@
// Event Model related classes
#include "eicd/CalorimeterHitCollection.h"
#include "eicd/VectorPolar.h"
using namespace Gaudi::Units;
typedef ROOT::Math::XYZPoint Point3D;
......@@ -97,9 +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()) {
Point3D p(h.position().x, h.position().y, h.position().z);
auto bins = std::make_pair(static_cast<int64_t>(pos2bin(p.eta(), gridSizes[0], 0.)),
static_cast<int64_t>(pos2bin(p.phi(), gridSizes[1], 0.)));
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.)));
merged_hits[bins].push_back(h);
}
......@@ -111,14 +111,12 @@ namespace Jug::Reco {
hit.layerID(ref.layerID());
// TODO, we can do timing cut to reject noises
hit.time(ref.time());
double r = std::hypot(ref.position().x, ref.position().y, ref.position().z);
double r = ref.position().mag();
double eta = bin2pos(bins.first, gridSizes[0], 0.);
double phi = bin2pos(bins.second, gridSizes[1], 1.);
double theta = std::atan(std::exp(-eta))*2.;
ROOT::Math::Polar3DVector p(r, theta, phi);
// TODO, we can do global to local conversion, which gives a better estimate
hit.local(ref.local());
hit.position({p.x(), p.y(), p.z()});
hit.position(eic::VectorPolar(r, theta, phi));
hit.dimension({gridSizes[0], gridSizes[1], 0.});
// merge energy
hit.energy(0.);
......
JugReco/src/components/CalorimeterIslandCluster.cpp
View file @ 2289a999
......@@ -14,8 +14,6 @@
#include <tuple>
#include "fmt/format.h"
#include "Math/Point2D.h"
#include "Math/Point3D.h"
#include "Gaudi/Property.h"
#include "GaudiAlg/GaudiAlgorithm.h"
......@@ -36,38 +34,37 @@
// Event Model related classes
#include "eicd/CalorimeterHitCollection.h"
#include "eicd/ClusterCollection.h"
#include "eicd/VectorXY.h"
#include "eicd/VectorXYZ.h"
using namespace Gaudi::Units;
typedef ROOT::Math::XYPoint Point;
typedef ROOT::Math::XYZPoint Point3D;
namespace Jug::Reco {
// helper functions to get distance between hits
static Point localDistXY(eic::ConstCalorimeterHit h1, eic::ConstCalorimeterHit h2) {
return Point(h1.local().local_x - h2.local().local_x, h1.local().local_y - h2.local().local_y);
static eic::VectorXY localDistXY(eic::ConstCalorimeterHit h1, eic::ConstCalorimeterHit h2) {
return {h1.local().x - h2.local().x, h1.local().y - h2.local().y};
}
static Point localDistXZ(eic::ConstCalorimeterHit h1, eic::ConstCalorimeterHit h2) {
return Point(h1.local().local_x - h2.local().local_x, h1.local().local_z - h2.local().local_z);
static eic::VectorXY localDistXZ(eic::ConstCalorimeterHit h1, eic::ConstCalorimeterHit h2) {
return {h1.local().x - h2.local().x, h1.local().z - h2.local().z};
}
static Point localDistYZ(eic::ConstCalorimeterHit h1, eic::ConstCalorimeterHit h2) {
return Point(h1.local().local_y - h2.local().local_y, h1.local().local_z - h2.local().local_z);
static eic::VectorXY localDistYZ(eic::ConstCalorimeterHit h1, eic::ConstCalorimeterHit h2) {
return {h1.local().y - h2.local().y, h1.local().z - h2.local().z};
}
static Point dimScaledLocalDistXY(eic::ConstCalorimeterHit h1, eic::ConstCalorimeterHit h2) {
return Point(2.*(h1.local().local_x - h2.local().local_x)/(h1.dimension().dim_x + h2.dimension().dim_x),
2.*(h1.local().local_y - h2.local().local_y)/(h1.dimension().dim_y + h2.dimension().dim_y));
static eic::VectorXY dimScaledLocalDistXY(eic::ConstCalorimeterHit h1, eic::ConstCalorimeterHit h2) {
return {2.*(h1.local().x - h2.local().x)/(h1.dimension().x + h2.dimension().x),
2.*(h1.local().y - h2.local().y)/(h1.dimension().y + h2.dimension().y)};
}
static Point globalDistRPhi(eic::ConstCalorimeterHit h1, eic::ConstCalorimeterHit h2) {
Point3D p1(h1.position().x, h1.position().y, h1.position().z), p2(h2.position().x, h2.position().y, h2.position().z);
return Point(p1.r() - p2.r(), p1.phi() - p2.phi());
static eic::VectorXY globalDistRPhi(eic::ConstCalorimeterHit h1, eic::ConstCalorimeterHit h2) {
return {h1.position().r() - h2.position().r(), h1.position().phi() - h2.position().phi()};
}
static Point globalDistEtaPhi(eic::ConstCalorimeterHit h1, eic::ConstCalorimeterHit h2) {
Point3D p1(h1.position().x, h1.position().y, h1.position().z), p2(h2.position().x, h2.position().y, h2.position().z);
return Point(p1.eta() - p2.eta(), p1.phi() - p2.phi());
static eic::VectorXY globalDistEtaPhi(eic::ConstCalorimeterHit h1, eic::ConstCalorimeterHit h2) {
return {h1.position().eta() - h2.position().eta(), h1.position().phi() - h2.position().phi()};
}
// name: {method, units}
static std::map< std::string, std::tuple<
std::function<Point(eic::ConstCalorimeterHit, eic::ConstCalorimeterHit)>,
std::vector<double> > > distMethods {
static std::map<std::string,
std::tuple<
std::function<eic::VectorXY(eic::ConstCalorimeterHit, eic::ConstCalorimeterHit)>,
std::vector<double>>> distMethods {
{"localDistXY", {localDistXY, {mm, mm}}},
{"localDistXZ", {localDistXZ, {mm, mm}}},
{"localDistYZ", {localDistYZ, {mm, mm}}},
......@@ -108,7 +105,7 @@ namespace Jug::Reco {
Gaudi::Property<std::vector<double>> u_globalDistEtaPhi{this, "globalDistEtaPhi", {}};
Gaudi::Property<std::vector<double>> u_dimScaledLocalDistXY{this, "dimScaledLocalDistXY", {1.8, 1.8}};
// neightbor checking function
std::function<Point(eic::ConstCalorimeterHit, eic::ConstCalorimeterHit)> hitsDist;
std::function<eic::VectorXY(eic::ConstCalorimeterHit, eic::ConstCalorimeterHit)> hitsDist;
// unitless counterparts of the input parameters
double minClusterHitEdep, minClusterCenterEdep, sectorDist;
......@@ -187,7 +184,7 @@ namespace Jug::Reco {
for (size_t i = 0; i < hits.size(); ++i) {
debug() << fmt::format("hit {:d}: energy = {:.4f} MeV, local = ({:.4f}, {:.4f}) mm, "
"global=({:.4f}, {:.4f}, {:.4f}) mm, layer = {:d}, sector = {:d}.",
i, hits[i].energy()*1000., hits[i].local().local_x, hits[i].local().local_y,
i, hits[i].energy()*1000., hits[i].local().x, hits[i].local().y,
hits[i].position().x, hits[i].position().y, hits[i].position().z,
hits[i].layerID(), hits[i].sectorID()) << endmsg;
// already in a group
......@@ -215,22 +212,17 @@ namespace Jug::Reco {
}
private:
template <typename T>
inline T dist3d(T t1, T t2, T t3) const
{
return std::sqrt(t1 * t1 + t2 * t2 + t3 * t3);
}
// helper function to group hits
inline bool is_neighbour(const eic::ConstCalorimeterHit& h1, const eic::ConstCalorimeterHit& h2) const
{
// in the same sector
if (h1.sectorID() == h2.sectorID()) {
auto dist = hitsDist(h1, h2);
return (dist.x() <= neighbourDist[0]) && (dist.y() <= neighbourDist[1]);
return (dist.x <= neighbourDist[0]) && (dist.y <= neighbourDist[1]);
// different sector, local coordinates do not work, using global coordinates
} else {
// sector may have rotation (barrel), so z is included
return dist3d(h1.position().x - h2.position().x, h1.position().y - h2.position().y, h1.position().z - h2.position().z) <= sectorDist;
return (h1.position().subtract(h2.position())).mag() <= sectorDist;
}
}
......@@ -347,9 +339,9 @@ namespace Jug::Reco {
size_t j = 0;
// calculate weights for local maxima
for (auto cit = maxima.begin(); cit != maxima.end(); ++cit, ++j) {
double dist_ref = cit->dimension().dim_x;
double dist_ref = cit->dimension().x;
double energy = cit->energy();
double dist = hitsDist(*cit, *it).r();
double dist = hitsDist(*cit, *it).mag();
weights[j] = std::exp(-dist / dist_ref) * energy;
}
......
JugReco/src/components/ClusterRecoCoG.cpp
View file @ 2289a999
......@@ -188,25 +188,19 @@ namespace Jug::Reco {
res.energy(totalE);
// center of gravity with logarithmic weighting
float tw = 0., x = 0., y = 0., z = 0.;
float tw = 0.;
eic::VectorXYZ v;
for (auto& hit : hits) {
// suppress low energy contributions
// info() << std::log(hit.energy()/totalE) << endmsg;
float w = weightFunc(hit.energy(), totalE, m_logWeightBase.value(), 0);
tw += w;
x += hit.position().x * w;
y += hit.position().y * w;
z += hit.position().z * w;
/*
debug() << hit.cellID() << ": (" << hit.local_x() << ", " << hit.local_y() << ", "
<< hit.local_z() << "), "
<< "(" << hit.position().x << ", " << hit.position().y << ", " << hit.position().z << "), " << endmsg;
*/
v = v.add(hit.position().scale(w));
}
if (tw == 0.) {
warning() << "zero total weights encountered, you may want to adjust your weighting parameter." << endmsg;
}
res.position({x / tw, y / tw, z / tw});
res.position(v.scale(1/tw));
// convert global position to local position, use the cell with max edep as a reference
const auto volman = m_geoSvc->detector()->volumeManager();
const auto alignment = volman.lookupDetElement(centerID).nominal();
......
JugReco/src/components/ImagingClusterReco.cpp
View file @ 2289a999
......@@ -185,7 +185,7 @@ namespace Jug::Reco {
double meta = 0.;
double mphi = 0.;
double edep = 0.;
double r = 9999 * cm;
float r = 9999 * cm;
for (const auto& hit : hits) {
meta += hit.eta() * hit.edep();
mphi += hit.polar().phi * hit.edep();
......
JugReco/src/components/ImagingTopoCluster.cpp
View file @ 2289a999
......@@ -135,7 +135,7 @@ namespace Jug::Reco {
/*
debug() << fmt::format("hit {:d}: local position = ({}, {}, {}), global position = ({}, {}, {})",
i + 1,
hits[i].local().local_x, hits[i].local().local_y, hits[i].local_z(),
hits[i].local().x, hits[i].local().y, hits[i].z(),
hits[i].position().x, hits[i].position().y, hits[i].position().z)
<< endmsg;
*/
......@@ -197,8 +197,8 @@ namespace Jug::Reco {
int ldiff = std::abs(h1.layerID() - h2.layerID());
// same layer, check local positions
if (!ldiff) {
return (std::abs(h1.local().local_x - h2.local().local_x) <= localDistXY[0]) &&
(std::abs(h1.local().local_y - h2.local().local_y) <= localDistXY[1]);
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.eta() - h2.eta()) <= layerDistEtaPhi[0]) &&
(std::abs(h1.polar().phi - h2.polar().phi) <= layerDistEtaPhi[1]);
......
JugReco/src/components/PhotoRingClusters.cpp
View file @ 2289a999
......@@ -87,7 +87,7 @@ namespace Jug::Reco {
MatrixXd data(rawhits.size(), 2);
for (int i = 0; i < data.rows(); ++i) {
if (rawhits[i].npe() > m_minNpe) {
data.row(i) << rawhits[i].local().local_x, rawhits[i].local().local_y;
data.row(i) << rawhits[i].local().x, rawhits[i].local().y;
}
}
......@@ -97,7 +97,7 @@ namespace Jug::Reco {
// local position
for (int i = 0; i < res.rows(); ++i) {
auto cl = clusters.create();
cl.position({res(i, 0), res(i, 1)});
cl.position({res(i, 0), res(i, 1), 0});
cl.radius(res(i, 2));
}
......
JugReco/src/components/TopologicalCellCluster.cpp
View file @ 2289a999
......@@ -198,8 +198,8 @@ namespace Jug::Reco {
int ldiff = std::abs(l1 - l2);
// same layer, check local positions
if (!ldiff) {
return (std::abs(h1.local().local_x - h2.local().local_x) <= u_localRanges[0]) &&
(std::abs(h1.local().local_y - h2.local().local_y) <= u_localRanges[1]);
return (std::abs(h1.local().x - h2.local().x) <= u_localRanges[0]) &&
(std::abs(h1.local().y - h2.local().y) <= u_localRanges[1]);
} else if (ldiff <= m_adjLayerDiff) {
double eta1, phi1, r1, eta2, phi2, r2;
calc_eta_phi_r(h1.position().x, h1.position().y, h1.position().z, eta1, phi1, r1);
......
JugTrack/src/components/GenFitTrackFitter.cpp
View file @ 2289a999
......@@ -85,7 +85,7 @@ namespace Jug::Reco {
// Read input data
const eic::TrackerHitCollection* hits = m_inputHitCollection.get();
const TrackParametersContainer* initialParameters = m_initialTrackParameters.get();
const ProtoTrackContainer* protoTracks = m_inputProtoTracks.get();
//const ProtoTrackContainer* protoTracks = m_inputProtoTracks.get();
const auto& single_track_param = (*initialParameters)[0];
// TrajectoryContainer trajectories;
......@@ -127,8 +127,8 @@ namespace Jug::Reco {
auto vol_id = vol_ctx->identifier;
TMatrixDSym hitCov(2);
hitCov.UnitMatrix();
hitCov(0, 0) = ahit.covMatrix().covsym_xx * ahit.covMatrix().covsym_xx;
hitCov(1, 1) = ahit.covMatrix().covsym_yy * ahit.covMatrix().covsym_yy;
hitCov(0, 0) = ahit.covMatrix().xx * ahit.covMatrix().xx;
hitCov(1, 1) = ahit.covMatrix().yy * ahit.covMatrix().yy;
TVector3 point = {ahit.position().x, ahit.position().y, ahit.position().z};
TVector3 u_dir = {1, 0, 0};
......
JugTrack/src/components/ParticlesFromTrackFit.cpp
View file @ 2289a999
......@@ -83,8 +83,8 @@ namespace Jug::Reco {
// Collect the trajectory summary info
auto trajState = Acts::MultiTrajectoryHelpers::trajectoryState(mj, trackTip);
int m_nMeasurements = trajState.nMeasurements;
int m_nStates = trajState.nStates;
//int m_nMeasurements = trajState.nMeasurements;
//int m_nStates = trajState.nStates;
// Get the fitted track parameter
bool m_hasFittedParams = false;
......@@ -148,9 +148,11 @@ namespace Jug::Reco {
return;
}
eic::Particle p({params[Acts::eBoundPhi], params[Acts::eBoundTheta],
1.0 / std::abs(params[Acts::eBoundQOverP]), 0.000511},
eic::Particle p({params[Acts::eBoundPhi],
params[Acts::eBoundTheta],
1.0 / std::abs(params[Acts::eBoundQOverP])},
{0.0, 0.0, 0.0, params[Acts::eBoundTime]},
0.000511,
(long long)11 * params[Acts::eBoundQOverP] /
std::abs(params[Acts::eBoundQOverP]),
0);
......
JugTrack/src/components/SingleTrackSourceLinker.cpp
View file @ 2289a999
......@@ -130,9 +130,9 @@ namespace Jug::Reco {
// construct the covariance matrix
Acts::SymMatrix2 cov = Acts::SymMatrix2::Zero();
cov(0, 0) =
ahit.covMatrix().covsym_xx * Acts::UnitConstants::mm * ahit.covMatrix().covsym_xx * Acts::UnitConstants::mm;
ahit.covMatrix().xx * Acts::UnitConstants::mm * ahit.covMatrix().xx * Acts::UnitConstants::mm;
cov(1, 1) =
ahit.covMatrix().covsym_yy * Acts::UnitConstants::mm * ahit.covMatrix().covsym_yy * Acts::UnitConstants::mm;
ahit.covMatrix().yy * Acts::UnitConstants::mm * ahit.covMatrix().yy * Acts::UnitConstants::mm;
// Above we only consider the two position coordinates the comment below shows how to add time
// which we will probably want to try later.
......@@ -141,7 +141,6 @@ namespace Jug::Reco {
// cov << 0.05, 0., 0., 0., 0.05, 0., 0., 0., 900. * Acts::UnitConstants::ps * Acts::UnitConstants::ps;
// Acts::Vector3 par(localX, localY, simHit.time());
Index hitIdx = ihit;
IndexSourceLink sourceLink(surface->geometryId(), ihit);
auto meas =
Acts::makeMeasurement(sourceLink, pos, cov, Acts::eBoundLoc0, Acts::eBoundLoc1); //, Acts::eBoundTime);
......
JugTrack/src/components/Tracker2SourceLinker.cpp
View file @ 2289a999
......@@ -126,8 +126,8 @@ namespace Jug::Reco {
allHits->push_back(ahit.clone());
Acts::SymMatrix2 cov = Acts::SymMatrix2::Zero();
cov(0, 0) = ahit.covMatrix().covsym_xx * Acts::UnitConstants::mm * ahit.covMatrix().covsym_xx * Acts::UnitConstants::mm;
cov(1, 1) = ahit.covMatrix().covsym_yy * Acts::UnitConstants::mm * ahit.covMatrix().covsym_yy * Acts::UnitConstants::mm;
cov(0, 0) = ahit.covMatrix().xx * Acts::UnitConstants::mm * ahit.covMatrix().xx * Acts::UnitConstants::mm;
cov(1, 1) = ahit.covMatrix().yy * Acts::UnitConstants::mm * ahit.covMatrix().yy * Acts::UnitConstants::mm;
auto vol_ctx = m_geoSvc->cellIDPositionConverter()->findContext(ahit.cellID());
auto vol_id = vol_ctx->identifier;
......@@ -143,8 +143,8 @@ namespace Jug::Reco {
// //allHits->push_back(ahit.clone());
//
// Acts::SymMatrix2 cov = Acts::SymMatrix2::Zero();
// cov(0, 0) = ahit.covsym_xx() * Acts::UnitConstants::mm * ahit.covsym_xx() * Acts::UnitConstants::mm;
// cov(1, 1) = ahit.covsym_yy() * Acts::UnitConstants::mm * ahit.covsym_yy() * Acts::UnitConstants::mm;
// cov(0, 0) = ahit.xx() * Acts::UnitConstants::mm * ahit.xx() * Acts::UnitConstants::mm;
// cov(1, 1) = ahit.yy() * Acts::UnitConstants::mm * ahit.yy() * Acts::UnitConstants::mm;
//
// auto vol_ctx = m_geoSvc->cellIDPositionConverter()->findContext(ahit.cellID());
// auto vol_id = vol_ctx->identifier;
......@@ -183,7 +183,6 @@ namespace Jug::Reco {
// the measurement container is unordered and the index under which the
// measurement will be stored is known before adding it.
Index hitIdx = measurements->size();
IndexSourceLink sourceLink(vol_id, ihit);
auto meas = Acts::makeMeasurement(sourceLink, loc, cov, Acts::eBoundLoc0, Acts::eBoundLoc1);
......
JugTrack/src/components/TrackerSourceLinker.cpp
View file @ 2289a999
......@@ -82,8 +82,8 @@ namespace Jug::Reco {
for(const auto& ahit : *hits) {
Acts::SymMatrix2 cov = Acts::SymMatrix2::Zero();
cov(0,0) = ahit.covMatrix().covsym_xx*Acts::UnitConstants::mm;//*ahit.covsym_xx()*Acts::UnitConstants::mm;
cov(1,1) = ahit.covMatrix().covsym_yy*Acts::UnitConstants::mm;//*ahit.covsym_yy()*Acts::UnitConstants::mm;
cov(0,0) = ahit.covMatrix().xx*Acts::UnitConstants::mm;//*ahit.xx()*Acts::UnitConstants::mm;
cov(1,1) = ahit.covMatrix().yy*Acts::UnitConstants::mm;//*ahit.yy()*Acts::UnitConstants::mm;
auto vol_ctx = m_geoSvc->cellIDPositionConverter()->findContext(ahit.cellID());
auto vol_id = vol_ctx->identifier;
......
JugTrack/src/components/TrackerSourcesLinker.cpp
View file @ 2289a999
......@@ -90,8 +90,8 @@ namespace Jug::Reco {
for(const auto col : u_hitCollections) {
for (const auto& ahit : *col->get()) {
Acts::SymMatrix2 cov = Acts::SymMatrix2::Zero();
cov(0, 0) = ahit.covMatrix().covsym_xx * Acts::UnitConstants::mm * ahit.covMatrix().covsym_xx * Acts::UnitConstants::mm;
cov(1, 1) = ahit.covMatrix().covsym_yy * Acts::UnitConstants::mm * ahit.covMatrix().covsym_yy * Acts::UnitConstants::mm;
cov(0, 0) = ahit.covMatrix().xx * Acts::UnitConstants::mm * ahit.covMatrix().xx * Acts::UnitConstants::mm;
cov(1, 1) = ahit.covMatrix().yy * Acts::UnitConstants::mm * ahit.covMatrix().yy * Acts::UnitConstants::mm;
auto vol_ctx = m_geoSvc->cellIDPositionConverter()->findContext(ahit.cellID());
auto vol_id = vol_ctx->identifier;
......
JugTrack/src/components/TrackingHitsSourceLinker.cpp
View file @ 2289a999
......@@ -109,9 +109,9 @@ namespace Jug::Reco {
// setup local covariance
Acts::BoundMatrix cov = Acts::BoundMatrix::Zero();
cov(Acts::eBoundLoc0, Acts::eBoundLoc0) =
ahit.covsym_xx() * Acts::UnitConstants::mm * ahit.covsym_xx() * Acts::UnitConstants::mm;
ahit.xx() * Acts::UnitConstants::mm * ahit.xx() * Acts::UnitConstants::mm;
cov(Acts::eBoundLoc1, Acts::eBoundLoc1) =
ahit.covsym_yy() * Acts::UnitConstants::mm * ahit.covsym_yy() * Acts::UnitConstants::mm;
ahit.yy() * Acts::UnitConstants::mm * ahit.yy() * Acts::UnitConstants::mm;
auto vol_ctx = m_geoSvc->cellIDPositionConverter()->findContext(ahit.cellID());
auto vol_id = vol_ctx->identifier;
......
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment