diff --git a/CMakeLists.txt b/CMakeLists.txt index 8cddf3ce165af7227ab0076db9155126f4eb4d8b..501ed418abd13941b00ef3c26cbca85cbf1e2c52 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -130,6 +130,7 @@ add_library(DetectorLib SHARED Detector/LHCb/src/DeLHCbHandles.cpp Detector/LHCb/src/InteractionRegion.cpp Detector/LHCb/src/VPGlobalCoordinateTransform.cpp + Detector/LHCb/src/OfflineMomentumScale.cpp Detector/LHCb/src/SMOGInfo.cpp Detector/LHCb/src/Tell40Links.cpp Detector/LHCb/src/LHCInfo.cpp diff --git a/Detector/LHCb/include/Detector/LHCb/DeLHCb.h b/Detector/LHCb/include/Detector/LHCb/DeLHCb.h index a4631024795411c525a50e1403d08127a9ade8e8..9fc0d58a589855e55d04496556b44dda182f413b 100644 --- a/Detector/LHCb/include/Detector/LHCb/DeLHCb.h +++ b/Detector/LHCb/include/Detector/LHCb/DeLHCb.h @@ -15,6 +15,7 @@ #include #include #include +#include #include #include #include @@ -37,6 +38,7 @@ namespace LHCb::Detector { std::optional m_lhcinfo; std::optional m_SMOG; std::optional m_vpCoordinateTransform; + std::optional m_momentumScale; }; // Utility method to lookup DeIOV object from the condition slice @@ -63,6 +65,9 @@ namespace LHCb::Detector { std::optional vpGlobalCoordinateTransformation() const { return this->access()->m_vpCoordinateTransform; } + + // Momentum scale, also only relevant offline + std::optional offlineMomentumScale() const { return this->access()->m_momentumScale; } }; // Utility method to setup DeLHCb diff --git a/Detector/LHCb/include/Detector/LHCb/OfflineMomentumScale.h b/Detector/LHCb/include/Detector/LHCb/OfflineMomentumScale.h new file mode 100644 index 0000000000000000000000000000000000000000..4fa23a2c7eab7685dc0320cc7804ce2f31375222 --- /dev/null +++ b/Detector/LHCb/include/Detector/LHCb/OfflineMomentumScale.h @@ -0,0 +1,131 @@ +/*****************************************************************************\ +* (c) Copyright 2025 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 + +#include +#include +#include + +#include +#include + +/** + * This condition is relevant for the offline processing of events, + * and is used to scale the values of the long-track states (in particular, + * the momentum) depending on their angles. + * + * The condition itself contains: + * - a boolean which indicates whether corrections are determined at all + * + * - a list of bins in tx, ty for the fine-grained correction factors + * - a per-run, or per-fill correction factor that is applied on top of these + * corrections in (tx,ty); this time a single number + * + * - a boolean which indicates that the corrections were determined + * after tracks have been refitted. This can be used to + * initate a warning in case the user applies these corrections on + * non-refitted data + * + * for 2 ranges of momenta, [0, [m_high_momentum_threshold] ], and anything above, + * the following information: + * - a map in terms of (tx, ty) of correction factors, in units of 1.e-4 + * which use the bins of (tx, ty) as shown above + * note that the (tx, ty) binning of high and low momentum tracks is the same + * + * The main spot in which this is used is in TrackScaleState, an algorithm + * in Rec. + */ +namespace LHCb::Detector { + struct OfflineMomentumScale { + + OfflineMomentumScale() = default; + OfflineMomentumScale( const nlohmann::json& obj ); + + /** When true, sensible values for the momentum-scale calibration are expected + to be present */ + bool enabled() const { return m_enabled; } + + /** When true, the corrections have been determined using tracks which have been + * refitted (see the RecoConf.track_refitting module). It is not adviced to use these + * to tracks which did not undergo the same treatment + */ + bool optimised_for_refit() const { return m_optimised_for_refit; } + + /** Retrieve the run-block-specific scaling factor for a given run number. + If the run is out of bounds for the known scaling factors, a nullopt is returned.*/ + std::optional getAlpha( unsigned int runNumber ) const; + + /** Returns the bin edges of the run numbers for which run-block-specific scaling factors + are applied */ + const std::set& binsRunNumber() const { return m_alpha_run_bins; } + + /** Momentum threshold at which the higher momentum correction map has been determined */ + float highMomentumThreshold() const { return m_high_momentum_threshold; } + + /** Returns the binning in tx used for the (tx, ty)-dependent correction map */ + const std::set& binsTx() const { return m_bins_tx; } + + /** Returns the binning in ty used for the (tx, ty)-dependent correction map */ + const std::set& binsTy() const { return m_bins_ty; } + + /** Returns the full correction map as a vector, of which the bin indices + * correspond to sequential (tx, ty) bins (see `getSequentialIndex` in this same + * class to have an accessor in terms of (tx, ty)). + * + * When highMomentum is true, the correction map is returned for tracks which have + * a momentum above the `highMomentumThreshold()`. + * Depending on q/p, positiveCharge should be set to true. + **/ + const std::vector& correctionMap( bool highMomentum, bool positiveCharge ) const; + + /** Retrieves the index of the correction map for a given tx, ty. + This can then be used in conjunction with the correctionMap() vector + to obtain the correction value. + + When the tx,ty vakues are out-of-bounds, a nullopt is returned. */ + std::optional getSequentialIndex( float tx, float ty ) const; + + private: + bool m_enabled = false; + bool m_optimised_for_refit = false; + + std::vector m_alpha = {}; + std::set m_alpha_run_bins = {}; + + float m_high_momentum_threshold = 20000.; + std::vector m_correction_map_plus_high_p = {}; + std::vector m_correction_map_minus_high_p = {}; + std::vector m_correction_map_plus_low_p = {}; + std::vector m_correction_map_minus_low_p = {}; + std::set m_bins_tx = {}; + std::set m_bins_ty = {}; + + /** + * Returns the index for a value within a range of bin edges. + * + * The lower bin edge is included: it acts as [begin, end] + * + * Note: bin id starts at 0. + * + * When it's out of bounds, the returned optional is not filled. + **/ + template + std::optional getBinId_1d( T val, std::set const& bins ) const { + if ( val < *bins.begin() || val > *bins.rbegin() ) return std::nullopt; + + auto lower = std::upper_bound( bins.begin(), bins.end(), val, std::less_equal{} ); + + return std::optional( std::distance( bins.begin(), lower ) - 1 ); + } + }; +} // namespace LHCb::Detector diff --git a/Detector/LHCb/src/DeLHCb.cpp b/Detector/LHCb/src/DeLHCb.cpp index 1f0c821857e10ece29231685ab1af295d35965b3..97f3f2a86278eb3bedfefee79e1fa0f81646bd60 100644 --- a/Detector/LHCb/src/DeLHCb.cpp +++ b/Detector/LHCb/src/DeLHCb.cpp @@ -49,6 +49,7 @@ LHCb::Detector::detail::DeLHCbObject::DeLHCbObject( const dd4hep::DetElement& if ( !SMOG.is_null() ) { m_SMOG = SMOG; } } } + { auto cond = ctxt.condition( hash_key( de, "OfflineVPGlobalCoordinateTransform" ), false ); if ( cond.isValid() ) { @@ -56,6 +57,15 @@ LHCb::Detector::detail::DeLHCbObject::DeLHCbObject( const dd4hep::DetElement& if ( !coordinateTransformData.is_null() ) { m_vpCoordinateTransform = coordinateTransformData; } } } + + { + auto cond = ctxt.condition( hash_key( de, "MomentumScale" ), false ); + if ( cond.isValid() ) { + auto momentumScaleData = cond.get(); + + if ( !momentumScaleData.is_null() ) { m_momentumScale = momentumScaleData; } + } + } } void LHCb::Detector::detail::DeLHCbObject::applyToAllChildren( @@ -154,6 +164,14 @@ void LHCb::Detector::setup_DeLHCb_callback( dd4hep::Detector& description ) { depbuilder.add( hash_key( de, "OfflineVPGlobalCoordinateTransform" ) ); } + // read possible offline momentum-scale calibration, if needed + if ( !schema || schema->has( "Conditions/LHCb/OfflineCalib/MomentumScale.yml", "MomentumScale" ) ) { + ( *requests ) + ->addLocation( de, LHCb::Detector::item_key( "MomentumScale" ), + "Conditions/LHCb/OfflineCalib/MomentumScale.yml", "MomentumScale" ); + depbuilder.add( hash_key( de, "MomentumScale" ) ); + } + if ( !schema || schema->has( "Conditions/LHCb/Online/LHC.yml", "LHC" ) ) { ( *requests )->addLocation( de, LHCb::Detector::item_key( "LHC" ), "Conditions/LHCb/Online/LHC.yml", "LHC" ); depbuilder.add( hash_key( de, "LHC" ) ); diff --git a/Detector/LHCb/src/OfflineMomentumScale.cpp b/Detector/LHCb/src/OfflineMomentumScale.cpp new file mode 100644 index 0000000000000000000000000000000000000000..92e27a7e541c129b668ba86615222029d5a5a0d7 --- /dev/null +++ b/Detector/LHCb/src/OfflineMomentumScale.cpp @@ -0,0 +1,82 @@ +/*****************************************************************************\ +* (c) Copyright 2025 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 + +#include + +#include +#include + +#include + +LHCb::Detector::OfflineMomentumScale::OfflineMomentumScale( const nlohmann::json& obj ) { + auto& enabledCondValue = obj.at( "enabled" ); + m_enabled = enabledCondValue.is_boolean() ? enabledCondValue.get() : bool( enabledCondValue.get() ); + + // if not enabled explicitly, don't even bother + // with setting the values and parsing the conditions database + if ( m_enabled ) { + auto& forRefitCondValue = obj.at( "optimised_for_refit" ); + m_optimised_for_refit = + forRefitCondValue.is_boolean() ? forRefitCondValue.get() : bool( forRefitCondValue.get() ); + + m_alpha = obj.at( "alpha" ).get>(); + m_alpha_run_bins = obj.at( "alpha_run_bins" ).get>(); + + m_high_momentum_threshold = obj.at( "high_momentum_threshold" ).get(); + + m_bins_tx = obj.at( "bins_tx" ).get>(); + m_bins_ty = obj.at( "bins_ty" ).get>(); + + if ( m_enabled && ( m_bins_tx.size() < 2 || m_bins_ty.size() < 2 ) ) { + throw std::runtime_error{ "The momentum-scale bin edges are incomplete, while the scaling is enabled." }; + } + + const auto expectedSizeOfCorrections = ( m_bins_tx.size() - 1 ) * ( m_bins_ty.size() - 1 ); + + m_correction_map_plus_high_p = obj.at( "correction_map_plus_high_p" ).get>(); + m_correction_map_minus_high_p = obj.at( "correction_map_minus_high_p" ).get>(); + m_correction_map_plus_low_p = obj.at( "correction_map_plus_low_p" ).get>(); + m_correction_map_minus_low_p = obj.at( "correction_map_minus_low_p" ).get>(); + + if ( expectedSizeOfCorrections != m_correction_map_plus_high_p.size() || + expectedSizeOfCorrections != m_correction_map_minus_high_p.size() || + expectedSizeOfCorrections != m_correction_map_plus_low_p.size() || + expectedSizeOfCorrections != m_correction_map_minus_low_p.size() ) { + throw std::runtime_error{ "The momentum-scale correction maps do not match the binning provided." }; + } + } +} + +const std::vector& LHCb::Detector::OfflineMomentumScale::correctionMap( bool highMomentum, + bool positiveCharge ) const { + if ( highMomentum ) return positiveCharge ? m_correction_map_plus_high_p : m_correction_map_minus_high_p; + + return positiveCharge ? m_correction_map_plus_low_p : m_correction_map_minus_low_p; +} + +std::optional LHCb::Detector::OfflineMomentumScale::getSequentialIndex( float tx, float ty ) const { + auto x_bin = getBinId_1d( tx, binsTx() ); + if ( !x_bin ) return std::nullopt; + + auto y_bin = getBinId_1d( ty, binsTy() ); + if ( !y_bin ) return std::nullopt; + + return *y_bin * ( m_bins_tx.size() - 1 ) + *x_bin; +} + +std::optional LHCb::Detector::OfflineMomentumScale::getAlpha( unsigned int runNumber ) const { + auto binId = getBinId_1d( runNumber, binsRunNumber() ); + + if ( !binId ) return std::nullopt; + + return m_alpha[*binId]; +} diff --git a/Detector/LHCb/src/VPGlobalCoordinateTransform.cpp b/Detector/LHCb/src/VPGlobalCoordinateTransform.cpp index 4631326a4c85f0fe44c6785b94cf7d150adc1e70..e39825cf0541fee9616df4b31332a9d409901e2c 100644 --- a/Detector/LHCb/src/VPGlobalCoordinateTransform.cpp +++ b/Detector/LHCb/src/VPGlobalCoordinateTransform.cpp @@ -42,10 +42,6 @@ LHCb::Detector::VPGlobalCoordinateTransform::VPGlobalCoordinateTransform( const "Could not find old position & rotation information - are you running the latest version of the software?" }; } - if ( !obj.contains( "newPosition" ) || !obj.contains( "newRotation" ) ) { - throw std::runtime_error{ "VPGlobalCoordinateTrransform condition is empty" }; - } - constexpr std::size_t posSize{ 3 }; constexpr std::size_t rotationSize{ 3 };