src/todo/CalorimeterHitReco.cpp

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 "fmt/format.h"
0009 #include "fmt/ranges.h"
0010 #include <algorithm>
0011 #include <bitset>
0012
0013 #include "Gaudi/Property.h"
0014 #include "GaudiAlg/GaudiAlgorithm.h"
0015 #include "GaudiAlg/GaudiTool.h"
0016 #include "GaudiAlg/Transformer.h"
0017 #include "GaudiKernel/PhysicalConstants.h"
0018 #include "GaudiKernel/RndmGenerators.h"
0019 #include "GaudiKernel/ToolHandle.h"
0020
0021 #include "DDRec/CellIDPositionConverter.h"
0022 #include "DDRec/Surface.h"
0023 #include "DDRec/SurfaceManager.h"
0024
0025 #include "JugBase/DataHandle.h"
0026 #include "JugBase/IGeoSvc.h"
0027
0028 // Event Model related classes
0029 #include "edm4eic/CalorimeterHitCollection.h"
0030 #include "edm4eic/RawCalorimeterHitCollection.h"
0031
0032 using namespace Gaudi::Units;
0033
0034 namespace Jug::Reco {
0035
0036 /** Calorimeter hit reconstruction.
0037  *
0038  * Reconstruct digitized outputs, paired with Jug::Digi::CalorimeterHitDigi
0039  * \ingroup reco
0040  */
0041 class CalorimeterHitReco : public GaudiAlgorithm {
0042 private:
0043   // length unit from dd4hep, should be fixed
0044   Gaudi::Property<double> m_lUnit{this, "lengthUnit", dd4hep::mm};
0045
0046   // digitization settings, must be consistent with digi class
0047   Gaudi::Property<unsigned int> m_capADC{this, "capacityADC", 8096};
0048   Gaudi::Property<double> m_dyRangeADC{this, "dynamicRangeADC", 100. * MeV};
0049   Gaudi::Property<unsigned int> m_pedMeanADC{this, "pedestalMean", 400};
0050   Gaudi::Property<double> m_pedSigmaADC{this, "pedestalSigma", 3.2};
0051   Gaudi::Property<double> m_resolutionTDC{this, "resolutionTDC", 10 * ps};
0052
0053   // zero suppression values
0054   Gaudi::Property<double> m_thresholdFactor{this, "thresholdFactor", 0.0};
0055   Gaudi::Property<double> m_thresholdValue{this, "thresholdValue", 0.0};
0056
0057   // energy correction with sampling fraction
0058   Gaudi::Property<double> m_sampFrac{this, "samplingFraction", 1.0};
0059
0060   // unitless counterparts of the input parameters
0061   double dyRangeADC{0};
0062   double thresholdADC{0};
0063   double stepTDC{0};
0064
0065   DataHandle<edm4eic::RawCalorimeterHitCollection> m_inputHitCollection{"inputHitCollection", Gaudi::DataHandle::Reader,
0066                                                                     this};
0067   DataHandle<edm4eic::CalorimeterHitCollection> m_outputHitCollection{"outputHitCollection", Gaudi::DataHandle::Writer,
0068                                                                   this};
0069
0070   // geometry service to get ids, ignored if no names provided
0071   Gaudi::Property<std::string> m_geoSvcName{this, "geoServiceName", "GeoSvc"};
0072   Gaudi::Property<std::string> m_readout{this, "readoutClass", ""};
0073   Gaudi::Property<std::string> m_layerField{this, "layerField", ""};
0074   Gaudi::Property<std::string> m_sectorField{this, "sectorField", ""};
0075   SmartIF<IGeoSvc> m_geoSvc;
0076   dd4hep::BitFieldCoder* id_dec = nullptr;
0077   size_t sector_idx{0}, layer_idx{0};
0078
0079   // name of detelment or fields to find the local detector (for global->local transform)
0080   // if nothing is provided, the lowest level DetElement (from cellID) will be used
0081   Gaudi::Property<std::string> m_localDetElement{this, "localDetElement", ""};
0082   Gaudi::Property<std::vector<std::string>> u_localDetFields{this, "localDetFields", {}};
0083   dd4hep::DetElement local;
0084   size_t local_mask = ~0;
0085
0086 public:
0087   CalorimeterHitReco(const std::string& name, ISvcLocator* svcLoc) : GaudiAlgorithm(name, svcLoc) {
0088     declareProperty("inputHitCollection", m_inputHitCollection, "");
0089     declareProperty("outputHitCollection", m_outputHitCollection, "");
0090   }
0091
0092   StatusCode initialize() override {
0093     if (GaudiAlgorithm::initialize().isFailure()) {
0094       return StatusCode::FAILURE;
0095     }
0096
0097     m_geoSvc = service(m_geoSvcName);
0098     if (!m_geoSvc) {
0099       error() << "Unable to locate Geometry Service. "
0100               << "Make sure you have GeoSvc and SimSvc in the right order in the configuration." << endmsg;
0101       return StatusCode::FAILURE;
0102     }
0103
0104     // unitless conversion
0105     dyRangeADC = m_dyRangeADC.value() / GeV;
0106
0107     // threshold for firing
0108     thresholdADC = m_thresholdFactor.value() * m_pedSigmaADC.value() + m_thresholdValue.value();
0109
0110     // TDC channels to timing conversion
0111     stepTDC = ns / m_resolutionTDC.value();
0112
0113     // do not get the layer/sector ID if no readout class provided
0114     if (m_readout.value().empty()) {
0115       return StatusCode::SUCCESS;
0116     }
0117
0118     auto id_spec = m_geoSvc->detector()->readout(m_readout).idSpec();
0119     try {
0120       id_dec = id_spec.decoder();
0121       if (!m_sectorField.value().empty()) {
0122         sector_idx = id_dec->index(m_sectorField);
0123         info() << "Find sector field " << m_sectorField.value() << ", index = " << sector_idx << endmsg;
0124       }
0125       if (!m_layerField.value().empty()) {
0126         layer_idx = id_dec->index(m_layerField);
0127         info() << "Find layer field " << m_layerField.value() << ", index = " << sector_idx << endmsg;
0128       }
0129     } catch (...) {
0130       error() << "Failed to load ID decoder for " << m_readout << endmsg;
0131       return StatusCode::FAILURE;
0132     }
0133
0134     // local detector name has higher priority
0135     if (!m_localDetElement.value().empty()) {
0136       try {
0137         local = m_geoSvc->detector()->detector(m_localDetElement.value());
0138         info() << "Local coordinate system from DetElement " << m_localDetElement.value() << endmsg;
0139       } catch (...) {
0140         error() << "Failed to locate local coordinate system from DetElement " << m_localDetElement.value() << endmsg;
0141         return StatusCode::FAILURE;
0142       }
0143       // or get from fields
0144     } else {
0145       std::vector<std::pair<std::string, int>> fields;
0146       for (auto& f : u_localDetFields.value()) {
0147         fields.emplace_back(f, 0);
0148       }
0149       local_mask = id_spec.get_mask(fields);
0150       // use all fields if nothing provided
0151       if (fields.empty()) {
0152         local_mask = ~0;
0153       }
0154       info() << fmt::format("Local DetElement mask {:#064b} from fields [{}]", local_mask, fmt::join(fields, ", "))
0155              << endmsg;
0156     }
0157
0158     return StatusCode::SUCCESS;
0159   }
0160
0161   StatusCode execute() override {
0162     // input collections
0163     const auto& rawhits = *m_inputHitCollection.get();
0164     // create output collections
0165     auto& hits     = *m_outputHitCollection.createAndPut();
0166     auto converter = m_geoSvc->cellIDPositionConverter();
0167
0168     // energy time reconstruction
0169     for (const auto& rh : rawhits) {
0170
0171       #pragma GCC diagnostic push
0172       #pragma GCC diagnostic error "-Wsign-conversion"
0173
0174       // did not pass the zero-suppression threshold
0175       if (rh.getAmplitude() < m_pedMeanADC + thresholdADC) {
0176         continue;
0177       }
0178
0179       // convert ADC -> energy
0180       const float energy =
0181         (((signed)rh.getAmplitude() - (signed)m_pedMeanADC)) / static_cast<float>(m_capADC.value()) * dyRangeADC / m_sampFrac;
0182       const float time  = rh.getTimeStamp() / stepTDC;
0183
0184       #pragma GCC diagnostic pop
0185
0186       const auto cellID = rh.getCellID();
0187       const int lid = id_dec != nullptr && !m_layerField.value().empty() ? static_cast<int>(id_dec->get(cellID, layer_idx)) : -1;
0188       const int sid = id_dec != nullptr && !m_sectorField.value().empty() ? static_cast<int>(id_dec->get(cellID, sector_idx)) : -1;
0189       // global positions
0190       const auto gpos = converter->position(cellID);
0191       // local positions
0192       if (m_localDetElement.value().empty()) {
0193         auto volman = m_geoSvc->detector()->volumeManager();
0194         local       = volman.lookupDetElement(cellID & local_mask);
0195         }
0196         const auto pos = local.nominal().worldToLocal(dd4hep::Position(gpos.x(), gpos.y(), gpos.z()));
0197         // cell dimension
0198         std::vector<double> cdim;
0199         // get segmentation dimensions
0200         if (converter->findReadout(local).segmentation().type() != "NoSegmentation") {
0201         cdim = converter->cellDimensions(cellID);
0202         // get volume dimensions (multiply by two to get fullsize)
0203         } else {
0204         // Using bounding box instead of actual solid so the dimensions are always in dim_x, dim_y, dim_z
0205         cdim = converter->findContext(cellID)->volumePlacement().volume().boundingBox().dimensions();
0206         std::transform(cdim.begin(), cdim.end(), cdim.begin(),
0207                        std::bind(std::multiplies<double>(), std::placeholders::_1, 2));
0208         }
0209
0210         // create const vectors for passing to hit initializer list
0211         const decltype(edm4eic::CalorimeterHitData::position) position(
0212           gpos.x() / m_lUnit, gpos.y() / m_lUnit, gpos.z() / m_lUnit
0213         );
0214         const decltype(edm4eic::CalorimeterHitData::dimension) dimension(
0215           cdim[0] / m_lUnit, cdim[1] / m_lUnit, cdim[2] / m_lUnit
0216         );
0217         const decltype(edm4eic::CalorimeterHitData::local) local_position(
0218           pos.x() / m_lUnit, pos.y() / m_lUnit, pos.z() / m_lUnit
0219         );
0220
0221         hits.push_back({
0222             rh.getCellID(), // cellID
0223             energy,         // energy
0224             0,              // @TODO: energy error
0225             time,           // time
0226             0,              // time error FIXME should be configurable
0227             position,       // global pos
0228             dimension,
0229             // Local hit info
0230             sid,
0231             lid,
0232             local_position, // local pos
0233         });
0234     }
0235
0236     return StatusCode::SUCCESS;
0237   }
0238
0239 }; // class CalorimeterHitReco
0240
0241 // NOLINTNEXTLINE(cppcoreguidelines-avoid-non-const-global-variables)
0242 DECLARE_COMPONENT(CalorimeterHitReco)
0243
0244 } // namespace Jug::Reco