Skip to content

GitLab

Commit 0ff24fa1 authored by Sylvester Joosten's avatar Sylvester Joosten
Browse files

Remove out VectorXY? --> Vector2/3f/d 

Swap out Covariance matrix structures

remove eic::Weight
parent 3ecf1202
Hide whitespace changes
Inline Side-by-side
JugDigi/src/components/CalorimeterHitDigi.cpp
View file @ 0ff24fa1
......@@ -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();
......
JugFast/src/components/ParticlesWithTruthPID.cpp
View file @ 0ff24fa1
......@@ -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;
}
}
......
JugFast/src/components/SmearedFarForwardParticles.cpp
View file @ 0ff24fa1
......@@ -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;
}
}
......
JugReco/src/components/CalorimeterHitsEtaPhiProjector.cpp
View file @ 0ff24fa1
......@@ -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) {
......
JugReco/src/components/CalorimeterIslandCluster.cpp
View file @ 0ff24fa1
......@@ -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;
}
......
JugReco/src/components/ClusterRecoCoG.cpp
View file @ 0ff24fa1
......@@ -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
......
JugReco/src/components/EnergyPositionClusterMerger.cpp
View file @ 0ff24fa1
......@@ -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;
}
}
......
JugReco/src/components/FarForwardParticles.cpp
View file @ 0ff24fa1
......@@ -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});