Back to home page

EIC code displayed by LXR
 
 

    

File indexing completed on 2025-12-16 09:25:04 

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 <boost/test/unit_test.hpp>
0010 
0011 #include "Acts/Seeding/CombinatorialSeedSolver.hpp"
0012 #include "Acts/Surfaces/detail/LineHelper.hpp"
0013 #include "Acts/Surfaces/detail/PlanarHelper.hpp"
0014 #include "Acts/Utilities/detail/Polynomials.hpp"
0015 
0016 #include <random>
0017 
0018 #include "StrawHitGeneratorHelper.hpp"
0019 
0020 using namespace Acts;
0021 using namespace Acts::detail;
0022 using namespace Acts::Experimental;
0023 using namespace Acts::Experimental::detail;
0024 using namespace Acts::UnitLiterals;
0025 using namespace Acts::PlanarHelper;
0026 using namespace Acts::detail::LineHelper;
0027 using namespace Acts::VectorHelpers;
0028 
0029 using Line_t = CompSpacePointAuxiliaries::Line_t;
0030 using Config_t = CompSpacePointAuxiliaries::Config;
0031 using ParIdx = CompSpacePointAuxiliaries::FitParIndex;
0032 using ResidualIdx = CompSpacePointAuxiliaries::ResidualIdx;
0033 
0034 using Vector = Line_t::Vector;
0035 using Pars_t = Line_t::ParamVector;
0036 
0037 constexpr auto logLvl = Logging::Level::INFO;
0038 constexpr std::size_t nEvents = 100;
0039 constexpr double tolerance = 1.e-3;
0040 
0041 ACTS_LOCAL_LOGGER(getDefaultLogger("StrawLineResidualTest", logLvl));
0042 
0043 namespace ActsTests {
0044 
0045 template <typename T>
0046 constexpr T absMax(const T& a, const T& b) {
0047   return std::max(Acts::abs(a), Acts::abs(b));
0048 }
0049 template <typename T, typename... argsT>
0050 constexpr T absMax(const T& a, const argsT... args) {
0051   return std::max(Acts::abs(a), absMax(args...));
0052 }
0053 double angle(const Vector& v1, const Vector& v2) {
0054   const double dots = Acts::abs(v1.dot(v2));
0055   return std::acos(std::clamp(dots, 0., 1.));
0056 }
0057 
0058 #define CHECK_CLOSE(a, b, tol) \
0059   BOOST_CHECK_LE(std::abs(a - b) / absMax(a, b, 1.), tol);
0060 
0061 #define COMPARE_VECTORS(v1, v2, tol)                             \
0062   {                                                              \
0063     CHECK_CLOSE(v1[toUnderlying(ResidualIdx::bending)],          \
0064                 v2[toUnderlying(ResidualIdx::bending)], tol);    \
0065     CHECK_CLOSE(v1[toUnderlying(ResidualIdx::nonBending)],       \
0066                 v2[toUnderlying(ResidualIdx::nonBending)], tol); \
0067     CHECK_CLOSE(v1[toUnderlying(ResidualIdx::time)],             \
0068                 v2[toUnderlying(ResidualIdx::time)], tol);       \
0069   }
0070 
0071 BOOST_AUTO_TEST_SUITE(SeedingSuite)
0072 
0073 void testResidual(const Pars_t& linePars, const FitTestSpacePoint& testPoint) {
0074   Config_t resCfg{};
0075   resCfg.useHessian = true;
0076   resCfg.calcAlongStrip = false;
0077   resCfg.parsToUse = {ParIdx::x0, ParIdx::y0, ParIdx::phi, ParIdx::theta};
0078   Line_t line{linePars};
0079 
0080   ACTS_INFO(__func__ << "() - " << __LINE__
0081                      << ":\n##############################################\n###"
0082                         "###########################################\n"
0083                      << "Residual test w.r.t. test line: "
0084                      << toString(line.position()) << ", "
0085                      << toString(line.direction())
0086                      << ", theta: " << (theta(line.direction()) / 1._degree)
0087                      << ", phi: " << (phi(line.direction())) / 1._degree);
0088 
0089   CompSpacePointAuxiliaries resCalc{resCfg,
0090                                     Acts::getDefaultLogger("testRes", logLvl)};
0091   resCalc.updateSpatialResidual(line, testPoint);
0092   ACTS_INFO(__func__ << "() - " << __LINE__ << ": Test residual w.r.t. "
0093                      << toString(testPoint));
0094   if (testPoint.isStraw()) {
0095     const double lineDist =
0096         signedDistance(testPoint.localPosition(), testPoint.sensorDirection(),
0097                        line.position(), line.direction());
0098 
0099     ACTS_INFO(__func__ << "() - " << __LINE__
0100                        << ": Residual: " << toString(resCalc.residual())
0101                        << ", line distance: " << lineDist
0102                        << ",  analog residual: "
0103                        << (lineDist - testPoint.driftRadius()));
0104     ///
0105     CHECK_CLOSE(resCalc.residual()[toUnderlying(ResidualIdx::bending)],
0106                 (lineDist - testPoint.driftRadius()), 1.e-8);
0107     if (testPoint.measuresLoc0()) {
0108       auto closePoint =
0109           lineIntersect(line.position(), line.direction(),
0110                         testPoint.localPosition(), testPoint.sensorDirection());
0111 
0112       ACTS_INFO(__func__
0113                 << "() - " << __LINE__ << ": Distance along wire "
0114                 << closePoint.pathLength() << ", residual: "
0115                 << resCalc.residual()[toUnderlying(ResidualIdx::nonBending)]);
0116 
0117       CHECK_CLOSE(resCalc.residual()[toUnderlying(ResidualIdx::nonBending)],
0118                   closePoint.pathLength(), 1.e-9);
0119     }
0120   } else {
0121     const Vector3& n{testPoint.planeNormal()};
0122     const Vector3 planeIsect = intersectPlane(line.position(), line.direction(),
0123                                               n, testPoint.localPosition())
0124                                    .position();
0125     const Vector3 delta = (planeIsect - testPoint.localPosition());
0126 
0127     ACTS_DEBUG(__func__ << "() - " << __LINE__
0128                         << ": Residual: " << toString(resCalc.residual())
0129                         << ", delta: " << toString(delta)
0130                         << ", <delta, n> =" << delta.dot(n));
0131 
0132     Acts::ActsSquareMatrix<3> coordTrf{Acts::ActsSquareMatrix<3>::Zero()};
0133     coordTrf.col(toUnderlying(ResidualIdx::bending)) =
0134         testPoint.measuresLoc1() ? testPoint.toNextSensor()
0135                                  : testPoint.sensorDirection();
0136     coordTrf.col(toUnderlying(ResidualIdx::nonBending)) =
0137         testPoint.measuresLoc1() ? testPoint.sensorDirection()
0138                                  : testPoint.toNextSensor();
0139     coordTrf.col(2) = n;
0140     ACTS_DEBUG(__func__ << "() - " << __LINE__ << ": Coordinate trf: \n"
0141                         << toString(coordTrf));
0142 
0143     const Vector3 trfDelta = coordTrf.inverse() * delta;
0144     ACTS_DEBUG(__func__ << "() - " << __LINE__
0145                         << ": Transformed delta: " << toString(trfDelta));
0146 
0147     if (testPoint.measuresLoc1()) {
0148       CHECK_CLOSE(trfDelta[toUnderlying(ResidualIdx::bending)],
0149                   resCalc.residual()[toUnderlying(ResidualIdx::bending)],
0150                   1.e-10);
0151     } else {
0152       CHECK_CLOSE(trfDelta[toUnderlying(ResidualIdx::bending)] *
0153                       static_cast<double>(resCfg.calcAlongStrip),
0154                   resCalc.residual()[toUnderlying(ResidualIdx::bending)],
0155                   1.e-10);
0156     }
0157     if (testPoint.measuresLoc0()) {
0158       CHECK_CLOSE(trfDelta[toUnderlying(ResidualIdx::nonBending)],
0159                   resCalc.residual()[toUnderlying(ResidualIdx::nonBending)],
0160                   1.e-10);
0161     } else {
0162       CHECK_CLOSE(trfDelta[toUnderlying(ResidualIdx::nonBending)] *
0163                       static_cast<double>(resCfg.calcAlongStrip),
0164                   resCalc.residual()[toUnderlying(ResidualIdx::nonBending)],
0165                   1.e-10);
0166     }
0167   }
0168 
0169   constexpr double h = 2.e-7;
0170   for (auto par : resCfg.parsToUse) {
0171     Pars_t lineParsUp{linePars};
0172     Pars_t lineParsDn{linePars};
0173     lineParsUp[toUnderlying(par)] += h;
0174     lineParsDn[toUnderlying(par)] -= h;
0175     const Line_t lineUp{lineParsUp};
0176     const Line_t lineDn{lineParsDn};
0177 
0178     CompSpacePointAuxiliaries resCalcUp{
0179         resCfg, Acts::getDefaultLogger("testResUp", logLvl)};
0180     CompSpacePointAuxiliaries resCalcDn{
0181         resCfg, Acts::getDefaultLogger("testResDn", logLvl)};
0182 
0183     resCalcUp.updateSpatialResidual(lineUp, testPoint);
0184     resCalcDn.updateSpatialResidual(lineDn, testPoint);
0185 
0186     const Vector numDeriv =
0187         (resCalcUp.residual() - resCalcDn.residual()) / (2. * h);
0188 
0189     ACTS_DEBUG(__func__ << "() - " << __LINE__ << ": Derivative test: "
0190                         << CompSpacePointAuxiliaries::parName(par)
0191                         << ", derivative: " << toString(resCalc.gradient(par))
0192                         << ",  numerical: " << toString(numDeriv));
0193 
0194     COMPARE_VECTORS(numDeriv, resCalc.gradient(par), tolerance);
0195     /// Next attempt to calculate the second derivative
0196     for (auto par1 : resCfg.parsToUse) {
0197       if (par1 > par) {
0198         break;
0199       }
0200       const Vector numDeriv1{
0201           (resCalcUp.gradient(par1) - resCalcDn.gradient(par1)) / (2. * h)};
0202       ACTS_DEBUG(__func__ << "() - " << __LINE__ << ": Second deriv ("
0203                           << CompSpacePointAuxiliaries::parName(par) << ", "
0204                           << CompSpacePointAuxiliaries::parName(par1)
0205                           << ") -- numerical: " << toString(numDeriv1)
0206                           << ", analytic: "
0207                           << toString(resCalc.hessian(par, par1)));
0208       COMPARE_VECTORS(numDeriv1, resCalc.hessian(par, par1), tolerance);
0209     }
0210   }
0211 }
0212 
0213 void testSeed(
0214     const std::array<std::shared_ptr<FitTestSpacePoint>, 4>& spacePoints,
0215     const std::array<double, 4>& truthDistances, const Vector3& truthPosZ0,
0216     const Vector3& truthDir) {
0217   const SquareMatrix2 bMatrix =
0218       CombinatorialSeedSolver::betaMatrix(spacePoints);
0219   if (std::abs(bMatrix.determinant()) < std::numeric_limits<float>::epsilon()) {
0220     ACTS_VERBOSE(__func__ << "() " << __LINE__ << ": Beta Matrix \n"
0221                           << toString(bMatrix)
0222                           << "\nhas zero determinant - skip the combination");
0223     return;
0224   }
0225   const std::array<double, 4> distsAlongStrip =
0226       CombinatorialSeedSolver::defineParameters(bMatrix, spacePoints);
0227   const auto [seedPos, seedDir] =
0228       CombinatorialSeedSolver::seedSolution(spacePoints, distsAlongStrip);
0229 
0230   ACTS_DEBUG(__func__ << "() " << __LINE__ << ": Reconstructed position "
0231                       << toString(seedPos) << " and direction "
0232                       << toString(seedDir));
0233   ACTS_DEBUG(__func__ << "() " << __LINE__ << ": Truth position "
0234                       << toString(truthPosZ0) << " and truth direction "
0235                       << toString(truthDir));
0236 
0237   // check the distances along the strips if they are the same
0238   for (std::size_t i = 0; i < truthDistances.size(); ++i) {
0239     CHECK_CLOSE(Acts::abs(distsAlongStrip[i] + truthDistances[i]), 0.,
0240                 tolerance);
0241   }
0242 
0243   COMPARE_VECTORS(seedPos, truthPosZ0, tolerance);
0244   // check the direction
0245   CHECK_CLOSE(std::abs(seedDir.dot(truthDir)), 1.,
0246               std::numeric_limits<float>::epsilon());
0247 }
0248 
0249 void timeStripResidualTest(const Pars_t& linePars, const double timeT0,
0250                            const FitTestSpacePoint& sp,
0251                            const Acts::Transform3& locToGlob) {
0252   using ResidualIdx = CompSpacePointAuxiliaries::ResidualIdx;
0253   Config_t resCfg{};
0254   resCfg.localToGlobal = locToGlob;
0255   resCfg.useHessian = true;
0256   resCfg.calcAlongStrip = true;
0257   resCfg.parsToUse = {ParIdx::x0, ParIdx::y0, ParIdx::phi, ParIdx::theta,
0258                       ParIdx::t0};
0259   Line_t line{linePars};
0260 
0261   ACTS_DEBUG(__func__ << "() - " << __LINE__ << ": Residual test for line: "
0262                       << toString(line.position()) << ", "
0263                       << toString(line.direction()) << " with t0: " << timeT0
0264                       << ".");
0265 
0266   CompSpacePointAuxiliaries resCalc{resCfg,
0267                                     Acts::getDefaultLogger("timeRes", logLvl)};
0268   resCalc.updateFullResidual(line, timeT0, sp);
0269 
0270   const Vector3 planeIsect =
0271       intersectPlane(line.position(), line.direction(), sp.planeNormal(),
0272                      sp.localPosition())
0273           .position();
0274   const double ToF =
0275       (locToGlob * planeIsect).norm() / Acts::PhysicalConstants::c + timeT0;
0276 
0277   CHECK_CLOSE(resCalc.residual()[toUnderlying(ResidualIdx::time)],
0278               (sp.time() - ToF), 1.e-10);
0279   ACTS_DEBUG(__func__ << "() - " << __LINE__ << ": Time of flight: " << ToF
0280                       << ", measured time : " << sp.time() << " --> residual: "
0281                       << (sp.time() - ToF) << " vs. from calculator: "
0282                       << resCalc.residual()[toUnderlying(ResidualIdx::time)]
0283                       << ".");
0284   constexpr double h = 5.e-9;
0285 
0286   Line_t lineUp{}, lineDn{};
0287   for (const auto partial : resCfg.parsToUse) {
0288     CompSpacePointAuxiliaries resCalcUp{
0289         resCfg, Acts::getDefaultLogger("timeResUp", logLvl)};
0290     CompSpacePointAuxiliaries resCalcDn{
0291         resCfg, Acts::getDefaultLogger("timeResDn", logLvl)};
0292 
0293     Pars_t lineParsUp{linePars}, lineParsDn{linePars};
0294 
0295     if (partial != ParIdx::t0) {
0296       lineParsUp[toUnderlying(partial)] += h;
0297       lineParsDn[toUnderlying(partial)] -= h;
0298       lineUp.updateParameters(lineParsUp);
0299       lineDn.updateParameters(lineParsDn);
0300       resCalcUp.updateFullResidual(lineUp, timeT0, sp);
0301       resCalcDn.updateFullResidual(lineDn, timeT0, sp);
0302     } else {
0303       lineUp.updateParameters(lineParsUp);
0304       lineDn.updateParameters(lineParsDn);
0305       resCalcUp.updateFullResidual(line, timeT0 + h, sp);
0306       resCalcDn.updateFullResidual(line, timeT0 - h, sp);
0307     }
0308 
0309     const Vector numDeriv =
0310         (resCalcUp.residual() - resCalcDn.residual()) / (2. * h);
0311 
0312     ACTS_DEBUG(__func__ << "() - " << __LINE__ << ": Derivative test: "
0313                         << CompSpacePointAuxiliaries::parName(partial)
0314                         << ", derivative: "
0315                         << toString(resCalc.gradient(partial))
0316                         << ",  numerical: " << toString(numDeriv));
0317     COMPARE_VECTORS(numDeriv, resCalc.gradient(partial), tolerance);
0318     for (const auto partial2 : resCfg.parsToUse) {
0319       const Vector numDeriv1{
0320           (resCalcUp.gradient(partial2) - resCalcDn.gradient(partial2)) /
0321           (2. * h)};
0322       ACTS_DEBUG(__func__ << "() - " << __LINE__ << ": Second deriv ("
0323                           << CompSpacePointAuxiliaries::parName(partial) << ", "
0324                           << CompSpacePointAuxiliaries::parName(partial2)
0325                           << ") -- numerical: " << toString(numDeriv1)
0326                           << ", analytic: "
0327                           << toString(resCalc.hessian(partial, partial2)));
0328       COMPARE_VECTORS(numDeriv1, resCalc.hessian(partial, partial2), tolerance);
0329     }
0330   }
0331 }
0332 BOOST_AUTO_TEST_CASE(StrawDriftTimeCase) {
0333   ACTS_INFO("Calibration straw point test ");
0334 
0335   constexpr double recordedT = 15._ns;
0336   constexpr std::array<double, 3> rtCoeffs{0., 75. / 1_ns,
0337                                            5000. / Acts::pow(1_ns, 2)};
0338 
0339   Config_t resCfg{};
0340   resCfg.useHessian = true;
0341   resCfg.calcAlongStrip = true;
0342   resCfg.parsToUse = {ParIdx::x0, ParIdx::phi, ParIdx::y0, ParIdx::theta,
0343                       ParIdx::t0};
0344 
0345   auto makeCalculator = [&rtCoeffs, &resCfg](
0346                             const std::string& calcName, const Pars_t& linePars,
0347                             const Vector& pos, const Vector& dir,
0348                             const double t0) {
0349     Line_t line{linePars};
0350     ACTS_DEBUG(__func__ << "() " << __LINE__ << ": Calculate residual w.r.t. "
0351                         << toString(line.position()) << ", "
0352                         << toString(line.direction()));
0353     auto isectP = lineIntersect(pos, dir, line.position(), line.direction());
0354 
0355     const double driftT = recordedT -
0356                           (resCfg.localToGlobal * isectP.position()).norm() /
0357                               PhysicalConstants::c -
0358                           t0;
0359 
0360     const double driftR = polynomialSum(driftT, rtCoeffs);
0361     const double driftV = derivativeSum<1>(driftT, rtCoeffs);
0362     const double driftA = derivativeSum<2>(driftT, rtCoeffs);
0363     ACTS_DEBUG(__func__ << "() " << __LINE__ << ": Create new space point @ "
0364                         << toString(pos) << ", " << toString(dir)
0365                         << ", using t0: " << t0 / 1._ns << " --> drfitT: "
0366                         << driftT / 1._ns << " --> driftR: " << driftR
0367                         << ", driftV: " << driftV << ", "
0368                         << "driftA: " << driftA);
0369 
0370     CompSpacePointAuxiliaries resCalc{resCfg,
0371                                       Acts::getDefaultLogger(calcName, logLvl)};
0372 
0373     resCalc.updateFullResidual(
0374         line, t0,
0375         FitTestSpacePoint{pos, dir, driftR, SpCalibrator::driftUncert(driftR)},
0376         driftV, driftA);
0377     return resCalc;
0378   };
0379 
0380   auto testTimingResidual = [&makeCalculator, &resCfg](
0381                                 const Pars_t& linePars, const Vector& wPos,
0382                                 const Vector& wDir, const double t0) {
0383     auto resCalc = makeCalculator("strawT0Res", linePars, wPos, wDir, t0);
0384     ACTS_DEBUG(__func__ << "() - " << __LINE__
0385                         << ": Residual: " << toString(resCalc.residual()));
0386     for (const auto partial : resCfg.parsToUse) {
0387       ACTS_DEBUG(__func__ << "() - " << __LINE__ << ": *** partial: "
0388                           << CompSpacePointAuxiliaries::parName(partial) << " "
0389                           << toString(resCalc.gradient(partial)));
0390     }
0391     constexpr double h = 1.e-7;
0392     for (const auto partial : resCfg.parsToUse) {
0393       Pars_t lineParsUp{linePars}, lineParsDn{linePars};
0394       double t0Up{t0}, t0Dn{t0};
0395       if (partial != ParIdx::t0) {
0396         lineParsUp[toUnderlying(partial)] += h;
0397         lineParsDn[toUnderlying(partial)] -= h;
0398       } else {
0399         t0Up += h;
0400         t0Dn -= h;
0401       }
0402       auto resCalcUp =
0403           makeCalculator("strawT0ResUp", lineParsUp, wPos, wDir, t0Up);
0404       auto resCalcDn =
0405           makeCalculator("strawT0ResDn", lineParsDn, wPos, wDir, t0Dn);
0406 
0407       const Vector numDeriv =
0408           (resCalcUp.residual() - resCalcDn.residual()) / (2. * h);
0409       ACTS_DEBUG(__func__ << "() - " << __LINE__ << ": Partial "
0410                           << CompSpacePointAuxiliaries::parName(partial)
0411                           << " --  numerical: " << toString(numDeriv)
0412                           << ", analytical: "
0413                           << toString(resCalc.gradient(partial)));
0414 
0415       COMPARE_VECTORS(numDeriv, resCalc.gradient(partial), tolerance);
0416       for (const auto partial2 : resCfg.parsToUse) {
0417         const Vector numDeriv1{
0418             (resCalcUp.gradient(partial2) - resCalcDn.gradient(partial2)) /
0419             (2. * h)};
0420         ACTS_DEBUG(__func__ << "() - " << __LINE__ << ": Second deriv ("
0421                             << CompSpacePointAuxiliaries::parName(partial)
0422                             << ", "
0423                             << CompSpacePointAuxiliaries::parName(partial2)
0424                             << ") -- numerical: " << toString(numDeriv1)
0425                             << ", analytic: "
0426                             << toString(resCalc.hessian(partial, partial2)));
0427         COMPARE_VECTORS(numDeriv1, resCalc.hessian(partial, partial2),
0428                         tolerance);
0429       }
0430     }
0431   };
0432 
0433   RandomEngine rndEngine{4711};
0434 
0435   const Vector wirePos{100._cm, 50._cm, 30._cm};
0436   for (std::size_t e = 0; e < nEvents; ++e) {
0437     break;
0438     ACTS_DEBUG(__func__ << "() - " << __LINE__ << ": Run test event: " << e);
0439     resCfg.localToGlobal = Acts::Transform3::Identity();
0440     Line_t line{generateLine(rndEngine, logger())};
0441     const double t0 = uniform{0_ns, 50._ns}(rndEngine);
0442 
0443     testTimingResidual(line.parameters(), wirePos, Vector::UnitX(), t0);
0444     testTimingResidual(line.parameters(), wirePos,
0445                        Vector{1., 1., 0.}.normalized(), t0);
0446     resCfg.localToGlobal.translation() =
0447         Vector{uniform{-10._cm, 10._cm}(rndEngine),
0448                uniform{-20._cm, 20._cm}(rndEngine),
0449                uniform{-30._cm, 30._cm}(rndEngine)};
0450     testTimingResidual(line.parameters(), wirePos, Vector::UnitX(), t0);
0451     testTimingResidual(line.parameters(), wirePos,
0452                        Vector{1., 1., 0.}.normalized(), t0);
0453 
0454     //// Next test the displacement
0455     resCfg.localToGlobal *=
0456         Acts::AngleAxis3{uniform{-30._degree, 30._degree}(rndEngine),
0457                          makeDirectionFromPhiTheta(
0458                              uniform{-30._degree, 30._degree}(rndEngine),
0459                              uniform{-45._degree, -35._degree}(rndEngine))};
0460     testTimingResidual(line.parameters(), wirePos, Vector::UnitX(), t0);
0461     testTimingResidual(line.parameters(), wirePos,
0462                        Vector{1., 1., 0.}.normalized(), t0);
0463   }
0464 }
0465 BOOST_AUTO_TEST_CASE(WireResidualTest) {
0466   RandomEngine rndEngine{2525};
0467   ACTS_INFO("Run WireResidualTest");
0468   const Vector wirePos{100._cm, 50._cm, 30._cm};
0469   const Vector wireDir1{Vector::UnitX()};
0470   const Vector wireDir2{makeDirectionFromPhiTheta(10._degree, 30._degree)};
0471 
0472   for (std::size_t e = 0; e < nEvents; ++e) {
0473     ACTS_DEBUG(__func__ << "() - " << __LINE__ << ": Run test event: " << e);
0474     const Line_t line{generateLine(rndEngine, logger())};
0475     // Generate the first test measurement
0476     constexpr double R = 10._cm;
0477     const double dR = SpCalibrator::driftUncert(R);
0478     constexpr double uncertWire = 1._cm;
0479     testResidual(line.parameters(),
0480                  FitTestSpacePoint{wirePos, wireDir1, R, dR});
0481 
0482     testResidual(line.parameters(),
0483                  FitTestSpacePoint{wirePos, wireDir2, R, dR});
0484 
0485     auto isect =
0486         lineIntersect(line.position(), line.direction(), wirePos, wireDir1);
0487 
0488     testResidual(line.parameters(),
0489                  FitTestSpacePoint{isect.position() + uniform{-50._cm, 50._cm}(
0490                                                           rndEngine)*wireDir1,
0491                                    wireDir1, R, dR, uncertWire});
0492 
0493     isect = lineIntersect(line.position(), line.direction(), wirePos, wireDir2);
0494     testResidual(line.parameters(),
0495                  FitTestSpacePoint{isect.position() + uniform{-50._cm, 50._cm}(
0496                                                           rndEngine)*wireDir2,
0497                                    wireDir2, R, dR, uncertWire});
0498   }
0499 }
0500 BOOST_AUTO_TEST_CASE(StripResidual) {
0501   ACTS_INFO("Run StripResidualTest");
0502   RandomEngine rndEngine{2505};
0503   const Vector stripPos{75._cm, -75._cm, 100._cm};
0504   const Vector b1{makeDirectionFromPhiTheta(90._degree, 90._degree)};
0505   const Vector b2{makeDirectionFromPhiTheta(0._degree, 90_degree)};
0506   for (std::size_t e = 0; e < nEvents; ++e) {
0507     ACTS_DEBUG(__func__ << "() - " << __LINE__ << ": Run test event: " << e);
0508     Pars_t linePars{generateLine(rndEngine, logger()).parameters()};
0509 
0510     /// Test solely the residual along the bending direction
0511     testResidual(linePars, FitTestSpacePoint{stripPos, b1, b2, 0._cm, 1._cm});
0512     /// Test solely the residual along the non-bending direction
0513     testResidual(linePars, FitTestSpacePoint{stripPos, b2, b1, 1._cm, 0._cm});
0514     /// Test the combined residual
0515     testResidual(linePars, FitTestSpacePoint{stripPos, b1, b2, 1._cm, 1._cm});
0516     const Vector b3 = makeDirectionFromPhiTheta(30._degree, 90._degree);
0517     const Vector b4 = makeDirectionFromPhiTheta(60._degree, 90._degree);
0518 
0519     testResidual(linePars, FitTestSpacePoint{stripPos, b3, b4, 1._cm, 1._cm});
0520   }
0521 }
0522 
0523 BOOST_AUTO_TEST_CASE(TimeStripResidual) {
0524   Pars_t linePars{};
0525   linePars[toUnderlying(ParIdx::phi)] = 60._degree;
0526   linePars[toUnderlying(ParIdx::theta)] = 45_degree;
0527 
0528   Acts::Transform3 locToGlob{Acts::Transform3::Identity()};
0529 
0530   const Vector pos{75._cm, -75._cm, 100._cm};
0531   const Vector pos1{75._cm, -75._cm, 200._cm};
0532 
0533   const Vector b1{makeDirectionFromPhiTheta(30._degree, 90._degree)};
0534   const Vector b2{makeDirectionFromPhiTheta(60._degree, 90._degree)};
0535   const Vector b3{makeDirectionFromPhiTheta(00._degree, 90._degree)};
0536   const Vector b4{makeDirectionFromPhiTheta(60._degree, 75._degree)};
0537   const std::array cov{Acts::square(10._cm), Acts::square(10._cm),
0538                        Acts::square(1._ns)};
0539   FitTestSpacePoint p1{pos, b1, b2, 15._ns, cov};
0540 
0541   FitTestSpacePoint p2{pos1, b3, b4, 15._ns, cov};
0542 
0543   timeStripResidualTest(linePars, 10., p1, locToGlob);
0544 
0545   timeStripResidualTest(linePars, 10., p2, locToGlob);
0546 
0547   locToGlob.translation() = Vector{75._cm, -75._cm, -35._cm};
0548   timeStripResidualTest(linePars, 10., p1, locToGlob);
0549 
0550   timeStripResidualTest(linePars, 10., p2, locToGlob);
0551   locToGlob *= Acts::AngleAxis3{
0552       30_degree, makeDirectionFromPhiTheta(30_degree, -45_degree)};
0553   timeStripResidualTest(linePars, 10., p1, locToGlob);
0554 
0555   timeStripResidualTest(linePars, 10., p2, locToGlob);
0556 }
0557 
0558 BOOST_AUTO_TEST_CASE(ChiSqEvaluation) {
0559   Pars_t linePars{};
0560   linePars[toUnderlying(ParIdx::phi)] = 30._degree;
0561   linePars[toUnderlying(ParIdx::theta)] = 75_degree;
0562   linePars[toUnderlying(ParIdx::x0)] = 10._cm;
0563   linePars[toUnderlying(ParIdx::y0)] = -10_cm;
0564 
0565   Config_t resCfg{};
0566   resCfg.useHessian = true;
0567   resCfg.calcAlongStrip = true;
0568   resCfg.parsToUse = {ParIdx::x0, ParIdx::phi, ParIdx::y0, ParIdx::theta,
0569                       ParIdx::t0};
0570 
0571   Line_t line{linePars};
0572 
0573   CompSpacePointAuxiliaries resCalc{resCfg,
0574                                     Acts::getDefaultLogger("testRes", logLvl)};
0575 
0576   const double t0{1._ns};
0577 
0578   auto testChi2 = [&line, &t0, &resCalc,
0579                    &resCfg](const FitTestSpacePoint& testMe) {
0580     using ChiSq_t = CompSpacePointAuxiliaries::ChiSqWithDerivatives;
0581     ChiSq_t chi2{};
0582     ACTS_DEBUG(__func__ << "() " << __LINE__ << ": Test space point "
0583                         << toString(testMe));
0584 
0585     resCalc.updateFullResidual(line, t0, testMe);
0586 
0587     resCalc.updateChiSq(chi2, testMe.covariance());
0588     resCalc.symmetrizeHessian(chi2);
0589     constexpr auto N = 5u;
0590     for (std::size_t i = 1; i < N; ++i) {
0591       for (std::size_t j = 0; j <= i; ++j) {
0592         BOOST_CHECK_EQUAL(chi2.hessian(i, j), chi2.hessian(j, i));
0593       }
0594     }
0595 
0596     ACTS_DEBUG(__func__ << "() - " << __LINE__
0597                         << ": Calculated chi2: " << chi2.chi2 << ", Gradient: "
0598                         << toString(chi2.gradient) << ", Hessian: \n"
0599                         << toString(chi2.hessian)
0600                         << "\ndeterminant: " << chi2.hessian.determinant());
0601     CHECK_CLOSE(CompSpacePointAuxiliaries::chi2Term(line, resCfg.localToGlobal,
0602                                                     t0, testMe),
0603                 chi2.chi2, 1.e-10);
0604 
0605     constexpr double h = 1.e-7;
0606     for (const auto par : resCfg.parsToUse) {
0607       Pars_t lineParsUp{line.parameters()};
0608       Pars_t lineParsDn{line.parameters()};
0609 
0610       const auto dIdx = toUnderlying(par);
0611       double t0Up{t0}, t0Dn{t0};
0612       if (par != ParIdx::t0) {
0613         lineParsUp[dIdx] += h;
0614         lineParsDn[dIdx] -= h;
0615       } else {
0616         t0Up += h;
0617         t0Dn -= h;
0618       }
0619       ChiSq_t chi2Up{}, chi2Dn{};
0620       /// Calculate the updated chi2
0621       line.updateParameters(lineParsUp);
0622       resCalc.updateFullResidual(line, t0Up, testMe);
0623       resCalc.updateChiSq(chi2Up, testMe.covariance());
0624 
0625       line.updateParameters(lineParsDn);
0626       resCalc.updateFullResidual(line, t0Dn, testMe);
0627       resCalc.updateChiSq(chi2Dn, testMe.covariance());
0628 
0629       const double anaDeriv = chi2.gradient[dIdx];
0630       const double numDeriv = (chi2Up.chi2 - chi2Dn.chi2) / (2. * h);
0631 
0632       ACTS_DEBUG(__func__ << "() - " << __LINE__ << ": Partial "
0633                           << CompSpacePointAuxiliaries::parName(par)
0634                           << " --  numerical: " << numDeriv
0635                           << ", analytical: " << anaDeriv);
0636       CHECK_CLOSE(numDeriv, anaDeriv, tolerance);
0637       for (const auto par2 : resCfg.parsToUse) {
0638         if (par2 > par) {
0639           break;
0640         }
0641         const auto dIdx2 = toUnderlying(par2);
0642         const double anaHess = chi2.hessian(dIdx, dIdx2);
0643         const double numHess =
0644             (chi2Up.gradient[dIdx2] - chi2Dn.gradient[dIdx2]) / (2. * h);
0645 
0646         ACTS_DEBUG(__func__ << "() - " << __LINE__ << ": Second deriv ("
0647                             << CompSpacePointAuxiliaries::parName(par) << ", "
0648                             << CompSpacePointAuxiliaries::parName(par2)
0649                             << ") -- numerical: " << numHess
0650                             << ", analytic: " << anaHess);
0651         CHECK_CLOSE(anaHess, numHess, tolerance);
0652       }
0653     }
0654   };
0655 
0656   /// Test orthogonal strips
0657   testChi2(FitTestSpacePoint{
0658       line.point(20._cm) + 5._cm * line.direction().cross(Vector::UnitX()),
0659       makeDirectionFromPhiTheta(0_degree, 90._degree),
0660       makeDirectionFromPhiTheta(90._degree, 0._degree),
0661       15._ns,
0662       {Acts::pow(5._cm, 2), Acts::pow(10._cm, 2), Acts::pow(1._ns, 2)}});
0663 
0664   /// Test strips with stereo
0665   testChi2(FitTestSpacePoint{
0666       line.point(20._cm) + 5._cm * line.direction().cross(Vector::UnitX()),
0667       makeDirectionFromPhiTheta(0_degree, 45._degree),
0668       makeDirectionFromPhiTheta(60._degree, 0._degree),
0669       15._ns,
0670       {Acts::pow(5._cm, 2), Acts::pow(10._cm, 2), Acts::pow(1._ns, 2)}});
0671   //// Test ordinary straws
0672   testChi2(FitTestSpacePoint{
0673       line.point(20._cm) + 5._cm * line.direction().cross(Vector::UnitX()),
0674       Vector3::UnitX(), 5._cm, 10._mm});
0675 
0676   /// Test straws with information on the position along the
0677   /// tube
0678   testChi2(FitTestSpacePoint{
0679       line.point(20._cm) + 5._cm * line.direction().cross(Vector::UnitX()) +
0680           4._cm * Vector::UnitX(),
0681       Vector3::UnitX(), 5._cm, 0.5_cm});
0682 }
0683 
0684 BOOST_AUTO_TEST_CASE(CombinatorialSeedSolverStripsTest) {
0685   ACTS_INFO("Run Combinatorial seed test");
0686   RandomEngine rndEngine{23568};
0687   constexpr std::size_t nStrips = 8;
0688 
0689   const std::array<Vector3, nStrips> stripDirections{
0690       Vector3::UnitX(),
0691       Vector3::UnitX(),
0692       makeDirectionFromPhiTheta(1.5_degree, 90_degree),
0693       makeDirectionFromPhiTheta(-1.5_degree, 90_degree),
0694       makeDirectionFromPhiTheta(1.5_degree, 90_degree),
0695       makeDirectionFromPhiTheta(-1.5_degree, 90_degree),
0696       Vector3::UnitX(),
0697       Vector3::UnitX()};
0698 
0699   constexpr std::array<double, nStrips> distancesZ{-20., -264., 0.,  300.,
0700                                                    350,  370.,  400, 450.};
0701 
0702   std::array<double, nStrips> parameters{};
0703   std::array<Vector3, nStrips> intersections{};
0704 
0705   // pseudo track initialization
0706 
0707   for (std::size_t i = 0; i < nEvents; ++i) {
0708     ACTS_DEBUG(__func__ << "() " << __LINE__ << " - Combinatorial Seed test - "
0709                         << "Processing Event: " << i);
0710     // update pseudo track parameters with random values
0711     Line_t line{generateLine(rndEngine, logger())};
0712 
0713     const Vector3& muonPos = line.position();
0714     const Vector3& muonDir = line.direction();
0715 
0716     std::array<std::shared_ptr<FitTestSpacePoint>, nStrips> spacePoints{};
0717     for (std::size_t layer = 0; layer < nStrips; ++layer) {
0718       auto intersection = PlanarHelper::intersectPlane(
0719           muonPos, muonDir, Vector3::UnitZ(), distancesZ[layer]);
0720       intersections[layer] = intersection.position();
0721       // where the hit is along the strip
0722       auto alongStrip = PlanarHelper::intersectPlane(
0723           intersections[layer], stripDirections[layer], Vector3::UnitX(), 0.);
0724 
0725       parameters[layer] = alongStrip.pathLength();
0726 
0727       Vector3 stripPos =
0728           intersections[layer] + parameters[layer] * stripDirections[layer];
0729 
0730       spacePoints[layer] = std::make_shared<FitTestSpacePoint>(
0731           stripPos, stripDirections[layer], Vector3::UnitZ(), 1._cm, 0._cm);
0732     }
0733 
0734     // let's pick four strips on fours different layers
0735     // check combinatorics
0736     std::array<std::shared_ptr<FitTestSpacePoint>, 4> seedSpacePoints{};
0737     for (unsigned int l = 0; l < distancesZ.size() - 3; ++l) {
0738       seedSpacePoints[0] = spacePoints[l];
0739       for (unsigned int k = l + 1; k < distancesZ.size() - 2; ++k) {
0740         seedSpacePoints[1] = spacePoints[k];
0741         for (unsigned int m = k + 1; m < distancesZ.size() - 1; ++m) {
0742           seedSpacePoints[2] = spacePoints[m];
0743           for (unsigned int n = m + 1; n < distancesZ.size(); ++n) {
0744             seedSpacePoints[3] = spacePoints[n];
0745 
0746             const std::array<double, 4> truthDistances = {
0747                 parameters[l], parameters[k], parameters[m], parameters[n]};
0748 
0749             testSeed(seedSpacePoints, truthDistances, muonPos, muonDir);
0750           }
0751         }
0752       }
0753     }
0754   }
0755 }
0756 
0757 }  // namespace ActsTests
0758 BOOST_AUTO_TEST_SUITE_END()
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated by the 2.3.7 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!