algorithms/digi/SiliconChargeSharing.cc

0001 // SPDX-License-Identifier: LGPL-3.0-or-later
0002 // Copyright (C) 2024-2025 Chun Yuen Tsang, Prithwish Tribedy, Simon Gardner
0003 //
0004 // Spread energy deposition from one strip to neighboring strips within sensor boundaries
0005
0006 #include <DD4hep/Alignments.h>
0007 #include <DD4hep/DetElement.h>
0008 #include <DD4hep/Handle.h>
0009 #include <DD4hep/Objects.h>
0010 #include <DD4hep/Readout.h>
0011 #include <DD4hep/Segmentations.h>
0012 #include <DD4hep/Shapes.h>
0013 #include <DD4hep/VolumeManager.h>
0014 #include <DD4hep/detail/SegmentationsInterna.h>
0015 #include <DDSegmentation/MultiSegmentation.h>
0016 #include <DDSegmentation/Segmentation.h>
0017 #include <Evaluator/DD4hepUnits.h>
0018 #include <Math/GenVector/Cartesian3D.h>
0019 #include <Math/GenVector/DisplacementVector3D.h>
0020 #include <RtypesCore.h>
0021 #include <TGeoBBox.h>
0022 #include <TGeoMatrix.h>
0023 #include <algorithms/geo.h>
0024 #include <edm4hep/Vector3d.h>
0025 #include <fmt/core.h>
0026 #include <cmath>
0027 #include <gsl/pointers>
0028 #include <numbers>
0029 #include <set>
0030 #include <stdexcept>
0031 #include <typeinfo>
0032 #include <utility>
0033
0034 #include "DD4hep/Detector.h"
0035 #include "SiliconChargeSharing.h"
0036 #include "algorithms/digi/SiliconChargeSharingConfig.h"
0037
0038 namespace eicrecon {
0039
0040 void SiliconChargeSharing::init() {
0041   m_converter = algorithms::GeoSvc::instance().cellIDPositionConverter();
0042   m_seg       = algorithms::GeoSvc::instance().detector()->readout(m_cfg.readout).segmentation();
0043 }
0044
0045 void SiliconChargeSharing::process(const SiliconChargeSharing::Input& input,
0046                                    const SiliconChargeSharing::Output& output) const {
0047   const auto [simhits] = input;
0048   auto [sharedHits]    = output;
0049
0050   for (const auto& hit : *simhits) {
0051
0052     auto cellID = hit.getCellID();
0053
0054     const auto* element = &m_converter->findContext(cellID)->element; // volume context
0055     // ToDo: Move this to init() and set it once for every detelement associated with the readout
0056     // Set transformation map if it hasn't already been set
0057     auto [transformIt, transformInserted] =
0058         m_transform_map.try_emplace(element, &element->nominal().worldTransformation());
0059     auto [segmentationIt, segmentationInserted] =
0060         m_segmentation_map.try_emplace(element, getLocalSegmentation(cellID));
0061
0062     // Try and get a box of the detectorElement solid requiring segmentation
0063     // to be a CartesianGridXY, throwing exception in getLocalSegmentation
0064     try {
0065       dd4hep::Box box = element->solid();
0066       m_xy_range_map.try_emplace(element, box->GetDX(), box->GetDY());
0067     } catch (const std::bad_cast& e) {
0068       error("Failed to cast solid to Box: {}", e.what());
0069     }
0070
0071     auto edep         = hit.getEDep();
0072     auto globalHitPos = hit.getPosition();
0073     auto hitPos =
0074         global2Local(dd4hep::Position(globalHitPos.x * dd4hep::mm, globalHitPos.y * dd4hep::mm,
0075                                       globalHitPos.z * dd4hep::mm),
0076                      transformIt->second);
0077
0078     // therefore, we search neighbors within the segmentation of the same volume
0079     // to find the cell ID that correspond to globalHitPos.
0080     // Precise reason unknown, but we suspect it's cause by steps in Geant4
0081     // Perhaps position is the average of all steps in volume while cellID is just the first cell the track hits
0082     // They disagree when there are multiple step and scattering inside the volume
0083     const dd4hep::Position dummy;
0084     cellID = segmentationIt->second->cellID(hitPos, dummy, cellID);
0085
0086     std::unordered_set<dd4hep::rec::CellID> tested_cells;
0087     std::unordered_map<dd4hep::rec::CellID, float> cell_charge;
0088
0089     // Warning: This function is recursive, it stops shen it finds the edge of a detector element
0090     // or when the energy deposited in a cell is below the configured threshold
0091     findAllNeighborsInSensor(cellID, tested_cells, edep, hitPos, segmentationIt->second,
0092                              m_xy_range_map[element], hit, sharedHits);
0093
0094   } // for simhits
0095 } // SiliconChargeSharing:process
0096
0097 // Recursively find neighbors where a charge is deposited
0098 void SiliconChargeSharing::findAllNeighborsInSensor(
0099     const dd4hep::rec::CellID testCellID, std::unordered_set<dd4hep::rec::CellID>& tested_cells,
0100     const float edep, const dd4hep::Position hitPos,
0101     const dd4hep::DDSegmentation::CartesianGridXY* segmentation,
0102     const std::pair<double, double>& xy_range, const edm4hep::SimTrackerHit& hit,
0103     edm4hep::SimTrackerHitCollection* sharedHits) const {
0104
0105   // Tag cell as tested
0106   tested_cells.insert(testCellID);
0107
0108   auto localPos = cell2LocalPosition(testCellID);
0109
0110   // Check if the cell is within the sensor boundaries
0111   if (std::abs(localPos.x()) > xy_range.first || std::abs(localPos.y()) > xy_range.second) {
0112     return;
0113   }
0114
0115   // cout the local position and hit position
0116   auto xDimension = segmentation->gridSizeX();
0117   auto yDimension = segmentation->gridSizeY();
0118   // Calculate deposited energy in cell
0119   float edepCell = energyAtCell(xDimension, yDimension, localPos, hitPos, edep);
0120   if (edepCell <= m_cfg.min_edep) {
0121     return;
0122   }
0123
0124   // Create a new simhit for cell with deposited energy
0125   auto globalCellPos = m_converter->position(testCellID);
0126
0127   edm4hep::MutableSimTrackerHit shared_hit = hit.clone();
0128   shared_hit.setCellID(testCellID);
0129   shared_hit.setEDep(edepCell);
0130   shared_hit.setPosition({globalCellPos.x() / dd4hep::mm, globalCellPos.y() / dd4hep::mm,
0131                           globalCellPos.z() / dd4hep::mm});
0132   sharedHits->push_back(shared_hit);
0133
0134   // As there is charge in the cell, test the neighbors too
0135   std::set<dd4hep::rec::CellID> testCellNeighbours;
0136   segmentation->neighbours(testCellID, testCellNeighbours);
0137
0138   for (const auto& neighbourCell : testCellNeighbours) {
0139     if (tested_cells.find(neighbourCell) == tested_cells.end()) {
0140       findAllNeighborsInSensor(neighbourCell, tested_cells, edep, hitPos, segmentation, xy_range,
0141                                hit, sharedHits);
0142     }
0143   }
0144 }
0145
0146 // Calculate integral of Gaussian distribution
0147 float SiliconChargeSharing::integralGaus(float mean, float sd, float low_lim, float up_lim) {
0148   // return integral Gauss(mean, sd) dx from x = low_lim to x = up_lim
0149   // default value is set when sd = 0
0150   float up  = mean > up_lim ? -0.5 : 0.5;
0151   float low = mean > low_lim ? -0.5 : 0.5;
0152   if (sd > 0) {
0153     up  = -0.5 * std::erf(std::numbers::sqrt2 * (mean - up_lim) / sd);
0154     low = -0.5 * std::erf(std::numbers::sqrt2 * (mean - low_lim) / sd);
0155   }
0156   return up - low;
0157 }
0158
0159 // Convert cellID to local position
0160 dd4hep::Position SiliconChargeSharing::cell2LocalPosition(const dd4hep::rec::CellID& cell) const {
0161   auto position = m_seg->position(cell); // local position
0162   return position;
0163 }
0164
0165 // Convert global position to local position
0166 dd4hep::Position SiliconChargeSharing::global2Local(const dd4hep::Position& globalPosition,
0167                                                     const TGeoHMatrix* transform) {
0168
0169   double g[3];
0170   double l[3];
0171
0172   globalPosition.GetCoordinates(static_cast<Double_t*>(g));
0173   transform->MasterToLocal(static_cast<const Double_t*>(g), static_cast<Double_t*>(l));
0174   dd4hep::Position localPosition;
0175   localPosition.SetCoordinates(static_cast<const Double_t*>(l));
0176   return localPosition;
0177 }
0178
0179 // Calculate energy deposition in a cell relative to the hit position
0180 float SiliconChargeSharing::energyAtCell(const double xDimension, const double yDimension,
0181                                          const dd4hep::Position localPos,
0182                                          const dd4hep::Position hitPos, const float edep) const {
0183   auto sigma_sharingx = m_cfg.sigma_sharingx;
0184   auto sigma_sharingy = m_cfg.sigma_sharingy;
0185   if (m_cfg.sigma_mode == SiliconChargeSharingConfig::ESigmaMode::rel) {
0186     sigma_sharingx *= xDimension;
0187     sigma_sharingy *= yDimension;
0188   }
0189   float energy = edep *
0190                  integralGaus(hitPos.x(), sigma_sharingx, localPos.x() - 0.5 * xDimension,
0191                               localPos.x() + 0.5 * xDimension) *
0192                  integralGaus(hitPos.y(), sigma_sharingy, localPos.y() - 0.5 * yDimension,
0193                               localPos.y() + 0.5 * yDimension);
0194   return energy;
0195 }
0196
0197 // Get the segmentation relevant to a cellID
0198 const dd4hep::DDSegmentation::CartesianGridXY*
0199 SiliconChargeSharing::getLocalSegmentation(const dd4hep::rec::CellID& cellID) const {
0200   // Get the segmentation type
0201   auto segmentation_type                                   = m_seg.type();
0202   const dd4hep::DDSegmentation::Segmentation* segmentation = m_seg.segmentation();
0203   // Check if the segmentation is a multi-segmentation
0204   while (segmentation_type == "MultiSegmentation") {
0205     const auto* multi_segmentation =
0206         dynamic_cast<const dd4hep::DDSegmentation::MultiSegmentation*>(segmentation);
0207     segmentation      = &multi_segmentation->subsegmentation(cellID);
0208     segmentation_type = segmentation->type();
0209   }
0210
0211   // Try to cast the segmentation to CartesianGridXY
0212   const auto* cartesianGrid =
0213       dynamic_cast<const dd4hep::DDSegmentation::CartesianGridXY*>(segmentation);
0214   if (cartesianGrid == nullptr) {
0215     throw std::runtime_error("Segmentation is not of type CartesianGridXY");
0216   }
0217
0218   return cartesianGrid;
0219 }
0220
0221 } // namespace eicrecon