Multi adaptive vertex fitter

Implements the multi adaptive vertex fitter and all related required tools.

Tests/Core/Vertexing/IterativeVertexFinderTests.cpp

@@ -24,7 +24,7 @@
 #include "Acts/Propagator/EigenStepper.hpp"
 
 #include "Acts/Vertexing/VertexFitterTools/FullBilloirVertexFitter.hpp"
-
+#include "Acts/Vertexing/VertexEventData/Vertex.hpp"
 #include "Acts/Vertexing/VertexFinderTools/IterativeVertexFinder.hpp"
 #include "Acts/Vertexing/VertexFinderUtils/FsmwMode1dFinder.hpp"
 #include "Acts/Vertexing/VertexFinderUtils/TrackToVertexIPEstimator.hpp"
@@ -57,22 +57,27 @@ namespace Test {
   std::uniform_real_distribution<> resAngDist(0., 0.1);
   // Track q/p resolution distribution
   std::uniform_real_distribution<> resQoPDist(-0.01, 0.01);
+  // Number of vertices per test event distribution
+  std::uniform_int_distribution<> nVertexDist(1, 6);
+  // Number of tracks per vertex distribution
+  std::uniform_int_distribution<> nTracksDist(5, 15);
+
 
   ///
   /// @brief Unit test for IterativeVertexFinder
   ///
   BOOST_AUTO_TEST_CASE(iterative_finder_test)
   {
-    unsigned int nTests = 1;
+  	bool debug = false;
 
-    for (unsigned int iTest = 0; iTest < nTests; ++iTest) {
+  	// Set up RNG
+    std::random_device rd;
+    std::mt19937       gen(rd());
 
-      // Number of tracks
-      unsigned int nTracks = 30;
+  	// Number of test events
+    unsigned int nEvents = 30; // = nTest
 
-      // Set up RNG
-      std::random_device rd;
-      std::mt19937       gen(rd());
+    for (unsigned int iEvent = 0; iEvent < nEvents; ++iEvent) {
 
       // Set up constant B-Field
       ConstantBField bField(Vector3D(0., 0., 1.) * units::_T);
@@ -83,92 +88,162 @@ namespace Test {
       // Set up propagator with void navigator
       Propagator<EigenStepper<ConstantBField>> propagator(stepper);
 
-      // Set up Billoir Vertex Fitter
-      FullBilloirVertexFitter<ConstantBField,
-                              BoundParameters,
-                              Propagator<EigenStepper<ConstantBField>>>::Config
-          vertexFitterCfg(bField);
-
-      std::
-          unique_ptr<FullBilloirVertexFitter<ConstantBField,
-                                             BoundParameters,
-                                             Propagator<EigenStepper<ConstantBField>>>>
-              bFitterPtr
-
-          = std::
-              make_unique<FullBilloirVertexFitter<ConstantBField,
-                                                  BoundParameters,
-                                                  Propagator<EigenStepper<ConstantBField>>>>(
-                  vertexFitterCfg);
-
-      IterativeVertexFinder<ConstantBField,
-                            BoundParameters,
-                            Propagator<EigenStepper<ConstantBField>>>::Config
-          cfg(bField, std::move(bFitterPtr));
-      IterativeVertexFinder<ConstantBField,
-                            BoundParameters,
-                            Propagator<EigenStepper<ConstantBField>>>::State
-          state;
-
-      IterativeVertexFinder<ConstantBField,
-                            BoundParameters,
-                            Propagator<EigenStepper<ConstantBField>>>
-          finder(cfg);
-
-      // Create perigee surface
-      std::shared_ptr<PerigeeSurface> perigeeSurface
-          = Surface::makeShared<PerigeeSurface>(Vector3D(0., 0., 0.));
-
-      // Create position of vertex and perigee surface
-      double x = vXYDist(gen);
-      double y = vXYDist(gen);
-      double z = vZDist(gen);
-
-      // Calculate d0 and z0 corresponding to vertex position
-      double d0_v = sqrt(x * x + y * y);
-      double z0_v = z;
-
-      // Start constructing nTracks tracks in the following
-      std::vector<BoundParameters> tracks;
-
-      // Construct random track emerging from vicinity of vertex position
-      // Vector to store track objects used for vertex fit
-      for (unsigned int iTrack = 0; iTrack < nTracks; iTrack++) {
-        // Construct positive or negative charge randomly
-        double q = qDist(gen) < 0 ? -1. : 1.;
-
-        // Construct random track parameters
-        TrackParametersBase::ParVector_t paramVec;
-        double                           z0track = z0_v + z0Dist(gen);
-        paramVec << d0_v + d0Dist(gen), z0track, phiDist(gen), thetaDist(gen),
-            q / pTDist(gen);
-
-        // Fill vector of track objects with simple covariance matrix
-        std::unique_ptr<ActsSymMatrixD<5>> covMat
-            = std::make_unique<ActsSymMatrixD<5>>();
-
-        // Resolutions
-        double res_d0 = resIPDist(gen);
-        double res_z0 = resIPDist(gen);
-        double res_ph = resAngDist(gen);
-        double res_th = resAngDist(gen);
-        double res_qp = resQoPDist(gen);
-
-        (*covMat) << res_d0 * res_d0, 0., 0., 0., 0., 0., res_z0 * res_z0, 0.,
-            0., 0., 0., 0., res_ph * res_ph, 0., 0., 0., 0., 0.,
-            res_th * res_th, 0., 0., 0., 0., 0., res_qp * res_qp;
-        tracks.push_back(
-            BoundParameters(std::move(covMat), paramVec, perigeeSurface));
-      }
+      typedef FullBilloirVertexFitter<ConstantBField,
+                                      BoundParameters,
+                                      Propagator<EigenStepper<ConstantBField>>>
+          BilloirFitter;
 
+      // Set up Billoir Vertex Fitter
+      BilloirFitter::Config vertexFitterCfg(bField);
+
+      std::unique_ptr<BilloirFitter> bFitterPtr
+
+          = std::make_unique<BilloirFitter>(vertexFitterCfg);
+
+      // Vertex Finder
+      typedef IterativeVertexFinder<ConstantBField,
+                                    BoundParameters,
+                                    Propagator<EigenStepper<ConstantBField>>>
+                           VertexFinder;
+      VertexFinder::Config cfg(bField, std::move(bFitterPtr));
+      cfg.reassignTracksAfterFirstFit = true;
+      VertexFinder::State  state;
+
+      VertexFinder finder(cfg);
+
+      // Vector to be filled with all tracks in current event
+	    std::vector<BoundParameters> tracks;
+
+	    // Vector to be filled with truth vertices for later comparison
+	    std::vector<Vertex<BoundParameters>> trueVertices;
+
+      // start creating event with nVertices vertices
+      unsigned int nVertices = nVertexDist(gen);
+      for(unsigned int iVertex = 0; iVertex < nVertices; ++iVertex){
+
+      	// Number of tracks
+     	 	unsigned int nTracks = nTracksDist(gen);
+
+     	 	if(debug){
+     	 		std::cout << "Event " << iEvent << ", Vertex " << iVertex << "/" << nVertices 
+     	 		<< " with " << nTracks << " tracks." << std::endl;
+     	 	}
+	      // Create perigee surface
+	      std::shared_ptr<PerigeeSurface> perigeeSurface
+	          = Surface::makeShared<PerigeeSurface>(Vector3D(0., 0., 0.));
+
+	      // Create position of vertex and perigee surface
+	      double x = vXYDist(gen);
+	      double y = vXYDist(gen);
+	      double z = vZDist(gen);
+
+	
+				// True vertex
+				Vertex<BoundParameters> trueV(Vector3D(x,y,z));
+				std::vector<TrackAtVertex<BoundParameters>> tracksAtTrueVtx;
+
+	      // Calculate d0 and z0 corresponding to vertex position
+	      double d0_v = sqrt(x * x + y * y);
+	      double z0_v = z;
+
+	      // Construct random track emerging from vicinity of vertex position
+	      // Vector to store track objects used for vertex fit
+	      for (unsigned int iTrack = 0; iTrack < nTracks; iTrack++) {
+	        // Construct positive or negative charge randomly
+	        double q = qDist(gen) < 0 ? -1. : 1.;
+
+	        // Construct random track parameters
+	        TrackParametersBase::ParVector_t paramVec;
+	        double                           z0track = z0_v + z0Dist(gen);
+	        paramVec << d0_v + d0Dist(gen), z0track, phiDist(gen), thetaDist(gen),
+	            q / pTDist(gen);
+
+	        // Fill vector of track objects with simple covariance matrix
+	        std::unique_ptr<ActsSymMatrixD<5>> covMat
+	            = std::make_unique<ActsSymMatrixD<5>>();
+
+	        // Resolutions
+	        double res_d0 = resIPDist(gen);
+	        double res_z0 = resIPDist(gen);
+	        double res_ph = resAngDist(gen);
+	        double res_th = resAngDist(gen);
+	        double res_qp = resQoPDist(gen);
+
+	        (*covMat) << res_d0 * res_d0, 0., 0., 0., 0., 0., res_z0 * res_z0, 0.,
+	            0., 0., 0., 0., res_ph * res_ph, 0., 0., 0., 0., 0.,
+	            res_th * res_th, 0., 0., 0., 0., 0., res_qp * res_qp;
+	        auto params = BoundParameters(std::move(covMat), paramVec, perigeeSurface);
+	        tracks.push_back(params);
+
+	        TrackAtVertex<BoundParameters> trAtVt(0., params, params);
+	        tracksAtTrueVtx.push_back(trAtVt);
+	      }
+
+	      trueV.setTracksAtVertex(tracksAtTrueVtx);
+	      trueVertices.push_back(trueV);
+
+    } // end loop over vertices
+
+    	// shuffle list of tracks
+    	std::shuffle(std::begin(tracks), std::end(tracks), gen);
+
+    	// find vertices
       VertexingStatus status = finder.find(tracks, state, propagator);
 
+      // check status
       bool statusBool(status == VertexingStatus::SUCCESS);
       BOOST_TEST(statusBool);
 
-      Vector3D result = state.vertexCollection[0].position();
-
-      // CHECK_CLOSE_ABS(result[eZ], z, 1 * units::_mm);
+      // check if same amount of vertices has been found with tolerance of 2
+	    CHECK_CLOSE_ABS(state.vertexCollection.size(), nVertices, 3);
+
+     
+      if(debug){
+      	std::cout << "########## RESULT: ########## Event " << iEvent << std::endl;
+      	std::cout << "Number of true vertices: " << nVertices << std::endl;
+      	std::cout << "Number of reco vertices: " << state.vertexCollection.size() << std::endl;
+
+	      int count = 1;
+	      std::cout << "----- True vertices -----" << std::endl;
+	      for(const auto& vertex : trueVertices){
+	      	Vector3D pos = vertex.position();
+	      	std::cout << count << ". True Vertex:\t Position:" 
+	      	<<  "(" << pos.x()<<","<<pos.y()<<","<<pos.z()<<")" << std::endl;
+	      	std::cout << "Number of tracks: " << vertex.tracks().size() << std::endl << std::endl;
+	      	count ++;
+	      }
+	      std::cout << "----- Reco vertices -----" << std::endl;
+	      count = 1;
+	      for(const auto& vertex : state.vertexCollection){
+	      	Vector3D pos = vertex.position();
+	      	std::cout << count << ". Reco Vertex:\t Position:" 
+	      	<<  "(" << pos.x()<<","<<pos.y()<<","<<pos.z()<<")" << std::endl;
+	      	std::cout << "Number of tracks: " << vertex.tracks().size() << std::endl << std::endl;
+	      	count ++;
+	      }
+	    }
+
+	    // Check if all vertices have been found with close z-values
+	    bool allVerticesFound = true;
+	    for(const auto& trueVertex : trueVertices){
+	    	Vector3D truePos = trueVertex.position();
+	    	bool currentVertexFound = false;
+	    	for(const auto& recoVertex : state.vertexCollection){
+	    		Vector3D recoPos = recoVertex.position();
+	    			// check only for close z distance
+	    			double zDistance = abs(truePos.z() - recoPos.z());
+	    			if(zDistance <  2 * units::_mm){
+	    				currentVertexFound = true;
+	    			}
+	    	}
+	    	if(!currentVertexFound){
+	    		allVerticesFound = false;
+	    	}
+	    }
+
+	    // check if found vertices have compatible z values
+	    BOOST_TEST(allVerticesFound);
+	    
     }
   }