AverageSimHits.hpp

Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file AverageSimHits.hpp

// This file is part of the Acts project.
//
// Copyright (C) 2020 CERN for the benefit of the Acts project
//
// This Source Code Form is subject to the terms of the Mozilla Public
// License, v. 2.0. If a copy of the MPL was not distributed with this
// file, You can obtain one at http://mozilla.org/MPL/2.0/.

#pragma once

#include "Acts/Definitions/Algebra.hpp"
#include "Acts/Definitions/Units.hpp"
#include "Acts/Surfaces/Surface.hpp"
#include "Acts/Utilities/Logger.hpp"
#include "ActsExamples/EventData/Index.hpp"
#include "ActsExamples/EventData/SimHit.hpp"
#include "ActsExamples/Utilities/Range.hpp"

#include <tuple>

namespace ActsExamples {

using HitSimHitsRange = Range<IndexMultimap<Index>::const_iterator>;

inline std::tuple<Acts::Vector2, Acts::Vector4, Acts::Vector3> averageSimHits(
    const Acts::GeometryContext& gCtx, const Acts::Surface& surface,
    const SimHitContainer& simHits, const HitSimHitsRange& hitSimHitsRange,
    const Acts::Logger& logger) {
  using namespace Acts::UnitLiterals;

  Acts::Vector2 avgLocal = Acts::Vector2::Zero();
  Acts::Vector4 avgPos4 = Acts::Vector4::Zero();
  Acts::Vector3 avgDir = Acts::Vector3::Zero();

  size_t n = 0u;
  for (auto [_, simHitIdx] : hitSimHitsRange) {
    n += 1u;

    // we assume that the indices are within valid ranges so we do not need to
    // check their validity again.
    const auto& simHit = *simHits.nth(simHitIdx);

    // transforming first to local positions and average that ensures that the
    // averaged position is still on the surface. the averaged global position
    // might not be on the surface anymore.
    auto result = surface.globalToLocal(gCtx, simHit.position(),
                                        simHit.direction(), 0.5_um);
    if (result.ok()) {
      avgLocal += result.value();
    } else {
      ACTS_WARNING("While averaging simhit, hit "
                   << simHitIdx << " is not on the corresponding surface "
                   << surface.geometryId() << "; use [0,0] as local position");
    }
    // global position should already be at the intersection. no need to perform
    // an additional intersection call.
    avgPos4 += simHit.fourPosition();
    avgDir += simHit.direction();
  }

  // only need to average if there are at least two inputs
  if (2u <= n) {
    double scale = 1.0 / n;
    avgLocal *= scale;
    avgPos4 *= scale;
    avgDir.normalize();
  }

  return {avgLocal, avgPos4, avgDir};
}

}  // namespace ActsExamples