TrackFinder.cxx

Go to the documentation of this file.

// Copyright 2019-2020 CERN and copyright holders of ALICE O2.
// See https://alice-o2.web.cern.ch/copyright for details of the copyright holders.
// All rights not expressly granted are reserved.
//
// This software is distributed under the terms of the GNU General Public
// License v3 (GPL Version 3), copied verbatim in the file "COPYING".
//
// In applying this license CERN does not waive the privileges and immunities
// granted to it by virtue of its status as an Intergovernmental Organization
// or submit itself to any jurisdiction.
 
 
#include "MCHTracking/TrackFinder.h"
 
#include <cassert>
#include <cmath>
#include <iostream>
#include <stdexcept>
 
#include <TGeoGlobalMagField.h>
#include <TMatrixD.h>
#include <TMath.h>
 
#include "Field/MagneticField.h"
#include "MCHBase/Error.h"
#include "MCHBase/TrackerParam.h"
#include "MCHTracking/TrackExtrap.h"
 
namespace o2
{
namespace mch
{
 
using namespace std;
 
constexpr double TrackFinder::SDefaultChamberZ[10];
constexpr double TrackFinder::SChamberThicknessInX0[10];
constexpr int TrackFinder::SNDE[10];
 
//_________________________________________________________________________________________________
void TrackFinder::init()
{
 
  // Set the parameters used for fitting the tracks during the tracking
  const auto& trackerParam = TrackerParam::Instance();
  mTrackFitter.setBendingVertexDispersion(trackerParam.bendingVertexDispersion);
  mTrackFitter.setChamberResolution(trackerParam.chamberResolutionX, trackerParam.chamberResolutionY);
  mTrackFitter.smoothTracks(true);
 
  // use the Runge-Kutta extrapolation v2
  TrackExtrap::useExtrapV2();
 
  // Pre-compute some parameters used during the tracking
  mChamberResolutionX2 = trackerParam.chamberResolutionX * trackerParam.chamberResolutionX;
  mChamberResolutionY2 = trackerParam.chamberResolutionY * trackerParam.chamberResolutionY;
  mBendingVertexDispersion2 = trackerParam.bendingVertexDispersion * trackerParam.bendingVertexDispersion;
  mMaxChi2ForTracking = 2. * trackerParam.sigmaCutForTracking * trackerParam.sigmaCutForTracking;
  mMaxChi2ForImprovement = 2. * trackerParam.sigmaCutForImprovement * trackerParam.sigmaCutForImprovement;
 
  // set the maximum MCS angle in chamber from the minimum acceptable momentum
  TrackParam param{};
  double inverseBendingP = (SMinBendingMomentum > 0.) ? 1. / SMinBendingMomentum : 1.;
  param.setInverseBendingMomentum(inverseBendingP);
  for (int iCh = 0; iCh < 10; ++iCh) {
    mMaxMCSAngle2[iCh] = TrackExtrap::getMCSAngle2(param, SChamberThicknessInX0[iCh], 1.);
  }
 
  // prepare the internal array of list of vector
  // grouping DEs in z-planes (2 for chambers 1-4 and 4 for chambers 5-10)
  for (int iCh = 0; iCh < 4; ++iCh) {
    mClusters[2 * iCh].reserve(2);
    mClusters[2 * iCh].emplace_back(100 * (iCh + 1) + 1, nullptr);
    mClusters[2 * iCh].emplace_back(100 * (iCh + 1) + 3, nullptr);
    mClusters[2 * iCh + 1].reserve(2);
    mClusters[2 * iCh + 1].emplace_back(100 * (iCh + 1), nullptr);
    mClusters[2 * iCh + 1].emplace_back(100 * (iCh + 1) + 2, nullptr);
  }
  for (int iCh = 4; iCh < 6; ++iCh) {
    mClusters[8 + 4 * (iCh - 4)].reserve(5);
    mClusters[8 + 4 * (iCh - 4)].emplace_back(100 * (iCh + 1), nullptr);
    mClusters[8 + 4 * (iCh - 4)].emplace_back(100 * (iCh + 1) + 2, nullptr);
    mClusters[8 + 4 * (iCh - 4)].emplace_back(100 * (iCh + 1) + 4, nullptr);
    mClusters[8 + 4 * (iCh - 4)].emplace_back(100 * (iCh + 1) + 14, nullptr);
    mClusters[8 + 4 * (iCh - 4)].emplace_back(100 * (iCh + 1) + 16, nullptr);
    mClusters[8 + 4 * (iCh - 4) + 1].reserve(4);
    mClusters[8 + 4 * (iCh - 4) + 1].emplace_back(100 * (iCh + 1) + 1, nullptr);
    mClusters[8 + 4 * (iCh - 4) + 1].emplace_back(100 * (iCh + 1) + 3, nullptr);
    mClusters[8 + 4 * (iCh - 4) + 1].emplace_back(100 * (iCh + 1) + 15, nullptr);
    mClusters[8 + 4 * (iCh - 4) + 1].emplace_back(100 * (iCh + 1) + 17, nullptr);
    mClusters[8 + 4 * (iCh - 4) + 2].reserve(4);
    mClusters[8 + 4 * (iCh - 4) + 2].emplace_back(100 * (iCh + 1) + 6, nullptr);
    mClusters[8 + 4 * (iCh - 4) + 2].emplace_back(100 * (iCh + 1) + 8, nullptr);
    mClusters[8 + 4 * (iCh - 4) + 2].emplace_back(100 * (iCh + 1) + 10, nullptr);
    mClusters[8 + 4 * (iCh - 4) + 2].emplace_back(100 * (iCh + 1) + 12, nullptr);
    mClusters[8 + 4 * (iCh - 4) + 3].reserve(5);
    mClusters[8 + 4 * (iCh - 4) + 3].emplace_back(100 * (iCh + 1) + 5, nullptr);
    mClusters[8 + 4 * (iCh - 4) + 3].emplace_back(100 * (iCh + 1) + 7, nullptr);
    mClusters[8 + 4 * (iCh - 4) + 3].emplace_back(100 * (iCh + 1) + 9, nullptr);
    mClusters[8 + 4 * (iCh - 4) + 3].emplace_back(100 * (iCh + 1) + 11, nullptr);
    mClusters[8 + 4 * (iCh - 4) + 3].emplace_back(100 * (iCh + 1) + 13, nullptr);
  }
  for (int iCh = 6; iCh < 10; ++iCh) {
    mClusters[8 + 4 * (iCh - 4)].reserve(7);
    mClusters[8 + 4 * (iCh - 4)].emplace_back(100 * (iCh + 1), nullptr);
    mClusters[8 + 4 * (iCh - 4)].emplace_back(100 * (iCh + 1) + 2, nullptr);
    mClusters[8 + 4 * (iCh - 4)].emplace_back(100 * (iCh + 1) + 4, nullptr);
    mClusters[8 + 4 * (iCh - 4)].emplace_back(100 * (iCh + 1) + 6, nullptr);
    mClusters[8 + 4 * (iCh - 4)].emplace_back(100 * (iCh + 1) + 20, nullptr);
    mClusters[8 + 4 * (iCh - 4)].emplace_back(100 * (iCh + 1) + 22, nullptr);
    mClusters[8 + 4 * (iCh - 4)].emplace_back(100 * (iCh + 1) + 24, nullptr);
    mClusters[8 + 4 * (iCh - 4) + 1].reserve(6);
    mClusters[8 + 4 * (iCh - 4) + 1].emplace_back(100 * (iCh + 1) + 1, nullptr);
    mClusters[8 + 4 * (iCh - 4) + 1].emplace_back(100 * (iCh + 1) + 3, nullptr);
    mClusters[8 + 4 * (iCh - 4) + 1].emplace_back(100 * (iCh + 1) + 5, nullptr);
    mClusters[8 + 4 * (iCh - 4) + 1].emplace_back(100 * (iCh + 1) + 21, nullptr);
    mClusters[8 + 4 * (iCh - 4) + 1].emplace_back(100 * (iCh + 1) + 23, nullptr);
    mClusters[8 + 4 * (iCh - 4) + 1].emplace_back(100 * (iCh + 1) + 25, nullptr);
    mClusters[8 + 4 * (iCh - 4) + 2].reserve(6);
    mClusters[8 + 4 * (iCh - 4) + 2].emplace_back(100 * (iCh + 1) + 8, nullptr);
    mClusters[8 + 4 * (iCh - 4) + 2].emplace_back(100 * (iCh + 1) + 10, nullptr);
    mClusters[8 + 4 * (iCh - 4) + 2].emplace_back(100 * (iCh + 1) + 12, nullptr);
    mClusters[8 + 4 * (iCh - 4) + 2].emplace_back(100 * (iCh + 1) + 14, nullptr);
    mClusters[8 + 4 * (iCh - 4) + 2].emplace_back(100 * (iCh + 1) + 16, nullptr);
    mClusters[8 + 4 * (iCh - 4) + 2].emplace_back(100 * (iCh + 1) + 18, nullptr);
    mClusters[8 + 4 * (iCh - 4) + 3].reserve(7);
    mClusters[8 + 4 * (iCh - 4) + 3].emplace_back(100 * (iCh + 1) + 7, nullptr);
    mClusters[8 + 4 * (iCh - 4) + 3].emplace_back(100 * (iCh + 1) + 9, nullptr);
    mClusters[8 + 4 * (iCh - 4) + 3].emplace_back(100 * (iCh + 1) + 11, nullptr);
    mClusters[8 + 4 * (iCh - 4) + 3].emplace_back(100 * (iCh + 1) + 13, nullptr);
    mClusters[8 + 4 * (iCh - 4) + 3].emplace_back(100 * (iCh + 1) + 15, nullptr);
    mClusters[8 + 4 * (iCh - 4) + 3].emplace_back(100 * (iCh + 1) + 17, nullptr);
    mClusters[8 + 4 * (iCh - 4) + 3].emplace_back(100 * (iCh + 1) + 19, nullptr);
  }
}
 
//_________________________________________________________________________________________________
void TrackFinder::initField(float l3Current, float dipoleCurrent)
{
  mTrackFitter.initField(l3Current, dipoleCurrent);
}
 
//_________________________________________________________________________________________________
const std::list<Track>& TrackFinder::findTracks(gsl::span<const Cluster> clusters)
{
  std::unordered_map<int, std::list<const Cluster*>> clustersPerDE{};
  for (const auto& cluster : clusters) {
    clustersPerDE[cluster.getDEId()].emplace_back(&cluster);
  }
  return findTracks(clustersPerDE);
}
 
//_________________________________________________________________________________________________
const std::list<Track>& TrackFinder::findTracks(const std::unordered_map<int, std::list<const Cluster*>>& clusters)
{
 
  mTracks.clear();
  mStartTime = std::chrono::steady_clock::now();
 
  // fill the internal array of pointers to the list of clusters per DE
  for (auto& plane : mClusters) {
    for (auto& de : plane) {
      auto itDE = clusters.find(de.first);
      if (itDE == clusters.end()) {
        de.second = nullptr;
      } else {
        de.second = &(itDE->second);
      }
    }
  }
 
  // use the chamber resolution when fitting the tracks during the tracking
  mTrackFitter.useChamberResolution();
 
  try {
 
    // find track candidates on stations 4 and 5
    auto tStart = std::chrono::high_resolution_clock::now();
    findTrackCandidates();
    auto tEnd = std::chrono::high_resolution_clock::now();
    mTimeFindCandidates += tEnd - tStart;
    if (TrackerParam::Instance().moreCandidates) {
      tStart = std::chrono::high_resolution_clock::now();
      findMoreTrackCandidates();
      tEnd = std::chrono::high_resolution_clock::now();
      mTimeFindMoreCandidates += tEnd - tStart;
    }
    mNCandidates += mTracks.size();
    print("------ list of track candidates ------");
    printTracks();
 
    // track each candidate down to chamber 1 and remove it
    tStart = std::chrono::high_resolution_clock::now();
    for (auto itTrack = mTracks.begin(); itTrack != mTracks.end();) {
      std::unordered_map<int, std::unordered_set<uint32_t>> excludedClusters{};
      followTrackInChamber(itTrack, 5, 0, false, excludedClusters);
      print("findTracks: removing candidate at position #", getTrackIndex(itTrack));
      itTrack = mTracks.erase(itTrack);
    }
    tEnd = std::chrono::high_resolution_clock::now();
    mTimeFollowTracks += tEnd - tStart;
    print("------ list of tracks before improvement and cleaning ------");
    printTracks();
 
  } catch (exception const& e) {
    LOG(warning) << e.what() << " --> abort";
    mTracks.clear();
    return mTracks;
  }
 
  // improve the reconstructed tracks
  auto tStart = std::chrono::high_resolution_clock::now();
  improveTracks();
  auto tEnd = std::chrono::high_resolution_clock::now();
  mTimeImproveTracks += tEnd - tStart;
 
  // remove connected tracks in stations(1..) 3, 4 and 5
  tStart = std::chrono::high_resolution_clock::now();
  removeConnectedTracks(2, 4);
  tEnd = std::chrono::high_resolution_clock::now();
  mTimeCleanTracks += tEnd - tStart;
 
  // refine the tracks using cluster resolution or just finalize them
  if (TrackerParam::Instance().refineTracks) {
    tStart = std::chrono::high_resolution_clock::now();
    mTrackFitter.useClusterResolution();
    refineTracks();
    tEnd = std::chrono::high_resolution_clock::now();
    mTimeRefineTracks += tEnd - tStart;
  } else {
    finalize();
  }
 
  return mTracks;
}
 
//_________________________________________________________________________________________________
void TrackFinder::findTrackCandidates()
{
 
  // start by looking for candidates on station 5
  findTrackCandidatesInSt5();
 
  for (auto itTrack = mTracks.begin(); itTrack != mTracks.end();) {
 
    // prepare backward tracking if not already done
    if (!itTrack->hasCurrentParam()) {
      prepareBackwardTracking(itTrack, false);
    }
 
    // look for compatible clusters on station 4
    std::unordered_map<int, std::unordered_set<uint32_t>> excludedClusters{};
    auto itNewTrack = followTrackInChamber(itTrack, 7, 6, false, excludedClusters);
 
    // keep the current candidate only if no compatible cluster is found and the station is not requested
    if (!TrackerParam::Instance().requestStation[3] && excludedClusters.empty() && itTrack->areCurrentParamValid()) {
      ++itTrack;
    } else {
      print("findTrackCandidates: removing candidate at position #", getTrackIndex(itTrack));
      itTrack = mTracks.erase(itTrack);
      // prepare backward tracking for the new tracks
      for (; itNewTrack != mTracks.end() && itNewTrack != itTrack; ++itNewTrack) {
        prepareBackwardTracking(itNewTrack, false);
      }
    }
  }
 
  // list the cluster combinations already used in stations 4 and 5
  std::vector<std::array<uint32_t, 8>> usedClusters(mTracks.size());
  int iTrack(0);
  for (const auto& track : mTracks) {
    for (const auto& param : track) {
      int iCl = 2 * (param.getClusterPtr()->getChamberId() - 6) + param.getClusterPtr()->getDEId() % 2;
      usedClusters[iTrack][iCl] = param.getClusterPtr()->uid;
    }
    ++iTrack;
  }
 
  auto itLastCandidateFromSt5 = mTracks.empty() ? mTracks.end() : std::prev(mTracks.end());
 
  // then look for candidates on station 4
  findTrackCandidatesInSt4();
 
  auto itFirstCandidateOnSt4 = (itLastCandidateFromSt5 == mTracks.end()) ? mTracks.begin() : std::next(itLastCandidateFromSt5);
  for (auto itTrack = itFirstCandidateOnSt4; itTrack != mTracks.end();) {
 
    // prepare forward tracking if not already done
    if (!itTrack->hasCurrentParam()) {
      try {
        prepareForwardTracking(itTrack, true);
      } catch (exception const&) {
        print("findTrackCandidates: removing candidate at position #", getTrackIndex(itTrack));
        itTrack = mTracks.erase(itTrack);
        continue;
      }
    }
 
    // look for compatible clusters on each chamber of station 5 separately,
    // exluding those already attached to an identical candidate on station 4
    // (cases where both chambers of station 5 are fired should have been found in the first step)
    std::unordered_map<int, std::unordered_set<uint32_t>> excludedClusters{};
    if (!usedClusters.empty()) {
      std::array<uint32_t, 4> currentClusters{};
      for (const auto& param : *itTrack) {
        int iCl = 2 * (param.getClusterPtr()->getChamberId() - 6) + param.getClusterPtr()->getDEId() % 2;
        currentClusters[iCl] = param.getClusterPtr()->uid;
      }
      excludeClustersFromIdenticalTracks(currentClusters, usedClusters, excludedClusters);
    }
    auto itFirstNewTrack = followTrackInChamber(itTrack, 8, 8, false, excludedClusters);
    auto itNewTrack = followTrackInChamber(itTrack, 9, 9, false, excludedClusters);
    if (itFirstNewTrack == mTracks.end()) {
      itFirstNewTrack = itNewTrack;
    }
 
    // keep the current candidate only if no compatible cluster is found and the station is not requested
    if (!TrackerParam::Instance().requestStation[4] && excludedClusters.empty()) {
      itFirstNewTrack = itTrack;
      ++itTrack;
    } else {
      print("findTrackCandidates: removing candidate at position #", getTrackIndex(itTrack));
      itTrack = mTracks.erase(itTrack);
    }
 
    // refit the track(s) and prepare to continue the tracking in the backward direction
    while (itFirstNewTrack != mTracks.end() && itFirstNewTrack != itTrack) {
      try {
        prepareBackwardTracking(itFirstNewTrack, true);
        ++itFirstNewTrack;
      } catch (exception const&) {
        print("findTrackCandidates: removing candidate at position #", getTrackIndex(itFirstNewTrack));
        itFirstNewTrack = mTracks.erase(itFirstNewTrack);
      }
    }
  }
}
 
//_________________________________________________________________________________________________
void TrackFinder::findTrackCandidatesInSt5()
{
 
  print("--- find candidates in station 5 ---");
 
  for (int iPlaneCh9 = 27; iPlaneCh9 > 23; --iPlaneCh9) {
 
    for (int iPlaneCh10 = 28; iPlaneCh10 < 32; ++iPlaneCh10) {
 
      // skip candidates that have already been found starting from a previous combination
      bool skipUsedPairs = (iPlaneCh9 == 24 || iPlaneCh9 == 26 || iPlaneCh10 == 29 || iPlaneCh10 == 31);
 
      // find all valid candidates between these 2 planes
      auto itTrack = findTrackCandidates(iPlaneCh9, iPlaneCh10, skipUsedPairs, mTracks.begin());
 
      // stop here if overlaps have already been checked on both chambers
      if ((iPlaneCh9 == 24 || iPlaneCh9 == 26) && (iPlaneCh10 == 29 || iPlaneCh10 == 31)) {
        continue;
      }
 
      while (itTrack != mTracks.end()) {
 
        auto itNextTrack = std::next(itTrack);
 
        if (iPlaneCh10 == 28 || iPlaneCh10 == 30) {
 
          // if not already done, look for compatible clusters in the overlapping regions of chamber 10
          try {
            prepareForwardTracking(itTrack, true);
          } catch (exception const&) {
            print("findTrackCandidatesInSt5: removing candidate at position #", getTrackIndex(itTrack));
            itTrack = mTracks.erase(itTrack);
            continue;
          }
          auto itNewTrack = followTrackInOverlapDE(itTrack, itTrack->last().getClusterPtr()->getDEId(), iPlaneCh10 + 1);
 
          if (itNewTrack != mTracks.end()) {
 
            // remove the initial candidate if compatible cluster(s) are found
            print("findTrackCandidatesInSt5: removing candidate at position #", getTrackIndex(itTrack));
            itTrack = mTracks.erase(itTrack);
 
            // refit the track(s) with new attached cluster(s) and prepare to continue the tracking in the backward direction
            bool stop(false);
            while (!stop) {
              itTrack = std::prev(itTrack);
              if (itTrack == itNewTrack) {
                stop = true;
              }
              try {
                prepareBackwardTracking(itTrack, true);
              } catch (exception const&) {
                print("findTrackCandidatesInSt5: removing candidate at position #", getTrackIndex(itTrack));
                itTrack = mTracks.erase(itTrack);
              }
            }
          } else {
            // prepare to continue the tracking in the backward direction with the initial candidate
            prepareBackwardTracking(itTrack, false);
          }
        }
 
        if (iPlaneCh9 == 25 || iPlaneCh9 == 27) {
 
          while (itTrack != itNextTrack) {
 
            // if not already done, look for compatible clusters in the overlapping regions of chamber 9
            if (!itTrack->hasCurrentParam()) {
              prepareBackwardTracking(itTrack, false);
            }
            auto itNewTrack = followTrackInOverlapDE(itTrack, itTrack->first().getClusterPtr()->getDEId(), iPlaneCh9 - 1);
 
            // keep the initial candidate only if no compatible cluster is found
            if (itNewTrack == mTracks.end()) {
              ++itTrack;
            } else {
              print("findTrackCandidatesInSt5: removing candidate at position #", getTrackIndex(itTrack));
              itTrack = mTracks.erase(itTrack);
            }
          }
        }
 
        itTrack = itNextTrack;
      }
    }
  }
 
  // remove tracks out of limits now that overlaps have been checked
  for (auto itTrack = mTracks.begin(); itTrack != mTracks.end();) {
    if (itTrack->isRemovable()) {
      itTrack = mTracks.erase(itTrack);
    } else {
      ++itTrack;
    }
  }
}
 
//_________________________________________________________________________________________________
void TrackFinder::findTrackCandidatesInSt4()
{
 
  print("--- find candidates in station 4 ---");
 
  auto itLastCandidateFromSt5 = mTracks.empty() ? mTracks.end() : std::prev(mTracks.end());
 
  for (int iPlaneCh8 = 20; iPlaneCh8 < 24; ++iPlaneCh8) {
 
    for (int iPlaneCh7 = 19; iPlaneCh7 > 15; --iPlaneCh7) {
 
      // skip candidates that have already been found starting from a previous combination
      bool skipUsedPairs = (iPlaneCh7 == 18 || iPlaneCh7 == 16 || iPlaneCh8 == 21 || iPlaneCh8 == 23);
 
      // find all valid candidates between these 2 planes
      auto itFirstCandidateOnSt4 = (itLastCandidateFromSt5 == mTracks.end()) ? mTracks.begin() : std::next(itLastCandidateFromSt5);
      auto itTrack = findTrackCandidates(iPlaneCh7, iPlaneCh8, skipUsedPairs, itFirstCandidateOnSt4);
 
      // stop here if overlaps have already been checked on both chambers
      if ((iPlaneCh7 == 18 || iPlaneCh7 == 16) && (iPlaneCh8 == 21 || iPlaneCh8 == 23)) {
        continue;
      }
 
      while (itTrack != mTracks.end()) {
 
        auto itNextTrack = std::next(itTrack);
 
        if (iPlaneCh7 == 19 || iPlaneCh7 == 17) {
 
          // if not already done, look for compatible clusters in the overlapping regions of chamber 7
          prepareBackwardTracking(itTrack, false);
          auto itNewTrack = followTrackInOverlapDE(itTrack, itTrack->first().getClusterPtr()->getDEId(), iPlaneCh7 - 1);
 
          // keep the initial candidate only if no compatible cluster is found
          if (itNewTrack != mTracks.end()) {
            print("findTrackCandidatesInSt4: removing candidate at position #", getTrackIndex(itTrack));
            mTracks.erase(itTrack);
            itTrack = itNewTrack;
          }
        }
 
        while (itTrack != itNextTrack) {
 
          // for every tracks, prepare to continue the tracking in the forward direction
          try {
            prepareForwardTracking(itTrack, true);
          } catch (exception const&) {
            print("findTrackCandidatesInSt4: removing candidate at position #", getTrackIndex(itTrack));
            itTrack = mTracks.erase(itTrack);
            continue;
          }
 
          if (iPlaneCh8 == 20 || iPlaneCh8 == 22) {
 
            // if not already done, look for compatible clusters in the overlapping regions of chamber 8
            auto itNewTrack = followTrackInOverlapDE(itTrack, itTrack->last().getClusterPtr()->getDEId(), iPlaneCh8 + 1);
 
            // keep the initial candidate only if no compatible cluster is found
            if (itNewTrack == mTracks.end()) {
              ++itTrack;
            } else {
              // prepare to continue the tracking of the new tracks from the last attached cluster
              for (; itNewTrack != itTrack; ++itNewTrack) {
                prepareForwardTracking(itNewTrack, false);
              }
              print("findTrackCandidatesInSt4: removing candidate at position #", getTrackIndex(itTrack));
              itTrack = mTracks.erase(itTrack);
            }
          } else {
            ++itTrack;
          }
        }
      }
    }
  }
 
  // remove tracks out of limits now that overlaps have been checked
  auto itTrack = (itLastCandidateFromSt5 == mTracks.end()) ? mTracks.begin() : ++itLastCandidateFromSt5;
  while (itTrack != mTracks.end()) {
    if (itTrack->isRemovable()) {
      itTrack = mTracks.erase(itTrack);
    } else {
      ++itTrack;
    }
  }
}
 
//_________________________________________________________________________________________________
void TrackFinder::findMoreTrackCandidates()
{
 
  print("--- find more candidates ---");
 
  auto itLastCandidate = mTracks.empty() ? mTracks.end() : std::prev(mTracks.end());
 
  for (int iPlaneSt4 = 23; iPlaneSt4 > 15; --iPlaneSt4) {
 
    for (int iPlaneSt5 = 24; iPlaneSt5 < 32; ++iPlaneSt5) {
 
      // find all valid candidates between these 2 planes
      auto itTrack = findTrackCandidates(iPlaneSt4, iPlaneSt5, true, mTracks.begin());
 
      // stop here if overlaps have already been checked on both chambers
      if ((iPlaneSt4 % 2 == 0) && (iPlaneSt5 % 2 == 1)) {
        continue;
      }
 
      while (itTrack != mTracks.end()) {
 
        auto itNextTrack = std::next(itTrack);
 
        if (iPlaneSt5 % 2 == 0) {
 
          // if not already done, look for compatible clusters in the overlapping regions of that chamber in station 5
          try {
            prepareForwardTracking(itTrack, true);
          } catch (exception const&) {
            print("findMoreTrackCandidates: removing candidate at position #", getTrackIndex(itTrack));
            itTrack = mTracks.erase(itTrack);
            continue;
          }
          auto itNewTrack = followTrackInOverlapDE(itTrack, itTrack->last().getClusterPtr()->getDEId(), iPlaneSt5 + 1);
 
          if (itNewTrack != mTracks.end()) {
 
            // remove the initial candidate if compatible cluster(s) are found
            print("findMoreTrackCandidates: removing candidate at position #", getTrackIndex(itTrack));
            itTrack = mTracks.erase(itTrack);
 
            // refit the track(s) with new cluster(s) and prepare to continue the tracking in the backward direction
            bool stop(false);
            while (!stop) {
              itTrack = std::prev(itTrack);
              if (itTrack == itNewTrack) {
                stop = true;
              }
              try {
                prepareBackwardTracking(itTrack, true);
              } catch (exception const&) {
                print("findMoreTrackCandidates: removing candidate at position #", getTrackIndex(itTrack));
                itTrack = mTracks.erase(itTrack);
              }
            }
          } else {
            // prepare to continue the tracking in the backward direction with the initial candidate
            prepareBackwardTracking(itTrack, false);
          }
        }
 
        if (iPlaneSt4 % 2 == 1) {
 
          while (itTrack != itNextTrack) {
 
            // if not already done, look for compatible clusters in the overlapping regions of that chamber in station 4
            if (!itTrack->hasCurrentParam()) {
              prepareBackwardTracking(itTrack, false);
            }
            auto itNewTrack = followTrackInOverlapDE(itTrack, itTrack->first().getClusterPtr()->getDEId(), iPlaneSt4 - 1);
 
            // keep the initial candidate only if no compatible cluster is found
            if (itNewTrack == mTracks.end()) {
              ++itTrack;
            } else {
              print("findMoreTrackCandidates: removing candidate at position #", getTrackIndex(itTrack));
              itTrack = mTracks.erase(itTrack);
            }
          }
        }
 
        itTrack = itNextTrack;
      }
    }
  }
 
  // remove tracks out of limits now that overlaps have been checked
  // and make sure every new tracks are prepared to continue the tracking in the backward direction
  auto itTrack = (itLastCandidate == mTracks.end()) ? mTracks.begin() : ++itLastCandidate;
  while (itTrack != mTracks.end()) {
    if (itTrack->isRemovable()) {
      itTrack = mTracks.erase(itTrack);
    } else {
      if (!itTrack->hasCurrentParam()) {
        prepareBackwardTracking(itTrack, false);
      }
      ++itTrack;
    }
  }
}
 
//_________________________________________________________________________________________________
std::list<Track>::iterator TrackFinder::findTrackCandidates(int plane1, int plane2, bool skipUsedPairs, const std::list<Track>::iterator& itFirstTrack)
{
 
  const auto& trackerParam = TrackerParam::Instance();
 
  // maximum impact parameter dispersion**2 due to MCS in chambers
  double impactMCS2(0.);
  int chamber1 = getChamberId(plane1);
  for (int iCh = 0; iCh <= chamber1; ++iCh) {
    impactMCS2 += SDefaultChamberZ[iCh] * SDefaultChamberZ[iCh] * mMaxMCSAngle2[iCh];
  }
 
  // create an iterator to the last track of the list before adding new ones
  auto itTrack = mTracks.empty() ? mTracks.end() : std::prev(mTracks.end());
 
  // list the cluster combinations used on the chambers corresponding to plane1 and plane2
  int chamber2 = getChamberId(plane2);
  std::vector<std::array<uint32_t, 4>> usedClusters{};
  if (skipUsedPairs) {
    usedClusters.reserve(mTracks.size());
    for (auto itTrk = itFirstTrack; itTrk != mTracks.end(); ++itTrk) {
      usedClusters.push_back({0, 0, 0, 0});
      for (const auto& param : *itTrk) {
        int ch = param.getClusterPtr()->getChamberId();
        if (ch == chamber1) {
          usedClusters.back()[param.getClusterPtr()->getDEId() % 2] = param.getClusterPtr()->uid;
        } else if (ch == chamber2) {
          usedClusters.back()[2 + param.getClusterPtr()->getDEId() % 2] = param.getClusterPtr()->uid;
        }
      }
    }
  }
 
  for (auto& de1 : mClusters[plane1]) {
 
    // skip DE without cluster
    if (de1.second == nullptr) {
      continue;
    }
 
    for (const auto cluster1 : *de1.second) {
 
      double z1 = cluster1->getZ();
 
      for (auto& de2 : mClusters[plane2]) {
 
        // skip DE without cluster
        if (de2.second == nullptr) {
          continue;
        }
 
        for (const auto cluster2 : *de2.second) {
 
          // skip combinations of clusters already part of a track if requested
          if (skipUsedPairs && areUsed(*cluster1, *cluster2, usedClusters)) {
            continue;
          }
 
          double z2 = cluster2->getZ();
          double dZ = z1 - z2;
 
          // check if non bending impact parameter is within tolerances
          double nonBendingSlope = (cluster1->getX() - cluster2->getX()) / dZ;
          double nonBendingImpactParam = TMath::Abs(cluster1->getX() - cluster1->getZ() * nonBendingSlope);
          double nonBendingImpactParamErr = TMath::Sqrt((z1 * z1 * mChamberResolutionX2 + z2 * z2 * mChamberResolutionX2) / dZ / dZ + impactMCS2);
          if ((nonBendingImpactParam - trackerParam.sigmaCutForTracking * nonBendingImpactParamErr) > (3. * trackerParam.nonBendingVertexDispersion)) {
            continue;
          }
 
          double bendingSlope = (cluster1->getY() - cluster2->getY()) / dZ;
          if (TrackExtrap::isFieldON()) { // depending whether the field is ON or OFF
            // check if bending momentum is within tolerances
            double bendingImpactParam = cluster1->getY() - cluster1->getZ() * bendingSlope;
            double bendingImpactParamErr2 = (z1 * z1 * mChamberResolutionY2 + z2 * z2 * mChamberResolutionY2) / dZ / dZ + impactMCS2;
            double bendingMomentum = TMath::Abs(TrackExtrap::getBendingMomentumFromImpactParam(bendingImpactParam));
            double bendingMomentumErr = TMath::Sqrt((mBendingVertexDispersion2 + bendingImpactParamErr2) / bendingImpactParam / bendingImpactParam + 0.01) * bendingMomentum;
            if ((bendingMomentum + 3. * bendingMomentumErr) < SMinBendingMomentum) {
              continue;
            }
          } else {
            // or check if bending impact parameter is within tolerances
            double bendingImpactParam = TMath::Abs(cluster1->getY() - cluster1->getZ() * bendingSlope);
            double bendingImpactParamErr = TMath::Sqrt((z1 * z1 * mChamberResolutionY2 + z2 * z2 * mChamberResolutionY2) / dZ / dZ + impactMCS2);
            if ((bendingImpactParam - trackerParam.sigmaCutForTracking * bendingImpactParamErr) > (3. * trackerParam.bendingVertexDispersion)) {
              continue;
            }
          }
 
          // create a new track candidate
          createTrack(*cluster1, *cluster2);
        }
      }
    }
  }
 
  return (itTrack == mTracks.end()) ? mTracks.begin() : ++itTrack;
}
 
//_________________________________________________________________________________________________
std::list<Track>::iterator TrackFinder::followTrackInOverlapDE(const std::list<Track>::iterator& itTrack, int currentDE, int plane)
{
 
  print("followTrackInOverlapDE: follow track #", getTrackIndex(itTrack), " currently at DE ", currentDE, " to plane ", plane);
  printTrack(*itTrack);
 
  // the current track parameters must be set
  assert(itTrack->hasCurrentParam());
 
  auto itNewTrack(itTrack);
 
  const TrackParam& currentParam = itTrack->getCurrentParam();
  int currentChamber = itTrack->getCurrentChamber();
 
  // loop over all DEs of plane
  TrackParam paramAtCluster{};
  for (auto& de : mClusters[plane]) {
 
    // skip DE without cluster
    if (de.second == nullptr) {
      continue;
    }
 
    // skip DE that do not overlap with the current DE
    if (de.first % 100 != (currentDE % 100 + 1) % SNDE[currentChamber] && de.first % 100 != (currentDE % 100 - 1 + SNDE[currentChamber]) % SNDE[currentChamber]) {
      continue;
    }
 
    // look for cluster candidate in this DE
    for (const auto cluster : *de.second) {
 
      // try to add the current cluster
      if (!isCompatible(currentParam, *cluster, paramAtCluster)) {
        continue;
      }
 
      // duplicate the track and add the new cluster
      itNewTrack = addTrack(itNewTrack, *itTrack);
      print("followTrackInOverlapDE: duplicating candidate at position #", getTrackIndex(itNewTrack), " to add cluster ", cluster->getIdAsString());
      itNewTrack->addParamAtCluster(paramAtCluster);
 
      // tag the track as removable (if it is not already the case) if it is out of limits
      if (!itNewTrack->isRemovable() && !isAcceptable(paramAtCluster)) {
        itNewTrack->removable();
      }
    }
  }
 
  return (itNewTrack == itTrack) ? mTracks.end() : itNewTrack;
}
 
//_________________________________________________________________________________________________
std::list<Track>::iterator TrackFinder::followTrackInChamber(std::list<Track>::iterator& itTrack,
                                                             int chamber, int lastChamber, bool canSkip,
                                                             std::unordered_map<int, std::unordered_set<uint32_t>>& excludedClusters)
{
 
  // list of (half-)planes, 2 or 4 per chamber, ordered according to the direction of propagation,
  // which is forward when going to station 5 and backward otherwise with the present algorithm
  static constexpr int plane[10][4] = {{1, 0, -1, -1}, {3, 2, -1, -1}, {5, 4, -1, -1}, {7, 6, -1, -1}, {11, 10, 9, 8}, {15, 14, 13, 12}, {19, 18, 17, 16}, {23, 22, 21, 20}, {24, 25, 26, 27}, {28, 29, 30, 31}};
 
  print("followTrackInChamber: follow track #", getTrackIndex(itTrack), " to chamber ", chamber + 1, " up to chamber ", lastChamber + 1);
 
  // the current track parameters must be set at a different chamber and valid
  if (!itTrack->areCurrentParamValid() || chamber == itTrack->getCurrentChamber()) {
    return mTracks.end();
  }
 
  std::chrono::duration<double> currentTrackingDuration = std::chrono::steady_clock::now() - mStartTime;
  if (currentTrackingDuration.count() > TrackerParam::Instance().maxTrackingDuration) {
    mErrorMap.add(ErrorType::Tracking_TooLong, 0, 0);
    throw length_error(string("Tracking is taking too long (") + std::round(currentTrackingDuration.count()) + " s)");
  }
 
  // determine whether the chamber is the first one reached on the station
  int currentChamber = itTrack->getCurrentChamber();
  bool isFirstOnStation = ((chamber < currentChamber && chamber % 2 == 1) || (chamber > currentChamber && chamber % 2 == 0));
 
  // follow the track in the 2 planes or 4 half-planes of the chamber
  auto itFirstNewTrack = followTrackInChamber(itTrack, plane[chamber][0], plane[chamber][1], lastChamber, excludedClusters);
  if (chamber > 3) {
    auto itNewTrack = followTrackInChamber(itTrack, plane[chamber][2], plane[chamber][3], lastChamber, excludedClusters);
    if (itFirstNewTrack == mTracks.end()) {
      itFirstNewTrack = itNewTrack;
    }
  }
 
  // add MCS effects in that chamber before going further with this track or stop here if the track could not reach that chamber
  if (itTrack->areCurrentParamValid()) {
    TrackExtrap::addMCSEffect(itTrack->getCurrentParam(), SChamberThicknessInX0[chamber], -1.);
  } else {
    return itFirstNewTrack;
  }
 
  if (chamber != lastChamber) {
 
    // save the current track parameters before going to the next chamber
    TrackParam currentParam = itTrack->getCurrentParam();
 
    // consider the possibility to skip the chamber if it is the first one of the station or if we know we can skip it,
    // i.e. if a compatible cluster has been found on the first chamber and none has been found on the second
    if (isFirstOnStation || (canSkip && excludedClusters.empty())) {
      int nextChamber = (chamber > lastChamber) ? chamber - 1 : chamber + 1;
      auto itNewTrack = followTrackInChamber(itTrack, nextChamber, lastChamber, false, excludedClusters);
      if (itFirstNewTrack == mTracks.end()) {
        itFirstNewTrack = itNewTrack;
      }
    }
 
    // consider the possibility to skip the entire station if not requested and not the last one
    if (isFirstOnStation && !TrackerParam::Instance().requestStation[chamber / 2] && chamber / 2 != lastChamber / 2) {
      int nextChamber = (chamber > lastChamber) ? chamber - 2 : chamber + 2;
      auto itNewTrack = followTrackInChamber(itTrack, nextChamber, lastChamber, false, excludedClusters);
      if (itFirstNewTrack == mTracks.end()) {
        itFirstNewTrack = itNewTrack;
      }
    }
 
    // reset the current track parameters to the ones at that chamber if needed
    // (not sure it is needed at all but that way it is clear what the current track parameters are at the end of this function)
    if (itTrack->getCurrentChamber() != chamber) {
      setCurrentParam(*itTrack, currentParam, chamber);
    }
  } else {
 
    // add a new track if a cluster has been found on the first chamber of the station but not on the second and last one
    // or if one reaches station 1 and it is not requested, whether a cluster has been found on it or not
    if ((!isFirstOnStation && canSkip && excludedClusters.empty()) ||
        (chamber / 2 == 0 && !TrackerParam::Instance().requestStation[0] && (isFirstOnStation || !canSkip))) {
      auto itNewTrack = addTrack(itTrack, *itTrack);
      if (itFirstNewTrack == mTracks.end()) {
        itFirstNewTrack = itNewTrack;
      }
      print("followTrackInChamber: duplicating candidate at position #", getTrackIndex(itFirstNewTrack));
    }
  }
 
  return itFirstNewTrack;
}
 
//_________________________________________________________________________________________________
std::list<Track>::iterator TrackFinder::followTrackInChamber(std::list<Track>::iterator& itTrack,
                                                             int plane1, int plane2, int lastChamber,
                                                             std::unordered_map<int, std::unordered_set<uint32_t>>& excludedClusters)
{
 
  print("followTrackInChamber: follow track #", getTrackIndex(itTrack), " to planes ", plane1, " and ", plane2, " up to chamber ", lastChamber + 1);
  printTrack(*itTrack);
 
  // the current track parameters must be set and valid
  if (!itTrack->areCurrentParamValid()) {
    return mTracks.end();
  }
 
  auto itFirstNewTrack(mTracks.end());
 
  // add MCS effects in the missing chambers if any. Update the current parameters in the process
  int chamber = getChamberId(plane1);
  if ((chamber < itTrack->getCurrentChamber() - 1 || chamber > itTrack->getCurrentChamber() + 1) &&
      !propagateCurrentParam(*itTrack, (chamber < itTrack->getCurrentChamber()) ? chamber + 1 : chamber - 1)) {
    return mTracks.end();
  }
 
  // extrapolate the candidate to the chamber if not already there
  TrackParam paramAtChamber = itTrack->getCurrentParam();
  if (itTrack->getCurrentChamber() != chamber && !TrackExtrap::extrapToZCov(paramAtChamber, SDefaultChamberZ[chamber], true)) {
    itTrack->invalidateCurrentParam();
    return mTracks.end();
  }
 
  // determine the next chamber to go to, if lastChamber is not yet reached
  int nextChamber(-1);
  if (chamber > lastChamber) {
    nextChamber = chamber - 1;
  } else if (chamber < lastChamber) {
    nextChamber = chamber + 1;
  }
 
  // loop over all DEs of plane1
  TrackParam paramAtCluster1{};
  TrackParam currentParamAtCluster1{};
  TrackParam paramAtCluster2{};
  std::unordered_map<int, std::unordered_set<uint32_t>> newExcludedClusters{};
  for (auto& de1 : mClusters[plane1]) {
 
    // skip DE without cluster
    if (de1.second == nullptr) {
      continue;
    }
 
    // get the list of excluded clusters for the DE
    auto itExcludedClusters = excludedClusters.find(de1.first);
    bool hasExcludedClusters = (itExcludedClusters != excludedClusters.end());
 
    // look for cluster candidate in this DE
    for (const auto cluster1 : *de1.second) {
 
      // skip excluded clusters
      if (hasExcludedClusters && itExcludedClusters->second.count(cluster1->uid) > 0) {
        continue;
      }
 
      // try to add the current cluster
      if (!isCompatible(paramAtChamber, *cluster1, paramAtCluster1)) {
        continue;
      }
 
      // add it to the list of excluded clusters for this candidate
      excludedClusters[de1.first].emplace(cluster1->uid);
 
      // skip tracks out of limits, but after checking for overlaps
      bool isAcceptableAtCluster1 = isAcceptable(paramAtCluster1);
 
      // save the current parameters at cluster1, reset the propagator and add MCS effects before going to plane2
      currentParamAtCluster1 = paramAtCluster1;
      currentParamAtCluster1.resetPropagator();
      TrackExtrap::addMCSEffect(currentParamAtCluster1, SChamberThicknessInX0[chamber], -1.);
 
      // loop over all DEs of plane2
      bool cluster2Found(false);
      for (auto& de2 : mClusters[plane2]) {
 
        // skip DE without cluster
        if (de2.second == nullptr) {
          continue;
        }
 
        // skip DE that do not overlap with the DE of plane1
        if (de2.first % 100 != (de1.first % 100 + 1) % SNDE[chamber] && de2.first % 100 != (de1.first % 100 - 1 + SNDE[chamber]) % SNDE[chamber]) {
          continue;
        }
 
        // look for cluster candidate in this DE
        for (const auto cluster2 : *de2.second) {
 
          // try to add the current cluster
          if (!isCompatible(currentParamAtCluster1, *cluster2, paramAtCluster2)) {
            continue;
          }
 
          cluster2Found = true;
 
          // add it to the list of excluded clusters for this candidate
          excludedClusters[de2.first].emplace(cluster2->uid);
 
          // skip tracks out of limits
          if (!isAcceptableAtCluster1 || !isAcceptable(paramAtCluster2)) {
            continue;
          }
 
          // continue the tracking to the next chambers and attach the 2 clusters to the new tracks if any
          auto itNewTrack = addClustersAndFollowTrack(itTrack, paramAtCluster1, &paramAtCluster2, nextChamber, lastChamber, newExcludedClusters);
          if (itFirstNewTrack == mTracks.end()) {
            itFirstNewTrack = itNewTrack;
          }
 
          // transfert the list of new excluded clusters to the full list for the initial candidate
          moveClusters(newExcludedClusters, excludedClusters);
        }
      }
 
      if (!cluster2Found && isAcceptableAtCluster1) {
 
        // continue the tracking with only cluster1 if no compatible cluster is found on plane2 and the track stays within limits
        auto itNewTrack = addClustersAndFollowTrack(itTrack, paramAtCluster1, nullptr, nextChamber, lastChamber, newExcludedClusters);
        if (itFirstNewTrack == mTracks.end()) {
          itFirstNewTrack = itNewTrack;
        }
 
        // transfert the list of new excluded clusters to the full list for the initial candidate
        moveClusters(newExcludedClusters, excludedClusters);
      }
    }
  }
 
  // loop over all DEs of plane2
  for (auto& de2 : mClusters[plane2]) {
 
    // skip DE without cluster
    if (de2.second == nullptr) {
      continue;
    }
 
    // get the list of excluded clusters for the DE
    auto itExcludedClusters = excludedClusters.find(de2.first);
    bool hasExcludedClusters = (itExcludedClusters != excludedClusters.end());
 
    // look for cluster candidate in this DE
    for (const auto cluster2 : *de2.second) {
 
      // skip excluded clusters (in particular the ones already attached together with a cluster on plane1)
      if (hasExcludedClusters && itExcludedClusters->second.count(cluster2->uid) > 0) {
        continue;
      }
 
      // try to add the current cluster
      if (!isCompatible(paramAtChamber, *cluster2, paramAtCluster2)) {
        continue;
      }
 
      // add it to the list of excluded clusters for this candidate
      excludedClusters[de2.first].emplace(cluster2->uid);
 
      // skip tracks out of limits
      if (!isAcceptable(paramAtCluster2)) {
        continue;
      }
 
      // continue the tracking to the next chambers and attach the cluster to the new tracks if any
      auto itNewTrack = addClustersAndFollowTrack(itTrack, paramAtCluster2, nullptr, nextChamber, lastChamber, newExcludedClusters);
      if (itFirstNewTrack == mTracks.end()) {
        itFirstNewTrack = itNewTrack;
      }
 
      // transfert the list of new excluded clusters to the full list for the initial candidate
      moveClusters(newExcludedClusters, excludedClusters);
    }
  }
 
  // reset the current parameters to the ones at that chamber if needed, not adding MCS effects yet
  if (itTrack->getCurrentChamber() != chamber) {
    setCurrentParam(*itTrack, paramAtChamber, chamber);
  }
 
  return itFirstNewTrack;
}
 
//_________________________________________________________________________________________________
std::list<Track>::iterator TrackFinder::addClustersAndFollowTrack(std::list<Track>::iterator& itTrack, const TrackParam& paramAtCluster1,
                                                                  const TrackParam* paramAtCluster2, int nextChamber, int lastChamber,
                                                                  std::unordered_map<int, std::unordered_set<uint32_t>>& excludedClusters)
{
 
  // the list of excluded clusters must be empty here as new cluster(s) are being attached to the candidate
  assert(excludedClusters.empty());
 
  auto itFirstNewTrack(mTracks.end());
 
  if (nextChamber >= 0) {
 
    // the tracking continues from paramAtCluster2, if any, or from paramAtCluster1
    if (paramAtCluster2) {
      print("addClustersAndFollowTrack: 2 clusters found (", paramAtCluster1.getClusterPtr()->getIdAsString(), " and ",
            paramAtCluster2->getClusterPtr()->getIdAsString(), "). Continuing the tracking of candidate #", getTrackIndex(itTrack));
      setCurrentParam(*itTrack, *paramAtCluster2, paramAtCluster2->getClusterPtr()->getChamberId());
    } else {
      print("addClustersAndFollowTrack: 1 cluster found (", paramAtCluster1.getClusterPtr()->getIdAsString(),
            "). Continuing the tracking of candidate #", getTrackIndex(itTrack));
      setCurrentParam(*itTrack, paramAtCluster1, paramAtCluster1.getClusterPtr()->getChamberId());
    }
 
    // follow the track to the next chamber, which can be skipped if it is on the same station
    bool canSkip = (nextChamber / 2 == paramAtCluster1.getClusterPtr()->getChamberId() / 2);
    auto itNewTrack = followTrackInChamber(itTrack, nextChamber, lastChamber, canSkip, excludedClusters);
    itFirstNewTrack = itNewTrack;
 
    // attach the current cluster(s) to every new tracks found
    if (itNewTrack != mTracks.end()) {
      while (itNewTrack != itTrack) {
        itNewTrack->addParamAtCluster(paramAtCluster1);
        if (paramAtCluster2) {
          itNewTrack->addParamAtCluster(*paramAtCluster2);
          print("addClustersAndFollowTrack: add to the candidate at position #", getTrackIndex(itNewTrack),
                " clusters ", paramAtCluster1.getClusterPtr()->getIdAsString(), " and ", paramAtCluster2->getClusterPtr()->getIdAsString());
        } else {
          print("addClustersAndFollowTrack: add to the candidate at position #", getTrackIndex(itNewTrack),
                " cluster ", paramAtCluster1.getClusterPtr()->getIdAsString());
        }
        ++itNewTrack;
      }
    }
 
  } else {
 
    // or duplicate the track and add the new cluster(s)
    itFirstNewTrack = addTrack(itTrack, *itTrack);
    itFirstNewTrack->addParamAtCluster(paramAtCluster1);
    if (paramAtCluster2) {
      itFirstNewTrack->addParamAtCluster(*paramAtCluster2);
      print("addClustersAndFollowTrack: duplicating candidate at position #", getTrackIndex(itFirstNewTrack), " to add 2 clusters (",
            paramAtCluster1.getClusterPtr()->getIdAsString(), " and ", paramAtCluster2->getClusterPtr()->getIdAsString(), ")");
    } else {
      print("addClustersAndFollowTrack: duplicating candidate at position #", getTrackIndex(itFirstNewTrack),
            " to add 1 cluster (", paramAtCluster1.getClusterPtr()->getIdAsString(), ")");
    }
  }
 
  return itFirstNewTrack;
}
 
//_________________________________________________________________________________________________
void TrackFinder::improveTracks()
{
 
  for (auto itTrack = mTracks.begin(); itTrack != mTracks.end();) {
 
    bool removeTrack(false);
 
    // At the first step, only run the smoother
    auto itStartingParam = std::prev(itTrack->rend());
 
    while (true) {
 
      // Refit the part of the track affected by the cluster removal, run the smoother, but do not finalize
      try {
        mTrackFitter.fit(*itTrack, true, false, (itStartingParam == itTrack->rbegin()) ? nullptr : &itStartingParam);
      } catch (exception const&) {
        removeTrack = true;
        break;
      }
 
      // Identify removable clusters
      itTrack->tagRemovableClusters(requestedStationMask(), !TrackerParam::Instance().moreCandidates);
 
      // Look for the cluster with the worst local chi2
      double worstLocalChi2(-1.);
      auto itWorstParam(itTrack->end());
      for (auto itParam = itTrack->begin(); itParam != itTrack->end(); ++itParam) {
        if (itParam->getLocalChi2() > worstLocalChi2) {
          worstLocalChi2 = itParam->getLocalChi2();
          itWorstParam = itParam;
        }
      }
 
      // If the worst chi2 is under requirement then the track is improved
      if (worstLocalChi2 < mMaxChi2ForImprovement) {
        break;
      }
 
      // If the worst cluster is not removable then the track cannot be improved
      if (!itWorstParam->isRemovable()) {
        removeTrack = true;
        break;
      }
 
      // Remove the worst cluster
      auto itNextParam = itTrack->removeParamAtCluster(itWorstParam);
 
      // Decide from where to refit the track: from the cluster next the one suppressed or
      // from scratch if the removed cluster was used to compute the tracking seed
      itStartingParam = itTrack->rbegin();
      auto itNextToNextParam = (itNextParam == itTrack->end()) ? itNextParam : std::next(itNextParam);
      while (itNextToNextParam != itTrack->end()) {
        if (itNextToNextParam->getClusterPtr()->getChamberId() != itNextParam->getClusterPtr()->getChamberId()) {
          itStartingParam = std::make_reverse_iterator(++itNextParam);
          break;
        }
        ++itNextToNextParam;
      }
    }
 
    // Remove the track if it couldn't be improved
    if (removeTrack) {
      print("improveTracks: removing candidate at position #", getTrackIndex(itTrack));
      itTrack = mTracks.erase(itTrack);
    } else {
      ++itTrack;
    }
  }
}
 
//_________________________________________________________________________________________________
void TrackFinder::removeConnectedTracks(int stMin, int stMax)
{
 
  if (mTracks.size() < 2) {
    return;
  }
 
  int chMin = 2 * stMin;
  int chMax = 2 * stMax + 1;
  int nPlane = 2 * (chMax - chMin + 1);
 
  // first loop to fill the arrays of cluster Ids, number of fired chambers and normalized chi2
  std::vector<uint32_t> ClIds(nPlane * mTracks.size());
  std::vector<uint8_t> nFiredCh(mTracks.size());
  std::vector<double> nChi2(mTracks.size());
  int previousCh(-1);
  int iTrack(0);
  for (auto itTrack = mTracks.begin(); itTrack != mTracks.end(); ++itTrack, ++iTrack) {
    for (auto itParam = itTrack->rbegin(); itParam != itTrack->rend(); ++itParam) {
      int ch = itParam->getClusterPtr()->getChamberId();
      if (ch != previousCh) {
        ++nFiredCh[iTrack];
        previousCh = ch;
      }
      if (ch >= chMin && ch <= chMax) {
        ClIds[nPlane * iTrack + 2 * (ch - chMin) + itParam->getClusterPtr()->getDEId() % 2] = itParam->getClusterPtr()->uid;
      }
    }
    nChi2[iTrack] = itTrack->first().getTrackChi2() / (itTrack->getNDF() - 1);
  }
 
  // second loop to tag the tracks to remove
  std::vector<bool> remove(mTracks.size(), false);
  int iindex = ClIds.size() - 1;
  for (int iTrack1 = mTracks.size() - 1; iTrack1 > -1; --iTrack1, iindex -= nPlane) {
    int jindex = iindex - nPlane;
    for (int iTrack2 = iTrack1 - 1; iTrack2 > -1; --iTrack2) {
      for (int iPlane = nPlane; iPlane > 0; --iPlane) {
        if (ClIds[iindex] == ClIds[jindex] && ClIds[iindex] > 0) {
          if (nFiredCh[iTrack2] > nFiredCh[iTrack1] ||
              (nFiredCh[iTrack2] == nFiredCh[iTrack1] && nChi2[iTrack2] < nChi2[iTrack1])) {
            remove[iTrack1] = true;
          } else {
            remove[iTrack2] = true;
          }
          iindex -= iPlane;
          jindex -= iPlane;
          break;
        }
        --iindex;
        --jindex;
      }
      iindex += nPlane;
    }
  }
 
  // third loop to remove them. That way all combinations are tested.
  iTrack = 0;
  for (auto itTrack = mTracks.begin(); itTrack != mTracks.end(); ++iTrack) {
    if (remove[iTrack]) {
      print("removeConnectedTracks: removing candidate at position #", getTrackIndex(itTrack));
      itTrack = mTracks.erase(itTrack);
    } else {
      ++itTrack;
    }
  }
}
 
//_________________________________________________________________________________________________
void TrackFinder::refineTracks()
{
 
  for (auto itTrack = mTracks.begin(); itTrack != mTracks.end();) {
    try {
      mTrackFitter.fit(*itTrack);
      ++itTrack;
    } catch (exception const&) {
      print("refineTracks: removing candidate at position #", getTrackIndex(itTrack));
      itTrack = mTracks.erase(itTrack);
    }
  }
}
 
//_________________________________________________________________________________________________
void TrackFinder::finalize()
{
  for (auto& track : mTracks) {
    for (auto& param : track) {
      param.setParameters(param.getSmoothParameters());
      param.setCovariances(param.getSmoothCovariances());
    }
  }
}
 
//_________________________________________________________________________________________________
void TrackFinder::createTrack(const Cluster& cl1, const Cluster& cl2)
{
 
  if (mTracks.size() >= TrackerParam::Instance().maxCandidates) {
    mErrorMap.add(ErrorType::Tracking_TooManyCandidates, 0, 0);
    throw length_error(string("Too many track candidates (") + mTracks.size() + ")");
  }
 
  // create the track and the trackParam at each cluster
  Track& track = mTracks.emplace_back();
  track.createParamAtCluster(cl2);
  track.createParamAtCluster(cl1);
  print("createTrack: creating candidate at position #", getTrackIndex(std::prev(mTracks.end())),
        " with clusters ", cl1.getIdAsString(), " and ", cl2.getIdAsString());
 
  // fit the track using the Kalman filter
  try {
    mTrackFitter.fit(track, false);
  } catch (exception const&) {
    print("... fit failed --> removing it");
    mTracks.erase(std::prev(mTracks.end()));
  }
}
 
//_________________________________________________________________________________________________
std::list<Track>::iterator TrackFinder::addTrack(const std::list<Track>::iterator& pos, const Track& track)
{
  if (mTracks.size() >= TrackerParam::Instance().maxCandidates) {
    mErrorMap.add(ErrorType::Tracking_TooManyCandidates, 0, 0);
    throw length_error(string("Too many track candidates (") + mTracks.size() + ")");
  }
  return mTracks.emplace(pos, track);
}
 
//_________________________________________________________________________________________________
bool TrackFinder::isAcceptable(const TrackParam& param) const
{
 
  // impact parameter dispersion**2 due to MCS in chambers
  int chamber = param.getClusterPtr()->getChamberId();
  double impactMCS2(0.);
  if (TrackExtrap::isFieldON() && chamber < 6) {
    // track momentum is known
    for (int iCh = 0; iCh <= chamber; ++iCh) {
      impactMCS2 += SDefaultChamberZ[iCh] * SDefaultChamberZ[iCh] * TrackExtrap::getMCSAngle2(param, SChamberThicknessInX0[iCh], 1.);
    }
  } else {
    // track momentum is unknown
    for (Int_t iCh = 0; iCh <= chamber; ++iCh) {
      impactMCS2 += SDefaultChamberZ[iCh] * SDefaultChamberZ[iCh] * mMaxMCSAngle2[iCh];
    }
  }
 
  const auto& trackerParam = TrackerParam::Instance();
  const TMatrixD& paramCov = param.getCovariances();
  double z = param.getZ();
 
  // check if non bending impact parameter is within tolerances
  double nonBendingImpactParam = TMath::Abs(param.getNonBendingCoor() - z * param.getNonBendingSlope());
  double nonBendingImpactParamErr = TMath::Sqrt(paramCov(0, 0) + z * z * paramCov(1, 1) - 2. * z * paramCov(0, 1) + impactMCS2);
  if ((nonBendingImpactParam - trackerParam.sigmaCutForTracking * nonBendingImpactParamErr) > (3. * trackerParam.nonBendingVertexDispersion)) {
    return false;
  }
 
  if (TrackExtrap::isFieldON()) { // depending whether the field is ON or OFF
    // check if bending momentum is within tolerances
    double bendingMomentum = TMath::Abs(1. / param.getInverseBendingMomentum());
    double bendingMomentumErr = TMath::Sqrt(paramCov(4, 4)) * bendingMomentum * bendingMomentum;
    if (chamber < 6 && (bendingMomentum + trackerParam.sigmaCutForTracking * bendingMomentumErr) < SMinBendingMomentum) {
      return false;
    } else if ((bendingMomentum + 3. * bendingMomentumErr) < SMinBendingMomentum) {
      return false;
    }
  } else {
    // or check if bending impact parameter is within tolerances
    double bendingImpactParam = TMath::Abs(param.getBendingCoor() - z * param.getBendingSlope());
    double bendingImpactParamErr = TMath::Sqrt(paramCov(2, 2) + z * z * paramCov(3, 3) - 2. * z * paramCov(2, 3) + impactMCS2);
    if ((bendingImpactParam - trackerParam.sigmaCutForTracking * bendingImpactParamErr) > (3. * trackerParam.bendingVertexDispersion)) {
      return false;
    }
  }
 
  return true;
}
 
//_________________________________________________________________________________________________
void TrackFinder::prepareForwardTracking(std::list<Track>::iterator& itTrack, bool runSmoother)
{
 
  if (runSmoother) {
    auto itStartingParam = std::prev(itTrack->rend());
    mTrackFitter.fit(*itTrack, true, false, &itStartingParam, true);
  }
 
  setCurrentParam(*itTrack, itTrack->last(), itTrack->last().getClusterPtr()->getChamberId(), runSmoother);
}
 
//_________________________________________________________________________________________________
void TrackFinder::prepareBackwardTracking(std::list<Track>::iterator& itTrack, bool refit)
{
 
  if (refit) {
    mTrackFitter.fit(*itTrack, false);
  }
 
  setCurrentParam(*itTrack, itTrack->first(), itTrack->first().getClusterPtr()->getChamberId());
}
 
//_________________________________________________________________________________________________
void TrackFinder::setCurrentParam(Track& track, const TrackParam& param, int chamber, bool smoothed)
{
 
  track.setCurrentParam(param, chamber);
 
  if (smoothed) {
    TrackParam& currentParam = track.getCurrentParam();
    currentParam.setParameters(param.getSmoothParameters());
    currentParam.setCovariances(param.getSmoothCovariances());
  }
 
  if (param.getClusterPtr()) {
    TrackParam& currentParam = track.getCurrentParam();
    currentParam.resetPropagator();
    TrackExtrap::addMCSEffect(currentParam, SChamberThicknessInX0[chamber], -1.);
  }
}
 
//_________________________________________________________________________________________________
bool TrackFinder::propagateCurrentParam(Track& track, int chamber)
{
 
  TrackParam& currentParam = track.getCurrentParam();
  int& currentChamber = track.getCurrentChamber();
  while (currentChamber != chamber) {
 
    currentChamber += (chamber < currentChamber) ? -1 : 1;
 
    if (!TrackExtrap::extrapToZCov(currentParam, SDefaultChamberZ[currentChamber], true)) {
      track.invalidateCurrentParam();
      return false;
    }
 
    TrackExtrap::addMCSEffect(currentParam, SChamberThicknessInX0[currentChamber], -1.);
  }
 
  return true;
}
 
//_________________________________________________________________________________________________
bool TrackFinder::areUsed(const Cluster& cl1, const Cluster& cl2, const std::vector<std::array<uint32_t, 4>>& usedClusters)
{
  int iCl1 = cl1.getDEId() % 2;
  int iCl2 = 2 + cl2.getDEId() % 2;
  for (const auto& clusters : usedClusters) {
    if (clusters[iCl1] == cl1.uid && clusters[iCl2] == cl2.uid) {
      return true;
    }
  }
  return false;
}
 
//_________________________________________________________________________________________________
void TrackFinder::excludeClustersFromIdenticalTracks(const std::array<uint32_t, 4>& currentClusters,
                                                     const std::vector<std::array<uint32_t, 8>>& usedClusters,
                                                     std::unordered_map<int, std::unordered_set<uint32_t>>& excludedClusters)
{
 
  for (const auto& clusters : usedClusters) {
 
    bool identicalTrack(true);
    for (int iCl = 0; iCl < 4; ++iCl) {
      if (clusters[iCl] != currentClusters[iCl] && currentClusters[iCl] > 0) {
        identicalTrack = false;
        break;
      }
    }
 
    if (identicalTrack) {
      for (int iCl = 4; iCl < 8; ++iCl) {
        if (clusters[iCl] > 0) {
          excludedClusters[Cluster::getDEId(clusters[iCl])].emplace(clusters[iCl]);
        }
      }
    }
  }
}
 
//_________________________________________________________________________________________________
void TrackFinder::moveClusters(std::unordered_map<int, std::unordered_set<uint32_t>>& source, std::unordered_map<int, std::unordered_set<uint32_t>>& destination)
{
  for (auto& sourceDE : source) {
    destination[sourceDE.first].insert(sourceDE.second.begin(), sourceDE.second.end());
  }
  source.clear();
}
 
//_________________________________________________________________________________________________
bool TrackFinder::isCompatible(const TrackParam& param, const Cluster& cluster, TrackParam& paramAtCluster)
{
 
  // fast try to add the current cluster
  if (!tryOneClusterFast(param, cluster)) {
    return false;
  }
 
  // try to add the current cluster accurately
  if (tryOneCluster(param, cluster, paramAtCluster) >= mMaxChi2ForTracking) {
    return false;
  }
 
  // save the extrapolated parameters and covariances for the smoother
  paramAtCluster.setExtrapParameters(paramAtCluster.getParameters());
  paramAtCluster.setExtrapCovariances(paramAtCluster.getCovariances());
 
  // compute the new track parameters including the cluster using the Kalman filter
  try {
    mTrackFitter.runKalmanFilter(paramAtCluster);
  } catch (exception const&) {
    return false;
  }
 
  return true;
}
 
//_________________________________________________________________________________________________
bool TrackFinder::tryOneClusterFast(const TrackParam& param, const Cluster& cluster)
{
 
  ++mNCallTryOneClusterFast;
 
  double dZ = cluster.getZ() - param.getZ();
  double dX = cluster.getX() - (param.getNonBendingCoor() + param.getNonBendingSlope() * dZ);
  double dY = cluster.getY() - (param.getBendingCoor() + param.getBendingSlope() * dZ);
  const TMatrixD& paramCov = param.getCovariances();
  double errX2 = paramCov(0, 0) + dZ * dZ * paramCov(1, 1) + 2. * dZ * paramCov(0, 1) + mChamberResolutionX2;
  double errY2 = paramCov(2, 2) + dZ * dZ * paramCov(3, 3) + 2. * dZ * paramCov(2, 3) + mChamberResolutionY2;
 
  double dXmax = TrackerParam::Instance().sigmaCutForTracking * TMath::Sqrt(2. * errX2) + SMaxNonBendingDistanceToTrack;
  double dYmax = TrackerParam::Instance().sigmaCutForTracking * TMath::Sqrt(2. * errY2) + SMaxBendingDistanceToTrack;
 
  if (TMath::Abs(dX) > dXmax || TMath::Abs(dY) > dYmax) {
    return false;
  }
  return true;
}
 
//_________________________________________________________________________________________________
double TrackFinder::tryOneCluster(const TrackParam& param, const Cluster& cluster, TrackParam& paramAtCluster)
{
 
  ++mNCallTryOneCluster;
 
  // Extrapolate the track parameters and covariances at the z position of the cluster
  paramAtCluster = param;
  paramAtCluster.setClusterPtr(&cluster);
  if (!TrackExtrap::extrapToZCov(paramAtCluster, cluster.getZ(), true)) {
    return mTrackFitter.getMaxChi2();
  }
 
  // Compute the cluster-track residuals in bending and non bending directions
  double dX = cluster.getX() - paramAtCluster.getNonBendingCoor();
  double dY = cluster.getY() - paramAtCluster.getBendingCoor();
 
  // Combine the cluster and track resolutions and covariances
  const TMatrixD& paramCov = paramAtCluster.getCovariances();
  double sigmaX2 = paramCov(0, 0) + mChamberResolutionX2;
  double sigmaY2 = paramCov(2, 2) + mChamberResolutionY2;
  double covXY = paramCov(0, 2);
  double det = sigmaX2 * sigmaY2 - covXY * covXY;
 
  // Compute and return the matching chi2
  if (det == 0.) {
    return mTrackFitter.getMaxChi2();
  }
  return (dX * dX * sigmaY2 + dY * dY * sigmaX2 - 2. * dX * dY * covXY) / det;
}
 
//_________________________________________________________________________________________________
uint8_t TrackFinder::requestedStationMask() const
{
  uint8_t mask(0);
  for (int i = 0; i < 5; ++i) {
    if (TrackerParam::Instance().requestStation[i]) {
      mask |= (1 << i);
    }
  }
  return mask;
}
 
//_________________________________________________________________________________________________
int TrackFinder::getTrackIndex(const std::list<Track>::iterator& itCurrentTrack) const
{
 
  if (mDebugLevel < 1) {
    return -3;
  }
 
  if (itCurrentTrack == mTracks.end()) {
    return -1;
  }
 
  int index(0);
  for (auto itTrack = mTracks.begin(); itTrack != itCurrentTrack; ++itTrack) {
    if (itTrack == mTracks.end()) {
      return -2;
    }
    ++index;
  }
 
  return index;
}
 
//_________________________________________________________________________________________________
void TrackFinder::printTracks() const
{
  if (mDebugLevel > 1) {
    for (const auto& track : mTracks) {
      track.print();
    }
  }
}
 
//_________________________________________________________________________________________________
void TrackFinder::printTrack(const Track& track) const
{
  if (mDebugLevel > 1) {
    track.print();
  }
}
 
//_________________________________________________________________________________________________
void TrackFinder::printTrackParam(const TrackParam& trackParam) const
{
  if (mDebugLevel > 1) {
    trackParam.print();
  }
}
 
//_________________________________________________________________________________________________
template <class... Args>
void TrackFinder::print(Args... args) const
{
  if (mDebugLevel > 0) {
    (cout << ... << args) << "\n";
  }
}
 
//_________________________________________________________________________________________________
void TrackFinder::printStats() const
{
  LOG(info) << "number of candidates tracked = " << mNCandidates;
  TrackExtrap::printNCalls();
  LOG(info) << "number of times tryOneClusterFast() is called = " << mNCallTryOneClusterFast;
  LOG(info) << "number of times tryOneCluster() is called = " << mNCallTryOneCluster;
}
 
//_________________________________________________________________________________________________
void TrackFinder::printTimers() const
{
  LOG(info) << "findTrackCandidates duration = " << mTimeFindCandidates.count() << " s";
  LOG(info) << "findMoreTrackCandidates duration = " << mTimeFindMoreCandidates.count() << " s";
  LOG(info) << "followTracks duration = " << mTimeFollowTracks.count() << " s";
  LOG(info) << "improveTracks duration = " << mTimeImproveTracks.count() << " s";
  LOG(info) << "removeConnectedTracks duration = " << mTimeCleanTracks.count() << " s";
  LOG(info) << "refineTracks duration = " << mTimeRefineTracks.count() << " s";
}
 
} // namespace mch
} // namespace o2