101 check that the vertexing works

CMakeLists.txt

@@ -91,7 +91,7 @@ target_link_libraries(Digitize Struct SANDGeoManager Utils)
 
 # Creates Reconstruct executable.
 add_executable(Reconstruct src/reconstruction.cpp)
-target_link_libraries(Reconstruct PUBLIC Struct Utils SANDGeoManager SANDClustering EDEPTree TrackletFinder SANDTrackerCluster SANDTrackerDigit SANDTrackerUtils SANDKalmanFilter ROOT::EG)
+target_link_libraries(Reconstruct PUBLIC Struct Utils SANDGeoManager SANDClustering EDEPTree TrackletFinder SANDTrackerCluster SANDTrackerDigit SANDTrackerUtils SANDKalmanFilter SANDTrackerVertexing ROOT::EG)
 
 # Creates Analyze executable.
 add_executable(Analyze src/analysis.cpp)
@@ -153,7 +153,7 @@ endif()
 # Copy setup.sh configuration file
 configure_file(setup.sh "${CMAKE_INSTALL_PREFIX}/setup.sh" COPYONLY)
 
- install(TARGETS Utils Struct SANDEventDisplay SANDGeoManager Digitize Reconstruct Display FastCheck eventDisplay TestKF TestSeed TestVertexing SANDClustering SANDECALClustering NoiseMatrixComputation TrackletFinder SANDTrackerCluster SANDTrackerDigit SANDKalmanFilter SANDKFTrack SANDTrackerUtils
+ install(TARGETS Utils Struct SANDEventDisplay SANDGeoManager Digitize Reconstruct Display FastCheck eventDisplay TestKF TestSeed TestVertexing SANDClustering SANDECALClustering NoiseMatrixComputation TrackletFinder SANDTrackerCluster SANDTrackerDigit SANDKalmanFilter SANDKFTrack SANDTrackerUtils SANDTrackerVertexing
         EXPORT SandRecoTargets
         RUNTIME DESTINATION "${CMAKE_INSTALL_PREFIX}/bin"
         LIBRARY DESTINATION "${CMAKE_INSTALL_PREFIX}/lib"

include/SANDTrackerPlane.h

@@ -74,9 +74,14 @@ class Plane
   void addCell(const double transverse_coordinate, Cell c)
   {
     if(coord_to_id_map_.find(transverse_coordinate) == coord_to_id_map_.end()) {
-      c.setPlane(this);
       coord_to_id_map_.insert({transverse_coordinate, c.getId()});
       id_to_cell_map_.insert({c.getId(), c});
+      id_to_cell_map_[c.getId()].setPlane(this);
+    }
+  }
+  void updateCells() {
+    for (auto& c:id_to_cell_map_) {
+      c.second.setPlane(this);
     }
   }
         std::map<CellID, Cell>& getIdToCellMap()       {return id_to_cell_map_;};

include/SANDTrackerVertexing.h

@@ -51,8 +51,11 @@ class Vertex {
 class TrackerVertexing {
 
   public:
-    void setParameters(double dz, double ip, double merging_radius, std::vector<Track> tracks);
+    void setParameters(double dz = 15, double ip = 15, double merging_radius = 15);
+    void setTracks(std::vector<track> tracks);
+    void setTracks(std::vector<Track> tracks);
     int run();
+    const std::vector<Vertex>& getVertices() {return vertices_;};
 
   private:
 
@@ -65,7 +68,7 @@ class TrackerVertexing {
     int mergeVertex();
     void refineVertexPosition(double stepSize = 0.001, int nSteps = 50);
 
-    int id_vert_;
+    int id_vert_ = 0;
     double dz_;
     double ip_;
     double merging_radius_;

include/StructLinkDef.h

@@ -10,6 +10,7 @@
 #pragma link C++ class std::vector < dg_ps> + ;
 #pragma link C++ class std::vector < dg_cell> + ;
 #pragma link C++ class std::vector < reco_cell > +;
+#pragma link C++ class std::vector < vertex > +;
 //#pragma link C++ class std::map < std::string, std::vector < hit>> + ;
 #pragma link C++ class std::vector < dg_wire > +;
 #pragma link C++ class std::vector < track > +;
@@ -27,4 +28,5 @@
 #pragma link C++ class track + ;
 // #pragma link C++ class particle + ;
 #pragma link C++ class event + ;
+#pragma link C++ class vertex + ;
 #endif

include/struct.h

@@ -244,6 +244,15 @@ struct particle
   std::vector<particle> daughters;
 };
 
+struct vertex 
+{
+  int id;
+  double x;
+  double y;
+  double z;
+  std::vector<int> track_ids;
+};
+
 struct event
 {
   double x;

src/SANDDigitizationEDEPSIM.cpp

@@ -3,6 +3,7 @@
 
 #include <iomanip>
 #include <iostream>
+#include <limits>
 
 #include <iomanip>
 
@@ -98,15 +99,21 @@ void CreateDigitsFromHits(const SANDGeoManager& geo,
        it != hits2cell.end(); ++it)  // run over wires
   {
     long did = it->first();  // wire unique id
-    const sand_geometry::tracker::WireInfo& wire_info = geo.getCellInfo(it->first())->second.getWire();
-    double wire_time = 999.;
-    double drift_time = 999.;
-    double signal_time = 999.;
-    double t_hit = 999.;
+    const auto& cell_info = geo.getCellInfo(it->first())->second;
+    const sand_geometry::tracker::WireInfo& wire_info = cell_info.getWire();
+    double wire_time = std::numeric_limits<double>::max();
+    double drift_time = std::numeric_limits<double>::max();
+    double signal_time = std::numeric_limits<double>::max();
+    double t_hit = std::numeric_limits<double>::max();
 
     dg_wire d;
     d.det = it->second[0].det;
     d.did = did;
+    if(!cell_info.getPlane()) {
+      std::cerr << "Plane not found. Skipping hit.." << std::endl;
+      continue;
+    }
+    d.hor = (cell_info.getPlane()->getRotation() == 0) ? true : false;
     d.de = 0;
     // To Do: what point do we want to save? 
     // Center or one of the attachment points?

src/SANDGeoManager.cpp

@@ -1516,8 +1516,13 @@ void SANDGeoManager::rearrangePlanes()
 
   for (auto it = planes_.begin();
             it != planes_.end(); ++it) {
-    id_to_plane_[it->uId()] = it;
-  }
+              id_to_plane_[it->uId()] = it;
+            }
+
+  for (auto it = planes_.begin();
+            it != planes_.end(); ++it) {
+              it->updateCells();
+            }
 }
 
 void SANDGeoManager::setTrackerInfo()

src/SANDTrackerVertexing.cpp

@@ -180,7 +180,7 @@ double pointTo3DlineDistance(double x0, double y0, double z0, Track tr) {
 int TrackerVertexing::mergeVertex() {
   int nVertex = static_cast<int>(vertices_list_.size());
   if (nVertex == 0) {
-    std::cout << "No vertex to merge" << std::endl;
+    // std::cout << "No vertex to merge" << std::endl;
     return 0;
   }
   for (int i = 0; i < nVertex; i++) {
@@ -321,25 +321,48 @@ void TrackerVertexing::flagVertex() {
   }
 }
 
-void TrackerVertexing::setParameters(double dz, double ip, double merging_radius, std::vector<Track> tracks) {
+// vertexing tracks
+void TrackerVertexing::setTracks(std::vector<Track> tracks) {
+  clearAll();
+  tracks_ = tracks;
+}
+
+// struct tracks
+void TrackerVertexing::setTracks(std::vector<track> tracks) {
+  clearAll();
+  for (const auto& t : tracks) {
+    Track converted_track;
+     converted_track.id = t.tid;
+     converted_track.x = t.x0;
+     converted_track.y = t.y0; 
+     converted_track.z = t.z0; 
+     converted_track.tx = t.b;
+
+     double m_yz = (t.y0 - t.yc) / (t.z0 - t.zc);
+     converted_track.ty = m_yz;
+
+    tracks_.push_back(converted_track);
+  }
+}
+
+void TrackerVertexing::setParameters(double dz, double ip, double merging_radius) {
   dz_ = dz;
   ip_ = ip;
   merging_radius_ = merging_radius;
-  tracks_ = tracks;
 
-  std::cout << tracks_.size() << " tracks read." << std::endl;
-  std::cout << "DZ: " << dz_ << "\nIP: " << ip_
-            << "\nMerging Radius: " << merging_radius << std::endl;
+  // std::cout << tracks_.size() << " tracks read." << std::endl;
+  // std::cout << "DZ: " << dz_ << "\nIP: " << ip_
+  //           << "\nMerging Radius: " << merging_radius << std::endl;
 }
 
 int TrackerVertexing::run() {
   int pairVertexes = doVertex();
-  dumpVertex("2-Prong.txt");
-  std::cout << "Start selecting neighboor vertexes (" << pairVertexes
-            << " - 2 Prong )" << std::endl;
+  // dumpVertex("2-Prong.txt");
+  // std::cout << "Start selecting neighboor vertexes (" << pairVertexes
+            // << " - 2 Prong )" << std::endl;
   while (static_cast<int>(vertices_.size()) > 0) selectVertex();
 
-  std::cout << "Start merging..." << std::endl;
+  // std::cout << "Start merging..." << std::endl;
   mergeVertex();
 
   refineVertexPosition();
@@ -353,16 +376,15 @@ int TrackerVertexing::run() {
     vertices_.push_back(vertices_2_prong_.at(j));
   }
 
-  std::cout << "Flag vertexes\n";
+  // std::cout << "Flag vertexes\n";
   flagVertex();
 
-  std::cout << "\n\n======== Results ========" << std::endl;
-  std::cout << "2-Prong: " << vertices_2_prong_.size() << std::endl;
-  std::cout << "Multi-Prong: " << vertices_multi_prong_.size() << std::endl;
-  std::cout << "=========================" << std::endl;
-  std::cout << "\nDump vertexes\n";
-  dumpVertex("vertices.txt");
+  // std::cout << "\n\n======== Results ========" << std::endl;
+  // std::cout << "2-Prong: " << vertices_2_prong_.size() << std::endl;
+  // std::cout << "Multi-Prong: " << vertices_multi_prong_.size() << std::endl;
+  // std::cout << "=========================" << std::endl;
+  // std::cout << "\nDump vertexes\n";
+  // dumpVertex("vertices.txt");
 
-  clearAll();
   return 0;
 }