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