diff --git a/Tr/TrackMonitors/CMakeLists.txt b/Tr/TrackMonitors/CMakeLists.txt index 2fb7950cae1583770429380a2c4e51e045749b42..fd9bf5b5cdadda50d4750cb918afb55316caa988 100644 --- a/Tr/TrackMonitors/CMakeLists.txt +++ b/Tr/TrackMonitors/CMakeLists.txt @@ -37,6 +37,7 @@ gaudi_add_module(TrackMonitors src/UTTrackMonitor.cpp src/VPTrackMonitor.cpp src/VertexCompare.cpp + src/BeamSpotMonitor.cpp LINK AIDA::aida Gaudi::GaudiAlgLib diff --git a/Tr/TrackMonitors/src/BeamSpotMonitor.cpp b/Tr/TrackMonitors/src/BeamSpotMonitor.cpp new file mode 100644 index 0000000000000000000000000000000000000000..440b07bc4b328991ba9b9e0008e27b3183c52d01 --- /dev/null +++ b/Tr/TrackMonitors/src/BeamSpotMonitor.cpp @@ -0,0 +1,522 @@ +/*****************************************************************************\ +* (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. * +\*****************************************************************************/ +#include "Event/ODIN.h" +#include "Event/RecVertex.h" + +#include "LHCbAlgs/Consumer.h" +#include <Gaudi/Accumulators/Histogram.h> + +#include <filesystem> +#include <mutex> +#include <regex> +#include <sstream> +#include <string> +#include <yaml-cpp/yaml.h> + +namespace { + using PVView = LHCb::RecVertex::Range; + + /// Unbiased sample covariance calculator for off-diagonals + template <typename X, typename XX> + auto calculate_spread_offdiag( X const& x, X const& y, XX const& xy ) { + return ( xy.sum() - xy.nEntries() * x.mean() * y.mean() ) / ( xy.nEntries() - 1 ); + } + + /// Format position condition "<value>*mm" + template <typename T> + inline std::string format_position( T x ) { + return std::to_string( x / Gaudi::Units::mm ) + "*mm"; + } + + /// Format spread condition "<value>*mm2" + template <typename T> + inline std::string format_spread( T x ) { + return std::to_string( x / Gaudi::Units::mm2 ) + "*mm2"; + } + + /// Find largest version number in a directory + /// Adapted from AlignOnlineYMLCopier::AlignOnlineYMLCopier + template <typename T> + int find_max_version_number( T dir ) { + int max_version_nr = -1; + std::regex regex_version_file( "^v([0-9]+)$" ); + for ( const auto& entry : std::filesystem::directory_iterator( dir ) ) { + auto version_file_name = entry.path().stem().string(); + std::smatch matches; + if ( std::regex_search( version_file_name, matches, regex_version_file ) && matches.size() > 1 && + matches[1].matched ) { + int version_nr = std::stoi( matches[1] ); + if ( version_nr > max_version_nr ) max_version_nr = version_nr; + } + } + return max_version_nr; + } +} // namespace + +class BeamSpotMonitor : public LHCb::Algorithm::Consumer<void( LHCb::ODIN const&, PVView const& )> { + +public: + /// Standard constructor + BeamSpotMonitor( const std::string& name, ISvcLocator* pSvcLocator ); + + /// Algorithm execute + void operator()( LHCb::ODIN const&, PVView const& pvcontainer ) const override; + + /// Initialization + StatusCode initialize() override { + return Consumer::initialize().andThen( [&] { + // Conditions writing path safety imitating RICH RefIndexCalib + const auto curPath = std::filesystem::current_path(); + if ( m_writeDBFiles ) { + try { + if ( !std::filesystem::exists( m_conditionsDbPath.value() ) ) { + std::filesystem::create_directories( m_conditionsDbPath.value() ); + if ( msgLevel( MSG::DEBUG ) ) debug() << "Created " << m_conditionsDbPath.value() << endmsg; + } + } catch ( const std::filesystem::filesystem_error& expt ) { + warning() << "Cannot write to '" << m_conditionsDbPath.value() << "' -> Resetting Db path to '" + << curPath.string() << "/conditions'" << endmsg; + m_conditionsDbPath = curPath.string() + "/conditions"; + } + + const std::filesystem::path conditions_dir = + std::filesystem::path( m_conditionsDbPath.value() ) / std::filesystem::path( m_conditionsPathInDb.value() ); + try { + if ( !std::filesystem::exists( conditions_dir ) ) { + std::filesystem::create_directories( conditions_dir ); + if ( msgLevel( MSG::DEBUG ) ) debug() << "Created " << conditions_dir.string() << endmsg; + } + } catch ( const std::filesystem::filesystem_error& expt ) { + // TODO: consider whether this is a recoverable condition. + throw GaudiException( "Conditions directory does not exist and cannot be created", name(), + StatusCode::FAILURE ); + } + } + } ); + } + +private: + // Methods + //--------------------------------------------------------------------------- + /// Non-thread-safe part of the event loop --- resetting and publishing + void check_and_publish_reset( unsigned ) const; + + /// Conditions for and implementation of copying accumulating ctrs to cache + bool check_cache_counters() const; + void cache_counters() const; + + /// Conditions for and implementation of publication of cached conditions + bool check_publish() const; + bool ymlFormatter( std::ostringstream& ) const; + bool ymlFormatter( YAML::Emitter& ) const; + YAML::Node yamlNode() const; + bool ymlWriter() const; + + /// Conditions for value drifts + bool is_accum_delta_pos_over_thresh() const { + return ( std::abs( m_pvXPosCtr.mean() - m_c_pvXPosCtr.mean() ) > m_maxDeltaX ) || + ( std::abs( m_pvYPosCtr.mean() - m_c_pvYPosCtr.mean() ) > m_maxDeltaY ) || + ( std::abs( m_pvZPosCtr.mean() - m_c_pvZPosCtr.mean() ) > m_maxDeltaZ ); + } + bool is_accum_delta_spread_over_thresh() const; + + /// Reset the accumulators + bool check_reset_accumulators( const unsigned ) const; + void reset_accumulators( const unsigned ) const; + + // Properties + //--------------------------------------------------------------------------- + + // Output configuration properties + //........................................................................... + /// Report conditions to log (INFO) + Gaudi::Property<bool> m_condToLog{this, "LogConditions", false, + "Write conditions to logfile with level INFO when updating"}; + + /// Create direct run conditions + Gaudi::Property<bool> m_writeDBFiles{this, "MakeDBRunFile", false}; + + /// LHCB DB to write to path + Gaudi::Property<std::string> m_conditionsDbPath{this, "conditionsDbPath", "/group/online/hlt/conditions.run3"}; + Gaudi::Property<std::string> m_conditionsPathInDb{ + this, "conditionsPathInDb", "lhcb-conditions-database/Conditions/LHCb/Online/InteractionRegion.yml/.pool"}; + + // Accumulation configuration properties + //........................................................................... + Gaudi::Property<long unsigned> m_minPVsForCalib{this, "MinPVsForCalib", 100ul, + "Minumum number of accumulated PVs for a calibration"}; + + /// Thresholds triggering updated publication + // TODO: Fix the implementation---histograms must be defined after config + Gaudi::Property<double> m_maxDeltaX{this, "DeltaXMax", 2 * Gaudi::Units::mm, + "Maximum allowed absolute difference in mean X from cached value"}; + Gaudi::Property<double> m_maxDeltaY{this, "DeltaYMax", 2 * Gaudi::Units::mm, + "Maximum allowed absolute difference in mean Y from cached value"}; + Gaudi::Property<double> m_maxDeltaZ{this, "DeltaZMax", 20 * Gaudi::Units::mm, + "Maximum allowed absolute difference in mean Z from cached value"}; + + Gaudi::Property<double> m_maxDeltaXX{this, "DeltaXXMax", 2 * Gaudi::Units::mm2, + "Maximum allowed absolute difference in spread-matrix X^2 from cached value"}; + Gaudi::Property<double> m_maxDeltaXY{this, "DeltaXYMax", 2 * Gaudi::Units::mm2, + "Maximum allowed absolute difference in spread-matrix XY from cached value"}; + Gaudi::Property<double> m_maxDeltaYY{this, "DeltaYYMax", 2 * Gaudi::Units::mm2, + "Maximum allowed absolute difference in spread-matrix Y^2 from cached value"}; + Gaudi::Property<double> m_maxDeltaZX{this, "DeltaZXMax", 2 * Gaudi::Units::mm2, + "Maximum allowed absolute difference in spread-matrix ZX from cached value"}; + Gaudi::Property<double> m_maxDeltaYZ{this, "DeltaYZMax", 2 * Gaudi::Units::mm2, + "Maximum allowed absolute difference in spread-matrix YZ from cached value"}; + Gaudi::Property<double> m_maxDeltaZZ{this, "DeltaZZMax", 2 * Gaudi::Units::mm2, + "Maximum allowed absolute difference in spread-matrix Z^2 from cached value"}; + + /// Histogram limits + Gaudi::Property<double> m_histMaxXPV{this, "HistMaxXPV", 2 * Gaudi::Units::mm}; + Gaudi::Property<double> m_histMaxYPV{this, "HistMaxYPV", 2 * Gaudi::Units::mm}; + Gaudi::Property<double> m_histMinZPV{this, "HistMinZPV", -200 * Gaudi::Units::mm}; + Gaudi::Property<double> m_histMaxZPV{this, "HistMaxZPV", 200 * Gaudi::Units::mm}; + Gaudi::Property<double> m_histMinZPV_wide{this, "HistMinZPV_Wide", -1.5e3 * Gaudi::Units::mm, + "Wide z window for PV plot"}; + Gaudi::Property<double> m_histMaxZPV_wide{this, "HistMaxZPV_Wide", 1.5e3 * Gaudi::Units::mm, + "Wide z window for PV plot"}; + Gaudi::Property<double> m_histMaxXYPV{this, "HistMaxXYPV", 0.1 * Gaudi::Units::mm2}; + Gaudi::Property<double> m_histMaxYZPV{this, "HistMaxYZPV", 20 * Gaudi::Units::mm2}; + Gaudi::Property<double> m_histMaxZXPV{this, "HistMaxZXPV", 10.0 * Gaudi::Units::mm2}; + + // Accumulating statistics counters + //--------------------------------------------------------------------------- + mutable unsigned m_accRunNumber{0}; // Run number of last processed event + mutable Gaudi::Accumulators::SigmaCounter<> m_pvXPosCtr{this, "PV x position"}; + mutable Gaudi::Accumulators::SigmaCounter<> m_pvYPosCtr{this, "PV y position"}; + mutable Gaudi::Accumulators::SigmaCounter<> m_pvZPosCtr{this, "PV z position"}; + mutable Gaudi::Accumulators::AveragingCounter<> m_pvXYProdCtr{this, "PV x * y for covariance"}; + mutable Gaudi::Accumulators::AveragingCounter<> m_pvYZProdCtr{this, "PV y * z for covariance"}; + mutable Gaudi::Accumulators::AveragingCounter<> m_pvZXProdCtr{this, "PV z * x for covariance"}; + + // Cached statistics counters and variables + //--------------------------------------------------------------------------- + mutable unsigned m_c_accRunNumber{0}; // Run number of last cached event + mutable Gaudi::Accumulators::SigmaCounter<> m_c_pvXPosCtr{this, "Cached PV x position"}; + mutable Gaudi::Accumulators::SigmaCounter<> m_c_pvYPosCtr{this, "Cached PV y position"}; + mutable Gaudi::Accumulators::SigmaCounter<> m_c_pvZPosCtr{this, "Cached PV z position"}; + mutable Gaudi::Accumulators::AveragingCounter<> m_c_pvXYProdCtr{this, "Cached PV x * y for covariance"}; + mutable Gaudi::Accumulators::AveragingCounter<> m_c_pvYZProdCtr{this, "Cached PV y * z for covariance"}; + mutable Gaudi::Accumulators::AveragingCounter<> m_c_pvZXProdCtr{this, "Cached PV z * x for covariance"}; + + // Histograms + //--------------------------------------------------------------------------- + mutable Gaudi::Accumulators::Histogram<1> m_numPrimaryVertices{ + this, "NumPrimaryVertices", "Primary vertices per event;Number of PVs;Entries", {11, -0.5, 10.5}}; + mutable Gaudi::Accumulators::Histogram<1> m_pvNTracks{ + this, "NumTracks", "Number of tracks per primary vertex;Number of tracks;Entries", {101, -0.5, 100.5}}; + mutable Gaudi::Accumulators::Histogram<1> m_pvXPosition{ + this, "PVPosX", "Primary vertex X position;PV X [mm];Entries", {200, -m_histMaxXPV, m_histMaxXPV}}; + mutable Gaudi::Accumulators::Histogram<1> m_pvYPosition{ + this, "PVPosY", "Primary vertex Y position;PV Y [mm];Entries", {200, -m_histMaxYPV, m_histMaxYPV}}; + mutable Gaudi::Accumulators::Histogram<1> m_pvZPosition{ + this, "PVPosZ", "Primary vertex Z position;PV Z [mm];Entries", {200, m_histMinZPV, m_histMaxZPV}}; + mutable Gaudi::Accumulators::Histogram<1> m_pvZPositionWide{this, + "PVPosZWide", + "Primary vertex Z position (wide);PV Z [mm];Entries", + {200, m_histMinZPV_wide, m_histMaxZPV_wide}}; + mutable Gaudi::Accumulators::Histogram<1> m_pvXYPosCov{ + this, + "PVCovXY", + "Interaction region covariance X-Y;PV (X - X_mean)*(Y - Y_mean) [mm2];Entries", + {200, -m_histMaxXYPV, m_histMaxXYPV}}; + mutable Gaudi::Accumulators::Histogram<1> m_pvYZPosCov{ + this, + "PVCovYZ", + "Interaction region covariance Y-Z;PV (Y - Y_mean)*(Z - Z_mean) [mm2];Entries", + {200, -m_histMaxYZPV, m_histMaxYZPV}}; + mutable Gaudi::Accumulators::Histogram<1> m_pvZXPosCov{ + this, + "PVCovZX", + "Interaction region covariance Z-X;PV (Z - Z_mean)*(X - X_mean) [mm2];Entries", + {200, -m_histMaxZXPV, m_histMaxZXPV}}; + + // Warning counters + //--------------------------------------------------------------------------- + mutable Gaudi::Accumulators::MsgCounter<MSG::ERROR> m_badPVCount{this, "Too few PVs to produce IR conditions"}; + + // Other data members + //--------------------------------------------------------------------------- + /// Lock for main event processing + mutable std::mutex m_mutex; +}; + +// Declaration of the Algorithm Factory +DECLARE_COMPONENT( BeamSpotMonitor ) + +//============================================================================= +// Standard constructor, initializes variables +//============================================================================= +BeamSpotMonitor::BeamSpotMonitor( const std::string& name, ISvcLocator* pSvcLocator ) + : Consumer{name, + pSvcLocator, + {KeyValue{"ODINLocation", LHCb::ODINLocation::Default}, + KeyValue{"PVContainer", LHCb::RecVertexLocation::Primary}}} {} + +//============================================================================= +// Algorithm execution +//============================================================================= +void BeamSpotMonitor::operator()( LHCb::ODIN const& odin, PVView const& pvcontainer ) const { + const auto curRunNumber = odin.runNumber(); + check_and_publish_reset( curRunNumber ); + + // Skip this event if it is from a run prior to the one we are accumulating. + if ( curRunNumber < m_accRunNumber ) return; + + // number of primary vertices + ++m_numPrimaryVertices[pvcontainer.size()]; + + for ( const auto& pv : pvcontainer ) { + m_pvXPosCtr += pv->position().x(); + m_pvYPosCtr += pv->position().y(); + m_pvZPosCtr += pv->position().z(); + m_pvXYProdCtr += pv->position().x() * pv->position().y(); + m_pvYZProdCtr += pv->position().y() * pv->position().z(); + m_pvZXProdCtr += pv->position().z() * pv->position().x(); + + ++m_pvNTracks[pv->tracks().size()]; + ++m_pvXPosition[pv->position().x()]; + ++m_pvYPosition[pv->position().y()]; + ++m_pvZPosition[pv->position().z()]; + ++m_pvZPositionWide[pv->position().z()]; + // Must be incremented after accumulators to avoid excursion on first PV. + ++m_pvXYPosCov[( pv->position().x() - m_pvXPosCtr.mean() ) * ( pv->position().y() - m_pvYPosCtr.mean() )]; + ++m_pvYZPosCov[( pv->position().y() - m_pvYPosCtr.mean() ) * ( pv->position().z() - m_pvZPosCtr.mean() )]; + ++m_pvZXPosCov[( pv->position().z() - m_pvZPosCtr.mean() ) * ( pv->position().x() - m_pvXPosCtr.mean() )]; + } +} + +//============================================================================= +// Non-thread-safe resetting and publishing +//============================================================================= +void BeamSpotMonitor::check_and_publish_reset( unsigned curRunNumber ) const { + std::scoped_lock lock{m_mutex}; + + // Handle cacheing and publication before processing the event. + if ( check_reset_accumulators( curRunNumber ) ) { + if ( check_cache_counters() ) { + // Copy values of accumulating counters to cache + cache_counters(); + + if ( check_publish() ) { ymlWriter(); } + } + // Reset accumulators + reset_accumulators( curRunNumber ); + } +} + +//============================================================================= +// Check for cacheing the counters +//============================================================================= +bool BeamSpotMonitor::check_cache_counters() const { + return ( m_pvXPosCtr.nEntries() >= m_minPVsForCalib ) && + ( m_c_pvXPosCtr.nEntries() == 0 || is_accum_delta_pos_over_thresh() || is_accum_delta_spread_over_thresh() ); +} + +//============================================================================= +// Check whether to publish the cached counters +//============================================================================= +bool BeamSpotMonitor::check_publish() const { return m_c_pvXPosCtr.nEntries() >= m_minPVsForCalib; } + +//============================================================================= +// Check whether to reset the accumulators +// Reset when +// - Accumulated at least the target number of PVs +// - Processing an event from a new, subsequent run. +//============================================================================= +bool BeamSpotMonitor::check_reset_accumulators( const unsigned curRunNumber ) const { + return ( m_pvXPosCtr.nEntries() >= m_minPVsForCalib ) || ( m_accRunNumber < curRunNumber ); +} + +//============================================================================= +// Cache the counters +//============================================================================= +void BeamSpotMonitor::cache_counters() const { + m_c_pvXPosCtr.reset(); + m_c_pvYPosCtr.reset(); + m_c_pvZPosCtr.reset(); + m_c_pvXYProdCtr.reset(); + m_c_pvYZProdCtr.reset(); + m_c_pvZXProdCtr.reset(); + + m_c_pvXPosCtr.mergeAndReset( m_pvXPosCtr ); + m_c_pvYPosCtr.mergeAndReset( m_pvYPosCtr ); + m_c_pvZPosCtr.mergeAndReset( m_pvZPosCtr ); + m_c_pvXYProdCtr.mergeAndReset( m_pvXYProdCtr ); + m_c_pvYZProdCtr.mergeAndReset( m_pvYZProdCtr ); + m_c_pvZXProdCtr.mergeAndReset( m_pvZXProdCtr ); + m_c_accRunNumber = m_accRunNumber; +} + +//============================================================================= +// Check whether deltas of the spread matrix have exceeded thresholds. +//============================================================================= +bool BeamSpotMonitor::is_accum_delta_spread_over_thresh() const { + return ( std::abs( m_c_pvXPosCtr.unbiased_sample_variance() - m_pvXPosCtr.unbiased_sample_variance() ) > + m_maxDeltaXX ) || + ( std::abs( m_c_pvYPosCtr.unbiased_sample_variance() - m_pvYPosCtr.unbiased_sample_variance() ) > + m_maxDeltaYY ) || + ( std::abs( m_c_pvZPosCtr.unbiased_sample_variance() - m_pvZPosCtr.unbiased_sample_variance() ) > + m_maxDeltaZZ ) || + ( std::abs( calculate_spread_offdiag( m_c_pvXPosCtr, m_c_pvYPosCtr, m_c_pvXYProdCtr ) - + calculate_spread_offdiag( m_pvXPosCtr, m_pvYPosCtr, m_pvXYProdCtr ) ) > m_maxDeltaXY ) || + ( std::abs( calculate_spread_offdiag( m_c_pvZPosCtr, m_c_pvXPosCtr, m_c_pvZXProdCtr ) - + calculate_spread_offdiag( m_pvZPosCtr, m_pvXPosCtr, m_pvZXProdCtr ) ) > m_maxDeltaZX ) || + ( std::abs( calculate_spread_offdiag( m_c_pvYPosCtr, m_c_pvZPosCtr, m_c_pvYZProdCtr ) - + calculate_spread_offdiag( m_pvYPosCtr, m_pvZPosCtr, m_pvYZProdCtr ) ) > m_maxDeltaYZ ); +} + +//============================================================================= +// Reset the accumulators +//============================================================================= +void BeamSpotMonitor::reset_accumulators( const unsigned curRunNumber ) const { + m_pvXPosCtr.reset(); + m_pvYPosCtr.reset(); + m_pvZPosCtr.reset(); + m_pvXYProdCtr.reset(); + m_pvYZProdCtr.reset(); + m_pvZXProdCtr.reset(); + m_accRunNumber = curRunNumber; +} + +//============================================================================= +// YML formatter +// +// Format of example condition +// lhcb-conditions-database/Conditions/LHCb/Online/InteractionRegion.yml/0/200 +// +// From its comments, the ordering of elements should be +// InteractionRegion: +// position: <avg_position_x_y_z: [float x3]> +// spread: <spread_matrix_xx_xy_yy_xz_yz_zz: [float x6]> +// +// The example conditions are +// InteractionRegion: +// position: [ "0.0*mm", "0.0*mm", "0.0*mm" ] +// spread: [ "0.0064*mm2", "0.0*mm2", "0.0064*mm2", "0.0*mm2", "0.0*mm2", "2809.0*mm2" ] +// +// Note that the numeric values and their units are enclosed in double quotes. +//============================================================================= +bool BeamSpotMonitor::ymlFormatter( std::ostringstream& ymlbuff ) const { + // + if ( m_c_pvXPosCtr.nEntries() < m_minPVsForCalib ) { + ++m_badPVCount; + return false; + } + + ymlbuff << "InteractionRegion:\n"; + ymlbuff << " position: [ \"" << format_position( m_c_pvXPosCtr.mean() ) << "\", \"" + << format_position( m_c_pvYPosCtr.mean() ) << "\", \"" << format_position( m_c_pvZPosCtr.mean() ) << "\" ]\n"; + ymlbuff << " spread: [ \"" << format_spread( m_c_pvXPosCtr.unbiased_sample_variance() ) << "\", \"" + << format_spread( calculate_spread_offdiag( m_c_pvXPosCtr, m_c_pvYPosCtr, m_c_pvXYProdCtr ) ) << "\", \"" + << format_spread( m_c_pvYPosCtr.unbiased_sample_variance() ) << "\", \"" + << format_spread( calculate_spread_offdiag( m_c_pvZPosCtr, m_c_pvXPosCtr, m_c_pvZXProdCtr ) ) << "\", \"" + << format_spread( calculate_spread_offdiag( m_c_pvYPosCtr, m_c_pvZPosCtr, m_c_pvYZProdCtr ) ) << "\", \"" + << format_spread( m_c_pvZPosCtr.unbiased_sample_variance() ) << "\" ]\n"; + + return true; +} + +bool BeamSpotMonitor::ymlFormatter( YAML::Emitter& yout ) const { + // + if ( m_c_pvXPosCtr.nEntries() < m_minPVsForCalib ) { + ++m_badPVCount; + return false; + } + + yout.SetSeqFormat( YAML::Flow ); + + yout << YAML::BeginMap; + yout << YAML::Key << "InteractionRegion"; + yout << YAML::Value << YAML::BeginMap; + yout << YAML::Key << "position"; + yout << YAML::Value << YAML::BeginSeq; + yout << YAML::DoubleQuoted << format_position( m_c_pvXPosCtr.mean() ) << YAML::DoubleQuoted + << format_position( m_c_pvYPosCtr.mean() ) << YAML::DoubleQuoted << format_position( m_c_pvZPosCtr.mean() ) + << YAML::EndSeq; + yout << YAML::Key << "spread"; + yout << YAML::Value << YAML::BeginSeq; + yout << YAML::DoubleQuoted << format_spread( m_c_pvXPosCtr.unbiased_sample_variance() ) << YAML::DoubleQuoted + << format_spread( calculate_spread_offdiag( m_c_pvXPosCtr, m_c_pvYPosCtr, m_c_pvXYProdCtr ) ) + << YAML::DoubleQuoted << format_spread( m_c_pvYPosCtr.unbiased_sample_variance() ) << YAML::DoubleQuoted + << format_spread( calculate_spread_offdiag( m_c_pvZPosCtr, m_c_pvXPosCtr, m_c_pvZXProdCtr ) ) + << YAML::DoubleQuoted + << format_spread( calculate_spread_offdiag( m_c_pvYPosCtr, m_c_pvZPosCtr, m_c_pvYZProdCtr ) ) + << YAML::DoubleQuoted << format_spread( m_c_pvZPosCtr.unbiased_sample_variance() ) << YAML::EndSeq; + yout << YAML::EndMap << YAML::EndMap; + + return true; +} + +//============================================================================= +// Experimental construction of YAML::Node +// Ultimately, the format of InteractionRegion.yml, with double quoted +// string values but unquoted string keys, requires a custom emitter method +// using the current yaml-cpp interface. +// +// This can be a more robuse way of construction the YAML if the formatting +// complications are resolved or determined to be irrelevant. +//============================================================================= +YAML::Node BeamSpotMonitor::yamlNode() const { + const std::string irkey = "InteractionRegion"; + const std::string pkey = "position"; + const std::string skey = "spread"; + + YAML::Node ir; + ir[irkey][pkey].push_back( format_position( m_c_pvXPosCtr.mean() ) ); + ir[irkey][pkey].push_back( format_position( m_c_pvYPosCtr.mean() ) ); + ir[irkey][pkey].push_back( format_position( m_c_pvZPosCtr.mean() ) ); + + ir[irkey][skey].push_back( format_spread( m_c_pvXPosCtr.unbiased_sample_variance() ) ); + ir[irkey][skey].push_back( + format_spread( calculate_spread_offdiag( m_c_pvXPosCtr, m_c_pvYPosCtr, m_c_pvXYProdCtr ) ) ); + ir[irkey][skey].push_back( format_spread( m_c_pvYPosCtr.unbiased_sample_variance() ) ); + ir[irkey][skey].push_back( + format_spread( calculate_spread_offdiag( m_c_pvZPosCtr, m_c_pvXPosCtr, m_c_pvZXProdCtr ) ) ); + ir[irkey][skey].push_back( + format_spread( calculate_spread_offdiag( m_c_pvYPosCtr, m_c_pvZPosCtr, m_c_pvYZProdCtr ) ) ); + ir[irkey][skey].push_back( format_spread( m_c_pvZPosCtr.unbiased_sample_variance() ) ); + + ir[irkey][pkey].SetStyle( YAML::EmitterStyle::Flow ); + ir[irkey][skey].SetStyle( YAML::EmitterStyle::Flow ); + + return ir; +} + +bool BeamSpotMonitor::ymlWriter() const { + std::ostringstream ymlBuff; + if ( !ymlFormatter( ymlBuff ) ) { + throw GaudiException( "YAML formatting of conditions failed - unable to write new conditions", name(), + StatusCode::FAILURE ); + } + + if ( m_condToLog ) info() << ymlBuff.str() << endmsg; + + if ( m_writeDBFiles ) { + // Adapting file handling from RICH RefIndexCalib + const std::filesystem::path conditions_dir = + std::filesystem::path( m_conditionsDbPath.value() ) / std::filesystem::path( m_conditionsPathInDb.value() ); + + const auto max_version_nr = find_max_version_number( conditions_dir ); + + const std::filesystem::path conditions_file( "v" + std::to_string( max_version_nr + 1 ) ); + const std::filesystem::path full_conditions_path = conditions_dir / conditions_file; + + // Outfile. + std::ofstream logging( full_conditions_path.string().c_str() ); + logging << ymlBuff.str(); + logging.close(); + } + + return true; +}