TrackFitter.cxx

Go to the documentation of this file.

 
// Copyright (C) 2002-2025 CERN for the benefit of the ATLAS collaboration
 
#include <algorithm>
#include <iostream>
#include <fstream>
#include <cmath>
 
#include "FPGATrackSimObjects/FPGATrackSimMultiTruth.h"
#include "GaudiKernel/MsgStream.h"
#include "AthenaKernel/getMessageSvc.h"
 
 
#include "FPGATrackSimAlgorithms/TrackFitter.h"
 
 
std::vector<FPGATrackSimTrack>::const_iterator getBestChi2(std::vector<FPGATrackSimTrack> const & tracks);
bool hasGoodFit(std::vector<FPGATrackSimTrack> const & track_cands, float minchi2);
 

// Constructor
 
TrackFitter::TrackFitter(const FPGATrackSimFitConstantBank *nominalbank,
             const std::vector<const FPGATrackSimFitConstantBank*>& droppedLayerbanks, bool guessingHits) :
  AthMessaging (Athena::getMessageSvc(), "FPGATrackSimTrackFitter") 
{
  m_guessinghits = guessingHits;
  m_nominalBank = nominalbank;
 
  if (!m_nominalBank) ATH_MSG_FATAL("Constant bank not set");
  if (!m_guessinghits)
    for (auto bank : droppedLayerbanks) 
      m_droppedLayerBanks.push_back(bank);
}
 
 
void TrackFitter::resetCounters()
{
    m_idbase = 0;
    m_nfits = 0;
    m_nfits_maj = 0;
    m_nfits_maj_pix = 0;
    m_nfits_maj_SCT = 0;
    m_nfits_rec = 0;
    m_nfits_addrec = 0;
}
 

