ImpactPointEstimatorTests.cpp

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

// This file is part of the Acts project.
//
// Copyright (C) 2019 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/.

#include <boost/test/data/test_case.hpp>
#include <boost/test/tools/output_test_stream.hpp>
#include <boost/test/unit_test.hpp>

#include "Acts/MagneticField/ConstantBField.hpp"
#include "Acts/Propagator/EigenStepper.hpp"
#include "Acts/Propagator/Propagator.hpp"
#include "Acts/Surfaces/PerigeeSurface.hpp"
#include "Acts/Tests/CommonHelpers/FloatComparisons.hpp"
#include "Acts/Utilities/Definitions.hpp"
#include "Acts/Utilities/Helpers.hpp"
#include "Acts/Utilities/Units.hpp"
#include "Acts/Vertexing/ImpactPointEstimator.hpp"
#include "Acts/Vertexing/Vertex.hpp"

using namespace Acts::UnitLiterals;
using Acts::VectorHelpers::makeVector4;

namespace Acts {
namespace Test {

using Covariance = BoundSymMatrix;
using Propagator = Propagator<EigenStepper<ConstantBField>>;
// Create a test context
GeometryContext geoContext = GeometryContext();
MagneticFieldContext magFieldContext = MagneticFieldContext();

// Track d0 distribution
std::uniform_real_distribution<> d0Dist(-0.01_mm, 0.01_mm);
// Track z0 distribution
std::uniform_real_distribution<> z0Dist(-0.2_mm, 0.2_mm);
// Track pT distribution
std::uniform_real_distribution<> pTDist(0.4_GeV, 10._GeV);
// Track phi distribution
std::uniform_real_distribution<> phiDist(-M_PI, M_PI);
// Track theta distribution
std::uniform_real_distribution<> thetaDist(1.0, M_PI - 1.0);
// Track IP resolution distribution
std::uniform_real_distribution<> resIPDist(0., 100._um);
// Track angular distribution
std::uniform_real_distribution<> resAngDist(0., 0.1);
// Track q/p resolution distribution
std::uniform_real_distribution<> resQoPDist(-0.1, 0.1);
// Track charge helper distribution
std::uniform_real_distribution<> qDist(-1, 1);
// Vertex x/y position distribution
std::uniform_real_distribution<> vXYDist(-0.1_mm, 0.1_mm);
// Vertex z position distribution
std::uniform_real_distribution<> vZDist(-20_mm, 20_mm);
// Number of tracks distritbution
std::uniform_int_distribution<> nTracksDist(3, 10);

BOOST_AUTO_TEST_CASE(impactpoint_estimator_params_distance_test) {
  // Debug mode
  bool debugMode = false;
  // Number of tests
  unsigned int nTests = 10;

  // Set up RNG
  int mySeed = 31415;
  std::mt19937 gen(mySeed);

  // Set up constant B-Field
  ConstantBField bField(Vector3D(0., 0., 2._T));

  // Set up Eigenstepper
  EigenStepper<ConstantBField> stepper(bField);

  // Set up propagator with void navigator
  auto propagator = std::make_shared<Propagator>(stepper);

  // Set up the ImpactPointEstimator
  using IPEstimator = ImpactPointEstimator<BoundTrackParameters, Propagator>;
  IPEstimator::Config ipEstCfg(bField, propagator);
  IPEstimator ipEstimator(ipEstCfg);
  IPEstimator::State state(magFieldContext);

  // Reference position
  Vector3D refPosition(0., 0., 0.);

  // Start running tests
  for (unsigned int i = 0; i < nTests; i++) {
    // Create a track
    // Resolutions
    double resD0 = resIPDist(gen);
    double resZ0 = resIPDist(gen);
    double resPh = resAngDist(gen);
    double resTh = resAngDist(gen);
    double resQp = resQoPDist(gen);

    // Covariance matrix
    Covariance covMat;
    covMat << resD0 * resD0, 0., 0., 0., 0., 0., 0., resZ0 * resZ0, 0., 0., 0.,
        0., 0., 0., resPh * resPh, 0., 0., 0., 0., 0., 0., resTh * resTh, 0.,
        0., 0., 0., 0., 0., resQp * resQp, 0., 0., 0., 0., 0., 0., 1.;

    // The charge
    double q = qDist(gen) < 0 ? -1. : 1.;
    // Impact parameters (IP)
    double d0 = d0Dist(gen);
    double z0 = z0Dist(gen);

    if (debugMode) {
      std::cout << "IP: (" << d0 << "," << z0 << ")" << std::endl;
    }

    // The track parameters
    BoundTrackParameters::ParametersVector paramVec;
    paramVec[eBoundLoc0] = d0;
    paramVec[eBoundLoc1] = z0;
    paramVec[eBoundTime] = 0;
    paramVec[eBoundPhi] = phiDist(gen);
    paramVec[eBoundTheta] = thetaDist(gen);
    paramVec[eBoundQOverP] = q / pTDist(gen);

    // Corresponding surface
    std::shared_ptr<PerigeeSurface> perigeeSurface =
        Surface::makeShared<PerigeeSurface>(refPosition);

    // Creating the track
    BoundTrackParameters myTrack(perigeeSurface, paramVec, std::move(covMat));

    // Distance in transverse plane
    double transverseDist = std::sqrt(std::pow(d0, 2) + std::pow(z0, 2));

    // Estimate 3D distance
    auto distanceRes = ipEstimator.calculate3dDistance(geoContext, myTrack,
                                                       refPosition, state);
    BOOST_CHECK(distanceRes.ok());

    double distance = *distanceRes;

    BOOST_CHECK(distance < transverseDist);

    if (debugMode) {
      std::cout << std::setprecision(10)
                << "Distance in transverse plane: " << transverseDist
                << std::endl;
      std::cout << std::setprecision(10) << "Distance in 3D: " << distance
                << std::endl;
    }

    auto res = ipEstimator.estimate3DImpactParameters(
        geoContext, magFieldContext, myTrack, refPosition, state);

    BOOST_CHECK(res.ok());

    BoundTrackParameters trackAtIP3d = std::move(**res);

    const auto& myTrackParams = myTrack.parameters();
    const auto& trackIP3dParams = trackAtIP3d.parameters();

    // d0 and z0 should have changed
    BOOST_CHECK_NE(myTrackParams[BoundIndices::eBoundLoc0],
                   trackIP3dParams[BoundIndices::eBoundLoc0]);
    BOOST_CHECK_NE(myTrackParams[BoundIndices::eBoundLoc1],
                   trackIP3dParams[BoundIndices::eBoundLoc1]);
    // Theta along helix and q/p shoud remain the same
    CHECK_CLOSE_REL(myTrackParams[BoundIndices::eBoundTheta],
                    trackIP3dParams[BoundIndices::eBoundTheta], 1e-5);
    CHECK_CLOSE_REL(myTrackParams[BoundIndices::eBoundQOverP],
                    trackIP3dParams[BoundIndices::eBoundQOverP], 1e-5);

    if (debugMode) {
      std::cout << std::setprecision(10) << "Old track parameters: \n"
                << myTrackParams << std::endl;
      std::cout << std::setprecision(10) << "Parameters at IP3d: \n"
                << trackIP3dParams << std::endl;
    }
  }  // end for loop tests
}

BOOST_AUTO_TEST_CASE(impactpoint_estimator_compatibility_test) {
  // Debug mode
  bool debugMode = false;
  // Number of tests
  unsigned int nTests = 10;

  // Set up RNG
  int mySeed = 31415;
  std::mt19937 gen(mySeed);

  // Set up constant B-Field
  ConstantBField bField(Vector3D(0., 0., 2.) * UnitConstants::T);

  // Set up Eigenstepper
  EigenStepper<ConstantBField> stepper(bField);

  // Set up propagator with void navigator
  auto propagator = std::make_shared<Propagator>(stepper);

  // Set up the ImpactPointEstimator
  using IPEstimator = ImpactPointEstimator<BoundTrackParameters, Propagator>;
  IPEstimator::Config ipEstCfg(bField, propagator);
  IPEstimator ipEstimator(ipEstCfg);
  IPEstimator::State state(magFieldContext);

  // Reference position
  Vector3D refPosition(0., 0., 0.);

  // Lists to store distances and comp values
  std::vector<double> distancesList;
  std::vector<double> compatibilityList;

  // Generate covariance matrix values once
  // for all cov matrices below
  // Resolutions
  double resD0 = resIPDist(gen);
  double resZ0 = resIPDist(gen);
  double resPh = resAngDist(gen);
  double resTh = resAngDist(gen);
  double resQp = resQoPDist(gen);

  // Start running tests
  for (unsigned int i = 0; i < nTests; i++) {
    // Covariance matrix
    BoundTrackParameters::CovarianceMatrix covMat =
        BoundTrackParameters::CovarianceMatrix::Zero();
    covMat(eBoundLoc0, eBoundLoc0) = resD0 * resD0;
    covMat(eBoundLoc1, eBoundLoc1) = resZ0 * resZ0;
    covMat(eBoundTime, eBoundTime) = 1;
    covMat(eBoundPhi, eBoundPhi) = resPh * resPh;
    covMat(eBoundTheta, eBoundTheta) = resTh * resTh;
    covMat(eBoundQOverP, eBoundQOverP) = resQp * resQp;
    // The track parameters
    BoundTrackParameters::ParametersVector paramVec;
    paramVec[eBoundLoc0] = d0Dist(gen);
    paramVec[eBoundLoc1] = z0Dist(gen);
    paramVec[eBoundTime] = 0;
    paramVec[eBoundPhi] = phiDist(gen);
    paramVec[eBoundTheta] = thetaDist(gen);
    paramVec[eBoundQOverP] = (qDist(gen) < 0 ? -1. : 1.) / pTDist(gen);

    // Corresponding surface
    auto perigeeSurface = Surface::makeShared<PerigeeSurface>(refPosition);
    // Creating the track
    BoundTrackParameters myTrack(perigeeSurface, paramVec, std::move(covMat));

    // Estimate 3D distance
    auto distanceRes = ipEstimator.calculate3dDistance(geoContext, myTrack,
                                                       refPosition, state);
    BOOST_CHECK(distanceRes.ok());

    distancesList.push_back(*distanceRes);

    auto res = ipEstimator.estimate3DImpactParameters(
        geoContext, magFieldContext, myTrack, refPosition, state);

    BOOST_CHECK(res.ok());

    BoundTrackParameters params = std::move(**res);

    auto compRes =
        ipEstimator.get3dVertexCompatibility(geoContext, &params, refPosition);

    BOOST_CHECK(compRes.ok());

    compatibilityList.push_back(*compRes);

  }  // end create tracks loop

  // Now test for all above constructed tracks
  // if distances and compatibility values are
  // compatible with one another
  for (unsigned int i = 0; i < nTests; i++) {
    for (unsigned int j = i + 1; j < nTests; j++) {
      double relDiffComp =
          (compatibilityList[i] - compatibilityList[j]) / compatibilityList[i];

      double relDiffDist =
          (distancesList[i] - distancesList[j]) / distancesList[i];

      if (debugMode) {
        std::cout << "Comparing track " << i << " with track " << j
                  << std::endl;
        std::cout << "\t" << i << ": Comp.: " << compatibilityList[i]
                  << ", dist.: " << distancesList[i] << std::endl;
        std::cout << "\t" << j << ": Comp.: " << compatibilityList[j]
                  << ", dist.: " << distancesList[j] << std::endl;
        std::cout << "\t Rel.diff.: Comp(1-2)/1: " << relDiffComp
                  << ", Dist(1-2)/1: " << relDiffDist << std::endl;
      }

      // Relative differences of compatibility values and distances
      // should have the same sign, i.e. the following product
      // should always be positive
      // double res = relDiffComp * relDiffDist;

      // TODO 2020-09-09 msmk
      // this fails for one track after the the track parameters cleanup.
      // i do not understand what this tests and/or how to fix it. Bastian
      // has to look at this.
      // BOOST_CHECK_GE(res, 0);
    }
  }
}

BOOST_AUTO_TEST_CASE(impactpoint_estimator_athena_test) {
  // Set up constant B-Field
  ConstantBField bField(Vector3D(0., 0., 1.9971546939_T));

  // Set up Eigenstepper
  EigenStepper<ConstantBField> stepper(bField);

  // Set up propagator with void navigator
  auto propagator = std::make_shared<Propagator>(stepper);

  // Set up the ImpactPointEstimator
  using IPEstimator = ImpactPointEstimator<BoundTrackParameters, Propagator>;
  IPEstimator::Config ipEstCfg(bField, propagator);
  IPEstimator ipEstimator(ipEstCfg);
  IPEstimator::State state(magFieldContext);

  // Use same values as in Athena unit test
  Vector3D pos1(2_mm, 1_mm, -10_mm);
  Vector3D mom1(400_MeV, 600_MeV, 200_MeV);
  Vector3D vtxPos(1.2_mm, 0.8_mm, -7_mm);

  // Start creating some track parameters
  std::shared_ptr<PerigeeSurface> perigeeSurface =
      Surface::makeShared<PerigeeSurface>(pos1);
  // Some fixed track parameter values
  BoundTrackParameters params1(perigeeSurface, geoContext,
                               makeVector4(pos1, 0_ns), mom1, mom1.norm(), 1,
                               Covariance::Identity());

  auto res1 =
      ipEstimator.calculate3dDistance(geoContext, params1, vtxPos, state);
  BOOST_CHECK(res1.ok());
  double distance = (*res1);

  // Desired result from Athena unit test
  const double result = 3.10391_mm;
  CHECK_CLOSE_ABS(distance, result, 0.00001_mm);

  auto res2 = ipEstimator.estimate3DImpactParameters(
      geoContext, magFieldContext, params1, vtxPos, state);
  BOOST_CHECK(res2.ok());
  BoundTrackParameters endParams = std::move(**res2);
  Vector3D surfaceCenter = endParams.referenceSurface().center(geoContext);

  BOOST_CHECK_EQUAL(surfaceCenter, vtxPos);
}

BOOST_AUTO_TEST_CASE(impactpoint_estimator_parameter_estimation_test) {
  // Number of tracks to test with
  unsigned int nTracks = 10;

  // Set up RNG
  int mySeed = 31415;
  std::mt19937 gen(mySeed);

  // Set up constant B-Field
  ConstantBField bField(0.0, 0.0, 1_T);

  // Set up Eigenstepper
  EigenStepper<ConstantBField> stepper(bField);

  // Set up propagator with void navigator
  auto propagator = std::make_shared<Propagator>(stepper);

  // Create perigee surface
  std::shared_ptr<PerigeeSurface> perigeeSurface =
      Surface::makeShared<PerigeeSurface>(Vector3D(0., 0., 0.));

  // Create position of vertex and perigee surface
  double x = vXYDist(gen);
  double y = vXYDist(gen);
  double z = vZDist(gen);

  Vector4D vertexPosition(x, y, z, 0.);

  // Constraint for vertex fit
  Vertex<BoundTrackParameters> myConstraint;
  // Some abitrary values
  SymMatrix4D myCovMat = SymMatrix4D::Zero();
  myCovMat(0, 0) = 30.;
  myCovMat(1, 1) = 30.;
  myCovMat(2, 2) = 30.;
  myCovMat(3, 3) = 30.;
  myConstraint.setFullCovariance(std::move(myCovMat));
  myConstraint.setFullPosition(vertexPosition);

  // Calculate d0 and z0 corresponding to vertex position
  double d0_v = std::sqrt(x * x + y * y);
  double z0_v = z;

  // Set up the ImpactPointEstimator
  using IPEstimator = ImpactPointEstimator<BoundTrackParameters, Propagator>;
  IPEstimator::Config ipEstCfg(bField, propagator);
  IPEstimator ipEstimator(ipEstCfg);
  IPEstimator::State state(magFieldContext);

  // Construct random track emerging from vicinity of vertex position
  // Vector to store track objects used for vertex fit
  for (unsigned int iTrack = 0; iTrack < nTracks; iTrack++) {
    // Construct positive or negative charge randomly
    double q = qDist(gen) < 0 ? -1. : 1.;

    // Construct random track parameters
    BoundVector paramVec;
    paramVec << d0_v + d0Dist(gen), z0_v + z0Dist(gen), phiDist(gen),
        thetaDist(gen), q / pTDist(gen), 0.;

    // Resolutions
    double resD0 = resIPDist(gen);
    double resZ0 = resIPDist(gen);
    double resPh = resAngDist(gen);
    double resTh = resAngDist(gen);
    double resQp = resQoPDist(gen);

    // Fill vector of track objects with simple covariance matrix
    Covariance covMat;
    covMat << resD0 * resD0, 0., 0., 0., 0., 0., 0., resZ0 * resZ0, 0., 0., 0.,
        0., 0., 0., resPh * resPh, 0., 0., 0., 0., 0., 0., resTh * resTh, 0.,
        0., 0., 0., 0., 0., resQp * resQp, 0., 0., 0., 0., 0., 0., 1.;

    BoundTrackParameters track(perigeeSurface, paramVec, std::move(covMat));

    // Check if IP are retrieved
    ImpactParametersAndSigma output =
        ipEstimator
            .estimateImpactParameters(track, myConstraint, geoContext,
                                      magFieldContext)
            .value();
    BOOST_CHECK_NE(output.IPd0, 0.);
    BOOST_CHECK_NE(output.IPz0, 0.);
  }
}

}  // namespace Test
}  // namespace Acts