src/Vertexing/FsmwMode1dFinder.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/FsmwMode1dFinder.hpp"
0010
0011 #include "Acts/Vertexing/VertexingError.hpp"
0012
0013 #include <algorithm>
0014 #include <cmath>
0015 #include <limits>
0016
0017 Acts::FsmwMode1dFinder::FsmwMode1dFinder(double firstFraction, double fraction)
0018     : m_firstFraction(firstFraction), m_fraction(fraction) {}
0019
0020 Acts::Result<double> Acts::FsmwMode1dFinder::getMode(
0021     std::vector<std::pair<double, double>> inputVector) const {
0022   if (inputVector.empty()) {
0023     return VertexingError::EmptyInput;
0024   }
0025   if (inputVector.size() == 1) {
0026     return inputVector.begin()->first;
0027   }
0028
0029   // first of all order the vector according to the double value
0030
0031   std::ranges::sort(inputVector, {}, [](const auto& i) { return i.first; });
0032
0033   // begin to consider a certain number of elements according to the fraction
0034   auto begin = inputVector.begin();
0035   auto end = inputVector.end();
0036
0037   double overallweight(0.);
0038   auto best_begin = begin;
0039   auto best_end = end;
0040
0041   double last_value = std::numeric_limits<double>::max();
0042
0043   bool isthelast = false;
0044
0045   int counter = 0;
0046   double fraction = m_firstFraction;
0047   while (!isthelast) {
0048     counter += 1;
0049     if (counter == 2) {
0050       fraction = m_fraction;
0051     }
0052     int step = static_cast<int>(std::floor(fraction * (end - begin + 1)));
0053     overallweight = 0.;
0054     {
0055       auto i = begin;
0056       if (step > 0) {
0057         auto j_end = i + step - 1;
0058         for (auto j = i; j != j_end; j++) {
0059           overallweight += j->second;
0060         }
0061       }
0062     }
0063     auto i_last = begin + step - 1;
0064
0065     for (auto i = begin; i != (end - step + 1); ++i, ++i_last) {
0066       // calculate the weight the interval should be divided into
0067       overallweight += i_last->second;
0068
0069       double new_value = ((i + step - 1)->first - i->first) / overallweight;
0070       if (new_value < last_value) {
0071         last_value = ((i + step - 1)->first - i->first) / overallweight;
0072         best_begin = i;
0073         best_end = i + step - 1;
0074       }
0075       overallweight -= i->second;
0076     }
0077
0078     // assign the new begin and end...
0079     begin = best_begin;
0080     end = best_end;
0081     last_value = std::numeric_limits<double>::max();
0082
0083     // Now it should have returned the value with smallest (x2-x1)/weight
0084     if (best_end - best_begin <= 2) {
0085       isthelast = true;
0086     }
0087   }
0088
0089   if (best_end - best_begin == 2) {
0090     auto medium = begin;
0091     medium++;
0092     return (begin->first * begin->second + medium->first * medium->second +
0093             end->first * end->second) /
0094            (begin->second + medium->second + end->second);
0095   }
0096
0097   return (begin->first * begin->second + end->first * end->second) /
0098          (begin->second + end->second);
0099 }