// fitTracks
 
 
int TrackFitter::fitTracks(const std::vector<std::shared_ptr<const FPGATrackSimRoad>>& roads, std::vector<FPGATrackSimTrack>& tracks) {
    resetCounters();
    for (const std::shared_ptr<const FPGATrackSimRoad>& cur_road : roads) {
      std::vector<FPGATrackSimTrack> t;
      int isOK = fitTracks(cur_road, t);
      if (isOK != FITTRACKS_OK) return isOK;
      else tracks.insert(tracks.end(), t.begin(), t.end());
    }
    return FITTRACKS_OK;
}
 

 int TrackFitter::fitTracks(const std::shared_ptr<const FPGATrackSimRoad> &road, std::vector<FPGATrackSimTrack>& tracks)
{
    if (not road){
      ATH_MSG_WARNING("road pointer is null in TrackFitter::fitTracks");
      return FITTRACKS_BAD;
    }
    
    m_tracks_missinghits_track.clear();
 
    double y = 0.0;
    double x = 0.0;
    if (m_IdealCoordFitType != TrackCorrType::None ) {
      y = road->getY();
      x = road->getX();
      ATH_MSG_DEBUG("Attempting to fit Hough road with y = " << y << ", x = " << x << ", sector = " << road->getSector() << "and nhits = " << road->getNHits());
    }
 
    int sector = 0;
    if (!m_fitFromRoad) {
      // Error checking
      sector = road->getSector();
      if (sector < 0) {
        ATH_MSG_DEBUG("Bad sector " << sector);
        return FITTRACKS_OK;
      }
      else if (sector >= m_nominalBank->getNSectors()) {
        ATH_MSG_WARNING("Constants for sector " << sector << " don't exist");
        return FITTRACKS_BAD;
      }
      else if (!m_nominalBank->getIsGood(sector)) {
    ATH_MSG_WARNING("Constants for sector " << sector << " are not valid");
    return FITTRACKS_BAD;
      }
    }
 
    // Get info on layers with missing hits
    int nMissing;
    bool missPixel;
    bool missStrip;
    layer_bitmask_t missing_mask;
    layer_bitmask_t norecovery_mask; // mask to prevent majority in planes with multiple hits
    getMissingInfo(*road, nMissing, missPixel, missStrip, missing_mask, norecovery_mask);
    // Create a template track with common parameters filled already for initializing below
    FPGATrackSimTrack temp;
    if(!m_do2ndStage){
      temp.setTrackStage(TrackStage::FIRST);
      if (!m_fitFromRoad) temp.setFirstSectorID(road->getSector());
    }
    else{
      temp.setTrackStage(TrackStage::SECOND);
      if (!m_fitFromRoad) temp.setSecondSectorID(road->getSector());
    }
    temp.setNLayers(m_pmap->getNLogiLayers());
    temp.setBankID(-1); // TODO
    temp.setPatternID(road->getPID());
    temp.setHitMap(missing_mask);
    temp.setNMissing(nMissing);
    temp.setHoughX(x);
    temp.setHoughY(y);
    temp.setQOverPt(y);
    temp.setTrackCorrType(m_IdealCoordFitType);
    temp.setDoDeltaGPhis(m_doDeltaGPhis);
 
    temp.setSubRegion(road->getSubRegion());
    temp.setHoughXBin(road->getXBin());
    temp.setHoughYBin(road->getYBin());
 
    temp.setBinIdx(road->getBinIdx());
 
    // Create a list of track candidates by taking all possible combinations of hits in road.
    std::vector<FPGATrackSimTrack> track_cands;
 
    // This may not even need to be a class member. Create the vector of possible combinations.
    std::vector<std::vector<int>> comboIndices = getComboIndices(road->getNHits_layer());
    ATH_MSG_DEBUG("There are theoretically " << comboIndices.size() << " combinations to process for this road");
    size_t nFits = 0;
    for (size_t icomb = 0; icomb < comboIndices.size(); icomb++) {
        FPGATrackSimTrack track_cand = makeTrackCandidate(*road, temp, comboIndices[icomb]);
 
        // Before we start, make sure this track candidate has not been marked as invalid
        // due to combinatorics issues with spacepoints.
        bool ok = track_cand.isValidCand();
        if (!ok) {
            continue;
        }
 
        // If the "fit from road" flag is set, then just assign track chi2 and parameters from that.
        if (!m_do2ndStage && m_fitFromRoad) {
            // Then actually do it using the road values.
            track_cand.setChi2(road->getFitChi2());
            track_cand.setChi2Phi(road->getFitChi2Phi());
            track_cand.setChi2Eta(road->getFitChi2Eta());
            track_cand.setPars(road->getFitParams());
            ATH_MSG_DEBUG("Assigned chi2 = " << track_cand.getChi2() << " and parameters from genscan tool");
            ATH_MSG_DEBUG("Set q/pt = " << track_cand.getQOverPt());
            ATH_MSG_DEBUG("Set d0 = " << track_cand.getD0());
            ATH_MSG_DEBUG("Set z0 = " << track_cand.getZ0());
            ATH_MSG_DEBUG("Set eta = " << track_cand.getEta());
            ATH_MSG_DEBUG("Set phi = " << track_cand.getPhi());
            tracks.push_back(track_cand);
            continue;
        } else {
 
        if (nMissing == 0 || m_guessinghits)
        {
            ok = m_nominalBank->linfit(sector, track_cand, m_do2ndStage);
        }
        else
        {
            unsigned int ic = 0; // this will be the leading coordinate with a missing hit
            for (unsigned int ic = 0; ic < m_pmap->getNCoords(); ic++) {
                if (!( (missing_mask >> ic) & 0x1)) break; // if we found one, stop
            }
            // banks are indexed by plane, not by coordinate
            ok = m_droppedLayerBanks[m_pmap->getCoordLayer(ic)]->linfit(sector, track_cand, m_do2ndStage);
        }
 
        if (!ok) {
            continue;
        }
        track_cand.setOrigChi2(track_cand.getChi2());
        if (getDoMissingHitsCheck()) {
            if (track_cand.getNMissing() == 0) { // missing 0 hits, drop hits one at a time and refit!
                for (unsigned icoord = 0; icoord < m_pmap->getNCoords(); icoord++) { // 9 coordinates to refit, nominally
 
                    FPGATrackSimTrack newtrack = track_cand;
                    newtrack.setNMissing(1);
                    unsigned int bitmask = ( 1 << m_pmap->getNCoords() ) - 1; // all bits set to 1
                    bitmask &= ~(1 << icoord); // remove this coordinate
                    if (m_pmap->getDim(icoord) == 2) { // if we need to do two coordinates at once
                        bitmask &= ~(1 << (icoord+1)); // remove the next coordinate, too
                        newtrack.setNMissing(2);
                    }
                    newtrack.setHitMap(bitmask); // update bitmask
                    m_nominalBank->linfit(sector, newtrack, m_do2ndStage);
                    m_tracks_missinghits_track.push_back(newtrack);
                    if (m_pmap->getDim(icoord) == 2) {
                        icoord++; // skip 2nd of pixel coordinates so we don't do them twice
                    }
                }
            }
        }
        }
        tracks.push_back(track_cand);
 
        // Enforce m_max_ncomb here, but only for the tracks that we actually fit.
        // Let's see how often this warning fires.
        nFits += 1;
        if (nFits > (size_t)m_max_ncomb) {
            ATH_MSG_WARNING("Reached " << m_max_ncomb << " fits when processing track combinations for single road, breaking");
        }
    }
    ATH_MSG_DEBUG("Processed " << tracks.size() << " actual/valid tracks for this road");
 
    // Increment counters
    m_nfits += nFits;
    if (nMissing > 0)
    {
        m_nfits_maj += nFits;
        if (missPixel) m_nfits_maj_pix += nFits;
        if (missStrip) m_nfits_maj_SCT += nFits;
    }
 
    // Do recovery fits
    if ((nMissing == 0)&&(!m_fitFromRoad)) {
        // In the case of m_do_majority > 1, we only do majority fits if ALL full fits fail the chi2 cut
        if (m_do_majority == 1 || (m_do_majority > 1 && !hasGoodFit(tracks, m_Chi2Dof_recovery_min)))
        {
            for (FPGATrackSimTrack & t : tracks)
          if (t.getChi2ndof() > m_Chi2Dof_recovery_min && t.getChi2ndof() < m_Chi2Dof_recovery_max){
        double y(0);
        if (road != nullptr && m_IdealCoordFitType != TrackCorrType::None)
          y = road->getY();
        t = recoverTrack(t, norecovery_mask, sector, y);
 
          }
        }
    }
 
    // Add truth info
    for (FPGATrackSimTrack & t : tracks) {
        compute_truth(t); // match the track to a geant particle using the channel-level geant info in the hit data.
    }
 
    return FITTRACKS_OK;
}
 

