Back to home page

sPhenix code displayed by LXR
 
 

    

File indexing completed on 2026-07-16 08:08:51 

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/Definitions/Units.hpp"
0012 #include "Acts/Seeding/detail/FastStrawLineFitter.hpp"
0013 #include "Acts/Surfaces/detail/LineHelper.hpp"
0014 #include "Acts/Surfaces/detail/PlanarHelper.hpp"
0015 #include "Acts/Utilities/StringHelpers.hpp"
0016 
0017 #include <format>
0018 #include <random>
0019 
0020 #include "TFile.h"
0021 #include "TTree.h"
0022 
0023 using namespace Acts;
0024 using namespace Acts::Experimental;
0025 using namespace Acts::Experimental::detail;
0026 using namespace Acts::UnitLiterals;
0027 using RandomEngine = std::mt19937;
0028 using uniform_t = std::uniform_real_distribution<double>;
0029 using gauss_t = std::normal_distribution<double>;
0030 
0031 constexpr std::size_t nTrials = 1;
0032 constexpr auto logLvl = Logging::Level::INFO;
0033 
0034 namespace ActsTests {
0035 
0036 ACTS_LOCAL_LOGGER(getDefaultLogger("FastStrawLineFitTests", logLvl));
0037 
0038 class StrawTestPoint;
0039 using TestStrawCont_t = std::vector<std::unique_ptr<StrawTestPoint>>;
0040 using Line_t = CompSpacePointAuxiliaries::Line_t;
0041 using ResidualIdx = FastStrawLineFitter::ResidualIdx;
0042 
0043 /// @brief Outstream operator for a vector
0044 template <typename T>
0045 std::ostream& operator<<(std::ostream& ostr, const std::vector<T>& v) {
0046   ostr << "[";
0047   for (std::size_t i = 0; i < v.size(); ++i) {
0048     ostr << v[i];
0049     if (i + 1 != v.size()) {
0050       ostr << ", ";
0051     }
0052   }
0053   ostr << "]";
0054   return ostr;
0055 }
0056 /// @brief Converts an ACTS time into SI [ns]
0057 constexpr double inNanoS(const double x) {
0058   return x / 1._ns;
0059 }
0060 
0061 class StrawTestPoint {
0062  public:
0063   StrawTestPoint(const Vector3& pos, const double driftR,
0064                  const double driftRUncert)
0065       : m_pos{pos}, m_driftR{abs(driftR)} {
0066     m_cov[toUnderlying(ResidualIdx::bending)] = square(driftRUncert);
0067   }
0068 
0069   StrawTestPoint(const Vector3& pos, const double uncert)
0070       : m_pos{pos}, m_isStraw{false} {
0071     m_cov[toUnderlying(ResidualIdx::bending)] = square(uncert);
0072   }
0073   /// @brief Straw tube's direction
0074   const Vector3& localPosition() const { return m_pos; }
0075   /// @brief Wire direction
0076   const Vector3& sensorDirection() const { return m_wireDir; }
0077   /// @brief To next sensor in the plane
0078   const Vector3& toNextSensor() const { return m_toNext; }
0079   /// @brief To next straw layer
0080   const Vector3& planeNormal() const { return m_planeNorm; }
0081   /// @brief Measurement's radius
0082   double driftRadius() const { return m_driftR; }
0083   /// @brief Measurement's covariance
0084   const std::array<double, 3>& covariance() const { return m_cov; }
0085   /// @brief Measurement's drift unceratinty
0086   double driftUncert() const {
0087     return std::sqrt(m_cov[toUnderlying(ResidualIdx::bending)]);
0088   }
0089   /// @brief Time of record
0090   double time() const { return m_drifT; }
0091   /// @brief All measurements are straws
0092   bool isStraw() const { return m_isStraw; }
0093   /// @brief Dummy return not used in test
0094   bool hasTime() const { return false; }
0095   /// @brief Dummy return not used in test
0096   bool measuresLoc0() const { return false; }
0097   /// @brief Dummy return not used in test
0098   bool measuresLoc1() const { return true; }
0099   void setRadius(const double r, const double uncertR) {
0100     m_driftR = abs(r);
0101     m_cov[toUnderlying(ResidualIdx::bending)] = square(uncertR);
0102   }
0103   void setTimeRecord(const double t) { m_drifT = t; }
0104 
0105  private:
0106   Vector3 m_pos{Vector3::Zero()};
0107   Vector3 m_wireDir{Vector3::UnitX()};
0108   Vector3 m_toNext{Vector3::UnitY()};
0109   Vector3 m_planeNorm{Vector3::UnitZ()};
0110   double m_driftR{0.};
0111   std::array<double, 3> m_cov{filledArray<double, 3>(0.)};
0112   double m_drifT{0.};
0113   bool m_isStraw{true};
0114 };
0115 static_assert(CompositeSpacePoint<StrawTestPoint>);
0116 
0117 class StrawTestCalibrator {
0118  public:
0119   /// @brief Choose the coefficient to arrive at a drift time of 750 ns
0120   ///        for 15 mm
0121   inline static const double CoeffRtoT = 750._ns / std::pow(15._mm, 1. / 3.);
0122   inline static const double CoeffTtoR = 1 / std::pow(CoeffRtoT, 3);
0123 
0124   static double calcDriftUncert(const double driftR) {
0125     return 0.1_mm + 0.15_mm * std::pow(1._mm + std::abs(driftR), -2.);
0126   }
0127   static double driftTime(const double r) {
0128     return CoeffRtoT * std::pow(r, 1. / 3.);
0129   }
0130   static double driftRadius(const double t) {
0131     return CoeffTtoR * std::pow(t, 3);
0132   }
0133 
0134   static double driftRadius(const CalibrationContext& /*ctx*/,
0135                             const StrawTestPoint& straw, const double t0) {
0136     return driftRadius(straw.time() - t0);
0137   }
0138   static double driftVelocity(const CalibrationContext& /*ctx*/,
0139                               const StrawTestPoint& straw, const double t0) {
0140     return 3 * CoeffTtoR * std::pow(straw.time() - t0, 2);
0141   }
0142   static double driftAcceleration(const Acts::CalibrationContext& /*ctx*/,
0143                                   const StrawTestPoint& straw,
0144                                   const double t0) {
0145     return 6 * CoeffTtoR * (straw.time() - t0);
0146   }
0147 };
0148 static_assert(
0149     CompositeSpacePointFastCalibrator<StrawTestCalibrator, StrawTestPoint>);
0150 
0151 /// @brief Generate a random straight track with a flat distribution in theta & y0
0152 ///        Range in theta is [0, 180] degrees in steps of 0.1 excluding the
0153 ///        vicinity around 90 degrees.
0154 ///          In y0, the generated range is [-500, 500] mm
0155 Line_t generateLine(RandomEngine& engine) {
0156   using ParIndex = Line_t::ParIndex;
0157   Line_t::ParamVector linePars{};
0158   linePars[toUnderlying(ParIndex::x0)] = 0.;
0159   linePars[toUnderlying(ParIndex::phi)] = 90._degree;
0160   linePars[toUnderlying(ParIndex::y0)] = uniform_t{-5000., 5000.}(engine);
0161   linePars[toUnderlying(ParIndex::theta)] =
0162       uniform_t{0.1_degree, 179.9_degree}(engine);
0163   if (abs(linePars[toUnderlying(ParIndex::theta)] - 90._degree) < 0.2_degree) {
0164     return generateLine(engine);
0165   }
0166   Line_t line{linePars};
0167   ACTS_DEBUG("Generated parameters theta: "
0168              << (linePars[toUnderlying(ParIndex::theta)] / 1._degree)
0169              << ", y0: " << linePars[toUnderlying(ParIndex::y0)] << " - "
0170              << toString(line.position()) << " + "
0171              << toString(line.direction()));
0172   return line;
0173 }
0174 /// @brief Extrapolate the straight line track through the straw layers to
0175 ///        evaluate which tubes were crossed by the track. The straw layers are
0176 ///        staggered in the z-direction and each layer expands in y. To
0177 ///        estimate, which tubes were crossed, the track is extrapolated to the
0178 ///        z-plane below & above the straw wires. Optionally, the true
0179 ///        drift-radius can be smeared assuming a Gaussian with a drift-radius
0180 ///        dependent uncertainty.
0181 /// @param trajLine: The track to extrapolate
0182 /// @param engine: Random number generator to smear the drift radius
0183 /// @param smearRadius: If true, the drift radius is smeared with a Gaussian
0184 TestStrawCont_t generateStrawCircles(const Line_t& trajLine,
0185                                      RandomEngine& engine, bool smearRadius) {
0186   const Vector3 posStaggering{0., std::cos(60._degree), std::sin(60._degree)};
0187   const Vector3 negStaggering{0., -std::cos(60._degree), std::sin(60._degree)};
0188   /// Number of tube layers per multilayer
0189   constexpr std::size_t nLayersPerMl = 8;
0190   /// Number of overall tubelayers
0191   constexpr std::size_t nTubeLayers = nLayersPerMl * 2;
0192   /// Radius of each straw
0193   constexpr double tubeRadius = 15._mm;
0194   /// Distance between the first <nLayersPerMl> layers and the second pack
0195   constexpr double tubeLayerDist = 1.2_m;
0196 
0197   std::array<Vector3, nTubeLayers> tubePositions{
0198       filledArray<Vector3, nTubeLayers>(Vector3{0., tubeRadius, tubeRadius})};
0199   /// Fill the positions of the reference tubes 1
0200   for (std::size_t l = 1; l < nTubeLayers; ++l) {
0201     const Vector3& layStag{l % 2 == 1 ? posStaggering : negStaggering};
0202     tubePositions[l] = tubePositions[l - 1] + 2. * tubeRadius * layStag;
0203 
0204     if (l == nLayersPerMl) {
0205       tubePositions[l] += tubeLayerDist * Vector3::UnitZ();
0206     }
0207   }
0208   /// Print the staggering
0209   ACTS_DEBUG("##############################################");
0210 
0211   for (std::size_t l = 0; l < nTubeLayers; ++l) {
0212     ACTS_DEBUG("  *** " << (l + 1) << " - " << toString(tubePositions[l]));
0213   }
0214   ACTS_DEBUG("##############################################");
0215 
0216   TestStrawCont_t circles{};
0217   /// Extrapolate the track to the z-planes of the tubes and determine which
0218   /// tubes were actually hit
0219   for (const auto& stag : tubePositions) {
0220     auto planeExtpLow =
0221         PlanarHelper::intersectPlane(trajLine.position(), trajLine.direction(),
0222                                      Vector3::UnitZ(), stag.z() - tubeRadius);
0223     auto planeExtpHigh =
0224         PlanarHelper::intersectPlane(trajLine.position(), trajLine.direction(),
0225                                      Vector3::UnitZ(), stag.z() + tubeRadius);
0226 
0227     ACTS_DEBUG("Extrapolated to plane " << toString(planeExtpLow.position())
0228                                         << " "
0229                                         << toString(planeExtpHigh.position()));
0230 
0231     const auto dToFirstLow = static_cast<int>(std::ceil(
0232         (planeExtpLow.position().y() - stag.y()) / (2. * tubeRadius)));
0233     const auto dToFirstHigh = static_cast<int>(std::ceil(
0234         (planeExtpHigh.position().y() - stag.y()) / (2. * tubeRadius)));
0235 
0236     ACTS_DEBUG("Extrapolated to plane "
0237                << toString(planeExtpLow.position()) << " "
0238                << toString(planeExtpHigh.position())
0239                << " Hit tubes: " << dToFirstLow << " " << dToFirstHigh);
0240 
0241     /// Does the track go from left to right or right to left?
0242     const int dT = dToFirstHigh > dToFirstLow ? 1 : -1;
0243     /// Loop over the candidate tubes and check each one whether the track
0244     /// actually crossed them. Then generate the circle and optionally smear the
0245     /// radius
0246     for (int tN = dToFirstLow - dT; tN != dToFirstHigh + 2 * dT; tN += dT) {
0247       const Vector3 tube = stag + 2. * tN * tubeRadius * Vector3::UnitY();
0248       const double rad = Acts::detail::LineHelper::signedDistance(
0249           tube, Vector3::UnitX(), trajLine.position(), trajLine.direction());
0250 
0251       if (std::abs(rad) > tubeRadius) {
0252         continue;
0253       }
0254       gauss_t dist{rad, StrawTestCalibrator::calcDriftUncert(rad)};
0255       const double smearedR = smearRadius ? std::abs(dist(engine)) : rad;
0256       if (smearedR > tubeRadius) {
0257         continue;
0258       }
0259       circles.emplace_back(std::make_unique<StrawTestPoint>(
0260           tube, smearedR, StrawTestCalibrator::calcDriftUncert(smearedR)));
0261       ACTS_DEBUG("Tube position: "
0262                  << toString(tube) << ", signedRadius: " << rad
0263                  << ", smearedRadius: " << smearedR << ", uncer: "
0264                  << StrawTestCalibrator::calcDriftUncert(smearedR));
0265     }
0266   }
0267   ACTS_DEBUG("Track hit in total " << circles.size() << " tubes ");
0268   return circles;
0269 }
0270 
0271 TestStrawCont_t generateStrips(const Line_t& trajLine, RandomEngine& engine) {
0272   TestStrawCont_t strips{};
0273   constexpr std::array<double, 16> stripZ{-500._mm, -545_mm, -540_mm, -535._mm,
0274                                           -250._mm, -245_mm, -240_mm, -235._mm,
0275                                           335._mm,  340_mm,  345_mm,  350._mm,
0276                                           435._mm,  440_mm,  445_mm,  400._mm};
0277   constexpr double stripPitch = 1._cm;
0278   for (const auto z : stripZ) {
0279     auto planeExtp = PlanarHelper::intersectPlane(
0280         trajLine.position(), trajLine.direction(), Vector3::UnitZ(), z);
0281     Vector3 stripPos = planeExtp.position();
0282     stripPos[eY] = gauss_t{stripPos[eY], stripPitch}(engine);
0283     strips.emplace_back(std::make_unique<StrawTestPoint>(stripPos, stripPitch));
0284   }
0285   return strips;
0286 }
0287 /// @brief Calculate the overall chi2 of the measurements to the track
0288 /// @param measurements: List of candidate straw measurements
0289 /// @param track: Straight line track
0290 double calcChi2(const TestStrawCont_t& measurements, const Line_t& track) {
0291   double chi2{0.};
0292   for (const auto& meas : measurements) {
0293     const double dist = Acts::detail::LineHelper::signedDistance(
0294         meas->localPosition(), meas->sensorDirection(), track.position(),
0295         track.direction());
0296     ACTS_DEBUG("Distance straw: " << toString(meas->localPosition())
0297                                   << ",  r: " << meas->driftRadius()
0298                                   << " - to track: " << abs(dist));
0299 
0300     chi2 += square((abs(dist) - meas->driftRadius()) / meas->driftUncert());
0301   }
0302   return chi2;
0303 }
0304 
0305 #define DECLARE_BRANCH(dTYPE, NAME) \
0306   dTYPE NAME{};                     \
0307   outTree->Branch(#NAME, &NAME);
0308 
0309 BOOST_AUTO_TEST_SUITE(SeedingSuite)
0310 
0311 void testSimpleStrawFit(RandomEngine& engine, TFile& outFile) {
0312   auto outTree = std::make_unique<TTree>("StrawFitTree", "FastFitTree");
0313 
0314   DECLARE_BRANCH(double, trueY0);
0315   DECLARE_BRANCH(double, trueTheta);
0316   DECLARE_BRANCH(double, recoY0);
0317   DECLARE_BRANCH(double, recoTheta);
0318   DECLARE_BRANCH(double, uncertY0);
0319   DECLARE_BRANCH(double, uncertTheta);
0320   DECLARE_BRANCH(double, chi2);
0321   DECLARE_BRANCH(std::size_t, nDoF);
0322   DECLARE_BRANCH(std::size_t, nIter);
0323 
0324   FastStrawLineFitter::Config cfg{};
0325   FastStrawLineFitter fastFitter{cfg, getDefaultLogger("FitterNoT0", logLvl)};
0326   ACTS_INFO("Start simple straw fit test.");
0327   for (std::size_t n = 0; n < nTrials; ++n) {
0328     if ((n + 1) % 1000 == 0) {
0329       ACTS_INFO(" -- processed " << (n + 1) << "/" << nTrials << " events");
0330     }
0331     auto track = generateLine(engine);
0332     auto strawPoints = generateStrawCircles(track, engine, true);
0333     if (strawPoints.size() < 3) {
0334       ACTS_WARNING(__func__ << "() - " << __LINE__ << ": -- event: " << n
0335                             << ", track " << toString(track.position()) << " + "
0336                             << toString(track.direction())
0337                             << " did not lead to any valid measurement ");
0338       continue;
0339     }
0340     const std::vector<int> trueDriftSigns =
0341         CompSpacePointAuxiliaries::strawSigns(track, strawPoints);
0342 
0343     BOOST_CHECK_LE(calcChi2(generateStrawCircles(track, engine, false), track),
0344                    1.e-12);
0345     ACTS_DEBUG("True drift signs: " << trueDriftSigns << ", chi2: " << chi2);
0346 
0347     auto fitResult = fastFitter.fit(strawPoints, trueDriftSigns);
0348     if (!fitResult) {
0349       continue;
0350     }
0351     auto trackPars = track.parameters();
0352 
0353     trueY0 = trackPars[toUnderlying(Line_t::ParIndex::y0)];
0354     trueTheta = trackPars[toUnderlying(Line_t::ParIndex::theta)];
0355     /// Calculate the chi2 again
0356     trackPars[toUnderlying(Line_t::ParIndex::theta)] = (*fitResult).theta;
0357     trackPars[toUnderlying(Line_t::ParIndex::y0)] = (*fitResult).y0;
0358     trackPars[toUnderlying(Line_t::ParIndex::phi)] = 90._degree;
0359     track.updateParameters(trackPars);
0360     ACTS_DEBUG("Updated parameters: "
0361                << (trackPars[toUnderlying(Line_t::ParIndex::theta)] / 1._degree)
0362                << ", y0: " << trackPars[toUnderlying(Line_t::ParIndex::y0)]
0363                << " -- " << toString(track.position()) << " + "
0364                << toString(track.direction()));
0365 
0366     const double testChi2 = calcChi2(strawPoints, track);
0367     ACTS_DEBUG("testChi2: " << testChi2 << ", fit:" << (*fitResult).chi2);
0368 
0369     BOOST_CHECK_LE(abs(testChi2 - (*fitResult).chi2), 1.e-9);
0370     recoTheta = (*fitResult).theta;
0371     recoY0 = (*fitResult).y0;
0372     uncertTheta = (*fitResult).dTheta;
0373     uncertY0 = (*fitResult).dY0;
0374     nDoF = (*fitResult).nDoF;
0375     chi2 = (*fitResult).chi2;
0376     nIter = (*fitResult).nIter;
0377     outTree->Fill();
0378   }
0379   outFile.WriteObject(outTree.get(), outTree->GetName());
0380 }
0381 
0382 void testFitWithT0(RandomEngine& engine, TFile& outFile) {
0383   auto outTree = std::make_unique<TTree>("StrawFitTreeT0", "FastFitTree");
0384 
0385   DECLARE_BRANCH(double, trueY0);
0386   DECLARE_BRANCH(double, trueTheta);
0387   DECLARE_BRANCH(double, trueT0);
0388   DECLARE_BRANCH(double, recoY0);
0389   DECLARE_BRANCH(double, recoTheta);
0390   DECLARE_BRANCH(double, recoT0);
0391   DECLARE_BRANCH(double, uncertY0);
0392   DECLARE_BRANCH(double, uncertTheta);
0393   DECLARE_BRANCH(double, uncertT0);
0394   DECLARE_BRANCH(double, chi2);
0395   DECLARE_BRANCH(double, meanSign);
0396   DECLARE_BRANCH(std::size_t, nDoF);
0397   DECLARE_BRANCH(std::size_t, nIter);
0398 
0399   FastStrawLineFitter::Config cfg{};
0400   cfg.maxIter = 500;
0401   FastStrawLineFitter fastFitter{cfg, getDefaultLogger("FitterWithT0", logLvl)};
0402   StrawTestCalibrator calibrator{};
0403   CalibrationContext cctx{};
0404   ACTS_INFO("Start straw fit with t0 test.");
0405   for (std::size_t n = 0; n < nTrials; ++n) {
0406     if ((n + 1) % 1000 == 0) {
0407       ACTS_INFO(" -- processed " << (n + 1) << "/" << nTrials << " events");
0408     }
0409     auto track = generateLine(engine);
0410     const double timeOffSet = uniform_t{-50._ns, 50._ns}(engine);
0411 
0412     ACTS_DEBUG("Generated time offset: " << inNanoS(timeOffSet) << " [ns]");
0413     auto strawPoints = generateStrawCircles(track, engine, true);
0414 
0415     if (strawPoints.size() < 4) {
0416       ACTS_WARNING(__func__ << "() - " << __LINE__ << ": -- event: " << n
0417                             << ", track " << toString(track.position()) << " + "
0418                             << toString(track.direction())
0419                             << " did not lead to any valid measurement ");
0420       continue;
0421     }
0422     BOOST_CHECK_LE(calcChi2(generateStrawCircles(track, engine, false), track),
0423                    1.e-12);
0424 
0425     /// Fold-in the general offset
0426     const std::vector<int> trueDriftSigns =
0427         CompSpacePointAuxiliaries::strawSigns(track, strawPoints);
0428 
0429     ACTS_DEBUG("Straw signs: " << trueDriftSigns);
0430 
0431     for (auto& meas : strawPoints) {
0432       const double dTime = StrawTestCalibrator::driftTime(meas->driftRadius());
0433       BOOST_CHECK_CLOSE(StrawTestCalibrator::driftRadius(dTime),
0434                         meas->driftRadius(), 1.e-12);
0435       meas->setTimeRecord(dTime + timeOffSet);
0436 
0437       const double updatedR =
0438           StrawTestCalibrator::driftRadius(dTime + timeOffSet);
0439 
0440       BOOST_CHECK_CLOSE(StrawTestCalibrator::driftRadius(dTime),
0441                         calibrator.driftRadius(cctx, *meas, timeOffSet), 1.e-3);
0442 
0443       ACTS_DEBUG("Update drift radius of tube "
0444                  << toString(meas->localPosition()) << " from "
0445                  << meas->driftRadius() << " to " << updatedR
0446                  << ", dTime: " << inNanoS(dTime));
0447       meas->setRadius(updatedR, StrawTestCalibrator::calcDriftUncert(updatedR));
0448 
0449       /// Calculate the numerical derivatives
0450       constexpr double h = 1.e-8_ns;
0451       const double numV =
0452           -(calibrator.driftRadius(cctx, *meas, timeOffSet + h) -
0453             calibrator.driftRadius(cctx, *meas, timeOffSet - h)) /
0454           (2. * h);
0455       BOOST_CHECK_LE(
0456           abs(numV - calibrator.driftVelocity(cctx, *meas, timeOffSet)) /
0457               std::max(numV, 1.),
0458           1.e-3);
0459     }
0460 
0461     auto result = fastFitter.fit(cctx, calibrator, strawPoints, trueDriftSigns,
0462                                  std::nullopt);
0463 
0464     if (!result) {
0465       continue;
0466     }
0467     auto linePars = track.parameters();
0468     /// True parameters
0469     trueY0 = linePars[toUnderlying(Line_t::ParIndex::y0)];
0470     trueTheta = linePars[toUnderlying(Line_t::ParIndex::theta)];
0471     trueT0 = inNanoS(timeOffSet);
0472     /// Fit parameters
0473     recoY0 = (*result).y0;
0474     recoTheta = (*result).theta;
0475     recoT0 = inNanoS((*result).t0);
0476     uncertY0 = (*result).dY0;
0477     uncertTheta = (*result).dTheta;
0478     uncertT0 = inNanoS((*result).dT0);
0479     nDoF = (*result).nDoF;
0480     chi2 = (*result).chi2;
0481     nIter = (*result).nIter;
0482     meanSign = {0.};
0483     std::ranges::for_each(trueDriftSigns,
0484                           [&meanSign](const int sign) { meanSign += sign; });
0485     meanSign /= static_cast<double>(trueDriftSigns.size());
0486     outTree->Fill();
0487   }
0488   outFile.WriteObject(outTree.get(), outTree->GetName());
0489 }
0490 
0491 void testStripFit(RandomEngine& engine, TFile& outFile) {
0492   auto outTree = std::make_unique<TTree>("StripFitTree", "FastFitTree");
0493 
0494   DECLARE_BRANCH(double, trueY0);
0495   DECLARE_BRANCH(double, trueTheta);
0496   DECLARE_BRANCH(double, recoY0);
0497   DECLARE_BRANCH(double, recoTheta);
0498   DECLARE_BRANCH(double, uncertY0);
0499   DECLARE_BRANCH(double, uncertTheta);
0500   DECLARE_BRANCH(double, chi2);
0501   DECLARE_BRANCH(std::size_t, nDoF);
0502   DECLARE_BRANCH(std::size_t, nIter);
0503 
0504   FastStrawLineFitter::Config cfg{};
0505   FastStrawLineFitter fastFitter{cfg, getDefaultLogger("StripFitter", logLvl)};
0506   ACTS_INFO("Start strip fit test.");
0507   for (std::size_t n = 0; n < nTrials; ++n) {
0508     if ((n + 1) % 1000 == 0) {
0509       ACTS_INFO(" -- processed " << (n + 1) << "/" << nTrials << " events");
0510     }
0511     auto track = generateLine(engine);
0512     auto linePars = track.parameters();
0513     /// True parameters
0514     trueY0 = linePars[toUnderlying(Line_t::ParIndex::y0)];
0515     trueTheta = linePars[toUnderlying(Line_t::ParIndex::theta)];
0516 
0517     auto stripPoints = generateStrips(track, engine);
0518 
0519     auto result = fastFitter.fit(stripPoints, ResidualIdx::bending);
0520     if (!result) {
0521       continue;
0522     }
0523     /// Fit parameters
0524     recoY0 = (*result).y0;
0525     recoTheta = (*result).theta;
0526     uncertY0 = (*result).dY0;
0527     uncertTheta = (*result).dTheta;
0528     nDoF = (*result).nDoF;
0529     nIter = (*result).nIter;
0530     chi2 = 0.;
0531     auto pos = recoY0 * Vector3::UnitY();
0532     auto dir = makeDirectionFromPhiTheta(90._degree, recoTheta);
0533     for (const auto& strip : stripPoints) {
0534       chi2 += CompSpacePointAuxiliaries::chi2Term(pos, dir, *strip);
0535     }
0536 
0537     outTree->Fill();
0538   }
0539 
0540   outFile.WriteObject(outTree.get(), outTree->GetName());
0541 }
0542 
0543 BOOST_AUTO_TEST_CASE(FitterTests) {
0544   RandomEngine engine{1800};
0545 
0546   std::unique_ptr<TFile> outFile{
0547       TFile::Open("FastStrawLineFitTest.root", "RECREATE")};
0548 
0549   BOOST_CHECK_EQUAL(outFile->IsZombie(), false);
0550 
0551   testSimpleStrawFit(engine, *outFile);
0552   testStripFit(engine, *outFile);
0553   testFitWithT0(engine, *outFile);
0554 }
0555 
0556 BOOST_AUTO_TEST_SUITE_END()
0557 }  // namespace ActsTests
[ 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!