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