File indexing completed on 2024-09-27 07:03:48
0001
0002
0003
0004
0005
0006
0007
0008
0009
0010
0011
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
0034 #include "edm4hep/SimCalorimeterHitCollection.h"
0035 #include "edm4hep/RawCalorimeterHitCollection.h"
0036
0037
0038 using namespace Gaudi::Units;
0039
0040 namespace Jug::Digi {
0041
0042
0043
0044
0045
0046
0047 class CalorimeterHitDigi : public GaudiAlgorithm {
0048 public:
0049
0050 Gaudi::Property<std::vector<double>> u_eRes{this, "energyResolutions", {}};
0051 Gaudi::Property<double> m_tRes{this, "timeResolution", 0.0 * ns};
0052
0053 Gaudi::Property<double> m_threshold{this, "threshold", 1. * keV};
0054
0055
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
0064
0065
0066
0067
0068
0069
0070
0071 Gaudi::Property<std::vector<std::string>> u_fields{this, "signalSumFields", {}};
0072
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
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
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
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
0107 for (size_t i = 0; i < u_eRes.size() && i < 3; ++i) {
0108 eRes[i] = u_eRes[i];
0109 }
0110
0111
0112 dyRangeADC = m_dyRangeADC.value() / GeV;
0113 tRes = m_tRes.value() / ns;
0114 stepTDC = ns / m_resolutionTDC.value();
0115
0116
0117 if (!u_fields.value().empty()) {
0118 m_geoSvc = service(m_geoSvcName);
0119
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
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
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
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
0170 const auto* const simhits = m_inputHitCollection.get();
0171
0172 auto* rawhits = m_outputHitCollection.createAndPut();
0173 for (const auto& ahit : *simhits) {
0174
0175 const double eDep = ahit.getEnergy();
0176
0177
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
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
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
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
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
0260 DECLARE_COMPONENT(CalorimeterHitDigi)
0261
0262 }