diff --git a/Pr/PrAlgorithms/src/PrForwardTracking.cpp b/Pr/PrAlgorithms/src/PrForwardTracking.cpp index 5138182846080177c331245e3af661ec6ccb3b84..c3939db7c9a822709de75aaf8ba021c792b4a7b0 100644 --- a/Pr/PrAlgorithms/src/PrForwardTracking.cpp +++ b/Pr/PrAlgorithms/src/PrForwardTracking.cpp @@ -45,9 +45,9 @@ // local #include "HoughSearch.h" +#include "PrKernel/FTGeometryCache.h" #include "PrKernel/IPrAddUTHitsTool.h" #include "PrKernel/IPrDebugTrackingTool.h" -#include "PrKernel/PrFTZoneHandler.h" #include "PrTrackModel.h" // NN for ghostprob @@ -185,7 +185,6 @@ namespace LHCb::Pr::Forward { using TracksVP = Velo::Tracks; using scalar = SIMDWrapper::scalar::types; using simd = SIMDWrapper::best::types; - using ZoneCache = FTZoneCache::ZoneCache; using PrForwardTracks = std::vector<PrForwardTrack>; using StereoSearch = ::Hough::HoughSearch<int, Hough::DEPTH, Hough::MAXCAND, Detector::FT::nUVLayersTotal, Hough::NBINS>; @@ -330,18 +329,12 @@ namespace LHCb::Pr::Forward { /** * @brief initialises parameters of y projection parameterisation */ - void initYFitParameters( XCandidate& trackCandidate, const VeloSeedExtended& veloSeed, const ZoneCache& cache ) { + 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}}; - auto& betterDz = trackCandidate.getBetterDz(); - for ( auto zone : Detector::FT::xZonesLower ) { - zone += veloSeed.upperHalfTrack; - const auto dz = cache.zoneZPos[zone] - zReference; - betterDz[zone / 4u] = dz + cache.zoneDzDy[zone] * trackCandidate.y( dz ); - } } /** @@ -623,9 +616,10 @@ namespace LHCb::Pr::Forward { while ( fit ) { float s0{0.f}, sz{0.f}, sz2{0.f}, sd{0.f}, sdz{0.f}; for ( auto iHit : coordsToFit ) { - const auto dz = veloSeed.betterZ[SciFiHits.planeCode( iHit )] - zReference; - const auto d = SciFiHits.x( iHit ) - protoCand.x( dz ); - const auto w = SciFiHits.w( iHit ); + 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; @@ -643,8 +637,9 @@ namespace LHCb::Pr::Forward { auto maxChi2{0.f}; const bool notMultiple = protoCand.nDifferentPlanes() == coordsToFit.size(); for ( auto itH = coordsToFit.begin(); itH != itEnd; ++itH ) { - const auto d = - SciFiHits.x( *itH ) - protoCand.x( veloSeed.betterZ[SciFiHits.planeCode( *itH )] - zReference ); + 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; @@ -779,15 +774,14 @@ namespace LHCb::Pr::Forward { auto fitXProjection( XCandidate& track, const PrParametersX& pars, const FT::Hits& SciFiHits ) { auto& coordsToFit = track.getCoordsToFit(); const auto minHits = pars.minXHits; - - bool fit{true}; + 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 ); + const auto dzNoScale = track.calcBetterDz( iHit ); const auto d = track.distanceXHit( iHit, dzNoScale ); const auto w = SciFiHits.w( iHit ); const auto dz = .001f * dzNoScale; @@ -821,7 +815,9 @@ namespace LHCb::Pr::Forward { const auto itEnd = coordsToFit.end(); auto worst = itEnd; for ( auto itH = coordsToFit.begin(); itH != itEnd; ++itH ) { - const auto chi2 = track.chi2XHit( *itH, track.getBetterDz( *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; @@ -973,12 +969,14 @@ namespace LHCb::Pr::Forward { const auto& uvZones = veloSeed.upperHalfTrack ? Detector::FT::uvZonesUpper : Detector::FT::uvZonesLower; for ( auto zoneNumber : uvZones ) { - const auto zZone = cache.zoneZPos[zoneNumber]; + 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.zoneDzDy[zoneNumber]; + const auto betterZ = zZone + yInZone * cache.dzdy( zoneNumber, side ); const auto xPred = track.x( betterZ - zReference ); - const auto dxDy = cache.zoneDxDy[zoneNumber]; + const auto dxDy = cache.dxdy( zoneNumber, side ); const auto xPredShifted = xPred - yInZone * dxDy; const auto dxTol = pars.tolY + @@ -988,11 +986,12 @@ namespace LHCb::Pr::Forward { 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 ); - const auto signedPred = xPredShifted * dxDySign; + 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 signedDx = SciFiHits.x( iUV ) * dxDySign - signedPred; + 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 ); } @@ -1200,11 +1199,12 @@ namespace LHCb::Pr::Forward { if ( track.isPlaneUsed( iPlane ) ) continue; const auto pc = Detector::FT::stereoLayers[iPlane]; const auto zoneNumber = 2 * pc + veloSeed.upperHalfTrack; - const auto zZone = cache.zoneZPos[zoneNumber]; - const auto betterZ = zZone + track.y( zZone - zReference ) * cache.zoneDzDy[zoneNumber]; + const auto side = track.x( cache.z( pc ) - zReference ) > 0.f ? ZoneCache::Side::A : ZoneCache::Side::C; + const auto zZone = cache.z( zoneNumber, side ); + const auto betterZ = zZone + track.y( zZone - zReference ) * cache.dzdy( zoneNumber, side ); const auto yInZone = track.y( betterZ - zReference ); const auto xPred = track.x( betterZ - zReference ); - const auto xPredShifted = xPred - yInZone * cache.zoneDxDy[zoneNumber]; + const auto xPredShifted = xPred - yInZone * cache.dxdy( zoneNumber, side ); // TODO: here and in collection step is it better to use xPredShifted? const auto dxTol = pars.tolY + pars.tolYSlope * ( std::abs( xPred - veloSeed.xStraightInZone[pc] ) + std::abs( yInZone ) ); @@ -1215,7 +1215,9 @@ namespace LHCb::Pr::Forward { auto bestChi2{pars.maxChi2StereoAdd}; auto best{0}; for ( auto iUV{iUVStart}; SciFiHits.x( iUV ) <= xMax; ++iUV ) { - const auto d = SciFiHits.x( iUV ) - xPredShifted; + const auto dz = SciFiHits.z( iUV ) + yInZone * SciFiHits.dzDy( iUV ) - zReference; + const auto predShifted = track.x( dz ) - yInZone * SciFiHits.dxDy( iUV ); + const auto d = SciFiHits.x( iUV ) - predShifted; if ( const auto chi2 = d * d * SciFiHits.w( iUV ); chi2 < bestChi2 ) { bestChi2 = chi2; best = iUV; @@ -1442,7 +1444,7 @@ namespace LHCb::Pr::Forward { for ( const auto iLayer : Detector::FT::stereoLayers ) { const auto zZone = veloSeed.betterZ[iLayer]; // dxdy is what makes stereo hits shifty - const auto dxDy = cache.zoneDxDy[2 * iLayer + veloSeed.upperHalfTrack]; + const auto dxDy = cache.dxdy( iLayer ); const auto xShiftAtY = veloSeed.seed.y( zZone ) * dxDy; const auto xStraightInZone = veloSeed.xStraightInZone[iLayer]; const auto step = dzInvRef * ( zZone - veloSeed.zMag ); @@ -1655,7 +1657,7 @@ namespace LHCb::Pr::Forward { : base_class_t( name, pSvcLocator, {typename base_class_t::KeyValue{"SciFiHits", ""}, typename base_class_t::KeyValue{"InputTracks", ""}, - typename base_class_t::KeyValue{"FTZoneCache", FTZoneCache::Location + name}, + typename base_class_t::KeyValue{"FTZoneCache", std::string{ZoneCache::Location} + name}, typename base_class_t::KeyValue{"AddUTHitsToolName", "PrAddUTHitsTool"}, typename base_class_t::KeyValue{"Magnet", LHCb::Det::Magnet::det_path}}, typename base_class_t::KeyValue{"OutputTracks", ""} ) @@ -1682,7 +1684,7 @@ namespace LHCb::Pr::Forward { template <bool secondLoop> void selectXCandidates( PrForwardTracks&, const VeloSeedExtended&, ModSciFiHits::ModPrSciFiHitsSOA&, - const FT::Hits&, const PrParametersX&, const ZoneCache& ) const; + const FT::Hits&, const PrParametersX& ) const; template <bool STORE_DATA = false> auto calcMVAInput( const VeloSeedExtended&, PrForwardTrack&, float ) const; @@ -1719,6 +1721,7 @@ namespace LHCb::Pr::Forward { Gaudi::Property<int> m_minPlanes{this, "MinPlanes", 4}; Gaudi::Property<int> m_minPlanesComplX{this, "MinPlanesComplX", 4}; Gaudi::Property<int> m_minPlanesComplUV{this, "MinPlanesComplUV", 4}; + Gaudi::Property<int> m_maxLinearYEndT{this, "MaxLinearYEndT", 2600 * Gaudi::Units::mm}; // first loop Hough Cluster search // unsigned to avoid warnings @@ -1888,6 +1891,15 @@ namespace LHCb::Pr::Forward { return std::tuple{endv_pos, endv_dir, nanMomentum}; } }(); + + // remove tracks that are so steep in y that they go out of acceptance + if ( const auto yStraightEndT = + endv_pos.y() + abs( endv_dir.y() ) * ( static_cast<float>( Z( State::Location::EndT ) ) - endv_pos.z() ); + yStraightEndT.cast() > m_maxLinearYEndT ) { + direct_debug( "yStraightEndT =", yStraightEndT, "not in acceptance ... skip" ); + continue; + } + const VeloSeed seed{endv_pos.x().cast(), endv_pos.y().cast(), endv_pos.z().cast(), endv_dir.x().cast(), endv_dir.y().cast(), qOverP, magScaleFactor}; @@ -1908,12 +1920,7 @@ namespace LHCb::Pr::Forward { if ( p < m_minP || pt < m_minPT ) continue; } } - // remove tracks that are so steep in y that they go out of acceptance - if ( const float yLastStation = seed.y( cache.zoneZPos.back() ); - !cache.lowerLastZone.isInsideY( yLastStation ) && !cache.upperLastZone.isInsideY( yLastStation ) ) { - direct_debug( "yLastStation =", yLastStation, "not in last layer ... skip" ); - continue; - } + ++acceptedInputTracks; const VeloSeedExtended veloSeed{iTrack, seed, cache}; @@ -1957,7 +1964,7 @@ namespace LHCb::Pr::Forward { }() ); candidates.clear(); - selectXCandidates<false>( candidates, veloSeed, allXHits, SciFiHits, pars, cache ); + selectXCandidates<false>( candidates, veloSeed, allXHits, SciFiHits, pars ); const auto candStart2nd = candidates.size(); xCandidates1Cnt += candStart2nd; auto [good, ok] = @@ -1968,7 +1975,7 @@ namespace LHCb::Pr::Forward { direct_debug( "Now second loop" ); direct_debug( "" ); const auto candStart2nd = candidates.size(); - selectXCandidates<true>( candidates, veloSeed, allXHits, SciFiHits, pars2ndLoop, cache ); + selectXCandidates<true>( candidates, veloSeed, allXHits, SciFiHits, pars2ndLoop ); xCandidates2Cnt += candidates.size(); auto candidates2ndLoop = LHCb::make_span( candidates.begin() + candStart2nd, candidates.end() ); const auto ok2nd = selectFullCandidates<true>( candidates2ndLoop, pars2ndLoop, cache, veloSeed, SciFiHits, @@ -2013,8 +2020,7 @@ namespace LHCb::Pr::Forward { template <bool secondLoop> void PrForwardTracking<TrackType>::selectXCandidates( PrForwardTracks& candidates, const VeloSeedExtended& veloSeed, ModSciFiHits::ModPrSciFiHitsSOA& allXHits, - const FT::Hits& SciFiHits, const PrParametersX& pars, - const ZoneCache& cache ) const { + const FT::Hits& SciFiHits, const PrParametersX& pars ) const { if constexpr ( secondLoop ) direct_debug( "------------selectXCandidates------------secondLoop" ); else @@ -2080,7 +2086,7 @@ namespace LHCb::Pr::Forward { direct_debug( "xParams after linear fit =", protoCand.getXParams() ); auto ok = protoCand.nDifferentPlanes() >= pars.minXHits; if ( ok ) { - initYFitParameters( protoCand, veloSeed, cache ); + initYFitParameters( protoCand, veloSeed ); ok = fitXProjection( protoCand, pars, SciFiHits ); if ( ok ) { if ( fillEmptyXPlanes<secondLoop>( protoCand, pars, veloSeed, SciFiHits, m_maxChi2XAddFull ) ) { diff --git a/Pr/PrAlgorithms/src/PrHybridSeeding.cpp b/Pr/PrAlgorithms/src/PrHybridSeeding.cpp index cd57cc9ee24cd460101a144e9aeb58dd84cae5e5..7c9574c0c2337d1c456433932adc434d52805d4d 100644 --- a/Pr/PrAlgorithms/src/PrHybridSeeding.cpp +++ b/Pr/PrAlgorithms/src/PrHybridSeeding.cpp @@ -32,8 +32,8 @@ #include "Magnet/DeMagnet.h" #include "Math/CholeskyDecomp.h" #include "PrHybridSeedTrack.h" +#include "PrKernel/FTGeometryCache.h" #include "PrKernel/PrFTHitHandler.h" -#include "PrKernel/PrFTZoneHandler.h" #include "PrTrackFitterXYZ.h" #include "PrTrackFitterXZ.h" #include "PrTrackFitterYZ.h" @@ -45,12 +45,13 @@ namespace LHCb::Pr { namespace { constexpr unsigned int nParts = 2; - using SmallDexes = std::vector<size_t>; - using Tracks = Seeding::Tracks; - using SeedTag = Seeding::Tag; - using HitIter = PrFTHitHandler<ModPrHit>::HitIter; - using HitRevIter = PrFTHitHandler<ModPrHit>::HitRevIter; - using FTZoneCache::ZoneCache; + static_assert( nParts <= 2 ); + using SmallDexes = std::vector<size_t>; + using Tracks = Seeding::Tracks; + using SeedTag = Seeding::Tag; + using HitIter = PrFTHitHandler<ModPrHit>::HitIter; + using HitRevIter = PrFTHitHandler<ModPrHit>::HitRevIter; + using ZoneCache = Detector::FT::Cache::GeometryCache; struct SearchWindow { HitIter begin; HitIter end; @@ -925,7 +926,7 @@ namespace LHCb::Pr { HybridSeeding::HybridSeeding( const std::string& name, ISvcLocator* pSvcLocator ) : Transformer( name, pSvcLocator, {KeyValue{"FTHitsLocation", PrFTInfo::SciFiHitsLocation}, - KeyValue{"ZoneCache", FTZoneCache::Location + name}, + KeyValue{"ZoneCache", std::string{ZoneCache::Location} + name}, KeyValue{"Magnet", LHCb::Det::Magnet::det_path}}, KeyValue{"OutputName", LHCb::TrackLocation::Seed} ) {} @@ -1118,9 +1119,11 @@ namespace LHCb::Pr { const ZoneCache& zoneCache, unsigned int part, const Pr::Hybrid::SeedTrackX& xProje, const ZoneLimitsUV& zoneLimits, SearchWindowsUV& searchWindows ) const noexcept { - auto calculate_factor = [zoneCache, shift = 0.f, - partSign = float( 1 - static_cast<int>( part ) * 2 )]( int layer ) -> float { - return partSign / ( ( zoneCache.zoneZPos[layer] - shift ) * zoneCache.zoneDxDy[layer] ); + auto calculate_factor = [&]( int zone ) -> float { + const auto layer = zone / Detector::FT::nHalfLayers; + const auto side = xProje.x( zoneCache.z( layer ) ) > 0.f ? ZoneCache::Side::A : ZoneCache::Side::C; + const auto partSign = part < 1 ? 1.f : -1.f; + return partSign / ( zoneCache.z( zone, side ) * zoneCache.dxdy( zone, side ) ); }; auto calculate_factors = [&]( auto... layers ) { return std::array{calculate_factor( layers )...}; }; @@ -1144,16 +1147,18 @@ namespace LHCb::Pr { template <int CASE> template <int UV> // 0 if U, 1 if V - void HybridSeeding::addStereo<CASE>::CollectLayerUV( const ZoneCache& zoneCache, unsigned int part, - unsigned int layer, const Pr::Hybrid::SeedTrackX& xProje, + void HybridSeeding::addStereo<CASE>::CollectLayerUV( const ZoneCache& zoneCache, unsigned int part, unsigned int zone, + const Pr::Hybrid::SeedTrackX& xProje, const ZoneLimitUV& zoneLimits, SearchWindowUV& searchWindows, HoughSearch& hough ) const noexcept { static_assert( UV == 0 || UV == 1, "UV must be 0 (U) or 1 (V)!" ); - float dxDy = zoneCache.zoneDxDy[layer]; // 0 if U, 1 if V - float xPred = xProje.xFromDz( zoneCache.zoneZPos[layer] - Pr::Hybrid::zReference ); + const auto side = + xProje.x( zoneCache.z( zone / Detector::FT::nHalfLayers ) ) > 0.f ? ZoneCache::Side::A : ZoneCache::Side::C; + const auto dxDy = zoneCache.dxdy( zone, side ); + const auto xPred = xProje.xFromDz( zoneCache.z( zone, side ) - Pr::Hybrid::zReference ); // Note: UV layers are 1, 3, 11, 13, 19, 21 - unsigned int iStation = layer / 8; - unsigned int iLayer = layer / 4; + unsigned int iStation = zone / ( Detector::FT::nLayers * Detector::FT::nHalfLayers ); + unsigned int iStereoLayer = zone / ( Detector::FT::nHalfLayers * 2 ); // only half of the layers are stereo // Part 0 std::array<float, 2> yMinMax{m_hybridSeeding.m_yMins[part][0], m_hybridSeeding.m_yMaxs[part][0]}; @@ -1165,7 +1170,7 @@ namespace LHCb::Pr { searchWindows[iStation][0].end != zoneLimits[iStation][0].end; ++searchWindows[iStation][0].end ) { if ( searchWindows[iStation][0].end->coord < xMinMax[0] ) continue; if ( searchWindows[iStation][0].end->coord > xMinMax[1] ) break; - AddUVHit( iLayer, searchWindows[iStation][0].end, xPred, hough ); + AddUVHit( iStereoLayer, searchWindows[iStation][0].end, xPred, hough ); } // Part 1 yMinMax = {m_hybridSeeding.m_yMins[part][1], m_hybridSeeding.m_yMaxs[part][1]}; @@ -1176,7 +1181,7 @@ namespace LHCb::Pr { searchWindows[iStation][1].end != zoneLimits[iStation][1].end; ++searchWindows[iStation][1].end ) { if ( searchWindows[iStation][1].end->coord < xMinMax[0] ) continue; if ( searchWindows[iStation][1].end->coord > xMinMax[1] ) break; - AddUVHit( iLayer, searchWindows[iStation][1].end, xPred, hough ); + AddUVHit( iStereoLayer, searchWindows[iStation][1].end, xPred, hough ); } } @@ -1773,23 +1778,23 @@ namespace LHCb::Pr { // Array[0] = first layer in T1 for 2-hit combo xZones.zones[0] = firstZoneId; - xZones.zLays[0] = zoneCache.zoneZPos[firstZoneId]; + xZones.zLays[0] = zoneCache.z( firstZoneId / 2 ); xZones.dzLays[0] = xZones.zLays[0] - Pr::Hybrid::zReference; xZones.dz2Lays[0] = xZones.dzLays[0] * xZones.dzLays[0] * ( 1.f + m_dRatio * xZones.dzLays[0] ); // Array[1] = last layer in T3 for 2-hit combo xZones.zones[1] = lastZoneId; - xZones.zLays[1] = zoneCache.zoneZPos[lastZoneId]; + xZones.zLays[1] = zoneCache.z( lastZoneId / 2 ); xZones.dzLays[1] = xZones.zLays[1] - Pr::Hybrid::zReference; xZones.dz2Lays[1] = xZones.dzLays[1] * xZones.dzLays[1] * ( 1.f + m_dRatio * xZones.dzLays[1] ); // Array[2] = T2-1st x-layers xZones.zones[2] = LHCb::Detector::FT::UpperZones::T2X1 - part; - xZones.zLays[2] = zoneCache.zoneZPos[LHCb::Detector::FT::UpperZones::T2X1 - part]; + xZones.zLays[2] = zoneCache.z( LHCb::Detector::FT::UpperZones::T2X1 / 2 ); xZones.dzLays[2] = xZones.zLays[2] - Pr::Hybrid::zReference; xZones.dz2Lays[2] = xZones.dzLays[2] * xZones.dzLays[2] * ( 1.f + m_dRatio * xZones.dzLays[2] ); // Array[3] = T2-2nd x-layers xZones.zones[3] = LHCb::Detector::FT::UpperZones::T2X2 - part; - xZones.zLays[3] = zoneCache.zoneZPos[LHCb::Detector::FT::UpperZones::T2X2 - part]; + xZones.zLays[3] = zoneCache.z( LHCb::Detector::FT::UpperZones::T2X2 / 2 ); xZones.dzLays[3] = xZones.zLays[3] - Pr::Hybrid::zReference; xZones.dz2Lays[3] = xZones.dzLays[3] * xZones.dzLays[3] * ( 1.f + m_dRatio * xZones.dzLays[3] ); @@ -1800,7 +1805,7 @@ namespace LHCb::Pr { xZoneId = xZoneId - part; if ( xZoneId != firstZoneId && xZoneId != lastZoneId ) { xZones.zones[i] = xZoneId; - xZones.zLays[i] = zoneCache.zoneZPos[xZoneId]; + xZones.zLays[i] = zoneCache.z( xZoneId / 2 ); xZones.dzLays[i] = xZones.zLays[i] - Pr::Hybrid::zReference; xZones.dz2Lays[i] = xZones.dzLays[i] * xZones.dzLays[i] * ( 1.f + m_dRatio * xZones.dzLays[i] ); ++i; diff --git a/Pr/PrAlgorithms/src/PrStoreSciFiHits.cpp b/Pr/PrAlgorithms/src/PrStoreSciFiHits.cpp index f7f3ff0b98f148c4cb0177e3c2521cb8d48d9a1a..701b2e3d84cbbe329c48682c24af7340940bcded 100644 --- a/Pr/PrAlgorithms/src/PrStoreSciFiHits.cpp +++ b/Pr/PrAlgorithms/src/PrStoreSciFiHits.cpp @@ -18,7 +18,6 @@ #include "LHCbAlgs/Transformer.h" #include "LHCbMath/bit_cast.h" #include "PrKernel/FTMatsCache.h" -#include "PrKernel/PrFTZoneHandler.h" #include <algorithm> #include <array> @@ -36,7 +35,7 @@ namespace LHCb::Pr::FT { using FTLiteClusters = FTLiteCluster::FTLiteClusters; - using MatsCache = FTMatsCache::MatsCache; + using MatsCache = Detector::FT::Cache::MatsCache; // TODO: get this from a tool? constexpr auto invClusRes2 = [] { @@ -54,7 +53,7 @@ namespace LHCb::Pr::FT { StoreHits( const std::string& name, ISvcLocator* pSvcLocator ) : Transformer( name, pSvcLocator, {KeyValue{"HitsLocation", FTLiteClusterLocation::Default}, - KeyValue{"FTMatsCache", FTZoneCache::MatLocation + name}}, + KeyValue{"FTMatsCache", std::string{MatsCache::Location} + name}}, KeyValue{"Output", PrFTInfo::SciFiHitsLocation} ) {} StatusCode initialize() override { diff --git a/Pr/PrAlgorithms/src/PrTrackModel.h b/Pr/PrAlgorithms/src/PrTrackModel.h index cea65bfbe64933bc39c4a6252349f46b155e1c32..d5fde98bd70da8cffcf23fd7c29ee4fbd9a51396 100644 --- a/Pr/PrAlgorithms/src/PrTrackModel.h +++ b/Pr/PrAlgorithms/src/PrTrackModel.h @@ -15,7 +15,7 @@ #include "Event/PrHits.h" #include "Event/SOACollection.h" #include "LHCbMath/SIMDWrapper.h" -#include "PrKernel/PrFTZoneHandler.h" +#include "PrKernel/FTGeometryCache.h" #include "boost/container/static_vector.hpp" #include <array> #include <bitset> @@ -116,7 +116,7 @@ namespace LHCb::Pr { namespace Forward { using boost::container::static_vector; - + using ZoneCache = Detector::FT::Cache::GeometryCache; inline constexpr auto nan = std::numeric_limits<float>::signaling_NaN(); inline constexpr auto nanMomentum = std::numeric_limits<float>::quiet_NaN(); @@ -288,7 +288,7 @@ namespace LHCb::Pr { * */ struct VeloSeedExtended { - VeloSeedExtended( int iTrack, const VeloSeed& veloseed, const FTZoneCache::ZoneCache& cache ) + VeloSeedExtended( int iTrack, const VeloSeed& veloseed, const ZoneCache& cache ) : veloSciFiMatch{veloseed.tx, veloseed.ty, veloseed.tx2, veloseed.ty2} , seed{veloseed} , xStraightAtRef{veloseed.x( zReference )} @@ -308,7 +308,8 @@ namespace LHCb::Pr { const auto ty3tx = ty3 * tx; const auto ty3tx2 = ty3 * tx2; - upperHalfTrack = yStraightAtRef > 0; + upperHalfTrack = yStraightAtRef > 0.f; + pointingSide = xStraightAtRef > 0.f ? ZoneCache::Side::A : ZoneCache::Side::C; zMagTerm = zMagnetParamsRef[2] * tx; cxTerm = cxParams[0] + cxParams[1] * tx + cxParams[2] * ty + cxParams[3] * tx2 + cxParams[4] * txty + cxParams[5] * ty2; @@ -331,10 +332,11 @@ namespace LHCb::Pr { yCorrTermAbsLayer[l] = yCorrParamsLayers[l][1] * ty + yCorrParamsLayers[l][3] * ty3 + yCorrParamsLayers[l][4] * tytx2 + yCorrParamsLayers[l][7] * ty3tx2; } - const auto iZone = 2 * iLayer + upperHalfTrack; - const auto yStraightInZone = seed.y( cache.zoneZPos[iZone] ); - betterZ[iLayer] = cache.zoneZPos[iZone] + yStraightInZone * cache.zoneDzDy[iZone]; - xStraightInZone[iLayer] = seed.x( betterZ[iLayer] ); + const auto zLayer = cache.z( iLayer ); + const auto yStraightInLayer = seed.y( zLayer ); + betterZ[iLayer] = zLayer + yStraightInLayer * cache.dzdy( iLayer ); + xStraightInZone[iLayer] = seed.x( betterZ[iLayer] ); + yStraightInZone[iLayer] = seed.y( betterZ[iLayer] ); } } @@ -460,6 +462,7 @@ namespace LHCb::Pr { std::array<float, Detector::FT::nLayersTotal> betterZ; std::array<float, Detector::FT::nLayersTotal> xStraightInZone; + std::array<float, Detector::FT::nLayersTotal> yStraightInZone; std::array<float, Detector::FT::nUVLayersTotal> yCorrTermLayer; std::array<float, Detector::FT::nUVLayersTotal> yCorrTermAbsLayer; VeloSciFiMatch<float> veloSciFiMatch; @@ -477,6 +480,7 @@ namespace LHCb::Pr { float zMag; float yStraightAtRef; int iTrack; + ZoneCache::Side pointingSide; bool upperHalfTrack; // param[0] + param[1]*tx^2 + param[2]*tx dSlope_fringe + param[3]*ty^2 + @@ -598,8 +602,6 @@ namespace LHCb::Pr { auto& nDifferentPlanes() { return m_nDifferentPlanes; } auto& getCoordsToFit() { return m_coordsToFit; } auto getCoordsToFit() const { return m_coordsToFit; } - auto& getBetterDz() { return m_betterDz; } - auto getBetterDz( int fulldex ) { return m_betterDz[m_SciFiHits.planeCode( fulldex ) / 2u]; } auto& getXParams() { return m_xParams; } auto getXParams() const { return m_xParams; } auto& getYParams() { return m_yParams; } @@ -727,7 +729,6 @@ namespace LHCb::Pr { alignas( 64 ) std::array<int, Detector::FT::nXLayersTotal * planeMulti> m_idxOnPlane{}; alignas( 64 ) std::array<int, Detector::FT::nXLayersTotal> m_planeSize{}; alignas( 64 ) std::array<float, 5> m_s{nan, nan, nan, nan, nan}; - alignas( 64 ) std::array<float, Detector::FT::nXLayersTotal> m_betterDz{}; FitCoords<int, Detector::FT::nXLayersTotal * planeMulti> m_coordsToFit{}; // unsigned to avoid warnings diff --git a/Pr/PrFitParams/src/PrParameterisationData.cpp b/Pr/PrFitParams/src/PrParameterisationData.cpp index 0d47b8f2e1539c628bd84be4381ce90e61a28963..b91a90bca9d45279bd8ea1a2498b115f607b0f7b 100644 --- a/Pr/PrFitParams/src/PrParameterisationData.cpp +++ b/Pr/PrFitParams/src/PrParameterisationData.cpp @@ -21,8 +21,8 @@ #include "LHCbAlgs/Consumer.h" #include "PrFitParams/IPrFitTool.h" +#include "PrKernel/FTGeometryCache.h" #include "PrKernel/IPrDebugTrackingTool.h" -#include "PrKernel/PrFTZoneHandler.h" #include "TrackInterfaces/ITrackExtrapolator.h" #include <algorithm> @@ -40,8 +40,8 @@ */ namespace { - using FTZoneCache::ZoneCache; using namespace LHCb::Detector::FT; + using ZoneCache = Cache::GeometryCache; } // namespace namespace LHCb::Pr { @@ -234,11 +234,8 @@ namespace LHCb::Pr { std::array<std::pair<std::string, double>, 5 * LHCb::Detector::FT::nLayersTotal> layer_states; std::vector<Gaudi::XYZPoint> extrapolated_positions( LHCb::Detector::FT::nLayersTotal ); if ( [&] { - // the zone z positions are taken according to the half the track passes through - const auto upper_half_track = y_ref > 0; for ( unsigned int iLayer{0}; iLayer < LHCb::Detector::FT::nLayersTotal; ++iLayer ) { - const auto iZone = 2 * iLayer + upper_half_track; - const auto zZone = ft_cache.zoneZPos[iZone]; + const auto zZone = ft_cache.z( iLayer ); if ( m_extrap->propagate( track_state, zZone, *lhcb.geometry() ).isFailure() ) { return true; } layer_states[iLayer * 5] = {"x_l" + std::to_string( iLayer ), track_state.x()}; layer_states[iLayer * 5 + 1] = {"y_l" + std::to_string( iLayer ), track_state.y()}; diff --git a/Pr/PrKernel/include/PrKernel/FTGeometryCache.h b/Pr/PrKernel/include/PrKernel/FTGeometryCache.h new file mode 100644 index 0000000000000000000000000000000000000000..2828156876ccbc1926286bba2a12259638274bce --- /dev/null +++ b/Pr/PrKernel/include/PrKernel/FTGeometryCache.h @@ -0,0 +1,145 @@ +/*****************************************************************************\ +* (c) Copyright 2023 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. * +\*****************************************************************************/ +#pragma once + +#include "FTDAQ/FTInfo.h" +#include "FTDet/DeFTDetector.h" +#include <array> +#include <fmt/ostream.h> +#include <limits> +#include <string_view> + +namespace LHCb::Detector::FT::Cache { + + struct Plane { + float z{}; + float dxdy{}; + float dzdy{}; + }; + + /** + * @brief This object caches geometry information of the SciFi tracker as needed by the + * pattern recognition algorithms. + * @details The GeometryCache stores two kind of information: z, dxdy, dzdy of SciFi layers + * which comes directly from the geometry and is not subject to alignment, and the same + * quantities for SciFi quarters which are derived as the average of the respective module + * quantities and thus profit from alignment. + */ + struct GeometryCache { + enum class Side { A = 1, C = 0 }; + static constexpr std::string_view Location = "AlgorithmSpecific-FTGeometryCache"; + + GeometryCache(){}; + GeometryCache( const DeFT& ftDet ) { +#ifndef USE_DD4HEP + auto iLayer{0}; +#endif + ftDet.applyToAllLayers( [&]( const DeFTLayer& layer ) { +#ifdef USE_DD4HEP + const auto index = layer.layerIdx(); +#else + const auto index = iLayer++; +#endif + m_layers.at( index ) = {layer.globalZ(), layer.dxdy(), layer.dzdy()}; + } ); +#ifndef USE_DD4HEP + auto iQuarter{0}; +#endif + ftDet.applyToAllQuarters( [&]( const DeFTQuarter& quarter ) { +#ifdef USE_DD4HEP + const auto index = quarter.quarterIdx(); +#else + const auto index = iQuarter++; +#endif + m_quarters.at( index ) = {quarter.meanModuleZ(), quarter.meanModuleDxdy(), quarter.meanModuleDzdy()}; + } ); + }; + + /** + * @brief Get z position of A or C side of a zone (i.e. a quarter). + * + * @param iZone Index of the zone. + * @param side A or C side of the detector. + * @return This is the z position of all quarter modules averaged, i.e. changes with alignment. + */ + auto z( int iZone, Side side ) const { return m_quarters[quarter( iZone, side )].z; } + /** + * @brief Get z position of a layer. + * + * @param iLayer Index of the layer. + * @return This is the z position of the layer as defined in the geometry. + */ + auto z( int iLayer ) const { return m_layers[iLayer].z; } + /** + * @brief Get dzdy slope of A or C side of a zone (i.e. a quarter). + * + * @param iZone Index of the zone. + * @param side A or C side of the detector. + * @return This is the dzdy slope of all quarter modules averaged, i.e. changes with alignment. + */ + auto dzdy( int iZone, Side side ) const { return m_quarters[quarter( iZone, side )].dzdy; } + /** + * @brief Get dzdy slope of a layer. + * + * @param iLayer Index of the layer. + * @return This is the dzdy slope of the layer as defined in the geometry. + */ + auto dzdy( int iLayer ) const { return m_layers[iLayer].dzdy; } + /** + * @brief Get dxdy slope of A or C side of a zone (i.e. a quarter). + * + * @param iZone Index of the zone. + * @param side A or C side of the detector. + * @return This is the dxdy slope of all quarter modules averaged, i.e. changes with alignment. + */ + auto dxdy( int iZone, Side side ) const { return m_quarters[quarter( iZone, side )].dxdy; } + /** + * @brief Get dxdy slope of a layer. + * + * @param iLayer Index of the layer. + * @return This is the dxdy slope of the layer as defined in the geometry. + */ + auto dxdy( int iLayer ) const { return m_layers[iLayer].dxdy; } + friend std::ostream& operator<<( std::ostream& os, const GeometryCache& c ); + + private: + int quarter( int iZone, Side side ) const { return 2 * iZone + static_cast<int>( side ); } + std::array<Plane, nLayersTotal> m_layers{{{std::numeric_limits<float>::signaling_NaN()}}}; + std::array<Plane, nQuartersTotal> m_quarters{{{std::numeric_limits<float>::signaling_NaN()}}}; + }; + + inline std::ostream& operator<<( std::ostream& os, const GeometryCache& c ) { + fmt::print( os, "LHCb::Detector::FT::Cache::GeometryCache:\n" ); + for ( auto layer{0u}; layer < nLayersTotal; ++layer ) { + fmt::print( os, "\n=================Layer {:2}==================\n", layer ); + fmt::print( os, "(z = {:6.1f}, dxdy = {:7.5f}, dzdy = {:7.5f})\n", c.m_layers[layer].z, c.m_layers[layer].dxdy, + c.m_layers[layer].dzdy ); + fmt::print( os, "\n--------------------------------------------\n" ); + const auto q = nQuarters * layer; + fmt::print( os, + "z = {:11.5f} | z = {:11.5f}\ndxdy = {:11.5f} | dxdy = {:11.5f}\ndzdy = {:11.5f} | " + " dzdy = " + "{:11.5f}", + c.m_quarters[q + 3].z, c.m_quarters[q + 2].z, c.m_quarters[q + 3].dxdy, c.m_quarters[q + 2].dxdy, + c.m_quarters[q + 3].dzdy, c.m_quarters[q + 2].dzdy ); + fmt::print( os, "\n--------------------------------------------\n" ); + fmt::print( + os, + "z = {:11.5f} | z = {:11.5f}\ndxdy = {:11.5f} | dxdy = {:11.5f}\ndzdy = {:11.5f} | dzdy " + "= {:11.5f}", + c.m_quarters[q + 1].z, c.m_quarters[q + 0].z, c.m_quarters[q + 1].dxdy, c.m_quarters[q + 0].dxdy, + c.m_quarters[q + 1].dzdy, c.m_quarters[q + 0].dzdy ); + fmt::print( os, "\n--------------------------------------------\n" ); + fmt::print( os, "\n============================================\n" ); + } + return os; + } +} // namespace LHCb::Detector::FT::Cache \ No newline at end of file diff --git a/Pr/PrKernel/include/PrKernel/FTMatsCache.h b/Pr/PrKernel/include/PrKernel/FTMatsCache.h index 12f32f62f4b327bf73cf4ff740f4b44f81278cd7..06b7ae733595668012273fd77928557690a5d299 100644 --- a/Pr/PrKernel/include/PrKernel/FTMatsCache.h +++ b/Pr/PrKernel/include/PrKernel/FTMatsCache.h @@ -14,27 +14,23 @@ #include "FTDet/DeFTDetector.h" #include <array> -namespace FTMatsCache { - - const std::string Location = "AlgorithmSpecific-FTMatsCache"; - +namespace LHCb::Detector::FT::Cache { struct MatsCache { + static constexpr std::string_view Location = "AlgorithmSpecific-FTMatsCache"; /** * partial SoA cache for mats, reserve enough (here 4096 which is more than enough) * space for all mats ( all mats should be less than 2 * 8 mats * 12 modules * 12 layers) */ - std::array<float, LHCb::Detector::FT::maxNumberMats> dxdy{}; - std::array<float, LHCb::Detector::FT::maxNumberMats> dzdy{}; - std::array<float, LHCb::Detector::FT::maxNumberMats> globaldy{}; - std::array<ROOT::Math::XYZPointF, LHCb::Detector::FT::maxNumberMats> mirrorPoint{}; - std::array<ROOT::Math::XYZVectorF, LHCb::Detector::FT::maxNumberMats> ddx{}; - - float uBegin{}; - float halfChannelPitch{}; - float dieGap{}; - float sipmPitch{}; - - std::array<std::vector<double>, LHCb::Detector::FT::maxNumberMats> matContractionParameterVector{}; + std::array<float, maxNumberMats> dxdy{}; + std::array<float, maxNumberMats> dzdy{}; + std::array<float, maxNumberMats> globaldy{}; + std::array<ROOT::Math::XYZPointF, maxNumberMats> mirrorPoint{}; + std::array<ROOT::Math::XYZVectorF, maxNumberMats> ddx{}; + std::array<std::vector<double>, maxNumberMats> matContractionParameterVector{}; + float uBegin{}; + float halfChannelPitch{}; + float dieGap{}; + float sipmPitch{}; MatsCache(){}; // Needed for DD4HEP MatsCache( const DeFT& ftDet ) { @@ -68,4 +64,4 @@ namespace FTMatsCache { ftDet.applyToAllMats( func ); }; }; -} // namespace FTMatsCache +} // namespace LHCb::Detector::FT::Cache diff --git a/Pr/PrKernel/include/PrKernel/PrFTZoneHandler.h b/Pr/PrKernel/include/PrKernel/PrFTZoneHandler.h deleted file mode 100644 index aa88005173e03f09bad9a275de245486ef935a01..0000000000000000000000000000000000000000 --- a/Pr/PrKernel/include/PrKernel/PrFTZoneHandler.h +++ /dev/null @@ -1,162 +0,0 @@ -/*****************************************************************************\ -* (c) Copyright 2000-2018 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. * -\*****************************************************************************/ -#pragma once - -// Include files -#include "FTDAQ/FTInfo.h" -#include "FTDet/DeFTDetector.h" -#include "Kernel/DetectorSegment.h" -#include "PrKernel/PrHitZone.h" -#include <cassert> -#include <stdexcept> - -namespace FTZoneCache { - - const std::string Location = "AlgorithmSpecific-FTZoneCache"; - const std::string MatLocation = "AlgorithmSpecific-FTMatsCache"; - - /** @class PrFTZoneHandler PrFTZoneHandler.h - * Handlers of zones, the object is stored in the detector event store as a condition - * and each algorithms reads the object from there without calling the HitManagers (tools) - * @author Renato Quagliani - * @author Sebastien Ponce - */ - - class PrFTZoneHandler final { - - public: - PrFTZoneHandler( DeFT const& ftDet ) { -#ifdef USE_DD4HEP - auto func = [this]( const DeFTLayer& layer ) { - const auto id = layer.layerIdx(); - // fixme - DetectorSegment seg( 0, layer.globalZ(), layer.dxdy(), layer.dzdy(), 0., 0. ); - const auto xmax = 0.5f * layer.sizeX(); - const auto ymax = 0.5f * layer.sizeY(); - // The setGeometry defines the z at y=0, the dxDy and the dzDy, as well as the isX properties of the zone. - // This is important, since these are used in the following. - // They are set once for each zone in this method. - const bool ok = - ( this->MakeZone( 2 * id + 1, seg, -xmax, xmax, -25.f, ymax ) && // Small overlap (25 mm) for stereo layers - this->MakeZone( 2 * id, seg, -xmax, xmax, -ymax, 25.f ) ); // Small overlap (25 mm) for stereo layers - if ( !ok ) { throw std::runtime_error( "Failed to create DeFT Zones for ID = " + std::to_string( id ) ); } - }; - ftDet.applyToAllLayers( func ); -#else - for ( auto station : ftDet.stations() ) { - for ( auto layer : station->layers() ) { - const auto id = 4 * ( station->stationID() - 1 ) + layer->layerID(); - DetectorSegment seg( 0, layer->globalZ(), layer->dxdy(), layer->dzdy(), 0., 0. ); - const auto xmax = 0.5f * layer->sizeX(); - const auto ymax = 0.5f * layer->sizeY(); - // The setGeometry defines the z at y=0, the dxDy and the dzDy, as well as the isX properties of the zone. - // This is important, since these are used in the following. - // They are set once for each zone in this method. - const bool ok = - ( this->MakeZone( 2 * id + 1, seg, -xmax, xmax, -25.f, ymax ) && // Small overlap (25 mm) for stereo - // layers - this->MakeZone( 2 * id, seg, -xmax, xmax, -ymax, 25.f ) ); // Small overlap (25 mm) for stereo layers - if ( !ok ) { throw std::runtime_error( "Failed to create DeFT Zones for ID = " + std::to_string( id ) ); } - } - } -#endif - } - /// Standard constructor - PrFTZoneHandler() = default; - - bool MakeZone( unsigned int n, DetectorSegment& seg, float xMin, float xMax, float yMin, float yMax ) { - if ( n < m_zones.size() ) { - m_zones[n].setZone( n, seg, xMin, xMax, yMin, yMax ); - return true; - } else { - return false; - } - } - const PrHitZone& zone( unsigned int n ) const { - if ( n >= m_zones.size() ) { - throw std::runtime_error( "Zone index " + std::to_string( n ) + " is out-of-range" ); - } - return m_zones[n]; - } - - template <PrHitZone::Side SIDE> - static int getXZone( int layer ) { - if constexpr ( SIDE == PrHitZone::Side::Upper ) { - return LHCb::Detector::FT::xZonesUpper[layer]; - } else { - return LHCb::Detector::FT::xZonesLower[layer]; - } - } - - template <PrHitZone::Side SIDE> - static int getUVZone( int layer ) { - if constexpr ( SIDE == PrHitZone::Side::Upper ) { - return LHCb::Detector::FT::uvZonesUpper[layer]; - } else { - return LHCb::Detector::FT::uvZonesLower[layer]; - } - } - - template <PrHitZone::Side SIDE> - static int getTriangleZone( int layer ) { - if constexpr ( SIDE == PrHitZone::Side::Upper ) { - return LHCb::Detector::FT::uvZonesLower[layer]; - } else { - return LHCb::Detector::FT::uvZonesUpper[layer]; - } - } - - private: - // plain vector with indexing needed! - std::array<PrHitZone, LHCb::Detector::FT::nZonesTotal> m_zones; - }; - - /** - * @class ZoneCache PrFTZoneHandler.h - * @brief Caches derived conditions, e.g. positions of SciFi layers - * @author André Günther - */ - struct ZoneCache final { - PrFTZoneHandler handler; - PrHitZone lowerLastZone{}, upperLastZone{}; - std::array<float, LHCb::Detector::FT::nZonesTotal> zoneDxDy{std::numeric_limits<float>::signaling_NaN()}; - std::array<float, LHCb::Detector::FT::nZonesTotal> zoneDzDy{std::numeric_limits<float>::signaling_NaN()}; - std::array<float, LHCb::Detector::FT::nZonesTotal> zoneZPos{std::numeric_limits<float>::signaling_NaN()}; - ZoneCache( const DeFT& ftDet ) : handler( ftDet ) { - for ( unsigned int i{0}; i < LHCb::Detector::FT::nLayersTotal; ++i ) { - zoneZPos[2 * i] = handler.zone( 2 * i ).z(); - zoneZPos[2 * i + 1] = handler.zone( 2 * i + 1 ).z(); - zoneDxDy[2 * i] = handler.zone( 2 * i ).dxDy(); - zoneDxDy[2 * i + 1] = handler.zone( 2 * i + 1 ).dxDy(); - zoneDzDy[2 * i] = handler.zone( 2 * i ).dzDy(); - zoneDzDy[2 * i + 1] = handler.zone( 2 * i + 1 ).dzDy(); - } - lowerLastZone = handler.zone( LHCb::Detector::FT::nZonesTotal - 2 ); - upperLastZone = handler.zone( LHCb::Detector::FT::nZonesTotal - 1 ); - } - // for debugging - friend std::ostream& operator<<( std::ostream& os, const ZoneCache& c ) { - os << "ZoneCache:" << std::endl; - os << "zoneDxDy = [ "; - for ( auto z : c.zoneDxDy ) { os << z << " "; } - os << "]" << std::endl; - os << "zoneDzDy = [ "; - for ( auto z : c.zoneDzDy ) { os << z << " "; } - os << "]" << std::endl; - os << "zoneZPos = [ "; - for ( auto z : c.zoneZPos ) { os << z << " "; } - os << "]" << std::endl; - os << c.lowerLastZone << c.upperLastZone; - return os; - } - }; - -} // namespace FTZoneCache diff --git a/Pr/PrKernel/include/PrKernel/PrHitZone.h b/Pr/PrKernel/include/PrKernel/PrHitZone.h deleted file mode 100644 index d65b60a09d7e8fc854949b133bc3b0bd222487a2..0000000000000000000000000000000000000000 --- a/Pr/PrKernel/include/PrKernel/PrHitZone.h +++ /dev/null @@ -1,80 +0,0 @@ -/*****************************************************************************\ -* (c) Copyright 2000-2023 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. * -\*****************************************************************************/ - -#pragma once - -#include "Kernel/DetectorSegment.h" -#include <cmath> -#include <iostream> - -/** @class PrHitZone PrHitZone.h PrKernel/PrHitZone.h - * Store the information of a zone in the T stations - * A zone is a part of a layer with boundaries, so that all hits are measuring - * one coordinate in this zone. - * For FT, this is the top or bottom parts of each layer. - * The zone has a number, a plane code, slopes (dxDy, dzDy) and boundaries - * - * @author Olivier Callot - * @date 2012-03-13 - */ -class PrHitZone final { -public: - enum class Side { Lower, Upper }; - - PrHitZone() = default; - // All geometry informations casted here! - void setZone( unsigned int number, DetectorSegment& seg, float xMin, float xMax, float yMin, float yMax ) { - m_number = number; - m_planeCode = number / 2; - m_z = seg.z( 0.f ); - m_dxDy = ( seg.x( 100.f ) - seg.x( 0.f ) ) / 100.f; - m_dzDy = ( seg.z( 100.f ) - seg.z( 0.f ) ) / 100.f; - m_isX = std::abs( m_dxDy ) < 0.010f; - m_xMin = xMin; - m_xMax = xMax; - m_yMin = yMin; - m_yMax = yMax; - } - unsigned int number() const { return m_number; } - unsigned int planeCode() const { return m_planeCode; } - float z( float y = 0.f ) const { return m_z + m_dzDy * y; } - float dxDy() const { return m_dxDy; } - float dzDy() const { return m_dzDy; } - bool isX() const { return m_isX; } - - /// Is the point in the surrounding box? - bool isInside( float x, float y ) const { return y < m_yMax && y > m_yMin && x > m_xMin && x < m_xMax; } - - bool isInsideY( float y ) const { return y < m_yMax && y > m_yMin; } - - float dxOnAFibre() const { return ( m_yMax - m_yMin ) * m_dxDy; } - - // for debugging - friend std::ostream& operator<<( std::ostream& os, const PrHitZone& zone ) { - os << "PrHitZone " << zone.m_number << ": z=" << zone.m_z << ", dxDy=" << zone.m_dxDy << ", dzDy=" << zone.m_dzDy - << ", isX=" << zone.m_isX << std::endl; - os << "xEdges=[" << zone.m_xMin << ", " << zone.m_xMax << "] " - << "yEdges=[" << zone.m_yMin << ", " << zone.m_yMax << "]" << std::endl; - return os; - } - -private: - unsigned int m_number{}; - unsigned int m_planeCode{}; - float m_z{}; - float m_dxDy{}; - float m_dzDy{}; - bool m_isX{}; - float m_xMin{}; - float m_xMax{}; - float m_yMin{}; - float m_yMax{}; -}; diff --git a/Pr/SciFiTrackForwarding/src/SciFiTrackForwarding.cpp b/Pr/SciFiTrackForwarding/src/SciFiTrackForwarding.cpp index bcfa6089b6a5408d19bfb9b06497b0bbdbf86c18..3f1e43f9834a8f133744e6b4610fb0c69919d0d4 100644 --- a/Pr/SciFiTrackForwarding/src/SciFiTrackForwarding.cpp +++ b/Pr/SciFiTrackForwarding/src/SciFiTrackForwarding.cpp @@ -10,6 +10,7 @@ \*****************************************************************************/ #include "DetDesc/Condition.h" #include "DetDesc/GenericConditionAccessorHolder.h" +#include "Detector/FT/FTConstants.h" #include "LHCbAlgs/Transformer.h" #include "Event/PrLongTracks.h" @@ -18,7 +19,7 @@ #include "FTDet/DeFTDetector.h" #include "Magnet/DeMagnet.h" -#include "PrKernel/PrFTZoneHandler.h" +#include "PrKernel/FTGeometryCache.h" #include "PrKernel/PrSelection.h" #include "SciFiTrackForwardingHits.h" #include "vdt/log.h" @@ -115,29 +116,29 @@ namespace { // polynomial coefficients for extrapolations std::array<float, 48> ext_coeff{std::numeric_limits<float>::signaling_NaN()}; - FTZoneCache::PrFTZoneHandler zoneHandler; - GeomCache() = default; - GeomCache( DeFT const& ftDet ) : zoneHandler( ftDet ) { + GeomCache( DeFT const& ftDet ) { + const auto fullCache = LHCb::Detector::FT::Cache::GeometryCache( ftDet ); // get the layer z-positions and slopes to cache them for later usage - for ( int i{0}; i < 12; ++i ) { + for ( auto i{0u}; i < LHCb::Detector::FT::nLayersTotal; ++i ) { // the first 12 values are the bottom layers - LayerZPos[i] = zoneHandler.zone( 2 * i ).z(); - dZdYFiber[i] = zoneHandler.zone( 2 * i ).dzDy(); - + LayerZPos[i] = fullCache.z( i ); + dZdYFiber[i] = fullCache.dzdy( i ); + // NOTE + FIXME: these are the values from the geometry and thus not subject to alignment! + // if this algorithm is developed further at some point this should be fixed // the last 12 values are the top layers - LayerZPos[12 + i] = zoneHandler.zone( 2 * i + 1 ).z(); - dZdYFiber[12 + i] = zoneHandler.zone( 2 * i + 1 ).dzDy(); + LayerZPos[12 + i] = fullCache.z( i ); + dZdYFiber[12 + i] = fullCache.dzdy( i ); } // get the uv layer fiber slopes to cache them for later usage for ( int i{0}; i < 6; ++i ) { // the first 6 values are the bottom uv layers - dXdYFiber[i] = -zoneHandler.zone( LHCb::Detector::FT::stereoZones[2 * i] ).dxDy(); + dXdYFiber[i] = -fullCache.dxdy( LHCb::Detector::FT::stereoZones[2 * i] / 2 ); // the last 6 values are the top uv layers - dXdYFiber[6 + i] = -zoneHandler.zone( LHCb::Detector::FT::stereoZones[2 * i + 1] ).dxDy(); + dXdYFiber[6 + i] = -fullCache.dxdy( LHCb::Detector::FT::stereoZones[2 * i + 1] / 2 ); } // calculate some coefficients for the extrapolation into other stations and put them in a cache