diff --git a/JugDigi/src/components/CalorimeterHitDigi.cpp b/JugDigi/src/components/CalorimeterHitDigi.cpp
index f7ba9a1b9ea9e43f264140ee08cf548a30fa1b56..9ae76a5424dc334e96750a1ba042ea47832fa624 100644
--- a/JugDigi/src/components/CalorimeterHitDigi.cpp
+++ b/JugDigi/src/components/CalorimeterHitDigi.cpp
@@ -33,10 +33,6 @@ public:
Gaudi::Property<std::vector<double>> u_eRes{this, "energyResolutions", {}}; // a/sqrt(E/GeV) + b + c/(E/GeV)
Gaudi::Property<double> m_tRes{this, "timineResolution", 0.0*ns};
- // input units, should be fixed
- Gaudi::Property<double> m_eUnit{this, "inputEnergyUnit", GeV};
- Gaudi::Property<double> m_tUnit{this, "inputTimeUnit", ns};
-
// digitization settings
Gaudi::Property<int> m_capADC{this, "capacityADC", 8096};
Gaudi::Property<double> m_dyRangeADC{this, "dynamicRangeADC", 100*MeV};
@@ -48,7 +44,8 @@ public:
m_inputHitCollection{"inputHitCollection", Gaudi::DataHandle::Reader, this};
DataHandle<eic::RawCalorimeterHitCollection>
m_outputHitCollection{"outputHitCollection", Gaudi::DataHandle::Writer, this};
- double res[3] = {0., 0., 0.};
+ // unitless counterparts of inputs
+ double dyRangeADC, tRes, eRes[3] = {0., 0., 0.};
// ill-formed: using GaudiAlgorithm::GaudiAlgorithm;
CalorimeterHitDigi(const std::string& name, ISvcLocator* svcLoc) : GaudiAlgorithm(name, svcLoc)
@@ -70,8 +67,13 @@ public:
}
// set energy resolution numbers
for (size_t i = 0; i < u_eRes.size() && i < 3; ++i) {
- res[i] = u_eRes[i];
+ eRes[i] = u_eRes[i];
}
+
+ // using juggler internal units (GeV, mm, radian, ns)
+ dyRangeADC = m_dyRangeADC.value()/GeV;
+ tRes = m_tRes.value()/ns;
+
return StatusCode::SUCCESS;
}
@@ -82,15 +84,16 @@ public:
// Create output collections
auto rawhits = m_outputHitCollection.createAndPut();
for (const auto& ahit : *simhits) {
- double eres = std::sqrt(std::pow(m_normDist()*res[0] / sqrt(ahit.energyDeposit()*m_eUnit/GeV), 2)
- + std::pow(m_normDist()*res[1], 2)
- + std::pow(m_normDist()*res[2] / (ahit.energyDeposit()*m_eUnit/GeV), 2));
+ // Note: juggler internal unit of energy is GeV
+ double eResRel = std::sqrt(std::pow(m_normDist()*eRes[0] / sqrt(ahit.energyDeposit()), 2)
+ + std::pow(m_normDist()*eRes[1], 2)
+ + std::pow(m_normDist()*eRes[2] / (ahit.energyDeposit()), 2));
double ped = m_pedMeanADC + m_normDist()*m_pedSigmaADC;
- long long adc = std::llround(ped + ahit.energyDeposit()*(1. + eres) * m_eUnit/m_dyRangeADC*m_capADC);
+ long long adc = std::llround(ped + ahit.energyDeposit()*(1. + eResRel)/dyRangeADC*m_capADC);
eic::RawCalorimeterHit rawhit(
(long long) ahit.cellID(),
- (adc > m_capADC ? m_capADC.value() : adc),
- (double) ahit.truth().time*m_tUnit/ns + m_normDist()*m_tRes/ns
+ (adc > m_capADC.value() ? m_capADC.value() : adc),
+ (double) ahit.truth().time + m_normDist()*tRes
);
rawhits->push_back(rawhit);
}
diff --git a/JugReco/src/components/CalorimeterHitReco.cpp b/JugReco/src/components/CalorimeterHitReco.cpp
index 2bea3613338027097519e684e11c7592b0100edc..c5b849d8d12d8f65deb85a716d8736dfcb628b68 100644
--- a/JugReco/src/components/CalorimeterHitReco.cpp
+++ b/JugReco/src/components/CalorimeterHitReco.cpp
@@ -44,7 +44,9 @@ namespace Jug::Reco {
// zero suppression values
Gaudi::Property<double> m_thresholdADC{this, "thresholdFactor", 3.0};
- Gaudi::Property<double> m_minEdep{this, "minimumHitEdep", 0.};
+
+ // unitless counterparts of the input parameters
+ double dyRangeADC;
DataHandle<eic::RawCalorimeterHitCollection> m_inputHitCollection{
"inputHitCollection", Gaudi::DataHandle::Reader, this};
@@ -87,6 +89,9 @@ namespace Jug::Reco {
return StatusCode::FAILURE;
}
+ // unitless conversion
+ dyRangeADC = m_dyRangeADC.value()/GeV;
+
// do not get the layer/sector ID if no readout class provided
if (m_readout.value().empty()) {
return StatusCode::SUCCESS;
@@ -97,9 +102,11 @@ namespace Jug::Reco {
id_dec = id_spec.decoder();
if (m_sectorField.value().size()) {
sector_idx = id_dec->index(m_sectorField);
+ info() << "Find sector field " << m_sectorField.value() << ", index = " << sector_idx << endmsg;
}
if (m_layerField.value().size()) {
layer_idx = id_dec->index(m_layerField);
+ info() << "Find layer field " << m_layerField.value() << ", index = " << sector_idx << endmsg;
}
} catch (...) {
error() << "Failed to load ID decoder for " << m_readout << endmsg;
@@ -117,7 +124,7 @@ namespace Jug::Reco {
<< 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()) {
@@ -151,20 +158,16 @@ namespace Jug::Reco {
}
// convert ADC -> energy
- float energy = (rh.amplitude() - m_pedMeanADC) / (float)m_capADC * m_dyRangeADC;
- if (energy < m_minEdep) {
- continue;
- }
+ float energy = (rh.amplitude() - m_pedMeanADC) / static_cast<float>(m_capADC.value()) * dyRangeADC;
float time = rh.timeStamp(); // ns
auto id = rh.cellID();
- int lid = ((id_dec != nullptr) & m_layerField.value().size())
+ int lid = ((id_dec != nullptr) && m_layerField.value().size())
? static_cast<int>(id_dec->get(id, layer_idx))
: -1;
- int sid = ((id_dec != nullptr) & m_sectorField.value().size())
+ int sid = ((id_dec != nullptr) && m_sectorField.value().size())
? static_cast<int>(id_dec->get(id, sector_idx))
: -1;
-
// global positions
auto gpos = m_geoSvc->cellIDPositionConverter()->position(id);
// local positions
@@ -173,9 +176,8 @@ namespace Jug::Reco {
local = volman.lookupDetElement(id & local_mask).nominal();
}
auto pos = local.worldToLocal(dd4hep::Position(gpos.x(), gpos.y(), gpos.z()));
- // auto pos =
- // m_geoSvc->cellIDPositionConverter()->findContext(id)->volumePlacement().position(); cell
- // dimension
+ // auto pos = m_geoSvc->cellIDPositionConverter()->findContext(id)->volumePlacement().position();
+ // cell dimension
auto cdim = m_geoSvc->cellIDPositionConverter()->cellDimensions(id);
double dim[3] = {0., 0., 0.};
for (size_t i = 0; i < cdim.size() && i < 3; ++i) {
diff --git a/JugReco/src/components/CalorimeterHitsMerger.cpp b/JugReco/src/components/CalorimeterHitsMerger.cpp
index 0d65598c88f9393f94639c843137e1e20adf3e52..f6340b40d6c7a2559fa0303be313a89387e22e02 100644
--- a/JugReco/src/components/CalorimeterHitsMerger.cpp
+++ b/JugReco/src/components/CalorimeterHitsMerger.cpp
@@ -4,8 +4,9 @@
*
* Author: Chao Peng (ANL), 03/31/2021
*/
-#include <algorithm>
+#include <tuple>
#include <bitset>
+#include <algorithm>
#include <unordered_map>
#include "Gaudi/Property.h"
@@ -102,50 +103,60 @@ namespace Jug::Reco {
StatusCode execute() override
{
// input collections
- const auto& hits = *m_inputHitCollection.get();
+ const auto& inputs = *m_inputHitCollection.get();
// Create output collections
- auto& mhits = *m_outputHitCollection.createAndPut();
-
- // dd4hep decoders
- auto poscon = m_geoSvc->cellIDPositionConverter();
- auto volman = m_geoSvc->detector()->volumeManager();
+ auto& outputs = *m_outputHitCollection.createAndPut();
- // sum energies that has the same id
- std::unordered_map<long long, size_t> merge_map;
- for (const auto& h : hits) {
+ // find the hits that belong to the same group (for merging)
+ std::unordered_map<long long, std::vector<eic::ConstCalorimeterHit>> merge_map;
+ for (const auto& h : inputs) {
int64_t id = h.cellID() & id_mask;
// use the reference field position
- int64_t ref_id = id | ref_mask;
- // debug() << h.cellID() << " - " << std::bitset<64>(h.cellID()) << endmsg;
auto it = merge_map.find(id);
- debug() << fmt::format("{:#064b} - {:#064b}, ref: {:#064b}", h.cellID(), id, ref_id)
- << endmsg;
if (it == merge_map.end()) {
- merge_map[id] = mhits.size();
- // global positions
- auto gpos = poscon->position(ref_id);
- // local positions
- auto alignment = volman.lookupDetElement(ref_id).nominal();
- // debug() << volman.lookupDetElement(ref_id).path() << ", " <<
- // volman.lookupDetector(ref_id).path() << endmsg;
- auto pos = alignment.worldToLocal(dd4hep::Position(gpos.x(), gpos.y(), gpos.z()));
-
- mhits.push_back({ref_id,
- h.clusterID(),
- h.layerID(),
- h.sectorID(),
- h.energy(),
- h.time(),
- {gpos.x() / dd4hep::mm, gpos.y() / dd4hep::mm, gpos.z() / dd4hep::mm},
- {pos.x() / dd4hep::mm, pos.y() / dd4hep::mm, pos.z() / dd4hep::mm},
- h.dimension(),
- h.type()});
+ merge_map[id] = {h};
} else {
- mhits[it->second].energy(mhits[it->second].energy() + h.energy());
+ it->second.push_back(h);
+ }
+ }
+
+ // reconstruct info for merged hits
+ // dd4hep decoders
+ auto poscon = m_geoSvc->cellIDPositionConverter();
+ auto volman = m_geoSvc->detector()->volumeManager();
+
+ for (auto &[id, hits] : merge_map) {
+ // reference fields id
+ int64_t ref_id = id | ref_mask;
+ // global positions
+ auto gpos = poscon->position(ref_id);
+ // local positions
+ auto alignment = volman.lookupDetElement(ref_id).nominal();
+ auto pos = alignment.worldToLocal(dd4hep::Position(gpos.x(), gpos.y(), gpos.z()));
+ // debug() << volman.lookupDetElement(ref_id).path() << ", "
+ // << volman.lookupDetector(ref_id).path() << endmsg;
+ // sum energy
+ float energy = 0.;
+ for (auto &hit : hits) {
+ energy += hit.energy();
+ // debug() << fmt::format("{:#064b} - {:#064b}, ref: {:#064b}", hit.cellID(), id, ref_id)
+ // << endmsg;
}
+ const auto &href = hits.front();
+ outputs.push_back(eic::CalorimeterHit{
+ ref_id,
+ href.clusterID(),
+ href.layerID(),
+ href.sectorID(),
+ energy,
+ href.time(),
+ {gpos.x() / dd4hep::mm, gpos.y() / dd4hep::mm, gpos.z() / dd4hep::mm},
+ {pos.x() / dd4hep::mm, pos.y() / dd4hep::mm, pos.z() / dd4hep::mm},
+ href.dimension(),
+ href.type()});
}
- debug() << "Size before = " << hits.size() << ", after = " << mhits.size() << endmsg;
+ debug() << "Size before = " << inputs.size() << ", after = " << outputs.size() << endmsg;
return StatusCode::SUCCESS;
}
diff --git a/JugReco/src/components/CalorimeterIslandCluster.cpp b/JugReco/src/components/CalorimeterIslandCluster.cpp
index a6dc5b679455fcced5132d527583039a42f574f1..48ab9d2cd60423b664578095435a181830e6eb57 100644
--- a/JugReco/src/components/CalorimeterIslandCluster.cpp
+++ b/JugReco/src/components/CalorimeterIslandCluster.cpp
@@ -1,8 +1,8 @@
/*
* Island Clustering Algorithm for Calorimeter Blocks
* 1. group all the adjacent modules
- * 2. split the groups between their local maxima with the energy deposit above
- * <minClusterCenterEdep>
+ * 2. split the groups between their local maxima with the energy deposit above <minClusterCenterEdep>
+ * Output hits collection with clusterID
*
* Author: Chao Peng (ANL), 09/27/2020
* References:
@@ -11,6 +11,7 @@
*/
#include <algorithm>
+#include "fmt/format.h"
#include "Gaudi/Property.h"
#include "GaudiAlg/GaudiAlgorithm.h"
#include "GaudiAlg/GaudiTool.h"
@@ -37,14 +38,18 @@ namespace Jug::Reco {
class CalorimeterIslandCluster : public GaudiAlgorithm {
public:
- Gaudi::Property<bool> m_splitCluster{this, "splitCluster", true};
- Gaudi::Property<double> m_groupRange{this, "groupRange", 1.8};
- Gaudi::Property<std::vector<double>> u_groupRanges{this, "groupRanges", {}};
- Gaudi::Property<double> m_minClusterCenterEdep{this, "minClusterCenterEdep", 50.0 * MeV};
- DataHandle<eic::CalorimeterHitCollection> m_inputHitCollection{"inputHitCollection",
- Gaudi::DataHandle::Reader, this};
- DataHandle<eic::CalorimeterHitCollection> m_splitHitCollection{"outputHitCollection",
- Gaudi::DataHandle::Writer, this};
+ Gaudi::Property<bool> m_splitCluster{this, "splitCluster", true};
+ Gaudi::Property<double> m_dimScaledDist{this, "dimScaledDist", 1.8};
+ Gaudi::Property<std::vector<double>> u_localDistXY{this, "localDistXY", {}};
+ Gaudi::Property<double> m_minClusterHitEdep{this, "minClusterHitEdep", 0.};
+ Gaudi::Property<double> m_minClusterCenterEdep{this, "minClusterCenterEdep", 50.0 * MeV};
+ DataHandle<eic::CalorimeterHitCollection> m_inputHitCollection{"inputHitCollection",
+ Gaudi::DataHandle::Reader, this};
+ DataHandle<eic::CalorimeterHitCollection> m_splitHitCollection{"outputHitCollection",
+ Gaudi::DataHandle::Writer, this};
+
+ // unitless counterparts of the input parameters
+ double minClusterHitEdep, minClusterCenterEdep, localDistXY[2] = {0., 0.};
CalorimeterIslandCluster(const std::string& name, ISvcLocator* svcLoc)
: GaudiAlgorithm(name, svcLoc)
@@ -58,6 +63,21 @@ namespace Jug::Reco {
if (GaudiAlgorithm::initialize().isFailure()) {
return StatusCode::FAILURE;
}
+
+ // unitless conversion, keep consistency with juggler internal units (GeV, mm, ns, rad)
+ minClusterHitEdep = m_minClusterHitEdep.value() / GeV;
+ minClusterCenterEdep = m_minClusterCenterEdep.value() / GeV;
+ // sanity check
+ if (u_localDistXY.size() == 2) {
+ localDistXY[0] = u_localDistXY.value()[0] / mm;
+ localDistXY[1] = u_localDistXY.value()[1] / mm;
+ info() << fmt::format("Clustering using distance [x, y] <= [{:.4f} mm, {:.4f} mm]",
+ localDistXY[0], localDistXY[1]) << endmsg;
+ } else if (u_localDistXY.size() > 0) {
+ warning() << fmt::format("Expect two values (x, y) from localDistXY, got {}. "
+ "Using dimScaledDist = {} instead",
+ u_localDistXY.size(), m_dimScaledDist.value()) << endmsg;
+ }
return StatusCode::SUCCESS;
}
@@ -73,6 +93,11 @@ namespace Jug::Reco {
std::vector<bool> visits(hits.size(), false);
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_x(), hits[i].local_y(),
+ hits[i].x(), hits[i].y(), hits[i].z(),
+ hits[i].layerID(), hits[i].sectorID()) << endmsg;
// already in a group
if (visits[i]) {
continue;
@@ -81,14 +106,16 @@ namespace Jug::Reco {
// create a new group, and group all the neighboring hits
dfs_group(groups.back(), i, hits, visits);
}
- debug() << "we have " << groups.size() << " groups of hits" << endmsg;
size_t clusterID = 0;
for (auto& group : groups) {
- auto maxima = find_maxima(group, !m_splitCluster);
+ if (group.empty()) {
+ continue;
+ }
+ auto maxima = find_maxima(group, !m_splitCluster.value());
split_group(group, maxima, clusterID, clustered_hits);
debug() << "hits in a group: " << group.size() << ", "
- << "local maxima: " << maxima.hits_size() << endmsg;
+ << "local maxima: " << maxima.size() << endmsg;
debug() << "total number of clusters so far: " << clusterID << ", " << endmsg;
}
@@ -112,23 +139,23 @@ namespace Jug::Reco {
{
// in the same sector
if (h1.sectorID() == h2.sectorID()) {
- if (u_groupRanges.size() >= 2) {
- return (std::abs(h1.local_x() - h2.local_x()) <= u_groupRanges[0]) &&
- (std::abs(h1.local_y() - h2.local_y()) <= u_groupRanges[1]);
+ // use absolute value
+ if (u_localDistXY.size() == 2) {
+ return (std::abs(h1.local_x() - h2.local_x()) <= localDistXY[0]) &&
+ (std::abs(h1.local_y() - h2.local_y()) <= localDistXY[1]);
+ // use scaled value (module size as the scales)
} else {
return (std::abs(h1.local_x() - h2.local_x()) <=
- m_groupRange * (h1.dim_x() + h2.dim_x()) / 2.) &&
+ m_dimScaledDist * (h1.dim_x() + h2.dim_x()) / 2.) &&
(std::abs(h1.local_y() - h2.local_y()) <=
- m_groupRange * (h1.dim_y() + h2.dim_y()) / 2.);
+ m_dimScaledDist * (h1.dim_y() + h2.dim_y()) / 2.);
}
- // different sector
+ // different sector, local coordinates do not work, using global coordinates
} else {
double range =
- (u_groupRanges.size() >= 2)
- ? std::sqrt(std::accumulate(u_groupRanges.begin(), u_groupRanges.end(), 0.,
- [](double s, double x) { return s + x * x; }))
- : m_groupRange *
- dist2d((h1.dim_x() + h2.dim_x()) / 2., (h1.dim_y() + h2.dim_y()) / 2.);
+ (u_localDistXY.size() == 2)
+ ? dist2d(localDistXY[0], localDistXY[1])
+ : m_dimScaledDist * dist2d((h1.dim_x() + h2.dim_x()) / 2., (h1.dim_y() + h2.dim_y()) / 2.);
// sector may have rotation (barrel), so z is included
return dist3d(h1.x() - h2.x(), h1.y() - h2.y(), h1.z() - h2.z()) <= range;
}
@@ -138,6 +165,12 @@ namespace Jug::Reco {
void dfs_group(std::vector<eic::CalorimeterHit>& group, int idx,
const eic::CalorimeterHitCollection& hits, std::vector<bool>& visits) const
{
+ // not a qualified hit to particpate clustering, stop here
+ if (hits[idx].energy() < minClusterHitEdep) {
+ visits[idx] = true;
+ return;
+ }
+
eic::CalorimeterHit hit{hits[idx].cellID(), hits[idx].clusterID(),
hits[idx].layerID(), hits[idx].sectorID(),
hits[idx].energy(), hits[idx].time(),
@@ -154,10 +187,10 @@ namespace Jug::Reco {
}
// find local maxima that above a certain threshold
- eic::Cluster find_maxima(const std::vector<eic::CalorimeterHit>& group,
- bool global = false) const
+ std::vector<eic::ConstCalorimeterHit> find_maxima(const std::vector<eic::CalorimeterHit>& group,
+ bool global = false) const
{
- eic::Cluster maxima;
+ std::vector<eic::ConstCalorimeterHit> maxima;
if (group.empty()) {
return maxima;
}
@@ -169,13 +202,15 @@ namespace Jug::Reco {
mpos = i;
}
}
- maxima.addhits(group[mpos]);
+ if (group[mpos].energy() >= minClusterCenterEdep) {
+ maxima.push_back(group[mpos]);
+ }
return maxima;
}
for (auto& hit : group) {
// not a qualified center
- if (hit.energy() < m_minClusterCenterEdep / GeV) {
+ if (hit.energy() < minClusterCenterEdep) {
continue;
}
@@ -191,7 +226,7 @@ namespace Jug::Reco {
}
if (maximum) {
- maxima.addhits(hit);
+ maxima.push_back(hit);
}
}
@@ -212,13 +247,14 @@ namespace Jug::Reco {
// split a group of hits according to the local maxima
// split_hits is used to persistify the data
- void split_group(std::vector<eic::CalorimeterHit>& group, const eic::Cluster& maxima,
+ void split_group(std::vector<eic::CalorimeterHit>& group,
+ const std::vector<eic::ConstCalorimeterHit>& maxima,
size_t& clusterID, eic::CalorimeterHitCollection& clustered_hits) const
{
// special cases
- if (maxima.hits_size() == 0) {
+ if (maxima.size() == 0) {
return;
- } else if (maxima.hits_size() == 1) {
+ } else if (maxima.size() == 1) {
for (auto& hit : group) {
hit.clusterID(clusterID);
clustered_hits.push_back(hit);
@@ -228,15 +264,15 @@ namespace Jug::Reco {
}
// split between maxima
- std::vector<double> weights(maxima.hits_size());
- std::vector<eic::Cluster> splits(maxima.hits_size());
+ std::vector<double> weights(maxima.size());
+ std::vector<eic::Cluster> splits(maxima.size());
size_t n_clus = clusterID + 1;
size_t i = 0;
for (auto it = group.begin(); it != group.end(); ++it, ++i) {
auto hedep = it->energy();
size_t j = 0;
// calculate weights for local maxima
- for (auto cit = maxima.hits_begin(); cit != maxima.hits_end(); ++cit, ++j) {
+ for (auto cit = maxima.begin(); cit != maxima.end(); ++cit, ++j) {
double dist_ref = cit->dim_x();
double energy = cit->energy();
double dist = std::sqrt(std::pow(it->local_x() - cit->local_x(), 2) +
diff --git a/JugReco/src/components/ImagingClusterReco.cpp b/JugReco/src/components/ImagingClusterReco.cpp
index 59bc273240ca66c5d92583abc2c6a62b5b5b1109..afdc4f72f9df196f96a3f176d7603a7b2375425c 100644
--- a/JugReco/src/components/ImagingClusterReco.cpp
+++ b/JugReco/src/components/ImagingClusterReco.cpp
@@ -102,7 +102,7 @@ namespace Jug::Reco {
// debug output
for (const auto& cl : clusters) {
debug() << fmt::format("Cluster {:d}: Edep = {:.3f} MeV, Dir = ({:.3f}, {:.3f}) deg",
- cl.clusterID(), cl.edep() / MeV, cl.cl_theta() / M_PI * 180.,
+ cl.clusterID(), cl.edep()*1000., cl.cl_theta() / M_PI * 180.,
cl.cl_phi() / M_PI * 180.)
<< endmsg;
}
diff --git a/JugReco/src/components/ImagingPixelReco.cpp b/JugReco/src/components/ImagingPixelReco.cpp
index 6b411b37012ff6602876b1150e66e853b3d709a5..9dab6d552e873ea49118a6ce7d1e9be544f599f7 100644
--- a/JugReco/src/components/ImagingPixelReco.cpp
+++ b/JugReco/src/components/ImagingPixelReco.cpp
@@ -44,6 +44,10 @@ namespace Jug::Reco {
Gaudi::Property<double> m_dyRangeADC{this, "dynamicRangeADC", 100 * MeV};
Gaudi::Property<double> m_pedSigmaADC{this, "pedestalSigma", 3.2};
Gaudi::Property<double> m_thresholdADC{this, "thresholdFactor", 3.0};
+
+ // unitless counterparts for the input parameters
+ double dyRangeADC;
+
// hits containers
DataHandle<eic::RawCalorimeterHitCollection> m_inputHitCollection{
"inputHitCollection", Gaudi::DataHandle::Reader, this};
@@ -90,6 +94,9 @@ namespace Jug::Reco {
return StatusCode::FAILURE;
}
+ // unitless conversion
+ dyRangeADC = m_dyRangeADC.value()/GeV;
+
return StatusCode::SUCCESS;
}
@@ -106,8 +113,7 @@ namespace Jug::Reco {
if ((rh.amplitude() - m_pedMeanADC) < m_thresholdADC * m_pedSigmaADC) {
continue;
}
- double edep = (rh.amplitude() - m_pedMeanADC) / (double)m_capADC *
- m_dyRangeADC; // convert ADC -> energy
+ double edep = (rh.amplitude() - m_pedMeanADC) / (double)m_capADC * dyRangeADC; // convert ADC -> energy
double time = rh.timeStamp(); // ns
auto id = rh.cellID();
int lid = (int)id_dec->get(id, layer_idx);
diff --git a/JugReco/src/components/ImagingTopoCluster.cpp b/JugReco/src/components/ImagingTopoCluster.cpp
index afe8f12cc551efe71ff7d469076a2320e87cbd2c..4922d4571b543fe712250c55dbba52b59a607f93 100644
--- a/JugReco/src/components/ImagingTopoCluster.cpp
+++ b/JugReco/src/components/ImagingTopoCluster.cpp
@@ -37,20 +37,22 @@ namespace Jug::Reco {
class ImagingTopoCluster : public GaudiAlgorithm {
public:
// maximum difference in layer numbers that can be considered as neighbours
- Gaudi::Property<int> m_adjLayerDiff{this, "adjLayerDiff", 1};
+ 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_localRanges{this, "localRanges", {1.0 * mm, 1.0 * mm}};
+ 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_adjLayerRanges{
- this, "adjLayerRanges", {0.01 * M_PI, 0.01 * M_PI}};
+ 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_adjSectorDist{this, "adjSectorDist", 1.0 * cm};
+ 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_minEdep{this, "minClusterEdep", 0.5 * MeV};
+ Gaudi::Property<double> m_minClusterEdep{this, "minClusterEdep", 0.5 * MeV};
// minimum number of hits (to save this cluster)
- Gaudi::Property<int> m_minNhits{this, "minClusterNhits", 10};
+ Gaudi::Property<int> m_minClusterNhits{this, "minClusterNhits", 10};
// input hits collection
DataHandle<eic::ImagingPixelCollection> m_inputHitCollection{"inputHitCollection",
Gaudi::DataHandle::Reader, this};
@@ -58,6 +60,10 @@ namespace Jug::Reco {
DataHandle<eic::ImagingPixelCollection> m_outputHitCollection{"outputHitCollection",
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, "");
@@ -70,25 +76,37 @@ namespace Jug::Reco {
return StatusCode::FAILURE;
}
- // check local clustering range
- if (u_localRanges.size() < 2) {
- error() << "Need 2-dimensional ranges for same-layer clustering, only have "
- << u_localRanges.size() << endmsg;
+ // unitless conversion
+ // sanity checks
+ if (u_localDistXY.size() != 2) {
+ error() << "Expected 2 values (x_dist, y_dist) for localDistXY" << endmsg;
return StatusCode::FAILURE;
}
- info() << "Same layer hits group ranges"
- << " (" << u_localRanges[0] / mm << " mm, " << u_localRanges[1] / mm << " mm)"
- << endmsg;
-
- // check adjacent layer clustering range
- if (u_adjLayerRanges.size() < 2) {
- error() << "Need 2-dimensional ranges for adjacent-layer clustering, only have "
- << u_adjLayerRanges.size() << endmsg;
+ if (u_layerDistEtaPhi.size() != 2) {
+ error() << "Expected 2 values (eta_dist, phi_dist) for layerDistEtaPhi" << endmsg;
return StatusCode::FAILURE;
}
- info() << "Same layer hits group ranges"
- << " (" << u_adjLayerRanges[0] / M_PI * 1000. << " mrad, "
- << u_adjLayerRanges[1] / M_PI * 1000. << " mrad)" << 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;
}
@@ -105,7 +123,7 @@ namespace Jug::Reco {
std::vector<std::vector<eic::ImagingPixel>> groups;
for (size_t i = 0; i < hits.size(); ++i) {
// already in a group, or not energetic enough to form a cluster
- if (visits[i] || hits[i].edep() < m_minClusterCenterEdep) {
+ if (visits[i] || hits[i].edep() < minClusterCenterEdep) {
continue;
}
groups.emplace_back();
@@ -116,14 +134,14 @@ namespace Jug::Reco {
size_t clusterID = 0;
for (const auto& group : groups) {
- if ((int)group.size() < m_minNhits) {
+ if (static_cast<int>(group.size()) < m_minClusterNhits.value()) {
continue;
}
double edep = 0.;
for (const auto& hit : group) {
edep += hit.edep();
}
- if (edep < m_minEdep) {
+ if (edep < minClusterEdep) {
continue;
}
for (auto hit : group) {
@@ -148,18 +166,18 @@ namespace Jug::Reco {
// different sectors, simple distance check
if (h1.sectorID() != h2.sectorID()) {
return std::sqrt(pow2(h1.x() - h2.x()) + pow2(h1.y() - h2.y()) + pow2(h1.z() - h2.z())) <=
- m_adjSectorDist;
+ sectorDist;
}
// layer check
int ldiff = std::abs(h1.layerID() - h2.layerID());
// same layer, check local positions
if (!ldiff) {
- 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) {
- return (std::abs(h1.eta() - h2.eta()) <= u_adjLayerRanges[0]) &&
- (std::abs(h1.phi() - h2.phi()) <= u_adjLayerRanges[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.phi() - h2.phi()) <= layerDistEtaPhi[1]);
}
// not in adjacent layers
@@ -170,6 +188,12 @@ namespace Jug::Reco {
void dfs_group(std::vector<eic::ImagingPixel>& group, int idx,
const eic::ImagingPixelCollection& hits, std::vector<bool>& visits) const
{
+ // not a qualified hit to participate in clustering, stop here
+ if (hits[idx].edep() < minClusterHitEdep) {
+ visits[idx] = true;
+ return;
+ }
+
eic::ImagingPixel hit{hits[idx].clusterID(), hits[idx].layerID(), hits[idx].sectorID(),
hits[idx].hitID(), hits[idx].edep(), hits[idx].time(),
hits[idx].eta(), hits[idx].local(), hits[idx].position(),