Attempt to optimize so only tracks with a reasonable distance of the test...

Attempt to optimize so only tracks with a reasonable distance of the test point are used in the computation.


Former-commit-id: 061ed3f39614b5e3a5340889bcd52d5150a2281d

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.

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"};
 
   };
 }

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() 
+{ }
 
 }