algorithms/digi/SiliconTrackerDigi.cc

0001 // SPDX-License-Identifier: LGPL-3.0-or-later
0002 // Copyright (C) 2022 Whitney Armstrong, Wouter Deconinck, Sylvester Joosten, Dmitry Romanov
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     // Create random gauss function
0025     m_gauss = [&](){
0026         return m_random.Gaus(0, m_cfg.timeResolution);
0027         //return m_rng.gaussian<double>(0., m_cfg.timeResolution);
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     // A map of unique cellIDs with temporary structure RawHit
0040     std::unordered_map<std::uint64_t, edm4eic::MutableRawTrackerHit> cell_hit_map;
0041
0042
0043     for (const auto& sim_hit : *sim_hits) {
0044
0045         // time smearing
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             // This cell doesn't have hits
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  // ns->ps
0073             };
0074         } else {
0075             // There is previous values in the cell
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             // keep earliest time for hit
0080             auto time_stamp = hit.getTimeStamp();
0081             hit.setTimeStamp(std::min(hit_time_stamp, hit.getTimeStamp()));
0082
0083             // sum deposited energy
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             // set association
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 } // namespace eicrecon