tests/algorithms_test/calorimetry_HEXPLIT.cc

0001 // SPDX-License-Identifier: LGPL-3.0-or-later
0002 // Copyright (C) 2024, Sebouh Paul
0003
0004 #include <DD4hep/Detector.h>       // for Detector
0005 #include <DD4hep/IDDescriptor.h>   // for IDDescriptor
0006 #include <DD4hep/Readout.h>        // for Readout
0007 #include <Evaluator/DD4hepUnits.h> // for MeV, mm, keV, ns
0008 #include <algorithms/geo.h>
0009 #include <catch2/catch_test_macros.hpp> // for AssertionHandler, operator""_catch_sr, StringRef, REQUIRE, operator<, operator==, operator>, TEST_CASE
0010 #include <edm4eic/CalorimeterHitCollection.h> // for CalorimeterHitCollection, MutableCalorimeterHit, CalorimeterHitMutableCollectionIterator
0011 #include <edm4hep/Vector3f.h> // for Vector3f
0012 #include <spdlog/common.h>    // for level_enum
0013 #include <spdlog/logger.h>    // for logger
0014 #include <spdlog/spdlog.h>    // for default_logger
0015 #include <array>              // for array
0016 #include <cmath>              // for sqrt, abs
0017 #include <cstddef>
0018 #include <gsl/pointers>
0019 #include <memory> // for allocator, unique_ptr, make_unique, shared_ptr, __shared_ptr_access
0020 #include <numbers>
0021 #include <string>
0022 #include <utility> // for pair
0023 #include <vector>
0024
0025 #include "algorithms/calorimetry/HEXPLIT.h"       // for HEXPLIT
0026 #include "algorithms/calorimetry/HEXPLITConfig.h" // for HEXPLITConfig
0027
0028 using eicrecon::HEXPLIT;
0029 using eicrecon::HEXPLITConfig;
0030
0031 TEST_CASE("the subcell-splitting algorithm runs", "[HEXPLIT]") {
0032   HEXPLIT algo("HEXPLIT");
0033
0034   std::shared_ptr<spdlog::logger> logger = spdlog::default_logger()->clone("HEXPLIT");
0035   logger->set_level(spdlog::level::trace);
0036
0037   HEXPLITConfig cfg;
0038   cfg.MIP  = 472. * dd4hep::keV;
0039   cfg.tmax = 1000. * dd4hep::ns;
0040
0041   auto detector = algorithms::GeoSvc::instance().detector();
0042   auto id_desc  = detector->readout("MockCalorimeterHits").idSpec();
0043
0044   //create a geometry for the fake detector.
0045   double side_length   = 31.3 * dd4hep::mm;
0046   double layer_spacing = 25.1 * dd4hep::mm;
0047   double thickness     = 3 * dd4hep::mm;
0048
0049   //dimension of a cell
0050   auto dimension = edm4hep::Vector3f(2 * side_length, std::numbers::sqrt3 * side_length, thickness);
0051
0052   algo.applyConfig(cfg);
0053   algo.init();
0054
0055   edm4eic::CalorimeterHitCollection hits_coll;
0056
0057   //create a set of 5 hits in consecutive layers, all of which overlap in a single rhombus,
0058   // centered at (3/8, sqrt(3)/8)*side_length
0059   std::array<double, 5> layer = {0, 1, 2, 3, 4};
0060   std::array<double, 5> x     = {0, 0.75 * side_length, 0, 0.75 * side_length, 0};
0061   std::array<double, 5> y     = {
0062       std::numbers::sqrt3 / 2 * side_length, -0.25 * std::numbers::sqrt3 * side_length, 0,
0063       0.25 * std::numbers::sqrt3 * side_length, std::numbers::sqrt3 / 2 * side_length};
0064   std::array<double, 5> E = {50 * dd4hep::MeV, 50 * dd4hep::MeV, 50 * dd4hep::MeV, 50 * dd4hep::MeV,
0065                              50 * dd4hep::MeV};
0066   for (std::size_t i = 0; i < 5; i++) {
0067     hits_coll.create(
0068         id_desc.encode({{"system", 255}, {"x", 0}, {"y", 0}}),   // std::uint64_t cellID,
0069         E[i],                                                    // float energy,
0070         0.0,                                                     // float energyError,
0071         0.0,                                                     // float time,
0072         0.0,                                                     // float timeError,
0073         edm4hep::Vector3f(x[i], y[i], layer[i] * layer_spacing), // edm4hep::Vector3f position,
0074         dimension,                                               // edm4hep::Vector3f dimension,
0075         0,                                                       // std::int32_t sector,
0076         layer[i],                                                // std::int32_t layer,
0077         edm4hep::Vector3f(x[i], y[i], layer[i] * layer_spacing)  // edm4hep::Vector3f local
0078     );
0079   }
0080
0081   auto subcellhits_coll = std::make_unique<edm4eic::CalorimeterHitCollection>();
0082   algo.process({&hits_coll}, {subcellhits_coll.get()});
0083
0084   //the number of subcell hits should be equal to the
0085   //number of subcells per cell (12) times the number of cells (5)
0086   REQUIRE((*subcellhits_coll).size() == 60);
0087
0088   //next check that the sum of the hit energies equals the energy that I gave the hits
0089   double tol  = 0.001;
0090   double Esum = 0;
0091   int i       = 0;
0092   for (auto subcell : *subcellhits_coll) {
0093     Esum += subcell.getEnergy();
0094     i++;
0095     if (i % 12 == 0) {
0096       REQUIRE(std::abs(Esum - E[i / 12 - 1]) / E[i / 12 - 1] < tol);
0097       Esum = 0;
0098     }
0099   }
0100   // next check that almost all of the energy of the hit in the middle layer
0101   // is in the subcell where the other hits overlap
0102   REQUIRE((*subcellhits_coll)[35].getEnergy() / E[2] > 0.95);
0103 }