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