// fitTracks core helper functions
 
 
// Given road, populates the supplied variables with info on which layers missed hits
void TrackFitter::getMissingInfo(const FPGATrackSimRoad & road, int & nMissing, bool & missPixel, bool & missStrip,
        layer_bitmask_t & missing_mask, layer_bitmask_t & norecovery_mask)
{
    nMissing = m_pmap->getNCoords(); // init with nCoords and decrement as we find misses
    missPixel = false;
    missStrip = false;
    missing_mask = 0;
    norecovery_mask = 0;
    unsigned int wclayers = road.getWCLayers();
    for (unsigned layer = 0; layer < m_pmap->getNLogiLayers(); layer++)
    {
        int nHits = road.getHits(layer).size();
        if (nHits==0)
        {
            if (m_IdealCoordFitType == TrackCorrType::None && ((wclayers >> layer) & 1))
            {
                int ix = m_pmap->getCoordOffset(layer);
                int iy = ix + 1;
                if (m_pmap->isSCT(layer))
                {
                    missing_mask |= 1 << ix;
                    nMissing -= 1;
                }
                else
                {
                    missing_mask |= (1<<ix) | (1<<iy);
                    nMissing -= 2;
                }
            }
            else
            {
                if (m_pmap->isSCT(layer)) missStrip = true;
                else missPixel = true;
            }
        }
        else if (!((wclayers >> layer) & 1)) { // we have a hit
            int ix = m_pmap->getCoordOffset(layer);
            int iy = ix + 1;
            if (m_pmap->isSCT(layer))
          {
                missing_mask |= 1 << ix;
                nMissing -= 1;
          }
            else
          {
                missing_mask |= (1<<ix) | (1<<iy);
                nMissing -= 2;
          }
 
            // set the mask to avoid recovery in crowded planes
            if (m_norecovery_nhits != -1 && nHits > m_norecovery_nhits)
                norecovery_mask |= (1<<layer);
        }
    }
}
 
 

