algorithms/calorimetry/CalorimeterHitReco.cc

0001 // SPDX-License-Identifier: LGPL-3.0-or-later
0002 // Copyright (C) 2022 Chao Peng, Sylvester Joosten, Wouter Deconinck, Chao, Whitney Armstrong
0003
0004 // Reconstruct digitized outputs, paired with Jug::Digi::CalorimeterHitDigi
0005 // Author: Chao Peng
0006 // Date: 06/14/2021
0007
0008 #include "CalorimeterHitReco.h"
0009
0010 #include <DD4hep/Alignments.h>
0011 #include <DD4hep/IDDescriptor.h>
0012 #include <DD4hep/Objects.h>
0013 #include <DD4hep/Readout.h>
0014 #include <DD4hep/Segmentations.h>
0015 #include <DD4hep/Shapes.h>
0016 #include <DD4hep/VolumeManager.h>
0017 #include <DD4hep/Volumes.h>
0018 #include <DD4hep/detail/SegmentationsInterna.h>
0019 #include <DDSegmentation/BitFieldCoder.h>
0020 #include <DDSegmentation/MultiSegmentation.h>
0021 #include <DDSegmentation/Segmentation.h>
0022 #include <Evaluator/DD4hepUnits.h>
0023 #include <Math/GenVector/Cartesian3D.h>
0024 #include <Math/GenVector/DisplacementVector3D.h>
0025 #include <algorithms/service.h>
0026 #include <edm4hep/Vector3f.h>
0027 #include <fmt/ranges.h>
0028 #include <algorithm>
0029 #include <cctype>
0030 #include <iterator>
0031 #include <sstream>
0032 #include <string>
0033 #include <tuple>
0034 #include <unordered_map>
0035 #include <utility>
0036 #include <vector>
0037
0038 #include "algorithms/calorimetry/CalorimeterHitRecoConfig.h"
0039 #include "services/evaluator/EvaluatorSvc.h"
0040
0041 using namespace dd4hep;
0042
0043 namespace eicrecon {
0044
0045 void CalorimeterHitReco::init() {
0046
0047   // threshold for firing
0048   // Should set either m_cfg.thresholdFactor or m_cfg.thresholdValue, not both
0049   if (m_cfg.thresholdFactor * m_cfg.thresholdValue != 0) {
0050     error("thresholdFactor = {}, thresholdValue = {}. Only one of these should be non-zero.",
0051           m_cfg.thresholdFactor, m_cfg.thresholdValue);
0052     throw; // throw with an argument doesn't trigger abort
0053   }
0054   thresholdADC = m_cfg.thresholdFactor * m_cfg.pedSigmaADC + m_cfg.thresholdValue;
0055   // TDC channels to timing conversion
0056   stepTDC = dd4hep::ns / m_cfg.resolutionTDC;
0057
0058   // do not get the layer/sector ID if no readout class provided
0059   if (m_cfg.readout.empty()) {
0060     return;
0061   }
0062
0063   // First, try and get the IDDescriptor. This will throw an exception if it fails.
0064   try {
0065     id_spec = m_detector->readout(m_cfg.readout).idSpec();
0066   } catch (...) {
0067     warning("Failed to get idSpec for {}", m_cfg.readout);
0068     return;
0069   }
0070   // Next, try and get the readout fields. This will throw a different exception.
0071   try {
0072     id_dec = id_spec.decoder();
0073     if (!m_cfg.sectorField.empty()) {
0074       sector_idx = id_dec->index(m_cfg.sectorField);
0075       debug("Find sector field {}, index = {}", m_cfg.sectorField, sector_idx);
0076     }
0077     if (!m_cfg.layerField.empty()) {
0078       layer_idx = id_dec->index(m_cfg.layerField);
0079       debug("Find layer field {}, index = {}", m_cfg.layerField, sector_idx);
0080     }
0081     if (!m_cfg.maskPosFields.empty()) {
0082       std::size_t tmp_mask = 0;
0083       for (auto& field : m_cfg.maskPosFields) {
0084         tmp_mask |= id_spec.field(field)->mask();
0085       }
0086       // assign this mask if all fields succeed
0087       gpos_mask = tmp_mask;
0088     }
0089   } catch (...) {
0090     if (id_dec == nullptr) {
0091       warning("Failed to load ID decoder for {}", m_cfg.readout);
0092       std::stringstream readouts;
0093       for (auto r : m_detector->readouts()) {
0094         readouts << "\"" << r.first << "\", ";
0095       }
0096       warning("Available readouts: {}", readouts.str());
0097     } else {
0098       warning("Failed to find field index for {}.", m_cfg.readout);
0099       if (!m_cfg.sectorField.empty()) {
0100         warning(" -- looking for sector field \"{}\".", m_cfg.sectorField);
0101       }
0102       if (!m_cfg.layerField.empty()) {
0103         warning(" -- looking for layer field  \"{}\".", m_cfg.layerField);
0104       }
0105       if (!m_cfg.maskPosFields.empty()) {
0106         warning(" -- looking for masking fields  \"{}\".", fmt::join(m_cfg.maskPosFields, ", "));
0107       }
0108       std::stringstream fields;
0109       for (auto field : id_spec.decoder()->fields()) {
0110         fields << "\"" << field.name() << "\", ";
0111       }
0112       warning("Available fields: {}", fields.str());
0113       warning("n.b. The local position, sector id and layer id will not be correct for this.");
0114       warning("Position masking may not be applied.");
0115       warning("however, the position, energy, and time values should still be good.");
0116     }
0117
0118     return;
0119   }
0120
0121   id_spec = m_detector->readout(m_cfg.readout).idSpec();
0122
0123   std::function hit_to_map = [this](const edm4hep::RawCalorimeterHit& h) {
0124     std::unordered_map<std::string, double> params;
0125     for (const auto& p : id_spec.fields()) {
0126       const std::string& name                  = p.first;
0127       const dd4hep::IDDescriptor::Field* field = p.second;
0128       params.emplace(name, field->value(h.getCellID()));
0129       trace("{} = {}", name, field->value(h.getCellID()));
0130     }
0131     return params;
0132   };
0133
0134   auto& serviceSvc = algorithms::ServiceSvc::instance();
0135   sampFrac = serviceSvc.service<EvaluatorSvc>("EvaluatorSvc")->compile(m_cfg.sampFrac, hit_to_map);
0136
0137   // local detector name has higher priority
0138   if (!m_cfg.localDetElement.empty()) {
0139     try {
0140       m_local = m_detector->detector(m_cfg.localDetElement);
0141       info("local coordinate system from DetElement {}", m_cfg.localDetElement);
0142     } catch (...) {
0143       error("failed to load local coordinate system from DetElement {}", m_cfg.localDetElement);
0144       return;
0145     }
0146   } else {
0147     std::vector<std::pair<std::string, int>> fields;
0148     for (auto f : m_cfg.localDetFields) {
0149       fields.emplace_back(f, 0);
0150     }
0151     local_mask = id_spec.get_mask(fields);
0152     // use all fields if nothing provided
0153     if (fields.empty()) {
0154       local_mask = ~static_cast<decltype(local_mask)>(0);
0155     }
0156   }
0157 }
0158
0159 void CalorimeterHitReco::process(const CalorimeterHitReco::Input& input,
0160                                  const CalorimeterHitReco::Output& output) const {
0161
0162   const auto [rawhits] = input;
0163   auto [recohits]      = output;
0164
0165   // For some detectors, the cellID in the raw hits may be broken
0166   // (currently this is the HcalBarrel). In this case, dd4hep
0167   // prints an error message and throws an exception. We catch
0168   // the exception and handle it, but the screen gets flooded
0169   // with these messages. Keep a count of these and if a max
0170   // number is encountered disable this algorithm. A useful message
0171   // indicating what is going on is printed below where the
0172   // error is detector.
0173   if (NcellIDerrors >= MaxCellIDerrors) {
0174     return;
0175   }
0176
0177   for (const auto& rh : *rawhits) {
0178
0179     //did not pass the zero-suppresion threshold
0180     const auto cellID = rh.getCellID();
0181     if (rh.getAmplitude() < m_cfg.pedMeanADC + thresholdADC) {
0182       continue;
0183     }
0184
0185     if (rh.getAmplitude() > static_cast<int>(m_cfg.capADC)) {
0186       error("Encountered hit with amplitude {} outside of ADC capacity {}", rh.getAmplitude(),
0187             m_cfg.capADC);
0188       continue;
0189     }
0190
0191     // get layer and sector ID
0192     const int lid = id_dec != nullptr && !m_cfg.layerField.empty()
0193                         ? static_cast<int>(id_dec->get(cellID, layer_idx))
0194                         : -1;
0195     const int sid = id_dec != nullptr && !m_cfg.sectorField.empty()
0196                         ? static_cast<int>(id_dec->get(cellID, sector_idx))
0197                         : -1;
0198
0199     // convert ADC to energy
0200     float sampFrac_value = sampFrac(rh);
0201     float energy         = (((signed)rh.getAmplitude() - (signed)m_cfg.pedMeanADC)) /
0202                            static_cast<float>(m_cfg.capADC) * m_cfg.dyRangeADC / sampFrac_value;
0203
0204     const float time = rh.getTimeStamp() / stepTDC;
0205     trace("cellID {}, \t energy: {},  TDC: {}, time: {}, sampFrac: {}", cellID, energy,
0206           rh.getTimeStamp(), time, sampFrac_value);
0207
0208     dd4hep::DetElement local;
0209     dd4hep::Position gpos;
0210     try {
0211       // global positions
0212       gpos = m_converter->position(cellID);
0213
0214       // masked position (look for a mother volume)
0215       if (gpos_mask != 0) {
0216         auto mpos = m_converter->position(cellID & ~gpos_mask);
0217         // replace corresponding coords
0218         for (const char& c : m_cfg.maskPos) {
0219           switch (std::tolower(c)) {
0220           case 'x':
0221             gpos.SetX(mpos.X());
0222             break;
0223           case 'y':
0224             gpos.SetY(mpos.Y());
0225             break;
0226           case 'z':
0227             gpos.SetZ(mpos.Z());
0228             break;
0229           default:
0230             break;
0231           }
0232         }
0233       }
0234
0235       // local positions
0236       if (m_cfg.localDetElement.empty()) {
0237         auto volman = m_detector->volumeManager();
0238         local       = volman.lookupDetElement(cellID & local_mask);
0239       } else {
0240         local = m_local;
0241       }
0242     } catch (...) {
0243       // Error looking up cellID. Messages should already have been printed.
0244       // Also, see comment at top of this method.
0245       if (++NcellIDerrors >= MaxCellIDerrors) {
0246         error("Maximum number of errors reached: {}", MaxCellIDerrors);
0247         error("This is likely an issue with the cellID being unknown.");
0248         error("Note: local_mask={:X} example cellID={:x}", local_mask, cellID);
0249         error("Disabling this algorithm since it requires a valid cellID.");
0250         error("(See {}:{})", __FILE__, __LINE__);
0251       }
0252       continue;
0253     }
0254
0255     const auto pos = local.nominal().worldToLocal(gpos);
0256     std::vector<double> cdim;
0257     // get segmentation dimensions
0258
0259     const dd4hep::DDSegmentation::Segmentation* segmentation =
0260         m_converter->findReadout(local).segmentation()->segmentation;
0261     auto segmentation_type = segmentation->type();
0262
0263     while (segmentation_type == "MultiSegmentation") {
0264       const auto* multi_segmentation =
0265           dynamic_cast<const dd4hep::DDSegmentation::MultiSegmentation*>(segmentation);
0266       const dd4hep::DDSegmentation::Segmentation& sub_segmentation =
0267           multi_segmentation->subsegmentation(cellID);
0268
0269       segmentation      = &sub_segmentation;
0270       segmentation_type = segmentation->type();
0271     }
0272
0273     if (segmentation_type == "CartesianGridXY" || segmentation_type == "HexGridXY" ||
0274         segmentation_type == "CartesianGridXYStaggered") {
0275       auto cell_dim = m_converter->cellDimensions(cellID);
0276       cdim.resize(3);
0277       cdim[0] = cell_dim[0];
0278       cdim[1] = cell_dim[1];
0279       debug("Using segmentation for cell dimensions: {}", fmt::join(cdim, ", "));
0280     } else if (segmentation_type == "CartesianStripZ") {
0281       auto cell_dim = m_converter->cellDimensions(cellID);
0282       cdim.resize(3);
0283       cdim[2] = cell_dim[0];
0284       debug("Using segmentation for cell dimensions: {}", fmt::join(cdim, ", "));
0285     } else {
0286       if ((segmentation_type != "NoSegmentation") && (!warned_unsupported_segmentation)) {
0287         warning("Unsupported segmentation type \"{}\"", segmentation_type);
0288         warned_unsupported_segmentation = true;
0289       }
0290
0291       // Using bounding box instead of actual solid so the dimensions are always in dim_x, dim_y, dim_z
0292       cdim =
0293           m_converter->findContext(cellID)->volumePlacement().volume().boundingBox().dimensions();
0294       std::ranges::transform(cdim, cdim.begin(), [](auto&& PH1) {
0295         return std::multiplies<double>()(std::forward<decltype(PH1)>(PH1), 2);
0296       });
0297       debug("Using bounding box for cell dimensions: {}", fmt::join(cdim, ", "));
0298     }
0299
0300     //create constant vectors for passing to hit initializer list
0301     //FIXME: needs to come from the geometry service/converter
0302     const decltype(edm4eic::CalorimeterHitData::position) position(
0303         gpos.x() / dd4hep::mm, gpos.y() / dd4hep::mm, gpos.z() / dd4hep::mm);
0304     const decltype(edm4eic::CalorimeterHitData::dimension) dimension(
0305         cdim.at(0) / dd4hep::mm, cdim.at(1) / dd4hep::mm, cdim.at(2) / dd4hep::mm);
0306     const decltype(edm4eic::CalorimeterHitData::local) local_position(
0307         pos.x() / dd4hep::mm, pos.y() / dd4hep::mm, pos.z() / dd4hep::mm);
0308
0309     auto recohit = recohits->create(rh.getCellID(), energy, 0, time, 0, position, dimension, sid,
0310                                     lid, local_position);
0311     recohit.setRawHit(rh);
0312   }
0313 }
0314
0315 } // namespace eicrecon