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