Commit 18e6ae32 authored by Oana Vickey Boeriu's avatar Oana Vickey Boeriu
Browse files

Merge branch '21.2_MultiSecVertex_update' into '21.2'

Update of MultiSecVertexTool

See merge request !43277
parents c04a1565 a7ff5b2a
/*
Copyright (C) 2002-2021 CERN for the benefit of the ATLAS collaboration
*/
///////////////////////////////////////////////////////////////////
// AllBVertexAlg.h, (c) ATLAS Detector software
// author: Vadim Kostyukhin (vadim.kostyukhin@cern.ch)
///////////////////////////////////////////////////////////////////
/**
* Reconstruct all B vertices in events, including soft ones
*/
#ifndef VRTSECINCLUSIVE_AllBVertexAlg_H
#define VRTSECINCLUSIVE_AllBVertexAlg_H
#include <string>
#include <vector>
#include "AthenaBaseComps/AthAlgorithm.h"
#include "GaudiKernel/ToolHandle.h"
#include "VrtSecInclusive/IVrtInclusive.h"
namespace Rec {
class AllBVertexAlg : public AthAlgorithm {
public:
AllBVertexAlg( const std::string& name, ISvcLocator* pSvcLocator );
StatusCode initialize() override;
StatusCode finalize() override;
StatusCode execute() override;
private:
std::string m_bvertexContainerName;
std::string m_trContainerName;
std::string m_pvContainerName;
ToolHandle < Rec::IVrtInclusive > m_bvertextool;
double vrtVrtDist(const xAOD::Vertex & primVrt, const xAOD::Vertex & secVrt) const;
};
}
#endif // VRTSECINCLUSIVE_AllBVertexAlg_H
......@@ -146,6 +146,7 @@ namespace Rec {
long int m_cutSctHits{};
long int m_cutPixelHits{};
long int m_cutSiHits{};
long int m_cutTRTHits{};
long int m_cutBLayHits{};
long int m_cutSharedHits{};
double m_cutPt{};
......
# Copyright (C) 2002-2020 CERN for the benefit of the ATLAS collaboration
# Copyright (C) 2002-2021 CERN for the benefit of the ATLAS collaboration
# Author: Vadim Kostyukhin vadim.kostyukhin@cern.ch
from GaudiKernel.GaudiHandles import *
from GaudiKernel.Proxy.Configurable import *
......@@ -69,15 +69,16 @@ class InclusiveBFinderTool( Rec__MultiSecVertexTool ):
Rec__MultiSecVertexTool.__init__( self, name = name,
VertexFitterTool = SVertexFitterTool,
CutBLayHits = 0,
CutPixelHits = 2,
CutPixelHits = 1,
CutSiHits = 8,
CutTRTHits = 10,
useVertexCleaning = True,
MultiWithOneTrkVrt = True,
removeTrkMatSignif = -1., # No additional material rejection
AntiPileupSigRCut = 2.0,
TrkSigCut = 2.0,
SelVrtSigCut = 3.0,
v2tIniBDTCut =-0.5,
v2tIniBDTCut =-0.7,
v2tFinBDTCut = 0.0,
CutPt = 500
)
......@@ -183,7 +184,29 @@ class DVFinderTool( Rec__MultiSecVertexTool ):
SelVrtSigCut = 8.0,
v2tIniBDTCut =-0.5,
v2tFinBDTCut = 0.0,
VrtMassLimit = 50000.,
Vrt2TrMassLimit = 50000.,
VrtMassLimit = 1000000.,
Vrt2TrMassLimit = 100000.,
CutPt = 1000.
)
##########################################################################################################
# define the Test algorithms
from VrtSecInclusiveConf import Rec__AllBVertexAlg
class AllBVertexFinderAlg( Rec__AllBVertexAlg ):
def __init__(self, name = 'AllBVertexFinderAlg' ):
from AthenaCommon.AppMgr import ToolSvc
mlog = logging.getLogger( 'AllBVertexFinderAlg::__init__ ' )
mlog.info("entering")
BFinderTool = InclusiveBFinderTool()
ToolSvc += BFinderTool
#----------------------
# All B-hadron vertex finder itself
#
Rec__AllBVertexAlg.__init__( self, name = name,
BVertexTool = BFinderTool
)
/*
Copyright (C) 2002-2021 CERN for the benefit of the ATLAS collaboration
*/
///////////////////////////////////////////////////////////////////
// AllBVertexAlg.cxx, (c) ATLAS Detector software
// Author: Vadim Kostyukhin (vadim.kostyukhin@cern.ch)
///////////////////////////////////////////////////////////////////
#include "VrtSecInclusive/AllBVertexAlg.h"
#include "xAODTracking/TrackParticleContainer.h"
#include "xAODTracking/VertexContainer.h"
#include "xAODTracking/VertexAuxContainer.h"
#include "VrtSecInclusive/MultiSecVertexTool.h"
#include "TLorentzVector.h"
#include "CxxUtils/sincos.h"
namespace Rec {
static const SG::AuxElement::Decorator<float> bvrtM("bvrtM");
static const SG::AuxElement::Decorator<float> bvrtPt("bvrtPt");
static const SG::AuxElement::Decorator<float> bvrtPhi("bvrtPhi");
static const SG::AuxElement::Decorator<float> bvrtEta("bvrtEta");
static const SG::AuxElement::Decorator<float> bvrtMinTPt("bvrtMinTPt");
static const SG::AuxElement::Decorator<float> bvrtSig3D("bvrtSig3D");
AllBVertexAlg::AllBVertexAlg(const std::string& name, ISvcLocator* pSvcLocator) :
AthAlgorithm( name, pSvcLocator ),
m_bvertexContainerName("AllBVertices"),
m_trContainerName("InDetTrackParticles"),
m_pvContainerName("PrimaryVertices"),
m_bvertextool("Rec::MultiSecVertexTool/InclusivetBVertexTool",this)
{
declareProperty("BVertexContainerName", m_bvertexContainerName);
declareProperty("TrackContainerName", m_trContainerName);
declareProperty("PVContainerName",m_pvContainerName);
declareProperty("BVertexTool",m_bvertextool);
}
StatusCode AllBVertexAlg::initialize()
{
if (m_bvertexContainerName=="") {
ATH_MSG_ERROR("No SG name provided for the output of AllBVertexAlg!");
return StatusCode::FAILURE;
}
CHECK( m_bvertextool.retrieve() );
return StatusCode::SUCCESS;
}
StatusCode AllBVertexAlg::finalize()
{
return StatusCode::SUCCESS;
}
StatusCode AllBVertexAlg::execute()
{
// check that container we want to write in SG does not yet exist
if (evtStore()->contains<xAOD::VertexContainer>(m_bvertexContainerName)) {
ATH_MSG_ERROR("Algorithm is attempting to write a StoreGate key " << m_bvertexContainerName
<< " which already exists. Please use a different key");
return StatusCode::FAILURE;
}
const xAOD::VertexContainer* pv_cont = nullptr;
ATH_CHECK( evtStore()->retrieve( pv_cont, m_pvContainerName ) );
std::unique_ptr<xAOD::VertexContainer> bVertexContainer = std::make_unique<xAOD::VertexContainer>();
std::unique_ptr<xAOD::VertexAuxContainer> bVertexAuxContainer = std::make_unique<xAOD::VertexAuxContainer>();
bVertexContainer->setStore(bVertexAuxContainer.get());
const xAOD::Vertex* pv = nullptr;
for ( auto v : *pv_cont ) {
if (v->vertexType()==xAOD::VxType::PriVtx) { pv = v; break; }
}
if (!pv) {
ATH_MSG_WARNING("Primary vertex not found");
CHECK(evtStore()->record( bVertexContainer.release(), m_bvertexContainerName));
CHECK(evtStore()->record( bVertexAuxContainer.release(), m_bvertexContainerName+"Aux."));
return StatusCode::SUCCESS;
}
ATH_MSG_DEBUG("Found PV");
// retrieve particle collections
const xAOD::TrackParticleContainer* tp_cont = nullptr;
ATH_CHECK( evtStore()->retrieve( tp_cont, m_trContainerName ) );
std::vector<const xAOD::TrackParticle*> trkparticles(0);
for(auto tp : (*tp_cont)) trkparticles.push_back(tp);
std::unique_ptr<Trk::VxSecVertexInfo> foundVrts = m_bvertextool->findAllVertices(trkparticles,*pv);
if(foundVrts && foundVrts->vertices().size()){
const std::vector<xAOD::Vertex*> vtmp=foundVrts->vertices();
for(auto & iv : vtmp) {
bVertexContainer->push_back(iv);
std::vector< Trk::VxTrackAtVertex > & vtrk = iv->vxTrackAtVertex();
TLorentzVector VSUM(0.,0.,0.,0.);
TLorentzVector tmp;
double minTPt=1.e10;
for(auto & it : vtrk){
const Trk::Perigee* mPer = dynamic_cast<const Trk::Perigee*>(it.perigeeAtVertex());
CxxUtils::sincos phi(mPer->parameters()[Trk::phi]);
CxxUtils::sincos theta(mPer->parameters()[Trk::theta]);
double absP = 1./std::abs(mPer->parameters()[Trk::qOverP]);
tmp.SetXYZM( phi.cs*theta.sn*absP, phi.sn*theta.sn*absP, theta.cs*absP, Trk::ParticleMasses().mass[Trk::pion]);
minTPt=std::min(minTPt,tmp.Pt());
VSUM+=tmp;
}
bvrtM(*iv) =VSUM.M();
bvrtPt(*iv) =VSUM.Pt();
bvrtEta(*iv)=VSUM.Eta();
bvrtPhi(*iv)=VSUM.Phi();
bvrtMinTPt(*iv)=minTPt;
bvrtSig3D(*iv)=vrtVrtDist(*pv,*iv);
}
}
ATH_MSG_DEBUG("BVertex container size: " << bVertexContainer->size());
CHECK(evtStore()->record( bVertexContainer.release(), m_bvertexContainerName));
CHECK(evtStore()->record( bVertexAuxContainer.release(), m_bvertexContainerName+"Aux."));
return StatusCode::SUCCESS;
}
double AllBVertexAlg::vrtVrtDist(const xAOD::Vertex & primVrt, const xAOD::Vertex & secVrt)
const
{
double distx = secVrt.x()- primVrt.x();
double disty = secVrt.y()- primVrt.y();
double distz = secVrt.z()- primVrt.z();
AmgSymMatrix(3) covMtx=primVrt.covariancePosition()+secVrt.covariancePosition();
AmgSymMatrix(3) wgtMtx = covMtx.inverse();
if(wgtMtx(0,0)<=0. || wgtMtx(1,1)<=0. || wgtMtx(2,2)<=0.) return 1.e10; //rudimentary corrrectness check
double signif = distx*wgtMtx(0,0)*distx
+disty*wgtMtx(1,1)*disty
+distz*wgtMtx(2,2)*distz
+2.*distx*wgtMtx(0,1)*disty
+2.*distx*wgtMtx(0,2)*distz
+2.*disty*wgtMtx(1,2)*distz;
signif=std::sqrt(std::abs(signif));
if( signif!=signif ) signif = 0.;
return signif;
}
}
......@@ -42,6 +42,7 @@ MultiSecVertexTool::MultiSecVertexTool(const std::string& type,
m_cutSctHits(4),
m_cutPixelHits(2),
m_cutSiHits(8),
m_cutTRTHits(10),
m_cutBLayHits(0),
m_cutSharedHits(1),
m_cutPt(500.),
......@@ -79,6 +80,7 @@ MultiSecVertexTool::MultiSecVertexTool(const std::string& type,
declareProperty("CutSctHits", m_cutSctHits , "Remove track is it has less SCT hits");
declareProperty("CutPixelHits", m_cutPixelHits, "Remove track is it has less Pixel hits");
declareProperty("CutSiHits", m_cutSiHits, "Remove track is it has less Pixel+SCT hits");
declareProperty("CutTRTHits", m_cutTRTHits, "Remove track is it has less TRT hits");
declareProperty("CutBLayHits", m_cutBLayHits, "Remove track is it has less B-layer hits");
declareProperty("CutSharedHits", m_cutSharedHits,"Reject final 2tr vertices if tracks have shared hits");
......@@ -1074,10 +1076,11 @@ MultiSecVertexTool::MultiSecVertexTool(const std::string& type,
const Trk::Perigee mPer=(*i_ntrk)->perigeeParameters() ;
double CovTrkMtx00 = (*(mPer.covariance()))(0,0);
uint8_t PixelHits,SctHits,BLayHits;
uint8_t PixelHits,SctHits,BLayHits,TRTHits;
if( !((*i_ntrk)->summaryValue(PixelHits,xAOD::numberOfPixelHits)) ) continue; // Track is
if( !((*i_ntrk)->summaryValue( SctHits,xAOD::numberOfSCTHits)) ) continue; // definitely
if( SctHits<3 ) continue; // bad
if( !((*i_ntrk)->summaryValue( TRTHits,xAOD::numberOfTRTHits)) ) TRTHits=0;
if( !((*i_ntrk)->summaryValue(BLayHits,xAOD::numberOfInnermostPixelLayerHits))) BLayHits=0;
if(m_fillHist){ m_hb_trkSelect->Fill( 2., m_w_1);}
......@@ -1109,6 +1112,7 @@ MultiSecVertexTool::MultiSecVertexTool(const std::string& type,
if(SctHits < m_cutSctHits) continue;
if((PixelHits+SctHits) < m_cutSiHits) continue;
if(BLayHits < m_cutBLayHits) continue;
if(std::abs((*i_ntrk)->eta())<1.9 && TRTHits < m_cutTRTHits) continue; //TRT hits must be present inside TRT
if(m_fillHist){ m_hb_trkSelect->Fill( 5., m_w_1);}
orderedTrk.emplace(signifBeam,*i_ntrk);
selectedTracks.push_back(*i_ntrk);
......
......@@ -3,17 +3,20 @@
#include "VrtSecInclusive/IPAugmentor.h"
#include "VrtSecInclusive/TrackRandomizer.h"
#include "VrtSecInclusive/SoftBtagTrackSelector.h"
#include "VrtSecInclusive/AllBVertexAlg.h"
#include "GaudiKernel/DeclareFactoryEntries.h"
DECLARE_NAMESPACE_ALGORITHM_FACTORY( VKalVrtAthena, VrtSecInclusive )
DECLARE_NAMESPACE_ALGORITHM_FACTORY( VKalVrtAthena, IPAugmentor )
DECLARE_NAMESPACE_ALGORITHM_FACTORY( VKalVrtAthena, TrackRandomizer )
DECLARE_NAMESPACE_ALGORITHM_FACTORY( VKalVrtAthena, SoftBtagTrackSelector )
DECLARE_NAMESPACE_ALGORITHM_FACTORY( Rec, AllBVertexAlg )
DECLARE_NAMESPACE_TOOL_FACTORY( Rec, MultiSecVertexTool )
DECLARE_FACTORY_ENTRIES(VrtSecInclusive) {
DECLARE_NAMESPACE_ALGORITHM( VKalVrtAthena, VrtSecInclusive )
DECLARE_NAMESPACE_ALGORITHM( VKalVrtAthena, IPAugmentor )
DECLARE_NAMESPACE_ALGORITHM( VKalVrtAthena, TrackRandomizer )
DECLARE_NAMESPACE_ALGORITHM( VKalVrtAthena, SoftBtagTrackSelector )
DECLARE_NAMESPACE_ALGORITHM( Rec, AllBVertexAlg )
DECLARE_NAMESPACE_ALGTOOL( Rec, MultiSecVertexTool )
}
//Notes:
......
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment