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