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 <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}}), // std::uint64_t cellID,
0061         5.0, // float energy,
0062         0.0, // float energyError,
0063         0.0, // float time,
0064         0.0, // float timeError,
0065         edm4hep::Vector3f(0.0, 0.0, 0.0), // edm4hep::Vector3f position,
0066         edm4hep::Vector3f(0.0, 0.0, 0.0), // edm4hep::Vector3f dimension,
0067         0, // std::int32_t sector,
0068         0, // std::int32_t layer,
0069         edm4hep::Vector3f(0.0, 0.0, 0.0) // edm4hep::Vector3f local
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, // std::uint64_t cellID,
0083         5.0, // float energy,
0084         0.0, // float energyError,
0085         0.0, // float time,
0086         0.0, // float timeError,
0087         edm4hep::Vector3f(0.0, 0.0, 0.0), // edm4hep::Vector3f position,
0088         edm4hep::Vector3f(1.0, 1.0, 0.0), // edm4hep::Vector3f dimension,
0089         0, // std::int32_t sector,
0090         0, // std::int32_t layer,
0091         edm4hep::Vector3f(0.0, 0.0, 0.0) // edm4hep::Vector3f local
0092       );
0093       hits_coll.create(
0094         1, // std::uint64_t cellID,
0095         6.0, // float energy,
0096         0.0, // float energyError,
0097         0.0, // float time,
0098         0.0, // float timeError,
0099         edm4hep::Vector3f(0.0, 0.0, 0.0), // edm4hep::Vector3f position,
0100         edm4hep::Vector3f(1.0, 1.0, 0.0), // edm4hep::Vector3f dimension,
0101         0, // std::int32_t sector,
0102         0, // std::int32_t layer,
0103         edm4hep::Vector3f(1.1 /* mm */, 1.1 /* mm */, 0.0) // edm4hep::Vector3f local
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}}), // std::uint64_t cellID,
0119         5.0, // float energy,
0120         0.0, // float energyError,
0121         0.0, // float time,
0122         0.0, // float timeError,
0123         edm4hep::Vector3f(0.0, 0.0, 0.0), // edm4hep::Vector3f position,
0124         edm4hep::Vector3f(1.0, 1.0, 0.0), // edm4hep::Vector3f dimension,
0125         0, // std::int32_t sector,
0126         0, // std::int32_t layer,
0127         edm4hep::Vector3f(0.0, 0.0, 0.0) // edm4hep::Vector3f local
0128       );
0129       hits_coll.create(
0130         id_desc.encode({{"system", 255}, {"x", 1}, {"y", 0}}), // std::uint64_t cellID,
0131         6.0, // float energy,
0132         0.0, // float energyError,
0133         0.0, // float time,
0134         0.0, // float timeError,
0135         edm4hep::Vector3f(0.0, 0.0, 0.0), // edm4hep::Vector3f position,
0136         edm4hep::Vector3f(1.0, 1.0, 0.0), // edm4hep::Vector3f dimension,
0137         0, // std::int32_t sector,
0138         0, // std::int32_t layer,
0139         edm4hep::Vector3f(0.9 /* mm */, 0.9 /* mm */, 0.0) // edm4hep::Vector3f local
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}}), // std::uint64_t cellID,
0177       5.0, // float energy,
0178       0.0, // float energyError,
0179       0.0, // float time,
0180       0.0, // float timeError,
0181       edm4hep::Vector3f(0.0, 0.0, 0.0), // edm4hep::Vector3f position,
0182       edm4hep::Vector3f(1.0, 1.0, 0.0), // edm4hep::Vector3f dimension,
0183       0, // std::int32_t sector,
0184       0, // std::int32_t layer,
0185       edm4hep::Vector3f(0.0, 0.0, 0.0) // edm4hep::Vector3f local
0186     );
0187     hits_coll.create(
0188       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}, {"x", test_diagonal_cluster ? 1 : 2}, {"y", test_diagonal_cluster ? 1 : 0}}), // std::uint64_t cellID,
0210       6.0, // float energy,
0211       0.0, // float energyError,
0212       0.0, // float time,
0213       0.0, // float timeError,
0214       edm4hep::Vector3f(0.0, 0.0, 0.0), // edm4hep::Vector3f position,
0215       edm4hep::Vector3f(1.0, 1.0, 0.0), // edm4hep::Vector3f dimension,
0216       0, // std::int32_t sector,
0217       0, // std::int32_t layer,
0218       edm4hep::Vector3f(1.8 /* mm */, 1.8 /* mm */, 0.0) // edm4hep::Vector3f local
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 }