diff --git a/Tr/TrackMonitors/src/VertexCompare.cpp b/Tr/TrackMonitors/src/VertexCompare.cpp
index 26f052b7cff8e6775b25dffbb04628bc89ab0112..644ec904e7d3f7fd59ce10315e9e97b79d1c53d6 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"
@@ -22,6 +21,7 @@
#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 "GaudiUtils/HistoStats.h"
@@ -29,9 +29,10 @@
#include "MCInterfaces/IForcedBDecayTool.h"
#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;
@@ -110,7 +111,7 @@ 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;
@@ -151,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 );
@@ -171,6 +172,32 @@ VtxVariables SetVtxVariables( const RecPVInfo<VertexType>& vrtf, const int& coun
return vtx;
}
+// ============================================================================
+// get the error in kurtosis
+// ============================================================================
+
+template <typename StatAccumulator, typename FourthSumAccumulator>
+double kurtosis( StatAccumulator const& stAcc, FourthSumAccumulator const& stAcc_fourth_sum ) {
+ const auto n = stAcc.nEntries();
+ 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 );
+}
+
+// ============================================================================
+// get an error in the rms value
+// ============================================================================
+
+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_fourth_sum );
+ result += 2.0 / ( n - 1 );
+ result /= 4.0 * n;
+ return stAcc.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>> {
@@ -190,9 +217,77 @@ public:
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 ); }
+
+ 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};
+ Gaudi::Property<bool> m_monitoring{this, "monitoring", false};
+ 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::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::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::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::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::StatAccumulator<Gaudi::Accumulators::atomicity::full, double> m_stat_pullz;
+ mutable Gaudi::Accumulators::SumAccumulator<Gaudi::Accumulators::atomicity::full, double> m_stat_pullz_fourth_sum;
+
+ 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, Gaudi::Histo1DDef def,
+ 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",
+ {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_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;
+
+ 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"};
mutable Gaudi::Accumulators::Counter<> m_nRec2{this, "Number of vertices in input2"};
@@ -217,7 +312,26 @@ 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 );
+ // 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,
@@ -237,10 +351,15 @@ void VertexCompare::operator()( Vertices const& recoVtx1, Vertices const& recoVt
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;
+
for ( auto const& pv : recoVtx1 ) {
- RecPVInfo<VertexType> recinfo1;
++m_nRec1;
recinfo1.pRECPV = &pv;
recinfo1.position = pv.position();
@@ -259,7 +378,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();
@@ -283,9 +401,8 @@ 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;
- int size_diff = fullVrt1.size() - fullVrt2.size();
+ 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();
@@ -295,28 +412,28 @@ void VertexCompare::operator()( Vertices const& recoVtx1, Vertices const& recoVt
}
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 multiplicity = -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 ); }
@@ -365,54 +482,56 @@ 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_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 );
- }
- }
+ 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;
+ for ( size_t i = 1; i < ntrack_bins.size(); ++i ) {
+ if ( std::lround( ( ntracks1 + ntracks2 ) / 2 ) <= ntrack_bins[i] ) { break; }
+ binCount++;
}
+ 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[pullx];
+ ++monitoringHistos.m_histo_pully_Monitoring[pully];
+ ++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() ) {
+ ++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 );
+ m_stat_dy += dy;
+ m_stat_dy_fourth_sum += std::pow( dy, 4 );
+ m_stat_dz += dz;
+ m_stat_dz_fourth_sum += std::pow( dz, 4 );
+ m_stat_pullx += pullx;
+ m_stat_pullx_fourth_sum += std::pow( pullx, 4 );
+ m_stat_pully += pully;
+ m_stat_pully_fourth_sum += std::pow( pully, 4 );
+ m_stat_pullz += pullz;
+ m_stat_pullz_fourth_sum += std::pow( pullz, 4 );
+ } // else end
if ( m_produceNtuple.value() ) {
Tuple myTuple = nTuple( "PV_nTuple", "PV_nTuple", CLID_ColumnWiseTuple );
@@ -475,7 +594,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++;
@@ -522,7 +641,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();
@@ -570,7 +689,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();
@@ -616,74 +735,63 @@ 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: "
- << 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 ) )
+ if ( m_produceHistogram.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(),
+ m_stat_dx.standard_deviation(), rmsErr( m_stat_dx, m_stat_dx_fourth_sum ) )
<< endmsg;
- }
- 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() << "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_fourth_sum ) )
<< 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() << "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_fourth_sum ) )
<< endmsg;
- }
- info() << " ---------------------------------------" << endmsg;
- if ( pullx ) {
+ info() << " ---------------------------------------" << endmsg;
+
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 ) )
+ << 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_fourth_sum ) )
<< 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 ) )
+ << 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_fourth_sum ) )
<< 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 ) )
+ << 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_fourth_sum ) )
<< endmsg;
- }
- info() << " ---------------------------------------" << endmsg;
- if ( pullx ) {
+
+ info() << " ---------------------------------------" << endmsg;
+
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() ) )
+ << format( "%5.3f +/- %5.3f", std::sqrt( 1.0 + m_stat_pullx.standard_deviation() ) - 1.0,
+ rmsErr( m_stat_pullx, m_stat_pullx_fourth_sum ) * 0.5 /
+ std::sqrt( 1.0 + m_stat_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() ) )
+ << format( "%5.3f +/- %5.3f", std::sqrt( 1.0 + m_stat_pully.standard_deviation() ) - 1.0,
+ rmsErr( m_stat_pully, m_stat_pully_fourth_sum ) * 0.5 /
+ std::sqrt( 1.0 + m_stat_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() ) )
+ << format( "%5.3f +/- %5.3f", std::sqrt( 1.0 + m_stat_pullz.standard_deviation() ) - 1.0,
+ rmsErr( m_stat_pullz, m_stat_pullz_fourth_sum ) * 0.5 /
+ std::sqrt( 1.0 + m_stat_pullz.standard_deviation() ) )
<< endmsg;
+ info() << " ============================================" << endmsg;
}
- info() << " ============================================" << endmsg;
-
return Consumer::finalize(); // Must be called after all other actions
}
@@ -742,4 +850,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
+}
diff --git a/Tr/TrackUtils/CMakeLists.txt b/Tr/TrackUtils/CMakeLists.txt
index 35e937a382f3db5c596e961fef14ef23b47923df..bb5b59d5b91ea2fa42b35132af468e8f319fdae3 100644
--- a/Tr/TrackUtils/CMakeLists.txt
+++ b/Tr/TrackUtils/CMakeLists.txt
@@ -26,6 +26,8 @@ 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
src/TrackListMerger.cpp
diff --git a/Tr/TrackUtils/src/RandomTrackContainerSplitter.cpp b/Tr/TrackUtils/src/RandomTrackContainerSplitter.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..a40cc8844d7e1c7d5afde7c879229a96927c649e
--- /dev/null
+++ b/Tr/TrackUtils/src/RandomTrackContainerSplitter.cpp
@@ -0,0 +1,67 @@
+/*****************************************************************************\
+* (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
+ * @author Bogdan Kutsenko
+ * @date 27/05/2020
+ */
+
+#include "Event/Track.h"
+#include "Functors/Function.h"
+#include "Functors/with_functors.h"
+#include "LHCbAlgs/Transformer.h"
+
+namespace {
+ /// The output data
+ using OutConts = std::tuple<LHCb::Track::Selection, LHCb::Track::Selection>;
+
+ struct TrackPredicate {
+ using Signature = bool( const LHCb::Track& );
+ static constexpr auto PropertyName = "Code";
+ };
+} // namespace
+
+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::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::Track::Range& input_tracks ) const override {
+ OutConts out;
+ auto& [passed, rejected] = out;
+ m_inputs += input_tracks.size();
+ for ( auto& trk : input_tracks ) {
+ auto& container = ( makePesudoRandomDecision( *trk ) ? passed : rejected );
+ container.insert( 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 )
\ 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 0000000000000000000000000000000000000000..15de57ebb5ddc0017ab1b84d24666ad79192bd92
--- /dev/null
+++ b/Tr/TrackUtils/src/RandomVeloTrackContainerSplitter.cpp
@@ -0,0 +1,70 @@
+/*****************************************************************************\
+* (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/PrVeloTracks.h"
+#include "Functors/Function.h"
+#include "Functors/with_functors.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 makePseudoRandomDecision( const PrTrackProxy& track ) {
+ const auto lhcb_id = track.lhcbIDs()[0];
+ return static_cast<int>( lhcb_id.lhcbID() ) % 2;
+}
+
+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;
+ 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();
+ 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 )