FPGATrackSimClusteringTool.cxx

Go to the documentation of this file.

/*
   Copyright (C) 2002-2025 CERN for the benefit of the ATLAS collaboration
   */
 
#include "FPGATrackSimClusteringTool.h"
#include "AthenaBaseComps/AthMsgStreamMacros.h"
#include "FPGATrackSimObjects/FPGATrackSimMultiTruth.h"
#include "FPGATrackSimObjects/FPGATrackSimConstants.h"
#include "CxxUtils/trapping_fp.h"
#include <algorithm>
#include <cmath>
 
 
namespace{
    //For deciding eta || phi columns in modules
    constexpr unsigned int ETA = 1;
    constexpr unsigned int PHI = 0;
}
 
FPGATrackSimClusteringTool::FPGATrackSimClusteringTool(const std::string& algname, const std::string &name, const IInterface *ifc) :
    base_class(algname, name, ifc)
{
}
 
StatusCode FPGATrackSimClusteringTool::initialize() {
    ATH_CHECK(m_lorentzAngleTool.retrieve(EnableTool{m_LorentzAngleShift >=0}));
    return StatusCode::SUCCESS;
}
 
 
StatusCode FPGATrackSimClusteringTool::DoClustering(FPGATrackSimLogicalEventInputHeader &header, std::vector<FPGATrackSimCluster> &clusters) const
{
    for (int i = 0; i<header.nTowers(); i++)
    {
        // Retreive the hits from the tower
        FPGATrackSimTowerInputHeader& tower = *header.getTower(i);
        HitPtrCollection hits;
        hits.reserve(tower.hits().size());
        for (auto& hit : tower.hits()) {
            hits.push_back(std::make_unique<FPGATrackSimHit>(hit));
        }
 
        HitPtrContainer hitsPerModule;
        std::vector<FPGATrackSimCluster> towerClusters;
 
        if (m_reduceCoordPrecision) {
            for (auto &hit : hits)
                reduceGlobalCoordPrecision(*hit);
        }
 
        splitAndSortHits(std::move(hits), hitsPerModule);
        SortedClustering(std::move(hitsPerModule), towerClusters);
        normaliseClusters(towerClusters);
 
        //remove the old hits from the tower...
        tower.clearHits();
        tower.reserveHits(towerClusters.size());
        clusters.clear();
        clusters.reserve(towerClusters.size());
        if(i > 1)
            ATH_MSG_WARNING("more than one tower, m_clusters is only going to contain those from the last one");
 
        unsigned cluster_count = 0;
        unsigned int pixelCounter =  0;
        unsigned int stripCounter = 0;
        for ( auto &cluster: towerClusters){
            SetMinMaxIndicies(cluster);
            if (m_reduceCoordPrecision)
                reduceGlobalCoordPrecision(cluster);
 
            FPGATrackSimHit cluster_as_FPGATrackSimhit = cluster.getClusterEquiv();
            cluster_as_FPGATrackSimhit.setHitType(HitType::clustered);
            cluster_as_FPGATrackSimhit.setParentageMask(cluster_count); // making use of unused m_parentageMask to keep track of cluster index
 
            if(cluster_as_FPGATrackSimhit.getDetType() == SiliconTech::pixel)
            {
              cluster_as_FPGATrackSimhit.setCluster1ID(pixelCounter);
              pixelCounter++;
            }
            else if(cluster_as_FPGATrackSimhit.getDetType() == SiliconTech::strip)
            {
              cluster_as_FPGATrackSimhit.setCluster1ID(stripCounter);
              stripCounter++;
            }
            tower.addHit(cluster_as_FPGATrackSimhit);
 
            if(m_LorentzAngleShift>= 0){
                ATH_CHECK(m_lorentzAngleTool->updateHitPosition(cluster_as_FPGATrackSimhit, m_LorentzAngleShift));
            }
            cluster.setClusterEquiv(cluster_as_FPGATrackSimhit);
            //send back a copy for monitoring and to check when writing out hits in each road
            clusters.push_back(cluster);
            cluster_count++;
        }
    }
    ATH_MSG_DEBUG("Produced "<< clusters.size()<< " clusters");
    return StatusCode::SUCCESS;
}
 
void FPGATrackSimClusteringTool::SetMinMaxIndicies(FPGATrackSimCluster &cluster) const {
    //Set the min and max indices for the cluster equivalent
    FPGATrackSimHit clusterEquiv = cluster.getClusterEquiv();
    int maxPhiIdx{-std::numeric_limits<int>::max()}, maxEtaIdx{-std::numeric_limits<int>::max()};
    int minPhiIdx{std::numeric_limits<int>::max()}, minEtaIdx{std::numeric_limits<int>::max()};
    for (const FPGATrackSimHit& hit : cluster.getHitList()){
         maxPhiIdx = std::max(maxPhiIdx, static_cast<int>(hit.getPhiIndex()));
         minPhiIdx = std::min(minPhiIdx, static_cast<int>(hit.getPhiIndex()));
         maxEtaIdx = std::max(maxEtaIdx, static_cast<int>(hit.getEtaIndex()));
         minEtaIdx = std::min(minEtaIdx, static_cast<int>(hit.getEtaIndex()));
    }
    clusterEquiv.setMinPhiIndex(minPhiIdx);
    clusterEquiv.setMaxPhiIndex(maxPhiIdx);
    clusterEquiv.setMinEtaIndex(minEtaIdx);
    clusterEquiv.setMaxEtaIndex(maxEtaIdx);
    cluster.setClusterEquiv(clusterEquiv);
}
 
