From 845200c9900021360a29a3c11f9b3a9f0cf6f02f Mon Sep 17 00:00:00 2001
From: Dave Casper <dcasper@uci.edu>
Date: Sun, 4 Mar 2018 20:51:14 -0800
Subject: [PATCH] Add truth vertex analysis to test algorithm; fix bug in
density calculation; tested successful on full-truth MC.
Former-commit-id: 1d841f3b9826d4159dc87e0d78ddc9a8a2ecb096
---
.../GaussianTrackDensity.h | 2 +-
.../GaussianDensityTestAlg_jobOptions.py | 7 +-
.../src/GaussianDensityTestAlg.cxx | 160 +++++++++++++++++-
.../src/GaussianDensityTestAlg.h | 27 +++
.../src/GaussianTrackDensity.cxx | 4 +-
5 files changed, 192 insertions(+), 8 deletions(-)
diff --git a/Tracking/TrkVertexFitter/TrkVertexSeedFinderUtils/TrkVertexSeedFinderUtils/GaussianTrackDensity.h b/Tracking/TrkVertexFitter/TrkVertexSeedFinderUtils/TrkVertexSeedFinderUtils/GaussianTrackDensity.h
index f7b451e83d1..1bb7c927b7f 100644
--- a/Tracking/TrkVertexFitter/TrkVertexSeedFinderUtils/TrkVertexSeedFinderUtils/GaussianTrackDensity.h
+++ b/Tracking/TrkVertexFitter/TrkVertexSeedFinderUtils/TrkVertexSeedFinderUtils/GaussianTrackDensity.h
@@ -154,7 +154,7 @@ namespace Trk
"Maximum radial impact parameter significance to use track" };
Gaudi::Property<double> m_z0MaxSignificance { this,
"MaxZ0Significance",
- 5.0,
+ 12.0,
"Maximum longitudinal impact parameter significance to include track in weight"};
};
diff --git a/Tracking/TrkVertexFitter/TrkVertexSeedFinderUtils/share/GaussianDensityTestAlg_jobOptions.py b/Tracking/TrkVertexFitter/TrkVertexSeedFinderUtils/share/GaussianDensityTestAlg_jobOptions.py
index 60b52d66795..6784f0b8c1e 100644
--- a/Tracking/TrkVertexFitter/TrkVertexSeedFinderUtils/share/GaussianDensityTestAlg_jobOptions.py
+++ b/Tracking/TrkVertexFitter/TrkVertexSeedFinderUtils/share/GaussianDensityTestAlg_jobOptions.py
@@ -98,11 +98,16 @@ InDetTrackSelectorTool = InDet__InDetTrackSelectionTool(name = "InDetDetailedTra
Extrapolator = InDetExtrapolator)
ToolSvc += InDetTrackSelectorTool
+from TrkVertexFitterUtils.TrkVertexFitterUtilsConf import Trk__TrackToVertexIPEstimator
+IPTool = Trk__TrackToVertexIPEstimator( name = "IPEstimator",
+ Extrapolator = InDetExtrapolator )
+
from TrkVertexSeedFinderUtils.TrkVertexSeedFinderUtilsConf import Trk__GaussianDensityTestAlg
myAlg = Trk__GaussianDensityTestAlg( name = "TestAlg",
Estimator = GaussianDensity,
- TrackSelector = InDetTrackSelectorTool)
+ TrackSelector = InDetTrackSelectorTool,
+ IPEstimator = IPTool)
myAlg.OutputLevel = DEBUG
algSeq += myAlg
diff --git a/Tracking/TrkVertexFitter/TrkVertexSeedFinderUtils/src/GaussianDensityTestAlg.cxx b/Tracking/TrkVertexFitter/TrkVertexSeedFinderUtils/src/GaussianDensityTestAlg.cxx
index e0383c553e3..d6481036388 100644
--- a/Tracking/TrkVertexFitter/TrkVertexSeedFinderUtils/src/GaussianDensityTestAlg.cxx
+++ b/Tracking/TrkVertexFitter/TrkVertexSeedFinderUtils/src/GaussianDensityTestAlg.cxx
@@ -54,14 +54,21 @@ StatusCode GaussianDensityTestAlg::initialize()
ATH_MSG_INFO ("Initializing " << name() << "...");
ATH_CHECK( m_trackParticlesKey.initialize() );
+ ATH_CHECK( m_truthEventsKey.initialize() );
+ ATH_CHECK( m_pileupEventsKey.initialize() );
ATH_CHECK( m_estimator.retrieve() );
ATH_CHECK( m_trackFilter.retrieve() );
+ ATH_CHECK( m_ipEstimator.retrieve() );
// setup histograms/trees
m_h_density = new TH1F("Density", "Density", 800, -200.0, 200.0);
+ m_h_truthDensity = new TH1F("Truth", "Truth", 800, -200.0, 200.0);
+ m_h_truthVertices = new TH1F("TruthVertices", "TruthVertices", 800, -200.0, 200.0);
CHECK( m_iTHistSvc->regHist("/file1/h/density", m_h_density) );
+ CHECK( m_iTHistSvc->regHist("/file1/h/truth", m_h_truthDensity) );
+ CHECK( m_iTHistSvc->regHist("/file1/h/truthvertices", m_h_truthVertices) );
return StatusCode::SUCCESS;
}
@@ -78,24 +85,28 @@ StatusCode GaussianDensityTestAlg::execute()
SG::ReadHandle<xAOD::TrackParticleContainer> trackParticles(m_trackParticlesKey);
+ ATH_MSG_VERBOSE("Selecting tracks");
std::vector<Trk::ITrackLink*> trackVector;
selectTracks(trackParticles.cptr(), trackVector);
std::vector<const Trk::TrackParameters*> perigeeList;
analyzeTracks(trackVector, perigeeList);
+ ATH_MSG_VERBOSE("Using density estimator");
m_estimator->reset();
m_estimator->addTracks(perigeeList);
- //ATH_MSG_DEBUG("Added " << perigeeList.size() << " tracks to density estimator");
-
for (int i = 0; i < 800; i++)
{
double z = -200.0 + 0.25 + i*0.5;
double density = m_estimator->trackDensity(z);
- //ATH_MSG_ALWAYS(" z: " << z << " Density: " << density);
m_h_density->Fill((float) z, (float) density);
}
+ ATH_MSG_VERBOSE("Analyzing MC truth");
+ std::vector<Amg::Vector3D> truth;
+ ATH_CHECK( findTruth(trackVector, truth) );
+ ATH_MSG_VERBOSE("Filling truth vertex histogram");
+ for (auto& v : truth) m_h_truthVertices->Fill( v[2] );
return StatusCode::SUCCESS;
@@ -152,4 +163,147 @@ void GaussianDensityTestAlg::selectTracks(const xAOD::TrackParticleContainer* tr
///////////////////////////////////////////////////////////////////
// Const methods:
///////////////////////////////////////////////////////////////////
+ StatusCode GaussianDensityTestAlg::findTruth(const std::vector<Trk::ITrackLink*>& trackVector, std::vector<Amg::Vector3D>& truth) const
+ {
+ xAOD::TrackParticle::ConstAccessor<ElementLink<xAOD::TruthParticleContainer> > truthParticleAssoc("truthParticleLink");
+
+ SG::ReadHandle<xAOD::TruthEventContainer> signalEvents(m_truthEventsKey);
+ if ( signalEvents.isValid() )
+ {
+ ATH_MSG_VERBOSE("Found signalEvents");
+ for (const xAOD::TruthEventBase* evt : *signalEvents)
+ {
+ const xAOD::TruthVertex* vLink = *(evt->truthVertexLink(0));
+ if (vLink == nullptr) ATH_MSG_ERROR("Invalid truthVertexLink from signalEvents");
+ Amg::Vector3D vTruth(vLink->x(),vLink->y(),vLink->z());
+ int nGoodTracks = 0;
+ for (auto trk : trackVector)
+ {
+ Trk::LinkToXAODTrackParticle* lxtp = dynamic_cast<Trk::LinkToXAODTrackParticle*>(trk);
+ if (lxtp)
+ {
+ bool isAssoc = truthParticleAssoc(**(*lxtp)).isValid();
+ if (isAssoc)
+ {
+ auto assocParticle = truthParticleAssoc(**(*lxtp));
+ ATH_MSG_VERBOSE("Found associated truth particle");
+ for (auto truthParticle : evt->truthParticleLinks())
+ {
+ if (!truthParticle.isValid()) continue;
+ if (assocParticle == truthParticle)
+ {
+ ATH_MSG_VERBOSE("Calling ipSignificance");
+ double significance = ipSignificance(trk->parameters(), &vTruth);
+ ATH_MSG_VERBOSE("ipSignificance returned " << significance);
+ if (significance <= m_significanceTruthCut)
+ {
+ ATH_MSG_VERBOSE("Adding good track");
+ nGoodTracks++;
+ const Trk::Perigee* perigee = dynamic_cast<const Trk::Perigee*>(trk->parameters());
+ if (perigee == nullptr) ATH_MSG_ERROR("Invalid Perigee");
+ m_h_truthDensity->Fill(vLink->z());
+ ATH_MSG_VERBOSE("Filled truth density histogram");
+ }
+ break;
+ }
+ }
+ }
+ else
+ {
+ ATH_MSG_VERBOSE("No associated truth particle found");
+ }
+ }
+ }
+ if (nGoodTracks >= m_truthVertexTracks)
+ {
+ truth.push_back(vTruth);
+ }
+ }
+ }
+ else
+ {
+ ATH_MSG_WARNING("No TruthEventContainer found");
+ }
+
+ SG::ReadHandle<xAOD::TruthPileupEventContainer> pileupEvents(m_pileupEventsKey);
+ if ( pileupEvents.isValid() )
+ {
+ ATH_MSG_VERBOSE("Found pileupEvents");
+ for (const xAOD::TruthEventBase* evt : *pileupEvents)
+ {
+ const xAOD::TruthVertex* vLink = *(evt->truthVertexLink(0));
+ if (vLink == nullptr) ATH_MSG_ERROR("Invalid truthVertexLink from pileupEvents");
+ Amg::Vector3D vTruth(vLink->x(),vLink->y(),vLink->z());
+ int nGoodTracks = 0;
+ for (auto trk : trackVector)
+ {
+ Trk::LinkToXAODTrackParticle* lxtp = dynamic_cast<Trk::LinkToXAODTrackParticle*>(trk);
+ if (lxtp)
+ {
+ bool isAssoc = truthParticleAssoc(**(*lxtp)).isValid();
+ if (isAssoc)
+ {
+ auto assocParticle = truthParticleAssoc(**(*lxtp));
+ ATH_MSG_VERBOSE("Found associated truth particle");
+ for (auto truthParticle : evt->truthParticleLinks())
+ {
+ if (!truthParticle.isValid()) continue;
+ if (assocParticle == truthParticle)
+ {
+ ATH_MSG_VERBOSE("Calling ipSignificance");
+ double significance = ipSignificance(trk->parameters(), &vTruth);
+ ATH_MSG_VERBOSE("ipSignificance returned " << significance);
+ if (significance <= m_significanceTruthCut)
+ {
+ ATH_MSG_VERBOSE("Adding good track");
+ nGoodTracks++;
+ const Trk::Perigee* perigee = dynamic_cast<const Trk::Perigee*>(trk->parameters());
+ if (perigee == nullptr) ATH_MSG_ERROR("Invalid Perigee");
+ m_h_truthDensity->Fill(vLink->z());
+ ATH_MSG_VERBOSE("Filled truth density histogram");
+ }
+ break;
+ }
+ }
+ }
+ else
+ {
+ ATH_MSG_VERBOSE("No associated truth particle found");
+ }
+ }
+ }
+ if (nGoodTracks >= m_truthVertexTracks)
+ {
+ truth.push_back(vTruth);
+ }
+ }
+ }
+ else
+ {
+ ATH_MSG_WARNING("No TruthPileupEventContainer found");
+ }
+ return StatusCode::SUCCESS;
+ }
+
+ double GaussianDensityTestAlg::ipSignificance( const Trk::TrackParameters* params,
+ const Amg::Vector3D * vertex ) const
+ {
+ xAOD::Vertex v;
+ v.makePrivateStore();
+ v.setPosition(*vertex);
+ v.setCovariancePosition(AmgSymMatrix(3)::Zero(3,3));
+ v.setFitQuality(0., 0.);
+
+ double significance = 0.0;
+ const ImpactParametersAndSigma* ipas = m_ipEstimator->estimate( params, &v );
+ if ( ipas != nullptr )
+ {
+ if ( ipas->sigmad0 > 0 && ipas->sigmaz0 > 0)
+ {
+ significance = sqrt( pow(ipas->IPd0/ipas->sigmad0,2) + pow(ipas->IPz0/ipas->sigmaz0,2) );
+ }
+ delete ipas;
+ }
+ return significance;
+ }
}
diff --git a/Tracking/TrkVertexFitter/TrkVertexSeedFinderUtils/src/GaussianDensityTestAlg.h b/Tracking/TrkVertexFitter/TrkVertexSeedFinderUtils/src/GaussianDensityTestAlg.h
index f959c3304c3..605f64fe7ce 100644
--- a/Tracking/TrkVertexFitter/TrkVertexSeedFinderUtils/src/GaussianDensityTestAlg.h
+++ b/Tracking/TrkVertexFitter/TrkVertexSeedFinderUtils/src/GaussianDensityTestAlg.h
@@ -25,10 +25,13 @@
#include "TrkParticleBase/TrackParticleBaseCollection.h"
#include "TrkParameters/TrackParameters.h"
#include "xAODTracking/TrackParticleContainer.h"
+#include "xAODTruth/TruthEventContainer.h"
+#include "xAODTruth/TruthPileupEventContainer.h"
#include "InDetBeamSpotService/IBeamCondSvc.h"
#include "TrkVertexFitterInterfaces/IVertexTrackDensityEstimator.h"
#include "InDetTrackSelectionTool/IInDetTrackSelectionTool.h"
+#include "TrkVertexFitterInterfaces/ITrackToVertexIPEstimator.h"
//Amg
#include "GeoPrimitives/GeoPrimitives.h"
@@ -69,6 +72,9 @@ class GaussianDensityTestAlg
///////////////////////////////////////////////////////////////////
// Const methods:
///////////////////////////////////////////////////////////////////
+ double ipSignificance(const Trk::TrackParameters* params, const Amg::Vector3D * vertex) const;
+
+ StatusCode findTruth(const std::vector<Trk::ITrackLink*>& trackVector, std::vector<Amg::Vector3D>& truth) const;
///////////////////////////////////////////////////////////////////
// Non-const methods:
@@ -84,6 +90,16 @@ class GaussianDensityTestAlg
///////////////////////////////////////////////////////////////////
private:
// Properties
+ Gaudi::Property<double> m_significanceTruthCut { this,
+ "SignificanceTruthCut",
+ 3.0,
+ "Reco track must pass within this many sigma of pp vertex to be good" };
+
+ Gaudi::Property<int> m_truthVertexTracks { this,
+ "MinTruthVertexTracks",
+ 2,
+ "Minimum associated reconstructed tracks for vertex to be considered visible" };
+
bool m_useBeamConstraint;
// Tools
@@ -94,6 +110,9 @@ class GaussianDensityTestAlg
ToolHandle< Trk::IVertexTrackDensityEstimator > m_estimator { this, "Estimator", "Trk::GaussianTrackDensity",
"Track density function" };
+ ToolHandle< Trk::ITrackToVertexIPEstimator > m_ipEstimator { this, "IPEstimator", "Trk::TrackToVertexIPEstimator",
+ "Impact point estimator" };
+
// Non-property private data
ServiceHandle< IBeamCondSvc > m_iBeamCondSvc;
@@ -103,9 +122,17 @@ class GaussianDensityTestAlg
SG::ReadHandleKey<xAOD::TrackParticleContainer> m_trackParticlesKey { this, "TrackParticles", "InDetTrackParticles",
"Input track particle collection" };
+ SG::ReadHandleKey<xAOD::TruthEventContainer> m_truthEventsKey { this, "TruthEvents", "TruthEvents",
+ "Key for truth event collection" };
+
+ SG::ReadHandleKey<xAOD::TruthPileupEventContainer> m_pileupEventsKey { this, "TruthPileupEvents", "TruthPileupEvents",
+ "Key for truth pileup event collection" };
+
/// Histograms and trees
TH1* m_h_density;
+ TH1* m_h_truthDensity;
+ TH1* m_h_truthVertices;
}; // class
} // namespace
diff --git a/Tracking/TrkVertexFitter/TrkVertexSeedFinderUtils/src/GaussianTrackDensity.cxx b/Tracking/TrkVertexFitter/TrkVertexSeedFinderUtils/src/GaussianTrackDensity.cxx
index 2e4aaafbf3d..cdc4c70dce3 100644
--- a/Tracking/TrkVertexFitter/TrkVertexSeedFinderUtils/src/GaussianTrackDensity.cxx
+++ b/Tracking/TrkVertexFitter/TrkVertexSeedFinderUtils/src/GaussianTrackDensity.cxx
@@ -110,7 +110,6 @@ namespace Trk
double cov_zz = perigeeCov(Trk::z0, Trk::z0);
if (cov_zz <= 0 ) continue;
double cov_dz = perigeeCov(Trk::d0, Trk::z0);
- //ATH_MSG_DEBUG("z0:" << z0 << " d0: " << d0 << " covdd, covdz, covzz " << cov_dd << " " << cov_dz << " " << cov_zz);
double covDeterminant = cov_dd*cov_zz - cov_dz*cov_dz;
if ( covDeterminant <= 0 ) continue;
double constantTerm = -(d0*d0*cov_zz + z0*z0*cov_dd + 2*d0*z0*cov_dz) / (2*covDeterminant);
@@ -121,8 +120,7 @@ namespace Trk
discriminant = sqrt(discriminant);
double zMax = (-linearTerm - discriminant)/(2*quadraticTerm);
double zMin = (-linearTerm + discriminant)/(2*quadraticTerm);
- //ATH_MSG_DEBUG("zMin: " << zMin << " zMax: " << zMax);
- constantTerm -= log(2*Gaudi::Units::pi*covDeterminant);
+ constantTerm -= log(2*Gaudi::Units::pi*sqrt(covDeterminant));
m_trackMap.emplace(std::piecewise_construct,
std::forward_as_tuple(*itrk),
std::forward_as_tuple(constantTerm, linearTerm, quadraticTerm, zMin, zMax));
--
GitLab