File indexing completed on 2024-09-28 07:02:23
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/GaudiTool.h"
0018 #include "GaudiAlg/Transformer.h"
0019 #include "GaudiKernel/PhysicalConstants.h"
0020 #include "Gaudi/Property.h"
0021 #include "GaudiKernel/RndmGenerators.h"
0022
0023 #include "DDRec/CellIDPositionConverter.h"
0024 #include "DDSegmentation/BitFieldCoder.h"
0025
0026 #include "JugBase/IGeoSvc.h"
0027 #include "JugBase/DataHandle.h"
0028
0029 #include "fmt/format.h"
0030 #include "fmt/ranges.h"
0031
0032
0033 #include "edm4hep/SimCalorimeterHitCollection.h"
0034 #include "edm4eic/RawCalorimeterHitCollection.h"
0035 #include "edm4eic/RawCalorimeterHitData.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<edm4eic::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->detector()->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, 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 edm4eic::RawCalorimeterHit rawhit(
0198 ahit.getCellID(),
0199 (adc > m_capADC.value() ? m_capADC.value() : adc),
0200 tdc
0201 );
0202 rawhits->push_back(rawhit);
0203 }
0204 }
0205
0206 void signal_sum_digi() {
0207 const auto* const simhits = m_inputHitCollection.get();
0208 auto* rawhits = m_outputHitCollection.createAndPut();
0209
0210
0211 std::unordered_map<long long, std::vector<edm4hep::SimCalorimeterHit>> merge_map;
0212 for (const auto &ahit : *simhits) {
0213 int64_t hid = (ahit.getCellID() & id_mask) | ref_mask;
0214 auto it = merge_map.find(hid);
0215
0216 if (it == merge_map.end()) {
0217 merge_map[hid] = {ahit};
0218 } else {
0219 it->second.push_back(ahit);
0220 }
0221 }
0222
0223
0224 for (auto &[id, hits] : merge_map) {
0225 double edep = hits[0].getEnergy();
0226 double time = hits[0].getContributions(0).getTime();
0227 double max_edep = hits[0].getEnergy();
0228
0229 for (size_t i = 1; i < hits.size(); ++i) {
0230 edep += hits[i].getEnergy();
0231 if (hits[i].getEnergy() > max_edep) {
0232 max_edep = hits[i].getEnergy();
0233 for (const auto& c : hits[i].getContributions()) {
0234 if (c.getTime() <= time) {
0235 time = c.getTime();
0236 }
0237 }
0238 }
0239 }
0240
0241
0242 const double eResRel = (edep > m_threshold)
0243 ? m_normDist() * eRes[0] / std::sqrt(edep) +
0244 m_normDist() * eRes[1] +
0245 m_normDist() * eRes[2] / edep
0246 : 0;
0247
0248 double ped = m_pedMeanADC + m_normDist() * m_pedSigmaADC;
0249 unsigned long long adc = std::llround(ped + edep * (1. + eResRel) / dyRangeADC * m_capADC);
0250 unsigned long long tdc = std::llround((time + m_normDist() * tRes) * stepTDC);
0251
0252 edm4eic::RawCalorimeterHit rawhit(
0253 id,
0254 (adc > m_capADC.value() ? m_capADC.value() : adc),
0255 tdc
0256 );
0257 rawhits->push_back(rawhit);
0258 }
0259 }
0260 };
0261
0262 DECLARE_COMPONENT(CalorimeterHitDigi)
0263
0264 }