diff --git a/Phys/FunctorCore/include/Functors/Functional.h b/Phys/FunctorCore/include/Functors/Functional.h index fa32bbb3a431c47a040362dea638dc6657582811..3fc49e6faf7801c18ac5093c06b274211c492d8b 100644 --- a/Phys/FunctorCore/include/Functors/Functional.h +++ b/Phys/FunctorCore/include/Functors/Functional.h @@ -252,7 +252,7 @@ namespace Functors::Functional { TrivialFunctor{ "Min", []( Range const& range ) { using std::ranges::begin; using std::ranges::end; - // FIXME doc. + // reduce to smallest value return std::accumulate( begin( range ), end( range ), std::ranges::range_value_t{ std::numeric_limits::max() }, []( auto curr_min_value, auto value ) { @@ -285,7 +285,7 @@ namespace Functors::Functional { constexpr auto Max = TrivialFunctor{ "Max", []( Range const& range ) { using std::ranges::begin; using std::ranges::end; - // FIXME doc. + // reduce to largest value return std::accumulate( begin( range ), end( range ), std::ranges::range_value_t{ std::numeric_limits::lowest() }, @@ -295,6 +295,29 @@ namespace Functors::Functional { } ); } }; + /* * + * @brief return nth element of range + * + * @param index of element to return + * */ + struct ElementAtIndex final : public Function { + + constexpr ElementAtIndex( int idx ) : m_idx( idx ) {} + + auto name() const { return "ElementAtIndex( " + std::to_string( m_idx ) + " )"; } + + auto operator()( std::ranges::range auto const& range ) const { + using std::ranges::begin; + using std::ranges::end; + using std::ranges::next; + return ( m_idx < std::ranges::size( range ) ) ? Functors::Optional{ *next( begin( range ), m_idx ) } + : std::nullopt; + } + + private: + unsigned long m_idx{ 0 }; + }; + /* * * @brief return sum element of range * */ diff --git a/Phys/FunctorCore/python/Functors/__init__.py b/Phys/FunctorCore/python/Functors/__init__.py index be1e5d9635bf7a63ec128143ad35bb8a704d9559..9a6f0ea1bac8c1718b91bdbbd93255c05ee670e9 100644 --- a/Phys/FunctorCore/python/Functors/__init__.py +++ b/Phys/FunctorCore/python/Functors/__init__.py @@ -3625,6 +3625,13 @@ MIN_ELEMENT_NOTZERO = Functor( MAX_ELEMENT = Functor("MAX_ELEMENT", "Functional::Max", "Maximum value of range") +ELEMENT_AT_INDEX = Functor( + "ELEMENT_AT_INDEX", + "Functional::ElementAtIndex", + "Element at given index of range", + Params=[("Index", "The index of the element to return", int)], +) + ENTRY_WITH_MIN_REL_VALUE_OF = Functor( "ENTRY_WITH_MIN_REL_VALUE_OF", "Functional::MinElement", diff --git a/Phys/FunctorCore/tests/src/TestFunctors.cpp b/Phys/FunctorCore/tests/src/TestFunctors.cpp index 199edb2b7b37e15d68271146f86559bbd978017a..91764f5bb4f1648e134adb0cb9b4be1b7945e724 100644 --- a/Phys/FunctorCore/tests/src/TestFunctors.cpp +++ b/Phys/FunctorCore/tests/src/TestFunctors.cpp @@ -290,7 +290,7 @@ BOOST_AUTO_TEST_CASE( test_rich_dll ) { auto const RICH_DLL_BT = Functors::PID::RichDLLbt; // define particle and protoparticle - auto part = LHCb::Particle( LHCb::ParticleID{ 11 } ); + auto part = LHCb::Particle( LHCb::ParticleIDs::electron ); auto proto = LHCb::ProtoParticle(); // add DLL info @@ -641,7 +641,7 @@ BOOST_AUTO_TEST_CASE( test_brem_functors ) { eproto1.setTrack( &etrack1 ); eproto1.setBremInfo( &ebrem1 ); Gaudi::LorentzVector emom1 = { 0., 0., 10. * Gaudi::Units::GeV, 10. * Gaudi::Units::GeV }; - auto electron1 = LHCb::Particle( LHCb::ParticleID{ 11 } ); + auto electron1 = LHCb::Particle( LHCb::ParticleIDs::electron ); electron1.setMomentum( emom1 ); electron1.setProto( &eproto1 ); BOOST_CHECK_EQUAL( threeMomentum( BREM( electron1 ) ).z().cast(), 12. * Gaudi::Units::GeV ); @@ -655,11 +655,11 @@ BOOST_AUTO_TEST_CASE( test_brem_functors ) { eproto2.setTrack( &etrack2 ); eproto2.setBremInfo( &ebrem2 ); Gaudi::LorentzVector emom2 = { 0., 0., -9. * Gaudi::Units::GeV, 9. * Gaudi::Units::GeV }; - auto electron2 = LHCb::Particle( LHCb::ParticleID{ -11 } ); + auto electron2 = LHCb::Particle( LHCb::ParticleIDs::positron ); electron2.setMomentum( emom2 ); electron2.setProto( &eproto2 ); - auto dielectron = LHCb::Particle( LHCb::ParticleID{ 443 } ); + auto dielectron = LHCb::Particle( LHCb::ParticleIDs::jpsi ); dielectron.addToDaughters( &electron1 ); dielectron.addToDaughters( &electron2 ); dielectron.setMomentum( electron1.momentum() + electron2.momentum() ); @@ -687,7 +687,7 @@ BOOST_AUTO_TEST_CASE( test_muonPID_functors ) { BOOST_AUTO_TEST_CASE( test_reconstructible_functor ) { // Build MCProperty - const auto electron = LHCb::ParticleID{ 11 }; + const auto electron = LHCb::ParticleIDs::electron; LHCb::MCProperty mc_property{}; LHCb::MCParticles mc_particles{}; mc_particles.insert( new LHCb::MCParticle{}, 1 ); @@ -1180,13 +1180,13 @@ BOOST_AUTO_TEST_CASE( test_v2_particle_functors ) { // Create some particles. @todo Need to be adapted to v2 LHCb::Particle::Container make_parts() { LHCb::Particle::Container parts; - auto p0 = std::make_unique( LHCb::ParticleID( 431 ), 0 ); // creator with PID and key - auto p1 = std::make_unique( LHCb::ParticleID( 431 ), 1 ); // creator with PID and key - auto p2 = std::make_unique( LHCb::ParticleID( 431 ), 2 ); // creator with PID and key - p0->setMomentum( Gaudi::LorentzVector( 100., 100., 1000., 5000. ) ); // 141 MeV pt - p1->setMomentum( Gaudi::LorentzVector( 0., 0., 1000., 5000. ) ); // 0 MeV pt - p2->setMomentum( Gaudi::LorentzVector( 1000., 0., 1000., 5000. ) ); // 1000 MeV pt - parts.insert( p0.release(), 0 ); // creator with PID and key + auto p0 = std::make_unique( LHCb::ParticleIDs::d_s_plus, 0 ); // creator with PID and key + auto p1 = std::make_unique( LHCb::ParticleIDs::d_s_plus, 1 ); // creator with PID and key + auto p2 = std::make_unique( LHCb::ParticleIDs::d_s_plus, 2 ); // creator with PID and key + p0->setMomentum( Gaudi::LorentzVector( 100., 100., 1000., 5000. ) ); // 141 MeV pt + p1->setMomentum( Gaudi::LorentzVector( 0., 0., 1000., 5000. ) ); // 0 MeV pt + p2->setMomentum( Gaudi::LorentzVector( 1000., 0., 1000., 5000. ) ); // 1000 MeV pt + parts.insert( p0.release(), 0 ); // creator with PID and key parts.insert( p1.release(), 1 ); parts.insert( p2.release(), 2 ); return parts; @@ -1292,10 +1292,8 @@ BOOST_AUTO_TEST_CASE( test_relation_table_isolation ) { // Create vector of fake particles auto parts = make_parts(); - float fake_w0_min = 0.25; - float fake_w1_min = 1.15; - float fake_w0_max = 0.35; - float fake_w1_max = 1.25; + float exp_w0 = 0.3; + float exp_w1 = 1.2; float fake_sum_pt_min = 1100.f; float fake_sum_pt_max = 1800.f; float fake_max_pt_min = 980.f; @@ -1304,7 +1302,7 @@ BOOST_AUTO_TEST_CASE( test_relation_table_isolation ) { float fake_min_pt_max = 150.f; float fake_asym_pt = -1.; - auto ref_part = std::make_unique( LHCb::ParticleID( 431 ), 0 ); + auto ref_part = std::make_unique( LHCb::ParticleIDs::d_s_plus, 0 ); ref_part->setMomentum( Gaudi::LorentzVector( 0., 0., 1000., 5000. ) ); auto vertex = std::make_unique(); @@ -1331,10 +1329,11 @@ BOOST_AUTO_TEST_CASE( test_relation_table_isolation ) { auto w1 = back_weight_prepped( mask_arg, true, table, parts( 2 ) ); BOOST_CHECK( w0.has_value() ); BOOST_CHECK( w1.has_value() ); - BOOST_CHECK( LHCb::Utils::as_arithmetic( w0.value() ) > fake_w0_min && - LHCb::Utils::as_arithmetic( w0.value() ) < fake_w0_max ); - BOOST_CHECK( LHCb::Utils::as_arithmetic( w1.value() ) > fake_w1_min && - LHCb::Utils::as_arithmetic( w1.value() ) < fake_w1_max ); + BOOST_CHECK( LHCb::Utils::as_arithmetic( w0.value() ) > exp_w0 - 0.01 && + LHCb::Utils::as_arithmetic( w0.value() ) < exp_w0 + 0.01 ); + BOOST_CHECK( LHCb::Utils::as_arithmetic( w1.value() ) > exp_w1 - 0.01 && + LHCb::Utils::as_arithmetic( w1.value() ) < exp_w1 + 0.01 ); + // the relation table is already sorted by the weight auto const pt_to_min_weight_entry = Functors::chain( PT, Functors::Common::To, Functors::Functional::MinElement( Functors::Common::Weight ), @@ -1521,7 +1520,7 @@ BOOST_AUTO_TEST_CASE( test_cov_matrix_functors ) { Gaudi::LorentzVector v1_p4{ 1., 2., 3., 4. }; // Create v1 particle - LHCb::Particle v1_particle{ LHCb::ParticleID( 511 ), 0 }; + LHCb::Particle v1_particle{ LHCb::ParticleIDs::b_zero, 0 }; v1_particle.setMomCovMatrix( v1_momCovMatrix ); v1_particle.setPosCovMatrix( v1_posCovMatrix ); v1_particle.setPosMomCovMatrix( v1_posMomCovMatrix ); @@ -1586,9 +1585,9 @@ BOOST_AUTO_TEST_CASE( test_cov_matrix_functors ) { } BOOST_AUTO_TEST_CASE( test_child ) { - auto B0 = LHCb::Particle{ LHCb::ParticleID( 511 ), 0 }; - auto K1 = LHCb::Particle( LHCb::ParticleID( 321 ), 1 ); - auto K2 = LHCb::Particle( LHCb::ParticleID( 321 ), 2 ); + auto B0 = LHCb::Particle{ LHCb::ParticleIDs::b_zero, 0 }; + auto K1 = LHCb::Particle( LHCb::ParticleIDs::kaon_plus, 1 ); + auto K2 = LHCb::Particle( LHCb::ParticleIDs::kaon_minus, 2 ); SmartRefVector daugs; daugs.push_back( &K1 ); @@ -1605,19 +1604,19 @@ BOOST_AUTO_TEST_CASE( test_vertex_diff_functors ) { const double delta_z{ 5.0 }; // B0 particle - auto B0_part = LHCb::Particle{ LHCb::ParticleID( 511 ), 0 }; + auto B0_part = LHCb::Particle{ LHCb::ParticleIDs::b_zero, 0 }; // set endvertex of B0 auto B0_endvertex = std::make_unique(); B0_endvertex->setPosition( Gaudi::XYZPoint{ 0.0, 0.0, B0_vertex_z } ); B0_part.setEndVertex( B0_endvertex.get() ); // Ds- particle with vertex and a basic K+ - auto Ds_part = LHCb::Particle( LHCb::ParticleID( -432 ), 1 ); + auto Ds_part = LHCb::Particle( LHCb::ParticleIDs::d_s_minus, 1 ); // set endvertex of Ds auto Ds_endvertex = std::make_unique(); Ds_endvertex->setPosition( Gaudi::XYZPoint{ 0.0, 0.0, B0_vertex_z + delta_z } ); Ds_part.setEndVertex( Ds_endvertex.get() ); - auto K_part = LHCb::Particle( LHCb::ParticleID( 321 ), 2 ); + auto K_part = LHCb::Particle( LHCb::ParticleIDs::kaon_plus, 2 ); // set Ds+ and K- as daughters of B0 SmartRefVector daugs; @@ -1640,7 +1639,7 @@ BOOST_AUTO_TEST_CASE( test_vertex_diff_functors ) { } BOOST_AUTO_TEST_CASE( test_decaytime_functor ) { - auto part = LHCb::Particle{ LHCb::ParticleID( 511 ), 0 }; + auto part = LHCb::Particle{ LHCb::ParticleIDs::b_zero, 0 }; // set endvertex of B0 auto decayvertex = LHCb::Vertex{}; @@ -1700,12 +1699,12 @@ BOOST_AUTO_TEST_CASE( test_mc_prompt ) { } BOOST_AUTO_TEST_CASE( test_mcorr_functors ) { - auto B0_parr = LHCb::Particle{ LHCb::ParticleID( 511 ), 0 }; + auto B0_parr = LHCb::Particle{ LHCb::ParticleIDs::b_zero, 0 }; auto decayvertex = LHCb::Vertex{}; decayvertex.setPosition( Gaudi::XYZPoint{ 0.0, 0.0, 5.0 } ); B0_parr.setEndVertex( &decayvertex ); - auto B0_perp = LHCb::Particle{ LHCb::ParticleID( 511 ), 0 }; + auto B0_perp = LHCb::Particle{ LHCb::ParticleIDs::b_zero, 0 }; B0_perp.setEndVertex( &decayvertex ); const double pz = 1.e5; @@ -1777,9 +1776,9 @@ BOOST_AUTO_TEST_CASE( test_mc_functors ) { auto mc_vertices = LHCb::MCVertices{}; // define properties for an mcp - const auto pid = LHCb::ParticleID{ 211 }; - const auto electron = LHCb::ParticleID{ 11 }; - const auto muon = LHCb::ParticleID{ 12 }; + const auto pid = LHCb::ParticleIDs::pion_plus; + const auto electron = LHCb::ParticleIDs::electron; + const auto muon = LHCb::ParticleIDs::muon_minus; const auto mother_key = 1; const auto gd_mother_key = 2; const auto child_key = 42; @@ -2171,27 +2170,41 @@ BOOST_AUTO_TEST_CASE( test_functional_functors ) { BOOST_CHECK( ( first_ip > 0.81649 && first_ip < 0.81650 ) ); BOOST_CHECK_EQUAL( tmp[1].cast(), 0 ); - auto min = Functors::Functional::Min; - auto max = Functors::Functional::Max; + auto min = Functors::Functional::Min; + auto max = Functors::Functional::Max; + constexpr auto first = Functors::Functional::ElementAtIndex{ 0 }; + constexpr auto second = Functors::Functional::ElementAtIndex{ 1 }; + constexpr auto third = Functors::Functional::ElementAtIndex{ 2 }; + constexpr auto tenth = Functors::Functional::ElementAtIndex{ 9 }; // outside of range std::vector vec{ 1.1, 2.2, 3.3, 1.5, 0.9, 1.0 }; - float output_min = min( vec ); - float output_max = max( vec ); + float output_min = min( vec ); + float output_max = max( vec ); + float output_first = first( vec ).value(); + float output_second = second( vec ).value(); + float output_third = third( vec ).value(); - float desired_min = 0.9; - float desired_max = 3.3; + float desired_min = 0.9; + float desired_max = 3.3; + float desired_first = 1.1; + float desired_second = 2.2; + float desired_third = 3.3; BOOST_CHECK_EQUAL( output_min, desired_min ); BOOST_CHECK_EQUAL( output_max, desired_max ); + BOOST_CHECK_EQUAL( output_first, desired_first ); + BOOST_CHECK_EQUAL( output_second, desired_second ); + BOOST_CHECK_EQUAL( output_third, desired_third ); + BOOST_CHECK_EQUAL( tenth( vec ).has_value(), false ); } BOOST_AUTO_TEST_CASE( test_ft_functors ) { // define three different B particles - LHCb::Particle particleA( LHCb::ParticleID( 511 ), 0 ); // B0 - LHCb::Particle particleB( LHCb::ParticleID( -511 ), 0 ); // B0bar - LHCb::Particle particleC( LHCb::ParticleID( 521 ), 0 ); // B+ + LHCb::Particle particleA( LHCb::ParticleIDs::b_zero, 0 ); // B0 + LHCb::Particle particleB( LHCb::ParticleIDs::anti_b_zero, 0 ); // B0bar + LHCb::Particle particleC( LHCb::ParticleIDs::b_plus, 0 ); // B+ Gaudi::LorentzVector momentumA( 1000, 1000, 9000, 10530 ); Gaudi::XYZPoint referencePointA( 0.75, 0.07, 25 ); @@ -2371,16 +2384,16 @@ BOOST_AUTO_TEST_CASE( check_binary_functors ) { BOOST_AUTO_TEST_CASE( test_flatten_functor ) { // First, make a couple of particles - auto B_plus = LHCb::Particle( LHCb::ParticleID( 521 ), 0 ); - auto B0 = LHCb::Particle( LHCb::ParticleID( 511 ), 0 ); - auto pi_plus = LHCb::Particle( LHCb::ParticleID( 211 ), 0 ); - auto pi_minus = LHCb::Particle( LHCb::ParticleID( -211 ), 0 ); - auto K0 = LHCb::Particle( LHCb::ParticleID( 311 ), 0 ); - auto K_plus = LHCb::Particle( LHCb::ParticleID( 321 ), 0 ); - auto K_minus = LHCb::Particle( LHCb::ParticleID( -321 ), 0 ); - auto Bs0 = LHCb::Particle( LHCb::ParticleID( 531 ), 0 ); - auto Ds_plus = LHCb::Particle( LHCb::ParticleID( 431 ), 0 ); - auto D_minus = LHCb::Particle( LHCb::ParticleID( -411 ), 0 ); + auto B_plus = LHCb::Particle( LHCb::ParticleIDs::b_plus, 0 ); + auto B0 = LHCb::Particle( LHCb::ParticleIDs::b_zero, 0 ); + auto pi_plus = LHCb::Particle( LHCb::ParticleIDs::pion_plus, 0 ); + auto pi_minus = LHCb::Particle( LHCb::ParticleIDs::pion_minus, 0 ); + auto K0 = LHCb::Particle( LHCb::ParticleIDs::kaon_neutral, 0 ); + auto K_plus = LHCb::Particle( LHCb::ParticleIDs::kaon_plus, 0 ); + auto K_minus = LHCb::Particle( LHCb::ParticleIDs::kaon_minus, 0 ); + auto Bs0 = LHCb::Particle( LHCb::ParticleIDs::b_s_zero, 0 ); + auto Ds_plus = LHCb::Particle( LHCb::ParticleIDs::d_s_plus, 0 ); + auto D_minus = LHCb::Particle( LHCb::ParticleIDs::d_minus, 0 ); // First decay: Bs0 -> ( B+ -> K0 pi+ ) pi- (DUMMY) // set pi- and K+ as daughters of B0 @@ -2475,49 +2488,70 @@ BOOST_AUTO_TEST_CASE( test_flatten_functor ) { // The logic is as such: Look at the Bs0, is it basic? No -> Look at the first // daughter, which is the B0. Is it basic? No. Look at the first daughter of the B0, // which is the pi-. Is it basic? Yes, then put it into the vector. - BOOST_CHECK_EQUAL( bs_basics[0]->particleID().pid(), 211 ); // is the first entry a pi+ - BOOST_CHECK_EQUAL( bs_basics[1]->particleID().pid(), 311 ); // is the second entry a K0 - BOOST_CHECK_EQUAL( bs_basics[2]->particleID().pid(), -211 ); // is the third entry a pi- - BOOST_CHECK_EQUAL( bs_descendants[0]->particleID().pid(), 521 ); // is the first entry a B+ - BOOST_CHECK_EQUAL( bs_descendants[1]->particleID().pid(), 211 ); // is the second entry a pi+ - BOOST_CHECK_EQUAL( bs_descendants[2]->particleID().pid(), 311 ); // is the third entry a K0 - BOOST_CHECK_EQUAL( bs_descendants[3]->particleID().pid(), -211 ); // is the fourth entry a pi- - BOOST_CHECK_EQUAL( bs_children[0]->particleID().pid(), 521 ); // is the first entry a B+ - BOOST_CHECK_EQUAL( bs_children[1]->particleID().pid(), -211 ); // is the second entry a pi- - BOOST_CHECK_EQUAL( bs_grandchildren[0]->particleID().pid(), 211 ); // is the first entry a pi+ - BOOST_CHECK_EQUAL( bs_grandchildren[1]->particleID().pid(), 311 ); // is the second entry a K0 - - BOOST_CHECK_EQUAL( b0_basics[0]->particleID().pid(), 321 ); // is the first entry a K+ - BOOST_CHECK_EQUAL( b0_basics[1]->particleID().pid(), -321 ); // is the second entry a K- - BOOST_CHECK_EQUAL( b0_basics[2]->particleID().pid(), 211 ); // is the third entry a pi+ - BOOST_CHECK_EQUAL( b0_basics[3]->particleID().pid(), 321 ); // is the first entry a K+ - BOOST_CHECK_EQUAL( b0_basics[4]->particleID().pid(), -211 ); // is the second entry a pi- - BOOST_CHECK_EQUAL( b0_basics[5]->particleID().pid(), -211 ); // is the third entry a pi- - BOOST_CHECK_EQUAL( b0_descendants[0]->particleID().pid(), 431 ); // is the first entry a Ds - BOOST_CHECK_EQUAL( b0_descendants[1]->particleID().pid(), 321 ); // is the second entry a K+ - BOOST_CHECK_EQUAL( b0_descendants[2]->particleID().pid(), -321 ); // is the third entry a K- - BOOST_CHECK_EQUAL( b0_descendants[3]->particleID().pid(), 211 ); // is the fourth entry a pi+ - BOOST_CHECK_EQUAL( b0_descendants[4]->particleID().pid(), -411 ); // is the first entry a D- - BOOST_CHECK_EQUAL( b0_descendants[5]->particleID().pid(), 321 ); // is the second entry a K+ - BOOST_CHECK_EQUAL( b0_descendants[6]->particleID().pid(), -211 ); // is the third entry a pi- - BOOST_CHECK_EQUAL( b0_descendants[7]->particleID().pid(), -211 ); // is the fourth entry a pi- - BOOST_CHECK_EQUAL( b0_children[0]->particleID().pid(), 431 ); // is the first entry a Ds+ - BOOST_CHECK_EQUAL( b0_children[1]->particleID().pid(), -411 ); // is the second entry a D- + BOOST_CHECK_EQUAL( bs_basics[0]->particleID(), LHCb::ParticleIDs::pion_plus ); // is the first entry a pi+ + BOOST_CHECK_EQUAL( bs_basics[1]->particleID(), + LHCb::ParticleIDs::kaon_neutral ); // is the second entry a K0 + BOOST_CHECK_EQUAL( bs_basics[2]->particleID(), + LHCb::ParticleIDs::pion_minus ); // is the third entry a pi- + BOOST_CHECK_EQUAL( bs_descendants[0]->particleID(), + LHCb::ParticleIDs::b_plus ); // is the first entry a B+ + BOOST_CHECK_EQUAL( bs_descendants[1]->particleID(), + LHCb::ParticleIDs::pion_plus ); // is the second entry a pi+ + BOOST_CHECK_EQUAL( bs_descendants[2]->particleID(), + LHCb::ParticleIDs::kaon_neutral ); // is the third entry a K0 + BOOST_CHECK_EQUAL( bs_descendants[3]->particleID(), + LHCb::ParticleIDs::pion_minus ); // is the fourth entry a pi- + BOOST_CHECK_EQUAL( bs_children[0]->particleID(), LHCb::ParticleIDs::b_plus ); // is the first entry a B+ + BOOST_CHECK_EQUAL( bs_children[1]->particleID(), + LHCb::ParticleIDs::pion_minus ); // is the second entry a pi- + BOOST_CHECK_EQUAL( bs_grandchildren[0]->particleID(), + LHCb::ParticleIDs::pion_plus ); // is the first entry a pi+ + BOOST_CHECK_EQUAL( bs_grandchildren[1]->particleID(), + LHCb::ParticleIDs::kaon_neutral ); // is the second entry a K0 + + BOOST_CHECK_EQUAL( b0_basics[0]->particleID(), LHCb::ParticleIDs::kaon_plus ); // is the first entry a K+ + BOOST_CHECK_EQUAL( b0_basics[1]->particleID(), + LHCb::ParticleIDs::kaon_minus ); // is the second entry a K- + BOOST_CHECK_EQUAL( b0_basics[2]->particleID(), LHCb::ParticleIDs::pion_plus ); // is the third entry a pi+ + BOOST_CHECK_EQUAL( b0_basics[3]->particleID(), LHCb::ParticleIDs::kaon_plus ); // is the first entry a K+ + BOOST_CHECK_EQUAL( b0_basics[4]->particleID(), + LHCb::ParticleIDs::pion_minus ); // is the second entry a pi- + BOOST_CHECK_EQUAL( b0_basics[5]->particleID(), + LHCb::ParticleIDs::pion_minus ); // is the third entry a pi- + BOOST_CHECK_EQUAL( b0_descendants[0]->particleID(), + LHCb::ParticleIDs::d_s_plus ); // is the first entry a Ds + BOOST_CHECK_EQUAL( b0_descendants[1]->particleID(), + LHCb::ParticleIDs::kaon_plus ); // is the second entry a K+ + BOOST_CHECK_EQUAL( b0_descendants[2]->particleID(), + LHCb::ParticleIDs::kaon_minus ); // is the third entry a K- + BOOST_CHECK_EQUAL( b0_descendants[3]->particleID(), + LHCb::ParticleIDs::pion_plus ); // is the fourth entry a pi+ + BOOST_CHECK_EQUAL( b0_descendants[4]->particleID(), + LHCb::ParticleIDs::d_minus ); // is the first entry a D- + BOOST_CHECK_EQUAL( b0_descendants[5]->particleID(), + LHCb::ParticleIDs::kaon_plus ); // is the second entry a K+ + BOOST_CHECK_EQUAL( b0_descendants[6]->particleID(), + LHCb::ParticleIDs::pion_minus ); // is the third entry a pi- + BOOST_CHECK_EQUAL( b0_descendants[7]->particleID(), + LHCb::ParticleIDs::pion_minus ); // is the fourth entry a pi- + BOOST_CHECK_EQUAL( b0_children[0]->particleID(), + LHCb::ParticleIDs::d_s_plus ); // is the first entry a Ds+ + BOOST_CHECK_EQUAL( b0_children[1]->particleID(), LHCb::ParticleIDs::d_minus ); // is the second entry a D- } BOOST_AUTO_TEST_CASE( test_quark_content_functor ) { // test to check quark content of particles auto B_plus = LHCb::MCParticle(); - auto B_plus_id = LHCb::ParticleID( 521 ); + auto B_plus_id = LHCb::ParticleIDs::b_plus; B_plus.setParticleID( B_plus_id ); auto Ds_plus = LHCb::MCParticle(); - auto Ds_plus_id = LHCb::ParticleID( 431 ); + auto Ds_plus_id = LHCb::ParticleIDs::d_s_plus; Ds_plus.setParticleID( Ds_plus_id ); auto muon = LHCb::MCParticle(); - auto muon_id = LHCb::ParticleID( 13 ); + auto muon_id = LHCb::ParticleIDs::muon_minus; muon.setParticleID( muon_id ); auto proton = LHCb::MCParticle(); - auto proton_id = LHCb::ParticleID( 2212 ); + auto proton_id = LHCb::ParticleIDs::proton; proton.setParticleID( proton_id ); const auto HAS_BOTTOM = Functors::Simulation::HasQuark( "b" ); @@ -2529,7 +2563,7 @@ BOOST_AUTO_TEST_CASE( test_quark_content_functor ) { const auto PARTICLE_ID_OBJ = Functors::Simulation::Particle_Id_Obj; const auto B_plus_id_obj = PARTICLE_ID_OBJ( B_plus ); - BOOST_CHECK_EQUAL( B_plus_id_obj.pid(), 521 ); + BOOST_CHECK_EQUAL( B_plus_id_obj.pid(), B_plus_id.pid() ); BOOST_CHECK_EQUAL( HAS_BOTTOM( B_plus ), true ); BOOST_CHECK_EQUAL( HAS_CHARM( B_plus ), false ); @@ -3002,13 +3036,13 @@ BOOST_AUTO_TEST_CASE( test_tree_functors ) { Functors::TopLevelInfo const top_level; // First, make a couple of particles - auto B0 = LHCb::Particle{ LHCb::ParticleID( 511 ), 0 }; - auto pi_plus = LHCb::Particle( LHCb::ParticleID( 211 ), 0 ); - auto pi_minus = LHCb::Particle( LHCb::ParticleID( -211 ), 0 ); - auto K_plus = LHCb::Particle( LHCb::ParticleID( 321 ), 0 ); - auto K_minus = LHCb::Particle( LHCb::ParticleID( -321 ), 0 ); - auto Ds_plus = LHCb::Particle( LHCb::ParticleID( 431 ), 0 ); - auto D_minus = LHCb::Particle( LHCb::ParticleID( -411 ), 0 ); + auto B0 = LHCb::Particle{ LHCb::ParticleIDs::b_zero, 0 }; + auto pi_plus = LHCb::Particle( LHCb::ParticleIDs::pion_plus, 0 ); + auto pi_minus = LHCb::Particle( LHCb::ParticleIDs::pion_minus, 0 ); + auto K_plus = LHCb::Particle( LHCb::ParticleIDs::kaon_plus, 0 ); + auto K_minus = LHCb::Particle( LHCb::ParticleIDs::kaon_minus, 0 ); + auto Ds_plus = LHCb::Particle( LHCb::ParticleIDs::d_s_plus, 0 ); + auto D_minus = LHCb::Particle( LHCb::ParticleIDs::d_minus, 0 ); // Decay: B0 -> ( Ds+ -> K+ K- pi+ ) ( D- -> K+ pi- pi- ) @@ -3133,14 +3167,14 @@ BOOST_AUTO_TEST_CASE( test_tree_functors ) { auto children = CHILD( B0 ); auto grandchildren = GCHILD( B0 ); - BOOST_CHECK_EQUAL( children[0]->particleID().pid(), 431 ); // is the first entry a Ds+ - BOOST_CHECK_EQUAL( children[1]->particleID().pid(), -411 ); // is the second entry a D- - BOOST_CHECK_EQUAL( grandchildren[0]->particleID().pid(), 321 ); // is the first entry a K+ - BOOST_CHECK_EQUAL( grandchildren[1]->particleID().pid(), -321 ); // is the second entry a K- - BOOST_CHECK_EQUAL( grandchildren[2]->particleID().pid(), 211 ); // is the third entry a pi+ - BOOST_CHECK_EQUAL( grandchildren[3]->particleID().pid(), 321 ); // is the first entry a K+ - BOOST_CHECK_EQUAL( grandchildren[4]->particleID().pid(), -211 ); // is the second entry a pi- - BOOST_CHECK_EQUAL( grandchildren[5]->particleID().pid(), -211 ); // is the third entry a pi- + BOOST_CHECK_EQUAL( children[0]->particleID(), LHCb::ParticleIDs::d_s_plus ); // is the first entry a Ds+ + BOOST_CHECK_EQUAL( children[1]->particleID(), LHCb::ParticleIDs::d_minus ); // is the second entry a D- + BOOST_CHECK_EQUAL( grandchildren[0]->particleID(), LHCb::ParticleIDs::kaon_plus ); // is the first entry a K+ + BOOST_CHECK_EQUAL( grandchildren[1]->particleID(), LHCb::ParticleIDs::kaon_minus ); // is the second entry a K- + BOOST_CHECK_EQUAL( grandchildren[2]->particleID(), LHCb::ParticleIDs::pion_plus ); // is the third entry a pi+ + BOOST_CHECK_EQUAL( grandchildren[3]->particleID(), LHCb::ParticleIDs::kaon_plus ); // is the first entry a K+ + BOOST_CHECK_EQUAL( grandchildren[4]->particleID(), LHCb::ParticleIDs::pion_minus ); // is the second entry a pi- + BOOST_CHECK_EQUAL( grandchildren[5]->particleID(), LHCb::ParticleIDs::pion_minus ); // is the third entry a pi- auto const CHILD_NINGENERATION = Functors::chain( ::Functors::detail::SizeOf, @@ -3338,9 +3372,9 @@ BOOST_AUTO_TEST_CASE( test_lhcb_math_functors ) { BOOST_AUTO_TEST_CASE( test_require_close ) { // make some dummy particles - auto B0 = LHCb::Particle{ LHCb::ParticleID( 511 ), 0 }; - auto pi_plus = LHCb::Particle( LHCb::ParticleID( 211 ), 0 ); - auto pi_minus = LHCb::Particle( LHCb::ParticleID( -211 ), 0 ); + auto B0 = LHCb::Particle{ LHCb::ParticleIDs::b_zero, 0 }; + auto pi_plus = LHCb::Particle( LHCb::ParticleIDs::pion_plus, 0 ); + auto pi_minus = LHCb::Particle( LHCb::ParticleIDs::pion_minus, 0 ); pi_plus.setMomentum( Gaudi::LorentzVector( 0., 0., 1000., 5000. ) ); pi_minus.setMomentum( Gaudi::LorentzVector( 0., 0., 1000., 5000. ) ); // set B0 daughters diff --git a/Phys/RelatedInfoTools/src/WeightedRelationTable.cpp b/Phys/RelatedInfoTools/src/WeightedRelationTable.cpp index a3404c2546f439897e84b5bd3c0b84c679958a59..b4ce34aec13758491df05cc30aadcd63b54927aa 100644 --- a/Phys/RelatedInfoTools/src/WeightedRelationTable.cpp +++ b/Phys/RelatedInfoTools/src/WeightedRelationTable.cpp @@ -11,16 +11,21 @@ // Include files #include "Event/Particle.h" +#include "Event/ParticleUtils.h" #include "Functors/with_functors.h" #include "Kernel/IParticlePropertySvc.h" +#include "Kernel/ParticleIDs.h" #include "Kernel/ParticleProperty.h" #include "LHCbAlgs/Transformer.h" #include "Relations/Relation1D.h" #include "Relations/RelationWeighted1D.h" namespace { - using WeightedRelationTable = LHCb::RelationWeighted1D; - using RelationTable = LHCb::Relation1D; + template ::value || + std::is_same::value>::type> + using WeightedRelationTable = LHCb::RelationWeighted1D; + using RelationTable = LHCb::Relation1D; + using WeightedRelationTables = std::tuple, WeightedRelationTable>; using RelTableTransformer = LHCb::Algorithm::Transformer; @@ -28,19 +33,26 @@ namespace { using Transformer = LHCb::Algorithm::Transformer; using WeightedSelectionTransformer = - LHCb::Algorithm::Transformer; + LHCb::Algorithm::Transformer& )>; - using weight_t = bool; + using WeightedRelTableTransformer = LHCb::Algorithm::MultiTransformer; struct PredFunct { constexpr static auto PropertyName = "Cut"; - using Signature = weight_t( LHCb::Particle const&, LHCb::Particle const& ); + using Signature = bool( LHCb::Particle const&, LHCb::Particle const& ); }; struct Funct { constexpr static auto PropertyName = "Functor"; using Signature = const LHCb::Particle::ConstVector( LHCb::Particle const& ); }; + + struct MVAFunct { + constexpr static auto PropertyName = "MVAFunctor"; + using Signature = SIMDWrapper::scalar::float_v( LHCb::Particle const& ); + }; + const SmartRefVector s_empty_calohypos{}; } // namespace @@ -92,7 +104,7 @@ struct WeightedRelTableAlg final : with_functors auto const& fun = this->template getFunctor(); for ( const auto& cand : cands_in_cone ) { for ( const auto& ref_part : ref_parts ) { - const weight_t& weight = fun( *ref_part, *cand ); + const bool& weight = fun( *ref_part, *cand ); if ( weight ) { iso_table.relate( ref_part, cand ).ignore(); ++m_outCount; @@ -150,22 +162,32 @@ struct Projection { // calo-track or track-calo -- in which case we order tracks first return lhs_t != nullptr; } + + bool operator==( Projection const& rhs ) const { + return ( ( track() && ( track() == rhs.track() ) ) || ( neutralPID() && ( neutralPID() == rhs.neutralPID() ) ) ) && + pid() == rhs.pid(); + } }; namespace { LHCb::Particle::Selection selectionFromProjectionVector( std::vector& wparts ) { std::sort( wparts.begin(), wparts.end() ); - auto end = std::unique( wparts.begin(), wparts.end(), []( const auto& lhs, const auto& rhs ) { - return ( ( lhs.track() && ( lhs.track() == rhs.track() ) ) || - ( lhs.neutralPID() && ( lhs.neutralPID() == rhs.neutralPID() ) ) ) && - lhs.pid() == rhs.pid(); - } ); + auto end = std::unique( wparts.begin(), wparts.end() ); LHCb::Particle::Selection output; std::transform( wparts.begin(), end, std::back_inserter( output ), []( const auto& p ) { return p.p; } ); return output; } + + bool isSameProjectionVector( std::vector& lhs, std::vector& rhs ) { + if ( lhs.size() != rhs.size() ) return false; + + std::sort( lhs.begin(), lhs.end() ); + std::sort( rhs.begin(), rhs.end() ); + + return std::equal( lhs.begin(), lhs.end(), rhs.begin() ); + } } // namespace // ============================================================================ @@ -211,23 +233,27 @@ struct SelectionFromWeightedRelationTable final : WeightedSelectionTransformer { SelectionFromWeightedRelationTable( const std::string& name, ISvcLocator* svcLoc ) : Transformer( name, svcLoc, KeyValue( "InputRelations", "" ), KeyValue( "OutputLocation", "" ) ) {} - LHCb::Particle::Selection operator()( const WeightedRelationTable& input ) const override { + LHCb::Particle::Selection operator()( const WeightedRelationTable& input ) const override { std::vector wparts; const auto& rels = input.relations(); wparts.reserve( rels.size() ); - // possibly some transform_if would be useful, but this is probably fine. for ( const auto& rel : rels ) { - if ( rel.weight() > m_minWeight.value() && rel.weight() < m_maxWeight.value() ) wparts.emplace_back( rel.to() ); + if ( ( std::isnan( m_minWeight.value() ) || rel.weight() > m_minWeight.value() ) && + ( std::isnan( m_maxWeight.value() ) || rel.weight() < m_maxWeight.value() ) ) { + wparts.emplace_back( rel.to() ); + } } return selectionFromProjectionVector( wparts ); } - Gaudi::Property m_minWeight{ this, "MinimumWeight", -1.0, "Minimum Weight" }; + Gaudi::Property m_minWeight{ this, "MinimumWeight", std::numeric_limits::quiet_NaN(), + "Minimum Weight" }; Gaudi::Property m_maxWeight{ this, "MaximumWeight", 1000.0, "Maximum Weight" }; }; +; // ============================================================================ /** @class ThOrParticleSelection @@ -269,7 +295,87 @@ private: mutable Gaudi::Accumulators::SummingCounter<> m_outCount{ this, "#OutputParticles" }; }; +/** @class ParticleToParticleWeightedRelTableAlg + * Algorithm that relates signal candidates with extra particles associating as weight the outcome of the MVA functor + * + * @param SignalCandidates Location of signal particles + * @param ExtraParticles Location of extra particles + * @param MVAFunctor MVA functor + * @param Predicate Cut applied to the composite to select candidates + * @returns OutputRelations location of relation table between signal candidates and extra particles + */ + +struct ParticleToParticleWeightedRelTableAlg final : with_functors { + /** the standard constructor + * @param name algorithm instance name + * @param pSvc service locator + */ + ParticleToParticleWeightedRelTableAlg( const std::string& name, ISvcLocator* pSvc ) + : with_functors::with_functors( + name, pSvc, { KeyValue{ "CompositeSelection", "" }, KeyValue{ "SignalSelection", "" } }, + { KeyValue{ "BOutputRelations", "" }, KeyValue{ "BstOutputRelations", "" } } ) {} + + // ========================================================================== + /// the standard execution of the algorithm + // ========================================================================== + + WeightedRelationTables operator()( const LHCb::Particle::Range& parents, + const LHCb::Particle::Range& signals ) const override { + + WeightedRelationTable b_iso_table, bst_iso_table; + auto const& fun = this->template getFunctor(); + + for ( const auto& p : parents ) { + // Evaluate the MVA score for parent + float mva_output = LHCb::Utils::as_arithmetic( fun( *p ) ); + + const LHCb::Particle* b_child = nullptr; + const LHCb::Particle* extra_track_child = nullptr; + // Loop over children + for ( const auto* c : p->daughtersVector() ) { + if ( c->particleID().abspid() == LHCb::ParticleIDs::pion_plus.abspid() ) { // Get pions from Bst + bst_iso_table.relate( p, c, mva_output ).ignore(); + extra_track_child = c; + ++m_nBstRelations; + } else { + b_child = c; + } + } + // b_child and signal do not have the same key, but the basics are the same + // Comparing projections... + for ( const auto& signal : signals ) { + std::vector signal_projections, b_child_projections; + + // Collect projections for signal + for ( const auto* p : LHCb::Utils::Particle::walk( *signal ) ) { + if ( p->isBasicParticle() ) { signal_projections.emplace_back( Projection{ p } ); } + }; + + // Collect projections for b_child + for ( const auto* p : LHCb::Utils::Particle::walk( *b_child ) ) { + if ( p->isBasicParticle() ) { b_child_projections.emplace_back( Projection{ p } ); } + }; + + // Compare projections + if ( isSameProjectionVector( signal_projections, b_child_projections ) ) { + b_iso_table.relate( signal, extra_track_child, mva_output ).ignore(); + ++m_nBRelations; + break; // Ensure that we have only one B for each Bstar + } + } + } + + return std::tuple{ b_iso_table, bst_iso_table }; + } + +private: + mutable Gaudi::Accumulators::Counter<> m_nBRelations{ this, "Number of relations for signal candidates" }; + mutable Gaudi::Accumulators::Counter<> m_nBstRelations{ + this, "Number of relations for combinations of signal and extra pion" }; +}; + DECLARE_COMPONENT( WeightedRelTableAlg ) DECLARE_COMPONENT( SelectionFromRelationTable ) DECLARE_COMPONENT( SelectionFromWeightedRelationTable ) DECLARE_COMPONENT( ThOrParticleSelection ) +DECLARE_COMPONENT( ParticleToParticleWeightedRelTableAlg )