//Attempt to implement clustering using FPGATrackSim objects.
void FPGATrackSimClusteringTool::SortedClustering(HitPtrContainer&& sorted_hits, std::vector<FPGATrackSimCluster> &clusters) const {
    std::vector<FPGATrackSimCluster> moduleClusters;
    //Loop over the sorted modules that we have
    for( HitPtrCollection & moduleHits: sorted_hits){
        //Make the clusters for this module
        Clustering(std::move(moduleHits), moduleClusters);
        //Put these clusters into the output list
        clusters.insert(clusters.end(), moduleClusters.begin(), moduleClusters.end());
        //Clear the vector or this will get messy
        moduleClusters.clear();
    }
}
 
 
void FPGATrackSimClusteringTool::Clustering(HitPtrCollection &&moduleHits, std::vector<FPGATrackSimCluster> &moduleClusters) const {
    FPGATrackSimHit clusterEquiv;
    bool newHit;
 
    //To hold the current cluster vars for comparison
    //loop over the hits that we have been passed for this module
    for( auto& hit: moduleHits){
        bool is_clustered_hit = false;
        std::vector<FPGATrackSimCluster>::iterator it_added_clus;
 
        //Loop over the clusters we have already made, check if this hit should be added to them?
        for(std::vector<FPGATrackSimCluster>::iterator it = moduleClusters.begin(); it != moduleClusters.end(); ++it) {
            if(hit->isPixel()) {
                if (FPGATrackSimCLUSTERING::updatePixelCluster(*it, *hit, false, m_digitalClustering)) {
                    if (!is_clustered_hit) {
                        is_clustered_hit = true;
                        it_added_clus = it;
                    } else {
                        int cPhi = it->getClusterEquiv().getPhiIndex();
                        int cPhiWidth = it->getClusterEquiv().getPhiWidth();
                        int cEta = it->getClusterEquiv().getEtaIndex();
                        int cEtaWidth = it->getClusterEquiv().getEtaWidth();
                        int fCPhi = it_added_clus->getClusterEquiv().getPhiIndex();
                        int fCPhiWidth = it_added_clus->getClusterEquiv().getPhiWidth();
                        int fCEta = it_added_clus->getClusterEquiv().getEtaIndex();
                        int fCEtaWidth = it_added_clus->getClusterEquiv().getEtaWidth();
 
                        clusterEquiv = it_added_clus->getClusterEquiv();
 
                        // set new phi & phi width
                        if (cPhi < fCPhi) {
                            clusterEquiv.setPhiIndex(cPhi);
                            if (cPhi + cPhiWidth < fCPhi + fCPhiWidth)
                                clusterEquiv.setPhiWidth(fCPhiWidth + (fCPhi - cPhi));
                            else
                                clusterEquiv.setPhiWidth(cPhiWidth);
                        } else {
                            clusterEquiv.setPhiIndex(fCPhi);
                            if (!(cPhi + cPhiWidth < fCPhi + fCPhiWidth))
                                clusterEquiv.setPhiWidth(cPhiWidth + (cPhi - fCPhi));
                            else
                                clusterEquiv.setPhiWidth(fCPhiWidth);
                        }
 
                        // set new eta & eta width
                        if (cEta < fCEta) {
                            clusterEquiv.setEtaIndex(cEta);
                            if (cEta + cEtaWidth < fCEta + fCEtaWidth)
                                clusterEquiv.setEtaWidth(fCEtaWidth + (fCEta - cEta));
                            else
                                clusterEquiv.setEtaWidth(cEtaWidth);
                        } else {
                            clusterEquiv.setEtaIndex(fCEta);
                            if (!(cEta + cEtaWidth < fCEta + fCEtaWidth))
                                clusterEquiv.setEtaWidth(cEtaWidth + (cEta - fCEta));
                            else
                                clusterEquiv.setEtaWidth(fCEtaWidth);
                        }
 
                        it_added_clus->setClusterEquiv(clusterEquiv);
 
                        for (auto& hit : it->getHitList()) {
                            newHit = true;
                            for (auto& finalHit : it_added_clus->getHitList()) {
                                if (hit.getEtaIndex() == finalHit.getEtaIndex() &&
                                        hit.getPhiIndex() == finalHit.getPhiIndex())
                                    newHit = false;
                            }
                            if (newHit) {
                                CXXUTILS_TRAPPING_FP;
                                clusterEquiv = it_added_clus->getClusterEquiv();
                                float xOld = clusterEquiv.getX();
                                float yOld = clusterEquiv.getY();
                                float zOld = clusterEquiv.getZ();
                                float xPhiOld = clusterEquiv.getPhiCoord();
                                float xEtaOld = clusterEquiv.getEtaCoord();
                                float cPhiOld = clusterEquiv.getCentroidPhiIndex();
                                float cEtaOld = clusterEquiv.getCentroidEtaIndex();
                                float xNew = hit.getX();
                                float yNew = hit.getY();
                                float zNew = hit.getZ();
                                float xPhiNew = hit.getPhiCoord();
                                float xEtaNew = hit.getEtaCoord();
                                float cPhiNew = hit.getPhiIndex();
                                float cEtaNew = hit.getEtaIndex();
                                int tot = clusterEquiv.getToT();
                                int totNew = hit.getToT();
                                if (m_digitalClustering) {
                                    // n+1 because that is old + new now
                                    int n = it_added_clus->getHitList().size();
                                    clusterEquiv.setX((xOld*n + xNew) / (n+1));
                                    clusterEquiv.setY((yOld*n + yNew) / (n+1));
                                    clusterEquiv.setZ((zOld*n + zNew) / (n+1));
                                    clusterEquiv.setPhiCoord((xPhiOld*n + xPhiNew) / (n+1));
                                    clusterEquiv.setEtaCoord((xEtaOld*n + xEtaNew) / (n+1));
                                    clusterEquiv.setCentroidPhiIndex((cPhiOld*n + cPhiNew) / (n+1));
                                    clusterEquiv.setCentroidEtaIndex((cEtaOld*n + cEtaNew) / (n+1));
                                } else {
                                    clusterEquiv.setX((xOld*tot + xNew*totNew) / (tot+totNew));
                                    clusterEquiv.setY((yOld*tot + yNew*totNew) / (tot+totNew));
                                    clusterEquiv.setZ((zOld*tot + zNew*totNew) / (tot+totNew));
                                    clusterEquiv.setPhiCoord((xPhiOld*tot + xPhiNew*totNew) / (tot+totNew));
                                    clusterEquiv.setEtaCoord((xEtaOld*tot + xEtaNew*totNew) / (tot+totNew));
                                    clusterEquiv.setCentroidPhiIndex((cPhiOld*tot + cPhiNew*totNew) / (tot+totNew));
                                    clusterEquiv.setCentroidEtaIndex((cEtaOld*tot + cEtaNew*totNew) / (tot+totNew));
                                }
                                clusterEquiv.setToT(tot + totNew);
                                it_added_clus->setClusterEquiv(clusterEquiv);
                                it_added_clus->push_backHitList(hit);
                            }
                        }
 
                        // move the last cluster to the newly freed spot in the cluster array if the merged cluster was not
                        // already the last cluster, decrease loop counter by one, so that this cluster gets also checked
                        // against the current hit
                        if (it != moduleClusters.end() - 1) {
                            *it = moduleClusters.back();
                            moduleClusters.pop_back();
                            it -= 1;
                        } else {
                            moduleClusters.pop_back();
                            break;
            }
                    }
                }
            }
            if(hit->isStrip()){
                if (FPGATrackSimCLUSTERING::updateStripCluster(*it, *hit, false, m_digitalClustering))
                    is_clustered_hit = true;
            }
        }
 
        //If it is the first hit or a not clustered hit, then start a new cluster and add it to the output vector
        if((!is_clustered_hit) || (moduleClusters.size() == 0)){
            FPGATrackSimCluster cluster;
            if(hit->isPixel()){
                // No need to check the return code here
                FPGATrackSimCLUSTERING::updatePixelCluster(cluster, *hit, true, m_digitalClustering);
            } else if(hit->isStrip()){
                FPGATrackSimCLUSTERING::updateStripCluster(cluster, *hit, true, m_digitalClustering);
            }
            //Put this cluster into the output hits. Will update it in place.
            moduleClusters.push_back(cluster);
        }
    }
    moduleHits.clear();
}
 
