Add method to return global maximum of density, and a histogram to test it

Former-commit-id: c776e396b639144bc6ce359ed2aa6d70e3d2d50b

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.
         */

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

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'" ]

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

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

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;