Skip to content
Snippets Groups Projects
PrForwardTracking.cpp 124 KiB
Newer Older
Andre Gunther's avatar
Andre Gunther committed
/*****************************************************************************\
* (c) Copyright 2021 CERN for the benefit of the LHCb Collaboration           *
*                                                                             *
* This software is distributed under the terms of the GNU General Public      *
* Licence version 3 (GPL Version 3), copied verbatim in the file "COPYING".   *
*                                                                             *
* In applying this licence, 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.                                       *
\*****************************************************************************/
// standard
#include <algorithm>
#include <array>
#include <cmath>
#include <functional>
#include <memory>
#include <numeric>
#include <string_view>
#include <utility>
#include <vector>

// boost
#include "boost/container/static_vector.hpp"

// vdt
#include "vdt/log.h"

#include "DetDesc/GenericConditionAccessorHolder.h"
#include "Event/PrLongTracks.h"
#include "Event/PrSciFiHits.h"
#include "Event/StateParameters.h"
Olli Lupton's avatar
Olli Lupton committed
#include "Event/Track.h"
#include "FTDet/DeFTDetector.h"
#include "Kernel/CountIterator.h"
#include "LHCbMath/SIMDWrapper.h"
#include "Magnet/DeMagnet.h"
#include "PrKernel/FTGeometryCache.h"
#include "PrKernel/IPrAddUTHitsTool.h"
#include "PrKernel/IPrDebugTrackingTool.h"
#include "PrTrackModel.h"
#include "weights/TMVA_MLP_GhostNN_PrForwardTracking.h"
#include "weights/TMVA_MLP_GhostNN_PrForwardTrackingVelo.h"
Andre Gunther's avatar
Andre Gunther committed
  ++++++++++++++++++++++++++++++++ Forward Tracking +++++++++++++++++++++++++++++++
*/

/** +++++++++++++++++++++ Short Introduction in geometry ++++++++++++++++++++++++++:
 *
 * The SciFi Tracker Detector, or simple Fibre Tracker (FT) consits out of 3 stations.
 * Each station consists out of 4 planes/layers. Thus there are in total 12 layers,
 * in which a particle can leave a hit. The reasonable maximum number of hits a track
 * can have is thus also 12.
 *
 * Each layer consists out of several Fibre mats. A fibre has a diameter of below a mm.
 * Several fibres are glued alongside each other to form a mat.
 * A Scintilating Fibre produces light, if a charged particle traverses. This light is then
 * detected on the outside of the Fibre mat.
 *
 * Looking from the collision point, one (X-)layer looks like the following:
 *
 *                    y       6m
 *                    ^  ||||||||||||| Upper side
 *                    |  ||||||||||||| 2.5m
 *                    |  |||||||||||||
 *                   -|--||||||o||||||----> -x
 *                       |||||||||||||
 *                       ||||||||||||| Lower side
 *                       ||||||||||||| 2.5m
 *
 * All fibres are aranged parallel to the y-axis. There are three different
 * kinds of layers, denoted by X,U,V. The U/V layers are rotated with respect to
 * the X-layers by +/- 5 degrees, to also get a handle of the y position of the
 * particle. As due to the magnetic field particles are deflected in
 * x-direction (only small deflection in y-direction).
 * The layer structure in the FT is XUVX-XUVX-XUVX.
 *
 * The detector is divided into an upeer and a lower side (>/< y=0). As particles
 * are deflected in x direction there are only very(!) few particles that go
 * from the lower to the upper side, or vice versa. The reconstruction algorithm therefore
 * treats upper and lower side tracks independently such that only hits in the respective
 * detector half are considered.
 * Only for U/V layers this can be different because in these layers the
 * complete fibre modules are rotated, producing a zic-zac pattern at y=0, also
 * called  "the triangles". In these layers hits are also searched for on the "other side",
 * if the track is close to
 * y=0. Sketch (rotation exagerated):
 *                                          _.*
 *     y ^   _.*                         _.*
 *       | .*._      Upper side       _.*._
 *       |     *._                 _.*     *._
 *       |--------*._           _.*           *._----------------> x
 *       |           *._     _.*                 *._     _.*
 *                      *._.*       Lower side      *._.*
 *
 *
 *
 *
 *
 *       Zone ordering defined in Detector/FT/FTConstants.h
 *
 *     y ^
 *       |    1  3  5  7     9 11 13 15    17 19 21 23
 *       |    |  |  |  |     |  |  |  |     |  |  |  |
 *       |    x  u  v  x     x  u  v  x     x  u  v  x   <-- type of layer
 *       |    |  |  |  |     |  |  |  |     |  |  |  |
 *       |------------------------------------------------> z
 *       |    |  |  |  |     |  |  |  |     |  |  |  |
 *       |    |  |  |  |     |  |  |  |     |  |  |  |
 *       |    0  2  4  6     8 10 12 14    16 18 20 22
 *
Andre Gunther's avatar
Andre Gunther committed
 *
 *
 *  A detailed introduction to the basic ideas of the Forward tracking can be
 *   found here:
 *   (2002) http://cds.cern.ch/record/684710/files/lhcb-2002-008.pdf
 *   (2007) http://cds.cern.ch/record/1033584/files/lhcb-2007-015.pdf
 *  The most recent note contains (now slightly outdated) information about parameterisations
 *   (2014) http://cds.cern.ch/record/1641927/files/LHCb-PUB-2014-001.pdf
 *
 *
 * Method overview
 *
 * The track reconstruction is done in several steps:
 *
 * 1) Hough-like transformation of SciFi hits using input tracks
 *    -> projectXHitsToHoughSpace and projectStereoHitsToHoughSpace
 *       filling a histogram-like data structure
 * 2) Threshold scan over "histogram" selecting local accumulations of hits
 *    -> pickUpCandidateBins and sortAndCopyBinContents
 * 3) Detailed selection of sets of x hits matching the input track
 *    -> selextXCandidates
 * 4) Completion of candidates by finding best matching u/v hits
 *    -> selectFullCandidates
 *        also assigning a "quality" to each candidate indicating a ghost
 *        probability
 * 5) Removing duplicates among all found tracks
 *    -> removeDuplicates
 */

