src/components/CalorimeterHitDigi.cpp

0001 // SPDX-License-Identifier: LGPL-3.0-or-later
0002 // Copyright (C) 2022 Chao Peng, Wouter Deconinck, Sylvester Joosten
0003
0004 // A general digitization for CalorimeterHit from simulation
0005 // 1. Smear energy deposit with a/sqrt(E/GeV) + b + c/E or a/sqrt(E/GeV) (relative value)
0006 // 2. Digitize the energy with dynamic ADC range and add pedestal (mean +- sigma)
0007 // 3. Time conversion with smearing resolution (absolute value)
0008 // 4. Signal is summed if the SumFields are provided
0009 //
0010 // Author: Chao Peng
0011 // Date: 06/02/2021
0012
0013 #include <algorithm>
0014 #include <cmath>
0015 #include <unordered_map>
0016
0017 #include "GaudiAlg/GaudiAlgorithm.h"
0018 #include "GaudiAlg/GaudiTool.h"
0019 #include "GaudiAlg/Transformer.h"
0020 #include "GaudiKernel/PhysicalConstants.h"
0021 #include "Gaudi/Property.h"
0022 #include "GaudiKernel/RndmGenerators.h"
0023
0024 #include "DDRec/CellIDPositionConverter.h"
0025 #include "DDSegmentation/BitFieldCoder.h"
0026
0027 #include <k4Interface/IGeoSvc.h>
0028 #include <k4FWCore/DataHandle.h>
0029
0030 #include "fmt/format.h"
0031 #include "fmt/ranges.h"
0032
0033 // Event Model related classes
0034 #include "edm4hep/SimCalorimeterHitCollection.h"
0035 #include "edm4hep/RawCalorimeterHitCollection.h"
0036
0037
0038 using namespace Gaudi::Units;
0039
0040 namespace Jug::Digi {
0041
0042   /** Generic calorimeter hit digitiziation.
0043    *
0044    * \ingroup digi
0045    * \ingroup calorimetry
0046    */
0047   class CalorimeterHitDigi : public GaudiAlgorithm {
0048   public:
0049     // additional smearing resolutions
0050     Gaudi::Property<std::vector<double>> u_eRes{this, "energyResolutions", {}}; // a/sqrt(E/GeV) + b + c/(E/GeV)
0051     Gaudi::Property<double>              m_tRes{this, "timeResolution", 0.0 * ns};
0052     // single hit energy deposition threshold
0053     Gaudi::Property<double>              m_threshold{this, "threshold", 1. * keV};
0054
0055     // digitization settings
0056     Gaudi::Property<unsigned int>       m_capADC{this, "capacityADC", 8096};
0057     Gaudi::Property<double>             m_dyRangeADC{this, "dynamicRangeADC", 100 * MeV};
0058     Gaudi::Property<unsigned int>       m_pedMeanADC{this, "pedestalMean", 400};
0059     Gaudi::Property<double>             m_pedSigmaADC{this, "pedestalSigma", 3.2};
0060     Gaudi::Property<double>             m_resolutionTDC{this, "resolutionTDC", 10 * ps};
0061
0062     Gaudi::Property<double>             m_corrMeanScale{this, "scaleResponse", 1.0};
0063     // These are better variable names for the "energyResolutions" array which is a bit
0064     // magic @FIXME
0065     //Gaudi::Property<double>             m_corrSigmaCoeffE{this, "responseCorrectionSigmaCoeffE", 0.0};
0066     //Gaudi::Property<double>             m_corrSigmaCoeffSqrtE{this, "responseCorrectionSigmaCoeffSqrtE", 0.0};
0067
0068     // signal sums
0069     // @TODO: implement signal sums with timing
0070     // field names to generate id mask, the hits will be grouped by masking the field
0071     Gaudi::Property<std::vector<std::string>> u_fields{this, "signalSumFields", {}};
0072     // ref field ids are used for the merged hits, 0 is used if nothing provided
0073     Gaudi::Property<std::vector<int>>         u_refs{this, "fieldRefNumbers", {}};
0074     Gaudi::Property<std::string>              m_geoSvcName{this, "geoServiceName", "GeoSvc"};
0075     Gaudi::Property<std::string>              m_readout{this, "readoutClass", ""};
0076
0077     // unitless counterparts of inputs
0078     double           dyRangeADC{0}, stepTDC{0}, tRes{0}, eRes[3] = {0., 0., 0.};
0079     Rndm::Numbers    m_normDist;
0080     SmartIF<IGeoSvc> m_geoSvc;
0081     uint64_t         id_mask{0}, ref_mask{0};
0082
0083     DataHandle<edm4hep::SimCalorimeterHitCollection> m_inputHitCollection{
0084       "inputHitCollection", Gaudi::DataHandle::Reader, this};
0085     DataHandle<edm4hep::RawCalorimeterHitCollection> m_outputHitCollection{
0086       "outputHitCollection", Gaudi::DataHandle::Writer, this};
0087
0088     //  ill-formed: using GaudiAlgorithm::GaudiAlgorithm;
0089     CalorimeterHitDigi(const std::string& name, ISvcLocator* svcLoc)
0090       : GaudiAlgorithm(name, svcLoc) {
0091       declareProperty("inputHitCollection", m_inputHitCollection, "");
0092       declareProperty("outputHitCollection", m_outputHitCollection, "");
0093     }
0094
0095     StatusCode initialize() override
0096     {
0097       if (GaudiAlgorithm::initialize().isFailure()) {
0098         return StatusCode::FAILURE;
0099       }
0100       // random number generator from service
0101       auto randSvc = svc<IRndmGenSvc>("RndmGenSvc", true);
0102       auto sc      = m_normDist.initialize(randSvc, Rndm::Gauss(0.0, 1.0));
0103       if (!sc.isSuccess()) {
0104         return StatusCode::FAILURE;
0105       }
0106       // set energy resolution numbers
0107       for (size_t i = 0; i < u_eRes.size() && i < 3; ++i) {
0108         eRes[i] = u_eRes[i];
0109       }
0110
0111       // using juggler internal units (GeV, mm, radian, ns)
0112       dyRangeADC = m_dyRangeADC.value() / GeV;
0113       tRes       = m_tRes.value() / ns;
0114       stepTDC    = ns / m_resolutionTDC.value();
0115
0116       // need signal sum
0117       if (!u_fields.value().empty()) {
0118         m_geoSvc = service(m_geoSvcName);
0119         // sanity checks
0120         if (!m_geoSvc) {
0121           error() << "Unable to locate Geometry Service. "
0122                   << "Make sure you have GeoSvc and SimSvc in the right order in the configuration."
0123                   << endmsg;
0124           return StatusCode::FAILURE;
0125         }
0126         if (m_readout.value().empty()) {
0127           error() << "readoutClass is not provided, it is needed to know the fields in readout ids"
0128                   << endmsg;
0129           return StatusCode::FAILURE;
0130         }
0131
0132         // get decoders
0133         try {
0134           auto id_desc = m_geoSvc->getDetector()->readout(m_readout).idSpec();
0135           id_mask      = 0;
0136           std::vector<std::pair<std::string, int>> ref_fields;
0137           for (size_t i = 0; i < u_fields.size(); ++i) {
0138             id_mask |= id_desc.field(u_fields[i])->mask();
0139             // use the provided id number to find ref cell, or use 0
0140             int ref = i < u_refs.size() ? u_refs[i] : 0;
0141             ref_fields.emplace_back(u_fields[i], ref);
0142           }
0143           ref_mask = id_desc.encode(ref_fields);
0144           // debug() << fmt::format("Referece id mask for the fields {:#064b}", ref_mask) << endmsg;
0145         } catch (...) {
0146           error() << "Failed to load ID decoder for " << m_readout << endmsg;
0147           return StatusCode::FAILURE;
0148         }
0149         id_mask = ~id_mask;
0150         info() << fmt::format("ID mask in {:s}: {:#064b}", m_readout.value(), id_mask) << endmsg;
0151         return StatusCode::SUCCESS;
0152       }
0153
0154       return StatusCode::SUCCESS;
0155     }
0156
0157     StatusCode execute() override
0158     {
0159       if (!u_fields.value().empty()) {
0160         signal_sum_digi();
0161       } else {
0162         single_hits_digi();
0163       }
0164       return StatusCode::SUCCESS;
0165     }
0166
0167   private:
0168     void single_hits_digi() {
0169       // input collections
0170       const auto* const simhits = m_inputHitCollection.get();
0171       // Create output collections
0172       auto* rawhits = m_outputHitCollection.createAndPut();
0173       for (const auto& ahit : *simhits) {
0174         // Note: juggler internal unit of energy is GeV
0175         const double eDep    = ahit.getEnergy();
0176
0177         // apply additional calorimeter noise to corrected energy deposit
0178         const double eResRel = (eDep > m_threshold)
0179             ? m_normDist() * std::sqrt(
0180                   std::pow(eRes[0] / std::sqrt(eDep), 2) +
0181                   std::pow(eRes[1], 2) +
0182                   std::pow(eRes[2] / (eDep), 2)
0183               )
0184             : 0;
0185
0186         const double ped    = m_pedMeanADC + m_normDist() * m_pedSigmaADC;
0187         const long long adc = std::llround(ped +  eDep * (m_corrMeanScale + eResRel) / dyRangeADC * m_capADC);
0188
0189         double time = std::numeric_limits<double>::max();
0190         for (const auto& c : ahit.getContributions()) {
0191           if (c.getTime() <= time) {
0192             time = c.getTime();
0193           }
0194         }
0195         const long long tdc = std::llround((time + m_normDist() * tRes) * stepTDC);
0196
0197         rawhits->create(
0198           ahit.getCellID(),
0199           (adc > m_capADC.value() ? m_capADC.value() : adc),
0200           tdc
0201         );
0202       }
0203     }
0204
0205     void signal_sum_digi() {
0206       const auto* const simhits = m_inputHitCollection.get();
0207       auto* rawhits = m_outputHitCollection.createAndPut();
0208
0209       // find the hits that belong to the same group (for merging)
0210       std::unordered_map<long long, std::vector<edm4hep::SimCalorimeterHit>> merge_map;
0211       for (const auto &ahit : *simhits) {
0212         int64_t hid = (ahit.getCellID() & id_mask) | ref_mask;
0213         auto    it  = merge_map.find(hid);
0214
0215         if (it == merge_map.end()) {
0216           merge_map[hid] = {ahit};
0217         } else {
0218           it->second.push_back(ahit);
0219         }
0220       }
0221
0222       // signal sum
0223       for (auto &[id, hits] : merge_map) {
0224         double edep     = hits[0].getEnergy();
0225         double time     = hits[0].getContributions(0).getTime();
0226         double max_edep = hits[0].getEnergy();
0227         // sum energy, take time from the most energetic hit
0228         for (size_t i = 1; i < hits.size(); ++i) {
0229           edep += hits[i].getEnergy();
0230           if (hits[i].getEnergy() > max_edep) {
0231             max_edep = hits[i].getEnergy();
0232             for (const auto& c : hits[i].getContributions()) {
0233               if (c.getTime() <= time) {
0234                 time = c.getTime();
0235               }
0236             }
0237           }
0238         }
0239
0240         // safety check
0241         const double eResRel = (edep > m_threshold)
0242             ? m_normDist() * eRes[0] / std::sqrt(edep) +
0243               m_normDist() * eRes[1] +
0244               m_normDist() * eRes[2] / edep
0245             : 0;
0246
0247         double    ped     = m_pedMeanADC + m_normDist() * m_pedSigmaADC;
0248         unsigned long long adc     = std::llround(ped + edep * (1. + eResRel) / dyRangeADC * m_capADC);
0249         unsigned long long tdc     = std::llround((time + m_normDist() * tRes) * stepTDC);
0250
0251         rawhits->create(
0252           id,
0253           (adc > m_capADC.value() ? m_capADC.value() : adc),
0254           tdc
0255         );
0256       }
0257     }
0258   };
0259   // NOLINTNEXTLINE(cppcoreguidelines-avoid-non-const-global-variables)
0260   DECLARE_COMPONENT(CalorimeterHitDigi)
0261
0262 } // namespace Jug::Digi