void FPGATrackSimClusteringTool::splitAndSortHits(HitPtrCollection &&hits, HitPtrContainer &hitsPerModule, int &eta_phi) const {
    splitHitsToModules(std::move(hits), hitsPerModule);
    sortHitsOnModules(hitsPerModule, eta_phi);
}
 
 
void FPGATrackSimClusteringTool::splitAndSortHits(HitPtrCollection &&hits, HitPtrContainer &hitsPerModule) const{
    splitHitsToModules(std::move(hits), hitsPerModule);
    sortHitsOnModules(hitsPerModule);
}
 
/*Temporarilly sort the hits into module by module packets
*/
void FPGATrackSimClusteringTool::splitHitsToModules(HitPtrCollection &&hits, HitPtrContainer &hitsPerModule) const{
    //To hold the current module
    HitPtrCollection currentModule;
    uint hashing = 0;
    //Split the incoming hits into hits by module
    for ( auto& hit:hits){
        if(hashing == 0){
            hashing = hit->getIdentifierHash();
            currentModule.push_back(std::move(hit));
        } else if (hit->getIdentifierHash() == hashing) {
            currentModule.push_back(std::move(hit));
        } else {
            hitsPerModule.push_back(std::exchange(currentModule, HitPtrCollection{}));
            hashing = hit->getIdentifierHash();
            currentModule.push_back(std::move(hit));
        }
    }
 
    // Now push that last one
    if (currentModule.size() > 0) hitsPerModule.push_back(std::move(currentModule));
}
 
