algorithms/calorimetry/CalorimeterHitDigi.cc

0001 // SPDX-License-Identifier: LGPL-3.0-or-later
0002 // Copyright (C) 2022 Chao Peng, Wouter Deconinck, Sylvester Joosten, Barak Schmookler, David Lawrence
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
0014 #include "CalorimeterHitDigi.h"
0015
0016 #include <DD4hep/Detector.h>
0017 #include <DD4hep/IDDescriptor.h>
0018 #include <DD4hep/Readout.h>
0019 #include <DD4hep/config.h>
0020 #include <DDSegmentation/BitFieldCoder.h>
0021 #include <Evaluator/DD4hepUnits.h>
0022 #include <algorithms/service.h>
0023 #include <edm4eic/EDM4eicVersion.h>
0024 #if EDM4EIC_VERSION_MAJOR >= 7
0025 #include <edm4eic/MCRecoCalorimeterHitAssociationCollection.h>
0026 #endif
0027 #include <edm4hep/CaloHitContributionCollection.h>
0028 #include <fmt/core.h>
0029 #include <podio/RelationRange.h>
0030 #include <algorithm>
0031 #include <cmath>
0032 #include <cstddef>
0033 #include <gsl/pointers>
0034 #include <limits>
0035 #include <stdexcept>
0036 #include <string>
0037 #include <unordered_map>
0038 #include <utility>
0039 #include <vector>
0040
0041 #include "algorithms/calorimetry/CalorimeterHitDigiConfig.h"
0042 #include "services/evaluator/EvaluatorSvc.h"
0043
0044 using namespace dd4hep;
0045
0046 namespace eicrecon {
0047
0048 //
0049 // TODO:
0050 // - Array type configuration parameters are not yet supported in JANA (needs to be added)
0051 // - Random number service needs to bew resolved (on global scale)
0052 // - It is possible standard running of this with Gaudi relied on a number of parameters
0053 //   being set in the config. If that is the case, they should be moved into the default
0054 //   values here. This needs to be confirmed.
0055
0056
0057 void CalorimeterHitDigi::init() {
0058
0059     // Gaudi implements a random number generator service. It is not clear to me how this
0060     // can work. There are multiple race conditions that occur in parallel event processing:
0061     // 1. The exact same events processed by a given thread in one invocation will not
0062     //    necessarily be the combination of events any thread sees in a subsequent
0063     //    invocation. Thus, you can't rely on thread_local storage.
0064     // 2. Its possible for the factory execution order to be modified by the presence of
0065     //    a processor (e.g. monitoring plugin). This is not as serious since changing the
0066     //    command line should cause one not to expect reproducibility. Still, one may
0067     //    expect the inclusion of an "observer" plugin not to have such side affects.
0068     //
0069     // More information will be needed. In the meantime, we implement a local random number
0070     // generator. Ideally, this would be seeded with the run number+event number, but for
0071     // now, just use default values defined in header file.
0072
0073     // set energy resolution numbers
0074     if (m_cfg.eRes.empty()) {
0075       m_cfg.eRes.resize(3);
0076     } else if (m_cfg.eRes.size() != 3) {
0077       error("Invalid m_cfg.eRes.size()");
0078       throw std::runtime_error("Invalid m_cfg.eRes.size()");
0079     }
0080
0081     // using juggler internal units (GeV, mm, radian, ns)
0082     tRes       = m_cfg.tRes / dd4hep::ns;
0083     stepTDC    = dd4hep::ns / m_cfg.resolutionTDC;
0084
0085     // sanity checks
0086     if (m_cfg.readout.empty()) {
0087         error("readoutClass is not provided, it is needed to know the fields in readout ids");
0088         throw std::runtime_error("readoutClass is not provided");
0089     }
0090
0091     // get decoders
0092     try {
0093         id_spec = m_geo.detector()->readout(m_cfg.readout).idSpec();
0094     } catch (...) {
0095         // Can not be more verbose. In JANA2, this will be attempted at each event, which
0096         // pollutes output for geometries that are less than complete.
0097         // We could save an exception and throw it from process.
0098         debug("Failed to load ID decoder for {}", m_cfg.readout);
0099         throw std::runtime_error(fmt::format("Failed to load ID decoder for {}", m_cfg.readout));
0100     }
0101
0102     decltype(id_mask) id_inverse_mask = 0;
0103     // all these are for signal sum at digitization level
0104     if (!m_cfg.fields.empty()) {
0105         for (auto & field : m_cfg.fields) {
0106             id_inverse_mask |= id_spec.field(field)->mask();
0107         }
0108         debug("ID mask in {:s}: {:#064b}", m_cfg.readout, id_mask);
0109     }
0110     id_mask = ~id_inverse_mask;
0111
0112     std::function hit_to_map = [this](const edm4hep::SimCalorimeterHit &h) {
0113       std::unordered_map<std::string, double> params;
0114       for(const auto &p : id_spec.fields()) {
0115         const std::string &name = p.first;
0116         const dd4hep::IDDescriptor::Field* field = p.second;
0117         params.emplace(name, field->value(h.getCellID()));
0118         trace("{} = {}", name, field->value(h.getCellID()));
0119       }
0120       return params;
0121     };
0122
0123     auto& serviceSvc = algorithms::ServiceSvc::instance();
0124     corrMeanScale = serviceSvc.service<EvaluatorSvc>("EvaluatorSvc")->compile(m_cfg.corrMeanScale, hit_to_map);
0125 }
0126
0127
0128 void CalorimeterHitDigi::process(
0129       const CalorimeterHitDigi::Input& input,
0130       const CalorimeterHitDigi::Output& output) const {
0131
0132     const auto [simhits] = input;
0133 #if EDM4EIC_VERSION_MAJOR >= 7
0134     auto [rawhits, rawassocs] = output;
0135 #else
0136     auto [rawhits] = output;
0137 #endif
0138
0139     // find the hits that belong to the same group (for merging)
0140     std::unordered_map<uint64_t, std::vector<std::size_t>> merge_map;
0141     std::size_t ix = 0;
0142     for (const auto &ahit : *simhits) {
0143         uint64_t hid = ahit.getCellID() & id_mask;
0144
0145         trace("org cell ID in {:s}: {:#064b}", m_cfg.readout, ahit.getCellID());
0146         trace("new cell ID in {:s}: {:#064b}", m_cfg.readout, hid);
0147
0148         merge_map[hid].push_back(ix);
0149
0150         ix++;
0151     }
0152
0153     // signal sum
0154     // NOTE: we take the cellID of the most energetic hit in this group so it is a real cellID from an MC hit
0155     for (const auto &[id, ixs] : merge_map) {
0156
0157         // create hit and association in advance
0158         edm4hep::MutableRawCalorimeterHit rawhit;
0159 #if EDM4EIC_VERSION_MAJOR >= 7
0160         std::vector<edm4eic::MutableMCRecoCalorimeterHitAssociation> rawassocs_staging;
0161 #endif
0162
0163         double edep     = 0;
0164         double time     = std::numeric_limits<double>::max();
0165         double max_edep = 0;
0166         auto   leading_hit = (*simhits)[ixs[0]];
0167         // sum energy, take time from the most energetic hit
0168         for (size_t i = 0; i < ixs.size(); ++i) {
0169             auto hit = (*simhits)[ixs[i]];
0170
0171             double timeC = std::numeric_limits<double>::max();
0172             for (const auto& c : hit.getContributions()) {
0173                 if (c.getTime() <= timeC) {
0174                     timeC = c.getTime();
0175                 }
0176             }
0177             if (timeC > m_cfg.capTime) continue;
0178             edep += hit.getEnergy();
0179             trace("adding {} \t total: {}", hit.getEnergy(), edep);
0180
0181             // change maximum hit energy & time if necessary
0182             if (hit.getEnergy() > max_edep) {
0183                 max_edep = hit.getEnergy();
0184                 leading_hit = hit;
0185                 if (timeC <= time) {
0186                     time = timeC;
0187                 }
0188             }
0189
0190 #if EDM4EIC_VERSION_MAJOR >= 7
0191             edm4eic::MutableMCRecoCalorimeterHitAssociation assoc;
0192             assoc.setRawHit(rawhit);
0193             assoc.setSimHit(hit);
0194             assoc.setWeight(hit.getEnergy());
0195             rawassocs_staging.push_back(assoc);
0196 #endif
0197         }
0198         if (time > m_cfg.capTime) continue;
0199
0200         // safety check
0201         const double eResRel = (edep > m_cfg.threshold)
0202                 ? m_gaussian(m_generator) * std::sqrt(
0203                      std::pow(m_cfg.eRes[0] / std::sqrt(edep), 2) +
0204                      std::pow(m_cfg.eRes[1], 2) +
0205                      std::pow(m_cfg.eRes[2] / (edep), 2)
0206                   )
0207                 : 0;
0208         double corrMeanScale_value = corrMeanScale(leading_hit);
0209
0210         double ped = m_cfg.pedMeanADC + m_gaussian(m_generator) * m_cfg.pedSigmaADC;
0211
0212         // Note: both adc and tdc values must be positive numbers to avoid integer wraparound
0213         unsigned long long adc = std::max(std::llround(ped + edep * corrMeanScale_value * (1.0 + eResRel) / m_cfg.dyRangeADC * m_cfg.capADC), 0LL);
0214         unsigned long long tdc = std::llround((time + m_gaussian(m_generator) * tRes) * stepTDC);
0215
0216         if (edep> 1.e-3) trace("E sim {} \t adc: {} \t time: {}\t maxtime: {} \t tdc: {} \t corrMeanScale: {}", edep, adc, time, m_cfg.capTime, tdc, corrMeanScale_value);
0217
0218         rawhit.setCellID(leading_hit.getCellID());
0219         rawhit.setAmplitude(adc > m_cfg.capADC ? m_cfg.capADC : adc);
0220         rawhit.setTimeStamp(tdc);
0221         rawhits->push_back(rawhit);
0222
0223 #if EDM4EIC_VERSION_MAJOR >= 7
0224         for (auto& assoc : rawassocs_staging) {
0225             assoc.setWeight(assoc.getWeight() / edep);
0226             rawassocs->push_back(assoc);
0227         }
0228 #endif
0229     }
0230 }
0231
0232 } // namespace eicrecon