Draft: Add PV splitting for vertex resolution studies (primarily BGI)

Adds PVSplitter which takes as input a set of VELO tracks, reconstructs PVs in the event and then splits and re-fits each PV using a specified algorithm for dividing the tracks into two subsets. Also modifies LumiPVs_nobeamline to tuple split-vertex information. Addresses #183.

Example options to test in Moore:

from Moore import options
from RecoConf.decoders import default_VeloCluster_source
from RecoConf.legacy_rec_hlt1_tracking import all_velo_track_types
from RecoConf.config import Reconstruction, run_reconstruction
from PyConf.application import make_odin, default_raw_banks, default_raw_event
from PyConf.Algorithms import HltRoutingBitsFilter, LumiPVs_nobeamline, PVSplitter
from PyConf.reading import get_decreports

default_raw_event.global_bind(stream="beamgas")

def make_reconstruction():
    phys_filter = HltRoutingBitsFilter(
        RawBanks=default_raw_banks('HltRoutingBits'),
        RequireMask=(0x0, 0x0, 0x80000000)
        )
    algs = []
    for splitalg in ['random', 'byvelohalf', 'equalvelohalf', 'bytopbottom', 'equaltopbottom']:
        pvsplit = PVSplitter(
            InputVeloTracks=all_velo_track_types()['v1'],
            SplitMethod=splitalg
            )
        tuple_maker = LumiPVs_nobeamline(
            name=f'PVs_{splitalg}',
            ODIN=make_odin(),
            InputPVs=pvsplit.OutputPVs,
            InputSplitPVs1=pvsplit.OutputSplitPVs1,
            InputSplitPVs2=pvsplit.OutputSplitPVs2,
            InputHlt=get_decreports('Hlt1')
            )
        algs += [pvsplit, tuple_maker]
    return Reconstruction("BeamGasTuple", algs, [phys_filter])

options.input_files = ['root://xrootd.echo.stfc.ac.uk//lhcb:prod/lhcb/LHCb/Vdm24/BEAMGAS.DST/00242987/0000/00242987_00005402_1.beamgas.dst']
options.ntuple_file = 'split_pvs.root'
options.input_type = 'ROOT'
options.evt_max = 100
options.geometry_version = 'run3/2024.Q1.2-v00.00'
options.conditions_version = 'master'
options.simulation = False

with default_VeloCluster_source.bind(bank_type="VPRetinaCluster"):
    run_reconstruction(options, make_reconstruction)

To do before un-drafting:

• Get things building against master
• Test PVSplitter runs and produces some output
• Test PVSplitter produces sensible output. Probably the easiest way to do this is to run over a reasonable sample of BGI vertices, and check the resolution de-convolved bunch width is approximately stable along z/with changing ntracks
• Add reco options in separate MR to Moore, alongside/to replace pvreco_beamgasstream_tupling_vdm_2024.py?

cc Lumi conveners @efranzos @vazhovko @elniel @fferrari

Rec/RecAlgs/src/LumiPVs_nobeamline.cpp

@@ -21,64 +21,94 @@
 #include "DetDesc/DetectorElement.h"
 #include "TrackInterfaces/IPVOfflineTool.h"
 
-using InputVeloTracks = LHCb::Event::v1::Tracks;
-using InputHlt        = LHCb::HltDecReports;
-using RecPV           = LHCb::RecVertex;
+using InputPVs = LHCb::RecVertices;
+using InputHlt = LHCb::HltDecReports;
+using RecPV    = LHCb::RecVertex;
+
+namespace {
+  int count_backwards_tracks( RecPV* pv ) {
+    int nbw = 0;
+    for ( auto tr : pv->tracks() ) {
+      if ( tr->isVeloBackward() ) nbw++;
+    }
+    return nbw;
+  }
+
+  // convert from LHCb::RecVertices -> std::vector<LHCb::RecVertex*>
+  void convert( const InputPVs& oldvec, std::vector<RecPV*>& newvec ) {
+    std::copy( oldvec.begin(), oldvec.end(), std::back_inserter( newvec ) );
+  }
+} // namespace
+
 class LumiPVs_nobeamline
-    : public virtual LHCb::Algorithm::Consumer<
-          void( const LHCb::ODIN&, const InputVeloTracks&, const InputHlt&, const DetectorElement& ),
-          LHCb::Algorithm::Traits::usesBaseAndConditions<GaudiTupleAlg, DetectorElement>> {
+    : public virtual LHCb::Algorithm::Consumer<void( const LHCb::ODIN&, const InputPVs&, const InputPVs&,
+                                                     const InputPVs&, const InputHlt& ),
+                                               LHCb::Algorithm::Traits::usesBaseAndConditions<GaudiTupleAlg>> {
+
 private:
   Gaudi::Property<std::string> m_tupleName{ this, "BeamGasTuple", "BeamGasTuple" };
-  ToolHandle<IPVOfflineTool>   m_toolpv{ this, "PVOfflineTool", "PVOfflineTool" };
-  LHCb::RecVertices            reconstructMultiPVFromTracks( std::vector<const LHCb::Event::v1::Track*>& tracks2use,
-                                                             Gaudi::XYZPoint const&                      beamSpot ) const;
 
 public:
-  virtual void operator()( const LHCb::ODIN&, const InputVeloTracks&, const InputHlt&,
-                           const DetectorElement& ) const override;
+  virtual void operator()( const LHCb::ODIN&, const InputPVs&, const InputPVs&, const InputPVs&,
+                           const InputHlt& ) const override;
   LumiPVs_nobeamline( const std::string& name, ISvcLocator* pSvcLocator );
 };
 
 DECLARE_COMPONENT_WITH_ID( LumiPVs_nobeamline, "LumiPVs_nobeamline" );
 LumiPVs_nobeamline::LumiPVs_nobeamline( const std::string& name, ISvcLocator* pSvcLocator )
     : Consumer( name, pSvcLocator,
-                { KeyValue{ "ODIN", "" }, KeyValue{ "InputVeloTracks", "" }, KeyValue{ "InputHlt", "" },
-                  KeyValue{ "StandardGeometryTop", LHCb::standard_geometry_top } } ) {}
-
-void LumiPVs_nobeamline::operator()( const LHCb::ODIN& odin, const InputVeloTracks& tracks, const InputHlt& decreport,
-                                     const DetectorElement& lhcb ) const {
+                { KeyValue{ "ODIN", "" }, KeyValue{ "InputPVs", "" }, KeyValue{ "InputSplitPVs1", "" },
+                  KeyValue{ "InputSplitPVs2", "" }, KeyValue{ "InputHlt", "" } } ) {}
 
+void LumiPVs_nobeamline::operator()( const LHCb::ODIN& odin, const InputPVs& _pvs, const InputPVs& _split1,
+                                     const InputPVs& _split2, const InputHlt& decreport ) const {
   auto tuple = nTuple( m_tupleName );
 
-  std::vector<const LHCb::Track*> selectedtracks;
-  std::vector<RecPV>              vtxlist;
-  std::copy_if( tracks.begin(), tracks.end(), std::back_inserter( selectedtracks ),
-                []( const LHCb::Track* tr ) { return tr->hasVelo(); } );
+  // convert to vectors, as I can't find a way to simultaneously iterate over KeyedContainers
+  std::vector<RecPV*> pvs, split1, split2;
+  convert( _pvs, pvs );
+  convert( _split1, split1 );
+  convert( _split2, split2 );
+
+  StatusCode  sc;
+  std::string suffix;
+  RecPV*      pv;
+  for ( int i = 0; i < static_cast<int>( pvs.size() ); i++ ) {
+    sc = tuple->column( "runNumber", odin.runNumber() );
+    sc &= tuple->column( "gpsTime", static_cast<unsigned long long>( odin.gpsTime() ) );
+    sc &= tuple->column( "eventNumber", static_cast<unsigned long long>( odin.eventNumber() ) );
+    sc &= tuple->column( "eventType", odin.eventType() );
+    sc &= tuple->column( "triggerType", odin.triggerType() );
+    sc &= tuple->column( "orbitNumber", odin.orbitNumber() );
+    sc &= tuple->column( "bxType", static_cast<std::uint16_t>( odin.bunchCrossingType() ) );
+    sc &= tuple->column( "bxId", odin.bunchId() );
+    sc &= tuple->column( "step", odin.calibrationStep() );
+    for ( const auto& kvp : decreport ) { sc &= tuple->column( kvp.first, kvp.second.decision() ); }
 
-  m_toolpv->reconstructMultiPVFromTracks( selectedtracks, vtxlist, *lhcb.geometry() ).ignore();
+    for ( int j = 0; j < 3; j++ ) {
+      if ( j == 0 ) {
+        suffix = "";
+        pv     = pvs.at( i );
+      } else if ( j == 1 ) {
+        suffix = "_1";
+        pv     = split1.at( i );
+      } else {
+        suffix = "_2";
+        pv     = split2.at( i );
+      }
 
-  auto sc = tuple->column( "runNumber", odin.runNumber() );
-  sc      = tuple->column( "gpsTime", static_cast<unsigned long long>( odin.gpsTime() ) );
-  sc      = tuple->column( "eventNumber", static_cast<unsigned long long>( odin.eventNumber() ) );
-  sc      = tuple->column( "eventType", odin.eventType() );
-  sc      = tuple->column( "triggerType", odin.triggerType() );
-  sc      = tuple->column( "orbitNumber", odin.orbitNumber() );
-  sc      = tuple->column( "bxType", static_cast<std::uint16_t>( odin.bunchCrossingType() ) );
-  sc      = tuple->column( "bxId", odin.bunchId() );
-  sc      = tuple->column( "step", odin.calibrationStep() );
-  sc      = tuple->column( "nVeloTracks", static_cast<unsigned long long>( tracks.size() ) );
-  for ( const auto& kvp : decreport ) { sc = tuple->column( kvp.first, kvp.second.decision() ); }
+      sc &= tuple->column( "PV_ntracks" + suffix, static_cast<unsigned long long>( pv->tracks().size() ) );
+      sc &= tuple->column( "PV_ntracksbw" + suffix, count_backwards_tracks( pv ) );
+      sc &= tuple->column( "PV_chi2" + suffix, pv->chi2() );
+      sc &= tuple->column( "PV_ndof" + suffix, pv->nDoF() );
+      sc &= tuple->column( "PVX" + suffix, pv->position().x() );
+      sc &= tuple->column( "PVY" + suffix, pv->position().y() );
+      sc &= tuple->column( "PVZ" + suffix, pv->position().z() );
+    }
 
-  sc = tuple->farray(
-      { { "PV_nTracks", +[]( const RecPV& pv ) -> double { return pv.tracks().size(); } },
-        { "PV_chi2ndof", +[]( const RecPV& pv ) { return pv.chi2() / pv.nDoF(); } },
-        { "PVX", +[]( const RecPV& pv ) { return pv.position().x(); } },
-        { "PVY", +[]( const RecPV& pv ) { return pv.position().y(); } },
-        { "PVZ", +[]( const RecPV& pv ) { return pv.position().z(); } },
-        { "PVR",
-          +[]( const RecPV& pv ) { return sqrt( pow( pv.position().x(), 2 ) + pow( pv.position().y(), 2 ) ); } } },
-      vtxlist.begin(), vtxlist.end(), "nPVs", 1000 );
+    if ( sc.isFailure() ) continue;
 
-  sc.andThen( [&] { return tuple->write(); } ).orThrow( "Failed to fill ntuple", "LumiPV_nobeamline Alg" );
+    sc = tuple->write();
+    if ( sc.isFailure() ) throw GaudiException( "Cannot fill ntuple", "LumiPVs_nobeamline", StatusCode::FAILURE );
+  }
 }