void FPGATrackSimClusteringTool::sortHitsOnModules(HitPtrContainer &hitsPerModule, int &eta_phi) const{
    //Loop over the module separated hits
    for ( auto& module:hitsPerModule){
        //Work out if columns are ETA (1) || PHI (0)
        if(etaOrPhi(*module.at(0)) == true){
            //Sort by ETA first
            eta_phi = ETA;
            if (module.size() > 1) {
                if (module.at(0)->isStrip())
                    std::stable_sort(module.begin(), module.end(), FPGATrackSimCLUSTERING::sortITkInputEta);
            }
            if (module.size() > 1) {
                if (module.at(0)->isStrip())
                    std::stable_sort(module.begin(), module.end(), FPGATrackSimCLUSTERING::sortITkInputPhi);
            }
        } else {
            //Sort by PHI first
            eta_phi = PHI;
            if (module.size() > 1) {
                if (module.at(0)->isStrip())
                    std::stable_sort(module.begin(), module.end(), FPGATrackSimCLUSTERING::sortITkInputPhi);
            }
            if (module.size() > 1) {
                if (module.at(0)->isStrip())
                    std::stable_sort(module.begin(), module.end(), FPGATrackSimCLUSTERING::sortITkInputEta);
            }
        }
    }
}
 
void FPGATrackSimClusteringTool::sortHitsOnModules(HitPtrContainer &hitsPerModule) const{
    //Loop over the module separated hits
    for ( auto& module:hitsPerModule){
        if (module.size() > 1) {
            if (module.at(0)->isStrip())
                std::stable_sort(module.begin(), module.end(), FPGATrackSimCLUSTERING::sortITkInputEta);
        }
        if (module.size() > 1) {
            if (module.at(0)->isStrip())
                std::stable_sort(module.begin(), module.end(), FPGATrackSimCLUSTERING::sortITkInputPhi);
        }
    }
}
 
 
 
//Need to remove the fpgatracksim::scaleHitFactor and normalise the widths
void FPGATrackSimClusteringTool::normaliseClusters(std::vector<FPGATrackSimCluster> &clusters) const {
    for( auto &cluster:clusters){
        //Grab the cluster equiv
        FPGATrackSimHit clusterEquiv = cluster.getClusterEquiv();
        //Update the clusterEquiv's position and width
        if(clusterEquiv.isStrip()){
            //Clear the groupsize, set this to be one as we are only clustering one row.
            clusterEquiv.setEtaWidth(1);
            clusterEquiv.setCentroidPhiIndex(clusterEquiv.getCentroidPhiIndex()/fpgatracksim::scaleHitFactor);
            clusterEquiv.setPhiWidth(clusterEquiv.getPhiWidth()+1);
        } else {
            // Nothing to normalise for pixel clusters
            continue;
        }
        cluster.setClusterEquiv(clusterEquiv);
    }
}
 
 
/*Function to work out if we need to sort by eta or phi first.
 *Depends on the position and orientation of the module in the detector.
 * Currently backwards engineered from the MC. Need to ask ITk people to confirm.
 */
