diff --git a/Control/AthenaBaseComps/AthenaBaseComps/AthReentrantAlgorithm.h b/Control/AthenaBaseComps/AthenaBaseComps/AthReentrantAlgorithm.h
index 21fb2f462383f848c26d6f5be329fbb7158dcedd..4784ffabce82c56c02a4eb27a532068688d3a828 100644
--- a/Control/AthenaBaseComps/AthenaBaseComps/AthReentrantAlgorithm.h
+++ b/Control/AthenaBaseComps/AthenaBaseComps/AthReentrantAlgorithm.h
@@ -136,6 +136,14 @@ class AthReentrantAlgorithm
*/
virtual const DataObjIDColl& extraOutputDeps() const override;
+ virtual bool filterPassed(const EventContext& ctx) const {
+ return execState( ctx ).filterPassed();
+ }
+
+ virtual void setFilterPassed( bool state, const EventContext& ctx ) const {
+ execState( ctx ).setFilterPassed( state );
+ }
+
private:
diff --git a/MuonSpectrometer/MuonCalib/NSWCalib/NSWCalibTools/src/NSWCalibTool.cxx b/MuonSpectrometer/MuonCalib/NSWCalib/NSWCalibTools/src/NSWCalibTool.cxx
index e062ad50ba94e0c86e717aa0fca9488b1e07ffd2..06069815b7a3fec0e39574c4c5d4ec516cae11bd 100644
--- a/MuonSpectrometer/MuonCalib/NSWCalib/NSWCalibTools/src/NSWCalibTool.cxx
+++ b/MuonSpectrometer/MuonCalib/NSWCalib/NSWCalibTools/src/NSWCalibTool.cxx
@@ -140,12 +140,9 @@ StatusCode Muon::NSWCalibTool::calibrateStrip(const Muon::MM_RawData* mmRawData,
{
Identifier rdoId = mmRawData->identify();
- // @TODO: for the shifter: this code will be used as soon as the NSW is built in the MuonDetectorCondAlg
- // for now, we will need to use the nominal manager from the detectorStore to run
- // SG::ReadCondHandle<MuonGM::MuonDetectorManager> muDetMgrHandle{m_muDetMgrKey};
- // const MuonGM::MuonDetectorManager* muDetMgr = muDetMgrHandle.cptr();
- const MuonGM::MuonDetectorManager* muDetMgr=nullptr;
- ATH_CHECK(detStore()->retrieve(muDetMgr));
+ // MuonDetectorManager from the conditions store
+ SG::ReadCondHandle<MuonGM::MuonDetectorManager> muDetMgrHandle{m_muDetMgrKey};
+ const MuonGM::MuonDetectorManager* muDetMgr = muDetMgrHandle.cptr();
//get globalPos
Amg::Vector3D globalPos;
diff --git a/MuonSpectrometer/MuonCnv/MuonSTGC_CnvTools/src/sTgcRdoToPrepDataToolCore.cxx b/MuonSpectrometer/MuonCnv/MuonSTGC_CnvTools/src/sTgcRdoToPrepDataToolCore.cxx
index 25924a048ded723f70e793f9a108882ddb3d4752..5166bfd20f33d8f0fd0720ebdfc0225a0d6b0e04 100644
--- a/MuonSpectrometer/MuonCnv/MuonSTGC_CnvTools/src/sTgcRdoToPrepDataToolCore.cxx
+++ b/MuonSpectrometer/MuonCnv/MuonSTGC_CnvTools/src/sTgcRdoToPrepDataToolCore.cxx
@@ -83,6 +83,15 @@ StatusCode Muon::sTgcRdoToPrepDataToolCore::processCollection(const STGC_RawData
std::vector<sTgcPrepData> sTgcPadPrds;
// convert the RDO collection to a PRD collection
STGC_RawDataCollection::const_iterator it = rdoColl->begin();
+
+ // MuonDetectorManager from the conditions store
+ SG::ReadCondHandle<MuonGM::MuonDetectorManager> detMgrHandle{m_muDetMgrKey};
+ const MuonGM::MuonDetectorManager* muonDetMgr = detMgrHandle.cptr();
+ if(!muonDetMgr){
+ ATH_MSG_ERROR("Null pointer to the read MuonDetectorManager conditions object");
+ return StatusCode::FAILURE;
+ }
+
for ( ; it != rdoColl->end() ; ++it ) {
ATH_MSG_DEBUG("Adding a new sTgc PrepRawData");
@@ -95,14 +104,9 @@ StatusCode Muon::sTgcRdoToPrepDataToolCore::processCollection(const STGC_RawData
}
std::vector<Identifier> rdoList;
rdoList.push_back(rdoId);
-
- // TODO: this needs to be replaced by SG::ReadCondHandle<MuonGM::MuonDetectorManager>
- // will do it in a follow-up MR, since for now, we need to get the Run2 detectors running, so skip sTGCs for now
- const MuonGM::MuonDetectorManager* muDetMgrNominal=nullptr;
- ATH_CHECK(detStore()->retrieve(muDetMgrNominal));
// get the local and global positions
- const MuonGM::sTgcReadoutElement* detEl = muDetMgrNominal->getsTgcReadoutElement(rdoId);
+ const MuonGM::sTgcReadoutElement* detEl = muonDetMgr->getsTgcReadoutElement(rdoId);
Amg::Vector2D localPos;
bool getLocalPos = detEl->stripPosition(rdoId,localPos);
diff --git a/MuonSpectrometer/MuonIdHelpersAlgs/MuonIdHelpersAlgs/TestMuonIdHelpers.h b/MuonSpectrometer/MuonIdHelpersAlgs/MuonIdHelpersAlgs/TestMuonIdHelpers.h
index f8e34a66ce252f7a5008dd2bd6b8f3e6c241277c..704ebace7eb1cfa2deafc13abd0b645401e6332f 100644
--- a/MuonSpectrometer/MuonIdHelpersAlgs/MuonIdHelpersAlgs/TestMuonIdHelpers.h
+++ b/MuonSpectrometer/MuonIdHelpersAlgs/MuonIdHelpersAlgs/TestMuonIdHelpers.h
@@ -5,16 +5,18 @@
#ifndef TestMuonIdHelpers_H
#define TestMuonIdHelpers_H
-#include <vector>
-
#include "AthenaBaseComps/AthAlgorithm.h"
#include "GaudiKernel/ServiceHandle.h"
-#include "MuonIdHelpers/IMuonIdHelperSvc.h"
+#include "StoreGate/ReadHandleKey.h"
-class StoreGateSvc;
-class ActiveStoreSvc;
+#include "MuonIdHelpers/IMuonIdHelperSvc.h"
+#include "MuonDigitContainer/MdtDigitContainer.h"
+#include "MuonDigitContainer/RpcDigitContainer.h"
+#include "MuonDigitContainer/CscDigitContainer.h"
+#include "MuonDigitContainer/TgcDigitContainer.h"
-/////////////////////////////////////////////////////////////////////////////
+#include <string>
+#include <vector>
class TestMuonIdHelpers : public AthAlgorithm {
@@ -44,9 +46,13 @@ private:
BooleanProperty m_testRPC;
BooleanProperty m_testTGC;
- ServiceHandle<ActiveStoreSvc> m_activeStore;
ServiceHandle<Muon::IMuonIdHelperSvc> m_idHelperSvc {this, "MuonIdHelperSvc", "Muon::MuonIdHelperSvc/MuonIdHelperSvc"};
+ SG::ReadHandleKey<MdtDigitContainer> m_mdtDigitContKey{this,"MdtDigits","MDT_DIGITS","ReadHandleKey for Input MdtDigitContainer"};
+ SG::ReadHandleKey<RpcDigitContainer> m_rpcDigitContKey{this,"RpcDigits","RPC_DIGITS","ReadHandleKey for Input RpcDigitContainer"};
+ SG::ReadHandleKey<CscDigitContainer> m_cscDigitContKey{this,"CscDigits","CSC_DIGITS","ReadHandleKey for Input CscDigitContainer"};
+ SG::ReadHandleKey<TgcDigitContainer> m_tgcDigitContKey{this,"TgcDigits","TGC_DIGITS","ReadHandleKey for Input TgcDigitContainer"};
+
mutable std::atomic<long long> m_deltaUser{0};
mutable std::atomic<long long> m_deltaKernel{0};
mutable std::atomic<long long> m_deltaElapsed{0};
diff --git a/MuonSpectrometer/MuonIdHelpersAlgs/src/TestMuonIdHelpers.cxx b/MuonSpectrometer/MuonIdHelpersAlgs/src/TestMuonIdHelpers.cxx
index 6cb950c02a130777da28783c103f11a567127d3c..d937dbe05b38ef13f2092ef76a89586ba35ebb4e 100644
--- a/MuonSpectrometer/MuonIdHelpersAlgs/src/TestMuonIdHelpers.cxx
+++ b/MuonSpectrometer/MuonIdHelpersAlgs/src/TestMuonIdHelpers.cxx
@@ -4,33 +4,11 @@
#include "MuonIdHelpersAlgs/TestMuonIdHelpers.h"
-#include "MuonDigitContainer/MdtDigitContainer.h"
-#include "MuonDigitContainer/MdtDigitCollection.h"
-#include "MuonDigitContainer/MdtDigit.h"
-
-#include "MuonDigitContainer/CscDigitContainer.h"
-#include "MuonDigitContainer/CscDigitCollection.h"
-#include "MuonDigitContainer/CscDigit.h"
-
-#include "MuonDigitContainer/RpcDigitContainer.h"
-#include "MuonDigitContainer/RpcDigitCollection.h"
-#include "MuonDigitContainer/RpcDigit.h"
-
-#include "MuonDigitContainer/TgcDigitContainer.h"
-#include "MuonDigitContainer/TgcDigitCollection.h"
-#include "MuonDigitContainer/TgcDigit.h"
-
-#include "MuonReadoutGeometry/MuonDetectorManager.h"
-
#include <algorithm>
#include <cmath>
-/////////////////////////////////////////////////////////////////////////////
-/////////////////////////////////////////////////////////////////////////////
-
TestMuonIdHelpers::TestMuonIdHelpers(const std::string& name, ISvcLocator* pSvcLocator) :
- AthAlgorithm(name, pSvcLocator),
- m_activeStore("ActiveStoreSvc", name) {
+ AthAlgorithm(name, pSvcLocator) {
m_testMDT = true;
m_testCSC = true;
@@ -51,10 +29,13 @@ StatusCode TestMuonIdHelpers::initialize(){
ATH_MSG_DEBUG(" in initialize()");
- ATH_CHECK(m_activeStore.retrieve());
-
ATH_CHECK(m_idHelperSvc.retrieve());
+ ATH_CHECK(m_mdtDigitContKey.initialize());
+ ATH_CHECK(m_rpcDigitContKey.initialize());
+ ATH_CHECK(m_cscDigitContKey.initialize());
+ ATH_CHECK(m_tgcDigitContKey.initialize());
+
return StatusCode::SUCCESS;
}
@@ -104,10 +85,12 @@ StatusCode TestMuonIdHelpers::testMdtIdHelper() const {
// retrieve data from StoreGate
typedef MdtDigitContainer::const_iterator collection_iterator;
typedef MdtDigitCollection::const_iterator digit_iterator;
-
- std::string key = "MDT_DIGITS";
- const DataHandle <MdtDigitContainer> container;
- ATH_CHECK( (*m_activeStore)->retrieve(container,key) );
+
+ SG::ReadHandle<MdtDigitContainer> container(m_mdtDigitContKey);
+ if (!container.isValid()) {
+ ATH_MSG_ERROR("Could not find MdtDigitContainer called " << container.name() << " in store " << container.store());
+ return StatusCode::FAILURE;
+ }
collection_iterator it1_coll= container->begin();
collection_iterator it2_coll= container->end();
@@ -121,7 +104,8 @@ StatusCode TestMuonIdHelpers::testMdtIdHelper() const {
++nc1;
Identifier moduleId = mdtCollection->identify();
if (!m_idHelperSvc->mdtIdHelper().validElement(moduleId)) {
- ATH_MSG_ERROR( "Invalid MDT module identifier " << m_idHelperSvc->mdtIdHelper().show_to_string(moduleId) );
+ ATH_MSG_ERROR( "Invalid MDT module identifier " << m_idHelperSvc->mdtIdHelper().show_to_string(moduleId));
+ return StatusCode::FAILURE;
}
digit_iterator it1_digit = mdtCollection->begin();
digit_iterator it2_digit = mdtCollection->end();
@@ -138,9 +122,11 @@ StatusCode TestMuonIdHelpers::testMdtIdHelper() const {
}
ATH_MSG_DEBUG(" Number of MDT Collections Accessed " << nc1 );
- if(error) ATH_MSG_ERROR( "Check of MDT ids FAILED " );
- else ATH_MSG_INFO( "Check of MDT ids OK " );
-
+ if(error) {
+ ATH_MSG_ERROR( "Check of MDT ids FAILED " );
+ return StatusCode::FAILURE;
+ }
+ else ATH_MSG_INFO( "Check of MDT ids OK " );
return StatusCode::SUCCESS;
}
@@ -158,10 +144,12 @@ StatusCode TestMuonIdHelpers::testCscIdHelper() const {
// retrieve data from StoreGate
typedef CscDigitContainer::const_iterator collection_iterator;
typedef CscDigitCollection::const_iterator digit_iterator;
-
- std::string key = "CSC_DIGITS";
- const DataHandle <CscDigitContainer> container;
- ATH_CHECK( (*m_activeStore)->retrieve(container,key) );
+
+ SG::ReadHandle<CscDigitContainer> container(m_cscDigitContKey);
+ if (!container.isValid()) {
+ ATH_MSG_ERROR("Could not find CscDigitContainer called " << container.name() << " in store " << container.store());
+ return StatusCode::FAILURE;
+ }
collection_iterator it1_coll= container->begin();
collection_iterator it2_coll= container->end();
@@ -175,7 +163,8 @@ StatusCode TestMuonIdHelpers::testCscIdHelper() const {
++nc1;
Identifier moduleId = cscCollection->identify();
if (!m_idHelperSvc->cscIdHelper().validElement(moduleId)) {
- ATH_MSG_ERROR( "Invalid CSC module identifier " << m_idHelperSvc->cscIdHelper().show_to_string(moduleId) );
+ ATH_MSG_ERROR( "Invalid CSC module identifier " << m_idHelperSvc->cscIdHelper().show_to_string(moduleId));
+ return StatusCode::FAILURE;
}
digit_iterator it1_digit = cscCollection->begin();
digit_iterator it2_digit = cscCollection->end();
@@ -192,9 +181,10 @@ StatusCode TestMuonIdHelpers::testCscIdHelper() const {
}
ATH_MSG_DEBUG("Number of CSC Collections Accessed " << nc1 );
- if(error) ATH_MSG_ERROR( "Check of CSC ids FAILED " );
- else ATH_MSG_INFO( "Check of CSC ids OK " );
-
+ if(error) {
+ ATH_MSG_ERROR( "Check of CSC ids FAILED " );
+ return StatusCode::FAILURE;
+ } else ATH_MSG_INFO( "Check of CSC ids OK " );
return StatusCode::SUCCESS;
}
@@ -211,10 +201,12 @@ StatusCode TestMuonIdHelpers::testRpcIdHelper() const {
// retrieve data from StoreGate
typedef RpcDigitContainer::const_iterator collection_iterator;
typedef RpcDigitCollection::const_iterator digit_iterator;
-
- std::string key = "RPC_DIGITS";
- const DataHandle <RpcDigitContainer> container;
- ATH_CHECK( (*m_activeStore)->retrieve(container,key) );
+
+ SG::ReadHandle<RpcDigitContainer> container(m_rpcDigitContKey);
+ if (!container.isValid()) {
+ ATH_MSG_ERROR("Could not find RpcDigitContainer called " << container.name() << " in store " << container.store());
+ return StatusCode::FAILURE;
+ }
collection_iterator it1_coll= container->begin();
collection_iterator it2_coll= container->end();
@@ -228,7 +220,8 @@ StatusCode TestMuonIdHelpers::testRpcIdHelper() const {
++nc1;
Identifier moduleId = rpcCollection->identify();
if (!m_idHelperSvc->rpcIdHelper().validElement(moduleId)) {
- ATH_MSG_ERROR( "Invalid RPC module identifier " << m_idHelperSvc->rpcIdHelper().show_to_string(moduleId) );
+ ATH_MSG_ERROR( "Invalid RPC module identifier " << m_idHelperSvc->rpcIdHelper().show_to_string(moduleId));
+ return StatusCode::FAILURE;
}
digit_iterator it1_digit = rpcCollection->begin();
digit_iterator it2_digit = rpcCollection->end();
@@ -245,8 +238,11 @@ StatusCode TestMuonIdHelpers::testRpcIdHelper() const {
}
ATH_MSG_DEBUG(" Number of RPC Collections Accessed " << nc1 );
- if(error) ATH_MSG_ERROR( "Check of RPC ids FAILED " );
- else ATH_MSG_INFO( "Check of RPC ids OK " );
+ if(error) {
+ ATH_MSG_ERROR("Check of RPC ids FAILED");
+ return StatusCode::FAILURE;
+ }
+ else ATH_MSG_INFO( "Check of RPC ids OK " );
return StatusCode::SUCCESS;
}
@@ -264,10 +260,12 @@ StatusCode TestMuonIdHelpers::testTgcIdHelper() const {
// retrieve data from StoreGate
typedef TgcDigitContainer::const_iterator collection_iterator;
typedef TgcDigitCollection::const_iterator digit_iterator;
-
- std::string key = "TGC_DIGITS";
- const DataHandle <TgcDigitContainer> container;
- ATH_CHECK( (*m_activeStore)->retrieve(container,key) );
+
+ SG::ReadHandle<TgcDigitContainer> container(m_tgcDigitContKey);
+ if (!container.isValid()) {
+ ATH_MSG_ERROR("Could not find TgcDigitContainer called " << container.name() << " in store " << container.store());
+ return StatusCode::FAILURE;
+ }
collection_iterator it1_coll= container->begin();
collection_iterator it2_coll= container->end();
@@ -281,7 +279,8 @@ StatusCode TestMuonIdHelpers::testTgcIdHelper() const {
++nc1;
Identifier moduleId = tgcCollection->identify();
if (!m_idHelperSvc->tgcIdHelper().validElement(moduleId)) {
- ATH_MSG_ERROR( "Invalid TGC module identifier " << m_idHelperSvc->tgcIdHelper().show_to_string(moduleId) );
+ ATH_MSG_ERROR( "Invalid TGC module identifier " << m_idHelperSvc->tgcIdHelper().show_to_string(moduleId));
+ return StatusCode::FAILURE;
}
digit_iterator it1_digit = tgcCollection->begin();
digit_iterator it2_digit = tgcCollection->end();
@@ -298,8 +297,10 @@ StatusCode TestMuonIdHelpers::testTgcIdHelper() const {
}
ATH_MSG_DEBUG("Number of TGC Collections Accessed " << nc1 );
- if(error) ATH_MSG_ERROR( "Check of TGC ids FAILED " );
- else ATH_MSG_INFO( "Check of TGC ids OK " );
+ if(error) {
+ ATH_MSG_ERROR( "Check of TGC ids FAILED ");
+ return StatusCode::FAILURE;
+ } else ATH_MSG_INFO( "Check of TGC ids OK ");
return StatusCode::SUCCESS;
}
diff --git a/PhysicsAnalysis/CommonTools/ExpressionEvaluation/ExpressionEvaluation/StackElement.icc b/PhysicsAnalysis/CommonTools/ExpressionEvaluation/ExpressionEvaluation/StackElement.icc
index 309b6b000b72c8373cfb6300906fe45026cda77a..9f728b7e38a8f67361a789e43829f2a88216c2ae 100644
--- a/PhysicsAnalysis/CommonTools/ExpressionEvaluation/ExpressionEvaluation/StackElement.icc
+++ b/PhysicsAnalysis/CommonTools/ExpressionEvaluation/ExpressionEvaluation/StackElement.icc
@@ -37,12 +37,12 @@ namespace ExpressionParsing {
m_variableType(static_cast<IProxyLoader::VariableType>(a.m_variableType)),
m_determinedVariableType(a.m_determinedVariableType ? true : false)
{
- if (a.m_moved) throw std::logic_error("Content already moved"); \
+ if (a.m_moved) throw std::logic_error("Content already moved");
}
inline StackElement &StackElement::operator=(StackElement &&a) {
if (&a != this) {
- if (a.m_moved) throw std::logic_error("Content already moved"); \
+ if (a.m_moved) throw std::logic_error("Content already moved");
m_type = a.m_type;
m_intVal = a.m_intVal;
m_doubleVal = a.m_doubleVal;
diff --git a/Reconstruction/VKalVrt/NewVrtSecInclusiveTool/CMakeLists.txt b/Reconstruction/VKalVrt/NewVrtSecInclusiveTool/CMakeLists.txt
new file mode 100644
index 0000000000000000000000000000000000000000..74c6928b2aa9f43cfaa6dc09848413b590ef39e8
--- /dev/null
+++ b/Reconstruction/VKalVrt/NewVrtSecInclusiveTool/CMakeLists.txt
@@ -0,0 +1,40 @@
+################################################################################
+# Package: NewVrtSecInclusiveTool
+################################################################################
+
+# Declare the package name:
+atlas_subdir( NewVrtSecInclusiveTool )
+
+# Declare the package's dependencies:
+atlas_depends_on_subdirs( PUBLIC
+ Control/AthenaBaseComps
+ Event/xAOD/xAODTracking
+ GaudiKernel
+ PRIVATE
+ Tools/PathResolver
+ DetectorDescription/GeoPrimitives
+ Tracking/TrkExtrapolation/TrkExInterfaces
+ Reconstruction/MVAUtils
+ Tracking/TrkEvent/VxSecVertex
+ Tracking/TrkVertexFitter/TrkVKalVrtFitter
+ InnerDetector/InDetConditions/InDetBeamSpotService)
+
+# External dependencies:
+find_package( Boost COMPONENTS filesystem thread system )
+find_package( Eigen )
+find_package( ROOT COMPONENTS Core Tree MathCore Hist RIO pthread TMVA )
+
+# Component(s) in the package:
+atlas_add_component( NewVrtSecInclusiveTool
+ src/*.cxx src/components/*.cxx
+ INCLUDE_DIRS ${ROOT_INCLUDE_DIRS} ${Boost_INCLUDE_DIRS} ${EIGEN_INCLUDE_DIRS}
+ LINK_LIBRARIES ${ROOT_LIBRARIES} ${Boost_LIBRARIES} ${EIGEN_LIBRARIES} AthenaBaseComps xAODTracking PathResolver GaudiKernel
+ Particle TrkVKalVrtFitterLib GeoPrimitives MVAUtils AnalysisUtilsLib TrkParticleBase
+ TrkTrackSummary TrkExInterfaces VxSecVertex )
+
+# Install files from the package:
+atlas_install_headers( NewVrtSecInclusiveTool )
+atlas_install_python_modules( python/*.py )
+#atlas_install_joboptions( share/*.txt )
+#atlas_install_data( share/*.root )
+
diff --git a/Reconstruction/VKalVrt/NewVrtSecInclusiveTool/NewVrtSecInclusiveTool/ATLAS_CHECK_THREAD_SAFETY b/Reconstruction/VKalVrt/NewVrtSecInclusiveTool/NewVrtSecInclusiveTool/ATLAS_CHECK_THREAD_SAFETY
new file mode 100644
index 0000000000000000000000000000000000000000..0c630dd36b6cb9c47f9ca3a056eb224244c94ba0
--- /dev/null
+++ b/Reconstruction/VKalVrt/NewVrtSecInclusiveTool/NewVrtSecInclusiveTool/ATLAS_CHECK_THREAD_SAFETY
@@ -0,0 +1 @@
+Reconstruction/VKalVrt/NewVrtSecInclusiveTool
diff --git a/Reconstruction/VKalVrt/NewVrtSecInclusiveTool/NewVrtSecInclusiveTool/IVrtInclusive.h b/Reconstruction/VKalVrt/NewVrtSecInclusiveTool/NewVrtSecInclusiveTool/IVrtInclusive.h
new file mode 100755
index 0000000000000000000000000000000000000000..3ab765839cbe0654b055558d41923b2854aa5297
--- /dev/null
+++ b/Reconstruction/VKalVrt/NewVrtSecInclusiveTool/NewVrtSecInclusiveTool/IVrtInclusive.h
@@ -0,0 +1,62 @@
+/*
+ Copyright (C) 2002-2020 CERN for the benefit of the ATLAS collaboration
+*/
+//
+// IVrtInclusive.h - Description
+//
+/*
+ Interface for inclusive secondary vertex reconstruction
+ It returns a pointer to Trk::VxSecVertexInfo object which contains
+ vector of pointers to xAOD::Vertex's of found secondary verteces.
+ In case of failure pointer to Trk::VxSecVertexInfo is 0.
+
+
+ Tool creates a derivative object VxSecVKalVertexInfo which contains also additional variables
+ see Tracking/TrkEvent/VxSecVertex/VxSecVertex/VxSecVKalVertexInfo.h
+
+ Author: Vadim Kostyukhin
+ e-mail: vadim.kostyukhin@cern.ch
+
+-----------------------------------------------------------------------------*/
+
+
+
+#ifndef _Rec_IVrtSecInclusive_H
+#define _Rec_IVrtSecInclusive_H
+// Normal STL and physical vectors
+#include <vector>
+// Gaudi includes
+#include "AthenaBaseComps/AthAlgTool.h"
+#include "xAODTracking/TrackParticleContainer.h"
+#include "xAODTracking/VertexContainer.h"
+#include "VxSecVertex/VxSecVertexInfo.h"
+
+
+//------------------------------------------------------------------------
+namespace Rec {
+
+//------------------------------------------------------------------------
+ static const InterfaceID IID_IVrtInclusive("IVrtInclusive", 1, 0);
+
+ class IVrtInclusive : virtual public IAlgTool {
+ public:
+ static const InterfaceID& interfaceID() { return IID_IVrtInclusive;}
+//---------------------------------------------------------------------------
+
+ /** @class IVrtInclusive
+
+ Interface class for inclusive vertex reconstruction.
+ All secondary vertices produced in a phase space around Primary Vertex(PV) and resulted in secondary tracks
+ provided by inputTracks collection are returned in Trk::VxSecVertexInfo.
+ Interface assumes that necessary TrackParticles can be pre-selected by user before calling the actual tool implementation
+ with default track quality cuts.
+
+ */
+ virtual std::unique_ptr<Trk::VxSecVertexInfo> findAllVertices( const std::vector<const xAOD::TrackParticle*> & inputTracks,
+ const xAOD::Vertex & PV) const =0;
+
+ };
+
+} //end namespace
+
+#endif
diff --git a/Reconstruction/VKalVrt/NewVrtSecInclusiveTool/NewVrtSecInclusiveTool/NewVrtSecInclusiveTool.h b/Reconstruction/VKalVrt/NewVrtSecInclusiveTool/NewVrtSecInclusiveTool/NewVrtSecInclusiveTool.h
new file mode 100755
index 0000000000000000000000000000000000000000..2749aa5c1aea70ac5329211bbd2e6b5917bb9c22
--- /dev/null
+++ b/Reconstruction/VKalVrt/NewVrtSecInclusiveTool/NewVrtSecInclusiveTool/NewVrtSecInclusiveTool.h
@@ -0,0 +1,396 @@
+/*
+ Copyright (C) 2002-2020 CERN for the benefit of the ATLAS collaboration
+*/
+
+//
+// NewVrtSecInclusiveTool.h - Description
+//
+/*
+ Tool for inclusive secondary vertex reconstruction
+ It returns a pointer to Trk::VxSecVertexInfo object which contains
+ vector of pointers to xAOD::Vertex's of found secondary verteces.
+ In case of failure pointer to Trk::VxSecVertexInfo is 0.
+
+
+ Tool creates a derivative object VxSecVKalVertexInfo which contains also additional variables
+ see Tracking/TrkEvent/VxSecVertex/VxSecVertex/VxSecVKalVertexInfo.h
+
+
+ Author: Vadim Kostyukhin
+ e-mail: vadim.kostyukhin@cern.ch
+
+-----------------------------------------------------------------------------*/
+
+
+
+#ifndef _VKalVrt_NewVrtSecInclusiveTool_H
+#define _VKalVrt_NewVrtSecInclusiveTool_H
+// Normal STL and physical vectors
+#include <vector>
+// Gaudi includes
+#include "AthenaBaseComps/AthAlgTool.h"
+#include "GaudiKernel/ToolHandle.h"
+#include "GaudiKernel/ServiceHandle.h"
+#include "boost/graph/adjacency_list.hpp"
+//
+#include "xAODTruth/TruthEventContainer.h"
+#include "TrkExInterfaces/IExtrapolator.h"
+//
+#include "VxSecVertex/VxSecVertexInfo.h"
+#include "NewVrtSecInclusiveTool/IVrtInclusive.h"
+
+class TH1D;
+class TH2D;
+class TH1F;
+class TProfile;
+class TTree;
+class IBeamCondSvc;
+
+namespace Trk{
+ class TrkVKalVrtFitter;
+ class IVertexFitter;
+ class IVKalState;
+}
+
+namespace MVAUtils { class BDT; }
+
+typedef std::vector<double> dvect;
+
+
+//------------------------------------------------------------------------
+namespace Rec {
+
+ struct workVectorArrxAOD{
+ std::vector<const xAOD::TrackParticle*> listSelTracks; // Selected tracks after quality cuts
+ std::vector<const xAOD::TrackParticle*> tmpListTracks;
+ std::vector<const xAOD::TrackParticle*> inpTrk; // All tracks provided to tool
+ double beamX=0.;
+ double beamY=0.;
+ double beamZ=0.;
+ double tanBeamTiltX=0.;
+ double tanBeamTiltY=0.;
+ };
+
+ class NewVrtSecInclusiveTool : public AthAlgTool, virtual public IVrtInclusive
+ {
+
+
+ public:
+ /* Constructor */
+ NewVrtSecInclusiveTool(const std::string& type, const std::string& name, const IInterface* parent);
+ /* Destructor */
+ virtual ~NewVrtSecInclusiveTool();
+
+
+ StatusCode initialize();
+ StatusCode finalize();
+
+
+
+ std::unique_ptr<Trk::VxSecVertexInfo> findAllVertices(const std::vector<const xAOD::TrackParticle*> & inputTracks,
+ const xAOD::Vertex & primaryVertex) const final;
+//------------------------------------------------------------------------------------------------------------------
+// Private data and functions
+//
+
+ private:
+
+ double m_w_1{};
+ TTree* m_tuple{};
+ TH1D* m_hb_massPiPi{};
+ TH1D* m_hb_massPiPi1{};
+ TH1D* m_hb_massPPi{};
+ TH1D* m_hb_massEE{};
+ TH1D* m_hb_nvrt2{};
+ TH1D* m_hb_ratio{};
+ TH1D* m_hb_totmass{};
+ TH1D* m_hb_impact{};
+ TH1D* m_hb_impactR{};
+ TH2D* m_hb_impactRZ{};
+ TH1D* m_hb_pileupRat{};
+ TH1D* m_hb_trkD0{};
+ TH1D* m_hb_trkZ{};
+ TH1F* m_hb_ntrksel{};
+ TH1F* m_hb_ntrkInput{};
+ TH1F* m_hb_trkSelect{};
+ TH1D* m_hb_impactZ{};
+ TH1D* m_hb_r2d{};
+ TH1D* m_hb_signif3D{};
+ TH1D* m_hb_impV0{};
+ TH1D* m_hb_sig3DTot{};
+ TH1F* m_hb_goodvrtN{};
+ TH1F* m_hb_goodvrt1N{};
+ TH1D* m_hb_distVV{};
+ TH1D* m_hb_diffPS{};
+ TH1D* m_hb_sig3D1tr{};
+ TH1D* m_hb_sig3D2tr{};
+ TH1D* m_hb_sig3DNtr{};
+ TH1F* m_hb_rawVrtN{};
+ TH1F* m_hb_cosSVMom{};
+ TH1F* m_hb_etaSV{};
+ TH1F* m_hb_fakeSVBDT{};
+//--
+
+ long int m_cutSctHits{};
+ long int m_cutPixelHits{};
+ long int m_cutSiHits{};
+ long int m_cutBLayHits{};
+ long int m_cutSharedHits{};
+ double m_cutPt{};
+ double m_cutZVrt{};
+ double m_cutA0{};
+ double m_cutChi2{};
+ double m_sel2VrtProbCut{};
+ double m_globVrtProbCut{};
+ double m_maxSVRadiusCut{};
+ double m_selVrtSigCut{};
+ double m_trkSigCut{};
+ float m_a0TrkErrorCut{};
+ float m_zTrkErrorCut{};
+ float m_VrtMassLimit{};
+ float m_Vrt2TrMassLimit{};
+ float m_Vrt2TrPtLimit{};
+ float m_antiPileupSigRCut{};
+ float m_dRdZRatioCut{};
+ float m_v2tIniBDTCut{};
+ float m_v2tFinBDTCut{};
+ float m_vertexMergeCut{};
+ float m_trackDetachCut{};
+ float m_beampipeR{};
+ float m_removeTrkMatSignif{};
+
+ bool m_fillHist{};
+ bool m_useVertexCleaning{};
+ bool m_multiWithOneTrkVrt{};
+ std::string m_calibFileName;
+
+ std::unique_ptr<MVAUtils::BDT> m_SV2T_BDT;
+
+ ServiceHandle< IBeamCondSvc > m_beamService;
+ //ToolHandle<Trk::IVertexFitter> m_fitterSvc;
+ ToolHandle<Trk::IExtrapolator> m_extrapolator{this,"ExtrapolatorName","Trk::Extrapolator/Extrapolator"};
+ ToolHandle<Trk::TrkVKalVrtFitter> m_fitSvc;
+ //Trk::TrkVKalVrtFitter* m_fitSvc{};
+
+ double m_massPi {};
+ double m_massP {};
+ double m_massE{};
+ double m_massK0{};
+ double m_massLam{};
+ std::string m_instanceName;
+
+//=======================================================================================
+// Functions and structure below are for algorithm development, debugging and calibration
+// NOT USED IN PRODUCTION!
+
+ int notFromBC(int PDGID) const;
+ const xAOD::TruthParticle * getPreviousParent(const xAOD::TruthParticle * child, int & ParentPDG) const;
+ int getIdHF(const xAOD::TrackParticle* TP ) const;
+ int getG4Inter( const xAOD::TrackParticle* TP ) const;
+ int getMCPileup(const xAOD::TrackParticle* TP ) const;
+
+ struct DevTuple
+ {
+ static const int maxNTrk=100;
+ static const int maxNVrt=100;
+ int nTrk;
+ float pttrk[maxNTrk];
+ float Sig3D[maxNTrk];
+ int n2Vrt;
+ int VrtTrkHF[maxNVrt];
+ int VrtTrkI[maxNVrt];
+ int VrtCh[maxNVrt];
+ int VrtIBL[maxNVrt];
+ int VrtBL[maxNVrt];
+ int VrtDisk[maxNVrt];
+ float VrtDist2D[maxNVrt];
+ float VrtSig3D[maxNVrt];
+ float VrtSig2D[maxNVrt];
+ float VrtM[maxNVrt];
+ float VrtZ[maxNVrt];
+ float VrtPt[maxNVrt];
+ float VrtEta[maxNVrt];
+ float VrtBDT[maxNVrt];
+ float VrtProb[maxNVrt];
+ float VrtHR1[maxNVrt];
+ float VrtHR2[maxNVrt];
+ float VrtSinSPM[maxNVrt];
+ float VMinPtT[maxNVrt];
+ float VMinS3DT[maxNVrt];
+ float VMaxS3DT[maxNVrt];
+ float VSigMat[maxNVrt];
+ //---
+ int nNVrt;
+ int NVrtTrk[maxNVrt];
+ int NVrtTrkHF[maxNVrt];
+ int NVrtTrkI[maxNVrt];
+ int NVrtCh[maxNVrt];
+ int NVrtIBL[maxNVrt];
+ int NVrtBL[maxNVrt];
+ float NVrtM[maxNVrt];
+ float NVrtPt[maxNVrt];
+ float NVrtEta[maxNVrt];
+ float NVrtSinSPM[maxNVrt];
+ float NVMinPtT[maxNVrt];
+ float NVMinS3DT[maxNVrt];
+ float NVMaxS3DT[maxNVrt];
+ float NVrtDist2D[maxNVrt];
+ float NVrtSig3D[maxNVrt];
+ float NVrtSig2D[maxNVrt];
+ float NVrtProb[maxNVrt];
+ float NVrtBDT[maxNVrt];
+ };
+ DevTuple* m_curTup;
+//
+// End of development stuff
+//============================================================
+
+
+ struct Vrt2Tr
+ {
+ int badVrt=0;
+ Amg::Vector3D fitVertex;
+ TLorentzVector momentum;
+ long int vertexCharge;
+ std::vector<double> errorMatrix;
+ std::vector<double> chi2PerTrk;
+ std::vector< std::vector<double> > trkAtVrt;
+ double chi2=0.;
+ double dRSVPV=-1.;
+ };
+
+
+// For multivertex version only
+
+ boost::adjacency_list<boost::listS, boost::vecS, boost::undirectedS> *m_compatibilityGraph{};
+ float m_chiScale[11]{};
+ struct WrkVrt
+ { bool Good=true;
+ std::deque<long int> selTrk;
+ Amg::Vector3D vertex;
+ TLorentzVector vertexMom;
+ long int vertexCharge{};
+ std::vector<double> vertexCov;
+ std::vector<double> chi2PerTrk;
+ std::vector< std::vector<double> > trkAtVrt;
+ double chi2{};
+ int nCloseVrt=0;
+ double dCloseVrt=1000000.;
+ double projectedVrt=0.;
+ int detachedTrack=-1;
+ };
+
+
+// Private technical functions
+//
+//
+ std::vector<xAOD::Vertex*> getVrtSecMulti( workVectorArrxAOD * inpParticlesxAOD, const xAOD::Vertex & primVrt) const;
+
+ void trackClassification(std::vector<WrkVrt> *WrkVrtSet,
+ std::vector< std::deque<long int> > *TrkInVrt) const;
+
+ double maxOfShared(std::vector<WrkVrt> *WrkVrtSet,
+ std::vector< std::deque<long int> > *TrkInVrt,
+ long int & selectedTrack,
+ long int & selectedVertex) const;
+ void removeTrackFromVertex(std::vector<WrkVrt> *WrkVrtSet,
+ std::vector< std::deque<long int> > *TrkInVrt,
+ long int selectedTrack,
+ long int selectedVertex) const;
+//
+//
+
+ void printWrkSet(const std::vector<WrkVrt> * WrkSet, const std::string &name ) const;
+
+//
+// Gives correct mass assignment in case of nonequal masses
+ double massV0(const std::vector< std::vector<double> >& TrkAtVrt, double massP, double massPi ) const;
+
+
+ TLorentzVector MomAtVrt(const std::vector<double>& inpTrk) const;
+ double VrtRadiusError(const Amg::Vector3D & secVrt, const std::vector<double> & vrtErr) const;
+
+ int nTrkCommon( std::vector<WrkVrt> *WrkVrtSet, int indexV1, int indexV2) const;
+ double minVrtVrtDist( std::vector<WrkVrt> *WrkVrtSet, int & indexV1, int & indexV2) const;
+ double minVrtVrtDistNext( std::vector<WrkVrt> *WrkVrtSet, int & indexV1, int & indexV2) const;
+ bool isPart( std::deque<long int> test, std::deque<long int> base) const;
+ void Clean1TrVertexSet(std::vector<WrkVrt> *WrkVrtSet) const;
+
+ double VrtVrtDist(const xAOD::Vertex & PrimVrt, const Amg::Vector3D & SecVrt,
+ const std::vector<double> VrtErr,double& Signif ) const;
+ double VrtVrtDist2D(const xAOD::Vertex & PrimVrt, const Amg::Vector3D & SecVrt,
+ const std::vector<double> VrtErr,double& Signif ) const;
+ double VrtVrtDist(const Amg::Vector3D & Vrt1, const std::vector<double>& VrtErr1,
+ const Amg::Vector3D & Vrt2, const std::vector<double>& VrtErr2) const;
+ double VrtVrtDist(const xAOD::Vertex & PrimVrt, const Amg::Vector3D & SecVrt,
+ const std::vector<double> SecVrtErr, const TLorentzVector & Dir) const;
+ double PntPntDist(const Amg::Vector3D & Vrt1, const Amg::Vector3D & Vrt2) const;
+
+ void DisassembleVertex(std::vector<WrkVrt> *WrkVrtSet, int iv,
+ std::vector<const xAOD::TrackParticle*> AllTracks,
+ Trk::IVKalState& istate) const;
+
+ double ProjSV_PV(const Amg::Vector3D & SV, const xAOD::Vertex & PV, const TLorentzVector & Direction) const;
+ double MomProjDist(const Amg::Vector3D & SV, const xAOD::Vertex & PV, const TLorentzVector & Direction) const;
+
+
+ double distToMatLayerSignificance(Vrt2Tr & Vrt) const;
+
+ StatusCode refitVertex( std::vector<WrkVrt> *WrkVrtSet, int SelectedVertex,
+ std::vector<const xAOD::TrackParticle*> & SelectedTracks,
+ Trk::IVKalState& istate,
+ bool ifCovV0) const;
+
+ double refitVertex( WrkVrt &Vrt,std::vector<const xAOD::TrackParticle*> & SelectedTracks,
+ Trk::IVKalState& istate,
+ bool ifCovV0) const;
+
+ int mostHeavyTrk(WrkVrt V, std::vector<const xAOD::TrackParticle*> AllTracks) const;
+
+ double mergeAndRefitVertices( std::vector<WrkVrt> *WrkVrtSet, int V1, int V2, WrkVrt & newvrt,
+ std::vector<const xAOD::TrackParticle*> & AllTrackList,
+ Trk::IVKalState& istate) const;
+ void mergeAndRefitOverlapVertices( std::vector<WrkVrt> *WrkVrtSet, int V1, int V2,
+ std::vector<const xAOD::TrackParticle*> & AllTrackLis,
+ Trk::IVKalState& istate) const;
+
+ double improveVertexChi2( std::vector<WrkVrt> *WrkVrtSet, int V, std::vector<const xAOD::TrackParticle*> & AllTracks,
+ Trk::IVKalState& istate,
+ bool ifCovV0) const;
+
+ void selGoodTrkParticle( workVectorArrxAOD * xAODwrk,
+ const xAOD::Vertex & PrimVrt) const;
+
+
+
+ void select2TrVrt(std::vector<const xAOD::TrackParticle*> & SelectedTracks,
+ const xAOD::Vertex & PrimVrt) const;
+
+
+ void getPixelDiscs (const xAOD::TrackParticle* Part, int &d0Hit, int &d1Hit, int &d2Hit) const;
+ int getIBLHit(const xAOD::TrackParticle* Part) const;
+ int getBLHit(const xAOD::TrackParticle* Part) const;
+
+ };
+
+
+ struct clique_visitor
+ {
+ clique_visitor(std::vector< std::vector<int> > & input): m_allCliques(input){ input.clear();}
+
+ template <typename Clique, typename Graph>
+ void clique(const Clique& clq, Graph& )
+ {
+ std::vector<int> new_clique(0);
+ for(auto i = clq.begin(); i != clq.end(); ++i) new_clique.push_back(*i);
+ m_allCliques.push_back(new_clique);
+ }
+
+ std::vector< std::vector<int> > & m_allCliques;
+
+ };
+
+
+} //end namespace
+
+#endif
diff --git a/Reconstruction/VKalVrt/NewVrtSecInclusiveTool/README b/Reconstruction/VKalVrt/NewVrtSecInclusiveTool/README
new file mode 100644
index 0000000000000000000000000000000000000000..fd65f5134f3eafe16e49c039be2f51d59b6c11f0
--- /dev/null
+++ b/Reconstruction/VKalVrt/NewVrtSecInclusiveTool/README
@@ -0,0 +1,28 @@
+A tool for inclusive secondary vertex reconstruction in ATLAS ID.
+
+This is an improved version of the VrtSecInclusive algorithm.
+
+Changes with respect to original version:
+
+1) Outdated Graph(CERNLIB) package is replaced by Bron&Kerbosch algorithm from BOOST
+2) Fake(fragmentation+pileup) and hadronic interaction vertices are removed using a dedicated MVA
+3) Additional material interaction rejection is based on the TrackingGeometry material layers
+4) Many improvements in vertex topology resolution
+5) One-track vertices are allowed (e.g. 3 track system, where one track crosses 2 others in different places,
+ so that 3-track vertex has inacceptable Chi2, will be resolved as 2-track + 1-track vertices).
+ This allows to get 2 vertices in case of e.g. cascade B->D decay with 3 reconstracted tracks in total.
+6) Algorithm is converted to Tool for convenience of usage in other Algorithms.
+7) Several tool configurations suitable for different problems are provided, e.g.:
+ SoftBFinderTool - targeting Soft b-tagging for SUSY
+ InclusiveBFinderTool - inclusive B-hadron decay reconstruction (a la CMS/D0) for b-tagging
+ DVFinderTool - search for displaced vertices produced by exotic particles
+
+
+Configured finders:
+1) SoftBFinderTool - looks for all low-pt (>1GeV) B-hadron vertices in event. No jet presence is assumed
+2) InclusiveBFinderTool - looks for all B-hadron vertices in event with low fake rate.
+ Suitable for ttbar and other processes with similar pt-range.
+3) HighPtBFinderTool - looks for all high-pt B-hadron vertices in event. No jet presence is assumed.
+4) MaterialSVFinderTool - looks for hadronic interaction vertices in detector material.
+5) DVFinderTool - looks for heavy displaced vertices in ID volume.
+
diff --git a/Reconstruction/VKalVrt/NewVrtSecInclusiveTool/python/NewVrtSecInclusive.py b/Reconstruction/VKalVrt/NewVrtSecInclusiveTool/python/NewVrtSecInclusive.py
new file mode 100644
index 0000000000000000000000000000000000000000..7d3c3cebf3f28fb7f46810638800bbaac1ed6fe0
--- /dev/null
+++ b/Reconstruction/VKalVrt/NewVrtSecInclusiveTool/python/NewVrtSecInclusive.py
@@ -0,0 +1,189 @@
+# Copyright (C) 2002-2020 CERN for the benefit of the ATLAS collaboration
+# Author: Vadim Kostyukhin vadim.kostyukhin@cern.ch
+from GaudiKernel.GaudiHandles import *
+from GaudiKernel.Proxy.Configurable import *
+from AthenaCommon.Include import Include, IncludeError, include
+from AthenaCommon.Logging import logging
+from AthenaCommon.AppMgr import ToolSvc
+from AthenaCommon import CfgMgr
+
+from NewVrtSecInclusiveTool.NewVrtSecInclusiveToolConf import Rec__NewVrtSecInclusiveTool
+
+# define the class
+# Search for soft B-hadron vertices.
+#------------------------------------
+class SoftBFinderTool( Rec__NewVrtSecInclusiveTool ):
+
+ def __init__(self, name = 'SoftBFinderTool' ):
+
+ from __main__ import ToolSvc
+ mlog = logging.getLogger( 'SoftBFinderTool::__init__ ' )
+ mlog.info("entering")
+ #----------------------
+ # VKalVrt vertex fitter
+ #
+ from TrkVKalVrtFitter.TrkVKalVrtFitterConf import Trk__TrkVKalVrtFitter
+ SVertexFitterTool = Trk__TrkVKalVrtFitter(name="SoftBVertexFitterTool",
+ Extrapolator="Trk::Extrapolator/AtlasExtrapolator")
+ ToolSvc += SVertexFitterTool
+ #----------------------
+ # Soft B-hadron vertex finder itself
+ #
+ Rec__NewVrtSecInclusiveTool.__init__( self, name = name,
+ VertexFitterTool = SVertexFitterTool,
+ CutBLayHits = 1,
+ CutPixelHits = 3,
+ CutSiHits = 8,
+ useVertexCleaning = True,
+ MultiWithOneTrkVrt = True,
+ removeTrkMatSignif = -1., # No additional material rejection
+ AntiPileupSigRCut = 2.0,
+ TrkSigCut = 2.0,
+ SelVrtSigCut = 3.0,
+ v2tIniBDTCut =-0.6,
+ v2tFinBDTCut = 0.2,
+ CutPt = 500,
+ MaxSVRadiusCut = 50 # Inside Pixel BL volume
+ )
+ mlog = logging.getLogger( 'SoftBFinderTool::__configured__ ' )
+
+##########################################################################################################
+# define the class
+class InclusiveBFinderTool( Rec__NewVrtSecInclusiveTool ):
+
+ def __init__(self, name = 'InclusiveBFinderTool' ):
+
+ from __main__ import ToolSvc
+ mlog = logging.getLogger( 'InclusiveBFinderTool::__init__ ' )
+ mlog.info("entering")
+ #----------------------
+ # VKalVrt vertex fitter
+ #
+ from TrkVKalVrtFitter.TrkVKalVrtFitterConf import Trk__TrkVKalVrtFitter
+ SVertexFitterTool = Trk__TrkVKalVrtFitter(name="InclsusiveBVertexFitterTool",
+ Extrapolator="Trk::Extrapolator/AtlasExtrapolator")
+ ToolSvc += SVertexFitterTool
+ #----------------------
+ # Soft B-hadron vertex finder itself
+ #
+ Rec__NewVrtSecInclusiveTool.__init__( self, name = name,
+ VertexFitterTool = SVertexFitterTool,
+ CutBLayHits = 0,
+ CutPixelHits = 2,
+ CutSiHits = 8,
+ useVertexCleaning = True,
+ MultiWithOneTrkVrt = True,
+ removeTrkMatSignif = -1., # No additional material rejection
+ AntiPileupSigRCut = 2.0,
+ TrkSigCut = 2.0,
+ SelVrtSigCut = 3.0,
+ v2tIniBDTCut =-0.5,
+ v2tFinBDTCut = 0.0,
+ CutPt = 500
+ )
+##########################################################################################################
+# define the class
+class HighPtBFinderTool( Rec__NewVrtSecInclusiveTool ):
+
+ def __init__(self, name = 'HighPtBFinderTool' ):
+
+ from __main__ import ToolSvc
+ mlog = logging.getLogger( 'HighPtBFinderTool::__init__ ' )
+ mlog.info("entering")
+ #----------------------
+ # VKalVrt vertex fitter
+ #
+ from TrkVKalVrtFitter.TrkVKalVrtFitterConf import Trk__TrkVKalVrtFitter
+ SVertexFitterTool = Trk__TrkVKalVrtFitter(name="HighPtBVertexFitterTool",
+ Extrapolator="Trk::Extrapolator/AtlasExtrapolator")
+ ToolSvc += SVertexFitterTool
+ #----------------------
+ # Soft B-hadron vertex finder itself
+ #
+ Rec__NewVrtSecInclusiveTool.__init__( self, name = name,
+ VertexFitterTool = SVertexFitterTool,
+ CutBLayHits = 0,
+ CutPixelHits = 2,
+ CutSiHits = 8,
+ useVertexCleaning = True,
+ MultiWithOneTrkVrt = True,
+ removeTrkMatSignif = -1., # No additional material rejection
+ AntiPileupSigRCut = 2.0,
+ TrkSigCut = 2.0,
+ SelVrtSigCut = 3.0,
+ v2tIniBDTCut =-0.6,
+ v2tFinBDTCut = 0.2,
+ CutPt = 1000
+ )
+##########################################################################################################
+# define the class
+class MaterialSVFinderTool( Rec__NewVrtSecInclusiveTool ):
+
+ def __init__(self, name = 'MaterialSVFinderTool' ):
+
+ from __main__ import ToolSvc
+ mlog = logging.getLogger( 'MaterialSVFinderTool::__init__ ' )
+ mlog.info("entering")
+ #----------------------
+ # VKalVrt vertex fitter
+ #
+ from TrkVKalVrtFitter.TrkVKalVrtFitterConf import Trk__TrkVKalVrtFitter
+ SVertexFitterTool = Trk__TrkVKalVrtFitter(name="MaterialSVVertexFitterTool",
+ Extrapolator="Trk::Extrapolator/AtlasExtrapolator")
+ ToolSvc += SVertexFitterTool
+ #----------------------
+ # Soft B-hadron vertex finder itself
+ #
+ Rec__NewVrtSecInclusiveTool.__init__( self, name = name,
+ VertexFitterTool = SVertexFitterTool,
+ CutBLayHits = 0,
+ CutPixelHits = 1,
+ CutSiHits = 8,
+ useVertexCleaning = False,
+ MultiWithOneTrkVrt = False,
+ removeTrkMatSignif = -1., # No additional material rejection
+ AntiPileupSigRCut = 5.0,
+ TrkSigCut = 5.0,
+ SelVrtSigCut = 10.0,
+ v2tIniBDTCut =-10., #Effectively remove MVA selection
+ v2tFinBDTCut =-10., #Effectively remove MVA selection
+ VrtMassLimit = 8000.,
+ Vrt2TrMassLimit = 8000.,
+ CutPt = 500.
+ )
+##########################################################################################################
+# define the class
+class DVFinderTool( Rec__NewVrtSecInclusiveTool ):
+
+ def __init__(self, name = 'DVFinderTool' ):
+
+ from __main__ import ToolSvc
+ mlog = logging.getLogger( 'DVFinderTool::__init__ ' )
+ mlog.info("entering")
+ #----------------------
+ # VKalVrt vertex fitter
+ #
+ from TrkVKalVrtFitter.TrkVKalVrtFitterConf import Trk__TrkVKalVrtFitter
+ DVertexFitterTool = Trk__TrkVKalVrtFitter(name="DVertexFitterTool",
+ Extrapolator="Trk::Extrapolator/AtlasExtrapolator")
+ ToolSvc += DVertexFitterTool
+ #----------------------
+ # Soft B-hadron vertex finder itself
+ #
+ Rec__NewVrtSecInclusiveTool.__init__( self, name = name,
+ VertexFitterTool = SVertexFitterTool,
+ CutBLayHits = 0,
+ CutPixelHits = 0,
+ CutSiHits = 0,
+ useVertexCleaning = False,
+ MultiWithOneTrkVrt = False,
+ removeTrkMatSignif = 2.,
+ AntiPileupSigRCut = 6.0,
+ TrkSigCut = 6.0,
+ SelVrtSigCut = 8.0,
+ v2tIniBDTCut =-0.5,
+ v2tFinBDTCut = 0.0,
+ VrtMassLimit = 50000.,
+ Vrt2TrMassLimit = 50000.,
+ CutPt = 1000.
+ )
diff --git a/Reconstruction/VKalVrt/NewVrtSecInclusiveTool/src/CutTrk.cxx b/Reconstruction/VKalVrt/NewVrtSecInclusiveTool/src/CutTrk.cxx
new file mode 100755
index 0000000000000000000000000000000000000000..122b6f0574f8e4cfeaa8a00b0afa3438b6661852
--- /dev/null
+++ b/Reconstruction/VKalVrt/NewVrtSecInclusiveTool/src/CutTrk.cxx
@@ -0,0 +1,97 @@
+/*
+ Copyright (C) 2002-2020 CERN for the benefit of the ATLAS collaboration
+*/
+///
+/// @author V.Kostykhin <Vadim.Kostyukhin@cern.ch>
+///
+
+// Header include
+#include "NewVrtSecInclusiveTool/NewVrtSecInclusiveTool.h"
+#include "AnalysisUtils/AnalysisMisc.h"
+#include "TrkVKalVrtFitter/TrkVKalVrtFitter.h"
+
+#include "TH1.h"
+
+//-------------------------------------------------
+namespace Rec{
+
+
+//==============================================================================================================
+// xAOD based stuff
+//
+ void NewVrtSecInclusiveTool::selGoodTrkParticle( workVectorArrxAOD * xAODwrk,
+ const xAOD::Vertex & PrimVrt)
+ const
+ {
+ std::vector<const xAOD::TrackParticle*>& inpTrk = xAODwrk->inpTrk;
+ std::vector<const xAOD::TrackParticle*>& selectedTracks = xAODwrk->listSelTracks;
+ std::vector<const xAOD::TrackParticle*>::const_iterator i_ntrk;
+ std::vector<double> impact,impactError;
+ std::multimap<double,const xAOD::TrackParticle*> orderedTrk;
+ if(m_fillHist){ m_hb_ntrkInput->Fill( inpTrk.size()+0.1, m_w_1);}
+ for (i_ntrk = inpTrk.begin(); i_ntrk < inpTrk.end(); ++i_ntrk) {
+//
+//-- Perigee in TrackParticle
+//
+ if(m_fillHist){ m_hb_trkSelect->Fill( 0., m_w_1);}
+ if((*i_ntrk)->numberDoF() == 0) continue; //Protection
+ double trkChi2 = (*i_ntrk)->chiSquared() / (*i_ntrk)->numberDoF();
+ if(trkChi2 > m_cutChi2) continue;
+ if( (*i_ntrk)->pt() < m_cutPt) continue;
+ if(m_fillHist){ m_hb_trkSelect->Fill( 1., m_w_1);}
+
+ const Trk::Perigee mPer=(*i_ntrk)->perigeeParameters() ;
+ double CovTrkMtx00 = (*(mPer.covariance()))(0,0);
+
+ uint8_t PixelHits,SctHits,BLayHits;
+ 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(BLayHits,xAOD::numberOfInnermostPixelLayerHits))) BLayHits=0;
+ if(m_fillHist){ m_hb_trkSelect->Fill( 2., m_w_1);}
+
+ uint8_t splSCTHits,outSCTHits,splPixHits,outPixHits;
+ if( !((*i_ntrk)->summaryValue(splSCTHits,xAOD::numberOfSCTSpoiltHits))) splSCTHits=0;
+ if( !((*i_ntrk)->summaryValue(outSCTHits,xAOD::numberOfSCTOutliers))) outSCTHits=0;
+ if( !((*i_ntrk)->summaryValue(splPixHits,xAOD::numberOfPixelSpoiltHits)))splPixHits=0;
+ if( !((*i_ntrk)->summaryValue(outPixHits,xAOD::numberOfPixelOutliers))) outPixHits=0;
+
+ Amg::Vector3D perigeePos=mPer.position();
+ double impactA0=sqrt( (perigeePos.x()-PrimVrt.x())*(perigeePos.x()-PrimVrt.x())
+ +(perigeePos.y()-PrimVrt.y())*(perigeePos.y()-PrimVrt.y()) );
+ double impactZ=perigeePos.z()-PrimVrt.z();
+ if(m_fillHist){
+ m_hb_trkD0->Fill( impactA0, m_w_1);
+ m_hb_trkZ ->Fill( impactZ, m_w_1);
+ }
+ if(fabs(impactZ)>m_cutZVrt) continue;
+ if(impactA0>m_cutA0) continue;
+ if(m_fillHist){ m_hb_trkSelect->Fill( 3., m_w_1);}
+
+ double bX=xAODwrk->beamX + (perigeePos.z()-xAODwrk->beamZ)*xAODwrk->tanBeamTiltX;
+ double bY=xAODwrk->beamY + (perigeePos.z()-xAODwrk->beamZ)*xAODwrk->tanBeamTiltY;
+ double impactBeam=sqrt( (perigeePos.x()-bX)*(perigeePos.x()-bX) + (perigeePos.y()-bY)*(perigeePos.y()-bY));
+//----Anti-pileup
+ double signifBeam = impactBeam / sqrt(CovTrkMtx00);
+ if(signifBeam < m_antiPileupSigRCut) continue; // cut against tracks from pileup vertices
+ if(m_fillHist){ m_hb_trkSelect->Fill( 4., m_w_1);}
+//----
+ if(PixelHits < m_cutPixelHits) continue;
+ if(SctHits < m_cutSctHits) continue;
+ if((PixelHits+SctHits) < m_cutSiHits) continue;
+ if(BLayHits < m_cutBLayHits) continue;
+ if(m_fillHist){ m_hb_trkSelect->Fill( 5., m_w_1);}
+//----
+ orderedTrk.emplace(signifBeam,*i_ntrk);
+ selectedTracks.push_back(*i_ntrk);
+ }
+//---- Order tracks according to pt
+// AnalysisUtils::Sort::pT (selectedTracks)
+//---- Order tracks according to ranks
+ std::map<double,const xAOD::TrackParticle*>::reverse_iterator rt=orderedTrk.rbegin();
+ selectedTracks.resize(orderedTrk.size());
+ for ( int cntt=0; rt!=orderedTrk.rend(); ++rt,++cntt) {selectedTracks[cntt]=(*rt).second;}
+ return;
+ }
+
+}//end namespace
diff --git a/Reconstruction/VKalVrt/NewVrtSecInclusiveTool/src/MultiUtilities.cxx b/Reconstruction/VKalVrt/NewVrtSecInclusiveTool/src/MultiUtilities.cxx
new file mode 100755
index 0000000000000000000000000000000000000000..7d7378aa8d226712c63b87a7d4da00dafdf0fdb5
--- /dev/null
+++ b/Reconstruction/VKalVrt/NewVrtSecInclusiveTool/src/MultiUtilities.cxx
@@ -0,0 +1,620 @@
+/*
+ Copyright (C) 2002-2020 CERN for the benefit of the ATLAS collaboration
+*/
+///
+/// @author Vadim Kostyukhin <vadim.kostyukhin@cern.ch>
+///
+// Header include
+#include "NewVrtSecInclusiveTool/NewVrtSecInclusiveTool.h"
+//-------------------------------------------------
+// Other stuff
+#include "AnalysisUtils/AnalysisMisc.h"
+#include "GeoPrimitives/GeoPrimitivesHelpers.h"
+#include "TrkVKalVrtFitter/TrkVKalVrtFitter.h"
+#include "TMath.h"
+#include <algorithm>
+//
+
+
+namespace Rec{
+
+
+ /* Service routines*/
+
+//-----------------------------------------------------------------------------------
+// Find track contributing most to the vertex invariant mass
+ int NewVrtSecInclusiveTool::mostHeavyTrk(WrkVrt V, std::vector<const xAOD::TrackParticle*> AllTracks)
+ const
+ {
+ int NTrk=V.selTrk.size(), SelT=-1;
+ if(NTrk<3)return -1;
+//--------------------------------------------------
+ TLorentzVector TSum;
+ for(int i=0; i<NTrk; i++){
+ TSum += AllTracks.at(V.selTrk[i])->p4();
+ }
+ double massMin=1.e9;
+ for(int i=0; i<NTrk; i++){
+ TLorentzVector TSum_m1 = TSum - AllTracks[V.selTrk[i]]->p4();
+ if(TSum_m1.M()<massMin){massMin=TSum_m1.M(); SelT=i;}
+ }
+ return SelT;
+ }
+//-----------------------------------------------------------------------------------
+// Detach the worst Chi2 track from the vertex and join it (if possible) with other track
+// from the vertex into additional 2-track vertices
+//
+ void NewVrtSecInclusiveTool::DisassembleVertex(std::vector<WrkVrt> *WrkVrtSet, int iv,
+ std::vector<const xAOD::TrackParticle*> AllTracks,
+ Trk::IVKalState& istate)
+ const
+ {
+ WrkVrt newvrt; newvrt.Good=true;
+ int NTrk=(*WrkVrtSet)[iv].selTrk.size(), SelT=-1;
+ if(NTrk<3)return;
+
+ std::vector<const xAOD::NeutralParticle*> neutralPartDummy(0);
+ std::vector<const xAOD::TrackParticle*> ListBaseTracks;
+ StatusCode sc; sc.setChecked();
+//=== To get robust definition of most bad outlier
+ m_fitSvc->setRobustness(5, istate);
+ sc = refitVertex( WrkVrtSet, iv, AllTracks, istate, false);
+ if(sc.isFailure()){ (*WrkVrtSet)[iv].Good=false; return; }
+ m_fitSvc->setRobustness(0, istate);
+//--------------------------------------------------
+ double Chi2Max=0.;
+ for(int i=0; i<NTrk; i++){
+ if( (*WrkVrtSet)[iv].chi2PerTrk[i]>Chi2Max) { Chi2Max=(*WrkVrtSet)[iv].chi2PerTrk[i]; SelT=i;}
+ }
+ int NVrtCur=WrkVrtSet->size();
+ for(int i=0; i<NTrk; i++){
+ if(i==SelT)continue;
+ ListBaseTracks.clear();
+ ListBaseTracks.push_back( AllTracks[(*WrkVrtSet)[iv].selTrk[i]] );
+ ListBaseTracks.push_back( AllTracks[(*WrkVrtSet)[iv].selTrk[SelT]] );
+ newvrt.selTrk.resize(2);
+ newvrt.selTrk[0]=(*WrkVrtSet)[iv].selTrk[i];
+ newvrt.selTrk[1]=(*WrkVrtSet)[iv].selTrk[SelT];
+ sc = m_fitSvc->VKalVrtFitFast(ListBaseTracks,newvrt.vertex,istate); /* Fast crude estimation*/
+ if( sc.isFailure() ) continue;
+ if( newvrt.vertex.perp() > m_maxSVRadiusCut ) newvrt.vertex=Amg::Vector3D(0.,0.,0.);
+ m_fitSvc->setApproximateVertex(newvrt.vertex[0],newvrt.vertex[1],newvrt.vertex[2],istate);
+ sc=m_fitSvc->VKalVrtFit(ListBaseTracks, neutralPartDummy,
+ newvrt.vertex,
+ newvrt.vertexMom,
+ newvrt.vertexCharge,
+ newvrt.vertexCov,
+ newvrt.chi2PerTrk,
+ newvrt.trkAtVrt,
+ newvrt.chi2,
+ istate, false);
+ if( sc.isFailure() ) continue;
+ if( newvrt.chi2>10.) continue; // Too bad 2-track vertex fit
+ newvrt.chi2PerTrk[0]=newvrt.chi2PerTrk[1]=newvrt.chi2/2.;
+ newvrt.nCloseVrt = 0;
+ newvrt.dCloseVrt = 1000000.;
+ newvrt.projectedVrt = 0.9999;
+ if((int)WrkVrtSet->size()==NVrtCur) { WrkVrtSet->push_back(newvrt); continue;} // just the first added vertex
+ if( (*WrkVrtSet).at(NVrtCur).chi2<newvrt.chi2 ) continue; // previously added 2tr vertex was better
+ WrkVrtSet->pop_back();
+ WrkVrtSet->push_back(newvrt);
+ }
+ (*WrkVrtSet)[iv].selTrk.erase( (*WrkVrtSet)[iv].selTrk.begin() + SelT ); //remove track
+ sc = refitVertex( WrkVrtSet, iv, AllTracks, istate, false);
+ if( sc.isFailure() ) {(*WrkVrtSet)[iv].Good=false; /*std::cout<<" Wrong vertex"<<'\n';*/}
+ }
+
+
+ void NewVrtSecInclusiveTool::Clean1TrVertexSet(std::vector<WrkVrt> *WrkVrtSet)
+ const
+ {
+ std::vector<int> countVT(WrkVrtSet->size(),0);
+ std::vector<int> linkedVrt(WrkVrtSet->size(),0);
+//--- Mark as bad the 1track vertex if the detached track is NOT present in any remaining good vertex (>=2tr)
+ for(int i1tv=0; i1tv<(int)WrkVrtSet->size(); i1tv++) {
+ if( (*WrkVrtSet)[i1tv].selTrk.size()!=1) continue;
+ if(!(*WrkVrtSet)[i1tv].Good) continue;
+ int Trk1=(*WrkVrtSet)[i1tv].detachedTrack;
+ int foundInGoodVrt=0;
+ for(int mtv=0; mtv<(int)WrkVrtSet->size(); mtv++) { //cycle over good vertices with many tracks
+ if( (*WrkVrtSet)[mtv].selTrk.size()<2) continue;
+ if(!(*WrkVrtSet)[mtv].Good) continue;
+ if( std::find((*WrkVrtSet)[mtv].selTrk.begin(),(*WrkVrtSet)[mtv].selTrk.end(), Trk1) != (*WrkVrtSet)[mtv].selTrk.end()){
+ //double m2Vrt=((*WrkVrtSet)[mtv].vertexMom+(*WrkVrtSet)[i1tv].vertexMom).M(); //VK Commented. M cut in other places
+ //if(m2Vrt>m_VrtMassLimit){ (*WrkVrtSet)[i1tv].Good=false; break; } //1Tr + manyTr system is too heavy
+ foundInGoodVrt++; countVT[mtv]++; linkedVrt[i1tv]=mtv; //Linked vertex found
+ }
+ }
+ if(!foundInGoodVrt)(*WrkVrtSet)[i1tv].Good=false; // Make the vertex bad
+ }
+//
+//---Select SINGLE 1tr-vertex from many pointing to one multi-track vertex
+ for(int mtv=0; mtv<(int)WrkVrtSet->size(); mtv++) {
+ if( (*WrkVrtSet)[mtv].selTrk.size()<2 ) continue;
+ if(!(*WrkVrtSet)[mtv].Good ) continue;
+ if( countVT[mtv] < 1 ) continue;
+ double distM=1.e9;
+ int best1TVrt=-1;
+ for(int i1tv=0; i1tv<(int)WrkVrtSet->size(); i1tv++) {
+ if( (*WrkVrtSet)[i1tv].selTrk.size()!=1 ) continue;
+ if(!(*WrkVrtSet)[i1tv].Good ) continue;
+ if( linkedVrt[i1tv] != mtv ) continue;
+ //double dist=((*WrkVrtSet)[mtv].vertex - (*WrkVrtSet)[i1tv].vertex).mag();
+ double dist=((*WrkVrtSet)[mtv].vertexMom+(*WrkVrtSet)[i1tv].vertexMom).M(); //Use
+ if( dist < distM ) {distM=dist; best1TVrt=i1tv;}
+ (*WrkVrtSet)[i1tv].Good=false;
+ }
+ if(best1TVrt>-1 && distM<m_VrtMassLimit)(*WrkVrtSet)[best1TVrt].Good=true;
+ }
+ }
+
+
+ void NewVrtSecInclusiveTool::trackClassification(std::vector<WrkVrt> *WrkVrtSet,
+ std::vector< std::deque<long int> > *TrkInVrt)
+ const
+ {
+ int NSet=WrkVrtSet->size();
+ for(int iv=0; iv<NSet; iv++) {
+ if(!(*WrkVrtSet)[iv].Good) continue;
+ int NTrkAtVrt=(*WrkVrtSet)[iv].selTrk.size();
+ if(NTrkAtVrt<2) continue;
+ for(int jt=0; jt<NTrkAtVrt; jt++){
+ int tracknum=(*WrkVrtSet)[iv].selTrk[jt];
+ (*TrkInVrt).at(tracknum).push_back(iv);
+ }
+ }
+ }
+
+
+ double NewVrtSecInclusiveTool::maxOfShared(std::vector<WrkVrt> *WrkVrtSet,
+ std::vector< std::deque<long int> > *TrkInVrt,
+ long int & SelectedTrack,
+ long int & SelectedVertex)
+ const
+ {
+ long int NTrack=TrkInVrt->size();
+ int it, jv, itmp, NVrt, VertexNumber;
+
+ double MaxOf=-999999999, Chi2Red=0., SelectedProb=1.;
+//
+ int NShareMax=0;
+ for( it=0; it<NTrack; it++) {
+ NVrt=(*TrkInVrt)[it].size(); /*Number of vertices for this track*/
+ if( NVrt > NShareMax ) NShareMax=NVrt;
+ }
+ if(NShareMax<=1)return MaxOf; /* No shared tracks */
+//
+ for( it=0; it<NTrack; it++) {
+ NVrt=(*TrkInVrt)[it].size(); /*Number of vertices for this track*/
+ if( NVrt <= 1 ) continue; /*Should ALWAYS be checked for safety*/
+ if( NVrt < NShareMax ) continue; /*Not a shared track with maximal sharing*/
+ double maxRadius=0;
+ for(auto &vrtn :(*TrkInVrt)[it] ){ maxRadius=std::max((*WrkVrtSet).at(vrtn).vertex.perp(),maxRadius); } //count number of 2-track vertices
+ for( jv=0; jv<NVrt; jv++ ){
+ VertexNumber=(*TrkInVrt)[it][jv];
+ if(!(*WrkVrtSet).at(VertexNumber).Good)continue;
+ int NTrkInVrt=(*WrkVrtSet)[VertexNumber].selTrk.size();
+ if( NTrkInVrt <= 1) continue; // one track vertex - nothing to do
+ ////if( N2trVrt>0 && N2trVrt<NShareMax && NTrkInVrt>2) continue; // Mixture of multi-track and 2-track vrts. Skip multi-track then.??? Bad!!!
+ for( itmp=0; itmp<NTrkInVrt; itmp++){
+ if( (*WrkVrtSet)[VertexNumber].selTrk[itmp] == it ) { /* Track found*/
+ //Chi2Red=(*WrkVrtSet)[VertexNumber].chi2PerTrk.at(itmp); // Normal Chi2
+ Chi2Red=(*WrkVrtSet)[VertexNumber].chi2PerTrk.at(itmp)*(*WrkVrtSet)[VertexNumber].vertex.perp()/maxRadius; // Works much better
+ double prob_vrt = TMath::Prob( (*WrkVrtSet)[VertexNumber].chi2, 2*(*WrkVrtSet)[VertexNumber].selTrk.size()-3);
+ if( MaxOf < Chi2Red ){
+ if(MaxOf>0 && prob_vrt>0.01 && SelectedProb<0.01 ) continue; // Don't disassemble good vertices if a bad one is present
+ MaxOf = Chi2Red;
+ SelectedTrack=it; SelectedVertex=VertexNumber;
+ SelectedProb = prob_vrt;
+ }
+ }
+ }
+ }
+ }
+//-- Additional check for a common track in 2tr-2tr/Ntr vertex topology
+ if( (*TrkInVrt)[SelectedTrack].size() == 2){
+ int v1=(*TrkInVrt)[SelectedTrack][0]; int v2=(*TrkInVrt)[SelectedTrack][1];
+ double prb1=TMath::Prob( (*WrkVrtSet)[v1].chi2, 2*(*WrkVrtSet)[v1].selTrk.size()-3),
+ prb2=TMath::Prob( (*WrkVrtSet)[v2].chi2, 2*(*WrkVrtSet)[v2].selTrk.size()-3);
+ double dst1=(*WrkVrtSet)[v1].projectedVrt, dst2=(*WrkVrtSet)[v2].projectedVrt;
+ if(prb1>0.05 && prb2>0.05){
+ if( (*WrkVrtSet)[v1].selTrk.size()==2 && (*WrkVrtSet)[v2].selTrk.size()==2){
+ if (SelectedVertex==v1 && dst2<dst1) SelectedVertex=v2; // Swap to remove the closest to PV vertex
+ else if(SelectedVertex==v2 && dst1<dst2) SelectedVertex=v1; // Swap to remove the closest to PV vertex
+ double M1=(*WrkVrtSet)[v1].vertexMom.M(); double M2=(*WrkVrtSet)[v2].vertexMom.M();
+ if( M1>m_VrtMassLimit && M2<m_VrtMassLimit ) SelectedVertex=v1;
+ if( M1<m_VrtMassLimit && M2>m_VrtMassLimit ) SelectedVertex=v2;
+ }
+ if( (*WrkVrtSet)[v1].selTrk.size()+(*WrkVrtSet)[v2].selTrk.size() > 4){
+ if( (*WrkVrtSet)[v1].selTrk.size()==2 && dst1>dst2) SelectedVertex=v2;
+ if( (*WrkVrtSet)[v2].selTrk.size()==2 && dst2>dst1) SelectedVertex=v1;
+ } }
+ }
+//
+ return MaxOf;
+ }
+
+
+ void NewVrtSecInclusiveTool::removeTrackFromVertex(std::vector<WrkVrt> *WrkVrtSet,
+ std::vector< std::deque<long int> > *TrkInVrt,
+ long int selectedTrack,
+ long int selectedVertex)
+ const
+ {
+ int posInVrtFit=0; //Position of SelectedTrack in vertex fit track list.
+ std::deque<long int>::iterator it;
+ WrkVrt & CVRT=(*WrkVrtSet).at(selectedVertex);
+//std::cout<<" In Found ="<<selectedTrack<<", "<<selectedVertex<<'\n';
+ for(it=CVRT.selTrk.begin();
+ it!=CVRT.selTrk.end(); it++) {
+ if( (*it) == selectedTrack ) {
+ CVRT.selTrk.erase(it); break;
+ }
+ posInVrtFit++;
+ }
+
+ for(it=(*TrkInVrt).at(selectedTrack).begin();
+ it!=(*TrkInVrt)[selectedTrack].end(); it++) {
+ if( (*it) == selectedVertex ) {
+ (*TrkInVrt)[selectedTrack].erase(it); break;
+ }
+ }
+ (*WrkVrtSet)[selectedVertex].detachedTrack=selectedTrack;
+//Check if track is removed from 2tr vertex => then sharing of track left should also be decreased
+ if( CVRT.selTrk.size() == 1){
+ long int LeftTrack=CVRT.selTrk[0]; // track left in 1tr vertex
+ for(it=(*TrkInVrt).at(LeftTrack).begin(); it!=(*TrkInVrt)[LeftTrack].end(); it++) {
+ if( (*it) == selectedVertex ) {
+ (*TrkInVrt)[LeftTrack].erase(it); break;
+ }
+ }
+ if( TMath::Prob( CVRT.chi2, 1) < 0.05 ) CVRT.Good=false; // Not really good Chi2 for one-track vertex
+ if( CVRT.vertexMom.M()>m_VrtMassLimit)CVRT.Good=false; // Vertex is too heavy
+ int ipos=0; if(posInVrtFit==0)ipos=1; // Position of remaining track in previous 2track vertex fit
+ CVRT.vertexMom=MomAtVrt(CVRT.trkAtVrt[ipos]); //Redefine vertexMom using remaining track
+ if((*TrkInVrt)[LeftTrack].size()>0)CVRT.Good=false; //Vertex is declared false only if remaining track
+ // is still linked to another vertex
+ if(CVRT.vertexMom.Pt()<2.*m_cutPt) CVRT.Good=false; //1track vertex should have twice momentum of track pt cut
+ CVRT.trkAtVrt.erase((*WrkVrtSet)[selectedVertex].trkAtVrt.begin()+posInVrtFit); //Remove also TrkAtVrt
+ }
+
+ }
+//
+// Number of common tracks for 2 vertices
+//
+ int NewVrtSecInclusiveTool::nTrkCommon( std::vector<WrkVrt> *WrkVrtSet,
+ int V1, int V2)
+ const
+ {
+ int NTrk_V1 = (*WrkVrtSet).at(V1).selTrk.size(); if( NTrk_V1< 2) return 0; /* Bad vertex */
+ int NTrk_V2 = (*WrkVrtSet).at(V2).selTrk.size(); if( NTrk_V2< 2) return 0; /* Bad vertex */
+ int nTrkCom=0;
+ if(NTrk_V1 < NTrk_V2){
+ for(int i=0; i<NTrk_V1; i++){
+ int trk=(*WrkVrtSet)[V1].selTrk[i];
+ if( std::find((*WrkVrtSet)[V2].selTrk.begin(),(*WrkVrtSet)[V2].selTrk.end(),trk) != (*WrkVrtSet)[V2].selTrk.end()) nTrkCom++;
+ }
+ }else{
+ for(int i=0; i<NTrk_V2; i++){
+ int trk=(*WrkVrtSet)[V2].selTrk[i];
+ if( std::find((*WrkVrtSet)[V1].selTrk.begin(),(*WrkVrtSet)[V1].selTrk.end(),trk) != (*WrkVrtSet)[V1].selTrk.end()) nTrkCom++;
+ }
+ }
+ return nTrkCom;
+ }
+//
+// Minimal normalized vertex-vertex distance
+//
+ double NewVrtSecInclusiveTool::minVrtVrtDist( std::vector<WrkVrt> *WrkVrtSet, int & V1, int & V2)
+ const
+ {
+ V1=V2=0;
+ for(int iv=0; iv<(int)WrkVrtSet->size(); iv++) { (*WrkVrtSet).at(iv).dCloseVrt=1000000.; (*WrkVrtSet)[iv].nCloseVrt=-1;}
+ double foundMinVrtDst=1000000.;
+
+ for(int iv=0; iv<(int)WrkVrtSet->size()-1; iv++) {
+ if( (*WrkVrtSet)[iv].selTrk.size()< 2) continue; /* Bad vertices */
+ if(!(*WrkVrtSet)[iv].Good ) continue; /* Bad vertices */
+ for(int jv=iv+1; jv<(int)WrkVrtSet->size(); jv++) {
+ if( (*WrkVrtSet)[jv].selTrk.size()< 2) continue; /* Bad vertices */
+ if(!(*WrkVrtSet)[jv].Good ) continue; /* Bad vertices */
+ //if(!nTrkCommon(WrkVrtSet, V1, V2)) continue; /* No common tracks*/
+ double tmp= fabs((*WrkVrtSet)[iv].vertex.x()-(*WrkVrtSet)[jv].vertex.x())
+ +fabs((*WrkVrtSet)[iv].vertex.y()-(*WrkVrtSet)[jv].vertex.y())
+ +fabs((*WrkVrtSet)[iv].vertex.z()-(*WrkVrtSet)[jv].vertex.z());
+ if(tmp>20.) continue;
+ double tmpDst = VrtVrtDist((*WrkVrtSet)[iv].vertex,(*WrkVrtSet)[iv].vertexCov,
+ (*WrkVrtSet)[jv].vertex,(*WrkVrtSet)[jv].vertexCov);
+ if( tmpDst < foundMinVrtDst){foundMinVrtDst = tmpDst; V1=iv; V2=jv;}
+ if( tmpDst < (*WrkVrtSet)[iv].dCloseVrt ) {(*WrkVrtSet)[iv].dCloseVrt=tmpDst; (*WrkVrtSet)[iv].nCloseVrt=jv;}
+ if( tmpDst < (*WrkVrtSet)[jv].dCloseVrt ) {(*WrkVrtSet)[jv].dCloseVrt=tmpDst; (*WrkVrtSet)[jv].nCloseVrt=iv;}
+ }
+ }
+ return foundMinVrtDst;
+ }
+//
+// Give minimal distance between nonmodifed yet vertices
+//
+ double NewVrtSecInclusiveTool::minVrtVrtDistNext( std::vector<WrkVrt> *WrkVrtSet, int & V1, int & V2)
+ const
+ {
+ V1=0; V2=0;
+ double foundMinVrtDst=1000000.;
+ for(int iv=0; iv<(int)WrkVrtSet->size()-1; iv++) {
+ if( (*WrkVrtSet)[iv].selTrk.size()< 2) continue; /* Bad vertex */
+ if( (*WrkVrtSet)[iv].nCloseVrt==-1) continue; /* Used vertex */
+ if( (*WrkVrtSet)[iv].dCloseVrt <foundMinVrtDst ) {
+ int vtmp=(*WrkVrtSet)[iv].nCloseVrt;
+ if((*WrkVrtSet)[vtmp].nCloseVrt==-1) { continue;} // Close vertex to given [iv] one is modified already
+ foundMinVrtDst=(*WrkVrtSet)[iv].dCloseVrt;
+ V1=iv;
+ V2=vtmp;
+ }
+ }
+ return foundMinVrtDst;
+ }
+
+//
+// Check two close vertices by explicit vertex fit and create combined vertex if successful
+//
+ double NewVrtSecInclusiveTool::mergeAndRefitVertices( std::vector<WrkVrt> *WrkVrtSet, int V1, int V2, WrkVrt & newvrt,
+ std::vector<const xAOD::TrackParticle*> & AllTrackList,
+ Trk::IVKalState& istate)
+ const
+ {
+ if(!(*WrkVrtSet).at(V1).Good)return -1.; //bad vertex
+ if(!(*WrkVrtSet).at(V2).Good)return -1.; //bad vertex
+ std::vector<const xAOD::NeutralParticle*> neutralPartDummy(0);
+ newvrt.Good=true;
+ int NTrk_V1=(*WrkVrtSet)[V1].selTrk.size();
+ int NTrk_V2=(*WrkVrtSet)[V2].selTrk.size();
+ newvrt.selTrk.resize(NTrk_V1+NTrk_V2);
+ std::copy((*WrkVrtSet)[V1].selTrk.begin(),(*WrkVrtSet)[V1].selTrk.end(), newvrt.selTrk.begin());
+ std::copy((*WrkVrtSet)[V2].selTrk.begin(),(*WrkVrtSet)[V2].selTrk.end(), newvrt.selTrk.begin()+NTrk_V1);
+
+ std::deque<long int>::iterator TransfEnd ;
+ sort( newvrt.selTrk.begin(), newvrt.selTrk.end() );
+ TransfEnd = unique( newvrt.selTrk.begin(), newvrt.selTrk.end() );
+ newvrt.selTrk.erase( TransfEnd, newvrt.selTrk.end());
+ std::vector<const xAOD::TrackParticle*> fitTrackList(0);
+ for(int it=0; it<(int)newvrt.selTrk.size(); it++)fitTrackList.push_back( AllTrackList[newvrt.selTrk[it]] );
+ m_fitSvc->setApproximateVertex( 0.5*((*WrkVrtSet)[V1].vertex[0]+(*WrkVrtSet)[V2].vertex[0]),
+ 0.5*((*WrkVrtSet)[V1].vertex[1]+(*WrkVrtSet)[V2].vertex[1]),
+ 0.5*((*WrkVrtSet)[V1].vertex[2]+(*WrkVrtSet)[V2].vertex[2]),
+ istate);
+
+ StatusCode sc=m_fitSvc->VKalVrtFit(fitTrackList,neutralPartDummy,
+ newvrt.vertex, newvrt.vertexMom, newvrt.vertexCharge, newvrt.vertexCov,
+ newvrt.chi2PerTrk, newvrt.trkAtVrt, newvrt.chi2,
+ istate, false);
+ if( sc.isFailure() ) return -1.;
+ if( newvrt.chi2>500. ) return -1.; //VK protection
+ if( newvrt.chi2PerTrk.size()==2) newvrt.chi2PerTrk[0]=newvrt.chi2PerTrk[1]=newvrt.chi2/2.;
+ return TMath::Prob( newvrt.chi2, 2*newvrt.selTrk.size()-3);
+ }
+ //
+ //
+ //================== Intelligent merge of multitrack vertices with 2 and more common tracks
+ void NewVrtSecInclusiveTool::mergeAndRefitOverlapVertices( std::vector<WrkVrt> *WrkVrtSet, int V1, int V2,
+ std::vector<const xAOD::TrackParticle*> & AllTrackList,
+ Trk::IVKalState& istate) const
+ {
+ if(!(*WrkVrtSet).at(V1).Good)return; //bad vertex
+ if(!(*WrkVrtSet).at(V2).Good)return; //bad vertex
+ std::vector<const xAOD::NeutralParticle*> neutralPartDummy(0);
+ WrkVrt newvrt;
+ newvrt.Good=true;
+ if( nTrkCommon( WrkVrtSet, V1, V2)<2 )return; //No overlap
+ //- V1 should become ref. vertex. Another Vrt tracks will be added to it.
+ if( (*WrkVrtSet)[V1].selTrk.size() < (*WrkVrtSet)[V2].selTrk.size() )
+ {int itmp=V2; V2=V1; V1=itmp;} //Vertex with NTrk=max is chosen
+ else if( (*WrkVrtSet)[V1].selTrk.size() == (*WrkVrtSet)[V2].selTrk.size() )
+ { if( (*WrkVrtSet)[V1].chi2 > (*WrkVrtSet)[V2].chi2 ) {int itmp=V2; V2=V1; V1=itmp;} } // Vertex with minimal Chi2 is chosen
+ //- Fill base track list for new vertex
+ newvrt.selTrk.resize( (*WrkVrtSet)[V1].selTrk.size() );
+ std::copy((*WrkVrtSet)[V1].selTrk.begin(),(*WrkVrtSet)[V1].selTrk.end(), newvrt.selTrk.begin());
+ //- Identify non-common tracks in second vertex
+ std::vector<int> noncommonTrk(0);
+ for(auto &it : (*WrkVrtSet)[V2].selTrk){
+ if( std::find((*WrkVrtSet)[V1].selTrk.begin(), (*WrkVrtSet)[V1].selTrk.end(), it) == (*WrkVrtSet)[V1].selTrk.end() )
+ noncommonTrk.push_back(it);
+ }
+ //
+ // Try to add non-common tracks one by one
+ std::vector<const xAOD::TrackParticle*> fitTrackList(0);
+ std::vector<int> detachedTrk(0);
+ StatusCode sc; sc.setChecked();
+ WrkVrt bestVrt;
+ bool foundMerged=false;
+ for( auto nct : noncommonTrk){
+ fitTrackList.clear();
+ for(int it=0; it<(int)newvrt.selTrk.size(); it++)fitTrackList.push_back( AllTrackList[newvrt.selTrk[it]] );
+ fitTrackList.push_back( AllTrackList.at(nct) );
+ m_fitSvc->setApproximateVertex( (*WrkVrtSet)[V1].vertex[0],(*WrkVrtSet)[V1].vertex[1],(*WrkVrtSet)[V1].vertex[2],istate);
+ sc=m_fitSvc->VKalVrtFit(fitTrackList, neutralPartDummy, newvrt.vertex, newvrt.vertexMom, newvrt.vertexCharge, newvrt.vertexCov,
+ newvrt.chi2PerTrk, newvrt.trkAtVrt, newvrt.chi2,
+ istate, false);
+ if( sc.isFailure() || TMath::Prob( newvrt.chi2, 2*newvrt.selTrk.size()+2-3)<m_globVrtProbCut ) {
+ detachedTrk.push_back(nct);
+ continue;
+ }
+ newvrt.selTrk.push_back(nct); // Compatible track. Add to common vertex.
+ bestVrt=newvrt;
+ foundMerged=true;
+ }
+ if(foundMerged)(*WrkVrtSet)[V1]=bestVrt;
+ (*WrkVrtSet)[V2].Good=false;
+ //
+ // Now detached tracks
+ if(detachedTrk.size()>1){
+ WrkVrt nVrt;
+ fitTrackList.clear(); nVrt.selTrk.clear();
+ for(auto nct : detachedTrk){ fitTrackList.push_back( AllTrackList[nct] ); nVrt.selTrk.push_back(nct); }
+ m_fitSvc->setApproximateVertex( (*WrkVrtSet)[V2].vertex[0],(*WrkVrtSet)[V2].vertex[1],(*WrkVrtSet)[V2].vertex[2],istate);
+ sc=m_fitSvc->VKalVrtFit(fitTrackList, neutralPartDummy, nVrt.vertex, nVrt.vertexMom, nVrt.vertexCharge, nVrt.vertexCov,
+ nVrt.chi2PerTrk, nVrt.trkAtVrt, nVrt.chi2,
+ istate, false);
+ if(sc.isSuccess()) (*WrkVrtSet).push_back(nVrt);
+ } else if( detachedTrk.size()==1 ){
+ bool tFound=false;
+ for( auto &vrt : (*WrkVrtSet) ){
+ if( !vrt.Good || vrt.selTrk.size()<2 )continue;
+ if( std::find(vrt.selTrk.begin(), vrt.selTrk.end(), detachedTrk[0]) != vrt.selTrk.end() ) { tFound=true; break; }
+ }
+ if( !tFound ) { //Track is not present in other vertices.
+ double Chi2min=1.e9; int selectedTrk=-1;
+ WrkVrt saveVrt;
+ fitTrackList.resize(2);
+ fitTrackList[0]= AllTrackList[detachedTrk[0]];
+ for(auto trk : (*WrkVrtSet)[V2].selTrk){
+ if(trk==detachedTrk[0])continue;
+ WrkVrt nVrt; nVrt.Good=true;
+ fitTrackList[1]=AllTrackList[trk];
+ m_fitSvc->setApproximateVertex( (*WrkVrtSet)[V2].vertex[0],(*WrkVrtSet)[V2].vertex[1],(*WrkVrtSet)[V2].vertex[2],istate);
+ sc=m_fitSvc->VKalVrtFit(fitTrackList, neutralPartDummy, nVrt.vertex, nVrt.vertexMom, nVrt.vertexCharge, nVrt.vertexCov,
+ nVrt.chi2PerTrk, nVrt.trkAtVrt, nVrt.chi2,
+ istate, false);
+ if(sc.isSuccess() && nVrt.chi2<Chi2min) {Chi2min=nVrt.chi2; saveVrt=nVrt; selectedTrk=trk; }
+ }
+ if(Chi2min<1.e9){
+ saveVrt.selTrk.resize(1); saveVrt.selTrk[0]=detachedTrk[0];
+ saveVrt.detachedTrack=selectedTrk;
+ saveVrt.vertexMom=MomAtVrt(saveVrt.trkAtVrt[0]); //redefine vertex momentum
+ (*WrkVrtSet).push_back(saveVrt);
+ }
+ }
+ }
+ return ;
+ }
+
+//
+// Iterate track removal until vertex get good Chi2
+//
+
+ double NewVrtSecInclusiveTool::improveVertexChi2( std::vector<WrkVrt> *WrkVrtSet, int V, std::vector<const xAOD::TrackParticle*> & AllTrackList,
+ Trk::IVKalState& istate,
+ bool ifCovV0)
+ const
+ {
+
+ int NTrk=(*WrkVrtSet)[V].selTrk.size();
+ if( NTrk<2 )return 0.;
+ double Prob=TMath::Prob( (*WrkVrtSet)[V].chi2, 2*NTrk-3);
+ //
+ //----Start track removal iterations
+ while(Prob<0.01){
+ NTrk=(*WrkVrtSet)[V].selTrk.size();
+ if(NTrk==2)return Prob;
+ int SelT=-1; double Chi2Max=0.;
+ for(int i=0; i<NTrk; i++){
+ if( (*WrkVrtSet)[V].chi2PerTrk[i]>Chi2Max) {
+ Chi2Max=(*WrkVrtSet)[V].chi2PerTrk[i];
+ SelT=i;
+ }
+ }
+ if (SelT<0) return 0;
+ (*WrkVrtSet)[V].detachedTrack=(*WrkVrtSet)[V].selTrk[SelT];
+ (*WrkVrtSet)[V].selTrk.erase( (*WrkVrtSet)[V].selTrk.begin() + SelT ); //remove track
+ StatusCode sc = refitVertex( WrkVrtSet, V, AllTrackList, istate, ifCovV0);
+ if(sc.isFailure())return 0.;
+ Prob=TMath::Prob( (*WrkVrtSet)[V].chi2, 2*(NTrk-1)-3);
+ }
+ return Prob;
+ }
+
+
+
+
+ StatusCode NewVrtSecInclusiveTool::refitVertex( std::vector<WrkVrt> *WrkVrtSet,
+ int SelectedVertex,
+ std::vector<const xAOD::TrackParticle*> & selectedTracks,
+ Trk::IVKalState& istate,
+ bool ifCovV0) const
+ {
+ int nth = (*WrkVrtSet)[SelectedVertex].selTrk.size();
+
+ if(nth<2) return StatusCode::SUCCESS;
+
+ std::vector<const xAOD::NeutralParticle*> neutralPartDummy(0);
+ std::vector<const xAOD::TrackParticle*> ListTracks(0);
+ for(int i=0;i<nth;i++) {
+ int j=(*WrkVrtSet)[SelectedVertex].selTrk[i]; /*Track number*/
+ ListTracks.push_back( selectedTracks[j] );
+ }
+ (*WrkVrtSet)[SelectedVertex].Good = false;
+ (*WrkVrtSet)[SelectedVertex].chi2PerTrk.resize(nth);
+ for(int i=0;i<nth;i++)(*WrkVrtSet)[SelectedVertex].chi2PerTrk[i]=100000.+i; //VK safety
+
+ m_fitSvc->setApproximateVertex((*WrkVrtSet)[SelectedVertex].vertex[0],
+ (*WrkVrtSet)[SelectedVertex].vertex[1],
+ (*WrkVrtSet)[SelectedVertex].vertex[2],
+ istate);
+ StatusCode sc=m_fitSvc->VKalVrtFit(ListTracks, neutralPartDummy,
+ (*WrkVrtSet)[SelectedVertex].vertex,
+ (*WrkVrtSet)[SelectedVertex].vertexMom,
+ (*WrkVrtSet)[SelectedVertex].vertexCharge,
+ (*WrkVrtSet)[SelectedVertex].vertexCov,
+ (*WrkVrtSet)[SelectedVertex].chi2PerTrk,
+ (*WrkVrtSet)[SelectedVertex].trkAtVrt,
+ (*WrkVrtSet)[SelectedVertex].chi2,
+ istate, ifCovV0);
+ if(sc.isSuccess())(*WrkVrtSet)[SelectedVertex].Good = true;
+ if((*WrkVrtSet)[SelectedVertex].chi2PerTrk.size()==2)
+ (*WrkVrtSet)[SelectedVertex].chi2PerTrk[0]=
+ (*WrkVrtSet)[SelectedVertex].chi2PerTrk[1]=(*WrkVrtSet)[SelectedVertex].chi2/2.;
+ return sc;
+ }
+
+
+
+ double NewVrtSecInclusiveTool::refitVertex(WrkVrt &Vrt,
+ std::vector<const xAOD::TrackParticle*> & selectedTracks,
+ Trk::IVKalState& istate,
+ bool ifCovV0) const
+ {
+ int i,j;
+ int nth = Vrt.selTrk.size();
+
+ if(nth<2) return -1.;
+
+ std::vector<const xAOD::NeutralParticle*> neutralPartDummy(0);
+ std::vector<const xAOD::TrackParticle*> ListTracks(0);
+ for(i=0;i<nth;i++) {
+ j=Vrt.selTrk[i]; /*Track number*/
+ ListTracks.push_back( selectedTracks[j] );
+ }
+ Vrt.Good = false;
+ Vrt.chi2PerTrk.resize(nth);
+ for(i=0;i<nth;i++)Vrt.chi2PerTrk[i]=100000.+i; //VK safety
+
+ m_fitSvc->setApproximateVertex(Vrt.vertex[0],Vrt.vertex[1],Vrt.vertex[2],istate);
+ StatusCode sc=m_fitSvc->VKalVrtFit(ListTracks,neutralPartDummy,Vrt.vertex,Vrt.vertexMom,Vrt.vertexCharge,
+ Vrt.vertexCov,Vrt.chi2PerTrk, Vrt.trkAtVrt,Vrt.chi2,
+ istate, ifCovV0);
+ if(sc.isSuccess())Vrt.Good = true;
+ else {Vrt.Good = false; return -1.;}
+ if(Vrt.chi2PerTrk.size()==2) Vrt.chi2PerTrk[0]=Vrt.chi2PerTrk[1]=Vrt.chi2/2.;
+ return TMath::Prob( Vrt.chi2, 2*nth-1);
+ }
+
+
+ bool NewVrtSecInclusiveTool::isPart( std::deque<long int> test, std::deque<long int> base)
+ const
+ {
+ std::deque<long int>::iterator trk=test.begin();
+ for(trk=test.begin(); trk!=test.end(); trk++)
+ if(std::find(base.begin(), base.end(), (*trk)) == base.end()) return false; //element not found => test is not part of base
+ return true;
+ }
+
+ double NewVrtSecInclusiveTool::MomProjDist(const Amg::Vector3D & SecVrt,const xAOD::Vertex &PrimVrt,const TLorentzVector &Mom)
+ const
+ {
+ Amg::Vector3D vv=SecVrt-PrimVrt.position();
+ return ( vv.x()*Mom.X()+vv.y()*Mom.Y()+vv.z()*Mom.Z() )/ Mom.P();
+ }
+
+} //end namespace
+
+
+
diff --git a/Reconstruction/VKalVrt/NewVrtSecInclusiveTool/src/NewVrtSecInclusiveTool.cxx b/Reconstruction/VKalVrt/NewVrtSecInclusiveTool/src/NewVrtSecInclusiveTool.cxx
new file mode 100644
index 0000000000000000000000000000000000000000..5b18338fc9941b66692c1097857d6e8ade31b29a
--- /dev/null
+++ b/Reconstruction/VKalVrt/NewVrtSecInclusiveTool/src/NewVrtSecInclusiveTool.cxx
@@ -0,0 +1,345 @@
+/*
+ Copyright (C) 2002-2020 CERN for the benefit of the ATLAS collaboration
+*/
+///
+/// @author Vadim Kostyukhin <vadim.kostyukhin@cern.ch>
+///
+
+
+// Header include
+#include "NewVrtSecInclusiveTool/NewVrtSecInclusiveTool.h"
+#include "VxSecVertex/VxSecVertexInfo.h"
+#include "TrkVKalVrtFitter/TrkVKalVrtFitter.h"
+#include "InDetBeamSpotService/IBeamCondSvc.h"
+#include "PathResolver/PathResolver.h"
+
+#include "GaudiKernel/ITHistSvc.h"
+#include "TH1.h"
+#include "TH2.h"
+#include "TTree.h"
+#include "TMath.h"
+#include "TFile.h"
+#include "MVAUtils/BDT.h"
+//
+
+
+namespace Rec {
+
+//
+//Constructor--------------------------------------------------------------
+NewVrtSecInclusiveTool::NewVrtSecInclusiveTool(const std::string& type,
+ const std::string& name,
+ const IInterface* parent):
+ AthAlgTool(type,name,parent),
+ m_cutSctHits(4),
+ m_cutPixelHits(2),
+ m_cutSiHits(8),
+ m_cutBLayHits(0),
+ m_cutSharedHits(1),
+ m_cutPt(500.),
+ m_cutZVrt(15.),
+ m_cutA0(10.),
+ m_cutChi2(5.),
+ m_sel2VrtProbCut(0.02),
+ m_globVrtProbCut(0.005),
+ m_maxSVRadiusCut(140.),
+ m_selVrtSigCut(3.0),
+ m_trkSigCut(2.0),
+ m_a0TrkErrorCut(1.0),
+ m_zTrkErrorCut(5.0),
+ m_VrtMassLimit(5500.),
+ m_Vrt2TrMassLimit(4000.),
+ m_Vrt2TrPtLimit(5.e5),
+ m_antiPileupSigRCut(2.0),
+ m_dRdZRatioCut(0.25),
+ m_v2tIniBDTCut(-0.6),
+ m_v2tFinBDTCut(0.),
+ m_vertexMergeCut(3.),
+ m_trackDetachCut(6.),
+ m_beampipeR(24.3),
+ m_removeTrkMatSignif(0.),
+ m_fillHist(false),
+ m_useVertexCleaning(true),
+ m_multiWithOneTrkVrt(true),
+ m_calibFileName("Fake2TrVertexReject.MVA.v01.root"),
+ m_SV2T_BDT(nullptr),
+ m_beamService("BeamCondSvc",name),
+ m_fitSvc("Trk::TrkVKalVrtFitter/VertexFitterTool",this)
+ {
+//
+// Declare additional interface
+//
+ declareInterface< IVrtInclusive >(this);
+// Properties
+//
+//
+ 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("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" );
+
+ declareProperty("CutPt", m_cutPt, "Track Pt selection cut" );
+ declareProperty("CutA0", m_cutA0, "Track A0 selection cut" );
+ declareProperty("CutZVrt", m_cutZVrt, "Track Z impact selection cut");
+ declareProperty("CutChi2", m_cutChi2, "Track Chi2 selection cut" );
+ declareProperty("TrkSigCut", m_trkSigCut, "Track 3D impact significance w/r primary vertex. Should be >=AntiPileupSigRCut" );
+
+ declareProperty("A0TrkErrorCut", m_a0TrkErrorCut, "Track A0 error cut" );
+ declareProperty("ZTrkErrorCut", m_zTrkErrorCut, "Track Z impact error cut" );
+ declareProperty("VrtMassLimit", m_VrtMassLimit, "Maximal allowed mass for found vertices" );
+ declareProperty("Vrt2TrMassLimit",m_Vrt2TrMassLimit,"Maximal allowed mass for 2-track vertices" );
+ declareProperty("Vrt2TrPtLimit", m_Vrt2TrPtLimit, "Maximal allowed Pt for 2-track vertices. Calibration limit" );
+
+ declareProperty("Sel2VrtProbCut", m_sel2VrtProbCut, "Cut on probability of 2-track vertex for initial selection" );
+ declareProperty("GlobVrtProbCut", m_globVrtProbCut, "Cut on probability of any vertex for final selection" );
+ declareProperty("MaxSVRadiusCut", m_maxSVRadiusCut, "Cut on maximal radius of SV (def = Pixel detector size)" );
+ declareProperty("SelVrtSigCut", m_selVrtSigCut, "Cut on significance of 3D distance between vertex and PV" );
+ declareProperty("AntiPileupSigRCut", m_antiPileupSigRCut, "Upper cut on significance of 2D distance between beam and perigee" );
+ declareProperty("dRdZRatioCut", m_dRdZRatioCut, "Cut on dR/dZ ratio to remove pileup tracks" );
+ declareProperty("v2tIniBDTCut", m_v2tIniBDTCut, "Initial BDT cut for 2track vertices selection " );
+ declareProperty("v2tFinBDTCut", m_v2tFinBDTCut, "Final BDT cut for 2track vertices selection " );
+
+ declareProperty("FillHist", m_fillHist, "Fill technical histograms" );
+
+
+ declareProperty("useVertexCleaning", m_useVertexCleaning, "Clean vertices by requiring pixel hit presence according to vertex position" );
+
+ declareProperty("MultiWithOneTrkVrt", m_multiWithOneTrkVrt,"Allow one-track-vertex addition to already found secondary vertices");
+
+ declareProperty("VertexMergeCut", m_vertexMergeCut, "To allow vertex merging for MultiVertex Finder" );
+ declareProperty("TrackDetachCut", m_trackDetachCut, "To allow track from vertex detachment for MultiVertex Finder" );
+
+ declareProperty("beampipeR", m_beampipeR, "Radius of the beampipe material" );
+ declareProperty("removeTrkMatSignif", m_removeTrkMatSignif, "Significance of Vertex-TrackingMaterial distance for removal. No removal if <=0." );
+
+ declareProperty("calibFileName", m_calibFileName, " MVA calibration file for 2-track fake vertices removal" );
+
+ declareProperty("BeamSpotSvc", m_beamService, "Name of the BeamSpot service");
+ declareProperty("VertexFitterTool", m_fitSvc, "Name of the Vertex Fitter tool");
+//
+ m_massPi = Trk::ParticleMasses().mass[Trk::pion];
+ m_massP = Trk::ParticleMasses().mass[Trk::proton];
+ m_massE = Trk::ParticleMasses().mass[Trk::electron];
+ m_massK0 = Trk::ParticleMasses().mass[Trk::k0];
+ m_massLam = 1115.683 ;
+ m_compatibilityGraph = nullptr;
+ m_instanceName=name;
+
+ }
+
+//Destructor---------------------------------------------------------------
+ NewVrtSecInclusiveTool::~NewVrtSecInclusiveTool(){
+ ATH_MSG_DEBUG("NewVrtSecInclusiveTool destructor called");
+ }
+
+//Initialize---------------------------------------------------------------
+ StatusCode NewVrtSecInclusiveTool::initialize(){
+ ATH_MSG_DEBUG( "Initialising NewVrtSecInclusiveTool- Package version: " << PACKAGE_VERSION );
+ m_compatibilityGraph = new boost::adjacency_list<boost::listS, boost::vecS, boost::undirectedS>();
+
+ ATH_CHECK( m_beamService.retrieve() );
+ ATH_CHECK( m_extrapolator.retrieve() );
+
+ ATH_CHECK( m_fitSvc.retrieve() );
+ ATH_MSG_DEBUG("NewVrtSecInclusiveTool TrkVKalVrtFitter found");
+
+//------------------------------------------
+//
+ ITHistSvc* hist_root=0;
+ if(m_fillHist){
+
+ StatusCode sc = service( "THistSvc", hist_root);
+ if( sc.isFailure() ) ATH_MSG_DEBUG("Could not find THistSvc service");
+ else ATH_MSG_DEBUG("NewVrtSecInclusiveTool Histograms found");
+
+ m_hb_massPiPi = new TH1D("massPiPi"," mass PiPi",200,0., 4000.);
+ m_hb_massPiPi1 = new TH1D("massPiPi1"," mass PiPi",200,0., 4000.);
+ m_hb_massPPi = new TH1D("massPPi"," massPPi", 100,1000., 1250.);
+ m_hb_massEE = new TH1D("massEE"," massEE", 100,0., 200.);
+ m_hb_nvrt2 = new TH1D("nvrt2"," vertices2", 50,0., 50.);
+ m_hb_ratio = new TH1D("ratio"," ratio", 51,0., 1.02);
+ m_hb_totmass = new TH1D("totmass"," totmass", 250,0., 10000.);
+ m_hb_impact = new TH1D("impact", " impact", 100,0., 20.);
+ m_hb_impactR = new TH1D("impactR"," impactR", 400,-30., 70.);
+ m_hb_impactZ = new TH1D("impactZ"," impactZ", 100,-30., 70.);
+ m_hb_impactRZ = new TH2D("impactRZ"," impactRZ", 40,-10., 10., 60, -30.,30. );
+ m_hb_pileupRat = new TH1D("pileupRatio"," dR/dZ ratio", 100, 0., 1.);
+ m_hb_trkD0 = new TH1D("trkD0"," d0 of tracks", 100, 0., 10.);
+ m_hb_trkZ = new TH1D("trkZ"," Z of tracks", 120,-60., 60.);
+ m_hb_r2d = new TH1D("r2interact","Interaction radius 2tr selected", 150,0., 150.);
+ m_hb_ntrksel = new TH1F("NTrkSel","Number of selected tracks", 200,0., 200.);
+ m_hb_ntrkInput = new TH1F("NTrkInput","Number of provided tracks", 200,0., 1000.);
+ m_hb_trkSelect = new TH1F("TrkSelect","Track selection efficiency", 15,0., 15.);
+ m_hb_signif3D = new TH1D("signif3D"," Signif3D for initial 2tr vertex", 140,-20., 50.);
+ m_hb_sig3DTot = new TH1D("sig3dcommon"," Signif3D for common vertex", 140,-20., 50.);
+ m_hb_sig3D1tr = new TH1D("sig3D1tr","Signif3D for 1tr vertices", 140,-20., 50.);
+ m_hb_sig3D2tr = new TH1D("sig3D2tr","Signif3D for 2tr single vertex", 140,-20., 50.);
+ m_hb_sig3DNtr = new TH1D("sig3DNtr","Signif3D for many-tr single vertex", 140,-20., 50.);
+ m_hb_goodvrtN = new TH1F("goodvrtN","Number of good vertices", 20,0., 20.);
+ m_hb_goodvrt1N = new TH1F("goodvrt1N","Number of good 1-track vertices", 20,0., 20.);
+ m_hb_distVV = new TH1D("distvv","Vertex-Vertex dist", 100,0., 20.);
+ m_hb_diffPS = new TH1D("diffPS","Primary-Secondary assoc", 200,-20., 20.);
+ m_hb_rawVrtN = new TH1F("rawVrtN","Number of raw vertices multivertex case", 20, 0., 20.);
+ m_hb_cosSVMom = new TH1F("cosSVMom","SV-PV vs SV momentum ", 100, 0., 1.);
+ m_hb_etaSV = new TH1F("etaSV"," Eta of SV-PV ", 100, -5., 5.);
+ m_hb_fakeSVBDT = new TH1F("fakeSVBDT"," BDT for fake SV rejection", 100, -1., 1.);
+//---
+ std::string histDir;
+ histDir="/file1/stat/MultiSVrt"+m_instanceName+"/";
+ sc = hist_root->regHist(histDir+"massPiPi", m_hb_massPiPi);
+ sc = hist_root->regHist(histDir+"massPiPi1", m_hb_massPiPi1);
+ sc = hist_root->regHist(histDir+"massPPi", m_hb_massPPi);
+ sc = hist_root->regHist(histDir+"massEE", m_hb_massEE );
+ sc = hist_root->regHist(histDir+"nvrt2", m_hb_nvrt2);
+ sc = hist_root->regHist(histDir+"ratio", m_hb_ratio);
+ sc = hist_root->regHist(histDir+"totmass", m_hb_totmass);
+ sc = hist_root->regHist(histDir+"impact", m_hb_impact);
+ sc = hist_root->regHist(histDir+"impactR", m_hb_impactR);
+ sc = hist_root->regHist(histDir+"impactZ", m_hb_impactZ);
+ sc = hist_root->regHist(histDir+"impactRZ", m_hb_impactRZ);
+ sc = hist_root->regHist(histDir+"pileupRatio", m_hb_pileupRat);
+ sc = hist_root->regHist(histDir+"trkD0", m_hb_trkD0);
+ sc = hist_root->regHist(histDir+"trkZ", m_hb_trkZ);
+ sc = hist_root->regHist(histDir+"r2interact",m_hb_r2d);
+ sc = hist_root->regHist(histDir+"NTrkSel", m_hb_ntrksel);
+ sc = hist_root->regHist(histDir+"NTrkInput", m_hb_ntrkInput);
+ sc = hist_root->regHist(histDir+"TrkSelect", m_hb_trkSelect);
+ sc = hist_root->regHist(histDir+"signif3D", m_hb_signif3D);
+ sc = hist_root->regHist(histDir+"sig3dcommon", m_hb_sig3DTot);
+ sc = hist_root->regHist(histDir+"sig3D1tr", m_hb_sig3D1tr);
+ sc = hist_root->regHist(histDir+"sig3D2tr", m_hb_sig3D2tr);
+ sc = hist_root->regHist(histDir+"sig3DNtr", m_hb_sig3DNtr);
+ sc = hist_root->regHist(histDir+"goodvrtN", m_hb_goodvrtN);
+ sc = hist_root->regHist(histDir+"goodvrt1N", m_hb_goodvrt1N);
+ sc = hist_root->regHist(histDir+"distVV", m_hb_distVV);
+ sc = hist_root->regHist(histDir+"diffPS", m_hb_diffPS);
+ sc = hist_root->regHist(histDir+"rawVrtN", m_hb_rawVrtN);
+ sc = hist_root->regHist(histDir+"cosSVMom", m_hb_cosSVMom);
+ sc = hist_root->regHist(histDir+"etaSV", m_hb_etaSV);
+ sc = hist_root->regHist(histDir+"fakeSVBDT", m_hb_fakeSVBDT);
+ if( sc.isFailure() ) ATH_MSG_INFO( "NewVrtSecInclusive Histogram registration failure!!!");
+ m_w_1 = 1.;
+ //------
+ m_curTup=new DevTuple();
+ m_tuple = new TTree("Vertices","Vertices");
+ std::string TreeDir;
+ TreeDir="/file1/stat/MultiSVrt"+m_instanceName+"/";
+ sc = hist_root->regTree(TreeDir,m_tuple);
+ if (sc.isSuccess()) {
+ m_tuple->Branch("ntrk", &m_curTup->nTrk, "ntrk/I");
+ m_tuple->Branch("pttrk", &m_curTup->pttrk, "pttrk[ntrk]/F");
+ m_tuple->Branch("Sig3D", &m_curTup->Sig3D, "Sig3D[ntrk]/F");
+
+ m_tuple->Branch("n2Vrt", &m_curTup->n2Vrt, "n2Vrt/I");
+ m_tuple->Branch("VrtTrkHF", &m_curTup->VrtTrkHF, "VrtTrkHF[n2Vrt]/I");
+ m_tuple->Branch("VrtTrkI", &m_curTup->VrtTrkI, "VrtTrkI[n2Vrt]/I");
+ m_tuple->Branch("VrtCh", &m_curTup->VrtCh, "VrtCh[n2Vrt]/I");
+ m_tuple->Branch("VrtDist2D", &m_curTup->VrtDist2D, "VrtDist2D[n2Vrt]/F");
+ m_tuple->Branch("VrtSig3D", &m_curTup->VrtSig3D, "VrtSig3D[n2Vrt]/F");
+ m_tuple->Branch("VrtSig2D", &m_curTup->VrtSig2D, "VrtSig2D[n2Vrt]/F");
+ m_tuple->Branch("VrtM", &m_curTup->VrtM, "VrtM[n2Vrt]/F");
+ m_tuple->Branch("VrtZ", &m_curTup->VrtZ, "VrtZ[n2Vrt]/F");
+ m_tuple->Branch("VrtPt", &m_curTup->VrtPt, "VrtPt[n2Vrt]/F");
+ m_tuple->Branch("VrtEta", &m_curTup->VrtEta, "VrtEta[n2Vrt]/F");
+ m_tuple->Branch("VrtIBL", &m_curTup->VrtIBL, "VrtIBL[n2Vrt]/I");
+ m_tuple->Branch("VrtBL", &m_curTup->VrtBL, "VrtBL[n2Vrt]/I");
+ m_tuple->Branch("VrtSinSPM", &m_curTup->VrtSinSPM, "VrtSinSPM[n2Vrt]/F");
+ m_tuple->Branch("VMinPtT", &m_curTup->VMinPtT, "VMinPtT[n2Vrt]/F");
+ m_tuple->Branch("VMinS3DT", &m_curTup->VMinS3DT, "VMinS3DT[n2Vrt]/F");
+ m_tuple->Branch("VMaxS3DT", &m_curTup->VMaxS3DT, "VMaxS3DT[n2Vrt]/F");
+ m_tuple->Branch("VrtProb", &m_curTup->VrtProb, "VrtProb[n2Vrt]/F");
+ m_tuple->Branch("VrtHR1", &m_curTup->VrtHR1, "VrtHR1[n2Vrt]/F");
+ m_tuple->Branch("VrtHR2", &m_curTup->VrtHR2, "VrtHR2[n2Vrt]/F");
+ m_tuple->Branch("VrtBDT", &m_curTup->VrtBDT, "VrtBDT[n2Vrt]/F");
+ m_tuple->Branch("VrtDisk", &m_curTup->VrtDisk, "VrtDisk[n2Vrt]/I");
+ m_tuple->Branch("VSigMat", &m_curTup->VSigMat, "VSigMat[n2Vrt]/F");
+
+ m_tuple->Branch("nNVrt", &m_curTup->nNVrt, "nNVrt/I");
+ m_tuple->Branch("NVrtTrk", &m_curTup->NVrtTrk, "NVrtTrk[nNVrt]/I");
+ m_tuple->Branch("NVrtTrkHF", &m_curTup->NVrtTrkHF, "NVrtTrkHF[nNVrt]/I");
+ m_tuple->Branch("NVrtTrkI", &m_curTup->NVrtTrkI, "NVrtTrkI[nNVrt]/I");
+ m_tuple->Branch("NVrtCh", &m_curTup->NVrtCh, "NVrtCh[nNVrt]/I");
+ m_tuple->Branch("NVrtDist2D", &m_curTup->NVrtDist2D, "NVrtDist2D[nNVrt]/F");
+ m_tuple->Branch("NVrtSig3D", &m_curTup->NVrtSig3D, "NVrtSig3D[nNVrt]/F");
+ m_tuple->Branch("NVrtSig2D", &m_curTup->NVrtSig2D, "NVrtSig2D[nNVrt]/F");
+ m_tuple->Branch("NVrtM", &m_curTup->NVrtM, "NVrtM[nNVrt]/F");
+ m_tuple->Branch("NVrtPt", &m_curTup->NVrtPt, "NVrtPt[nNVrt]/F");
+ m_tuple->Branch("NVrtEta", &m_curTup->NVrtEta, "NVrtEta[nNVrt]/F");
+ m_tuple->Branch("NVrtIBL", &m_curTup->NVrtIBL, "NVrtIBL[nNVrt]/I");
+ m_tuple->Branch("NVrtBL", &m_curTup->NVrtBL, "NVrtBL[nNVrt]/I");
+ m_tuple->Branch("NVrtSinSPM", &m_curTup->NVrtSinSPM, "NVrtSinSPM[nNVrt]/F");
+ m_tuple->Branch("NVMinPtT", &m_curTup->NVMinPtT, "NVMinPtT[nNVrt]/F");
+ m_tuple->Branch("NVMinS3DT", &m_curTup->NVMinS3DT, "NVMinS3DT[nNVrt]/F");
+ m_tuple->Branch("NVrtProb", &m_curTup->NVrtProb, "NVrtProb[nNVrt]/F");
+ m_tuple->Branch("NVrtBDT", &m_curTup->NVrtBDT, "NVrtBDT[nNVrt]/F");
+ }
+ }
+
+//--------------------------------------------------------
+ //std::string fileName="NewVrtSecInclusiveTool/Fake2TrVertexReject.MVA.v01.root"; ///For local calibration file
+ //std::string rootFilePath = PathResolver::find_file(fileName, "DATAPATH"); ///
+ std::string rootFilePath = PathResolver::find_calib_file("NewVrtSecInclusiveTool/"+m_calibFileName);
+ TFile* rootFile = TFile::Open(rootFilePath.c_str(), "READ");
+ if (!rootFile) {
+ ATH_MSG_FATAL("Could not retrieve root file: " << m_calibFileName);
+ return StatusCode::FAILURE;
+ }
+ TTree * training = (TTree*)rootFile->Get("BDT");
+ m_SV2T_BDT = std::make_unique<MVAUtils::BDT>(training);
+//--------------------------------------------------------
+ return StatusCode::SUCCESS;
+
+ }
+
+
+
+
+ StatusCode NewVrtSecInclusiveTool::finalize()
+ {
+ if(m_compatibilityGraph)delete m_compatibilityGraph;
+ ATH_MSG_DEBUG("NewVrtSecInclusiveTool finalize()");
+ return StatusCode::SUCCESS;
+ }
+
+
+
+
+ std::unique_ptr<Trk::VxSecVertexInfo> NewVrtSecInclusiveTool::findAllVertices (
+ const std::vector<const xAOD::TrackParticle*> & inpTrk,
+ const xAOD::Vertex & primVrt ) const
+ {
+ std::vector<xAOD::Vertex*> listVrtSec(0);
+
+ if(m_fillHist && m_curTup){
+ m_curTup->nTrk=0;
+ m_curTup->n2Vrt=0;
+ m_curTup->nNVrt=0;
+ };
+
+ workVectorArrxAOD * tmpVectxAOD=new workVectorArrxAOD();
+ tmpVectxAOD->inpTrk.resize(inpTrk.size());
+ std::copy(inpTrk.begin(),inpTrk.end(), tmpVectxAOD->inpTrk.begin());
+ tmpVectxAOD->beamX=m_beamService->beamPos().x();
+ tmpVectxAOD->beamY=m_beamService->beamPos().y();
+ tmpVectxAOD->beamZ=m_beamService->beamPos().z();
+ tmpVectxAOD->tanBeamTiltX=tan(m_beamService->beamTilt(0));
+ tmpVectxAOD->tanBeamTiltY=tan(m_beamService->beamTilt(1));
+
+
+ listVrtSec = getVrtSecMulti(tmpVectxAOD,primVrt);
+ delete tmpVectxAOD;
+
+
+
+ std::vector<const xAOD::IParticle*> iparTrkV0(0);
+ std::unique_ptr<Trk::VxSecVertexInfo> res = std::make_unique<Trk::VxSecVertexInfo>(Trk::VxSecVertexInfo(listVrtSec));
+
+ if(m_fillHist && m_tuple){ m_tuple->Fill(); };
+
+ m_compatibilityGraph->clear();
+ return res;
+ }
+
+
+} // end Rec namespace
diff --git a/Reconstruction/VKalVrt/NewVrtSecInclusiveTool/src/Sel2TrkVertices.cxx b/Reconstruction/VKalVrt/NewVrtSecInclusiveTool/src/Sel2TrkVertices.cxx
new file mode 100755
index 0000000000000000000000000000000000000000..1f9e9b53f9a5e904cfc76b0485dc23978d13e497
--- /dev/null
+++ b/Reconstruction/VKalVrt/NewVrtSecInclusiveTool/src/Sel2TrkVertices.cxx
@@ -0,0 +1,215 @@
+/*
+ Copyright (C) 2002-2020 CERN for the benefit of the ATLAS collaboration
+*/
+///
+/// @author Vadim Kostyukhin <vadim.kostyukhin@cern.ch>
+///
+// Header include
+#include "NewVrtSecInclusiveTool/NewVrtSecInclusiveTool.h"
+#include "GeoPrimitives/GeoPrimitivesHelpers.h"
+#include "MVAUtils/BDT.h"
+//-------------------------------------------------
+// Other stuff
+#include "AnalysisUtils/AnalysisMisc.h"
+#include "TrkVKalVrtFitter/TrkVKalVrtFitter.h"
+
+#include "TMath.h"
+#include "TH1.h"
+#include "TH2D.h"
+//
+//#include<iostream>
+
+
+namespace Rec{
+
+
+//
+//
+//--------------------------------------------------------
+// Template routine for 2track secondary vertices selection
+//
+
+ void NewVrtSecInclusiveTool::select2TrVrt(std::vector<const xAOD::TrackParticle*> & selectedTracks,
+ const xAOD::Vertex & PrimVrt)
+ const
+ {
+ std::vector<const xAOD::NeutralParticle*> neutralPartDummy(0);
+ std::vector<const xAOD::TrackParticle*> tracksForFit(2,0);
+ std::vector<double> impact,impactError;
+ std::vector<double> inpMass(2,m_massPi);
+ long int Charge;
+ int i,j;
+ StatusCode sc; sc.setChecked();
+ Vrt2Tr tmpVrt;
+ std::vector<Vrt2Tr> all2TrVrt(0);
+ TLorentzVector PSum2T;
+ Amg::Vector3D iniVrt(0.,0.,0.);
+ //
+ int NTracks = (int) (selectedTracks.size());
+//
+// impact parameters with sign calculations
+//
+ std::vector<float> covPV=PrimVrt.covariance();
+ double signifR=0.,signifZ=0.;
+ std::vector<int> hitIBL(NTracks,0), hitBL(NTracks,0);
+ std::vector<double> trackSignif(NTracks),dRdZratio(NTracks);
+ for (i=0; i<NTracks; i++) {
+ m_fitSvc->VKalGetImpact(selectedTracks[i], PrimVrt.position(), 1, impact, impactError);
+ signifR = impact[0]/ sqrt(impactError[0]);
+ signifZ = impact[1]/ sqrt(impactError[2]);
+ trackSignif[i] = sqrt( signifR*signifR + signifZ*signifZ);
+ dRdZratio[i] = fabs(signifR/signifZ);
+ if(m_fillHist){
+ m_hb_impactR->Fill( signifR, m_w_1);
+ m_hb_impactZ->Fill( signifZ, m_w_1);
+ m_hb_impactRZ->Fill(signifR, signifZ, m_w_1);
+ m_hb_impact->Fill( trackSignif[i], m_w_1);
+ if(trackSignif[i]>2.) m_hb_pileupRat->Fill(dRdZratio[i],1.);
+ if( i<DevTuple::maxNTrk && m_curTup){
+ m_curTup->pttrk[i]=selectedTracks[i]->pt();
+ m_curTup->Sig3D[i]=trackSignif[i];
+ }
+ }
+ }
+
+ if( m_fillHist && m_curTup ){
+ m_curTup->nTrk=std::min(NTracks,DevTuple::maxNTrk);
+ m_curTup->n2Vrt=0;
+ }
+
+ for (i=0; i<NTracks-1; i++) {
+ if(trackSignif[i]<m_trkSigCut || dRdZratio[i]<m_dRdZRatioCut )continue;
+ for (j=i+1; j<NTracks; j++) {
+ if(trackSignif[j]<m_trkSigCut || dRdZratio[j]<m_dRdZRatioCut )continue;
+ PSum2T=selectedTracks[i]->p4()+selectedTracks[j]->p4();
+ if(PSum2T.M()>m_Vrt2TrMassLimit)continue;
+
+ std::unique_ptr<Trk::IVKalState> state = m_fitSvc->makeState();
+ m_fitSvc->setMassInputParticles( inpMass, *state ); // Use pion masses for fit
+ tracksForFit[0]=selectedTracks[i];
+ tracksForFit[1]=selectedTracks[j];
+ sc=m_fitSvc->VKalVrtFitFast(tracksForFit,tmpVrt.fitVertex,*state); /* Fast crude estimation*/
+
+ if( sc.isFailure() ) { /* No initial estimation */
+ iniVrt=PrimVrt.position();
+ } else {
+ double PMomVrtDir = ProjSV_PV(tmpVrt.fitVertex,PrimVrt,PSum2T);
+ if( PMomVrtDir>0. ) iniVrt=tmpVrt.fitVertex; /* Good initial estimation */
+ else iniVrt=PrimVrt.position();
+ }
+ m_fitSvc->setApproximateVertex(iniVrt.x(), iniVrt.y(), iniVrt.z(),*state);
+ sc=m_fitSvc->VKalVrtFit(tracksForFit, neutralPartDummy, tmpVrt.fitVertex, tmpVrt.momentum, Charge,
+ tmpVrt.errorMatrix, tmpVrt.chi2PerTrk, tmpVrt.trkAtVrt, tmpVrt.chi2, *state, true );
+ if(sc.isFailure()) continue; /* No fit */
+ double Prob2v=TMath::Prob(tmpVrt.chi2,1);
+ if( Prob2v < m_sel2VrtProbCut ) continue;
+ if( tmpVrt.momentum.M()> m_Vrt2TrMassLimit ) continue;
+ if( ProjSV_PV(tmpVrt.fitVertex, PrimVrt, tmpVrt.momentum) < 0.) continue; // SV-PV don't coincide with summary momentum direction
+ if( tmpVrt.fitVertex.perp() > m_maxSVRadiusCut) continue; // Too far from interaction point
+ double cosSVPV=ProjSV_PV(tmpVrt.fitVertex, PrimVrt, tmpVrt.momentum);
+ TLorentzVector SVPV(tmpVrt.fitVertex.x()-PrimVrt.x(),
+ tmpVrt.fitVertex.y()-PrimVrt.y(),
+ tmpVrt.fitVertex.z()-PrimVrt.z(), 10.);
+ if(m_fillHist){
+ if(Charge==0){m_hb_massPiPi->Fill(tmpVrt.momentum.M(),1.);}
+ m_hb_cosSVMom->Fill(cosSVPV,1.);
+ m_hb_etaSV->Fill(SVPV.Eta(),1.);
+ }
+ if(cosSVPV<0.8)continue;
+ if(tmpVrt.momentum.Pt()<1000.)continue;
+ float vrtR=tmpVrt.fitVertex.perp();
+//Check close material layer
+ double dstMatSignif=1.e4;
+ if(m_removeTrkMatSignif>0. && vrtR>20.){
+ if(vrtR<30.){ dstMatSignif=fabs(vrtR-m_beampipeR)/VrtRadiusError(tmpVrt.fitVertex,tmpVrt.errorMatrix );} //beampipe
+ else { dstMatSignif=distToMatLayerSignificance(tmpVrt);} //Material in Pixel volume
+ if(dstMatSignif<m_removeTrkMatSignif)continue;
+ }
+//
+// Check pixel hits vs vertex positions.
+ int ihitIBL = getIBLHit(selectedTracks[i]);
+ int jhitIBL = getIBLHit(selectedTracks[j]);
+ if( ihitIBL<0 && jhitIBL<0)continue;
+ float ihitR = selectedTracks[i]->radiusOfFirstHit();
+ float jhitR = selectedTracks[j]->radiusOfFirstHit();
+ if(m_useVertexCleaning){
+ if(fabs(ihitR-jhitR)>15.) continue;
+ if( std::min(ihitR,jhitR)-vrtR > 36.) continue; // Too big dR between vertex and hit in pixel
+ // Should be another layer in between
+ if( std::max(ihitR,jhitR)-vrtR <-10.) continue; // Vertex is behind hit in pixel
+ }
+//
+// Debugging and BDT
+ double minPtT = std::min(tracksForFit[0]->pt(),tracksForFit[1]->pt());
+ double Sig3D=0.,Sig2D=0., Dist2D=0.;
+ if( m_fillHist && m_curTup ){
+ int ihitBL = getBLHit (selectedTracks[i]);
+ int jhitBL = getBLHit (selectedTracks[j]);
+ int idisk1=0,idisk2=0,idisk3=0,jdisk1=0,jdisk2=0,jdisk3=0;
+ int sumIBLHits = std::max(ihitIBL,0)+std::max(jhitIBL,0);
+ int sumBLHits = std::max(ihitBL, 0)+std::max(jhitBL, 0);
+ getPixelDiscs(selectedTracks[i],idisk1,idisk2,idisk3);
+ getPixelDiscs(selectedTracks[j],jdisk1,jdisk2,jdisk3);
+ VrtVrtDist(PrimVrt, tmpVrt.fitVertex, tmpVrt.errorMatrix, Sig3D);
+ Dist2D=VrtVrtDist2D(PrimVrt, tmpVrt.fitVertex, tmpVrt.errorMatrix, Sig2D);
+ m_hb_signif3D->Fill(Sig3D,1.);
+ m_curTup->VrtTrkHF [m_curTup->n2Vrt] = getIdHF(tracksForFit[0])+ getIdHF(tracksForFit[1]);
+ m_curTup->VrtTrkI [m_curTup->n2Vrt] = getG4Inter(tracksForFit[0])+ getG4Inter(tracksForFit[1]);
+ m_curTup->VrtCh [m_curTup->n2Vrt] = Charge;
+ m_curTup->VrtProb [m_curTup->n2Vrt] = Prob2v;
+ m_curTup->VrtSig3D [m_curTup->n2Vrt] = Sig3D;
+ m_curTup->VrtSig2D [m_curTup->n2Vrt] = Sig2D;
+ m_curTup->VrtDist2D[m_curTup->n2Vrt] = vrtR<20. ? Dist2D : vrtR;
+ m_curTup->VrtM [m_curTup->n2Vrt] = tmpVrt.momentum.M();
+ m_curTup->VrtZ [m_curTup->n2Vrt] = tmpVrt.fitVertex.z();
+ m_curTup->VrtPt [m_curTup->n2Vrt] = tmpVrt.momentum.Pt();
+ m_curTup->VrtEta [m_curTup->n2Vrt] = SVPV.Eta();
+ m_curTup->VrtIBL [m_curTup->n2Vrt] = sumIBLHits;
+ m_curTup->VrtBL [m_curTup->n2Vrt] = sumBLHits;
+ m_curTup->VrtSinSPM[m_curTup->n2Vrt] = sqrt(1.-cosSVPV*cosSVPV);
+ m_curTup->VMinPtT [m_curTup->n2Vrt] = minPtT;
+ m_curTup->VMinS3DT [m_curTup->n2Vrt] = std::min(trackSignif[i],trackSignif[j]);
+ m_curTup->VMaxS3DT [m_curTup->n2Vrt] = std::max(trackSignif[i],trackSignif[j]);
+ m_curTup->VrtBDT [m_curTup->n2Vrt] = 1.1;
+ m_curTup->VrtHR1 [m_curTup->n2Vrt] = ihitR;
+ m_curTup->VrtHR2 [m_curTup->n2Vrt] = jhitR;
+ m_curTup->VrtDisk [m_curTup->n2Vrt] = idisk1+10*idisk2+20*idisk3+30*jdisk1+40*jdisk2+50*jdisk3;
+ m_curTup->VSigMat [m_curTup->n2Vrt] = dstMatSignif;
+ if(m_curTup->n2Vrt<DevTuple::maxNVrt-1)m_curTup->n2Vrt++;
+ }
+//-------------------BDT based rejection
+ if(tmpVrt.momentum.Pt() > m_Vrt2TrPtLimit) continue;
+ VrtVrtDist(PrimVrt, tmpVrt.fitVertex, tmpVrt.errorMatrix, Sig3D);
+ Dist2D=VrtVrtDist2D(PrimVrt, tmpVrt.fitVertex, tmpVrt.errorMatrix, Sig2D);
+ std::vector<float> VARS(10);
+ VARS[0]=Prob2v;
+ VARS[1]=log(tmpVrt.momentum.Pt());
+ VARS[2]=log(std::max(minPtT,m_cutPt));
+ VARS[3]=log(vrtR<20. ? Dist2D : vrtR);
+ VARS[4]=log(std::max(std::min(trackSignif[i],trackSignif[j]),m_trkSigCut));
+ VARS[5]=log(std::max(trackSignif[i],trackSignif[j]));
+ VARS[6]=tmpVrt.momentum.M();
+ VARS[7]=sqrt(fabs(1.-cosSVPV*cosSVPV));
+ VARS[8]=SVPV.Eta();
+ VARS[9]=std::max(ihitR,jhitR);
+ //VARS[9]=sumIBLHits;
+ //VARS[10]=sumBLHits;
+ //VARS[11]=std::max(Sig3D,m_selVrtSigCut);
+ float wgtSelect=m_SV2T_BDT->GetGradBoostMVA(VARS);
+ //std::vector<float> weights=m_SV2T_BDT->GetMultiResponse(VARS,3);
+ //float wgtSelect=weights[0];
+ if( m_fillHist && m_curTup ) m_curTup->VrtBDT[m_curTup->n2Vrt-1] = wgtSelect;
+ if(wgtSelect<m_v2tIniBDTCut) continue;
+ //if(weights[2]>weights[0] && weights[2]>weights[1])continue;
+//
+//--- Save good candidate for multi-vertex fit
+//
+ add_edge(i,j,*m_compatibilityGraph);
+ }
+ }
+
+ return;
+ }
+
+
+} //end of namespace
diff --git a/Reconstruction/VKalVrt/NewVrtSecInclusiveTool/src/Utilities.cxx b/Reconstruction/VKalVrt/NewVrtSecInclusiveTool/src/Utilities.cxx
new file mode 100755
index 0000000000000000000000000000000000000000..9377d1a92fb868ac24e8a363d891d905cd168be8
--- /dev/null
+++ b/Reconstruction/VKalVrt/NewVrtSecInclusiveTool/src/Utilities.cxx
@@ -0,0 +1,388 @@
+/*
+ Copyright (C) 2002-2019 CERN for the benefit of the ATLAS collaboration
+*/
+///
+/// @author Vadim Kostyukhin <vadim.kostyukhin@cern.ch>
+///
+// Header include
+#include "NewVrtSecInclusiveTool/NewVrtSecInclusiveTool.h"
+#include "TrkNeutralParameters/NeutralParameters.h"
+#include "TrkTrackSummary/TrackSummary.h"
+#include "TrkVKalVrtFitter/TrkVKalVrtFitter.h"
+//-------------------------------------------------
+#include "TrkGeometry/TrackingGeometry.h"
+#include "TrkGeometry/TrackingVolume.h"
+#include "TrkGeometry/Layer.h"
+ // Other stuff
+#include <cmath>
+
+
+namespace Rec{
+
+
+// void NewVrtSecInclusiveTool::printWrkSet(const std::vector<WrkVrt> *, const & std::string ) const {
+ void NewVrtSecInclusiveTool::printWrkSet(const std::vector<WrkVrt> *WrkVrtSet, const std::string &name) const {
+ int nGoodV=0;
+ for(int iv=0; iv<(int)WrkVrtSet->size(); iv++) {
+ std::cout<<name
+ <<"= "<<(*WrkVrtSet)[iv].vertex[0]
+ <<", "<<(*WrkVrtSet)[iv].vertex[1]
+ <<", "<<(*WrkVrtSet)[iv].vertex[2]
+ <<" NTrk="<<(*WrkVrtSet)[iv].selTrk.size()
+ <<" is good="<<std::boolalpha<<(*WrkVrtSet)[iv].Good<<std::noboolalpha
+ <<" Chi2="<<(*WrkVrtSet)[iv].chi2
+ <<" Mass="<<(*WrkVrtSet)[iv].vertexMom.M()
+ <<" detached="<<(*WrkVrtSet)[iv].detachedTrack
+ <<" proj.dist="<<(*WrkVrtSet)[iv].projectedVrt
+ <<" trk=";
+ for(int kk=0; kk<(int)(*WrkVrtSet)[iv].selTrk.size(); kk++) {
+ std::cout<<", "<<(*WrkVrtSet)[iv].selTrk[kk];}
+ //for(int kk=0; kk<(int)(*WrkVrtSet)[iv].selTrk.size(); kk++) {
+ // std::cout<<", "<<MomAtVrt((*WrkVrtSet)[iv].trkAtVrt[kk]).Perp();}
+ std::cout<<'\n';
+ if((*WrkVrtSet)[iv].Good)nGoodV++;
+ }
+ std::cout<<name<<" N="<<nGoodV<<'\n';
+ }
+
+ /* Technicalities */
+ double NewVrtSecInclusiveTool::ProjSV_PV(const Amg::Vector3D & SV, const xAOD::Vertex & PV, const TLorentzVector & Direction) const
+ {
+ TVector3 SV_PV( SV.x()-PV.x(), SV.y()-PV.y(), SV.z()-PV.z() );
+ return Direction.Vect().Unit()*SV_PV.Unit();
+ }
+
+
+ double NewVrtSecInclusiveTool::VrtVrtDist(const xAOD::Vertex & PrimVrt, const Amg::Vector3D & SecVrt,
+ const std::vector<double> SecVrtErr, double& signif)
+ const
+ {
+ double distx = PrimVrt.x()- SecVrt.x();
+ double disty = PrimVrt.y()- SecVrt.y();
+ double distz = PrimVrt.z()- SecVrt.z();
+
+
+ AmgSymMatrix(3) PrimCovMtx=PrimVrt.covariancePosition(); //Create
+ PrimCovMtx(0,0) += SecVrtErr[0];
+ PrimCovMtx(0,1) += SecVrtErr[1];
+ PrimCovMtx(1,0) += SecVrtErr[1];
+ PrimCovMtx(1,1) += SecVrtErr[2];
+ PrimCovMtx(0,2) += SecVrtErr[3];
+ PrimCovMtx(2,0) += SecVrtErr[3];
+ PrimCovMtx(1,2) += SecVrtErr[4];
+ PrimCovMtx(2,1) += SecVrtErr[4];
+ PrimCovMtx(2,2) += SecVrtErr[5];
+
+ AmgSymMatrix(3) WgtMtx = PrimCovMtx.inverse();
+
+ 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=sqrt(fabs(signif));
+ if( signif!=signif ) signif = 0.;
+ return sqrt(distx*distx+disty*disty+distz*distz);
+ }
+
+ double NewVrtSecInclusiveTool::VrtVrtDist2D(const xAOD::Vertex & PrimVrt, const Amg::Vector3D & SecVrt,
+ const std::vector<double> SecVrtErr, double& signif)
+ const
+ {
+ double distx = PrimVrt.x()- SecVrt.x();
+ double disty = PrimVrt.y()- SecVrt.y();
+
+
+ AmgSymMatrix(3) PrimCovMtx=PrimVrt.covariancePosition(); //Create
+ AmgSymMatrix(2) CovMtx;
+ CovMtx(0,0) = PrimCovMtx(0,0) + SecVrtErr[0];
+ CovMtx(0,1) = PrimCovMtx(0,1) + SecVrtErr[1];
+ CovMtx(1,0) = PrimCovMtx(1,0) + SecVrtErr[1];
+ CovMtx(1,1) = PrimCovMtx(1,1) + SecVrtErr[2];
+
+ AmgSymMatrix(2) WgtMtx = CovMtx.inverse();
+
+ signif = distx*WgtMtx(0,0)*distx
+ +disty*WgtMtx(1,1)*disty
+ +2.*distx*WgtMtx(0,1)*disty;
+ signif=sqrt(fabs(signif));
+ if( signif!=signif ) signif = 0.;
+ return sqrt(distx*distx+disty*disty);
+ }
+
+
+ double NewVrtSecInclusiveTool::VrtVrtDist(const Amg::Vector3D & Vrt1, const std::vector<double> & VrtErr1,
+ const Amg::Vector3D & Vrt2, const std::vector<double> & VrtErr2)
+ const
+ {
+ double distx = Vrt1.x()- Vrt2.x();
+ double disty = Vrt1.y()- Vrt2.y();
+ double distz = Vrt1.z()- Vrt2.z();
+
+ AmgSymMatrix(3) PrimCovMtx; //Create
+ PrimCovMtx(0,0) = VrtErr1[0]+VrtErr2[0];
+ PrimCovMtx(0,1) = PrimCovMtx(1,0) = VrtErr1[1]+VrtErr2[1];
+ PrimCovMtx(1,1) = VrtErr1[2]+VrtErr2[2];
+ PrimCovMtx(0,2) = PrimCovMtx(2,0) = VrtErr1[3]+VrtErr2[3];
+ PrimCovMtx(1,2) = PrimCovMtx(2,1) = VrtErr1[4]+VrtErr2[4];
+ PrimCovMtx(2,2) = VrtErr1[5]+VrtErr2[5];
+
+ AmgSymMatrix(3) WgtMtx = PrimCovMtx.inverse();
+
+
+ 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=sqrt(fabs(signif));
+ if(signif != signif) signif = 0.;
+ return signif;
+ }
+//
+ double NewVrtSecInclusiveTool::PntPntDist(const Amg::Vector3D & Vrt1, const Amg::Vector3D & Vrt2) const
+ {
+ double dx = Vrt1.x()- Vrt2.x();
+ double dy = Vrt1.y()- Vrt2.y();
+ double dz = Vrt1.z()- Vrt2.z();
+ return sqrt(dx*dx+dy*dy*dz*dz);
+ }
+//--------------------------------------------------
+// significance along some direction
+//--------------------------------------------------
+double NewVrtSecInclusiveTool::VrtVrtDist(const xAOD::Vertex & PrimVrt, const Amg::Vector3D & SecVrt,
+ const std::vector<double> SecVrtErr, const TLorentzVector & Dir)
+ const
+ {
+ Amg::Vector3D dir(Dir.Vect().Unit().X(), Dir.Vect().Unit().Y(), Dir.Vect().Unit().Z());
+ double projDist=(SecVrt-PrimVrt.position()).dot(dir);
+ double distx = dir.x()*projDist;
+ double disty = dir.y()*projDist;
+ double distz = dir.z()*projDist;
+
+ AmgSymMatrix(3) PrimCovMtx=PrimVrt.covariancePosition(); //Create
+ PrimCovMtx(0,0) += SecVrtErr[0];
+ PrimCovMtx(0,1) += SecVrtErr[1];
+ PrimCovMtx(1,0) += SecVrtErr[1];
+ PrimCovMtx(1,1) += SecVrtErr[2];
+ PrimCovMtx(0,2) += SecVrtErr[3];
+ PrimCovMtx(2,0) += SecVrtErr[3];
+ PrimCovMtx(1,2) += SecVrtErr[4];
+ PrimCovMtx(2,1) += SecVrtErr[4];
+ PrimCovMtx(2,2) += SecVrtErr[5];
+
+ AmgSymMatrix(3) WgtMtx = PrimCovMtx.inverse();
+
+ 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=sqrt(fabs(signif));
+ if( signif!=signif ) signif = 0.;
+ if(projDist<0)signif=-signif;
+ return signif;
+ }
+
+//----------------------------
+// Vertex error along radius
+//----------------------------
+ double NewVrtSecInclusiveTool::VrtRadiusError(const Amg::Vector3D & SecVrt, const std::vector<double> & VrtErr) const
+ {
+ double DirX=SecVrt.x(), DirY=SecVrt.y();
+ double Covar = DirX*VrtErr[0]*DirX
+ +2.*DirX*VrtErr[1]*DirY
+ +DirY*VrtErr[2]*DirY;
+ Covar /= DirX*DirX + DirY*DirY;
+ Covar=sqrt(fabs(Covar));
+ if(Covar != Covar) Covar = 0.;
+ return Covar;
+ }
+
+
+ /* Invariant mass calculation for V0 decays*/
+ /* Gives correct mass assignment in case of nonequal masses*/
+
+
+ double NewVrtSecInclusiveTool::massV0(const std::vector< std::vector<double> >& TrkAtVrt,
+ double massP, double massPi )
+ const
+ {
+ double ap1i=fabs(TrkAtVrt[0][2]); double ap2i=fabs(TrkAtVrt[1][2]);
+ double px = cos(TrkAtVrt[0][0])*sin(TrkAtVrt[0][1])*ap1i
+ + cos(TrkAtVrt[1][0])*sin(TrkAtVrt[1][1])*ap2i;
+ double py = sin(TrkAtVrt[0][0])*sin(TrkAtVrt[0][1])*ap1i
+ + sin(TrkAtVrt[1][0])*sin(TrkAtVrt[1][1])*ap2i;
+ double pz = cos(TrkAtVrt[0][1])*ap1i
+ + cos(TrkAtVrt[1][1])*ap2i;
+ double ee= (ap1i > ap2i) ?
+ (sqrt(ap1i*ap1i+massP*massP)+sqrt(ap2i*ap2i+massPi*massPi)):
+ (sqrt(ap2i*ap2i+massP*massP)+sqrt(ap1i*ap1i+massPi*massPi));
+ double test=(ee-pz)*(ee+pz)-px*px-py*py;
+ return test>0 ? sqrt(test) : 0.;
+ }
+
+
+
+ TLorentzVector NewVrtSecInclusiveTool::MomAtVrt(const std::vector< double >& inpTrk)
+ const
+ {
+ double api=1./fabs(inpTrk[2]);
+ double px = cos ( inpTrk[0]) * sin(inpTrk[1])*api;
+ double py = sin ( inpTrk[0]) * sin(inpTrk[1])*api;
+ double pz = cos(inpTrk[1])*api;
+ double ee = sqrt( px*px + py*py + pz*pz + m_massPi*m_massPi);
+ return TLorentzVector(px,py,pz,ee);
+ }
+
+
+/*************************************************************************************************************/
+ int NewVrtSecInclusiveTool::getIBLHit(const xAOD::TrackParticle* Part) const
+ {
+ uint8_t IBLhit,IBLexp;
+ if(!Part->summaryValue( IBLexp, xAOD::expectInnermostPixelLayerHit) ) IBLexp = 0;
+ if( IBLexp==0 ) return -1;
+ if(!Part->summaryValue( IBLhit, xAOD::numberOfInnermostPixelLayerHits) ) IBLhit = 0;
+ if(IBLhit) return 1;
+ else return 0;
+ }
+ int NewVrtSecInclusiveTool::getBLHit(const xAOD::TrackParticle* Part) const
+ {
+ uint8_t BLhit,BLexp;
+ if(!Part->summaryValue( BLexp, xAOD::expectNextToInnermostPixelLayerHit) ) BLexp = 0;
+ if( BLexp==0 ) return -1;
+ if(!Part->summaryValue( BLhit, xAOD::numberOfNextToInnermostPixelLayerHits) ) BLhit = 0;
+ if(BLhit) return 1;
+ else return 0;
+ }
+
+ void NewVrtSecInclusiveTool::getPixelDiscs(const xAOD::TrackParticle* Part, int &d0Hit, int &d1Hit, int &d2Hit) const
+ {
+ uint32_t HitPattern=Part->hitPattern();
+ d0Hit=0; if( HitPattern&((1<<Trk::pixelEndCap0)) ) d0Hit=1;
+ d1Hit=0; if( HitPattern&((1<<Trk::pixelEndCap1)) ) d1Hit=1;
+ d2Hit=0; if( HitPattern&((1<<Trk::pixelEndCap2)) ) d2Hit=1;
+ }
+/*************************************************************************************************************/
+
+ int NewVrtSecInclusiveTool::getIdHF(const xAOD::TrackParticle* TP ) const {
+ if( TP->isAvailable< ElementLink< xAOD::TruthParticleContainer> >( "truthParticleLink") ) {
+ const ElementLink<xAOD::TruthParticleContainer>& tplink =
+ TP->auxdata< ElementLink< xAOD::TruthParticleContainer > >("truthParticleLink");
+ if( !tplink.isValid() ) return 0;
+ if( TP->auxdata< float >( "truthMatchProbability" ) < 0.75 ) return 0;
+ if( (*tplink)->barcode() > 200000) return 0;
+ if( (*tplink)->hasProdVtx()){
+ if( (*tplink)->prodVtx()->nIncomingParticles()==1){
+ int PDGID1=0, PDGID2=0, PDGID3=0;
+ const xAOD::TruthParticle * parTP1=getPreviousParent(*tplink, PDGID1);
+ const xAOD::TruthParticle * parTP2=0;
+ int noBC1=notFromBC(PDGID1);
+ if(noBC1) parTP2 = getPreviousParent(parTP1, PDGID2);
+ int noBC2=notFromBC(PDGID2);
+ if(noBC2 && parTP2) getPreviousParent(parTP2, PDGID3);
+ int noBC3=notFromBC(PDGID3);
+ if(noBC1 && noBC2 && noBC3)return 0;
+ return 1; //This is a reconstructed track from B/C decays
+ } } }
+ return 0;
+ }
+
+ int NewVrtSecInclusiveTool::notFromBC(int PDGID) const {
+ int noBC=0;
+ if(PDGID<=0)return 1;
+ if(PDGID>600 && PDGID<4000)noBC=1;
+ if(PDGID<400 || PDGID>5600)noBC=1;
+ if(PDGID==513 || PDGID==523 || PDGID==533 || PDGID==543)noBC=1; //Remove tracks from B* (they are in PV)
+ if(PDGID==5114 || PDGID==5214 || PDGID==5224 || PDGID==5314 || PDGID==5324)noBC=1; //Remove tracks from B_Barions* (they are in PV)
+ //if(PDGID==413 || PDGID==423 || PDGID==433 )continue; //Keep tracks from D* (they are from B vertex)
+ //if(PDGID==4114 || PDGID==4214 || PDGID==4224 || PDGID==4314 || PDGID==4324)continue;
+ return noBC;
+ }
+ const xAOD::TruthParticle * NewVrtSecInclusiveTool::getPreviousParent(const xAOD::TruthParticle * child, int & ParentPDG) const {
+ ParentPDG=0;
+ if( child->hasProdVtx() ){
+ if( child->prodVtx()->nIncomingParticles()==1 ){
+ ParentPDG = abs((*(child->prodVtx()->incomingParticleLinks())[0])->pdgId());
+ return *(child->prodVtx()->incomingParticleLinks())[0];
+ }
+ }
+ return 0;
+ }
+
+
+ int NewVrtSecInclusiveTool::getG4Inter(const xAOD::TrackParticle* TP ) const {
+ if( TP->isAvailable< ElementLink< xAOD::TruthParticleContainer> >( "truthParticleLink") ) {
+ const ElementLink<xAOD::TruthParticleContainer>& tplink =
+ TP->auxdata< ElementLink< xAOD::TruthParticleContainer > >("truthParticleLink");
+ if( tplink.isValid() && (*tplink)->barcode()>200000) return 1;
+ }
+ return 0;
+ }
+ int NewVrtSecInclusiveTool::getMCPileup(const xAOD::TrackParticle* TP ) const {
+ if( TP->isAvailable< ElementLink< xAOD::TruthParticleContainer> >( "truthParticleLink") ) {
+ const ElementLink<xAOD::TruthParticleContainer>& tplink =
+ TP->auxdata< ElementLink< xAOD::TruthParticleContainer > >("truthParticleLink");
+ if( !tplink.isValid() ) return 1;
+ } else { return 1; }
+ return 0;
+ }
+
+ double NewVrtSecInclusiveTool::distToMatLayerSignificance(Vrt2Tr & Vrt) const
+ {
+ if(Vrt.fitVertex.perp()<20.) return 1.e9;
+ double normP=1./Vrt.momentum.P();
+ Amg::Vector3D momentumP(Vrt.momentum.Px()*normP,Vrt.momentum.Py()*normP,Vrt.momentum.Pz()*normP);
+ Amg::Vector3D momentumN=-momentumP;
+
+ const Trk::Layer * someLayer = 0;
+ const Trk::Layer * nextLayerP = 0;
+ const Trk::Layer * nextLayerN = 0;
+ auto volume = m_extrapolator->trackingGeometry()->lowestTrackingVolume(Vrt.fitVertex);
+ someLayer = volume->associatedLayer(Vrt.fitVertex);
+ auto material = someLayer->layerMaterialProperties();
+ if(material){
+ nextLayerP=someLayer;
+ } else {
+ nextLayerP = someLayer->nextLayer(Vrt.fitVertex,momentumP);
+ if(nextLayerP){ if(!nextLayerP->layerMaterialProperties())nextLayerP=0; }
+ nextLayerN = someLayer->nextLayer(Vrt.fitVertex,momentumN);
+ if(nextLayerN){ if(!nextLayerN->layerMaterialProperties())nextLayerN=0; }
+ }
+ momentumP *= 1.e5; //100GeV to have straight trajectory
+ double charge = 1.;
+ const Trk::Perigee pseudoVrtPart(Vrt.fitVertex, momentumP, charge, Vrt.fitVertex);
+
+ const Trk::TrackParameters * extrapParP=0; //along momentum
+ const Trk::TrackParameters * extrapParN=0; //backward
+ if(nextLayerP){ extrapParP = m_extrapolator->extrapolate(pseudoVrtPart,
+ nextLayerP->surfaceRepresentation(), Trk::anyDirection, false, Trk::nonInteractingMuon) ;}
+ if(nextLayerN){ extrapParN = m_extrapolator->extrapolate(pseudoVrtPart,
+ nextLayerN->surfaceRepresentation(), Trk::anyDirection, false, Trk::nonInteractingMuon) ;}
+
+ float distanceP=1.e9, distanceN=1.e9;
+ if(extrapParP)distanceP=PntPntDist(extrapParP->position(), Vrt.fitVertex);
+ if(extrapParN)distanceN=PntPntDist(extrapParN->position(), Vrt.fitVertex);
+ if(distanceP==1.e9 && distanceN==1.e9) return 1.e9;
+
+ //std::pair<const Trk::TrackParameters*,const Trk::Layer*> next=
+ // m_extrapolator->extrapolateToNextActiveLayer(pseudoVrtPart,Trk::anyDirection,true,Trk::pion) ;
+
+ double signif=1.e9;
+ std::vector<double> pntCovar={1.e-2,0.,1.e-2,0.,0.,4.e-2};
+ if(distanceP<distanceN)signif=VrtVrtDist(Vrt.fitVertex, Vrt.errorMatrix, extrapParP->position(), pntCovar);
+ else signif=VrtVrtDist(Vrt.fitVertex, Vrt.errorMatrix, extrapParN->position(), pntCovar);
+ delete extrapParP;
+ delete extrapParN;
+ return signif;
+ }
+
+
+
+
+
+} //end namespace
diff --git a/Reconstruction/VKalVrt/NewVrtSecInclusiveTool/src/VrtSecMulti.cxx b/Reconstruction/VKalVrt/NewVrtSecInclusiveTool/src/VrtSecMulti.cxx
new file mode 100755
index 0000000000000000000000000000000000000000..48342e3225eb6ec03c6d5f77817c7469fa764a60
--- /dev/null
+++ b/Reconstruction/VKalVrt/NewVrtSecInclusiveTool/src/VrtSecMulti.cxx
@@ -0,0 +1,645 @@
+/*
+ Copyright (C) 2002-2020 CERN for the benefit of the ATLAS collaboration
+*/
+///
+/// @author Vadim Kostyukhin <vadim.kostyukhin@cern.ch>
+///
+//
+// Header include
+#include "NewVrtSecInclusiveTool/NewVrtSecInclusiveTool.h"
+//-------------------------------------------------
+// Other stuff
+#include "AnalysisUtils/AnalysisMisc.h"
+#include "GeoPrimitives/GeoPrimitivesHelpers.h"
+#include "TrkVKalVrtFitter/TrkVKalVrtFitter.h"
+#include "MVAUtils/BDT.h"
+
+#include "boost/graph/bron_kerbosch_all_cliques.hpp"
+#include "TMath.h"
+#include "TH1.h"
+
+#include <algorithm>
+//
+
+
+namespace Rec{
+
+
+ std::vector<xAOD::Vertex*> NewVrtSecInclusiveTool::getVrtSecMulti( workVectorArrxAOD * xAODwrk,
+ const xAOD::Vertex & PrimVrt )
+ const
+ {
+
+ const double probVrtMergeLimit=0.01;
+
+ int inpNPart=0;
+ int i,j;
+ inpNPart=xAODwrk->inpTrk.size();
+ std::vector<const xAOD::NeutralParticle*> neutralPartDummy(0);
+ ATH_MSG_DEBUG( "getVrtSecMulti() called with NPart=" << inpNPart);
+
+ std::vector<xAOD::Vertex*> finalVertices(0);
+
+ if( inpNPart < 2 ) { return finalVertices;} // 0,1 track => nothing to do!
+
+
+ long int NTracks = 0;
+ selGoodTrkParticle( xAODwrk, PrimVrt);
+ NTracks = xAODwrk->listSelTracks.size();
+
+ if( NTracks < 2 ) { return finalVertices;} // 0,1 selected track => nothing to do!
+
+ ATH_MSG_DEBUG( "Number of selected tracks = " <<NTracks);
+
+ if(m_fillHist){ m_hb_ntrksel->Fill( (double) NTracks, m_w_1); }
+
+//
+// inpTrk[] - input track list
+// listSelTracks[] - list of good tracks in jet for vertex search
+//------------------------------------------------------------
+// Initial track list ready
+// Find 2track vertices
+//
+
+ select2TrVrt(xAODwrk->listSelTracks, PrimVrt);
+
+//---
+ ATH_MSG_DEBUG(" Defined edges in the graph="<< num_edges(*m_compatibilityGraph));
+
+//
+// m_Incomp[] - main vector of pointers for multivertex search
+//-----------------------------------------------------------------------------------------------------
+// Secondary track list is ready
+// Creation of initial vertex set
+//
+
+
+ std::vector<WrkVrt> *WrkVrtSet= new std::vector<WrkVrt>;
+ WrkVrt newvrt; newvrt.Good=true;
+ std::unique_ptr<Trk::IVKalState> state = m_fitSvc->makeState();
+ StatusCode sc; sc.setChecked();
+ long int NPTR=0, nth=2; // VK nth=2 to speed up PGRAPH when it's used
+
+
+ std::vector< std::vector<int> > allCliques;
+ bron_kerbosch_all_cliques(*m_compatibilityGraph, clique_visitor(allCliques));
+ for(int cq=0; cq<(int)allCliques.size();cq++){
+ newvrt.selTrk.clear();
+ NPTR=allCliques[cq].size();
+ for(i=0;i<NPTR;i++){ newvrt.selTrk.push_back(allCliques[cq][i]);}
+//==================================================
+ xAODwrk->tmpListTracks.clear();
+ TLorentzVector vpsum;
+ for(i=0;i<NPTR;i++) {
+ xAODwrk->tmpListTracks.push_back( xAODwrk->listSelTracks.at(newvrt.selTrk[i]) );
+ vpsum += xAODwrk->listSelTracks.at(newvrt.selTrk[i])->p4();
+ }
+ sc = m_fitSvc->VKalVrtFitFast(xAODwrk->tmpListTracks,newvrt.vertex,*state); /* Fast crude estimation*/
+
+ if( sc.isFailure() ) { /* No initial estimation */
+ m_fitSvc->setApproximateVertex(PrimVrt.x(), PrimVrt.y(), PrimVrt.z(),*state); /*Use as starting point*/
+ } else {
+ Amg::Vector3D vDist=newvrt.vertex-PrimVrt.position();
+ double MomVrtDir = vpsum.Px()*vDist.x() + vpsum.Py()*vDist.y() + vpsum.Pz()*vDist.z();
+ if( MomVrtDir>0. ) { /* Good initial estimation */
+ m_fitSvc->setApproximateVertex(newvrt.vertex.x(),newvrt.vertex.y(),newvrt.vertex.z(),*state); /*Use as starting point*/
+ }else{
+ m_fitSvc->setApproximateVertex(PrimVrt.x(), PrimVrt.y(), PrimVrt.z(),*state);
+ }
+ }
+// std::cout<<"FoundAppVrt="<<newvrt.vertex[0]<<", "<<newvrt.vertex[1]<<", "<<newvrt.vertex[2]<<'\n';
+ sc = StatusCode::FAILURE;
+ sc=m_fitSvc->VKalVrtFit(xAODwrk->tmpListTracks, neutralPartDummy,
+ newvrt.vertex, newvrt.vertexMom, newvrt.vertexCharge, newvrt.vertexCov,
+ newvrt.chi2PerTrk, newvrt.trkAtVrt, newvrt.chi2,
+ *state, false);
+
+// std::cout << "Res="<<newvrt.chi2<<", "<<NPTR<<", "<<newvrt.selTrk[0]<<", "<<newvrt.selTrk[1]<<'\n';
+ if( sc.isFailure() ) continue; /* Bad fit - goto next solution */
+ if(NPTR==2 && newvrt.chi2>10.) continue; /* Bad 2track vertex */
+ if(newvrt.chi2PerTrk.size()==2) newvrt.chi2PerTrk[0]=newvrt.chi2PerTrk[1]=newvrt.chi2/2.;
+ newvrt.Good = true;
+ newvrt.nCloseVrt = 0;
+ newvrt.dCloseVrt = 1000000.;
+ newvrt.projectedVrt=MomProjDist(newvrt.vertex, PrimVrt, newvrt.vertexMom); //3D SV-PV distance
+ WrkVrtSet->push_back(newvrt);
+ }
+//==================================================================================
+// boost::adjacency_list<boost::listS, boost::vecS, boost::undirectedS>::vertex_iterator vertexIt, vertexEnd;
+// boost::adjacency_list<boost::listS, boost::vecS, boost::undirectedS>::adjacency_iterator neighbourIt, neighbourEnd;
+// tie(vertexIt, vertexEnd) = vertices(*m_compatibilityGraph);
+// for (; vertexIt != vertexEnd; ++vertexIt) { std::cout << *vertexIt << " is connected with ";
+// tie(neighbourIt, neighbourEnd) = adjacent_vertices(*vertexIt, *m_compatibilityGraph);
+// for (; neighbourIt != neighbourEnd; ++neighbourIt) std::cout << *neighbourIt << " "; std::cout << "\n"; }
+//==================================================================================
+//========= Initial cleaning of solutions
+//-Remove worst track from vertices with very bad Chi2
+ bool disassembled=false;
+ int NSoluI=0;
+ do{
+ disassembled=false;
+ NSoluI=(*WrkVrtSet).size();
+ for(int iv=0; iv<NSoluI; iv++){
+ if(!(*WrkVrtSet)[iv].Good || (*WrkVrtSet)[iv].selTrk.size() == 2 ) continue;
+ if( TMath::Prob( (*WrkVrtSet)[iv].chi2, 2*(*WrkVrtSet)[iv].selTrk.size()-3) <1.e-3){
+ //printWrkSet(WrkVrtSet,"BigChi2Vertex present");
+ DisassembleVertex(WrkVrtSet,iv,xAODwrk->listSelTracks,*state);
+ disassembled=true;
+ } }
+ }while(disassembled);
+//-Modify too heavy vertices
+ for(auto & iv : (*WrkVrtSet)){
+ if( iv.vertexMom.M()>m_VrtMassLimit ) {
+ int it_bad=mostHeavyTrk(iv,xAODwrk->listSelTracks);
+ if(it_bad>-1){
+ iv.selTrk.erase( iv.selTrk.begin() + it_bad );
+ refitVertex(iv, xAODwrk->listSelTracks, *state, false);
+ } } }
+//-Remove vertices fully contained in other vertices
+ while( (*WrkVrtSet).size()>1 ){
+ int tmpN=(*WrkVrtSet).size(); int iv=0;
+ for(; iv<tmpN-1; iv++){ int jv=iv+1;
+ if(!(*WrkVrtSet)[iv].Good ) continue;
+ for(; jv<tmpN; jv++){
+ if(!(*WrkVrtSet)[jv].Good ) continue;
+ int nTCom=nTrkCommon( WrkVrtSet, iv, jv);
+ if( nTCom==(int)(*WrkVrtSet)[iv].selTrk.size()){ (*WrkVrtSet).erase((*WrkVrtSet).begin()+iv); break; }
+ else if( nTCom==(int)(*WrkVrtSet)[jv].selTrk.size()){ (*WrkVrtSet).erase((*WrkVrtSet).begin()+jv); break; }
+ } if(jv!=tmpN) break; // One vertex is erased. Restart check
+ } if(iv==tmpN-1) break; // No vertex deleted
+ }
+//
+//- Try to merge all vertices with common tracks
+ std::multimap<int,std::pair<int,int>> vrtWithCommonTrk;
+ std::multimap<int,std::pair<int,int>>::reverse_iterator icvrt;
+ do{
+ NSoluI=(*WrkVrtSet).size();
+ vrtWithCommonTrk.clear();
+ for(int iv=0; iv<NSoluI-1; iv++ ){ for(int jv=iv+1; jv<NSoluI; jv++){
+ if(!(*WrkVrtSet)[iv].Good || !(*WrkVrtSet)[jv].Good) continue;
+ int nTCom=nTrkCommon( WrkVrtSet, iv, jv); if(!nTCom)continue;
+ vrtWithCommonTrk.emplace(nTCom,std::make_pair(iv,jv));
+ } }
+ //============================== DEBUG output
+ //printWrkSet(WrkVrtSet,"InitialVrts");
+ //for(auto ku : vrtWithCommonTrk)std::cout<<" nCom="<<ku.first<<" v1="<<ku.second.first<<" v2="<<ku.second.second<<'\n';
+ //===========================================
+ for(icvrt=vrtWithCommonTrk.rbegin(); icvrt!=vrtWithCommonTrk.rend(); icvrt++){
+ int nTCom=(*icvrt).first;
+ int iv=(*icvrt).second.first;
+ int jv=(*icvrt).second.second;
+ int nTrkI=(*WrkVrtSet)[iv].selTrk.size();
+ int nTrkJ=(*WrkVrtSet)[jv].selTrk.size();
+ double probV=-1.;
+ probV=mergeAndRefitVertices( WrkVrtSet, iv, jv, newvrt, xAODwrk->listSelTracks, *state);
+ ////std::cout<<" ntcommon="<<(*icvrt).first<<" prb="<<probV<<" itrk="<<nTrkI<<" jtrk="<<nTrkJ<<'\n';
+ if(newvrt.vertexMom.M()>m_VrtMassLimit) continue; //Do not merge vertices if summary one is too heavy
+ if(probV<probVrtMergeLimit){
+ if(nTrkI==2 || nTrkJ==2 || nTCom<2) continue;
+ if((nTCom>nTrkI-nTCom || nTCom>nTrkJ-nTCom)){ //2 and more common tracks for NTr>=3 vertices. Merge anyway.
+ mergeAndRefitOverlapVertices( WrkVrtSet, iv, jv, xAODwrk->listSelTracks, *state);
+ break; //Vertex list is changed. Restart merging from scratch.
+ }
+ continue; //Continue merging loop
+ }
+ newvrt.Good = true;
+ (*WrkVrtSet).push_back(newvrt);
+ (*WrkVrtSet)[iv].Good=false;
+ (*WrkVrtSet)[jv].Good=false;
+ break; //Merging successful. Break merging loop and remake list of connected vertices
+ }
+ } while( icvrt != vrtWithCommonTrk.rend() );
+ if(m_fillHist){ int cvgood=0; for(int iv=0; iv<(int)(*WrkVrtSet).size(); iv++) if((*WrkVrtSet)[iv].Good)cvgood++;
+ m_hb_rawVrtN->Fill( (float)cvgood, m_w_1); }
+//
+//-Remove all bad vertices from the working set
+ int tmpV=0; while( tmpV<(int)(*WrkVrtSet).size() )if( !(*WrkVrtSet)[tmpV].Good ) { (*WrkVrtSet).erase((*WrkVrtSet).begin()+tmpV);} else {tmpV++;}
+ if((*WrkVrtSet).size()==0){ // No vertices at all
+ delete WrkVrtSet;
+ return finalVertices;
+ }
+ //------
+ //printWrkSet(WrkVrtSet,"Interm");
+ //------
+ for( auto &tmpV : (*WrkVrtSet) ) tmpV.projectedVrt=MomProjDist(tmpV.vertex, PrimVrt, tmpV.vertexMom ); //Setup ProjectedVrt
+//----------------------------------------------------------------------------
+// Here we have the overlapping solutions.
+// Vertices may have only 1 common track.
+// Now solution cleaning
+
+ int nGoodVertices=0; // Final number of good vertices
+ int n2trVrt=0; // N of vertices with 2 tracks
+ int nNtrVrt=0; // N vertices with 3 and more tracks
+ std::vector< std::deque<long int> > *TrkInVrt =new std::vector< std::deque<long int> >(NTracks);
+ double foundMaxT; long int SelectedTrack, SelectedVertex;
+ int foundV1, foundV2;
+
+ trackClassification( WrkVrtSet, TrkInVrt);
+
+ state = m_fitSvc->makeState();
+ while((foundMaxT=maxOfShared( WrkVrtSet, TrkInVrt, SelectedTrack, SelectedVertex))>0) {
+ //std::cout << "MAX="<<foundMaxT<<", "<<SelectedTrack<<", "<<SelectedVertex<<'\n';
+ //std::cout << "VRT="<<minVrtVrtDist( WrkVrtSet, foundV1, foundV2)<<", "<<foundV1<<", "<<foundV2<<'\n';
+ //printWrkSet(WrkVrtSet,"Iterat");
+
+ double foundMinVrtDst = 1000000.;
+ if(foundMaxT<m_trackDetachCut) foundMinVrtDst = minVrtVrtDist( WrkVrtSet, foundV1, foundV2);
+
+//Choice of action
+ if( foundMaxT<m_trackDetachCut && foundMinVrtDst<m_vertexMergeCut && nTrkCommon( WrkVrtSet, foundV1, foundV2)){
+ //if( foundMaxT<m_trackDetachCut && foundMinVrtDst<m_vertexMergeCut){
+ bool vrtMerged=false; //to check whether something is really merged
+ while(foundMinVrtDst<m_vertexMergeCut){
+ if(foundV1<foundV2) { int tmp=foundV1; foundV1=foundV2; foundV2=tmp;} /*Always drop vertex with smallest number*/
+ double probV=0.;
+ probV=mergeAndRefitVertices( WrkVrtSet, foundV1, foundV2, newvrt, xAODwrk->listSelTracks, *state);
+ if(probV>probVrtMergeLimit && newvrt.vertexMom.M()<m_VrtMassLimit){ // Good merged vertex found
+ double tstDst=MomProjDist(newvrt.vertex, PrimVrt, newvrt.vertexMom);
+ if(tstDst>0.){ // only positive vertex directions are accepted as merging result
+ std::swap((*WrkVrtSet)[foundV1],newvrt);
+ (*WrkVrtSet)[foundV1].projectedVrt=tstDst;
+ (*WrkVrtSet)[foundV2].Good=false; //Drop vertex
+ (*WrkVrtSet)[foundV2].selTrk.clear(); //Clean dropped vertex
+ vrtMerged=true;
+ }
+ }
+ (*WrkVrtSet)[foundV1].nCloseVrt=-1; (*WrkVrtSet)[foundV1].dCloseVrt=1000000.; //For minVrtVrtDistNext optimisation(!)
+ (*WrkVrtSet)[foundV2].nCloseVrt=-1; (*WrkVrtSet)[foundV2].dCloseVrt=1000000.; //Exclude given pair
+ foundMinVrtDst=minVrtVrtDistNext( WrkVrtSet, foundV1, foundV2); //Check other vertices
+ }
+ if(vrtMerged){
+ delete TrkInVrt;
+ TrkInVrt = new std::vector< std::deque<long int> >(NTracks);
+ trackClassification( WrkVrtSet, TrkInVrt);
+ continue; // Something was merged => goto next cycle. Otherwise break the found track-vertex link
+ }
+ }
+ removeTrackFromVertex(WrkVrtSet, TrkInVrt, SelectedTrack, SelectedVertex);
+ sc = refitVertex( WrkVrtSet, SelectedVertex, xAODwrk->listSelTracks, *state, false);
+ if( sc.isFailure() ) { (*WrkVrtSet)[SelectedVertex].Good=false; continue; } //bad vertex
+ (*WrkVrtSet)[SelectedVertex].projectedVrt=MomProjDist((*WrkVrtSet)[SelectedVertex].vertex, PrimVrt, (*WrkVrtSet)[SelectedVertex].vertexMom);
+ if( (*WrkVrtSet)[SelectedVertex].projectedVrt<0. && (*WrkVrtSet)[SelectedVertex].selTrk.size()==2 )
+ (*WrkVrtSet)[SelectedVertex].Good=false; // 2track vertex migrates behind PV - drop it.
+ //std::cout<<"Dropped v="<<SelectedVertex<<" dst="<<(*WrkVrtSet)[SelectedVertex].projectedVrt<<'\n';
+ }
+//
+// Final check/merge for close vertices
+//
+ double minDistVV=minVrtVrtDist( WrkVrtSet, foundV1, foundV2); //recalculate VV distances
+ while ( minDistVV < m_vertexMergeCut) {
+ if(foundV1<foundV2) { int tmp=foundV1; foundV1=foundV2; foundV2=tmp;}
+ double probV=0.;
+ probV=mergeAndRefitVertices( WrkVrtSet, foundV1, foundV2, newvrt, xAODwrk->listSelTracks, *state);
+ if(probV>probVrtMergeLimit && newvrt.vertexMom.M()<m_VrtMassLimit){ // Good merged vertex found
+ double tstDst=MomProjDist(newvrt.vertex, PrimVrt, newvrt.vertexMom);
+ if(tstDst>0.){ // only positive vertex directions are accepted as merging result
+ std::swap((*WrkVrtSet)[foundV1],newvrt);
+ (*WrkVrtSet)[foundV1].projectedVrt=tstDst;
+ (*WrkVrtSet)[foundV2].Good=false; //Drop vertex
+ (*WrkVrtSet)[foundV2].selTrk.clear(); //Clean dropped vertex
+ }
+ }
+ (*WrkVrtSet)[foundV1].nCloseVrt=-1; (*WrkVrtSet)[foundV1].dCloseVrt=1000000.; //For minVrtVrtDistNext optimisation(!)
+ (*WrkVrtSet)[foundV2].nCloseVrt=-1; (*WrkVrtSet)[foundV2].dCloseVrt=1000000.; //Exclude given pair
+ minDistVV=minVrtVrtDistNext( WrkVrtSet, foundV1, foundV2);
+ }
+//
+// Try to improve vertices with big Chi2 if something went wrong. Just precaution.
+ for(int iv=0; iv<(int)WrkVrtSet->size(); iv++) {
+ if(!(*WrkVrtSet)[iv].Good ) continue; //don't work on vertex which is already bad
+ if( (*WrkVrtSet)[iv].selTrk.size()<3 ) continue;
+ double tmpProb=TMath::Prob( (*WrkVrtSet)[iv].chi2, 2*(*WrkVrtSet)[iv].selTrk.size()-3 ); //Chi2 of the original vertex
+ if(tmpProb<m_globVrtProbCut){
+ tmpProb=improveVertexChi2( WrkVrtSet, iv, xAODwrk->listSelTracks, *state, false);
+ if(tmpProb<m_globVrtProbCut)(*WrkVrtSet)[iv].Good=false;
+ (*WrkVrtSet)[iv].projectedVrt=MomProjDist((*WrkVrtSet)[iv].vertex, PrimVrt, (*WrkVrtSet)[iv].vertexMom);
+ }
+ }
+ //
+ //printWrkSet(WrkVrtSet,"Joined");
+//--------------------------------------------------------------------------------------------------------
+// Final vertex selection/cleaning
+//
+ state = m_fitSvc->makeState();
+ double signif3D=0., signif2D=0., vProb=0.;
+ //double Dist3D=0;
+
+//--------- Start with 1-track vertices
+//=First check if the track was detached from a multitrack vertex. If so - reattach.
+ for(auto &ntrVrt : (*WrkVrtSet)){ if(!ntrVrt.Good || ntrVrt.selTrk.size()<=1) continue;
+ for(auto &onetVrt : (*WrkVrtSet)){ if(!onetVrt.Good || onetVrt.selTrk.size()!=1) continue;
+ if(ntrVrt.detachedTrack==onetVrt.selTrk[0]){
+ WrkVrt newV(ntrVrt); newV.selTrk.push_back(ntrVrt.detachedTrack);
+ vProb = refitVertex( newV, xAODwrk->listSelTracks, *state, true);
+ if( vProb>probVrtMergeLimit && newV.vertexMom.M()<m_VrtMassLimit ){ onetVrt.Good=false; ntrVrt=newV; ntrVrt.detachedTrack=-1;}
+ break;
+ } } }
+//=Recheck if the track was detached from a 2track vertex. If so - reattach.
+ for(auto &iVrt : (*WrkVrtSet)){ if(!iVrt.Good || iVrt.selTrk.size()!=1) continue;
+ for(auto &jVrt : (*WrkVrtSet)){ if(!jVrt.Good || jVrt.selTrk.size()!=1) continue;
+ if(iVrt.detachedTrack==jVrt.selTrk[0]){
+ WrkVrt newV(iVrt); newV.selTrk.push_back(iVrt.detachedTrack);
+ vProb = refitVertex( newV, xAODwrk->listSelTracks, *state, true);
+ if( vProb>probVrtMergeLimit && newV.vertexMom.M()<m_VrtMassLimit ){ jVrt.Good=false; iVrt=newV; iVrt.detachedTrack=-1;}
+ break;
+ } } }
+ for(auto &ntrVrt : (*WrkVrtSet)){ if(!ntrVrt.Good || ntrVrt.selTrk.size()<=1) continue;
+ for(auto tr : ntrVrt.selTrk){
+ for(auto &onetVrt : (*WrkVrtSet)){ if(!onetVrt.Good || onetVrt.selTrk.size()!=1) continue;
+ if(onetVrt.detachedTrack==tr){
+ WrkVrt newV(ntrVrt); newV.selTrk.push_back(onetVrt.selTrk[0]);
+ vProb = refitVertex( newV, xAODwrk->listSelTracks, *state, true);
+ if( vProb>probVrtMergeLimit && newV.vertexMom.M()<m_VrtMassLimit ){ onetVrt.Good=false; ntrVrt=newV;}
+ } } } }
+ for(auto & curVrt : (*WrkVrtSet) ) {
+ if(!curVrt.Good ) continue; //don't work on vertex which is already bad
+ if(curVrt.selTrk.size() != 1) continue;
+ curVrt.Good=false; // Make them bad by default
+ if(m_multiWithOneTrkVrt){ /* 1track vertices left unassigned from good 2tr vertices */
+ VrtVrtDist(PrimVrt,curVrt.vertex, curVrt.vertexCov, signif3D); //VK non-projected signif3D is worse
+ double tmpProb=TMath::Prob( curVrt.chi2, 1); //Chi2 of the original 2tr vertex
+ if(tmpProb>0.01){ /* accept only good tracks coming from good 2tr vertex*/
+ std::vector<double> impact,impactError; double signif3DP = 0;
+ signif3DP=m_fitSvc->VKalGetImpact(xAODwrk->listSelTracks[curVrt.selTrk[0]],PrimVrt.position(), 1, impact, impactError, *state);
+ if(m_fillHist){
+ m_hb_diffPS->Fill( signif3DP, m_w_1);
+ m_hb_sig3D1tr->Fill( signif3D, m_w_1);
+ }
+ if( signif3DP>2.*m_trkSigCut && signif3D>m_selVrtSigCut) curVrt.Good=true; // accept only tracks which are far from PV
+ }
+ }
+ }
+//----- Vertices with >1 tracks
+ for(int iv=0; iv<(int)WrkVrtSet->size(); iv++) {
+ WrkVrt & curVrt=(*WrkVrtSet)[iv];
+ if(!curVrt.Good ) continue; //don't work on vertex which is already bad
+ nth=(*WrkVrtSet)[iv].selTrk.size(); if(nth == 1) continue; // 1track vertices are treated already
+ VrtVrtDist(PrimVrt,curVrt.vertex, curVrt.vertexCov, signif3D);
+ xAODwrk->tmpListTracks.resize(nth);
+ for(i=0;i<nth;i++) {
+ xAODwrk->tmpListTracks[i]=xAODwrk->listSelTracks[curVrt.selTrk[i]];
+ }
+ (*WrkVrtSet)[iv].Good = false; /* Make all vertices bad for futher selection */
+ if(nth <= 1) continue; /* Definitely bad vertices */
+ if((*WrkVrtSet)[iv].projectedVrt<0.) continue; /* Remove vertices behind primary one */
+ if( TMath::Prob( curVrt.chi2, 2*nth-3)<m_globVrtProbCut) continue; /* Bad Chi2 of refitted vertex */
+//-----------------------------------------------------------------------------------------
+ if(m_fillHist){
+ if(nth==2 && curVrt.vertexCharge==0) m_hb_massPiPi1->Fill(curVrt.vertexMom.M(), m_w_1);
+ m_hb_sig3DTot->Fill( signif3D, m_w_1);
+ if(nth==2)m_hb_sig3D2tr->Fill( signif3D, m_w_1);
+ if(nth >2)m_hb_sig3DNtr->Fill( signif3D, m_w_1);
+ }
+//
+//--- Check V0s and conversions ???
+/* if(nth==2 && curVrt.vertexCharge==0 && curVrt.detachedTrack<0){
+ double mass_PiPi = curVrt.vertexMom.M();
+ double mass_PPi = massV0(curVrt.trkAtVrt,m_massP,m_massPi);
+ double mass_EE = massV0(curVrt.trkAtVrt,m_massE,m_massE);
+ if(m_fillHist){ m_hb_massPiPi->Fill( mass_PiPi, m_w_1);
+ m_hb_massPPi ->Fill( mass_PPi, m_w_1);
+ if( curVrt.vertex.perp()>20.)m_hb_massEE ->Fill( mass_EE, m_w_1); }
+ if( fabs(mass_PiPi-m_massK0) < 22.) continue;
+ if( fabs(mass_PPi-m_massLam) < 8.) continue;
+ if( mass_EE < 60. && curVrt.vertex.perp() > 20.) continue;
+ }
+*/
+//---
+ if(signif3D<m_selVrtSigCut+1.) continue; //Main PV-SV distance quality cut
+ if(curVrt.vertex.perp() > m_maxSVRadiusCut) continue; // Too far from interaction point
+//---
+ curVrt.Good = true; /* Vertex is absolutely good */
+ }
+//-------------------------------------------
+// Debug ntuple filling and BDT application
+//
+ std::vector<double> impact,impactError;
+ for(int iv=0; iv<(int)WrkVrtSet->size(); iv++) {
+ WrkVrt & curVrt=(*WrkVrtSet)[iv];
+ nth=curVrt.selTrk.size();
+ if(!curVrt.Good ) continue; //don't work on vertex which is already bad
+ double minPtT=1.e6, minSig3DT=1.e6, maxSig3DT=0.;
+ int ntrkBC=0,ntrkI=0,sumIBLHits=0,sumBLHits=0;
+ for(i=0;i<nth;i++) {
+ j=curVrt.selTrk[i]; /*Track number*/
+ minPtT=std::min( minPtT, xAODwrk->listSelTracks[j]->pt());
+ m_fitSvc->VKalGetImpact(xAODwrk->listSelTracks[j], PrimVrt.position(), 1, impact, impactError);
+ double SigR2 = impact[0]*impact[0]/impactError[0];
+ double SigZ2 = impact[1]*impact[1]/impactError[2];
+ minSig3DT=std::min( minSig3DT, sqrt( SigR2 + SigZ2) );
+ maxSig3DT=std::max( maxSig3DT, sqrt( SigR2 + SigZ2) );
+ sumIBLHits += std::max(getIBLHit(xAODwrk->listSelTracks[j]),0);
+ sumBLHits += std::max(getBLHit(xAODwrk->listSelTracks[j]),0);
+ if(m_fillHist && m_curTup) {
+ ntrkBC += getIdHF(xAODwrk->listSelTracks[j]);
+ ntrkI += getG4Inter(xAODwrk->listSelTracks[j]);
+ }
+ }
+ float vProb=TMath::Prob(curVrt.chi2, 2*nth-3);
+ float cosSVPVM=ProjSV_PV(curVrt.vertex, PrimVrt, curVrt.vertexMom);
+ float vrtR=curVrt.vertex.perp();
+ TLorentzVector SVPV(curVrt.vertex.x()-PrimVrt.x(),curVrt.vertex.y()-PrimVrt.y(),curVrt.vertex.z()-PrimVrt.z(), 10.);
+ if(m_fillHist){
+ if( m_curTup && nth>1 ){
+ VrtVrtDist(PrimVrt,curVrt.vertex, curVrt.vertexCov, signif3D);
+ float Dist2D=VrtVrtDist2D(PrimVrt,curVrt.vertex, curVrt.vertexCov, signif2D);
+ m_curTup->NVrtTrk [m_curTup->nNVrt] = nth;
+ m_curTup->NVrtTrkHF [m_curTup->nNVrt] = ntrkBC;
+ m_curTup->NVrtTrkI [m_curTup->nNVrt] = ntrkI;
+ m_curTup->NVrtProb [m_curTup->nNVrt] = vProb;
+ m_curTup->NVrtSig3D [m_curTup->nNVrt] = signif3D;
+ m_curTup->NVrtSig2D [m_curTup->nNVrt] = signif2D;
+ m_curTup->NVrtDist2D[m_curTup->nNVrt] = vrtR<20. ? Dist2D : vrtR;
+ m_curTup->NVrtM [m_curTup->nNVrt] = curVrt.vertexMom.M();
+ m_curTup->NVrtPt [m_curTup->nNVrt] = curVrt.vertexMom.Pt();
+ m_curTup->NVrtEta [m_curTup->nNVrt] = SVPV.Eta();
+ m_curTup->NVrtIBL [m_curTup->nNVrt] = sumIBLHits;
+ m_curTup->NVrtBL [m_curTup->nNVrt] = sumBLHits;
+ m_curTup->NVrtSinSPM[m_curTup->nNVrt] = sqrt(1.-cosSVPVM*cosSVPVM);
+ m_curTup->NVrtCh [m_curTup->nNVrt] = curVrt.vertexCharge;
+ m_curTup->NVMinPtT [m_curTup->nNVrt] = minPtT;
+ m_curTup->NVMinS3DT [m_curTup->nNVrt] = minSig3DT;
+ m_curTup->NVrtBDT [m_curTup->nNVrt] = 1.1;
+ if( m_curTup->nNVrt < DevTuple::maxNVrt-1 )m_curTup->nNVrt++;
+ }
+ }
+//-------------------BDT based rejection
+ if(nth==2){
+ if(curVrt.vertexMom.Pt() > m_Vrt2TrPtLimit){ curVrt.Good = false; continue; }
+ float rhit0=xAODwrk->listSelTracks[curVrt.selTrk[0]]->radiusOfFirstHit();
+ float rhit1=xAODwrk->listSelTracks[curVrt.selTrk[1]]->radiusOfFirstHit();
+ VrtVrtDist(PrimVrt,curVrt.vertex, curVrt.vertexCov, signif3D);
+ float Dist2D=VrtVrtDist2D(PrimVrt,curVrt.vertex, curVrt.vertexCov, signif2D);
+ std::vector<float> VARS(10);
+ VARS[0]=vProb;
+ VARS[1]=log(curVrt.vertexMom.Pt());
+ VARS[2]=log(std::max(minPtT,m_cutPt));
+ VARS[3]=log(vrtR<20. ? Dist2D : vrtR);
+ VARS[4]=log(std::max(minSig3DT,m_trkSigCut));
+ VARS[5]=log(maxSig3DT);
+ VARS[6]=curVrt.vertexMom.M();
+ VARS[7]=sqrt(fabs(1.-cosSVPVM*cosSVPVM));
+ VARS[8]=SVPV.Eta();
+ VARS[9]=std::max(rhit0,rhit1);
+ //VARS[9]=sumIBLHits;
+ //VARS[10]=sumBLHits;
+ //VARS[4]=std::max(signif3D,m_selVrtSigCut);
+ float wgtSelect=m_SV2T_BDT->GetGradBoostMVA(VARS);
+ //std::vector<float> weights=m_SV2T_BDT->GetMultiResponse(VARS,3);
+ //float wgtSelect=weights[0];
+ if(m_fillHist){
+ m_hb_fakeSVBDT->Fill(wgtSelect,1.);
+ if( m_curTup ) m_curTup->NVrtBDT[m_curTup->nNVrt-1] = wgtSelect;
+ }
+ if(wgtSelect<m_v2tFinBDTCut) curVrt.Good = false;
+ }
+ }
+//
+//--Final cleaning of the 1-track vertices set. Must be behind all other cleanings.
+//-- Debug ntuple for 1track vertex is filled here as well
+//
+ if(m_multiWithOneTrkVrt) Clean1TrVertexSet(WrkVrtSet);
+ if(m_fillHist && m_curTup && m_multiWithOneTrkVrt){
+ for(auto & V : (*WrkVrtSet)) {
+ nth=V.selTrk.size();
+ if( !V.Good || nth != 1 ) continue; // Good 1track vertices
+ m_fitSvc->VKalGetImpact(xAODwrk->listSelTracks[V.selTrk[0]], PrimVrt.position(), 1, impact, impactError);
+ int ntrkBC = getIdHF(xAODwrk->listSelTracks[V.selTrk[0]]);
+ double SigR2 = impact[0]*impact[0]/impactError[0];
+ double SigZ2 = impact[1]*impact[1]/impactError[2];
+ float Dist2D=VrtVrtDist2D(PrimVrt,V.vertex, V.vertexCov, signif2D);
+ m_curTup->NVrtTrk [m_curTup->nNVrt] = nth;
+ m_curTup->NVrtTrkHF [m_curTup->nNVrt] = ntrkBC;
+ m_curTup->NVrtProb [m_curTup->nNVrt] = TMath::Prob(V.chi2, 1);
+ m_curTup->NVrtSig3D [m_curTup->nNVrt] = 0.;
+ m_curTup->NVrtSig2D [m_curTup->nNVrt] = signif2D;
+ m_curTup->NVrtDist2D[m_curTup->nNVrt] = Dist2D;
+ m_curTup->NVrtM [m_curTup->nNVrt] = V.vertexMom.M();
+ m_curTup->NVrtPt [m_curTup->nNVrt] = V.vertexMom.Pt();
+ m_curTup->NVrtSinSPM[m_curTup->nNVrt] = 0.;
+ m_curTup->NVrtCh [m_curTup->nNVrt] = V.vertexCharge;
+ m_curTup->NVMinPtT [m_curTup->nNVrt] = xAODwrk->listSelTracks[V.selTrk[0]]->pt();
+ m_curTup->NVMinS3DT [m_curTup->nNVrt] = sqrt(SigR2 + SigZ2);
+ m_curTup->NVrtBDT [m_curTup->nNVrt] = 1.1;
+ if( m_curTup->nNVrt < DevTuple::maxNVrt-1 )m_curTup->nNVrt++;
+ } }
+//-------------------------------------------
+//Checks
+ std::vector<WrkVrt> GoodVertices(0);
+ nGoodVertices=0; // Final number of good vertices
+ n2trVrt=nNtrVrt=0; // N of vertices with different amount of tracks
+ for(auto & iv : (*WrkVrtSet) ) {
+ nth=iv.selTrk.size(); if(nth == 0) continue; /* Definitely bad vertices */
+ if( iv.Good) {
+ nGoodVertices++;
+ GoodVertices.emplace_back(iv); /* add it */
+ if(nth==2)n2trVrt++;
+ if(nth>=3)nNtrVrt++;
+ }
+ }
+ if(nGoodVertices == 0 || (n2trVrt+nNtrVrt)==0){ // No good vertices at all
+ delete WrkVrtSet; delete TrkInVrt;
+ if(m_fillHist && m_curTup ) m_curTup->nNVrt=0;
+ return finalVertices;
+ }
+//
+//sorting
+ while(1){ bool swapDone=false; // Sort on XY distance from (0.,0.)
+ for(int iv=0; iv<(int)GoodVertices.size()-1; iv++) {
+ if( GoodVertices[iv].vertex.perp() > GoodVertices[iv+1].vertex.perp()){
+ std::swap( GoodVertices[iv], GoodVertices[iv+1] );
+ swapDone=true;
+ }
+ } if(!swapDone) break;
+ }
+ while(1){ bool swapDone=false; // Then sort on Projected dist if R<20.
+ for(int iv=0; iv<(int)GoodVertices.size()-1; iv++) {
+ if(GoodVertices[iv+1].vertex.perp() > 400.) continue;
+ if(GoodVertices[iv].projectedVrt > GoodVertices[iv+1].projectedVrt){
+ std::swap( GoodVertices[iv], GoodVertices[iv+1] );
+ swapDone=true;
+ }
+ } if(!swapDone) break;
+ }
+ if(nGoodVertices>1){
+ if( GoodVertices[1].vertexMom.M()-GoodVertices[0].vertexMom.M() > 5000.) std::swap( GoodVertices[0], GoodVertices[1] );
+ }
+ if(m_fillHist){m_hb_distVV->Fill( minVrtVrtDist( WrkVrtSet, foundV1, foundV2), m_w_1); }
+//-------------------------------------------
+// Saving and final cleaning
+//
+ nGoodVertices=0; // Final number of good vertices
+ int n1trVrt=0; // Final number of good 1-track vertices
+ TLorentzVector VertexMom;
+ for(int iv=0; iv<(int)GoodVertices.size(); iv++) {
+ nth=GoodVertices[iv].selTrk.size();
+ xAODwrk->tmpListTracks.clear();
+ for(i=0;i<nth;i++) {
+ j=GoodVertices[iv].selTrk[i]; /*Track number*/
+ xAODwrk->tmpListTracks.push_back( xAODwrk->listSelTracks[j] );
+ }
+ if(m_fillHist){
+ m_hb_totmass->Fill(GoodVertices[iv].vertexMom.M(), m_w_1);
+ m_hb_r2d->Fill( GoodVertices[iv].vertex.perp(), m_w_1); }
+// xAOD::Vertex creation-----------------------------
+ xAOD::Vertex * tmpVertex=new xAOD::Vertex();
+ tmpVertex->makePrivateStore();
+ tmpVertex->setPosition(GoodVertices[iv].vertex);
+ int NERMS=GoodVertices[iv].vertexCov.size();
+ NERMS=6;
+ std::vector<float> floatErrMtx; floatErrMtx.resize(NERMS);
+ for(int i=0; i<NERMS; i++) floatErrMtx[i]=GoodVertices[iv].vertexCov[i];
+ tmpVertex->setCovariance(floatErrMtx);
+ tmpVertex->setFitQuality(GoodVertices[iv].chi2, (float)(nth*2-3));
+
+ std::vector<Trk::VxTrackAtVertex> & tmpVTAV=tmpVertex->vxTrackAtVertex(); tmpVTAV.clear();
+ for(int ii=0; ii<nth; ii++) {
+ AmgSymMatrix(5) *CovMtxP=new AmgSymMatrix(5); (*CovMtxP).setIdentity();
+ Trk::Perigee * tmpMeasPer = new Trk::Perigee( 0.,0., GoodVertices[iv].trkAtVrt[ii][0],
+ GoodVertices[iv].trkAtVrt[ii][1],
+ GoodVertices[iv].trkAtVrt[ii][2],
+ Trk::PerigeeSurface(GoodVertices[iv].vertex), CovMtxP );
+ tmpVTAV.push_back( Trk::VxTrackAtVertex( 1., tmpMeasPer) );
+ if(xAODwrk){ //No way to store a link to Rec::TrackParticleContainer
+ ElementLink<xAOD::TrackParticleContainer> TEL; TEL.setElement( xAODwrk->tmpListTracks[ii] );
+ const xAOD::TrackParticleContainer* cont = (const xAOD::TrackParticleContainer* ) (xAODwrk->tmpListTracks[ii]->container() );
+ TEL.setStorableObject(*cont);
+ tmpVertex->addTrackAtVertex(TEL,1.);
+ }
+ }
+ finalVertices.push_back(tmpVertex);
+//-----
+ nGoodVertices++;
+ if(nth==1)n1trVrt++;
+//-----
+ VertexMom += GoodVertices[iv].vertexMom;
+ }
+ if(m_fillHist){
+ m_hb_goodvrtN->Fill( nGoodVertices+0.1, m_w_1);
+ m_hb_goodvrt1N->Fill( n1trVrt+0.1, m_w_1);
+ }
+
+
+ if(nGoodVertices == 0){
+ delete WrkVrtSet;
+ delete TrkInVrt;
+ return finalVertices;
+ }
+
+//-----------------------------------------------------------------------------------
+// Saving of results
+//
+//
+
+
+ delete WrkVrtSet; delete TrkInVrt;
+
+ return finalVertices;
+
+
+ }
+
+
+} //end namespace
+
+
diff --git a/Reconstruction/VKalVrt/NewVrtSecInclusiveTool/src/components/NewVrtSecInclusiveTool_entries.cxx b/Reconstruction/VKalVrt/NewVrtSecInclusiveTool/src/components/NewVrtSecInclusiveTool_entries.cxx
new file mode 100644
index 0000000000000000000000000000000000000000..245d17b9d69d448ffb1fbb235aebb787298a3f1f
--- /dev/null
+++ b/Reconstruction/VKalVrt/NewVrtSecInclusiveTool/src/components/NewVrtSecInclusiveTool_entries.cxx
@@ -0,0 +1,7 @@
+#include "NewVrtSecInclusiveTool/NewVrtSecInclusiveTool.h"
+
+using namespace Rec;
+
+DECLARE_COMPONENT( NewVrtSecInclusiveTool )
+
+
diff --git a/Trigger/TrigMonitoring/TrigIDtrkMonitoring/python/TrigIDPhysValMonitoringConfig.py b/Trigger/TrigMonitoring/TrigIDtrkMonitoring/python/TrigIDPhysValMonitoringConfig.py
index f08a8b013d43287ba3f9a28a48ea6aaf427a776d..1fc3ac1d3b0ff4064af328641f301137e1c5c23e 100644
--- a/Trigger/TrigMonitoring/TrigIDtrkMonitoring/python/TrigIDPhysValMonitoringConfig.py
+++ b/Trigger/TrigMonitoring/TrigIDtrkMonitoring/python/TrigIDPhysValMonitoringConfig.py
@@ -123,7 +123,9 @@ def TrigIDPhysValMonitoringTool( legacy_monitoring=False ):
if mt_chains:
chainnames = [
"HLT_mu.*idperf.*:key=HLT_IDTrack_Muon_FTF:roi=HLT_Roi_L2SAMuon",
- "HLT_mu.*idperf.*:key=HLT_IDTrack_Muon_IDTrig:roi=HLT_Roi_L2SAMuon"
+ "HLT_mu.*idperf.*:key=HLT_IDTrack_Muon_FTF:roi=HLT_Roi_L2SAMuonForEF",
+ "HLT_mu.*idperf.*:key=HLT_IDTrack_Muon_IDTrig:roi=HLT_Roi_L2SAMuon",
+ "HLT_mu.*idperf.*:key=HLT_IDTrack_Muon_IDTrig:roi=HLT_Roi_L2SAMuonForEF"
]
else:
chainnames = [
diff --git a/Trigger/TrigMonitoring/TrigIDtrkMonitoring/python/TrigIDtrkMonitoringConfig.py b/Trigger/TrigMonitoring/TrigIDtrkMonitoring/python/TrigIDtrkMonitoringConfig.py
index 48ebb4505f53451c9e6c1566dc954ed22271c012..55bd062918373294cd0de0d2bd72cab802146d23 100644
--- a/Trigger/TrigMonitoring/TrigIDtrkMonitoring/python/TrigIDtrkMonitoringConfig.py
+++ b/Trigger/TrigMonitoring/TrigIDtrkMonitoring/python/TrigIDtrkMonitoringConfig.py
@@ -178,10 +178,9 @@ def TrigIDtrkMonitoringTool( legacy_monitoring=False ):
tidamuon.ntupleChainNames += [
"Offline",
"HLT_mu.*idperf.*:key=HLT_IDTrack_Muon_FTF:roi=HLT_Roi_L2SAMuon",
+ "HLT_mu.*idperf.*:key=HLT_IDTrack_Muon_FTF:roi=HLT_Roi_L2SAMuonForEF",
"HLT_mu.*idperf.*:key=HLT_IDTrack_Muon_IDTrig:roi=HLT_Roi_L2SAMuon",
- "HLT_mu.*idperf.*:key=HLT_IDTrack_MuonIso_FTF",
- "HLT_mu.*idperf.*:key=HLT_IDTrack_MuonIso_IDTrig",
- "HLT_mu.*idperf.*:key=HLT_IDTrack_MuonLate_IDTrig"
+ "HLT_mu.*idperf.*:key=HLT_IDTrack_Muon_IDTrig:roi=HLT_Roi_L2SAMuonForEF"
]
else:
tidamuon.ntupleChainNames += [
diff --git a/Trigger/TrigSteer/DecisionHandling/src/RoRSeqFilter.cxx b/Trigger/TrigSteer/DecisionHandling/src/RoRSeqFilter.cxx
index a5775365246c8f1726071837f5f7c752df58f1ce..ecb012901a29d84399e670888d8512770cb15304 100644
--- a/Trigger/TrigSteer/DecisionHandling/src/RoRSeqFilter.cxx
+++ b/Trigger/TrigSteer/DecisionHandling/src/RoRSeqFilter.cxx
@@ -18,7 +18,7 @@ using TrigCompositeUtils::newDecisionIn;
RoRSeqFilter::RoRSeqFilter( const std::string& name,
ISvcLocator* pSvcLocator ) :
- ::AthAlgorithm( name, pSvcLocator )
+ ::AthReentrantAlgorithm( name, pSvcLocator )
{}
@@ -71,15 +71,17 @@ StatusCode RoRSeqFilter::initialize()
}
-StatusCode RoRSeqFilter::execute() {
+StatusCode RoRSeqFilter::execute(const EventContext& ctx) const {
ATH_MSG_DEBUG ( "Executing " << name() << "..." );
- auto inputHandles = m_inputKeys.makeHandles();
- auto outputHandles = m_outputKeys.makeHandles();
+ auto inputHandles = m_inputKeys.makeHandles(ctx);
+ auto outputHandles = m_outputKeys.makeHandles(ctx);
std::vector<std::string> inputNames({"exec", "anyvalid"});
std::vector<bool> inputStats({true, false}); // position 0 for number of execs, always true, bool at position 1 is set later
+ inputNames.reserve(inputHandles.size() + 2);
+ inputStats.reserve(inputHandles.size() + 2);
bool validInputs = false;
- for ( auto inputHandle: inputHandles ) {
+ for ( auto& inputHandle: inputHandles ) {
inputNames.push_back(inputHandle.name());
if( inputHandle.isValid() ) {// this is because input is implicit
validInputs = true;
@@ -94,7 +96,7 @@ StatusCode RoRSeqFilter::execute() {
Monitored::Group( m_monTool, inputStat, inputName );
if (!validInputs) {
- setFilterPassed(false);
+ setFilterPassed(false, ctx);
ATH_MSG_DEBUG ( "No valid inputs found, filter failed. Return...." );
return StatusCode::SUCCESS;
}
@@ -108,25 +110,23 @@ StatusCode RoRSeqFilter::execute() {
ATH_MSG_DEBUG( "Recording " << m_outputKeys[ 0 ].key() );
createAndStore(outputHandles[0]);
DecisionContainer* output = outputHandles[0].ptr();
- for ( auto inputKey: m_inputKeys ) {
- auto inputHandle = SG::makeHandle( inputKey );
- if( inputHandle.isValid() )
+ for ( auto& inputHandle: inputHandles) {
+ if( inputHandle.isValid() ) {
passCounter += copyPassing( *inputHandle, *output );
+ }
}
outputIndex++;
} else { // Not merging inputs
- for ( auto inputKey: m_inputKeys ) {
- // already made handles, so this code could be simplified to a loop over inputHandles.
- auto inputHandle = SG::makeHandle( inputKey );
+ for ( auto& inputHandle: inputHandles ) {
if( not inputHandle.isValid() ) {
- ATH_MSG_DEBUG( "InputHandle "<< inputKey.key() <<" not present" );
+ ATH_MSG_DEBUG( "InputHandle "<< inputHandle.key() <<" not present" );
} else if ( inputHandle->empty() ) {
- ATH_MSG_DEBUG( "InputHandle "<< inputKey.key() <<" contains all rejected decisions, skipping" );
+ ATH_MSG_DEBUG( "InputHandle "<< inputHandle.key() <<" contains all rejected decisions, skipping" );
} else {
- ATH_MSG_DEBUG( "Checking inputHandle: "<< inputKey.key() <<" has " << inputHandle->size() <<" elements");
+ ATH_MSG_DEBUG( "Checking inputHandle: "<< inputHandle.key() <<" has " << inputHandle->size() <<" elements");
createAndStore(outputHandles[outputIndex]);
DecisionContainer* output = outputHandles[outputIndex].ptr();
passCounter += copyPassing( *inputHandle, *output );
@@ -134,13 +134,12 @@ StatusCode RoRSeqFilter::execute() {
}
outputIndex++; // Keep the mapping of inputKey<->outputKey correct
}
-
}
- setFilterPassed( passCounter != 0 );
- ATH_MSG_DEBUG( "Filter " << ( filterPassed() ? "passed" : "rejected") <<"; creating "<< outputIndex<<" valid outDecisions DH:");
+ setFilterPassed( passCounter != 0, ctx );
+ ATH_MSG_DEBUG( "Filter " << ( filterPassed(ctx) ? "passed" : "rejected") <<"; creating "<< outputIndex<<" valid outDecisions DH:");
if (msgLvl(MSG::DEBUG)){
- for (auto output: outputHandles){
+ for (auto& output: outputHandles){
if( output.isValid() ) ATH_MSG_DEBUG(" "<<output.key());
}
}
diff --git a/Trigger/TrigSteer/DecisionHandling/src/RoRSeqFilter.h b/Trigger/TrigSteer/DecisionHandling/src/RoRSeqFilter.h
index cfb712abc0bcbcac2fd0c344e0d5ee5f2cc685ed..ac24edb149317fb1dd4986f82cfd09575b60de35 100644
--- a/Trigger/TrigSteer/DecisionHandling/src/RoRSeqFilter.h
+++ b/Trigger/TrigSteer/DecisionHandling/src/RoRSeqFilter.h
@@ -7,7 +7,7 @@
#include <string>
#include <set>
-#include "AthenaBaseComps/AthAlgorithm.h"
+#include "AthenaBaseComps/AthReentrantAlgorithm.h"
#include "AthContainers/ConstDataVector.h"
#include "TrigCompositeUtils/TrigCompositeUtils.h"
#include "TrigCompositeUtils/HLTIdentifier.h"
@@ -44,7 +44,7 @@
**/
class RoRSeqFilter
- : public ::AthAlgorithm
+ : public ::AthReentrantAlgorithm
{
public:
RoRSeqFilter( const std::string& name, ISvcLocator* pSvcLocator );
@@ -58,7 +58,7 @@ class RoRSeqFilter
* @brief Apply this filter in-between Steps of trigger execution. Fully implicit inputs, requires Control Flow to unlock.
* will signal a negative filter result to the Scheduler if zero chains remain active upon termination.
**/
- virtual StatusCode execute() override;
+ virtual StatusCode execute(const EventContext& ctx) const override;
private:
RoRSeqFilter();
diff --git a/Trigger/TrigValidation/TrigInDetValidation/test/test_trigID_all_ttbar_pu80_mp_grid.py b/Trigger/TrigValidation/TrigInDetValidation/test/test_trigID_all_ttbar_pu80_mt_grid.py
similarity index 97%
rename from Trigger/TrigValidation/TrigInDetValidation/test/test_trigID_all_ttbar_pu80_mp_grid.py
rename to Trigger/TrigValidation/TrigInDetValidation/test/test_trigID_all_ttbar_pu80_mt_grid.py
index e417f55a46a1345fbc3f1c68199573342220fb0a..19ba137045c02c905eac2614fc66a40403af595d 100755
--- a/Trigger/TrigValidation/TrigInDetValidation/test/test_trigID_all_ttbar_pu80_mp_grid.py
+++ b/Trigger/TrigValidation/TrigInDetValidation/test/test_trigID_all_ttbar_pu80_mt_grid.py
@@ -92,9 +92,8 @@ preexec_all = ';'.join([
rdo2aod = ExecStep.ExecStep()
rdo2aod.type = 'Reco_tf'
rdo2aod.max_events = 1000 # TODO: 2000 events
-#rdo2aod.threads = 1 # TODO: change to 4
-#rdo2aod.concurrent_events = 1 # TODO: change to 4
-rdo2aod.forks = 4
+rdo2aod.threads = 1 # TODO: change to 4
+rdo2aod.concurrent_events = 4
rdo2aod.perfmon = False
rdo2aod.args = '--outputAODFile=AOD.pool.root --steering="doRDO_TRIG" '
if local:
diff --git a/Trigger/TriggerCommon/TriggerJobOpts/python/Lvl1SimulationConfig.py b/Trigger/TriggerCommon/TriggerJobOpts/python/Lvl1SimulationConfig.py
index beb53063b7a6694bdafa103a80888da9a0ac69f4..6794393247cb09762dded97b37abc6c358f9cb6f 100644
--- a/Trigger/TriggerCommon/TriggerJobOpts/python/Lvl1SimulationConfig.py
+++ b/Trigger/TriggerCommon/TriggerJobOpts/python/Lvl1SimulationConfig.py
@@ -139,6 +139,9 @@ def Lvl1SimulationSequence( flags = None ):
muctpi = L1Muctpi()
muctpi.LVL1ConfigSvc = svcMgr.LVL1ConfigSvc
+ # Sets up and configures the muon alignment and detector manager
+ from MuonRecExample import MuonAlignConfig # noqa: F401
+
l1MuonSim = seqAND("l1MuonSim", [
MuonRdoToMuonDigit( "MuonRdoToMuonDigit",