diff --git a/Pr/PrPixel/src/VeloKalman.cpp b/Pr/PrPixel/src/VeloKalman.cpp
index c9850c37f668e792806ef1a671ec9dc4d4ed68da..566f92871fd2e965cd62b1ec6082a6659b98764d 100644
--- a/Pr/PrPixel/src/VeloKalman.cpp
+++ b/Pr/PrPixel/src/VeloKalman.cpp
@@ -22,6 +22,7 @@
#include "Event/PrHits.h"
#include "Event/PrLongTracks.h"
#include "Event/PrVeloTracks.h"
+#include "Event/Track_v1.h"
#include "Event/Track_v3.h"
#include "VeloKalmanHelpers.h"
@@ -130,11 +131,93 @@ namespace LHCb::Pr::Velo {
}
return out;
- };
+ }
+
+ private:
+ mutable Gaudi::Accumulators::SummingCounter<> m_nbTracksCounter{this, "Nb of Produced Tracks"};
+ };
+
+ class KalmanTrackV1 : public Algorithm::Transformer<Event::v1::Tracks( const LHCb::Track::Range&, const DeVP&,
+ const Conditions::InteractionRegion& ),
+ DetDesc::usesConditions<DeVP, Conditions::InteractionRegion>> {
+ public:
+ KalmanTrackV1( const std::string& name, ISvcLocator* pSvcLocator )
+ : Transformer( name, pSvcLocator,
+ {KeyValue{"TracksLocation", "Rec/Track/VP"}, KeyValue{"DEVP", Det::VP::det_path},
+ KeyValue{"InteractionRegionCache", "AlgorithmSpecific-" + name + "-InteractionRegion"}},
+ KeyValue{"OutputTracksLocation", "Rec/Track/Fit"} ) {}
+
+ StatusCode initialize() override {
+ return Transformer::initialize().andThen( [&]() {
+ // This is only needed to have a fallback in case the IR condition does not exist. In that case, the information
+ // is taken from DeVP and thus the Velo motion system which is not exactly the same.
+ Conditions::InteractionRegion::addConditionDerivation( this, inputLocation<Conditions::InteractionRegion>() );
+ } );
+ }
+
+ Event::v1::Tracks operator()( const LHCb::Track::Range& input_tracks, const DeVP& deVP,
+ const LHCb::Conditions::InteractionRegion& region ) const override {
+ auto out_tracks = LHCb::Event::v1::Tracks{};
+ out_tracks.reserve( input_tracks.size() );
+
+ std::vector<LHCb::LinAlg::Vec<float, 3>> hits;
+ hits.reserve( 50 );
+
+ for ( auto const* track : input_tracks ) {
+ auto new_track = new Event::v1::Track( *track );
+
+ hits.clear();
+ for ( const LHCb::VPMicroCluster& cluster : new_track->vpClusters() ) {
+ auto chanId = cluster.channelID();
+ auto fx = static_cast<float>( cluster.xfraction() ) / 255.f;
+ auto fy = static_cast<float>( cluster.yfraction() ) / 255.f;
+ auto ltg = deVP.ltg( chanId.sensor() );
+ float cx = static_cast<float>( chanId.scol() );
+ float cy = static_cast<float>( chanId.row() );
+ const float local_x = deVP.local_x( cx ) + fx * deVP.x_pitch( cx );
+ const float local_y = ( 0.5f + cy + fy ) * deVP.pixel_size();
+ const float gx = ( ltg[0] * local_x + ltg[1] * local_y + ltg[9] );
+ const float gy = ( ltg[3] * local_x + ltg[4] * local_y + ltg[10] );
+ const float gz = ( ltg[6] * local_x + ltg[7] * local_y + ltg[11] );
+ hits.emplace_back( gx, gy, gz );
+ }
+
+ auto* state = new_track->stateAt( LHCb::State::Location::ClosestToBeam );
+ auto tx = state->tx();
+ auto ty = state->ty();
+
+ if ( m_resetFirstState.value() ) {
+ // Seed with a simple fit, similar to what is done in VeloClusterTracking
+ const auto& p1 = hits[hits.size() - 3];
+ const auto& p2 = hits[hits.size() - 1];
+ auto d = p1 - p2;
+ tx = d.x() / d.z();
+ ty = d.y() / d.z();
+ }
+
+ LHCb::LinAlg::Vec<float, 3> dir{tx, ty, 1.f};
+ if ( new_track->isVeloBackward() ) {
+ auto s = fitBackwardV1( hits, dir );
+ s.transportTo( s.zBeam( region.avgPosition.x(), region.avgPosition.y() ) );
+ state->setState( s.x, s.y, s.z, s.tx, s.ty, state->qOverP() );
+ } else {
+ auto s = fitForwardV1( hits, dir );
+ s.transportTo( s.zBeam( region.avgPosition.x(), region.avgPosition.y() ) );
+ state->setState( s.x, s.y, s.z, s.tx, s.ty, state->qOverP() );
+ }
+
+ out_tracks.insert( new_track );
+ }
+
+ m_nbTracksCounter += out_tracks.size();
+ return out_tracks;
+ }
private:
mutable Gaudi::Accumulators::SummingCounter<> m_nbTracksCounter{this, "Nb of Produced Tracks"};
+ Gaudi::Property<bool> m_resetFirstState{this, "ResetFirstState", true};
};
} // namespace LHCb::Pr::Velo
DECLARE_COMPONENT_WITH_ID( LHCb::Pr::Velo::Kalman, "VeloKalman" )
+DECLARE_COMPONENT_WITH_ID( LHCb::Pr::Velo::KalmanTrackV1, "VeloKalmanTrackV1" )
diff --git a/Pr/PrPixel/src/VeloKalmanHelpers.h b/Pr/PrPixel/src/VeloKalmanHelpers.h
index 5df693d5148203ec8d91625441ddec396803bef8..9a7664487677e49df1bdd7640ccc7ab9001e53ad 100644
--- a/Pr/PrPixel/src/VeloKalmanHelpers.h
+++ b/Pr/PrPixel/src/VeloKalmanHelpers.h
@@ -64,7 +64,12 @@ public:
const T x0 = x - z * tx;
const T y0 = y - z * ty;
T denom = tx * tx + ty * ty;
- return select( denom < 0.001f * 0.001f, z, ( ( beamspot_x - x0 ) * tx + ( beamspot_y - y0 ) * ty ) / denom );
+ if constexpr ( std::is_same_v<T, float> ) {
+ if ( denom < 0.001f * 0.001f ) { return z; }
+ return ( ( beamspot_x - x0 ) * tx + ( beamspot_y - y0 ) * ty ) / denom;
+ } else {
+ return select( denom < 0.001f * 0.001f, z, ( ( beamspot_x - x0 ) * tx + ( beamspot_y - y0 ) * ty ) / denom );
+ }
}
inline void transportTo( const T& toZ ) {
@@ -80,6 +85,11 @@ public:
covYY = covYY + dz2 * covTyTy + 2.f * dz * covYTy;
covYTy = covYTy + dz * covTyTy;
}
+
+ inline T noise2PerLayer() const {
+ // Parameters for kalmanfit scattering. calibrated on MC, shamelessly hardcoded:
+ return 1e-8f + 7e-6f * ( tx * tx + ty * ty );
+ }
};
template <typename M, typename F>
@@ -122,8 +132,7 @@ fitBackward( const M track_mask, I& nHits, const LHCb::Pr::VP::Hits& hits, LHCb:
FittedState<F> s = FittedState<F>( pos, dir, 100.0f, 0.f, 0.01f, 100.0f, 0.f, 0.01f );
- // Parameters for kalmanfit scattering. calibrated on MC, shamelessly hardcoded:
- const F noise2PerLayer = 1e-8f + 7e-6f * ( s.tx * s.tx + s.ty * s.ty );
+ const F noise2PerLayer = s.noise2PerLayer();
for ( int i = 1; i < maxHits; i++ ) {
auto mask = track_mask && ( I( i ) < nHits );
@@ -156,8 +165,7 @@ fitForward( const M track_mask, I& nHits, const LHCb::Pr::VP::Hits& hits, LHCb::
FittedState<F> s = FittedState<F>( pos, dir, 100.0f, 0.f, 0.01f, 100.0f, 0.f, 0.01f );
- // Parameters for kalmanfit scattering. calibrated on MC, shamelessly hardcoded:
- const F noise2PerLayer = 1e-8f + 7e-6f * ( s.tx * s.tx + s.ty * s.ty );
+ const F noise2PerLayer = s.noise2PerLayer();
for ( int i = maxHits - 2; i >= 0; i-- ) {
I idxHit = vp_index[i];
@@ -182,6 +190,73 @@ fitForward( const M track_mask, I& nHits, const LHCb::Pr::VP::Hits& hits, LHCb::
return s;
}
+inline __attribute__( ( always_inline ) ) void filterV1( const float z, float& x, float& tx, float& covXX,
+ float& covXTx, float& covTxTx, const float zhit,
+ const float xhit, const float winv ) {
+ // compute prediction
+ const float dz = zhit - z;
+ const float predx = x + dz * tx;
+
+ const float dz_t_covTxTx = dz * covTxTx;
+ const float predcovXTx = covXTx + dz_t_covTxTx;
+ const float dz_t_covXTx = dz * covXTx;
+
+ const float predcovXX = covXX + 2.f * dz_t_covXTx + dz * dz_t_covTxTx;
+
+ // compute the gain matrix
+ const float R = 1.0f / ( winv + predcovXX );
+ const float Kx = predcovXX * R;
+ const float KTx = predcovXTx * R;
+
+ // update the state vector
+ const float r = xhit - predx;
+ x = predx + Kx * r;
+ tx = tx + KTx * r;
+
+ // update the covariance matrix
+ covTxTx = R * ( covTxTx * ( winv + covXX ) - covXTx * covXTx );
+ covXTx = winv * KTx;
+ covXX = winv * Kx;
+}
+
+FittedState<float> inline __attribute__( ( always_inline ) )
+fitBackwardV1( std::vector<LHCb::LinAlg::Vec<float, 3>>& hits, LHCb::LinAlg::Vec<float, 3>& dir ) {
+ FittedState<float> s = FittedState<float>( hits[0], dir, 100.0f, 0.f, 0.01f, 100.0f, 0.f, 0.01f );
+
+ const float noise2PerLayer = s.noise2PerLayer();
+ for ( int i = 1; i < (int)hits.size(); i++ ) {
+ auto hit = hits[i];
+ s.covTxTx = s.covTxTx + noise2PerLayer;
+ s.covTyTy = s.covTyTy + noise2PerLayer;
+
+ filterV1( s.z, s.x, s.tx, s.covXX, s.covXTx, s.covTxTx, hit.z(), hit.x(), VeloKalmanParam::wx );
+ filterV1( s.z, s.y, s.ty, s.covYY, s.covYTy, s.covTyTy, hit.z(), hit.y(), VeloKalmanParam::wy );
+ s.z = hit.z();
+ }
+ s.covTxTx = s.covTxTx + noise2PerLayer;
+ s.covTyTy = s.covTyTy + noise2PerLayer;
+ return s;
+}
+
+FittedState<float> inline __attribute__( ( always_inline ) )
+fitForwardV1( std::vector<LHCb::LinAlg::Vec<float, 3>>& hits, LHCb::LinAlg::Vec<float, 3>& dir ) {
+ FittedState<float> s = FittedState<float>( hits[hits.size() - 1], dir, 100.0f, 0.f, 0.01f, 100.0f, 0.f, 0.01f );
+
+ const float noise2PerLayer = s.noise2PerLayer();
+ for ( int i = hits.size() - 2; i >= 0; i-- ) {
+ auto hit = hits[i];
+ s.covTxTx = s.covTxTx + noise2PerLayer;
+ s.covTyTy = s.covTyTy + noise2PerLayer;
+
+ filterV1( s.z, s.x, s.tx, s.covXX, s.covXTx, s.covTxTx, hit.z(), hit.x(), VeloKalmanParam::wx );
+ filterV1( s.z, s.y, s.ty, s.covYY, s.covYTy, s.covTyTy, hit.z(), hit.y(), VeloKalmanParam::wy );
+ s.z = hit.z();
+ }
+ s.covTxTx = s.covTxTx + noise2PerLayer;
+ s.covTyTy = s.covTyTy + noise2PerLayer;
+ return s;
+}
+
template <typename M, typename F>
inline F filterWithMomentum( const M mask, const F z, F& x, F& tx, F& covXX, F& covXTx, F& covTxTx, const F zhit,
const F xhit, const F winv, const F qop ) {