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