From 725a4c220316b2ae604b58a2eab097d2aa5c7926 Mon Sep 17 00:00:00 2001
From: Dave Casper <dcasper@uci.edu>
Date: Thu, 8 Mar 2018 08:16:08 -0800
Subject: [PATCH] Add method to return global maximum of density, and a
histogram to test it
Former-commit-id: c776e396b639144bc6ce359ed2aa6d70e3d2d50b
---
.../IVertexTrackDensityEstimator.h | 6 ++-
.../GaussianTrackDensity.h | 20 +++++++++
.../GaussianDensityTestAlg_jobOptions.py | 3 +-
.../src/GaussianDensityTestAlg.cxx | 25 +++++++----
.../src/GaussianDensityTestAlg.h | 1 +
.../src/GaussianTrackDensity.cxx | 41 +++++++++++++++++--
6 files changed, 83 insertions(+), 13 deletions(-)
diff --git a/Tracking/TrkVertexFitter/TrkVertexFitterInterfaces/TrkVertexFitterInterfaces/IVertexTrackDensityEstimator.h b/Tracking/TrkVertexFitter/TrkVertexFitterInterfaces/TrkVertexFitterInterfaces/IVertexTrackDensityEstimator.h
index f6d60e480cb..09e109b66d7 100755
--- a/Tracking/TrkVertexFitter/TrkVertexFitterInterfaces/TrkVertexFitterInterfaces/IVertexTrackDensityEstimator.h
+++ b/Tracking/TrkVertexFitter/TrkVertexFitterInterfaces/TrkVertexFitterInterfaces/IVertexTrackDensityEstimator.h
@@ -44,7 +44,6 @@ namespace Trk
*/
static const InterfaceID& interfaceID() { return IID_IVertexTrackDensityEstimator; };
-
/**
* Adds a list of tracks, whose impact parameters will contribute to the density function.
*/
@@ -77,6 +76,11 @@ namespace Trk
*/
virtual void trackDensity(double z, double& density, double& firstDerivative, double& secondDerivative) const = 0;
+ /*
+ * Find position of global maximum for density function
+ */
+ virtual double globalMaximum() const = 0;
+
/**
* Resets the internal state of the tool, forgetting all tracks previously added.
*/
diff --git a/Tracking/TrkVertexFitter/TrkVertexSeedFinderUtils/TrkVertexSeedFinderUtils/GaussianTrackDensity.h b/Tracking/TrkVertexFitter/TrkVertexSeedFinderUtils/TrkVertexSeedFinderUtils/GaussianTrackDensity.h
index 1bb7c927b7f..9d64a2d92bf 100644
--- a/Tracking/TrkVertexFitter/TrkVertexSeedFinderUtils/TrkVertexSeedFinderUtils/GaussianTrackDensity.h
+++ b/Tracking/TrkVertexFitter/TrkVertexSeedFinderUtils/TrkVertexSeedFinderUtils/GaussianTrackDensity.h
@@ -72,6 +72,11 @@ namespace Trk
*/
virtual void trackDensity(double z, double& density, double& firstDerivative, double& secondDerivativec) const;
+ /**
+ * Find location of global maximum for density function
+ */
+ virtual double globalMaximum() const;
+
/**
* Resets the internal state of the tool, forgetting all tracks previously added.
*/
@@ -139,6 +144,13 @@ namespace Trk
GaussianTrackDensity::pred_entry_by_min >::const_iterator upperMapIterator;
private:
+ inline void updateMaximum(double trialZ, double trialValue, double& maxZ, double& maxValue) const
+ { if (trialValue > maxValue) { maxZ = trialZ; maxValue = trialValue; } }
+
+
+ inline double stepSize(double y, double dy, double ddy) const
+ { return ( m_gaussStep ? 1/(dy/y-ddy/dy) : -dy/ddy ); }
+
// Cache for track information
trackMap m_trackMap;
@@ -152,11 +164,19 @@ namespace Trk
"MaxD0Significance",
3.5,
"Maximum radial impact parameter significance to use track" };
+
+ // Tracks within this many sigma(z) are added to weight; increasing cut trades CPU efficiency for improved smoothness in tails
Gaudi::Property<double> m_z0MaxSignificance { this,
"MaxZ0Significance",
12.0,
"Maximum longitudinal impact parameter significance to include track in weight"};
+ // Assumed shape of density function near peak; Gaussian (true) or parabolic (false)
+ Gaudi::Property<bool> m_gaussStep { this,
+ "GaussianStep",
+ true,
+ "Peak search: True means assume Gaussian behavior, False means Newton/parabolic" };
+
};
}
#endif
diff --git a/Tracking/TrkVertexFitter/TrkVertexSeedFinderUtils/share/GaussianDensityTestAlg_jobOptions.py b/Tracking/TrkVertexFitter/TrkVertexSeedFinderUtils/share/GaussianDensityTestAlg_jobOptions.py
index 6784f0b8c1e..553052c06c2 100644
--- a/Tracking/TrkVertexFitter/TrkVertexSeedFinderUtils/share/GaussianDensityTestAlg_jobOptions.py
+++ b/Tracking/TrkVertexFitter/TrkVertexSeedFinderUtils/share/GaussianDensityTestAlg_jobOptions.py
@@ -7,7 +7,8 @@ import AthenaPython.ConfigLib as apcl
cfg = apcl.AutoCfg(name = 'InDetRecExampleAutoConfig', input_files=athenaCommonFlags.FilesInput())
cfg.configure_job()
-theApp.EvtMax = 1
+theApp.EvtMax = -1
+#athenaCommonFlags.SkipEvents = 1
svcMgr += CfgMgr.THistSvc()
svcMgr.THistSvc.Output = [ "file1 DATAFILE='gaussianDensity.root' OPT='RECREATE'" ]
diff --git a/Tracking/TrkVertexFitter/TrkVertexSeedFinderUtils/src/GaussianDensityTestAlg.cxx b/Tracking/TrkVertexFitter/TrkVertexSeedFinderUtils/src/GaussianDensityTestAlg.cxx
index d6481036388..f7940784e21 100644
--- a/Tracking/TrkVertexFitter/TrkVertexSeedFinderUtils/src/GaussianDensityTestAlg.cxx
+++ b/Tracking/TrkVertexFitter/TrkVertexSeedFinderUtils/src/GaussianDensityTestAlg.cxx
@@ -65,10 +65,12 @@ StatusCode GaussianDensityTestAlg::initialize()
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);
+ m_h_modeCheck = new TH1F("ModeOffset", "ModeOffset", 100, -10.0, 10.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) );
+ CHECK( m_iTHistSvc->regHist("/file1/h/modeoffset", m_h_modeCheck) );
return StatusCode::SUCCESS;
}
@@ -96,17 +98,17 @@ StatusCode GaussianDensityTestAlg::execute()
m_estimator->reset();
m_estimator->addTracks(perigeeList);
- for (int i = 0; i < 800; i++)
- {
- double z = -200.0 + 0.25 + i*0.5;
- double density = m_estimator->trackDensity(z);
- m_h_density->Fill((float) z, (float) density);
- }
+ // for (int i = 0; i < 800; i++)
+ // {
+ // double z = -200.0 + 0.25 + i*0.5;
+ // double density = m_estimator->trackDensity(z);
+ // 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] );
+ //for (auto& v : truth) m_h_truthVertices->Fill( v[2] );
return StatusCode::SUCCESS;
@@ -165,6 +167,9 @@ void GaussianDensityTestAlg::selectTracks(const xAOD::TrackParticleContainer* tr
///////////////////////////////////////////////////////////////////
StatusCode GaussianDensityTestAlg::findTruth(const std::vector<Trk::ITrackLink*>& trackVector, std::vector<Amg::Vector3D>& truth) const
{
+ double modeClosestDistance = std::numeric_limits<double>::max();
+ double mode = m_estimator->globalMaximum();
+
xAOD::TrackParticle::ConstAccessor<ElementLink<xAOD::TruthParticleContainer> > truthParticleAssoc("truthParticleLink");
SG::ReadHandle<xAOD::TruthEventContainer> signalEvents(m_truthEventsKey);
@@ -201,7 +206,7 @@ void GaussianDensityTestAlg::selectTracks(const xAOD::TrackParticleContainer* tr
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());
+ //m_h_truthDensity->Fill(vLink->z());
ATH_MSG_VERBOSE("Filled truth density histogram");
}
break;
@@ -217,6 +222,7 @@ void GaussianDensityTestAlg::selectTracks(const xAOD::TrackParticleContainer* tr
if (nGoodTracks >= m_truthVertexTracks)
{
truth.push_back(vTruth);
+ modeClosestDistance = std::min(modeClosestDistance, mode - vTruth[2]);
}
}
}
@@ -275,6 +281,7 @@ void GaussianDensityTestAlg::selectTracks(const xAOD::TrackParticleContainer* tr
if (nGoodTracks >= m_truthVertexTracks)
{
truth.push_back(vTruth);
+ modeClosestDistance = std::min(modeClosestDistance, mode - vTruth[2]);
}
}
}
@@ -282,6 +289,8 @@ void GaussianDensityTestAlg::selectTracks(const xAOD::TrackParticleContainer* tr
{
ATH_MSG_WARNING("No TruthPileupEventContainer found");
}
+
+ m_h_modeCheck->Fill( modeClosestDistance );
return StatusCode::SUCCESS;
}
diff --git a/Tracking/TrkVertexFitter/TrkVertexSeedFinderUtils/src/GaussianDensityTestAlg.h b/Tracking/TrkVertexFitter/TrkVertexSeedFinderUtils/src/GaussianDensityTestAlg.h
index 605f64fe7ce..2ef16e3ff4c 100644
--- a/Tracking/TrkVertexFitter/TrkVertexSeedFinderUtils/src/GaussianDensityTestAlg.h
+++ b/Tracking/TrkVertexFitter/TrkVertexSeedFinderUtils/src/GaussianDensityTestAlg.h
@@ -133,6 +133,7 @@ class GaussianDensityTestAlg
TH1* m_h_density;
TH1* m_h_truthDensity;
TH1* m_h_truthVertices;
+ TH1* m_h_modeCheck;
}; // class
} // namespace
diff --git a/Tracking/TrkVertexFitter/TrkVertexSeedFinderUtils/src/GaussianTrackDensity.cxx b/Tracking/TrkVertexFitter/TrkVertexSeedFinderUtils/src/GaussianTrackDensity.cxx
index cdc4c70dce3..0d7f3e3bfdb 100644
--- a/Tracking/TrkVertexFitter/TrkVertexSeedFinderUtils/src/GaussianTrackDensity.cxx
+++ b/Tracking/TrkVertexFitter/TrkVertexSeedFinderUtils/src/GaussianTrackDensity.cxx
@@ -61,14 +61,49 @@ namespace Trk
if (firstLoop->second.lowerBound > z || finalLoop->second.upperBound < z) return sum;
for (auto itrk = firstLoop; itrk != finalLoop && itrk != m_trackMap.end(); itrk++)
{
- //ATH_MSG_DEBUG("@z=" << z << " adding contrib from z0=" << itrk->first.parameters()[Trk::z0] <<
- // "c_0, c_1, c2 = " << itrk->second.c_0 << " " << itrk->second.c_1 << " " <<
- // itrk->second.c_2 << " arg=" << (itrk->second.c_0+z*(itrk->second.c_1 + z*itrk->second.c_2)));
sum += exp(itrk->second.c_0+z*(itrk->second.c_1 + z*itrk->second.c_2));
}
return sum;
}
+ double GaussianTrackDensity::globalMaximum() const
+ {
+ // strategy:
+ // the global maximum must be somewhere near a track...
+ // since we can calculate the first and second derivatives, at each point we can determine
+ // a) whether the function is curved up (minimum) or down (maximum)
+ // b) the distance to nearest maximum, assuming either Newton (parabolic) or Gaussian local behavior
+ //
+ // For each track where the second derivative is negative, find step to nearest maximum
+ // Take that step, and then do one final refinement
+ // The largest density encountered in this procedure (after checking all tracks) is considered the maximum
+ //
+ double maximumPosition = 0.0;
+ double maximumDensity = 0.0;
+
+ for (const auto& entry : m_trackMap)
+ {
+ double trialZ = entry.first.parameters()[Trk::z0];
+ double density = 0.0;
+ double slope = 0.0;
+ double curvature = 0.0;
+ trackDensity( trialZ, density, slope, curvature );
+ if ( curvature >= 0.0 || density <= 0.0 ) continue;
+ updateMaximum( trialZ, density, maximumPosition, maximumDensity );
+ trialZ += stepSize( density, slope, curvature );
+ trackDensity( trialZ, density, slope, curvature );
+ if ( curvature >= 0.0 || density <= 0.0 ) continue;
+ updateMaximum( trialZ, density, maximumPosition, maximumDensity );
+ trialZ += stepSize( density, slope, curvature );
+ trackDensity( trialZ, density, slope, curvature );
+ if ( curvature >= 0.0 || density <= 0.0) continue;
+ updateMaximum( trialZ, density, maximumPosition, maximumDensity );
+ }
+ if ( maximumDensity <= 0 ) ATH_MSG_WARNING("Global maximum at density of 0; track map contains " <<
+ m_trackMap.size() << " tracks");
+ return maximumPosition;
+ }
+
void GaussianTrackDensity::addTracks(const std::vector<const Track*>& vectorTrk)
{
std::vector<const TrackParameters*> perigeeList;
--
GitLab