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 <set>
0029 #include <stdexcept>
0030 #include <typeinfo>
0031 #include <utility>
0032
0033 #include "DD4hep/Detector.h"
0034 #include "SiliconChargeSharing.h"
0035 #include "algorithms/digi/SiliconChargeSharingConfig.h"
0036
0037 namespace eicrecon {
0038
0039 void SiliconChargeSharing::init() {
0040   m_converter = algorithms::GeoSvc::instance().cellIDPositionConverter();
0041   m_seg       = algorithms::GeoSvc::instance().detector()->readout(m_cfg.readout).segmentation();
0042 }
0043
0044 void SiliconChargeSharing::process(const SiliconChargeSharing::Input& input,
0045                                    const SiliconChargeSharing::Output& output) const {
0046   const auto [simhits] = input;
0047   auto [sharedHits]    = output;
0048
0049   for (const auto& hit : *simhits) {
0050
0051     auto cellID = hit.getCellID();
0052
0053     const auto* element = &m_converter->findContext(cellID)->element; // volume context
0054     // ToDo: Move this to init() and set it once for every detelement associated with the readout
0055     // Set transformation map if it hasn't already been set
0056     auto [transformIt, transformInserted] =
0057         m_transform_map.try_emplace(element, &element->nominal().worldTransformation());
0058     auto [segmentationIt, segmentationInserted] =
0059         m_segmentation_map.try_emplace(element, getLocalSegmentation(cellID));
0060
0061     // Try and get a box of the detectorElement solid requiring segmentation
0062     // to be a CartesianGridXY, throwing exception in getLocalSegmentation
0063     try {
0064       dd4hep::Box box = element->solid();
0065       m_xy_range_map.try_emplace(element, box->GetDX(), box->GetDY());
0066     } catch (const std::bad_cast& e) {
0067       error("Failed to cast solid to Box: {}", e.what());
0068     }
0069
0070     auto edep         = hit.getEDep();
0071     auto globalHitPos = hit.getPosition();
0072     auto hitPos =
0073         global2Local(dd4hep::Position(globalHitPos.x * dd4hep::mm, globalHitPos.y * dd4hep::mm,
0074                                       globalHitPos.z * dd4hep::mm),
0075                      transformIt->second);
0076
0077     std::unordered_set<dd4hep::rec::CellID> tested_cells;
0078     std::unordered_map<dd4hep::rec::CellID, float> cell_charge;
0079
0080     // Warning: This function is recursive, it stops shen it finds the edge of a detector element
0081     // or when the energy deposited in a cell is below the configured threshold
0082     findAllNeighborsInSensor(cellID, tested_cells, edep, hitPos, segmentationIt->second,
0083                              m_xy_range_map[element], hit, sharedHits);
0084
0085   } // for simhits
0086 } // SiliconChargeSharing:process
0087
0088 // Recursively find neighbors where a charge is deposited
0089 void SiliconChargeSharing::findAllNeighborsInSensor(
0090     const dd4hep::rec::CellID testCellID, std::unordered_set<dd4hep::rec::CellID>& tested_cells,
0091     const float edep, const dd4hep::Position hitPos,
0092     const dd4hep::DDSegmentation::CartesianGridXY* segmentation,
0093     const std::pair<double, double>& xy_range, const edm4hep::SimTrackerHit& hit,
0094     edm4hep::SimTrackerHitCollection* sharedHits) const {
0095
0096   // Tag cell as tested
0097   tested_cells.insert(testCellID);
0098
0099   auto localPos = cell2LocalPosition(testCellID);
0100
0101   // Check if the cell is within the sensor boundaries
0102   if (std::abs(localPos.x()) > xy_range.first || std::abs(localPos.y()) > xy_range.second) {
0103     return;
0104   }
0105
0106   // cout the local position and hit position
0107   auto xDimension = segmentation->gridSizeX();
0108   auto yDimension = segmentation->gridSizeY();
0109   // Calculate deposited energy in cell
0110   float edepCell = energyAtCell(xDimension, yDimension, localPos, hitPos, edep);
0111   if (edepCell <= m_cfg.min_edep) {
0112     return;
0113   }
0114
0115   // Create a new simhit for cell with deposited energy
0116   auto globalCellPos = m_converter->position(testCellID);
0117
0118   edm4hep::MutableSimTrackerHit shared_hit = hit.clone();
0119   shared_hit.setCellID(testCellID);
0120   shared_hit.setEDep(edepCell);
0121   shared_hit.setPosition({globalCellPos.x() / dd4hep::mm, globalCellPos.y() / dd4hep::mm,
0122                           globalCellPos.z() / dd4hep::mm});
0123   sharedHits->push_back(shared_hit);
0124
0125   // As there is charge in the cell, test the neighbors too
0126   std::set<dd4hep::rec::CellID> testCellNeighbours;
0127   segmentation->neighbours(testCellID, testCellNeighbours);
0128
0129   for (const auto& neighbourCell : testCellNeighbours) {
0130     if (tested_cells.find(neighbourCell) == tested_cells.end()) {
0131       findAllNeighborsInSensor(neighbourCell, tested_cells, edep, hitPos, segmentation, xy_range,
0132                                hit, sharedHits);
0133     }
0134   }
0135 }
0136
0137 // Calculate integral of Gaussian distribution
0138 float SiliconChargeSharing::integralGaus(float mean, float sd, float low_lim, float up_lim) {
0139   // return integral Gauss(mean, sd) dx from x = low_lim to x = up_lim
0140   // default value is set when sd = 0
0141   float up  = mean > up_lim ? -0.5 : 0.5;
0142   float low = mean > low_lim ? -0.5 : 0.5;
0143   if (sd > 0) {
0144     up  = -0.5 * std::erf(std::sqrt(2) * (mean - up_lim) / sd);
0145     low = -0.5 * std::erf(std::sqrt(2) * (mean - low_lim) / sd);
0146   }
0147   return up - low;
0148 }
0149
0150 // Convert cellID to local position
0151 dd4hep::Position SiliconChargeSharing::cell2LocalPosition(const dd4hep::rec::CellID& cell) const {
0152   auto position = m_seg->position(cell); // local position
0153   return position;
0154 }
0155
0156 // Convert global position to local position
0157 dd4hep::Position SiliconChargeSharing::global2Local(const dd4hep::Position& globalPosition,
0158                                                     const TGeoHMatrix* transform) {
0159
0160   double g[3];
0161   double l[3];
0162
0163   globalPosition.GetCoordinates(static_cast<Double_t*>(g));
0164   transform->MasterToLocal(static_cast<const Double_t*>(g), static_cast<Double_t*>(l));
0165   dd4hep::Position localPosition;
0166   localPosition.SetCoordinates(static_cast<const Double_t*>(l));
0167   return localPosition;
0168 }
0169
0170 // Calculate energy deposition in a cell relative to the hit position
0171 float SiliconChargeSharing::energyAtCell(const double xDimension, const double yDimension,
0172                                          const dd4hep::Position localPos,
0173                                          const dd4hep::Position hitPos, const float edep) const {
0174   float energy = edep *
0175                  integralGaus(hitPos.x(), m_cfg.sigma_sharingx, localPos.x() - 0.5 * xDimension,
0176                               localPos.x() + 0.5 * xDimension) *
0177                  integralGaus(hitPos.y(), m_cfg.sigma_sharingy, localPos.y() - 0.5 * yDimension,
0178                               localPos.y() + 0.5 * yDimension);
0179   return energy;
0180 }
0181
0182 // Get the segmentation relevant to a cellID
0183 const dd4hep::DDSegmentation::CartesianGridXY*
0184 SiliconChargeSharing::getLocalSegmentation(const dd4hep::rec::CellID& cellID) const {
0185   // Get the segmentation type
0186   auto segmentation_type                                   = m_seg.type();
0187   const dd4hep::DDSegmentation::Segmentation* segmentation = m_seg.segmentation();
0188   // Check if the segmentation is a multi-segmentation
0189   while (segmentation_type == "MultiSegmentation") {
0190     const auto* multi_segmentation =
0191         dynamic_cast<const dd4hep::DDSegmentation::MultiSegmentation*>(segmentation);
0192     segmentation      = &multi_segmentation->subsegmentation(cellID);
0193     segmentation_type = segmentation->type();
0194   }
0195
0196   // Try to cast the segmentation to CartesianGridXY
0197   const auto* cartesianGrid =
0198       dynamic_cast<const dd4hep::DDSegmentation::CartesianGridXY*>(segmentation);
0199   if (cartesianGrid == nullptr) {
0200     throw std::runtime_error("Segmentation is not of type CartesianGridXY");
0201   }
0202
0203   return cartesianGrid;
0204 }
0205
0206 } // namespace eicrecon