From c628db3d00442877681732dcd5c7634cfaadcbb2 Mon Sep 17 00:00:00 2001
From: Bogdan Kutsenko <bokutsen@n4050101.lbdaq.cern.ch>
Date: Mon, 22 Apr 2024 22:41:55 +0200
Subject: [PATCH 01/21] Modifications for vertex compare to work with RecoMon
---
Tr/TrackMonitors/src/VertexCompare.cpp | 53 ++++++++++++++++++++------
1 file changed, 42 insertions(+), 11 deletions(-)
diff --git a/Tr/TrackMonitors/src/VertexCompare.cpp b/Tr/TrackMonitors/src/VertexCompare.cpp
index 26f052b7cff..a9d2c209720 100644
--- a/Tr/TrackMonitors/src/VertexCompare.cpp
+++ b/Tr/TrackMonitors/src/VertexCompare.cpp
@@ -21,6 +21,7 @@
#include "Event/Track_v2.h"
#include "GaudiAlg/GaudiAlgorithm.h"
#include "GaudiAlg/GaudiTupleAlg.h"
+#include <Gaudi/Accumulators/Histogram.h>
#include "GaudiAlg/Tuples.h"
#include "GaudiKernel/PhysicalConstants.h"
#include "GaudiKernel/SystemOfUnits.h"
@@ -191,8 +192,16 @@ public:
private:
bool debugLevel() const { return msgLevel( MSG::DEBUG ) || msgLevel( MSG::VERBOSE ); }
- Gaudi::Property<bool> m_produceHistogram{this, "produceHistogram", true};
- Gaudi::Property<bool> m_produceNtuple{this, "produceNtuple", true};
+
+ mutable std::optional<Gaudi::Accumulators::Histogram<1>> m_histo_dx;
+ mutable std::optional<Gaudi::Accumulators::Histogram<1>> m_histo_dy;
+ mutable std::optional<Gaudi::Accumulators::Histogram<1>> m_histo_dz;
+ mutable std::optional<Gaudi::Accumulators::Histogram<1>> m_histo_pullx;
+ mutable std::optional<Gaudi::Accumulators::Histogram<1>> m_histo_pully;
+ mutable std::optional<Gaudi::Accumulators::Histogram<1>> m_histo_pullz;
+
+ Gaudi::Property<bool> m_produceHistogram{this, "produceHistogram", false};
+ Gaudi::Property<bool> m_produceNtuple{this, "produceNtuple", false};
mutable Gaudi::Accumulators::Counter<> m_nVtx{this, "Number of pairs of vertices in processed events"};
mutable Gaudi::Accumulators::Counter<> m_nRec1{this, "Number of vertices in input1"};
mutable Gaudi::Accumulators::Counter<> m_nRec2{this, "Number of vertices in input2"};
@@ -217,7 +226,17 @@ DECLARE_COMPONENT_WITH_ID( VertexCompare, "VertexCompare" )
// Initialization
//=============================================================================
StatusCode VertexCompare::initialize() {
+
+
return Consumer::initialize().andThen( [&]() {
+ using axis1D = Gaudi::Accumulators::Axis<decltype( m_histo_dx )::value_type::AxisArithmeticType>;
+ m_histo_dx.emplace( this, "dx", "dx, mm", axis1D{50, -0.25, 0.25} );
+ m_histo_dy.emplace( this, "dy", "dx, mm", axis1D{50, -0.25, 0.25} );
+ m_histo_dz.emplace( this, "dz", "dx, mm", axis1D{50, -1.5, 1.5} );
+ m_histo_pullx.emplace( this, "pullx", "pull x", axis1D{20, -2, 2} );
+ m_histo_pully.emplace( this, "pully", "pull y", axis1D{20, -2, 2} );
+ m_histo_pullz.emplace( this, "pullz", "pull z", axis1D{20, -2, 2} );
+
m_nVtx.reset();
m_nEvt.reset();
LHCb::Conditions::InteractionRegion::addConditionDerivation( this,
@@ -370,7 +389,16 @@ void VertexCompare::operator()( Vertices const& recoVtx1, Vertices const& recoVt
multiplicity = link.at( oIt ) + 1;
}
+
+ ++m_histo_dx.value()[dx];
+ ++m_histo_dy.value()[dy];
+ ++m_histo_dz.value()[dz];
+ ++m_histo_pullx.value()[dx / errx];
+ ++m_histo_pully.value()[dy / erry];
+ ++m_histo_pullz.value()[dz / errz];
+
if ( m_produceHistogram.value() ) {
+
plot( dx, 1021, "dx", -0.25, 0.25, 50 );
plot( dy, 1022, "dy", -0.25, 0.25, 50 );
plot( dz, 1023, "dz", -1.5, 1.5, 60 );
@@ -615,13 +643,16 @@ StatusCode VertexCompare::finalize() {
<< ( (double)m_nVtx.nEntries() - (double)m_nVtx_geq10tr_matched2.nEntries() ) /
( (double)m_nRec2.nEntries() - (double)m_nRec_geq10tr_2.nEntries() ) * 100
<< " % " << endmsg;
-
- const AIDA::IHistogram1D* dx = histo( HistoID( 1021 ) );
- const AIDA::IHistogram1D* pullx = histo( HistoID( 1031 ) );
- const AIDA::IHistogram1D* dy = histo( HistoID( 1022 ) );
- const AIDA::IHistogram1D* pully = histo( HistoID( 1032 ) );
- const AIDA::IHistogram1D* dz = histo( HistoID( 1023 ) );
- const AIDA::IHistogram1D* pullz = histo( HistoID( 1033 ) );
+
+ //const AIDA::IHistogram1D* dx = histo( HistoID( 1021 ) );
+ //const AIDA::IHistogram1D* pullx = histo( HistoID( 1031 ) );
+ //const AIDA::IHistogram1D* dy = histo( HistoID( 1022 ) );
+ //const AIDA::IHistogram1D* pully = histo( HistoID( 1032 ) );
+ //const AIDA::IHistogram1D* dz = histo( HistoID( 1023 ) );
+ //const AIDA::IHistogram1D* pullz = histo( HistoID( 1033 ) );
+
+
+ /*
if ( dx ) {
info() << " ---------------------------------------" << endmsg;
info() << "dx: "
@@ -683,7 +714,7 @@ StatusCode VertexCompare::finalize() {
<< endmsg;
}
info() << " ============================================" << endmsg;
-
+ */
return Consumer::finalize(); // Must be called after all other actions
}
@@ -742,4 +773,4 @@ void VertexCompare::printRat( std::string mes, int a, int b ) {
pmes += " : ";
info() << pmes << format( " %6.3f ( %7d / %8d )", rat, a, b ) << endmsg;
-}
\ No newline at end of file
+}
--
GitLab
From d979458039f2004deca3090afbdb9b28eb55c4d2 Mon Sep 17 00:00:00 2001
From: Bogdan Kutsenko <bokutsen@n4050101.lbdaq.cern.ch>
Date: Thu, 25 Apr 2024 15:58:09 +0200
Subject: [PATCH 02/21] Selections for monitoring added
---
Tr/TrackMonitors/src/VertexCompare.cpp | 62 +++++++++++++-------------
1 file changed, 30 insertions(+), 32 deletions(-)
diff --git a/Tr/TrackMonitors/src/VertexCompare.cpp b/Tr/TrackMonitors/src/VertexCompare.cpp
index a9d2c209720..149d1635a29 100644
--- a/Tr/TrackMonitors/src/VertexCompare.cpp
+++ b/Tr/TrackMonitors/src/VertexCompare.cpp
@@ -192,16 +192,19 @@ public:
private:
bool debugLevel() const { return msgLevel( MSG::DEBUG ) || msgLevel( MSG::VERBOSE ); }
-
- mutable std::optional<Gaudi::Accumulators::Histogram<1>> m_histo_dx;
- mutable std::optional<Gaudi::Accumulators::Histogram<1>> m_histo_dy;
- mutable std::optional<Gaudi::Accumulators::Histogram<1>> m_histo_dz;
- mutable std::optional<Gaudi::Accumulators::Histogram<1>> m_histo_pullx;
- mutable std::optional<Gaudi::Accumulators::Histogram<1>> m_histo_pully;
- mutable std::optional<Gaudi::Accumulators::Histogram<1>> m_histo_pullz;
+
+ using axis1D = Gaudi::Accumulators::Axis<Gaudi::Accumulators::Histogram<1>::AxisArithmeticType>;
+ mutable Gaudi::Accumulators::Histogram<1> m_histo_dx{ this, "dx", "dx, mm", axis1D{50, -0.25, 0.25} };
+ mutable Gaudi::Accumulators::Histogram<1> m_histo_dy{ this, "dy", "dy, mm", axis1D{50, -0.25, 0.25} };
+ mutable Gaudi::Accumulators::Histogram<1> m_histo_dz{ this, "dz", "dz, mm", axis1D{50, -1.5, 1.5} };
+ mutable Gaudi::Accumulators::Histogram<1> m_histo_pullx{ this, "pullx", "pull x", axis1D{20, -2, 2} };
+ mutable Gaudi::Accumulators::Histogram<1> m_histo_pully{ this, "pully", "pull y", axis1D{20, -2, 2} };
+ mutable Gaudi::Accumulators::Histogram<1> m_histo_pullz{ this, "pullz", "pull z", axis1D{20, -2, 2} };
Gaudi::Property<bool> m_produceHistogram{this, "produceHistogram", false};
Gaudi::Property<bool> m_produceNtuple{this, "produceNtuple", false};
+ Gaudi::Property<bool> m_monitoring{this, "monitoring", true};
+
mutable Gaudi::Accumulators::Counter<> m_nVtx{this, "Number of pairs of vertices in processed events"};
mutable Gaudi::Accumulators::Counter<> m_nRec1{this, "Number of vertices in input1"};
mutable Gaudi::Accumulators::Counter<> m_nRec2{this, "Number of vertices in input2"};
@@ -229,14 +232,6 @@ StatusCode VertexCompare::initialize() {
return Consumer::initialize().andThen( [&]() {
- using axis1D = Gaudi::Accumulators::Axis<decltype( m_histo_dx )::value_type::AxisArithmeticType>;
- m_histo_dx.emplace( this, "dx", "dx, mm", axis1D{50, -0.25, 0.25} );
- m_histo_dy.emplace( this, "dy", "dx, mm", axis1D{50, -0.25, 0.25} );
- m_histo_dz.emplace( this, "dz", "dx, mm", axis1D{50, -1.5, 1.5} );
- m_histo_pullx.emplace( this, "pullx", "pull x", axis1D{20, -2, 2} );
- m_histo_pully.emplace( this, "pully", "pull y", axis1D{20, -2, 2} );
- m_histo_pullz.emplace( this, "pullz", "pull z", axis1D{20, -2, 2} );
-
m_nVtx.reset();
m_nEvt.reset();
LHCb::Conditions::InteractionRegion::addConditionDerivation( this,
@@ -389,16 +384,18 @@ void VertexCompare::operator()( Vertices const& recoVtx1, Vertices const& recoVt
multiplicity = link.at( oIt ) + 1;
}
+ if ( m_monitoring.value() ) {
+ if(dz<2 && size_diff < 3 && multiplicity < 2){
+ ++m_histo_dx[dx];
+ ++m_histo_dy[dy];
+ ++m_histo_dz[dz];
+ ++m_histo_pullx[dx / errx];
+ ++m_histo_pully[dy / erry];
+ ++m_histo_pullz[dz / errz];
+ }
+ }
- ++m_histo_dx.value()[dx];
- ++m_histo_dy.value()[dy];
- ++m_histo_dz.value()[dz];
- ++m_histo_pullx.value()[dx / errx];
- ++m_histo_pully.value()[dy / erry];
- ++m_histo_pullz.value()[dz / errz];
-
if ( m_produceHistogram.value() ) {
-
plot( dx, 1021, "dx", -0.25, 0.25, 50 );
plot( dy, 1022, "dy", -0.25, 0.25, 50 );
plot( dz, 1023, "dz", -1.5, 1.5, 60 );
@@ -643,16 +640,16 @@ StatusCode VertexCompare::finalize() {
<< ( (double)m_nVtx.nEntries() - (double)m_nVtx_geq10tr_matched2.nEntries() ) /
( (double)m_nRec2.nEntries() - (double)m_nRec_geq10tr_2.nEntries() ) * 100
<< " % " << endmsg;
-
- //const AIDA::IHistogram1D* dx = histo( HistoID( 1021 ) );
- //const AIDA::IHistogram1D* pullx = histo( HistoID( 1031 ) );
- //const AIDA::IHistogram1D* dy = histo( HistoID( 1022 ) );
- //const AIDA::IHistogram1D* pully = histo( HistoID( 1032 ) );
- //const AIDA::IHistogram1D* dz = histo( HistoID( 1023 ) );
- //const AIDA::IHistogram1D* pullz = histo( HistoID( 1033 ) );
+
+ const AIDA::IHistogram1D* dx = histo( HistoID( 1021 ) );
+ const AIDA::IHistogram1D* pullx = histo( HistoID( 1031 ) );
+ const AIDA::IHistogram1D* dy = histo( HistoID( 1022 ) );
+ const AIDA::IHistogram1D* pully = histo( HistoID( 1032 ) );
+ const AIDA::IHistogram1D* dz = histo( HistoID( 1023 ) );
+ const AIDA::IHistogram1D* pullz = histo( HistoID( 1033 ) );
+
- /*
if ( dx ) {
info() << " ---------------------------------------" << endmsg;
info() << "dx: "
@@ -660,6 +657,7 @@ StatusCode VertexCompare::finalize() {
Gaudi::Utils::HistoStats::meanErr( dx ), dx->rms(), Gaudi::Utils::HistoStats::rmsErr( dx ) )
<< endmsg;
}
+
if ( dy ) {
info() << "dy: "
<< format( "mean = %5.3f +/- %5.3f, RMS = %5.3f +/- %5.3f", dy->mean(),
@@ -714,7 +712,7 @@ StatusCode VertexCompare::finalize() {
<< endmsg;
}
info() << " ============================================" << endmsg;
- */
+
return Consumer::finalize(); // Must be called after all other actions
}
--
GitLab
From 80a8f59c174ab534036681b5b763178543f6342c Mon Sep 17 00:00:00 2001
From: Bogdan Kutsenko <bogdan.kutsenko@cern.ch>
Date: Thu, 25 Apr 2024 19:14:41 +0200
Subject: [PATCH 03/21] Add ntracks binning
---
Tr/TrackMonitors/src/VertexCompare.cpp | 62 +++++++++++++++++++-------
1 file changed, 47 insertions(+), 15 deletions(-)
diff --git a/Tr/TrackMonitors/src/VertexCompare.cpp b/Tr/TrackMonitors/src/VertexCompare.cpp
index 149d1635a29..619954b8bae 100644
--- a/Tr/TrackMonitors/src/VertexCompare.cpp
+++ b/Tr/TrackMonitors/src/VertexCompare.cpp
@@ -184,7 +184,7 @@ public:
{KeyValue{"inputVerticesName1", LHCb::Event::PV::DefaultLocation},
KeyValue{"inputVerticesName2", LHCb::Event::PV::DefaultLocation},
KeyValue{"InteractionRegionCache", "AlgorithmSpecific-" + name + "-InteractionRegion"}}} {}
-
+
StatusCode initialize() override; ///< Algorithm initialization
void operator()( Vertices const& recoVtx1, Vertices const& recoVtx2,
LHCb::Conditions::InteractionRegion const& ) const override; ///< Algorithm execution
@@ -193,13 +193,36 @@ public:
private:
bool debugLevel() const { return msgLevel( MSG::DEBUG ) || msgLevel( MSG::VERBOSE ); }
- using axis1D = Gaudi::Accumulators::Axis<Gaudi::Accumulators::Histogram<1>::AxisArithmeticType>;
- mutable Gaudi::Accumulators::Histogram<1> m_histo_dx{ this, "dx", "dx, mm", axis1D{50, -0.25, 0.25} };
- mutable Gaudi::Accumulators::Histogram<1> m_histo_dy{ this, "dy", "dy, mm", axis1D{50, -0.25, 0.25} };
- mutable Gaudi::Accumulators::Histogram<1> m_histo_dz{ this, "dz", "dz, mm", axis1D{50, -1.5, 1.5} };
- mutable Gaudi::Accumulators::Histogram<1> m_histo_pullx{ this, "pullx", "pull x", axis1D{20, -2, 2} };
- mutable Gaudi::Accumulators::Histogram<1> m_histo_pully{ this, "pully", "pull y", axis1D{20, -2, 2} };
- mutable Gaudi::Accumulators::Histogram<1> m_histo_pullz{ this, "pullz", "pull z", axis1D{20, -2, 2} };
+ std::vector<double> ntrack_bins = {4, 17.2, 30.4, 43.6, 56.8, 70};
+
+ struct monitoringHistos {
+ mutable std::array<Gaudi::Accumulators::Histogram<1>, 5> m_histo_dx;
+ mutable std::array<Gaudi::Accumulators::Histogram<1>, 5> m_histo_dy;
+ mutable std::array<Gaudi::Accumulators::Histogram<1>, 5> m_histo_dz;
+ mutable Gaudi::Accumulators::Histogram<1> m_histo_pullx;
+ mutable Gaudi::Accumulators::Histogram<1> m_histo_pully;
+ mutable Gaudi::Accumulators::Histogram<1> m_histo_pullz;
+
+ template <std::size_t... IDXs>
+ static std::array<Gaudi::Accumulators::Histogram<1>, sizeof...( IDXs )>
+ histo1DArrayBuilder( const VertexCompare* owner, const std::string& name, const std::string& title, std::tuple<unsigned, double, double> xbins, std::index_sequence<IDXs...> ) {
+ return {{{owner,
+ name + std::to_string( IDXs ),
+ title + std::to_string( IDXs ),
+ {std::get<0>( xbins ), std::get<1>( xbins ), std::get<2>( xbins )}}...}};
+ }
+ monitoringHistos( const VertexCompare* owner):
+ m_histo_dx{histo1DArrayBuilder( owner, "dx","dx, mm", {50, -0.25,0.25}, std::make_index_sequence<5>())},
+ m_histo_dy{histo1DArrayBuilder( owner, "dy","dy, mm", {50, -0.25,0.25},std::make_index_sequence<5>())},
+ m_histo_dz{histo1DArrayBuilder( owner, "dz","dz, mm", {50, -1.5,1.5},std::make_index_sequence<5>())},
+ m_histo_pullx{ owner, "pullx","pull x", {20, -2, 2}},
+ m_histo_pully{ owner, "pully","pull y", {20, -2, 2}},
+ m_histo_pullz{ owner, "pulz","pull z", {20, -2, 2}}
+ {}
+
+ };
+
+ std::unique_ptr<monitoringHistos> m_monitoringHistos;
Gaudi::Property<bool> m_produceHistogram{this, "produceHistogram", false};
Gaudi::Property<bool> m_produceNtuple{this, "produceNtuple", false};
@@ -232,6 +255,7 @@ StatusCode VertexCompare::initialize() {
return Consumer::initialize().andThen( [&]() {
+ if(m_monitoring.value()) m_monitoringHistos = std::make_unique<monitoringHistos>(this);
m_nVtx.reset();
m_nEvt.reset();
LHCb::Conditions::InteractionRegion::addConditionDerivation( this,
@@ -383,15 +407,23 @@ void VertexCompare::operator()( Vertices const& recoVtx1, Vertices const& recoVt
} else if ( ntracks2 > 1 ) {
multiplicity = link.at( oIt ) + 1;
}
-
+
if ( m_monitoring.value() ) {
if(dz<2 && size_diff < 3 && multiplicity < 2){
- ++m_histo_dx[dx];
- ++m_histo_dy[dy];
- ++m_histo_dz[dz];
- ++m_histo_pullx[dx / errx];
- ++m_histo_pully[dy / erry];
- ++m_histo_pullz[dz / errz];
+ int binCount = 0;
+ for (size_t i = 0; i < ntrack_bins.size(); ++i) {
+ if ((ntracks1+ntracks2)/2 <= ntrack_bins[i]) {
+ break;
+ }
+ binCount++;
+ }
+ auto& histos = *m_monitoringHistos.get();
+ ++histos.m_histo_dx[binCount][dx];
+ ++histos.m_histo_dy[binCount][dy];
+ ++histos.m_histo_dz[binCount][dz];
+ ++histos.m_histo_pullx[dx / errx];
+ ++histos.m_histo_pully[dy / erry];
+ ++histos.m_histo_pullz[dz / errz];
}
}
--
GitLab
From d381572c3669d66966f47cfb067cbfcef7625e37 Mon Sep 17 00:00:00 2001
From: gitlabCI <gitlab-ci@cern.ch>
Date: Fri, 26 Apr 2024 14:41:36 +0200
Subject: [PATCH 04/21] adding 'requireSingle' optional selection
---
Tr/TrackMonitors/src/VertexCompare.cpp | 258 +++++++++++++++----------
1 file changed, 154 insertions(+), 104 deletions(-)
diff --git a/Tr/TrackMonitors/src/VertexCompare.cpp b/Tr/TrackMonitors/src/VertexCompare.cpp
index 619954b8bae..79a8f8de179 100644
--- a/Tr/TrackMonitors/src/VertexCompare.cpp
+++ b/Tr/TrackMonitors/src/VertexCompare.cpp
@@ -21,7 +21,6 @@
#include "Event/Track_v2.h"
#include "GaudiAlg/GaudiAlgorithm.h"
#include "GaudiAlg/GaudiTupleAlg.h"
-#include <Gaudi/Accumulators/Histogram.h>
#include "GaudiAlg/Tuples.h"
#include "GaudiKernel/PhysicalConstants.h"
#include "GaudiKernel/SystemOfUnits.h"
@@ -30,6 +29,7 @@
#include "MCInterfaces/IForcedBDecayTool.h"
#include "VPDet/DeVP.h"
#include <Event/MCTrackInfo.h>
+#include <Gaudi/Accumulators/Histogram.h>
#include <LHCbDet/InteractionRegion.h>
#include <Linker/LinkedTo.h>
@@ -184,49 +184,49 @@ public:
{KeyValue{"inputVerticesName1", LHCb::Event::PV::DefaultLocation},
KeyValue{"inputVerticesName2", LHCb::Event::PV::DefaultLocation},
KeyValue{"InteractionRegionCache", "AlgorithmSpecific-" + name + "-InteractionRegion"}}} {}
-
+
StatusCode initialize() override; ///< Algorithm initialization
void operator()( Vertices const& recoVtx1, Vertices const& recoVtx2,
LHCb::Conditions::InteractionRegion const& ) const override; ///< Algorithm execution
StatusCode finalize() override; ///< Algorithm finalization
private:
- bool debugLevel() const { return msgLevel( MSG::DEBUG ) || msgLevel( MSG::VERBOSE ); }
+ bool debugLevel() const { return msgLevel( MSG::DEBUG ) || msgLevel( MSG::VERBOSE ); }
std::vector<double> ntrack_bins = {4, 17.2, 30.4, 43.6, 56.8, 70};
-
- struct monitoringHistos {
- mutable std::array<Gaudi::Accumulators::Histogram<1>, 5> m_histo_dx;
- mutable std::array<Gaudi::Accumulators::Histogram<1>, 5> m_histo_dy;
- mutable std::array<Gaudi::Accumulators::Histogram<1>, 5> m_histo_dz;
- mutable Gaudi::Accumulators::Histogram<1> m_histo_pullx;
- mutable Gaudi::Accumulators::Histogram<1> m_histo_pully;
- mutable Gaudi::Accumulators::Histogram<1> m_histo_pullz;
-
- template <std::size_t... IDXs>
- static std::array<Gaudi::Accumulators::Histogram<1>, sizeof...( IDXs )>
- histo1DArrayBuilder( const VertexCompare* owner, const std::string& name, const std::string& title, std::tuple<unsigned, double, double> xbins, std::index_sequence<IDXs...> ) {
- return {{{owner,
- name + std::to_string( IDXs ),
- title + std::to_string( IDXs ),
- {std::get<0>( xbins ), std::get<1>( xbins ), std::get<2>( xbins )}}...}};
- }
- monitoringHistos( const VertexCompare* owner):
- m_histo_dx{histo1DArrayBuilder( owner, "dx","dx, mm", {50, -0.25,0.25}, std::make_index_sequence<5>())},
- m_histo_dy{histo1DArrayBuilder( owner, "dy","dy, mm", {50, -0.25,0.25},std::make_index_sequence<5>())},
- m_histo_dz{histo1DArrayBuilder( owner, "dz","dz, mm", {50, -1.5,1.5},std::make_index_sequence<5>())},
- m_histo_pullx{ owner, "pullx","pull x", {20, -2, 2}},
- m_histo_pully{ owner, "pully","pull y", {20, -2, 2}},
- m_histo_pullz{ owner, "pulz","pull z", {20, -2, 2}}
- {}
-
+
+ struct monitoringHistos {
+ mutable std::array<Gaudi::Accumulators::Histogram<1>, 5> m_histo_dx;
+ mutable std::array<Gaudi::Accumulators::Histogram<1>, 5> m_histo_dy;
+ mutable std::array<Gaudi::Accumulators::Histogram<1>, 5> m_histo_dz;
+ mutable Gaudi::Accumulators::Histogram<1> m_histo_pullx;
+ mutable Gaudi::Accumulators::Histogram<1> m_histo_pully;
+ mutable Gaudi::Accumulators::Histogram<1> m_histo_pullz;
+
+ template <std::size_t... IDXs>
+ static std::array<Gaudi::Accumulators::Histogram<1>, sizeof...( IDXs )>
+ histo1DArrayBuilder( const VertexCompare* owner, const std::string& name, const std::string& title,
+ std::tuple<unsigned, double, double> xbins, std::index_sequence<IDXs...> ) {
+ return {{{owner,
+ name + std::to_string( IDXs ),
+ title + std::to_string( IDXs ),
+ {std::get<0>( xbins ), std::get<1>( xbins ), std::get<2>( xbins )}}...}};
+ }
+ monitoringHistos( const VertexCompare* owner )
+ : m_histo_dx{histo1DArrayBuilder( owner, "dx", "dx, mm", {50, -0.25, 0.25}, std::make_index_sequence<5>() )}
+ , m_histo_dy{histo1DArrayBuilder( owner, "dy", "dy, mm", {50, -0.25, 0.25}, std::make_index_sequence<5>() )}
+ , m_histo_dz{histo1DArrayBuilder( owner, "dz", "dz, mm", {50, -1.5, 1.5}, std::make_index_sequence<5>() )}
+ , m_histo_pullx{owner, "pullx", "pull x", {20, -2, 2}}
+ , m_histo_pully{owner, "pully", "pull y", {20, -2, 2}}
+ , m_histo_pullz{owner, "pulz", "pull z", {20, -2, 2}} {}
};
std::unique_ptr<monitoringHistos> m_monitoringHistos;
-
+
Gaudi::Property<bool> m_produceHistogram{this, "produceHistogram", false};
Gaudi::Property<bool> m_produceNtuple{this, "produceNtuple", false};
Gaudi::Property<bool> m_monitoring{this, "monitoring", true};
+ Gaudi::Property<bool> m_requireSingle{this, "requireSingle", false};
mutable Gaudi::Accumulators::Counter<> m_nVtx{this, "Number of pairs of vertices in processed events"};
mutable Gaudi::Accumulators::Counter<> m_nRec1{this, "Number of vertices in input1"};
@@ -252,10 +252,9 @@ DECLARE_COMPONENT_WITH_ID( VertexCompare, "VertexCompare" )
// Initialization
//=============================================================================
StatusCode VertexCompare::initialize() {
-
return Consumer::initialize().andThen( [&]() {
- if(m_monitoring.value()) m_monitoringHistos = std::make_unique<monitoringHistos>(this);
+ if ( m_monitoring.value() ) m_monitoringHistos = std::make_unique<monitoringHistos>( this );
m_nVtx.reset();
m_nEvt.reset();
LHCb::Conditions::InteractionRegion::addConditionDerivation( this,
@@ -277,59 +276,114 @@ void VertexCompare::operator()( Vertices const& recoVtx1, Vertices const& recoVt
if ( debugLevel() ) debug() << " Vtx Properities x y z chi2/ndof ntracks" << endmsg;
- for ( auto const& pv : recoVtx1 ) {
- RecPVInfo<VertexType> recinfo1;
- ++m_nRec1;
- recinfo1.pRECPV = &pv;
- recinfo1.position = pv.position();
- recinfo1.covPV = pv.covMatrix();
- recinfo1.positionSigma = Gaudi::XYZPoint{std::sqrt( recinfo1.covPV( 0, 0 ) ), std::sqrt( recinfo1.covPV( 1, 1 ) ),
- std::sqrt( recinfo1.covPV( 2, 2 ) )};
- recinfo1.ntracks = pv.nTracks();
- if ( recinfo1.ntracks >= 10 ) { ++m_nRec_geq10tr_1; }
- recinfo1.chi2 = pv.chi2();
- recinfo1.nDoF = pv.nDoF();
+ int size_diff = -9999;
+
+ if ( m_requireSingle == false ) {
+ for ( auto const& pv : recoVtx1 ) {
+ RecPVInfo<VertexType> recinfo1;
+ ++m_nRec1;
+ recinfo1.pRECPV = &pv;
+ recinfo1.position = pv.position();
+ recinfo1.covPV = pv.covMatrix();
+ recinfo1.positionSigma = Gaudi::XYZPoint{std::sqrt( recinfo1.covPV( 0, 0 ) ), std::sqrt( recinfo1.covPV( 1, 1 ) ),
+ std::sqrt( recinfo1.covPV( 2, 2 ) )};
+ recinfo1.ntracks = pv.nTracks();
+ if ( recinfo1.ntracks >= 10 ) { ++m_nRec_geq10tr_1; }
+ recinfo1.chi2 = pv.chi2();
+ recinfo1.nDoF = pv.nDoF();
+
+ if ( debugLevel() )
+ debug() << " " << pv.position().x() << " " << pv.position().y() << " " << pv.position().z()
+ << " " << pv.chi2PerDoF() << " " << pv.nTracks() << endmsg;
+ fullVrt1.push_back( recinfo1 );
+ } // end of loop over vertices1
+
+ for ( auto const& pv : recoVtx2 ) {
+ RecPVInfo<VertexType> recinfo2;
+ ++m_nRec2;
+ recinfo2.pRECPV = &pv;
+ recinfo2.position = pv.position();
+ recinfo2.covPV = pv.covMatrix();
+ recinfo2.positionSigma = Gaudi::XYZPoint{std::sqrt( recinfo2.covPV( 0, 0 ) ), std::sqrt( recinfo2.covPV( 1, 1 ) ),
+ std::sqrt( recinfo2.covPV( 2, 2 ) )};
+ recinfo2.ntracks = pv.nTracks();
+ if ( recinfo2.ntracks >= 10 ) { ++m_nRec_geq10tr_2; }
+ recinfo2.chi2 = pv.chi2();
+ recinfo2.nDoF = pv.nDoF();
+
+ if ( debugLevel() )
+ debug() << " " << pv.position().x() << " " << pv.position().y() << " " << pv.position().z()
+ << " " << pv.chi2PerDoF() << " " << pv.nTracks() << endmsg;
+ fullVrt2.push_back( recinfo2 );
+ } // end of loop over vertices2
+
+ // added sorting to mirror the 'multirec' variable implementation in PVChecker
+ std::sort( fullVrt1.begin(), fullVrt1.end(), trackcomp );
+ std::sort( fullVrt2.begin(), fullVrt2.end(), trackcomp );
+
+ if ( debugLevel() ) debug() << "fullVrt1 size " << fullVrt1.size() << endmsg;
+ if ( debugLevel() ) debug() << "fullVrt2 size " << fullVrt2.size() << endmsg;
+ size_diff = fullVrt1.size() - fullVrt2.size();
+ } else if ( m_requireSingle == true && recoVtx1.size() == 1 && recoVtx2.size() == 1 ) {
+ for ( auto const& pv : recoVtx1 ) {
+ RecPVInfo<VertexType> recinfo1;
+ ++m_nRec1;
+ recinfo1.pRECPV = &pv;
+ recinfo1.position = pv.position();
+ recinfo1.covPV = pv.covMatrix();
+ recinfo1.positionSigma = Gaudi::XYZPoint{std::sqrt( recinfo1.covPV( 0, 0 ) ), std::sqrt( recinfo1.covPV( 1, 1 ) ),
+ std::sqrt( recinfo1.covPV( 2, 2 ) )};
+ recinfo1.ntracks = pv.nTracks();
+ if ( recinfo1.ntracks >= 10 ) { ++m_nRec_geq10tr_1; }
+ recinfo1.chi2 = pv.chi2();
+ recinfo1.nDoF = pv.nDoF();
+
+ if ( debugLevel() )
+ debug() << " " << pv.position().x() << " " << pv.position().y() << " " << pv.position().z()
+ << " " << pv.chi2PerDoF() << " " << pv.nTracks() << endmsg;
+ fullVrt1.push_back( recinfo1 );
+ } // end of loop over vertices1
+
+ for ( auto const& pv : recoVtx2 ) {
+ RecPVInfo<VertexType> recinfo2;
+ ++m_nRec2;
+ recinfo2.pRECPV = &pv;
+ recinfo2.position = pv.position();
+ recinfo2.covPV = pv.covMatrix();
+ recinfo2.positionSigma = Gaudi::XYZPoint{std::sqrt( recinfo2.covPV( 0, 0 ) ), std::sqrt( recinfo2.covPV( 1, 1 ) ),
+ std::sqrt( recinfo2.covPV( 2, 2 ) )};
+ recinfo2.ntracks = pv.nTracks();
+ if ( recinfo2.ntracks >= 10 ) { ++m_nRec_geq10tr_2; }
+ recinfo2.chi2 = pv.chi2();
+ recinfo2.nDoF = pv.nDoF();
+
+ if ( debugLevel() )
+ debug() << " " << pv.position().x() << " " << pv.position().y() << " " << pv.position().z()
+ << " " << pv.chi2PerDoF() << " " << pv.nTracks() << endmsg;
+ fullVrt2.push_back( recinfo2 );
+ } // end of loop over vertices2
+
+ // sorting is unnecessary for singular vertices
if ( debugLevel() )
- debug() << " " << pv.position().x() << " " << pv.position().y() << " " << pv.position().z()
- << " " << pv.chi2PerDoF() << " " << pv.nTracks() << endmsg;
- fullVrt1.push_back( recinfo1 );
- } // end of loop over vertices1
-
- for ( auto const& pv : recoVtx2 ) {
- RecPVInfo<VertexType> recinfo2;
- ++m_nRec2;
- recinfo2.pRECPV = &pv;
- recinfo2.position = pv.position();
- recinfo2.covPV = pv.covMatrix();
- recinfo2.positionSigma = Gaudi::XYZPoint{std::sqrt( recinfo2.covPV( 0, 0 ) ), std::sqrt( recinfo2.covPV( 1, 1 ) ),
- std::sqrt( recinfo2.covPV( 2, 2 ) )};
- recinfo2.ntracks = pv.nTracks();
- if ( recinfo2.ntracks >= 10 ) { ++m_nRec_geq10tr_2; }
- recinfo2.chi2 = pv.chi2();
- recinfo2.nDoF = pv.nDoF();
-
+ debug() << "fullVrt1 size "
+ << "1" << endmsg;
if ( debugLevel() )
- debug() << " " << pv.position().x() << " " << pv.position().y() << " " << pv.position().z()
- << " " << pv.chi2PerDoF() << " " << pv.nTracks() << endmsg;
- fullVrt2.push_back( recinfo2 );
- } // end of loop over vertices2
-
- // added sorting to mirror the 'multirec' variable implementation in PVChecker
- std::sort( fullVrt1.begin(), fullVrt1.end(), trackcomp );
- std::sort( fullVrt2.begin(), fullVrt2.end(), trackcomp );
-
- if ( debugLevel() ) debug() << "fullVrt1 size " << fullVrt1.size() << endmsg;
- if ( debugLevel() ) debug() << "fullVrt2 size " << fullVrt2.size() << endmsg;
- int size_diff = fullVrt1.size() - fullVrt2.size();
-
- if ( m_produceHistogram.value() ) { plot( double( size_diff ), 1001, "size_diff", -5.5, 5.5, 11 ); }
- if ( m_produceNtuple.value() ) {
- Tuple myTuple_evt = nTuple( "PV_nTuple_evt", "PV_nTuple_evt", CLID_ColumnWiseTuple );
- myTuple_evt->column( "size_diff", double( size_diff ) ).ignore();
- myTuple_evt->column( "size_1", double( fullVrt1.size() ) ).ignore();
- myTuple_evt->column( "size_2", double( fullVrt2.size() ) ).ignore();
- myTuple_evt->write().ignore();
+ debug() << "fullVrt2 size "
+ << "1" << endmsg;
+ size_diff = 0;
+ }
+
+ // this check takes care of events which evade the selection
+ if ( !( m_requireSingle == true && ( recoVtx1.size() != 1 || recoVtx2.size() != 1 ) ) ) {
+ if ( m_produceHistogram.value() ) { plot( double( size_diff ), 1001, "size_diff", -5.5, 5.5, 11 ); }
+ if ( m_produceNtuple.value() ) {
+ Tuple myTuple_evt = nTuple( "PV_nTuple_evt", "PV_nTuple_evt", CLID_ColumnWiseTuple );
+ myTuple_evt->column( "size_diff", double( size_diff ) ).ignore();
+ myTuple_evt->column( "size_1", double( fullVrt1.size() ) ).ignore();
+ myTuple_evt->column( "size_2", double( fullVrt2.size() ) ).ignore();
+ myTuple_evt->write().ignore();
+ }
}
std::vector<int> link;
@@ -407,26 +461,24 @@ void VertexCompare::operator()( Vertices const& recoVtx1, Vertices const& recoVt
} else if ( ntracks2 > 1 ) {
multiplicity = link.at( oIt ) + 1;
}
-
+
if ( m_monitoring.value() ) {
- if(dz<2 && size_diff < 3 && multiplicity < 2){
- int binCount = 0;
- for (size_t i = 0; i < ntrack_bins.size(); ++i) {
- if ((ntracks1+ntracks2)/2 <= ntrack_bins[i]) {
- break;
- }
- binCount++;
+ if ( dz < 2 && size_diff < 3 && multiplicity < 2 ) {
+ int binCount = 0;
+ for ( size_t i = 0; i < ntrack_bins.size(); ++i ) {
+ if ( ( ntracks1 + ntracks2 ) / 2 <= ntrack_bins[i] ) { break; }
+ binCount++;
}
- auto& histos = *m_monitoringHistos.get();
- ++histos.m_histo_dx[binCount][dx];
- ++histos.m_histo_dy[binCount][dy];
- ++histos.m_histo_dz[binCount][dz];
- ++histos.m_histo_pullx[dx / errx];
- ++histos.m_histo_pully[dy / erry];
- ++histos.m_histo_pullz[dz / errz];
- }
+ auto& histos = *m_monitoringHistos.get();
+ ++histos.m_histo_dx[binCount][dx];
+ ++histos.m_histo_dy[binCount][dy];
+ ++histos.m_histo_dz[binCount][dz];
+ ++histos.m_histo_pullx[dx / errx];
+ ++histos.m_histo_pully[dy / erry];
+ ++histos.m_histo_pullz[dz / errz];
+ }
}
-
+
if ( m_produceHistogram.value() ) {
plot( dx, 1021, "dx", -0.25, 0.25, 50 );
plot( dy, 1022, "dy", -0.25, 0.25, 50 );
@@ -673,15 +725,13 @@ StatusCode VertexCompare::finalize() {
( (double)m_nRec2.nEntries() - (double)m_nRec_geq10tr_2.nEntries() ) * 100
<< " % " << endmsg;
-
const AIDA::IHistogram1D* dx = histo( HistoID( 1021 ) );
const AIDA::IHistogram1D* pullx = histo( HistoID( 1031 ) );
const AIDA::IHistogram1D* dy = histo( HistoID( 1022 ) );
const AIDA::IHistogram1D* pully = histo( HistoID( 1032 ) );
const AIDA::IHistogram1D* dz = histo( HistoID( 1023 ) );
const AIDA::IHistogram1D* pullz = histo( HistoID( 1033 ) );
-
-
+
if ( dx ) {
info() << " ---------------------------------------" << endmsg;
info() << "dx: "
@@ -689,7 +739,7 @@ StatusCode VertexCompare::finalize() {
Gaudi::Utils::HistoStats::meanErr( dx ), dx->rms(), Gaudi::Utils::HistoStats::rmsErr( dx ) )
<< endmsg;
}
-
+
if ( dy ) {
info() << "dy: "
<< format( "mean = %5.3f +/- %5.3f, RMS = %5.3f +/- %5.3f", dy->mean(),
@@ -744,7 +794,7 @@ StatusCode VertexCompare::finalize() {
<< endmsg;
}
info() << " ============================================" << endmsg;
-
+
return Consumer::finalize(); // Must be called after all other actions
}
--
GitLab
From 26d523164ba5d201cd1ff6b991c2e9c3227abaf0 Mon Sep 17 00:00:00 2001
From: Gerhard Raven <gerhard.raven@nikhef.nl>
Date: Fri, 26 Apr 2024 14:55:12 +0200
Subject: [PATCH 05/21] Apply 1 suggestion(s) to 1 file(s)
---
Tr/TrackMonitors/src/VertexCompare.cpp | 2 +-
1 file changed, 1 insertion(+), 1 deletion(-)
diff --git a/Tr/TrackMonitors/src/VertexCompare.cpp b/Tr/TrackMonitors/src/VertexCompare.cpp
index 79a8f8de179..853c5dde0e2 100644
--- a/Tr/TrackMonitors/src/VertexCompare.cpp
+++ b/Tr/TrackMonitors/src/VertexCompare.cpp
@@ -193,7 +193,7 @@ public:
private:
bool debugLevel() const { return msgLevel( MSG::DEBUG ) || msgLevel( MSG::VERBOSE ); }
- std::vector<double> ntrack_bins = {4, 17.2, 30.4, 43.6, 56.8, 70};
+ static constexpr auto ntrack_bins = std::array{4., 17.2, 30.4, 43.6, 56.8, 70.};
struct monitoringHistos {
mutable std::array<Gaudi::Accumulators::Histogram<1>, 5> m_histo_dx;
--
GitLab
From 6f935d6cd80c0aed789e711f15d06993386ab4da Mon Sep 17 00:00:00 2001
From: Bogdan Kutsenko <bogdan.kutsenko@cern.ch>
Date: Mon, 29 Apr 2024 13:49:17 +0200
Subject: [PATCH 06/21] Vertex compare cleaned
---
Tr/TrackMonitors/src/VertexCompare.cpp | 10 +++++-----
1 file changed, 5 insertions(+), 5 deletions(-)
diff --git a/Tr/TrackMonitors/src/VertexCompare.cpp b/Tr/TrackMonitors/src/VertexCompare.cpp
index 853c5dde0e2..a02f9028f5f 100644
--- a/Tr/TrackMonitors/src/VertexCompare.cpp
+++ b/Tr/TrackMonitors/src/VertexCompare.cpp
@@ -223,9 +223,9 @@ private:
std::unique_ptr<monitoringHistos> m_monitoringHistos;
- Gaudi::Property<bool> m_produceHistogram{this, "produceHistogram", false};
- Gaudi::Property<bool> m_produceNtuple{this, "produceNtuple", false};
- Gaudi::Property<bool> m_monitoring{this, "monitoring", true};
+ Gaudi::Property<bool> m_produceHistogram{this, "produceHistogram", true};
+ Gaudi::Property<bool> m_produceNtuple{this, "produceNtuple", true};
+ Gaudi::Property<bool> m_monitoring{this, "monitoring", false};
Gaudi::Property<bool> m_requireSingle{this, "requireSingle", false};
mutable Gaudi::Accumulators::Counter<> m_nVtx{this, "Number of pairs of vertices in processed events"};
@@ -277,10 +277,11 @@ void VertexCompare::operator()( Vertices const& recoVtx1, Vertices const& recoVt
if ( debugLevel() ) debug() << " Vtx Properities x y z chi2/ndof ntracks" << endmsg;
int size_diff = -9999;
+ RecPVInfo<VertexType> recinfo1;
+ RecPVInfo<VertexType> recinfo2;
if ( m_requireSingle == false ) {
for ( auto const& pv : recoVtx1 ) {
- RecPVInfo<VertexType> recinfo1;
++m_nRec1;
recinfo1.pRECPV = &pv;
recinfo1.position = pv.position();
@@ -299,7 +300,6 @@ void VertexCompare::operator()( Vertices const& recoVtx1, Vertices const& recoVt
} // end of loop over vertices1
for ( auto const& pv : recoVtx2 ) {
- RecPVInfo<VertexType> recinfo2;
++m_nRec2;
recinfo2.pRECPV = &pv;
recinfo2.position = pv.position();
--
GitLab
From 4f324fb85261f6f75300348057e0051fbb4339cb Mon Sep 17 00:00:00 2001
From: Bogdan Kutsenko <bogdan.kutsenko@cern.ch>
Date: Tue, 7 May 2024 10:24:31 +0200
Subject: [PATCH 07/21] VertexCompare update
---
Tr/TrackMonitors/src/VertexCompare.cpp | 8 ++++----
1 file changed, 4 insertions(+), 4 deletions(-)
diff --git a/Tr/TrackMonitors/src/VertexCompare.cpp b/Tr/TrackMonitors/src/VertexCompare.cpp
index a02f9028f5f..d71efaee88e 100644
--- a/Tr/TrackMonitors/src/VertexCompare.cpp
+++ b/Tr/TrackMonitors/src/VertexCompare.cpp
@@ -193,7 +193,7 @@ public:
private:
bool debugLevel() const { return msgLevel( MSG::DEBUG ) || msgLevel( MSG::VERBOSE ); }
- static constexpr auto ntrack_bins = std::array{4., 17.2, 30.4, 43.6, 56.8, 70.};
+ static constexpr auto ntrack_bins = std::array{4. , 7.2, 10.4, 13.6, 16.8, 20.};
struct monitoringHistos {
mutable std::array<Gaudi::Accumulators::Histogram<1>, 5> m_histo_dx;
@@ -463,10 +463,10 @@ void VertexCompare::operator()( Vertices const& recoVtx1, Vertices const& recoVt
}
if ( m_monitoring.value() ) {
- if ( dz < 2 && size_diff < 3 && multiplicity < 2 ) {
+ //if ( dz < 2 && size_diff < 3) {
int binCount = 0;
for ( size_t i = 0; i < ntrack_bins.size(); ++i ) {
- if ( ( ntracks1 + ntracks2 ) / 2 <= ntrack_bins[i] ) { break; }
+ if ( std::lround(( ntracks1 + ntracks2 ) / 2) <= ntrack_bins[i] ) { break; }
binCount++;
}
auto& histos = *m_monitoringHistos.get();
@@ -476,7 +476,7 @@ void VertexCompare::operator()( Vertices const& recoVtx1, Vertices const& recoVt
++histos.m_histo_pullx[dx / errx];
++histos.m_histo_pully[dy / erry];
++histos.m_histo_pullz[dz / errz];
- }
+ // } if end
}
if ( m_produceHistogram.value() ) {
--
GitLab
From f37fdc1d32317021141961fd09fb9fea1ec639c5 Mon Sep 17 00:00:00 2001
From: gitlabCI <gitlab-ci@cern.ch>
Date: Thu, 4 Apr 2024 20:08:27 +0200
Subject: [PATCH 08/21] random track splitter
---
Tr/TrackUtils/CMakeLists.txt | 1 +
.../src/RandomTrackContainerSplitter.cpp | 103 ++++++++++++++++++
2 files changed, 104 insertions(+)
create mode 100644 Tr/TrackUtils/src/RandomTrackContainerSplitter.cpp
diff --git a/Tr/TrackUtils/CMakeLists.txt b/Tr/TrackUtils/CMakeLists.txt
index 35e937a382f..2e523f5f873 100644
--- a/Tr/TrackUtils/CMakeLists.txt
+++ b/Tr/TrackUtils/CMakeLists.txt
@@ -26,6 +26,7 @@ gaudi_add_module(TrackUtils
src/TrackCompetition.cpp
src/TrackContainerCopy.cpp
src/TrackContainerSplitter.cpp
+ src/RandomTrackContainerSplitter.cpp
src/TrackContainersMerger.cpp
src/TrackListFilter.cpp
src/TrackListMerger.cpp
diff --git a/Tr/TrackUtils/src/RandomTrackContainerSplitter.cpp b/Tr/TrackUtils/src/RandomTrackContainerSplitter.cpp
new file mode 100644
index 00000000000..248e4ab35ed
--- /dev/null
+++ b/Tr/TrackUtils/src/RandomTrackContainerSplitter.cpp
@@ -0,0 +1,103 @@
+/*****************************************************************************\
+* (c) Copyright 2000-2020 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. *
+\*****************************************************************************/
+
+/** @class RandomTrackContainerSplitter RandomTrackContainerSplitter.h
+ *
+ * Split container to two (passed and rejected) wrt to given selection criteria
+ *
+ * @author Andrii Usachov
+ * @date 30/04/2020
+ */
+
+#include "Event/Track.h"
+#include "Functors/Function.h"
+#include "Functors/with_functors.h"
+#include "GaudiKernel/IEventTimeDecoder.h"
+#include "GaudiKernel/IRndmEngine.h"
+#include "GaudiKernel/IRndmGenSvc.h"
+#include "GaudiKernel/RndmGenerators.h"
+#include "GaudiKernel/Service.h" // Include the header file for Service
+#include "GaudiKernel/ServiceHandle.h"
+#include "LHCbAlgs/Transformer.h"
+
+namespace {
+ /// The output data
+ using OutConts = std::tuple<LHCb::Tracks, LHCb::Tracks>;
+
+ struct TrackPredicate {
+ using Signature = bool( const LHCb::Track& );
+ static constexpr auto PropertyName = "Code";
+ };
+} // namespace
+
+// bool my_splitter(auto trk){
+// Rndm::Numbers m_rnd;
+// IRndmGenSvc* rndm = svc<IRndmGenSvc>( "RndmGenSvc", true );
+// StatusCode sc = m_rnd.initialize( rndm, Rndm::Flat( 0.0, 100 ) );
+// if ( sc.isFailure() )
+// throw GaudiException{"Unable to initialize Flat distribution", "addConditionDerivation", StatusCode::FAILURE};
+// if ( m_rnd < 0.5 )
+// {
+// return false;
+// }
+// else
+// {
+// return true;
+// }
+// }
+
+bool makeRandomDecision( const LHCb::Track& track, ServiceHandle<IRndmGenSvc>& rndmSvc ) {
+ if ( !rndmSvc.retrieve().isSuccess() ) {
+ std::cerr << "Error: Failed to retrieve IRndmGenSvc service." << std::endl;
+ return false;
+ }
+
+ IRndmGenSvc* rndmGenSvcPtr = rndmSvc.operator->();
+
+ if ( !rndmGenSvcPtr ) {
+ std::cerr << "Error: IRndmGenSvc pointer is null." << std::endl;
+ return false;
+ }
+
+ Rndm::Numbers rndm( rndmGenSvcPtr, Rndm::Flat( 0., 1.0 ) );
+ double randomNumber = rndm();
+
+ return randomNumber < 0.5;
+}
+
+class RandomTrackContainerSplitter final
+ : public with_functors<LHCb::Algorithm::MultiTransformer<OutConts( const LHCb::Tracks& )>, TrackPredicate> {
+public:
+ RandomTrackContainerSplitter( const std::string& name, ISvcLocator* pSvcLocator )
+ : with_functors( name, pSvcLocator, {KeyValue{"TracksInContainer", ""}},
+ {KeyValue{"PassedContainer", ""}, KeyValue{"RejectedContainer", ""}} ) {}
+
+ OutConts operator()( const LHCb::Tracks& input_tracks ) const override {
+ OutConts out;
+ ServiceHandle<IRndmGenSvc> rndmSvc( "RndmGenSvc", "RandomTrackContainerSplitter" );
+ auto& [passed, rejected] = out;
+ m_inputs += input_tracks.size();
+ // auto const& pred = my_splitter<TrackPredicate>();
+ for ( auto& trk : input_tracks ) {
+ auto& container = ( makeRandomDecision( *trk, rndmSvc ) ? passed : rejected );
+ container.insert( new LHCb::Track( *trk ) );
+ }
+ m_passed += passed.size();
+
+ return out;
+ }
+
+private:
+ mutable Gaudi::Accumulators::StatCounter<> m_inputs{this, "#inputs"};
+ mutable Gaudi::Accumulators::StatCounter<> m_passed{this, "#passed"};
+};
+
+DECLARE_COMPONENT( RandomTrackContainerSplitter )
--
GitLab
From 03520ad17e1de23a60da0d2b78bd32c140c37ca9 Mon Sep 17 00:00:00 2001
From: Bogdan Kutsenko <bogdan.kutsenko@cern.ch>
Date: Fri, 24 May 2024 16:14:37 +0200
Subject: [PATCH 09/21] Transition to Gaudi::Accumulator::Histogram
---
Tr/TrackMonitors/src/VertexCompare.cpp | 408 +++++++++++++------------
1 file changed, 209 insertions(+), 199 deletions(-)
diff --git a/Tr/TrackMonitors/src/VertexCompare.cpp b/Tr/TrackMonitors/src/VertexCompare.cpp
index d71efaee88e..1e42189e092 100644
--- a/Tr/TrackMonitors/src/VertexCompare.cpp
+++ b/Tr/TrackMonitors/src/VertexCompare.cpp
@@ -8,7 +8,6 @@
* granted to it by virtue of its status as an Intergovernmental Organization *
* or submit itself to any jurisdiction. *
\*****************************************************************************/
-#include "AIDA/IHistogram1D.h"
#include "Event/MCHeader.h"
#include "Event/MCParticle.h"
#include "Event/MCProperty.h"
@@ -30,9 +29,9 @@
#include "VPDet/DeVP.h"
#include <Event/MCTrackInfo.h>
#include <Gaudi/Accumulators/Histogram.h>
+#include <Gaudi/Accumulators/RootHistogram.h>
#include <LHCbDet/InteractionRegion.h>
#include <Linker/LinkedTo.h>
-
using Vertices = LHCb::Event::PV::PrimaryVertexContainer;
using VertexType = Vertices::value_type;
@@ -111,10 +110,11 @@ struct VtxVariables {
double nDoF = -9999.;
int dsize = -9999;
int match = -9999;
- int multiplicity = -9999;
+ int pv_rank = -9999;
bool equal_sizes = false;
bool single = false;
bool opposite_container = false;
+ std::vector<double> dx_vector{}, dy_vector{}, dz_vector{}, pullx_vector{}, pully_vector{}, pullz_vector{};
};
VtxVariables SetVtxVariables( const RecPVInfo<VertexType>& vrtf, const int& counter, const int& size1, const int& size2,
@@ -152,9 +152,9 @@ VtxVariables SetVtxVariables( const RecPVInfo<VertexType>& vrtf, const int& coun
vtx.ntracks = vrtf.ntracks;
if ( vtx.ntracks > 1 ) {
- vtx.multiplicity = counter + 1;
+ vtx.pv_rank = counter + 1;
} else {
- vtx.multiplicity = -99999;
+ vtx.pv_rank = -99999;
}
vtx.dsize = size1 - size2;
vtx.equal_sizes = ( size1 == size2 );
@@ -172,6 +172,76 @@ VtxVariables SetVtxVariables( const RecPVInfo<VertexType>& vrtf, const int& coun
return vtx;
}
+// ============================================================================
+// get an error in the mean value
+// ============================================================================
+
+template <typename Histogram>
+double meanErr( Histogram const & h ) {
+ const auto n = h.nEntries();
+ return ( 0 >= n ? 0.0 : h.standard_deviation() / std::sqrt( n ) );
+}
+
+// ============================================================================
+// Fourth momentum calculation
+// ============================================================================
+template <typename Histogram>
+double fourthMomentum( Histogram const & h ) {
+ const auto n = h.nEntries();
+ if ( 0 >= n ) { return 0.0; } // RETURN
+ // get the mean
+ const auto h_mean = h.mean();
+ // number of bins
+ const auto& axis = h.axis();
+ const int nBins = h.totNBins();
+ double result{ 0 }, weight{ 0 };
+
+ // loop over all bins
+ double minVal = axis[0].minValue;
+ double maxVal = axis[0].maxValue;
+ double binWidth = (maxVal-minVal)/nBins;
+
+ std::vector<double> binCenters;
+ for (int i = 0; i < nBins; ++i) {
+ double center = minVal + (i + 0.5) * binWidth;
+ binCenters.push_back(center);
+ }
+ for (int i = 0; i < nBins; ++i) {
+ const auto yBin = h.binValue( i ); // bin weight
+ weight += yBin;
+ result += yBin * std::pow( binCenters[i] - h_mean, 4 );
+ }
+ if ( weight > std::numeric_limits<double>::epsilon() ) { result /= weight; }
+
+ return result;
+}
+
+// ============================================================================
+// get the error in kurtosis
+// ============================================================================
+
+template <typename Histogram>
+double kurtosis( Histogram const & h ) {
+ const auto mu4 = fourthMomentum( h );
+ const auto s4 = std::pow( h.standard_deviation() , 4 );
+ return ( std::fabs( s4 ) > 0 ? mu4 / s4 - 3.0 : 0.0 );
+}
+
+// ============================================================================
+// get an error in the rms value
+// ============================================================================
+
+ template <typename Histogram>
+ double rmsErr( Histogram const & h ) {
+ const auto n = h.nEntries();
+ if ( 1 >= n ) { return 0.0; }
+ auto result = 2.0 + kurtosis( h );
+ result += 2.0 / ( n - 1 );
+ result /= 4.0 * n;
+ return h.standard_deviation() * std::sqrt( std::max( result, 0.0 ) );
+ }
+
+
class VertexCompare : public LHCb::Algorithm::Consumer<
void( Vertices const&, Vertices const&, LHCb::Conditions::InteractionRegion const& ),
LHCb::DetDesc::usesBaseAndConditions<GaudiTupleAlg, LHCb::Conditions::InteractionRegion>> {
@@ -193,16 +263,24 @@ public:
private:
bool debugLevel() const { return msgLevel( MSG::DEBUG ) || msgLevel( MSG::VERBOSE ); }
- static constexpr auto ntrack_bins = std::array{4. , 7.2, 10.4, 13.6, 16.8, 20.};
+ static constexpr auto ntrack_bins = std::array{ 3.5 , 16.5 , 29.5 , 43.5 , 56.5 , 69.5 };
- struct monitoringHistos {
- mutable std::array<Gaudi::Accumulators::Histogram<1>, 5> m_histo_dx;
- mutable std::array<Gaudi::Accumulators::Histogram<1>, 5> m_histo_dy;
- mutable std::array<Gaudi::Accumulators::Histogram<1>, 5> m_histo_dz;
- mutable Gaudi::Accumulators::Histogram<1> m_histo_pullx;
- mutable Gaudi::Accumulators::Histogram<1> m_histo_pully;
- mutable Gaudi::Accumulators::Histogram<1> m_histo_pullz;
+ Gaudi::Property<bool> m_produceHistogram{this, "produceHistogram", true};
+ Gaudi::Property<bool> m_produceNtuple{this, "produceNtuple", true};
+ Gaudi::Property<bool> m_monitoring{this, "monitoring", false};
+ Gaudi::Property<bool> m_requireSingle{this, "requireSingle", false};
+ mutable std::vector<double> dx_vector{}, dy_vector{}, dz_vector{}, pullx_vector{}, pully_vector{}, pullz_vector{};
+
+ struct monitoringHistos {
+ mutable std::array<Gaudi::Accumulators::Histogram<1>, 5> m_histo_nTracksBins_dx;
+ mutable std::array<Gaudi::Accumulators::Histogram<1>, 5> m_histo_nTracksBins_dy;
+ mutable std::array<Gaudi::Accumulators::Histogram<1>, 5> m_histo_nTracksBins_dz;
+ mutable Gaudi::Accumulators::Histogram<1> m_histo_pullx_Monitoring;
+ mutable Gaudi::Accumulators::Histogram<1> m_histo_pully_Monitoring;
+ mutable Gaudi::Accumulators::Histogram<1> m_histo_pullz_Monitoring;
+
+
template <std::size_t... IDXs>
static std::array<Gaudi::Accumulators::Histogram<1>, sizeof...( IDXs )>
histo1DArrayBuilder( const VertexCompare* owner, const std::string& name, const std::string& title,
@@ -212,22 +290,28 @@ private:
title + std::to_string( IDXs ),
{std::get<0>( xbins ), std::get<1>( xbins ), std::get<2>( xbins )}}...}};
}
- monitoringHistos( const VertexCompare* owner )
- : m_histo_dx{histo1DArrayBuilder( owner, "dx", "dx, mm", {50, -0.25, 0.25}, std::make_index_sequence<5>() )}
- , m_histo_dy{histo1DArrayBuilder( owner, "dy", "dy, mm", {50, -0.25, 0.25}, std::make_index_sequence<5>() )}
- , m_histo_dz{histo1DArrayBuilder( owner, "dz", "dz, mm", {50, -1.5, 1.5}, std::make_index_sequence<5>() )}
- , m_histo_pullx{owner, "pullx", "pull x", {20, -2, 2}}
- , m_histo_pully{owner, "pully", "pull y", {20, -2, 2}}
- , m_histo_pullz{owner, "pulz", "pull z", {20, -2, 2}} {}
+ monitoringHistos( const VertexCompare* owner )
+ : m_histo_nTracksBins_dx{histo1DArrayBuilder( owner, "dx_Monitoring_", "dx, mm", {50, -0.15, 0.15}, std::make_index_sequence<5>() )}
+ , m_histo_nTracksBins_dy{histo1DArrayBuilder( owner, "dy_Monitoring_", "dy, mm", {50, -0.15, 0.15}, std::make_index_sequence<5>() )}
+ , m_histo_nTracksBins_dz{histo1DArrayBuilder( owner, "dz_Monitoring_", "dz, mm", {50, -1.5, 1.5}, std::make_index_sequence<5>() )}
+ , m_histo_pullx_Monitoring{owner, "pullx_Monitoring", "pull x", {20, -5, 5}}
+ , m_histo_pully_Monitoring{owner, "pully_Monitoring", "pull y", {20, -5, 5}}
+ , m_histo_pullz_Monitoring{owner, "pullz_Monitoring", "pull z", {20, -5, 5}} {}
};
std::unique_ptr<monitoringHistos> m_monitoringHistos;
- Gaudi::Property<bool> m_produceHistogram{this, "produceHistogram", true};
- Gaudi::Property<bool> m_produceNtuple{this, "produceNtuple", true};
- Gaudi::Property<bool> m_monitoring{this, "monitoring", false};
- Gaudi::Property<bool> m_requireSingle{this, "requireSingle", false};
-
+ using axis1D = Gaudi::Accumulators::Axis<Gaudi::Accumulators::Histogram<1>::AxisArithmeticType>;
+ mutable Gaudi::Accumulators::RootHistogram<1> m_histo_dx{this, "dx", "dx, mm", axis1D{50, -0.15, 0.15} };
+ mutable Gaudi::Accumulators::RootHistogram<1> m_histo_dy{ this, "dy", "dy, mm", axis1D{50, -0.15, 0.15} };
+ mutable Gaudi::Accumulators::RootHistogram<1> m_histo_dz{ this, "dz", "dz, mm", axis1D{50, -1.5, 1.5} };
+ mutable Gaudi::Accumulators::RootHistogram<1> m_histo_pullx{ this, "pullx", "pull x", axis1D{20, -5, 5} };
+ mutable Gaudi::Accumulators::RootHistogram<1> m_histo_pully{ this, "pully", "pull y", axis1D{20, -5, 5} };
+ mutable Gaudi::Accumulators::RootHistogram<1> m_histo_pullz{ this, "pullz", "pull z", axis1D{20, -5, 5} };
+ mutable Gaudi::Accumulators::Histogram<1> m_nTracks1{ this, "ntracks1", "Number of tracks in PV set 1", axis1D{50, 0., 150.} };
+ mutable Gaudi::Accumulators::Histogram<1> m_nTracks2{ this, "ntracks2", "Number of tracks in PV set 2", axis1D{ 50, 0., 150.}};
+ mutable Gaudi::Accumulators::Histogram<1> m_nTracks_dif{ this, "dtracks", "Difference in the number of tracks in two sets of PVs", axis1D{ 50, 0., 150.} };
+
mutable Gaudi::Accumulators::Counter<> m_nVtx{this, "Number of pairs of vertices in processed events"};
mutable Gaudi::Accumulators::Counter<> m_nRec1{this, "Number of vertices in input1"};
mutable Gaudi::Accumulators::Counter<> m_nRec2{this, "Number of vertices in input2"};
@@ -267,20 +351,21 @@ StatusCode VertexCompare::initialize() {
//=============================================================================
void VertexCompare::operator()( Vertices const& recoVtx1, Vertices const& recoVtx2,
LHCb::Conditions::InteractionRegion const& region ) const {
- if ( msgLevel( MSG::DEBUG ) ) debug() << "==> Execute" << endmsg;
+ if ( msgLevel( MSG::DEBUG ) ) debug() << "==> Execute" << endmsg;
const auto beamspot = region.avgPosition;
std::vector<RecPVInfo<VertexType>> fullVrt1, fullVrt2;
fullVrt1.reserve( recoVtx1.size() );
fullVrt2.reserve( recoVtx2.size() );
+ if(m_requireSingle == true && !(recoVtx1.size() == 1 && recoVtx2.size() == 1)) return;
+
if ( debugLevel() ) debug() << " Vtx Properities x y z chi2/ndof ntracks" << endmsg;
int size_diff = -9999;
RecPVInfo<VertexType> recinfo1;
RecPVInfo<VertexType> recinfo2;
- if ( m_requireSingle == false ) {
for ( auto const& pv : recoVtx1 ) {
++m_nRec1;
recinfo1.pRECPV = &pv;
@@ -324,59 +409,7 @@ void VertexCompare::operator()( Vertices const& recoVtx1, Vertices const& recoVt
if ( debugLevel() ) debug() << "fullVrt1 size " << fullVrt1.size() << endmsg;
if ( debugLevel() ) debug() << "fullVrt2 size " << fullVrt2.size() << endmsg;
size_diff = fullVrt1.size() - fullVrt2.size();
- } else if ( m_requireSingle == true && recoVtx1.size() == 1 && recoVtx2.size() == 1 ) {
- for ( auto const& pv : recoVtx1 ) {
- RecPVInfo<VertexType> recinfo1;
- ++m_nRec1;
- recinfo1.pRECPV = &pv;
- recinfo1.position = pv.position();
- recinfo1.covPV = pv.covMatrix();
- recinfo1.positionSigma = Gaudi::XYZPoint{std::sqrt( recinfo1.covPV( 0, 0 ) ), std::sqrt( recinfo1.covPV( 1, 1 ) ),
- std::sqrt( recinfo1.covPV( 2, 2 ) )};
- recinfo1.ntracks = pv.nTracks();
- if ( recinfo1.ntracks >= 10 ) { ++m_nRec_geq10tr_1; }
- recinfo1.chi2 = pv.chi2();
- recinfo1.nDoF = pv.nDoF();
-
- if ( debugLevel() )
- debug() << " " << pv.position().x() << " " << pv.position().y() << " " << pv.position().z()
- << " " << pv.chi2PerDoF() << " " << pv.nTracks() << endmsg;
- fullVrt1.push_back( recinfo1 );
- } // end of loop over vertices1
-
- for ( auto const& pv : recoVtx2 ) {
- RecPVInfo<VertexType> recinfo2;
- ++m_nRec2;
- recinfo2.pRECPV = &pv;
- recinfo2.position = pv.position();
- recinfo2.covPV = pv.covMatrix();
- recinfo2.positionSigma = Gaudi::XYZPoint{std::sqrt( recinfo2.covPV( 0, 0 ) ), std::sqrt( recinfo2.covPV( 1, 1 ) ),
- std::sqrt( recinfo2.covPV( 2, 2 ) )};
- recinfo2.ntracks = pv.nTracks();
- if ( recinfo2.ntracks >= 10 ) { ++m_nRec_geq10tr_2; }
- recinfo2.chi2 = pv.chi2();
- recinfo2.nDoF = pv.nDoF();
-
- if ( debugLevel() )
- debug() << " " << pv.position().x() << " " << pv.position().y() << " " << pv.position().z()
- << " " << pv.chi2PerDoF() << " " << pv.nTracks() << endmsg;
- fullVrt2.push_back( recinfo2 );
- } // end of loop over vertices2
-
- // sorting is unnecessary for singular vertices
-
- if ( debugLevel() )
- debug() << "fullVrt1 size "
- << "1" << endmsg;
- if ( debugLevel() )
- debug() << "fullVrt2 size "
- << "1" << endmsg;
- size_diff = 0;
- }
-
- // this check takes care of events which evade the selection
- if ( !( m_requireSingle == true && ( recoVtx1.size() != 1 || recoVtx2.size() != 1 ) ) ) {
- if ( m_produceHistogram.value() ) { plot( double( size_diff ), 1001, "size_diff", -5.5, 5.5, 11 ); }
+
if ( m_produceNtuple.value() ) {
Tuple myTuple_evt = nTuple( "PV_nTuple_evt", "PV_nTuple_evt", CLID_ColumnWiseTuple );
myTuple_evt->column( "size_diff", double( size_diff ) ).ignore();
@@ -384,7 +417,6 @@ void VertexCompare::operator()( Vertices const& recoVtx1, Vertices const& recoVt
myTuple_evt->column( "size_2", double( fullVrt2.size() ) ).ignore();
myTuple_evt->write().ignore();
}
- }
std::vector<int> link;
int ntracks1 = 0;
@@ -408,7 +440,7 @@ void VertexCompare::operator()( Vertices const& recoVtx1, Vertices const& recoVt
double dy = -99999.;
double dz = -99999.;
int oIt = 0;
- int multiplicity = -99999;
+ int pv_rank = -99999;
if ( fullVrt1.size() != 0 && fullVrt2.size() != 0 ) { matchByDistance( fullVrt1, fullVrt2, link ); }
@@ -457,71 +489,44 @@ void VertexCompare::operator()( Vertices const& recoVtx1, Vertices const& recoVt
// no pass is ntracks == 1
// we have reco cut on ntracks >=4
if ( ntracks1 > 1 ) {
- multiplicity = oIt + 1;
+ pv_rank = oIt + 1;
} else if ( ntracks2 > 1 ) {
- multiplicity = link.at( oIt ) + 1;
+ pv_rank = link.at( oIt ) + 1;
}
if ( m_monitoring.value() ) {
- //if ( dz < 2 && size_diff < 3) {
+ if ( dz < 2 && size_diff < 3 && pv_rank < 2 && ( ntracks1 + ntracks2 ) / 2 > ntrack_bins[0]) {
int binCount = 0;
- for ( size_t i = 0; i < ntrack_bins.size(); ++i ) {
+ for ( size_t i = 1; i < ntrack_bins.size(); ++i ) {
if ( std::lround(( ntracks1 + ntracks2 ) / 2) <= ntrack_bins[i] ) { break; }
binCount++;
}
- auto& histos = *m_monitoringHistos.get();
- ++histos.m_histo_dx[binCount][dx];
- ++histos.m_histo_dy[binCount][dy];
- ++histos.m_histo_dz[binCount][dz];
- ++histos.m_histo_pullx[dx / errx];
- ++histos.m_histo_pully[dy / erry];
- ++histos.m_histo_pullz[dz / errz];
- // } if end
- }
-
- if ( m_produceHistogram.value() ) {
- plot( dx, 1021, "dx", -0.25, 0.25, 50 );
- plot( dy, 1022, "dy", -0.25, 0.25, 50 );
- plot( dz, 1023, "dz", -1.5, 1.5, 60 );
- plot( x1, 2021, "x1", -0.25, 0.25, 50 );
- plot( y1, 2022, "y1", -0.25, 0.25, 50 );
- plot( z1, 2023, "z1", -100., 100., 50 );
- plot( x2, 3021, "x2", -0.25, 0.25, 50 );
- plot( y2, 3022, "y2", -0.25, 0.25, 50 );
- plot( z2, 3023, "z2", -100., 100., 50 );
- plot( std::sqrt( sigx_part1 ), 4011, "x err 1", 0., .1, 50 );
- plot( std::sqrt( sigy_part1 ), 4012, "y err 1", 0., .1, 50 );
- plot( std::sqrt( sigz_part1 ), 4013, "z err 1", 0., .5, 50 );
- plot( std::sqrt( sigx_part2 ), 4021, "x err 2", 0., .1, 50 );
- plot( std::sqrt( sigy_part2 ), 4022, "y err 2", 0., .1, 50 );
- plot( std::sqrt( sigz_part2 ), 4023, "z err 2", 0., .5, 50 );
- plot( dx / errx, 1031, "pullx", -2., 2., 20 );
- plot( dy / erry, 1032, "pully", -2., 2., 20 );
- plot( dz / errz, 1033, "pullz", -2., 2., 20 );
- plot( double( ntracks1 ), 1041, "ntracks1", 0., 150., 50 );
- plot( double( ntracks2 ), 1042, "ntracks2", 0., 150., 50 );
- plot( double( dtracks ), 1043, "dtracks", 0., 150., 50 );
- if ( recoVtx1.size() > 1 ) {
- for ( unsigned int i1 = 0; i1 < recoVtx1.size(); i1++ ) {
- for ( unsigned int i2 = 0; i2 < recoVtx1.size(); i2++ ) {
- if ( i2 != i1 ) {
- double vdz = recoVtx1[i1].position().z() - recoVtx1[i2].position().z();
- plot( vdz, 1201, "dz vertices 1", -150., 150., 100 );
- }
- }
- }
- }
- if ( recoVtx2.size() > 1 ) {
- for ( unsigned int i1 = 0; i1 < recoVtx2.size(); i1++ ) {
- for ( unsigned int i2 = 0; i2 < recoVtx2.size(); i2++ ) {
- if ( i2 != i1 ) {
- double vdz = recoVtx2[i1].position().z() - recoVtx2[i2].position().z();
- plot( vdz, 1202, "dz vertices 2", -150., 150., 100 );
- }
- }
- }
- }
+ auto& monitoringHistos = *m_monitoringHistos.get();
+ ++monitoringHistos.m_histo_nTracksBins_dx[binCount][dx];
+ ++monitoringHistos.m_histo_nTracksBins_dy[binCount][dy];
+ ++monitoringHistos.m_histo_nTracksBins_dz[binCount][dz];
+ ++monitoringHistos.m_histo_pullx_Monitoring[dx / errx];
+ ++monitoringHistos.m_histo_pully_Monitoring[dy / erry];
+ ++monitoringHistos.m_histo_pullz_Monitoring[dz / errz];
+ }
}
+ if (m_produceHistogram.value()) {
+ dx_vector.push_back(dx);
+ dy_vector.push_back(dy);
+ dz_vector.push_back(dz);
+ pullx_vector.push_back(dx / errx);
+ pully_vector.push_back(dy / erry);
+ pullz_vector.push_back(dz / errz);
+ ++m_histo_dx[dx];
+ ++m_histo_dy[dy];
+ ++m_histo_dz[dz];
+ ++m_histo_pullx[dx / errx];
+ ++m_histo_pully[dy / erry];
+ ++m_histo_pullz[dz / errz];
+ ++m_nTracks1[ntracks1];
+ ++m_nTracks2[ntracks2];
+ ++m_nTracks_dif[ntracks2-ntracks1];
+ } //else end
if ( m_produceNtuple.value() ) {
Tuple myTuple = nTuple( "PV_nTuple", "PV_nTuple", CLID_ColumnWiseTuple );
@@ -584,7 +589,7 @@ void VertexCompare::operator()( Vertices const& recoVtx1, Vertices const& recoVt
myTuple->column( "match", set1.match ).ignore();
myTuple->column( "isopposite", set1.opposite_container ).ignore();
myTuple->column( "evt", int( m_nEvt.nEntries() ) ).ignore();
- myTuple->column( "multiplicity", int( multiplicity ) ).ignore();
+ myTuple->column( "pv_rank", int( pv_rank ) ).ignore();
myTuple->write().ignore();
}
oIt++;
@@ -631,7 +636,7 @@ void VertexCompare::operator()( Vertices const& recoVtx1, Vertices const& recoVt
myTuple->column( "equal_sizes", set1.equal_sizes ).ignore();
myTuple->column( "count", int( counter1 ) ).ignore();
myTuple->column( "match", set1.match ).ignore();
- myTuple->column( "multiplicity", set1.multiplicity ).ignore();
+ myTuple->column( "pv_rank", set1.pv_rank ).ignore();
myTuple->column( "isopposite", set1.opposite_container ).ignore();
myTuple->column( "evt", int( m_nEvt.nEntries() ) ).ignore();
myTuple->write().ignore();
@@ -679,7 +684,7 @@ void VertexCompare::operator()( Vertices const& recoVtx1, Vertices const& recoVt
myTuple->column( "equal_sizes", set2.equal_sizes ).ignore();
myTuple->column( "count", int( counter2 ) ).ignore();
myTuple->column( "match", set2.match ).ignore();
- myTuple->column( "multiplicity", set2.multiplicity ).ignore();
+ myTuple->column( "pv_rank", set2.pv_rank ).ignore();
myTuple->column( "isopposite", set2.opposite_container ).ignore();
myTuple->column( "evt", int( m_nEvt.nEntries() ) ).ignore();
myTuple->write().ignore();
@@ -687,7 +692,7 @@ void VertexCompare::operator()( Vertices const& recoVtx1, Vertices const& recoVt
oIt++;
m_nEvt += 1;
}
-}
+ }
//=============================================================================
// Finalize
@@ -724,80 +729,85 @@ StatusCode VertexCompare::finalize() {
<< ( (double)m_nVtx.nEntries() - (double)m_nVtx_geq10tr_matched2.nEntries() ) /
( (double)m_nRec2.nEntries() - (double)m_nRec_geq10tr_2.nEntries() ) * 100
<< " % " << endmsg;
+
+ if(m_produceHistogram.value()){
- const AIDA::IHistogram1D* dx = histo( HistoID( 1021 ) );
- const AIDA::IHistogram1D* pullx = histo( HistoID( 1031 ) );
- const AIDA::IHistogram1D* dy = histo( HistoID( 1022 ) );
- const AIDA::IHistogram1D* pully = histo( HistoID( 1032 ) );
- const AIDA::IHistogram1D* dz = histo( HistoID( 1023 ) );
- const AIDA::IHistogram1D* pullz = histo( HistoID( 1033 ) );
-
- if ( dx ) {
- info() << " ---------------------------------------" << endmsg;
- info() << "dx: "
- << format( "mean = %5.3f +/- %5.3f, RMS = %5.3f +/- %5.3f", dx->mean(),
- Gaudi::Utils::HistoStats::meanErr( dx ), dx->rms(), Gaudi::Utils::HistoStats::rmsErr( dx ) )
- << endmsg;
+ for (double value : dx_vector) {
+ ++m_histo_dx[value];
+ }
+ for (double value : dy_vector) {
+ ++m_histo_dy[value];
+ }
+ for (double value : dz_vector) {
+ ++m_histo_dz[value];
+ }
+ for (double value : pullx_vector) {
+ ++m_histo_pullx[value];
+ }
+ for (double value : pully_vector) {
+ ++m_histo_pully[value];
+ }
+ for (double value : pullz_vector) {
+ ++m_histo_pullz[value];
}
- if ( dy ) {
- info() << "dy: "
- << format( "mean = %5.3f +/- %5.3f, RMS = %5.3f +/- %5.3f", dy->mean(),
- Gaudi::Utils::HistoStats::meanErr( dy ), dy->rms(), Gaudi::Utils::HistoStats::rmsErr( dy ) )
+ info() << "dx: "
+ << format( "mean = %5.3f +/- %5.3f, RMS = %5.3f +/- %5.3f", m_histo_dx.mean(),
+ meanErr( m_histo_dx ), m_histo_dx.standard_deviation(), rmsErr( m_histo_dx) )
<< endmsg;
- }
- if ( dz ) {
- info() << "dz: "
- << format( "mean = %5.3f +/- %5.3f, RMS = %5.3f +/- %5.3f", dz->mean(),
- Gaudi::Utils::HistoStats::meanErr( dz ), dz->rms(), Gaudi::Utils::HistoStats::rmsErr( dz ) )
+
+
+ info() << "dy: "
+ << format( "mean = %5.3f +/- %5.3f, RMS = %5.3f +/- %5.3f", m_histo_dy.mean(),
+ meanErr( m_histo_dy ), m_histo_dy.standard_deviation(), rmsErr( m_histo_dy) )
<< endmsg;
- }
+
+
+ info() << "dz: "
+ << format( "mean = %5.3f +/- %5.3f, RMS = %5.3f +/- %5.3f", m_histo_dz.mean(),
+ meanErr( m_histo_dz ), m_histo_dz.standard_deviation(), rmsErr( m_histo_dz) )
+ << endmsg;
+
info() << " ---------------------------------------" << endmsg;
- if ( pullx ) {
- info() << "pullx: "
- << format( "mean = %5.3f +/- %5.3f, RMS = %5.3f +/- %5.3f", pullx->mean(),
- Gaudi::Utils::HistoStats::meanErr( pullx ), pullx->rms(),
- Gaudi::Utils::HistoStats::rmsErr( pullx ) )
+
+ info() << "pullx: "
+ << format( "mean = %5.3f +/- %5.3f, RMS = %5.3f +/- %5.3f", m_histo_pullx.mean(),
+ meanErr( m_histo_pullx ), m_histo_pullx.standard_deviation(), rmsErr( m_histo_pullx ) )
<< endmsg;
- }
- if ( pully ) {
- info() << "pully: "
- << format( "mean = %5.3f +/- %5.3f, RMS = %5.3f +/- %5.3f", pully->mean(),
- Gaudi::Utils::HistoStats::meanErr( pully ), pully->rms(),
- Gaudi::Utils::HistoStats::rmsErr( pully ) )
+
+ info() << "pully: "
+ << format( "mean = %5.3f +/- %5.3f, RMS = %5.3f +/- %5.3f", m_histo_pully.mean(),
+ meanErr( m_histo_pully ), m_histo_pully.standard_deviation(), rmsErr( m_histo_pully) )
<< endmsg;
- }
- if ( pullz ) {
- info() << "pullz: "
- << format( "mean = %5.3f +/- %5.3f, RMS = %5.3f +/- %5.3f", pullz->mean(),
- Gaudi::Utils::HistoStats::meanErr( pullz ), pullz->rms(),
- Gaudi::Utils::HistoStats::rmsErr( pullz ) )
+
+ info() << "pullz: "
+ << format( "mean = %5.3f +/- %5.3f, RMS = %5.3f +/- %5.3f", m_histo_pullz.mean(),
+ meanErr( m_histo_pullz ), m_histo_pullz.standard_deviation(), rmsErr( m_histo_pullz) )
<< endmsg;
- }
+
info() << " ---------------------------------------" << endmsg;
- if ( pullx ) {
- info() << "diff in x: "
- << format( "%5.3f +/- %5.3f", std::sqrt( 1.0 + pullx->rms() ) - 1.0,
- Gaudi::Utils::HistoStats::rmsErr( pullx ) * 0.5 / std::sqrt( 1.0 + pullx->rms() ) )
+
+ info() << "diff in x: "
+ << format( "%5.3f +/- %5.3f", std::sqrt( 1.0 + m_histo_pullx.standard_deviation() ) - 1.0,
+ rmsErr( m_histo_pullx ) * 0.5 / std::sqrt( 1.0 + m_histo_pullx.standard_deviation() ) )
<< endmsg;
- }
- if ( pully ) {
- info() << "diff in y: "
- << format( "%5.3f +/- %5.3f", std::sqrt( 1.0 + pully->rms() ) - 1.0,
- Gaudi::Utils::HistoStats::rmsErr( pully ) * 0.5 / std::sqrt( 1.0 + pully->rms() ) )
+
+ info() << "diff in y: "
+ << format( "%5.3f +/- %5.3f", std::sqrt( 1.0 + m_histo_pully.standard_deviation() ) - 1.0,
+ rmsErr( m_histo_pully ) * 0.5 / std::sqrt( 1.0 + m_histo_pully.standard_deviation() ) )
<< endmsg;
- }
- if ( pullz ) {
- info() << "diff in z: "
- << format( "%5.3f +/- %5.3f", std::sqrt( 1.0 + pullz->rms() ) - 1.0,
- Gaudi::Utils::HistoStats::rmsErr( pullz ) * 0.5 / std::sqrt( 1.0 + pullz->rms() ) )
+
+ info() << "diff in z: "
+ << format( "%5.3f +/- %5.3f", std::sqrt( 1.0 + m_histo_pullz.standard_deviation() ) - 1.0,
+ rmsErr( m_histo_pullz ) * 0.5 / std::sqrt( 1.0 + m_histo_pullz.standard_deviation() ) )
<< endmsg;
- }
- info() << " ============================================" << endmsg;
+ info() << " ============================================" << endmsg;
+ }
return Consumer::finalize(); // Must be called after all other actions
}
+
//=============================================================================
// Match vertices by distance
//=============================================================================
--
GitLab
From d8b30a4c834a0a4918683137acdfb1570dc5fe72 Mon Sep 17 00:00:00 2001
From: Gitlab CI <noreply@cern.ch>
Date: Fri, 24 May 2024 14:15:47 +0000
Subject: [PATCH 10/21] Fixed formatting
patch generated by https://gitlab.cern.ch/lhcb/Rec/-/jobs/39274748
---
Tr/TrackMonitors/src/VertexCompare.cpp | 409 ++++++++++++-------------
1 file changed, 199 insertions(+), 210 deletions(-)
diff --git a/Tr/TrackMonitors/src/VertexCompare.cpp b/Tr/TrackMonitors/src/VertexCompare.cpp
index 1e42189e092..b8294a5d5dc 100644
--- a/Tr/TrackMonitors/src/VertexCompare.cpp
+++ b/Tr/TrackMonitors/src/VertexCompare.cpp
@@ -89,31 +89,31 @@ bool trackcomp( const RecPVInfo<VertexType>& first, const RecPVInfo<VertexType>&
}
struct VtxVariables {
- int ntracks = -9999;
- double x = -9999.;
- double y = -9999.;
- double z = -9999.;
- double dxr = -9999.;
- double dyr = -9999.;
- double r = -9999.;
- double errx = -9999.;
- double erry = -9999.;
- double errz = -9999.;
- double err_r = -9999.;
- double covxx = -9999.;
- double covyy = -9999.;
- double covzz = -9999.;
- double covxy = -9999.;
- double covxz = -9999.;
- double covyz = -9999.;
- double chi2 = -9999.;
- double nDoF = -9999.;
- int dsize = -9999;
- int match = -9999;
- int pv_rank = -9999;
- bool equal_sizes = false;
- bool single = false;
- bool opposite_container = false;
+ int ntracks = -9999;
+ double x = -9999.;
+ double y = -9999.;
+ double z = -9999.;
+ double dxr = -9999.;
+ double dyr = -9999.;
+ double r = -9999.;
+ double errx = -9999.;
+ double erry = -9999.;
+ double errz = -9999.;
+ double err_r = -9999.;
+ double covxx = -9999.;
+ double covyy = -9999.;
+ double covzz = -9999.;
+ double covxy = -9999.;
+ double covxz = -9999.;
+ double covyz = -9999.;
+ double chi2 = -9999.;
+ double nDoF = -9999.;
+ int dsize = -9999;
+ int match = -9999;
+ int pv_rank = -9999;
+ bool equal_sizes = false;
+ bool single = false;
+ bool opposite_container = false;
std::vector<double> dx_vector{}, dy_vector{}, dz_vector{}, pullx_vector{}, pully_vector{}, pullz_vector{};
};
@@ -177,7 +177,7 @@ VtxVariables SetVtxVariables( const RecPVInfo<VertexType>& vrtf, const int& coun
// ============================================================================
template <typename Histogram>
-double meanErr( Histogram const & h ) {
+double meanErr( Histogram const& h ) {
const auto n = h.nEntries();
return ( 0 >= n ? 0.0 : h.standard_deviation() / std::sqrt( n ) );
}
@@ -186,28 +186,28 @@ double meanErr( Histogram const & h ) {
// Fourth momentum calculation
// ============================================================================
template <typename Histogram>
-double fourthMomentum( Histogram const & h ) {
+double fourthMomentum( Histogram const& h ) {
const auto n = h.nEntries();
if ( 0 >= n ) { return 0.0; } // RETURN
// get the mean
const auto h_mean = h.mean();
// number of bins
- const auto& axis = h.axis();
- const int nBins = h.totNBins();
- double result{ 0 }, weight{ 0 };
+ const auto& axis = h.axis();
+ const int nBins = h.totNBins();
+ double result{0}, weight{0};
// loop over all bins
- double minVal = axis[0].minValue;
- double maxVal = axis[0].maxValue;
- double binWidth = (maxVal-minVal)/nBins;
+ double minVal = axis[0].minValue;
+ double maxVal = axis[0].maxValue;
+ double binWidth = ( maxVal - minVal ) / nBins;
std::vector<double> binCenters;
- for (int i = 0; i < nBins; ++i) {
- double center = minVal + (i + 0.5) * binWidth;
- binCenters.push_back(center);
- }
- for (int i = 0; i < nBins; ++i) {
- const auto yBin = h.binValue( i ); // bin weight
+ for ( int i = 0; i < nBins; ++i ) {
+ double center = minVal + ( i + 0.5 ) * binWidth;
+ binCenters.push_back( center );
+ }
+ for ( int i = 0; i < nBins; ++i ) {
+ const auto yBin = h.binValue( i ); // bin weight
weight += yBin;
result += yBin * std::pow( binCenters[i] - h_mean, 4 );
}
@@ -221,26 +221,25 @@ double fourthMomentum( Histogram const & h ) {
// ============================================================================
template <typename Histogram>
-double kurtosis( Histogram const & h ) {
+double kurtosis( Histogram const& h ) {
const auto mu4 = fourthMomentum( h );
- const auto s4 = std::pow( h.standard_deviation() , 4 );
+ const auto s4 = std::pow( h.standard_deviation(), 4 );
return ( std::fabs( s4 ) > 0 ? mu4 / s4 - 3.0 : 0.0 );
}
// ============================================================================
// get an error in the rms value
// ============================================================================
-
- template <typename Histogram>
- double rmsErr( Histogram const & h ) {
- const auto n = h.nEntries();
- if ( 1 >= n ) { return 0.0; }
- auto result = 2.0 + kurtosis( h );
- result += 2.0 / ( n - 1 );
- result /= 4.0 * n;
- return h.standard_deviation() * std::sqrt( std::max( result, 0.0 ) );
- }
+template <typename Histogram>
+double rmsErr( Histogram const& h ) {
+ const auto n = h.nEntries();
+ if ( 1 >= n ) { return 0.0; }
+ auto result = 2.0 + kurtosis( h );
+ result += 2.0 / ( n - 1 );
+ result /= 4.0 * n;
+ return h.standard_deviation() * std::sqrt( std::max( result, 0.0 ) );
+}
class VertexCompare : public LHCb::Algorithm::Consumer<
void( Vertices const&, Vertices const&, LHCb::Conditions::InteractionRegion const& ),
@@ -263,7 +262,7 @@ public:
private:
bool debugLevel() const { return msgLevel( MSG::DEBUG ) || msgLevel( MSG::VERBOSE ); }
- static constexpr auto ntrack_bins = std::array{ 3.5 , 16.5 , 29.5 , 43.5 , 56.5 , 69.5 };
+ static constexpr auto ntrack_bins = std::array{3.5, 16.5, 29.5, 43.5, 56.5, 69.5};
Gaudi::Property<bool> m_produceHistogram{this, "produceHistogram", true};
Gaudi::Property<bool> m_produceNtuple{this, "produceNtuple", true};
@@ -271,7 +270,7 @@ private:
Gaudi::Property<bool> m_requireSingle{this, "requireSingle", false};
mutable std::vector<double> dx_vector{}, dy_vector{}, dz_vector{}, pullx_vector{}, pully_vector{}, pullz_vector{};
-
+
struct monitoringHistos {
mutable std::array<Gaudi::Accumulators::Histogram<1>, 5> m_histo_nTracksBins_dx;
mutable std::array<Gaudi::Accumulators::Histogram<1>, 5> m_histo_nTracksBins_dy;
@@ -280,7 +279,6 @@ private:
mutable Gaudi::Accumulators::Histogram<1> m_histo_pully_Monitoring;
mutable Gaudi::Accumulators::Histogram<1> m_histo_pullz_Monitoring;
-
template <std::size_t... IDXs>
static std::array<Gaudi::Accumulators::Histogram<1>, sizeof...( IDXs )>
histo1DArrayBuilder( const VertexCompare* owner, const std::string& name, const std::string& title,
@@ -290,10 +288,13 @@ private:
title + std::to_string( IDXs ),
{std::get<0>( xbins ), std::get<1>( xbins ), std::get<2>( xbins )}}...}};
}
- monitoringHistos( const VertexCompare* owner )
- : m_histo_nTracksBins_dx{histo1DArrayBuilder( owner, "dx_Monitoring_", "dx, mm", {50, -0.15, 0.15}, std::make_index_sequence<5>() )}
- , m_histo_nTracksBins_dy{histo1DArrayBuilder( owner, "dy_Monitoring_", "dy, mm", {50, -0.15, 0.15}, std::make_index_sequence<5>() )}
- , m_histo_nTracksBins_dz{histo1DArrayBuilder( owner, "dz_Monitoring_", "dz, mm", {50, -1.5, 1.5}, std::make_index_sequence<5>() )}
+ monitoringHistos( const VertexCompare* owner )
+ : m_histo_nTracksBins_dx{histo1DArrayBuilder( owner, "dx_Monitoring_", "dx, mm", {50, -0.15, 0.15},
+ std::make_index_sequence<5>() )}
+ , m_histo_nTracksBins_dy{histo1DArrayBuilder( owner, "dy_Monitoring_", "dy, mm", {50, -0.15, 0.15},
+ std::make_index_sequence<5>() )}
+ , m_histo_nTracksBins_dz{histo1DArrayBuilder( owner, "dz_Monitoring_", "dz, mm", {50, -1.5, 1.5},
+ std::make_index_sequence<5>() )}
, m_histo_pullx_Monitoring{owner, "pullx_Monitoring", "pull x", {20, -5, 5}}
, m_histo_pully_Monitoring{owner, "pully_Monitoring", "pull y", {20, -5, 5}}
, m_histo_pullz_Monitoring{owner, "pullz_Monitoring", "pull z", {20, -5, 5}} {}
@@ -302,16 +303,19 @@ private:
std::unique_ptr<monitoringHistos> m_monitoringHistos;
using axis1D = Gaudi::Accumulators::Axis<Gaudi::Accumulators::Histogram<1>::AxisArithmeticType>;
- mutable Gaudi::Accumulators::RootHistogram<1> m_histo_dx{this, "dx", "dx, mm", axis1D{50, -0.15, 0.15} };
- mutable Gaudi::Accumulators::RootHistogram<1> m_histo_dy{ this, "dy", "dy, mm", axis1D{50, -0.15, 0.15} };
- mutable Gaudi::Accumulators::RootHistogram<1> m_histo_dz{ this, "dz", "dz, mm", axis1D{50, -1.5, 1.5} };
- mutable Gaudi::Accumulators::RootHistogram<1> m_histo_pullx{ this, "pullx", "pull x", axis1D{20, -5, 5} };
- mutable Gaudi::Accumulators::RootHistogram<1> m_histo_pully{ this, "pully", "pull y", axis1D{20, -5, 5} };
- mutable Gaudi::Accumulators::RootHistogram<1> m_histo_pullz{ this, "pullz", "pull z", axis1D{20, -5, 5} };
- mutable Gaudi::Accumulators::Histogram<1> m_nTracks1{ this, "ntracks1", "Number of tracks in PV set 1", axis1D{50, 0., 150.} };
- mutable Gaudi::Accumulators::Histogram<1> m_nTracks2{ this, "ntracks2", "Number of tracks in PV set 2", axis1D{ 50, 0., 150.}};
- mutable Gaudi::Accumulators::Histogram<1> m_nTracks_dif{ this, "dtracks", "Difference in the number of tracks in two sets of PVs", axis1D{ 50, 0., 150.} };
-
+ mutable Gaudi::Accumulators::RootHistogram<1> m_histo_dx{this, "dx", "dx, mm", axis1D{50, -0.15, 0.15}};
+ mutable Gaudi::Accumulators::RootHistogram<1> m_histo_dy{this, "dy", "dy, mm", axis1D{50, -0.15, 0.15}};
+ mutable Gaudi::Accumulators::RootHistogram<1> m_histo_dz{this, "dz", "dz, mm", axis1D{50, -1.5, 1.5}};
+ mutable Gaudi::Accumulators::RootHistogram<1> m_histo_pullx{this, "pullx", "pull x", axis1D{20, -5, 5}};
+ mutable Gaudi::Accumulators::RootHistogram<1> m_histo_pully{this, "pully", "pull y", axis1D{20, -5, 5}};
+ mutable Gaudi::Accumulators::RootHistogram<1> m_histo_pullz{this, "pullz", "pull z", axis1D{20, -5, 5}};
+ mutable Gaudi::Accumulators::Histogram<1> m_nTracks1{this, "ntracks1", "Number of tracks in PV set 1",
+ axis1D{50, 0., 150.}};
+ mutable Gaudi::Accumulators::Histogram<1> m_nTracks2{this, "ntracks2", "Number of tracks in PV set 2",
+ axis1D{50, 0., 150.}};
+ mutable Gaudi::Accumulators::Histogram<1> m_nTracks_dif{
+ this, "dtracks", "Difference in the number of tracks in two sets of PVs", axis1D{50, 0., 150.}};
+
mutable Gaudi::Accumulators::Counter<> m_nVtx{this, "Number of pairs of vertices in processed events"};
mutable Gaudi::Accumulators::Counter<> m_nRec1{this, "Number of vertices in input1"};
mutable Gaudi::Accumulators::Counter<> m_nRec2{this, "Number of vertices in input2"};
@@ -351,96 +355,96 @@ StatusCode VertexCompare::initialize() {
//=============================================================================
void VertexCompare::operator()( Vertices const& recoVtx1, Vertices const& recoVtx2,
LHCb::Conditions::InteractionRegion const& region ) const {
- if ( msgLevel( MSG::DEBUG ) ) debug() << "==> Execute" << endmsg;
+ if ( msgLevel( MSG::DEBUG ) ) debug() << "==> Execute" << endmsg;
const auto beamspot = region.avgPosition;
std::vector<RecPVInfo<VertexType>> fullVrt1, fullVrt2;
fullVrt1.reserve( recoVtx1.size() );
fullVrt2.reserve( recoVtx2.size() );
- if(m_requireSingle == true && !(recoVtx1.size() == 1 && recoVtx2.size() == 1)) return;
-
+ if ( m_requireSingle == true && !( recoVtx1.size() == 1 && recoVtx2.size() == 1 ) ) return;
+
if ( debugLevel() ) debug() << " Vtx Properities x y z chi2/ndof ntracks" << endmsg;
- int size_diff = -9999;
+ int size_diff = -9999;
RecPVInfo<VertexType> recinfo1;
RecPVInfo<VertexType> recinfo2;
- for ( auto const& pv : recoVtx1 ) {
- ++m_nRec1;
- recinfo1.pRECPV = &pv;
- recinfo1.position = pv.position();
- recinfo1.covPV = pv.covMatrix();
- recinfo1.positionSigma = Gaudi::XYZPoint{std::sqrt( recinfo1.covPV( 0, 0 ) ), std::sqrt( recinfo1.covPV( 1, 1 ) ),
- std::sqrt( recinfo1.covPV( 2, 2 ) )};
- recinfo1.ntracks = pv.nTracks();
- if ( recinfo1.ntracks >= 10 ) { ++m_nRec_geq10tr_1; }
- recinfo1.chi2 = pv.chi2();
- recinfo1.nDoF = pv.nDoF();
-
- if ( debugLevel() )
- debug() << " " << pv.position().x() << " " << pv.position().y() << " " << pv.position().z()
- << " " << pv.chi2PerDoF() << " " << pv.nTracks() << endmsg;
- fullVrt1.push_back( recinfo1 );
- } // end of loop over vertices1
-
- for ( auto const& pv : recoVtx2 ) {
- ++m_nRec2;
- recinfo2.pRECPV = &pv;
- recinfo2.position = pv.position();
- recinfo2.covPV = pv.covMatrix();
- recinfo2.positionSigma = Gaudi::XYZPoint{std::sqrt( recinfo2.covPV( 0, 0 ) ), std::sqrt( recinfo2.covPV( 1, 1 ) ),
- std::sqrt( recinfo2.covPV( 2, 2 ) )};
- recinfo2.ntracks = pv.nTracks();
- if ( recinfo2.ntracks >= 10 ) { ++m_nRec_geq10tr_2; }
- recinfo2.chi2 = pv.chi2();
- recinfo2.nDoF = pv.nDoF();
-
- if ( debugLevel() )
- debug() << " " << pv.position().x() << " " << pv.position().y() << " " << pv.position().z()
- << " " << pv.chi2PerDoF() << " " << pv.nTracks() << endmsg;
- fullVrt2.push_back( recinfo2 );
- } // end of loop over vertices2
-
- // added sorting to mirror the 'multirec' variable implementation in PVChecker
- std::sort( fullVrt1.begin(), fullVrt1.end(), trackcomp );
- std::sort( fullVrt2.begin(), fullVrt2.end(), trackcomp );
-
- if ( debugLevel() ) debug() << "fullVrt1 size " << fullVrt1.size() << endmsg;
- if ( debugLevel() ) debug() << "fullVrt2 size " << fullVrt2.size() << endmsg;
- size_diff = fullVrt1.size() - fullVrt2.size();
-
- if ( m_produceNtuple.value() ) {
- Tuple myTuple_evt = nTuple( "PV_nTuple_evt", "PV_nTuple_evt", CLID_ColumnWiseTuple );
- myTuple_evt->column( "size_diff", double( size_diff ) ).ignore();
- myTuple_evt->column( "size_1", double( fullVrt1.size() ) ).ignore();
- myTuple_evt->column( "size_2", double( fullVrt2.size() ) ).ignore();
- myTuple_evt->write().ignore();
- }
+ for ( auto const& pv : recoVtx1 ) {
+ ++m_nRec1;
+ recinfo1.pRECPV = &pv;
+ recinfo1.position = pv.position();
+ recinfo1.covPV = pv.covMatrix();
+ recinfo1.positionSigma = Gaudi::XYZPoint{std::sqrt( recinfo1.covPV( 0, 0 ) ), std::sqrt( recinfo1.covPV( 1, 1 ) ),
+ std::sqrt( recinfo1.covPV( 2, 2 ) )};
+ recinfo1.ntracks = pv.nTracks();
+ if ( recinfo1.ntracks >= 10 ) { ++m_nRec_geq10tr_1; }
+ recinfo1.chi2 = pv.chi2();
+ recinfo1.nDoF = pv.nDoF();
+
+ if ( debugLevel() )
+ debug() << " " << pv.position().x() << " " << pv.position().y() << " " << pv.position().z()
+ << " " << pv.chi2PerDoF() << " " << pv.nTracks() << endmsg;
+ fullVrt1.push_back( recinfo1 );
+ } // end of loop over vertices1
+
+ for ( auto const& pv : recoVtx2 ) {
+ ++m_nRec2;
+ recinfo2.pRECPV = &pv;
+ recinfo2.position = pv.position();
+ recinfo2.covPV = pv.covMatrix();
+ recinfo2.positionSigma = Gaudi::XYZPoint{std::sqrt( recinfo2.covPV( 0, 0 ) ), std::sqrt( recinfo2.covPV( 1, 1 ) ),
+ std::sqrt( recinfo2.covPV( 2, 2 ) )};
+ recinfo2.ntracks = pv.nTracks();
+ if ( recinfo2.ntracks >= 10 ) { ++m_nRec_geq10tr_2; }
+ recinfo2.chi2 = pv.chi2();
+ recinfo2.nDoF = pv.nDoF();
+
+ if ( debugLevel() )
+ debug() << " " << pv.position().x() << " " << pv.position().y() << " " << pv.position().z()
+ << " " << pv.chi2PerDoF() << " " << pv.nTracks() << endmsg;
+ fullVrt2.push_back( recinfo2 );
+ } // end of loop over vertices2
+
+ // added sorting to mirror the 'multirec' variable implementation in PVChecker
+ std::sort( fullVrt1.begin(), fullVrt1.end(), trackcomp );
+ std::sort( fullVrt2.begin(), fullVrt2.end(), trackcomp );
+
+ if ( debugLevel() ) debug() << "fullVrt1 size " << fullVrt1.size() << endmsg;
+ if ( debugLevel() ) debug() << "fullVrt2 size " << fullVrt2.size() << endmsg;
+ size_diff = fullVrt1.size() - fullVrt2.size();
+
+ if ( m_produceNtuple.value() ) {
+ Tuple myTuple_evt = nTuple( "PV_nTuple_evt", "PV_nTuple_evt", CLID_ColumnWiseTuple );
+ myTuple_evt->column( "size_diff", double( size_diff ) ).ignore();
+ myTuple_evt->column( "size_1", double( fullVrt1.size() ) ).ignore();
+ myTuple_evt->column( "size_2", double( fullVrt2.size() ) ).ignore();
+ myTuple_evt->write().ignore();
+ }
std::vector<int> link;
- int ntracks1 = 0;
- int ntracks2 = 0;
- int dtracks = 0;
- int size1 = 0;
- int size2 = 0;
- double sigx_part1 = -99999.;
- double sigy_part1 = -99999.;
- double sigz_part1 = -99999.;
- double sigx_part2 = -99999.;
- double sigy_part2 = -99999.;
- double sigz_part2 = -99999.;
- double x1 = -99999.;
- double y1 = -99999.;
- double z1 = -99999.;
- double x2 = -99999.;
- double y2 = -99999.;
- double z2 = -99999.;
- double dx = -99999.;
- double dy = -99999.;
- double dz = -99999.;
- int oIt = 0;
- int pv_rank = -99999;
+ int ntracks1 = 0;
+ int ntracks2 = 0;
+ int dtracks = 0;
+ int size1 = 0;
+ int size2 = 0;
+ double sigx_part1 = -99999.;
+ double sigy_part1 = -99999.;
+ double sigz_part1 = -99999.;
+ double sigx_part2 = -99999.;
+ double sigy_part2 = -99999.;
+ double sigz_part2 = -99999.;
+ double x1 = -99999.;
+ double y1 = -99999.;
+ double z1 = -99999.;
+ double x2 = -99999.;
+ double y2 = -99999.;
+ double z2 = -99999.;
+ double dx = -99999.;
+ double dy = -99999.;
+ double dz = -99999.;
+ int oIt = 0;
+ int pv_rank = -99999;
if ( fullVrt1.size() != 0 && fullVrt2.size() != 0 ) { matchByDistance( fullVrt1, fullVrt2, link ); }
@@ -495,10 +499,10 @@ void VertexCompare::operator()( Vertices const& recoVtx1, Vertices const& recoVt
}
if ( m_monitoring.value() ) {
- if ( dz < 2 && size_diff < 3 && pv_rank < 2 && ( ntracks1 + ntracks2 ) / 2 > ntrack_bins[0]) {
+ if ( dz < 2 && size_diff < 3 && pv_rank < 2 && ( ntracks1 + ntracks2 ) / 2 > ntrack_bins[0] ) {
int binCount = 0;
for ( size_t i = 1; i < ntrack_bins.size(); ++i ) {
- if ( std::lround(( ntracks1 + ntracks2 ) / 2) <= ntrack_bins[i] ) { break; }
+ if ( std::lround( ( ntracks1 + ntracks2 ) / 2 ) <= ntrack_bins[i] ) { break; }
binCount++;
}
auto& monitoringHistos = *m_monitoringHistos.get();
@@ -508,15 +512,15 @@ void VertexCompare::operator()( Vertices const& recoVtx1, Vertices const& recoVt
++monitoringHistos.m_histo_pullx_Monitoring[dx / errx];
++monitoringHistos.m_histo_pully_Monitoring[dy / erry];
++monitoringHistos.m_histo_pullz_Monitoring[dz / errz];
- }
+ }
}
- if (m_produceHistogram.value()) {
- dx_vector.push_back(dx);
- dy_vector.push_back(dy);
- dz_vector.push_back(dz);
- pullx_vector.push_back(dx / errx);
- pully_vector.push_back(dy / erry);
- pullz_vector.push_back(dz / errz);
+ if ( m_produceHistogram.value() ) {
+ dx_vector.push_back( dx );
+ dy_vector.push_back( dy );
+ dz_vector.push_back( dz );
+ pullx_vector.push_back( dx / errx );
+ pully_vector.push_back( dy / erry );
+ pullz_vector.push_back( dz / errz );
++m_histo_dx[dx];
++m_histo_dy[dy];
++m_histo_dz[dz];
@@ -525,8 +529,8 @@ void VertexCompare::operator()( Vertices const& recoVtx1, Vertices const& recoVt
++m_histo_pullz[dz / errz];
++m_nTracks1[ntracks1];
++m_nTracks2[ntracks2];
- ++m_nTracks_dif[ntracks2-ntracks1];
- } //else end
+ ++m_nTracks_dif[ntracks2 - ntracks1];
+ } // else end
if ( m_produceNtuple.value() ) {
Tuple myTuple = nTuple( "PV_nTuple", "PV_nTuple", CLID_ColumnWiseTuple );
@@ -692,7 +696,7 @@ void VertexCompare::operator()( Vertices const& recoVtx1, Vertices const& recoVt
oIt++;
m_nEvt += 1;
}
- }
+}
//=============================================================================
// Finalize
@@ -729,85 +733,70 @@ StatusCode VertexCompare::finalize() {
<< ( (double)m_nVtx.nEntries() - (double)m_nVtx_geq10tr_matched2.nEntries() ) /
( (double)m_nRec2.nEntries() - (double)m_nRec_geq10tr_2.nEntries() ) * 100
<< " % " << endmsg;
-
- if(m_produceHistogram.value()){
- for (double value : dx_vector) {
- ++m_histo_dx[value];
- }
- for (double value : dy_vector) {
- ++m_histo_dy[value];
- }
- for (double value : dz_vector) {
- ++m_histo_dz[value];
- }
- for (double value : pullx_vector) {
- ++m_histo_pullx[value];
- }
- for (double value : pully_vector) {
- ++m_histo_pully[value];
- }
- for (double value : pullz_vector) {
- ++m_histo_pullz[value];
- }
+ if ( m_produceHistogram.value() ) {
- info() << "dx: "
- << format( "mean = %5.3f +/- %5.3f, RMS = %5.3f +/- %5.3f", m_histo_dx.mean(),
- meanErr( m_histo_dx ), m_histo_dx.standard_deviation(), rmsErr( m_histo_dx) )
+ for ( double value : dx_vector ) { ++m_histo_dx[value]; }
+ for ( double value : dy_vector ) { ++m_histo_dy[value]; }
+ for ( double value : dz_vector ) { ++m_histo_dz[value]; }
+ for ( double value : pullx_vector ) { ++m_histo_pullx[value]; }
+ for ( double value : pully_vector ) { ++m_histo_pully[value]; }
+ for ( double value : pullz_vector ) { ++m_histo_pullz[value]; }
+
+ info() << "dx: "
+ << format( "mean = %5.3f +/- %5.3f, RMS = %5.3f +/- %5.3f", m_histo_dx.mean(), meanErr( m_histo_dx ),
+ m_histo_dx.standard_deviation(), rmsErr( m_histo_dx ) )
<< endmsg;
-
-
- info() << "dy: "
- << format( "mean = %5.3f +/- %5.3f, RMS = %5.3f +/- %5.3f", m_histo_dy.mean(),
- meanErr( m_histo_dy ), m_histo_dy.standard_deviation(), rmsErr( m_histo_dy) )
+
+ info() << "dy: "
+ << format( "mean = %5.3f +/- %5.3f, RMS = %5.3f +/- %5.3f", m_histo_dy.mean(), meanErr( m_histo_dy ),
+ m_histo_dy.standard_deviation(), rmsErr( m_histo_dy ) )
<< endmsg;
-
- info() << "dz: "
- << format( "mean = %5.3f +/- %5.3f, RMS = %5.3f +/- %5.3f", m_histo_dz.mean(),
- meanErr( m_histo_dz ), m_histo_dz.standard_deviation(), rmsErr( m_histo_dz) )
+ info() << "dz: "
+ << format( "mean = %5.3f +/- %5.3f, RMS = %5.3f +/- %5.3f", m_histo_dz.mean(), meanErr( m_histo_dz ),
+ m_histo_dz.standard_deviation(), rmsErr( m_histo_dz ) )
<< endmsg;
-
- info() << " ---------------------------------------" << endmsg;
- info() << "pullx: "
- << format( "mean = %5.3f +/- %5.3f, RMS = %5.3f +/- %5.3f", m_histo_pullx.mean(),
- meanErr( m_histo_pullx ), m_histo_pullx.standard_deviation(), rmsErr( m_histo_pullx ) )
+ info() << " ---------------------------------------" << endmsg;
+
+ info() << "pullx: "
+ << format( "mean = %5.3f +/- %5.3f, RMS = %5.3f +/- %5.3f", m_histo_pullx.mean(), meanErr( m_histo_pullx ),
+ m_histo_pullx.standard_deviation(), rmsErr( m_histo_pullx ) )
<< endmsg;
-
- info() << "pully: "
- << format( "mean = %5.3f +/- %5.3f, RMS = %5.3f +/- %5.3f", m_histo_pully.mean(),
- meanErr( m_histo_pully ), m_histo_pully.standard_deviation(), rmsErr( m_histo_pully) )
+
+ info() << "pully: "
+ << format( "mean = %5.3f +/- %5.3f, RMS = %5.3f +/- %5.3f", m_histo_pully.mean(), meanErr( m_histo_pully ),
+ m_histo_pully.standard_deviation(), rmsErr( m_histo_pully ) )
<< endmsg;
-
- info() << "pullz: "
- << format( "mean = %5.3f +/- %5.3f, RMS = %5.3f +/- %5.3f", m_histo_pullz.mean(),
- meanErr( m_histo_pullz ), m_histo_pullz.standard_deviation(), rmsErr( m_histo_pullz) )
+
+ info() << "pullz: "
+ << format( "mean = %5.3f +/- %5.3f, RMS = %5.3f +/- %5.3f", m_histo_pullz.mean(), meanErr( m_histo_pullz ),
+ m_histo_pullz.standard_deviation(), rmsErr( m_histo_pullz ) )
<< endmsg;
-
- info() << " ---------------------------------------" << endmsg;
- info() << "diff in x: "
+ info() << " ---------------------------------------" << endmsg;
+
+ info() << "diff in x: "
<< format( "%5.3f +/- %5.3f", std::sqrt( 1.0 + m_histo_pullx.standard_deviation() ) - 1.0,
rmsErr( m_histo_pullx ) * 0.5 / std::sqrt( 1.0 + m_histo_pullx.standard_deviation() ) )
<< endmsg;
- info() << "diff in y: "
+ info() << "diff in y: "
<< format( "%5.3f +/- %5.3f", std::sqrt( 1.0 + m_histo_pully.standard_deviation() ) - 1.0,
rmsErr( m_histo_pully ) * 0.5 / std::sqrt( 1.0 + m_histo_pully.standard_deviation() ) )
<< endmsg;
- info() << "diff in z: "
+ info() << "diff in z: "
<< format( "%5.3f +/- %5.3f", std::sqrt( 1.0 + m_histo_pullz.standard_deviation() ) - 1.0,
rmsErr( m_histo_pullz ) * 0.5 / std::sqrt( 1.0 + m_histo_pullz.standard_deviation() ) )
<< endmsg;
- info() << " ============================================" << endmsg;
+ info() << " ============================================" << endmsg;
}
return Consumer::finalize(); // Must be called after all other actions
}
-
//=============================================================================
// Match vertices by distance
//=============================================================================
--
GitLab
From 4a7a8085c7747dfe26dbfa6f4561736ee057a990 Mon Sep 17 00:00:00 2001
From: Bogdan Kutsenko <bogdan.kutsenko@cern.ch>
Date: Sat, 25 May 2024 09:45:59 +0200
Subject: [PATCH 11/21] Random split of velo track for Moore TBLV resolution
monitor
---
Tr/TrackMonitors/src/VertexCompare.cpp | 411 +++++++++---------
Tr/TrackUtils/CMakeLists.txt | 1 +
.../src/RandomTrackContainerSplitter.cpp | 19 +-
.../src/RandomVeloTrackContainerSplitter.cpp | 93 ++++
4 files changed, 308 insertions(+), 216 deletions(-)
create mode 100644 Tr/TrackUtils/src/RandomVeloTrackContainerSplitter.cpp
diff --git a/Tr/TrackMonitors/src/VertexCompare.cpp b/Tr/TrackMonitors/src/VertexCompare.cpp
index b8294a5d5dc..360d24fc961 100644
--- a/Tr/TrackMonitors/src/VertexCompare.cpp
+++ b/Tr/TrackMonitors/src/VertexCompare.cpp
@@ -89,31 +89,31 @@ bool trackcomp( const RecPVInfo<VertexType>& first, const RecPVInfo<VertexType>&
}
struct VtxVariables {
- int ntracks = -9999;
- double x = -9999.;
- double y = -9999.;
- double z = -9999.;
- double dxr = -9999.;
- double dyr = -9999.;
- double r = -9999.;
- double errx = -9999.;
- double erry = -9999.;
- double errz = -9999.;
- double err_r = -9999.;
- double covxx = -9999.;
- double covyy = -9999.;
- double covzz = -9999.;
- double covxy = -9999.;
- double covxz = -9999.;
- double covyz = -9999.;
- double chi2 = -9999.;
- double nDoF = -9999.;
- int dsize = -9999;
- int match = -9999;
- int pv_rank = -9999;
- bool equal_sizes = false;
- bool single = false;
- bool opposite_container = false;
+ int ntracks = -9999;
+ double x = -9999.;
+ double y = -9999.;
+ double z = -9999.;
+ double dxr = -9999.;
+ double dyr = -9999.;
+ double r = -9999.;
+ double errx = -9999.;
+ double erry = -9999.;
+ double errz = -9999.;
+ double err_r = -9999.;
+ double covxx = -9999.;
+ double covyy = -9999.;
+ double covzz = -9999.;
+ double covxy = -9999.;
+ double covxz = -9999.;
+ double covyz = -9999.;
+ double chi2 = -9999.;
+ double nDoF = -9999.;
+ int dsize = -9999;
+ int match = -9999;
+ int pv_rank = -9999;
+ bool equal_sizes = false;
+ bool single = false;
+ bool opposite_container = false;
std::vector<double> dx_vector{}, dy_vector{}, dz_vector{}, pullx_vector{}, pully_vector{}, pullz_vector{};
};
@@ -177,7 +177,7 @@ VtxVariables SetVtxVariables( const RecPVInfo<VertexType>& vrtf, const int& coun
// ============================================================================
template <typename Histogram>
-double meanErr( Histogram const& h ) {
+double meanErr( Histogram const & h ) {
const auto n = h.nEntries();
return ( 0 >= n ? 0.0 : h.standard_deviation() / std::sqrt( n ) );
}
@@ -186,28 +186,28 @@ double meanErr( Histogram const& h ) {
// Fourth momentum calculation
// ============================================================================
template <typename Histogram>
-double fourthMomentum( Histogram const& h ) {
+double fourthMomentum( Histogram const & h ) {
const auto n = h.nEntries();
if ( 0 >= n ) { return 0.0; } // RETURN
// get the mean
const auto h_mean = h.mean();
// number of bins
- const auto& axis = h.axis();
- const int nBins = h.totNBins();
- double result{0}, weight{0};
+ const auto& axis = h.axis();
+ const int nBins = h.totNBins();
+ double result{ 0 }, weight{ 0 };
// loop over all bins
- double minVal = axis[0].minValue;
- double maxVal = axis[0].maxValue;
- double binWidth = ( maxVal - minVal ) / nBins;
+ double minVal = axis[0].minValue;
+ double maxVal = axis[0].maxValue;
+ double binWidth = (maxVal-minVal)/nBins;
std::vector<double> binCenters;
- for ( int i = 0; i < nBins; ++i ) {
- double center = minVal + ( i + 0.5 ) * binWidth;
- binCenters.push_back( center );
- }
- for ( int i = 0; i < nBins; ++i ) {
- const auto yBin = h.binValue( i ); // bin weight
+ for (int i = 0; i < nBins; ++i) {
+ double center = minVal + (i + 0.5) * binWidth;
+ binCenters.push_back(center);
+ }
+ for (int i = 0; i < nBins; ++i) {
+ const auto yBin = h.binValue( i ); // bin weight
weight += yBin;
result += yBin * std::pow( binCenters[i] - h_mean, 4 );
}
@@ -221,25 +221,26 @@ double fourthMomentum( Histogram const& h ) {
// ============================================================================
template <typename Histogram>
-double kurtosis( Histogram const& h ) {
+double kurtosis( Histogram const & h ) {
const auto mu4 = fourthMomentum( h );
- const auto s4 = std::pow( h.standard_deviation(), 4 );
+ const auto s4 = std::pow( h.standard_deviation() , 4 );
return ( std::fabs( s4 ) > 0 ? mu4 / s4 - 3.0 : 0.0 );
}
// ============================================================================
// get an error in the rms value
// ============================================================================
+
+ template <typename Histogram>
+ double rmsErr( Histogram const & h ) {
+ const auto n = h.nEntries();
+ if ( 1 >= n ) { return 0.0; }
+ auto result = 2.0 + kurtosis( h );
+ result += 2.0 / ( n - 1 );
+ result /= 4.0 * n;
+ return h.standard_deviation() * std::sqrt( std::max( result, 0.0 ) );
+ }
-template <typename Histogram>
-double rmsErr( Histogram const& h ) {
- const auto n = h.nEntries();
- if ( 1 >= n ) { return 0.0; }
- auto result = 2.0 + kurtosis( h );
- result += 2.0 / ( n - 1 );
- result /= 4.0 * n;
- return h.standard_deviation() * std::sqrt( std::max( result, 0.0 ) );
-}
class VertexCompare : public LHCb::Algorithm::Consumer<
void( Vertices const&, Vertices const&, LHCb::Conditions::InteractionRegion const& ),
@@ -262,7 +263,7 @@ public:
private:
bool debugLevel() const { return msgLevel( MSG::DEBUG ) || msgLevel( MSG::VERBOSE ); }
- static constexpr auto ntrack_bins = std::array{3.5, 16.5, 29.5, 43.5, 56.5, 69.5};
+ static constexpr auto ntrack_bins = std::array{ 3.5 , 16.5 , 29.5 , 43.5 , 56.5 , 69.5 };
Gaudi::Property<bool> m_produceHistogram{this, "produceHistogram", true};
Gaudi::Property<bool> m_produceNtuple{this, "produceNtuple", true};
@@ -270,7 +271,7 @@ private:
Gaudi::Property<bool> m_requireSingle{this, "requireSingle", false};
mutable std::vector<double> dx_vector{}, dy_vector{}, dz_vector{}, pullx_vector{}, pully_vector{}, pullz_vector{};
-
+
struct monitoringHistos {
mutable std::array<Gaudi::Accumulators::Histogram<1>, 5> m_histo_nTracksBins_dx;
mutable std::array<Gaudi::Accumulators::Histogram<1>, 5> m_histo_nTracksBins_dy;
@@ -279,6 +280,7 @@ private:
mutable Gaudi::Accumulators::Histogram<1> m_histo_pully_Monitoring;
mutable Gaudi::Accumulators::Histogram<1> m_histo_pullz_Monitoring;
+
template <std::size_t... IDXs>
static std::array<Gaudi::Accumulators::Histogram<1>, sizeof...( IDXs )>
histo1DArrayBuilder( const VertexCompare* owner, const std::string& name, const std::string& title,
@@ -288,34 +290,32 @@ private:
title + std::to_string( IDXs ),
{std::get<0>( xbins ), std::get<1>( xbins ), std::get<2>( xbins )}}...}};
}
+ Gaudi::Histo1DDef m_histodx{"dx, mm", -0.15, 0.15, 50};
monitoringHistos( const VertexCompare* owner )
- : m_histo_nTracksBins_dx{histo1DArrayBuilder( owner, "dx_Monitoring_", "dx, mm", {50, -0.15, 0.15},
- std::make_index_sequence<5>() )}
- , m_histo_nTracksBins_dy{histo1DArrayBuilder( owner, "dy_Monitoring_", "dy, mm", {50, -0.15, 0.15},
- std::make_index_sequence<5>() )}
- , m_histo_nTracksBins_dz{histo1DArrayBuilder( owner, "dz_Monitoring_", "dz, mm", {50, -1.5, 1.5},
- std::make_index_sequence<5>() )}
+ : m_histo_nTracksBins_dx{histo1DArrayBuilder( owner, "dx_Monitoring_", "dx, mm", {50, -0.15, 0.15}, std::make_index_sequence<5>() )}
+ , m_histo_nTracksBins_dy{histo1DArrayBuilder( owner, "dy_Monitoring_", "dy, mm", {50, -0.15, 0.15}, std::make_index_sequence<5>() )}
+ , m_histo_nTracksBins_dz{histo1DArrayBuilder( owner, "dz_Monitoring_", "dz, mm", {50, -1.5, 1.5}, std::make_index_sequence<5>() )}
, m_histo_pullx_Monitoring{owner, "pullx_Monitoring", "pull x", {20, -5, 5}}
, m_histo_pully_Monitoring{owner, "pully_Monitoring", "pull y", {20, -5, 5}}
, m_histo_pullz_Monitoring{owner, "pullz_Monitoring", "pull z", {20, -5, 5}} {}
+
+
+
};
std::unique_ptr<monitoringHistos> m_monitoringHistos;
using axis1D = Gaudi::Accumulators::Axis<Gaudi::Accumulators::Histogram<1>::AxisArithmeticType>;
- mutable Gaudi::Accumulators::RootHistogram<1> m_histo_dx{this, "dx", "dx, mm", axis1D{50, -0.15, 0.15}};
- mutable Gaudi::Accumulators::RootHistogram<1> m_histo_dy{this, "dy", "dy, mm", axis1D{50, -0.15, 0.15}};
- mutable Gaudi::Accumulators::RootHistogram<1> m_histo_dz{this, "dz", "dz, mm", axis1D{50, -1.5, 1.5}};
- mutable Gaudi::Accumulators::RootHistogram<1> m_histo_pullx{this, "pullx", "pull x", axis1D{20, -5, 5}};
- mutable Gaudi::Accumulators::RootHistogram<1> m_histo_pully{this, "pully", "pull y", axis1D{20, -5, 5}};
- mutable Gaudi::Accumulators::RootHistogram<1> m_histo_pullz{this, "pullz", "pull z", axis1D{20, -5, 5}};
- mutable Gaudi::Accumulators::Histogram<1> m_nTracks1{this, "ntracks1", "Number of tracks in PV set 1",
- axis1D{50, 0., 150.}};
- mutable Gaudi::Accumulators::Histogram<1> m_nTracks2{this, "ntracks2", "Number of tracks in PV set 2",
- axis1D{50, 0., 150.}};
- mutable Gaudi::Accumulators::Histogram<1> m_nTracks_dif{
- this, "dtracks", "Difference in the number of tracks in two sets of PVs", axis1D{50, 0., 150.}};
-
+ mutable Gaudi::Accumulators::RootHistogram<1> m_histo_dx{this, "dx", "dx, mm", axis1D{50, -0.15, 0.15} };
+ mutable Gaudi::Accumulators::RootHistogram<1> m_histo_dy{ this, "dy", "dy, mm", axis1D{50, -0.15, 0.15} };
+ mutable Gaudi::Accumulators::RootHistogram<1> m_histo_dz{ this, "dz", "dz, mm", axis1D{50, -1.5, 1.5} };
+ mutable Gaudi::Accumulators::RootHistogram<1> m_histo_pullx{ this, "pullx", "pull x", axis1D{20, -5, 5} };
+ mutable Gaudi::Accumulators::RootHistogram<1> m_histo_pully{ this, "pully", "pull y", axis1D{20, -5, 5} };
+ mutable Gaudi::Accumulators::RootHistogram<1> m_histo_pullz{ this, "pullz", "pull z", axis1D{20, -5, 5} };
+ mutable Gaudi::Accumulators::Histogram<1> m_nTracks1{ this, "ntracks1", "Number of tracks in PV set 1", axis1D{50, 0., 150.} };
+ mutable Gaudi::Accumulators::Histogram<1> m_nTracks2{ this, "ntracks2", "Number of tracks in PV set 2", axis1D{ 50, 0., 150.}};
+ mutable Gaudi::Accumulators::Histogram<1> m_nTracks_dif{ this, "dtracks", "Difference in the number of tracks in two sets of PVs", axis1D{ 50, 0., 150.} };
+
mutable Gaudi::Accumulators::Counter<> m_nVtx{this, "Number of pairs of vertices in processed events"};
mutable Gaudi::Accumulators::Counter<> m_nRec1{this, "Number of vertices in input1"};
mutable Gaudi::Accumulators::Counter<> m_nRec2{this, "Number of vertices in input2"};
@@ -355,96 +355,96 @@ StatusCode VertexCompare::initialize() {
//=============================================================================
void VertexCompare::operator()( Vertices const& recoVtx1, Vertices const& recoVtx2,
LHCb::Conditions::InteractionRegion const& region ) const {
- if ( msgLevel( MSG::DEBUG ) ) debug() << "==> Execute" << endmsg;
+ if ( msgLevel( MSG::DEBUG ) ) debug() << "==> Execute" << endmsg;
const auto beamspot = region.avgPosition;
std::vector<RecPVInfo<VertexType>> fullVrt1, fullVrt2;
fullVrt1.reserve( recoVtx1.size() );
fullVrt2.reserve( recoVtx2.size() );
- if ( m_requireSingle == true && !( recoVtx1.size() == 1 && recoVtx2.size() == 1 ) ) return;
-
+ if(m_requireSingle == true && !(recoVtx1.size() == 1 && recoVtx2.size() == 1)) return;
+
if ( debugLevel() ) debug() << " Vtx Properities x y z chi2/ndof ntracks" << endmsg;
- int size_diff = -9999;
+ int size_diff = -9999;
RecPVInfo<VertexType> recinfo1;
RecPVInfo<VertexType> recinfo2;
- for ( auto const& pv : recoVtx1 ) {
- ++m_nRec1;
- recinfo1.pRECPV = &pv;
- recinfo1.position = pv.position();
- recinfo1.covPV = pv.covMatrix();
- recinfo1.positionSigma = Gaudi::XYZPoint{std::sqrt( recinfo1.covPV( 0, 0 ) ), std::sqrt( recinfo1.covPV( 1, 1 ) ),
- std::sqrt( recinfo1.covPV( 2, 2 ) )};
- recinfo1.ntracks = pv.nTracks();
- if ( recinfo1.ntracks >= 10 ) { ++m_nRec_geq10tr_1; }
- recinfo1.chi2 = pv.chi2();
- recinfo1.nDoF = pv.nDoF();
-
- if ( debugLevel() )
- debug() << " " << pv.position().x() << " " << pv.position().y() << " " << pv.position().z()
- << " " << pv.chi2PerDoF() << " " << pv.nTracks() << endmsg;
- fullVrt1.push_back( recinfo1 );
- } // end of loop over vertices1
-
- for ( auto const& pv : recoVtx2 ) {
- ++m_nRec2;
- recinfo2.pRECPV = &pv;
- recinfo2.position = pv.position();
- recinfo2.covPV = pv.covMatrix();
- recinfo2.positionSigma = Gaudi::XYZPoint{std::sqrt( recinfo2.covPV( 0, 0 ) ), std::sqrt( recinfo2.covPV( 1, 1 ) ),
- std::sqrt( recinfo2.covPV( 2, 2 ) )};
- recinfo2.ntracks = pv.nTracks();
- if ( recinfo2.ntracks >= 10 ) { ++m_nRec_geq10tr_2; }
- recinfo2.chi2 = pv.chi2();
- recinfo2.nDoF = pv.nDoF();
-
- if ( debugLevel() )
- debug() << " " << pv.position().x() << " " << pv.position().y() << " " << pv.position().z()
- << " " << pv.chi2PerDoF() << " " << pv.nTracks() << endmsg;
- fullVrt2.push_back( recinfo2 );
- } // end of loop over vertices2
-
- // added sorting to mirror the 'multirec' variable implementation in PVChecker
- std::sort( fullVrt1.begin(), fullVrt1.end(), trackcomp );
- std::sort( fullVrt2.begin(), fullVrt2.end(), trackcomp );
-
- if ( debugLevel() ) debug() << "fullVrt1 size " << fullVrt1.size() << endmsg;
- if ( debugLevel() ) debug() << "fullVrt2 size " << fullVrt2.size() << endmsg;
- size_diff = fullVrt1.size() - fullVrt2.size();
-
- if ( m_produceNtuple.value() ) {
- Tuple myTuple_evt = nTuple( "PV_nTuple_evt", "PV_nTuple_evt", CLID_ColumnWiseTuple );
- myTuple_evt->column( "size_diff", double( size_diff ) ).ignore();
- myTuple_evt->column( "size_1", double( fullVrt1.size() ) ).ignore();
- myTuple_evt->column( "size_2", double( fullVrt2.size() ) ).ignore();
- myTuple_evt->write().ignore();
- }
+ for ( auto const& pv : recoVtx1 ) {
+ ++m_nRec1;
+ recinfo1.pRECPV = &pv;
+ recinfo1.position = pv.position();
+ recinfo1.covPV = pv.covMatrix();
+ recinfo1.positionSigma = Gaudi::XYZPoint{std::sqrt( recinfo1.covPV( 0, 0 ) ), std::sqrt( recinfo1.covPV( 1, 1 ) ),
+ std::sqrt( recinfo1.covPV( 2, 2 ) )};
+ recinfo1.ntracks = pv.nTracks();
+ if ( recinfo1.ntracks >= 10 ) { ++m_nRec_geq10tr_1; }
+ recinfo1.chi2 = pv.chi2();
+ recinfo1.nDoF = pv.nDoF();
+
+ if ( debugLevel() )
+ debug() << " " << pv.position().x() << " " << pv.position().y() << " " << pv.position().z()
+ << " " << pv.chi2PerDoF() << " " << pv.nTracks() << endmsg;
+ fullVrt1.push_back( recinfo1 );
+ } // end of loop over vertices1
+
+ for ( auto const& pv : recoVtx2 ) {
+ ++m_nRec2;
+ recinfo2.pRECPV = &pv;
+ recinfo2.position = pv.position();
+ recinfo2.covPV = pv.covMatrix();
+ recinfo2.positionSigma = Gaudi::XYZPoint{std::sqrt( recinfo2.covPV( 0, 0 ) ), std::sqrt( recinfo2.covPV( 1, 1 ) ),
+ std::sqrt( recinfo2.covPV( 2, 2 ) )};
+ recinfo2.ntracks = pv.nTracks();
+ if ( recinfo2.ntracks >= 10 ) { ++m_nRec_geq10tr_2; }
+ recinfo2.chi2 = pv.chi2();
+ recinfo2.nDoF = pv.nDoF();
+
+ if ( debugLevel() )
+ debug() << " " << pv.position().x() << " " << pv.position().y() << " " << pv.position().z()
+ << " " << pv.chi2PerDoF() << " " << pv.nTracks() << endmsg;
+ fullVrt2.push_back( recinfo2 );
+ } // end of loop over vertices2
+
+ // added sorting to mirror the 'multirec' variable implementation in PVChecker
+ std::sort( fullVrt1.begin(), fullVrt1.end(), trackcomp );
+ std::sort( fullVrt2.begin(), fullVrt2.end(), trackcomp );
+
+ if ( debugLevel() ) debug() << "fullVrt1 size " << fullVrt1.size() << endmsg;
+ if ( debugLevel() ) debug() << "fullVrt2 size " << fullVrt2.size() << endmsg;
+ size_diff = fullVrt1.size() - fullVrt2.size();
+
+ if ( m_produceNtuple.value() ) {
+ Tuple myTuple_evt = nTuple( "PV_nTuple_evt", "PV_nTuple_evt", CLID_ColumnWiseTuple );
+ myTuple_evt->column( "size_diff", double( size_diff ) ).ignore();
+ myTuple_evt->column( "size_1", double( fullVrt1.size() ) ).ignore();
+ myTuple_evt->column( "size_2", double( fullVrt2.size() ) ).ignore();
+ myTuple_evt->write().ignore();
+ }
std::vector<int> link;
- int ntracks1 = 0;
- int ntracks2 = 0;
- int dtracks = 0;
- int size1 = 0;
- int size2 = 0;
- double sigx_part1 = -99999.;
- double sigy_part1 = -99999.;
- double sigz_part1 = -99999.;
- double sigx_part2 = -99999.;
- double sigy_part2 = -99999.;
- double sigz_part2 = -99999.;
- double x1 = -99999.;
- double y1 = -99999.;
- double z1 = -99999.;
- double x2 = -99999.;
- double y2 = -99999.;
- double z2 = -99999.;
- double dx = -99999.;
- double dy = -99999.;
- double dz = -99999.;
- int oIt = 0;
- int pv_rank = -99999;
+ int ntracks1 = 0;
+ int ntracks2 = 0;
+ int dtracks = 0;
+ int size1 = 0;
+ int size2 = 0;
+ double sigx_part1 = -99999.;
+ double sigy_part1 = -99999.;
+ double sigz_part1 = -99999.;
+ double sigx_part2 = -99999.;
+ double sigy_part2 = -99999.;
+ double sigz_part2 = -99999.;
+ double x1 = -99999.;
+ double y1 = -99999.;
+ double z1 = -99999.;
+ double x2 = -99999.;
+ double y2 = -99999.;
+ double z2 = -99999.;
+ double dx = -99999.;
+ double dy = -99999.;
+ double dz = -99999.;
+ int oIt = 0;
+ int pv_rank = -99999;
if ( fullVrt1.size() != 0 && fullVrt2.size() != 0 ) { matchByDistance( fullVrt1, fullVrt2, link ); }
@@ -499,10 +499,10 @@ void VertexCompare::operator()( Vertices const& recoVtx1, Vertices const& recoVt
}
if ( m_monitoring.value() ) {
- if ( dz < 2 && size_diff < 3 && pv_rank < 2 && ( ntracks1 + ntracks2 ) / 2 > ntrack_bins[0] ) {
+ if ( dz < 2 && size_diff < 3 && pv_rank < 2 && ( ntracks1 + ntracks2 ) / 2 > ntrack_bins[0]) {
int binCount = 0;
for ( size_t i = 1; i < ntrack_bins.size(); ++i ) {
- if ( std::lround( ( ntracks1 + ntracks2 ) / 2 ) <= ntrack_bins[i] ) { break; }
+ if ( std::lround(( ntracks1 + ntracks2 ) / 2) <= ntrack_bins[i] ) { break; }
binCount++;
}
auto& monitoringHistos = *m_monitoringHistos.get();
@@ -512,15 +512,15 @@ void VertexCompare::operator()( Vertices const& recoVtx1, Vertices const& recoVt
++monitoringHistos.m_histo_pullx_Monitoring[dx / errx];
++monitoringHistos.m_histo_pully_Monitoring[dy / erry];
++monitoringHistos.m_histo_pullz_Monitoring[dz / errz];
- }
+ }
}
- if ( m_produceHistogram.value() ) {
- dx_vector.push_back( dx );
- dy_vector.push_back( dy );
- dz_vector.push_back( dz );
- pullx_vector.push_back( dx / errx );
- pully_vector.push_back( dy / erry );
- pullz_vector.push_back( dz / errz );
+ if (m_produceHistogram.value()) {
+ dx_vector.push_back(dx);
+ dy_vector.push_back(dy);
+ dz_vector.push_back(dz);
+ pullx_vector.push_back(dx / errx);
+ pully_vector.push_back(dy / erry);
+ pullz_vector.push_back(dz / errz);
++m_histo_dx[dx];
++m_histo_dy[dy];
++m_histo_dz[dz];
@@ -529,8 +529,8 @@ void VertexCompare::operator()( Vertices const& recoVtx1, Vertices const& recoVt
++m_histo_pullz[dz / errz];
++m_nTracks1[ntracks1];
++m_nTracks2[ntracks2];
- ++m_nTracks_dif[ntracks2 - ntracks1];
- } // else end
+ ++m_nTracks_dif[ntracks2-ntracks1];
+ } //else end
if ( m_produceNtuple.value() ) {
Tuple myTuple = nTuple( "PV_nTuple", "PV_nTuple", CLID_ColumnWiseTuple );
@@ -696,7 +696,7 @@ void VertexCompare::operator()( Vertices const& recoVtx1, Vertices const& recoVt
oIt++;
m_nEvt += 1;
}
-}
+ }
//=============================================================================
// Finalize
@@ -733,70 +733,85 @@ StatusCode VertexCompare::finalize() {
<< ( (double)m_nVtx.nEntries() - (double)m_nVtx_geq10tr_matched2.nEntries() ) /
( (double)m_nRec2.nEntries() - (double)m_nRec_geq10tr_2.nEntries() ) * 100
<< " % " << endmsg;
+
+ if(m_produceHistogram.value()){
- if ( m_produceHistogram.value() ) {
-
- for ( double value : dx_vector ) { ++m_histo_dx[value]; }
- for ( double value : dy_vector ) { ++m_histo_dy[value]; }
- for ( double value : dz_vector ) { ++m_histo_dz[value]; }
- for ( double value : pullx_vector ) { ++m_histo_pullx[value]; }
- for ( double value : pully_vector ) { ++m_histo_pully[value]; }
- for ( double value : pullz_vector ) { ++m_histo_pullz[value]; }
+ for (double value : dx_vector) {
+ ++m_histo_dx[value];
+ }
+ for (double value : dy_vector) {
+ ++m_histo_dy[value];
+ }
+ for (double value : dz_vector) {
+ ++m_histo_dz[value];
+ }
+ for (double value : pullx_vector) {
+ ++m_histo_pullx[value];
+ }
+ for (double value : pully_vector) {
+ ++m_histo_pully[value];
+ }
+ for (double value : pullz_vector) {
+ ++m_histo_pullz[value];
+ }
- info() << "dx: "
- << format( "mean = %5.3f +/- %5.3f, RMS = %5.3f +/- %5.3f", m_histo_dx.mean(), meanErr( m_histo_dx ),
- m_histo_dx.standard_deviation(), rmsErr( m_histo_dx ) )
+ info() << "dx: "
+ << format( "mean = %5.3f +/- %5.3f, RMS = %5.3f +/- %5.3f", m_histo_dx.mean(),
+ meanErr( m_histo_dx ), m_histo_dx.standard_deviation(), rmsErr( m_histo_dx) )
<< endmsg;
-
- info() << "dy: "
- << format( "mean = %5.3f +/- %5.3f, RMS = %5.3f +/- %5.3f", m_histo_dy.mean(), meanErr( m_histo_dy ),
- m_histo_dy.standard_deviation(), rmsErr( m_histo_dy ) )
+
+
+ info() << "dy: "
+ << format( "mean = %5.3f +/- %5.3f, RMS = %5.3f +/- %5.3f", m_histo_dy.mean(),
+ meanErr( m_histo_dy ), m_histo_dy.standard_deviation(), rmsErr( m_histo_dy) )
<< endmsg;
+
- info() << "dz: "
- << format( "mean = %5.3f +/- %5.3f, RMS = %5.3f +/- %5.3f", m_histo_dz.mean(), meanErr( m_histo_dz ),
- m_histo_dz.standard_deviation(), rmsErr( m_histo_dz ) )
+ info() << "dz: "
+ << format( "mean = %5.3f +/- %5.3f, RMS = %5.3f +/- %5.3f", m_histo_dz.mean(),
+ meanErr( m_histo_dz ), m_histo_dz.standard_deviation(), rmsErr( m_histo_dz) )
<< endmsg;
+
+ info() << " ---------------------------------------" << endmsg;
- info() << " ---------------------------------------" << endmsg;
-
- info() << "pullx: "
- << format( "mean = %5.3f +/- %5.3f, RMS = %5.3f +/- %5.3f", m_histo_pullx.mean(), meanErr( m_histo_pullx ),
- m_histo_pullx.standard_deviation(), rmsErr( m_histo_pullx ) )
+ info() << "pullx: "
+ << format( "mean = %5.3f +/- %5.3f, RMS = %5.3f +/- %5.3f", m_histo_pullx.mean(),
+ meanErr( m_histo_pullx ), m_histo_pullx.standard_deviation(), rmsErr( m_histo_pullx ) )
<< endmsg;
-
- info() << "pully: "
- << format( "mean = %5.3f +/- %5.3f, RMS = %5.3f +/- %5.3f", m_histo_pully.mean(), meanErr( m_histo_pully ),
- m_histo_pully.standard_deviation(), rmsErr( m_histo_pully ) )
+
+ info() << "pully: "
+ << format( "mean = %5.3f +/- %5.3f, RMS = %5.3f +/- %5.3f", m_histo_pully.mean(),
+ meanErr( m_histo_pully ), m_histo_pully.standard_deviation(), rmsErr( m_histo_pully) )
<< endmsg;
-
- info() << "pullz: "
- << format( "mean = %5.3f +/- %5.3f, RMS = %5.3f +/- %5.3f", m_histo_pullz.mean(), meanErr( m_histo_pullz ),
- m_histo_pullz.standard_deviation(), rmsErr( m_histo_pullz ) )
+
+ info() << "pullz: "
+ << format( "mean = %5.3f +/- %5.3f, RMS = %5.3f +/- %5.3f", m_histo_pullz.mean(),
+ meanErr( m_histo_pullz ), m_histo_pullz.standard_deviation(), rmsErr( m_histo_pullz) )
<< endmsg;
+
+ info() << " ---------------------------------------" << endmsg;
- info() << " ---------------------------------------" << endmsg;
-
- info() << "diff in x: "
+ info() << "diff in x: "
<< format( "%5.3f +/- %5.3f", std::sqrt( 1.0 + m_histo_pullx.standard_deviation() ) - 1.0,
rmsErr( m_histo_pullx ) * 0.5 / std::sqrt( 1.0 + m_histo_pullx.standard_deviation() ) )
<< endmsg;
- info() << "diff in y: "
+ info() << "diff in y: "
<< format( "%5.3f +/- %5.3f", std::sqrt( 1.0 + m_histo_pully.standard_deviation() ) - 1.0,
rmsErr( m_histo_pully ) * 0.5 / std::sqrt( 1.0 + m_histo_pully.standard_deviation() ) )
<< endmsg;
- info() << "diff in z: "
+ info() << "diff in z: "
<< format( "%5.3f +/- %5.3f", std::sqrt( 1.0 + m_histo_pullz.standard_deviation() ) - 1.0,
rmsErr( m_histo_pullz ) * 0.5 / std::sqrt( 1.0 + m_histo_pullz.standard_deviation() ) )
<< endmsg;
- info() << " ============================================" << endmsg;
+ info() << " ============================================" << endmsg;
}
return Consumer::finalize(); // Must be called after all other actions
}
+
//=============================================================================
// Match vertices by distance
//=============================================================================
diff --git a/Tr/TrackUtils/CMakeLists.txt b/Tr/TrackUtils/CMakeLists.txt
index 2e523f5f873..bb5b59d5b91 100644
--- a/Tr/TrackUtils/CMakeLists.txt
+++ b/Tr/TrackUtils/CMakeLists.txt
@@ -26,6 +26,7 @@ gaudi_add_module(TrackUtils
src/TrackCompetition.cpp
src/TrackContainerCopy.cpp
src/TrackContainerSplitter.cpp
+ src/RandomVeloTrackContainerSplitter.cpp
src/RandomTrackContainerSplitter.cpp
src/TrackContainersMerger.cpp
src/TrackListFilter.cpp
diff --git a/Tr/TrackUtils/src/RandomTrackContainerSplitter.cpp b/Tr/TrackUtils/src/RandomTrackContainerSplitter.cpp
index 248e4ab35ed..f9a9ac691f4 100644
--- a/Tr/TrackUtils/src/RandomTrackContainerSplitter.cpp
+++ b/Tr/TrackUtils/src/RandomTrackContainerSplitter.cpp
@@ -38,22 +38,6 @@ namespace {
};
} // namespace
-// bool my_splitter(auto trk){
-// Rndm::Numbers m_rnd;
-// IRndmGenSvc* rndm = svc<IRndmGenSvc>( "RndmGenSvc", true );
-// StatusCode sc = m_rnd.initialize( rndm, Rndm::Flat( 0.0, 100 ) );
-// if ( sc.isFailure() )
-// throw GaudiException{"Unable to initialize Flat distribution", "addConditionDerivation", StatusCode::FAILURE};
-// if ( m_rnd < 0.5 )
-// {
-// return false;
-// }
-// else
-// {
-// return true;
-// }
-// }
-
bool makeRandomDecision( const LHCb::Track& track, ServiceHandle<IRndmGenSvc>& rndmSvc ) {
if ( !rndmSvc.retrieve().isSuccess() ) {
std::cerr << "Error: Failed to retrieve IRndmGenSvc service." << std::endl;
@@ -85,7 +69,6 @@ public:
ServiceHandle<IRndmGenSvc> rndmSvc( "RndmGenSvc", "RandomTrackContainerSplitter" );
auto& [passed, rejected] = out;
m_inputs += input_tracks.size();
- // auto const& pred = my_splitter<TrackPredicate>();
for ( auto& trk : input_tracks ) {
auto& container = ( makeRandomDecision( *trk, rndmSvc ) ? passed : rejected );
container.insert( new LHCb::Track( *trk ) );
@@ -100,4 +83,4 @@ private:
mutable Gaudi::Accumulators::StatCounter<> m_passed{this, "#passed"};
};
-DECLARE_COMPONENT( RandomTrackContainerSplitter )
+DECLARE_COMPONENT( RandomTrackContainerSplitter )
\ No newline at end of file
diff --git a/Tr/TrackUtils/src/RandomVeloTrackContainerSplitter.cpp b/Tr/TrackUtils/src/RandomVeloTrackContainerSplitter.cpp
new file mode 100644
index 00000000000..8ec152d1fba
--- /dev/null
+++ b/Tr/TrackUtils/src/RandomVeloTrackContainerSplitter.cpp
@@ -0,0 +1,93 @@
+/*****************************************************************************\
+* (c) Copyright 2000-2020 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. *
+\*****************************************************************************/
+
+/** @class RandomVeloTrackContainerSplitter RandomVeloTrackContainerSplitter.h
+ *
+ * Split Velo container to two (passed and rejected) wrt to given selection criteria
+ *
+ * @author Andrii Usachov
+ * @date 30/04/2020
+ * @author Bogdan Kutsenko
+ * @date 25/05/2024
+ */
+
+#include "Event/Track.h"
+#include "Functors/Function.h"
+#include "Functors/with_functors.h"
+#include "GaudiKernel/IEventTimeDecoder.h"
+#include "GaudiKernel/IRndmEngine.h"
+#include "GaudiKernel/IRndmGenSvc.h"
+#include "GaudiKernel/RndmGenerators.h"
+#include "GaudiKernel/Service.h" // Include the header file for Service
+#include "GaudiKernel/ServiceHandle.h"
+#include "LHCbAlgs/Transformer.h"
+
+namespace {
+ /// The output data
+ using OutConts = std::tuple<LHCb::Pr::Velo::Tracks, LHCb::Pr::Velo::Tracks>;
+
+ struct TrackPredicate {
+ using Signature = bool( const LHCb::Pr::Velo::Tracks& );
+ static constexpr auto PropertyName = "Code";
+ };
+} // namespace
+
+bool makeRandomDecision( const LHCb::Pr::Velo::Tracks&, ServiceHandle<IRndmGenSvc>& rndmSvc ) {
+ if ( !rndmSvc.retrieve().isSuccess() ) {
+ std::cerr << "Error: Failed to retrieve IRndmGenSvc service." << std::endl;
+ return false;
+ }
+
+ IRndmGenSvc* rndmGenSvcPtr = rndmSvc.operator->();
+
+ if ( !rndmGenSvcPtr ) {
+ std::cerr << "Error: IRndmGenSvc pointer is null." << std::endl;
+ return false;
+ }
+
+ Rndm::Numbers rndm( rndmGenSvcPtr, Rndm::Flat( 0., 1.0 ) );
+ double randomNumber = rndm();
+
+ return randomNumber < 0.5;
+}
+
+class RandomVeloTrackContainerSplitter final
+ : public with_functors<LHCb::Algorithm::MultiTransformer<OutConts( const LHCb::Pr::Velo::Tracks& )>, TrackPredicate> {
+public:
+ RandomVeloTrackContainerSplitter( const std::string& name, ISvcLocator* pSvcLocator )
+ : with_functors( name, pSvcLocator, {KeyValue{"TracksInContainer", ""}},
+ {KeyValue{"PassedContainer", ""}, KeyValue{"RejectedContainer", ""}} ) {}
+
+ OutConts operator()( const LHCb::Pr::Velo::Tracks& input_tracks ) const override {
+ OutConts out;
+ ServiceHandle<IRndmGenSvc> rndmSvc( "RndmGenSvc", "RandomVeloTrackContainerSplitter" );
+ auto& [passed, rejected] = out;
+ m_inputs += input_tracks.size();
+
+ for ( auto const& trk : input_tracks.scalar() ) {
+ const int t = trk.offset();
+ LHCb::Pr::Velo::Tracks tmpTrk;
+ tmpTrk.reserve( 1 );
+ tmpTrk.copy_back<SIMDWrapper::scalar::types>(input_tracks, t, true);
+ auto& container = ( makeRandomDecision( tmpTrk, rndmSvc ) ? passed : rejected );
+ container.copy_back<SIMDWrapper::scalar::types>(input_tracks, t, true);
+ }
+ m_passed += passed.size();
+
+ return out;
+ }
+
+private:
+ mutable Gaudi::Accumulators::StatCounter<> m_inputs{this, "#inputs"};
+ mutable Gaudi::Accumulators::StatCounter<> m_passed{this, "#passed"};
+};
+
+DECLARE_COMPONENT( RandomVeloTrackContainerSplitter )
--
GitLab
From 980eb3668206a55ef0872cc38d7108341cba9b9f Mon Sep 17 00:00:00 2001
From: Gitlab CI <noreply@cern.ch>
Date: Sat, 25 May 2024 07:47:57 +0000
Subject: [PATCH 12/21] Fixed formatting
patch generated by https://gitlab.cern.ch/lhcb/Rec/-/jobs/39299105
---
Tr/TrackMonitors/src/VertexCompare.cpp | 412 +++++++++---------
.../src/RandomVeloTrackContainerSplitter.cpp | 9 +-
2 files changed, 204 insertions(+), 217 deletions(-)
diff --git a/Tr/TrackMonitors/src/VertexCompare.cpp b/Tr/TrackMonitors/src/VertexCompare.cpp
index 360d24fc961..d432c6679e4 100644
--- a/Tr/TrackMonitors/src/VertexCompare.cpp
+++ b/Tr/TrackMonitors/src/VertexCompare.cpp
@@ -89,31 +89,31 @@ bool trackcomp( const RecPVInfo<VertexType>& first, const RecPVInfo<VertexType>&
}
struct VtxVariables {
- int ntracks = -9999;
- double x = -9999.;
- double y = -9999.;
- double z = -9999.;
- double dxr = -9999.;
- double dyr = -9999.;
- double r = -9999.;
- double errx = -9999.;
- double erry = -9999.;
- double errz = -9999.;
- double err_r = -9999.;
- double covxx = -9999.;
- double covyy = -9999.;
- double covzz = -9999.;
- double covxy = -9999.;
- double covxz = -9999.;
- double covyz = -9999.;
- double chi2 = -9999.;
- double nDoF = -9999.;
- int dsize = -9999;
- int match = -9999;
- int pv_rank = -9999;
- bool equal_sizes = false;
- bool single = false;
- bool opposite_container = false;
+ int ntracks = -9999;
+ double x = -9999.;
+ double y = -9999.;
+ double z = -9999.;
+ double dxr = -9999.;
+ double dyr = -9999.;
+ double r = -9999.;
+ double errx = -9999.;
+ double erry = -9999.;
+ double errz = -9999.;
+ double err_r = -9999.;
+ double covxx = -9999.;
+ double covyy = -9999.;
+ double covzz = -9999.;
+ double covxy = -9999.;
+ double covxz = -9999.;
+ double covyz = -9999.;
+ double chi2 = -9999.;
+ double nDoF = -9999.;
+ int dsize = -9999;
+ int match = -9999;
+ int pv_rank = -9999;
+ bool equal_sizes = false;
+ bool single = false;
+ bool opposite_container = false;
std::vector<double> dx_vector{}, dy_vector{}, dz_vector{}, pullx_vector{}, pully_vector{}, pullz_vector{};
};
@@ -177,7 +177,7 @@ VtxVariables SetVtxVariables( const RecPVInfo<VertexType>& vrtf, const int& coun
// ============================================================================
template <typename Histogram>
-double meanErr( Histogram const & h ) {
+double meanErr( Histogram const& h ) {
const auto n = h.nEntries();
return ( 0 >= n ? 0.0 : h.standard_deviation() / std::sqrt( n ) );
}
@@ -186,28 +186,28 @@ double meanErr( Histogram const & h ) {
// Fourth momentum calculation
// ============================================================================
template <typename Histogram>
-double fourthMomentum( Histogram const & h ) {
+double fourthMomentum( Histogram const& h ) {
const auto n = h.nEntries();
if ( 0 >= n ) { return 0.0; } // RETURN
// get the mean
const auto h_mean = h.mean();
// number of bins
- const auto& axis = h.axis();
- const int nBins = h.totNBins();
- double result{ 0 }, weight{ 0 };
+ const auto& axis = h.axis();
+ const int nBins = h.totNBins();
+ double result{0}, weight{0};
// loop over all bins
- double minVal = axis[0].minValue;
- double maxVal = axis[0].maxValue;
- double binWidth = (maxVal-minVal)/nBins;
+ double minVal = axis[0].minValue;
+ double maxVal = axis[0].maxValue;
+ double binWidth = ( maxVal - minVal ) / nBins;
std::vector<double> binCenters;
- for (int i = 0; i < nBins; ++i) {
- double center = minVal + (i + 0.5) * binWidth;
- binCenters.push_back(center);
- }
- for (int i = 0; i < nBins; ++i) {
- const auto yBin = h.binValue( i ); // bin weight
+ for ( int i = 0; i < nBins; ++i ) {
+ double center = minVal + ( i + 0.5 ) * binWidth;
+ binCenters.push_back( center );
+ }
+ for ( int i = 0; i < nBins; ++i ) {
+ const auto yBin = h.binValue( i ); // bin weight
weight += yBin;
result += yBin * std::pow( binCenters[i] - h_mean, 4 );
}
@@ -221,26 +221,25 @@ double fourthMomentum( Histogram const & h ) {
// ============================================================================
template <typename Histogram>
-double kurtosis( Histogram const & h ) {
+double kurtosis( Histogram const& h ) {
const auto mu4 = fourthMomentum( h );
- const auto s4 = std::pow( h.standard_deviation() , 4 );
+ const auto s4 = std::pow( h.standard_deviation(), 4 );
return ( std::fabs( s4 ) > 0 ? mu4 / s4 - 3.0 : 0.0 );
}
// ============================================================================
// get an error in the rms value
// ============================================================================
-
- template <typename Histogram>
- double rmsErr( Histogram const & h ) {
- const auto n = h.nEntries();
- if ( 1 >= n ) { return 0.0; }
- auto result = 2.0 + kurtosis( h );
- result += 2.0 / ( n - 1 );
- result /= 4.0 * n;
- return h.standard_deviation() * std::sqrt( std::max( result, 0.0 ) );
- }
+template <typename Histogram>
+double rmsErr( Histogram const& h ) {
+ const auto n = h.nEntries();
+ if ( 1 >= n ) { return 0.0; }
+ auto result = 2.0 + kurtosis( h );
+ result += 2.0 / ( n - 1 );
+ result /= 4.0 * n;
+ return h.standard_deviation() * std::sqrt( std::max( result, 0.0 ) );
+}
class VertexCompare : public LHCb::Algorithm::Consumer<
void( Vertices const&, Vertices const&, LHCb::Conditions::InteractionRegion const& ),
@@ -263,7 +262,7 @@ public:
private:
bool debugLevel() const { return msgLevel( MSG::DEBUG ) || msgLevel( MSG::VERBOSE ); }
- static constexpr auto ntrack_bins = std::array{ 3.5 , 16.5 , 29.5 , 43.5 , 56.5 , 69.5 };
+ static constexpr auto ntrack_bins = std::array{3.5, 16.5, 29.5, 43.5, 56.5, 69.5};
Gaudi::Property<bool> m_produceHistogram{this, "produceHistogram", true};
Gaudi::Property<bool> m_produceNtuple{this, "produceNtuple", true};
@@ -271,7 +270,7 @@ private:
Gaudi::Property<bool> m_requireSingle{this, "requireSingle", false};
mutable std::vector<double> dx_vector{}, dy_vector{}, dz_vector{}, pullx_vector{}, pully_vector{}, pullz_vector{};
-
+
struct monitoringHistos {
mutable std::array<Gaudi::Accumulators::Histogram<1>, 5> m_histo_nTracksBins_dx;
mutable std::array<Gaudi::Accumulators::Histogram<1>, 5> m_histo_nTracksBins_dy;
@@ -280,7 +279,6 @@ private:
mutable Gaudi::Accumulators::Histogram<1> m_histo_pully_Monitoring;
mutable Gaudi::Accumulators::Histogram<1> m_histo_pullz_Monitoring;
-
template <std::size_t... IDXs>
static std::array<Gaudi::Accumulators::Histogram<1>, sizeof...( IDXs )>
histo1DArrayBuilder( const VertexCompare* owner, const std::string& name, const std::string& title,
@@ -290,32 +288,35 @@ private:
title + std::to_string( IDXs ),
{std::get<0>( xbins ), std::get<1>( xbins ), std::get<2>( xbins )}}...}};
}
- Gaudi::Histo1DDef m_histodx{"dx, mm", -0.15, 0.15, 50};
+ Gaudi::Histo1DDef m_histodx{"dx, mm", -0.15, 0.15, 50};
monitoringHistos( const VertexCompare* owner )
- : m_histo_nTracksBins_dx{histo1DArrayBuilder( owner, "dx_Monitoring_", "dx, mm", {50, -0.15, 0.15}, std::make_index_sequence<5>() )}
- , m_histo_nTracksBins_dy{histo1DArrayBuilder( owner, "dy_Monitoring_", "dy, mm", {50, -0.15, 0.15}, std::make_index_sequence<5>() )}
- , m_histo_nTracksBins_dz{histo1DArrayBuilder( owner, "dz_Monitoring_", "dz, mm", {50, -1.5, 1.5}, std::make_index_sequence<5>() )}
+ : m_histo_nTracksBins_dx{histo1DArrayBuilder( owner, "dx_Monitoring_", "dx, mm", {50, -0.15, 0.15},
+ std::make_index_sequence<5>() )}
+ , m_histo_nTracksBins_dy{histo1DArrayBuilder( owner, "dy_Monitoring_", "dy, mm", {50, -0.15, 0.15},
+ std::make_index_sequence<5>() )}
+ , m_histo_nTracksBins_dz{histo1DArrayBuilder( owner, "dz_Monitoring_", "dz, mm", {50, -1.5, 1.5},
+ std::make_index_sequence<5>() )}
, m_histo_pullx_Monitoring{owner, "pullx_Monitoring", "pull x", {20, -5, 5}}
, m_histo_pully_Monitoring{owner, "pully_Monitoring", "pull y", {20, -5, 5}}
, m_histo_pullz_Monitoring{owner, "pullz_Monitoring", "pull z", {20, -5, 5}} {}
-
-
-
};
std::unique_ptr<monitoringHistos> m_monitoringHistos;
using axis1D = Gaudi::Accumulators::Axis<Gaudi::Accumulators::Histogram<1>::AxisArithmeticType>;
- mutable Gaudi::Accumulators::RootHistogram<1> m_histo_dx{this, "dx", "dx, mm", axis1D{50, -0.15, 0.15} };
- mutable Gaudi::Accumulators::RootHistogram<1> m_histo_dy{ this, "dy", "dy, mm", axis1D{50, -0.15, 0.15} };
- mutable Gaudi::Accumulators::RootHistogram<1> m_histo_dz{ this, "dz", "dz, mm", axis1D{50, -1.5, 1.5} };
- mutable Gaudi::Accumulators::RootHistogram<1> m_histo_pullx{ this, "pullx", "pull x", axis1D{20, -5, 5} };
- mutable Gaudi::Accumulators::RootHistogram<1> m_histo_pully{ this, "pully", "pull y", axis1D{20, -5, 5} };
- mutable Gaudi::Accumulators::RootHistogram<1> m_histo_pullz{ this, "pullz", "pull z", axis1D{20, -5, 5} };
- mutable Gaudi::Accumulators::Histogram<1> m_nTracks1{ this, "ntracks1", "Number of tracks in PV set 1", axis1D{50, 0., 150.} };
- mutable Gaudi::Accumulators::Histogram<1> m_nTracks2{ this, "ntracks2", "Number of tracks in PV set 2", axis1D{ 50, 0., 150.}};
- mutable Gaudi::Accumulators::Histogram<1> m_nTracks_dif{ this, "dtracks", "Difference in the number of tracks in two sets of PVs", axis1D{ 50, 0., 150.} };
-
+ mutable Gaudi::Accumulators::RootHistogram<1> m_histo_dx{this, "dx", "dx, mm", axis1D{50, -0.15, 0.15}};
+ mutable Gaudi::Accumulators::RootHistogram<1> m_histo_dy{this, "dy", "dy, mm", axis1D{50, -0.15, 0.15}};
+ mutable Gaudi::Accumulators::RootHistogram<1> m_histo_dz{this, "dz", "dz, mm", axis1D{50, -1.5, 1.5}};
+ mutable Gaudi::Accumulators::RootHistogram<1> m_histo_pullx{this, "pullx", "pull x", axis1D{20, -5, 5}};
+ mutable Gaudi::Accumulators::RootHistogram<1> m_histo_pully{this, "pully", "pull y", axis1D{20, -5, 5}};
+ mutable Gaudi::Accumulators::RootHistogram<1> m_histo_pullz{this, "pullz", "pull z", axis1D{20, -5, 5}};
+ mutable Gaudi::Accumulators::Histogram<1> m_nTracks1{this, "ntracks1", "Number of tracks in PV set 1",
+ axis1D{50, 0., 150.}};
+ mutable Gaudi::Accumulators::Histogram<1> m_nTracks2{this, "ntracks2", "Number of tracks in PV set 2",
+ axis1D{50, 0., 150.}};
+ mutable Gaudi::Accumulators::Histogram<1> m_nTracks_dif{
+ this, "dtracks", "Difference in the number of tracks in two sets of PVs", axis1D{50, 0., 150.}};
+
mutable Gaudi::Accumulators::Counter<> m_nVtx{this, "Number of pairs of vertices in processed events"};
mutable Gaudi::Accumulators::Counter<> m_nRec1{this, "Number of vertices in input1"};
mutable Gaudi::Accumulators::Counter<> m_nRec2{this, "Number of vertices in input2"};
@@ -355,96 +356,96 @@ StatusCode VertexCompare::initialize() {
//=============================================================================
void VertexCompare::operator()( Vertices const& recoVtx1, Vertices const& recoVtx2,
LHCb::Conditions::InteractionRegion const& region ) const {
- if ( msgLevel( MSG::DEBUG ) ) debug() << "==> Execute" << endmsg;
+ if ( msgLevel( MSG::DEBUG ) ) debug() << "==> Execute" << endmsg;
const auto beamspot = region.avgPosition;
std::vector<RecPVInfo<VertexType>> fullVrt1, fullVrt2;
fullVrt1.reserve( recoVtx1.size() );
fullVrt2.reserve( recoVtx2.size() );
- if(m_requireSingle == true && !(recoVtx1.size() == 1 && recoVtx2.size() == 1)) return;
-
+ if ( m_requireSingle == true && !( recoVtx1.size() == 1 && recoVtx2.size() == 1 ) ) return;
+
if ( debugLevel() ) debug() << " Vtx Properities x y z chi2/ndof ntracks" << endmsg;
- int size_diff = -9999;
+ int size_diff = -9999;
RecPVInfo<VertexType> recinfo1;
RecPVInfo<VertexType> recinfo2;
- for ( auto const& pv : recoVtx1 ) {
- ++m_nRec1;
- recinfo1.pRECPV = &pv;
- recinfo1.position = pv.position();
- recinfo1.covPV = pv.covMatrix();
- recinfo1.positionSigma = Gaudi::XYZPoint{std::sqrt( recinfo1.covPV( 0, 0 ) ), std::sqrt( recinfo1.covPV( 1, 1 ) ),
- std::sqrt( recinfo1.covPV( 2, 2 ) )};
- recinfo1.ntracks = pv.nTracks();
- if ( recinfo1.ntracks >= 10 ) { ++m_nRec_geq10tr_1; }
- recinfo1.chi2 = pv.chi2();
- recinfo1.nDoF = pv.nDoF();
-
- if ( debugLevel() )
- debug() << " " << pv.position().x() << " " << pv.position().y() << " " << pv.position().z()
- << " " << pv.chi2PerDoF() << " " << pv.nTracks() << endmsg;
- fullVrt1.push_back( recinfo1 );
- } // end of loop over vertices1
-
- for ( auto const& pv : recoVtx2 ) {
- ++m_nRec2;
- recinfo2.pRECPV = &pv;
- recinfo2.position = pv.position();
- recinfo2.covPV = pv.covMatrix();
- recinfo2.positionSigma = Gaudi::XYZPoint{std::sqrt( recinfo2.covPV( 0, 0 ) ), std::sqrt( recinfo2.covPV( 1, 1 ) ),
- std::sqrt( recinfo2.covPV( 2, 2 ) )};
- recinfo2.ntracks = pv.nTracks();
- if ( recinfo2.ntracks >= 10 ) { ++m_nRec_geq10tr_2; }
- recinfo2.chi2 = pv.chi2();
- recinfo2.nDoF = pv.nDoF();
-
- if ( debugLevel() )
- debug() << " " << pv.position().x() << " " << pv.position().y() << " " << pv.position().z()
- << " " << pv.chi2PerDoF() << " " << pv.nTracks() << endmsg;
- fullVrt2.push_back( recinfo2 );
- } // end of loop over vertices2
-
- // added sorting to mirror the 'multirec' variable implementation in PVChecker
- std::sort( fullVrt1.begin(), fullVrt1.end(), trackcomp );
- std::sort( fullVrt2.begin(), fullVrt2.end(), trackcomp );
-
- if ( debugLevel() ) debug() << "fullVrt1 size " << fullVrt1.size() << endmsg;
- if ( debugLevel() ) debug() << "fullVrt2 size " << fullVrt2.size() << endmsg;
- size_diff = fullVrt1.size() - fullVrt2.size();
-
- if ( m_produceNtuple.value() ) {
- Tuple myTuple_evt = nTuple( "PV_nTuple_evt", "PV_nTuple_evt", CLID_ColumnWiseTuple );
- myTuple_evt->column( "size_diff", double( size_diff ) ).ignore();
- myTuple_evt->column( "size_1", double( fullVrt1.size() ) ).ignore();
- myTuple_evt->column( "size_2", double( fullVrt2.size() ) ).ignore();
- myTuple_evt->write().ignore();
- }
+ for ( auto const& pv : recoVtx1 ) {
+ ++m_nRec1;
+ recinfo1.pRECPV = &pv;
+ recinfo1.position = pv.position();
+ recinfo1.covPV = pv.covMatrix();
+ recinfo1.positionSigma = Gaudi::XYZPoint{std::sqrt( recinfo1.covPV( 0, 0 ) ), std::sqrt( recinfo1.covPV( 1, 1 ) ),
+ std::sqrt( recinfo1.covPV( 2, 2 ) )};
+ recinfo1.ntracks = pv.nTracks();
+ if ( recinfo1.ntracks >= 10 ) { ++m_nRec_geq10tr_1; }
+ recinfo1.chi2 = pv.chi2();
+ recinfo1.nDoF = pv.nDoF();
+
+ if ( debugLevel() )
+ debug() << " " << pv.position().x() << " " << pv.position().y() << " " << pv.position().z()
+ << " " << pv.chi2PerDoF() << " " << pv.nTracks() << endmsg;
+ fullVrt1.push_back( recinfo1 );
+ } // end of loop over vertices1
+
+ for ( auto const& pv : recoVtx2 ) {
+ ++m_nRec2;
+ recinfo2.pRECPV = &pv;
+ recinfo2.position = pv.position();
+ recinfo2.covPV = pv.covMatrix();
+ recinfo2.positionSigma = Gaudi::XYZPoint{std::sqrt( recinfo2.covPV( 0, 0 ) ), std::sqrt( recinfo2.covPV( 1, 1 ) ),
+ std::sqrt( recinfo2.covPV( 2, 2 ) )};
+ recinfo2.ntracks = pv.nTracks();
+ if ( recinfo2.ntracks >= 10 ) { ++m_nRec_geq10tr_2; }
+ recinfo2.chi2 = pv.chi2();
+ recinfo2.nDoF = pv.nDoF();
+
+ if ( debugLevel() )
+ debug() << " " << pv.position().x() << " " << pv.position().y() << " " << pv.position().z()
+ << " " << pv.chi2PerDoF() << " " << pv.nTracks() << endmsg;
+ fullVrt2.push_back( recinfo2 );
+ } // end of loop over vertices2
+
+ // added sorting to mirror the 'multirec' variable implementation in PVChecker
+ std::sort( fullVrt1.begin(), fullVrt1.end(), trackcomp );
+ std::sort( fullVrt2.begin(), fullVrt2.end(), trackcomp );
+
+ if ( debugLevel() ) debug() << "fullVrt1 size " << fullVrt1.size() << endmsg;
+ if ( debugLevel() ) debug() << "fullVrt2 size " << fullVrt2.size() << endmsg;
+ size_diff = fullVrt1.size() - fullVrt2.size();
+
+ if ( m_produceNtuple.value() ) {
+ Tuple myTuple_evt = nTuple( "PV_nTuple_evt", "PV_nTuple_evt", CLID_ColumnWiseTuple );
+ myTuple_evt->column( "size_diff", double( size_diff ) ).ignore();
+ myTuple_evt->column( "size_1", double( fullVrt1.size() ) ).ignore();
+ myTuple_evt->column( "size_2", double( fullVrt2.size() ) ).ignore();
+ myTuple_evt->write().ignore();
+ }
std::vector<int> link;
- int ntracks1 = 0;
- int ntracks2 = 0;
- int dtracks = 0;
- int size1 = 0;
- int size2 = 0;
- double sigx_part1 = -99999.;
- double sigy_part1 = -99999.;
- double sigz_part1 = -99999.;
- double sigx_part2 = -99999.;
- double sigy_part2 = -99999.;
- double sigz_part2 = -99999.;
- double x1 = -99999.;
- double y1 = -99999.;
- double z1 = -99999.;
- double x2 = -99999.;
- double y2 = -99999.;
- double z2 = -99999.;
- double dx = -99999.;
- double dy = -99999.;
- double dz = -99999.;
- int oIt = 0;
- int pv_rank = -99999;
+ int ntracks1 = 0;
+ int ntracks2 = 0;
+ int dtracks = 0;
+ int size1 = 0;
+ int size2 = 0;
+ double sigx_part1 = -99999.;
+ double sigy_part1 = -99999.;
+ double sigz_part1 = -99999.;
+ double sigx_part2 = -99999.;
+ double sigy_part2 = -99999.;
+ double sigz_part2 = -99999.;
+ double x1 = -99999.;
+ double y1 = -99999.;
+ double z1 = -99999.;
+ double x2 = -99999.;
+ double y2 = -99999.;
+ double z2 = -99999.;
+ double dx = -99999.;
+ double dy = -99999.;
+ double dz = -99999.;
+ int oIt = 0;
+ int pv_rank = -99999;
if ( fullVrt1.size() != 0 && fullVrt2.size() != 0 ) { matchByDistance( fullVrt1, fullVrt2, link ); }
@@ -499,10 +500,10 @@ void VertexCompare::operator()( Vertices const& recoVtx1, Vertices const& recoVt
}
if ( m_monitoring.value() ) {
- if ( dz < 2 && size_diff < 3 && pv_rank < 2 && ( ntracks1 + ntracks2 ) / 2 > ntrack_bins[0]) {
+ if ( dz < 2 && size_diff < 3 && pv_rank < 2 && ( ntracks1 + ntracks2 ) / 2 > ntrack_bins[0] ) {
int binCount = 0;
for ( size_t i = 1; i < ntrack_bins.size(); ++i ) {
- if ( std::lround(( ntracks1 + ntracks2 ) / 2) <= ntrack_bins[i] ) { break; }
+ if ( std::lround( ( ntracks1 + ntracks2 ) / 2 ) <= ntrack_bins[i] ) { break; }
binCount++;
}
auto& monitoringHistos = *m_monitoringHistos.get();
@@ -512,15 +513,15 @@ void VertexCompare::operator()( Vertices const& recoVtx1, Vertices const& recoVt
++monitoringHistos.m_histo_pullx_Monitoring[dx / errx];
++monitoringHistos.m_histo_pully_Monitoring[dy / erry];
++monitoringHistos.m_histo_pullz_Monitoring[dz / errz];
- }
+ }
}
- if (m_produceHistogram.value()) {
- dx_vector.push_back(dx);
- dy_vector.push_back(dy);
- dz_vector.push_back(dz);
- pullx_vector.push_back(dx / errx);
- pully_vector.push_back(dy / erry);
- pullz_vector.push_back(dz / errz);
+ if ( m_produceHistogram.value() ) {
+ dx_vector.push_back( dx );
+ dy_vector.push_back( dy );
+ dz_vector.push_back( dz );
+ pullx_vector.push_back( dx / errx );
+ pully_vector.push_back( dy / erry );
+ pullz_vector.push_back( dz / errz );
++m_histo_dx[dx];
++m_histo_dy[dy];
++m_histo_dz[dz];
@@ -529,8 +530,8 @@ void VertexCompare::operator()( Vertices const& recoVtx1, Vertices const& recoVt
++m_histo_pullz[dz / errz];
++m_nTracks1[ntracks1];
++m_nTracks2[ntracks2];
- ++m_nTracks_dif[ntracks2-ntracks1];
- } //else end
+ ++m_nTracks_dif[ntracks2 - ntracks1];
+ } // else end
if ( m_produceNtuple.value() ) {
Tuple myTuple = nTuple( "PV_nTuple", "PV_nTuple", CLID_ColumnWiseTuple );
@@ -696,7 +697,7 @@ void VertexCompare::operator()( Vertices const& recoVtx1, Vertices const& recoVt
oIt++;
m_nEvt += 1;
}
- }
+}
//=============================================================================
// Finalize
@@ -733,85 +734,70 @@ StatusCode VertexCompare::finalize() {
<< ( (double)m_nVtx.nEntries() - (double)m_nVtx_geq10tr_matched2.nEntries() ) /
( (double)m_nRec2.nEntries() - (double)m_nRec_geq10tr_2.nEntries() ) * 100
<< " % " << endmsg;
-
- if(m_produceHistogram.value()){
- for (double value : dx_vector) {
- ++m_histo_dx[value];
- }
- for (double value : dy_vector) {
- ++m_histo_dy[value];
- }
- for (double value : dz_vector) {
- ++m_histo_dz[value];
- }
- for (double value : pullx_vector) {
- ++m_histo_pullx[value];
- }
- for (double value : pully_vector) {
- ++m_histo_pully[value];
- }
- for (double value : pullz_vector) {
- ++m_histo_pullz[value];
- }
+ if ( m_produceHistogram.value() ) {
+
+ for ( double value : dx_vector ) { ++m_histo_dx[value]; }
+ for ( double value : dy_vector ) { ++m_histo_dy[value]; }
+ for ( double value : dz_vector ) { ++m_histo_dz[value]; }
+ for ( double value : pullx_vector ) { ++m_histo_pullx[value]; }
+ for ( double value : pully_vector ) { ++m_histo_pully[value]; }
+ for ( double value : pullz_vector ) { ++m_histo_pullz[value]; }
- info() << "dx: "
- << format( "mean = %5.3f +/- %5.3f, RMS = %5.3f +/- %5.3f", m_histo_dx.mean(),
- meanErr( m_histo_dx ), m_histo_dx.standard_deviation(), rmsErr( m_histo_dx) )
+ info() << "dx: "
+ << format( "mean = %5.3f +/- %5.3f, RMS = %5.3f +/- %5.3f", m_histo_dx.mean(), meanErr( m_histo_dx ),
+ m_histo_dx.standard_deviation(), rmsErr( m_histo_dx ) )
<< endmsg;
-
-
- info() << "dy: "
- << format( "mean = %5.3f +/- %5.3f, RMS = %5.3f +/- %5.3f", m_histo_dy.mean(),
- meanErr( m_histo_dy ), m_histo_dy.standard_deviation(), rmsErr( m_histo_dy) )
+
+ info() << "dy: "
+ << format( "mean = %5.3f +/- %5.3f, RMS = %5.3f +/- %5.3f", m_histo_dy.mean(), meanErr( m_histo_dy ),
+ m_histo_dy.standard_deviation(), rmsErr( m_histo_dy ) )
<< endmsg;
-
- info() << "dz: "
- << format( "mean = %5.3f +/- %5.3f, RMS = %5.3f +/- %5.3f", m_histo_dz.mean(),
- meanErr( m_histo_dz ), m_histo_dz.standard_deviation(), rmsErr( m_histo_dz) )
+ info() << "dz: "
+ << format( "mean = %5.3f +/- %5.3f, RMS = %5.3f +/- %5.3f", m_histo_dz.mean(), meanErr( m_histo_dz ),
+ m_histo_dz.standard_deviation(), rmsErr( m_histo_dz ) )
<< endmsg;
-
- info() << " ---------------------------------------" << endmsg;
- info() << "pullx: "
- << format( "mean = %5.3f +/- %5.3f, RMS = %5.3f +/- %5.3f", m_histo_pullx.mean(),
- meanErr( m_histo_pullx ), m_histo_pullx.standard_deviation(), rmsErr( m_histo_pullx ) )
+ info() << " ---------------------------------------" << endmsg;
+
+ info() << "pullx: "
+ << format( "mean = %5.3f +/- %5.3f, RMS = %5.3f +/- %5.3f", m_histo_pullx.mean(), meanErr( m_histo_pullx ),
+ m_histo_pullx.standard_deviation(), rmsErr( m_histo_pullx ) )
<< endmsg;
-
- info() << "pully: "
- << format( "mean = %5.3f +/- %5.3f, RMS = %5.3f +/- %5.3f", m_histo_pully.mean(),
- meanErr( m_histo_pully ), m_histo_pully.standard_deviation(), rmsErr( m_histo_pully) )
+
+ info() << "pully: "
+ << format( "mean = %5.3f +/- %5.3f, RMS = %5.3f +/- %5.3f", m_histo_pully.mean(), meanErr( m_histo_pully ),
+ m_histo_pully.standard_deviation(), rmsErr( m_histo_pully ) )
<< endmsg;
-
- info() << "pullz: "
- << format( "mean = %5.3f +/- %5.3f, RMS = %5.3f +/- %5.3f", m_histo_pullz.mean(),
- meanErr( m_histo_pullz ), m_histo_pullz.standard_deviation(), rmsErr( m_histo_pullz) )
+
+ info() << "pullz: "
+ << format( "mean = %5.3f +/- %5.3f, RMS = %5.3f +/- %5.3f", m_histo_pullz.mean(), meanErr( m_histo_pullz ),
+ m_histo_pullz.standard_deviation(), rmsErr( m_histo_pullz ) )
<< endmsg;
-
- info() << " ---------------------------------------" << endmsg;
- info() << "diff in x: "
+ info() << " ---------------------------------------" << endmsg;
+
+ info() << "diff in x: "
<< format( "%5.3f +/- %5.3f", std::sqrt( 1.0 + m_histo_pullx.standard_deviation() ) - 1.0,
rmsErr( m_histo_pullx ) * 0.5 / std::sqrt( 1.0 + m_histo_pullx.standard_deviation() ) )
<< endmsg;
- info() << "diff in y: "
+ info() << "diff in y: "
<< format( "%5.3f +/- %5.3f", std::sqrt( 1.0 + m_histo_pully.standard_deviation() ) - 1.0,
rmsErr( m_histo_pully ) * 0.5 / std::sqrt( 1.0 + m_histo_pully.standard_deviation() ) )
<< endmsg;
- info() << "diff in z: "
+ info() << "diff in z: "
<< format( "%5.3f +/- %5.3f", std::sqrt( 1.0 + m_histo_pullz.standard_deviation() ) - 1.0,
rmsErr( m_histo_pullz ) * 0.5 / std::sqrt( 1.0 + m_histo_pullz.standard_deviation() ) )
<< endmsg;
- info() << " ============================================" << endmsg;
+ info() << " ============================================" << endmsg;
}
return Consumer::finalize(); // Must be called after all other actions
}
-
//=============================================================================
// Match vertices by distance
//=============================================================================
diff --git a/Tr/TrackUtils/src/RandomVeloTrackContainerSplitter.cpp b/Tr/TrackUtils/src/RandomVeloTrackContainerSplitter.cpp
index 8ec152d1fba..7188900e766 100644
--- a/Tr/TrackUtils/src/RandomVeloTrackContainerSplitter.cpp
+++ b/Tr/TrackUtils/src/RandomVeloTrackContainerSplitter.cpp
@@ -60,7 +60,8 @@ bool makeRandomDecision( const LHCb::Pr::Velo::Tracks&, ServiceHandle<IRndmGenSv
}
class RandomVeloTrackContainerSplitter final
- : public with_functors<LHCb::Algorithm::MultiTransformer<OutConts( const LHCb::Pr::Velo::Tracks& )>, TrackPredicate> {
+ : public with_functors<LHCb::Algorithm::MultiTransformer<OutConts( const LHCb::Pr::Velo::Tracks& )>,
+ TrackPredicate> {
public:
RandomVeloTrackContainerSplitter( const std::string& name, ISvcLocator* pSvcLocator )
: with_functors( name, pSvcLocator, {KeyValue{"TracksInContainer", ""}},
@@ -73,12 +74,12 @@ public:
m_inputs += input_tracks.size();
for ( auto const& trk : input_tracks.scalar() ) {
- const int t = trk.offset();
+ const int t = trk.offset();
LHCb::Pr::Velo::Tracks tmpTrk;
tmpTrk.reserve( 1 );
- tmpTrk.copy_back<SIMDWrapper::scalar::types>(input_tracks, t, true);
+ tmpTrk.copy_back<SIMDWrapper::scalar::types>( input_tracks, t, true );
auto& container = ( makeRandomDecision( tmpTrk, rndmSvc ) ? passed : rejected );
- container.copy_back<SIMDWrapper::scalar::types>(input_tracks, t, true);
+ container.copy_back<SIMDWrapper::scalar::types>( input_tracks, t, true );
}
m_passed += passed.size();
--
GitLab
From e2436e5df7a9d7904bb7855b5105958a0a7827b2 Mon Sep 17 00:00:00 2001
From: Bogdan Kutsenko <bogdan.kutsenko@cern.ch>
Date: Sat, 25 May 2024 11:17:16 +0200
Subject: [PATCH 13/21] Histo stats calculation simplified
---
Tr/TrackMonitors/src/VertexCompare.cpp | 15 ++++++---------
1 file changed, 6 insertions(+), 9 deletions(-)
diff --git a/Tr/TrackMonitors/src/VertexCompare.cpp b/Tr/TrackMonitors/src/VertexCompare.cpp
index d432c6679e4..72215c7f2bb 100644
--- a/Tr/TrackMonitors/src/VertexCompare.cpp
+++ b/Tr/TrackMonitors/src/VertexCompare.cpp
@@ -186,7 +186,7 @@ double meanErr( Histogram const& h ) {
// Fourth momentum calculation
// ============================================================================
template <typename Histogram>
-double fourthMomentum( Histogram const& h ) {
+double fourthMoment( Histogram const& h ) {
const auto n = h.nEntries();
if ( 0 >= n ) { return 0.0; } // RETURN
// get the mean
@@ -201,16 +201,13 @@ double fourthMomentum( Histogram const& h ) {
double maxVal = axis[0].maxValue;
double binWidth = ( maxVal - minVal ) / nBins;
- std::vector<double> binCenters;
- for ( int i = 0; i < nBins; ++i ) {
- double center = minVal + ( i + 0.5 ) * binWidth;
- binCenters.push_back( center );
- }
- for ( int i = 0; i < nBins; ++i ) {
+ for ( int i = 0; i < nBins; ++i ) {
const auto yBin = h.binValue( i ); // bin weight
weight += yBin;
- result += yBin * std::pow( binCenters[i] - h_mean, 4 );
+ double center = minVal + ( i + 0.5 ) * binWidth;
+ result += yBin * std::pow( center - h_mean, 4 );
}
+
if ( weight > std::numeric_limits<double>::epsilon() ) { result /= weight; }
return result;
@@ -222,7 +219,7 @@ double fourthMomentum( Histogram const& h ) {
template <typename Histogram>
double kurtosis( Histogram const& h ) {
- const auto mu4 = fourthMomentum( h );
+ const auto mu4 = fourthMoment( h );
const auto s4 = std::pow( h.standard_deviation(), 4 );
return ( std::fabs( s4 ) > 0 ? mu4 / s4 - 3.0 : 0.0 );
}
--
GitLab
From 6636bb3a49a1070cd602b5e68bc7bf63ff12cd70 Mon Sep 17 00:00:00 2001
From: Gitlab CI <noreply@cern.ch>
Date: Sat, 25 May 2024 09:18:07 +0000
Subject: [PATCH 14/21] Fixed formatting
patch generated by https://gitlab.cern.ch/lhcb/Rec/-/jobs/39300421
---
Tr/TrackMonitors/src/VertexCompare.cpp | 2 +-
1 file changed, 1 insertion(+), 1 deletion(-)
diff --git a/Tr/TrackMonitors/src/VertexCompare.cpp b/Tr/TrackMonitors/src/VertexCompare.cpp
index 72215c7f2bb..0657c7d314e 100644
--- a/Tr/TrackMonitors/src/VertexCompare.cpp
+++ b/Tr/TrackMonitors/src/VertexCompare.cpp
@@ -201,7 +201,7 @@ double fourthMoment( Histogram const& h ) {
double maxVal = axis[0].maxValue;
double binWidth = ( maxVal - minVal ) / nBins;
- for ( int i = 0; i < nBins; ++i ) {
+ for ( int i = 0; i < nBins; ++i ) {
const auto yBin = h.binValue( i ); // bin weight
weight += yBin;
double center = minVal + ( i + 0.5 ) * binWidth;
--
GitLab
From 2474b06ef4c25c69f92606111fb4f9b78bd361bd Mon Sep 17 00:00:00 2001
From: Bogdan Kutsenko <bogdan.kutsenko@cern.ch>
Date: Mon, 27 May 2024 18:29:23 +0200
Subject: [PATCH 15/21] Track Splitter and vertex compare - optimized
---
Tr/TrackMonitors/src/VertexCompare.cpp | 178 +++++++++---------
.../src/RandomTrackContainerSplitter.cpp | 39 +---
.../src/RandomVeloTrackContainerSplitter.cpp | 34 +---
3 files changed, 99 insertions(+), 152 deletions(-)
diff --git a/Tr/TrackMonitors/src/VertexCompare.cpp b/Tr/TrackMonitors/src/VertexCompare.cpp
index 0657c7d314e..0f93b548b61 100644
--- a/Tr/TrackMonitors/src/VertexCompare.cpp
+++ b/Tr/TrackMonitors/src/VertexCompare.cpp
@@ -23,6 +23,7 @@
#include "GaudiAlg/Tuples.h"
#include "GaudiKernel/PhysicalConstants.h"
#include "GaudiKernel/SystemOfUnits.h"
+#include "GaudiKernel/HistoDef.h"
#include "GaudiUtils/HistoStats.h"
#include "LHCbAlgs/Consumer.h"
#include "MCInterfaces/IForcedBDecayTool.h"
@@ -114,7 +115,6 @@ struct VtxVariables {
bool equal_sizes = false;
bool single = false;
bool opposite_container = false;
- std::vector<double> dx_vector{}, dy_vector{}, dz_vector{}, pullx_vector{}, pully_vector{}, pullz_vector{};
};
VtxVariables SetVtxVariables( const RecPVInfo<VertexType>& vrtf, const int& counter, const int& size1, const int& size2,
@@ -172,55 +172,15 @@ VtxVariables SetVtxVariables( const RecPVInfo<VertexType>& vrtf, const int& coun
return vtx;
}
-// ============================================================================
-// get an error in the mean value
-// ============================================================================
-
-template <typename Histogram>
-double meanErr( Histogram const& h ) {
- const auto n = h.nEntries();
- return ( 0 >= n ? 0.0 : h.standard_deviation() / std::sqrt( n ) );
-}
-
-// ============================================================================
-// Fourth momentum calculation
-// ============================================================================
-template <typename Histogram>
-double fourthMoment( Histogram const& h ) {
- const auto n = h.nEntries();
- if ( 0 >= n ) { return 0.0; } // RETURN
- // get the mean
- const auto h_mean = h.mean();
- // number of bins
- const auto& axis = h.axis();
- const int nBins = h.totNBins();
- double result{0}, weight{0};
-
- // loop over all bins
- double minVal = axis[0].minValue;
- double maxVal = axis[0].maxValue;
- double binWidth = ( maxVal - minVal ) / nBins;
-
- for ( int i = 0; i < nBins; ++i ) {
- const auto yBin = h.binValue( i ); // bin weight
- weight += yBin;
- double center = minVal + ( i + 0.5 ) * binWidth;
- result += yBin * std::pow( center - h_mean, 4 );
- }
-
- if ( weight > std::numeric_limits<double>::epsilon() ) { result /= weight; }
-
- return result;
-}
-
// ============================================================================
// get the error in kurtosis
// ============================================================================
-template <typename Histogram>
-double kurtosis( Histogram const& h ) {
- const auto mu4 = fourthMoment( h );
- const auto s4 = std::pow( h.standard_deviation(), 4 );
+template <typename StatAccumulator>
+double kurtosis( StatAccumulator const& stAcc , StatAccumulator const& stAcc_sq ) {
+ const auto n = stAcc.nEntries();
+ const auto mu4 = stAcc_sq.sum2()/n; // Calculate a 4t moment using statAccumulater of squared variables
+ const auto s4 = std::pow( stAcc.standard_deviation(), 4 );
return ( std::fabs( s4 ) > 0 ? mu4 / s4 - 3.0 : 0.0 );
}
@@ -228,14 +188,14 @@ double kurtosis( Histogram const& h ) {
// get an error in the rms value
// ============================================================================
-template <typename Histogram>
-double rmsErr( Histogram const& h ) {
- const auto n = h.nEntries();
+template <typename StatAccumulator>
+double rmsErr( StatAccumulator const& stAcc, StatAccumulator const& stAcc_sq ) {
+ const auto n = stAcc.nEntries();
if ( 1 >= n ) { return 0.0; }
- auto result = 2.0 + kurtosis( h );
+ auto result = 2.0 + kurtosis( stAcc , stAcc_sq );
result += 2.0 / ( n - 1 );
result /= 4.0 * n;
- return h.standard_deviation() * std::sqrt( std::max( result, 0.0 ) );
+ return stAcc.standard_deviation() * std::sqrt( std::max( result, 0.0 ) );
}
class VertexCompare : public LHCb::Algorithm::Consumer<
@@ -266,8 +226,25 @@ private:
Gaudi::Property<bool> m_monitoring{this, "monitoring", false};
Gaudi::Property<bool> m_requireSingle{this, "requireSingle", false};
- mutable std::vector<double> dx_vector{}, dy_vector{}, dz_vector{}, pullx_vector{}, pully_vector{}, pullz_vector{};
+ mutable Gaudi::Accumulators::StatAccumulator<Gaudi::Accumulators::atomicity::full, double> m_stat_dx;
+ mutable Gaudi::Accumulators::StatAccumulator<Gaudi::Accumulators::atomicity::full, double> m_stat_dx_sq;
+
+ mutable Gaudi::Accumulators::StatAccumulator<Gaudi::Accumulators::atomicity::full, double> m_stat_dy;
+ mutable Gaudi::Accumulators::StatAccumulator<Gaudi::Accumulators::atomicity::full, double> m_stat_dy_sq;
+ mutable Gaudi::Accumulators::StatAccumulator<Gaudi::Accumulators::atomicity::full, double> m_stat_dz;
+ mutable Gaudi::Accumulators::StatAccumulator<Gaudi::Accumulators::atomicity::full, double> m_stat_dz_sq;
+
+ mutable Gaudi::Accumulators::StatAccumulator<Gaudi::Accumulators::atomicity::full, double> m_stat_pullx;
+ mutable Gaudi::Accumulators::StatAccumulator<Gaudi::Accumulators::atomicity::full, double> m_stat_pullx_sq;
+
+ mutable Gaudi::Accumulators::StatAccumulator<Gaudi::Accumulators::atomicity::full, double> m_stat_pully;
+ mutable Gaudi::Accumulators::StatAccumulator<Gaudi::Accumulators::atomicity::full, double> m_stat_pully_sq;
+
+ mutable Gaudi::Accumulators::StatAccumulator<Gaudi::Accumulators::atomicity::full, double> m_stat_pullz;
+ mutable Gaudi::Accumulators::StatAccumulator<Gaudi::Accumulators::atomicity::full, double> m_stat_pullz_sq;
+
+ mutable std::vector<double> dx_vector{}, dy_vector{}, dz_vector{}, pullx_vector{}, pully_vector{}, pullz_vector{};
struct monitoringHistos {
mutable std::array<Gaudi::Accumulators::Histogram<1>, 5> m_histo_nTracksBins_dx;
mutable std::array<Gaudi::Accumulators::Histogram<1>, 5> m_histo_nTracksBins_dy;
@@ -276,28 +253,26 @@ private:
mutable Gaudi::Accumulators::Histogram<1> m_histo_pully_Monitoring;
mutable Gaudi::Accumulators::Histogram<1> m_histo_pullz_Monitoring;
- template <std::size_t... IDXs>
+ template <std::size_t... IDXs>
static std::array<Gaudi::Accumulators::Histogram<1>, sizeof...( IDXs )>
- histo1DArrayBuilder( const VertexCompare* owner, const std::string& name, const std::string& title,
- std::tuple<unsigned, double, double> xbins, std::index_sequence<IDXs...> ) {
+ histo1DArrayBuilder( const VertexCompare* owner, const std::string& name, Gaudi::Histo1DDef def, std::index_sequence<IDXs...> ) {
return {{{owner,
name + std::to_string( IDXs ),
- title + std::to_string( IDXs ),
- {std::get<0>( xbins ), std::get<1>( xbins ), std::get<2>( xbins )}}...}};
+ def.title() + std::to_string( IDXs ),
+ {static_cast<unsigned int>(def.bins()), def.lowEdge(), def.highEdge()}}...}};
+
}
- Gaudi::Histo1DDef m_histodx{"dx, mm", -0.15, 0.15, 50};
monitoringHistos( const VertexCompare* owner )
- : m_histo_nTracksBins_dx{histo1DArrayBuilder( owner, "dx_Monitoring_", "dx, mm", {50, -0.15, 0.15},
+ : m_histo_nTracksBins_dx{histo1DArrayBuilder( owner, "dx_Monitoring_", Gaudi::Histo1DDef{ "dx, mm", -0.15, 0.15, 50},
std::make_index_sequence<5>() )}
- , m_histo_nTracksBins_dy{histo1DArrayBuilder( owner, "dy_Monitoring_", "dy, mm", {50, -0.15, 0.15},
+ , m_histo_nTracksBins_dy{histo1DArrayBuilder( owner, "dy_Monitoring_", Gaudi::Histo1DDef{ "dy, mm", -0.15, 0.15, 50},
std::make_index_sequence<5>() )}
- , m_histo_nTracksBins_dz{histo1DArrayBuilder( owner, "dz_Monitoring_", "dz, mm", {50, -1.5, 1.5},
+ , m_histo_nTracksBins_dz{histo1DArrayBuilder( owner, "dz_Monitoring_", Gaudi::Histo1DDef{ "dz, mm", -1.5, 1.5, 50},
std::make_index_sequence<5>() )}
, m_histo_pullx_Monitoring{owner, "pullx_Monitoring", "pull x", {20, -5, 5}}
, m_histo_pully_Monitoring{owner, "pully_Monitoring", "pull y", {20, -5, 5}}
, m_histo_pullz_Monitoring{owner, "pullz_Monitoring", "pull z", {20, -5, 5}} {}
};
-
std::unique_ptr<monitoringHistos> m_monitoringHistos;
using axis1D = Gaudi::Accumulators::Axis<Gaudi::Accumulators::Histogram<1>::AxisArithmeticType>;
@@ -496,6 +471,9 @@ void VertexCompare::operator()( Vertices const& recoVtx1, Vertices const& recoVt
pv_rank = link.at( oIt ) + 1;
}
+ double pullx = dx / errx;
+ double pully = dy / erry;
+ double pullz = dz / errz;
if ( m_monitoring.value() ) {
if ( dz < 2 && size_diff < 3 && pv_rank < 2 && ( ntracks1 + ntracks2 ) / 2 > ntrack_bins[0] ) {
int binCount = 0;
@@ -507,27 +485,40 @@ void VertexCompare::operator()( Vertices const& recoVtx1, Vertices const& recoVt
++monitoringHistos.m_histo_nTracksBins_dx[binCount][dx];
++monitoringHistos.m_histo_nTracksBins_dy[binCount][dy];
++monitoringHistos.m_histo_nTracksBins_dz[binCount][dz];
- ++monitoringHistos.m_histo_pullx_Monitoring[dx / errx];
- ++monitoringHistos.m_histo_pully_Monitoring[dy / erry];
- ++monitoringHistos.m_histo_pullz_Monitoring[dz / errz];
+ ++monitoringHistos.m_histo_pullx_Monitoring[pullx];
+ ++monitoringHistos.m_histo_pully_Monitoring[pully];
+ ++monitoringHistos.m_histo_pullz_Monitoring[pullz];
}
}
if ( m_produceHistogram.value() ) {
dx_vector.push_back( dx );
dy_vector.push_back( dy );
dz_vector.push_back( dz );
- pullx_vector.push_back( dx / errx );
- pully_vector.push_back( dy / erry );
- pullz_vector.push_back( dz / errz );
+ pullx_vector.push_back( pullx );
+ pully_vector.push_back( pully );
+ pullz_vector.push_back( pullz );
++m_histo_dx[dx];
++m_histo_dy[dy];
++m_histo_dz[dz];
- ++m_histo_pullx[dx / errx];
- ++m_histo_pully[dy / erry];
- ++m_histo_pullz[dz / errz];
+ ++m_histo_pullx[pullx];
+ ++m_histo_pully[pully];
+ ++m_histo_pullz[pullz];
++m_nTracks1[ntracks1];
++m_nTracks2[ntracks2];
++m_nTracks_dif[ntracks2 - ntracks1];
+
+ m_stat_dx += dx;
+ m_stat_dx_sq += dx*dx;
+ m_stat_dy += dy;
+ m_stat_dy_sq += dy*dy;
+ m_stat_dz += dz;
+ m_stat_dz_sq += dz*dz;
+ m_stat_pullx += pullx;
+ m_stat_pullx_sq += pullx * pullx ;
+ m_stat_pully += pully;
+ m_stat_pully_sq += pully * pully ;
+ m_stat_pullz += pullz;
+ m_stat_pullz_sq += pullz * pullz ;
} // else end
if ( m_produceNtuple.value() ) {
@@ -741,55 +732,54 @@ StatusCode VertexCompare::finalize() {
for ( double value : pully_vector ) { ++m_histo_pully[value]; }
for ( double value : pullz_vector ) { ++m_histo_pullz[value]; }
+
+ info() << " stat = " << m_stat_dx.standard_deviation() << endmsg;
+ info() << " stat = " << m_stat_dx_sq.sum2() << endmsg;
info() << "dx: "
- << format( "mean = %5.3f +/- %5.3f, RMS = %5.3f +/- %5.3f", m_histo_dx.mean(), meanErr( m_histo_dx ),
- m_histo_dx.standard_deviation(), rmsErr( m_histo_dx ) )
+ << format( "mean = %5.3f +/- %5.3f, RMS = %5.3f +/- %5.3f", static_cast<double>(m_stat_dx.sum())/m_stat_dx.nEntries(), m_stat_dx.meanErr(),
+ m_stat_dx.standard_deviation(), rmsErr( m_stat_dx , m_stat_dx_sq) )
<< endmsg;
-
info() << "dy: "
- << format( "mean = %5.3f +/- %5.3f, RMS = %5.3f +/- %5.3f", m_histo_dy.mean(), meanErr( m_histo_dy ),
- m_histo_dy.standard_deviation(), rmsErr( m_histo_dy ) )
+ << format( "mean = %5.3f +/- %5.3f, RMS = %5.3f +/- %5.3f", static_cast<double>(m_stat_dy.sum())/m_stat_dy.nEntries(), m_stat_dy.meanErr(),
+ m_stat_dy.standard_deviation(), rmsErr( m_stat_dy , m_stat_dy_sq))
<< endmsg;
-
info() << "dz: "
- << format( "mean = %5.3f +/- %5.3f, RMS = %5.3f +/- %5.3f", m_histo_dz.mean(), meanErr( m_histo_dz ),
- m_histo_dz.standard_deviation(), rmsErr( m_histo_dz ) )
- << endmsg;
-
+ << format( "mean = %5.3f +/- %5.3f, RMS = %5.3f +/- %5.3f", static_cast<double>(m_stat_dz.sum())/m_stat_dz.nEntries(), m_stat_dz.meanErr(),
+ m_stat_dz.standard_deviation(), rmsErr( m_stat_dz , m_stat_dz_sq) )
+ << endmsg;
info() << " ---------------------------------------" << endmsg;
info() << "pullx: "
- << format( "mean = %5.3f +/- %5.3f, RMS = %5.3f +/- %5.3f", m_histo_pullx.mean(), meanErr( m_histo_pullx ),
- m_histo_pullx.standard_deviation(), rmsErr( m_histo_pullx ) )
+ << format( "mean = %5.3f +/- %5.3f, RMS = %5.3f +/- %5.3f", static_cast<double>(m_stat_pullx.sum())/m_stat_pullx.nEntries(), m_stat_pullx.meanErr(),
+ m_stat_pullx.standard_deviation(), rmsErr( m_stat_pullx , m_stat_pullx_sq) )
<< endmsg;
info() << "pully: "
- << format( "mean = %5.3f +/- %5.3f, RMS = %5.3f +/- %5.3f", m_histo_pully.mean(), meanErr( m_histo_pully ),
- m_histo_pully.standard_deviation(), rmsErr( m_histo_pully ) )
+ << format( "mean = %5.3f +/- %5.3f, RMS = %5.3f +/- %5.3f", static_cast<double>(m_stat_pully.sum())/m_stat_pully.nEntries(), m_stat_pully.meanErr(),
+ m_stat_pully.standard_deviation(), rmsErr( m_stat_pully , m_stat_pully_sq) )
<< endmsg;
info() << "pullz: "
- << format( "mean = %5.3f +/- %5.3f, RMS = %5.3f +/- %5.3f", m_histo_pullz.mean(), meanErr( m_histo_pullz ),
- m_histo_pullz.standard_deviation(), rmsErr( m_histo_pullz ) )
+ << format( "mean = %5.3f +/- %5.3f, RMS = %5.3f +/- %5.3f", static_cast<double>(m_stat_pullz.sum())/m_stat_pullz.nEntries(), m_stat_pullz.meanErr(),
+ m_stat_pullz.standard_deviation(), rmsErr( m_stat_pullz , m_stat_pullz_sq) )
<< endmsg;
info() << " ---------------------------------------" << endmsg;
info() << "diff in x: "
- << format( "%5.3f +/- %5.3f", std::sqrt( 1.0 + m_histo_pullx.standard_deviation() ) - 1.0,
- rmsErr( m_histo_pullx ) * 0.5 / std::sqrt( 1.0 + m_histo_pullx.standard_deviation() ) )
+ << format( "%5.3f +/- %5.3f", std::sqrt( 1.0 + m_stat_pullx.standard_deviation() ) - 1.0,
+ rmsErr( m_stat_pullx, m_stat_pullx_sq ) * 0.5 / std::sqrt( 1.0 + m_stat_pullx.standard_deviation() ) )
<< endmsg;
info() << "diff in y: "
- << format( "%5.3f +/- %5.3f", std::sqrt( 1.0 + m_histo_pully.standard_deviation() ) - 1.0,
- rmsErr( m_histo_pully ) * 0.5 / std::sqrt( 1.0 + m_histo_pully.standard_deviation() ) )
+ << format( "%5.3f +/- %5.3f", std::sqrt( 1.0 + m_stat_pully.standard_deviation() ) - 1.0,
+ rmsErr( m_stat_pully, m_stat_pully_sq ) * 0.5 / std::sqrt( 1.0 + m_stat_pully.standard_deviation() ) )
<< endmsg;
info() << "diff in z: "
- << format( "%5.3f +/- %5.3f", std::sqrt( 1.0 + m_histo_pullz.standard_deviation() ) - 1.0,
- rmsErr( m_histo_pullz ) * 0.5 / std::sqrt( 1.0 + m_histo_pullz.standard_deviation() ) )
+ << format( "%5.3f +/- %5.3f", std::sqrt( 1.0 + m_stat_pullz.standard_deviation() ) - 1.0,
+ rmsErr( m_stat_pullz, m_stat_pullz_sq ) * 0.5 / std::sqrt( 1.0 + m_stat_pullz.standard_deviation() ) )
<< endmsg;
-
info() << " ============================================" << endmsg;
}
return Consumer::finalize(); // Must be called after all other actions
diff --git a/Tr/TrackUtils/src/RandomTrackContainerSplitter.cpp b/Tr/TrackUtils/src/RandomTrackContainerSplitter.cpp
index f9a9ac691f4..5df05286cd8 100644
--- a/Tr/TrackUtils/src/RandomTrackContainerSplitter.cpp
+++ b/Tr/TrackUtils/src/RandomTrackContainerSplitter.cpp
@@ -15,22 +15,18 @@
*
* @author Andrii Usachov
* @date 30/04/2020
+ * @author Bogdan Kutsenko
+ * @date 27/05/2020
*/
#include "Event/Track.h"
#include "Functors/Function.h"
#include "Functors/with_functors.h"
-#include "GaudiKernel/IEventTimeDecoder.h"
-#include "GaudiKernel/IRndmEngine.h"
-#include "GaudiKernel/IRndmGenSvc.h"
-#include "GaudiKernel/RndmGenerators.h"
-#include "GaudiKernel/Service.h" // Include the header file for Service
-#include "GaudiKernel/ServiceHandle.h"
#include "LHCbAlgs/Transformer.h"
namespace {
/// The output data
- using OutConts = std::tuple<LHCb::Tracks, LHCb::Tracks>;
+ using OutConts = std::tuple<LHCb::Track::Selection, LHCb::Track::Selection>;
struct TrackPredicate {
using Signature = bool( const LHCb::Track& );
@@ -38,40 +34,25 @@ namespace {
};
} // namespace
-bool makeRandomDecision( const LHCb::Track& track, ServiceHandle<IRndmGenSvc>& rndmSvc ) {
- if ( !rndmSvc.retrieve().isSuccess() ) {
- std::cerr << "Error: Failed to retrieve IRndmGenSvc service." << std::endl;
- return false;
- }
-
- IRndmGenSvc* rndmGenSvcPtr = rndmSvc.operator->();
-
- if ( !rndmGenSvcPtr ) {
- std::cerr << "Error: IRndmGenSvc pointer is null." << std::endl;
- return false;
- }
-
- Rndm::Numbers rndm( rndmGenSvcPtr, Rndm::Flat( 0., 1.0 ) );
- double randomNumber = rndm();
-
- return randomNumber < 0.5;
+bool makePesudoRandomDecision( const LHCb::Track& track) {
+ const auto lhcb_id = track.lhcbIDs()[0];
+ return static_cast<int>( lhcb_id.lhcbID() ) % 2;
}
class RandomTrackContainerSplitter final
- : public with_functors<LHCb::Algorithm::MultiTransformer<OutConts( const LHCb::Tracks& )>, TrackPredicate> {
+ : public with_functors<LHCb::Algorithm::MultiTransformer<OutConts( const LHCb::Track::Range& )>, TrackPredicate> {
public:
RandomTrackContainerSplitter( const std::string& name, ISvcLocator* pSvcLocator )
: with_functors( name, pSvcLocator, {KeyValue{"TracksInContainer", ""}},
{KeyValue{"PassedContainer", ""}, KeyValue{"RejectedContainer", ""}} ) {}
- OutConts operator()( const LHCb::Tracks& input_tracks ) const override {
+ OutConts operator()( const LHCb::Track::Range& input_tracks ) const override {
OutConts out;
- ServiceHandle<IRndmGenSvc> rndmSvc( "RndmGenSvc", "RandomTrackContainerSplitter" );
auto& [passed, rejected] = out;
m_inputs += input_tracks.size();
for ( auto& trk : input_tracks ) {
- auto& container = ( makeRandomDecision( *trk, rndmSvc ) ? passed : rejected );
- container.insert( new LHCb::Track( *trk ) );
+ auto& container = ( makePesudoRandomDecision( *trk) ? passed : rejected );
+ container.insert( trk );
}
m_passed += passed.size();
diff --git a/Tr/TrackUtils/src/RandomVeloTrackContainerSplitter.cpp b/Tr/TrackUtils/src/RandomVeloTrackContainerSplitter.cpp
index 7188900e766..8d3b5a48049 100644
--- a/Tr/TrackUtils/src/RandomVeloTrackContainerSplitter.cpp
+++ b/Tr/TrackUtils/src/RandomVeloTrackContainerSplitter.cpp
@@ -22,41 +22,21 @@
#include "Event/Track.h"
#include "Functors/Function.h"
#include "Functors/with_functors.h"
-#include "GaudiKernel/IEventTimeDecoder.h"
-#include "GaudiKernel/IRndmEngine.h"
-#include "GaudiKernel/IRndmGenSvc.h"
-#include "GaudiKernel/RndmGenerators.h"
-#include "GaudiKernel/Service.h" // Include the header file for Service
-#include "GaudiKernel/ServiceHandle.h"
#include "LHCbAlgs/Transformer.h"
namespace {
/// The output data
using OutConts = std::tuple<LHCb::Pr::Velo::Tracks, LHCb::Pr::Velo::Tracks>;
-
+ using PrTrackProxy = decltype( *std::declval<const LHCb::Pr::Velo::Tracks>().scalar().begin() );
struct TrackPredicate {
using Signature = bool( const LHCb::Pr::Velo::Tracks& );
static constexpr auto PropertyName = "Code";
};
} // namespace
-bool makeRandomDecision( const LHCb::Pr::Velo::Tracks&, ServiceHandle<IRndmGenSvc>& rndmSvc ) {
- if ( !rndmSvc.retrieve().isSuccess() ) {
- std::cerr << "Error: Failed to retrieve IRndmGenSvc service." << std::endl;
- return false;
- }
-
- IRndmGenSvc* rndmGenSvcPtr = rndmSvc.operator->();
-
- if ( !rndmGenSvcPtr ) {
- std::cerr << "Error: IRndmGenSvc pointer is null." << std::endl;
- return false;
- }
-
- Rndm::Numbers rndm( rndmGenSvcPtr, Rndm::Flat( 0., 1.0 ) );
- double randomNumber = rndm();
-
- return randomNumber < 0.5;
+bool makePseudoRandomDecision( const PrTrackProxy& track) {
+ const auto lhcb_id = track.lhcbIDs()[0];
+ return static_cast<int>( lhcb_id.lhcbID() ) % 2;
}
class RandomVeloTrackContainerSplitter final
@@ -69,16 +49,12 @@ public:
OutConts operator()( const LHCb::Pr::Velo::Tracks& input_tracks ) const override {
OutConts out;
- ServiceHandle<IRndmGenSvc> rndmSvc( "RndmGenSvc", "RandomVeloTrackContainerSplitter" );
auto& [passed, rejected] = out;
m_inputs += input_tracks.size();
for ( auto const& trk : input_tracks.scalar() ) {
+ auto& container = ( makePseudoRandomDecision( trk ) ? passed : rejected );
const int t = trk.offset();
- LHCb::Pr::Velo::Tracks tmpTrk;
- tmpTrk.reserve( 1 );
- tmpTrk.copy_back<SIMDWrapper::scalar::types>( input_tracks, t, true );
- auto& container = ( makeRandomDecision( tmpTrk, rndmSvc ) ? passed : rejected );
container.copy_back<SIMDWrapper::scalar::types>( input_tracks, t, true );
}
m_passed += passed.size();
--
GitLab
From 59aa931da6f642a8556980492348afa6f408b5b4 Mon Sep 17 00:00:00 2001
From: Gitlab CI <noreply@cern.ch>
Date: Mon, 27 May 2024 16:30:21 +0000
Subject: [PATCH 16/21] Fixed formatting
patch generated by https://gitlab.cern.ch/lhcb/Rec/-/jobs/39371144
---
Tr/TrackMonitors/src/VertexCompare.cpp | 138 +++++++++---------
.../src/RandomTrackContainerSplitter.cpp | 6 +-
.../src/RandomVeloTrackContainerSplitter.cpp | 10 +-
3 files changed, 81 insertions(+), 73 deletions(-)
diff --git a/Tr/TrackMonitors/src/VertexCompare.cpp b/Tr/TrackMonitors/src/VertexCompare.cpp
index 0f93b548b61..dfef252ff00 100644
--- a/Tr/TrackMonitors/src/VertexCompare.cpp
+++ b/Tr/TrackMonitors/src/VertexCompare.cpp
@@ -21,9 +21,9 @@
#include "GaudiAlg/GaudiAlgorithm.h"
#include "GaudiAlg/GaudiTupleAlg.h"
#include "GaudiAlg/Tuples.h"
+#include "GaudiKernel/HistoDef.h"
#include "GaudiKernel/PhysicalConstants.h"
#include "GaudiKernel/SystemOfUnits.h"
-#include "GaudiKernel/HistoDef.h"
#include "GaudiUtils/HistoStats.h"
#include "LHCbAlgs/Consumer.h"
#include "MCInterfaces/IForcedBDecayTool.h"
@@ -90,31 +90,31 @@ bool trackcomp( const RecPVInfo<VertexType>& first, const RecPVInfo<VertexType>&
}
struct VtxVariables {
- int ntracks = -9999;
- double x = -9999.;
- double y = -9999.;
- double z = -9999.;
- double dxr = -9999.;
- double dyr = -9999.;
- double r = -9999.;
- double errx = -9999.;
- double erry = -9999.;
- double errz = -9999.;
- double err_r = -9999.;
- double covxx = -9999.;
- double covyy = -9999.;
- double covzz = -9999.;
- double covxy = -9999.;
- double covxz = -9999.;
- double covyz = -9999.;
- double chi2 = -9999.;
- double nDoF = -9999.;
- int dsize = -9999;
- int match = -9999;
- int pv_rank = -9999;
- bool equal_sizes = false;
- bool single = false;
- bool opposite_container = false;
+ int ntracks = -9999;
+ double x = -9999.;
+ double y = -9999.;
+ double z = -9999.;
+ double dxr = -9999.;
+ double dyr = -9999.;
+ double r = -9999.;
+ double errx = -9999.;
+ double erry = -9999.;
+ double errz = -9999.;
+ double err_r = -9999.;
+ double covxx = -9999.;
+ double covyy = -9999.;
+ double covzz = -9999.;
+ double covxy = -9999.;
+ double covxz = -9999.;
+ double covyz = -9999.;
+ double chi2 = -9999.;
+ double nDoF = -9999.;
+ int dsize = -9999;
+ int match = -9999;
+ int pv_rank = -9999;
+ bool equal_sizes = false;
+ bool single = false;
+ bool opposite_container = false;
};
VtxVariables SetVtxVariables( const RecPVInfo<VertexType>& vrtf, const int& counter, const int& size1, const int& size2,
@@ -177,9 +177,9 @@ VtxVariables SetVtxVariables( const RecPVInfo<VertexType>& vrtf, const int& coun
// ============================================================================
template <typename StatAccumulator>
-double kurtosis( StatAccumulator const& stAcc , StatAccumulator const& stAcc_sq ) {
- const auto n = stAcc.nEntries();
- const auto mu4 = stAcc_sq.sum2()/n; // Calculate a 4t moment using statAccumulater of squared variables
+double kurtosis( StatAccumulator const& stAcc, StatAccumulator const& stAcc_sq ) {
+ const auto n = stAcc.nEntries();
+ const auto mu4 = stAcc_sq.sum2() / n; // Calculate a 4t moment using statAccumulater of squared variables
const auto s4 = std::pow( stAcc.standard_deviation(), 4 );
return ( std::fabs( s4 ) > 0 ? mu4 / s4 - 3.0 : 0.0 );
}
@@ -189,10 +189,10 @@ double kurtosis( StatAccumulator const& stAcc , StatAccumulator const& stAcc_sq
// ============================================================================
template <typename StatAccumulator>
-double rmsErr( StatAccumulator const& stAcc, StatAccumulator const& stAcc_sq ) {
+double rmsErr( StatAccumulator const& stAcc, StatAccumulator const& stAcc_sq ) {
const auto n = stAcc.nEntries();
if ( 1 >= n ) { return 0.0; }
- auto result = 2.0 + kurtosis( stAcc , stAcc_sq );
+ auto result = 2.0 + kurtosis( stAcc, stAcc_sq );
result += 2.0 / ( n - 1 );
result /= 4.0 * n;
return stAcc.standard_deviation() * std::sqrt( std::max( result, 0.0 ) );
@@ -243,7 +243,7 @@ private:
mutable Gaudi::Accumulators::StatAccumulator<Gaudi::Accumulators::atomicity::full, double> m_stat_pullz;
mutable Gaudi::Accumulators::StatAccumulator<Gaudi::Accumulators::atomicity::full, double> m_stat_pullz_sq;
-
+
mutable std::vector<double> dx_vector{}, dy_vector{}, dz_vector{}, pullx_vector{}, pully_vector{}, pullz_vector{};
struct monitoringHistos {
mutable std::array<Gaudi::Accumulators::Histogram<1>, 5> m_histo_nTracksBins_dx;
@@ -253,22 +253,22 @@ private:
mutable Gaudi::Accumulators::Histogram<1> m_histo_pully_Monitoring;
mutable Gaudi::Accumulators::Histogram<1> m_histo_pullz_Monitoring;
- template <std::size_t... IDXs>
+ template <std::size_t... IDXs>
static std::array<Gaudi::Accumulators::Histogram<1>, sizeof...( IDXs )>
- histo1DArrayBuilder( const VertexCompare* owner, const std::string& name, Gaudi::Histo1DDef def, std::index_sequence<IDXs...> ) {
+ histo1DArrayBuilder( const VertexCompare* owner, const std::string& name, Gaudi::Histo1DDef def,
+ std::index_sequence<IDXs...> ) {
return {{{owner,
name + std::to_string( IDXs ),
def.title() + std::to_string( IDXs ),
- {static_cast<unsigned int>(def.bins()), def.lowEdge(), def.highEdge()}}...}};
-
+ {static_cast<unsigned int>( def.bins() ), def.lowEdge(), def.highEdge()}}...}};
}
monitoringHistos( const VertexCompare* owner )
- : m_histo_nTracksBins_dx{histo1DArrayBuilder( owner, "dx_Monitoring_", Gaudi::Histo1DDef{ "dx, mm", -0.15, 0.15, 50},
- std::make_index_sequence<5>() )}
- , m_histo_nTracksBins_dy{histo1DArrayBuilder( owner, "dy_Monitoring_", Gaudi::Histo1DDef{ "dy, mm", -0.15, 0.15, 50},
- std::make_index_sequence<5>() )}
- , m_histo_nTracksBins_dz{histo1DArrayBuilder( owner, "dz_Monitoring_", Gaudi::Histo1DDef{ "dz, mm", -1.5, 1.5, 50},
- std::make_index_sequence<5>() )}
+ : m_histo_nTracksBins_dx{histo1DArrayBuilder(
+ owner, "dx_Monitoring_", Gaudi::Histo1DDef{"dx, mm", -0.15, 0.15, 50}, std::make_index_sequence<5>() )}
+ , m_histo_nTracksBins_dy{histo1DArrayBuilder(
+ owner, "dy_Monitoring_", Gaudi::Histo1DDef{"dy, mm", -0.15, 0.15, 50}, std::make_index_sequence<5>() )}
+ , m_histo_nTracksBins_dz{histo1DArrayBuilder(
+ owner, "dz_Monitoring_", Gaudi::Histo1DDef{"dz, mm", -1.5, 1.5, 50}, std::make_index_sequence<5>() )}
, m_histo_pullx_Monitoring{owner, "pullx_Monitoring", "pull x", {20, -5, 5}}
, m_histo_pully_Monitoring{owner, "pully_Monitoring", "pull y", {20, -5, 5}}
, m_histo_pullz_Monitoring{owner, "pullz_Monitoring", "pull z", {20, -5, 5}} {}
@@ -508,17 +508,17 @@ void VertexCompare::operator()( Vertices const& recoVtx1, Vertices const& recoVt
++m_nTracks_dif[ntracks2 - ntracks1];
m_stat_dx += dx;
- m_stat_dx_sq += dx*dx;
+ m_stat_dx_sq += dx * dx;
m_stat_dy += dy;
- m_stat_dy_sq += dy*dy;
+ m_stat_dy_sq += dy * dy;
m_stat_dz += dz;
- m_stat_dz_sq += dz*dz;
+ m_stat_dz_sq += dz * dz;
m_stat_pullx += pullx;
- m_stat_pullx_sq += pullx * pullx ;
+ m_stat_pullx_sq += pullx * pullx;
m_stat_pully += pully;
- m_stat_pully_sq += pully * pully ;
+ m_stat_pully_sq += pully * pully;
m_stat_pullz += pullz;
- m_stat_pullz_sq += pullz * pullz ;
+ m_stat_pullz_sq += pullz * pullz;
} // else end
if ( m_produceNtuple.value() ) {
@@ -732,53 +732,61 @@ StatusCode VertexCompare::finalize() {
for ( double value : pully_vector ) { ++m_histo_pully[value]; }
for ( double value : pullz_vector ) { ++m_histo_pullz[value]; }
-
info() << " stat = " << m_stat_dx.standard_deviation() << endmsg;
info() << " stat = " << m_stat_dx_sq.sum2() << endmsg;
info() << "dx: "
- << format( "mean = %5.3f +/- %5.3f, RMS = %5.3f +/- %5.3f", static_cast<double>(m_stat_dx.sum())/m_stat_dx.nEntries(), m_stat_dx.meanErr(),
- m_stat_dx.standard_deviation(), rmsErr( m_stat_dx , m_stat_dx_sq) )
+ << format( "mean = %5.3f +/- %5.3f, RMS = %5.3f +/- %5.3f",
+ static_cast<double>( m_stat_dx.sum() ) / m_stat_dx.nEntries(), m_stat_dx.meanErr(),
+ m_stat_dx.standard_deviation(), rmsErr( m_stat_dx, m_stat_dx_sq ) )
<< endmsg;
info() << "dy: "
- << format( "mean = %5.3f +/- %5.3f, RMS = %5.3f +/- %5.3f", static_cast<double>(m_stat_dy.sum())/m_stat_dy.nEntries(), m_stat_dy.meanErr(),
- m_stat_dy.standard_deviation(), rmsErr( m_stat_dy , m_stat_dy_sq))
+ << format( "mean = %5.3f +/- %5.3f, RMS = %5.3f +/- %5.3f",
+ static_cast<double>( m_stat_dy.sum() ) / m_stat_dy.nEntries(), m_stat_dy.meanErr(),
+ m_stat_dy.standard_deviation(), rmsErr( m_stat_dy, m_stat_dy_sq ) )
<< endmsg;
info() << "dz: "
- << format( "mean = %5.3f +/- %5.3f, RMS = %5.3f +/- %5.3f", static_cast<double>(m_stat_dz.sum())/m_stat_dz.nEntries(), m_stat_dz.meanErr(),
- m_stat_dz.standard_deviation(), rmsErr( m_stat_dz , m_stat_dz_sq) )
- << endmsg;
+ << format( "mean = %5.3f +/- %5.3f, RMS = %5.3f +/- %5.3f",
+ static_cast<double>( m_stat_dz.sum() ) / m_stat_dz.nEntries(), m_stat_dz.meanErr(),
+ m_stat_dz.standard_deviation(), rmsErr( m_stat_dz, m_stat_dz_sq ) )
+ << endmsg;
info() << " ---------------------------------------" << endmsg;
info() << "pullx: "
- << format( "mean = %5.3f +/- %5.3f, RMS = %5.3f +/- %5.3f", static_cast<double>(m_stat_pullx.sum())/m_stat_pullx.nEntries(), m_stat_pullx.meanErr(),
- m_stat_pullx.standard_deviation(), rmsErr( m_stat_pullx , m_stat_pullx_sq) )
+ << format( "mean = %5.3f +/- %5.3f, RMS = %5.3f +/- %5.3f",
+ static_cast<double>( m_stat_pullx.sum() ) / m_stat_pullx.nEntries(), m_stat_pullx.meanErr(),
+ m_stat_pullx.standard_deviation(), rmsErr( m_stat_pullx, m_stat_pullx_sq ) )
<< endmsg;
info() << "pully: "
- << format( "mean = %5.3f +/- %5.3f, RMS = %5.3f +/- %5.3f", static_cast<double>(m_stat_pully.sum())/m_stat_pully.nEntries(), m_stat_pully.meanErr(),
- m_stat_pully.standard_deviation(), rmsErr( m_stat_pully , m_stat_pully_sq) )
+ << format( "mean = %5.3f +/- %5.3f, RMS = %5.3f +/- %5.3f",
+ static_cast<double>( m_stat_pully.sum() ) / m_stat_pully.nEntries(), m_stat_pully.meanErr(),
+ m_stat_pully.standard_deviation(), rmsErr( m_stat_pully, m_stat_pully_sq ) )
<< endmsg;
info() << "pullz: "
- << format( "mean = %5.3f +/- %5.3f, RMS = %5.3f +/- %5.3f", static_cast<double>(m_stat_pullz.sum())/m_stat_pullz.nEntries(), m_stat_pullz.meanErr(),
- m_stat_pullz.standard_deviation(), rmsErr( m_stat_pullz , m_stat_pullz_sq) )
+ << format( "mean = %5.3f +/- %5.3f, RMS = %5.3f +/- %5.3f",
+ static_cast<double>( m_stat_pullz.sum() ) / m_stat_pullz.nEntries(), m_stat_pullz.meanErr(),
+ m_stat_pullz.standard_deviation(), rmsErr( m_stat_pullz, m_stat_pullz_sq ) )
<< endmsg;
info() << " ---------------------------------------" << endmsg;
info() << "diff in x: "
<< format( "%5.3f +/- %5.3f", std::sqrt( 1.0 + m_stat_pullx.standard_deviation() ) - 1.0,
- rmsErr( m_stat_pullx, m_stat_pullx_sq ) * 0.5 / std::sqrt( 1.0 + m_stat_pullx.standard_deviation() ) )
+ rmsErr( m_stat_pullx, m_stat_pullx_sq ) * 0.5 /
+ std::sqrt( 1.0 + m_stat_pullx.standard_deviation() ) )
<< endmsg;
info() << "diff in y: "
<< format( "%5.3f +/- %5.3f", std::sqrt( 1.0 + m_stat_pully.standard_deviation() ) - 1.0,
- rmsErr( m_stat_pully, m_stat_pully_sq ) * 0.5 / std::sqrt( 1.0 + m_stat_pully.standard_deviation() ) )
+ rmsErr( m_stat_pully, m_stat_pully_sq ) * 0.5 /
+ std::sqrt( 1.0 + m_stat_pully.standard_deviation() ) )
<< endmsg;
info() << "diff in z: "
<< format( "%5.3f +/- %5.3f", std::sqrt( 1.0 + m_stat_pullz.standard_deviation() ) - 1.0,
- rmsErr( m_stat_pullz, m_stat_pullz_sq ) * 0.5 / std::sqrt( 1.0 + m_stat_pullz.standard_deviation() ) )
+ rmsErr( m_stat_pullz, m_stat_pullz_sq ) * 0.5 /
+ std::sqrt( 1.0 + m_stat_pullz.standard_deviation() ) )
<< endmsg;
info() << " ============================================" << endmsg;
}
diff --git a/Tr/TrackUtils/src/RandomTrackContainerSplitter.cpp b/Tr/TrackUtils/src/RandomTrackContainerSplitter.cpp
index 5df05286cd8..a40cc8844d7 100644
--- a/Tr/TrackUtils/src/RandomTrackContainerSplitter.cpp
+++ b/Tr/TrackUtils/src/RandomTrackContainerSplitter.cpp
@@ -34,7 +34,7 @@ namespace {
};
} // namespace
-bool makePesudoRandomDecision( const LHCb::Track& track) {
+bool makePesudoRandomDecision( const LHCb::Track& track ) {
const auto lhcb_id = track.lhcbIDs()[0];
return static_cast<int>( lhcb_id.lhcbID() ) % 2;
}
@@ -47,11 +47,11 @@ public:
{KeyValue{"PassedContainer", ""}, KeyValue{"RejectedContainer", ""}} ) {}
OutConts operator()( const LHCb::Track::Range& input_tracks ) const override {
- OutConts out;
+ OutConts out;
auto& [passed, rejected] = out;
m_inputs += input_tracks.size();
for ( auto& trk : input_tracks ) {
- auto& container = ( makePesudoRandomDecision( *trk) ? passed : rejected );
+ auto& container = ( makePesudoRandomDecision( *trk ) ? passed : rejected );
container.insert( trk );
}
m_passed += passed.size();
diff --git a/Tr/TrackUtils/src/RandomVeloTrackContainerSplitter.cpp b/Tr/TrackUtils/src/RandomVeloTrackContainerSplitter.cpp
index 8d3b5a48049..cb735e11bdf 100644
--- a/Tr/TrackUtils/src/RandomVeloTrackContainerSplitter.cpp
+++ b/Tr/TrackUtils/src/RandomVeloTrackContainerSplitter.cpp
@@ -26,7 +26,7 @@
namespace {
/// The output data
- using OutConts = std::tuple<LHCb::Pr::Velo::Tracks, LHCb::Pr::Velo::Tracks>;
+ using OutConts = std::tuple<LHCb::Pr::Velo::Tracks, LHCb::Pr::Velo::Tracks>;
using PrTrackProxy = decltype( *std::declval<const LHCb::Pr::Velo::Tracks>().scalar().begin() );
struct TrackPredicate {
using Signature = bool( const LHCb::Pr::Velo::Tracks& );
@@ -34,7 +34,7 @@ namespace {
};
} // namespace
-bool makePseudoRandomDecision( const PrTrackProxy& track) {
+bool makePseudoRandomDecision( const PrTrackProxy& track ) {
const auto lhcb_id = track.lhcbIDs()[0];
return static_cast<int>( lhcb_id.lhcbID() ) % 2;
}
@@ -48,13 +48,13 @@ public:
{KeyValue{"PassedContainer", ""}, KeyValue{"RejectedContainer", ""}} ) {}
OutConts operator()( const LHCb::Pr::Velo::Tracks& input_tracks ) const override {
- OutConts out;
+ OutConts out;
auto& [passed, rejected] = out;
m_inputs += input_tracks.size();
for ( auto const& trk : input_tracks.scalar() ) {
- auto& container = ( makePseudoRandomDecision( trk ) ? passed : rejected );
- const int t = trk.offset();
+ auto& container = ( makePseudoRandomDecision( trk ) ? passed : rejected );
+ const int t = trk.offset();
container.copy_back<SIMDWrapper::scalar::types>( input_tracks, t, true );
}
m_passed += passed.size();
--
GitLab
From a745a4148b118f75e4a4e0235bcb98693d4e1014 Mon Sep 17 00:00:00 2001
From: Bogdan Kutsenko <bogdan.kutsenko@cern.ch>
Date: Mon, 27 May 2024 20:18:28 +0200
Subject: [PATCH 17/21] Sum for rmsError accumulator is added
---
Tr/TrackMonitors/src/VertexCompare.cpp | 56 +++++++++++++-------------
1 file changed, 27 insertions(+), 29 deletions(-)
diff --git a/Tr/TrackMonitors/src/VertexCompare.cpp b/Tr/TrackMonitors/src/VertexCompare.cpp
index dfef252ff00..abf0da1513a 100644
--- a/Tr/TrackMonitors/src/VertexCompare.cpp
+++ b/Tr/TrackMonitors/src/VertexCompare.cpp
@@ -176,10 +176,10 @@ VtxVariables SetVtxVariables( const RecPVInfo<VertexType>& vrtf, const int& coun
// get the error in kurtosis
// ============================================================================
-template <typename StatAccumulator>
-double kurtosis( StatAccumulator const& stAcc, StatAccumulator const& stAcc_sq ) {
+template <typename StatAccumulator, typename FourthSumAccumulator>
+double kurtosis( StatAccumulator const& stAcc, FourthSumAccumulator const& stAcc_fourth_sum ) {
const auto n = stAcc.nEntries();
- const auto mu4 = stAcc_sq.sum2() / n; // Calculate a 4t moment using statAccumulater of squared variables
+ const auto mu4 = stAcc_fourth_sum.sum() / n; // Calculate a 4t moment using statAccumulater of squared variables
const auto s4 = std::pow( stAcc.standard_deviation(), 4 );
return ( std::fabs( s4 ) > 0 ? mu4 / s4 - 3.0 : 0.0 );
}
@@ -188,11 +188,11 @@ double kurtosis( StatAccumulator const& stAcc, StatAccumulator const& stAcc_sq )
// get an error in the rms value
// ============================================================================
-template <typename StatAccumulator>
-double rmsErr( StatAccumulator const& stAcc, StatAccumulator const& stAcc_sq ) {
+template <typename StatAccumulator, typename FourthSumAccumulator >
+double rmsErr( StatAccumulator const& stAcc, FourthSumAccumulator const& stAcc_fourth_sum ) {
const auto n = stAcc.nEntries();
if ( 1 >= n ) { return 0.0; }
- auto result = 2.0 + kurtosis( stAcc, stAcc_sq );
+ auto result = 2.0 + kurtosis( stAcc, stAcc_fourth_sum );
result += 2.0 / ( n - 1 );
result /= 4.0 * n;
return stAcc.standard_deviation() * std::sqrt( std::max( result, 0.0 ) );
@@ -227,22 +227,22 @@ private:
Gaudi::Property<bool> m_requireSingle{this, "requireSingle", false};
mutable Gaudi::Accumulators::StatAccumulator<Gaudi::Accumulators::atomicity::full, double> m_stat_dx;
- mutable Gaudi::Accumulators::StatAccumulator<Gaudi::Accumulators::atomicity::full, double> m_stat_dx_sq;
+ mutable Gaudi::Accumulators::SumAccumulator<Gaudi::Accumulators::atomicity::full, double> m_stat_dx_fourth_sum;
mutable Gaudi::Accumulators::StatAccumulator<Gaudi::Accumulators::atomicity::full, double> m_stat_dy;
- mutable Gaudi::Accumulators::StatAccumulator<Gaudi::Accumulators::atomicity::full, double> m_stat_dy_sq;
+ mutable Gaudi::Accumulators::SumAccumulator<Gaudi::Accumulators::atomicity::full, double> m_stat_dy_fourth_sum;
mutable Gaudi::Accumulators::StatAccumulator<Gaudi::Accumulators::atomicity::full, double> m_stat_dz;
- mutable Gaudi::Accumulators::StatAccumulator<Gaudi::Accumulators::atomicity::full, double> m_stat_dz_sq;
+ mutable Gaudi::Accumulators::SumAccumulator<Gaudi::Accumulators::atomicity::full, double> m_stat_dz_fourth_sum;
mutable Gaudi::Accumulators::StatAccumulator<Gaudi::Accumulators::atomicity::full, double> m_stat_pullx;
- mutable Gaudi::Accumulators::StatAccumulator<Gaudi::Accumulators::atomicity::full, double> m_stat_pullx_sq;
+ mutable Gaudi::Accumulators::SumAccumulator<Gaudi::Accumulators::atomicity::full, double> m_stat_pullx_fourth_sum;
mutable Gaudi::Accumulators::StatAccumulator<Gaudi::Accumulators::atomicity::full, double> m_stat_pully;
- mutable Gaudi::Accumulators::StatAccumulator<Gaudi::Accumulators::atomicity::full, double> m_stat_pully_sq;
+ mutable Gaudi::Accumulators::SumAccumulator<Gaudi::Accumulators::atomicity::full, double> m_stat_pully_fourth_sum;
mutable Gaudi::Accumulators::StatAccumulator<Gaudi::Accumulators::atomicity::full, double> m_stat_pullz;
- mutable Gaudi::Accumulators::StatAccumulator<Gaudi::Accumulators::atomicity::full, double> m_stat_pullz_sq;
+ mutable Gaudi::Accumulators::SumAccumulator<Gaudi::Accumulators::atomicity::full, double> m_stat_pullz_fourth_sum;
mutable std::vector<double> dx_vector{}, dy_vector{}, dz_vector{}, pullx_vector{}, pully_vector{}, pullz_vector{};
struct monitoringHistos {
@@ -508,17 +508,17 @@ void VertexCompare::operator()( Vertices const& recoVtx1, Vertices const& recoVt
++m_nTracks_dif[ntracks2 - ntracks1];
m_stat_dx += dx;
- m_stat_dx_sq += dx * dx;
+ m_stat_dx_fourth_sum += std::pow(dx, 4);
m_stat_dy += dy;
- m_stat_dy_sq += dy * dy;
+ m_stat_dy_fourth_sum += std::pow(dy, 4);
m_stat_dz += dz;
- m_stat_dz_sq += dz * dz;
+ m_stat_dz_fourth_sum += std::pow(dz, 4);
m_stat_pullx += pullx;
- m_stat_pullx_sq += pullx * pullx;
+ m_stat_pullx_fourth_sum += std::pow(pullx, 4);
m_stat_pully += pully;
- m_stat_pully_sq += pully * pully;
+ m_stat_pully_fourth_sum += std::pow(pully, 4);
m_stat_pullz += pullz;
- m_stat_pullz_sq += pullz * pullz;
+ m_stat_pullz_fourth_sum += std::pow(pullz, 4);
} // else end
if ( m_produceNtuple.value() ) {
@@ -732,60 +732,58 @@ StatusCode VertexCompare::finalize() {
for ( double value : pully_vector ) { ++m_histo_pully[value]; }
for ( double value : pullz_vector ) { ++m_histo_pullz[value]; }
- info() << " stat = " << m_stat_dx.standard_deviation() << endmsg;
- info() << " stat = " << m_stat_dx_sq.sum2() << endmsg;
info() << "dx: "
<< format( "mean = %5.3f +/- %5.3f, RMS = %5.3f +/- %5.3f",
static_cast<double>( m_stat_dx.sum() ) / m_stat_dx.nEntries(), m_stat_dx.meanErr(),
- m_stat_dx.standard_deviation(), rmsErr( m_stat_dx, m_stat_dx_sq ) )
+ m_stat_dx.standard_deviation(), rmsErr( m_stat_dx, m_stat_dx_fourth_sum ) )
<< endmsg;
info() << "dy: "
<< format( "mean = %5.3f +/- %5.3f, RMS = %5.3f +/- %5.3f",
static_cast<double>( m_stat_dy.sum() ) / m_stat_dy.nEntries(), m_stat_dy.meanErr(),
- m_stat_dy.standard_deviation(), rmsErr( m_stat_dy, m_stat_dy_sq ) )
+ m_stat_dy.standard_deviation(), rmsErr( m_stat_dy, m_stat_dy_fourth_sum ) )
<< endmsg;
info() << "dz: "
<< format( "mean = %5.3f +/- %5.3f, RMS = %5.3f +/- %5.3f",
static_cast<double>( m_stat_dz.sum() ) / m_stat_dz.nEntries(), m_stat_dz.meanErr(),
- m_stat_dz.standard_deviation(), rmsErr( m_stat_dz, m_stat_dz_sq ) )
+ m_stat_dz.standard_deviation(), rmsErr( m_stat_dz, m_stat_dz_fourth_sum ) )
<< endmsg;
info() << " ---------------------------------------" << endmsg;
info() << "pullx: "
<< format( "mean = %5.3f +/- %5.3f, RMS = %5.3f +/- %5.3f",
static_cast<double>( m_stat_pullx.sum() ) / m_stat_pullx.nEntries(), m_stat_pullx.meanErr(),
- m_stat_pullx.standard_deviation(), rmsErr( m_stat_pullx, m_stat_pullx_sq ) )
+ m_stat_pullx.standard_deviation(), rmsErr( m_stat_pullx, m_stat_pullx_fourth_sum ) )
<< endmsg;
info() << "pully: "
<< format( "mean = %5.3f +/- %5.3f, RMS = %5.3f +/- %5.3f",
static_cast<double>( m_stat_pully.sum() ) / m_stat_pully.nEntries(), m_stat_pully.meanErr(),
- m_stat_pully.standard_deviation(), rmsErr( m_stat_pully, m_stat_pully_sq ) )
+ m_stat_pully.standard_deviation(), rmsErr( m_stat_pully, m_stat_pully_fourth_sum ) )
<< endmsg;
info() << "pullz: "
<< format( "mean = %5.3f +/- %5.3f, RMS = %5.3f +/- %5.3f",
static_cast<double>( m_stat_pullz.sum() ) / m_stat_pullz.nEntries(), m_stat_pullz.meanErr(),
- m_stat_pullz.standard_deviation(), rmsErr( m_stat_pullz, m_stat_pullz_sq ) )
+ m_stat_pullz.standard_deviation(), rmsErr( m_stat_pullz, m_stat_pullz_fourth_sum ) )
<< endmsg;
info() << " ---------------------------------------" << endmsg;
info() << "diff in x: "
<< format( "%5.3f +/- %5.3f", std::sqrt( 1.0 + m_stat_pullx.standard_deviation() ) - 1.0,
- rmsErr( m_stat_pullx, m_stat_pullx_sq ) * 0.5 /
+ rmsErr( m_stat_pullx, m_stat_pullx_fourth_sum ) * 0.5 /
std::sqrt( 1.0 + m_stat_pullx.standard_deviation() ) )
<< endmsg;
info() << "diff in y: "
<< format( "%5.3f +/- %5.3f", std::sqrt( 1.0 + m_stat_pully.standard_deviation() ) - 1.0,
- rmsErr( m_stat_pully, m_stat_pully_sq ) * 0.5 /
+ rmsErr( m_stat_pully, m_stat_pully_fourth_sum ) * 0.5 /
std::sqrt( 1.0 + m_stat_pully.standard_deviation() ) )
<< endmsg;
info() << "diff in z: "
<< format( "%5.3f +/- %5.3f", std::sqrt( 1.0 + m_stat_pullz.standard_deviation() ) - 1.0,
- rmsErr( m_stat_pullz, m_stat_pullz_sq ) * 0.5 /
+ rmsErr( m_stat_pullz, m_stat_pullz_fourth_sum ) * 0.5 /
std::sqrt( 1.0 + m_stat_pullz.standard_deviation() ) )
<< endmsg;
info() << " ============================================" << endmsg;
--
GitLab
From 8fb0bf863b2267f398625945b6216e7d861ff532 Mon Sep 17 00:00:00 2001
From: Gitlab CI <noreply@cern.ch>
Date: Mon, 27 May 2024 18:19:22 +0000
Subject: [PATCH 18/21] Fixed formatting
patch generated by https://gitlab.cern.ch/lhcb/Rec/-/jobs/39374176
---
Tr/TrackMonitors/src/VertexCompare.cpp | 26 +++++++++++++-------------
1 file changed, 13 insertions(+), 13 deletions(-)
diff --git a/Tr/TrackMonitors/src/VertexCompare.cpp b/Tr/TrackMonitors/src/VertexCompare.cpp
index abf0da1513a..8d248bd774f 100644
--- a/Tr/TrackMonitors/src/VertexCompare.cpp
+++ b/Tr/TrackMonitors/src/VertexCompare.cpp
@@ -188,7 +188,7 @@ double kurtosis( StatAccumulator const& stAcc, FourthSumAccumulator const& stAcc
// get an error in the rms value
// ============================================================================
-template <typename StatAccumulator, typename FourthSumAccumulator >
+template <typename StatAccumulator, typename FourthSumAccumulator>
double rmsErr( StatAccumulator const& stAcc, FourthSumAccumulator const& stAcc_fourth_sum ) {
const auto n = stAcc.nEntries();
if ( 1 >= n ) { return 0.0; }
@@ -227,22 +227,22 @@ private:
Gaudi::Property<bool> m_requireSingle{this, "requireSingle", false};
mutable Gaudi::Accumulators::StatAccumulator<Gaudi::Accumulators::atomicity::full, double> m_stat_dx;
- mutable Gaudi::Accumulators::SumAccumulator<Gaudi::Accumulators::atomicity::full, double> m_stat_dx_fourth_sum;
+ mutable Gaudi::Accumulators::SumAccumulator<Gaudi::Accumulators::atomicity::full, double> m_stat_dx_fourth_sum;
mutable Gaudi::Accumulators::StatAccumulator<Gaudi::Accumulators::atomicity::full, double> m_stat_dy;
- mutable Gaudi::Accumulators::SumAccumulator<Gaudi::Accumulators::atomicity::full, double> m_stat_dy_fourth_sum;
+ mutable Gaudi::Accumulators::SumAccumulator<Gaudi::Accumulators::atomicity::full, double> m_stat_dy_fourth_sum;
mutable Gaudi::Accumulators::StatAccumulator<Gaudi::Accumulators::atomicity::full, double> m_stat_dz;
- mutable Gaudi::Accumulators::SumAccumulator<Gaudi::Accumulators::atomicity::full, double> m_stat_dz_fourth_sum;
+ mutable Gaudi::Accumulators::SumAccumulator<Gaudi::Accumulators::atomicity::full, double> m_stat_dz_fourth_sum;
mutable Gaudi::Accumulators::StatAccumulator<Gaudi::Accumulators::atomicity::full, double> m_stat_pullx;
- mutable Gaudi::Accumulators::SumAccumulator<Gaudi::Accumulators::atomicity::full, double> m_stat_pullx_fourth_sum;
+ mutable Gaudi::Accumulators::SumAccumulator<Gaudi::Accumulators::atomicity::full, double> m_stat_pullx_fourth_sum;
mutable Gaudi::Accumulators::StatAccumulator<Gaudi::Accumulators::atomicity::full, double> m_stat_pully;
- mutable Gaudi::Accumulators::SumAccumulator<Gaudi::Accumulators::atomicity::full, double> m_stat_pully_fourth_sum;
+ mutable Gaudi::Accumulators::SumAccumulator<Gaudi::Accumulators::atomicity::full, double> m_stat_pully_fourth_sum;
mutable Gaudi::Accumulators::StatAccumulator<Gaudi::Accumulators::atomicity::full, double> m_stat_pullz;
- mutable Gaudi::Accumulators::SumAccumulator<Gaudi::Accumulators::atomicity::full, double> m_stat_pullz_fourth_sum;
+ mutable Gaudi::Accumulators::SumAccumulator<Gaudi::Accumulators::atomicity::full, double> m_stat_pullz_fourth_sum;
mutable std::vector<double> dx_vector{}, dy_vector{}, dz_vector{}, pullx_vector{}, pully_vector{}, pullz_vector{};
struct monitoringHistos {
@@ -508,17 +508,17 @@ void VertexCompare::operator()( Vertices const& recoVtx1, Vertices const& recoVt
++m_nTracks_dif[ntracks2 - ntracks1];
m_stat_dx += dx;
- m_stat_dx_fourth_sum += std::pow(dx, 4);
+ m_stat_dx_fourth_sum += std::pow( dx, 4 );
m_stat_dy += dy;
- m_stat_dy_fourth_sum += std::pow(dy, 4);
+ m_stat_dy_fourth_sum += std::pow( dy, 4 );
m_stat_dz += dz;
- m_stat_dz_fourth_sum += std::pow(dz, 4);
+ m_stat_dz_fourth_sum += std::pow( dz, 4 );
m_stat_pullx += pullx;
- m_stat_pullx_fourth_sum += std::pow(pullx, 4);
+ m_stat_pullx_fourth_sum += std::pow( pullx, 4 );
m_stat_pully += pully;
- m_stat_pully_fourth_sum += std::pow(pully, 4);
+ m_stat_pully_fourth_sum += std::pow( pully, 4 );
m_stat_pullz += pullz;
- m_stat_pullz_fourth_sum += std::pow(pullz, 4);
+ m_stat_pullz_fourth_sum += std::pow( pullz, 4 );
} // else end
if ( m_produceNtuple.value() ) {
--
GitLab
From c4e100471e10d1894e9b6254210ecf67c613f444 Mon Sep 17 00:00:00 2001
From: Bogdan Kutsenko <bogdan.kutsenko@cern.ch>
Date: Thu, 30 May 2024 11:00:28 +0200
Subject: [PATCH 19/21] Clang compilation fix
---
Tr/TrackUtils/src/RandomVeloTrackContainerSplitter.cpp | 2 +-
1 file changed, 1 insertion(+), 1 deletion(-)
diff --git a/Tr/TrackUtils/src/RandomVeloTrackContainerSplitter.cpp b/Tr/TrackUtils/src/RandomVeloTrackContainerSplitter.cpp
index cb735e11bdf..15de57ebb5d 100644
--- a/Tr/TrackUtils/src/RandomVeloTrackContainerSplitter.cpp
+++ b/Tr/TrackUtils/src/RandomVeloTrackContainerSplitter.cpp
@@ -19,7 +19,7 @@
* @date 25/05/2024
*/
-#include "Event/Track.h"
+#include "Event/PrVeloTracks.h"
#include "Functors/Function.h"
#include "Functors/with_functors.h"
#include "LHCbAlgs/Transformer.h"
--
GitLab
From 4ec2660e8f3086afcfcb1a425fb73d7f1200d0d3 Mon Sep 17 00:00:00 2001
From: Bogdan Kutsenko <bogdan.kutsenko@cern.ch>
Date: Fri, 14 Jun 2024 18:52:59 +0200
Subject: [PATCH 20/21] Cleaned / titles changed
---
Tr/TrackMonitors/src/VertexCompare.cpp | 85 +++++++++++++-------------
1 file changed, 43 insertions(+), 42 deletions(-)
diff --git a/Tr/TrackMonitors/src/VertexCompare.cpp b/Tr/TrackMonitors/src/VertexCompare.cpp
index 8d248bd774f..5d7e5b8bd4d 100644
--- a/Tr/TrackMonitors/src/VertexCompare.cpp
+++ b/Tr/TrackMonitors/src/VertexCompare.cpp
@@ -219,7 +219,7 @@ public:
private:
bool debugLevel() const { return msgLevel( MSG::DEBUG ) || msgLevel( MSG::VERBOSE ); }
- static constexpr auto ntrack_bins = std::array{3.5, 16.5, 29.5, 43.5, 56.5, 69.5};
+ static constexpr auto ntrack_bins = std::array{3.5, 15.5, 30.5, 45.5, 58.5, 70.5};
Gaudi::Property<bool> m_produceHistogram{this, "produceHistogram", true};
Gaudi::Property<bool> m_produceNtuple{this, "produceNtuple", true};
@@ -244,7 +244,6 @@ private:
mutable Gaudi::Accumulators::StatAccumulator<Gaudi::Accumulators::atomicity::full, double> m_stat_pullz;
mutable Gaudi::Accumulators::SumAccumulator<Gaudi::Accumulators::atomicity::full, double> m_stat_pullz_fourth_sum;
- mutable std::vector<double> dx_vector{}, dy_vector{}, dz_vector{}, pullx_vector{}, pully_vector{}, pullz_vector{};
struct monitoringHistos {
mutable std::array<Gaudi::Accumulators::Histogram<1>, 5> m_histo_nTracksBins_dx;
mutable std::array<Gaudi::Accumulators::Histogram<1>, 5> m_histo_nTracksBins_dy;
@@ -259,35 +258,32 @@ private:
std::index_sequence<IDXs...> ) {
return {{{owner,
name + std::to_string( IDXs ),
- def.title() + std::to_string( IDXs ),
+ def.title() + ", ntracks > " + std::to_string( ntrack_bins[IDXs] ) + " & "+ std::to_string( ntrack_bins[IDXs+1]) +
+ " <ntracks",
{static_cast<unsigned int>( def.bins() ), def.lowEdge(), def.highEdge()}}...}};
}
monitoringHistos( const VertexCompare* owner )
: m_histo_nTracksBins_dx{histo1DArrayBuilder(
- owner, "dx_Monitoring_", Gaudi::Histo1DDef{"dx, mm", -0.15, 0.15, 50}, std::make_index_sequence<5>() )}
+ owner, "dx_Monitoring_ntracks_bin", Gaudi::Histo1DDef{"dx, mm", -0.15, 0.15, 50}, std::make_index_sequence<5>() )}
, m_histo_nTracksBins_dy{histo1DArrayBuilder(
- owner, "dy_Monitoring_", Gaudi::Histo1DDef{"dy, mm", -0.15, 0.15, 50}, std::make_index_sequence<5>() )}
+ owner, "dy_Monitoring_ntracks_bin", Gaudi::Histo1DDef{"dy, mm", -0.15, 0.15, 50}, std::make_index_sequence<5>() )}
, m_histo_nTracksBins_dz{histo1DArrayBuilder(
- owner, "dz_Monitoring_", Gaudi::Histo1DDef{"dz, mm", -1.5, 1.5, 50}, std::make_index_sequence<5>() )}
+ owner, "dz_Monitoring_ntracks_bin", Gaudi::Histo1DDef{"dz, mm", -1.5, 1.5, 50}, std::make_index_sequence<5>() )}
, m_histo_pullx_Monitoring{owner, "pullx_Monitoring", "pull x", {20, -5, 5}}
, m_histo_pully_Monitoring{owner, "pully_Monitoring", "pull y", {20, -5, 5}}
, m_histo_pullz_Monitoring{owner, "pullz_Monitoring", "pull z", {20, -5, 5}} {}
};
std::unique_ptr<monitoringHistos> m_monitoringHistos;
- using axis1D = Gaudi::Accumulators::Axis<Gaudi::Accumulators::Histogram<1>::AxisArithmeticType>;
- mutable Gaudi::Accumulators::RootHistogram<1> m_histo_dx{this, "dx", "dx, mm", axis1D{50, -0.15, 0.15}};
- mutable Gaudi::Accumulators::RootHistogram<1> m_histo_dy{this, "dy", "dy, mm", axis1D{50, -0.15, 0.15}};
- mutable Gaudi::Accumulators::RootHistogram<1> m_histo_dz{this, "dz", "dz, mm", axis1D{50, -1.5, 1.5}};
- mutable Gaudi::Accumulators::RootHistogram<1> m_histo_pullx{this, "pullx", "pull x", axis1D{20, -5, 5}};
- mutable Gaudi::Accumulators::RootHistogram<1> m_histo_pully{this, "pully", "pull y", axis1D{20, -5, 5}};
- mutable Gaudi::Accumulators::RootHistogram<1> m_histo_pullz{this, "pullz", "pull z", axis1D{20, -5, 5}};
- mutable Gaudi::Accumulators::Histogram<1> m_nTracks1{this, "ntracks1", "Number of tracks in PV set 1",
- axis1D{50, 0., 150.}};
- mutable Gaudi::Accumulators::Histogram<1> m_nTracks2{this, "ntracks2", "Number of tracks in PV set 2",
- axis1D{50, 0., 150.}};
- mutable Gaudi::Accumulators::Histogram<1> m_nTracks_dif{
- this, "dtracks", "Difference in the number of tracks in two sets of PVs", axis1D{50, 0., 150.}};
+ mutable std::optional<Gaudi::Accumulators::Histogram<1>> m_nTracks1;
+ mutable std::optional<Gaudi::Accumulators::Histogram<1>> m_nTracks2;
+ mutable std::optional<Gaudi::Accumulators::RootHistogram<1>> m_histo_dx;
+ mutable std::optional<Gaudi::Accumulators::RootHistogram<1>> m_histo_dy;
+ mutable std::optional<Gaudi::Accumulators::RootHistogram<1>> m_histo_dz;
+ mutable std::optional<Gaudi::Accumulators::RootHistogram<1>> m_histo_pullx;
+ mutable std::optional<Gaudi::Accumulators::RootHistogram<1>> m_histo_pully;
+ mutable std::optional<Gaudi::Accumulators::RootHistogram<1>> m_histo_pullz;
+ mutable std::optional<Gaudi::Accumulators::Histogram<1>> m_nTracks_dif;
mutable Gaudi::Accumulators::Counter<> m_nVtx{this, "Number of pairs of vertices in processed events"};
mutable Gaudi::Accumulators::Counter<> m_nRec1{this, "Number of vertices in input1"};
@@ -316,6 +312,23 @@ StatusCode VertexCompare::initialize() {
return Consumer::initialize().andThen( [&]() {
if ( m_monitoring.value() ) m_monitoringHistos = std::make_unique<monitoringHistos>( this );
+ // define range of new histograms from properties
+ using axis1D = Gaudi::Accumulators::Axis<Gaudi::Accumulators::Histogram<1>::AxisArithmeticType>;
+
+ if (m_produceHistogram.value() || m_monitoring.value() ){
+ m_nTracks1.emplace( this, "ntracks_set1", "Number of tracks in PV set 1", axis1D{50, 0., 150.} );
+ m_nTracks2.emplace( this, "ntracks_set2", "Number of tracks in PV set 2", axis1D{50, 0., 150.} );
+ }
+ if(m_produceHistogram.value() ){
+ m_histo_dx.emplace( this, "dx", "dx, mm", axis1D{50, -0.15, 0.15} );
+ m_histo_dy.emplace( this, "dy", "dy, mm", axis1D{50, -0.15, 0.15} );
+ m_histo_dz.emplace( this, "dz", "dz, mm", axis1D{50, -1.5, 1.5} );
+ m_histo_pullx.emplace( this, "pullx", "pull x", axis1D{20, -5, 5} );
+ m_histo_pully.emplace( this, "pully", "pull y", axis1D{20, -5, 5} );
+ m_histo_pullz.emplace( this, "pullz", "pull z", axis1D{20, -5, 5} );
+ m_nTracks_dif.emplace(
+ this, "dtracks", "Difference in the number of tracks in two sets of PVs", axis1D{50, 0., 150.});
+ }
m_nVtx.reset();
m_nEvt.reset();
LHCb::Conditions::InteractionRegion::addConditionDerivation( this,
@@ -490,22 +503,18 @@ void VertexCompare::operator()( Vertices const& recoVtx1, Vertices const& recoVt
++monitoringHistos.m_histo_pullz_Monitoring[pullz];
}
}
+ if(m_produceHistogram.value() || m_monitoring.value() ){
+ ++m_nTracks1.value()[ntracks1];
+ ++m_nTracks2.value()[ntracks2];
+ }
if ( m_produceHistogram.value() ) {
- dx_vector.push_back( dx );
- dy_vector.push_back( dy );
- dz_vector.push_back( dz );
- pullx_vector.push_back( pullx );
- pully_vector.push_back( pully );
- pullz_vector.push_back( pullz );
- ++m_histo_dx[dx];
- ++m_histo_dy[dy];
- ++m_histo_dz[dz];
- ++m_histo_pullx[pullx];
- ++m_histo_pully[pully];
- ++m_histo_pullz[pullz];
- ++m_nTracks1[ntracks1];
- ++m_nTracks2[ntracks2];
- ++m_nTracks_dif[ntracks2 - ntracks1];
+ ++m_histo_dx.value()[dx];
+ ++m_histo_dy.value()[dy];
+ ++m_histo_dz.value()[dz];
+ ++m_histo_pullx.value()[pullx];
+ ++m_histo_pully.value()[pully];
+ ++m_histo_pullz.value()[pullz];
+ ++m_nTracks_dif.value()[ntracks2 - ntracks1];
m_stat_dx += dx;
m_stat_dx_fourth_sum += std::pow( dx, 4 );
@@ -724,14 +733,6 @@ StatusCode VertexCompare::finalize() {
<< " % " << endmsg;
if ( m_produceHistogram.value() ) {
-
- for ( double value : dx_vector ) { ++m_histo_dx[value]; }
- for ( double value : dy_vector ) { ++m_histo_dy[value]; }
- for ( double value : dz_vector ) { ++m_histo_dz[value]; }
- for ( double value : pullx_vector ) { ++m_histo_pullx[value]; }
- for ( double value : pully_vector ) { ++m_histo_pully[value]; }
- for ( double value : pullz_vector ) { ++m_histo_pullz[value]; }
-
info() << "dx: "
<< format( "mean = %5.3f +/- %5.3f, RMS = %5.3f +/- %5.3f",
static_cast<double>( m_stat_dx.sum() ) / m_stat_dx.nEntries(), m_stat_dx.meanErr(),
--
GitLab
From ea803a373a785512fdfb53754221d44694b3d5c6 Mon Sep 17 00:00:00 2001
From: Gitlab CI <noreply@cern.ch>
Date: Fri, 14 Jun 2024 17:03:50 +0000
Subject: [PATCH 21/21] Fixed formatting
patch generated by https://gitlab.cern.ch/lhcb/Rec/-/jobs/40081826
---
Tr/TrackMonitors/src/VertexCompare.cpp | 43 ++++++++++++++------------
1 file changed, 23 insertions(+), 20 deletions(-)
diff --git a/Tr/TrackMonitors/src/VertexCompare.cpp b/Tr/TrackMonitors/src/VertexCompare.cpp
index 5d7e5b8bd4d..644ec904e7d 100644
--- a/Tr/TrackMonitors/src/VertexCompare.cpp
+++ b/Tr/TrackMonitors/src/VertexCompare.cpp
@@ -258,17 +258,20 @@ private:
std::index_sequence<IDXs...> ) {
return {{{owner,
name + std::to_string( IDXs ),
- def.title() + ", ntracks > " + std::to_string( ntrack_bins[IDXs] ) + " & "+ std::to_string( ntrack_bins[IDXs+1]) +
- " <ntracks",
+ def.title() + ", ntracks > " + std::to_string( ntrack_bins[IDXs] ) + " & " +
+ std::to_string( ntrack_bins[IDXs + 1] ) + " <ntracks",
{static_cast<unsigned int>( def.bins() ), def.lowEdge(), def.highEdge()}}...}};
}
monitoringHistos( const VertexCompare* owner )
- : m_histo_nTracksBins_dx{histo1DArrayBuilder(
- owner, "dx_Monitoring_ntracks_bin", Gaudi::Histo1DDef{"dx, mm", -0.15, 0.15, 50}, std::make_index_sequence<5>() )}
- , m_histo_nTracksBins_dy{histo1DArrayBuilder(
- owner, "dy_Monitoring_ntracks_bin", Gaudi::Histo1DDef{"dy, mm", -0.15, 0.15, 50}, std::make_index_sequence<5>() )}
- , m_histo_nTracksBins_dz{histo1DArrayBuilder(
- owner, "dz_Monitoring_ntracks_bin", Gaudi::Histo1DDef{"dz, mm", -1.5, 1.5, 50}, std::make_index_sequence<5>() )}
+ : m_histo_nTracksBins_dx{histo1DArrayBuilder( owner, "dx_Monitoring_ntracks_bin",
+ Gaudi::Histo1DDef{"dx, mm", -0.15, 0.15, 50},
+ std::make_index_sequence<5>() )}
+ , m_histo_nTracksBins_dy{histo1DArrayBuilder( owner, "dy_Monitoring_ntracks_bin",
+ Gaudi::Histo1DDef{"dy, mm", -0.15, 0.15, 50},
+ std::make_index_sequence<5>() )}
+ , m_histo_nTracksBins_dz{histo1DArrayBuilder( owner, "dz_Monitoring_ntracks_bin",
+ Gaudi::Histo1DDef{"dz, mm", -1.5, 1.5, 50},
+ std::make_index_sequence<5>() )}
, m_histo_pullx_Monitoring{owner, "pullx_Monitoring", "pull x", {20, -5, 5}}
, m_histo_pully_Monitoring{owner, "pully_Monitoring", "pull y", {20, -5, 5}}
, m_histo_pullz_Monitoring{owner, "pullz_Monitoring", "pull z", {20, -5, 5}} {}
@@ -315,20 +318,20 @@ StatusCode VertexCompare::initialize() {
// define range of new histograms from properties
using axis1D = Gaudi::Accumulators::Axis<Gaudi::Accumulators::Histogram<1>::AxisArithmeticType>;
- if (m_produceHistogram.value() || m_monitoring.value() ){
+ if ( m_produceHistogram.value() || m_monitoring.value() ) {
m_nTracks1.emplace( this, "ntracks_set1", "Number of tracks in PV set 1", axis1D{50, 0., 150.} );
m_nTracks2.emplace( this, "ntracks_set2", "Number of tracks in PV set 2", axis1D{50, 0., 150.} );
}
- if(m_produceHistogram.value() ){
- m_histo_dx.emplace( this, "dx", "dx, mm", axis1D{50, -0.15, 0.15} );
- m_histo_dy.emplace( this, "dy", "dy, mm", axis1D{50, -0.15, 0.15} );
- m_histo_dz.emplace( this, "dz", "dz, mm", axis1D{50, -1.5, 1.5} );
- m_histo_pullx.emplace( this, "pullx", "pull x", axis1D{20, -5, 5} );
- m_histo_pully.emplace( this, "pully", "pull y", axis1D{20, -5, 5} );
- m_histo_pullz.emplace( this, "pullz", "pull z", axis1D{20, -5, 5} );
- m_nTracks_dif.emplace(
- this, "dtracks", "Difference in the number of tracks in two sets of PVs", axis1D{50, 0., 150.});
- }
+ if ( m_produceHistogram.value() ) {
+ m_histo_dx.emplace( this, "dx", "dx, mm", axis1D{50, -0.15, 0.15} );
+ m_histo_dy.emplace( this, "dy", "dy, mm", axis1D{50, -0.15, 0.15} );
+ m_histo_dz.emplace( this, "dz", "dz, mm", axis1D{50, -1.5, 1.5} );
+ m_histo_pullx.emplace( this, "pullx", "pull x", axis1D{20, -5, 5} );
+ m_histo_pully.emplace( this, "pully", "pull y", axis1D{20, -5, 5} );
+ m_histo_pullz.emplace( this, "pullz", "pull z", axis1D{20, -5, 5} );
+ m_nTracks_dif.emplace( this, "dtracks", "Difference in the number of tracks in two sets of PVs",
+ axis1D{50, 0., 150.} );
+ }
m_nVtx.reset();
m_nEvt.reset();
LHCb::Conditions::InteractionRegion::addConditionDerivation( this,
@@ -503,7 +506,7 @@ void VertexCompare::operator()( Vertices const& recoVtx1, Vertices const& recoVt
++monitoringHistos.m_histo_pullz_Monitoring[pullz];
}
}
- if(m_produceHistogram.value() || m_monitoring.value() ){
+ if ( m_produceHistogram.value() || m_monitoring.value() ) {
++m_nTracks1.value()[ntracks1];
++m_nTracks2.value()[ntracks2];
}
--
GitLab