EIC code displayed by LXR master/​EICrecon/​src/​algorithms/​digi/​MPGDTrackerDigi.cc	Source navigation Diff markup Identifier search General search
	Version: master test Revert   Change

File indexing completed on 2025-02-22 09:38:59

 // SPDX-License-Identifier: LGPL-3.0-or-later
 // Copyright (C) 2022 - 2025 Whitney Armstrong, Wouter Deconinck, Sylvester Joosten, Dmitry Romanov, Yann Bedfer
 
 /*
   Digitization specific to MPGDs.
   - What's special in MPGDs is their 2D-strip readout: i.e. simultaneous
    registering along two sets of coordinate strips.
   - The "process" method involves a combination of segmentation, simulation and
    digitization.
   - The segmentation done when producing "sim_hits", and stored as a pixel
    cellID, is overwritten by strip cellIDs, via "dd4hep::MultiSegmentation".
   - The simulation will eventually involve simulating the amplitude and timing
    correlation of the two coordinates, and the spreading of the charge on
    adjacent strips producing multi-hit clusters.
     The preliminary version is rudimentary: single-hit clusters, with identical
    timing and uncorrelated amplitudes.
   - The digitization follows the standard steps, only that it needs to be
    convoluted with the simulation.
   NOTA BENE: The file could be simplified, see issue #1702.
  */
 
 
 #include "MPGDTrackerDigi.h"
 
 #include <DD4hep/Alignments.h>
 #include <DD4hep/DetElement.h>
 #include <DD4hep/Handle.h>
 #include <DD4hep/IDDescriptor.h>
 #include <DD4hep/Objects.h>
 #include <DD4hep/Readout.h>
 #include <DD4hep/VolumeManager.h>
 #include <DD4hep/config.h>
 #include <DD4hep/detail/SegmentationsInterna.h>
 #include <DDSegmentation/BitFieldCoder.h>
 #include <Evaluator/DD4hepUnits.h>
 #include <JANA/JException.h>
 #include <Math/GenVector/Cartesian3D.h>
 #include <Math/GenVector/DisplacementVector3D.h>
 #include <Parsers/Primitives.h>
 // Access "algorithms:GeoSvc"
 #include <algorithms/geo.h>
 #include <algorithms/logger.h>
 #include <edm4hep/EDM4hepVersion.h>
 #include <edm4eic/EDM4eicVersion.h>
 #include <edm4hep/MCParticleCollection.h>
 #include <edm4hep/Vector3d.h>
 #include <edm4hep/Vector3f.h>
 #include <fmt/core.h>
 #include <algorithm>
 #include <cmath>
 #include <cstdint>
 #include <gsl/pointers>
 #include <initializer_list>
 #include <unordered_map>
 #include <utility>
 #include <vector>
 
 #include "algorithms/digi/MPGDTrackerDigiConfig.h"
 
 using namespace dd4hep;
 
 namespace eicrecon {
 
 void MPGDTrackerDigi::init() {
     // Create random gauss function
     m_gauss = [&](){
         return m_random.Gaus(0, m_cfg.timeResolution);
         //return m_rng.gaussian<double>(0., m_cfg.timeResolution);
     };
 
     // Access id decoder
     m_detector = algorithms::GeoSvc::instance().detector();
     const dd4hep::BitFieldCoder* m_id_dec;
     if (m_cfg.readout.empty()) {
         throw JException("Readout is empty");
     }
     try {
         m_seg =    m_detector->readout(m_cfg.readout).segmentation();
         m_id_dec = m_detector->readout(m_cfg.readout).idSpec().decoder();
     } catch (...) {
         critical("Failed to load ID decoder for \"{}\" readout.", m_cfg.readout);
         throw JException("Failed to load ID decoder");
     }
     // Method "process" relies on a strict assumption on the IDDescriptor:
     // - Must have a "strip" field.
     // - That "strip" field includes bits 30|31.
     // Let's check.
     if (m_id_dec->get(((CellID)0x3)<<30,"strip") != 0x3) {
         critical("Missing or invalid \"strip\" field in IDDescriptor for \"{}\" readout.",
                  m_cfg.readout);
         throw JException("Invalid IDDescriptor");
     }
     debug("Find valid \"strip\" field in IDDescriptor for \"{}\" readout.",
           m_cfg.readout);
 }
 
 
 void MPGDTrackerDigi::process(
         const MPGDTrackerDigi::Input& input,
         const MPGDTrackerDigi::Output& output) const {
 
     // ********** SIMULATE THE 2D-strip READOUT of MPGDs.
     // - Overwrite and extend segmentation stored in "sim_hit", which is anyway
     //  expected to be along a single coordinate (this happens to allow one to
     //  reconstruct data w/ a segmentation differing from that used when
     //  generating the data).
     // - New segmentation is along two coordinates, described by two cellID's
     //  with each a distinctive "strip" field.
     //   N.B.: Assumptions on the IDDescriptor: the "strip" specification
     //  is fixed = cellID>>32&0x3.
     // - The simulation is simplistic: single-hit cluster per coordinate.
 
     const auto [sim_hits] = input;
     auto [raw_hits,associations] = output;
 
     // A map of unique cellIDs with temporary structure RawHit
     std::unordered_map<std::uint64_t, edm4eic::MutableRawTrackerHit> cell_hit_map;
     // Prepare for strip segmentation
     const Position dummy(0,0,0);
     const VolumeManager& volman = m_detector->volumeManager();
 
     using CellIDs = std::pair<CellID,CellID>;
     using Sim2IDs = std::vector<CellIDs>; Sim2IDs sim2IDs;
     for (const edm4hep::SimTrackerHit& sim_hit : *sim_hits) {
 
         // ***** TIME SMEARING
         // - Simplistic treatment.
         // - A more realistic one would have to distinguish a smearing common to
         //  both coordinates of the 2D-strip readout (due to the drifting of the
         //  leading primary electrons) from other smearing effects, specific to
         //  each coordinate.
         double time_smearing = m_gauss();
         double result_time = sim_hit.getTime() + time_smearing;
         auto hit_time_stamp = (std::int32_t) (result_time * 1e3);
 
         // ***** SEGMENTATION
         // - The two cellID's are encoded via a "dd4hep::MultiSegmentation"
         //  discriminating on the strip field, w/ "strip" setting of 0x1 (
         //  called 'p') and 0x2 (called 'n').
         // - They are evaluated based on "sim_hit" Cartesian coordinates
         //  positions
         //   Given that all the segmentation classes foreseen for MPGDs (
         //  "CartesianGrid.." for Outer and EndCaps, "CylindricalGridPhiZ" for
         //  "CyMBaL") disregard the _global_ position argument to
         //  "dd4hep::Segmentation::cellID", we need the _local_ position and
         //  only that.
         const edm4hep::Vector3d &pos = sim_hit.getPosition();
         using dd4hep::mm;
         Position gpos(pos.x*mm,pos.y*mm,pos.z*mm);
         CellID vID = // Note: Only the bits corresponding to the volumeID will
           // be used. The rest, encoding the segmentation stored in "sim_hit",
           // being disregared.
           sim_hit.getCellID();
         DetElement local = volman.lookupDetElement(vID);
         const auto lpos = local.nominal().worldToLocal(gpos);
         // p "strip"
         CellID stripBitp = ((CellID)0x1)<<30, vIDp = vID|stripBitp;
         CellID cIDp = m_seg->cellID(lpos,dummy,vIDp);
         // n "strip"
         CellID stripBitn = ((CellID)0x2)<<30, vIDn = vID|stripBitn;
         CellID cIDn = m_seg->cellID(lpos,dummy,vIDn);
 
         sim2IDs.push_back({cIDp,cIDn}); // Remember cellIDs.
         // ***** DEBUGGING INFO
         if(level() >= algorithms::LogLevel::kDebug) {
             CellID hIDp = cIDp>>32, sIDp = cIDp>>30&0x3;
             debug("--------------------");
             debug("Hit cellIDp  = 0x{:08x}, 0x{:08x} 0x{:02x}", hIDp, vIDp, sIDp);
             CellID hIDn = cIDn>>32, sIDn = cIDn>>30&0x3;
             debug("Hit cellIDn  = 0x{:08x}, 0x{:08x} 0x{:02x}", hIDn, vIDn, sIDn);
             debug("   position  = ({:.2f}, {:.2f}, {:.2f})", sim_hit.getPosition().x, sim_hit.getPosition().y, sim_hit.getPosition().z);
             debug("   xy_radius = {:.2f}", std::hypot(sim_hit.getPosition().x, sim_hit.getPosition().y));
             debug("   momentum  = ({:.2f}, {:.2f}, {:.2f})", sim_hit.getMomentum().x, sim_hit.getMomentum().y, sim_hit.getMomentum().z);
             debug("   edep = {:.2f}", sim_hit.getEDep());
             debug("   time = {:.4f}[ns]", sim_hit.getTime());
 #if EDM4HEP_BUILD_VERSION >= EDM4HEP_VERSION(0, 99, 0)
             debug("   particle time = {}[ns]", sim_hit.getParticle().getTime());
 #else
             debug("   particle time = {}[ns]", sim_hit.getMCParticle().getTime());
 #endif
             debug("   time smearing: {:.4f}, resulting time = {:.4f} [ns]", time_smearing, result_time);
             debug("   hit_time_stamp: {} [~ps]", hit_time_stamp);
         }
 
         // ***** APPLY THRESHOLD
         if (sim_hit.getEDep() < m_cfg.threshold) {
             debug("  edep is below threshold of {:.2f} [keV]", m_cfg.threshold / keV);
             continue;
         }
 
         // ***** HIT ACCUMULATION
         for (CellID cID : {cIDp,cIDn}) {
             if (cell_hit_map.count(cID) == 0) {
                 // This cell doesn't have hits
                 cell_hit_map[cID] = {
                     cID,
                     (std::int32_t) std::llround(sim_hit.getEDep() * 1e6),
                     hit_time_stamp  // ns->ps
                 };
             } else {
                 // There is previous values in the cell
                 auto& hit = cell_hit_map[cID];
                 debug("  Hit already exists in cell ID={}, prev. hit time: {}", cID, hit.getTimeStamp());
 
                 // keep earliest time for hit
                 auto time_stamp = hit.getTimeStamp();
                 hit.setTimeStamp(std::min(hit_time_stamp, hit.getTimeStamp()));
 
                 // sum deposited energy
                 auto charge = hit.getCharge();
                 hit.setCharge(charge + (std::int32_t) std::llround(sim_hit.getEDep() * 1e6));
             }
         }
     }
 
     // ***** raw_hit INSTANTIATION AND raw<-sim_hit's ASSOCIATION
     for (auto item : cell_hit_map) {
         raw_hits->push_back(item.second);
         Sim2IDs::const_iterator sim_it = sim2IDs.cbegin();
         for (const auto& sim_hit : *sim_hits) {
             CellIDs cIDs = *sim_it++;
             for (CellID cID : {cIDs.first,cIDs.second}) {
                 if (item.first == cID) {
                     // set association
                     auto hitassoc = associations->create();
                     hitassoc.setWeight(1.0);
                     hitassoc.setRawHit(item.second);
 #if EDM4EIC_VERSION_MAJOR >= 6
                     hitassoc.setSimHit(sim_hit);
 #else
                     hitassoc.addToSimHits(sim_hit);
 #endif
                 }
             }
         }
 
     }
 }
 
 } // namespace eicrecon

This page was automatically generated by the 2.3.7 LXR engine. The LXR team