bool FPGATrackSimClusteringTool::etaOrPhi(const FPGATrackSimHit& hit) const {
 
    /*This currently might get complicated as eta/phi ordering varies depending on the position in the detector.
     * For ITK-20-00-00 Step 2.2 inclined duals layout, worked out by Naoki.
     * Detector position               |  Module Type (# chip)  |  Column  |  Row  |
     * =============================================================================
     * Barrel Layer 0 Eta module 1-6   |  2           (double)  |  eta     |  phi  |
     * Barrel Layer 0 Eta module 7-22  |  1           (single)  |  phi     |  eta  |
     * Barrel Layer 1 Eta module 1-6   |  3           (quad)    |  eta     |  phi  |
     * Barrel Layer 1 Eta module 7-19  |  3                     |  phi     |  eta  |
     * Barrel Layer 2 Eta module 1-11  |  3                     |  eta     |  phi  |
     * Barrel Layer 2 Eta module 12-22 |  2                     |  phi     |  eta  |
     * Barrel Layer 3 Eta module 1-12  |  3                     |  eta     |  phi  |
     * Barrel Layer 3 Eta module 13-25 |  2                     |  phi     |  eta  |
     * Barrel Layer 4 Eta module 1-13  |  3                     |  eta     |  phi  |
     * Barrel Layer 4 Eta module 14-26 |  2                     |  phi     |  eta  |
     * All Endcap modules              |  3                     |  eta     |  phi  |
     * =============================================================================
     * Module Type 1 = Single, 2, Dual, 3 Quad
     * 328x400 blocks
     * Hit type is essentially isPixel
     * DetectorZone 0 = Barrel, -ive/+ive = endcaps
     */
 
    //Check if the two hits are from the same module
    //If it is not a barrel module then sort eta as column
    if(hit.getDetectorZone() != DetectorZone::barrel){
        return ETA;
    }
    //Otherwise it is a barrel module and now things get more complicated
    else {
        //Start by looking at what layer it is in
        if (hit.getPhysLayer() == 0 || hit.getPhysLayer() == 1) {
            if (hit.getEtaModule() <=6) {
                return ETA;
            } else {
                return PHI;
            }
        } else if (hit.getPhysLayer() == 2) {
            if (hit.getEtaModule() <=11) {
                return ETA;
            } else {
                return PHI;
            }
        } else if (hit.getPhysLayer() == 3) {
            if (hit.getEtaModule() <=12) {
                return ETA;
            } else {
                return PHI;
            }
        } else if (hit.getPhysLayer() == 4) {
            if (hit.getEtaModule() <=13) {
                return ETA;
            } else {
                return PHI;
            }
        }
    }
    //Default to ETA, but shouldn't reach here
    return ETA;
}
 
 
 
void FPGATrackSimCLUSTERING::attachTruth(std::vector<FPGATrackSimHit> &hits){
    for( auto& hit : hits) {
        FPGATrackSimMultiTruth mt;
        // record highest pt contribution to the combination (cluster
        if(!hit.getTruth().isEmpty()) {
            mt.add(hit.getTruth());
            hit.setTruth(mt);
        } else {
            FPGATrackSimMultiTruth::Barcode uniquecode(hit.getEventIndex(), hit.getBarcode());
            mt.maximize(uniquecode, hit.getBarcodePt());
            hit.setTruth(mt);
        }
    }
} //record truth for each raw channel in the cluster
 
/*
 * This function is used in the FPGATrackSimClusteringTools to see if a new hit should be added to the current cluster under construction.
 * It checks if the hit is in a number of positions w.r.t. the cluster being formed: up/right, down/right, above, right, or inside a cluster that has formed a horseshoe.
 */