void TrackFitter::makeTrackCandidates(const FPGATrackSimRoad & road, const FPGATrackSimTrack & temp, std::vector<FPGATrackSimTrack>& track_cands)
{
    std::vector<std::vector<int>> combs = getComboIndices(road.getNHits_layer());
    track_cands.resize(combs.size(), temp);
 
    //get the WC hits:
    layer_bitmask_t wcbits= road.getWCLayers();
    for (size_t icomb = 0; icomb < combs.size(); icomb++)
    {
      //Need to set the ID and the hits size of this track
      track_cands[icomb].setTrackID(m_idbase + icomb);
      track_cands[icomb].setNLayers(m_pmap->getNLogiLayers());
 
      // If this is an idealized coordinate fit; keep references to the idealized radii.
      track_cands[icomb].setIdealRadii(m_rmap->getAvgRadii(0));
 
      std::vector<int> const & hit_indices = m_comboIndices[icomb]; // size nLayers
        for (unsigned layer = 0; layer < m_pmap->getNLogiLayers(); layer++)
        {
            if (hit_indices[layer] < 0) // Set a dummy hit if road has no hits in this layer
            {
                FPGATrackSimHit newhit=FPGATrackSimHit();
                newhit.setLayer(layer);
                newhit.setSection(0);
                if (m_pmap->getDim(layer) == 2) newhit.setDetType(SiliconTech::pixel);
                else newhit.setDetType(SiliconTech::strip);
 
                if (wcbits & (1 << layer ) ) {
                    newhit.setHitType(HitType::wildcard);
            newhit.setLayer(layer);
        }
                
                track_cands[icomb].setFPGATrackSimHit(layer, newhit);
            }
            else
            {
                const std::shared_ptr<const FPGATrackSimHit> hit = road.getHits(layer)[hit_indices[layer]];
                // If this is an outer spacepoint, and it is not the same as the inner spacepoint, reject it.
                // Here we "reject" it by marking the candidate as "invalid", to be rejected later.
                // That require another field on the track object, but it avoids having to change the sizes
                // of arrays computed above.
                if (hit->getHitType() == HitType::spacepoint && hit->getSide() == 1 && layer > 0) {
                    const FPGATrackSimHit inner_hit = track_cands[icomb].getFPGATrackSimHits().at(layer - 1);//avoid negative index
                    if ((hit->getX() != inner_hit.getX()) || (hit->getY() != inner_hit.getY()) || (hit->getZ() != inner_hit.getZ())) {
                        track_cands[icomb].setValidCand(false);
                    }
                }
                track_cands[icomb].setFPGATrackSimHit(layer, *hit);
            }
        }
    }
 
    m_idbase += combs.size();
}
 
// This is a version of the above that produces a single track candidate on demand.
FPGATrackSimTrack TrackFitter::makeTrackCandidate(const FPGATrackSimRoad & road, const FPGATrackSimTrack & temp, const std::vector<int>& hit_indices)
{
    m_idbase += 1;
 
    FPGATrackSimTrack track_cand(temp);
 
    //get the WC hits:
    layer_bitmask_t wcbits= road.getWCLayers();
 
    //Need to set the ID and the hits size of this track
    track_cand.setTrackID(m_idbase);
    track_cand.setNLayers(m_pmap->getNLogiLayers());
 
    // If this is an idealized coordinate fit; keep references to the idealized radii.
    track_cand.setIdealRadii(m_rmap->getAvgRadii(0));
 
    for (unsigned layer = 0; layer < m_pmap->getNLogiLayers(); layer++)
    {
        if (hit_indices[layer] < 0) // Set a dummy hit if road has no hits in this layer
        {
            FPGATrackSimHit newhit=FPGATrackSimHit();
            newhit.setLayer(layer);
            newhit.setSection(0);
            if (m_pmap->getDim(layer) == 2) newhit.setDetType(SiliconTech::pixel);
            else newhit.setDetType(SiliconTech::strip);
 
            if (wcbits & (1 << layer ) ) {
                newhit.setHitType(HitType::wildcard);
                newhit.setLayer(layer);
            }
 
            track_cand.setFPGATrackSimHit(layer, newhit);
        }
        else
        {
            const std::shared_ptr<const FPGATrackSimHit> hit = road.getHits(layer)[hit_indices[layer]];
            // If this is an outer spacepoint, and it is not the same as the inner spacepoint, reject it.
            // Here we "reject" it by marking the candidate as "invalid", to be rejected later.
            // That require another field on the track object, but it avoids having to change the sizes
            // of arrays computed above.
            if (hit->getHitType() == HitType::spacepoint && hit->getSide() == 1 && layer > 0) {
                const FPGATrackSimHit inner_hit = track_cand.getFPGATrackSimHits().at(layer - 1);//avoid negative index
                if ((hit->getX() != inner_hit.getX()) || (hit->getY() != inner_hit.getY()) || (hit->getZ() != inner_hit.getZ())) {
                    track_cand.setValidCand(false);
                    break;
                }
            }
            track_cand.setFPGATrackSimHit(layer, *hit);
        }
    }
 
    return track_cand;
}