/** @class PrForwardTracking PrForwardTracking.cpp
 *
 *  @author Olivier Callot
 *  @date   2012-03-20
 *  @author Thomas Nikodem
 *  @date   2013-03-15
 *  @author Michel De Cian
 *  @date   2014-03-12 Changes with make code more standard and faster
 *  @author Sevda Esen
 *  @date   2015-02-13 additional search in the triangles by Marian Stahl
 *  @author Thomas Nikodem
 *  @date   2016-03-09 complete restructuring
 *  @author Olli Lupton
 *  @date   2018-11-07 Imported PrForwardTool into PrForwardTracking as a step towards making PrForwardTracking accept
 *                     selections
 *  @author André Günther
 *  @date 2019-12-03 adapted PrForwardTracking to new SoA hit class PrSciFiHits
 *  @date 2021-03-07 complete redesign
namespace LHCb::Pr::Forward {
  namespace Hough {
    constexpr int DEPTH   = 2;
    constexpr int MAXCAND = 2;
    constexpr int NBINS   = 15;
  } // namespace Hough
  namespace {
    using namespace FT;
    using TracksTag       = Long::Tag;
    using TracksFT        = Long::Tracks;
    using TracksUT        = Upstream::Tracks;
    using TracksVP        = Velo::Tracks;
    using scalar          = SIMDWrapper::scalar::types;
    using simd            = SIMDWrapper::best::types;
    using PrForwardTracks = std::vector<PrForwardTrack>;
Louis Henry's avatar
Louis Henry committed
        ::Hough::HoughSearch<int, Hough::DEPTH, Hough::MAXCAND, Detector::FT::nUVLayersTotal, Hough::NBINS>;
    // names of variables used in multilayer perceptron used to identify ghosts
    constexpr std::array<std::string_view, 9> NNVars{"redChi2",
                                                     "abs((x+(zMagMatch-770.0)*tx)-(xEndT+(zMagMatch-9410.0)*txEndT))",
                                                     "abs(ySeedMatch-yEndT)",
                                                     "abs(yParam0Final-yParam0Init)",
                                                     "abs(yParam1Final-yParam1Init)",
                                                     "abs(ty)",
                                                     "abs(qop)",
                                                     "abs(tx)",
                                                     "abs(xParam1Final-xParam1Init)"};

    constexpr std::array<std::string_view, 10> NNVarsVeloUT{
        "log(abs(1./qop-1./qopUT))",
        "redChi2",
        "abs((x+(zMagMatch-770.0)*tx)-(xEndT+(zMagMatch-9410.0)*txEndT))",
        "abs(ySeedMatch-yEndT)",
        "abs(yParam0Final-yParam0Init)",
        "abs(yParam1Final-yParam1Init)",
        "abs(ty)",
        "abs(qop)",
        "abs(tx)",
        "abs(xParam1Final-xParam1Init)"};
    /**
     * @return larger n such that it is divisible by simd::size and padded by one more simd::size
     */
    constexpr size_t alignAndPad( size_t n ) { return ( n / simd::size + 2 ) * simd::size; }

    /**
     * @class PrParametersX PrForwardTracking.cpp
     * @brief Bundles parameters used selection of x candidates
     */
    struct PrParametersX {
      unsigned minXHits{};
      float    maxXWindow{};
      float    maxXWindowSlope{};
      float    maxXGap{};
      unsigned minStereoHits{};
      float    maxChi2PerDoF{};
      float    maxChi2XProjection{};
      float    maxChi2LinearFit{};
    };
    /**
     * @class PrParametersY PrForwardTracking.cpp
     * @brief Bundles parameters used selection of stereo and full candidates
     */
    struct PrParametersY {
      float    maxTolY{};
      float    uvSearchBinWidth{};
      float    tolY{};
      float    tolYSlope{};
      float    yTolUVSearch{};
      float    tolYTriangleSearch{};
      unsigned minStereoHits{};
      float    maxChi2StereoLinear{};
      float    maxChi2Stereo{};
      float    maxChi2StereoAdd{};
      float    tolYMag{};
      float    tolYMagSlope{};
    };
    template <typename TrackType>
    auto make_TracksFT_from_ancestors( const TrackType& input_tracks ) {
      const TracksVP*                  velo_ancestors{nullptr};
      const TracksUT*                  upstream_ancestors{nullptr};
      const LHCb::Pr::Seeding::Tracks* seed_ancestors{nullptr};
      const auto                       history = Event::Enum::Track::History::PrForward;
      if constexpr ( std::is_same_v<TrackType, TracksUT> ) {
        velo_ancestors     = input_tracks.getVeloAncestors();
        upstream_ancestors = &input_tracks;
      } else {
        velo_ancestors = &input_tracks;
      return TracksFT{velo_ancestors,
                      upstream_ancestors,
                      seed_ancestors,
                      history,
                      Zipping::generateZipIdentifier(),
                      {input_tracks.get_allocator().resource()}};
    /**
     * @brief remove hits from range if found on track
     * @param idx1 first index of range
     * @param idxEnd end index of range
     * @param track track containing hits that are removed from range
     * @details Note that idxEnd is modified by this function because hits are removed by shuffling them to the
     * end of the range and simply moving the end index such that they are not contained anymore.
     * @return true if range is not empty after removing hits
     */
    auto removeUsedHits( int idx1, int& idxEnd, const XCandidate& track, ModSciFiHits::ModPrSciFiHitsSOA& allXHits ) {
      const auto& hitsOnTrack = track.getCoordsToFit();
      const auto  hits_simd   = allXHits.simd();
      for ( auto idx{idx1}; idx < idxEnd; idx += simd::size ) {
        const auto loopMask = simd::loop_mask( idx, idxEnd );
        auto       keepMask = loopMask;

        const auto hit_proxy = hits_simd[idx];
        const auto fulldexes = hit_proxy.get<ModSciFiHits::HitTag::fulldex>();
        const auto coords    = hit_proxy.get<ModSciFiHits::HitTag::coord>();

        for ( auto iHit : hitsOnTrack ) {
          const auto foundMask = fulldexes == iHit;
          keepMask             = keepMask && !foundMask;
        }

        allXHits.store<ModSciFiHits::HitTag::fulldex>( idx1, fulldexes, keepMask );
        allXHits.store<ModSciFiHits::HitTag::coord>( idx1, coords, keepMask );
        idx1 += popcount( keepMask );

        if ( any( !loopMask ) ) {
          const auto coords_new    = select( loopMask, hit_proxy.get<ModSciFiHits::HitTag::coord>(), coords );
          const auto fulldexes_new = select( loopMask, hit_proxy.get<ModSciFiHits::HitTag::fulldex>(), fulldexes );
          allXHits.store<ModSciFiHits::HitTag::fulldex>( idx, fulldexes_new );
          allXHits.store<ModSciFiHits::HitTag::coord>( idx, coords_new );
        }
      }
      const auto ok = idx1 != idxEnd;
      idxEnd        = idx1;
      return ok;
Olli Lupton's avatar
Olli Lupton committed
    }
    /**
     * @brief initialises parameters of x projection parameterisation
     */
    void initXFitParameters( XCandidate& trackCandidate, const VeloSeedExtended& veloSeed ) {
      const auto xAtRef           = trackCandidate.xAtRef();
      auto       dSlope           = ( xAtRef - veloSeed.xStraightAtRef ) / ( zReference - veloSeed.zMag );
      const auto zMag             = veloSeed.calcZMag( dSlope );
      const auto xMag             = veloSeed.seed.x( zMag );
      const auto slopeT           = ( xAtRef - xMag ) / ( zReference - zMag );
      dSlope                      = slopeT - veloSeed.seed.tx;
      const auto CX               = veloSeed.calcCX( dSlope );
      const auto DX               = veloSeed.calcDX( dSlope );
      trackCandidate.getXParams() = {{xAtRef, slopeT, CX, DX}};
     * @brief initialises parameters of y projection parameterisation
    void initYFitParameters( XCandidate& trackCandidate, const VeloSeedExtended& veloSeed ) {
      const auto dSlope           = trackCandidate.getXParams()[1] - veloSeed.seed.tx;
      const auto AY               = veloSeed.yStraightAtRef + veloSeed.calcYCorr( dSlope );
      const auto BY               = veloSeed.seed.ty + veloSeed.calcTyCorr( dSlope );
      const auto CY               = veloSeed.calcCY( dSlope );
      trackCandidate.getYParams() = {{AY, BY, CY}};
     * @class HoughTransformation PrForwardTracking.cpp
     * @brief Contains data and methods used by the Hough-transformation-like step in the Forward Tracking
     * Main methods are defined outside of class
     * Short methods are defined inline
    class HoughTransformation {
    public:
      static constexpr int      minBinOffset{2};
      static constexpr int      nBins{1150 + minBinOffset};
      static constexpr int      reservedBinContent{16};
      static constexpr float    rangeMax{3000.f * Gaudi::Units::mm};
      static constexpr float    rangeMin{-3000.f * Gaudi::Units::mm};
      static constexpr unsigned nDiffPlanesBits{8};
      static constexpr unsigned diffPlanesFlags{0xFF}; // 1111 1111
      static constexpr unsigned uvFlags{0x666};        // 0110 0110 0110
      static constexpr unsigned xFlags{0x999};         // 1001 1001 1001
      static_assert( sizeof( unsigned ) >= 4 );        // at the moment at least 32 bit datatype for plane encoding
      static_assert( ( reservedBinContent % simd::size == 0 ) &&
                     "Currently only multiples of avx2 vector width supported" );

      void projectStereoHitsToHoughSpace( const VeloSeedExtended&, const ZoneCache&, std::pair<float, float>,
                                          const FT::Hits& );
      void projectXHitsToHoughSpace( const VeloSeedExtended&, std::pair<float, float>, const FT::Hits& );
      auto pickUpCandidateBins( int, int, int, int );
      void sortAndCopyBinContents( int, ModSciFiHits::ModPrSciFiHitsSOA& );
      /**
       * @brief clears containers for size and therefore also content, plane counter and candidates
       */
      void clear() {
        m_candSize = 0;
        m_planeCounter.fill( 0 );
        m_binContentSize.fill( 0 );
      }
      /**
       * @brief Calculates bin numbers for x values in "Hough space".
       * @param xAtRef x values in "Hough space", i.e. on reference plane.
       * @param idx index of x hit
       * @param maxIdx last index/end of range
       * @return bin numbers
       * @details The binning is non-simple and follows the occupancy in the detector. Furthermore, only values within
       * a certain range are considered. x values outside of this range will cause a return of bin number 0, the
       * "garbage bin". The binning function is a fast sigmoid. The parameterisation is obtained using the
       * Reco-Parametersiation-Tuner.
       * https://gitlab.cern.ch/gunther/prforwardtracking-parametrisation-tuner/-/tree/master
       */
      auto calculateBinNumber( simd::float_v xAtRef, int idx, int maxIdx ) const {
        const auto boundaryMask = xAtRef < rangeMax && xAtRef > rangeMin && simd::loop_mask( idx, maxIdx );
        // p[0] + p[1] * x / (1 + p[2] * abs(x)) for nBins = 1152
        constexpr auto p                 = std::array{576.9713937732083f, 0.5780266207743059f, 0.0006728484590464921f};
        const auto     floatingBinNumber = p[0] + p[1] * xAtRef / ( 1.f + p[2] * abs( xAtRef ) );
        // it can happen that two almost equal xAtRef lead to floating bin numbers right on a bin edge, where one
        // falls into the left, the other into the right bin due to numerical differences coming from the calculations
        // above. This messes with the sorting that is only done per bin later on, but has no dramatic consequences
        // (the hits are so close to each other in that case that no deep logic breaks)
        const auto binNumber = static_cast<simd::int_v>( floatingBinNumber );
        return select( boundaryMask, binNumber, 0 );
      }
      /**
       * @brief Performs a compressed store of candidates and handles size incrementation.
       * @param mask Mask of candidates to store.
       * @param bins The bin number(s) qualifying as candidates.
       */
      void compressstoreCandidates( simd::mask_v mask, simd::int_v bins ) {
        bins.compressstore( mask, m_candidateBins.data() + m_candSize );
        m_candSize += popcount( mask );
Marco Clemencic's avatar
Marco Clemencic committed

      /**
       * @brief Decodes number of different planes encoded in number stored in plane counting array.
       * @param bin The bin number
       */
      template <typename I>
      auto decodeNbDifferentPlanes( int bin ) const {
        return I{m_planeCounter.data() + bin} & diffPlanesFlags;
      }

      /**
       * @brief return integer with bits encoding hit planes shifted to the beginning
       * @param bin The bin number
       */
      template <typename I>
      auto planeBits( int bin ) const {
        return I{m_planeCounter.data() + bin} >> nDiffPlanesBits;
      /**
       * @brief Make a span for one bin
       * @param iBin bin number
       * @return span for the coordinates in this bin
       */
      auto getBinContentCoord( int iBin ) {
        return LHCb::span{m_binContentCoord}.subspan( iBin * reservedBinContent, m_binContentSize[iBin] );
      }
      auto getBinContentCoord( int iBin ) const {
        return LHCb::span{m_binContentCoord}.subspan( iBin * reservedBinContent, m_binContentSize[iBin] );
      }

      /**
       * @brief Make a span for one bin
       * @param iBin bin number
       * @return span for the PrSciFiHits indices in this bin
       */
      auto getBinContentFulldex( int iBin ) {
        return LHCb::span{m_binContentFulldex}.subspan( iBin * reservedBinContent, m_binContentSize[iBin] );
      }
      auto getBinContentFulldex( int iBin ) const {
        return LHCb::span{m_binContentFulldex}.subspan( iBin * reservedBinContent, m_binContentSize[iBin] );
      }

      /**
       * @brief Sorts and then copies content of a bin to a new container
       * @param bin bin number
       * @param allXHits container to copy to (SoA)
       * @details The sorting is done by an insertion sort because the number of elements to sort is small by design.
       * The sorting is done by storing the permutations and then using a gather operation to shuffle according to the
       * permutations, followed by the storing.
       */
      [[gnu::always_inline]] void sortAndCopyBin( int bin, ModSciFiHits::ModPrSciFiHitsSOA& allXHits ) {
        if ( const auto binSize = m_binContentSize[bin]; binSize ) {
          std::iota( m_binPermutations.begin(), m_binPermutations.end(), 0 );
          const auto coordContent = getBinContentCoord( bin );
          for ( auto iCoord{1}; iCoord < binSize; ++iCoord ) {
            const auto insertVal = coordContent[iCoord];
            auto       insertPos = iCoord;
            for ( auto movePos = iCoord; movePos && insertVal < coordContent[m_binPermutations[--movePos]];
                  --insertPos ) {
              m_binPermutations[insertPos] = m_binPermutations[movePos];
            }
            m_binPermutations[insertPos] = iCoord;
          }
          const auto fulldexContent = getBinContentFulldex( bin );
          const auto max_size       = allXHits.size() + binSize;
          auto       i{0};
          do {
            const auto permut_v  = simd::int_v{m_binPermutations.data() + i};
            const auto fulldex_v = gather( fulldexContent.data(), permut_v );
            const auto coord_v   = gather( coordContent.data(), permut_v );
            auto       hit       = allXHits.template emplace_back<>( max_size );
            hit.template field<ModSciFiHits::HitTag::fulldex>().set( fulldex_v );
            hit.template field<ModSciFiHits::HitTag::coord>().set( coord_v );
            i += simd::size;
          } while ( i < binSize );
        }
      }

      // for debugging
      [[maybe_unused]] friend std::ostream& operator<<( std::ostream& os, const HoughTransformation& h ) {
        const auto candidateBins = LHCb::span( h.m_candidateBins.begin(), h.m_candSize );
        os << "Found " << h.m_candSize << " candidate bins:" << std::endl;
        for ( auto bin : candidateBins ) {
          os << "Bin " << bin << " x coords: [ ";
          const auto coords    = h.getBinContentCoord( bin );
          const auto fulldexes = h.getBinContentFulldex( bin );
          for ( auto i{0}; i < h.m_binContentSize[bin]; ++i ) {
            os << "(" << fulldexes[i] << ", " << coords[i] << ")";
            os << ( i != h.m_binContentSize[bin] - 1 ? " | " : "" );
          }
          os << "]"
             << " size=" << h.m_binContentSize[bin] << std::endl;
          const auto planes =
              ( h.m_planeCounter[bin] | h.m_planeCounter[bin - 1] | h.m_planeCounter[bin + 1] ) >> h.nDiffPlanesBits;
          os << "Used planes including neighbours: " << std::bitset<12>( planes ) << std::endl;
          os << std::endl;
        }
        return os;
      }

    private:
      /**
       * @property Contain the coordinates (x hits on reference plane) and PrSciFiHits indices in each bin
       */
      alignas( 64 ) std::array<float, alignAndPad( nBins ) * reservedBinContent> m_binContentCoord;
      alignas( 64 ) std::array<int, alignAndPad( nBins ) * reservedBinContent> m_binContentFulldex;
       * @property Contains the the number of unique planes within a bin encoded together with bit flags for each plane.
      alignas( 64 ) std::array<int, alignAndPad( nBins )> m_planeCounter{};
       * @property Contains the number of hits in each bin
      alignas( 64 ) std::array<int, alignAndPad( nBins )> m_binContentSize{};
      /**
       * @property Contains a selection of bins that qualify as part of a track candidates.
       * @details The number of candidates is tracked by candSize.
       */
      alignas( 64 ) std::array<int, alignAndPad( nBins )> m_candidateBins;
      /**
       * @property sstorage that contains temporary permutations from insertion sort
       */
      alignas( 64 ) std::array<int, reservedBinContent> m_binPermutations;
      int m_candSize{0};
    };
Marco Clemencic's avatar
Marco Clemencic committed

     * @brief separate hits that are alone on their SciFi plane from ones that have friends on the plane
     * @param otherPlanes Container keeping track of the plane numbers containing multiples
     * @param protoCand the candidate under consideration
     * @details single hits are directly put into coordsToFit container inside protoCand
    template <typename Container>
    [[gnu::always_inline]] inline void separateSingleHitsForFit( Container& otherPlanes, XCandidate& protoCand,
                                                                 const VeloSeedExtended& veloSeed ) {
      otherPlanes.clear();
Louis Henry's avatar
Louis Henry committed
      for ( unsigned int iPlane{0}; iPlane < Detector::FT::nXLayersTotal; ++iPlane ) {
        if ( protoCand.planeSize( iPlane ) == 1 ) {
          const auto idx = protoCand.getIdx( iPlane * protoCand.planeMulti );
          protoCand.addHitForLineFit( idx, veloSeed );
        } else if ( protoCand.planeSize( iPlane ) ) {
          otherPlanes.push_back( iPlane );
        }
      }
     * @brief add hits to candidate from planes that have multiple hits
     * @param otherPlanes Container keeping track of the plane numbers containing multiples
     * @param protoCand the candidate under consideration
     * @details only one hit per plane is allowed, selected by smallest chi2 calculated
     * using a straight line fit and putting them into coordsToFit container inside protoCand
    template <typename Container>
    [[gnu::always_inline]] inline void addBestOtherHits( const Container& otherPlanes, XCandidate& protoCand,
                                                         const VeloSeedExtended& veloSeed ) {
      for ( auto iPlane : otherPlanes ) {
        const auto idxSpan = protoCand.getIdxSpan( iPlane );
        assert( idxSpan.size() > 0 );
        const auto bestIdx = [&] {
          auto best{0};
          auto bestChi2{std::numeric_limits<float>::max()};
          for ( auto idx : idxSpan ) {
            const auto chi2 = protoCand.lineChi2( idx, veloSeed );
            if ( chi2 < bestChi2 ) {
              bestChi2 = chi2;
              best     = idx;
            }
          }
          return best;
        }();
        protoCand.addHitForLineFit( bestIdx, veloSeed );
        protoCand.setPlaneSize( iPlane, 1 );
        protoCand.solveLineFit();
      }
     * @brief store all hits present in candidate for fitting
    void prepareAllHitsForFit( const ModSciFiHits::ModPrSciFiHitsSOA& allXHits, XCandidate& protoCand ) {
Louis Henry's avatar
Louis Henry committed
      for ( unsigned int iPlane{0}; iPlane < Detector::FT::nXLayersTotal; ++iPlane ) {
        const auto idxSpan = protoCand.getIdxSpan( iPlane );
        std::transform( idxSpan.begin(), idxSpan.end(), std::back_inserter( protoCand.getCoordsToFit() ),
                        [&]( auto idx ) {
                          protoCand.xAtRef() += allXHits.coord( idx );
                          return allXHits.fulldex( idx );
                        } );
      }
      protoCand.xAtRef() /= static_cast<float>( protoCand.getCoordsToFit().size() );
     * @brief fit linear parameters of parameterisation in x and remove bad hits
     * @param protoCand candidate to fit
     * @param pars bundle of parameters that tune hit removal
     * @param SciFiHits PrSciFiHits container
     * @details although only the linear parameters are fitted, the full parameterisation is used meaning
     * that the quadratic and the cubic term are fixed to their initial values.
    template <bool secondLoop>
    void fitLinearXProjection( XCandidate& protoCand, const PrParametersX& pars, const FT::Hits& SciFiHits,
                               const VeloSeedExtended& veloSeed ) {

      auto& coordsToFit = protoCand.getCoordsToFit();
      bool  fit{true};
      while ( fit ) {
        float s0{0.f}, sz{0.f}, sz2{0.f}, sd{0.f}, sdz{0.f};
        for ( auto iHit : coordsToFit ) {
          auto dz = SciFiHits.z( iHit ) - zReference;
          dz += veloSeed.yStraightInZone[SciFiHits.planeCode( iHit )] * SciFiHits.dzDy( iHit );
          const auto d = SciFiHits.x( iHit ) - protoCand.x( dz );
          const auto w = SciFiHits.w( iHit );
          s0 += w;
          sz += w * dz;
          sz2 += w * dz * dz;
          sd += w * d;
          sdz += w * d * dz;
        }
        const auto den = sz * sz - s0 * sz2;
        if ( den == 0.f ) return;
        const auto da = ( sdz * sz - sd * sz2 ) / den;
        const auto db = ( sd * sz - s0 * sdz ) / den;
        protoCand.addXParams<2>( std::array{da, db} );

        const auto itEnd = coordsToFit.end();
        auto       worst = itEnd;
        auto       maxChi2{0.f};
        const bool notMultiple = protoCand.nDifferentPlanes() == coordsToFit.size();
        for ( auto itH = coordsToFit.begin(); itH != itEnd; ++itH ) {
          auto dz = SciFiHits.z( *itH ) - zReference;
          dz += veloSeed.yStraightInZone[SciFiHits.planeCode( *itH )] * SciFiHits.dzDy( *itH );
          const auto d    = SciFiHits.x( *itH ) - protoCand.x( dz );
          const auto chi2 = d * d * SciFiHits.w( *itH );
          if ( chi2 > maxChi2 && ( notMultiple || protoCand.nInSamePlane( *itH ) > 1 ) ) {
            maxChi2 = chi2;
            worst   = itH;
        fit = false;
        if ( const int ip = SciFiHits.planeCode( *worst ) / 2u;
             maxChi2 > pars.maxChi2LinearFit || protoCand.planeSize( ip ) > 1 ) {
          protoCand.removePlane( ip );
          std::iter_swap( worst, std::prev( itEnd ) );
          coordsToFit.pop_back();
          if ( coordsToFit.size() < pars.minXHits ) return;
          fit = true;
        }
Gitlab CI's avatar
Gitlab CI committed

     * @brief try to find an x hit on a still empty x plane that matches the track
     * @param protoCand track candidate
     * @param pars bundle of parameters that tune the window defining if an x hit matches or not
     * @param veloSeed extended velo track
     * @param SciFiHits PrSciFiHits container
     * @param maxChi2XAdd defines the maximum allowed chi2 deviation from track
     * @details Hits are only added if they are within a window around the extrapolated position on the plane
     * under consideration and do not deviate too much in chi2.
    template <bool secondLoop>
    auto fillEmptyXPlanes( XCandidate& protoCand, const PrParametersX& pars, const VeloSeedExtended& veloSeed,
                           const FT::Hits& SciFiHits, float maxChi2XAdd ) {
Louis Henry's avatar
Louis Henry committed
      if ( protoCand.nDifferentPlanes() == Detector::FT::nXLayersTotal ) return false;

      bool       added{false};
      const auto xAtRef  = protoCand.getXParams()[0];
      const auto xWindow = pars.maxXWindow +
                           ( std::abs( xAtRef ) + std::abs( xAtRef - veloSeed.xStraightAtRef ) ) * pars.maxXWindowSlope;

Louis Henry's avatar
Louis Henry committed
      for ( unsigned int iPlane{0}; iPlane < Detector::FT::nXLayersTotal; ++iPlane ) {
        if ( protoCand.planeSize( iPlane ) ) continue;
        const auto pc             = Detector::FT::xLayers[iPlane];
        const auto [xStart, xEnd] = SciFiHits.getZoneIndices( 2 * pc + veloSeed.upperHalfTrack );

        const auto xPred = protoCand.x( veloSeed.betterZ[pc] - zReference );

        const auto startwin = SciFiHits.get_lower_bound_fast<4>( xStart, xEnd, xPred - xWindow );
        auto       endwin{startwin};
        auto       bestChi2{maxChi2XAdd};
        auto       best{0};
        for ( const auto maxX = xPred + xWindow; SciFiHits.x( endwin ) <= maxX; ++endwin ) {
          const auto d = SciFiHits.x( endwin ) - xPred;
          if ( const auto chi2 = d * d * SciFiHits.w( endwin ); chi2 < bestChi2 ) {
            best     = endwin;
        if ( best ) {
          protoCand.getCoordsToFit().push_back( best );
          ++protoCand.nDifferentPlanes();
          ++protoCand.planeSize( iPlane );
          added = true;
      return added;
    /**
     * @brief try to find an x hit on a still empty x plane that matches the track
     * @param track track candidate
     * @param pars bundle of parameters that tune the window defining if an x hit matches or not
     * @param veloSeed extended velo track
     * @param SciFiHits PrSciFiHits container
     * @param maxChi2XAdd defines the maximum allowed chi2 deviation from track
     * @details Hits are only added if they are within a window around the extrapolated position on the plane
     * under consideration and do not deviate too much in chi2. The empty planes are determined first.
     */
    auto fillEmptyXPlanes( PrForwardTrack& track, const PrParametersX& pars, const VeloSeedExtended& veloSeed,
                           const FT::Hits& SciFiHits, float maxChi2XAdd ) {

Louis Henry's avatar
Louis Henry committed
      if ( track.size() == Detector::FT::nLayersTotal ) return false;
      auto&                                    coordsToFit = track.getCoordsToFit();
      std::bitset<Detector::FT::nXLayersTotal> hasXLayer{};
      unsigned int                             nX{0};
      for ( auto iHit : coordsToFit ) {
        if ( const auto pc = SciFiHits.planeCode( iHit ); Detector::FT::isXLayer[pc] ) {
          hasXLayer.set( pc / 2u );
          ++nX;
Louis Henry's avatar
Louis Henry committed
      if ( nX == Detector::FT::nXLayersTotal ) return false;
      bool       added{false};
      const auto xAtRef  = track.getXParams()[0];
      const auto xWindow = pars.maxXWindow +
                           ( std::abs( xAtRef ) + std::abs( xAtRef - veloSeed.xStraightAtRef ) ) * pars.maxXWindowSlope;

Louis Henry's avatar
Louis Henry committed
      for ( unsigned int iPlane{0}; iPlane < Detector::FT::nXLayersTotal; ++iPlane ) {
        if ( hasXLayer.test( iPlane ) ) continue;
        const auto pc             = Detector::FT::xLayers[iPlane];
        const auto [xStart, xEnd] = SciFiHits.getZoneIndices( 2 * pc + veloSeed.upperHalfTrack );

        const auto xPred = track.x( veloSeed.betterZ[pc] - zReference );

        const auto startwin = SciFiHits.get_lower_bound_fast<4>( xStart, xEnd, xPred - xWindow );
        auto       endwin{startwin};
        auto       bestChi2{maxChi2XAdd};
        auto       best{0};
        for ( const auto maxX = xPred + xWindow; SciFiHits.x( endwin ) <= maxX; ++endwin ) {
          const auto d = SciFiHits.x( endwin ) - xPred;
          if ( const auto chi2 = d * d * SciFiHits.w( endwin ); chi2 < bestChi2 ) {
            bestChi2 = chi2;
            best     = endwin;
          }
        }
        if ( best ) {
          track.getCoordsToFit().push_back( best );
          added = true;
      return added;
Olli Lupton's avatar
Olli Lupton committed

    /**
     * @brief fit up to quadratic parameters of parameterisation in x and remove bad hits
     * @param track candidate to fit
     * @param pars bundle of parameters that tune hit removal
     * @param SciFiHits PrSciFiHits container
     * @details although only the parameters up to the quadratic term are fitted, the full parameterisation is used
     * meaning that the cubic term is fixed to its initial value.
     */
    auto fitXProjection( XCandidate& track, const PrParametersX& pars, const FT::Hits& SciFiHits ) {
      auto&      coordsToFit = track.getCoordsToFit();
      const auto minHits     = pars.minXHits;
      while ( fit ) {
        // plus one because we are fitting all params but one
        const auto nDoF = coordsToFit.size() - track.getXParams().size() + 1;
        if ( nDoF < 1 ) return false;
        auto s0{0.f}, sz{0.f}, sz2{0.f}, sz3{0.f}, sz4{0.f}, sd{0.f}, sdz{0.f}, sdz2{0.f};
        for ( auto iHit : coordsToFit ) {
          const auto dzNoScale = track.calcBetterDz( iHit );
          const auto d         = track.distanceXHit( iHit, dzNoScale );
          const auto w         = SciFiHits.w( iHit );
          const auto dz        = .001f * dzNoScale;
          const auto dz2       = dz * dz;
          const auto wdz       = w * dz;
          s0 += w;
          sz += wdz;
          sz2 += wdz * dz;
          sz3 += wdz * dz2;
          sz4 += wdz * dz2 * dz;
          sd += w * d;
          sdz += wdz * d;
          sdz2 += w * d * dz2;
        const auto b1  = sz * sz - s0 * sz2;
        const auto c1  = sz2 * sz - s0 * sz3;
        const auto d1  = sd * sz - s0 * sdz;
        const auto b2  = sz2 * sz2 - sz * sz3;
        const auto c2  = sz3 * sz2 - sz * sz4;
        const auto d2  = sdz * sz2 - sz * sdz2;
        const auto den = b1 * c2 - b2 * c1;
        if ( den == 0.f ) return false;
        const auto db = ( d1 * c2 - d2 * c1 ) / den;
        const auto dc = ( d2 * b1 - d1 * b2 ) / den;
        const auto da = ( sd - db * sz - dc * sz2 ) / s0;
        track.addXParams<3>( std::array{da, db * 1.e-3f, dc * 1.e-6f} );

        auto maxChi2{0.f};
        auto totChi2{0.f};

        const auto itEnd = coordsToFit.end();
        auto       worst = itEnd;
        for ( auto itH = coordsToFit.begin(); itH != itEnd; ++itH ) {
          auto dz = SciFiHits.z( *itH ) - zReference;
          dz += track.y( dz ) * SciFiHits.dzDy( *itH );
          const auto chi2 = track.chi2XHit( *itH, dz );
          totChi2 += chi2;
          if ( chi2 > maxChi2 ) {
            maxChi2 = chi2;
            worst   = itH;
          }
        }
        track.setChi2NDoF( {totChi2, static_cast<float>( nDoF )} );
        if ( worst == itEnd ) return true;
        fit = false;
        if ( totChi2 > pars.maxChi2PerDoF * nDoF || maxChi2 > pars.maxChi2XProjection ) {
          track.removeFromPlane( *worst );
          std::iter_swap( worst, std::prev( itEnd ) );
          coordsToFit.pop_back();
          if ( coordsToFit.size() < minHits ) return false;
          fit = true;
        }
      return true;
    /**
     * @brief fit up to quadratic parameters of parameterisation in x and remove bad hits
     * @param track candidate to fit
     * @param pars bundle of parameters that tune hit removal
     * @param SciFiHits PrSciFiHits container
     * @details although only the parameters up to the quadratic term are fitted, the full parameterisation is used
     * meaning that the cubic term is fixed to its initial value.
     */
    auto fitXProjection( PrForwardTrack& track, const PrParametersX& pars, const FT::Hits& SciFiHits ) {
      auto&      coordsToFit = track.getCoordsToFit();
      const auto minHits     = pars.minXHits + pars.minStereoHits;

      bool fit{true};
      while ( fit ) {
        // plus one because we are fitting all params but one
        const auto nDoF = coordsToFit.size() - track.getXParams().size() + 1;
        if ( nDoF < 1 ) return false;
        auto s0{0.f}, sz{0.f}, sz2{0.f}, sz3{0.f}, sz4{0.f}, sd{0.f}, sdz{0.f}, sdz2{0.f};
        for ( auto iHit : coordsToFit ) {
          const auto dzNoScale = track.getBetterDz( iHit, SciFiHits.z( iHit ) - zReference, SciFiHits );
          const auto d         = track.distance( iHit, dzNoScale, SciFiHits );
          const auto w         = SciFiHits.w( iHit );
          const auto dz        = .001f * dzNoScale;
          const auto dz2       = dz * dz;
          const auto wdz       = w * dz;
          s0 += w;
          sz += wdz;
          sz2 += wdz * dz;
          sz3 += wdz * dz2;
          sz4 += wdz * dz2 * dz;
          sd += w * d;
          sdz += wdz * d;
          sdz2 += w * d * dz2;
        const auto b1  = sz * sz - s0 * sz2;
        const auto c1  = sz2 * sz - s0 * sz3;
        const auto d1  = sd * sz - s0 * sdz;
        const auto b2  = sz2 * sz2 - sz * sz3;
        const auto c2  = sz3 * sz2 - sz * sz4;
        const auto d2  = sdz * sz2 - sz * sdz2;
        const auto den = b1 * c2 - b2 * c1;
        if ( den == 0.f ) return false;
        const auto db = ( d1 * c2 - d2 * c1 ) / den;
        const auto dc = ( d2 * b1 - d1 * b2 ) / den;
        const auto da = ( sd - db * sz - dc * sz2 ) / s0;
        track.addXParams<3>( std::array{da, db * 1.e-3f, dc * 1.e-6f} );

        auto maxChi2{0.f};
        auto totChi2{0.f};

        const auto itEnd = coordsToFit.end();
        auto       worst = itEnd;
        for ( auto itH = coordsToFit.begin(); itH != itEnd; ++itH ) {
          const auto dz   = track.getBetterDz( *itH, SciFiHits.z( *itH ) - zReference, SciFiHits );
          const auto chi2 = track.chi2( *itH, dz, SciFiHits );
          totChi2 += chi2;
          if ( chi2 > maxChi2 ) {
            maxChi2 = chi2;
            worst   = itH;
          }
        }
        track.setChi2NDoF( {totChi2, static_cast<float>( nDoF )} );
        if ( worst == itEnd ) return true;
        fit = false;
        if ( totChi2 > pars.maxChi2PerDoF * nDoF || maxChi2 > pars.maxChi2XProjection ) {
          std::iter_swap( worst, std::prev( itEnd ) );
          coordsToFit.pop_back();
          if ( coordsToFit.size() < minHits ) return false;
          fit = true;
        }
      return true;
Andre Gunther's avatar
Andre Gunther committed
    }
     * @brief Collect the stereo layer hits that are on the other side because of triangle reaching into side.
     *
     * @param zoneNumberOS zoneNumber of the other side
     * @param xMin minimum coordinate to check on this layer
     * @param xMax maximum coordinate to check on this layer
     * @param dxDySign the sign of stereo rotation of layer
     * @param xPredShifted predicated coordinate corrected for shift due to stereo angle
     * @param yInZone y position of track in layer
     * @param SciFiHits hits in detector
     * @param pars y parameter pack
     * @param hough instance of HoughSerach
    [[gnu::noinline]] void collectTriangleHits( unsigned zoneNumberOS, float xMin, float xMax, float dxDySign,
                                                float xPredShifted, float yInZone, const FT::Hits& SciFiHits,
                                                const PrParametersY& pars, StereoSearch& hough ) {
      const auto [uvStart, uvEnd] = SciFiHits.getZoneIndices( zoneNumberOS );
      const auto iUVStart         = SciFiHits.get_lower_bound_fast<4>( uvStart, uvEnd, xMin );
      const auto yMin             = yInZone + pars.yTolUVSearch;
      const auto yMax             = yInZone - pars.yTolUVSearch;
      auto       triangleOK       = [&]( int index ) {
        const auto [yMinHit, yMaxHit] = SciFiHits.yEnd( index );
        return yMax < yMaxHit && yMin > yMinHit;
      };
      for ( auto iUV{iUVStart}; SciFiHits.x( iUV ) <= xMax; ++iUV ) {
        const auto signedDx = ( SciFiHits.x( iUV ) - xPredShifted ) * dxDySign;
        if ( triangleOK( iUV ) ) { hough.add( zoneNumberOS / 4u, signedDx, iUV ); }
     * @brief Collect stereo hits close to the x projection of the track in all layers.
     *
     * @param search_result this is the container where the candidates from the search are written to
     * @param track the track with mostly x information yet
     * @param SciFiHits hits in the detector
     * @param veloSeed velo track parameters
     * @param cache cache for position of scifi zones
     * @param pars parameter pack for stereo candidate
     * @return gsl::span span over candidates
    auto collectStereoHits( StereoSearch::result_type& search_result, const PrForwardTrack& track,
                            const FT::Hits& SciFiHits, const VeloSeedExtended& veloSeed, const ZoneCache& cache,
                            const PrParametersY& pars ) {

      const auto minTol = -pars.tolY - pars.tolYSlope * ( HoughTransformation::rangeMax + pars.maxTolY );
      // the first parameter controls the threshold which bins in the hough search have to have at least
Louis Henry's avatar
Louis Henry committed
      StereoSearch hough{Detector::FT::nUVLayersTotal / 2, minTol, pars.uvSearchBinWidth};
      direct_debug( "Collect stereo hits with maximal", std::abs( minTol ), "mm deviation from x track." );
      const auto& uvZones = veloSeed.upperHalfTrack ? Detector::FT::uvZonesUpper : Detector::FT::uvZonesLower;
      for ( auto zoneNumber : uvZones ) {
        const auto side =
            track.x( cache.z( zoneNumber / 2u ) - zReference ) > 0.f ? ZoneCache::Side::A : ZoneCache::Side::C;
        const auto zZone = cache.z( zoneNumber, side );
        // TODO: can we improve here by using the ShiftCalculator?
        const auto yInZone      = track.y( zZone - zReference );
        const auto betterZ      = zZone + yInZone * cache.dzdy( zoneNumber, side );
        const auto xPred        = track.x( betterZ - zReference );
        const auto dxDy         = cache.dxdy( zoneNumber, side );
        const auto xPredShifted = xPred - yInZone * dxDy;
        const auto dxTol =
            pars.tolY +
            pars.tolYSlope * ( std::abs( xPred - veloSeed.xStraightInZone[zoneNumber / 2u] ) + std::abs( yInZone ) );
        const auto [uvStart, uvEnd] = SciFiHits.getZoneIndices( zoneNumber );
        const auto xMin             = xPredShifted - dxTol;
        const auto xMax             = xPredShifted + dxTol;
        // the difference only is interepreation: currently: negative dx * sign means underestimated y position and vice
        // versa but it makes a difference for the sorting!
        const auto dxDySign = std::copysign( 1.f, dxDy );
        const auto iUVStart = SciFiHits.get_lower_bound_fast<4>( uvStart, uvEnd, xMin );
        for ( auto iUV{iUVStart}; SciFiHits.x( iUV ) <= xMax; ++iUV ) {
          const auto dz          = SciFiHits.z( iUV ) + yInZone * SciFiHits.dzDy( iUV ) - zReference;
          const auto predShifted = track.x( dz ) - yInZone * SciFiHits.dxDy( iUV );
          const auto signedDx    = ( SciFiHits.x( iUV ) - predShifted ) * dxDySign;
          // we actually only care about 0-5 for layers, so divide by 4
          hough.add( zoneNumber / 4u, signedDx, iUV );
        }
        if ( std::abs( yInZone ) < pars.tolYTriangleSearch ) {
          const auto zoneNumberOS = veloSeed.upperHalfTrack ? zoneNumber - 1 : zoneNumber + 1;
          collectTriangleHits( zoneNumberOS, xMin, xMax, dxDySign, xPredShifted, yInZone, SciFiHits, pars, hough );