File indexing completed on 2024-09-27 07:03:08
0001
0002
0003
0004 #include <DD4hep/Detector.h>
0005 #include <DD4hep/IDDescriptor.h>
0006 #include <DD4hep/Readout.h>
0007 #include <Evaluator/DD4hepUnits.h>
0008 #include <algorithms/geo.h>
0009 #include <catch2/catch_test_macros.hpp>
0010 #include <catch2/generators/catch_generators.hpp>
0011 #include <catch2/matchers/catch_matchers.hpp>
0012 #include <catch2/matchers/catch_matchers_floating_point.hpp>
0013 #include <edm4eic/CalorimeterHitCollection.h>
0014 #include <edm4eic/ProtoClusterCollection.h>
0015 #include <edm4hep/Vector3f.h>
0016 #include <podio/RelationRange.h>
0017 #include <spdlog/common.h>
0018 #include <spdlog/logger.h>
0019 #include <spdlog/spdlog.h>
0020 #include <gsl/pointers>
0021 #include <limits>
0022 #include <memory>
0023 #include <utility>
0024 #include <vector>
0025
0026 #include "algorithms/calorimetry/CalorimeterIslandCluster.h"
0027 #include "algorithms/calorimetry/CalorimeterIslandClusterConfig.h"
0028
0029 using eicrecon::CalorimeterIslandCluster;
0030 using eicrecon::CalorimeterIslandClusterConfig;
0031
0032 TEST_CASE( "the clustering algorithm runs", "[CalorimeterIslandCluster]" ) {
0033 CalorimeterIslandCluster algo("CalorimeterIslandCluster");
0034
0035 std::shared_ptr<spdlog::logger> logger = spdlog::default_logger()->clone("CalorimeterIslandCluster");
0036 logger->set_level(spdlog::level::trace);
0037
0038 CalorimeterIslandClusterConfig cfg;
0039 cfg.minClusterHitEdep = 0. * dd4hep::GeV;
0040 cfg.minClusterCenterEdep = 0. * dd4hep::GeV;
0041
0042 auto detector = algorithms::GeoSvc::instance().detector();
0043 auto id_desc = detector->readout("MockCalorimeterHits").idSpec();
0044
0045 SECTION( "without splitting" ) {
0046 bool use_adjacencyMatrix = GENERATE(false, true);
0047 cfg.splitCluster = false;
0048 if (use_adjacencyMatrix) {
0049 cfg.adjacencyMatrix = "abs(x_1 - x_2) + abs(y_1 - y_2) == 1";
0050 cfg.readout = "MockCalorimeterHits";
0051 } else {
0052 cfg.localDistXY = {1 * dd4hep::mm, 1 * dd4hep::mm};
0053 }
0054 algo.applyConfig(cfg);
0055 algo.init();
0056
0057 SECTION( "on a single cell" ) {
0058 edm4eic::CalorimeterHitCollection hits_coll;
0059 hits_coll.create(
0060 id_desc.encode({{"system", 255}, {"x", 0}, {"y", 0}}),
0061 5.0,
0062 0.0,
0063 0.0,
0064 0.0,
0065 edm4hep::Vector3f(0.0, 0.0, 0.0),
0066 edm4hep::Vector3f(0.0, 0.0, 0.0),
0067 0,
0068 0,
0069 edm4hep::Vector3f(0.0, 0.0, 0.0)
0070 );
0071 auto protoclust_coll = std::make_unique<edm4eic::ProtoClusterCollection>();
0072 algo.process({&hits_coll}, {protoclust_coll.get()});
0073
0074 REQUIRE( (*protoclust_coll).size() == 1 );
0075 REQUIRE( (*protoclust_coll)[0].hits_size() == 1 );
0076 REQUIRE( (*protoclust_coll)[0].weights_size() == 1 );
0077 }
0078
0079 SECTION( "on two separated cells" ) {
0080 edm4eic::CalorimeterHitCollection hits_coll;
0081 hits_coll.create(
0082 0,
0083 5.0,
0084 0.0,
0085 0.0,
0086 0.0,
0087 edm4hep::Vector3f(0.0, 0.0, 0.0),
0088 edm4hep::Vector3f(1.0, 1.0, 0.0),
0089 0,
0090 0,
0091 edm4hep::Vector3f(0.0, 0.0, 0.0)
0092 );
0093 hits_coll.create(
0094 1,
0095 6.0,
0096 0.0,
0097 0.0,
0098 0.0,
0099 edm4hep::Vector3f(0.0, 0.0, 0.0),
0100 edm4hep::Vector3f(1.0, 1.0, 0.0),
0101 0,
0102 0,
0103 edm4hep::Vector3f(1.1 , 1.1 , 0.0)
0104 );
0105 auto protoclust_coll = std::make_unique<edm4eic::ProtoClusterCollection>();
0106 algo.process({&hits_coll}, {protoclust_coll.get()});
0107
0108 REQUIRE( (*protoclust_coll).size() == 2 );
0109 REQUIRE( (*protoclust_coll)[0].hits_size() == 1 );
0110 REQUIRE( (*protoclust_coll)[0].weights_size() == 1 );
0111 REQUIRE( (*protoclust_coll)[1].hits_size() == 1 );
0112 REQUIRE( (*protoclust_coll)[1].weights_size() == 1 );
0113 }
0114
0115 SECTION( "on two adjacent cells" ) {
0116 edm4eic::CalorimeterHitCollection hits_coll;
0117 hits_coll.create(
0118 id_desc.encode({{"system", 255}, {"x", 0}, {"y", 0}}),
0119 5.0,
0120 0.0,
0121 0.0,
0122 0.0,
0123 edm4hep::Vector3f(0.0, 0.0, 0.0),
0124 edm4hep::Vector3f(1.0, 1.0, 0.0),
0125 0,
0126 0,
0127 edm4hep::Vector3f(0.0, 0.0, 0.0)
0128 );
0129 hits_coll.create(
0130 id_desc.encode({{"system", 255}, {"x", 1}, {"y", 0}}),
0131 6.0,
0132 0.0,
0133 0.0,
0134 0.0,
0135 edm4hep::Vector3f(0.0, 0.0, 0.0),
0136 edm4hep::Vector3f(1.0, 1.0, 0.0),
0137 0,
0138 0,
0139 edm4hep::Vector3f(0.9 , 0.9 , 0.0)
0140 );
0141 auto protoclust_coll = std::make_unique<edm4eic::ProtoClusterCollection>();
0142 algo.process({&hits_coll}, {protoclust_coll.get()});
0143
0144 REQUIRE( (*protoclust_coll).size() == 1 );
0145 REQUIRE( (*protoclust_coll)[0].hits_size() == 2 );
0146 REQUIRE( (*protoclust_coll)[0].weights_size() == 2 );
0147 }
0148 }
0149
0150 SECTION( "run on three adjacent cells" ) {
0151 bool use_adjacencyMatrix = GENERATE(false, true);
0152 if (use_adjacencyMatrix) {
0153 cfg.adjacencyMatrix = "abs(x_1 - x_2) + abs(y_1 - y_2) == 1";
0154 } else {
0155 cfg.localDistXY = {1 * dd4hep::mm, 1 * dd4hep::mm};
0156 }
0157 bool disalow_diagonal_peaks = GENERATE(false, true);
0158 if (disalow_diagonal_peaks) {
0159 cfg.peakNeighbourhoodMatrix = "max(abs(x_1 - x_2), abs(y_1 - y_2)) == 1";
0160 } else {
0161 cfg.peakNeighbourhoodMatrix = "abs(x_1 - x_2) + abs(y_1 - y_2) == 1";
0162 }
0163 cfg.readout = "MockCalorimeterHits";
0164
0165 cfg.splitCluster = GENERATE(false, true);
0166 if (cfg.splitCluster) {
0167 cfg.transverseEnergyProfileMetric = "localDistXY";
0168 cfg.transverseEnergyProfileScale = std::numeric_limits<decltype(cfg.transverseEnergyProfileScale)>::infinity();
0169 }
0170 cfg.localDistXY = {1 * dd4hep::mm, 1 * dd4hep::mm};
0171 algo.applyConfig(cfg);
0172 algo.init();
0173
0174 edm4eic::CalorimeterHitCollection hits_coll;
0175 hits_coll.create(
0176 id_desc.encode({{"system", 255}, {"x", 0}, {"y", 0}}),
0177 5.0,
0178 0.0,
0179 0.0,
0180 0.0,
0181 edm4hep::Vector3f(0.0, 0.0, 0.0),
0182 edm4hep::Vector3f(1.0, 1.0, 0.0),
0183 0,
0184 0,
0185 edm4hep::Vector3f(0.0, 0.0, 0.0)
0186 );
0187 hits_coll.create(
0188 id_desc.encode({{"system", 255}, {"x", 1}, {"y", 0}}),
0189 1.0,
0190 0.0,
0191 0.0,
0192 0.0,
0193 edm4hep::Vector3f(0.0, 0.0, 0.0),
0194 edm4hep::Vector3f(1.0, 1.0, 0.0),
0195 0,
0196 0,
0197 edm4hep::Vector3f(0.9 , 0.9 , 0.0)
0198 );
0199
0200 bool test_diagonal_cluster = GENERATE(false, true);
0201
0202
0203
0204
0205
0206
0207
0208 hits_coll.create(
0209 id_desc.encode({{"system", 255}, {"x", test_diagonal_cluster ? 1 : 2}, {"y", test_diagonal_cluster ? 1 : 0}}),
0210 6.0,
0211 0.0,
0212 0.0,
0213 0.0,
0214 edm4hep::Vector3f(0.0, 0.0, 0.0),
0215 edm4hep::Vector3f(1.0, 1.0, 0.0),
0216 0,
0217 0,
0218 edm4hep::Vector3f(1.8 , 1.8 , 0.0)
0219 );
0220 auto protoclust_coll = std::make_unique<edm4eic::ProtoClusterCollection>();
0221 algo.process({&hits_coll}, {protoclust_coll.get()});
0222
0223 bool expect_two_peaks = cfg.splitCluster;
0224 if (cfg.splitCluster && disalow_diagonal_peaks) {
0225 expect_two_peaks = not test_diagonal_cluster;
0226 }
0227 if (expect_two_peaks) {
0228 REQUIRE( (*protoclust_coll).size() == 2 );
0229 REQUIRE( (*protoclust_coll)[0].hits_size() == 3 );
0230 REQUIRE( (*protoclust_coll)[0].weights_size() == 3 );
0231 for (double weight : (*protoclust_coll)[0].getWeights()) {
0232 double energy_fraction = hits_coll[0].getEnergy() / (hits_coll[0].getEnergy() + hits_coll[2].getEnergy());
0233 REQUIRE_THAT( weight, Catch::Matchers::WithinAbs(energy_fraction, 1e-5) );
0234 }
0235 REQUIRE( (*protoclust_coll)[1].hits_size() == 3 );
0236 REQUIRE( (*protoclust_coll)[1].weights_size() == 3 );
0237 for (double weight : (*protoclust_coll)[1].getWeights()) {
0238 double energy_fraction = hits_coll[2].getEnergy() / (hits_coll[0].getEnergy() + hits_coll[2].getEnergy());
0239 REQUIRE_THAT( weight, Catch::Matchers::WithinAbs(energy_fraction, 1e-5) );
0240 }
0241 } else {
0242 REQUIRE( (*protoclust_coll).size() == 1 );
0243 REQUIRE( (*protoclust_coll)[0].hits_size() == 3 );
0244 REQUIRE( (*protoclust_coll)[0].weights_size() == 3 );
0245 }
0246 }
0247 }