diff --git a/.clang-format b/.clang-format
index 6fd6f902b661fe1114f04c8d482b7059b3c2a980..5737aca1ec7df306124ffc6d805befe3ff6ade12 100644
--- a/.clang-format
+++ b/.clang-format
@@ -1,112 +1,15 @@
---
-Language: Cpp
-BasedOnStyle: Chromium
-AccessModifierOffset: -2
-AlignAfterOpenBracket: Align
-AlignConsecutiveAssignments: true
-AlignConsecutiveDeclarations: true
-AlignEscapedNewlines: Right
-AlignOperands: true
-AlignTrailingComments: true
-AllowAllParametersOfDeclarationOnNextLine: true
-AllowShortBlocksOnASingleLine: false
-AllowShortCaseLabelsOnASingleLine: false
-AllowShortFunctionsOnASingleLine: All
-AllowShortIfStatementsOnASingleLine: false
-AllowShortLoopsOnASingleLine: false
-AlwaysBreakAfterDefinitionReturnType: None
-AlwaysBreakAfterReturnType: None
-AlwaysBreakBeforeMultilineStrings: false
-AlwaysBreakTemplateDeclarations: true
-BinPackArguments: true
-BinPackParameters: true
-BraceWrapping:
- AfterClass: false
- AfterControlStatement: false
- AfterEnum: false
- AfterFunction: true
- AfterNamespace: false
- AfterObjCDeclaration: false
- AfterStruct: false
- AfterUnion: false
- AfterExternBlock: false
- BeforeCatch: false
- BeforeElse: false
- IndentBraces: false
- SplitEmptyFunction: true
- SplitEmptyRecord: true
- SplitEmptyNamespace: true
-BreakBeforeBinaryOperators: None
-BreakBeforeBraces: Custom
-BreakBeforeInheritanceComma: false
-BreakBeforeTernaryOperators: true
-BreakConstructorInitializersBeforeComma: false
-BreakConstructorInitializers: BeforeColon
-BreakAfterJavaFieldAnnotations: false
-BreakStringLiterals: true
-ColumnLimit: 120
-CommentPragmas: '^ IWYU pragma:'
-CompactNamespaces: false
-ConstructorInitializerAllOnOneLineOrOnePerLine: false
-ConstructorInitializerIndentWidth: 4
-ContinuationIndentWidth: 4
+BasedOnStyle: LLVM
+BreakConstructorInitializersBeforeComma: true
+ConstructorInitializerAllOnOneLineOrOnePerLine: true
Cpp11BracedListStyle: true
-DerivePointerAlignment: false
-DisableFormat: false
-ExperimentalAutoDetectBinPacking: false
-FixNamespaceComments: true
-ForEachMacros:
- - foreach
- - Q_FOREACH
- - BOOST_FOREACH
+Standard: Cpp11
+#SpaceBeforeParens: ControlStatements
+SpaceAfterControlStatementKeyword: true
+PointerBindsToType: true
IncludeBlocks: Preserve
-IncludeCategories:
- - Regex: '^"(llvm|llvm-c|clang|clang-c)/'
- Priority: 2
- - Regex: '^(<|"(gtest|gmock|isl|json)/)'
- Priority: 3
- - Regex: '.*'
- Priority: 1
-IncludeIsMainRegex: '(Test)?$'
-IndentCaseLabels: false
-IndentPPDirectives: None
-IndentWidth: 2
-IndentWrappedFunctionNames: false
-JavaScriptQuotes: Leave
-JavaScriptWrapImports: true
-KeepEmptyLinesAtTheStartOfBlocks: true
-MacroBlockBegin: ''
-MacroBlockEnd: ''
-MaxEmptyLinesToKeep: 1
-NamespaceIndentation: All
-ObjCBlockIndentWidth: 2
-ObjCSpaceAfterProperty: false
-ObjCSpaceBeforeProtocolList: true
-PenaltyBreakAssignment: 2
-PenaltyBreakBeforeFirstCallParameter: 19
-PenaltyBreakComment: 300
-PenaltyBreakFirstLessLess: 120
-PenaltyBreakString: 1000
-PenaltyExcessCharacter: 1000000
-PenaltyReturnTypeOnItsOwnLine: 60
-PointerAlignment: Left
-ReflowComments: true
-SortIncludes: true
-SortUsingDeclarations: true
-SpaceAfterCStyleCast: false
-SpaceAfterTemplateKeyword: true
-SpaceBeforeAssignmentOperators: true
-SpaceBeforeParens: ControlStatements
-#SpaceBeforeRangeBasedForLoopColon: true
-SpaceInEmptyParentheses: false
-SpacesBeforeTrailingComments: 1
-SpacesInAngles: false
-SpacesInContainerLiterals: true
-SpacesInCStyleCastParentheses: false
-SpacesInParentheses: false
-SpacesInSquareBrackets: false
-Standard: Cpp11
-TabWidth: 4
UseTab: Never
+ColumnLimit: 120
+NamespaceIndentation: Inner
+AlignConsecutiveAssignments: true
...
-
diff --git a/JugDigi/src/components/CalorimeterHitDigi.cpp b/JugDigi/src/components/CalorimeterHitDigi.cpp
index 2dc231974219e6e45f6a3c945e0dadcabd8b3a0f..ec07b3892e21b71752ee757a6b79b9fb4ac128c4 100644
--- a/JugDigi/src/components/CalorimeterHitDigi.cpp
+++ b/JugDigi/src/components/CalorimeterHitDigi.cpp
@@ -88,6 +88,7 @@ namespace Jug::Digi {
const auto simhits = m_inputHitCollection.get();
// Create output collections
auto rawhits = m_outputHitCollection.createAndPut();
+ int nhits = 0;
for (const auto& ahit : *simhits) {
// Note: juggler internal unit of energy is GeV
double eResRel = std::sqrt(std::pow(m_normDist() * eRes[0] / sqrt(ahit.energyDeposit()), 2) +
@@ -95,8 +96,11 @@ namespace Jug::Digi {
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. + eResRel) / dyRangeADC * m_capADC);
- eic::RawCalorimeterHit rawhit((long long)ahit.cellID(), (adc > m_capADC.value() ? m_capADC.value() : adc),
- (double)ahit.truth().time + m_normDist() * tRes);
+ eic::RawCalorimeterHit rawhit(
+ (long long)ahit.cellID(),
+ (adc > m_capADC.value() ? m_capADC.value() : adc),
+ static_cast<int64_t>(1e6*(double)ahit.truth().time + m_normDist() * tRes), // @FIXME: this shouldn't be hardcoded, but should still be stored as an integer type
+ nhits++);
rawhits->push_back(rawhit);
}
return StatusCode::SUCCESS;
diff --git a/JugDigi/src/components/CrystalEndcapsDigi.cpp b/JugDigi/src/components/CrystalEndcapsDigi.cpp
index b6bff48bd83a682429dcd590764af70a30058eef..9f3abfe09fbf81e21072d07ac7c75d463094a7c7 100644
--- a/JugDigi/src/components/CrystalEndcapsDigi.cpp
+++ b/JugDigi/src/components/CrystalEndcapsDigi.cpp
@@ -57,13 +57,14 @@ namespace Jug {
// Create output collections
auto rawhits = m_outputHitCollection.createAndPut();
eic::RawCalorimeterHitCollection* rawHitCollection = new eic::RawCalorimeterHitCollection();
+ int nhits = 0;
for (const auto& ahit : *simhits) {
double res = m_gaussDist()/sqrt(ahit.energyDeposit());
eic::RawCalorimeterHit rawhit(
(long long) ahit.cellID(),
std::llround(ahit.energyDeposit() * (1. + res)*1.0e6), // convert to keV integer
- (double) ahit.truth().time);
- rawhits->push_back(rawhit);
+ (double) ahit.truth().time, nhits++);
+ rawhits->push_back(rawhit);
}
return StatusCode::SUCCESS;
}
diff --git a/JugDigi/src/components/EMCalorimeterDigi.cpp b/JugDigi/src/components/EMCalorimeterDigi.cpp
index 5a678d5c1e89f095d67afbb2cce870db21c268a3..0be573682c6c62b8c056ef5db9a2ff2867e94154 100644
--- a/JugDigi/src/components/EMCalorimeterDigi.cpp
+++ b/JugDigi/src/components/EMCalorimeterDigi.cpp
@@ -57,13 +57,14 @@ namespace Jug {
// Create output collections
auto rawhits = m_outputHitCollection.createAndPut();
eic::RawCalorimeterHitCollection* rawHitCollection = new eic::RawCalorimeterHitCollection();
+ int nhits = 0;
for (const auto& ahit : *simhits) {
// std::cout << ahit << "\n";
double sqrtE = std::sqrt(ahit.energyDeposit()) ;
double aterm = m_gaussDist()*sqrtE;
eic::RawCalorimeterHit rawhit((long long)ahit.cellID(),
std::llround((ahit.energyDeposit() + aterm) * 1e6),
- ahit.truth().time * 1e6);
+ ahit.truth().time * 1e6, nhits++);
rawhits->push_back(rawhit);
}
return StatusCode::SUCCESS;
diff --git a/JugDigi/src/components/EcalTungstenSamplingDigi.cpp b/JugDigi/src/components/EcalTungstenSamplingDigi.cpp
index 452634075465aeba6f8d919bf49b6ba8f80d549a..88d67871ed32ec1bf618ff51b9c33ace92b86f03 100644
--- a/JugDigi/src/components/EcalTungstenSamplingDigi.cpp
+++ b/JugDigi/src/components/EcalTungstenSamplingDigi.cpp
@@ -77,6 +77,7 @@ namespace Jug {
// Create output collections
auto rawhits = m_outputHitCollection.createAndPut();
eic::RawCalorimeterHitCollection* rawHitCollection = new eic::RawCalorimeterHitCollection();
+ int nhits = 0;
for (const auto& ahit : *simhits) {
double resval = std::pow(m_normDist()*res[0] / sqrt(ahit.energyDeposit()*m_eUnit/GeV), 2)
+ std::pow(m_normDist()*res[1], 2)
@@ -87,7 +88,8 @@ namespace Jug {
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);
+ (double)ahit.truth().time*m_tUnit/ns + m_normDist()*m_tRes/ns,
+ nhits++);
rawhits->push_back(rawhit);
}
return StatusCode::SUCCESS;
diff --git a/JugDigi/src/components/ExampleCaloDigi.cpp b/JugDigi/src/components/ExampleCaloDigi.cpp
index 7e248182e23084928ea02535504d524e855512c6..4152034a24e381f12f690b5f5e2974e5777f05c0 100644
--- a/JugDigi/src/components/ExampleCaloDigi.cpp
+++ b/JugDigi/src/components/ExampleCaloDigi.cpp
@@ -110,9 +110,10 @@ namespace Jug {
// Create output collections
auto rawhits = m_outputHitCollection.createAndPut();
eic::RawCalorimeterHitCollection* rawHitCollection = new eic::RawCalorimeterHitCollection();
+ int nhits = 0;
for(const auto& ahit : *simhits) {
//std::cout << ahit << "\n";
- eic::RawCalorimeterHit rawhit((long long)ahit.cellID(), std::llround(ahit.energyDeposit() * 100), 0);
+ eic::RawCalorimeterHit rawhit((long long)ahit.cellID(), std::llround(ahit.energyDeposit() * 100), 0, nhits++);
rawhits->push_back(rawhit);
}
return StatusCode::SUCCESS;
diff --git a/JugDigi/src/components/HadronicCalDigi.cpp b/JugDigi/src/components/HadronicCalDigi.cpp
index b0572782a809cc65181dfd63ef01eaae7253f130..a11a2ac95ef878d95450a6a9cf77125c462a6bc5 100644
--- a/JugDigi/src/components/HadronicCalDigi.cpp
+++ b/JugDigi/src/components/HadronicCalDigi.cpp
@@ -71,6 +71,7 @@ namespace Jug {
// Create output collections
auto rawhits = m_outputHitCollection.createAndPut();
eic::RawCalorimeterHitCollection* rawHitCollection = new eic::RawCalorimeterHitCollection();
+ int nhits = 0;
for (const auto& ahit : *simhits) {
// std::cout << ahit << "\n";
double sqrtE = std::sqrt(ahit.energyDeposit()) ;
@@ -78,7 +79,7 @@ namespace Jug {
double bterm = ahit.energyDeposit()*m_gaussDist_b();
// here 1000 is arbitrary scale factor
eic::RawCalorimeterHit rawhit((long long)ahit.cellID(),
- std::llround(ahit.energyDeposit() +aterm + bterm * 1000), 0);
+ std::llround(ahit.energyDeposit() +aterm + bterm * 1000), 0, nhits++);
rawhits->push_back(rawhit);
}
return StatusCode::SUCCESS;
diff --git a/JugReco/src/components/:wq b/JugReco/src/components/:wq
new file mode 100644
index 0000000000000000000000000000000000000000..7c16d5534de6a009709619bb8cf7790561a5fc4f
--- /dev/null
+++ b/JugReco/src/components/:wq
@@ -0,0 +1,229 @@
+/*
+ * Reconstruct the cluster with Center of Gravity method
+ * Logarithmic weighting is used for mimicing energy deposit in transverse direction
+ *
+ * Author: Chao Peng (ANL), 09/27/2020
+ */
+#include <algorithm>
+#include <cstring>
+#include <functional>
+#include <map>
+
+#include "boost/algorithm/string/join.hpp"
+#include "fmt/format.h"
+#include <boost/range/adaptor/map.hpp>
+
+#include "Gaudi/Property.h"
+#include "GaudiAlg/GaudiAlgorithm.h"
+#include "GaudiAlg/GaudiTool.h"
+#include "GaudiAlg/Transformer.h"
+#include "GaudiKernel/PhysicalConstants.h"
+#include "GaudiKernel/RndmGenerators.h"
+#include "GaudiKernel/ToolHandle.h"
+
+#include "DDRec/CellIDPositionConverter.h"
+#include "DDRec/Surface.h"
+#include "DDRec/SurfaceManager.h"
+
+// FCCSW
+#include "JugBase/DataHandle.h"
+#include "JugBase/IGeoSvc.h"
+
+// Event Model related classes
+#include "eicd/CalorimeterHitCollection.h"
+#include "eicd/ProtoClusterCollection.h"
+#include "eicd/ClusterCollection.h"
+#include "eicd/Cluster2DInfoCollection.h"
+#include "eicd/VectorPolar.h"
+
+using namespace Gaudi::Units;
+
+namespace Jug::Reco {
+
+ // weighting functions (with place holders for hit energy, total energy, one parameter and module
+ // type enum
+ static double constWeight(double /*E*/, double /*tE*/, double /*p*/, int /*type*/) { return 1.0; }
+ static double linearWeight(double E, double /*tE*/, double /*p*/, int /*type*/) { return E; }
+ static double logWeight(double E, double tE, double base, int /*type*/)
+ {
+ return std::max(0., base + std::log(E / tE));
+ }
+
+ static const std::map<std::string, std::function<double(double, double, double, int)>>
+ weightMethods {
+ {"none", constWeight},
+ {"linear", linearWeight},
+ {"log", logWeight},
+ };
+
+ /** Clustering with center of gravity method.
+ *
+ * Reconstruct the cluster with Center of Gravity method
+ * Logarithmic weighting is used for mimicing energy deposit in transverse direction
+ *
+ * \ingroup reco
+ */
+ class ClusterRecoCoG : public GaudiAlgorithm {
+ public:
+ Gaudi::Property<double> m_sampFrac{this, "samplingFraction", 1.0};
+ Gaudi::Property<double> m_logWeightBase{this, "logWeightBase", 3.6};
+ Gaudi::Property<double> m_depthCorrection{this, "depthCorrection", 0.0};
+ Gaudi::Property<std::string> m_energyWeight{this, "energyWeight", "log"};
+ Gaudi::Property<std::string> m_moduleDimZName{this, "moduleDimZName", ""};
+ DataHandle<eic::CalorimeterHitCollection> m_inputHits{"inputHitCollection",
+ Gaudi::DataHandle::Reader, this};
+ DataHandle<eic::ProtoClusterCollection> m_inputProto{"inputProtoClusterCollection",
+ Gaudi::DataHandle::Reader, this};
+ DataHandle<eic::ClusterCollection> m_outputClusters{"outputClusterCollection",
+ Gaudi::DataHandle::Writer, this};
+ DataHandle<eic::Cluster2DInfoCollection> m_outputInfo{"outputInfoCollection",
+ Gaudi::DataHandle::Writer, this};
+ // Pointer to the geometry service
+ SmartIF<IGeoSvc> m_geoSvc;
+ double m_depthCorr;
+ std::function<double(double, double, double, int)> weightFunc;
+
+ ClusterRecoCoG(const std::string& name, ISvcLocator* svcLoc) : GaudiAlgorithm(name, svcLoc)
+ {
+ declareProperty("inputHitCollection", m_inputHits, "");
+ declareProperty("inputProtoClusterCollection", m_inputProto, "");
+ declareProperty("outputClusterCollection", m_outputClusters, "");
+ declareProperty("outputInfoCollection", m_outputInfo, "");
+ }
+
+ StatusCode initialize() override
+ {
+ if (GaudiAlgorithm::initialize().isFailure()) {
+ return StatusCode::FAILURE;
+ }
+ m_geoSvc = service("GeoSvc");
+ if (!m_geoSvc) {
+ error() << "Unable to locate Geometry Service. "
+ << "Make sure you have GeoSvc and SimSvc in the right order in the configuration."
+ << endmsg;
+ return StatusCode::FAILURE;
+ }
+ // update depth correction if a name is provided
+ if (!m_moduleDimZName.value().empty()) {
+ m_depthCorrection = m_geoSvc->detector()->constantAsDouble(m_moduleDimZName);
+ }
+
+ // select weighting method
+ std::string ew = m_energyWeight.value();
+ // make it case-insensitive
+ std::transform(ew.begin(), ew.end(), ew.begin(), [](char s) { return std::tolower(s); });
+ auto it = weightMethods.find(ew);
+ if (it == weightMethods.end()) {
+ error() << fmt::format(
+ "Cannot find energy weighting method {}, choose one from [{}]",
+ m_energyWeight,
+ boost::algorithm::join(weightMethods | boost::adaptors::map_keys, ", "))
+ << endmsg;
+ return StatusCode::FAILURE;
+ }
+ weightFunc = it->second;
+ // info() << "z_length " << depth << endmsg;
+ return StatusCode::SUCCESS;
+ }
+
+ StatusCode execute() override
+ {
+ // input collections
+ const auto& hits = *m_inputHits.get();
+ const auto& proto = *m_inputProto.get();
+ auto& clusters = *m_outputClusters.createAndPut();
+ auto& info = *m_outputInfo.createAndPut();
+
+ // Create a map of clusterID --> associated hits by looping over our clustered hits
+ std::map<int, std::vector<std::pair<eic::ConstProtoCluster,
+ eic::ConstCalorimeterHit>>> cluster_map;
+ for (const auto& pc : proto) {
+ const size_t clusterID = pc.clusterID();
+ if (!cluster_map.count(clusterID)) {
+ cluster_map[clusterID] = {};
+ }
+ size_t idx;
+ for (idx = 0; idx < hits.size(); ++idx) {
+ if (hits[idx].ID() == pc.hitID()) {
+ break;
+ }
+ }
+ if (idx >= hits.size()) {
+ continue;
+ }
+ cluster_map[clusterID].push_back({pc, hits[idx]});
+ }
+
+ for (const auto& [idx, hit_info] : cluster_map) {
+ auto cl = reconstruct(hit_info, idx);
+
+ debug() << cl.nhits() << " hits: " << cl.energy() / GeV << " GeV, (" << cl.position().x / mm
+ << ", " << cl.position().y / mm << ", " << cl.position().z / mm << ")" << endmsg;
+ clusters.push_back(cl);
+ info.push_back({cl.ID(), cl.position(), cl.position().eta()});
+ }
+
+ return StatusCode::SUCCESS;
+ }
+
+ private:
+ template <typename T1>
+ eic::VectorPolar cart_to_polar(const T1& cart)
+ {
+ auto r = std::sqrt(cart.x * cart.x + cart.y * cart.y + cart.z * cart.z);
+ return eic::VectorPolar{r, std::acos(cart.z / r), std::atan2(cart.y, cart.x)};
+ }
+
+ eic::Cluster reconstruct(const std::vector<std::pair<eic::ConstProtoCluster,
+ eic::ConstCalorimeterHit>>& hit_info,
+ const int idx) const
+ {
+ eic::Cluster cl;
+ cl.ID(idx);
+ cl.nhits(hit_info.size());
+
+ // no hits
+ debug() << "hit size = " << hit_info.size() << endmsg;
+ if (hit_info.empty()) {
+ return cl;
+ }
+
+ // calculate total energy, find the cell with the maximum energy deposit
+ float totalE = 0.;
+ float maxE = 0.;
+ auto time = hit_info[0].second.time();
+ for (const auto& [proto, hit] : hit_info) {
+ debug() << "hit energy = " << hit.energy() << " hit weight: " << proto.weight() << endmsg;
+ auto energy = hit.energy() * proto.weight();
+ totalE += energy;
+ if (energy > maxE) {
+ maxE = energy;
+ time = hit.time();
+ }
+ }
+ cl.energy(totalE/m_sampFrac);
+ cl.energyError(0.);
+ cl.time(time);
+
+ // center of gravity with logarithmic weighting
+ float tw = 0.;
+ eic::VectorXYZ v;
+ for (const auto& [proto, hit] : hit_info) {
+ float w = weightFunc(hit.energy()*proto.weight(), totalE, m_logWeightBase.value(), 0);
+ tw += w;
+ v = v.add(hit.position().scale(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.positionError({}); // @TODO: Covariance matrix
+
+ return cl;
+ }
+ };
+
+ DECLARE_COMPONENT(ClusterRecoCoG)
+
+} // namespace Jug::Reco
+
diff --git a/JugReco/src/components/CalorimeterHitReco.cpp b/JugReco/src/components/CalorimeterHitReco.cpp
index 5efbc0b21bf618014ac2f5cc95cb4493c5c72203..5dc5cfc2323dbb46083e5611a1ff5f44a3a10642 100644
--- a/JugReco/src/components/CalorimeterHitReco.cpp
+++ b/JugReco/src/components/CalorimeterHitReco.cpp
@@ -51,6 +51,9 @@ namespace Jug::Reco {
// zero suppression values
Gaudi::Property<double> m_thresholdADC{this, "thresholdFactor", 3.0};
+ // Calibration!
+ Gaudi::Property<double> m_sampFrac{this, "samplingFraction", 1.0};
+
// unitless counterparts of the input parameters
double dyRangeADC;
@@ -165,35 +168,34 @@ namespace Jug::Reco {
}
// convert ADC -> energy
- float energy = (rh.amplitude() - m_pedMeanADC) / static_cast<float>(m_capADC.value()) * dyRangeADC;
+ float energy = (rh.amplitude() - m_pedMeanADC) / static_cast<float>(m_capADC.value()) * dyRangeADC / m_sampFrac;
- float time = rh.timeStamp(); // ns
- auto id = rh.cellID();
+ float time = rh.time()*1.e-6; // ns @FIXME: this should not be hardcoded
+ auto cellID = rh.cellID();
int lid = ((id_dec != nullptr) && m_layerField.value().size())
- ? static_cast<int>(id_dec->get(id, layer_idx))
+ ? static_cast<int>(id_dec->get(cellID, layer_idx))
: -1;
int sid = ((id_dec != nullptr) && m_sectorField.value().size())
- ? static_cast<int>(id_dec->get(id, sector_idx))
+ ? static_cast<int>(id_dec->get(cellID, sector_idx))
: -1;
// global positions
- auto gpos = converter->position(id);
+ auto gpos = converter->position(cellID);
// local positions
if (m_localDetElement.value().empty()) {
auto volman = m_geoSvc->detector()->volumeManager();
- local = volman.lookupDetElement(id & local_mask);
+ local = volman.lookupDetElement(cellID & local_mask);
}
auto pos = local.nominal().worldToLocal(dd4hep::Position(gpos.x(), gpos.y(), gpos.z()));
- // auto pos = m_geoSvc->cellIDPositionConverter()->findContext(id)->volumePlacement().position();
// cell dimension
std::vector<double> cdim;
// get segmentation dimensions
if (converter->findReadout(local).segmentation().type() != "NoSegmentation") {
- cdim = converter->cellDimensions(id);
+ cdim = converter->cellDimensions(cellID);
// get volume dimensions (multiply by two to get fullsize)
} else {
// cdim = converter->findContext(id)->volumePlacement().volume().solid().dimensions();
// Using bounding box instead of actual solid so the dimensions are always in dim_x, dim_y, dim_z
- cdim = converter->findContext(id)->volumePlacement().volume().boundingBox().dimensions();
+ cdim = converter->findContext(cellID)->volumePlacement().volume().boundingBox().dimensions();
std::transform(cdim.begin(), cdim.end(), cdim.begin(),
std::bind(std::multiplies<double>(), std::placeholders::_1, 2));
}
@@ -206,16 +208,17 @@ namespace Jug::Reco {
// m_geoSvc->cellIDPositionConverter()->findContext(id)->volumePlacement().volIDs().str()
// << endmsg;
hits.push_back({
- id,
- -1,
+ rh.cellID(),
+ rh.ID(),
lid,
- sid, // cell id, cluster id, layer id, sector id
+ sid, // cell id, cluster id, layer id, sector id
+ 0, // @TODO: hit type
energy,
+ 0, // @TODO: energy error
time, // energy, time
{gpos.x() / m_lUnit, gpos.y() / m_lUnit, gpos.z() / m_lUnit},
{pos.x() / m_lUnit, pos.y() / m_lUnit, pos.z() / m_lUnit},
- {dim[0], dim[1], dim[2]},
- 0 // @TODO: hit type
+ {dim[0], dim[1], dim[2]}
});
}
diff --git a/JugReco/src/components/CalorimeterHitsEtaPhiProjector.cpp b/JugReco/src/components/CalorimeterHitsEtaPhiProjector.cpp
index 71e86e4f80a24a4d732561664618444f69e0e551..5faaa418266ce4fc5cfed3dac4fa55355a26a638 100644
--- a/JugReco/src/components/CalorimeterHitsEtaPhiProjector.cpp
+++ b/JugReco/src/components/CalorimeterHitsEtaPhiProjector.cpp
@@ -107,8 +107,8 @@ namespace Jug::Reco {
const auto ref = hits.front();
eic::CalorimeterHit hit;
hit.cellID(ref.cellID());
- hit.sectorID(ref.sectorID());
- hit.layerID(ref.layerID());
+ hit.sector(ref.sector());
+ hit.layer(ref.layer());
// TODO, we can do timing cut to reject noises
hit.time(ref.time());
double r = ref.position().mag();
diff --git a/JugReco/src/components/CalorimeterHitsMerger.cpp b/JugReco/src/components/CalorimeterHitsMerger.cpp
index e2e95a4f48e8c00f691252c1b1f1f6f399c10190..5d18447ce02a4ec1be90e5c1f64c9a793a0d5dbc 100644
--- a/JugReco/src/components/CalorimeterHitsMerger.cpp
+++ b/JugReco/src/components/CalorimeterHitsMerger.cpp
@@ -132,6 +132,7 @@ namespace Jug::Reco {
auto poscon = m_geoSvc->cellIDPositionConverter();
auto volman = m_geoSvc->detector()->volumeManager();
+ int nresults = 0;
for (auto &[id, hits] : merge_map) {
// reference fields id
int64_t ref_id = id | ref_mask;
@@ -151,16 +152,17 @@ namespace Jug::Reco {
}
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()});
+ ref_id,
+ nresults++,
+ href.layer(),
+ href.sector(),
+ href.type(),
+ energy,
+ 0, //@TODO: energy uncertainty
+ 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()});
}
debug() << "Size before = " << inputs.size() << ", after = " << outputs.size() << endmsg;
diff --git a/JugReco/src/components/CalorimeterIslandCluster.cpp b/JugReco/src/components/CalorimeterIslandCluster.cpp
index b1930c49f9e3978e92780681aaee7fd2b967a614..7c05213b6eb3fb1a7d99a21139bd78115ef1715f 100644
--- a/JugReco/src/components/CalorimeterIslandCluster.cpp
+++ b/JugReco/src/components/CalorimeterIslandCluster.cpp
@@ -34,6 +34,7 @@
// Event Model related classes
#include "eicd/CalorimeterHitCollection.h"
#include "eicd/ClusterCollection.h"
+#include "eicd/ProtoClusterCollection.h"
#include "eicd/VectorXY.h"
#include "eicd/VectorXYZ.h"
@@ -93,7 +94,7 @@ namespace Jug::Reco {
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",
+ DataHandle<eic::ProtoClusterCollection> m_outputProtoCollection{"outputProtoClusterCollection",
Gaudi::DataHandle::Writer, this};
// neighbour checking distances
@@ -104,7 +105,7 @@ namespace Jug::Reco {
Gaudi::Property<std::vector<double>> u_globalDistRPhi{this, "globalDistRPhi", {}};
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
+ // neighbor checking function
std::function<eic::VectorXY(eic::ConstCalorimeterHit, eic::ConstCalorimeterHit)> hitsDist;
// unitless counterparts of the input parameters
@@ -115,7 +116,7 @@ namespace Jug::Reco {
: GaudiAlgorithm(name, svcLoc)
{
declareProperty("inputHitCollection", m_inputHitCollection, "");
- declareProperty("outputHitCollection", m_splitHitCollection, "");
+ declareProperty("outputProtoClusterCollection", m_outputProtoCollection, "");
}
StatusCode initialize() override
@@ -175,10 +176,10 @@ namespace Jug::Reco {
// input collections
const auto& hits = *m_inputHitCollection.get();
// Create output collections
- auto& clustered_hits = *m_splitHitCollection.createAndPut();
+ auto& proto = *m_outputProtoCollection.createAndPut();
// group neighboring hits
- std::vector<std::vector<eic::CalorimeterHit>> groups;
+ std::vector<std::vector<eic::ConstCalorimeterHit>> groups;
std::vector<bool> visits(hits.size(), false);
for (size_t i = 0; i < hits.size(); ++i) {
@@ -186,7 +187,7 @@ namespace Jug::Reco {
"global=({:.4f}, {:.4f}, {:.4f}) mm, layer = {:d}, sector = {:d}.",
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;
+ hits[i].layer(), hits[i].sector()) << endmsg;
// already in a group
if (visits[i]) {
continue;
@@ -202,7 +203,7 @@ namespace Jug::Reco {
continue;
}
auto maxima = find_maxima(group, !m_splitCluster.value());
- split_group(group, maxima, clusterID, clustered_hits);
+ split_group(group, maxima, clusterID, proto);
debug() << "hits in a group: " << group.size() << ", "
<< "local maxima: " << maxima.size() << endmsg;
debug() << "total number of clusters so far: " << clusterID << ", " << endmsg;
@@ -216,7 +217,7 @@ namespace Jug::Reco {
inline bool is_neighbour(const eic::ConstCalorimeterHit& h1, const eic::ConstCalorimeterHit& h2) const
{
// in the same sector
- if (h1.sectorID() == h2.sectorID()) {
+ if (h1.sector() == h2.sector()) {
auto dist = hitsDist(h1, h2);
return (dist.x <= neighbourDist[0]) && (dist.y <= neighbourDist[1]);
// different sector, local coordinates do not work, using global coordinates
@@ -227,7 +228,7 @@ namespace Jug::Reco {
}
// grouping function with Depth-First Search
- void dfs_group(std::vector<eic::CalorimeterHit>& group, int idx,
+ void dfs_group(std::vector<eic::ConstCalorimeterHit>& group, int idx,
const eic::CalorimeterHitCollection& hits, std::vector<bool>& visits) const
{
// not a qualified hit to particpate clustering, stop here
@@ -236,15 +237,10 @@ namespace Jug::Reco {
return;
}
- eic::CalorimeterHit hit{hits[idx].cellID(), hits[idx].clusterID(),
- hits[idx].layerID(), hits[idx].sectorID(),
- hits[idx].energy(), hits[idx].time(),
- hits[idx].position(), hits[idx].local(),
- hits[idx].dimension(), 1};
- group.push_back(hit);
+ group.push_back(hits[idx]);
visits[idx] = true;
for (size_t i = 0; i < hits.size(); ++i) {
- if (visits[i] || !is_neighbour(hit, hits[i])) {
+ if (visits[i] || !is_neighbour(hits[idx], hits[i])) {
continue;
}
dfs_group(group, i, hits, visits);
@@ -252,7 +248,7 @@ namespace Jug::Reco {
}
// find local maxima that above a certain threshold
- std::vector<eic::ConstCalorimeterHit> find_maxima(const std::vector<eic::CalorimeterHit>& group,
+ std::vector<eic::ConstCalorimeterHit> find_maxima(const std::vector<eic::ConstCalorimeterHit>& group,
bool global = false) const
{
std::vector<eic::ConstCalorimeterHit> maxima;
@@ -312,17 +308,16 @@ 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,
+ void split_group(std::vector<eic::ConstCalorimeterHit>& group,
const std::vector<eic::ConstCalorimeterHit>& maxima,
- size_t& clusterID, eic::CalorimeterHitCollection& clustered_hits) const
+ size_t& clusterID, eic::ProtoClusterCollection& proto) const
{
// special cases
if (maxima.size() == 0) {
return;
} else if (maxima.size() == 1) {
for (auto& hit : group) {
- hit.clusterID(clusterID);
- clustered_hits.push_back(hit);
+ proto.create(hit.ID(), clusterID, 1.);
}
clusterID += 1;
return;
@@ -331,11 +326,9 @@ namespace Jug::Reco {
// split between maxima
// TODO, here we can implement iterations with profile, or even ML for better splits
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.begin(); cit != maxima.end(); ++cit, ++j) {
@@ -357,19 +350,15 @@ namespace Jug::Reco {
vec_normalize(weights);
// split energy between local maxima
- for (size_t k = 0; k < weights.size(); ++k) {
+ for (size_t k = 0; k < maxima.size(); ++k) {
double weight = weights[k];
if (weight <= 1e-6) {
continue;
}
- auto hit = it->clone();
- hit.energy(hedep * weight);
- hit.type(1);
- hit.clusterID(n_clus + k);
- clustered_hits.push_back(hit);
+ proto.create(it->ID(), n_clus + k, weight);
}
}
- clusterID += splits.size();
+ clusterID += maxima.size();
return;
}
};
diff --git a/JugReco/src/components/ClusterRecoCoG.cpp b/JugReco/src/components/ClusterRecoCoG.cpp
index b3764a92977d308ff1c7ba2563a5b6f6593796c9..7c16d5534de6a009709619bb8cf7790561a5fc4f 100644
--- a/JugReco/src/components/ClusterRecoCoG.cpp
+++ b/JugReco/src/components/ClusterRecoCoG.cpp
@@ -31,7 +31,10 @@
// Event Model related classes
#include "eicd/CalorimeterHitCollection.h"
+#include "eicd/ProtoClusterCollection.h"
#include "eicd/ClusterCollection.h"
+#include "eicd/Cluster2DInfoCollection.h"
+#include "eicd/VectorPolar.h"
using namespace Gaudi::Units;
@@ -69,8 +72,12 @@ namespace Jug::Reco {
Gaudi::Property<std::string> m_moduleDimZName{this, "moduleDimZName", ""};
DataHandle<eic::CalorimeterHitCollection> m_inputHits{"inputHitCollection",
Gaudi::DataHandle::Reader, this};
+ DataHandle<eic::ProtoClusterCollection> m_inputProto{"inputProtoClusterCollection",
+ Gaudi::DataHandle::Reader, this};
DataHandle<eic::ClusterCollection> m_outputClusters{"outputClusterCollection",
- Gaudi::DataHandle::Writer, this};
+ Gaudi::DataHandle::Writer, this};
+ DataHandle<eic::Cluster2DInfoCollection> m_outputInfo{"outputInfoCollection",
+ Gaudi::DataHandle::Writer, this};
// Pointer to the geometry service
SmartIF<IGeoSvc> m_geoSvc;
double m_depthCorr;
@@ -79,7 +86,9 @@ namespace Jug::Reco {
ClusterRecoCoG(const std::string& name, ISvcLocator* svcLoc) : GaudiAlgorithm(name, svcLoc)
{
declareProperty("inputHitCollection", m_inputHits, "");
+ declareProperty("inputProtoClusterCollection", m_inputProto, "");
declareProperty("outputClusterCollection", m_outputClusters, "");
+ declareProperty("outputInfoCollection", m_outputInfo, "");
}
StatusCode initialize() override
@@ -121,33 +130,37 @@ namespace Jug::Reco {
{
// input collections
const auto& hits = *m_inputHits.get();
+ const auto& proto = *m_inputProto.get();
auto& clusters = *m_outputClusters.createAndPut();
+ auto& info = *m_outputInfo.createAndPut();
// Create a map of clusterID --> associated hits by looping over our clustered hits
- std::map<int, std::vector<eic::ConstCalorimeterHit>> cluster_map;
- for (const auto& hit : hits) {
- const size_t idx = hit.clusterID();
- if (!cluster_map.count(idx)) {
- cluster_map[idx] = {};
+ std::map<int, std::vector<std::pair<eic::ConstProtoCluster,
+ eic::ConstCalorimeterHit>>> cluster_map;
+ for (const auto& pc : proto) {
+ const size_t clusterID = pc.clusterID();
+ if (!cluster_map.count(clusterID)) {
+ cluster_map[clusterID] = {};
}
- cluster_map[idx].push_back(hit);
+ size_t idx;
+ for (idx = 0; idx < hits.size(); ++idx) {
+ if (hits[idx].ID() == pc.hitID()) {
+ break;
+ }
+ }
+ if (idx >= hits.size()) {
+ continue;
+ }
+ cluster_map[clusterID].push_back({pc, hits[idx]});
}
- for (const auto& [idx, clhits] : cluster_map) {
- auto clhit = reconstruct(clhits);
- eic::Cluster cl{idx,
- static_cast<float>(clhit.energy() / m_sampFrac),
- clhit.energy(),
- static_cast<int>(clhits.size()),
- clhit.position(),
- {},
- cart_to_polar(clhit.position()),
- 0.,
- 0.};
+ for (const auto& [idx, hit_info] : cluster_map) {
+ auto cl = reconstruct(hit_info, idx);
debug() << cl.nhits() << " hits: " << cl.energy() / GeV << " GeV, (" << cl.position().x / mm
<< ", " << cl.position().y / mm << ", " << cl.position().z / mm << ")" << endmsg;
clusters.push_back(cl);
+ info.push_back({cl.ID(), cl.position(), cl.position().eta()});
}
return StatusCode::SUCCESS;
@@ -161,54 +174,52 @@ namespace Jug::Reco {
return eic::VectorPolar{r, std::acos(cart.z / r), std::atan2(cart.y, cart.x)};
}
- eic::CalorimeterHit reconstruct(const std::vector<eic::ConstCalorimeterHit>& hits) const
+ eic::Cluster reconstruct(const std::vector<std::pair<eic::ConstProtoCluster,
+ eic::ConstCalorimeterHit>>& hit_info,
+ const int idx) const
{
- eic::CalorimeterHit res;
+ eic::Cluster cl;
+ cl.ID(idx);
+ cl.nhits(hit_info.size());
+
// no hits
- debug() << "hit size = " << hits.size() << endmsg;
- if (hits.empty()) {
- return res;
- ;
+ debug() << "hit size = " << hit_info.size() << endmsg;
+ if (hit_info.empty()) {
+ return cl;
}
// calculate total energy, find the cell with the maximum energy deposit
float totalE = 0.;
float maxE = 0.;
- auto centerID = hits.begin()->cellID();
- for (const auto& hit : hits) {
- debug() << "hit energy = " << hit.energy() << endmsg;
- auto energy = hit.energy();
+ auto time = hit_info[0].second.time();
+ for (const auto& [proto, hit] : hit_info) {
+ debug() << "hit energy = " << hit.energy() << " hit weight: " << proto.weight() << endmsg;
+ auto energy = hit.energy() * proto.weight();
totalE += energy;
if (energy > maxE) {
maxE = energy;
- centerID = hit.cellID();
+ time = hit.time();
}
}
- res.cellID(centerID);
- res.energy(totalE);
+ cl.energy(totalE/m_sampFrac);
+ cl.energyError(0.);
+ cl.time(time);
// center of gravity with logarithmic weighting
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);
+ for (const auto& [proto, hit] : hit_info) {
+ float w = weightFunc(hit.energy()*proto.weight(), totalE, m_logWeightBase.value(), 0);
tw += w;
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(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();
- const auto lpos = alignment.worldToLocal(dd4hep::Position(res.position().x, res.position().y, res.position().z));
-
- // TODO: may need convert back to have depthCorrection in global positions
- res.local({lpos.x(), lpos.y(), lpos.z() + m_depthCorrection});
- return res;
+ cl.position(v.scale(1/tw));
+ cl.positionError({}); // @TODO: Covariance matrix
+
+ return cl;
}
};
diff --git a/JugReco/src/components/CrystalEndcapsReco.cpp b/JugReco/src/components/CrystalEndcapsReco.cpp
index ac58cd6341915f1ac1c91b2b6a638c91e1ed0819..9c28f8a62d176e38776311153dcb45f22446357a 100644
--- a/JugReco/src/components/CrystalEndcapsReco.cpp
+++ b/JugReco/src/components/CrystalEndcapsReco.cpp
@@ -70,10 +70,11 @@ namespace Jug::Reco {
auto& hits = *m_outputHitCollection.createAndPut();
// energy time reconstruction
+ int nhits = 0;
for (const auto& rh : rawhits) {
float energy = rh.amplitude() / 1.0e6; // convert keV -> GeV
if (energy >= (m_minModuleEdep / GeV)) {
- float time = rh.timeStamp();
+ float time = rh.time();
auto id = rh.cellID();
// global positions
auto gpos = m_geoSvc->cellIDPositionConverter()->position(id);
@@ -82,16 +83,9 @@ namespace Jug::Reco {
m_geoSvc->cellIDPositionConverter()->findContext(id)->volumePlacement().position();
// cell dimension
auto dim = m_geoSvc->cellIDPositionConverter()->cellDimensions(id);
- hits.push_back(eic::CalorimeterHit{id,
- -1,
- -1,
- -1,
- energy,
- time,
- {gpos.x(), gpos.y(), gpos.z()},
- {pos.x(), pos.y(), pos.z()},
- {dim[0], dim[1], 0.0},
- 0});
+ hits.push_back(eic::CalorimeterHit{id, nhits++, -1, -1, 0, energy, 0., time,
+ {gpos.x(), gpos.y(), gpos.z()},
+ {pos.x(), pos.y(), pos.z()}, {dim[0], dim[1], 0.}});
}
}
diff --git a/JugReco/src/components/EMCalReconstruction.cpp b/JugReco/src/components/EMCalReconstruction.cpp
index 7396392c0acddfe6ad4799f17dce8b5195da93d6..b4ff87369f8f15ae7122a61917136a87b18a9a5c 100644
--- a/JugReco/src/components/EMCalReconstruction.cpp
+++ b/JugReco/src/components/EMCalReconstruction.cpp
@@ -76,10 +76,11 @@ namespace Jug::Reco {
auto& hits = *m_outputHitCollection.createAndPut();
// energy time reconstruction
+ int nhits = 0;
for (const auto& rh : rawhits) {
float energy = rh.amplitude() / 1e6; // GeV
if (energy >= m_minModuleEdep) {
- float time = rh.timeStamp() / 1e6; // ns
+ float time = rh.time() / 1e6; // ns
auto id = rh.cellID();
// global positions
auto gpos = m_geoSvc->cellIDPositionConverter()->position(id);
@@ -90,15 +91,16 @@ namespace Jug::Reco {
auto dim = m_geoSvc->cellIDPositionConverter()->cellDimensions(id);
hits.push_back(eic::CalorimeterHit{
id,
+ nhits++,
-1,
-1,
- -1,
- static_cast<float>(energy / m_samplingFraction),
+ 0,
+ static_cast<float>(energy / m_samplingFraction), // Why are we applying the sampling fraction already here?
+ 0,
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},
- {dim[0] / dd4hep::mm, dim[1] / dd4hep::mm, 0.0},
- 0});
+ {dim[0] / dd4hep::mm, dim[1] / dd4hep::mm, 0.0}});
}
}
diff --git a/JugReco/src/components/EcalTungstenSamplingCluster.cpp b/JugReco/src/components/EcalTungstenSamplingCluster.cpp
index 4a59df2cd86de8d2e7c2d2e9b3444e7d070bfc97..39d0dce573cf00b1864d73e84649a26f1b26eea2 100644
--- a/JugReco/src/components/EcalTungstenSamplingCluster.cpp
+++ b/JugReco/src/components/EcalTungstenSamplingCluster.cpp
@@ -96,7 +96,7 @@ namespace Jug::Reco {
std::vector<bool>& visits) const
{
visits[index] = true;
- auto tot_edep = hits[index].energy();
+ auto tot_energy = hits[index].energy();
double pos_x = hits[index].position().x;
double pos_y = hits[index].position().y;
double pos_z = hits[index].position().z;
@@ -110,14 +110,14 @@ namespace Jug::Reco {
}
// Add up energy deposit based on the same z position and above energy threshold
if ((double)hits[i].position().z == ref_pos_z && hits[i].energy() > m_minModuleEdep / GeV) {
- tot_edep += hits[i].energy();
+ tot_energy += hits[i].energy();
visits[i] = true;
}
}
// Save info as a cluster
// TODO: position x and y determination
eic::Cluster cl;
- cl.edep(tot_edep);
+ cl.energy(tot_energy);
cl.position({pos_x, pos_y, pos_z});
return cl;
}
diff --git a/JugReco/src/components/EcalTungstenSamplingReco.cpp b/JugReco/src/components/EcalTungstenSamplingReco.cpp
index c9bd7975a92687bb811e3e38e905251a7fc4d800..e565f683036a08d692d9e58a2e781225790710ff 100644
--- a/JugReco/src/components/EcalTungstenSamplingReco.cpp
+++ b/JugReco/src/components/EcalTungstenSamplingReco.cpp
@@ -106,6 +106,7 @@ namespace Jug::Reco {
auto& hits = *m_outputHitCollection.createAndPut();
// energy time reconstruction
+ int nhits = 0;
for (const auto& rh : rawhits) {
// did not pass the threshold
if ((rh.amplitude() - m_pedMeanADC) < m_thresholdADC * m_pedSigmaADC) {
@@ -113,7 +114,7 @@ namespace Jug::Reco {
}
float energy = (rh.amplitude() - m_pedMeanADC) / (float)m_capADC *
m_dyRangeADC; // convert ADC -> energy
- float time = rh.timeStamp(); // ns
+ float time = rh.time(); // ns
auto id = rh.cellID();
int lid = ((id_dec != nullptr) & m_layerField.value().size())
? static_cast<int>(id_dec->get(id, layer_idx))
@@ -142,15 +143,16 @@ namespace Jug::Reco {
// << endmsg;
hits.push_back(
eic::CalorimeterHit{id,
- -1,
+ nhits++,
lid,
sid,
+ 0,
energy,
+ 0.,
time,
{gpos.x() / m_lUnit, gpos.y() / m_lUnit, gpos.z() / m_lUnit},
{pos.x() / m_lUnit, pos.y() / m_lUnit, pos.z() / m_lUnit},
- {dim[0], dim[1], dim[2]},
- 0});
+ {dim[0], dim[1], dim[2]}});
}
return StatusCode::SUCCESS;
diff --git a/JugReco/src/components/HCalReconstruction.cpp b/JugReco/src/components/HCalReconstruction.cpp
index b85fd222caf4e238deb07dc4e0d82b9556bbb2b7..e834af473b2c23a9f0ab19ceef79ae398b95deb0 100644
--- a/JugReco/src/components/HCalReconstruction.cpp
+++ b/JugReco/src/components/HCalReconstruction.cpp
@@ -76,10 +76,11 @@ namespace Jug::Reco {
auto& hits = *m_outputHitCollection.createAndPut();
// energy time reconstruction
+ int nhits = 0;
for (const auto& rh : rawhits) {
float energy = rh.amplitude() / 1e6; // GeV
if (energy >= m_minModuleEdep) {
- float time = rh.timeStamp() / 1e6; // ns
+ float time = rh.time() / 1e6; // ns
auto id = rh.cellID();
// global positions
auto gpos = m_geoSvc->cellIDPositionConverter()->position(id);
@@ -90,15 +91,16 @@ namespace Jug::Reco {
auto dim = m_geoSvc->cellIDPositionConverter()->cellDimensions(id);
hits.push_back(eic::CalorimeterHit{
id,
+ nhits++,
-1,
-1,
- -1,
+ 0,
static_cast<float>(energy / m_samplingFraction),
+ 0,
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},
- {dim[0] / dd4hep::mm, dim[1] / dd4hep::mm, 0.0},
- 0});
+ {dim[0] / dd4hep::mm, dim[1] / dd4hep::mm, 0.0}});
}
}
diff --git a/JugReco/src/components/ImagingClusterReco.cpp b/JugReco/src/components/ImagingClusterReco.cpp
index 84dde0bb64b4056a829b16c6b77ac9365acee324..b563f271e7233a84cb35dc08c7a7d4874f98be95 100644
--- a/JugReco/src/components/ImagingClusterReco.cpp
+++ b/JugReco/src/components/ImagingClusterReco.cpp
@@ -26,239 +26,242 @@
#include "JugBase/Utilities/Utils.hpp"
// Event Model related classes
-#include "eicd/ImagingClusterCollection.h"
-#include "eicd/ImagingLayerCollection.h"
-#include "eicd/ImagingPixelCollection.h"
+#include "eicd/CalorimeterHitCollection.h"
+#include "eicd/Cluster3DInfoCollection.h"
+#include "eicd/ClusterCollection.h"
+#include "eicd/ClusterLayerCollection.h"
+#include "eicd/ProtoClusterCollection.h"
+#include "eicd/VectorPolar.h"
+#include "eicd/VectorXYZ.h"
using namespace Gaudi::Units;
using namespace Eigen;
namespace Jug::Reco {
- /** Imaging cluster reconstruction.
- *
- * Reconstruct the cluster/layer info for imaging calorimeter
- * Logarithmic weighting is used to describe energy deposit in transverse direction
- *
- * \ingroup reco
- */
- class ImagingClusterReco : public GaudiAlgorithm {
- public:
- Gaudi::Property<double> m_sampFrac{this, "samplingFraction", 1.0};
- Gaudi::Property<int> m_trackStopLayer{this, "trackStopLayer", 9};
-
- DataHandle<eic::ImagingPixelCollection> m_inputHitCollection{"inputHitCollection",
- Gaudi::DataHandle::Reader, this};
- DataHandle<eic::ImagingLayerCollection> m_outputLayerCollection{
- "outputLayerCollection", Gaudi::DataHandle::Writer, this};
- DataHandle<eic::ImagingClusterCollection> m_outputClusterCollection{
- "outputClusterCollection", Gaudi::DataHandle::Reader, this};
-
- ImagingClusterReco(const std::string& name, ISvcLocator* svcLoc) : GaudiAlgorithm(name, svcLoc)
- {
- declareProperty("inputHitCollection", m_inputHitCollection, "");
- declareProperty("outputLayerCollection", m_outputLayerCollection, "");
- declareProperty("outputClusterCollection", m_outputClusterCollection, "");
+/** Imaging cluster reconstruction.
+ *
+ * Reconstruct the cluster/layer info for imaging calorimeter
+ * Logarithmic weighting is used to describe energy deposit in transverse direction
+ *
+ * \ingroup reco
+ */
+class ImagingClusterReco : public GaudiAlgorithm {
+public:
+ Gaudi::Property<double> m_sampFrac{this, "samplingFraction", 1.0};
+ Gaudi::Property<int> m_trackStopLayer{this, "trackStopLayer", 9};
+
+ DataHandle<eic::ProtoClusterCollection> m_inputProtoClusterCollection{"inputProtoClusterCollection",
+ Gaudi::DataHandle::Reader, this};
+ DataHandle<eic::CalorimeterHitCollection> m_inputHitCollection{"inputHitCollection", Gaudi::DataHandle::Reader, this};
+ DataHandle<eic::ClusterLayerCollection> m_outputLayerCollection{"outputLayerCollection", Gaudi::DataHandle::Writer,
+ this};
+ DataHandle<eic::ClusterCollection> m_outputClusterCollection{"outputClusterCollection", Gaudi::DataHandle::Reader,
+ this};
+ DataHandle<eic::Cluster3DInfoCollection> m_outputInfoCollection{"outputInfoCollection", Gaudi::DataHandle::Reader,
+ this};
+
+ ImagingClusterReco(const std::string& name, ISvcLocator* svcLoc) : GaudiAlgorithm(name, svcLoc) {
+ declareProperty("inputProtoClusterCollection", m_inputProtoClusterCollection, "");
+ declareProperty("inputHitCollection", m_inputHitCollection, "");
+ declareProperty("outputLayerCollection", m_outputLayerCollection, "");
+ declareProperty("outputClusterCollection", m_outputClusterCollection, "");
+ declareProperty("outputInfoCollection", m_outputInfoCollection, "");
+ }
+
+ StatusCode initialize() override {
+ if (GaudiAlgorithm::initialize().isFailure()) {
+ return StatusCode::FAILURE;
}
- StatusCode initialize() override
- {
- if (GaudiAlgorithm::initialize().isFailure()) {
- return StatusCode::FAILURE;
+ return StatusCode::SUCCESS;
+ }
+
+ StatusCode execute() override {
+ // input collections
+ const auto& proto = *m_inputProtoClusterCollection.get();
+ const auto& hits = *m_inputHitCollection.get();
+ // output collections
+ auto& layers = *m_outputLayerCollection.createAndPut();
+ auto& clusters = *m_outputClusterCollection.createAndPut();
+ auto& info = *m_outputInfoCollection.createAndPut();
+
+ // Create a map of clusterID --> associated ProtoCluster by looping over our clustered hits
+ std::map<int, std::vector<std::pair<eic::ConstProtoCluster,
+ eic::ConstCalorimeterHit>>> cluster_map;
+ for (const auto& pc : proto) {
+ const size_t clusterID = pc.clusterID();
+ if (!cluster_map.count(clusterID)) {
+ cluster_map[clusterID] = {};
}
-
- return StatusCode::SUCCESS;
- }
-
- StatusCode execute() override
- {
- // input collections
- const auto& hits = *m_inputHitCollection.get();
- // output collections
- auto& layers = *m_outputLayerCollection.createAndPut();
- auto& clusters = *m_outputClusterCollection.createAndPut();
-
- // Create a map of clusterID --> associated hits by looping over our clustered hits
- std::map<int, std::vector<eic::ConstImagingPixel>> cluster_map;
- for (const auto& hit : hits) {
- const size_t cid = hit.clusterID();
- if (!cluster_map.count(cid)) {
- cluster_map[cid] = {};
+ size_t idx;
+ for (idx = 0; idx < hits.size(); ++idx) {
+ if (hits[idx].ID() == pc.hitID()) {
+ break;
}
- cluster_map[cid].push_back(hit);
}
+ cluster_map[clusterID].push_back({pc, hits[idx]});
+ }
- for (const auto& [cid, clhits] : cluster_map) {
- // get cluster and associated layers
- auto cl = reconstruct_cluster(clhits, cid);
- auto cl_layers = reconstruct_cluster_layers(clhits, cid);
-
- // reconstruct cluster direction
- const auto [cl_theta, cl_phi] = fit_track(cl_layers, m_trackStopLayer);
- cl.cl_theta(cl_theta);
- cl.cl_phi(cl_phi);
+ for (const auto& [cid, hit_info] : cluster_map) {
+ // get cluster and associated layers
+ auto cl = reconstruct_cluster(hit_info, cid);
+ auto cl_layers = reconstruct_cluster_layers(hit_info, cid);
- // store layer and clusters on the datastore
- for (auto& layer : cl_layers) {
- layers.push_back(layer);
- // cl.addlayers(layer); // deprectated
- }
- clusters.push_back(cl);
- }
+ // reconstruct cluster direction
+ eic::Cluster3DInfo cl_info{cl.ID(), cl.position(), cl.position().eta(), fit_track(cl_layers, m_trackStopLayer)};
- // debug output
- for (const auto& cl : clusters) {
- debug() << fmt::format("Cluster {:d}: Edep = {:.3f} MeV, Dir = ({:.3f}, {:.3f}) deg",
- cl.clusterID(), cl.edep()*1000., cl.cl_theta() / M_PI * 180.,
- cl.cl_phi() / M_PI * 180.)
- << endmsg;
+ // store layer and clusters on the datastore
+ for (auto& layer : cl_layers) {
+ layer.ID(layers.size()); // unique ID for this clusterlayer
+ layers.push_back(layer);
+ // cl.addlayers(layer); // deprectated
}
-
- return StatusCode::SUCCESS;
+ clusters.push_back(cl);
+ info.push_back(cl_info);
}
- private:
- template <typename T>
- static inline T pow2(const T& x)
- {
- return x * x;
+ // debug output
+ int idx = 0;
+ for (const auto& cl : clusters) {
+ debug() << fmt::format("Cluster {:d}: Edep = {:.3f} MeV, Dir = ({:.3f}, {:.3f}) deg", cl.ID(), cl.energy() * 1000.,
+ info[idx].direction().theta / M_PI * 180., info[idx].direction().phi / M_PI * 180.)
+ << endmsg;
+ idx += 1;
}
- std::vector<eic::ImagingLayer>
- reconstruct_cluster_layers(const std::vector<eic::ConstImagingPixel>& hits, const int cid) const
- {
- // using map to have hits sorted by layer
- std::map<int, std::vector<eic::ConstImagingPixel>> layer_map;
- for (const auto& hit : hits) {
- auto lid = hit.layerID();
- if (!layer_map.count(lid)) {
- layer_map[lid] = {};
- }
- layer_map[lid].push_back(hit);
- }
-
- // create layers
- std::vector<eic::ImagingLayer> cl_layers;
- for (const auto& [lid, lhits] : layer_map) {
- auto layer = reconstruct_layer(lhits, cid, lid);
- cl_layers.push_back(layer);
+ return StatusCode::SUCCESS;
+ }
+
+private:
+ template <typename T> static inline T pow2(const T& x) { return x * x; }
+
+ std::vector<eic::ClusterLayer>
+ reconstruct_cluster_layers(const std::vector<std::pair<eic::ConstProtoCluster, eic::ConstCalorimeterHit>>& hit_info,
+ const int cid) const {
+ // using map to have hits sorted by layer
+ std::map<int, std::vector<std::pair<eic::ConstProtoCluster,
+ eic::ConstCalorimeterHit>>> layer_map;
+ for (const auto& [proto, hit] : hit_info) {
+ auto lid = hit.layer();
+ if (!layer_map.count(lid)) {
+ layer_map[lid] = {};
}
- return cl_layers;
+ layer_map[lid].push_back({proto, hit});
}
- eic::ImagingLayer reconstruct_layer(const std::vector<eic::ConstImagingPixel>& hits,
- const int cid, const int lid) const
- {
- // use full members initialization here so it could catch changes in ecid
- eic::ImagingLayer layer{cid, lid, static_cast<int>(hits.size()), 0., 0., 0., 0., {}, {}};
-
- // mean position and total edep
- double mx = 0.;
- double my = 0.;
- double mz = 0.;
- double edep = 0.;
- for (const auto& hit : hits) {
- mx += hit.position().x;
- my += hit.position().y;
- mz += hit.position().z;
- edep += hit.edep();
- }
-
- layer.position({mx / layer.nhits(), my / layer.nhits(), mz / layer.nhits()});
- layer.edep(edep);
-
- double radius = 0.;
- for (const auto& hit : hits) {
- radius += std::sqrt(pow2(hit.position().x - layer.position().x) + pow2(hit.position().y - layer.position().y) +
- pow2(hit.position().z - layer.position().z));
- }
- layer.radius(radius / layer.nhits());
- return layer;
+ // create layers
+ std::vector<eic::ClusterLayer> cl_layers;
+ for (const auto& [lid, layer_hit_info] : layer_map) {
+ auto layer = reconstruct_layer(layer_hit_info, cid, lid);
+ cl_layers.push_back(layer);
+ }
+ return cl_layers;
+ }
+
+ eic::ClusterLayer reconstruct_layer(const std::vector<std::pair<eic::ConstProtoCluster, eic::ConstCalorimeterHit>>& hit_info,
+ const int cid, const int lid) const {
+ // use full members initialization here so it could catch changes in ecid
+ eic::ClusterLayer layer{-1, cid, lid, static_cast<int>(hit_info.size()), 0, 0., 0., 0., 0., {}};
+
+ // mean position and total energy
+ eic::VectorXYZ pos;
+ double energy = 0.;
+ for (const auto& [proto, hit] : hit_info) {
+ pos = pos.add(hit.position());
+ energy += hit.energy() * proto.weight();
}
- eic::ImagingCluster reconstruct_cluster(const std::vector<eic::ConstImagingPixel>& hits,
- const int cid) const
- {
- eic::ImagingCluster cluster;
- cluster.clusterID(cid);
- // eta, phi center, weighted by energy
- double meta = 0.;
- double mphi = 0.;
- double edep = 0.;
- float r = 9999 * cm;
- for (const auto& hit : hits) {
- meta += hit.eta() * hit.edep();
- mphi += hit.polar().phi * hit.edep();
- edep += hit.edep();
- r = std::min(hit.polar().r, r);
- }
- const double eta = meta / edep;
- const double phi = mphi / edep;
- const double theta = 2. * std::atan(std::exp(-eta));
- cluster.nhits(hits.size());
- cluster.edep(edep);
- cluster.energy(edep / m_sampFrac); // simple energy reconstruction
- cluster.eta(eta);
- cluster.polar({r, phi, theta});
- // cartesian coordinates
- ROOT::Math::Polar3DVectorD polar{r, theta, (phi > M_PI ? phi - M_PI : phi)};
- cluster.position({polar.x(), polar.y(),polar.z()});
-
- // shower radius estimate (eta-phi plane)
- double radius = 0.;
- for (auto hit : hits) {
- radius += std::sqrt(pow2(hit.eta() - cluster.eta()) + pow2(hit.polar().phi - cluster.polar().phi));
- }
- cluster.radius(radius / cluster.nhits());
+ pos = pos.scale(1 / layer.nhits());
+ layer.energy(energy);
- return cluster;
+ double radius = 0.;
+ for (const auto& [proto, hit] : hit_info) {
+ radius += std::sqrt(pow2(hit.position().x - layer.position().x) + pow2(hit.position().y - layer.position().y) +
+ pow2(hit.position().z - layer.position().z));
}
-
- std::pair<double /* theta */, double /* phi */>
- fit_track(const std::vector<eic::ImagingLayer>& layers, const int stop_layer) const
- {
- int nrows = 0;
- double mx = 0.;
- double my = 0.;
- double mz = 0.;
- for (const auto& layer : layers) {
- if ((layer.layerID() <= stop_layer) && (layer.nhits() > 0)) {
- mx += layer.position().x;
- my += layer.position().y;
- mz += layer.position().z;
- nrows += 1;
- }
- }
- // cannot fit
- if (nrows < 2) {
- return {};
+ layer.radius(radius / layer.nhits());
+ return layer;
+ }
+
+ eic::Cluster
+ reconstruct_cluster(const std::vector<std::pair<eic::ConstProtoCluster, eic::ConstCalorimeterHit>>& hit_info,
+ const int cid) const {
+ eic::Cluster cluster;
+ cluster.ID(cid);
+ // eta, phi center, weighted by energy
+ double meta = 0.;
+ double mphi = 0.;
+ double energy = 0.;
+ float r = 9999 * cm;
+ for (const auto& [proto, hit] : hit_info) {
+ meta += hit.position().eta() * hit.energy() * proto.weight();
+ mphi += hit.position().phi() * hit.energy() * proto.weight();
+ energy += hit.energy() * proto.weight();
+ r = std::min(hit.position().r(), r);
+ }
+ const double eta = meta / energy;
+ const double phi = mphi / energy;
+ const double theta = 2. * std::atan(std::exp(-eta));
+ cluster.nhits(hit_info.size());
+ cluster.energy(energy / m_sampFrac); // simple energy reconstruction //DEPRECATED
+ eic::VectorPolar polar{r, theta, phi > M_PI ? phi - M_PI : phi};
+ cluster.position(polar);
+
+ // shower radius estimate (eta-phi plane)
+ double radius = 0.;
+ for (const auto& [proto, hit] : hit_info) {
+ radius += std::sqrt(pow2(hit.position().eta() - cluster.position().eta()) +
+ pow2(hit.position().phi() - cluster.position().phi()));
+ }
+ cluster.radius(radius / cluster.nhits());
+
+ return cluster;
+ }
+
+ eic::Direction fit_track(const std::vector<eic::ClusterLayer>& layers, const int stop_layer) const {
+ int nrows = 0;
+ double mx = 0.;
+ double my = 0.;
+ double mz = 0.;
+ for (const auto& layer : layers) {
+ if ((layer.layer() <= stop_layer) && (layer.nhits() > 0)) {
+ mx += layer.position().x;
+ my += layer.position().y;
+ mz += layer.position().z;
+ nrows += 1;
}
+ }
+ // cannot fit
+ if (nrows < 2) {
+ return {};
+ }
- mx /= nrows;
- my /= nrows;
- mz /= nrows;
- // fill position data
- MatrixXd pos(nrows, 3);
- int ir = 0;
- for (const auto& layer : layers) {
- if ((layer.layerID() <= stop_layer) && (layer.nhits() > 0)) {
- pos(ir, 0) = layer.position().x - mx;
- pos(ir, 1) = layer.position().y - my;
- pos(ir, 2) = layer.position().z - mz;
- ir += 1;
- }
+ mx /= nrows;
+ my /= nrows;
+ mz /= nrows;
+ // fill position data
+ MatrixXd pos(nrows, 3);
+ int ir = 0;
+ for (const auto& layer : layers) {
+ if ((layer.layer() <= stop_layer) && (layer.nhits() > 0)) {
+ pos(ir, 0) = layer.position().x - mx;
+ pos(ir, 1) = layer.position().y - my;
+ pos(ir, 2) = layer.position().z - mz;
+ ir += 1;
}
-
- JacobiSVD<MatrixXd> svd(pos, ComputeThinU | ComputeThinV);
- // debug() << pos << endmsg;
- // debug() << svd.matrixV() << endmsg;
- const auto dir = svd.matrixV().col(0);
- // theta and phi
- return {std::acos(dir(2)), std::atan2(dir(1), dir(0))};
}
- };
- DECLARE_COMPONENT(ImagingClusterReco)
+ JacobiSVD<MatrixXd> svd(pos, ComputeThinU | ComputeThinV);
+ // debug() << pos << endmsg;
+ // debug() << svd.matrixV() << endmsg;
+ const auto dir = svd.matrixV().col(0);
+ // theta and phi
+ return {std::acos(dir(2)), std::atan2(dir(1), dir(0))};
+ }
+};
+
+DECLARE_COMPONENT(ImagingClusterReco)
} // namespace Jug::Reco
diff --git a/JugReco/src/components/ImagingPixelMerger.cpp b/JugReco/src/components/ImagingPixelMerger.cpp
index e616ea87c448b9cc6a64afbae3bcefddef47912b..a971c58fcc73c58a348057ae30f2b947bcc672ad 100644
--- a/JugReco/src/components/ImagingPixelMerger.cpp
+++ b/JugReco/src/components/ImagingPixelMerger.cpp
@@ -29,8 +29,10 @@
#include "JugBase/Utilities/Utils.hpp"
// Event Model related classes
-#include "eicd/ImagingPixel.h"
-#include "eicd/ImagingPixelCollection.h"
+#include "eicd/VectorPolar.h"
+#include "eicd/VectorXYZ.h"
+#include "eicd/CalorimeterHit.h"
+#include "eicd/CalorimeterHitCollection.h"
using namespace Gaudi::Units;
@@ -60,9 +62,9 @@ namespace Jug::Reco {
Gaudi::Property<int> m_nLayers{this, "numberOfLayers", 20};
Gaudi::Property<double> m_etaSize{this, "etaSize", 0.001};
Gaudi::Property<double> m_phiSize{this, "phiSize", 0.001};
- DataHandle<eic::ImagingPixelCollection> m_inputHitCollection{"inputHitCollection",
+ DataHandle<eic::CalorimeterHitCollection> m_inputHitCollection{"inputHitCollection",
Gaudi::DataHandle::Reader, this};
- DataHandle<eic::ImagingPixelCollection> m_outputHitCollection{"outputHitCollection",
+ DataHandle<eic::CalorimeterHitCollection> m_outputHitCollection{"outputHitCollection",
Gaudi::DataHandle::Writer, this};
ImagingPixelMerger(const std::string& name, ISvcLocator* svcLoc) : GaudiAlgorithm(name, svcLoc)
@@ -90,7 +92,7 @@ namespace Jug::Reco {
// group the hits by layer
std::vector<std::unordered_map<std::pair<int, int>, double, pair_hash>> group_hits(m_nLayers);
for (const auto& h : hits) {
- auto k = h.layerID();
+ auto k = h.layer();
if ((int)k >= m_nLayers) {
continue;
}
@@ -105,9 +107,9 @@ namespace Jug::Reco {
auto it = layer.find(g);
// merge energy
if (it != layer.end()) {
- it->second += h.edep();
+ it->second += h.energy();
} else {
- layer[g] = h.edep();
+ layer[g] = h.energy();
}
}
@@ -132,12 +134,13 @@ namespace Jug::Reco {
if (k < grids.size()) {
grid = grids[k];
}
+ eic::VectorXYZ pos {eic::VectorPolar(0., 0., grid.phi)};
+ // @TODO: This seems incomplete...
auto h = mhits.create();
- h.edep(grid.energy);
- h.eta(grid.eta);
- h.polar({0., 0., grid.phi});
- h.layerID(i);
- h.hitID(k);
+ h.ID(k);
+ h.layer(i);
+ h.energy(grid.energy);
+ h.position(pos);
}
}
diff --git a/JugReco/src/components/ImagingPixelReco.cpp b/JugReco/src/components/ImagingPixelReco.cpp
index 90a3de4add3089f2066b44ab975badc7b3c47f83..aa6730ede888f3baeb0d645de2336159163ee5ed 100644
--- a/JugReco/src/components/ImagingPixelReco.cpp
+++ b/JugReco/src/components/ImagingPixelReco.cpp
@@ -22,7 +22,7 @@
#include "JugBase/IGeoSvc.h"
// Event Model related classes
-#include "eicd/ImagingPixelCollection.h"
+#include "eicd/CalorimeterHitCollection.h"
#include "eicd/RawCalorimeterHitCollection.h"
using namespace Gaudi::Units;
@@ -51,6 +51,8 @@ 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};
+ // Calibration!
+ Gaudi::Property<double> m_sampFrac{this, "samplingFraction", 1.0};
// unitless counterparts for the input parameters
double dyRangeADC;
@@ -58,8 +60,8 @@ namespace Jug::Reco {
// hits containers
DataHandle<eic::RawCalorimeterHitCollection> m_inputHitCollection{
"inputHitCollection", Gaudi::DataHandle::Reader, this};
- DataHandle<eic::ImagingPixelCollection> m_outputHitCollection{"outputHitCollection",
- Gaudi::DataHandle::Writer, this};
+ DataHandle<eic::CalorimeterHitCollection> m_outputHitCollection{"outputHitCollection",
+ Gaudi::DataHandle::Writer, this};
// Pointer to the geometry service
SmartIF<IGeoSvc> m_geoSvc;
@@ -115,13 +117,14 @@ namespace Jug::Reco {
auto& hits = *m_outputHitCollection.createAndPut();
// energy time reconstruction
+ int nhits = 0;
for (const auto& rh : rawhits) {
// did not pass the threshold
if ((rh.amplitude() - m_pedMeanADC) < m_thresholdADC * m_pedSigmaADC) {
continue;
}
- double edep = (rh.amplitude() - m_pedMeanADC) / (double)m_capADC * dyRangeADC; // convert ADC -> energy
- double time = rh.timeStamp(); // ns
+ double energy = (rh.amplitude() - m_pedMeanADC) / (double)m_capADC * dyRangeADC / m_sampFrac; // convert ADC -> energy
+ double time = rh.time() * 1.e-6; // ns
auto id = rh.cellID();
int lid = (int)id_dec->get(id, layer_idx);
int sid = (int)id_dec->get(id, sector_idx);
@@ -132,24 +135,19 @@ namespace Jug::Reco {
auto volman = m_geoSvc->detector()->volumeManager();
auto alignment = volman.lookupDetElement(id).nominal();
auto pos = alignment.worldToLocal(dd4hep::Position(gpos.x(), gpos.y(), gpos.z()));
- // polar coordinates
- double r = std::sqrt(gpos.x() * gpos.x() + gpos.y() * gpos.y() + gpos.z() * gpos.z());
- double th = std::acos(gpos.z() / r);
- double eta = -std::log(std::tan(th / 2.));
- double phi = std::atan2(gpos.y(), gpos.x());
-
- hits.push_back(eic::ImagingPixel{
- -1,
+
+ hits.push_back(eic::CalorimeterHit{
+ id,
+ nhits++,
lid,
sid,
- -1, // cluster id, layer id, sector id, hit id
- edep,
- time,
- eta, // edep, time, pseudo-rapidity
+ 0,
+ static_cast<float>(energy),
+ 0,
+ static_cast<float>(time),
{pos.x() / m_lUnit, pos.y() / m_lUnit, pos.z() / m_lUnit}, // local pos
{gpos.x() / m_lUnit, gpos.y() / m_lUnit, gpos.z() / m_lUnit}, // global pos
- {r, th, phi} // polar global pos
- });
+ {0, 0, 0}}); // @TODO: add dimension
}
return StatusCode::SUCCESS;
}
diff --git a/JugReco/src/components/ImagingTopoCluster.cpp b/JugReco/src/components/ImagingTopoCluster.cpp
index 916f26030dfcba65f0bb197ed9838ba9bbce3f34..4c70715b04c5591deef29164676c6934478a196f 100644
--- a/JugReco/src/components/ImagingTopoCluster.cpp
+++ b/JugReco/src/components/ImagingTopoCluster.cpp
@@ -27,8 +27,8 @@
#include "JugBase/IGeoSvc.h"
// Event Model related classes
-#include "eicd/ImagingClusterCollection.h"
-#include "eicd/ImagingPixelCollection.h"
+#include "eicd/ProtoClusterCollection.h"
+#include "eicd/CalorimeterHitCollection.h"
using namespace Gaudi::Units;
@@ -64,11 +64,11 @@ namespace Jug::Reco {
// minimum number of hits (to save this cluster)
Gaudi::Property<int> m_minClusterNhits{this, "minClusterNhits", 10};
// input hits collection
- DataHandle<eic::ImagingPixelCollection> m_inputHitCollection{"inputHitCollection",
- Gaudi::DataHandle::Reader, this};
+ DataHandle<eic::CalorimeterHitCollection> m_inputHitCollection{"inputHitCollection",
+ Gaudi::DataHandle::Reader, this};
// output clustered hits
- DataHandle<eic::ImagingPixelCollection> m_outputHitCollection{"outputHitCollection",
- Gaudi::DataHandle::Writer, this};
+ DataHandle<eic::ProtoClusterCollection> m_outputProtoClusterCollection{"outputProtoClusterCollection",
+ Gaudi::DataHandle::Writer, this};
// unitless counterparts of the input parameters
double localDistXY[2], layerDistEtaPhi[2], sectorDist;
@@ -77,7 +77,7 @@ namespace Jug::Reco {
ImagingTopoCluster(const std::string& name, ISvcLocator* svcLoc) : GaudiAlgorithm(name, svcLoc)
{
declareProperty("inputHitCollection", m_inputHitCollection, "");
- declareProperty("outputHitCollection", m_outputHitCollection, "");
+ declareProperty("outputProtoClusterCollection", m_outputProtoClusterCollection, "");
}
StatusCode initialize() override
@@ -126,11 +126,11 @@ namespace Jug::Reco {
// input collections
const auto& hits = *m_inputHitCollection.get();
// Create output collections
- auto& clustered_hits = *m_outputHitCollection.createAndPut();
+ auto& proto = *m_outputProtoClusterCollection.createAndPut();
// group neighboring hits
std::vector<bool> visits(hits.size(), false);
- std::vector<std::vector<eic::ImagingPixel>> groups;
+ std::vector<std::vector<eic::ConstCalorimeterHit>> groups;
for (size_t i = 0; i < hits.size(); ++i) {
/*
debug() << fmt::format("hit {:d}: local position = ({}, {}, {}), global position = ({}, {}, {})",
@@ -140,7 +140,7 @@ namespace Jug::Reco {
<< endmsg;
*/
// already in a group, or not energetic enough to form a cluster
- if (visits[i] || hits[i].edep() < minClusterCenterEdep) {
+ if (visits[i] || hits[i].energy() < minClusterCenterEdep) {
continue;
}
groups.emplace_back();
@@ -160,16 +160,15 @@ namespace Jug::Reco {
if (static_cast<int>(group.size()) < m_minClusterNhits.value()) {
continue;
}
- double edep = 0.;
+ double energy = 0.;
for (const auto& hit : group) {
- edep += hit.edep();
+ energy += hit.energy();
}
- if (edep < minClusterEdep) {
+ if (energy < minClusterEdep) {
continue;
}
- for (auto hit : group) {
- hit.clusterID(clusterID);
- clustered_hits.push_back(hit);
+ for (const auto& hit : group) {
+ proto.create(hit.ID(), clusterID, 1.);
}
clusterID += 1;
}
@@ -185,23 +184,23 @@ namespace Jug::Reco {
}
// helper function to group hits
- bool is_neighbor(const eic::ConstImagingPixel& h1, const eic::ConstImagingPixel& h2) const
+ bool is_neighbor(const eic::ConstCalorimeterHit& h1, const eic::ConstCalorimeterHit& h2) const
{
// different sectors, simple distance check
- if (h1.sectorID() != h2.sectorID()) {
+ 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;
}
// layer check
- int ldiff = std::abs(h1.layerID() - h2.layerID());
+ 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.eta() - h2.eta()) <= layerDistEtaPhi[0]) &&
- (std::abs(h1.polar().phi - h2.polar().phi) <= layerDistEtaPhi[1]);
+ return (std::abs(h1.position().eta() - h2.position().eta()) <= layerDistEtaPhi[0]) &&
+ (std::abs(h1.position().phi() - h2.position().phi()) <= layerDistEtaPhi[1]);
}
// not in adjacent layers
@@ -209,24 +208,20 @@ namespace Jug::Reco {
}
// grouping function with Depth-First Search
- void dfs_group(std::vector<eic::ImagingPixel>& group, int idx,
- const eic::ImagingPixelCollection& hits, std::vector<bool>& visits) const
+ void dfs_group(std::vector<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].edep() < minClusterHitEdep) {
+ if (hits[idx].energy() < 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(),
- hits[idx].polar()};
- group.push_back(hit);
+ group.push_back(hits[idx]);
visits[idx] = true;
for (size_t i = 0; i < hits.size(); ++i) {
// visited, or not a neighbor
- if (visits[i] || !is_neighbor(hit, hits[i])) {
+ if (visits[i] || !is_neighbor(hits[idx], hits[i])) {
continue;
}
dfs_group(group, i, hits, visits);
diff --git a/JugReco/src/components/SimpleClustering.cpp b/JugReco/src/components/SimpleClustering.cpp
index 8863a9357012b460b460da36ff2e320e0bb85cb8..c20ec353765b0d9ec255d6c07a8f2ea38e16be89 100644
--- a/JugReco/src/components/SimpleClustering.cpp
+++ b/JugReco/src/components/SimpleClustering.cpp
@@ -19,6 +19,7 @@
// Event Model related classes
#include "eicd/CalorimeterHitCollection.h"
#include "eicd/ClusterCollection.h"
+#include "eicd/ProtoClusterCollection.h"
#include "eicd/RawCalorimeterHitCollection.h"
using namespace Gaudi::Units;
@@ -32,12 +33,15 @@ namespace Jug::Reco {
class SimpleClustering : public GaudiAlgorithm {
public:
using RecHits = eic::CalorimeterHitCollection;
+ using ProtoClusters = eic::ProtoClusterCollection;
using Clusters = eic::ClusterCollection;
- DataHandle<RecHits> m_inputHitCollection{"inputHitCollection", Gaudi::DataHandle::Reader, this};
- DataHandle<Clusters> m_outputClusters{"outputClusters", Gaudi::DataHandle::Writer, this};
- Gaudi::Property<double> m_minModuleEdep{this, "minModuleEdep", 5.0 * MeV};
- Gaudi::Property<double> m_maxDistance{this, "maxDistance", 20.0 * cm};
+ DataHandle<RecHits> m_inputHitCollection{"inputHitCollection", Gaudi::DataHandle::Reader, this};
+ DataHandle<ProtoClusters> m_outputProtoClusters{"outputProtoCluster", Gaudi::DataHandle::Writer, this};
+ DataHandle<Clusters> m_outputClusters{"outputClusterCollection", Gaudi::DataHandle::Writer, this};
+
+ Gaudi::Property<double> m_minModuleEdep{this, "minModuleEdep", 5.0 * MeV};
+ Gaudi::Property<double> m_maxDistance{this, "maxDistance", 20.0 * cm};
/// Pointer to the geometry service
SmartIF<IGeoSvc> m_geoSvc;
@@ -45,7 +49,8 @@ namespace Jug::Reco {
SimpleClustering(const std::string& name, ISvcLocator* svcLoc) : GaudiAlgorithm(name, svcLoc)
{
declareProperty("inputHitCollection", m_inputHitCollection, "");
- declareProperty("outputClusters", m_outputClusters, "Output clusters");
+ declareProperty("outputProtoClusterCollection", m_outputClusters, "Output proto clusters");
+ declareProperty("outputClusterCollection", m_outputClusters, "Output clusters");
}
StatusCode initialize() override
@@ -67,19 +72,19 @@ namespace Jug::Reco {
// input collections
const auto& hits = *m_inputHitCollection.get();
// Create output collections
+ auto& proto = *m_outputProtoClusters.createAndPut();
auto& clusters = *m_outputClusters.createAndPut();
- std::vector<eic::CalorimeterHit> the_hits;
- std::vector<eic::CalorimeterHit> remaining_hits;
+ std::vector<eic::ConstCalorimeterHit> the_hits;
+ std::vector<eic::ConstCalorimeterHit> remaining_hits;
double max_dist = m_maxDistance.value() / mm;
double min_energy = m_minModuleEdep.value() / GeV;
- eic::CalorimeterHit ref_hit;
+ eic::ConstCalorimeterHit ref_hit;
bool have_ref = false;
- for (const auto& ch : hits) {
- const eic::CalorimeterHit h{ch.cellID(), ch.clusterID(), ch.layerID(), ch.sectorID(), ch.energy(),
- ch.time(), ch.position(), ch.local(), ch.dimension(), 1};
+ for (const auto& h : hits) {
+ //const eic::CalorimeterHit h = ch.clone();
if (!have_ref || h.energy() > ref_hit.energy()) {
ref_hit = h;
have_ref = true;
@@ -91,31 +96,31 @@ namespace Jug::Reco {
while (have_ref && ref_hit.energy() > min_energy) {
- std::vector<eic::CalorimeterHit> cluster_hits;
+ std::vector<eic::ConstCalorimeterHit> cluster_hits;
for (const auto& h : the_hits) {
if (std::hypot(h.position().x - ref_hit.position().x, h.position().y - ref_hit.position().y,
h.position().z - ref_hit.position().z) < max_dist) {
cluster_hits.push_back(h);
} else {
- remaining_hits.push_back(h);
- }
+ remaining_hits.push_back(h); }
}
debug() << " cluster size " << cluster_hits.size() << endmsg;
double total_energy =
std::accumulate(std::begin(cluster_hits), std::end(cluster_hits), 0.0,
- [](double t, const eic::CalorimeterHit& h1) { return (t + h1.energy()); });
+ [](double t, const eic::ConstCalorimeterHit& h1) { return (t + h1.energy()); });
debug() << " total_energy = " << total_energy << endmsg;
- eic::Cluster cluster0;
+ eic::Cluster cl;
+ cl.ID(clusters.size());
+ cl.nhits(cluster_hits.size());
for (const auto& h : cluster_hits) {
- cluster0.energy(cluster0.energy() + h.energy());
- cluster0.position({cluster0.position().x + h.energy() * h.position().x / total_energy,
- cluster0.position().y + h.energy() * h.position().y / total_energy,
- cluster0.position().z + h.energy() * h.position().z / total_energy});
+ cl.energy(cl.energy() + h.energy());
+ cl.position(cl.position().add(h.position().scale(h.energy()/total_energy)));
+ proto.create(h.ID(), cl.ID());
}
- clusters.push_back(cluster0);
+ clusters.push_back(cl);
have_ref = false;
if ((remaining_hits.size() > 5) && (clusters.size() < 10)) {
diff --git a/JugReco/src/components/TopologicalCellCluster.cpp b/JugReco/src/components/TopologicalCellCluster.cpp
index fc858af64e5c9cda6b22592b6eaf001842c0f3d5..da8f1e96af9d2c43e24aa9a6d04cd8558cddf2bd 100644
--- a/JugReco/src/components/TopologicalCellCluster.cpp
+++ b/JugReco/src/components/TopologicalCellCluster.cpp
@@ -8,6 +8,9 @@
* References: https://arxiv.org/pdf/1603.02934.pdf
*
*/
+
+// @FIXME: Deprecated?
+
#include "fmt/format.h"
#include <algorithm>
@@ -64,8 +67,8 @@ namespace Jug::Reco {
// maximum global distance to be considered as neighbors in different sectors
Gaudi::Property<double> m_adjSectorDist{this, "adjSectorDist", 1.0 * cm};
// minimum cluster center energy (to be considered as a seed for cluster)
- // @TODO One can not simply find a center by edep with extremely fine granularity.
- // Projecting to (eta, phi) with crude pixel size may help determine the center edep, which can
+ // @TODO One can not simply find a center by energy with extremely fine granularity.
+ // Projecting to (eta, phi) with crude pixel size may help determine the center energy, which can
// happen in the following reconstruction step
Gaudi::Property<double> m_minClusterCenterEdep{this, "minClusterCenterEdep", 50 * keV};
// input hits collection
@@ -154,7 +157,7 @@ namespace Jug::Reco {
// group neighboring hits
std::vector<bool> visits(hits.size(), false);
- std::vector<std::vector<eic::CalorimeterHit>> groups;
+ std::vector<std::vector<eic::ConstCalorimeterHit>> 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].energy() < m_minClusterCenterEdep) {
@@ -222,19 +225,14 @@ namespace Jug::Reco {
}
// grouping function with Depth-First Search
- void dfs_group(std::vector<eic::CalorimeterHit>& group, int idx,
+ void dfs_group(std::vector<eic::ConstCalorimeterHit>& group, int idx,
const eic::CalorimeterHitCollection& hits, std::vector<bool>& visits) const
{
- const eic::CalorimeterHit hit{hits[idx].cellID(), hits[idx].clusterID(),
- hits[idx].layerID(), hits[idx].sectorID(),
- hits[idx].energy(), hits[idx].time(),
- hits[idx].position(), hits[idx].local(),
- hits[idx].dimension(), 1};
- group.push_back(hit);
+ group.push_back(hits[idx]);
visits[idx] = true;
for (size_t i = 0; i < hits.size(); ++i) {
// visited, or not a neighbor
- if (visits[i] || !is_neighbor(hit, hits[i])) {
+ if (visits[i] || !is_neighbor(hits[idx], hits[i])) {
continue;
}
dfs_group(group, i, hits, visits);
diff --git a/JugTrack/src/components/TrackParamImagingClusterInit.cpp b/JugTrack/src/components/TrackParamImagingClusterInit.cpp
index a0eb8d5c81b1a9604f144a110467f6be677831b8..d893cde9754b979968e1ff60142818e6bf28b549 100644
--- a/JugTrack/src/components/TrackParamImagingClusterInit.cpp
+++ b/JugTrack/src/components/TrackParamImagingClusterInit.cpp
@@ -14,7 +14,7 @@
#include "Acts/Definitions/Common.hpp"
#include "eicd/TrackerHitCollection.h"
-#include "eicd/ImagingClusterCollection.h"
+#include "eicd/ClusterCollection.h"
#include "Math/Vector3D.h"
#include "Acts/Surfaces/PerigeeSurface.hpp"
@@ -41,7 +41,7 @@ namespace Jug::Reco {
*/
class TrackParamImagingClusterInit : public GaudiAlgorithm {
public:
- using ImagingClusters = eic::ImagingClusterCollection;
+ using ImagingClusters = eic::ClusterCollection;
DataHandle<ImagingClusters> m_inputClusters{"inputClusters", Gaudi::DataHandle::Reader, this};
DataHandle<TrackParametersContainer> m_outputInitialTrackParameters{"outputInitialTrackParameters",
@@ -65,7 +65,7 @@ namespace Jug::Reco {
StatusCode execute() override {
// input collection
- const eic::ImagingClusterCollection* clusters = m_inputClusters.get();
+ const eic::ClusterCollection* clusters = m_inputClusters.get();
// Create output collections
auto init_trk_params = m_outputInitialTrackParameters.createAndPut();
@@ -77,10 +77,11 @@ namespace Jug::Reco {
using Acts::UnitConstants::ns;
double p = c.energy()*GeV;
+ // FIXME hardcoded value
if( p < 0.1*GeV) {
continue;
}
- double len = std::hypot( c.position().x , c.position().y , c.position().z );
+ 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