bool FPGATrackSimCLUSTERING::updatePixelCluster(FPGATrackSimCluster &currentCluster, FPGATrackSimHit &incomingHit, bool newCluster, bool digitalClustering){
 
    if(newCluster){
        FPGATrackSimHit newHit = incomingHit;
        newHit.setEtaIndex(incomingHit.getEtaIndex());
        newHit.setPhiIndex(incomingHit.getPhiIndex());
        newHit.setEtaCoord(incomingHit.getEtaCoord());
        newHit.setPhiCoord(incomingHit.getPhiCoord());
        newHit.setCentroidPhiIndex(incomingHit.getPhiIndex());
        newHit.setCentroidEtaIndex(incomingHit.getEtaIndex());
        newHit.setEtaWidth(1);
        newHit.setPhiWidth(1);
        //Set the initial clusterEquiv to be the incoming hit with double precision
        currentCluster.setClusterEquiv(newHit);
        //Add the current hit to the list of hits
        currentCluster.push_backHitList(incomingHit);
        //It doesn't really matter, as we will be at the end of the hit loop, but we did technically "cluster" this hit
        return true;
    } else {
        int hitCol = incomingHit.getEtaIndex();
        int hitRow = incomingHit.getPhiIndex();
 
        FPGATrackSimHit clusterEquiv = currentCluster.getClusterEquiv();
        int clusterCol = clusterEquiv.getEtaIndex();
        int clusterColWidth = clusterEquiv.getEtaWidth();
        int clusterRow = clusterEquiv.getPhiIndex();
        int clusterRowWidth = clusterEquiv.getPhiWidth();
 
        if ((hitCol == clusterCol + clusterColWidth) && (hitRow == clusterRow + clusterRowWidth)) {
            clusterColWidth++;
            clusterRowWidth++;
 
            return FPGATrackSimCLUSTERING::updateClusterContents(currentCluster, clusterRow, clusterRowWidth, clusterCol, clusterColWidth, incomingHit, digitalClustering);
        } else if ((hitCol == clusterCol + clusterColWidth) && (hitRow == clusterRow - 1)) {
            clusterColWidth++;
            clusterRow--;
            clusterRowWidth++;
 
            return FPGATrackSimCLUSTERING::updateClusterContents(currentCluster, clusterRow, clusterRowWidth, clusterCol, clusterColWidth, incomingHit, digitalClustering);
        } else if ((hitCol >= clusterCol) && (hitCol < clusterCol + clusterColWidth) && (hitRow == clusterRow + clusterRowWidth)) {
            clusterRowWidth++;
 
            return FPGATrackSimCLUSTERING::updateClusterContents(currentCluster, clusterRow, clusterRowWidth, clusterCol, clusterColWidth, incomingHit, digitalClustering);
        } else if ((hitCol == clusterCol + clusterColWidth) && (hitRow >= clusterRow) && (hitRow < clusterRow + clusterRowWidth)) {
            clusterColWidth++;
 
            return FPGATrackSimCLUSTERING::updateClusterContents(currentCluster, clusterRow, clusterRowWidth, clusterCol, clusterColWidth, incomingHit, digitalClustering);
        } else if ((hitCol >= clusterCol) && (hitCol < clusterCol + clusterColWidth) && (hitRow == clusterRow - 1)) {
            clusterRow--;
            clusterRowWidth++;
 
            return FPGATrackSimCLUSTERING::updateClusterContents(currentCluster, clusterRow, clusterRowWidth, clusterCol, clusterColWidth, incomingHit, digitalClustering);
        } else if ((hitCol == clusterCol - 1) && (hitRow == clusterRow - 1)) {
            clusterCol--;
            clusterColWidth++;
            clusterRow--;
            clusterRowWidth++;
 
            return FPGATrackSimCLUSTERING::updateClusterContents(currentCluster, clusterRow, clusterRowWidth, clusterCol, clusterColWidth, incomingHit, digitalClustering);
        } else if ((hitCol == clusterCol - 1) && (hitRow >= clusterRow) && (hitRow < clusterRow + clusterRowWidth)) {
            clusterCol--;
            clusterColWidth++;
 
            return FPGATrackSimCLUSTERING::updateClusterContents(currentCluster, clusterRow, clusterRowWidth, clusterCol, clusterColWidth, incomingHit, digitalClustering);
        } else if ((hitCol == clusterCol - 1) && (hitRow == clusterRow + clusterRowWidth)) {
            clusterCol--;
            clusterColWidth++;
            clusterRowWidth++;
 
            return FPGATrackSimCLUSTERING::updateClusterContents(currentCluster, clusterRow, clusterRowWidth, clusterCol, clusterColWidth, incomingHit, digitalClustering);
        } else if ((hitCol >= clusterCol) && (hitCol < clusterCol + clusterColWidth) && (hitRow >= clusterRow) && (hitRow < clusterRow + clusterRowWidth)) {
            return FPGATrackSimCLUSTERING::updateClusterContents(currentCluster, clusterRow, clusterRowWidth, clusterCol, clusterColWidth, incomingHit, digitalClustering);
        } else {
            return false;
        }
    }
}
 
/*
 * This function is used in the FPGATrackSimClusteringTools to see if a new hit should be added to the current cluster under construction. It assumes double precision hits.
 * It checks if the hit is in a number of positions w.r.t. the cluster being formed: up/right, down/right, above, right, or inside a cluster that has formed a horseshoe.
 */
