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