diff --git a/Tracking/TrkVertexFitter/TrkVertexFitterInterfaces/TrkVertexFitterInterfaces/IVertexTrackDensityEstimator.h b/Tracking/TrkVertexFitter/TrkVertexFitterInterfaces/TrkVertexFitterInterfaces/IVertexTrackDensityEstimator.h
index f6d60e480cb7da24a35b0ba2c80bc7d03f994884..f8b9487109a863b71c3b70e191e396a652ec64a2 100755
--- a/Tracking/TrkVertexFitter/TrkVertexFitterInterfaces/TrkVertexFitterInterfaces/IVertexTrackDensityEstimator.h
+++ b/Tracking/TrkVertexFitter/TrkVertexFitterInterfaces/TrkVertexFitterInterfaces/IVertexTrackDensityEstimator.h
@@ -70,12 +70,12 @@ namespace Trk
/**
* Evaluate the density function at the specified coordinate along the beam-line.
*/
- virtual double trackDensity(double z) const = 0;
+ virtual double trackDensity(double z) = 0;
/*
* Evaluate the density function and its first two derivatives at the specified coordinate.
*/
- virtual void trackDensity(double z, double& density, double& firstDerivative, double& secondDerivative) const = 0;
+ virtual void trackDensity(double z, double& density, double& firstDerivative, double& secondDerivative) = 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 43b555f3f12352d9fd11ce5d552463a64e868bfa..61f184bbf2b510851de9424eb483757eccee5b8e 100644
--- a/Tracking/TrkVertexFitter/TrkVertexSeedFinderUtils/TrkVertexSeedFinderUtils/GaussianTrackDensity.h
+++ b/Tracking/TrkVertexFitter/TrkVertexSeedFinderUtils/TrkVertexSeedFinderUtils/GaussianTrackDensity.h
@@ -8,12 +8,13 @@
#include "AthenaBaseComps/AthAlgTool.h"
#include "TrkVertexFitterInterfaces/IVertexTrackDensityEstimator.h"
-#include <unordered_map>
+#include <map>
namespace Trk
{
class Track;
+ class GaussianTrackDensity;
/**
@class GaussianTrackDensity
@@ -64,57 +65,62 @@ namespace Trk
/**
* Evaluate the density function at the specified coordinate along the beam-line.
*/
- virtual double trackDensity(double z) const;
+ virtual double trackDensity(double z);
/**
* Evaluate the density and its first two derivatives at the specified coordinate.
*/
- virtual void trackDensity(double z, double& density, double& firstDerivative, double& secondDerivativec) const;
+ virtual void trackDensity(double z, double& density, double& firstDerivative, double& secondDerivativec);
/**
* Resets the internal state of the tool, forgetting all tracks previously added.
*/
virtual void reset();
+ // functor to compare two unordered_map Key values for equality
+ struct pred_perigee {
+ bool operator()(const Trk::Perigee& left, const Trk::Perigee& right) const
+ {
+ return left.parameters()[Trk::z0] < right.parameters()[Trk::z0];
+ }
+ };
- private:
struct TrackEntry
{
- TrackEntry(double c0, double c1, double c2);
+ TrackEntry(double c0, double c1, double c2, double covz);
// Cached information for a single track
double constant; // z-independent factor
double c_1; // linear coefficient in exponent
double c_2; // quadratic coefficient in exponent
+ double cov_z; // z0 variance from track error matrix
+ std::map< Perigee, TrackEntry, pred_perigee >::const_iterator start;
+ // will point to the left-most TrackEntry close enough to affect the weight at this one
+ std::map< Perigee, TrackEntry, pred_perigee >::const_iterator finish;
+ // will point to the right-most TrackEntry close enough to affect the weight at this one
};
- // functor to hash key for unordered_map
- struct hash_perigee {
- size_t operator()(const Trk::Perigee& perigee) const
- {
- return
- std::hash<double>()(perigee.parameters()[Trk::d0]) ^
- std::hash<double>()(perigee.parameters()[Trk::z0]) ^
- std::hash<double>()(perigee.parameters()[Trk::phi]) ^
- std::hash<double>()(perigee.parameters()[Trk::theta]) ^
- std::hash<double>()(perigee.parameters()[Trk::qOverP]);
- }
- };
+ typedef std::map< Perigee,
+ GaussianTrackDensity::TrackEntry,
+ GaussianTrackDensity::pred_perigee > trackMap;
+ typedef std::map< Perigee,
+ GaussianTrackDensity::TrackEntry,
+ GaussianTrackDensity::pred_perigee >::const_iterator trackMapIterator;
+
+ private:
+
+ trackMapIterator findStart(double z);
+
+ trackMapIterator findFinish(double z);
- // functor to compare two unordered_map Key values for equality
- struct pred_perigee {
- bool operator()(const Trk::Perigee& left, const Trk::Perigee& right) const
- {
- return
- (left.parameters()[Trk::d0] == right.parameters()[Trk::d0]) &&
- (left.parameters()[Trk::z0] == right.parameters()[Trk::z0]) &&
- (left.parameters()[Trk::phi] == right.parameters()[Trk::phi]) &&
- (left.parameters()[Trk::theta] == right.parameters()[Trk::theta]) &&
- (left.parameters()[Trk::qOverP] == right.parameters()[Trk::qOverP]);
- }
- };
// Cache for track information
- std::unordered_map< Trk::Perigee, Trk::GaussianTrackDensity::TrackEntry, hash_perigee, pred_perigee> m_trackMap;
+ //std::unordered_map< Trk::Perigee, Trk::GaussianTrackDensity::TrackEntry, hash_perigee, pred_perigee> m_trackMap;
+ trackMap m_trackMap;
+
+ // Indicates whether track data needs to be refreshed before evaluation
+ bool m_dirty;
+
+ void prepareTracks();
// Cuts set by configurable properties
@@ -123,6 +129,10 @@ namespace Trk
"MaxD0Significance",
3.5,
"Maximum radial impact parameter significance to use track" };
+ Gaudi::Property<double> m_z0MaxSignificance { this,
+ "MaxZ0Significance",
+ 5.0,
+ "Maximum longitudinal impact parameter significance to include track in weight"};
};
}
diff --git a/Tracking/TrkVertexFitter/TrkVertexSeedFinderUtils/src/GaussianTrackDensity.cxx b/Tracking/TrkVertexFitter/TrkVertexSeedFinderUtils/src/GaussianTrackDensity.cxx
index 8b626a27a0b45b40916442b930c461434d0ea95c..7a16f5827c541e77772f230e9a8a1a0fbce9a049 100644
--- a/Tracking/TrkVertexFitter/TrkVertexSeedFinderUtils/src/GaussianTrackDensity.cxx
+++ b/Tracking/TrkVertexFitter/TrkVertexSeedFinderUtils/src/GaussianTrackDensity.cxx
@@ -13,7 +13,8 @@
namespace Trk
{
GaussianTrackDensity::GaussianTrackDensity(const std::string& t, const std::string& n, const IInterface* p):
- AthAlgTool(t, n, p)
+ AthAlgTool(t, n, p),
+ m_dirty(false)
{
declareInterface<IVertexTrackDensityEstimator>(this);
}
@@ -28,30 +29,55 @@ namespace Trk
return StatusCode::SUCCESS;
}
- double GaussianTrackDensity::trackDensity(double z) const
+ GaussianTrackDensity::trackMapIterator GaussianTrackDensity::findStart(double z)
{
- double sum = 0.0;
- for (auto& entry : m_trackMap)
+ trackMapIterator result = m_trackMap.cbegin();
+ while (result != m_trackMap.cend() && result->first.parameters()[Trk::z0] <= z)
{
- sum += entry.second.constant * exp(z*(entry.second.c_1 + z*entry.second.c_2));
+ result++;
}
- return sum;
+ // always goes one past the track immediately to left of z
+ return (--result)->second.start;
+ }
+
+ GaussianTrackDensity::trackMapIterator GaussianTrackDensity::findFinish(double z)
+ {
+ trackMapIterator result = m_trackMap.cend();
+ result--;
+ while (result != m_trackMap.cbegin() && result->first.parameters()[Trk::z0] > z)
+ {
+ result--;
+ }
+ if (result->first.parameters()[Trk::z0] < z) result++;
+ return result->second.finish;
}
- void GaussianTrackDensity::trackDensity(double z, double& density, double& firstDerivative, double& secondDerivative) const
+ void GaussianTrackDensity::trackDensity(double z, double& density, double& firstDerivative, double& secondDerivative)
{
+ if (m_dirty) prepareTracks();
density = 0.0;
firstDerivative = 0.0;
secondDerivative = 0.0;
- for (auto& entry : m_trackMap)
+ for (auto it = findStart(z); it != findFinish(z); ++it)
{
- double delta = entry.second.constant * exp(z*(entry.second.c_1 + z*entry.second.c_2));
+ double delta = it->second.constant * exp(z*(it->second.c_1 + z*it->second.c_2));
density += delta;
- double qPrime = entry.second.c_1 + 2*z*entry.second.c_2;
+ double qPrime = it->second.c_1 + 2*z*it->second.c_2;
double deltaPrime = delta * qPrime;
firstDerivative += deltaPrime;
- secondDerivative += 2*entry.second.c_2*delta + qPrime*deltaPrime;
+ secondDerivative += 2*it->second.c_2*delta + qPrime*deltaPrime;
+ }
+ }
+
+ double GaussianTrackDensity::trackDensity(double z)
+ {
+ if (m_dirty) prepareTracks();
+ double sum = 0.0;
+ for (auto it = findStart(z); it != findFinish(z); ++it)
+ {
+ sum += it->second.constant * exp(z*(it->second.c_1 + z*it->second.c_2));
}
+ return sum;
}
void GaussianTrackDensity::addTracks(const std::vector<const Trk::Track*>& vectorTrk)
@@ -92,6 +118,7 @@ namespace Trk
if ( cov_dd <= 0 ) continue;
if ( d0*d0/cov_dd > significanceCut ) continue;
double cov_zz = perigeeCov(Trk::z0, Trk::z0);
+ if (cov_zz <= 0 ) continue;
double cov_dz = perigeeCov(Trk::d0, Trk::z0);
double covDeterminant = cov_dd*cov_zz - cov_dz*cov_dz;
if ( covDeterminant <= 0 ) continue;
@@ -101,7 +128,8 @@ namespace Trk
double quadraticTerm = -cov_dd / (2*covDeterminant);
m_trackMap.emplace(std::piecewise_construct,
std::forward_as_tuple(*itrk),
- std::forward_as_tuple(constantFactor, linearTerm, quadraticTerm));
+ std::forward_as_tuple(constantFactor, linearTerm, quadraticTerm, cov_zz));
+ m_dirty = true;
}
}
}
@@ -115,15 +143,43 @@ namespace Trk
{
if (m_trackMap.count(*itrk) == 0) continue;
m_trackMap.erase(*itrk);
+ m_dirty = true;
+ }
+ }
+ }
+
+ void GaussianTrackDensity::prepareTracks()
+ {
+ trackMapIterator left = m_trackMap.cbegin();
+ double z0Cut = m_z0MaxSignificance * m_z0MaxSignificance;
+ for (auto& entry : m_trackMap)
+ {
+ double myZ0 = entry.first.parameters()[Trk::z0];
+ // left should advance to the left-most track to the left of "entry" that is within n sigma of it,
+ // or no track to the left of "entry" passes, it should stop at "entry" itself
+ while (pow(myZ0 - left->first.parameters()[Trk::z0],2) > z0Cut*left->second.cov_z &&
+ left->first.parameters()[Trk::z0] < myZ0)
+ left++;
+ // right should advance to the element after the right-most element within n-sigma of "entry"
+ trackMapIterator right = m_trackMap.cend();
+ while((--right)->first.parameters()[Trk::z0] > myZ0)
+ {
+ if (pow(myZ0 - right->first.parameters()[Trk::z0],2) < z0Cut*right->second.cov_z) break;
}
+ entry.second.start = left;
+ entry.second.finish = ++right;
}
+ m_dirty = false;
}
void GaussianTrackDensity::reset()
{
m_trackMap.clear();
+ m_dirty = false;
}
- GaussianTrackDensity::TrackEntry::TrackEntry(double c0, double c1, double c2) : constant(c0), c_1(c1), c_2(c2) { }
+ GaussianTrackDensity::TrackEntry::TrackEntry(double c0, double c1, double c2, double covz)
+ : constant(c0), c_1(c1), c_2(c2), cov_z(covz), start(), finish()
+{ }
}