algorithms/calorimetry/HEXPLIT.cc

0001 // Copyright (C) 2023 Sebouh Paul
0002 // SPDX-License-Identifier: LGPL-3.0-or-later
0003
0004 // References:
0005 //   https://arxiv.org/abs/2308.06939
0006
0007 #include <DD4hep/Alignments.h>
0008 #include <DD4hep/DetElement.h>
0009 #include <DD4hep/Objects.h>
0010 #include <DD4hep/VolumeManager.h>
0011 #include <Evaluator/DD4hepUnits.h>
0012 #include <Math/GenVector/Cartesian3D.h>
0013 #include <Math/GenVector/DisplacementVector3D.h>
0014 #include <edm4hep/Vector3f.h>
0015 #include <cstdlib>
0016 #include <algorithm>
0017 #include <cmath>
0018 #include <gsl/pointers> // for not_null
0019 #include <vector>
0020
0021 #include "HEXPLIT.h"
0022 #include "algorithms/calorimetry/HEXPLITConfig.h"
0023
0024 namespace eicrecon {
0025
0026 //positions where the overlapping cells are relative to a given cell (in units of hexagon side length)
0027 const std::vector<double> HEXPLIT::neighbor_offsets_x = []() {
0028   std::vector<double> x;
0029   double rs[2]      = {1.5, sqrt(3) / 2.};
0030   double offsets[2] = {0, M_PI / 2};
0031   for (int i = 0; i < 2; i++) {
0032     for (int j = 0; j < 6; j += 1) {
0033       x.push_back(rs[i] * cos(j * M_PI / 3 + offsets[i]));
0034     }
0035   }
0036   return x;
0037 }();
0038
0039 const std::vector<double> HEXPLIT::neighbor_offsets_y = []() {
0040   std::vector<double> y;
0041   double rs[2]      = {1.5, sqrt(3) / 2.};
0042   double offsets[2] = {0, M_PI / 2};
0043   for (int i = 0; i < 2; i++) {
0044     for (int j = 0; j < 6; j += 1) {
0045       y.push_back(rs[i] * sin(j * M_PI / 3 + offsets[i]));
0046     }
0047   }
0048   return y;
0049 }();
0050
0051 //indices of the neighboring cells which overlap to produce a given subcell
0052 const int HEXPLIT::neighbor_indices[SUBCELLS][OVERLAP] = {
0053     {0, 11, 10}, {1, 6, 11}, {2, 7, 6}, {3, 8, 7},  {4, 9, 8},   {5, 10, 9},
0054     {6, 11, 7},  {7, 6, 8},  {8, 7, 9}, {9, 8, 10}, {10, 9, 11}, {11, 10, 6}};
0055
0056 //positions of the centers of subcells
0057 const std::vector<double> HEXPLIT::subcell_offsets_x = []() {
0058   std::vector<double> x;
0059   double rs[2]      = {0.75, sqrt(3) / 4.};
0060   double offsets[2] = {0, M_PI / 2};
0061   for (int i = 0; i < 2; i++) {
0062     for (int j = 0; j < 6; j += 1) {
0063       x.push_back(rs[i] * cos(j * M_PI / 3 + offsets[i]));
0064     }
0065   }
0066   return x;
0067 }();
0068
0069 const std::vector<double> HEXPLIT::subcell_offsets_y = []() {
0070   std::vector<double> y;
0071   double rs[2]      = {0.75, sqrt(3) / 4.};
0072   double offsets[2] = {0, M_PI / 2};
0073   for (int i = 0; i < 2; i++) {
0074     for (int j = 0; j < 6; j += 1) {
0075       y.push_back(rs[i] * sin(j * M_PI / 3 + offsets[i]));
0076     }
0077   }
0078   return y;
0079 }();
0080
0081 void HEXPLIT::init() {}
0082
0083 void HEXPLIT::process(const HEXPLIT::Input& input, const HEXPLIT::Output& output) const {
0084
0085   const auto [hits]  = input;
0086   auto [subcellHits] = output;
0087
0088   double MIP   = m_cfg.MIP / dd4hep::GeV;
0089   double delta = m_cfg.delta_in_MIPs * MIP;
0090   double Emin  = m_cfg.Emin_in_MIPs * MIP;
0091   double tmax  = m_cfg.tmax / dd4hep::ns;
0092
0093   auto volman = m_detector->volumeManager();
0094
0095   for (const auto& hit : *hits) {
0096     //skip hits that do not pass E and t cuts
0097     if (hit.getEnergy() < Emin || hit.getTime() > tmax) {
0098       continue;
0099     }
0100
0101     //keep track of the energy in each neighboring cell
0102     std::vector<double> Eneighbors(NEIGHBORS, 0.0);
0103
0104     double sl = hit.getDimension().x / 2.;
0105     for (const auto& other_hit : *hits) {
0106       // maximum distance between where the neighboring cell is and where it should be
0107       // based on an ideal geometry using the staggered tessellation pattern.
0108       // Deviations could arise from rounding errors or from detector misalignment.
0109       double tol = 0.1; // in units of side lengths.
0110
0111       //only look at hits nearby within two layers of the current layer
0112       int dz = std::abs(hit.getLayer() - other_hit.getLayer());
0113       if (dz > 2 || dz == 0) {
0114         continue;
0115       }
0116       if (other_hit.getEnergy() < Emin || other_hit.getTime() > tmax) {
0117         continue;
0118       }
0119       //difference in transverse position (in units of side lengths)
0120       double dx = (other_hit.getLocal().x - hit.getLocal().x) / sl;
0121       double dy = (other_hit.getLocal().y - hit.getLocal().y) / sl;
0122       if (std::abs(dx) > 2 || std::abs(dy) > sqrt(3)) {
0123         continue;
0124       }
0125
0126       //loop over locations of the neighboring cells
0127       //and check if the jth hit matches this location
0128       for (int k = 0; k < NEIGHBORS; k++) {
0129         if (std::abs(dx - neighbor_offsets_x[k]) < tol &&
0130             std::abs(dy - neighbor_offsets_y[k]) < tol) {
0131           Eneighbors[k] += other_hit.getEnergy();
0132           break;
0133         }
0134       }
0135     }
0136     double weights[SUBCELLS];
0137     for (int k = 0; k < NEIGHBORS; k++) {
0138       Eneighbors[k] = std::max(Eneighbors[k], delta);
0139     }
0140     double sum_weights = 0;
0141     for (int k = 0; k < SUBCELLS; k++) {
0142       weights[k] = Eneighbors[neighbor_indices[k][0]] * Eneighbors[neighbor_indices[k][1]] *
0143                    Eneighbors[neighbor_indices[k][2]];
0144       sum_weights += weights[k];
0145     }
0146     for (int k = 0; k < SUBCELLS; k++) {
0147
0148       //create the subcell hits.  First determine their positions in local coordinates.
0149       const decltype(edm4eic::CalorimeterHitData::local) local(
0150           hit.getLocal().x + subcell_offsets_x[k] * sl,
0151           hit.getLocal().y + subcell_offsets_y[k] * sl, hit.getLocal().z);
0152
0153       //convert this to a position object so that the global position can be determined
0154       dd4hep::Position local_position;
0155       local_position.SetX(local.x * dd4hep::mm);
0156       local_position.SetY(local.y * dd4hep::mm);
0157       local_position.SetZ(local.z * dd4hep::mm);
0158
0159       dd4hep::Position global_position;
0160       try {
0161
0162         //also convert this to the detector's global coordinates.  To do: check if this is correct
0163         auto alignment = volman.lookupDetElement(hit.getCellID()).nominal();
0164
0165         global_position = alignment.localToWorld(local_position);
0166
0167       } catch (...) {
0168         // do this to prevent errors when running the test on the mock detector
0169         warning("Cannot find transformation from local to global coordinates.");
0170         global_position = local_position;
0171       }
0172
0173       //convert this from position object to a vector object
0174       const decltype(edm4eic::CalorimeterHitData::position) position = {
0175           static_cast<float>(global_position.X() / dd4hep::mm),
0176           static_cast<float>(global_position.Y() / dd4hep::mm),
0177           static_cast<float>(global_position.Z() / dd4hep::mm)};
0178
0179       //bounding box dimensions depend on the orientation of the rhombus
0180       int orientation = static_cast<int>(k % 3 == 0);
0181       const decltype(edm4eic::CalorimeterHitData::dimension) dimension(
0182           sl * (orientation != 0 ? 1 : 1.5), sl * sqrt(3) / 2. * (orientation != 0 ? 2 : 1),
0183           hit.getDimension()[2]);
0184
0185       subcellHits->create(hit.getCellID(), hit.getEnergy() * weights[k] / sum_weights, 0,
0186                           hit.getTime(), 0, position, dimension, hit.getSector(), hit.getLayer(),
0187                           local);
0188     }
0189   }
0190 }
0191
0192 } // namespace eicrecon