bool FPGATrackSimCLUSTERING::updateStripCluster(FPGATrackSimCluster &currentCluster, FPGATrackSimHit &incomingHit, bool newCluster, bool digitalClustering){
 
    CXXUTILS_TRAPPING_FP;
 
    // Shift initial widths 1->0, 2->2, 3->4, 4->6 etc...
    //The groupSize is stored in the EtaWidth
    int tempWidth = (incomingHit.getPhiWidth()*fpgatracksim::scaleHitFactor)-fpgatracksim::scaleHitFactor;
    // Now shift to pixel width equivalents, 0->0, 2->1, 4->2, 6->3 etc...
    if(tempWidth > 0) tempWidth = tempWidth/fpgatracksim::scaleHitFactor;
    if(newCluster){
        FPGATrackSimHit newHit = incomingHit;
        //Double the precision of the strip positions.
        int tempCentroid = incomingHit.getPhiIndex()*fpgatracksim::scaleHitFactor;
        // Now shift the centroid phi+phiWidth, and store the width (put it back in the PhiWidth)
        newHit.setPhiIndex(incomingHit.getPhiIndex());
        newHit.setCentroidPhiIndex(tempCentroid+tempWidth);
        newHit.setPhiWidth(tempWidth);
        //Set the initial clusterEquiv to be the incoming hit with double precision
        currentCluster.setClusterEquiv(newHit);
        //Add the current hit to the list of hits
        currentCluster.push_backHitList(incomingHit);
        //It doesn't really matter, as we will be at the end of the hit loop, but we did technically "cluster" this hit
        return true;
    } else {
        //Now get the --START-- of the new strip cluster
        int hitRow = incomingHit.getEtaIndex();
        int hitCol = incomingHit.getPhiIndex()*fpgatracksim::scaleHitFactor;
 
        FPGATrackSimHit clusterEquiv = currentCluster.getClusterEquiv();
        int clusterRow = clusterEquiv.getEtaIndex();
        int clusterRowWidth = clusterEquiv.getEtaWidth();
        int clusterCol = clusterEquiv.getCentroidPhiIndex();
        int clusterColWidth = clusterEquiv.getPhiWidth();
 
        //Looking for a neighbour to the right. i.e. find the end of the current cluster (Col+width) and look in the next cell (+2). Compare this to the start of the new cluster. This is unlikely/impossible(?) to happen due to preclustering.
        if(hitCol == clusterCol+clusterColWidth+fpgatracksim::scaleHitFactor && hitRow == clusterRow) {
            //The new centroid will be the original column position, minus its width, plus the new width
            //So subtract the original width...
            clusterCol = clusterCol - clusterColWidth;
            //The new width will be the combination of the current widths, ++
            clusterColWidth = clusterColWidth+tempWidth+1;
            //And add on the new width
            clusterCol = clusterCol + clusterColWidth;
            FPGATrackSimCLUSTERING::updateClusterContents(currentCluster, clusterRow, clusterRowWidth, clusterCol, clusterColWidth, incomingHit, digitalClustering);
            return true;
        } else return false;
    }
}
 
 
bool FPGATrackSimCLUSTERING::updateClusterContents(FPGATrackSimCluster &currentCluster, int &clusterRow, int &clusterRowWidth, int &clusterCol, int &clusterColWidth, FPGATrackSimHit &incomingHit, bool digitalClustering) {
    //Grab the cluster equiv
    FPGATrackSimHit clusterEquiv = currentCluster.getClusterEquiv();
    bool isConnected = false;
 
    //Check if connected to another hit in the cluster
    if(incomingHit.isPixel()){
        for (auto & hit : currentCluster.getHitList()) {
            auto hitEta = hit.getEtaIndex();
            auto hitPhi = hit.getPhiIndex();
            auto inHitEta = incomingHit.getEtaIndex();
            auto inHitPhi = incomingHit.getPhiIndex();
 
            if (((inHitEta == hitEta - 1) && (inHitPhi == hitPhi - 1)) ||
                    ((inHitEta == hitEta + 1) && (inHitPhi == hitPhi - 1)) ||
                    ((inHitEta == hitEta - 1) && (inHitPhi == hitPhi + 1)) ||
                    ((inHitEta == hitEta + 1) && (inHitPhi == hitPhi + 1)) ||
                    ((inHitEta == hitEta) && (inHitPhi == hitPhi - 1)) ||
                    ((inHitEta == hitEta) && (inHitPhi == hitPhi + 1)) ||
                    ((inHitEta == hitEta - 1) && (inHitPhi == hitPhi)) ||
                    ((inHitEta == hitEta + 1) && (inHitPhi == hitPhi))) {
                isConnected = true;
                break;
            }
        }
        if (!isConnected)
            return false;
    }
 
    //Update the clusterEquiv's position and width
    if (incomingHit.isPixel()) {
        clusterEquiv.setEtaIndex(clusterCol);
        clusterEquiv.setEtaWidth(clusterColWidth);
        clusterEquiv.setPhiIndex(clusterRow);
        clusterEquiv.setPhiWidth(clusterRowWidth);
    } else {
        clusterEquiv.setEtaIndex(clusterRow);
        clusterEquiv.setEtaWidth(clusterRowWidth);
        clusterEquiv.setCentroidPhiIndex(clusterCol);
        clusterEquiv.setPhiWidth(clusterColWidth);
    }
 
 
    float xOld = clusterEquiv.getX();
    float yOld = clusterEquiv.getY();
    float zOld = clusterEquiv.getZ();
    float xPhiOld = clusterEquiv.getPhiCoord();
    float xEtaOld = clusterEquiv.getEtaCoord();
    float cPhiOld = clusterEquiv.getCentroidPhiIndex();
    float cEtaOld = clusterEquiv.getCentroidEtaIndex();
    float xNew = incomingHit.getX();
    float yNew = incomingHit.getY();
    float zNew = incomingHit.getZ();
    float xPhiNew = incomingHit.getPhiCoord();
    float xEtaNew = incomingHit.getEtaCoord();
    float cPhiNew = incomingHit.getPhiIndex();
    float cEtaNew = incomingHit.getEtaIndex();
    int tot = clusterEquiv.getToT();
    int totNew = incomingHit.getToT();
    //As strips arrive pre-clustered, this is different for pixels/strips
    if(incomingHit.isPixel()){
        CXXUTILS_TRAPPING_FP;
        if (digitalClustering) {
            // n+1 because that is old + new now
            int n = currentCluster.getHitList().size();
            clusterEquiv.setX((xOld*n + xNew) / (n+1));
            clusterEquiv.setY((yOld*n + yNew) / (n+1));
            clusterEquiv.setZ((zOld*n + zNew) / (n+1));
            clusterEquiv.setPhiCoord((xPhiOld*n + xPhiNew) / (n+1));
            clusterEquiv.setEtaCoord((xEtaOld*n + xEtaNew) / (n+1));
            clusterEquiv.setCentroidPhiIndex((cPhiOld*n + cPhiNew) / (n+1));
            clusterEquiv.setCentroidEtaIndex((cEtaOld*n + cEtaNew) / (n+1));
        } else {
            clusterEquiv.setX((xOld*tot + xNew*totNew) / (tot+totNew));
            clusterEquiv.setY((yOld*tot + yNew*totNew) / (tot+totNew));
            clusterEquiv.setZ((zOld*tot + zNew*totNew) / (tot+totNew));
            clusterEquiv.setPhiCoord((xPhiOld*tot + xPhiNew*totNew) / (tot+totNew));
            clusterEquiv.setEtaCoord((xEtaOld*tot + xEtaNew*totNew) / (tot+totNew));
            clusterEquiv.setCentroidPhiIndex((cPhiOld*tot + cPhiNew*totNew) / (tot+totNew));
            clusterEquiv.setCentroidEtaIndex((cEtaOld*tot + cEtaNew*totNew) / (tot+totNew));
        }
    } else {
        //Phi width + 1 for the seed is the width of the current cluster
        int N = currentCluster.getClusterEquiv().getPhiWidth()+1;
        //Phi width of an incoming strip is the width of the cluster
        int newN = incomingHit.getPhiWidth();
        //Now as above, N+newN
        clusterEquiv.setPhiCoord((xPhiOld*N + xPhiNew*newN) / (N+newN));
        clusterEquiv.setX((xOld*N + xNew*newN) / (N+newN));
        clusterEquiv.setY((yOld*N + yNew*newN) / (N+newN));
        clusterEquiv.setZ((zOld*N + zNew*newN) / (N+newN));
    }
    clusterEquiv.setToT(tot + totNew);
 
    //Put it back
    currentCluster.setClusterEquiv(clusterEquiv);
 
    //Pushback the hit into the hitlist
    currentCluster.push_backHitList(incomingHit);
 
    return true;
}
 
