src/Vertexing/ZScanVertexFinder.cpp

0001 // This file is part of the ACTS project.
0002 //
0003 // Copyright (C) 2016 CERN for the benefit of the ACTS project
0004 //
0005 // This Source Code Form is subject to the terms of the Mozilla Public
0006 // License, v. 2.0. If a copy of the MPL was not distributed with this
0007 // file, You can obtain one at https://mozilla.org/MPL/2.0/.
0008
0009 #include "Acts/Vertexing/ZScanVertexFinder.hpp"
0010
0011 Acts::ZScanVertexFinder::ZScanVertexFinder(const Config& cfg,
0012                                            std::unique_ptr<const Logger> logger)
0013     : m_cfg(cfg), m_logger(std::move(logger)) {
0014   if (!m_cfg.extractParameters.connected()) {
0015     throw std::invalid_argument(
0016         "ZScanVertexFinder: "
0017         "No track parameter extractor provided.");
0018   }
0019 }
0020
0021 Acts::Result<std::vector<Acts::Vertex>> Acts::ZScanVertexFinder::find(
0022     const std::vector<InputTrack>& trackVector,
0023     const VertexingOptions& vertexingOptions,
0024     IVertexFinder::State& /*state*/) const {
0025   double ZResult = 0.;
0026   // Prepare the vector of points, on which the 3d mode has later to be
0027   // calculated
0028   std::vector<std::pair<double, double>> zPositions;
0029
0030   for (const auto& iTrk : trackVector) {
0031     // Extract BoundTrackParameters from InputTrack_t object
0032     const BoundTrackParameters& params = m_cfg.extractParameters(iTrk);
0033
0034     std::pair<double, double> z0AndWeight;
0035     ImpactParametersAndSigma ipas;
0036     if (vertexingOptions.useConstraintInFit &&
0037         vertexingOptions.constraint.covariance()(0, 0) != 0) {
0038       auto estRes = m_cfg.ipEstimator.getImpactParameters(
0039           params, vertexingOptions.constraint, vertexingOptions.geoContext,
0040           vertexingOptions.magFieldContext);
0041       if (estRes.ok()) {
0042         ipas = *estRes;
0043       } else {
0044         return estRes.error();
0045       }
0046     }
0047
0048     if (ipas.sigmaD0 > 0) {
0049       // calculate z0
0050       z0AndWeight.first = ipas.z0 + vertexingOptions.constraint.position().z();
0051
0052       // calculate chi2 of IP
0053       double chi2IP = std::pow(ipas.d0 / ipas.sigmaD0, 2);
0054
0055       if (!m_cfg.disableAllWeights) {
0056         z0AndWeight.second =
0057             1. / (1. + std::exp((chi2IP - m_cfg.constraintcutoff) /
0058                                 m_cfg.constrainttemp));
0059         // overflow protection
0060         if (!std::isnormal(z0AndWeight.second)) {
0061           z0AndWeight.second = 0.;
0062         }
0063       } else {
0064         z0AndWeight.second = 1.;
0065       }
0066     } else {
0067       ACTS_DEBUG(
0068           "Unable to compute IP significance. "
0069           "Setting IP weight to 1.");
0070
0071       z0AndWeight.first = params.position(vertexingOptions.geoContext)[eZ];
0072       z0AndWeight.second = 1.;
0073     }
0074
0075     // apply pT weighting as/if configured
0076     if (!m_cfg.disableAllWeights && (m_cfg.usePt || m_cfg.useLogPt)) {
0077       double Pt =
0078           std::abs(1. / params.parameters()[BoundIndices::eBoundQOverP]) *
0079           std::sin(params.parameters()[BoundIndices::eBoundTheta]);
0080       if (m_cfg.usePt) {
0081         z0AndWeight.second *= std::pow(Pt, m_cfg.expPt);
0082       } else {
0083         z0AndWeight.second *=
0084             Pt > m_cfg.minPt ? std::log(Pt / m_cfg.minPt) : 0.;
0085       }
0086     }
0087
0088     if (z0AndWeight.second >= m_cfg.minWeight) {
0089       zPositions.push_back(z0AndWeight);
0090     }
0091   }  // end of loop over perigeeList
0092
0093   if (!zPositions.empty()) {
0094     auto res = m_cfg.mode1dFinder.getMode(zPositions);
0095     if (res.ok()) {
0096       ZResult = *res;
0097     } else {
0098       return res.error();
0099     }
0100
0101     ACTS_DEBUG("Resulting mean Z position found: " << ZResult);
0102   }
0103
0104   // constraint x()/y() equals 0 if no constraint
0105   Vector4 output(vertexingOptions.constraint.position().x(),
0106                  vertexingOptions.constraint.position().y(), ZResult,
0107                  vertexingOptions.constraint.time());
0108   Vertex vtxResult = Vertex(output);
0109
0110   return std::vector<Vertex>{vtxResult};
0111 }