FPGATrackSimTrack TrackFitter::recoverTrack(FPGATrackSimTrack const & t, sector_t sector, layer_bitmask_t norecovery_mask, double qoverpt)
{
    m_nfits_rec += 1;
 
    std::vector<FPGATrackSimTrack> recovered_tracks(m_pmap->getNLogiLayers(),t); // start with the complete track in all layers
 
    float best_chi2ndof = t.getChi2ndof();
    int best_i = -1;
 
    // Remove a hit from each layer and do a fit
    for (unsigned layer = (m_require_first ? 1 : 0); layer < recovered_tracks.size(); ++layer)
    {
        // Skip planes with multiple hits
        if (norecovery_mask & (1<<layer)) continue;
        m_nfits_addrec += 1;
        // Zero the cluster for the missing layer, make sure to set the number of dimensions for it
        
        FPGATrackSimHit newhit=FPGATrackSimHit();
        newhit.setLayer(layer);
        newhit.setSection(0);
        newhit.setHitType(HitType::guessed);
        recovered_tracks[layer].setFPGATrackSimHit(layer,newhit);
        // Set the number of missing points and the related bitmask
        int ix = m_pmap->getCoordOffset(layer);
 
        unsigned int missing_mask = t.getHitMap();
        missing_mask &= ~(1 << ix);
 
        if (m_pmap->isPixel(layer)) missing_mask &= ~(1 << (ix + 1));
 
        recovered_tracks[layer].setHitMap(missing_mask);
        recovered_tracks[layer].setNMissing(t.getNMissing() + m_pmap->getDim(layer));
 
    // reset the pt
    recovered_tracks[layer].setQOverPt(qoverpt);
    recovered_tracks[layer].setTrackCorrType(m_IdealCoordFitType);
 
 
        // If not guessing do the fit using layer+1 since 1st set of constants are the 8/8
    if (m_guessinghits) m_nominalBank->linfit(sector, recovered_tracks[layer], m_do2ndStage);
    else m_droppedLayerBanks[layer]->linfit(sector, recovered_tracks[layer], m_do2ndStage);
 
 
        // Check if the chi2ndof is better
        if (recovered_tracks[layer].getChi2ndof() < best_chi2ndof)
        {
            best_chi2ndof = recovered_tracks[layer].getChi2ndof();
            best_i = layer;
        }
    }
 
    if (best_i >= 0)  {// a better track was found
        return recovered_tracks[best_i];
    } else
        return t;
}
 
 
 
 
// compute the best geant match, the barcode with the largest number of hits contributing to the track,
// and store in the track.
void TrackFitter::compute_truth(FPGATrackSimTrack & t) const
{
    std::vector<FPGATrackSimMultiTruth> mtv;
 
    for (unsigned layer = 0; layer < m_pmap->getNLogiLayers(); layer++)
    {
        if (t.getHitMap() & (1 << m_pmap->getCoordOffset(layer))) continue; // no hit in this plane
    //Sanity check that we have enough hits.
    if (layer < t.getFPGATrackSimHits().size()) mtv.push_back(t.getFPGATrackSimHits().at(layer).getTruth());
 
        // adjust weight for hits without (and also with) a truth match, so that each is counted with the same weight.
        mtv.back().assign_equal_normalization();
    }
 
    FPGATrackSimMultiTruth mt( std::accumulate(mtv.begin(), mtv.end(), FPGATrackSimMultiTruth(), FPGATrackSimMultiTruth::AddAccumulator()) );
    // frac is then the fraction of the total number of hits on the track attributed to the barcode.
 
    FPGATrackSimMultiTruth::Barcode tbarcode;
    FPGATrackSimMultiTruth::Weight tfrac;
    const bool ok = mt.best(tbarcode, tfrac);
    if (ok)
    {
        t.setEventIndex(tbarcode.first);
        t.setBarcode(tbarcode.second);
        t.setBarcodeFrac(tfrac);
    }
}
 

// findTracks utility helper functions
 
// Returns true if there is a fit in track_cands with chi2ndof < minchi2ndof
bool hasGoodFit(std::vector<FPGATrackSimTrack> const & track_cands, float minchi2ndof)
{   
    for (FPGATrackSimTrack const & t: track_cands)
        {
        if (t.getChi2ndof() > 0 && t.getChi2ndof() < minchi2ndof)
            return true;
        }   
    return false;
}