Back to home page

EIC code displayed by LXR
 
 

    

File indexing completed on 2025-07-06 07:55:53 

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 =
0036       spdlog::default_logger()->clone("CalorimeterIslandCluster");
0037   logger->set_level(spdlog::level::trace);
0038 
0039   CalorimeterIslandClusterConfig cfg;
0040   cfg.minClusterHitEdep    = 0. * dd4hep::GeV;
0041   cfg.minClusterCenterEdep = 0. * dd4hep::GeV;
0042 
0043   auto detector = algorithms::GeoSvc::instance().detector();
0044   auto id_desc  = detector->readout("MockCalorimeterHits").idSpec();
0045 
0046   SECTION("without splitting") {
0047     bool use_adjacencyMatrix = GENERATE(false, true);
0048     cfg.splitCluster         = false;
0049     if (use_adjacencyMatrix) {
0050       cfg.adjacencyMatrix = "abs(x_1 - x_2) + abs(y_1 - y_2) == 1";
0051       cfg.readout         = "MockCalorimeterHits";
0052     } else {
0053       cfg.localDistXY = {1 * dd4hep::mm, 1 * dd4hep::mm};
0054     }
0055     algo.applyConfig(cfg);
0056     algo.init();
0057 
0058     SECTION("on a single cell") {
0059       edm4eic::CalorimeterHitCollection hits_coll;
0060       hits_coll.create(
0061           id_desc.encode({{"system", 255}, {"x", 0}, {"y", 0}}), // std::uint64_t cellID,
0062           5.0,                                                   // float energy,
0063           0.0,                                                   // float energyError,
0064           0.0,                                                   // float time,
0065           0.0,                                                   // float timeError,
0066           edm4hep::Vector3f(0.0, 0.0, 0.0),                      // edm4hep::Vector3f position,
0067           edm4hep::Vector3f(0.0, 0.0, 0.0),                      // edm4hep::Vector3f dimension,
0068           0,                                                     // std::int32_t sector,
0069           0,                                                     // std::int32_t layer,
0070           edm4hep::Vector3f(0.0, 0.0, 0.0)                       // edm4hep::Vector3f local
0071       );
0072       auto protoclust_coll = std::make_unique<edm4eic::ProtoClusterCollection>();
0073       algo.process({&hits_coll}, {protoclust_coll.get()});
0074 
0075       REQUIRE((*protoclust_coll).size() == 1);
0076       REQUIRE((*protoclust_coll)[0].hits_size() == 1);
0077       REQUIRE((*protoclust_coll)[0].weights_size() == 1);
0078     }
0079 
0080     SECTION("on two separated cells") {
0081       edm4eic::CalorimeterHitCollection hits_coll;
0082       hits_coll.create(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(1,                                // std::uint64_t cellID,
0094                        6.0,                              // float energy,
0095                        0.0,                              // float energyError,
0096                        0.0,                              // float time,
0097                        0.0,                              // float timeError,
0098                        edm4hep::Vector3f(0.0, 0.0, 0.0), // edm4hep::Vector3f position,
0099                        edm4hep::Vector3f(1.0, 1.0, 0.0), // edm4hep::Vector3f dimension,
0100                        0,                                // std::int32_t sector,
0101                        0,                                // std::int32_t layer,
0102                        edm4hep::Vector3f(1.1 /* mm */, 1.1 /* mm */, 0.0) // edm4hep::Vector3f local
0103       );
0104       auto protoclust_coll = std::make_unique<edm4eic::ProtoClusterCollection>();
0105       algo.process({&hits_coll}, {protoclust_coll.get()});
0106 
0107       REQUIRE((*protoclust_coll).size() == 2);
0108       REQUIRE((*protoclust_coll)[0].hits_size() == 1);
0109       REQUIRE((*protoclust_coll)[0].weights_size() == 1);
0110       REQUIRE((*protoclust_coll)[1].hits_size() == 1);
0111       REQUIRE((*protoclust_coll)[1].weights_size() == 1);
0112     }
0113 
0114     SECTION("on two adjacent cells") {
0115       edm4eic::CalorimeterHitCollection hits_coll;
0116       hits_coll.create(
0117           id_desc.encode({{"system", 255}, {"x", 0}, {"y", 0}}), // std::uint64_t cellID,
0118           5.0,                                                   // float energy,
0119           0.0,                                                   // float energyError,
0120           0.0,                                                   // float time,
0121           0.0,                                                   // float timeError,
0122           edm4hep::Vector3f(0.0, 0.0, 0.0),                      // edm4hep::Vector3f position,
0123           edm4hep::Vector3f(1.0, 1.0, 0.0),                      // edm4hep::Vector3f dimension,
0124           0,                                                     // std::int32_t sector,
0125           0,                                                     // std::int32_t layer,
0126           edm4hep::Vector3f(0.0, 0.0, 0.0)                       // edm4hep::Vector3f local
0127       );
0128       hits_coll.create(
0129           id_desc.encode({{"system", 255}, {"x", 1}, {"y", 0}}), // std::uint64_t cellID,
0130           6.0,                                                   // float energy,
0131           0.0,                                                   // float energyError,
0132           0.0,                                                   // float time,
0133           0.0,                                                   // float timeError,
0134           edm4hep::Vector3f(0.0, 0.0, 0.0),                      // edm4hep::Vector3f position,
0135           edm4hep::Vector3f(1.0, 1.0, 0.0),                      // edm4hep::Vector3f dimension,
0136           0,                                                     // std::int32_t sector,
0137           0,                                                     // std::int32_t layer,
0138           edm4hep::Vector3f(0.9 /* mm */, 0.9 /* mm */, 0.0)     // edm4hep::Vector3f local
0139       );
0140       auto protoclust_coll = std::make_unique<edm4eic::ProtoClusterCollection>();
0141       algo.process({&hits_coll}, {protoclust_coll.get()});
0142 
0143       REQUIRE((*protoclust_coll).size() == 1);
0144       REQUIRE((*protoclust_coll)[0].hits_size() == 2);
0145       REQUIRE((*protoclust_coll)[0].weights_size() == 2);
0146     }
0147   }
0148 
0149   SECTION("run on three adjacent cells") {
0150     bool use_adjacencyMatrix = GENERATE(false, true);
0151     if (use_adjacencyMatrix) {
0152       cfg.adjacencyMatrix = "abs(x_1 - x_2) + abs(y_1 - y_2) == 1";
0153     } else {
0154       cfg.localDistXY = {1 * dd4hep::mm, 1 * dd4hep::mm};
0155     }
0156     bool disalow_diagonal_peaks = GENERATE(false, true);
0157     if (disalow_diagonal_peaks) {
0158       cfg.peakNeighbourhoodMatrix = "max(abs(x_1 - x_2), abs(y_1 - y_2)) == 1";
0159     } else {
0160       cfg.peakNeighbourhoodMatrix = "abs(x_1 - x_2) + abs(y_1 - y_2) == 1";
0161     }
0162     cfg.readout = "MockCalorimeterHits";
0163 
0164     cfg.splitCluster = GENERATE(false, true);
0165     if (cfg.splitCluster) {
0166       cfg.transverseEnergyProfileMetric = "localDistXY";
0167       cfg.transverseEnergyProfileScale =
0168           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(id_desc.encode({{"system", 255}, {"x", 0}, {"y", 0}}), // std::uint64_t cellID,
0176                      5.0,                                                   // float energy,
0177                      0.0,                                                   // float energyError,
0178                      0.0,                                                   // float time,
0179                      0.0,                                                   // float timeError,
0180                      edm4hep::Vector3f(0.0, 0.0, 0.0), // edm4hep::Vector3f position,
0181                      edm4hep::Vector3f(1.0, 1.0, 0.0), // edm4hep::Vector3f dimension,
0182                      0,                                // std::int32_t sector,
0183                      0,                                // std::int32_t layer,
0184                      edm4hep::Vector3f(0.0, 0.0, 0.0)  // edm4hep::Vector3f local
0185     );
0186     hits_coll.create(id_desc.encode({{"system", 255}, {"x", 1}, {"y", 0}}), // std::uint64_t cellID,
0187                      1.0,                                                   // float energy,
0188                      0.0,                                                   // float energyError,
0189                      0.0,                                                   // float time,
0190                      0.0,                                                   // float timeError,
0191                      edm4hep::Vector3f(0.0, 0.0, 0.0), // edm4hep::Vector3f position,
0192                      edm4hep::Vector3f(1.0, 1.0, 0.0), // edm4hep::Vector3f dimension,
0193                      0,                                // std::int32_t sector,
0194                      0,                                // std::int32_t layer,
0195                      edm4hep::Vector3f(0.9 /* mm */, 0.9 /* mm */, 0.0) // edm4hep::Vector3f local
0196     );
0197 
0198     bool test_diagonal_cluster = GENERATE(false, true);
0199     // If false, test a cluster with two maxima:
0200     //  XxX
0201     // If true, test a diagonal cluster:
0202     //  Xx
0203     //   X
0204     // The idea is to test whether peakNeighbourhoodMatrix allows increasing
0205     // peak resolution threshold while keeping the Island algorithm the same.
0206     hits_coll.create(
0207         id_desc.encode({{"system", 255},
0208                         {"x", test_diagonal_cluster ? 1 : 2},
0209                         {"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 =
0233             hits_coll[0].getEnergy() / (hits_coll[0].getEnergy() + hits_coll[2].getEnergy());
0234         REQUIRE_THAT(weight, Catch::Matchers::WithinAbs(energy_fraction, 1e-5));
0235       }
0236       REQUIRE((*protoclust_coll)[1].hits_size() == 3);
0237       REQUIRE((*protoclust_coll)[1].weights_size() == 3);
0238       for (double weight : (*protoclust_coll)[1].getWeights()) {
0239         double energy_fraction =
0240             hits_coll[2].getEnergy() / (hits_coll[0].getEnergy() + hits_coll[2].getEnergy());
0241         REQUIRE_THAT(weight, Catch::Matchers::WithinAbs(energy_fraction, 1e-5));
0242       }
0243     } else {
0244       REQUIRE((*protoclust_coll).size() == 1);
0245       REQUIRE((*protoclust_coll)[0].hits_size() == 3);
0246       REQUIRE((*protoclust_coll)[0].weights_size() == 3);
0247     }
0248   }
0249 }
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated by the 2.3.7 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!