src/Vertexing/FsmwMode1dFinder.cpp

0001 // This file is part of the Acts project.
0002 //
0003 // Copyright (C) 2019 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 http://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::sort(inputVector.begin(), inputVector.end(),
0032             [](std::pair<double, double> a, std::pair<double, double> b) {
0033               return a.first < b.first;
0034             });
0035
0036   // begin to consider a certain number of elements according to the fraction
0037   auto begin = inputVector.begin();
0038   auto end = inputVector.end();
0039
0040   double overallweight(0.);
0041   auto best_begin = begin;
0042   auto best_end = end;
0043
0044   double last_value = std::numeric_limits<double>::max();
0045
0046   bool isthelast = false;
0047
0048   int counter = 0;
0049   double fraction = m_firstFraction;
0050   while (!isthelast) {
0051     counter += 1;
0052     if (counter == 2) {
0053       fraction = m_fraction;
0054     }
0055     int step = (int)std::floor(fraction * (end - begin + 1));
0056     overallweight = 0.;
0057     {
0058       auto i = begin;
0059       if (step > 0) {
0060         auto j_end = i + step - 1;
0061         for (auto j = i; j != j_end; j++) {
0062           overallweight += j->second;
0063         }
0064       }
0065     }
0066     auto i_last = begin + step - 1;
0067
0068     for (auto i = begin; i != (end - step + 1); ++i, ++i_last) {
0069       // calculate the weight the interval should be divided into
0070       overallweight += i_last->second;
0071
0072       double new_value = ((i + step - 1)->first - i->first) / overallweight;
0073       if (new_value < last_value) {
0074         last_value = ((i + step - 1)->first - i->first) / overallweight;
0075         best_begin = i;
0076         best_end = i + step - 1;
0077       }
0078       overallweight -= i->second;
0079     }
0080
0081     // assign the new begin and end...
0082     begin = best_begin;
0083     end = best_end;
0084     last_value = std::numeric_limits<double>::max();
0085
0086     // Now it should have returned the value with smallest (x2-x1)/weight
0087     if (best_end - best_begin <= 2) {
0088       isthelast = true;
0089     }
0090   }
0091
0092   if (best_end - best_begin == 2) {
0093     auto medium = begin;
0094     medium++;
0095     return (begin->first * begin->second + medium->first * medium->second +
0096             end->first * end->second) /
0097            (begin->second + medium->second + end->second);
0098   }
0099
0100   return (begin->first * begin->second + end->first * end->second) /
0101          (begin->second + end->second);
0102 }