diff --git a/Pr/PrAlgorithms/src/PrDownTrack.h b/Pr/PrAlgorithms/src/PrDownTrack.h index c3ea7d2630c23a4f60da1a0c47d925da826cce4c..67a5b180bac4441d77abfb6dfdfd8a01d6b7108a 100644 --- a/Pr/PrAlgorithms/src/PrDownTrack.h +++ b/Pr/PrAlgorithms/src/PrDownTrack.h @@ -111,35 +111,35 @@ public: // using Hits = boost::container::static_vector; // Until we can put a bound on the number of hits, use a small_vector - PrDownTrack( Gaudi::TrackVector stateVector, double stateZ, double zUT, LHCb::span magnetParams, - LHCb::span yParams, LHCb::span momPar, double magnetScale ) + PrDownTrack( Gaudi::TrackVector stateVector, float stateZ, float zUT, LHCb::span magnetParams, + LHCb::span yParams, LHCb::span momPar, float magnetScale ) : m_stateVector( stateVector ), m_stateZ( stateZ ), m_zUT( zUT ) { const auto tx2 = stateTx() * stateTx(); const auto ty2 = stateTy() * stateTy(); m_momentumParam = ( momPar[0] + momPar[1] * tx2 + momPar[2] * ty2 ) * magnetScale; // -- See PrFitKsParams to see how these coefficients are derived. - const double zMagnet = magnetParams[0] + magnetParams[1] * ty2 + magnetParams[2] * tx2 + - magnetParams[3] * std::abs( stateQoP() ) + /// this is where the old one stopped. - magnetParams[4] * std::abs( stateX() ) + magnetParams[5] * std::abs( stateY() ) + - magnetParams[6] * std::abs( stateTy() ); + const float zMagnet = magnetParams[0] + magnetParams[1] * ty2 + magnetParams[2] * tx2 + + magnetParams[3] * std::abs( stateQoP() ) + /// this is where the old one stopped. + magnetParams[4] * std::abs( stateX() ) + magnetParams[5] * std::abs( stateY() ) + + magnetParams[6] * std::abs( stateTy() ); - const double dz = zMagnet - stateZ; - const double xMagnet = stateX() + dz * stateTx(); - m_slopeX = xMagnet / zMagnet; - const double dSlope = std::abs( m_slopeX - stateTx() ); - const double dSlope2 = dSlope * dSlope; + const float dz = zMagnet - stateZ; + const float xMagnet = stateX() + dz * stateTx(); + m_slopeX = xMagnet / zMagnet; + const float dSlope = std::abs( m_slopeX - stateTx() ); + const float dSlope2 = dSlope * dSlope; - const double by = stateY() / ( stateZ + ( yParams[0] * fabs( stateTy() ) * zMagnet + yParams[1] ) * dSlope2 ); - m_slopeY = by * ( 1. + yParams[0] * fabs( by ) * dSlope2 ); + const float by = stateY() / ( stateZ + ( yParams[0] * fabs( stateTy() ) * zMagnet + yParams[1] ) * dSlope2 ); + m_slopeY = by * ( 1. + yParams[0] * fabs( by ) * dSlope2 ); - const double yMagnet = stateY() + dz * by - yParams[1] * by * dSlope2; + const float yMagnet = stateY() + dz * by - yParams[1] * by * dSlope2; // -- These resolutions are semi-empirical and are obtained by fitting residuals // -- with MCHits and reconstructed tracks // -- See Tracking &Alignment meeting, 19.2.2015, for the idea - double errXMag = dSlope2 * 15.0 + dSlope * 15.0 + 3.0; - double errYMag = dSlope2 * 80.0 + dSlope * 10.0 + 4.0; + float errXMag = dSlope2 * 15.0 + dSlope * 15.0 + 3.0; + float errYMag = dSlope2 * 80.0 + dSlope * 10.0 + 4.0; // -- Assume better resolution for SciFi than for OT // -- obviously this should be properly tuned... @@ -163,59 +163,57 @@ public: } /// getters - double stateX() const noexcept { return m_stateVector[0]; } - double stateY() const noexcept { return m_stateVector[1]; } - double stateZ() const noexcept { return m_stateZ; } - double stateTx() const noexcept { return m_stateVector[2]; } - double stateTy() const noexcept { return m_stateVector[3]; } - double stateQoP() const noexcept { return m_stateVector[4]; } + float stateX() const noexcept { return m_stateVector[0]; } + float stateY() const noexcept { return m_stateVector[1]; } + float stateZ() const noexcept { return m_stateZ; } + float stateTx() const noexcept { return m_stateVector[2]; } + float stateTy() const noexcept { return m_stateVector[3]; } + float stateQoP() const noexcept { return m_stateVector[4]; } Hits& hits() noexcept { return m_hits; } const Hits& hits() const noexcept { return m_hits; } - double xMagnet() const noexcept { return m_magnet.x(); } - double yMagnet() const noexcept { return m_magnet.y(); } - double zMagnet() const noexcept { return m_magnet.z(); } - double slopeX() const noexcept { return m_slopeX; } - double slopeY() const noexcept { return m_slopeY; } - double weightXMag() const noexcept { return m_weightXMag; } - double weightYMag() const noexcept { return m_weightYMag; } - double chi2() const noexcept { return m_chi2; } - double zUT() const noexcept { return m_zUT; } + float xMagnet() const noexcept { return m_magnet.x(); } + float yMagnet() const noexcept { return m_magnet.y(); } + float zMagnet() const noexcept { return m_magnet.z(); } + float slopeX() const noexcept { return m_slopeX; } + float slopeY() const noexcept { return m_slopeY; } + float weightXMag() const noexcept { return m_weightXMag; } + float weightYMag() const noexcept { return m_weightYMag; } + float chi2() const noexcept { return m_chi2; } + float zUT() const noexcept { return m_zUT; } /// setters - void setSlopeX( double slopeX ) noexcept { m_slopeX = slopeX; } - void setChi2( double chi2 ) noexcept { m_chi2 = chi2; } + void setSlopeX( float slopeX ) noexcept { m_slopeX = slopeX; } + void setChi2( float chi2 ) noexcept { m_chi2 = chi2; } // functions - double xAtZ( double z ) const noexcept { - const double curvature = 1.6e-5 * ( stateTx() - m_slopeX ); + float xAtZ( float z ) const noexcept { + const float curvature = 1.6e-5 * ( stateTx() - m_slopeX ); return xMagnet() + ( z - zMagnet() ) * m_slopeX + curvature * ( z - m_zUT ) * ( z - m_zUT ); } - double yAtZ( double z ) const noexcept { return yMagnet() + m_displY + ( z - zMagnet() ) * slopeY(); } + float yAtZ( float z ) const noexcept { return yMagnet() + m_displY + ( z - zMagnet() ) * slopeY(); } - void updateX( double dx, double dsl ) noexcept { + void updateX( float dx, float dsl ) noexcept { m_displX += dx; m_magnet = Gaudi::XYZPoint( m_magnet.x() + dx, m_magnet.y(), m_magnet.z() ); m_slopeX += dsl; } - void updateY( double dy ) noexcept { m_displY += dy; } + void updateY( float dy ) noexcept { m_displY += dy; } - double dxMagnet() const noexcept { return -m_displX; } + float dxMagnet() const noexcept { return -m_displX; } - double initialChi2() const noexcept { - return m_displX * m_displX * m_weightXMag + m_displY * m_displY * m_weightYMag; - } + float initialChi2() const noexcept { return m_displX * m_displX * m_weightXMag + m_displY * m_displY * m_weightYMag; } - double momentum() const noexcept { return m_momentumParam / ( stateTx() - m_slopeX ); } + float momentum() const noexcept { return m_momentumParam / ( stateTx() - m_slopeX ); } - double pt() const noexcept { - const double tx2 = slopeX() * slopeX(); - const double ty2 = slopeY() * slopeY(); - const double sinTrack = sqrt( 1. - 1. / ( 1. + tx2 + ty2 ) ); + float pt() const noexcept { + const float tx2 = slopeX() * slopeX(); + const float ty2 = slopeY() * slopeY(); + const float sinTrack = sqrt( 1. - 1. / ( 1. + tx2 + ty2 ) ); return sinTrack * std::abs( momentum() ); } - double distance( const Downstream::Hit& hit ) const noexcept { return hit.xAt( yAtZ( hit.z ) ) - xAtZ( hit.z ); } + float distance( const Downstream::Hit& hit ) const noexcept { return hit.xAt( yAtZ( hit.z ) ) - xAtZ( hit.z ); } bool isYCompatible( const float tol ) const noexcept { return std ::all_of( m_hits.begin(), m_hits.end(), @@ -230,18 +228,18 @@ public: private: Gaudi::TrackVector m_stateVector; - double m_stateZ{ 0 }; + float m_stateZ{ 0 }; Gaudi::XYZPoint m_magnet; - double m_momentumParam{ 0 }; - double m_zUT{ 0 }; - double m_slopeX{ 0 }; - double m_slopeY{ 0 }; - double m_displX{ 0 }; - double m_displY{ 0 }; - double m_weightXMag{ 0 }; - double m_weightYMag{ 0 }; - double m_chi2{ 0 }; + float m_momentumParam{ 0 }; + float m_zUT{ 0 }; + float m_slopeX{ 0 }; + float m_slopeY{ 0 }; + float m_displX{ 0 }; + float m_displY{ 0 }; + float m_weightXMag{ 0 }; + float m_weightYMag{ 0 }; + float m_chi2{ 0 }; Hits m_hits; /// working list of hits on this track }; diff --git a/Pr/PrAlgorithms/src/PrLongLivedTracking.cpp b/Pr/PrAlgorithms/src/PrLongLivedTracking.cpp index 4387ed79f6cb03cbfd15eda3a72640d9a5794bdf..3361d140b2d198c54c25855871ccf1966e747136 100644 --- a/Pr/PrAlgorithms/src/PrLongLivedTracking.cpp +++ b/Pr/PrAlgorithms/src/PrLongLivedTracking.cpp @@ -121,7 +121,7 @@ public: DownstreamTracks finalTracks{ &InputTracks, Zipping::generateZipIdentifier(), LHCb::getMemResource( evtCtx ) }; finalTracks.reserve( InputTracks.size() ); - const double magScaleFactor = magnet.signedRelativeCurrent(); + const float magScaleFactor = magnet.signedRelativeCurrent(); bool magnetOff = std::abs( magScaleFactor ) > 1e-6 ? false : true; @@ -133,7 +133,7 @@ public: for ( const auto& _tr : InputTracks.scalar() ) { // -- simple Fisher discriminant to reject bad seed tracks // -- tune this! - // const double fisher = evaluateFisher( tr ); + // const float fisher = evaluateFisher( tr ); // if( fisher < m_seedCut ) continue; @@ -152,7 +152,7 @@ public: // -- tune this! // -- check for compatible momentum - const double deltaP = refTrack.momentum() * refTrack.stateQoP() - 1.; + const float deltaP = refTrack.momentum() * refTrack.stateQoP() - 1.; if ( maxDeltaP( refTrack ) < fabs( deltaP ) ) { if ( !magnetOff ) continue; } @@ -169,42 +169,17 @@ public: //============================================================== // Try to find a candidate: X first, then UV. //============================================================== - int myPlane = 0; - const int otherPlane = ( myPlane == 0 ) ? 3 : 0; - int nXTrack = 0; - for ( auto& myHit : preSelHits[myPlane] ) { - const double meanZ = myHit.z; - const double posX = myHit.x; - - PrDownTrack track( refTrack ); - - // -- Create track estimate with one x hit - const double slopeX = ( track.xMagnet() - posX ) / ( track.zMagnet() - meanZ ); - track.setSlopeX( slopeX ); - - // -- Fit x projection - findMatchingHits( track, preSelHits[otherPlane], matchingXHits ); - fitXProjection( track, myHit, matchingXHits, goodXTracks ); - - nXTrack += goodXTracks.size(); + nXTrack += createTrackCandidates( trackCandidates, refTrack, preSelHits, matchingXHits, uHitsTemp, goodXTracks, + goodXUTracks, { 0, 3, 1, 2 }, false ); - // -- Loop over good x tracks - for ( PrDownTrack& xTrack : goodXTracks ) { - // -- Take all xTracks into account whose chi2 is close to the best - // if ( xTrack.chi2() - m_maxChi2DistXTracks >= goodXTracks[0].chi2() ) break; + bool has4LayerTrack = false; + for ( PrDownTrack& track : trackCandidates ) { has4LayerTrack |= track.hits().size() == 4; } - addUHits( xTrack, preSelHits[1], uHitsTemp, goodXUTracks ); - - // -- Loop over good xu tracks - for ( PrDownTrack& xuTrack : goodXUTracks ) { - addVHits( xuTrack, preSelHits[2] ); - fitAndRemove( xuTrack ); - if ( xuTrack.chi2() > m_maxChi2 || !xuTrack.isYCompatible( m_yTol ) ) continue; - trackCandidates.push_back( std::move( xuTrack ) ); - } // Loop over good xu tracks - } // Loop over good x tracks + if ( !has4LayerTrack ) { + nXTrack += createTrackCandidates( trackCandidates, refTrack, preSelHits, matchingXHits, uHitsTemp, goodXTracks, + goodXUTracks, { 3, 0, 2, 1 }, true ); } m_nXcand += nXTrack; @@ -220,6 +195,7 @@ public: downstream_tracks_mlp_datas.reserve( 16 ); for ( PrDownTrack& track : trackCandidates ) { + if ( has4LayerTrack && track.hits().size() < 4 ) continue; addOverlapRegions( track, preSelHits ); @@ -259,54 +235,54 @@ private: // // - XPredTolConst: x-window for preselection is XPredTolConst/p + XPredTolOffset - Gaudi::Property m_xPredTolConst{ this, "XPredTolConst", 200. * Gaudi::Units::mm* Gaudi::Units::GeV }; + Gaudi::Property m_xPredTolConst{ this, "XPredTolConst", 200. * Gaudi::Units::mm* Gaudi::Units::GeV }; // - XPredTolOffset: x-window for preselection is XPredTolConst/p + XPredTolOffset - Gaudi::Property m_xPredTolOffset{ this, "XPredTolOffset", 6. * Gaudi::Units::mm }; + Gaudi::Property m_xPredTolOffset{ this, "XPredTolOffset", 6. * Gaudi::Units::mm }; // - TolMatchConst: x-window for matching x hits is TolMatchConst/p + TolMatchOffset - Gaudi::Property m_tolMatchConst{ this, "TolMatchConst", 20000. }; + Gaudi::Property m_tolMatchConst{ this, "TolMatchConst", 20000. }; // - TolMatchOffset: x-window for matching x hits is TolMatchConst/p + TolMatchOffset - Gaudi::Property m_tolMatchOffset{ this, "TolMatchOffset", 1.5 * Gaudi::Units::mm }; + Gaudi::Property m_tolMatchOffset{ this, "TolMatchOffset", 1.5 * Gaudi::Units::mm }; // - TolUConst: window for u hits is TolUConst/p + TolUOffset - Gaudi::Property m_tolUConst{ this, "TolUConst", 20000.0 }; + Gaudi::Property m_tolUConst{ this, "TolUConst", 20000.0 }; // - TolUOffset: window for u hits is TolUConst/p + TolUOffset - Gaudi::Property m_tolUOffset{ this, "TolUOffset", 2.5 }; + Gaudi::Property m_tolUOffset{ this, "TolUOffset", 2.5 }; // - TolVConst: window for v hits is TolVConst/p + TolVOffset - Gaudi::Property m_tolVConst{ this, "TolVConst", 2000.0 }; + Gaudi::Property m_tolVConst{ this, "TolVConst", 2000.0 }; // - TolVOffset: window for v hits is TolVConst/p + TolVOffset - Gaudi::Property m_tolVOffset{ this, "TolVOffset", 0.5 }; + Gaudi::Property m_tolVOffset{ this, "TolVOffset", 0.5 }; // - MaxWindowSize: maximum window size for matching x hits - Gaudi::Property m_maxWindow{ this, "MaxWindowSize", 10.0 * Gaudi::Units::mm }; + Gaudi::Property m_maxWindow{ this, "MaxWindowSize", 10.0 * Gaudi::Units::mm }; // - MaxChi2: Maximum chi2 for tracks with at least 4 hits - Gaudi::Property m_maxChi2{ this, "MaxChi2", 20. }; + Gaudi::Property m_maxChi2{ this, "MaxChi2", 20. }; // - MaxChi2ThreeHits: Maximum chi2 for tracks with 3 hits - Gaudi::Property m_maxChi2ThreeHits{ this, "MaxChi2ThreeHits", 10.0 }; + Gaudi::Property m_maxChi2ThreeHits{ this, "MaxChi2ThreeHits", 10.0 }; // - MinUTx: half-width of beampipe rectangle - Gaudi::Property m_minUTx{ this, "MinUTx", 25. * Gaudi::Units::mm }; + Gaudi::Property m_minUTx{ this, "MinUTx", 25. * Gaudi::Units::mm }; // - MinUTy: half-height of of beampipe rectangle - Gaudi::Property m_minUTy{ this, "MinUTy", 25. * Gaudi::Units::mm }; + Gaudi::Property m_minUTy{ this, "MinUTy", 25. * Gaudi::Units::mm }; // - MaxGhostProb: Maximum ghost probability for tracks - Gaudi::Property m_maxGhostProb{ this, "MaxGhostProb", 0.75 }; + Gaudi::Property m_maxGhostProb{ this, "MaxGhostProb", 0.75 }; // Define parameters for MC09 field, zState = 9410 // - ZMagnetParams: Parameters to determine the z-position of the magnet point. Tune with PrKsParams. - Gaudi::Property> m_zMagnetParams{ + Gaudi::Property> m_zMagnetParams{ this, "ZMagnetParams", { 5379.88, -2143.93, 366.124, 119074, -0.0100333, -0.146055, 1260.96 } }; // - MomentumParams: Parameters to determine the momentum. Tune with PrKsParams. - Gaudi::Property> m_momentumParams{ this, "MomentumParams", { 1217.77, 454.598, 3353.39 } }; + Gaudi::Property> m_momentumParams{ this, "MomentumParams", { 1217.77, 454.598, 3353.39 } }; // - YParams: Parameters to determine the bending in y. Tune with PrKsParams. - Gaudi::Property> m_yParams{ this, "YParams", { 5., 2000. } }; + Gaudi::Property> m_yParams{ this, "YParams", { 5.f, 2000.f } }; // - ZUT: z-position of middle of UT. - Gaudi::Property m_zUT{ this, "ZUT", 2485. * Gaudi::Units::mm }; + Gaudi::Property m_zUT{ this, "ZUT", 2485. * Gaudi::Units::mm }; // - ZUTa: z-position of first UT station - Gaudi::Property m_zUTa{ this, "ZUTa", 2350. * Gaudi::Units::mm }; + Gaudi::Property m_zUTa{ this, "ZUTa", 2350. * Gaudi::Units::mm }; // - InitialMinPt: Minimum pT of the track from initial estimate - Gaudi::Property m_initialMinPt{ this, "InitialMinPt", 0. * Gaudi::Units::MeV }; + Gaudi::Property m_initialMinPt{ this, "InitialMinPt", 0. * Gaudi::Units::MeV }; // - InitialMinMomentum: Minimum momentum of the track from initial estimate - Gaudi::Property m_initialMinMomentum{ this, "InitialMinMomentum", 1400. * Gaudi::Units::MeV }; + Gaudi::Property m_initialMinMomentum{ this, "InitialMinMomentum", 1400. * Gaudi::Units::MeV }; // - MinPt: Minimum pT of the track - Gaudi::Property m_minPt{ this, "MinPt", 0. * Gaudi::Units::MeV }; + Gaudi::Property m_minPt{ this, "MinPt", 0. * Gaudi::Units::MeV }; // - MinMomentum: Minimum momentum of the track - Gaudi::Property m_minMomentum{ this, "MinMomentum", 0. * Gaudi::Units::MeV }; + Gaudi::Property m_minMomentum{ this, "MinMomentum", 0. * Gaudi::Units::MeV }; // - MinNumUTHits: Minimum number of UT hits required Gaudi::Property m_MinNumUTHits{ this, "MinNumUTHits", 3 }; @@ -315,35 +291,35 @@ private: // -- Parameter to reject seed track which are likely ghosts // - FisherCut: Cut on Fisher-discriminant to reject bad seed tracks. - Gaudi::Property m_seedCut{ this, "FisherCut", -1.0 }; + Gaudi::Property m_seedCut{ this, "FisherCut", -1.0 }; // -- Parameters for the cut on deltaP (momentum estimate from Seeding and Downstream kink) // - MaxDeltaPConst: Window for deltaP is: MaxDeltaPConst/p + MaxDeltaPOffset - Gaudi::Property m_maxDeltaPConst{ this, "MaxDeltaPConst", 0.0 }; + Gaudi::Property m_maxDeltaPConst{ this, "MaxDeltaPConst", 0.0 }; // - MaxDeltaPOffset: Window for deltaP is: MaxDeltaPConst/p + MaxDeltaPOffset - Gaudi::Property m_maxDeltaPOffset{ this, "MaxDeltaPOffset", 0.25 }; + Gaudi::Property m_maxDeltaPOffset{ this, "MaxDeltaPOffset", 0.25 }; // -- Parameters for correcting the predicted position // - XCorrectionConst: Correction for x-position of search window in preselection is XCorrectionConst/p + // XCorrestionOffset - Gaudi::Property m_xCorrectionConst{ this, "XCorrectionConst", 23605.0 }; + Gaudi::Property m_xCorrectionConst{ this, "XCorrectionConst", 23605.0 }; // - XCorrectionOffset: Correction for x-position of search window in preselection is XCorrectionConst/p + // XCorrestionOffset - Gaudi::Property m_xCorrectionOffset{ this, "XCorrestionOffset", 0.4 }; + Gaudi::Property m_xCorrectionOffset{ this, "XCorrestionOffset", 0.4 }; // - MaxXTracks: Maximum number of x-tracklets to process further Gaudi::Property m_maxXTracks{ this, "MaxXTracks", 2 }; // - MaxChi2DistXTracks: Maximum chi2 difference to x-tracklet with best chi2 - Gaudi::Property m_maxChi2DistXTracks{ this, "MaxChi2DistXTracks", 0.2 }; + Gaudi::Property m_maxChi2DistXTracks{ this, "MaxChi2DistXTracks", 0.2 }; // - MaxXUTracks: Maximum number of xu-tracklets to process further Gaudi::Property m_maxXUTracks{ this, "MaxXUTracks", 3 }; - Gaudi::Property m_fitXProjChi2Offset{ this, "FitXProjChi2Offset", 4.5 }; - Gaudi::Property m_fitXProjChi2Const{ this, "FitXProjChi2Const", 35000.0 }; + Gaudi::Property m_fitXProjChi2Offset{ this, "FitXProjChi2Offset", 4.5 }; + Gaudi::Property m_fitXProjChi2Const{ this, "FitXProjChi2Const", 35000.0 }; // -- Tolerance for adding overlap hits // - OverlapTol: Tolerance for adding overlap hits - Gaudi::Property m_overlapTol{ this, "OverlapTol", 2.0 * Gaudi::Units::mm }; + Gaudi::Property m_overlapTol{ this, "OverlapTol", 2.0 * Gaudi::Units::mm }; // - YTol: YTolerance for adding / removing hits. - Gaudi::Property m_yTol{ this, "YTol", 2.0 * Gaudi::Units::mm }; + Gaudi::Property m_yTol{ this, "YTol", 2.0 * Gaudi::Units::mm }; // Change this in order to remove hits and T-tracks used for longtracks. // RemoveAll configures that everything is removed. // If false only hits and T-tracks from good longtracks are removed. @@ -353,7 +329,7 @@ private: // - RemoveAll: Remove seed tracks and used UT hits (withoug chi2-cut on long track)? Gaudi::Property m_removeAll{ this, "RemoveAll", false }; // - LongChi2: Chi2-cut for the removal - Gaudi::Property m_longChi2{ this, "LongChi2", 1.5 }; + Gaudi::Property m_longChi2{ this, "LongChi2", 1.5 }; // properties Gaudi::Property m_weightsfilename{ this, "WeightsFileName", @@ -370,6 +346,48 @@ private: // void prepareSeeds(const Tracks& inTracks, std::vector& myInTracks)const; ///< Tag already used T-Seeds + int createTrackCandidates( PrDownTracks& trackCandidates, const PrDownTrack& refTrack, + std::array& preSelHits, Downstream::Hits matchingXHits, + Downstream::Hits uHitsTemp, PrDownTracks goodXTracks, PrDownTracks goodXUTracks, + const std::array plane_indx, bool skip_second_layer = false ) const { + int nXTrack = 0; + + for ( auto& myHit : preSelHits[plane_indx[0]] ) { + const float meanZ = myHit.z; + const float posX = myHit.x; + + PrDownTrack track( refTrack ); + + // -- Create track estimate with one x hit + const float slopeX = ( track.xMagnet() - posX ) / ( track.zMagnet() - meanZ ); + track.setSlopeX( slopeX ); + + // -- Fit x projection + findMatchingHits( track, preSelHits[plane_indx[1]], matchingXHits ); + fitXProjection( track, myHit, matchingXHits, goodXTracks, skip_second_layer ); + + nXTrack += goodXTracks.size(); + + // -- Loop over good x tracks + for ( PrDownTrack& xTrack : goodXTracks ) { + // -- Take all xTracks into account whose chi2 is close to the best + // if ( xTrack.chi2() - m_maxChi2DistXTracks >= goodXTracks[0].chi2() ) break; + + addUHits( xTrack, preSelHits[plane_indx[2]], uHitsTemp, goodXUTracks ); + + // -- Loop over good xu tracks + for ( PrDownTrack& xuTrack : goodXUTracks ) { + addVHits( xuTrack, preSelHits[plane_indx[3]] ); + if ( xuTrack.hits().size() < 3 ) continue; + fitAndRemove( xuTrack ); + if ( xuTrack.chi2() > m_maxChi2 || !xuTrack.isYCompatible( m_yTol ) ) continue; + trackCandidates.push_back( std::move( xuTrack ) ); + } // Loop over good xu tracks + } // Loop over good x tracks + } + return nXTrack; + } + //========================================================================= // Get the PreSelection of hits around a first track estimate //========================================================================= @@ -377,21 +395,21 @@ private: const LHCb::Pr::Hits& hitHandler, const LHCb::UTDAQ::GeomCache& geom ) const { // - Max Pt around 100 MeV for strange particle decay -> maximum displacement is in 1/p. - double xPredTol = m_xPredTolOffset; + float xPredTol = m_xPredTolOffset; // P dependance + overal tol. auto p = std::abs( track.momentum() ); if ( p > 1e-6 ) xPredTol = m_xPredTolConst / p + m_xPredTolOffset; - const double yTol = xPredTol / 2.0 + 7.5; // this is a little vague and not the final word + const float yTol = xPredTol / 2.0 + 7.5; // this is a little vague and not the final word // -- a correction turns out to be beneficial // -- maybe to compensate tracks not coming from 0/0 (?) - const double correction = xPosCorrection( track ); + const float correction = xPosCorrection( track ); for ( auto& i : preSelHits ) i.clear(); - const double yTrack = track.yAtZ( 0. ); - const double tyTr = track.slopeY(); + const float yTrack = track.yAtZ( 0. ); + const float tyTr = track.slopeY(); boost::container::small_vector sectors; @@ -402,7 +420,7 @@ private: sectors.clear(); - const double zLayer = geom.layers[layer].z; + const float zLayer = geom.layers[layer].z; geom.findSectorsFullID( layer, track.xAtZ( zLayer ), track.yAtZ( zLayer ), xPredTol, yTol + 20.0, sectors ); std::sort( sectors.begin(), sectors.end() ); @@ -419,16 +437,16 @@ private: const auto myHs = hitHandler.scalar(); const auto firstHit = myHs[sector.first]; - const double dxDy = firstHit.get().cast(); - const double zLayer = firstHit.get().cast(); - const double yPredLay = track.yAtZ( zLayer ); - const double xPredLay = track.xAtZ( zLayer ); + const float dxDy = firstHit.get().cast(); + const float zLayer = firstHit.get().cast(); + const float yPredLay = track.yAtZ( zLayer ); + const float xPredLay = track.xAtZ( zLayer ); - const double pos = xPredLay - correction; - const auto y = yTrack + tyTr * zLayer; + const float pos = xPredLay - correction; + const auto y = yTrack + tyTr * zLayer; // this should sort of take the stereo angle and some tolerance into account. - const double lowerBoundX = xPredLay - xPredTol - dxDy * yPredLay - 2.0; + const float lowerBoundX = xPredLay - xPredTol - dxDy * yPredLay - 2.0; for ( auto itHit = sector.first; itHit < sector.second; ++itHit ) { @@ -470,7 +488,7 @@ private: again = false; //== Fit, using the magnet point as constraint. - double mat[6], rhs[3]; + float mat[6], rhs[3]; mat[0] = track.weightXMag(); mat[1] = 0.; mat[2] = 0.; @@ -484,10 +502,10 @@ private: std::array differentPlanes = { 0, 0, 0, 0 }; for ( auto& hit : track.hits() ) { - const double dz = 0.001 * ( hit.z - track.zMagnet() ); - const double dist = track.distance( hit ); - const double w = hit.weight(); - const double t = hit.sin(); + const float dz = 0.001 * ( hit.z - track.zMagnet() ); + const float dist = track.distance( hit ); + const float w = hit.weight(); + const float t = hit.sin(); mat[0] += w; mat[1] += w * dz; @@ -508,7 +526,7 @@ private: const int nbUV = differentPlanes[1] + differentPlanes[2]; // -- solve the equation and update the parameters of the track - CholeskyDecomp decomp( mat ); + CholeskyDecomp decomp( mat ); if ( !decomp ) { track.setChi2( 1e42 ); return; @@ -516,24 +534,24 @@ private: decomp.Solve( rhs ); } - const double dx = rhs[0]; - const double dsl = 0.001 * rhs[1]; - const double dy = rhs[2]; + const float dx = rhs[0]; + const float dsl = 0.001 * rhs[1]; + const float dy = rhs[2]; if ( nbUV < 4 ) track.updateX( dx, dsl ); track.updateY( dy ); //== Remove worst hit and retry, if too far. - double chi2 = track.initialChi2(); + float chi2 = track.initialChi2(); - double maxDist = -1.; + float maxDist = -1.; PrDownTrack::Hits::iterator worst; for ( auto itH = track.hits().begin(); itH != track.hits().end(); ++itH ) { Downstream::Hit& hit = *itH; if ( !onlyFit ) { - const double yTrackAtZ = track.yAtZ( hit.z ); + const float yTrackAtZ = track.yAtZ( hit.z ); if ( !hit.isYCompatible( yTrackAtZ, m_yTol ) ) { track.hits().erase( itH ); if ( 2 < track.hits().size() ) again = true; @@ -541,8 +559,8 @@ private: } } - const double dist = std::abs( track.distance( hit ) ); - hit.projection = dist; + const float dist = std::abs( track.distance( hit ) ); + hit.projection = dist; chi2 += dist * dist * hit.weight(); // -- only flag this hit as removable if it is not alone in a plane or there are 4 planes that fired if ( !onlyFit && maxDist < dist && ( 1 < differentPlanes[hit.planeCode()] || nDoF == track.hits().size() ) ) { @@ -573,12 +591,12 @@ private: matchingXHits.clear(); if ( preSelHits.empty() ) return; - const double tol = + const float tol = std::min( m_maxWindow.value(), m_tolMatchOffset + m_tolMatchConst * std::abs( track.stateQoP() ) ); - const double xPred = track.xAtZ( preSelHits.front().z ); + const float xPred = track.xAtZ( preSelHits.front().z ); for ( auto& hit : preSelHits ) { - const double adist = std::abs( hit.x - xPred ); + const float adist = std::abs( hit.x - xPred ); if ( adist <= tol ) { matchingXHits.push_back( hit ); } } @@ -635,15 +653,15 @@ private: uHitsTemp.clear(); - const double tol = m_tolUOffset + m_tolUConst / std::abs( track.momentum() ); + const float tol = m_tolUOffset + m_tolUConst / std::abs( track.momentum() ); // -- these numbers are a little arbitrary - double minChi2 = ( track.hits().size() == 1 ) ? 800 : 300; + float minChi2 = ( track.hits().size() == 1 ) ? 800 : 300; // -- first select all hits, and then // -- accept until over a tolerance for ( auto& hit : preSelHits ) { - const double dist = std::abs( track.distance( hit ) ); + const float dist = std::abs( track.distance( hit ) ); if ( dist > tol ) continue; hit.projection = dist; uHitsTemp.push_back( hit ); @@ -669,6 +687,9 @@ private: } if ( goodXUTracks.size() >= m_maxXUTracks ) { break; } } + + // 3-hit case handling + if ( ( goodXUTracks.size() == 0 ) && ( track.hits().size() == 2 ) ) { goodXUTracks.emplace_back( track ); } } //========================================================================= @@ -677,30 +698,30 @@ private: void addVHits( PrDownTrack& track, Downstream::Hits& preSelHits ) const { if ( preSelHits.empty() ) { return; } - auto p = std::abs( track.momentum() ); - double tol = ( track.hits().size() == 2 ) ? m_tolUOffset + m_tolUConst / p : m_tolVOffset + m_tolVConst / p; - - double minChi2 = 10000; + auto p = std::abs( track.momentum() ); + float tol = ( track.hits().size() == 2 ) ? m_tolUOffset + m_tolUConst / p : m_tolVOffset + m_tolVConst / p; - Downstream::Hit* bestHit = nullptr; + // Find the closest hit to the track + // It's enough to look at the closest hit only, since absolute majority of the tracks + // with not a closest hit get killed later anyways + // This also doesn't affect efficiency, as the addOverlapRegions function probes all close enough hits as well + float bestDist = 1e9f; + Downstream::Hit* bestHit = nullptr; for ( auto& hit : preSelHits ) { - const double adist = std::abs( track.distance( hit ) ); - - if ( adist < tol ) { - hit.projection = adist; - track.hits().push_back( hit ); - fitAndRemove( track ); - track.hits().pop_back(); - - if ( track.chi2() < minChi2 ) { - bestHit = &hit; - minChi2 = track.chi2(); - } + const float dist = std::abs( track.distance( hit ) ); + if ( dist > tol ) continue; + hit.projection = dist; + if ( dist < bestDist ) { + bestDist = dist; + bestHit = &hit; } } - if ( bestHit ) track.hits().push_back( *bestHit ); - + // Fit the track + if ( bestHit ) { + track.hits().push_back( *bestHit ); + fitAndRemove( track ); + } track.sortFinalHits(); } @@ -710,10 +731,10 @@ private: // Fit the projection in the zx plane, one hit in each x layer //============================================================================= void fitXProjection( const PrDownTrack& track, Downstream::Hit& firstHit, Downstream::Hits& matchingXHits, - PrDownTracks& goodXTracks ) const { + PrDownTracks& goodXTracks, const bool skip_second_layer ) const { goodXTracks.clear(); - const double maxChi2 = m_fitXProjChi2Offset + m_fitXProjChi2Const / std::abs( track.momentum() ); + const float maxChi2 = m_fitXProjChi2Offset + m_fitXProjChi2Const / std::abs( track.momentum() ); // Catch if there is no second hit in other station for ( const auto& hit : matchingXHits ) { @@ -729,6 +750,17 @@ private: if ( goodXTracks.size() >= 3 ) { break; } } + + if ( skip_second_layer && goodXTracks.size() != 0 ) { + goodXTracks.clear(); + return; + } + + // 3-hit case handling + if ( goodXTracks.size() == 0 ) { + auto& tr = goodXTracks.emplace_back( track ); + tr.hits().push_back( firstHit ); + } } //========================================================================= @@ -742,14 +774,14 @@ private: // -- use a tighter chi2 for 3 hit tracks // -- as they are more likely to be ghosts - const double maxChi2 = ( nbMeasureOK == 3 ) ? m_maxChi2ThreeHits : m_maxChi2; + const float maxChi2 = ( nbMeasureOK == 3 ) ? m_maxChi2ThreeHits : m_maxChi2; //== Good enough Chi2/ndof if ( maxChi2 < track.chi2() ) { return false; } //== Compatible momentum - const double p = track.momentum(); - const double deltaP = p * track.stateQoP() - 1.; + const float p = track.momentum(); + const float deltaP = p * track.stateQoP() - 1.; if ( maxDeltaP( track ) < fabs( deltaP ) && !magnetOff ) { return false; } if ( std::abs( p ) < m_minMomentum || track.pt() < m_minPt ) { return false; } @@ -773,7 +805,7 @@ private: for ( auto& hit : preSelHits[trackHit.planeCode()] ) { if ( m_overlapTol > std::abs( track.distance( hit ) ) ) { - double yTrack = track.yAtZ( hit.z ); + float yTrack = track.yAtZ( hit.z ); if ( !hit.isYCompatible( yTrack, m_yTol ) ) continue; @@ -795,12 +827,12 @@ private: //============================================================================= // Evaluate the Fisher discriminant for a preselection of seed tracks //============================================================================= - /* double evaluateFisher( const Track* track ) { + /* float evaluateFisher( const Track* track ) { const unsigned int nbIT = std::count_if( track->lhcbIDs().begin(), track->lhcbIDs().end(), [](const LHCb::LHCbID id){ return id.isIT();}); - double nbITD = static_cast(nbIT); - double lhcbIDSizeD = static_cast(track->lhcbIDs().size()); - std::array vals = { track->chi2PerDoF(), track->p(), track->pt(), nbITD, lhcbIDSizeD }; + float nbITD = static_cast(nbIT); + float lhcbIDSizeD = static_cast(track->lhcbIDs().size()); + std::array vals = { track->chi2PerDoF(), track->p(), track->pt(), nbITD, lhcbIDSizeD }; return getFisher( vals ); } */ @@ -809,7 +841,6 @@ private: // Fit hits in x layers, using the magnet point as constraint. //========================================================================= void xFit( PrDownTrack& track, const Downstream::Hit& hit1, const Downstream::Hit& hit2 ) const { - double mat[3], rhs[2]; const auto w1 = hit1.weight(); const auto w2 = hit2.weight(); @@ -820,24 +851,20 @@ private: const auto dz1 = hit1.z - track.zMagnet(); const auto dz2 = hit2.z - track.zMagnet(); - mat[0] = track.weightXMag() + w1 + w2; - mat[1] = w1 * dz1 + w2 * dz2; - mat[2] = w1 * dz1 * dz1 + w2 * dz2 * dz2; - - rhs[0] = w1 * d1 + w2 * d2; - rhs[1] = w1 * d1 * dz1 + w2 * d2 * dz2; + auto mat = std::array{ track.weightXMag() + w1 + w2, w1 * dz1 + w2 * dz2, w1 * dz1 * dz1 + w2 * dz2 * dz2 }; + auto rhs = std::array{ w1 * d1 + w2 * d2, w1 * d1 * dz1 + w2 * d2 * dz2 }; // Solve linear system - const double det = mat[0] * mat[2] - mat[1] * mat[1]; - const double dx = ( mat[2] * rhs[0] - mat[1] * rhs[1] ) / det; - const double dsl = ( mat[0] * rhs[1] - mat[1] * rhs[0] ) / det; + const float det = mat[0] * mat[2] - mat[1] * mat[1]; + const float dx = ( mat[2] * rhs[0] - mat[1] * rhs[1] ) / det; + const float dsl = ( mat[0] * rhs[1] - mat[1] * rhs[0] ) / det; track.updateX( dx, dsl ); d1 = track.distance( hit1 ); d2 = track.distance( hit2 ); - const double chi2 = track.initialChi2() + w1 * d1 * d1 + w2 * d2 * d2; + const float chi2 = track.initialChi2() + w1 * d1 * d1 + w2 * d2 * d2; track.setChi2( chi2 ); } @@ -847,22 +874,22 @@ private: } /// Helper to evaluate the Fisher discriminant - /* double getFisher(const std::array vals) { - double c_fishConst = -1.69860581797; - std::array c_fishCoefficients = {-0.241020410138, 3.03197732663e-07, -1.14400162824e-05, + /* float getFisher(const std::array vals) { + float c_fishConst = -1.69860581797; + std::array c_fishCoefficients = {-0.241020410138, 3.03197732663e-07, -1.14400162824e-05, 0.126857153245, 0.122359738469}; - double fishVal = c_fishConst; + float fishVal = c_fishConst; for(int i = 0; i < 5; i++) fishVal += vals[i] * c_fishCoefficients[i]; return fishVal; } */ /// Helper to evaluate the maximum discrepancy between momentum from kink and curvature in T-stations - double maxDeltaP( const PrDownTrack& track ) const { + float maxDeltaP( const PrDownTrack& track ) const { return m_maxDeltaPConst / std::abs( track.momentum() ) + m_maxDeltaPOffset; } /// Helper to evaluate the correction to the x position in the UT - double xPosCorrection( const PrDownTrack& track ) const { + float xPosCorrection( const PrDownTrack& track ) const { auto p = track.momentum(); return std::copysign( m_xCorrectionOffset + m_xCorrectionConst / std::abs( p ), p ); }