File indexing completed on 2024-09-27 07:02:58
0001
0002
0003
0004 #include "SiliconTrackerDigi.h"
0005
0006 #include <Evaluator/DD4hepUnits.h>
0007 #include <edm4eic/EDM4eicVersion.h>
0008 #include <edm4hep/MCParticleCollection.h>
0009 #include <edm4hep/Vector3d.h>
0010 #include <edm4hep/Vector3f.h>
0011 #include <fmt/core.h>
0012 #include <algorithm>
0013 #include <cmath>
0014 #include <cstdint>
0015 #include <gsl/pointers>
0016 #include <unordered_map>
0017 #include <utility>
0018
0019 #include "algorithms/digi/SiliconTrackerDigiConfig.h"
0020
0021 namespace eicrecon {
0022
0023 void SiliconTrackerDigi::init() {
0024
0025 m_gauss = [&](){
0026 return m_random.Gaus(0, m_cfg.timeResolution);
0027
0028 };
0029 }
0030
0031
0032 void SiliconTrackerDigi::process(
0033 const SiliconTrackerDigi::Input& input,
0034 const SiliconTrackerDigi::Output& output) const {
0035
0036 const auto [sim_hits] = input;
0037 auto [raw_hits,associations] = output;
0038
0039
0040 std::unordered_map<std::uint64_t, edm4eic::MutableRawTrackerHit> cell_hit_map;
0041
0042
0043 for (const auto& sim_hit : *sim_hits) {
0044
0045
0046 double time_smearing = m_gauss();
0047 double result_time = sim_hit.getTime() + time_smearing;
0048 auto hit_time_stamp = (std::int32_t) (result_time * 1e3);
0049
0050 debug("--------------------");
0051 debug("Hit cellID = {}", sim_hit.getCellID());
0052 debug(" position = ({:.2f}, {:.2f}, {:.2f})", sim_hit.getPosition().x, sim_hit.getPosition().y, sim_hit.getPosition().z);
0053 debug(" xy_radius = {:.2f}", std::hypot(sim_hit.getPosition().x, sim_hit.getPosition().y));
0054 debug(" momentum = ({:.2f}, {:.2f}, {:.2f})", sim_hit.getMomentum().x, sim_hit.getMomentum().y, sim_hit.getMomentum().z);
0055 debug(" edep = {:.2f}", sim_hit.getEDep());
0056 debug(" time = {:.4f}[ns]", sim_hit.getTime());
0057 debug(" particle time = {}[ns]", sim_hit.getMCParticle().getTime());
0058 debug(" time smearing: {:.4f}, resulting time = {:.4f} [ns]", time_smearing, result_time);
0059 debug(" hit_time_stamp: {} [~ps]", hit_time_stamp);
0060
0061
0062 if (sim_hit.getEDep() < m_cfg.threshold) {
0063 debug(" edep is below threshold of {:.2f} [keV]", m_cfg.threshold / dd4hep::keV);
0064 continue;
0065 }
0066
0067 if (cell_hit_map.count(sim_hit.getCellID()) == 0) {
0068
0069 cell_hit_map[sim_hit.getCellID()] = {
0070 sim_hit.getCellID(),
0071 (std::int32_t) std::llround(sim_hit.getEDep() * 1e6),
0072 hit_time_stamp
0073 };
0074 } else {
0075
0076 auto& hit = cell_hit_map[sim_hit.getCellID()];
0077 debug(" Hit already exists in cell ID={}, prev. hit time: {}", sim_hit.getCellID(), hit.getTimeStamp());
0078
0079
0080 auto time_stamp = hit.getTimeStamp();
0081 hit.setTimeStamp(std::min(hit_time_stamp, hit.getTimeStamp()));
0082
0083
0084 auto charge = hit.getCharge();
0085 hit.setCharge(charge + (std::int32_t) std::llround(sim_hit.getEDep() * 1e6));
0086 }
0087 }
0088
0089 for (auto item : cell_hit_map) {
0090 raw_hits->push_back(item.second);
0091
0092 for (const auto& sim_hit : *sim_hits) {
0093 if (item.first == sim_hit.getCellID()) {
0094
0095 auto hitassoc = associations->create();
0096 hitassoc.setWeight(1.0);
0097 hitassoc.setRawHit(item.second);
0098 #if EDM4EIC_VERSION_MAJOR >= 6
0099 hitassoc.setSimHit(sim_hit);
0100 #else
0101 hitassoc.addToSimHits(sim_hit);
0102 #endif
0103 }
0104 }
0105
0106 }
0107 }
0108
0109 }