TrackTimeDefAndQualityAlg.cxx

Go to the documentation of this file.

 
#include "TrackTimeDefAndQualityAlg.h"
 
#include "StoreGate/ReadDecorHandle.h"
#include "StoreGate/ReadHandle.h"
#include "StoreGate/WriteDecorHandle.h"
 
#include <optional>
 
#include "xAODInDetMeasurement/PixelCluster.h"
#include "xAODInDetMeasurement/StripCluster.h"
#include "Acts/Utilities/TrackHelpers.hpp"
 
#include "ActsEvent/TrackContainer.h"
#include "ActsCalibrators/xAODUncalibMeasCalibrator.h"
namespace HGTD {
 
TrackTimeDefAndQualityAlg::TrackTimeDefAndQualityAlg(const std::string& name,
                                                     ISvcLocator* pSvcLocator)
    : AthReentrantAlgorithm(name, pSvcLocator) {}
 
StatusCode TrackTimeDefAndQualityAlg::initialize() {
 
  ATH_CHECK(m_trackParticleContainerKey.initialize());
  ATH_CHECK(m_layerHasExtensionKey.initialize());
  ATH_CHECK(m_holesHGTDKey.initialize(!m_doActs)); //HGTD_holes not produced in ActsHGTDTrackExtensionAlg for now so make it optional
  ATH_CHECK(m_layerClusterTimeKey.initialize());
  ATH_CHECK(m_layerClusterTruthClassKey.initialize(m_doTruth)); //Only initialize if truth information is available
  ATH_CHECK(m_time_dec_key.initialize());
  ATH_CHECK(m_time_res_dec_key.initialize());
  ATH_CHECK(m_hasValidTime_dec_key.initialize());
  ATH_CHECK(m_summarypattern_dec_key.initialize());
 
  return StatusCode::SUCCESS;
}

 
StatusCode TrackTimeDefAndQualityAlg::execute(const EventContext& ctx) const {
 
  const xAOD::TrackParticleContainer* track_particles{nullptr};
  ATH_CHECK(SG::get(track_particles, m_trackParticleContainerKey, ctx));
 
  SG::WriteDecorHandle<xAOD::TrackParticleContainer, float> time_handle(
      m_time_dec_key, ctx);
  SG::WriteDecorHandle<xAOD::TrackParticleContainer, float> timeres_handle(
      m_time_res_dec_key, ctx);
  SG::WriteDecorHandle<xAOD::TrackParticleContainer, uint8_t> hasValidTime_handle(
      m_hasValidTime_dec_key, ctx);
  SG::WriteDecorHandle<xAOD::TrackParticleContainer, uint32_t> summary_handle(
      m_summarypattern_dec_key, ctx);
 
  SG::ReadDecorHandle<xAOD::TrackParticleContainer, std::vector<float>>
    layerClusterTimeHandle(m_layerClusterTimeKey, ctx);
  ATH_CHECK(layerClusterTimeHandle.isValid());
 
  SG::ReadDecorHandle<xAOD::TrackParticleContainer, std::vector<char>>
      layerHasExtensionHandle(m_layerHasExtensionKey, ctx);
  ATH_CHECK(layerHasExtensionHandle.isValid());
 
  static const std::vector<char> s_no_holes(s_hgtd_layers, false);
  std::optional<SG::ReadDecorHandle<xAOD::TrackParticleContainer, std::vector<char>>>
      holesHGTDHandle;
  if (!m_doActs) {
    holesHGTDHandle.emplace(m_holesHGTDKey, ctx);
    ATH_CHECK(holesHGTDHandle->isValid());
  }
 
  static const std::vector<int> s_no_truth(s_hgtd_layers, 0);
  std::optional<SG::ReadDecorHandle<xAOD::TrackParticleContainer, std::vector<int>>>
      layerClusterTruthClass;
  if (m_doTruth) {
    layerClusterTruthClass.emplace(m_layerClusterTruthClassKey, ctx);
    ATH_CHECK(layerClusterTruthClass->isValid());
  }
 
  for (const auto* track_ptkl : *track_particles) {
    // runs the time consistency checks
    // if no hits are found in HGTD, returns a default time
    const std::vector<float>& times = layerClusterTimeHandle(*track_ptkl);
    const std::vector<char>& has_clusters = layerHasExtensionHandle(*track_ptkl);
    const std::vector<int>& hit_classification = m_doTruth ? (*layerClusterTruthClass)(*track_ptkl) : s_no_truth;
    const std::vector<char>& holes_HGTD = m_doActs ? s_no_holes : (*holesHGTDHandle)(*track_ptkl);
 
    CleaningResult res = runTimeConsistencyCuts(times,
                        has_clusters,
                        hit_classification);
 
    // check if the last hit on track was within the predefined area
    if (lastHitIsOnLastSurface(*track_ptkl)) {
      res.m_field |= (0b0000 << m_holes_ptrn_sft);
    } else {
      res.m_field |= (0b0001 << m_holes_ptrn_sft);
    }
    
    // keep which of the hits associated in reco were primary hits (truth info!)
    short prime_pattern = 0x0;
    for (short i = 0; i < s_hgtd_layers; i++) {
      if (res.m_hits.at(i).m_isprime) {
        prime_pattern |= (1 << i);
      }
    }
    res.m_field |= (prime_pattern << m_primes_ptrn_sft);
 
    // expected pattern : 'on which HGTD layer a hit was expected?' which means extrapolation has
    // reached an active sensor, whether a matching cluster was found or not. 
    // So ‘expected = has_cluster OR HGTD_holes’. 
    short expected_pattern = 0x0;
    for (short i = 0; i < s_hgtd_layers; i++) {
      if (has_clusters.at(i) || holes_HGTD.at(i))  {
        expected_pattern |= (1 << i);
      }
    }
    res.m_field |= (expected_pattern << m_exp_ptrn_sft);
 
    // decorate the track again with this info
    time_handle(*track_ptkl) = res.m_time;
    timeres_handle(*track_ptkl) = res.m_resolution;
    hasValidTime_handle(*track_ptkl) = res.m_hasValidTime;
    summary_handle(*track_ptkl) = res.m_field;
  }
  return StatusCode::SUCCESS;
}

 
TrackTimeDefAndQualityAlg::CleaningResult
TrackTimeDefAndQualityAlg::runTimeConsistencyCuts(const std::vector<float>& times,
                          const std::vector<char>& has_clusters,
                          const std::vector<int>& hit_classification) const {
  // get all available hits (see the struct Hit) in a first step
  std::array<Hit, s_hgtd_layers> valid_hits = getValidHits(times,
                               has_clusters,
                               hit_classification);
 
  CleaningResult result;
  result.m_hits = valid_hits;
  result.m_field = 0x0;
  result.m_time = m_default_time;
  result.m_resolution = m_default_time_res;
  result.m_hasValidTime = 0;
 
  short recoed_pattern = getValidPattern(valid_hits);
  // stored the pattern of hits as retrieved from the iterative extension
  result.m_field |= (recoed_pattern << m_recoed_ptrn_sft);
 
  short nhits = std::count_if(valid_hits.begin(), valid_hits.end(),
                              [](const Hit& hit) { return hit.m_isvalid; });
  if (nhits < 2) {
    // fill the patern with the 1 hit (or none) and return
    result.m_field |= (recoed_pattern << m_comp_ptrn_sft);
    result.m_time = meanTime(valid_hits);
    result.m_resolution = trackTimeResolution(valid_hits);
    result.m_hasValidTime = recoed_pattern ? 1 : 0;
    return result;
  } else if (nhits == 2) {
    // if the deltaT cut is  passed, the pattern stays the same, otherwise set
    // to 0 as no hit passes
    // TODO: find better way to treat this!
    if (passesDeltaT(valid_hits)) {
      result.m_field |= (recoed_pattern << m_comp_ptrn_sft); // stays the same
      result.m_time = meanTime(valid_hits);
      result.m_resolution = trackTimeResolution(valid_hits);
      result.m_hasValidTime = 1;
      return result;
    } else {
      result.m_field |= (0b0000 << m_comp_ptrn_sft); // no hit passes
      result.m_time = m_default_time; // TODO should I just use the mean?
      result.m_resolution = m_default_time_res;
      return result;
    }
 
  } else {
    // for 3 or 4 hits, remove hit(s) with worst chi2 if needed
    float chi2 = calculateChi2(valid_hits);
    // if the chi2 is below the threshold, keep all hits
    bool searching = chi2 > m_chi2_threshold;
    while (searching) {
      short remove_layer = findLayerWithBadChi2(valid_hits);
      setLayerAsInvalid(valid_hits, remove_layer);
      float new_chi2 = calculateChi2(valid_hits);
      nhits = std::count_if(valid_hits.begin(), valid_hits.end(),
                            [](const Hit& hit) { return hit.m_isvalid; });
      if (new_chi2 <= m_chi2_threshold or nhits < 3) {
        searching = false;
      }
    } // while loop ended
 
    short chi2_rej_pattern = getValidPattern(valid_hits);
 
    if (nhits == 2) {
      if (passesDeltaT(valid_hits)) {
        result.m_field |= (chi2_rej_pattern << m_comp_ptrn_sft);
        result.m_time = meanTime(valid_hits);
        result.m_resolution = trackTimeResolution(valid_hits);
        result.m_hasValidTime = 1;
        return result;
      } else {
        result.m_field |= (0b0000 << m_comp_ptrn_sft); // no hit passes
        result.m_time = m_default_time; // TODO should I just use the mean?
        result.m_resolution = m_default_time_res;
        return result;
      }
    } else {
      // 3 or 4 hits, chi2 passed
      result.m_field |= (chi2_rej_pattern << m_comp_ptrn_sft);
      result.m_time = meanTime(valid_hits);
      result.m_resolution = trackTimeResolution(valid_hits);
      result.m_hasValidTime = 1;
      return result;
    }
  }
}
 
std::array<TrackTimeDefAndQualityAlg::Hit, s_hgtd_layers>
TrackTimeDefAndQualityAlg::getValidHits(const std::vector<float>& times,
                    const std::vector<char>& has_clusters,
                    const std::vector<int>& hit_classification) const {
  std::array<Hit, s_hgtd_layers> valid_hits {};
 
  for (size_t i = 0; i < s_hgtd_layers; i++) {
    Hit& newhit = valid_hits[i];
    if (has_clusters.at(i)) {
      newhit.m_time = times.at(i);
      newhit.m_isprime = hit_classification.at(i) == 1;
      newhit.m_isvalid = true;
    }
    newhit.m_layer = i;
  }
 
  return valid_hits;
}
 
short TrackTimeDefAndQualityAlg::getValidPattern(
    const std::array<TrackTimeDefAndQualityAlg::Hit, s_hgtd_layers>& hits)
    const {
  short pattern = 0x0;
  for (short i = 0; i < s_hgtd_layers; i++) {
    if (hits.at(i).m_isvalid) {
      pattern |= (1 << i);
    }
  }
  return pattern;
}
 
float TrackTimeDefAndQualityAlg::calculateChi2(
    const std::array<Hit, s_hgtd_layers>& hits) const {
 
  float mean = meanTime(hits);
 
  float chi2 = 0.;
  for (const auto& hit : hits) {
    if (hit.m_isvalid) {
      chi2 += (hit.m_time - mean) * (hit.m_time - mean) /
              (hit.m_resolution * hit.m_resolution);
    }
  }
  return chi2;
}
 
bool TrackTimeDefAndQualityAlg::passesDeltaT(
    const std::array<TrackTimeDefAndQualityAlg::Hit, s_hgtd_layers>& hits)
    const {
  // don't trust the user here.
  short n_valid = std::count_if(hits.begin(), hits.end(),
                                [](const Hit& hit) { return hit.m_isvalid; });
  if (n_valid != 2) {
    return false;
  }
  // FIXME this should be doable in a simpler manner...
  std::vector<float> times;
  std::vector<float> res;
  for (const auto& hit : hits) {
    if (hit.m_isvalid) {
      times.push_back(hit.m_time);
      res.push_back(hit.m_resolution);
    }
  }
  // pass if the distance in units of the resolution passes the cut
  return std::abs(times.at(0) - times.at(1)) <
         m_deltat_cut * hypot(res.at(0), res.at(1));
}
 
float TrackTimeDefAndQualityAlg::meanTime(
    const std::array<TrackTimeDefAndQualityAlg::Hit, s_hgtd_layers>& hits)
    const {
  float sum = 0.;
  short n = 0;
  for (const auto& hit : hits) {
    if (hit.m_isvalid) {
      sum += hit.m_time;
      n++;
    }
  }
  return n == 0 ? m_default_time.value() : sum / (float)n;
}
 
float TrackTimeDefAndQualityAlg::trackTimeResolution(
    const std::array<TrackTimeDefAndQualityAlg::Hit, s_hgtd_layers>& hits)
    const {
 
  float sum = 0.;
  for (const auto& hit : hits) {
    if (hit.m_isvalid) {
      sum += 1. / (hit.m_resolution * hit.m_resolution);
    }
  }
  return sum == 0. ? m_default_time_res.value()
                   : static_cast<float>(std::sqrt(1. / sum));
}
 
short TrackTimeDefAndQualityAlg::findLayerWithBadChi2(
    std::array<TrackTimeDefAndQualityAlg::Hit, s_hgtd_layers> hits) const {
  short remove_layer = -1;
  float local_min_chi2 = 999999;
  for (auto& hit : hits) {
    // "turn off" hits one after the other to test their impact on the chi2
    bool validbuff = hit.m_isvalid;
    hit.m_isvalid = false;
    float local_chi2 = calculateChi2(hits);
    hit.m_isvalid = validbuff;
    if (local_chi2 < local_min_chi2) {
      local_min_chi2 = local_chi2;
      remove_layer = hit.m_layer;
    }
  }
  return remove_layer;
}
 
void TrackTimeDefAndQualityAlg::setLayerAsInvalid(
    std::array<TrackTimeDefAndQualityAlg::Hit, s_hgtd_layers>& hits,
    short layer) const {
  for (auto& hit : hits) {
    if (hit.m_layer == layer) {
      hit.m_isvalid = false;
    }
  }
}
 
const Trk::TrackParameters*
TrackTimeDefAndQualityAlg::getLastHitOnTrack(const Trk::Track& track) const {
 
  const Trk::TrackStates* tsos =
      track.trackStateOnSurfaces();
  if (not tsos) {
    ATH_MSG_ERROR("Failed to retrieve track state on surfaces");
    return nullptr;
  }
  // loop over the associated hits in ITk in reverse order, since we want to
  // select the one closest to HGTD to start the extrapolation
  for (auto i = tsos->rbegin(); i != tsos->rend(); ++i) {
    const auto* curr_last_tsos = *i;
    if (not curr_last_tsos) {
      continue;
    }
    
    if (curr_last_tsos->type(Trk::TrackStateOnSurface::Measurement) and
        curr_last_tsos->trackParameters() and
        curr_last_tsos->measurementOnTrack()) {
      return curr_last_tsos->trackParameters();
    }
  }
  
  return nullptr;
}
 
std::pair<float, float> TrackTimeDefAndQualityAlg::getRadiusAndZ(const xAOD::TrackParticle& track_particle) const
{
  float radius = 0.f;
  float abs_z = 0.f;
 
  if (not m_doActs) {
    const Trk::Track* track = track_particle.track();
    if (not track) throw std::runtime_error("Cannot retrieve Trk track from Track Particle");
    const Trk::TrackParameters* last_hit_param = getLastHitOnTrack(*track);
    if (not last_hit_param) throw std::runtime_error("Cannot retrieve Trk track parameters from Trk track");
 
    radius = std::hypot(last_hit_param->position().x(),
            last_hit_param->position().y());
    abs_z = std::abs(last_hit_param->position().z());
  } else {
    // ACTS
    static const SG::ConstAccessor< ElementLink<ActsTrk::TrackContainer> > actsTrackLink("actsTrack");
    if (not actsTrackLink.isAvailable(track_particle)) throw std::runtime_error("Track particle does not have link to acts track");
 
    ElementLink<ActsTrk::TrackContainer> link_to_track = actsTrackLink(track_particle);
    if (not link_to_track.isValid()) throw std::runtime_error("Element link to acts track is not valid");
 
    std::optional<ActsTrk::TrackContainer::ConstTrackProxy> optional_track = *link_to_track;
    if (not optional_track.has_value()) throw std::runtime_error("Link to acts track has no value");
 
    const ActsTrk::TrackContainer::ConstTrackProxy& track = optional_track.value();
    const auto lastMeasurementState = Acts::findLastMeasurementState(track);
    const auto state = lastMeasurementState.value();
 
    const xAOD::UncalibratedMeasurement *cluster = ActsTrk::detail::xAODUncalibMeasCalibrator::unpack(state.getUncalibratedSourceLink());
    xAOD::UncalibMeasType clusterType = cluster->type();
 
    switch (clusterType) {
    case xAOD::UncalibMeasType::PixelClusterType:
      {
        auto glob = static_cast<const xAOD::PixelCluster*>(cluster)->globalPosition();
        radius = glob.perp();
        abs_z = std::abs(glob.z());
      }
      break;
    case xAOD::UncalibMeasType::StripClusterType:
      {
        auto glob = static_cast<const xAOD::StripCluster*>(cluster)->globalPosition();
        radius = glob.perp();
        abs_z = std::abs(glob.z());
      }
      break;
    case xAOD::UncalibMeasType::HGTDClusterType:
      // return some default numbers
      return std::make_pair(700, 3000);
    default:
      return std::make_pair(radius, abs_z);
    }; // switch
  } // acts
 
  return std::make_pair(radius, abs_z);
}
  
bool TrackTimeDefAndQualityAlg::lastHitIsOnLastSurface(
    const xAOD::TrackParticle& track_particle) const {
  auto [radius, abs_z] = getRadiusAndZ(track_particle);
 
  if (abs_z > 2700) {
    return true;
  }
  if (radius < 350 and abs_z > 2400) {
    return true;
  }
  // region 2
  if (radius > 205 and radius < 350 and abs_z > 2100) {
    return true;
  }
  // region 3
  if (radius < 220 and abs_z > 2200) {
    return true;
  }
 
  if (radius < 140 and abs_z > 1890) {
    return true;
  }
 
  return false;
}
 
} // namespace HGTD