Core/Seeding/StrawLineFitterTest.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 <boost/test/unit_test.hpp>
0010
0011 #include <chrono>
0012 #include <format>
0013 #include <future>
0014 #include <iostream>
0015 #include <random>
0016 #include <ranges>
0017 #include <span>
0018 #include <thread>
0019
0020 #include <TFile.h>
0021 #include <TH1D.h>
0022 #include <TTree.h>
0023
0024 #include "StrawHitGeneratorHelper.hpp"
0025
0026 using TimePoint_t = std::chrono::system_clock::time_point;
0027 using Fitter_t = CompositeSpacePointLineFitter;
0028
0029 constexpr auto logLvl = Acts::Logging::Level::INFO;
0030 constexpr std::size_t nEvents = 1;
0031 constexpr long int nThreads = 1;
0032 std::mutex writeMutex{};
0033
0034 ACTS_LOCAL_LOGGER(getDefaultLogger("StrawLineFitterTest", logLvl));
0035
0036 namespace ActsTests {
0037
0038 using GenCfg_t = MeasurementGenerator::Config;
0039
0040 #define DECLARE_BRANCH(dType, bName) \
0041   dType bName{};                     \
0042   outTree->Branch(#bName, &bName);
0043 // NOLINTBEGIN
0044 long int runFitTest(Fitter_t::Config fitCfg, GenCfg_t genCfg,
0045                     const std::string& testName, const unsigned seed,
0046                     TFile& outFile) {
0047   // NOLINTEND
0048   auto outTree = std::make_unique<TTree>(
0049       std::format("{:}Tree", testName).c_str(), "MonitorTree");
0050   outTree->SetDirectory(nullptr);
0051
0052   DECLARE_BRANCH(double, trueY0);
0053   DECLARE_BRANCH(double, trueX0);
0054   DECLARE_BRANCH(double, trueTheta);
0055   DECLARE_BRANCH(double, truePhi);
0056   DECLARE_BRANCH(double, trueT0);
0057   DECLARE_BRANCH(double, trueProjTheta);
0058   DECLARE_BRANCH(double, trueProjPhi);
0059
0060   DECLARE_BRANCH(double, recoY0);
0061   DECLARE_BRANCH(double, recoX0);
0062   DECLARE_BRANCH(double, recoTheta);
0063   DECLARE_BRANCH(double, recoPhi);
0064   DECLARE_BRANCH(double, recoT0);
0065   DECLARE_BRANCH(double, recoProjTheta);
0066   DECLARE_BRANCH(double, recoProjPhi);
0067
0068   DECLARE_BRANCH(double, sigmaY0);
0069   DECLARE_BRANCH(double, sigmaX0);
0070   DECLARE_BRANCH(double, sigmaTheta);
0071   DECLARE_BRANCH(double, sigmaPhi);
0072   DECLARE_BRANCH(double, sigmaT0);
0073
0074   DECLARE_BRANCH(double, chi2);
0075   DECLARE_BRANCH(unsigned, nIter);
0076   DECLARE_BRANCH(unsigned, nDoF);
0077   DECLARE_BRANCH(char, converged);
0078
0079   RandomEngine engine{seed};
0080   Fitter_t fitter{fitCfg, getDefaultLogger(
0081                               std::format("LineFitter_{:}", testName), logLvl)};
0082
0083   TimePoint_t start = std::chrono::system_clock::now();
0084
0085   ACTS_INFO("Start test " << testName << ".");
0086
0087   /// @brief Fill the parameter array to the tree variables
0088   /// @param pars: Parameter array to safe
0089   /// @param y0: Reference to the variable storing y0
0090   /// @param x0: Reference to the variable storing x0
0091   /// @param theta: Reference to the variable storing theta
0092   /// @param phi: Reference to the variable storing phi
0093   auto fillPars = [](const auto pars, double& y0, double& x0, double& theta,
0094                      double& phi) {
0095     y0 = pars[toUnderlying(FitParIndex::y0)];
0096     x0 = pars[toUnderlying(FitParIndex::x0)];
0097     Vector3 properDir =
0098         makeDirectionFromPhiTheta(pars[toUnderlying(FitParIndex::phi)],
0099                                   pars[toUnderlying(FitParIndex::theta)]);
0100     theta =
0101         VectorHelpers::theta(copySign(properDir, properDir.z())) / 1._degree;
0102     phi = VectorHelpers::phi(copySign(properDir, properDir.z())) / 1._degree;
0103   };
0104   /// @brief Fill the
0105   auto fillProjected = [](const auto pars, double& projTheta, double& projPhi) {
0106     auto dir =
0107         makeDirectionFromPhiTheta(pars[toUnderlying(FitParIndex::phi)],
0108                                   pars[toUnderlying(FitParIndex::theta)]);
0109     projTheta = std::atan(dir[ePos1] / dir[ePos2]) / 1._degree;
0110     projPhi = std::atan(dir[ePos0] / dir[ePos2]) / 1._degree;
0111   };
0112   // Pass a localToGlobal transform to the calibrator to proper handling the ToF
0113   auto calibrator = std::make_unique<SpCalibrator>();
0114   std::size_t goodFits{0};
0115   for (std::size_t evt = 0; evt < nEvents; ++evt) {
0116     const auto line = generateLine(engine, logger());
0117     fillPars(line.parameters(), trueY0, trueX0, trueTheta, truePhi);
0118     fillProjected(line.parameters(), trueProjTheta, trueProjPhi);
0119     const double t0 = uniform{-50._ns, 50._ns}(engine);
0120     trueT0 = t0 / 1._ns;
0121
0122     using FitOpts_t = Fitter_t::FitOptions<Container_t, SpCalibrator>;
0123
0124     FitOpts_t fitOpts{};
0125     fitOpts.calibrator = calibrator.get();
0126
0127     fitOpts.selector.connect<&isGoodHit>();
0128     fitOpts.measurements =
0129         MeasurementGenerator::spawn(line, t0, engine, genCfg, logger());
0130     fitOpts.startParameters = startParameters(line, fitOpts.measurements);
0131     fillPars(fitOpts.startParameters, recoY0, recoX0, recoTheta, recoPhi);
0132
0133     //
0134
0135     auto result = fitter.fit(std::move(fitOpts));
0136     if (!result.converged) {
0137       ACTS_DEBUG("Fit " << outTree->GetName() << " failed.");
0138       converged = 0;
0139       chi2 = -1.;
0140       nDoF = 1;
0141       nIter = fitter.config().maxIter;
0142       outTree->Fill();
0143       continue;
0144     }
0145
0146     ACTS_DEBUG("Fit Successful.");
0147     converged = 1;
0148     ++goodFits;
0149     fillPars(result.parameters, recoY0, recoX0, recoTheta, recoPhi);
0150     fillProjected(result.parameters, recoProjTheta, recoProjPhi);
0151
0152     recoT0 = result.parameters[toUnderlying(FitParIndex::t0)] / 1._ns;
0153
0154     auto extractUncert = [&result](const auto idx) {
0155       return std::sqrt(result.covariance(toUnderlying(idx), toUnderlying(idx)));
0156     };
0157     sigmaY0 = extractUncert(FitParIndex::y0);
0158     sigmaX0 = extractUncert(FitParIndex::x0);
0159     sigmaTheta = extractUncert(FitParIndex::theta) / 1._degree;
0160     sigmaPhi = extractUncert(FitParIndex::phi) / 1._degree;
0161     sigmaT0 = extractUncert(FitParIndex::t0) / 1._ns;
0162
0163     chi2 = result.chi2;
0164     nDoF = result.nDoF;
0165     nIter = result.nIter;
0166
0167     outTree->Fill();
0168     if ((evt + 1) % 10000 == 0u) {
0169       ACTS_INFO("Processed " << (evt + 1) << "/" << nEvents
0170                              << " events. Test: " << outTree->GetName());
0171     }
0172   }
0173
0174   TimePoint_t end = std::chrono::system_clock::now();  // timing: get end time
0175   auto diff = std::chrono::duration_cast<std::chrono::milliseconds>(end - start)
0176                   .count();
0177   std::unique_lock guard{writeMutex};
0178   outFile.WriteObject(outTree.get(), outTree->GetName());
0179   ACTS_INFO("Test " << outTree->GetName() << " finished. " << goodFits
0180                     << " tracks written. It took " << (diff / 1000)
0181                     << " seconds.");
0182   return diff;
0183 }
0184 #undef DECLARE_BRANCH
0185
0186 BOOST_AUTO_TEST_SUITE(StrawLineFitTestSuite)
0187
0188 BOOST_AUTO_TEST_CASE(SimpleLineFit) {
0189   using namespace std::chrono_literals;
0190
0191   auto outFile =
0192       std::make_unique<TFile>("StrawLineFitterTest.root", "RECREATE");
0193
0194   Fitter_t::Config fitCfg{};
0195   fitCfg.useHessian = false;
0196   fitCfg.calcAlongStraw = true;
0197   fitCfg.recalibrate = false;
0198   fitCfg.useFastFitter = false;
0199   fitCfg.ranges[toUnderlying(FitParIndex::theta)] =
0200       std::array{1._degree, 179._degree};
0201   fitCfg.ranges[toUnderlying(FitParIndex::phi)] =
0202       std::array{-179._degree, 179._degree};
0203   fitCfg.ranges[toUnderlying(FitParIndex::x0)] = std::array{-1000., 1000.};
0204   fitCfg.ranges[toUnderlying(FitParIndex::y0)] = std::array{-1000., 1000.};
0205   /// Configuration for fast pre-fit
0206   Fitter_t::Config fastPreCfg{fitCfg};
0207   fastPreCfg.useFastFitter = true;
0208   /// Configuration for fast only fit
0209   Fitter_t::Config fastCfg{fastPreCfg};
0210   fastCfg.fastPreFitter = false;
0211
0212   // 2D straw only test
0213   std::vector<std::pair<std::string, std::future<long int>>> timings{};
0214   using namespace std::chrono_literals;
0215
0216   auto sendSleep = [&timings]() {
0217     do {
0218       std::this_thread::sleep_for(100ms);
0219     } while (std::ranges::count_if(timings, [](const auto& timeObj) {
0220                return timeObj.second.wait_for(0ms) != std::future_status::ready;
0221              }) >= nThreads);
0222   };
0223   /// @brief Helper lambda to launch the test for a given measurement configuration.
0224   ///        Always three tests are launched:
0225   ///           Fast: Use only the fast fitter to obtain the results
0226   ///           FastPre: Use the fast fitter for a pre estimate of the
0227   ///           parameters
0228   ///                    followed by the full fitter
0229   ///           Full: Just use the full fitter to obtain the results
0230   /// @param testName: Name of the tests
0231   /// @param genCfg: Configuration object to generate the measurements
0232   /// @param seed: Seed number for the random number generator
0233   auto launchTest = [&](const std::string& testName, const GenCfg_t& genCfg,
0234                         const unsigned seed) {
0235     sendSleep();
0236     timings.emplace_back(
0237         "Fast" + testName, std::async(std::launch::async, [&]() {
0238           return runFitTest(fastCfg, genCfg, "Fast" + testName, seed, *outFile);
0239         }));
0240     sendSleep();
0241     timings.emplace_back(testName, std::async(std::launch::async, [&]() {
0242                            return runFitTest(fitCfg, genCfg, testName, seed,
0243                                              *outFile);
0244                          }));
0245
0246     sendSleep();
0247     timings.emplace_back(
0248         "FastPre" + testName, std::async(std::launch::async, [&]() {
0249           return runFitTest(fastPreCfg, genCfg, "FastPre" + testName, seed,
0250                             *outFile);
0251         }));
0252     sendSleep();
0253   };
0254   {
0255     GenCfg_t genCfg{};
0256     genCfg.twinStraw = false;
0257     genCfg.createStrips = false;
0258     launchTest("StrawOnlyTest", genCfg, 1602);
0259   }
0260   // 2D straws + twin measurement test
0261   {
0262     GenCfg_t genCfg{};
0263     genCfg.createStraws = true;
0264     genCfg.twinStraw = true;
0265     genCfg.createStrips = false;
0266
0267     launchTest("StrawAndTwinTest", genCfg, 1503);
0268   }
0269   // 1D straws + single strip measurements
0270   {
0271     GenCfg_t genCfg{};
0272     genCfg.createStrips = true;
0273     genCfg.twinStraw = false;
0274     genCfg.combineSpacePoints = false;
0275     genCfg.discretizeStrips = true;
0276     genCfg.createStraws = true;
0277     genCfg.stripPitchLoc1 = 2._cm;
0278     genCfg.stripPitchLoc0 = 3.5_cm;
0279     launchTest("StrawAndStripTest", genCfg, 1701);
0280   }
0281   // 1D straws + 2D strip measurements
0282   {
0283     GenCfg_t genCfg{};
0284
0285     genCfg.createStrips = true;
0286     genCfg.twinStraw = false;
0287     genCfg.combineSpacePoints = true;
0288     genCfg.discretizeStrips = true;
0289     genCfg.createStraws = true;
0290     genCfg.stripPitchLoc1 = 2._cm;
0291     genCfg.stripPitchLoc0 = 3.5_cm;
0292     launchTest("StrawAndStrip2DTest", genCfg, 1404);
0293   }
0294   // Strip only
0295   {
0296     GenCfg_t genCfg{};
0297     genCfg.createStrips = true;
0298     genCfg.twinStraw = false;
0299     genCfg.combineSpacePoints = false;
0300     genCfg.discretizeStrips = true;
0301     genCfg.createStraws = false;
0302     genCfg.stripPitchLoc1 = 500._um;
0303     genCfg.stripPitchLoc0 = 3._cm;
0304     launchTest("StripOnlyTest", genCfg, 2070);
0305   }
0306   // 2D Strip only
0307   {
0308     GenCfg_t genCfg{};
0309     genCfg.createStrips = true;
0310     genCfg.twinStraw = false;
0311     genCfg.combineSpacePoints = true;
0312     genCfg.discretizeStrips = true;
0313     genCfg.createStraws = false;
0314     genCfg.stripPitchLoc1 = 500._um;
0315     genCfg.stripPitchLoc0 = 3._cm;
0316     launchTest("Strip2DOnlyTest", genCfg, 2225);
0317   }
0318   // Strip stereo test
0319   {
0320     GenCfg_t genCfg{};
0321     genCfg.createStrips = true;
0322     genCfg.twinStraw = false;
0323     genCfg.combineSpacePoints = false;
0324     genCfg.discretizeStrips = true;
0325     genCfg.createStraws = false;
0326     genCfg.stripPitchLoc1 = 500._um;
0327     genCfg.stripDirLoc0.clear();
0328     genCfg.stripDirLoc1 = {
0329         makeDirectionFromPhiTheta(0._degree, 90._degree),
0330         makeDirectionFromPhiTheta(0._degree, 90._degree),
0331         makeDirectionFromPhiTheta(-1.5_degree, 90._degree),
0332         makeDirectionFromPhiTheta(1.5_degree, 90._degree),
0333         makeDirectionFromPhiTheta(-1.5_degree, 90._degree),
0334         makeDirectionFromPhiTheta(1.5_degree, 90._degree),
0335         makeDirectionFromPhiTheta(0._degree, 90._degree),
0336         makeDirectionFromPhiTheta(0._degree, 90._degree),
0337
0338     };
0339
0340     fitCfg.parsToUse = {FitParIndex::x0, FitParIndex::y0, FitParIndex::theta,
0341                         FitParIndex::phi};
0342     launchTest("StereoStripTest", genCfg, 1800);
0343   }
0344   /// Wait until all tests ar ecompleted
0345   while (std::ranges::any_of(timings, [](const auto& lblTh) {
0346     return lblTh.second.wait_for(0ms) != std::future_status::ready;
0347   })) {
0348     std::this_thread::sleep_for(10ms);
0349   }
0350   {
0351     auto timeHisto = std::make_unique<TH1D>("TestTimings", "timings",
0352                                             timings.size(), 0, timings.size());
0353     int bin{1};
0354     for (auto& [label, result] : timings) {
0355       timeHisto->GetXaxis()->SetBinLabel(bin, label.c_str());
0356       timeHisto->SetBinContent(bin, static_cast<double>(result.get()) /
0357                                         static_cast<double>(nEvents));
0358       ++bin;
0359     }
0360     outFile->WriteObject(timeHisto.get(), timeHisto->GetName());
0361   }
0362 }
0363
0364 BOOST_AUTO_TEST_SUITE_END()
0365 }  // namespace ActsTests