diff --git a/Reconstruction/VKalVrt/VrtSecInclusive/VrtSecInclusive/VrtSecInclusive.h b/Reconstruction/VKalVrt/VrtSecInclusive/VrtSecInclusive/VrtSecInclusive.h index d76bec130245bd08b013f38600c74aae5b77eb9b..8fd81e54bb03ab3ff097147cb060f2b40b40bcdf 100755 --- a/Reconstruction/VKalVrt/VrtSecInclusive/VrtSecInclusive/VrtSecInclusive.h +++ b/Reconstruction/VKalVrt/VrtSecInclusive/VrtSecInclusive/VrtSecInclusive.h @@ -466,8 +466,8 @@ namespace VKalVrtAthena { // // - StatusCode augmentDVimpactParametersToMuons(); - StatusCode augmentDVimpactParametersToElectrons(); + template<class LeptonFlavor> + StatusCode augmentDVimpactParametersToLeptons( const std::string& containerName ); }; diff --git a/Reconstruction/VKalVrt/VrtSecInclusive/VrtSecInclusive/details/Utilities.h b/Reconstruction/VKalVrt/VrtSecInclusive/VrtSecInclusive/details/Utilities.h index 71ffe475d2f484c8db4fc9d28195c76ac7d9d595..87d3c3d7acb851d725494b4228a853785c4aa949 100644 --- a/Reconstruction/VKalVrt/VrtSecInclusive/VrtSecInclusive/details/Utilities.h +++ b/Reconstruction/VKalVrt/VrtSecInclusive/VrtSecInclusive/details/Utilities.h @@ -6,6 +6,7 @@ #define _VrtSecInclusive_VrtSecInclusive_Utilities_H #include "VrtSecInclusive/IntersectionPos.h" +#include <numeric> namespace VKalVrtAthena { @@ -82,6 +83,124 @@ namespace VKalVrtAthena { delete Output; } + + //____________________________________________________________________________________________________ + template<class LeptonFlavor> + void genSequence( const LeptonFlavor*, std::vector<unsigned>& ) {} + + template<> void genSequence( const xAOD::Muon*, std::vector<unsigned>& trackTypes ); + template<> void genSequence( const xAOD::Electron* electron, std::vector<unsigned>& trackTypes ); + + //____________________________________________________________________________________________________ + template<class LeptonFlavor> + const xAOD::TrackParticle* getLeptonTrackParticle( const LeptonFlavor*, const unsigned& ) { return nullptr; } + + template<> const xAOD::TrackParticle* getLeptonTrackParticle( const xAOD::Muon* muon, const unsigned& trackType ); + template<> const xAOD::TrackParticle* getLeptonTrackParticle( const xAOD::Electron* electron, const unsigned& trackType ); + + //____________________________________________________________________________________________________ + template<class LeptonFlavor> + StatusCode VrtSecInclusive::augmentDVimpactParametersToLeptons( const std::string& containerName ) + { + + const xAOD::VertexContainer *secondaryVertexContainer( nullptr ); + ATH_CHECK( evtStore()->retrieve( secondaryVertexContainer, "VrtSecInclusive_" + m_jp.secondaryVerticesContainerName ) ); + + using LeptonContainer = DataVector<LeptonFlavor>; + + const LeptonContainer *leptonContainer( nullptr ); + ATH_CHECK( evtStore()->retrieve( leptonContainer, containerName ) ); + + using IPDecoratorType = SG::AuxElement::Decorator< std::vector< std::vector<float> > >; + static IPDecoratorType decor_d0wrtSV ( "d0_wrtSVs" ); + static IPDecoratorType decor_z0wrtSV ( "z0_wrtSVs" ); + static IPDecoratorType decor_ptwrtSV ( "pt_wrtSVs" ); + static IPDecoratorType decor_etawrtSV ( "eta_wrtSVs" ); + static IPDecoratorType decor_phiwrtSV ( "phi_wrtSVs" ); + static IPDecoratorType decor_d0errWrtSV ( "d0err_wrtSVs" ); + static IPDecoratorType decor_z0errWrtSV ( "z0err_wrtSVs" ); + + // Grouping decorators + std::vector< IPDecoratorType > decor_ipWrtSVs { decor_d0wrtSV, decor_z0wrtSV, decor_ptwrtSV, decor_etawrtSV, decor_phiwrtSV, decor_d0errWrtSV, decor_z0errWrtSV }; + enum { k_ip_d0, k_ip_z0, k_ip_pt, k_ip_eta, k_ip_phi, k_ip_d0err, k_ip_z0err }; + + static SG::AuxElement::Decorator< std::vector<ElementLink< xAOD::VertexContainer > > > decor_svLink("svLinks"); + + // Loop over leptons + for( const auto& lepton : *leptonContainer ) { + + std::vector< std::vector< std::vector<float> > > ip_wrtSVs( decor_ipWrtSVs.size() ); // triple nest of { ip parameters, tracks, DVs } + + bool linkFlag { false }; + + std::vector<unsigned> trackTypes; + genSequence<LeptonFlavor>( lepton, trackTypes ); + + // Loop over lepton types + for( auto& trackType : trackTypes ) { + + std::vector< std::vector<float> > ip_wrtSV( decor_ipWrtSVs.size() ); // nest of { tracks, DVs } + + const auto* trk = getLeptonTrackParticle<LeptonFlavor>( lepton, trackType ); + + if( !trk ) continue; + + std::map< const xAOD::Vertex*, std::vector<double> > distanceMap; + + std::vector<ElementLink< xAOD::VertexContainer > > links; + + // Loop over vertices + for( const auto& vtx : *secondaryVertexContainer ) { + + std::vector<double> impactParameters; + std::vector<double> impactParErrors; + + m_fitSvc->VKalGetImpact( trk, vtx->position(), static_cast<int>( lepton->charge() ), impactParameters, impactParErrors ); + + enum { k_d0, k_z0, k_theta, k_phi, k_qOverP }; // for the impact parameter + enum { k_d0d0, k_d0z0, k_z0z0 }; // for the par errors + + const auto& theta = impactParameters.at( k_theta ); + const auto& phi = impactParameters.at( k_phi ); + const auto p = fabs( 1.0 / impactParameters.at(k_qOverP) ); + const auto pt = fabs( p * sin( theta ) ); + const auto eta = -log( tan(theta/2.) ); + + // filling the parameters to the corresponding container + ip_wrtSV.at( k_ip_d0 ) .emplace_back( impactParameters.at(k_d0) ); + ip_wrtSV.at( k_ip_z0 ) .emplace_back( impactParameters.at(k_d0) ); + ip_wrtSV.at( k_ip_pt ) .emplace_back( pt ); + ip_wrtSV.at( k_ip_eta ) .emplace_back( eta ); + ip_wrtSV.at( k_ip_phi ) .emplace_back( phi ); + ip_wrtSV.at( k_ip_d0err ) .emplace_back( impactParErrors.at(k_d0d0) ); + ip_wrtSV.at( k_ip_z0err ) .emplace_back( impactParErrors.at(k_z0z0) ); + + if( !linkFlag ) { + + ElementLink<xAOD::VertexContainer> link_SV( *( dynamic_cast<const xAOD::VertexContainer*>( vtx->container() ) ), static_cast<size_t>( vtx->index() ) ); + links.emplace_back( link_SV ); + + } + + } // end of vertex loop + + // The linking to the vertices need to be done only once + if( !linkFlag ) { + decor_svLink ( *lepton ) = links; + linkFlag = true; + } + + for( size_t ipar = 0; ipar < ip_wrtSVs.size(); ipar++ ) ip_wrtSVs.at( ipar ).emplace_back( ip_wrtSV.at( ipar ) ); + + } // end of track type loop + + // decoration + for( size_t ipar = 0; ipar < decor_ipWrtSVs.size(); ipar++ ) decor_ipWrtSVs.at( ipar )( *lepton ) = ip_wrtSVs.at( ipar ); + + } // end of lepton container loop + + return StatusCode::SUCCESS; + } } diff --git a/Reconstruction/VKalVrt/VrtSecInclusive/share/VrtSecInclusive_DV_postInclude_ESDtoAOD.py b/Reconstruction/VKalVrt/VrtSecInclusive/share/VrtSecInclusive_DV_postInclude_ESDtoAOD.py new file mode 100644 index 0000000000000000000000000000000000000000..fc1d3d9eb10c343142236a85b61466ec318a2298 --- /dev/null +++ b/Reconstruction/VKalVrt/VrtSecInclusive/share/VrtSecInclusive_DV_postInclude_ESDtoAOD.py @@ -0,0 +1,6 @@ +StreamAOD.ItemList+=["xAOD::TrackParticleContainer#*"] +StreamAOD.ItemList+=["xAOD::TrackParticleAuxContainer#*"] +StreamAOD.ItemList+=["xAOD::VertexContainer#*"] +StreamAOD.ItemList+=["xAOD::VertexAuxContainer#*"] + + diff --git a/Reconstruction/VKalVrt/VrtSecInclusive/src/Utilities.cxx b/Reconstruction/VKalVrt/VrtSecInclusive/src/Utilities.cxx index cb44cdcb39e27584c87a90495384184c87023d8f..80725000ef81a69f610a02c0aa7e268e38b486fb 100644 --- a/Reconstruction/VKalVrt/VrtSecInclusive/src/Utilities.cxx +++ b/Reconstruction/VKalVrt/VrtSecInclusive/src/Utilities.cxx @@ -73,6 +73,34 @@ namespace VKalVrtAthena { } } + //____________________________________________________________________________________________________ + template<> + void genSequence( const xAOD::Muon*, std::vector<unsigned>& trackTypes ) { + trackTypes = { xAOD::Muon::Primary, + xAOD::Muon::InnerDetectorTrackParticle, + xAOD::Muon::MuonSpectrometerTrackParticle, + xAOD::Muon::CombinedTrackParticle, + xAOD::Muon::ExtrapolatedMuonSpectrometerTrackParticle + }; + } + + template<> + void genSequence( const xAOD::Electron* electron, std::vector<unsigned>& trackTypes ) { + for( size_t i=0; i<electron->nTrackParticles(); i++ ) trackTypes.emplace_back( i ); + } + + + //____________________________________________________________________________________________________ + template<> + const xAOD::TrackParticle* getLeptonTrackParticle( const xAOD::Muon* muon, const unsigned& trackType ) { + return muon->trackParticle( static_cast<xAOD::Muon::TrackParticleType>( trackType ) ); + } + + template<> + const xAOD::TrackParticle* getLeptonTrackParticle( const xAOD::Electron* electron, const unsigned& trackType ) { + return electron->trackParticle( trackType ); + } + //____________________________________________________________________________________________________ double VrtSecInclusive::distanceBetweenVertices( const WrkVrt& v1, const WrkVrt& v2 ) const { @@ -1556,120 +1584,7 @@ namespace VKalVrtAthena { } - //____________________________________________________________________________________________________ - StatusCode VrtSecInclusive::augmentDVimpactParametersToMuons() - { - - const xAOD::VertexContainer *secondaryVertexContainer( nullptr ); - ATH_CHECK( evtStore()->retrieve( secondaryVertexContainer, "VrtSecInclusive_" + m_jp.secondaryVerticesContainerName ) ); - - const xAOD::MuonContainer *muonContainer( nullptr ); - ATH_CHECK( evtStore()->retrieve( muonContainer, "Muons" ) ); - - static SG::AuxElement::Decorator< std::vector<float> > decor_d0wrtSV( "d0_wrtSVs" ); - static SG::AuxElement::Decorator< std::vector<float> > decor_z0wrtSV( "z0_wrtSVs" ); - static SG::AuxElement::Decorator< std::vector<ElementLink< xAOD::VertexContainer > > > decor_svLink("svLinks"); - - for( const auto& muon : *muonContainer ) { - // Loop over muons - - const auto& primaryTrackLink = muon->primaryTrackParticleLink(); - const auto* trk = *primaryTrackLink; - if( trk ) { - - std::map< const xAOD::Vertex*, std::vector<double> > distanceMap; - - std::vector<float> d0wrtSV; - std::vector<float> z0wrtSV; - std::vector<ElementLink< xAOD::VertexContainer > > links; - - for( const auto& vtx : *secondaryVertexContainer ) { - - std::vector<double> impactParameters; - std::vector<double> impactParErrors; - - m_fitSvc->VKalGetImpact( trk, vtx->position(), static_cast<int>( muon->charge() ), impactParameters, impactParErrors); - - enum { k_d0, k_z0, k_theta, k_phi, k_qOverP }; - - d0wrtSV.emplace_back( impactParameters.at(k_d0) ); - z0wrtSV.emplace_back( impactParameters.at(k_z0) ); - - ElementLink<xAOD::VertexContainer> link_SV( *( dynamic_cast<const xAOD::VertexContainer*>( vtx->container() ) ), static_cast<size_t>( vtx->index() ) ); - - links.emplace_back( link_SV ); - - } - - decor_d0wrtSV( *muon ) = d0wrtSV; - decor_z0wrtSV( *muon ) = z0wrtSV; - decor_svLink ( *muon ) = links; - - } - } - - return StatusCode::SUCCESS; - } - - - //____________________________________________________________________________________________________ - StatusCode VrtSecInclusive::augmentDVimpactParametersToElectrons() - { - - const xAOD::VertexContainer *secondaryVertexContainer( nullptr ); - ATH_CHECK( evtStore()->retrieve( secondaryVertexContainer, "VrtSecInclusive_" + m_jp.secondaryVerticesContainerName ) ); - - const xAOD::ElectronContainer *electronContainer( nullptr ); - ATH_CHECK( evtStore()->retrieve( electronContainer, "Electrons" ) ); - - static SG::AuxElement::Decorator< std::vector<float> > decor_d0wrtSV( "d0_wrtSVs" ); - static SG::AuxElement::Decorator< std::vector<float> > decor_z0wrtSV( "z0_wrtSVs" ); - static SG::AuxElement::Decorator< std::vector<ElementLink< xAOD::VertexContainer > > > decor_svLink("svLinks"); - - for( const auto& electron : *electronContainer ) { - // Loop over electrons - - if( 0 == electron->nTrackParticles() ) continue; - - // The first track is the best-matched track - const auto* trk = electron->trackParticle(0); - if( trk ) { - - std::map< const xAOD::Vertex*, std::vector<double> > distanceMap; - - std::vector<float> d0wrtSV; - std::vector<float> z0wrtSV; - std::vector<ElementLink< xAOD::VertexContainer > > links; - - for( const auto& vtx : *secondaryVertexContainer ) { - - std::vector<double> impactParameters; - std::vector<double> impactParErrors; - - m_fitSvc->VKalGetImpact( trk, vtx->position(), static_cast<int>( electron->charge() ), impactParameters, impactParErrors); - - enum { k_d0, k_z0, k_theta, k_phi, k_qOverP }; - - d0wrtSV.emplace_back( impactParameters.at(k_d0) ); - z0wrtSV.emplace_back( impactParameters.at(k_z0) ); - - ElementLink<xAOD::VertexContainer> link_SV( *( dynamic_cast<const xAOD::VertexContainer*>( vtx->container() ) ), static_cast<size_t>( vtx->index() ) ); - - links.emplace_back( link_SV ); - - } - - decor_d0wrtSV( *electron ) = d0wrtSV; - decor_z0wrtSV( *electron ) = z0wrtSV; - decor_svLink ( *electron ) = links; - - } - } - - return StatusCode::SUCCESS; - } - - + } // end of namespace VKalVrtAthena diff --git a/Reconstruction/VKalVrt/VrtSecInclusive/src/VertexingAlgs.cxx b/Reconstruction/VKalVrt/VrtSecInclusive/src/VertexingAlgs.cxx index 8051b38b30aad2832db7b725555050c385bb7361..a61d4791ed48a18805f71ab22d7694286766022c 100644 --- a/Reconstruction/VKalVrt/VrtSecInclusive/src/VertexingAlgs.cxx +++ b/Reconstruction/VKalVrt/VrtSecInclusive/src/VertexingAlgs.cxx @@ -1564,8 +1564,8 @@ namespace VKalVrtAthena { // Post process -- Additional augmentations - if( m_jp.doAugmentDVimpactParametersToMuons ) { ATH_CHECK( augmentDVimpactParametersToMuons() ); } - if( m_jp.doAugmentDVimpactParametersToElectrons ) { ATH_CHECK( augmentDVimpactParametersToElectrons() ); } + if( m_jp.doAugmentDVimpactParametersToMuons ) { ATH_CHECK( augmentDVimpactParametersToLeptons<xAOD::Muon> ( "Muons" ) ); } + if( m_jp.doAugmentDVimpactParametersToElectrons ) { ATH_CHECK( augmentDVimpactParametersToLeptons<xAOD::Electron>( "Electrons" ) ); } return StatusCode::SUCCESS; }