1) Add example algorithm for MultiSecVertexTool

2) Add TRT hit cut for MultiSecVertexTool

Reconstruction/VKalVrt/VrtSecInclusive/VrtSecInclusive/AllBVertexAlg.h

+/*
+  Copyright (C) 2002-2021 CERN for the benefit of the ATLAS collaboration
+*/
+
+///////////////////////////////////////////////////////////////////
+// AllBVertexAlg.h, (c) ATLAS Detector software
+// author: Vadim Kostyukhin (vadim.kostyukhin@cern.ch)
+///////////////////////////////////////////////////////////////////
+
+/**
+ * Reconstruct all B vertices in events, including soft ones
+*/
+
+#ifndef VRTSECINCLUSIVE_AllBVertexAlg_H
+#define VRTSECINCLUSIVE_AllBVertexAlg_H
+
+#include <string>
+#include <vector>
+
+#include "AthenaBaseComps/AthAlgorithm.h"
+#include "GaudiKernel/ToolHandle.h"
+#include "VrtSecInclusive/IVrtInclusive.h"
+
+
+namespace Rec {
+
+  class AllBVertexAlg : public AthAlgorithm {
+    public: 
+
+      AllBVertexAlg( const std::string& name, ISvcLocator* pSvcLocator );
+
+      StatusCode initialize() override;
+      StatusCode finalize() override;
+      StatusCode execute() override;
+
+    private:
+
+      std::string m_bvertexContainerName;
+      std::string m_trContainerName;
+      std::string m_pvContainerName;
+
+      ToolHandle < Rec::IVrtInclusive > m_bvertextool;
+      double vrtVrtDist(const xAOD::Vertex & primVrt, const xAOD::Vertex & secVrt) const;
+
+  };
+}
+
+#endif // VRTSECINCLUSIVE_AllBVertexAlg_H

Reconstruction/VKalVrt/VrtSecInclusive/VrtSecInclusive/MultiSecVertexTool.h

@@ -146,6 +146,7 @@ namespace Rec {
       long int m_cutSctHits{};
       long int m_cutPixelHits{};
       long int m_cutSiHits{};
+      long int m_cutTRTHits{};
       long int m_cutBLayHits{};
       long int m_cutSharedHits{};
       double m_cutPt{};

Reconstruction/VKalVrt/VrtSecInclusive/python/MultiSecVertexTool.py

-# Copyright (C) 2002-2020 CERN for the benefit of the ATLAS collaboration
+# Copyright (C) 2002-2021 CERN for the benefit of the ATLAS collaboration
 # Author: Vadim Kostyukhin vadim.kostyukhin@cern.ch
 from GaudiKernel.GaudiHandles import *
 from GaudiKernel.Proxy.Configurable import *
@@ -69,15 +69,16 @@ class InclusiveBFinderTool( Rec__MultiSecVertexTool ):
         Rec__MultiSecVertexTool.__init__( self, name = name,
                                              VertexFitterTool   = SVertexFitterTool,
 					     CutBLayHits  = 0,
-					     CutPixelHits = 2,
+					     CutPixelHits = 1,
 					     CutSiHits    = 8,
+					     CutTRTHits   = 10,
 					     useVertexCleaning  = True,
 					     MultiWithOneTrkVrt = True,
 					     removeTrkMatSignif = -1.,    # No additional material rejection
 					     AntiPileupSigRCut = 2.0,
 					     TrkSigCut         = 2.0,
 					     SelVrtSigCut      = 3.0,
-					     v2tIniBDTCut      =-0.5,
+					     v2tIniBDTCut      =-0.7,
 					     v2tFinBDTCut      = 0.0,
 					     CutPt             = 500
                                              )
@@ -183,7 +184,29 @@ class DVFinderTool( Rec__MultiSecVertexTool ):
 					     SelVrtSigCut      = 8.0,
 					     v2tIniBDTCut      =-0.5,
 					     v2tFinBDTCut      = 0.0,
-					     VrtMassLimit      = 50000.,
-					     Vrt2TrMassLimit   = 50000.,
+					     VrtMassLimit      = 1000000.,
+					     Vrt2TrMassLimit   = 100000.,
 					     CutPt             = 1000.
                                              )
+
+##########################################################################################################
+# define the Test algorithms
+
+from VrtSecInclusiveConf import Rec__AllBVertexAlg
+class AllBVertexFinderAlg( Rec__AllBVertexAlg ):
+
+    def __init__(self, name = 'AllBVertexFinderAlg'  ):
+
+        from AthenaCommon.AppMgr import ToolSvc
+        mlog = logging.getLogger( 'AllBVertexFinderAlg::__init__ ' )
+        mlog.info("entering")
+        
+        BFinderTool = InclusiveBFinderTool()
+        ToolSvc += BFinderTool
+        #----------------------
+        # All B-hadron vertex finder itself
+        #
+        Rec__AllBVertexAlg.__init__( self, name = name,
+                                     BVertexTool = BFinderTool
+                                   )
+ 

Reconstruction/VKalVrt/VrtSecInclusive/src/AllBVertexAlg.cxx

+/*
+  Copyright (C) 2002-2021 CERN for the benefit of the ATLAS collaboration
+*/
+
+///////////////////////////////////////////////////////////////////
+// AllBVertexAlg.cxx, (c) ATLAS Detector software
+// Author: Vadim Kostyukhin (vadim.kostyukhin@cern.ch)
+///////////////////////////////////////////////////////////////////
+
+#include "VrtSecInclusive/AllBVertexAlg.h"
+
+#include "xAODTracking/TrackParticleContainer.h"
+#include "xAODTracking/VertexContainer.h"
+#include "xAODTracking/VertexAuxContainer.h"
+#include "VrtSecInclusive/MultiSecVertexTool.h"
+
+#include "TLorentzVector.h"
+#include "CxxUtils/sincos.h"
+
+namespace Rec {
+
+  static const SG::AuxElement::Decorator<float> bvrtM("bvrtM");  
+  static const SG::AuxElement::Decorator<float> bvrtPt("bvrtPt");  
+  static const SG::AuxElement::Decorator<float> bvrtPhi("bvrtPhi");  
+  static const SG::AuxElement::Decorator<float> bvrtEta("bvrtEta");  
+  static const SG::AuxElement::Decorator<float> bvrtMinTPt("bvrtMinTPt");  
+  static const SG::AuxElement::Decorator<float> bvrtSig3D("bvrtSig3D");  
+
+   AllBVertexAlg::AllBVertexAlg(const std::string& name, ISvcLocator* pSvcLocator) :
+    AthAlgorithm( name, pSvcLocator ),
+    m_bvertexContainerName("AllBVertices"),
+    m_trContainerName("InDetTrackParticles"),
+    m_pvContainerName("PrimaryVertices"),
+    m_bvertextool("Rec::MultiSecVertexTool/InclusivetBVertexTool",this)
+  {
+    declareProperty("BVertexContainerName", m_bvertexContainerName);    
+    declareProperty("TrackContainerName", m_trContainerName);
+    declareProperty("PVContainerName",m_pvContainerName);
+    declareProperty("BVertexTool",m_bvertextool);
+  }
+
+  StatusCode AllBVertexAlg::initialize()
+  {
+    if (m_bvertexContainerName=="") {
+      ATH_MSG_ERROR("No SG name provided for the output of AllBVertexAlg!");
+      return StatusCode::FAILURE;
+    }
+
+    CHECK( m_bvertextool.retrieve() );
+
+    return StatusCode::SUCCESS;
+  }
+
+  StatusCode AllBVertexAlg::finalize()
+  {
+    return StatusCode::SUCCESS;
+  }
+
+  StatusCode AllBVertexAlg::execute()
+  {
+    // check that container we want to write in SG does not yet exist
+    if (evtStore()->contains<xAOD::VertexContainer>(m_bvertexContainerName)) {
+      ATH_MSG_ERROR("Algorithm is attempting to write a StoreGate key " << m_bvertexContainerName 
+                     << " which already exists. Please use a different key");
+      return StatusCode::FAILURE;
+    }
+
+    const xAOD::VertexContainer* pv_cont = nullptr;
+    ATH_CHECK( evtStore()->retrieve( pv_cont, m_pvContainerName ) );
+
+    std::unique_ptr<xAOD::VertexContainer> bVertexContainer = std::make_unique<xAOD::VertexContainer>();
+    std::unique_ptr<xAOD::VertexAuxContainer> bVertexAuxContainer = std::make_unique<xAOD::VertexAuxContainer>();
+    bVertexContainer->setStore(bVertexAuxContainer.get());
+
+    const xAOD::Vertex* pv = nullptr;
+    for ( auto v : *pv_cont ) {
+      if (v->vertexType()==xAOD::VxType::PriVtx) { pv = v;  break; }
+    }
+    if (!pv) {
+      ATH_MSG_WARNING("Primary vertex not found");
+      CHECK(evtStore()->record( bVertexContainer.release(), m_bvertexContainerName));
+      CHECK(evtStore()->record( bVertexAuxContainer.release(), m_bvertexContainerName+"Aux."));
+      return StatusCode::SUCCESS;
+    }
+    ATH_MSG_DEBUG("Found PV");
+
+    // retrieve particle collections
+    const xAOD::TrackParticleContainer* tp_cont = nullptr;
+    ATH_CHECK( evtStore()->retrieve( tp_cont, m_trContainerName ) );
+    std::vector<const xAOD::TrackParticle*> trkparticles(0);
+    for(auto tp : (*tp_cont)) trkparticles.push_back(tp);
+
+    std::unique_ptr<Trk::VxSecVertexInfo> foundVrts = m_bvertextool->findAllVertices(trkparticles,*pv);
+    if(foundVrts && foundVrts->vertices().size()){
+       const std::vector<xAOD::Vertex*> vtmp=foundVrts->vertices();
+       for(auto & iv :  vtmp) {
+          bVertexContainer->push_back(iv);
+          std::vector< Trk::VxTrackAtVertex > & vtrk = iv->vxTrackAtVertex();
+          TLorentzVector VSUM(0.,0.,0.,0.);
+          TLorentzVector tmp;
+          double minTPt=1.e10;
+          for(auto & it : vtrk){
+             const Trk::Perigee* mPer = dynamic_cast<const Trk::Perigee*>(it.perigeeAtVertex());
+             CxxUtils::sincos   phi(mPer->parameters()[Trk::phi]);
+             CxxUtils::sincos theta(mPer->parameters()[Trk::theta]);
+             double absP  =  1./std::abs(mPer->parameters()[Trk::qOverP]);
+             tmp.SetXYZM( phi.cs*theta.sn*absP, phi.sn*theta.sn*absP, theta.cs*absP, Trk::ParticleMasses().mass[Trk::pion]);
+             minTPt=std::min(minTPt,tmp.Pt());
+             VSUM+=tmp;
+          }
+          bvrtM(*iv)  =VSUM.M();
+          bvrtPt(*iv) =VSUM.Pt();
+          bvrtEta(*iv)=VSUM.Eta();
+          bvrtPhi(*iv)=VSUM.Phi();
+          bvrtMinTPt(*iv)=minTPt;
+          bvrtSig3D(*iv)=vrtVrtDist(*pv,*iv);
+       }
+    }
+
+    ATH_MSG_DEBUG("BVertex container size: " << bVertexContainer->size());
+    CHECK(evtStore()->record( bVertexContainer.release(), m_bvertexContainerName));
+    CHECK(evtStore()->record( bVertexAuxContainer.release(), m_bvertexContainerName+"Aux."));
+
+    return StatusCode::SUCCESS;
+  }
+
+  double AllBVertexAlg::vrtVrtDist(const xAOD::Vertex & primVrt, const xAOD::Vertex & secVrt)
+  const
+  {
+     double distx =  secVrt.x()- primVrt.x();
+     double disty =  secVrt.y()- primVrt.y();
+     double distz =  secVrt.z()- primVrt.z();
+
+     AmgSymMatrix(3)  covMtx=primVrt.covariancePosition()+secVrt.covariancePosition();
+ 
+     AmgSymMatrix(3)  wgtMtx = covMtx.inverse();
+     if(wgtMtx(0,0)<=0. || wgtMtx(1,1)<=0. || wgtMtx(2,2)<=0.) return 1.e10; //rudimentary corrrectness check
+
+     double signif = distx*wgtMtx(0,0)*distx
+                    +disty*wgtMtx(1,1)*disty
+                    +distz*wgtMtx(2,2)*distz
+                 +2.*distx*wgtMtx(0,1)*disty
+                 +2.*distx*wgtMtx(0,2)*distz
+                 +2.*disty*wgtMtx(1,2)*distz;
+     signif=std::sqrt(std::abs(signif));
+     if( signif!=signif ) signif = 0.;
+     return signif;
+  }
+
+}

Reconstruction/VKalVrt/VrtSecInclusive/src/MultiSecVertexTool.cxx

@@ -42,6 +42,7 @@ MultiSecVertexTool::MultiSecVertexTool(const std::string& type,
     m_cutSctHits(4),
     m_cutPixelHits(2),
     m_cutSiHits(8),
+    m_cutTRTHits(10),
     m_cutBLayHits(0),
     m_cutSharedHits(1),
     m_cutPt(500.),
@@ -79,6 +80,7 @@ MultiSecVertexTool::MultiSecVertexTool(const std::string& type,
     declareProperty("CutSctHits",    m_cutSctHits ,  "Remove track is it has less SCT hits");
     declareProperty("CutPixelHits",  m_cutPixelHits, "Remove track is it has less Pixel hits");
     declareProperty("CutSiHits",     m_cutSiHits,    "Remove track is it has less Pixel+SCT hits");
+    declareProperty("CutTRTHits",    m_cutTRTHits,   "Remove track is it has less TRT hits");
     declareProperty("CutBLayHits",   m_cutBLayHits,  "Remove track is it has less B-layer hits");
     declareProperty("CutSharedHits", m_cutSharedHits,"Reject final 2tr vertices if tracks have shared hits");
 
@@ -1074,10 +1076,11 @@ MultiSecVertexTool::MultiSecVertexTool(const std::string& type,
           const Trk::Perigee mPer=(*i_ntrk)->perigeeParameters() ;
           double CovTrkMtx00 = (*(mPer.covariance()))(0,0);
  
-          uint8_t PixelHits,SctHits,BLayHits;
+          uint8_t PixelHits,SctHits,BLayHits,TRTHits;
           if( !((*i_ntrk)->summaryValue(PixelHits,xAOD::numberOfPixelHits)) )   continue; // Track is 
           if( !((*i_ntrk)->summaryValue(  SctHits,xAOD::numberOfSCTHits))   )   continue; // definitely  
           if( SctHits<3 )                                                       continue; // bad
+          if( !((*i_ntrk)->summaryValue(  TRTHits,xAOD::numberOfTRTHits))   )              TRTHits=0;
           if( !((*i_ntrk)->summaryValue(BLayHits,xAOD::numberOfInnermostPixelLayerHits)))  BLayHits=0;
           if(m_fillHist){ m_hb_trkSelect->Fill( 2., m_w_1);}
 
@@ -1109,6 +1112,7 @@ MultiSecVertexTool::MultiSecVertexTool(const std::string& type,
           if(SctHits	    < m_cutSctHits) 		continue;
           if((PixelHits+SctHits) < m_cutSiHits) 	continue;
           if(BLayHits	    < m_cutBLayHits) 		continue;
+          if(std::abs((*i_ntrk)->eta())<1.9 && TRTHits < m_cutTRTHits) continue; //TRT hits must be present inside TRT
           if(m_fillHist){ m_hb_trkSelect->Fill( 5., m_w_1);}
           orderedTrk.emplace(signifBeam,*i_ntrk);
           selectedTracks.push_back(*i_ntrk);

Reconstruction/VKalVrt/VrtSecInclusive/src/components/VrtSecInclusive_entries.cxx

@@ -3,17 +3,20 @@
 #include "VrtSecInclusive/IPAugmentor.h"
 #include "VrtSecInclusive/TrackRandomizer.h"
 #include "VrtSecInclusive/SoftBtagTrackSelector.h"
+#include "VrtSecInclusive/AllBVertexAlg.h"
 #include "GaudiKernel/DeclareFactoryEntries.h"
 DECLARE_NAMESPACE_ALGORITHM_FACTORY( VKalVrtAthena, VrtSecInclusive )
 DECLARE_NAMESPACE_ALGORITHM_FACTORY( VKalVrtAthena, IPAugmentor )
 DECLARE_NAMESPACE_ALGORITHM_FACTORY( VKalVrtAthena, TrackRandomizer )
 DECLARE_NAMESPACE_ALGORITHM_FACTORY( VKalVrtAthena, SoftBtagTrackSelector )
+DECLARE_NAMESPACE_ALGORITHM_FACTORY( Rec, AllBVertexAlg )
 DECLARE_NAMESPACE_TOOL_FACTORY( Rec, MultiSecVertexTool )
 DECLARE_FACTORY_ENTRIES(VrtSecInclusive) {
     DECLARE_NAMESPACE_ALGORITHM( VKalVrtAthena, VrtSecInclusive )
     DECLARE_NAMESPACE_ALGORITHM( VKalVrtAthena, IPAugmentor )
     DECLARE_NAMESPACE_ALGORITHM( VKalVrtAthena, TrackRandomizer )
     DECLARE_NAMESPACE_ALGORITHM( VKalVrtAthena, SoftBtagTrackSelector )
+    DECLARE_NAMESPACE_ALGORITHM( Rec, AllBVertexAlg )
     DECLARE_NAMESPACE_ALGTOOL( Rec, MultiSecVertexTool )
 }
 //Notes: