Newer
Older
/*****************************************************************************\
* (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. *
\*****************************************************************************/
Andre Gunther
committed
// standard
#include <algorithm>
#include <array>
Andre Gunther
committed
#include <cassert>
Andre Gunther
committed
#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"
// from Gaudi
#include "DetDesc/GenericConditionAccessorHolder.h"
Andre Gunther
committed
#include "GaudiAlg/GaudiTupleTool.h"
#include "LHCbAlgs/Transformer.h"
// from LHCb
#include "Event/PrSciFiHits.h"
#include "Event/StateParameters.h"
#include "Event/StateVector.h"
#include "FTDet/DeFTDetector.h"
#include "Kernel/CountIterator.h"
#include "LHCbMath/SIMDWrapper.h"
// local
#include "HoughSearch.h"
#include "PrKernel/FTGeometryCache.h"
Andre Gunther
committed
#include "PrKernel/IPrDebugTrackingTool.h"
#include "PrTrackModel.h"
// NN for ghostprob
#include "weights/TMVA_MLP_GhostNN_PrForwardTracking.h"
#include "weights/TMVA_MLP_GhostNN_PrForwardTrackingVelo.h"
/**
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
++++++++++++++++++++++++++++++++ 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
*
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
*
*
* 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 Hough {
constexpr int DEPTH = 2;
constexpr int MAXCAND = 2;
constexpr int NBINS = 15;
} // namespace Hough
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>;
using StereoSearch =
::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(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{};
};
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 );
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;
/**
* @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}};
Andre Gunther
committed
* @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 );
/**
* @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;
Andre Gunther
committed
/**
* @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;
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;
* @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 ) {
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 ) {
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;
}
}
if ( worst == itEnd ) return;
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;
}
* @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 ) {
if ( protoCand.nDifferentPlanes() == Detector::FT::nXLayersTotal ) return false;
const auto xAtRef = protoCand.getXParams()[0];
const auto xWindow = pars.maxXWindow +
( std::abs( xAtRef ) + std::abs( xAtRef - veloSeed.xStraightAtRef ) ) * pars.maxXWindowSlope;
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 ) {
Andre Gunther
committed
bestChi2 = chi2;
Andre Gunther
committed
}
}
if ( best ) {
protoCand.getCoordsToFit().push_back( best );
++protoCand.nDifferentPlanes();
++protoCand.planeSize( iPlane );
added = true;
/**
* @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.
*/
template <bool secondLoop>
auto fillEmptyXPlanes( PrForwardTrack& track, const PrParametersX& pars, const VeloSeedExtended& veloSeed,
const FT::Hits& SciFiHits, float maxChi2XAdd ) {
if ( track.size() == Detector::FT::nLayersTotal ) return false;
auto& coordsToFit = track.getCoordsToFit();
std::bitset<Detector::FT::nXLayersTotal> hasXLayer{};
unsigned int nX{0};
if ( const auto pc = SciFiHits.planeCode( iHit ); Detector::FT::isXLayer[pc] ) {
}
}
if ( nX == Detector::FT::nXLayersTotal ) return false;
const auto xAtRef = track.getXParams()[0];
const auto xWindow = pars.maxXWindow +
( std::abs( xAtRef ) + std::abs( xAtRef - veloSeed.xStraightAtRef ) ) * pars.maxXWindowSlope;
for ( unsigned int iPlane{0}; iPlane < Detector::FT::nXLayersTotal; ++iPlane ) {
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;
}
}
/**
* @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 )} );
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;
}
/**
* @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;
877
878
879
880
881
882
883
884
885
886
887
888
889
890
891
892
893
894
895
896
897
898
899
900
901
902
903
904
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 )} );
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;
}
* @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
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;
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 );