Add truth vertex analysis to test algorithm; fix bug in density calculation;...

Add truth vertex analysis to test algorithm; fix bug in density calculation; tested successful on full-truth MC.


Former-commit-id: 1d841f3b9826d4159dc87e0d78ddc9a8a2ecb096

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

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
 

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

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

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