algorithms/fardetectors/FarDetectorTrackerCluster.cc

0001 // SPDX-License-Identifier: LGPL-3.0-or-later
0002 // Copyright (C) 2023 - 2024, Simon Gardner
0003
0004 #include <DD4hep/Handle.h>
0005 #include <DD4hep/IDDescriptor.h>
0006 #include <DD4hep/Objects.h>
0007 #include <DD4hep/Readout.h>
0008 #include <DD4hep/VolumeManager.h>
0009 #include <DD4hep/detail/SegmentationsInterna.h>
0010 #include <DDSegmentation/BitFieldCoder.h>
0011 #include <JANA/JException.h>
0012 #include <Math/GenVector/Cartesian3D.h>
0013 #include <Math/GenVector/DisplacementVector3D.h>
0014 #include <ROOT/RVec.hxx>
0015 #include <algorithms/geo.h>
0016 #include <edm4hep/Vector3d.h>
0017 #include <fmt/core.h>
0018 #include <sys/types.h>
0019 #include <gsl/pointers>
0020
0021 #include "algorithms/fardetectors/FarDetectorTrackerCluster.h"
0022 #include "algorithms/fardetectors/FarDetectorTrackerClusterConfig.h"
0023
0024 namespace eicrecon {
0025
0026   void FarDetectorTrackerCluster::init(std::shared_ptr<spdlog::logger>& log) {
0027
0028     m_log = log;
0029     m_detector = algorithms::GeoSvc::instance().detector();
0030     m_cellid_converter = algorithms::GeoSvc::instance().cellIDPositionConverter();
0031
0032     if (m_cfg.readout.empty()) {
0033       throw JException("Readout is empty");
0034     }
0035     try {
0036       m_seg    = m_detector->readout(m_cfg.readout).segmentation();
0037       m_id_dec = m_detector->readout(m_cfg.readout).idSpec().decoder();
0038       if (!m_cfg.x_field.empty()) {
0039         m_x_idx = m_id_dec->index(m_cfg.x_field);
0040         m_log->debug("Find layer field {}, index = {}",  m_cfg.x_field, m_x_idx);
0041       }
0042       if (!m_cfg.y_field.empty()) {
0043         m_y_idx = m_id_dec->index(m_cfg.y_field);
0044         m_log->debug("Find layer field {}, index = {}", m_cfg.y_field, m_y_idx);
0045       }
0046     } catch (...) {
0047       m_log->error("Failed to load ID decoder for {}", m_cfg.readout);
0048       throw JException("Failed to load ID decoder");
0049     }
0050
0051   }
0052
0053   void FarDetectorTrackerCluster::process(
0054       const FarDetectorTrackerCluster::Input& input,
0055       const FarDetectorTrackerCluster::Output& output) const {
0056
0057     const auto [inputHitsCollections] = input;
0058     auto [outputClustersCollection]  = output;
0059
0060     //Loop over input and output collections - Any collection should only contain hits from a single surface
0061     for(size_t i=0; i<inputHitsCollections.size(); i++){
0062       auto inputHits = inputHitsCollections[i];
0063       if(inputHits->size() == 0) continue;
0064       auto outputClusters = outputClustersCollection[i];
0065
0066       // Make clusters
0067       auto clusters = ClusterHits(*inputHits);
0068
0069       // Create TrackerHits from 2D cluster positions
0070       ConvertClusters(clusters, *outputClusters);
0071
0072     }
0073   }
0074
0075   // Create vector of FDTrackerCluster from list of hits
0076   std::vector<FDTrackerCluster> FarDetectorTrackerCluster::ClusterHits(const edm4eic::RawTrackerHitCollection& inputHits) const {
0077
0078     std::vector<FDTrackerCluster> clusters;
0079
0080     ROOT::VecOps::RVec<unsigned long>  id;
0081     ROOT::VecOps::RVec<int>   x;
0082     ROOT::VecOps::RVec<int>   y;
0083     ROOT::VecOps::RVec<float> e;
0084     ROOT::VecOps::RVec<float> t;
0085
0086     // Gather detector id positions
0087     for(const auto& hit: inputHits){
0088       auto cellID = hit.getCellID();
0089       id.push_back(cellID);
0090       x.push_back (m_id_dec->get( cellID, m_x_idx      ));
0091       y.push_back (m_id_dec->get( cellID, m_y_idx      ));
0092       e.push_back (hit.getCharge());
0093       t.push_back (hit.getTimeStamp());
0094     }
0095
0096     // Set up clustering variables
0097     ROOT::VecOps::RVec<bool> available(id.size(), 1);
0098     auto indices = Enumerate(id);
0099
0100     // Loop while there are unclustered hits
0101     while(ROOT::VecOps::Any(available)){
0102
0103       dd4hep::Position localPos = {0,0,0};
0104       float  weightSum = 0;
0105
0106       float esum   = 0;
0107       float t0     = 0;
0108       float tError = 0;
0109       auto maxIndex = ROOT::VecOps::ArgMax(e*available);
0110
0111       available[maxIndex] = 0;
0112
0113       ROOT::VecOps::RVec<unsigned long> clusterList = {maxIndex};
0114       ROOT::VecOps::RVec<float> clusterT;
0115       std::vector<podio::ObjectID> clusterHits;
0116
0117       // Loop over hits, adding neighbouring hits as relevant
0118       while(clusterList.size()){
0119
0120         // Takes first remaining hit in cluster list
0121         auto index  = clusterList[0];
0122
0123         // Finds neighbours of cluster within time limit
0124         auto filter = available*(abs(x-x[index])<=1)*(abs(y-y[index])<=1)*(abs(t-t[index])<m_cfg.hit_time_limit);
0125
0126         // Adds the found hits to the cluster
0127         clusterList = Concatenate(clusterList,indices[filter]);
0128
0129         // Removes the found hits from the list of still available hits
0130         available = available*(!filter);
0131
0132         // Removes current hit from remaining found cluster hits
0133         clusterList.erase(clusterList.begin());
0134
0135         // Adds raw hit to TrackerHit contribution
0136         clusterHits.push_back((inputHits)[index].getObjectID());
0137
0138         // Energy
0139         auto hitE = e[index];
0140         esum += hitE;
0141         // TODO - See if now a single detector element is expected a better function is available.
0142         auto pos = m_seg->position(id[index]);
0143
0144         //Weighted position
0145         float weight = hitE; // TODO - Calculate appropriate weighting based on sensor charge sharing
0146         weightSum += weight;
0147         localPos  += pos*weight;
0148
0149         //Time
0150         clusterT.push_back(t[index]);
0151
0152       }
0153
0154       // Finalise position
0155       localPos/=weightSum;
0156
0157       // Finalise time
0158       t0      = Mean(clusterT);
0159       tError  = StdDev(clusterT); // TODO fold detector timing resolution into error
0160
0161       // Create cluster
0162       clusters.push_back(FDTrackerCluster{
0163         .cellID    = id[maxIndex],
0164         .x         = localPos.x(),
0165         .y         = localPos.y(),
0166         .energy    = esum,
0167         .time      = t0,
0168         .timeError = tError,
0169         .rawHits   = clusterHits
0170       });
0171
0172
0173     }
0174
0175     return clusters;
0176
0177   }
0178
0179   // Convert to global coordinates and create TrackerHits
0180   void FarDetectorTrackerCluster::ConvertClusters(const std::vector<FDTrackerCluster>& clusters, edm4hep::TrackerHitCollection& outputClusters) const {
0181
0182     // Get context of first hit
0183     const dd4hep::VolumeManagerContext* context = m_cellid_converter->findContext(clusters[0].cellID);
0184
0185     for(auto cluster: clusters){
0186       auto hitPos = outputClusters.create();
0187
0188       auto globalPos = context->localToWorld({cluster.x,cluster.y,0});
0189
0190       // Set cluster members
0191       hitPos.setCellID  (cluster.cellID);
0192       hitPos.setPosition(edm4hep::Vector3d(globalPos.x()/dd4hep::mm,globalPos.y()/dd4hep::mm,globalPos.z()/dd4hep::mm));
0193       hitPos.setEDep    (cluster.energy);
0194       hitPos.setTime    (cluster.time);
0195
0196       // Add raw hits to cluster
0197       for(auto hit: cluster.rawHits){
0198         hitPos.addToRawHits(hit);
0199       }
0200
0201     }
0202
0203   }
0204
0205
0206 }