Commit f0c8b3fa authored by Vakhtang Tsulaia's avatar Vakhtang Tsulaia
Browse files

Merge branch 'master-vertex-tracks' into 'master'

Add track information to evaluation of vertex finding performance

See merge request atlas/athena!45574
parents 2d811b64 bb9ef5e3
......@@ -69,6 +69,9 @@ public:
double chi2() const { return m_chi2; }
int ndof() const { return m_ndof; }
protected:
void setNtracks( int nTracks ) { m_Ntracks = nTracks; }
private:
......
/* emacs: this is -*- c++ -*- */
/**
** @file TIDAVertexNew.h
**
** @author emil haines
** @date Wed 28 Jul 2021
**
** Copyright (C) 2002-2021 CERN for the benefit of the ATLAS collaboration
**/
#ifndef TIDAVERTEXNEW_H
#define TIDAVERTEXNEW_H
#include <iostream>
#include <vector>
#include "TrigInDetAnalysis/TIDAVertex.h"
#include "TrigInDetAnalysis/TrackTrigObject.h"
#include "TrigInDetAnalysis/Track.h"
namespace TIDA {
class VertexNew : public TIDA::Vertex {
public:
VertexNew( const TIDA::Vertex& v ) : Vertex(v) {}
VertexNew( const TIDA::Vertex& v,
const TrackTrigObject& tobject,
const std::vector<TIDA::Track*>* refTracks=0 );
// copy constructor
VertexNew( const TIDA::VertexNew& v,
const std::vector<TIDA::Track*>* refTracks=0,
const TrackTrigObject* tobject = 0 );
virtual ~VertexNew() {}
void selectTracks( const std::vector<TIDA::Track*>* refTracks,
const TrackTrigObject* tobject = 0 );
const std::vector<TIDA::Track*>& tracks() const { return m_tracks; }
const std::vector<unsigned long>& ids() const { return m_ids; }
private:
void selectTracksInternal( const std::vector<TIDA::Track*>* refTracks );
void addTrack(TIDA::Track* trk) { m_tracks.push_back(trk); }
std::vector<TIDA::Track*> m_tracks;
std::vector<unsigned long> m_ids;
};
}
#endif // TIDAVERTEXNEW.H
\ No newline at end of file
/**
** @file TIDAVertexNew.cxx
**
** @author emil haines
** @date Wed 28 Jul 2021
**
** Copyright (C) 2002-2021 CERN for the benefit of the ATLAS collaboration
**/
#include "TrigInDetAnalysis/TIDAVertexNew.h"
TIDA::VertexNew::VertexNew(
const TIDA::Vertex& v,
const TrackTrigObject& tobject,
const std::vector<TIDA::Track*>* refTracks )
: TIDA::Vertex(v) {
const std::vector<unsigned long>& ids = tobject.children();
m_ids = ids;
if ( refTracks ) selectTracks( refTracks );
}
TIDA::VertexNew::VertexNew(
const TIDA::VertexNew& v,
const std::vector<TIDA::Track*>* refTracks,
const TrackTrigObject* tobject )
: TIDA::Vertex( v ),
m_tracks( v.tracks() ),
m_ids( v.ids() ) {
if ( tobject ) {
const std::vector<unsigned long>& ids = tobject->children();
m_ids = ids;
}
if ( refTracks ) selectTracks( refTracks );
}
void TIDA::VertexNew::selectTracks(
const std::vector<TIDA::Track*>* refTracks,
const TrackTrigObject* tobject ) {
if ( tobject ) {
const std::vector<unsigned long>& ids = tobject->children();
m_ids = ids;
}
if ( m_ids.size() == 0 ) {
return;
}
m_tracks.clear();
selectTracksInternal( refTracks );
}
void TIDA::VertexNew::selectTracksInternal(
const std::vector<TIDA::Track*>* refTracks ) {
int trackcount = 0;
for ( size_t iid = 0; iid < m_ids.size(); iid++ ) {
for (unsigned itr = 0; itr < refTracks->size(); itr++) {
TIDA::Track* tr = (*refTracks)[itr];
if ( tr->id() == m_ids[iid] ) {
addTrack( tr );
trackcount++;
break;
}
}
}
setNtracks( trackcount );
}
\ No newline at end of file
......@@ -474,6 +474,10 @@ void AnalysisConfigMT_Ntuple::loop() {
std::vector<TIDA::Vertex> vertices;
// store offline vertex track ids along with vertices
std::vector<TrackTrigObject> offVertexTracks;
// std::vector<TIDA::Vertex> vertices;
m_provider->msg(MSG::VERBOSE) << "fetching AOD Primary vertex container" << endmsg;
......@@ -496,18 +500,35 @@ void AnalysisConfigMT_Ntuple::loop() {
// std::cout << "SUTT xAOD::Vertex::type() " << (*vtxitr)->type() << "\tvtxtype " << (*vtxitr)->vertexType() << "\tntrax " << (*vtxitr)->nTrackParticles() << std::endl;
if ( (*vtxitr)->nTrackParticles()>0 && (*vtxitr)->vertexType()!=0 ) {
vertices.push_back( TIDA::Vertex( (*vtxitr)->x(),
(*vtxitr)->y(),
(*vtxitr)->z(),
/// variances
(*vtxitr)->covariancePosition()(Trk::x,Trk::x),
(*vtxitr)->covariancePosition()(Trk::y,Trk::y),
(*vtxitr)->covariancePosition()(Trk::z,Trk::z),
(*vtxitr)->nTrackParticles(),
/// quality
(*vtxitr)->chiSquared(),
(*vtxitr)->numberDoF() ) );
}
vertices.push_back( TIDA::Vertex( (*vtxitr)->x(),
(*vtxitr)->y(),
(*vtxitr)->z(),
/// variances
(*vtxitr)->covariancePosition()(Trk::x,Trk::x),
(*vtxitr)->covariancePosition()(Trk::y,Trk::y),
(*vtxitr)->covariancePosition()(Trk::z,Trk::z),
(*vtxitr)->nTrackParticles(),
/// quality
(*vtxitr)->chiSquared(),
(*vtxitr)->numberDoF() ) );
// get tracks associated to vertex
const std::vector< ElementLink< xAOD::TrackParticleContainer > >& tracks = (*vtxitr)->trackParticleLinks();
// convert from xAOD into TIDA::Track
TrigTrackSelector selector( &filter_etaPT ); // not sure about the filter, copied line 147
selector.selectTracks( tracks );
const std::vector<TIDA::Track*>& tidatracks = selector.tracks();
// Store ids of tracks belonging to vertex in TrackTrigObject
TrackTrigObject vertexTracks = TrackTrigObject();
for ( auto trkitr = tidatracks.begin(); trkitr != tidatracks.end(); ++trkitr ) {
vertexTracks.addChild( (*trkitr)->id() );
}
offVertexTracks.push_back( vertexTracks );
}
}
}
......@@ -552,6 +573,7 @@ void AnalysisConfigMT_Ntuple::loop() {
m_event->addChain( "Vertex" );
m_event->back().addRoi(TIDARoiDescriptor(true));
m_event->back().back().addVertices( vertices );
m_event->back().back().addObjects( offVertexTracks );
}
......@@ -1049,6 +1071,10 @@ void AnalysisConfigMT_Ntuple::loop() {
/// fetch vertices if available ...
std::vector<TIDA::Vertex> tidavertices;
// store trigger vertex track ids along with vertices
std::vector<TrackTrigObject> tidaVertexTracks;
if ( vtx_name!="" ) {
......@@ -1069,7 +1095,7 @@ void AnalysisConfigMT_Ntuple::loop() {
xAOD::VertexContainer::const_iterator vtxitr = vtx_itrpair.first;
for ( ; vtxitr!=vtx_itrpair.second ; vtxitr++ ) {
for ( ; vtxitr!=vtx_itrpair.second ; vtxitr++ ) {
/// leave this code commented so that we have a record of the change - as soon as we can
/// fix the missing track multiplicity from the vertex this will need to go back
......@@ -1086,12 +1112,29 @@ void AnalysisConfigMT_Ntuple::loop() {
/// quality
(*vtxitr)->chiSquared(),
(*vtxitr)->numberDoF() ) );
// get tracks associated to vertex
const std::vector< ElementLink< xAOD::TrackParticleContainer > >& tracks = (*vtxitr)->trackParticleLinks();
// covert from xAOD into TIDA::Track
TrigTrackSelector selector( &filter ); // not sure about the filter, copied line 148
selector.selectTracks( tracks );
const std::vector<TIDA::Track*>& tidatracks = selector.tracks();
// Store ids of tracks belonging to vertex in TrackTrigObject
TrackTrigObject vertexTracks = TrackTrigObject();
for ( auto trkitr = tidatracks.begin(); trkitr != tidatracks.end(); ++trkitr ) {
vertexTracks.addChild( (*trkitr)->id() );
}
tidaVertexTracks.push_back( vertexTracks );
}
}
}
}
}
}
}
#if 0
//// not yet ready to get the jet yet - this can come once everything else is working
......@@ -1118,10 +1161,10 @@ void AnalysisConfigMT_Ntuple::loop() {
}
chain.addRoi( *roi_tmp );
chain.back().addTracks(testTracks);
chain.back().addVertices(tidavertices);
chain.back().addObjects(tidaVertexTracks);
#if 0
/// jets can't be added yet
......
......@@ -11,6 +11,7 @@
#include "ConfVtxAnalysis.h"
#include "TrigInDetAnalysisUtils/VertexMatcher.h"
#include "TrigInDetAnalysisUtils/VertexNewMatcher.h"
#include "TrigInDetAnalysis/TIDAEvent.h"
......@@ -164,119 +165,152 @@ void ConfVtxAnalysis::execute( const std::vector<TIDA::Vertex*>& vtx0,
std::cout << "\ttest: " << vtx1 << std::endl;
#endif
VertexMatcher m( "vtx_matcher", 3 );
execute_internal<TIDA::Vertex, VertexMatcher>(vtx0, vtx1, m, tevt);
VertexMatcher m( "vtx_matcher", 3 );
}
m.match( vtx0, vtx1 );
hnvtx->Fill( vtx0.size() );
hnvtx_rec->Fill( vtx1.size() );
rnvtxrec_nvtx->Fill( vtx0.size(), vtx1.size() );
void ConfVtxAnalysis::execute( const std::vector<TIDA::VertexNew*>& vtx0,
const std::vector<TIDA::VertexNew*>& vtx1,
const TIDA::Event* tevt ) {
// std::cout << "gevent " << gevent << std::endl;
if ( !m_initialised ) return;
/// pass in a parameter now, rather than using a global
// double mu = gevent->mu();
// double lb = gevent->lumi_block();
double mu = tevt->mu();
double lb = tevt->lumi_block();
#if 0
std::cout << "ConfVtxAnalysis::execute() " << name()
<< "\tvtx0.size() " << vtx0.size()
<< "\tvtx1.size() " << vtx1.size()
<< std::endl;
for ( unsigned i=0 ; i<vtx0.size() ; i++ ) {
std::cout << "\tref: " << vtx0 << std::endl;
std::cout << "\ttest: " << vtx1 << std::endl;
#endif
/// reject vertices with no tracks in the Roi ...
if ( vtx0[i]->Ntracks() == 0 ) continue;
// use new matcher if tracks included
if ( vtx0[0]->tracks().size() > 0 ) {
VertexNewMatcher m( "vtx_matcher", 0.5 );
execute_internal<TIDA::VertexNew, VertexNewMatcher>(vtx0, vtx1, m, tevt);
}
else {
VertexMatcher m("vtx_matcher", 3 );
// m.match() requires std::vector<TIDA::Vertex*>
std::vector<TIDA::Vertex*> vtx0_(vtx0.begin(), vtx0.end());
std::vector<TIDA::Vertex*> vtx1_(vtx1.begin(), vtx1.end());
execute_internal<TIDA::Vertex, VertexMatcher>(vtx0_, vtx1_, m, tevt);
}
}
hzed->Fill( vtx0[i]->z() );
hntrax->Fill( vtx0[i]->Ntracks() );
template<typename Vertex, typename Matcher>
void ConfVtxAnalysis::execute_internal( const std::vector<Vertex*>& vtx0,
const std::vector<Vertex*>& vtx1,
Matcher& m,
const TIDA::Event* tevt ) {
m.match( vtx0, vtx1 );
hnvtx->Fill( vtx0.size() );
hnvtx_rec->Fill( vtx1.size() );
rnvtxrec_nvtx->Fill( vtx0.size(), vtx1.size() );
hlb->Fill( lb );
hmu->Fill( mu );
const TIDA::Vertex* mv = m.matched( vtx0[i] );
// std::cout << "\tvtx match: " << i << " " << mv << std::endl;
// keep alternate functionality commented
// std::cout << "gevent " << gevent << std::endl;
if ( mv ) {
// std::cout << "\ttest z " << mv->z() << " : delta z " << (mv->z()-vtx0[i]->z()) << std::endl;
/// ah ha ! can fill some silly old histograms here
/// ...
/// pass in a parameter now, rather than using a global
// double mu = gevent->mu();
// double lb = gevent->lumi_block();
double mu = tevt->mu();
double lb = tevt->lumi_block();
hzed_rec->Fill( mv->z() );
hntrax_rec->Fill( mv->Ntracks() );
for ( unsigned i=0 ; i<vtx0.size() ; i++ ) {
hzed_res->Fill( mv->z() - vtx0[i]->z() );
/// reject vertices with no tracks in the Roi ...
if ( vtx0[i]->Ntracks() == 0 ) continue;
rdz_vs_zed->Fill( vtx0[i]->z(), mv->z() - vtx0[i]->z() );
rdz_vs_ntrax->Fill( vtx0[i]->Ntracks(), mv->z() - vtx0[i]->z() );
rdz_vs_nvtx->Fill( vtx0.size(), mv->z() - vtx0[i]->z() ); /// this isn't really legitimate
rdz_vs_mu->Fill( mu, mv->z() - vtx0[i]->z() ); /// this isn't really legitimate
eff_zed->Fill( vtx0[i]->z() );
eff_ntrax->Fill( vtx0[i]->Ntracks() );
eff_nvtx->Fill( vtx0.size() );
hzed->Fill( vtx0[i]->z() );
hntrax->Fill( vtx0[i]->Ntracks() );
eff_mu->Fill( mu );
eff_lb->Fill( lb );
// std::cout << "found vtx ref vertex size " << vtx0.size() << "\tonline " << vtx1.size() << std::endl;
// std::cout << "\tref: " << *vtx0[i] << std::endl;
// for ( unsigned iv=0 ; iv<vtx1.size() ; iv++ ) if ( vtx1[iv] ) std::cout << "\t" << iv << " : " << *vtx1[iv] << std::endl;
/// what about beam tilts etc? where are these defined with respect to ?
rdz_vs_lb->Fill( lb, mv->z() - vtx0[i]->z() );
rdx_vs_lb->Fill( lb, mv->x() - vtx0[i]->x() );
rdy_vs_lb->Fill( lb, mv->y() - vtx0[i]->y() );
hlb->Fill( lb );
hmu->Fill( mu );
const Vertex* mv = m.matched( vtx0[i] );
// std::cout << "\tvtx match: " << i << " " << mv << std::endl;
if ( mv ) {
// std::cout << "\ttest z " << mv->z() << " : delta z " << (mv->z()-vtx0[i]->z()) << std::endl;
/// ah ha ! can fill some silly old histograms here
/// ...
hzed_rec->Fill( mv->z() );
hntrax_rec->Fill( mv->Ntracks() );
hzed_res->Fill( mv->z() - vtx0[i]->z() );
rdz_vs_zed->Fill( vtx0[i]->z(), mv->z() - vtx0[i]->z() );
rdz_vs_ntrax->Fill( vtx0[i]->Ntracks(), mv->z() - vtx0[i]->z() );
rdz_vs_nvtx->Fill( vtx0.size(), mv->z() - vtx0[i]->z() ); /// this isn't really legitimate
rdz_vs_mu->Fill( mu, mv->z() - vtx0[i]->z() ); /// this isn't really legitimate
eff_zed->Fill( vtx0[i]->z() );
eff_ntrax->Fill( vtx0[i]->Ntracks() );
eff_nvtx->Fill( vtx0.size() );
eff_mu->Fill( mu );
eff_lb->Fill( lb );
// std::cout << "found vtx ref vertex size " << vtx0.size() << "\tonline " << vtx1.size() << std::endl;
// std::cout << "\tref: " << *vtx0[i] << std::endl;
// for ( unsigned iv=0 ; iv<vtx1.size() ; iv++ ) if ( vtx1[iv] ) std::cout << "\t" << iv << " : " << *vtx1[iv] << std::endl;
/// what about beam tilts etc? where are these defined with respect to ?
rdz_vs_lb->Fill( lb, mv->z() - vtx0[i]->z() );
rdx_vs_lb->Fill( lb, mv->x() - vtx0[i]->x() );
rdy_vs_lb->Fill( lb, mv->y() - vtx0[i]->y() );
}
else {
// std::cout << "\t" << "------" << std::endl;
// std::cout << "missing vtx ref vertex size " << vtx0.size() << "\tonline " << vtx1.size() << std::endl;
// std::cout << "\tref: " << *vtx0[i] << std::endl;
// for ( unsigned iv=0 ; iv<vtx1.size() ; iv++ ) if ( vtx1[iv] ) std::cout << "\t" << iv << " : " << *vtx1[iv] << std::endl;
} else {
// std::cout << "\t" << "------" << std::endl;
// std::cout << "missing vtx ref vertex size " << vtx0.size() << "\tonline " << vtx1.size() << std::endl;
// std::cout << "\tref: " << *vtx0[i] << std::endl;
// for ( unsigned iv=0 ; iv<vtx1.size() ; iv++ ) if ( vtx1[iv] ) std::cout << "\t" << iv << " : " << *vtx1[iv] << std::endl;
#if 0
static int vtxcount=0;
if ( vtxcount<100 ) {
std::cout << "ConfVtxAnalysis::execute() " << name() << "\tnomatch\n"
<< "\tvtx0.size() " << vtx0.size()
<< "\tvtx1.size() " << vtx1.size()
<< std::endl;
std::cout << "\tref: " << vtx0 << std::endl;
std::cout << "\ttest: " << vtx1 << std::endl;
vtxcount++;
// std::cout << *gevent << std::endl;
}
static int vtxcount=0;
if ( vtxcount<100 ) {
std::cout << "ConfVtxAnalysis::execute() " << name() << "\tnomatch\n"
<< "\tvtx0.size() " << vtx0.size()
<< "\tvtx1.size() " << vtx1.size()
<< std::endl;
std::cout << "\tref: " << vtx0 << std::endl;
std::cout << "\ttest: " << vtx1 << std::endl;
vtxcount++;
// std::cout << *gevent << std::endl;
}
#endif
eff_zed->FillDenom( vtx0[i]->z() );
eff_ntrax->FillDenom( vtx0[i]->Ntracks() );
eff_nvtx->FillDenom( vtx0.size() );
eff_zed->FillDenom( vtx0[i]->z() );
eff_ntrax->FillDenom( vtx0[i]->Ntracks() );
eff_nvtx->FillDenom( vtx0.size() );
eff_mu->FillDenom( mu );
eff_lb->FillDenom( lb );
eff_mu->FillDenom( mu );
eff_lb->FillDenom( lb );
}
}
}
}
void ConfVtxAnalysis::finalise() {
if ( !m_initialised ) return;
......
......@@ -15,6 +15,7 @@
#include <iostream>
#include "TrigInDetAnalysis/VertexAnalysis.h"
#include "TrigInDetAnalysis/TIDAVertexNew.h"
#include "TrigInDetAnalysis/TIDDirectory.h"
#include "TrigInDetAnalysis/Efficiency.h"
......@@ -31,13 +32,23 @@ public:
void initialise();
void execute(const std::vector<TIDA::Vertex*>& vtx0,
const std::vector<TIDA::Vertex*>& vtx1,
const TIDA::Event* tevt=0 );
const std::vector<TIDA::Vertex*>& vtx1,
const TIDA::Event* tevt=0 );
void execute(const std::vector<TIDA::VertexNew*>& vtx0,
const std::vector<TIDA::VertexNew*>& vtx1,
const TIDA::Event* tevt=0 );
void finalise();
private:
template<typename Vertex, typename Matcher>
void execute_internal(const std::vector<Vertex*>& vtx0,
const std::vector<Vertex*>& vtx1,
Matcher& m,
const TIDA::Event* tevt=0);
bool m_initialised;
bool m_finalised;
......
......@@ -31,6 +31,7 @@
#include "TrigInDetAnalysis/Track.h"
#include "TrigInDetAnalysis/TIDAEvent.h"
#include "TrigInDetAnalysis/TrackSelector.h"
#include "TrigInDetAnalysis/TIDAVertexNew.h"
#include "TrigInDetAnalysisUtils/Associator_BestMatch.h"
#include "TrigInDetAnalysisUtils/Filters.h"
......@@ -575,7 +576,6 @@ int main(int argc, char** argv)
//bool printflag = false; // JK removed (unused)
bool rotate_testtracks = false;
if ( inputdata.isTagDefined("RotateTestTracks") ) rotate_testtracks = ( inputdata.GetValue("RotateTestTracks") ? true : false );
......@@ -860,7 +860,9 @@ int main(int argc, char** argv)