diff --git a/Pr/PrConverters/src/VPMicroClustersToVPLightClustersConverter.cpp b/Pr/PrConverters/src/VPMicroClustersToVPLightClustersConverter.cpp index 930db47fc36546ebf5b3d0ed29633a40b23b2010..0549dc4d725ed9a0845dcd95eb9d103269dff398 100644 --- a/Pr/PrConverters/src/VPMicroClustersToVPLightClustersConverter.cpp +++ b/Pr/PrConverters/src/VPMicroClustersToVPLightClustersConverter.cpp @@ -52,21 +52,9 @@ LHCb::VPLightClusters LHCb::VPMicroClustersToVPLightClustersConverter::operator( clusters.reserve( dst.size() ); for ( const auto& vp_micro_cluster : dst.range() ) { - const auto& vpID = vp_micro_cluster.channelID(); - - const float fx = vp_micro_cluster.xfraction() / 255.; - const float fy = vp_micro_cluster.yfraction() / 255.; - const uint32_t cx = vpID.scol(); - const float local_x = devp.local_x( cx ) + fx * devp.x_pitch( cx ); - const float local_y = ( 0.5f + fy + to_unsigned( vpID.row() ) ) * devp.pixel_size(); - - const auto& ltg = devp.ltg( vpID.sensor() ); - - 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] ); - - clusters.emplace_back( vp_micro_cluster.xfraction(), vp_micro_cluster.yfraction(), gx, gy, gz, vpID ); + const auto pos = vp_micro_cluster.globalPosition( devp ); + clusters.emplace_back( vp_micro_cluster.xfraction(), vp_micro_cluster.yfraction(), pos.x(), pos.y(), pos.z(), + vp_micro_cluster.channelID() ); } sortClusterContainer( clusters ); diff --git a/Pr/PrKernel/include/PrKernel/FTMatsCache.h b/Pr/PrKernel/include/PrKernel/FTMatsCache.h index 6e8d2deb7bc502a4c754db02aef52c698c46d0a9..0ec4ff3dcf25856d822039652cc8a416e4ac90d4 100644 --- a/Pr/PrKernel/include/PrKernel/FTMatsCache.h +++ b/Pr/PrKernel/include/PrKernel/FTMatsCache.h @@ -59,7 +59,7 @@ namespace LHCb::Detector::FT::Cache { this->uBegin = first_mat->uBegin(); this->halfChannelPitch = first_mat->halfChannelPitch(); #endif - auto func = [this, applyMatContractionCalibration]( const DeFTMat& mat ) { + auto func = [this, applyMatContractionCalibration, parent]( const DeFTMat& mat ) { assert( LHCb::essentiallyEqual( this->dieGap, mat.dieGap() ) && "Unexpected difference in dieGap" ); assert( LHCb::essentiallyEqual( this->sipmPitch, mat.sipmPitch() ) && "Unexpected difference in sipmPitch" ); assert( LHCb::essentiallyEqual( this->uBegin, mat.uBegin() ) && "Unexpected difference in uBegin" ); @@ -79,16 +79,15 @@ namespace LHCb::Detector::FT::Cache { if ( applyMatContractionCalibration ) { // compute all the calibrated channel positions and ensure always x_i-1 // < x_i - // check the readout direction - const bool normalReadoutDir = matId.normalReadoutDir(); - float prevX = ( normalReadoutDir ? -1 : 1 ) * 10e3; - - bool alwaysIncreasing = true; + // check the readout direction + const int readoutdir = matId.normalReadoutDir() ? 1 : -1; + bool alwaysIncreasing = true; // loop over SiPMs for ( size_t sipm = 0; sipm < FTConstants::nSiPM; sipm++ ) { // loop over channels - for ( size_t chan = 1; chan < FTConstants::nChannels; chan++ ) { + float prevX{ 0 }; + for ( size_t chan = 0; chan < FTConstants::nChannels; chan++ ) { const auto id = LHCb::Detector::FTChannelID( mat.stationID(), mat.layerID(), mat.quarterID(), mat.moduleID(), mat.matID(), sipm, chan ); @@ -97,21 +96,29 @@ namespace LHCb::Detector::FT::Cache { if ( !this->matContractionParameterVector[index].empty() ) { // is applyToAllMats SIMD? Is this safe? x0 = calculateCalibratedX( id, x0_orig ); } - - const bool increasingNormal = ( x0 > prevX ) && normalReadoutDir; - const bool decreasingReverse = ( x0 < prevX ) && !normalReadoutDir; - alwaysIncreasing = alwaysIncreasing && ( increasingNormal || decreasingReverse ); - prevX = x0; + if ( chan > 0 ) { + const bool increasing = ( x0 - prevX ) * readoutdir > 0; + alwaysIncreasing = alwaysIncreasing && increasing; + if ( !increasing ) { + parent->debug() << "Invalid channel ordering for mat: " + << " station=" << int( mat.stationID() ) << " layer=" << int( mat.layerID() ) + << " quarter=" << int( mat.quarterID() ) << " module=" << int( mat.moduleID() ) + << " mat=" << int( mat.matID() ) << " simp=" << sipm << " channel=" << chan + << " xprev=" << prevX << " x=" << x0 << " dir=" << readoutdir + << " delta: " << this->matContractionParameterVector[index][sipm * 128 + chan] << " " + << this->matContractionParameterVector[index][sipm * 128 + chan - 1] << " " << endmsg; + } + } + prevX = x0; } } this->nInvalid += !alwaysIncreasing; - - std::transform( this->matContractionParameterVector[index].begin(), - this->matContractionParameterVector[index].end(), - this->matContractionParameterVector[index].begin(), [alwaysIncreasing]( double value ) { - return alwaysIncreasing * value; // if the mat is not always increasing then - // all the values are overwritten with a 0 for that mat - } ); + std::transform( + this->matContractionParameterVector[index].begin(), this->matContractionParameterVector[index].end(), + this->matContractionParameterVector[index].begin(), [alwaysIncreasing]( double value ) { + return alwaysIncreasing ? value : 0.; // if the mat is not always increasing then + // all the values are overwritten with a 0 for that mat + } ); } }; ftDet.applyToAllMats( func ); diff --git a/Tr/PrKalmanFilter/include/PrKalmanFilter/KF.h b/Tr/PrKalmanFilter/include/PrKalmanFilter/KF.h index 9c3da4547b30981bf99d6428a2f7335a86239340..95f56c15c32128de4f44b07ea14d759f94f72b62 100644 --- a/Tr/PrKalmanFilter/include/PrKalmanFilter/KF.h +++ b/Tr/PrKalmanFilter/include/PrKalmanFilter/KF.h @@ -640,7 +640,7 @@ namespace LHCb::Pr::Tracks::Fit { for ( auto& node : fitnodes ) node.ref_vector() = node.final_state_vec; } - template + template void classical_smoother_iteration_for_alignment( LHCb::span fitnodes, LHCb::span gain_matrices ) { // classical smoother does NOT rely on any information from the backward pass. @@ -655,20 +655,14 @@ namespace LHCb::Pr::Tracks::Fit { // start at the last node in forward direction int idx = fitnodes.size() - 1; // we start smoothing at the first node that has upstream info + // note added March 2025: because we skip nodes with outliers, we may miss some gain matrices at the end of the + // track. these are not used for alignment so we leave it like this. do { auto& node = fitnodes[idx]; node.final_state_vec = node.filtered_state_vec[Node::forward]; node.final_state_cov = node.filtered_state_cov[Node::forward]; - if constexpr ( UpdateReference ) { - // FIXME should we set ref vector and update ref-residual? - // I think so right? best info for alignment would be the updated info? - node.ref_vector() = node.filtered_state_vec[Node::forward]; - // updates H and the ref residual based on the new ref_vector - node.project_reference(); - // - } --idx; - } while ( idx > 0 && !fitnodes[idx].isHitOnTrack() ); + } while ( !fitnodes[idx + 1].isHitOnTrack() && idx > 0 && !fitnodes[idx].isHitOnTrack() ); // all remaining nodes need to be smoothed for ( ; idx >= 0; --idx ) { @@ -694,13 +688,19 @@ namespace LHCb::Pr::Tracks::Fit { // FIXME throw std::runtime_error( "inversion in classical smoother failed" ); } - // calculate gain matrix A. we can make this quicker by epxloiting that F is empty - A = node.final_state_cov * Transpose( F ) * inv_nextnode_cov; - Gaudi::TrackSymMatrix FCFinv = LHCb::Math::Similarity( F, node.final_state_cov ); // is also nextNodeC - Q + A = node.final_state_cov * Transpose( F ) * inv_nextnode_cov; +#ifdef COMMONEXPRESSION + Gaudi::TrackSymMatrix delta = nextnode.final_state_cov - nextnode.predicted_state_cov[Node::forward]; + LHCb::Math::Similarity( A, delta, node.final_state_cov ); + node.final_state_cov += node.filtered_state_cov[Node::forward]; +#else + Gaudi::TrackSymMatrix FCFinv = + LHCb::Math::Similarity( F, node.filtered_state_cov[Node::forward] ); // is also nextNodeC - Q FCFinv.InvertChol(); Gaudi::TrackSymMatrix sum = nextnode.final_state_cov + Q + LHCb::Math::Similarity( Q, FCFinv ); LHCb::Math::Similarity( A, sum, node.final_state_cov ); +#endif } else { // if there is no noise, the gain matrix is just the inverse of // the transport matrix @@ -710,14 +710,15 @@ namespace LHCb::Pr::Tracks::Fit { LHCb::Math::Similarity( A, nextnode.final_state_cov, node.final_state_cov ); } node.final_state_vec += A * ( nextnode.final_state_vec - nextnode.predicted_state_vec[Node::forward] ); - if constexpr ( UpdateReference ) { + } + if constexpr ( UpdateReference ) + for ( auto& node : fitnodes ) { // FIXME should we set ref vector and update ref-residual? // I think so right? best info for alignment would be the updated info? node.ref_vector() = node.final_state_vec; // updates H and the ref residual based on the new ref_vector node.project_reference(); } - } } inline Gaudi::Matrix1x6 alignmentDerivatives( const Node& node, const Gaudi::XYZPoint& pivot ) { diff --git a/Tr/TrackMonitors/CMakeLists.txt b/Tr/TrackMonitors/CMakeLists.txt index 609ad38ea7c0b99004fc59d6bb16f6b639785109..9bbd9354ea695233d32d8e1c127fdfcd61361cf3 100644 --- a/Tr/TrackMonitors/CMakeLists.txt +++ b/Tr/TrackMonitors/CMakeLists.txt @@ -16,7 +16,6 @@ Tr/TrackMonitors gaudi_add_module(TrackMonitors SOURCES src/FTTrackMonitor.cpp - src/FTMatCalibrationMonitor.cpp src/FTHitEfficiencyMonitor.cpp src/FTPseudoHitEfficiencyMonitor.cpp src/PrSciFiHitsMonitor.cpp diff --git a/Tr/TrackMonitors/src/FTMatCalibrationMonitor.cpp b/Tr/TrackMonitors/src/FTMatCalibrationMonitor.cpp deleted file mode 100644 index b50f464a229096186ed0bd8a9ffa6e4696bc59e3..0000000000000000000000000000000000000000 --- a/Tr/TrackMonitors/src/FTMatCalibrationMonitor.cpp +++ /dev/null @@ -1,321 +0,0 @@ -/*****************************************************************************\ -* (c) Copyright 2024 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. * -\*****************************************************************************/ -#include "DetDesc/DetectorElement.h" -#include "DetDesc/GenericConditionAccessorHolder.h" -#include "Detector/FT/FTChannelID.h" -#include "Detector/FT/FTConstants.h" -#include "Event/FitNode.h" -#include "Event/KalmanFitResult.h" -#include "Event/Measurement.h" -#include "Event/PrFitNode.h" -#include "Event/PrKalmanFitResult.h" -#include "Event/Track.h" -#include "Event/TrackParameters.h" -#include "FTDet/DeFTDetector.h" -#include "FTDet/DeFTMat.h" -#include "GaudiAlg/GaudiTupleAlg.h" -#include "GaudiKernel/SystemOfUnits.h" -#include "GaudiKernel/Vector3DTypes.h" -#include "Kernel/LHCbID.h" -#include "LHCbAlgs/Consumer.h" -#include "Map.h" -#include "TrackKernel/TrackFunctors.h" -#include "TrackMonitorTupleBase.h" -#include -#include -#include - -namespace FTConstants = LHCb::Detector::FT; -// backward compatibility with Gaudi < v39 where StaticProfileHistogram is just the plain ProfileHistogram -#if GAUDI_MAJOR_VERSION < 39 -namespace Gaudi::Accumulators { - template - using StaticProfileHistogram = ProfileHistogram; -} // namespace Gaudi::Accumulators -#endif - -class FTMatCalibrationMonitor - : public LHCb::Algorithm::Consumer> { - -public: - FTMatCalibrationMonitor( const std::string& name, ISvcLocator* pSvcLocator ) - : Consumer( name, pSvcLocator, - { KeyValue{ "TracksInContainer", LHCb::TrackLocation::Default }, - KeyValue{ "FTDetectorLocation", DeFTDetectorLocation::Default } } ){}; - - StatusCode initialize() override { return Consumer::initialize(); } - - StatusCode finalize() override { - return Consumer::finalize().andThen( [&] { - if ( m_produceYamlFile ) { - YAML::Node yamlNode; - // function that subtracts the average from a sequence of optionals - auto subtractaverage = []( auto begin, auto end ) { - size_t n{ 0 }; - double sum{ 0 }; - std::for_each( begin, end, [&n, &sum]( const auto& i ) { - if ( i ) { - ++n; - sum += *i; - }; - } ); - if ( n > 0 ) { - const double av = sum / n; - std::for_each( begin, end, [&av]( auto& i ) { - if ( i ) { ( *i ) -= av; }; - } ); - } - }; - - // Loop over the stations, layers, quarters, modules and mats to make the names of the conditions as strings - constexpr auto nChannelsPerMat = FTConstants::nChannels * FTConstants::nSiPM; - constexpr auto nChannelsPerModule = nChannelsPerMat * FTConstants::nMats; - for ( const auto stationID : LHCb::Detector::FTChannelID::allStationIDs ) { - for ( const auto layerID : LHCb::Detector::FTChannelID::allLayerIDs ) { - for ( const auto quarterID : LHCb::Detector::FTChannelID::allQuarterIDs ) { - for ( const auto moduleID : LHCb::Detector::FTChannelID::allModuleIDsInStation( stationID ) ) { - // compute the deltas per module such that we can make a half-module average. - // use std::optional such that we can keep track of which values were set. - std::array, nChannelsPerModule> deltas; - deltas.fill( std::optional{} ); // does default initialization work? - - for ( const auto matID : LHCb::Detector::FTChannelID::allMatIDs ) { - LHCb::Detector::FTChannelID chanID_mat( stationID, layerID, quarterID, moduleID, matID, 0, 0 ); - const auto matoffset = nChannelsPerMat * to_unsigned( matID ); //---LoH: TODO: Detector!442 - // the local and global numbering can be different depending on the quarter - - // get the histogram for this mat - const auto& matHisto = m_histograms_mat_wrt0.at( chanID_mat.globalMatIdx() ); - - if ( m_calibmodel == Model::PerChannel ) { - for ( std::size_t ichan = 1; ichan < nChannelsPerMat; ++ichan ) { - const auto& bincontents = matHisto.UnbiasedresidualPerChannel.binValue( ichan + 1 ); - const float sumOfSqs = std::get<1>( bincontents ); - const float sum = std::get<1>( std::get<0>( bincontents ) ); - const auto nEntries = std::get<0>( std::get<0>( bincontents ) ); - const float val = sum / nEntries; - const auto variance = sumOfSqs / nEntries - val * val; - const auto error = sqrt( variance ); - // it is not necessary to use the variance: to measure an mean with a precision of 1/sqrt(N) of - // the variance of the parent distribution, you need at least N entries. so here we'd like N to - // be at least 100 or so. - if ( nEntries >= m_minEntries && val / error > m_significanceThreshold ) { - deltas[ichan + matoffset] = -1 * val; - // initialized[ichan + matoffset ] = true; - } - } - } else { - for ( unsigned iasic = 0; iasic < FTConstants::nSiPM; ++iasic ) { - // perform a fit to the model y = [0] + [1]*x - using Vector = ROOT::Math::SVector; - using SymMatrix = ROOT::Math::SMatrix>; - Vector halfDChi2DX; - SymMatrix halfD2Chi2DX2; - for ( unsigned ichan = 0; ichan < FTConstants::nChannels; ++ichan ) { - const auto ibin = iasic * FTConstants::nChannels + ichan + 1; - const auto& bincontents = matHisto.UnbiasedresidualPerChannel.binValue( ibin ); - const double x = ichan; - const double wy = std::get<1>( std::get<0>( bincontents ) ); - const auto w = std::get<0>( std::get<0>( bincontents ) ); - halfDChi2DX( 0 ) += -wy; - halfDChi2DX( 1 ) += -wy * x; - halfD2Chi2DX2( 0, 0 ) += w; - halfD2Chi2DX2( 0, 1 ) += w * x; - halfD2Chi2DX2( 1, 1 ) += w * x * x; - } - if ( halfD2Chi2DX2( 0, 0 ) > m_minEntries ) { - SymMatrix C = halfD2Chi2DX2; - const auto ok = C.Invert(); - if ( ok ) { - Vector par = -C * halfDChi2DX; - // compute the channel deltas from the model - for ( unsigned ichan = 0; ichan < FTConstants::nChannels; ++ichan ) - deltas[ichan + iasic * FTConstants::nChannels + matoffset] = -( par[0] + ichan * par[1] ); - } - } - } - } - // subtract the average per mat (only if we run separate alignment for mat Tx) - if ( m_subtractMatAverage ) - subtractaverage( deltas.begin() + matoffset, deltas.begin() + matoffset + nChannelsPerMat ); - } // end loop over mats - - // subtract the average per module - subtractaverage( deltas.begin(), deltas.end() ); - - // for each mat write the results to yaml - for ( const auto matID : LHCb::Detector::FTChannelID::allMatIDs ) { - const auto matoffset = nChannelsPerMat * to_unsigned( matID ); - const std::string conditionName = - fmt::format( "matContraction{matName}", - fmt::arg( "matName", LHCb::Detector::FTChannelID::matStringName( - stationID, layerID, quarterID, moduleID, matID ) ) ); - yamlNode["matContraction"][conditionName] = - YAML::Load( "[]" ); // is there a better way to make it appear flow style? - auto thisnode = yamlNode["matContraction"][conditionName]; - std::for_each( deltas.begin() + matoffset, deltas.begin() + matoffset + nChannelsPerMat, - [&thisnode]( const auto& d ) { thisnode.push_back( d ? *d : 0.0 ); } ); - } // end loop over mats - } // end loop over modules - } // end loop over quarters - } // end loop over layers - } // end loop over stations - std::ofstream fileOut( m_outputFileName ); - fileOut << yamlNode; - fileOut.close(); - } // end if - } ); - } - - void operator()( const LHCb::Track::Range& tracks, const DeFT& ftDet ) const override; - -private: - enum Model { PerChannel = 0, PerAsic = 1 }; - Gaudi::Property m_produceYamlFile{ this, "ProduceYamlFile", false }; - Gaudi::Property m_outputFileName{ this, "OutputFileName", "matContraction.yaml" }; - Gaudi::Property m_significanceThreshold{ this, "SignificanceThreshold", -100. }; - Gaudi::Property m_minEntries{ this, "MinEntriesPerChannel", 100 }; - Gaudi::Property m_calibmodel{ this, "Model", Model::PerAsic }; - Gaudi::Property m_subtractMatAverage{ this, "SubtractMatAverage", false }; - - template - void fillContainers( const TFitResult& track, const DeFT& ftDet ) const; - -private: - // histogram to look at the unbiased residual, considering the current calibration values - struct DifferencePerChannelHistogram { - mutable Gaudi::Accumulators::StaticProfileHistogram<1> UnbiasedresidualPerChannel; - - DifferencePerChannelHistogram( const FTMatCalibrationMonitor* owner, std::string const& mat ) - : UnbiasedresidualPerChannel{ - owner, - "DeltaUnbiasedresidualPerChannel" + mat, - "Difference in unbiased residual per Channel, " + mat, - { FTConstants::nSiPM * FTConstants::nChannels, 0, FTConstants::nSiPM * FTConstants::nChannels } } {} - }; - - // histogram to look at the unbiased residual, subtracting the current calibration values (i.e. residual wrt no - // calibration) - struct Ref0PerChannelHistogram { - mutable Gaudi::Accumulators::StaticProfileHistogram<1> UnbiasedresidualPerChannel; - - Ref0PerChannelHistogram( const FTMatCalibrationMonitor* owner, std::string const& mat ) - : UnbiasedresidualPerChannel{ - owner, - "Ref0UnbiasedresidualPerChannel" + mat, - "Unbiased residual per Channel wrt 0, " + mat, - { FTConstants::nSiPM * FTConstants::nChannels, 0, FTConstants::nSiPM * FTConstants::nChannels } } {} - }; - - std::map m_histograms_mat_wrt_current_cond = - []( const FTMatCalibrationMonitor* parent ) { - std::map map; - for ( size_t globalMatIdx = 0; globalMatIdx < LHCb::Detector::FT::nMatsTotal; globalMatIdx++ ) { - // todo: make naming consistent with the condition file - map.emplace( std::piecewise_construct, std::forward_as_tuple( globalMatIdx ), - std::forward_as_tuple( parent, "GlobalMat" + std::to_string( globalMatIdx ) ) ); - } - return map; - }( this ); - - std::map m_histograms_mat_wrt0 = []( const FTMatCalibrationMonitor* parent ) { - std::map map; - for ( size_t globalMatNumber = 0; globalMatNumber < LHCb::Detector::FT::nMatsTotal; globalMatNumber++ ) { - // todo: make naming consistent with the condition file - map.emplace( std::piecewise_construct, std::forward_as_tuple( globalMatNumber ), - std::forward_as_tuple( parent, "GlobalMat" + std::to_string( globalMatNumber ) ) ); - } - return map; - }( this ); -}; - -DECLARE_COMPONENT_WITH_ID( FTMatCalibrationMonitor, "FTMatCalibrationMonitor" ); - -void FTMatCalibrationMonitor::operator()( const LHCb::Track::Range& tracks, const DeFT& ftDet ) const { - - for ( const LHCb::Track* track : tracks ) { - const auto type = track->type(); - if ( !LHCb::Event::Enum::Track::hasT( type ) ) { continue; } - if ( track->checkFitStatus( LHCb::Track::FitStatus::Fitted ) ) { - const auto fr = track->fitResult(); - const auto prkfr = dynamic_cast( fr ); - if ( prkfr ) { - fillContainers( *prkfr, ftDet ); - } else { - const auto tmffr = dynamic_cast( fr ); - if ( tmffr ) { - fillContainers( *tmffr, ftDet ); - } else { - throw GaudiException( "Fit result is NULL - track was not fitted or uses wrong fit result type", name(), - StatusCode::FAILURE ); - } - } - } - } -} - -namespace { - Gaudi::XYZPoint mpos( const LHCb::FitNode& node ) { - const LHCb::Measurement& meas = node.measurement(); - LHCb::LineTraj mtraj = - meas.visit( [&]( const auto& m ) { return m.trajectory; }, - [&]( const LHCb::Measurement::VP2D& ) { - throw std::runtime_error( "Alignment derivatives not implemented for 2D measurements." ); - LHCb::LineTraj dummy; - return dummy; - } ); - return mtraj.beginPoint(); - } - Gaudi::XYZPoint mpos( const LHCb::Pr::Tracks::Fit::Node& node ) { - return Gaudi::XYZPoint{ node.measurement_pos[0], node.measurement_pos[1], node.measurement_pos[2] }; - } -} // namespace - -template -void FTMatCalibrationMonitor::fillContainers( const TFitResult& fitResult, const DeFT& ftDet ) const { - - for ( const auto& node : nodes( fitResult ) ) { - if ( !node.hasMeasurement() ) continue; - if ( node.isHitOnTrack() && node.isFT() ) { - LHCb::LHCbID lhcbID = id( node ); - assert( lhcbID.isFT() ); - - LHCb::Detector::FTChannelID chan = lhcbID.ftID(); - const auto globalMatIdx = chan.globalMatIdx(); - - // get the histograms for this mat - auto& histos_mats_wrt_current_cond = m_histograms_mat_wrt_current_cond.at( globalMatIdx ); - auto& histos_mats_wrt0 = m_histograms_mat_wrt0.at( globalMatIdx ); - const auto& ftMat = ftDet.findMat( chan ); - const auto ichan = - chan.channel() + chan.sipm() * FTConstants::nChannels; //---LoH: This is in local orientation (not CtoE) - - // get the unbiased track state - const auto trkstate = node.unbiasedState(); - // rotate position and direction to the local frame - const auto localtrkpos = ftMat->toLocal( trkstate.position() ); - const auto localtrkdir = ftMat->toLocal( trkstate.slopes() ); - // compute the local x coordinate at local z = 0 - const auto localtrkx = localtrkpos.x() - localtrkpos.z() * localtrkdir.x() / localtrkdir.z(); - // now do the same for the nodes measurement position: this should come out with z=0 - const auto localmpos = ftMat->toLocal( mpos( node ) ); - // the difference between the two is the residual that we need - const auto residual = localmpos.x() - localtrkx; - // You could check that the local z position is (very close) to zero: - // std::cout << "Local mpos: " << localmpos << std::endl ; - - histos_mats_wrt_current_cond.UnbiasedresidualPerChannel[ichan] += residual; - // double negative - histos_mats_wrt0.UnbiasedresidualPerChannel[ichan] += residual - ftMat->getmatContractionParameterVector()[ichan]; - } - } -} diff --git a/Tr/TrackMonitors/src/TrackFitMatchMonitor.cpp b/Tr/TrackMonitors/src/TrackFitMatchMonitor.cpp index 21eac5f44b9aa34efc77ead13a4c601ceb668429..accb9872167e22ffee045c76f0a629a2d25e34ff 100644 --- a/Tr/TrackMonitors/src/TrackFitMatchMonitor.cpp +++ b/Tr/TrackMonitors/src/TrackFitMatchMonitor.cpp @@ -89,8 +89,12 @@ private: mutable Gaudi::Accumulators::Histogram<1> m_curvatureRatioVeloUTToLongPullH1{ this, "curvatureRatioVeloUTToLongPull", "curvature ratio Velo-UT to Long pull", { 40, -5, 5 } }; mutable Gaudi::Accumulators::Histogram<1> m_kickZH1{ this, "kickZ", "Z position of magnet kick", { 40, 4900, 5400 } }; + mutable Gaudi::Accumulators::Histogram<1> m_kickZerrH1{ + this, "kickZerr", "error on Z position of magnet kick", { 40, 0, 50 } }; mutable Gaudi::Accumulators::ProfileHistogram<1> m_kickZVsXPr{ - this, "kickZVsXPr", "Z position of magnet kick versus x", { 40, -1500, 1500 } }; + this, "kickZVsXPr", "Z position of magnet kick versus x", { 40, -1500, 1500, "x [mm]" } }; + mutable Gaudi::Accumulators::ProfileHistogram<1> m_kickZVsQoPPr{ + this, "kickZVsQoPPr", "Z position of magnet kick versus q/p", { 50, -0.25, 0.25, "q/p [1/GeV]" } }; //============================================================================= // Define struct to store delta plot histograms @@ -293,6 +297,10 @@ void TrackFitMatchMonitor::makePlots( Trackers tracker, Gaudi::TrackVector delta histos->dty_qop_profile[state[4] * Gaudi::Units::GeV] += deltac( 3 ); histos->dty_tx_profile[state[2]] += deltac( 3 ); histos->dty_ty_profile[state[3]] += deltac( 3 ); + + histos->dty_pull_qop_profile[state[4] * Gaudi::Units::GeV] += deltacpull( 3 ); + histos->dty_pull_tx_profile[state[2]] += deltacpull( 3 ); + histos->dty_pull_ty_profile[state[3]] += deltacpull( 3 ); } } @@ -495,22 +503,46 @@ void TrackFitMatchMonitor::fill( const TFitResult& fitResult, const LHCb::Track& // compute the (x,z) point of the intersection of the 2 segments for linear propagation // FIXME: it must be better to take a fixed z position in T. if ( 1 / std::abs( qop ) > 5 * Gaudi::Units::GeV ) { + const LHCb::State stateVelo = + upstream ? filteredStateBackward( *lastVelo ) : filteredStateForward( *lastVelo ); + const auto xT = stateT.x(); const auto txT = stateT.tx(); const auto zT = stateT.z(); - const auto xVeloUT = stateVeloUT.x(); - const auto txVeloUT = stateVeloUT.tx(); - const auto zVeloUT = stateVeloUT.z(); - - const double zkick = ( zVeloUT * txVeloUT - xVeloUT + xT - zT * txT ) / ( txVeloUT - txT ); - const double xkick = xT + ( zkick - zT ) * txT; - // double xkickprime = stateVeloUT.x() + (zkick - stateVeloUT.z()) * stateVeloUT.tx() ; - ++m_kickZH1[zkick]; - // when rejecting reference nodes zkick distribution peaks around 5200mm - // and has a longer tail on left side, - // for profile plot remove outliers - if ( ( ZMagnetKickMean - 200 ) < zkick && zkick < ( ZMagnetKickMean + 200 ) ) m_kickZVsXPr[xkick] += zkick; + const auto xVelo = stateVelo.x(); + const auto txVelo = stateVelo.tx(); + const auto zVelo = stateVelo.z(); + + const double zkick = ( zVelo * txVelo - xVelo + xT - zT * txT ) / ( txVelo - txT ); + // we also want to have the error: formulas comes from the fast vertex fit + auto stateTatZkick = stateT; + stateTatZkick.linearTransportTo( zkick ); + auto stateVeloatZkick = stateVelo; + stateVeloatZkick.linearTransportTo( zkick ); + const Gaudi::SymMatrix2x2 cov = stateTatZkick.covariance().Sub( 0, 0 ) + + stateVeloatZkick.covariance().Sub( 0, 0 ); + Gaudi::SymMatrix2x2 W = cov; + const auto ok = W.InvertChol(); + if ( ok ) { + const Gaudi::Vector2 dtx = stateTatZkick.stateVector().Sub( 2 ) - + stateVeloatZkick.stateVector().Sub( 2 ); + const double Wpoca = ROOT::Math::Similarity( dtx, W ); + const double zkickerr = 1 / std::sqrt( Wpoca ); + ++m_kickZerrH1[zkickerr]; + if ( zkickerr < 50 * Gaudi::Units::mm ) { + const double xkick = xT + ( zkick - zT ) * txT; + // double xkickprime = stateVelo.x() + (zkick - stateVelo.z()) * stateVelo.tx() ; + ++m_kickZH1[zkick]; + // when rejecting reference nodes zkick distribution peaks around 5200mm + // and has a longer tail on left side, + // for profile plot remove outliers + if ( ( ZMagnetKickMean - 200 ) < zkick && zkick < ( ZMagnetKickMean + 200 ) ) { + m_kickZVsXPr[xkick] += zkick; + m_kickZVsQoPPr[stateT.qOverP() * Gaudi::Units::GeV] += zkick; + } + } + } } } } diff --git a/Tr/TrackMonitors/src/TrackMonitor.cpp b/Tr/TrackMonitors/src/TrackMonitor.cpp index 920c166d40d28febbc0721b5dd3106cb012aa4ab..497f6bc17efa53033745838007536720c168ec92 100644 --- a/Tr/TrackMonitors/src/TrackMonitor.cpp +++ b/Tr/TrackMonitors/src/TrackMonitor.cpp @@ -149,6 +149,8 @@ private: Gaudi::Accumulators::Histogram<1> nVPHits; Gaudi::Accumulators::Histogram<1> nFTHits; Gaudi::Accumulators::Histogram<1> NumOutliers; + Gaudi::Accumulators::Histogram<1> outlierType; + Gaudi::Accumulators::Histogram<1> outlierLayer; Gaudi::Accumulators::Histogram<1> nVeloALayers; Gaudi::Accumulators::Histogram<1> nVeloCLayers; Gaudi::Accumulators::Histogram<1> nVeloOverlapLayers; @@ -231,6 +233,8 @@ private: , nVPHits{ owner, prefix + "nVPHits", "# VP hits", { 27, -0.5, 26.5 } } , nFTHits{ owner, prefix + "nFTHits", "# FT hits", { 16, -0.5, 15.5 } } , NumOutliers{ owner, prefix + "noutliers", "#outliers", { 11, -0.5, 10.5 } } + , outlierType{ owner, prefix + "outliertype", "outlier type", { NHitTypes + 1, -0.5, int( NHitTypes ) + 0.5 } } + , outlierLayer{ owner, prefix + "outlierlayer", "outlier layer", { 51, -0.5, 50.5 } } , nVeloALayers{ owner, prefix + "nVeloAHitLayers", "# Velo-A hit layers", { 26, -0.5, 25.5 } } , nVeloCLayers{ owner, prefix + "nVeloCHitLayers", "# Velo-C hit layers", { 26, -0.5, 25.5 } } , nVeloOverlapLayers{ owner, prefix + "nVeloOverlapLayers", "# Velo overlap layers", { 26, -0.5, 25.5 } } @@ -491,15 +495,28 @@ void TrackMonitor::fillFitResultHistograms( const Track& track, const FitResultT // discard extremely small fraction of hits with zero error // on residual. (e.g. a downstream track with only one // active TT hit) - if ( node.isHitOnTrack() && node.errResidual2() > TrackParameters::lowTolerance && + if ( ( node.isHitOnTrack() || node.isOutlier() ) && node.errResidual2() > TrackParameters::lowTolerance && ( mtype = hittypemap( node ) ) != NHitTypes ) { - // factor for unbiasing the rms (not the mean!) - double f = std::sqrt( node.errMeasure2() / node.errResidual2() ); - auto& namedhistos = histos.hithistograms[mtype]; - ++namedhistos->Residual[f * node.residual()]; - ++namedhistos->ResidualPull[node.residual() / node.errResidual()]; - } else if ( node.isOutlier() ) { - ++numoutliers; + if ( node.isHitOnTrack() ) { + // factor for unbiasing the rms (not the mean!) + double f = std::sqrt( node.errMeasure2() / node.errResidual2() ); + auto& namedhistos = histos.hithistograms[mtype]; + ++namedhistos->Residual[f * node.residual()]; + ++namedhistos->ResidualPull[node.residual() / node.errResidual()]; + } else if ( node.isOutlier() ) { + ++numoutliers; + ++histos.outlierType[mtype]; + // We would like to know in which 'layer' the outlier is: construt a global layer from the LHCb ID + unsigned layer = 50; + const auto lhcbid = id( node ); + if ( lhcbid.isVP() ) + layer = lhcbid.vpID().station(); + else if ( lhcbid.isUT() ) + layer = lhcbid.utID().layer() + 30; + else if ( lhcbid.isFT() ) + layer = lhcbid.ftID().globalLayerIdx() + 38; + ++histos.outlierLayer[layer]; + } } } diff --git a/Tr/TrackMonitors/src/TrackParticleMonitor.cpp b/Tr/TrackMonitors/src/TrackParticleMonitor.cpp index d64c242d7e0661e38b42ef2a8209588543d4dafa..7a5eda088f90f9305ab8fce12fd4f26f1df30a72 100644 --- a/Tr/TrackMonitors/src/TrackParticleMonitor.cpp +++ b/Tr/TrackMonitors/src/TrackParticleMonitor.cpp @@ -13,14 +13,18 @@ #include "DetDesc/DetectorElement.h" #include "DetDesc/GenericConditionAccessorHolder.h" #include "Event/ChiSquare.h" +#include "Event/FunctorDefaults.h" #include "Event/Particle.h" #include "Event/ParticleUtils.h" +#include "Event/PrKalmanFitResult.h" +#include "Event/TrackFitResult.h" #include "Gaudi/Accumulators/Histogram.h" #include "GaudiKernel/PhysicalConstants.h" #include "GaudiKernel/ToolHandle.h" #include "Kernel/IParticlePropertySvc.h" #include "Kernel/ParticleProperty.h" #include "LHCbAlgs/Consumer.h" +#include "LHCbMath/StateVertexUtils.h" #include "Magnet/DeMagnet.h" #include "TrackInterfaces/ITrackStateProvider.h" #include "TrackKernel/TrackStateVertex.h" @@ -108,23 +112,43 @@ private: // Track and vertex histograms struct TrackVertexHistos { - mutable Gaudi::Accumulators::Histogram<1> track_chi2PerDoF; - mutable Gaudi::Accumulators::Histogram<1> track_pt; - mutable Gaudi::Accumulators::Histogram<1> track_p; - mutable Gaudi::Accumulators::Histogram<1> vertex_chi2; - mutable Gaudi::Accumulators::Histogram<1> vertex_chi2prob; - mutable Gaudi::Accumulators::Histogram<1> vertex_x_pos; - mutable Gaudi::Accumulators::Histogram<1> vertex_y_pos; - mutable Gaudi::Accumulators::Histogram<1> vertex_z_pos; - mutable Gaudi::Accumulators::Histogram<1> vertex_x_error; - mutable Gaudi::Accumulators::Histogram<1> vertex_y_error; - mutable Gaudi::Accumulators::Histogram<1> vertex_z_error; + mutable Gaudi::Accumulators::Histogram<1> track_chi2PerDoF; + mutable Gaudi::Accumulators::Histogram<1> track_pt; + mutable Gaudi::Accumulators::Histogram<1> track_p; + mutable Gaudi::Accumulators::Histogram<1> track_chi2Prob; + mutable Gaudi::Accumulators::Histogram<1> track_chi2ProbDownstream; + mutable Gaudi::Accumulators::Histogram<1> track_chi2ProbVelo; + mutable Gaudi::Accumulators::Histogram<1> track_chi2ProbMatch; + mutable Gaudi::Accumulators::ProfileHistogram<1> track_chi2ProbDownstreamVsPhi; + mutable Gaudi::Accumulators::Histogram<1> track_ipX; + mutable Gaudi::Accumulators::Histogram<1> track_ipY; + mutable Gaudi::Accumulators::Histogram<1> vertex_chi2; + mutable Gaudi::Accumulators::Histogram<1> vertex_chi2prob; + mutable Gaudi::Accumulators::Histogram<1> vertex_x_pos; + mutable Gaudi::Accumulators::Histogram<1> vertex_y_pos; + mutable Gaudi::Accumulators::Histogram<1> vertex_z_pos; + mutable Gaudi::Accumulators::Histogram<1> vertex_x_error; + mutable Gaudi::Accumulators::Histogram<1> vertex_y_error; + mutable Gaudi::Accumulators::Histogram<1> vertex_z_error; TrackVertexHistos( const TrackParticleMonitor* owner, float const& x_range, float const& y_range, float const& z_range ) : track_chi2PerDoF{ owner, "trackChi2", "Track #chi^2 per dof", { 100, 0, 5, "#chi^2 / nDoF" } } , track_pt{ owner, "trackPt", "Track p_{T}", { 100, 0, 20, "p_{T} [GeV]" } } , track_p{ owner, "trackP", "Track momentum", { 100, 0, 100, "p [GeV]" } } + , track_chi2Prob{ owner, "trackChi2Prob", "Track #chi^2 probability", { 50, 0, 1 } } + , track_chi2ProbDownstream{ owner, + "trackChi2ProbDownstream", + "Track downstream (T) #chi^2 probability", + { 50, 0, 1 } } + , track_chi2ProbVelo{ owner, "trackChi2ProbVelo", "Track velo #chi^2 probability", { 50, 0, 1 } } + , track_chi2ProbMatch{ owner, "trackChi2ProbMatch", "Track match #chi^2 probability", { 50, 0, 1 } } + , track_chi2ProbDownstreamVsPhi{ owner, + "trackChi2ProbDownstreamVsPhi", + "Track downstream #chi^2 probability vs phi", + { 50, -M_PI, M_PI } } + , track_ipX{ owner, "trackIPx", "Track IPx", { 100, -0.1, 0.1, "IPx [mm]" } } + , track_ipY{ owner, "trackIPy", "Track IPy", { 100, -0.1, 0.1, "IPy [mm]" } } , vertex_chi2{ owner, "chi2", "vertex chi2", { 100, 0, 5, "#chi^2" } } , vertex_chi2prob{ owner, "chi2prob", "vertex chi2prob", { 100, 0, 1, "#chi^2 probability" } } , vertex_x_pos{ owner, "vtxx", "x position of vertex", { 100, -x_range, x_range, "x [mm]" } } @@ -162,6 +186,8 @@ private: mutable Gaudi::Accumulators::Histogram<1> phi; mutable Gaudi::Accumulators::Histogram<1> matt; mutable Gaudi::Accumulators::Histogram<1> openingAngle; + mutable Gaudi::Accumulators::Histogram<1> ipchi2; + mutable Gaudi::Accumulators::Histogram<1> dls; BiasHistos( const TrackParticleMonitor* owner, float const& min_mass, float const& max_mass, float const& pdiff_range, unsigned int const& mass_bins, unsigned int const& pdiff_bins, @@ -201,7 +227,9 @@ private: "phimatt", "Angle between the decay plane and the y axis", { bins_angles * 5, 0, M_PI, "phiMatt [rad]" } } - , openingAngle{ owner, "openingangle", "Opening angle", { bins_angles * 5, 0, 0.3, "Opening angle [rad]" } } {} + , openingAngle{ owner, "openingangle", "Opening angle", { bins_angles * 5, 0, 0.3, "Opening angle [rad]" } } + , ipchi2{ owner, "ipchi2", "IP chi2", { 50, 0, 20 } } + , dls{ owner, "dls", "decay length significance", { 50, -5, 5 } } {} }; // Profile histograms of mass vs bias variables @@ -341,6 +369,20 @@ TrackParticleMonitor::TrackParticleMonitor( const std::string& name, ISvcLocator { KeyValue{ "InputLocation", "" }, KeyValue{ "StandardGeometryTop", LHCb::standard_geometry_top }, KeyValue{ "Magnet", LHCb::Det::Magnet::det_path } } } {} +namespace { + template + void fillFitResultHistograms( const LHCb::Track& track, const TrackFitResult& fitResult, TrackVertexHistos& histos ) { + if ( const auto tmp = fitResult.chi2(); tmp.nDoF() > 0 ) ++histos.track_chi2Prob[tmp.prob()]; + if ( const auto tmp = fitResult.chi2Velo(); tmp.nDoF() > 0 ) ++histos.track_chi2ProbVelo[tmp.prob()]; + if ( const auto tmp = fitResult.chi2Downstream(); tmp.nDoF() > 0 ) { + const auto pr = tmp.prob(); + ++histos.track_chi2ProbDownstream[pr]; + histos.track_chi2ProbDownstreamVsPhi[track.phi()] += pr; + } + if ( const auto tmp = fitResult.chi2Match(); tmp.nDoF() > 0 ) ++histos.track_chi2ProbMatch[tmp.prob()]; + } +} // namespace + // Algorithm void TrackParticleMonitor::operator()( LHCb::Particle::Range const& particles, DetectorElement const& lhcb, DeMagnet const& magnet ) const { @@ -360,11 +402,30 @@ void TrackParticleMonitor::operator()( LHCb::Particle::Range const& particles, D ++track_vertex_histos.track_pt[track->pt() / Gaudi::Units::GeV]; ++track_vertex_histos.track_p[track->p() / Gaudi::Units::GeV]; + // fill some histograms using the chi2 values from the kalmanfitresult + auto prFitResult = dynamic_cast( track->fitResult() ); + if ( prFitResult ) + fillFitResultHistograms( *track, *prFitResult, track_vertex_histos ); + else { + auto masterFitResult = dynamic_cast( track->fitResult() ); + if ( masterFitResult ) fillFitResultHistograms( *track, *masterFitResult, track_vertex_histos ); + } + LHCb::State* state = new LHCb::State(); m_stateprovider->stateFromTrajectory( *state, *track, z, *lhcb.geometry() ).ignore(); states.push_back( state ); } + // This is easier with the particles + for ( const auto& daughter : particle->daughters() ) + if ( daughter->proto() && daughter->proto()->track() && daughter->pv() ) { + const auto& state = daughter->proto()->track()->firstState(); + const auto& pvpos = daughter->pv()->position(); + const auto dz = pvpos.z() - state.z(); + ++track_vertex_histos.track_ipX[state.x() + dz * state.tx() - pvpos.x()]; + ++track_vertex_histos.track_ipY[state.y() + dz * state.ty() - pvpos.y()]; + } + // Obtain masses from PIDs double pdgmass( 0 ), pdgwidth( 0 ); const LHCb::ParticleProperty* prop = m_propertysvc->find( particle->particleID() ); @@ -463,7 +524,16 @@ void TrackParticleMonitor::operator()( LHCb::Particle::Range const& particles, D ++bias.phi[phiangle]; ++bias.matt[phimatt]; ++bias.openingAngle[openangle]; - + if ( particle->pv() ) { + // FIXME: what a mess this has become. I also don't manage to call it anymore with the output of TrackStateVertex + auto [chi2, decaylength, decaylength_err] = LHCb::StateVertexUtils::computeChiSquare( + referencePoint( *particle ), threeMomentum( *particle ), LHCb::Event::covMatrix( *particle ), + endVertexPos( *( particle->pv() ) ), posCovMatrix( *( particle->pv() ) ) ); + if ( decaylength_err > 0 ) { + ++bias.ipchi2[chi2.cast()]; + ++bias.dls[decaylength.cast() / decaylength_err.cast()]; + } + } // Profile histograms of mass vs bias variables auto& mass_vs_bias = m_histograms_mass_vs_bias.at( 0 ); mass_vs_bias.mass_vs_p[mom] += ( mass - pdgmass ) * Gaudi::Units::GeV; diff --git a/Tr/TrackMonitors/src/TrackVPOverlapMonitor.cpp b/Tr/TrackMonitors/src/TrackVPOverlapMonitor.cpp index c8d60cb4950f5a4ca870adf45b810d78eef92b4d..6e39e6792ffb8eacf89aadf3ed09b7dd0760f465 100644 --- a/Tr/TrackMonitors/src/TrackVPOverlapMonitor.cpp +++ b/Tr/TrackMonitors/src/TrackVPOverlapMonitor.cpp @@ -63,8 +63,16 @@ namespace LHCb::Tr::Monitor { LHCb::Detector::VPChannelID vpid; double residualX{ 0 }; double residualY{ 0 }; + const VPMicroCluster* cluster{ 0 }; }; + const VPMicroCluster* findcluster( const std::vector& clusters, + const LHCb::Detector::VPChannelID& vpid ) { + auto it = std::find_if( clusters.begin(), clusters.end(), + [vpid]( const auto& clus ) { return clus.channelID() == vpid; } ); + return it != clusters.end() ? &( *it ) : nullptr; + } + } // namespace struct TrackVPOverlapMonitor final @@ -225,8 +233,9 @@ namespace LHCb::Tr::Monitor { if ( node.isVP() ) { const auto vpid = id( node ).vpID(); if ( !currentnode.xnode || vpid != currentnode.vpid ) { - currentnode.vpid = vpid; - currentnode.xnode = &node; + currentnode.vpid = vpid; + currentnode.xnode = &node; + currentnode.cluster = findcluster( track.vpClusters(), vpid ); } else { currentnode.ynode = &node; xynodes.emplace_back( currentnode ); @@ -315,14 +324,32 @@ namespace LHCb::Tr::Monitor { ++x_y_tileoverlap_histo[{ globalpos.x(), globalpos.y() }]; // Obtain information about sensors and the position of the hits in global and local coordinates - const DeVPSensor& sens1 = det.sensor( node1.vpid ); - const DeVPSensor& sens2 = det.sensor( node2.vpid ); - const auto nodeGlobal1 = state( *( node1.xnode ) ).position(); - const auto nodeGlobal2 = state( *( node2.xnode ) ).position(); + const auto nodeGlobal1 = state( *( node1.xnode ) ).position(); + const auto nodeGlobal2 = state( *( node2.xnode ) ).position(); + // Compute the overlap residual in the local frame of the N-side sensor. Node 2 is the n-side node. + const DeVPSensor& sens1 = det.sensor( node1.vpid ); + const DeVPSensor& sens2 = det.sensor( node2.vpid ); const auto nodeLocal1 = sens1.globalToLocal( Gaudi::XYZPoint{ nodeGlobal1.x(), nodeGlobal1.y(), nodeGlobal1.z() } ); const auto nodeLocal2 = sens2.globalToLocal( Gaudi::XYZPoint{ nodeGlobal2.x(), nodeGlobal2.y(), nodeGlobal2.z() } ); + Gaudi::XYZVector localresidual; + // If we have the clusters, then derive the residual from the cluster positions and the local track slope. + if ( node1.cluster && node2.cluster ) { + // Translate both cluster positions and the track directions to the local frame of the n-side sensor + const auto localpos1 = sens2.globalToLocal( Gaudi::XYZPoint{ node1.cluster->globalPosition( det ) } ); + const auto localpos2 = sens2.globalToLocal( Gaudi::XYZPoint{ node2.cluster->globalPosition( det ) } ); + const auto localdir = sens2.globalToLocal( state( *node2.xnode ).slopes() ); + const auto dz = localpos1.z() - localpos2.z(); + const auto dx = localpos1.x() - ( localpos2.x() + dz * localdir.x() / localdir.z() ); + const auto dy = localpos1.y() - ( localpos2.y() + dz * localdir.y() / localdir.z() ); + localresidual = Gaudi::XYZVector{ dx, dy, 0 }; + } else { + // transform the residuals to a local frame (approximately) + Gaudi::XYZVector globalresidual{ node1.residualX - node2.residualX, node1.residualY - node2.residualY, + 0 }; + localresidual = sens2.globalToLocal( globalresidual ); + } // Implement cuts removing hits on the last row or column of pixels of the sensors and require 3 hits // before and after the hit of interest @@ -344,16 +371,16 @@ namespace LHCb::Tr::Monitor { } if ( sensor1 == 0 && sensor2 == 1 ) { // CLI-NLO - ++m_overlap_residual_x_CLI_NLO[{ module, node1.residualX - node2.residualX }]; - ++m_overlap_residual_y_CLI_NLO[{ module, node1.residualY - node2.residualY }]; + ++m_overlap_residual_x_CLI_NLO[{ module, localresidual.x() }]; + ++m_overlap_residual_y_CLI_NLO[{ module, localresidual.y() }]; ++module_CLI_NLO_histo[module]; } else if ( sensor1 == 0 && sensor2 == 2 ) { // CLI-NSI - ++m_overlap_residual_x_CLI_NSI[{ module, node1.residualX - node2.residualX }]; - ++m_overlap_residual_y_CLI_NSI[{ module, node1.residualY - node2.residualY }]; + ++m_overlap_residual_x_CLI_NSI[{ module, localresidual.x() }]; + ++m_overlap_residual_y_CLI_NSI[{ module, localresidual.y() }]; ++module_CLI_NSI_histo[module]; } else if ( sensor1 == 2 && sensor2 == 3 ) { // NSI-CSO - ++m_overlap_residual_x_CSO_NSI[{ module, node1.residualX - node2.residualX }]; - ++m_overlap_residual_y_CSO_NSI[{ module, node1.residualY - node2.residualY }]; + ++m_overlap_residual_x_CSO_NSI[{ module, localresidual.x() }]; + ++m_overlap_residual_y_CSO_NSI[{ module, localresidual.y() }]; ++module_CSO_NSI_histo[module]; } else { warning() << "How can these sensors overlap?: " << sensor1 << " " << sensor2 << endmsg; diff --git a/Tr/TrackMonitors/src/VPTrackMonitor.cpp b/Tr/TrackMonitors/src/VPTrackMonitor.cpp index c5ddc699bc39b0143f0d77da719c337b80d433e9..8c5ea845a99499dfa8127c05d13fa03d929cc190 100644 --- a/Tr/TrackMonitors/src/VPTrackMonitor.cpp +++ b/Tr/TrackMonitors/src/VPTrackMonitor.cpp @@ -423,10 +423,17 @@ namespace LHCb::Tr::Monitor { size_t tracknumber = 0; - auto get_cluster_ptr = [&clusters]( const auto& node ) { - auto i = std::find_if( begin( clusters ), end( clusters ), - [id = id( node ).vpID()]( const auto& c ) { return c.channelID() == id; } ); - return i != end( clusters ) ? &*i : nullptr; + auto get_cluster_position = [&clusters, &det]( const LHCb::Track& track, const Detector::VPChannelID vpid ) { + // use the micro-clusters if they are on the track. that saves a lot of time. + auto it = std::find_if( track.vpClusters().begin(), track.vpClusters().end(), + [vpid]( const auto& clus ) { return clus.channelID() == vpid; } ); + if ( it != track.vpClusters().end() ) { + return Gaudi::XYZPoint{ it->globalPosition( det ) }; + } else { + auto it = std::find_if( begin( clusters ), end( clusters ), + [vpid]( const auto& c ) { return c.channelID() == vpid; } ); + return it != clusters.end() ? Gaudi::XYZPoint{ it->x(), it->y(), it->z() } : Gaudi::XYZPoint{}; + } }; std::array, 8> N_exp{ { { 0 }, { 0 } } }; @@ -521,7 +528,6 @@ namespace LHCb::Tr::Monitor { // Skip non-VP measurements. if ( !node.isVP() ) continue; - auto iclus = get_cluster_ptr( node ); auto isXRes = hittypemap( node ) == VPX; auto isYRes = hittypemap( node ) == VPY; @@ -530,10 +536,9 @@ namespace LHCb::Tr::Monitor { // Get the sensor. const DeVPSensor& sens = det.sensor( chan ); - const auto corner = sens.localToGlobal( origin ); - const Gaudi::XYZPoint cluGlobal( iclus->x(), iclus->y(), iclus->z() ); - const Gaudi::XYZPoint nodeGlobal( state( node ).position().x(), state( node ).position().y(), - state( node ).position().z() ); + const auto corner = sens.localToGlobal( origin ); + const auto cluGlobal = get_cluster_position( *track, chan ); + const auto nodeGlobal = state( node ).position(); LHCb::State unbiasedState; if constexpr ( isTrackFitResult ) { diff --git a/Tr/TrackTools/src/TrackSelector.cpp b/Tr/TrackTools/src/TrackSelector.cpp index 1b123fb1dc2e50d11e8138a0c691a15da9ac425e..c2cb3ff414232eb2e8b19efca019a361ac31242d 100644 --- a/Tr/TrackTools/src/TrackSelector.cpp +++ b/Tr/TrackTools/src/TrackSelector.cpp @@ -33,6 +33,12 @@ using namespace LHCb; DECLARE_COMPONENT( TrackSelector ) //----------------------------------------------------------------------------- +namespace { + auto nOutliers( const PrKalmanFitResult& fr ) { + return std::count_if( fr.fitnodes.begin(), fr.fitnodes.end(), + [&]( const Pr::Tracks::Fit::Node& node ) { return node.isOutlier(); } ); + } +} // namespace bool TrackSelector::accept( const Track& aTrack ) const { @@ -63,7 +69,8 @@ bool TrackSelector::accept( const Track& aTrack ) const { const auto* fit = dynamic_cast( aTrack.fitResult() ); const auto* prfit = dynamic_cast( aTrack.fitResult() ); - if ( ( m_maxChi2Velo > 0 ) || ( m_maxChi2Upstream > 0 ) || ( m_maxChi2Downstream > 0 ) || ( m_maxChi2Match > 0 ) ) { + if ( ( m_maxChi2Velo > 0 ) || ( m_maxChi2Upstream > 0 ) || ( m_maxChi2Downstream > 0 ) || ( m_maxChi2Match > 0 ) || + m_maxNOutlier < 999 ) { if ( !fit && !prfit ) throw GaudiException( "TrackSelector", "Unknown or empty TrackFitResult", StatusCode::FAILURE ); } @@ -88,6 +95,13 @@ bool TrackSelector::accept( const Track& aTrack ) const { if ( msgLevel( MSG::VERBOSE ) ) verbose() << " -> Match Chi^2 " << chi2 << " failed cut" << endmsg; return false; } + if ( m_maxNOutlier < 999 ) { + const auto noutlier = fit ? fit->nOutliers() : nOutliers( *prfit ); + if ( noutlier > m_maxNOutlier ) { + if ( msgLevel( MSG::VERBOSE ) ) verbose() << " -> MaxNOutlier " << noutlier << " failed cut" << endmsg; + return false; + } + } // cut p const double p = aTrack.p(); diff --git a/Tr/TrackTools/src/TrackSelector.h b/Tr/TrackTools/src/TrackSelector.h index d5e09c86f99538df66f6f489398aaf56e7187407..05cbf6d556e9ec629fe52eebdec9d797271437af 100644 --- a/Tr/TrackTools/src/TrackSelector.h +++ b/Tr/TrackTools/src/TrackSelector.h @@ -112,6 +112,7 @@ private: Gaudi::Property m_minNUTLayers{ this, "MinNUTLayers", 0 }; Gaudi::Property m_maxNVeloHoles{ this, "MaxNVeloHoles", 999, "Maximum number of holes in Velo segment" }; Gaudi::Property m_maxNTHoles{ this, "MaxNTHoles", 999, "Maximum number of holes in T segment" }; + Gaudi::Property m_maxNOutlier{ this, "MaxNOutlier", 999, "Maximum number of outliers" }; }; #endif // TRACKTOOLS_TrackSelector_H