/* Sort for the ordering of ITk modules: Sort by ETA.
*/
bool FPGATrackSimCLUSTERING::sortITkInputEta(const std::unique_ptr<FPGATrackSimHit>& hitA, const std::unique_ptr<FPGATrackSimHit>& hitB)
{
    if (hitA->getIdentifierHash() != hitB->getIdentifierHash())
        return hitA->getIdentifierHash() < hitB->getIdentifierHash();
    return hitA->getEtaIndex() < hitB->getEtaIndex();
}
 
/* Sort for the ordering of ITk modules: Sort by PHI.
*/
bool FPGATrackSimCLUSTERING::sortITkInputPhi(const std::unique_ptr<FPGATrackSimHit>& hitA, const std::unique_ptr<FPGATrackSimHit>& hitB)
{
    if (hitA->getIdentifierHash() != hitB->getIdentifierHash())
        return hitA->getIdentifierHash() < hitB->getIdentifierHash();
    return hitA->getPhiIndex() < hitB->getPhiIndex();
}
 
/* Cap precision of the global coordinates in r, phi, z */
void FPGATrackSimClusteringTool::reduceGlobalCoordPrecision(FPGATrackSimCluster &cluster) const {
    FPGATrackSimHit clusterEquiv = cluster.getClusterEquiv();
    float pos[3] = { clusterEquiv.getR(), clusterEquiv.getGPhi(), clusterEquiv.getZ() };
 
    pos[0] = std::trunc(pos[0] / m_coordRPrecision) * m_coordRPrecision;
    pos[1] = std::trunc(pos[1] / m_coordPhiPrecision) * m_coordPhiPrecision;
    pos[2] = std::trunc(pos[2] / m_coordZPrecision) * m_coordZPrecision;
 
    clusterEquiv.setX(pos[0] * std::cos(pos[1]));
    clusterEquiv.setY(pos[0] * std::sin(pos[1]));
    clusterEquiv.setZ(pos[2]);
 
    cluster.setClusterEquiv(clusterEquiv);
}
 
void FPGATrackSimClusteringTool::reduceGlobalCoordPrecision(FPGATrackSimHit &hit) const {
    float pos[3] = { hit.getR(), hit.getGPhi(), hit.getZ() };
 
    pos[0] = std::trunc(pos[0] / m_coordRPrecision) * m_coordRPrecision;
    pos[1] = std::trunc(pos[1] / m_coordPhiPrecision) * m_coordPhiPrecision;
    pos[2] = std::trunc(pos[2] / m_coordZPrecision) * m_coordZPrecision;
 
    hit.setX(pos[0] * std::cos(pos[1]));
    hit.setY(pos[0] * std::sin(pos[1]));
    hit.setZ(pos[2]);
}