tests/algorithms_test/calorimetry_CalorimeterIslandCluster.cc

0001 // SPDX-License-Identifier: LGPL-3.0-or-later
0002 // Copyright (C) 2023, Dmitry Kalinkin
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 <string>
0024 #include <utility>
0025 #include <variant>
0026 #include <vector>
0027
0028 #include "algorithms/calorimetry/CalorimeterIslandCluster.h"
0029 #include "algorithms/calorimetry/CalorimeterIslandClusterConfig.h"
0030
0031 using eicrecon::CalorimeterIslandCluster;
0032 using eicrecon::CalorimeterIslandClusterConfig;
0033
0034 TEST_CASE("the clustering algorithm runs", "[CalorimeterIslandCluster]") {
0035   CalorimeterIslandCluster algo("CalorimeterIslandCluster");
0036
0037   std::shared_ptr<spdlog::logger> logger =
0038       spdlog::default_logger()->clone("CalorimeterIslandCluster");
0039   logger->set_level(spdlog::level::trace);
0040
0041   CalorimeterIslandClusterConfig cfg;
0042   cfg.minClusterHitEdep    = 0. * dd4hep::GeV;
0043   cfg.minClusterCenterEdep = 0. * dd4hep::GeV;
0044
0045   auto detector = algorithms::GeoSvc::instance().detector();
0046   auto id_desc  = detector->readout("MockCalorimeterHits").idSpec();
0047
0048   SECTION("without splitting") {
0049     bool use_adjacencyMatrix = GENERATE(false, true);
0050     cfg.splitCluster         = false;
0051     if (use_adjacencyMatrix) {
0052       cfg.adjacencyMatrix = "abs(x_1 - x_2) + abs(y_1 - y_2) == 1";
0053       cfg.readout         = "MockCalorimeterHits";
0054     } else {
0055       cfg.localDistXY = {1 * dd4hep::mm, 1 * dd4hep::mm};
0056     }
0057     algo.applyConfig(cfg);
0058     algo.init();
0059
0060     SECTION("on a single cell") {
0061       edm4eic::CalorimeterHitCollection hits_coll;
0062       hits_coll.create(
0063           id_desc.encode({{"system", 255}, {"x", 0}, {"y", 0}}), // std::uint64_t cellID,
0064           5.0,                                                   // float energy,
0065           0.0,                                                   // float energyError,
0066           0.0,                                                   // float time,
0067           0.0,                                                   // float timeError,
0068           edm4hep::Vector3f(0.0, 0.0, 0.0),                      // edm4hep::Vector3f position,
0069           edm4hep::Vector3f(0.0, 0.0, 0.0),                      // edm4hep::Vector3f dimension,
0070           0,                                                     // std::int32_t sector,
0071           0,                                                     // std::int32_t layer,
0072           edm4hep::Vector3f(0.0, 0.0, 0.0)                       // edm4hep::Vector3f local
0073       );
0074       auto protoclust_coll = std::make_unique<edm4eic::ProtoClusterCollection>();
0075       algo.process({&hits_coll}, {protoclust_coll.get()});
0076
0077       REQUIRE((*protoclust_coll).size() == 1);
0078       REQUIRE((*protoclust_coll)[0].hits_size() == 1);
0079       REQUIRE((*protoclust_coll)[0].weights_size() == 1);
0080     }
0081
0082     SECTION("on two separated cells") {
0083       edm4eic::CalorimeterHitCollection hits_coll;
0084       hits_coll.create(0,                                // std::uint64_t cellID,
0085                        5.0,                              // float energy,
0086                        0.0,                              // float energyError,
0087                        0.0,                              // float time,
0088                        0.0,                              // float timeError,
0089                        edm4hep::Vector3f(0.0, 0.0, 0.0), // edm4hep::Vector3f position,
0090                        edm4hep::Vector3f(1.0, 1.0, 0.0), // edm4hep::Vector3f dimension,
0091                        0,                                // std::int32_t sector,
0092                        0,                                // std::int32_t layer,
0093                        edm4hep::Vector3f(0.0, 0.0, 0.0)  // edm4hep::Vector3f local
0094       );
0095       hits_coll.create(1,                                // std::uint64_t cellID,
0096                        6.0,                              // float energy,
0097                        0.0,                              // float energyError,
0098                        0.0,                              // float time,
0099                        0.0,                              // float timeError,
0100                        edm4hep::Vector3f(0.0, 0.0, 0.0), // edm4hep::Vector3f position,
0101                        edm4hep::Vector3f(1.0, 1.0, 0.0), // edm4hep::Vector3f dimension,
0102                        0,                                // std::int32_t sector,
0103                        0,                                // std::int32_t layer,
0104                        edm4hep::Vector3f(1.1 /* mm */, 1.1 /* mm */, 0.0) // edm4hep::Vector3f local
0105       );
0106       auto protoclust_coll = std::make_unique<edm4eic::ProtoClusterCollection>();
0107       algo.process({&hits_coll}, {protoclust_coll.get()});
0108
0109       REQUIRE((*protoclust_coll).size() == 2);
0110       REQUIRE((*protoclust_coll)[0].hits_size() == 1);
0111       REQUIRE((*protoclust_coll)[0].weights_size() == 1);
0112       REQUIRE((*protoclust_coll)[1].hits_size() == 1);
0113       REQUIRE((*protoclust_coll)[1].weights_size() == 1);
0114     }
0115
0116     SECTION("on two adjacent cells") {
0117       edm4eic::CalorimeterHitCollection hits_coll;
0118       hits_coll.create(
0119           id_desc.encode({{"system", 255}, {"x", 0}, {"y", 0}}), // std::uint64_t cellID,
0120           5.0,                                                   // float energy,
0121           0.0,                                                   // float energyError,
0122           0.0,                                                   // float time,
0123           0.0,                                                   // float timeError,
0124           edm4hep::Vector3f(0.0, 0.0, 0.0),                      // edm4hep::Vector3f position,
0125           edm4hep::Vector3f(1.0, 1.0, 0.0),                      // edm4hep::Vector3f dimension,
0126           0,                                                     // std::int32_t sector,
0127           0,                                                     // std::int32_t layer,
0128           edm4hep::Vector3f(0.0, 0.0, 0.0)                       // edm4hep::Vector3f local
0129       );
0130       hits_coll.create(
0131           id_desc.encode({{"system", 255}, {"x", 1}, {"y", 0}}), // std::uint64_t cellID,
0132           6.0,                                                   // float energy,
0133           0.0,                                                   // float energyError,
0134           0.0,                                                   // float time,
0135           0.0,                                                   // float timeError,
0136           edm4hep::Vector3f(0.0, 0.0, 0.0),                      // edm4hep::Vector3f position,
0137           edm4hep::Vector3f(1.0, 1.0, 0.0),                      // edm4hep::Vector3f dimension,
0138           0,                                                     // std::int32_t sector,
0139           0,                                                     // std::int32_t layer,
0140           edm4hep::Vector3f(0.9 /* mm */, 0.9 /* mm */, 0.0)     // edm4hep::Vector3f local
0141       );
0142       auto protoclust_coll = std::make_unique<edm4eic::ProtoClusterCollection>();
0143       algo.process({&hits_coll}, {protoclust_coll.get()});
0144
0145       REQUIRE((*protoclust_coll).size() == 1);
0146       REQUIRE((*protoclust_coll)[0].hits_size() == 2);
0147       REQUIRE((*protoclust_coll)[0].weights_size() == 2);
0148     }
0149   }
0150
0151   SECTION("run on three adjacent cells") {
0152     bool use_adjacencyMatrix = GENERATE(false, true);
0153     if (use_adjacencyMatrix) {
0154       cfg.adjacencyMatrix = "abs(x_1 - x_2) + abs(y_1 - y_2) == 1";
0155     } else {
0156       cfg.localDistXY = {1 * dd4hep::mm, 1 * dd4hep::mm};
0157     }
0158     bool disalow_diagonal_peaks = GENERATE(false, true);
0159     if (disalow_diagonal_peaks) {
0160       cfg.peakNeighbourhoodMatrix = "max(abs(x_1 - x_2), abs(y_1 - y_2)) == 1";
0161     } else {
0162       cfg.peakNeighbourhoodMatrix = "abs(x_1 - x_2) + abs(y_1 - y_2) == 1";
0163     }
0164     cfg.readout = "MockCalorimeterHits";
0165
0166     cfg.splitCluster = GENERATE(false, true);
0167     if (cfg.splitCluster) {
0168       cfg.transverseEnergyProfileMetric = "localDistXY";
0169       cfg.transverseEnergyProfileScale =
0170           std::numeric_limits<decltype(cfg.transverseEnergyProfileScale)>::infinity();
0171     }
0172     cfg.localDistXY = {1 * dd4hep::mm, 1 * dd4hep::mm};
0173     algo.applyConfig(cfg);
0174     algo.init();
0175
0176     edm4eic::CalorimeterHitCollection hits_coll;
0177     hits_coll.create(id_desc.encode({{"system", 255}, {"x", 0}, {"y", 0}}), // std::uint64_t cellID,
0178                      5.0,                                                   // float energy,
0179                      0.0,                                                   // float energyError,
0180                      0.0,                                                   // float time,
0181                      0.0,                                                   // float timeError,
0182                      edm4hep::Vector3f(0.0, 0.0, 0.0), // edm4hep::Vector3f position,
0183                      edm4hep::Vector3f(1.0, 1.0, 0.0), // edm4hep::Vector3f dimension,
0184                      0,                                // std::int32_t sector,
0185                      0,                                // std::int32_t layer,
0186                      edm4hep::Vector3f(0.0, 0.0, 0.0)  // edm4hep::Vector3f local
0187     );
0188     hits_coll.create(id_desc.encode({{"system", 255}, {"x", 1}, {"y", 0}}), // std::uint64_t cellID,
0189                      1.0,                                                   // float energy,
0190                      0.0,                                                   // float energyError,
0191                      0.0,                                                   // float time,
0192                      0.0,                                                   // float timeError,
0193                      edm4hep::Vector3f(0.0, 0.0, 0.0), // edm4hep::Vector3f position,
0194                      edm4hep::Vector3f(1.0, 1.0, 0.0), // edm4hep::Vector3f dimension,
0195                      0,                                // std::int32_t sector,
0196                      0,                                // std::int32_t layer,
0197                      edm4hep::Vector3f(0.9 /* mm */, 0.9 /* mm */, 0.0) // edm4hep::Vector3f local
0198     );
0199
0200     bool test_diagonal_cluster = GENERATE(false, true);
0201     // If false, test a cluster with two maxima:
0202     //  XxX
0203     // If true, test a diagonal cluster:
0204     //  Xx
0205     //   X
0206     // The idea is to test whether peakNeighbourhoodMatrix allows increasing
0207     // peak resolution threshold while keeping the Island algorithm the same.
0208     hits_coll.create(
0209         id_desc.encode({{"system", 255},
0210                         {"x", test_diagonal_cluster ? 1 : 2},
0211                         {"y", test_diagonal_cluster ? 1 : 0}}), // std::uint64_t cellID,
0212         6.0,                                                    // float energy,
0213         0.0,                                                    // float energyError,
0214         0.0,                                                    // float time,
0215         0.0,                                                    // float timeError,
0216         edm4hep::Vector3f(0.0, 0.0, 0.0),                       // edm4hep::Vector3f position,
0217         edm4hep::Vector3f(1.0, 1.0, 0.0),                       // edm4hep::Vector3f dimension,
0218         0,                                                      // std::int32_t sector,
0219         0,                                                      // std::int32_t layer,
0220         edm4hep::Vector3f(1.8 /* mm */, 1.8 /* mm */, 0.0)      // edm4hep::Vector3f local
0221     );
0222     auto protoclust_coll = std::make_unique<edm4eic::ProtoClusterCollection>();
0223     algo.process({&hits_coll}, {protoclust_coll.get()});
0224
0225     bool expect_two_peaks = cfg.splitCluster;
0226     if (cfg.splitCluster && disalow_diagonal_peaks) {
0227       expect_two_peaks = not test_diagonal_cluster;
0228     }
0229     if (expect_two_peaks) {
0230       REQUIRE((*protoclust_coll).size() == 2);
0231       REQUIRE((*protoclust_coll)[0].hits_size() == 3);
0232       REQUIRE((*protoclust_coll)[0].weights_size() == 3);
0233       for (double weight : (*protoclust_coll)[0].getWeights()) {
0234         double energy_fraction =
0235             hits_coll[0].getEnergy() / (hits_coll[0].getEnergy() + hits_coll[2].getEnergy());
0236         REQUIRE_THAT(weight, Catch::Matchers::WithinAbs(energy_fraction, 1e-5));
0237       }
0238       REQUIRE((*protoclust_coll)[1].hits_size() == 3);
0239       REQUIRE((*protoclust_coll)[1].weights_size() == 3);
0240       for (double weight : (*protoclust_coll)[1].getWeights()) {
0241         double energy_fraction =
0242             hits_coll[2].getEnergy() / (hits_coll[0].getEnergy() + hits_coll[2].getEnergy());
0243         REQUIRE_THAT(weight, Catch::Matchers::WithinAbs(energy_fraction, 1e-5));
0244       }
0245     } else {
0246       REQUIRE((*protoclust_coll).size() == 1);
0247       REQUIRE((*protoclust_coll)[0].hits_size() == 3);
0248       REQUIRE((*protoclust_coll)[0].weights_size() == 3);
0249     }
0250   }
0251 }